Changed the estimateEP_Brinkley convergence criteria to more stringent
than the main convergence criteria in equilibrium(). This was necessary to avoid another nonconvergence case. You don't want equilibrate() routine's matrix algorithm to handle ill-conditioned matrices. Changed the storage of element potentials. Now, dimensionless element potentials are storred, instead of dimensional one. There is less dependence on temperature. All debugging print statements now use Cantera's writelog() utility.
This commit is contained in:
Harry Moffat
parent
e433bed576
commit
6d85b768d8
14 changed files with 439 additions and 299 deletions
|
|
@ -13,7 +13,10 @@
|
|||
using namespace Cantera;
|
||||
using namespace std;
|
||||
#ifdef DEBUG_HKM
|
||||
int debug_print_lvl = 0;
|
||||
namespace Cantera {
|
||||
int Cantera::BasisOptimize_print_lvl = 0;
|
||||
static char sbuf[1024];
|
||||
}
|
||||
static void print_stringTrunc(const char *str, int space, int alignment);
|
||||
#endif
|
||||
static int amax(double *x, int j, int n);
|
||||
|
|
@ -111,38 +114,38 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn,
|
|||
|
||||
#ifdef DEBUG_HKM
|
||||
double molSave = 0.0;
|
||||
if (debug_print_lvl >= 1) {
|
||||
printf(" "); for(i=0; i<77; i++) printf("-"); printf("\n");
|
||||
printf(" --- Subroutine BASOPT called to ");
|
||||
printf("calculate the number of components and ");
|
||||
printf("evaluate the formation matrix\n");
|
||||
if (debug_print_lvl > 0) {
|
||||
printf(" ---\n");
|
||||
if (BasisOptimize_print_lvl >= 1) {
|
||||
writelog(" "); for(i=0; i<77; i++) writelog("-"); writelog("\n");
|
||||
writelog(" --- Subroutine BASOPT called to ");
|
||||
writelog("calculate the number of components and ");
|
||||
writelog("evaluate the formation matrix\n");
|
||||
if (BasisOptimize_print_lvl > 0) {
|
||||
writelog(" ---\n");
|
||||
|
||||
printf(" --- Formula Matrix used in BASOPT calculation\n");
|
||||
printf(" --- Species | Order | ");
|
||||
writelog(" --- Formula Matrix used in BASOPT calculation\n");
|
||||
writelog(" --- Species | Order | ");
|
||||
for (j = 0; j < ne; j++) {
|
||||
jj = orderVectorElements[j];
|
||||
printf(" ");
|
||||
writelog(" ");
|
||||
ename = mphase->elementName(jj);
|
||||
print_stringTrunc(ename.c_str(), 4, 1);
|
||||
printf("(%1d)", j);
|
||||
sprintf(sbuf,"(%1d)", j); writelog(sbuf);
|
||||
}
|
||||
printf("\n");
|
||||
writelog("\n");
|
||||
for (k = 0; k < nspecies; k++) {
|
||||
kk = orderVectorSpecies[k];
|
||||
printf(" --- ");
|
||||
writelog(" --- ");
|
||||
sname = mphase->speciesName(kk);
|
||||
print_stringTrunc(sname.c_str(), 11, 1);
|
||||
printf(" | %4d |", k);
|
||||
sprintf(sbuf," | %4d |", k); writelog(sbuf);
|
||||
for (j = 0; j < ne; j++) {
|
||||
jj = orderVectorElements[j];
|
||||
double num = mphase->nAtoms(kk,jj);
|
||||
printf("%6.1g ", num);
|
||||
sprintf(sbuf,"%6.1g ", num); writelog(sbuf);
|
||||
}
|
||||
printf("\n");
|
||||
writelog("\n");
|
||||
}
|
||||
printf(" --- \n");
|
||||
writelog(" --- \n");
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -283,14 +286,16 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn,
|
|||
/* ****************************************** */
|
||||
if (jr != k) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (debug_print_lvl >= 1) {
|
||||
if (BasisOptimize_print_lvl >= 1) {
|
||||
kk = orderVectorSpecies[k];
|
||||
sname = mphase->speciesName(kk);
|
||||
printf(" --- %-12.12s", sname.c_str());
|
||||
sprintf(sbuf," --- %-12.12s", sname.c_str()); writelog(sbuf);
|
||||
jj = orderVectorSpecies[jr];
|
||||
ename = mphase->speciesName(jj);
|
||||
printf("(%9.2g) replaces %-12.12s", molSave, ename.c_str());
|
||||
printf("(%9.2g) as component %3d\n", molNum[jj], jr);
|
||||
sprintf(sbuf,"(%9.2g) replaces %-12.12s", molSave, ename.c_str());
|
||||
writelog(sbuf);
|
||||
sprintf(sbuf,"(%9.2g) as component %3d\n", molNum[jj], jr);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
switch_pos(orderVectorSpecies, jr, k);
|
||||
|
|
@ -363,49 +368,51 @@ int Cantera::BasisOptimize(int *usedZeroedSpecies, bool doFormRxn,
|
|||
*/
|
||||
j = mlequ(DATA_PTR(sm), ne, nComponents, DATA_PTR(formRxnMatrix), nNonComponents);
|
||||
if (j == 1) {
|
||||
printf("ERROR: mlequ returned an error condition\n");
|
||||
writelog("ERROR: mlequ returned an error condition\n");
|
||||
throw CanteraError("basopt", "mlequ returned an error condition");
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM
|
||||
if (debug_print_lvl >= 1) {
|
||||
printf(" ---\n");
|
||||
printf(" --- Number of Components = %d\n", nComponents);
|
||||
printf(" --- Formula Matrix:\n");
|
||||
printf(" --- Components: ");
|
||||
if (Cantera::BasisOptimize_print_lvl >= 1) {
|
||||
writelog(" ---\n");
|
||||
sprintf(sbuf," --- Number of Components = %d\n", nComponents);
|
||||
writelog(sbuf);
|
||||
writelog(" --- Formula Matrix:\n");
|
||||
writelog(" --- Components: ");
|
||||
for (k = 0; k < nComponents; k++) {
|
||||
kk = orderVectorSpecies[k];
|
||||
printf(" %3d (%3d) ", k, kk);
|
||||
sprintf(sbuf," %3d (%3d) ", k, kk); writelog(sbuf);
|
||||
}
|
||||
printf("\n --- Components Moles: ");
|
||||
writelog("\n --- Components Moles: ");
|
||||
for (k = 0; k < nComponents; k++) {
|
||||
kk = orderVectorSpecies[k];
|
||||
printf("%-11.3g", molNumBase[kk]);
|
||||
sprintf(sbuf,"%-11.3g", molNumBase[kk]); writelog(sbuf);
|
||||
}
|
||||
printf("\n --- NonComponent | Moles | ");
|
||||
writelog("\n --- NonComponent | Moles | ");
|
||||
for (i = 0; i < nComponents; i++) {
|
||||
kk = orderVectorSpecies[i];
|
||||
sname = mphase->speciesName(kk);
|
||||
printf("%-11.10s", sname.c_str());
|
||||
sprintf(sbuf,"%-11.10s", sname.c_str()); writelog(sbuf);
|
||||
}
|
||||
printf("\n");
|
||||
writelog("\n");
|
||||
|
||||
for (i = 0; i < nNonComponents; i++) {
|
||||
k = i + nComponents;
|
||||
kk = orderVectorSpecies[k];
|
||||
printf(" --- %3d (%3d) ", k, kk);
|
||||
sprintf(sbuf," --- %3d (%3d) ", k, kk); writelog(sbuf);
|
||||
sname = mphase->speciesName(kk);
|
||||
printf("%-10.10s", sname.c_str());
|
||||
printf("|%10.3g|", molNumBase[kk]);
|
||||
sprintf(sbuf,"%-10.10s", sname.c_str()); writelog(sbuf);
|
||||
sprintf(sbuf,"|%10.3g|", molNumBase[kk]); writelog(sbuf);
|
||||
/*
|
||||
* Print the negative of formRxnMatrix[]; it's easier to interpret.
|
||||
*/
|
||||
for (j = 0; j < nComponents; j++) {
|
||||
printf(" %6.2f", - formRxnMatrix[j + i * ne]);
|
||||
sprintf(sbuf," %6.2f", - formRxnMatrix[j + i * ne]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
printf("\n");
|
||||
writelog("\n");
|
||||
}
|
||||
printf(" "); for (i=0; i<77; i++) printf("-"); printf("\n");
|
||||
writelog(" "); for (i=0; i<77; i++) writelog("-"); writelog("\n");
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -435,7 +442,7 @@ static void print_stringTrunc(const char *str, int space, int alignment)
|
|||
int len = strlen(str);
|
||||
if ((len) >= space) {
|
||||
for (i = 0; i < space; i++) {
|
||||
printf("%c", str[i]);
|
||||
sprintf(sbuf,"%c", str[i]); writelog(sbuf);
|
||||
}
|
||||
} else {
|
||||
if (alignment == 1) {
|
||||
|
|
@ -447,11 +454,11 @@ static void print_stringTrunc(const char *str, int space, int alignment)
|
|||
rs = space - len - ls;
|
||||
}
|
||||
if (ls != 0) {
|
||||
for (i = 0; i < ls; i++) printf(" ");
|
||||
for (i = 0; i < ls; i++) writelog(" ");
|
||||
}
|
||||
printf("%s", str);
|
||||
sprintf(sbuf,"%s", str); writelog(sbuf);
|
||||
if (rs != 0) {
|
||||
for (i = 0; i < rs; i++) printf(" ");
|
||||
for (i = 0; i < rs; i++) writelog(" ");
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -532,7 +539,8 @@ static int amax(double *x, int j, int n) {
|
|||
for (k = i + 1; k < n; ++k) {
|
||||
if (c[k + i * idem] != 0.0) goto FOUND_PIVOT;
|
||||
}
|
||||
printf("vcs_mlequ ERROR: Encountered a zero column: %d\n", i);
|
||||
sprintf(sbuf,"vcs_mlequ ERROR: Encountered a zero column: %d\n", i);
|
||||
writelog(sbuf);
|
||||
return 1;
|
||||
FOUND_PIVOT: ;
|
||||
for (j = 0; j < n; ++j) c[i + j * idem] += c[k + j * idem];
|
||||
|
|
@ -612,11 +620,11 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances,
|
|||
|
||||
double test = -1.0E10;
|
||||
#ifdef DEBUG_HKM
|
||||
if (debug_print_lvl > 0) {
|
||||
printf(" "); for(i=0; i<77; i++) printf("-"); printf("\n");
|
||||
printf(" --- Subroutine ElemRearrange() called to ");
|
||||
printf("check stoich. coefficent matrix\n");
|
||||
printf(" --- and to rearrange the element ordering once\n");
|
||||
if (BasisOptimize_print_lvl > 0) {
|
||||
writelog(" "); for(i=0; i<77; i++) writelog("-"); writelog("\n");
|
||||
writelog(" --- Subroutine ElemRearrange() called to ");
|
||||
writelog("check stoich. coefficent matrix\n");
|
||||
writelog(" --- and to rearrange the element ordering once\n");
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -704,8 +712,9 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances,
|
|||
// We haven't found the number of elements necessary.
|
||||
// This is signalled by returning jr != nComponents.
|
||||
#ifdef DEBUG_HKM
|
||||
if (debug_print_lvl > 0) {
|
||||
printf("Error exit: returning with nComponents = %d\n", jr);
|
||||
if (BasisOptimize_print_lvl > 0) {
|
||||
sprintf(sbuf,"Error exit: returning with nComponents = %d\n", jr);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
return jr;
|
||||
|
|
@ -782,15 +791,16 @@ int Cantera::ElemRearrange(int nComponents, const vector_fp & elementAbundances,
|
|||
/* ****************************************** */
|
||||
if (jr != k) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (debug_print_lvl > 0) {
|
||||
if (BasisOptimize_print_lvl > 0) {
|
||||
kk = orderVectorElements[k];
|
||||
ename = mphase->elementName(kk);
|
||||
printf(" --- "); printf("%-2.2s", ename.c_str());
|
||||
printf("replaces ");
|
||||
writelog(" --- ");
|
||||
sprintf(sbuf,"%-2.2s", ename.c_str()); writelog(sbuf);
|
||||
writelog("replaces ");
|
||||
kk = orderVectorElements[jr];
|
||||
ename = mphase->elementName(kk);
|
||||
printf("%-2.2s", ename.c_str());
|
||||
printf(" as element %3d\n", jr);
|
||||
sprintf(sbuf,"%-2.2s", ename.c_str()); writelog(sbuf);
|
||||
sprintf(sbuf," as element %3d\n", jr); writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
switch_pos(orderVectorElements, jr, k);
|
||||
|
|
|
|||
|
|
@ -30,9 +30,10 @@ using namespace std;
|
|||
#include "stringUtils.h"
|
||||
#include "MultiPhase.h"
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
#ifdef DEBUG_HKM
|
||||
#include "stdio.h"
|
||||
int debug_prnt_lvl = 0;
|
||||
int Cantera::ChemEquil_print_lvl = 0;
|
||||
static char sbuf[1024];
|
||||
#endif
|
||||
#ifndef MIN
|
||||
#define MIN(x,y) (( (x) < (y) ) ? (x) : (y))
|
||||
|
|
@ -65,6 +66,22 @@ namespace Cantera {
|
|||
m_doResPerturb(false)
|
||||
{}
|
||||
|
||||
//! Constructor combined with the initialization function
|
||||
/*!
|
||||
* This constructor initializes the ChemEquil object with everything it
|
||||
* needs to start solving equilibrium problems.
|
||||
* @param s ThermoPhase object that will be used in the equilibrium calls.
|
||||
*/
|
||||
ChemEquil::ChemEquil(thermo_t& s) :
|
||||
m_skip(-1), m_p1(0), m_p2(0),
|
||||
m_elementTotalSum(1.0),
|
||||
m_p0(OneAtm), m_eloc(-1),
|
||||
m_elemFracCutoff(1.0E-100),
|
||||
m_doResPerturb(false)
|
||||
{
|
||||
initialize(s);
|
||||
}
|
||||
|
||||
/// Destructor
|
||||
ChemEquil::~ChemEquil(){
|
||||
if (m_p1) {
|
||||
|
|
@ -251,23 +268,22 @@ namespace Cantera {
|
|||
*/
|
||||
update(s);
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("setInitialMoles: Estimated Mole Fractions\n");
|
||||
printf(" Temperature = %g\n", s.temperature());
|
||||
printf(" Pressure = %g\n", s.pressure());
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog("setInitialMoles: Estimated Mole Fractions\n");
|
||||
sprintf(sbuf," Temperature = %g\n", s.temperature()); writelog(sbuf);
|
||||
sprintf(sbuf," Pressure = %g\n", s.pressure()); writelog(sbuf);
|
||||
for (int k = 0; k < m_kk; k++) {
|
||||
string nnn = s.speciesName(k);
|
||||
double mf = s.moleFraction(k);
|
||||
printf(" %-12s % -10.5g\n", nnn.c_str(), mf);
|
||||
sprintf(sbuf," %-12s % -10.5g\n", nnn.c_str(), mf); writelog(sbuf);
|
||||
}
|
||||
printf(" Element_Name ElementGoal ElementMF\n");
|
||||
writelog(" Element_Name ElementGoal ElementMF\n");
|
||||
for (int m = 0; m < m_mm; m++) {
|
||||
string nnn = s.elementName(m);
|
||||
printf(" %-12s % -10.5g% -10.5g\n",
|
||||
nnn.c_str(), elMoleGoal[m], m_elementmolefracs[m]);
|
||||
sprintf(sbuf," %-12s % -10.5g% -10.5g\n",
|
||||
nnn.c_str(), elMoleGoal[m], m_elementmolefracs[m]); writelog(sbuf);
|
||||
}
|
||||
|
||||
}
|
||||
#endif
|
||||
|
||||
|
|
@ -347,24 +363,24 @@ namespace Cantera {
|
|||
doublereal rrt = 1.0/(GasConstant* s.temperature());
|
||||
scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt);
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
for (m = 0; m < m_nComponents; m++) {
|
||||
int isp = m_component[m];
|
||||
string nnn = s.speciesName(isp);
|
||||
printf("isp = %d, %s\n", isp, nnn.c_str());
|
||||
sprintf(sbuf,"isp = %d, %s\n", isp, nnn.c_str());
|
||||
writelog(sbuf);
|
||||
}
|
||||
double pres = s.pressure();
|
||||
double temp = s.temperature();
|
||||
printf("Pressure = %g\n", pres);
|
||||
printf("Temperature = %g\n", temp);
|
||||
printf(" id Name MF mu/RT \n");
|
||||
|
||||
|
||||
sprintf(sbuf,"Pressure = %g\n", pres); writelog(sbuf);
|
||||
sprintf(sbuf,"Temperature = %g\n", temp); writelog(sbuf);
|
||||
writelog(" id Name MF mu/RT \n");
|
||||
for (n = 0; n < s.nSpecies(); n++) {
|
||||
string nnn = s.speciesName(n);
|
||||
printf("%10d %15s %10.5g %10.5g\n",
|
||||
sprintf(sbuf,"%10d %15s %10.5g %10.5g\n",
|
||||
n, nnn.c_str(), xMF_est[n], mu_RT[n]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -396,9 +412,9 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf(" id CompSpecies ChemPot EstChemPot Diff\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog(" id CompSpecies ChemPot EstChemPot Diff\n");
|
||||
for (m = 0; m < m_nComponents; m++) {
|
||||
int isp = m_component[m];
|
||||
double tmp = 0.0;
|
||||
|
|
@ -406,14 +422,16 @@ namespace Cantera {
|
|||
for (n = 0; n < m_mm; n++) {
|
||||
tmp += nAtoms(isp, n) * lambda_RT[n];
|
||||
}
|
||||
printf("%3d %16s %10.5g %10.5g %10.5g\n",
|
||||
sprintf(sbuf,"%3d %16s %10.5g %10.5g %10.5g\n",
|
||||
m, sname.c_str(), mu_RT[isp], tmp, tmp - mu_RT[isp]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
|
||||
printf(" id ElName Lambda_RT\n");
|
||||
writelog(" id ElName Lambda_RT\n");
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
string ename = s.elementName(m);
|
||||
printf(" %3d %6s %10.5g\n", m, ename.c_str(), lambda_RT[m]);
|
||||
sprintf(sbuf," %3d %6s %10.5g\n", m, ename.c_str(), lambda_RT[m]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -471,7 +489,15 @@ namespace Cantera {
|
|||
vector_fp state;
|
||||
s.saveState(state);
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
/*
|
||||
* Check Compatibility
|
||||
*/
|
||||
if (m_mm != s.nElements() || m_kk != s.nSpecies()) {
|
||||
throw CanteraError("ChemEquil::equilibrate ERROR",
|
||||
"Input ThermoPhase is incompatible with initialization");
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM
|
||||
int n;
|
||||
const vector<string>& eNames = s.elementNames();
|
||||
#endif
|
||||
|
|
@ -654,7 +680,7 @@ namespace Cantera {
|
|||
if (useThermoPhaseElementPotentials) {
|
||||
bool haveEm = s.getElementPotentials(DATA_PTR(x));
|
||||
if (haveEm) {
|
||||
doublereal rt = GasConstant * m_phase->temperature();
|
||||
doublereal rt = GasConstant * s.temperature();
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
x[m] /= rt;
|
||||
}
|
||||
|
|
@ -692,7 +718,7 @@ namespace Cantera {
|
|||
addLogEntry("estimateEP_Brinkley had a singular Jacobian. Continuing anyway");
|
||||
}
|
||||
} else {
|
||||
setToEquilState(s, x, m_phase->temperature());
|
||||
setToEquilState(s, x, s.temperature());
|
||||
// Tempting -> However, nonideal is a problem. Turn on if not worried
|
||||
// about nonideality and you are having problems with the main
|
||||
// algorithm.
|
||||
|
|
@ -706,7 +732,7 @@ namespace Cantera {
|
|||
* Install the log(temp) into the last solution unknown
|
||||
* slot.
|
||||
*/
|
||||
x[m_mm] = log(m_phase->temperature());
|
||||
x[m_mm] = log(s.temperature());
|
||||
|
||||
/*
|
||||
* Setting the max and min values for x[]. Also, if element
|
||||
|
|
@ -720,8 +746,8 @@ namespace Cantera {
|
|||
below[m] = -2000.0;
|
||||
if (elMolesGoal[m] < m_elemFracCutoff && m != m_eloc) x[m] = -1000.0;
|
||||
}
|
||||
above[mm] = log(m_phase->maxTemp() + 1.0);
|
||||
below[mm] = log(m_phase->minTemp() - 1.0);
|
||||
above[mm] = log(s.maxTemp() + 1.0);
|
||||
below[mm] = log(s.minTemp() - 1.0);
|
||||
|
||||
vector_fp grad(nvar, 0.0); // gradient of f = F*F/2
|
||||
vector_fp oldx(nvar, 0.0); // old solution
|
||||
|
|
@ -747,15 +773,15 @@ namespace Cantera {
|
|||
// Compute the Jacobian matrix
|
||||
equilJacobian(s, x, elMolesGoal, jac, xval, yval);
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("Jacobian matrix %d:\n", iter);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf,"Jacobian matrix %d:\n", iter); writelog(sbuf);
|
||||
for (m = 0; m <= m_mm; m++) {
|
||||
printf(" [ ");
|
||||
writelog(" [ ");
|
||||
for (n = 0; n <= m_mm; n++) {
|
||||
printf("%10.5g ", jac(m,n));
|
||||
sprintf(sbuf,"%10.5g ", jac(m,n)); writelog(sbuf);
|
||||
}
|
||||
printf(" ]");
|
||||
writelog(" ]");
|
||||
char xName[32];
|
||||
if (m < m_mm) {
|
||||
string nnn = eNames[m];
|
||||
|
|
@ -769,8 +795,9 @@ namespace Cantera {
|
|||
if (m == m_skip) {
|
||||
sprintf(xName, "x_YY");
|
||||
}
|
||||
printf("%-12s", xName);
|
||||
printf(" = - (%10.5g)\n", res_trial[m]);
|
||||
sprintf(sbuf,"%-12s", xName); writelog(sbuf);
|
||||
sprintf(sbuf, " = - (%10.5g)\n", res_trial[m]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -873,10 +900,10 @@ namespace Cantera {
|
|||
addLogEntry("Relative error in "+m_p2->symbol(),deltay);
|
||||
addLogEntry("Max residual",rmax);
|
||||
beginLogGroup("Element potentials");
|
||||
doublereal rt = GasConstant*m_phase->temperature();
|
||||
doublereal rt = GasConstant* s.temperature();
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
m_lambda[m] = x[m]*rt;
|
||||
addLogEntry("element "+m_phase->elementName(m), fp2str(x[m]));
|
||||
addLogEntry("element "+ s.elementName(m), fp2str(x[m]));
|
||||
}
|
||||
|
||||
if (m_eloc >= 0) {
|
||||
|
|
@ -890,12 +917,12 @@ namespace Cantera {
|
|||
addLogEntry("Saving Element Potentials to ThermoPhase Object");
|
||||
endLogGroup("Element potentials");
|
||||
|
||||
if (m_phase->temperature() > m_phase->maxTemp() + 1.0 ||
|
||||
m_phase->temperature() < m_phase->minTemp() - 1.0 ) {
|
||||
if (s.temperature() > s.maxTemp() + 1.0 ||
|
||||
s.temperature() < s.minTemp() - 1.0 ) {
|
||||
writelog("Warning: Temperature ("
|
||||
+fp2str(m_phase->temperature())+" K) outside "
|
||||
"valid range of "+fp2str(m_phase->minTemp())+" K to "
|
||||
+fp2str(m_phase->maxTemp())+" K\n");
|
||||
+fp2str(s.temperature())+" K) outside "
|
||||
"valid range of "+fp2str(s.minTemp())+" K to "
|
||||
+fp2str(s.maxTemp())+" K\n");
|
||||
}
|
||||
endLogGroup("Converged solution");
|
||||
endLogGroup("ChemEquil::equilibrate");
|
||||
|
|
@ -967,12 +994,13 @@ namespace Cantera {
|
|||
for (m = 0; m < nvar; m++) {
|
||||
x[m] = oldx[m] + damp * step[m];
|
||||
}
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("Solution Unknowns: damp = %g\n", damp);
|
||||
printf(" X_new X_old Step\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf, "Solution Unknowns: damp = %g\n", damp); writelog(sbuf);
|
||||
writelog(" X_new X_old Step\n");
|
||||
for (m = 0; m < nvar; m++) {
|
||||
printf(" %10.5g %10.5g %10.5g\n", x[m], oldx[m], step[m]);
|
||||
sprintf(sbuf," % -10.5g % -10.5g % -10.5g\n", x[m], oldx[m], step[m]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -983,7 +1011,7 @@ namespace Cantera {
|
|||
/**
|
||||
* Evaluates the residual vector F, of length mm
|
||||
*/
|
||||
void ChemEquil::equilResidual(thermo_t& mix, const vector_fp& x,
|
||||
void ChemEquil::equilResidual(thermo_t& s, const vector_fp& x,
|
||||
const vector_fp& elmFracGoal, vector_fp& resid,
|
||||
doublereal xval, doublereal yval)
|
||||
{
|
||||
|
|
@ -991,7 +1019,7 @@ namespace Cantera {
|
|||
int n, m;
|
||||
doublereal xx, yy;
|
||||
doublereal temp = exp(x[m_mm]);
|
||||
setToEquilState(mix, x, temp);
|
||||
setToEquilState(s, x, temp);
|
||||
|
||||
// residuals are the total element moles
|
||||
vector_fp& elmFrac = m_elementmolefracs;
|
||||
|
|
@ -1014,22 +1042,23 @@ namespace Cantera {
|
|||
resid[m] = log( (1.0 + elmFracGoal[m]) / (1.0 + elmFrac[m]) );
|
||||
}
|
||||
}
|
||||
addLogEntry(m_phase->elementName(m),fp2str(elmFrac[m])+" ("
|
||||
addLogEntry(s.elementName(m),fp2str(elmFrac[m])+" ("
|
||||
+fp2str(elmFracGoal[m])+")");
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0 && !m_doResPerturb) {
|
||||
printf("Residual: ElFracGoal ElFracCurrent Resid\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0 && !m_doResPerturb) {
|
||||
writelog("Residual: ElFracGoal ElFracCurrent Resid\n");
|
||||
for (n = 0; n < m_mm; n++) {
|
||||
printf(" % -14.7E % -14.7E % -10.5E\n",
|
||||
sprintf(sbuf," % -14.7E % -14.7E % -10.5E\n",
|
||||
elmFracGoal[n], elmFrac[n], resid[n]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
xx = m_p1->value(mix);
|
||||
yy = m_p2->value(mix);
|
||||
xx = m_p1->value(s);
|
||||
yy = m_p2->value(s);
|
||||
resid[m_mm] = xx/xval - 1.0;
|
||||
resid[m_skip] = yy/yval - 1.0;
|
||||
string xstr = fp2str(xx)+" ("+fp2str(xval)+")";
|
||||
|
|
@ -1038,11 +1067,13 @@ namespace Cantera {
|
|||
addLogEntry(m_p2->symbol(), ystr);
|
||||
endLogGroup("ChemEquil::equilResidual");
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0 && !m_doResPerturb) {
|
||||
printf(" Goal Xvalue Resid\n");
|
||||
printf(" XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]);
|
||||
printf(" YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0 && !m_doResPerturb) {
|
||||
writelog(" Goal Xvalue Resid\n");
|
||||
sprintf(sbuf," XX : % -14.7E % -14.7E % -10.5E\n", xval, xx, resid[m_mm]);
|
||||
writelog(sbuf);
|
||||
sprintf(sbuf," YY(%1d): % -14.7E % -14.7E % -10.5E\n", m_skip, yval, yy, resid[m_skip]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
|
@ -1264,27 +1295,29 @@ namespace Cantera {
|
|||
}
|
||||
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
#ifdef DEBUG_HKM
|
||||
const vector<string>& eNames = s.elementNames();
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("estimateEP_Brinkley::\n\n");
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog("estimateEP_Brinkley::\n\n");
|
||||
double temp = s.temperature();
|
||||
double pres = s.pressure();
|
||||
printf("temp = %g\n", temp);
|
||||
printf("pres = %g\n", pres);
|
||||
printf("Initial mole numbers and mu_SS:\n");
|
||||
printf(" Name MoleNum mu_SS actCoeff\n");
|
||||
sprintf(sbuf, "temp = %g\n", temp); writelog(sbuf);
|
||||
sprintf(sbuf, "pres = %g\n", pres); writelog(sbuf);
|
||||
writelog("Initial mole numbers and mu_SS:\n");
|
||||
writelog(" Name MoleNum mu_SS actCoeff\n");
|
||||
for (k = 0; k < m_kk; k++) {
|
||||
string nnn = s.speciesName(k);
|
||||
printf("%15s %13.5g %13.5g %13.5g\n",
|
||||
sprintf(sbuf,"%15s %13.5g %13.5g %13.5g\n",
|
||||
nnn.c_str(), n_i[k], m_muSS_RT[k], actCoeff[k]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
printf("Initial n_t = %10.5g\n", n_t);
|
||||
printf("Comparison of Goal Element Abundance with Initial Guess:\n");
|
||||
printf(" eName eCurrent eGoal\n");
|
||||
sprintf(sbuf,"Initial n_t = %10.5g\n", n_t); writelog(sbuf);
|
||||
writelog("Comparison of Goal Element Abundance with Initial Guess:\n");
|
||||
writelog(" eName eCurrent eGoal\n");
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
string nnn = s.elementName(m);
|
||||
printf("%5s %13.5g %13.5g\n",nnn.c_str(), eMolesFix[m], elMoles[m]);
|
||||
sprintf(sbuf,"%5s %13.5g %13.5g\n",nnn.c_str(), eMolesFix[m], elMoles[m]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1310,9 +1343,9 @@ namespace Cantera {
|
|||
/*
|
||||
* Calculate the mole numbers of species
|
||||
*/
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("START ITERATION %d:\n", iter);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf, "START ITERATION %d:\n", iter); writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
/*
|
||||
|
|
@ -1324,18 +1357,22 @@ namespace Cantera {
|
|||
Xmol_i_calc[k] = n_i_calc[k]/n_t_calc;
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf(" Species: Calculated_Moles Calculated_Mole_Fraction\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog(" Species: Calculated_Moles Calculated_Mole_Fraction\n");
|
||||
for (k = 0; k < m_kk; k++) {
|
||||
string nnn = s.speciesName(k);
|
||||
printf("%15s: %10.5g %10.5g\n", nnn.c_str(), n_i_calc[k], Xmol_i_calc[k]);
|
||||
sprintf(sbuf,"%15s: %10.5g %10.5g\n", nnn.c_str(), n_i_calc[k], Xmol_i_calc[k]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
printf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc);
|
||||
printf("(iter %d) element moles bal: Goal Calculated\n", iter);
|
||||
sprintf(sbuf,"%15s: %10.5g\n", "Total Molar Sum", n_t_calc);
|
||||
writelog(sbuf);
|
||||
sprintf(sbuf,"(iter %d) element moles bal: Goal Calculated\n", iter);
|
||||
writelog(sbuf);
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
string nnn = eNames[m];
|
||||
printf(" %8s: %10.5g %10.5g \n", nnn.c_str(), elMoles[m], eMolesCalc[m]);
|
||||
sprintf(sbuf," %8s: %10.5g %10.5g \n", nnn.c_str(), elMoles[m], eMolesCalc[m]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1359,10 +1396,11 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
if (!normalStep) {
|
||||
printf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM);
|
||||
sprintf(sbuf," NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1427,9 +1465,9 @@ namespace Cantera {
|
|||
}
|
||||
|
||||
nCutoff = 1.0E-9 * n_t_calc;
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf(" Lump Sum Elements Calculation: \n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog(" Lump Sum Elements Calculation: \n");
|
||||
}
|
||||
#endif
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
|
|
@ -1455,11 +1493,11 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
string nnn = eNames[m];
|
||||
printf(" %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2);
|
||||
|
||||
sprintf(sbuf," %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
|
|
@ -1514,9 +1552,10 @@ namespace Cantera {
|
|||
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
if (a1(m,m) < 1.0E-50) {
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf(" NOTE: Diagonalizing the analytical Jac row %d\n", m);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf," NOTE: Diagonalizing the analytical Jac row %d\n", m);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
for (n = 0; n < m_mm; n++) {
|
||||
|
|
@ -1536,32 +1575,39 @@ namespace Cantera {
|
|||
|
||||
resid[m_mm] = n_t - n_t_calc;
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("Matrix:\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog("Matrix:\n");
|
||||
for (m = 0; m <= m_mm; m++) {
|
||||
printf(" [");
|
||||
writelog(" [");
|
||||
for (n = 0; n <= m_mm; n++) {
|
||||
printf(" %10.5g", a1(m,n));
|
||||
sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf);
|
||||
}
|
||||
printf("] = %10.5g\n", resid[m]);
|
||||
sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
||||
tmp = resid[m_mm] /(n_t + 1.0E-15);
|
||||
sum += tmp * tmp;
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("(it %d) Convergence = %g\n", iter, sum);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf,"(it %d) Convergence = %g\n", iter, sum);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
if (sum < 100. * options.relTolerance) {
|
||||
/*
|
||||
* Insist on 20x accuracy compared to the top routine.
|
||||
* There are instances, for ill-conditioned or
|
||||
* singular matrices where this is needed to move
|
||||
* the system to a point where the matrices aren't
|
||||
* singular.
|
||||
*/
|
||||
if (sum < 0.05 * options.relTolerance) {
|
||||
retn = 0;
|
||||
goto exit;
|
||||
}
|
||||
|
||||
|
||||
/*
|
||||
* Row Sum scaling
|
||||
*/
|
||||
|
|
@ -1571,9 +1617,10 @@ namespace Cantera {
|
|||
tmp += fabs(a1(m,n));
|
||||
}
|
||||
if (m < m_mm && tmp < 1.0E-30) {
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf(" NOTE: Diagonalizing row %d\n", m);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf," NOTE: Diagonalizing row %d\n", m);
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
for (n = 0; n <= m_mm; n++) {
|
||||
|
|
@ -1590,15 +1637,15 @@ namespace Cantera {
|
|||
resid[m] *= tmp;
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("Row Summed Matrix:\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog("Row Summed Matrix:\n");
|
||||
for (m = 0; m <= m_mm; m++) {
|
||||
printf(" [");
|
||||
writelog(" [");
|
||||
for (n = 0; n <= m_mm; n++) {
|
||||
printf(" %10.5g", a1(m,n));
|
||||
sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf);
|
||||
}
|
||||
printf("] = %10.5g\n", resid[m]);
|
||||
sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1640,12 +1687,14 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
if (sameAsRow >= 0 || lumpSum[m]) {
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
if (lumpSum[m]) {
|
||||
printf("Lump summing row %d, due to rank deficiency analysis\n", m);
|
||||
sprintf(sbuf,"Lump summing row %d, due to rank deficiency analysis\n", m);
|
||||
writelog(sbuf);
|
||||
} else if (sameAsRow >= 0) {
|
||||
printf("Identified that rows %d and %d are the same\n", m, sameAsRow);
|
||||
sprintf(sbuf,"Identified that rows %d and %d are the same\n", m, sameAsRow);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1659,15 +1708,15 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0 && modifiedMatrix) {
|
||||
printf("Row Summed, MODIFIED Matrix:\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0 && modifiedMatrix) {
|
||||
writelog("Row Summed, MODIFIED Matrix:\n");
|
||||
for (m = 0; m <= m_mm; m++) {
|
||||
printf(" [");
|
||||
writelog(" [");
|
||||
for (n = 0; n <= m_mm; n++) {
|
||||
printf(" %10.5g", a1(m,n));
|
||||
sprintf(sbuf," %10.5g", a1(m,n)); writelog(sbuf);
|
||||
}
|
||||
printf("] = %10.5g\n", resid[m]);
|
||||
sprintf(sbuf,"] = %10.5g\n", resid[m]); writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1677,9 +1726,9 @@ namespace Cantera {
|
|||
}
|
||||
catch (CanteraError) {
|
||||
addLogEntry("estimateEP_Brinkley:Jacobian is singular.");
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("Matrix is SINGULAR.ERROR\n");
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
writelog("Matrix is SINGULAR.ERROR\n");
|
||||
}
|
||||
#endif
|
||||
s.restoreState(state);
|
||||
|
|
@ -1710,10 +1759,10 @@ namespace Cantera {
|
|||
}
|
||||
}
|
||||
}
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
if (beta != 1.0) {
|
||||
printf("(it %d) Beta = %g\n", iter, beta);
|
||||
sprintf(sbuf,"(it %d) Beta = %g\n", iter, beta); writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1728,29 +1777,34 @@ namespace Cantera {
|
|||
n_t *= exp(beta * resid[m_mm]);
|
||||
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
printf("(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
sprintf(sbuf,"(it %d) OLD_SOLUTION NEW SOLUTION (undamped updated)\n", iter);
|
||||
writelog(sbuf);
|
||||
for (m = 0; m < m_mm; m++) {
|
||||
string eee = eNames[m];
|
||||
printf(" %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]);
|
||||
sprintf(sbuf," %5s %10.5g %10.5g %10.5g\n", eee.c_str(), x_old[m], x[m], resid[m]);
|
||||
writelog(sbuf);
|
||||
}
|
||||
printf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm]));
|
||||
sprintf(sbuf," n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm]));
|
||||
writelog(sbuf);
|
||||
}
|
||||
#endif
|
||||
}
|
||||
exit:
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
double temp = s.temperature();
|
||||
double pres = s.pressure();
|
||||
|
||||
if (retn == 0) {
|
||||
printf(" ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n",
|
||||
sprintf(sbuf," ChemEquil::estimateEP_Brinkley() SUCCESS: equilibrium found at T = %g, Pres = %g\n",
|
||||
temp, pres);
|
||||
writelog(sbuf);
|
||||
} else {
|
||||
printf(" ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n",
|
||||
sprintf(sbuf," ChemEquil::estimateEP_Brinkley() FAILURE: equilibrium not found at T = %g, Pres = %g\n",
|
||||
temp, pres);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1767,12 +1821,12 @@ namespace Cantera {
|
|||
s.getMoleFractions(DATA_PTR(m_molefractions));
|
||||
int k;
|
||||
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
#ifdef DEBUG_HKM
|
||||
int maxPosEloc = -1;
|
||||
int maxNegEloc = -1;
|
||||
double maxPosVal = -1.0;
|
||||
double maxNegVal = -1.0;
|
||||
if (debug_prnt_lvl > 0) {
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
for (k = 0; k < m_kk; k++) {
|
||||
if (nAtoms(k,m_eloc) > 0.0) {
|
||||
if (m_molefractions[k] > maxPosVal && m_molefractions[k] > 0.0) {
|
||||
|
|
@ -1805,13 +1859,14 @@ namespace Cantera {
|
|||
if (sumPos >= sumNeg) {
|
||||
if ( sumPos <= 0.0) return;
|
||||
double factor = (elMolesGoal[m_eloc] + sumNeg) / sumPos;
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
if (factor < 0.9999999999) {
|
||||
string nnn = s.speciesName(maxPosEloc);
|
||||
printf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
|
||||
sprintf(sbuf,"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
|
||||
nnn.c_str(),
|
||||
m_molefractions[maxPosEloc], m_molefractions[maxPosEloc]*factor);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
@ -1822,13 +1877,14 @@ namespace Cantera {
|
|||
}
|
||||
} else {
|
||||
double factor = (-elMolesGoal[m_eloc] + sumPos) / sumNeg;
|
||||
#ifdef DEBUG_HKM_EPEQUIL
|
||||
if (debug_prnt_lvl > 0) {
|
||||
#ifdef DEBUG_HKM
|
||||
if (ChemEquil_print_lvl > 0) {
|
||||
if (factor < 0.9999999999) {
|
||||
string nnn = s.speciesName(maxNegEloc);
|
||||
printf("adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
|
||||
sprintf(sbuf,"adjustEloc: adjusted %s and friends from %g to %g to ensure neutrality condition\n",
|
||||
nnn.c_str(),
|
||||
m_molefractions[maxNegEloc], m_molefractions[maxNegEloc]*factor);
|
||||
writelog(sbuf);
|
||||
}
|
||||
}
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -98,7 +98,17 @@ namespace Cantera {
|
|||
class ChemEquil {
|
||||
|
||||
public:
|
||||
//! Default Constructor
|
||||
ChemEquil();
|
||||
|
||||
//! Constructor combined with the initialization function
|
||||
/*!
|
||||
* This constructor initializes the ChemEquil object with everything it
|
||||
* needs to start solving equilibrium problems.
|
||||
* @param s ThermoPhase object that will be used in the equilibrium calls.
|
||||
*/
|
||||
ChemEquil(thermo_t& s);
|
||||
|
||||
virtual ~ChemEquil();
|
||||
|
||||
int equilibrate(thermo_t& s, const char* XY,
|
||||
|
|
@ -116,6 +126,14 @@ namespace Cantera {
|
|||
|
||||
protected:
|
||||
|
||||
//! Pointer to the %ThermoPhase object used to initialize this object.
|
||||
|
||||
/*!
|
||||
* This %ThermoPhase object must be compatible with the %ThermoPhase
|
||||
* objects input from the equilibrate function. Currently, this
|
||||
* means that the 2 %ThermoPhases have to have consist of the same
|
||||
* species and elements.
|
||||
*/
|
||||
thermo_t* m_phase;
|
||||
|
||||
/// number of atoms of element m in species k.
|
||||
|
|
@ -226,9 +244,14 @@ namespace Cantera {
|
|||
|
||||
vector_int m_orderVectorElements;
|
||||
vector_int m_orderVectorSpecies;
|
||||
|
||||
|
||||
};
|
||||
|
||||
#ifdef DEBUG_HKM
|
||||
extern int ChemEquil_print_lvl;
|
||||
#endif
|
||||
|
||||
}
|
||||
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -308,6 +308,9 @@ namespace Cantera {
|
|||
MultiPhase *mphase,
|
||||
vector_int & orderVectorSpecies,
|
||||
vector_int & orderVectorElements);
|
||||
#ifdef DEBUG_HKM
|
||||
extern int BasisOptimize_print_lvl;
|
||||
#endif
|
||||
}
|
||||
|
||||
#endif
|
||||
|
|
|
|||
|
|
@ -77,7 +77,7 @@ namespace Cantera {
|
|||
|
||||
m_index = right.m_index;
|
||||
m_phi = right.m_phi;
|
||||
m_lambda = right.m_lambda;
|
||||
m_lambdaRRT = right.m_lambdaRRT;
|
||||
m_hasElementPotentials = right.m_hasElementPotentials;
|
||||
|
||||
return *this;
|
||||
|
|
@ -379,31 +379,75 @@ namespace Cantera {
|
|||
|
||||
}
|
||||
|
||||
/**
|
||||
* Set the thermodynamic state.
|
||||
*/
|
||||
/**
|
||||
* Set the thermodynamic state.
|
||||
*/
|
||||
void ThermoPhase::setStateFromXML(const XML_Node& state) {
|
||||
|
||||
string comp = getString(state,"moleFractions");
|
||||
if (comp != "")
|
||||
setMoleFractionsByName(comp);
|
||||
else {
|
||||
comp = getString(state,"massFractions");
|
||||
if (comp != "")
|
||||
setMassFractionsByName(comp);
|
||||
}
|
||||
if (state.hasChild("temperature")) {
|
||||
double t = getFloat(state, "temperature", "temperature");
|
||||
setTemperature(t);
|
||||
}
|
||||
if (state.hasChild("pressure")) {
|
||||
double p = getFloat(state, "pressure", "pressure");
|
||||
setPressure(p);
|
||||
}
|
||||
if (state.hasChild("density")) {
|
||||
double rho = getFloat(state, "density", "density");
|
||||
setDensity(rho);
|
||||
}
|
||||
string comp = getString(state,"moleFractions");
|
||||
if (comp != "")
|
||||
setMoleFractionsByName(comp);
|
||||
else {
|
||||
comp = getString(state,"massFractions");
|
||||
if (comp != "")
|
||||
setMassFractionsByName(comp);
|
||||
}
|
||||
if (state.hasChild("temperature")) {
|
||||
double t = getFloat(state, "temperature", "temperature");
|
||||
setTemperature(t);
|
||||
}
|
||||
if (state.hasChild("pressure")) {
|
||||
double p = getFloat(state, "pressure", "pressure");
|
||||
setPressure(p);
|
||||
}
|
||||
if (state.hasChild("density")) {
|
||||
double rho = getFloat(state, "density", "density");
|
||||
setDensity(rho);
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
/*
|
||||
* Called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
* the element potentials to this object after every successful
|
||||
* equilibration routine.
|
||||
* The element potentials are storred in their dimensionless
|
||||
* forms, calculated by dividing by RT.
|
||||
* @param lambda vector containing the element potentials.
|
||||
* Length = nElements. Units are Joules/kmol.
|
||||
*/
|
||||
void ThermoPhase::setElementPotentials(const vector_fp& lambda) {
|
||||
doublereal rrt = 1.0/(GasConstant* temperature());
|
||||
int mm = nElements();
|
||||
if (lambda.size() < (size_t) mm) {
|
||||
throw CanteraError("setElementPotentials", "lambda too small");
|
||||
}
|
||||
if (!m_hasElementPotentials) {
|
||||
m_lambdaRRT.resize(mm);
|
||||
}
|
||||
for (int m = 0; m < mm; m++) {
|
||||
m_lambdaRRT[m] = lambda[m] * rrt;
|
||||
}
|
||||
m_hasElementPotentials = true;
|
||||
}
|
||||
|
||||
/*
|
||||
* Returns the storred element potentials.
|
||||
* The element potentials are retrieved from their storred
|
||||
* dimensionless forms by multiplying by RT.
|
||||
* @param lambda Vector containing the element potentials.
|
||||
* Length = nElements. Units are Joules/kmol.
|
||||
*/
|
||||
bool ThermoPhase::getElementPotentials(doublereal* lambda) const {
|
||||
doublereal rt = GasConstant* temperature();
|
||||
int mm = nElements();
|
||||
if (m_hasElementPotentials) {
|
||||
for (int m = 0; m < mm; m++) {
|
||||
lambda[m] = m_lambdaRRT[m] * rt;
|
||||
}
|
||||
}
|
||||
return (m_hasElementPotentials);
|
||||
}
|
||||
|
||||
}
|
||||
|
|
|
|||
|
|
@ -751,19 +751,32 @@ namespace Cantera {
|
|||
err("setToEquilState");
|
||||
}
|
||||
|
||||
// Called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
// the element potentials to this object after every successful
|
||||
// equilibration routine.
|
||||
void setElementPotentials(const vector_fp& lambda) {
|
||||
m_lambda = lambda;
|
||||
m_hasElementPotentials = true;
|
||||
}
|
||||
//! Stores the element potentials in the ThermoPhase object
|
||||
/*!
|
||||
* Called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
* the element potentials to this object after every successful
|
||||
* equilibration routine.
|
||||
* The element potentials are storred in their dimensionless
|
||||
* forms, calculated by dividing by RT.
|
||||
* @param lambda Input vector containing the element potentials.
|
||||
* Length = nElements. Units are Joules/kmol.
|
||||
*/
|
||||
void setElementPotentials(const vector_fp& lambda);
|
||||
|
||||
|
||||
bool getElementPotentials(doublereal* lambda) {
|
||||
if (m_hasElementPotentials)
|
||||
std::copy(m_lambda.begin(), m_lambda.end(), lambda);
|
||||
return (m_hasElementPotentials);
|
||||
}
|
||||
//! Returns the element potentials storred in the ThermoPhase object
|
||||
/*!
|
||||
* Returns the storred element potentials.
|
||||
* The element potentials are retrieved from their storred
|
||||
* dimensionless forms by multiplying by RT.
|
||||
* @param lambda Output vector containing the element potentials.
|
||||
* Length = nElements. Units are Joules/kmol.
|
||||
* @return bool indicating whether thare are any valid storred element
|
||||
* potentials. The calling routine should check this
|
||||
* bool. In the case that there aren't any, lambda is not
|
||||
* touched.
|
||||
*/
|
||||
bool getElementPotentials(doublereal* lambda) const;
|
||||
|
||||
//@}
|
||||
|
||||
|
|
@ -1017,7 +1030,7 @@ namespace Cantera {
|
|||
doublereal m_phi;
|
||||
/// Vector of element potentials.
|
||||
/// -> length equal to number of elements
|
||||
vector_fp m_lambda;
|
||||
vector_fp m_lambdaRRT;
|
||||
bool m_hasElementPotentials;
|
||||
|
||||
private:
|
||||
|
|
|
|||
|
|
@ -1,14 +1,13 @@
|
|||
/**
|
||||
* @file sort.cpp
|
||||
*
|
||||
* $Id$
|
||||
*/
|
||||
|
||||
|
||||
#ifdef WIN32
|
||||
#pragma warning(disable:4786)
|
||||
#endif
|
||||
|
||||
|
||||
#include "sort.h"
|
||||
|
||||
namespace Cantera {
|
||||
|
|
@ -17,6 +16,7 @@ namespace Cantera {
|
|||
|
||||
void heapsort(vector_fp& x, vector_int& y) {
|
||||
int n = x.size();
|
||||
if (n < 2) return;
|
||||
doublereal rra;
|
||||
integer rrb;
|
||||
int ll = n/2;
|
||||
|
|
@ -66,6 +66,7 @@ namespace Cantera {
|
|||
|
||||
void heapsort(vector_fp& x, vector_fp& y) {
|
||||
int n = x.size();
|
||||
if (n < 2) return;
|
||||
doublereal rra;
|
||||
doublereal rrb;
|
||||
int ll = n/2;
|
||||
|
|
@ -115,8 +116,3 @@ namespace Cantera {
|
|||
|
||||
}
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -25,8 +25,21 @@ namespace Cantera {
|
|||
*/
|
||||
std::string fp2str(double x, std::string fmt) {
|
||||
char buf[30];
|
||||
sprintf(buf, fmt.c_str(), x);
|
||||
return std::string(buf);
|
||||
int n = snprintf(buf, 30, fmt.c_str(), x);
|
||||
if (n > 0) {
|
||||
buf[29] = '\0';
|
||||
return std::string(buf);
|
||||
}
|
||||
return std::string(" ");
|
||||
}
|
||||
std::string fp2str(double x) {
|
||||
char buf[30];
|
||||
int n = snprintf(buf, 30, "%g" , x);
|
||||
if (n > 0) {
|
||||
buf[29] = '\0';
|
||||
return std::string(buf);
|
||||
}
|
||||
return std::string(" ");
|
||||
}
|
||||
|
||||
/**
|
||||
|
|
@ -34,8 +47,24 @@ namespace Cantera {
|
|||
*/
|
||||
std::string int2str(int n, std::string fmt) {
|
||||
char buf[30];
|
||||
int m = snprintf(buf, 30, fmt.c_str(), n);
|
||||
sprintf(buf, fmt.c_str(), n);
|
||||
return std::string(buf);
|
||||
if (m > 0) {
|
||||
buf[29] = '\0';
|
||||
return std::string(buf);
|
||||
}
|
||||
return std::string(" ");
|
||||
}
|
||||
|
||||
|
||||
std::string int2str(int n) {
|
||||
char buf[30];
|
||||
int m = snprintf(buf, 30, "%d", n);
|
||||
if (m > 0) {
|
||||
buf[29] = '\0';
|
||||
return std::string(buf);
|
||||
}
|
||||
return std::string(" ");
|
||||
}
|
||||
|
||||
std::string lowercase(std::string s) {
|
||||
|
|
|
|||
|
|
@ -17,8 +17,10 @@ namespace Cantera {
|
|||
class Phase;
|
||||
class ThermoPhase;
|
||||
|
||||
std::string fp2str(double x, std::string fmt = "%g");
|
||||
std::string int2str(int n, std::string fmt = "%d");
|
||||
std::string fp2str(double x, std::string fmt);
|
||||
std::string fp2str(double x);
|
||||
std::string int2str(int n, std::string fmt);
|
||||
std::string int2str(int n);
|
||||
std::string stripws(std::string s);
|
||||
std::string stripnonprint(std::string s);
|
||||
std::string lowercase(std::string s);
|
||||
|
|
|
|||
|
|
@ -677,15 +677,6 @@ namespace Cantera {
|
|||
err("setToEquilState");
|
||||
}
|
||||
|
||||
// called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
// the element potentials to this object
|
||||
void setElementPotentials(const vector_fp& lambda) {
|
||||
m_lambda = lambda;
|
||||
}
|
||||
|
||||
void getElementPotentials(doublereal* lambda) {
|
||||
copy(m_lambda.begin(), m_lambda.end(), lambda);
|
||||
}
|
||||
|
||||
//@}
|
||||
|
||||
|
|
|
|||
|
|
@ -684,16 +684,6 @@ namespace Cantera {
|
|||
err("setToEquilState");
|
||||
}
|
||||
|
||||
// called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
// the element potentials to this object
|
||||
void setElementPotentials(const vector_fp& lambda) {
|
||||
m_lambda = lambda;
|
||||
}
|
||||
|
||||
void getElementPotentials(doublereal* lambda) {
|
||||
copy(m_lambda.begin(), m_lambda.end(), lambda);
|
||||
}
|
||||
|
||||
//@}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -582,16 +582,6 @@ namespace Cantera {
|
|||
err("setToEquilState");
|
||||
}
|
||||
|
||||
// called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
// the element potentials to this object
|
||||
void setElementPotentials(const vector_fp& lambda) {
|
||||
m_lambda = lambda;
|
||||
}
|
||||
|
||||
void getElementPotentials(doublereal* lambda) {
|
||||
copy(m_lambda.begin(), m_lambda.end(), lambda);
|
||||
}
|
||||
|
||||
//@}
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -365,15 +365,7 @@ namespace Cantera {
|
|||
err("setToEquilState");
|
||||
}
|
||||
|
||||
// called by function 'equilibrate' in ChemEquil.h to transfer
|
||||
// the element potentials to this object
|
||||
void setElementPotentials(const vector_fp& lambda) {
|
||||
m_lambda = lambda;
|
||||
}
|
||||
|
||||
void getElementPotentials(doublereal* lambda) {
|
||||
copy(m_lambda.begin(), m_lambda.end(), lambda);
|
||||
}
|
||||
|
||||
//@}
|
||||
|
||||
|
|
|
|||
|
|
@ -1,3 +1,4 @@
|
|||
Makefile
|
||||
SunWS_cache
|
||||
.depends
|
||||
*.d
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue