research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

enhanced and cleaned up equilibrium code

Browse source
This commit is contained in:
Dave Goodwin 2005-07-22 03:27:12 +00:00
parent 2f1d6f3179
commit 44a4226481
8 changed files with 508 additions and 413 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

362
Cantera/src/ChemEquil.cpp
View file

@ -2,8 +2,8 @@
*
* @file ChemEquil.cpp
*
* Chemical equilibrium.
* Implementation file for class ChemEquil
* Chemical equilibrium. Implementation file for class
* ChemEquil.
*
* Copyright 2001 California Institute of Technology
*
@ -19,16 +19,17 @@ using namespace std;
#include "ChemEquil.h"
#include "DenseMatrix.h"
#include "recipes.h"
#include "sort.h"
#include "PropertyCalculator.h"
#include "ctexceptions.h"
#include "vec_functions.h"
#include "stringUtils.h"
#include "MultiPhase.h"
namespace Cantera {
/// map property strings to integers
int _equilflag(const char* xy) {
string flag = string(xy);
if (flag == "TP") return TP;
@ -61,13 +62,13 @@ namespace Cantera {
/**
* Prepare for equilibrium calculations with a specified
* mixture.
* @param s mixture
* Prepare for equilibrium calculations.
* @param s object representing the solution phase.
*/
void ChemEquil::initialize(thermo_t& s)
{
// store a pointer to s and some of its properties locally
// store a pointer to s and some of its properties locally.
// Note: the use of two pointers is a historical artifact.
m_thermo = &s;
m_phase = &s;
@ -89,6 +90,7 @@ namespace Cantera {
m_startSoln.resize(m_mm+1);
m_grt.resize(m_kk);
m_mu_RT.resize(m_kk);
m_component.resize(m_mm,-2);
// set up elemental composition matrix
int m, k, mneg = -1;
@ -97,8 +99,9 @@ namespace Cantera {
for (k = 0; k < m_kk; k++) {
na = m_phase->nAtoms(k,m);
// handle the case of negative atom numbers (used to
// represent positive ions)
// handle the case of negative atom numbers (used to
// represent positive ions, where the 'element' is an
// electron
if (na < 0.0) {
// if negative atom numbers have already been specified
@ -145,7 +148,7 @@ namespace Cantera {
void ChemEquil::setToEquilState(thermo_t& s,
const vector_fp& lambda_RT, doublereal t)
{
// compute the chemical potentials by summing element potentials
// construct the chemical potentials by summing element potentials
fill(m_mu_RT.begin(), m_mu_RT.end(), 0.0);
for (int k = 0; k < m_kk; k++)
for (int m = 0; m < m_mm; m++)
@ -193,75 +196,43 @@ namespace Cantera {
for (m = 0; m < m_mm; m++) m_elementmolefracs[m] /= sum;
}
/**
* Estimate the initial mole fractions. Uses the Simplex method
* to estimate the initial number of moles of each species. The
* linear Gibbs minimization problem is solved, neglecting the
* free energy of mixing terms. This procedure produces a good
* estimate of the low-temperature equilibrium composition.
*
* @param s phase object
* @param elementMoles vector of elemental moles
*/
int ChemEquil::setInitialMoles(thermo_t& s,
vector_fp& elementMoles)
{
int m, n;
double pres = s.pressure();
double lp = log(pres/m_p0);
integer mm = m_phase->nElements();
integer kksp = m_phase->nSpecies();
DenseMatrix aa(mm+2, kksp+1, 0.0);
// first column contains fixed element moles
for (m = 0; m < mm; m++) {
aa(m+1,0) = elementMoles[m]; // + 0.01;
}
/// Estimate the initial mole numbers. This version borrows from the
/// MultiPhaseEquil solver.
int ChemEquil::setInitialMoles(thermo_t& s) {
MultiPhase* mp = 0;
MultiPhaseEquil* e = 0;
int iok = 0;
beginLogGroup("ChemEquil::setInitialMoles");
try {
mp = new MultiPhase;
mp->addPhase(&s, 1.0);
mp->init();
e = new MultiPhaseEquil(mp, true);
e->setInitialMixMoles();
// get the array of non-dimensional Gibbs functions for the pure
// species
s.getGibbs_RT(m_grt.begin());
int kpp = 0;
for (int k = 0; k < kksp; k++) {
kpp++;
aa(0, kpp) = -m_grt[k];
aa(0, kpp) -= lp; // ideal gas
for (int q = 0; q < mm; q++)
aa(q+1, kpp) = -nAtoms(k, q);
}
integer mp = mm+2; // parameters for SIMPLX
integer np = kksp+1;
integer m1 = 0;
integer m2 = 0;
integer m3 = mm;
integer icase=0;
vector_int iposv(mm);
vector_int izrov(kksp);
// solve the linear programming problem
simplx_(&aa(0,0), &mm, &kksp, &mp, &np, &m1, &m2, &m3,
&icase, izrov.begin(), iposv.begin());
fill(m_molefractions.begin(), m_molefractions.end(), 0.0);
for (n = 0; n < mm; n++) {
int ksp = 0;
int ip = iposv[n] - 1;
for (int k = 0; k < kksp; k++) {
if (ip == ksp) {
m_molefractions[k] = aa(n+1, 0);
}
ksp++;
// store component indices
for (int m = 0; m < m_mm; m++) {
m_component[m] = e->componentIndex(m);
}
for (int k = 0; k < m_kk; k++) {
if (m_phase->moleFraction(k) > 0.0) {
addLogEntry(m_phase->speciesName(k),
m_phase->moleFraction(k));
}
}
update(s);
delete e;
delete mp;
iok = 0;
}
s.setState_PX(pres, m_molefractions.begin());
update(s);
return icase;
catch (CanteraError) {
delete e;
delete mp;
iok = -1;
}
endLogGroup();
return iok;
}
@ -271,60 +242,29 @@ namespace Cantera {
int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda)
{
int k, ksp, m, n;
for (k = 0; k < m_kk; k++) {
if (m_molefractions[k] > 0.0) {
m_molefractions[k] = fmaxx(m_molefractions[k], 0.05);
}
//else
// m_molefractions[k] = 0.001;
}
s.setState_PX(s.pressure(), m_molefractions.begin());
beginLogGroup("estimateElementPotentials");
//for (k = 0; k < m_kk; k++) {
// if (m_molefractions[k] > 0.0) {
// m_molefractions[k] = fmaxx(m_molefractions[k], 0.05);
// }
//}
//s.setState_PX(s.pressure(), m_molefractions.begin());
// sort mole fractions
vector_fp mol(m_kk, 0.0);
vector_int index(m_kk, 0);
for (k = 0; k < m_kk; k++) {
mol[k] = m_molefractions[k];
index[k] = k;
}
heapsort(mol, index);
DenseMatrix aa(m_mm, m_mm, 0.0);
vector_fp b(m_mm, -999.0);
vector_fp ipvt(m_mm, 0);
// find a set of constituents
vector_int kc(m_mm, 0);
vector_fp tmp(m_mm, 0.0);
vector_fp mu_RT(m_kk, 0.0);
s.getChemPotentials(mu_RT.begin());
doublereal rrt = 1.0/(GasConstant*m_phase->temperature());
scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt);
int j = 0;
for (k = m_kk - 1; k >= 0; k--) {
ksp = index[k];
if ( mol[k] > 0.0 ) {
kc[j] = ksp;
j++;
if (j == m_mm) break;
}
}
//if (j < m_mm)
// return -1;
//throw CanteraError("estimateElementPotentials",
// "too few species (" + int2str(j) + ").");
for (m = 0; m < j; m++) {
for (m = 0; m < m_mm; m++) {
for (n = 0; n < m_mm; n++) {
aa(m,n) = nAtoms(kc[m], n);
aa(m,n) = nAtoms(m_component[m], n);
}
b[m] = mu_RT[kc[m]];
}
for (m = j+1; m < m_mm; m++) {
aa(m,m) = 1.0;
b[m] = lambda[m];
b[m] = mu_RT[m_component[m]];
}
int info;
@ -332,23 +272,27 @@ namespace Cantera {
info = solve(aa, b.begin());
}
catch (CanteraError) {
return -2;
addLogEntry("failed to estimate initial element potentials.");
info = -2;
}
if (info == 0) {
for (m = 0; m < m_mm; m++) {
lambda[m] = b[m];
addLogEntry(m_phase->elementName(m),b[m]);
}
}
endLogGroup();
return info;
}
/**
* Equilibrate a phase, holding the elemental composition fixed
* at the initial vaollue.
*/
int ChemEquil::equilibrate(thermo_t& s, int XY) {
int ChemEquil::equilibrate(thermo_t& s, const char* XY) {
vector_fp emol(s.nElements());
initialize(s);
update(s);
@ -356,13 +300,13 @@ namespace Cantera {
emol.begin());
return equilibrate(s, XY, emol);
}
/**
* compute the equilibrium composition for 2 specified
* properties and specified element moles.
*/
int ChemEquil::equilibrate(thermo_t& s, int XY, vector_fp& elMoles)
int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp& elMoles)
{
doublereal xval, yval;
int fail = 0;
@ -370,35 +314,51 @@ namespace Cantera {
delete m_p1;
delete m_p2;
bool tempFixed = true;
int XY = _equilflag(XYstr);
vector_fp state;
s.saveState(state);
beginLogGroup("ChemEquil::equilibrate");
initialize(s);
update(s);
switch (XY) {
case TP: case PT:
m_p1 = new TemperatureCalculator<thermo_t>;
m_p2 = new PressureCalculator<thermo_t>; break;
m_p2 = new PressureCalculator<thermo_t>;
break;
case HP: case PH:
tempFixed = false;
m_p1 = new EnthalpyCalculator<thermo_t>;
m_p2 = new PressureCalculator<thermo_t>; break;
m_p2 = new PressureCalculator<thermo_t>;
break;
case SP: case PS:
tempFixed = false;
m_p1 = new EntropyCalculator<thermo_t>;
m_p2 = new PressureCalculator<thermo_t>; break;
m_p2 = new PressureCalculator<thermo_t>;
break;
case SV: case VS:
tempFixed = false;
m_p1 = new EntropyCalculator<thermo_t>;
m_p2 = new DensityCalculator<thermo_t>; break;
m_p2 = new DensityCalculator<thermo_t>;
break;
case TV: case VT:
m_p1 = new TemperatureCalculator<thermo_t>;
m_p2 = new DensityCalculator<thermo_t>; break;
m_p2 = new DensityCalculator<thermo_t>;
break;
case UV: case VU:
tempFixed = false;
m_p1 = new IntEnergyCalculator<thermo_t>;
m_p2 = new DensityCalculator<thermo_t>; break;
m_p2 = new DensityCalculator<thermo_t>;
break;
default:
throw CanteraError("equilibrate","illegal property pair.");
}
addLogEntry("Problem type","fixed "+m_p1->symbol()+", "+m_p2->symbol());
addLogEntry(m_p1->symbol(), m_p1->value(s));
addLogEntry(m_p2->symbol(), m_p2->value(s));
// If the temperature is one of the specified variables, and
// it is outside the valid range, throw an exception.
if (tempFixed) {
@ -421,41 +381,44 @@ namespace Cantera {
vector_fp x(nvar, -102.0); // solution vector
vector_fp res_trial(nvar);
vector_fp elementMol(mm, 0.0);
double perturb;
for (m = 0; m < mm; m++) {
if (m_skip < 0 && elMoles[m] > 0.0 ) m_skip = m;
#define PERTURB_ELEMENT_MOLES
#ifdef PERTURB_ELEMENT_MOLES
perturb = Cutoff*(1.0 + rand());
#else
perturb = 0.0;
#endif
elementMol[m] = elMoles[m] + perturb;
}
// start with a composition with everything non-zero. Note
// that since we have already save the target element moles,
// changing the composition at this point only affects the
// starting point, not the final solution.
vector_fp xmm(m_kk,0.0);
for (int k = 0; k < m_kk; k++) {
xmm[k] = m_phase->moleFraction(k) + Cutoff;
}
m_phase->setMoleFractions(xmm.begin());
update(s);
// loop to estimate T
if (!tempFixed) {
beginLogGroup("Initial T Estimate");
doublereal tmax = m_thermo->maxTemp();
doublereal tmin = m_thermo->minTemp();
doublereal slope, phigh, plow, pval, dt;
// first get the property values at the upper and lower
// temperature limits. Since p1 (h, s, or u) is monotonic
// in T, these values determine the upper and lower
// bounnds (phigh, plow) for p1.
m_phase->setTemperature(tmax);
setInitialMoles(s, elementMol);
setInitialMoles(s);
phigh = m_p1->value(s);
m_phase->setTemperature(tmin);
setInitialMoles(s, elementMol);
setInitialMoles(s);
plow = m_p1->value(s);
// start with T at the midpoint of the range
@ -466,7 +429,7 @@ namespace Cantera {
for (int it = 0; it < 5; it++) {
// set the composition and get p1
setInitialMoles(s, elementMol);
setInitialMoles(s);
pval = m_p1->value(s);
@ -493,19 +456,22 @@ namespace Cantera {
// update the T estimate
t0 = tmin + dt;
addLogEntry("new T estimate", t0);
m_phase->setTemperature(t0);
}
endLogGroup(); // initial T estimate
}
if (m_lambda[0] == -100.0) {
setInitialMoles(s, elementMol);
for (int ii = 0; ii < m_mm; ii++) x[ii] = -101.0;
estimateElementPotentials(s, x);
}
else {
doublereal rt = GasConstant * m_phase->temperature();
for (int ii = 0; ii < m_mm; ii++) x[ii] = m_lambda[ii]/rt;
}
//if (m_lambda[0] == -100.0) {
setInitialMoles(s);
for (int ii = 0; ii < m_mm; ii++) x[ii] = -101.0;
estimateElementPotentials(s, x);
//}
//else {
// doublereal rt = GasConstant * m_phase->temperature();
// for (int ii = 0; ii < m_mm; ii++) x[ii] = m_lambda[ii]/rt;
//}
x[m_mm] = log(m_phase->temperature());
@ -513,17 +479,16 @@ namespace Cantera {
vector_fp below(nvar);
for (m = 0; m < mm; m++) {
above[m] = 200.0; // 30.0;
above[m] = 200.0;
below[m] = -2000.0;
if (elMoles[m] < Cutoff && m != m_eloc) x[m] = -1000.0;
//if (m == m_eloc) x[m] = -10.0;
}
above[mm] = log(m_thermo->maxTemp() + 1.0); //log(1.e4);
below[mm] = log(m_thermo->minTemp() - 1.0); //log(10.0);
above[mm] = log(m_thermo->maxTemp() + 1.0);
below[mm] = log(m_thermo->minTemp() - 1.0);
vector_fp grad(nvar, 0.0); // gradient of f = F*F/2
vector_fp oldx(nvar, 0.0); // old solution
vector_fp prevx(nvar, 0.0); // old solution
//vector_fp prevx(nvar, 0.0); // old solution
vector_fp oldresid(nvar, 0.0);
doublereal f, oldf;
@ -549,30 +514,33 @@ namespace Cantera {
iter++;
if (iter > 1) endLogGroup(); // iteration
beginLogGroup("Iteration "+int2str(iter));
// compute the residual and the jacobian using the current
// solution vector
equilResidual(s, x, elMoles, res_trial, XY, xval, yval);
equilResidual(s, x, elMoles, res_trial, xval, yval);
f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
equilJacobian(s, x, elMoles, jac, XY, xval, yval);
addLogEntry("Residual norm", f);
equilJacobian(s, x, elMoles, jac, xval, yval);
// compute grad f = F*J
jac.leftMult(res_trial.begin(), grad.begin());
copy(x.begin(), x.end(), oldx.begin());
copy(oldx.begin(), oldx.end(), prevx.begin());
oldf = f;
scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0);
try {
info = solve(jac, res_trial.begin());
}
catch (CanteraError) {
cout << x << endl;
//cout << res_trial << endl;
//cout << grad << endl;
cout << elMoles << endl;
//cout << jac << endl;
addLogEntry("Jacobian is singular.");
endLogGroup(); // iteration
endLogGroup(); // equilibrate
s.restoreState(state);
//cout << "m_skip = " << m_skip << endl;
throw CanteraError("equilibrate",
"Jacobian is singular. \nTry adding more species, "
"changing the elemental composition slightly, \nor removing "
@ -593,15 +561,19 @@ namespace Cantera {
fctr = fminn(fctr, 0.8*(x[m] - below[m])/(x[m] - newval));
}
}
if (fctr != 1.0) addLogEntry("factor to keep solution in bounds",
fctr);
// multiply the step by the scaing factor
// multiply the step by the scaling factor
scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr);
if (!dampStep(s, oldx, oldf, grad, res_trial,
x, f, elMoles , XY, xval, yval))
x, f, elMoles , xval, yval))
{
fail++;
if (fail > 3) {
addLogEntry("dampStep","Failed 3 times. Giving up.");
s.restoreState(state);
throw CanteraError("equilibrate",
"Cannot find an acceptable Newton damping coefficient.");
return -4;
@ -613,24 +585,32 @@ converge:
// check for convergence.
equilResidual(s, x, elMoles, res_trial, XY, xval, yval);
equilResidual(s, x, elMoles, res_trial, xval, yval);
f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin());
doublereal xx, yy, deltax, deltay;
xx = m_p1->value(s);
yy = m_p2->value(s);
deltax = (xx - xval)/xval;
deltay = (yy - yval)/yval;
if (absmax(res_trial.begin(), res_trial.end()) < options.relTolerance
doublereal rmax = absmax(res_trial.begin(), res_trial.end());
if (iter > 0 && rmax < options.relTolerance
&& fabs(deltax) < options.relTolerance
&& fabs(deltay) < options.relTolerance) {
options.iterations = iter;
endLogGroup(); // iteration
m_lambda.resize(m_mm);
beginLogGroup("Converged solution");
addLogEntry("Iterations",iter);
addLogEntry("Relative error in "+m_p1->symbol(),deltax);
addLogEntry("Relative error in "+m_p2->symbol(),deltay);
addLogEntry("Max residual",rmax);
beginLogGroup("Element potentials");
doublereal rt = GasConstant*m_thermo->temperature();
for (m = 0; m < m_mm; m++) {
m_lambda[m] = x[m]*rt;
addLogEntry("element "+m_phase->elementName(m), fp2str(x[m]));
}
endLogGroup(); // element potentials
if (m_thermo->temperature() > m_thermo->maxTemp() + 1.0 ||
m_thermo->temperature() < m_thermo->minTemp() - 1.0 ) {
@ -639,13 +619,19 @@ converge:
"valid range of "+fp2str(m_thermo->minTemp())+" K to "
+fp2str(m_thermo->maxTemp())+" K\n");
}
endLogGroup(); // converged solution
endLogGroup(); // equilibrate
return 0;
}
// no convergence
if (iter > options.maxIterations) {
addLogEntry("equilibrate","no convergence");
endLogGroup(); // iteration
endLogGroup(); // equilibrate
s.restoreState(state);
throw CanteraError("equilibrate",
"no convergence in "+int2str(options.maxIterations)
+" iterations.");
@ -657,7 +643,7 @@ converge:
int ChemEquil::dampStep(thermo_t& mix, vector_fp& oldx,
double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
double& f, vector_fp& elmols, int XY, double xval, double yval )
double& f, vector_fp& elmols, double xval, double yval )
{
int nvar = x.size();
@ -696,7 +682,7 @@ converge:
add_each(x, oldx);
equilResidual(mix, x, elmols, res_new, XY, xval, yval);
equilResidual(mix, x, elmols, res_new, xval, yval);
//f = 0.5*(res_new*res_new);
f = 0.5*dot(res_new.begin(), res_new.end(), res_new.begin());
if (damp < minDamp && damp < 1.0)
@ -748,8 +734,9 @@ converge:
*/
void ChemEquil::equilResidual(thermo_t& mix, const vector_fp& x,
const vector_fp& elmtotal, vector_fp& resid,
int XY, doublereal xval, doublereal yval)
doublereal xval, doublereal yval)
{
beginLogGroup("ChemEquil::equilResidual");
int n;
doublereal xx, yy;
doublereal temp = exp(x[m_mm]);
@ -763,8 +750,9 @@ converge:
if (elmtotal[n] < Cutoff && n != m_eloc)
resid[n] = x[n] + 1000.0;
else
// resid[n] = elmtotal[n] - elm[n]; // log( (1.0 + elmtotal[n]) / (1.0 + elm[n]) );
resid[n] = log( (1.0 + elmtotal[n]) / (1.0 + elm[n]) );
addLogEntry(m_phase->elementName(n),fp2str(elm[n])+" ("
+fp2str(elmtotal[n])+")");
}
if (m_eloc >= 0) {
doublereal chrg, sumnet = 0.0, sumabs = 0.0;
@ -773,12 +761,16 @@ converge:
sumnet += chrg;
sumabs += fabs(chrg);
}
addLogEntry("net charge",sumnet);
resid[m_eloc] = sumnet/m_abscharge; // log((1.0 + sumnet/sumabs));
}
xx = m_p1->value(mix);
yy = m_p2->value(mix);
resid[m_mm] = xx/xval - 1.0;
resid[m_skip] = yy/yval - 1.0;
resid[m_skip] = yy/yval - 1.0;
addLogEntry(m_p1->symbol(), fp2str(xx)+" ("+fp2str(xval)+")");
addLogEntry(m_p2->symbol(), fp2str(yy)+" ("+fp2str(yval)+")");
endLogGroup();
}
@ -786,8 +778,10 @@ converge:
void ChemEquil::equilJacobian(thermo_t& mix, vector_fp& x,
const vector_fp& elmols, DenseMatrix& jac,
int XY, doublereal xval, doublereal yval)
doublereal xval, doublereal yval)
{
beginLogGroup("equilJacobian",0);
int len = x.size();
vector_fp& r0 = m_jwork1;
vector_fp& r1 = m_jwork2;
@ -798,7 +792,7 @@ converge:
doublereal rdx, dx, xsave;
doublereal atol = 1.e-10;
equilResidual(mix, x, elmols, r0, XY, xval, yval);
equilResidual(mix, x, elmols, r0, xval, yval);
for (n = 0; n < len; n++)
{
@ -812,7 +806,7 @@ converge:
// calculate perturbed residual
equilResidual(mix, x, elmols, r1, XY, xval, yval);
equilResidual(mix, x, elmols, r1, xval, yval);
// compute nth column of Jacobian
@ -821,9 +815,13 @@ converge:
}
x[n] = xsave;
}
endLogGroup();
}
} // namespace
// $Log: ChemEquil.cpp,v

77
Cantera/src/ChemEquil.h
View file

@ -57,7 +57,7 @@ namespace Cantera {
* Continuation flag. Set true if the calculation should be
* initialized from the last calculation. Otherwise, the
* calculation will be started from scratch and the initial
* composition and element potentials estimated. (Not Implemented.)
* composition and element potentials estimated.
*/
bool contin;
};
@ -66,17 +66,35 @@ namespace Cantera {
class PropertyCalculator;
/**
* Chemical equilibrium processor. Sets a mixture to a state of
* chemical equilibrium.
*/
* Class ChemEquil implements a chemical equilibrium solver for
* single-phase solutions. It is a "non-stoichiometric" solver in
* the terminology of Smith and Missen, meaning that every
* intermediate state is a valid chemical equilibrium state, but
* does not necessarily satisfy the element constraints. In
* contrast, the solver implemented in class MultiPhaseEquil uses
* a "stoichiometric" algorithm, in which each intermediate state
* satisfies the element constraints but is not a state of
* chemical equilibrium. Non-stoichiometric methods are faster
* when they converge, but stoichiometric ones tend to be more
* robust and can be used also for problems with multiple
* condensed phases. As expected, the ChemEquil solver is faster
* than MultiPhaseEquil for many single-phase equilibrium
* problems (particularly if there are only a few elements but
* vvery many species), but can be less stable. Problem
* situations include low temperatures where only a few species
* have non-zero mole fractions, precisely stoichiometric
* compositions (e.g. 2 H2 + O2). In general, if speed is
* important, this solver should be tried first, and if it fails
* then use MultiPhaseEquil.
*/
class ChemEquil {
public:
ChemEquil();
virtual ~ChemEquil();
int equilibrate(thermo_t& s, int XY = 0);
int equilibrate(thermo_t& s, int XY, vector_fp& elMoles);
int equilibrate(thermo_t& s, const char* XY);
int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles);
const vector_fp& elementPotentials() const { return m_lambda; }
/**
@ -99,22 +117,22 @@ namespace Cantera {
void setToEquilState(thermo_t& s,
const vector_fp& x, doublereal t);
int setInitialMoles(thermo_t& s, vector_fp& elementMoles);
int setInitialMoles(thermo_t& s);
int estimateElementPotentials(thermo_t& s,
vector_fp& lambda);
int dampStep(thermo_t& s, vector_fp& oldx,
double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
double& f, vector_fp& elmols, int XY, double xval, double yval );
double& f, vector_fp& elmols, double xval, double yval );
void equilResidual(thermo_t& s, const vector_fp& x,
const vector_fp& elmtotal, vector_fp& resid,
int XY, double xval, double yval);
double xval, double yval);
void equilJacobian(thermo_t& s, vector_fp& x,
const vector_fp& elmols, DenseMatrix& jac,
int XY, double xval, double yval);
double xval, double yval);
void update(const thermo_t& s);
@ -141,46 +159,9 @@ namespace Cantera {
vector_fp m_grt;
vector_fp m_mu_RT;
vector_int m_component;
};
//-----------------------------------------------------------
// convenience functions
//-----------------------------------------------------------
/**
* Set a mixture to a state of chemical equilibrium. The flag 'XY'
* determines the two properties that will be held fixed in the
* calculation.
*/
inline void equilibrate(thermo_t& s, int XY, int solver = 0,
doublereal rtol = 1.0e-9, int maxsteps = 1000,
int loglevel = 0) {
if (solver > 0) {
MultiPhase mix;
mix.addPhase(&s, 1.0);
mix.init();
mix.setTemperature(s.temperature());
mix.setPressure(s.pressure());
equilibrate(mix, XY, rtol, maxsteps, loglevel);
}
else {
ChemEquil e;
e.equilibrate(s,XY);
s.setElementPotentials(e.elementPotentials());
}
}
/**
* Set a mixture to a state of chemical equilibrium. The flag 'XY'
* determines the two properties that will be held fixed in the
* calculation.
*/
inline void equilibrate(thermo_t& s, const char* XY,
int solver = 0) {
equilibrate(s,_equilflag(XY), solver);
}
}

161
Cantera/src/MultiPhase.cpp
View file

@ -412,23 +412,18 @@ namespace Cantera {
MultiPhaseEquil* e = 0;
if (!m_init) init();
if (loglevel > 0) {
beginLogGroup("MultiPhase::equilibrate");
}
beginLogGroup("MultiPhase::equilibrate", loglevel);
if (XY == TP) {
if (loglevel > 0) {
addLogEntry("problem type","fixed T,P");
}
addLogEntry("problem type","fixed T,P");
// create an equilibrium manager
e = new MultiPhaseEquil(this);
try {
error = e->equilibrate(XY, err, maxsteps, loglevel-1);
error = e->equilibrate(XY, err, maxsteps);
}
catch (CanteraError err) {
if (loglevel > 0) {
endLogGroup();
//write_logfile("equil_err.html");
}
endLogGroup();
delete e;
e = 0;
throw err;
@ -440,12 +435,9 @@ namespace Cantera {
h0 = enthalpy();
Tlow = 0.5*m_Tmin; // lower bound on T
Thigh = 2.0*m_Tmax; // upper bound on T
if (loglevel > 0) {
addLogEntry("problem type","fixed H,P");
addLogEntry("H target",fp2str(h0));
//addLogEntry("min T",fp2str(Tlow));
//addLogEntry("max T",fp2str(Thigh));
}
addLogEntry("problem type","fixed H,P");
addLogEntry("H target",fp2str(h0));
for (n = 0; n < maxiter; n++) {
// if 'strt' is false, the current composition will be used as
@ -453,11 +445,10 @@ namespace Cantera {
e = new MultiPhaseEquil(this, strt);
// start with a loose error tolerance, but tighten it as we get
// close to the final temperature
if (loglevel > 0) {
beginLogGroup("iteration "+int2str(n));
}
beginLogGroup("iteration "+int2str(n));
try {
error = e->equilibrate(TP, err, maxsteps, loglevel-1);
error = e->equilibrate(TP, err, maxsteps);
hnow = enthalpy();
// the equilibrium enthalpy monotonically increases with T;
// if the current value is below the target, the we know the
@ -490,22 +481,18 @@ namespace Cantera {
}
herr = fabs((h0 - hnow)/h0);
if (loglevel > 0) {
addLogEntry("T",fp2str(temperature()));
addLogEntry("H",fp2str(hnow));
addLogEntry("H rel error",fp2str(herr));
addLogEntry("lower T bound",fp2str(Tlow));
addLogEntry("upper T bound",fp2str(Thigh));
endLogGroup();
}
addLogEntry("T",fp2str(temperature()));
addLogEntry("H",fp2str(hnow));
addLogEntry("H rel error",fp2str(herr));
addLogEntry("lower T bound",fp2str(Tlow));
addLogEntry("upper T bound",fp2str(Thigh));
endLogGroup(); // iteration
if (herr < err) { // || dta < 1.0e-4) {
if (loglevel > 0) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("H rel error",fp2str(herr));
}
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("H rel error",fp2str(herr));
goto done;
}
tnew = m_temp + dt;
@ -521,27 +508,24 @@ namespace Cantera {
catch (CanteraError err) {
if (!strt) {
if (loglevel > 0)
addLogEntry("no convergence",
"try estimating starting composition");
addLogEntry("no convergence",
"try estimating starting composition");
strt = true;
}
else {
tnew = 0.5*(Tlow + Thigh);
tnew = 0.5*(m_temp + Thigh);
if (fabs(tnew - m_temp) < 1.0) tnew = m_temp + 1.0;
setTemperature(tnew);
if (loglevel > 0)
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
}
endLogGroup();
}
}
delete e;
e = 0;
if (loglevel > 0) {
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
}
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
throw CanteraError("MultiPhase::equilibrate",
"No convergence for T");
}
@ -551,22 +535,20 @@ namespace Cantera {
start = true;
Tlow = 1.0; // m_Tmin; // lower bound on T
Thigh = 1.0e6; // m_Tmax; // upper bound on T
if (loglevel > 0) {
addLogEntry("problem type","fixed S,P");
addLogEntry("S target",fp2str(s0));
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
}
addLogEntry("problem type","fixed S,P");
addLogEntry("S target",fp2str(s0));
addLogEntry("min T",fp2str(Tlow));
addLogEntry("max T",fp2str(Thigh));
for (n = 0; n < maxiter; n++) {
e = new MultiPhaseEquil(this, strt);
ferr = 0.1;
if (fabs(dt) < 1.0) ferr = err;
//start = false;
if (loglevel > 0) {
beginLogGroup("iteration "+int2str(n));
}
beginLogGroup("iteration "+int2str(n));
try {
error = e->equilibrate(TP, err, maxsteps, loglevel-1);
error = e->equilibrate(TP, err, maxsteps);
snow = entropy();
if (snow < s0) {
if (m_temp > Tlow) Tlow = m_temp;
@ -575,23 +557,20 @@ namespace Cantera {
if (m_temp < Thigh) Thigh = m_temp;
}
serr = fabs((s0 - snow)/s0);
if (loglevel > 0) {
addLogEntry("T",fp2str(temperature()));
addLogEntry("S",fp2str(snow));
addLogEntry("S rel error",fp2str(serr));
endLogGroup();
}
addLogEntry("T",fp2str(temperature()));
addLogEntry("S",fp2str(snow));
addLogEntry("S rel error",fp2str(serr));
endLogGroup();
dt = (s0 - snow)*m_temp/cp();
dtmax = 0.5*fabs(Thigh - Tlow);
dtmax = (dtmax > 500.0 ? 500.0 : dtmax);
dta = fabs(dt);
if (dta > dtmax) dt *= dtmax/dta;
if (herr < err || dta < 1.0e-4) {
if (loglevel > 0) {
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("S rel error",fp2str(serr));
}
addLogEntry("T iterations",int2str(n));
addLogEntry("Final T",fp2str(temperature()));
addLogEntry("S rel error",fp2str(serr));
goto done;
}
tnew = m_temp + dt;
@ -604,17 +583,15 @@ namespace Cantera {
catch (CanteraError err) {
if (!strt) {
if (loglevel > 0)
addLogEntry("no convergence",
"setting strt to True");
addLogEntry("no convergence",
"setting strt to True");
strt = true;
}
else {
tnew = 0.5*(m_temp + Thigh);
setTemperature(tnew);
if (loglevel > 0)
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
addLogEntry("no convergence",
"trying T = "+fp2str(m_temp));
}
endLogGroup();
@ -622,10 +599,8 @@ namespace Cantera {
}
delete e;
e = 0;
if (loglevel > 0) {
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
}
addLogEntry("reached max number of T iterations",int2str(maxiter));
endLogGroup();
throw CanteraError("MultiPhase::equilibrate",
"No convergence for T");
}
@ -691,9 +666,7 @@ namespace Cantera {
// "No convergence for T");
// }
else if (XY == TV) {
if (loglevel > 0) {
addLogEntry("problem type","fixed T, V");
}
addLogEntry("problem type","fixed T, V");
doublereal dt = 1.0e3;
doublereal v0 = volume();
doublereal dVdP;
@ -704,23 +677,19 @@ namespace Cantera {
pnow = pressure();
MultiPhaseEquil e(this, start);
start = false;
if (loglevel > 1) {
beginLogGroup("iteration "+int2str(n));
}
error = e.equilibrate(TP, err, maxsteps, loglevel-1);
beginLogGroup("iteration "+int2str(n));
error = e.equilibrate(TP, err, maxsteps);
vnow = volume();
verr = fabs((v0 - vnow)/v0);
if (loglevel > 1) {
addLogEntry("P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
endLogGroup();
}
addLogEntry("P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
endLogGroup();
if (verr < err) {
if (loglevel > 0) {
addLogEntry("P iterations",int2str(n));
addLogEntry("Final P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
}
addLogEntry("P iterations",int2str(n));
addLogEntry("Final P",fp2str(pressure()));
addLogEntry("V rel error",fp2str(verr));
goto done;
}
// find dV/dP
@ -731,16 +700,14 @@ namespace Cantera {
}
else {
if (loglevel > 0) endLogGroup();
endLogGroup();
throw CanteraError("MultiPhase::equilibrate","unknown option");
}
return -1.0;
done:
delete e;
e = 0;
if (loglevel > 0) {
endLogGroup();
}
endLogGroup();
return err;
}

2
Cantera/src/MultiPhase.h
View file

@ -149,7 +149,7 @@ namespace Cantera {
/// @param loglevel Level of diagnostic output, written to a
/// file in HTML format.
doublereal equilibrate(int XY, doublereal err = 1.0e-9,
int maxsteps = 1000, int maxiter = 200, int loglevel = 0);
int maxsteps = 1000, int maxiter = 200, int loglevel = -99);
/// Set the temperature [K].

104
Cantera/src/MultiPhaseEquil.cpp
View file

@ -64,7 +64,9 @@ namespace Cantera {
m_incl_element.resize(m_nel_mix,1);
for (m = 0; m < m_nel_mix; m++) {
string enm = mix->elementName(m);
if (enm == "E" || enm == "e") m_eloc = m;
if (enm == "E" || enm == "e") {
m_eloc = m;
}
if (m_mix->elementMoles(m) <= 0.0) {
if (m != m_eloc) {
m_incl_element[m] = 0;
@ -163,36 +165,29 @@ namespace Cantera {
int maxsteps, int loglevel) {
int i;
m_iter = 0;
if (loglevel > 0)
beginLogGroup("MultiPhaseEquil::equilibrate");
string iterstr;
beginLogGroup("MultiPhaseEquil::equilibrate", loglevel);
for (i = 0; i < maxsteps; i++) {
stepComposition(loglevel - 1);
if (loglevel > 1) {
beginLogGroup("iteration "+int2str(i));
addLogEntry("error",fp2str(error()));
endLogGroup();
}
if (loglevel > 2) printInfo();
iterstr = "iteration "+int2str(i);
beginLogGroup(iterstr);
stepComposition();
addLogEntry("error",fp2str(error()));
endLogGroup(iterstr);
if (error() < err) break;
}
if (i >= maxsteps) {
if (loglevel > 0) {
addLogEntry("Error","no convergence in "+int2str(maxsteps)
+" iterations");
if (loglevel > 2) printInfo();
endLogGroup();
}
addLogEntry("Error","no convergence in "+int2str(maxsteps)
+" iterations");
endLogGroup("MultiPhaseEquil::equilibrate");
throw CanteraError("MultiPhaseEquil::equilibrate",
"no convergence in " + int2str(maxsteps) +
" iterations. Error = " + fp2str(error()));
}
if (loglevel > 0) {
addLogEntry("iterations",int2str(iterations()));
addLogEntry("error tolerance",fp2str(err));
addLogEntry("error",fp2str(error()));
endLogGroup();
}
addLogEntry("iterations",int2str(iterations()));
addLogEntry("error tolerance",fp2str(err));
addLogEntry("error",fp2str(error()));
endLogGroup("MultiPhaseEquil::equilibrate");
finish();
return error();
}
@ -282,29 +277,30 @@ namespace Cantera {
// set the moles of the phase objects to match
updateMixMoles();
}
return 0;
}
/// This method finds a set of constituent species and a complete
/// set of formation reactions for the non-constituents in terms
/// of the constituents. Note that in most cases, many different
/// constituent sets are possible, and therefore neither the
/// constituents returned by this method nor the formation
/// This method finds a set of component species and a complete
/// set of formation reactions for the non-components in terms of
/// the components. Note that in most cases, many different
/// component sets are possible, and therefore neither the
/// components returned by this method nor the formation
/// reactions are unique. The algorithm used here is described in
/// Smith and Missen, Chemical Reaction Equilibrium Analysis.
///
/// The constituent species are taken to be the first M species
/// The component species are taken to be the first M species
/// in array 'species' that have linearly-independent compositions.
///
/// @param order On entry, vector \a order should contain species
/// index numbers in the order of decreasing desirability as a
/// constituent. For example, if it is desired to choose the
/// constituents from among the major species, this array might
/// component. For example, if it is desired to choose the
/// components from among the major species, this array might
/// list species index numbers in decreasing order of mole
/// fraction. If array 'species' does not have length =
/// nSpecies(), then the species will be considered as candidates
/// to be constituents in declaration order, beginning with the
/// to be components in declaration order, beginning with the
/// first phase added.
///
void MultiPhaseEquil::getComponents(const vector_int& order) {
@ -312,7 +308,7 @@ namespace Cantera {
int n;
// if the input species array has the wrong size, ignore it
// and consider the species for constituents in declarationi order.
// and consider the species for components in declarationi order.
if (order.size() != m_nsp) {
for (k = 0; k < m_nsp; k++) m_order[k] = k;
}
@ -335,8 +331,20 @@ namespace Cantera {
// Do Gauss elimination
for (m = 0; m < nRows; m++) {
// if a pivot is zero, exchange columns
// If a pivot is zero, exchange columns. This occurs when
// a species has an elemental composition that is not
// linearly independent of the component species that have
// already been assigned
if (m_A(m,m) == 0.0) {
// First, we need to find a good candidate for a
// component species to swap in for the one that has
// zero pivot. It must contain element m, be linearly
// independent of the components processed so far
// (m_A(m,k) != 0), and should be a major species if
// possible. We'll choose the species with greatest
// mole fraction that satisfies these criteria.
doublereal maxmoles = -999.0;
index_t kmax = 0;
for (k = m+1; k < nColumns; k++) {
@ -347,11 +355,15 @@ namespace Cantera {
}
}
}
// Now exchange the column with zero pivot with the
// column for this major species
for (n = 0; n < int(nRows); n++) {
tmp = m_A(n,m);
m_A(n, m) = m_A(n, kmax);
m_A(n, kmax) = tmp;
}
}
// exchange the species labels on the columns
itmp = m_order[m];
m_order[m] = m_order[kmax];
@ -365,8 +377,8 @@ namespace Cantera {
m_A(m,k) *= fctr;
}
// subtract A(n,m)/A(m,m) * (row m) from row n, so that
// A(n,m) = 0.
// For all rows below the diagonal, subtract A(n,m)/A(m,m)
// * (row m) from row n, so that A(n,m) = 0.
for (n = int(m+1); n < int(m_nel); n++) {
fctr = m_A(n,m)/m_A(m,m);
for (k = 0; k < m_nsp; k++) {
@ -376,8 +388,8 @@ namespace Cantera {
}
// The left m_nel columns of A are now upper-diagonal.
// Now reduce it to diagonal form by back-solving
// The left m_nel columns of A are now upper-diagonal. Now
// reduce the m_nel columns to diagonal form by back-solving
for (m = nRows-1; m > 0; m--) {
for (n = m-1; n>= 0; n--) {
if (m_A(n,m) != 0.0) {
@ -411,6 +423,9 @@ namespace Cantera {
}
}
/// Re-arrange a vector of species properties in sorted form
/// (components first) into unsorted, sequential form.
void MultiPhaseEquil::unsort(vector_fp& x) {
@ -421,6 +436,7 @@ namespace Cantera {
}
}
void MultiPhaseEquil::printInfo() {
index_t m, ik, k;
beginLogGroup("info");
@ -498,9 +514,9 @@ namespace Cantera {
/// Take one step in composition, given the gradient of G at the
/// starting point, and a vector of reaction steps dxi.
doublereal MultiPhaseEquil::
stepComposition(int loglevel) {
stepComposition() {
if (loglevel > 0) beginLogGroup("MultiPhaseEquil::stepComposition");
beginLogGroup("MultiPhaseEquil::stepComposition");
m_iter++;
index_t ik, j, k = 0;
@ -555,7 +571,7 @@ namespace Cantera {
}
}
}
if (loglevel > 0 && m_moles[k] < -Tiny) {
if (m_moles[k] < -Tiny) {
addLogEntry("Negative moles for "
+m_mix->speciesName(m_species[k]), fp2str(m_moles[k]));
}
@ -564,7 +580,7 @@ namespace Cantera {
}
// now take a step with this scaled omega
if (loglevel > 0) addLogEntry("Stepping by ", fp2str(omegamax));
addLogEntry("Stepping by ", fp2str(omegamax));
step(omegamax, m_work);
// compute the gradient of G at this new position in the
@ -581,10 +597,10 @@ namespace Cantera {
if (grad1 > 0.0) {
omega *= fabs(grad0) / (grad1 + fabs(grad0));
for (k = 0; k < m_nsp; k++) m_moles[k] = m_lastmoles[k];
if (loglevel > 0) addLogEntry("Stepped over minimum. Take smaller step ", fp2str(omega));
addLogEntry("Stepped over minimum. Take smaller step ", fp2str(omega));
step(omega, m_work);
}
if (loglevel > 0) endLogGroup();
endLogGroup("MultiPhaseEquil::stepComposition");
return omega;
}

50
Cantera/src/MultiPhaseEquil.h
View file

@ -6,8 +6,6 @@
namespace Cantera {
int _equilflag(const char* xy);
class MultiPhaseEquil {
public:
@ -37,19 +35,25 @@ namespace Cantera {
int iterations() { return m_iter; }
doublereal equilibrate(int XY, doublereal err = 1.0e-9,
int maxsteps = 1000, int loglevel=0);
int maxsteps = 1000, int loglevel=-99);
string reactionString(index_t j);
doublereal error();
void printInfo();
void setInitialMixMoles() {
setInitialMoles();
finish();
}
index_t componentIndex(index_t n) { return m_species[m_order[n]]; }
protected:
void getComponents(const vector_int& order);
int setInitialMoles();
int setInitialMoles2();
void computeN();
doublereal stepComposition(int loglevel);
doublereal stepComposition();
//void sort(vector_fp& x);
void unsort(vector_fp& x);
void step(doublereal omega, vector_fp& deltaN);
@ -92,42 +96,6 @@ namespace Cantera {
bool m_force;
};
//-----------------------------------------------------------
// convenience functions
//-----------------------------------------------------------
/**
* Set a mixture to a state of chemical equilibrium. The flag 'XY'
* determines the two properties that will be held fixed in the
* calculation.
*/
inline doublereal equilibrate(MultiPhase& s, int XY,
doublereal tol = 1.0e-9, int maxsteps = 1000, int loglevel = 0) {
s.init();
writelog("in equilibrate(MultiPhase, ...)\n");
if (XY == TP || XY == HP || XY == SP || XY == TV) {
double err = s.equilibrate(XY, tol, maxsteps, maxsteps, loglevel);
if (loglevel > 0) {
writelog("writing log file\n");
write_logfile("equilibrate.html");
}
return err;
}
else {
throw CanteraError("equilibrate","unsupported option");
return -1.0;
}
}
/**
* Set a mixture to a state of chemical equilibrium. The flag 'XY'
* determines the two properties that will be held fixed in the
* calculation.
*/
inline doublereal equilibrate(MultiPhase& s, const char* XY,
doublereal tol = 1.0e-9, int maxsteps = 1000, int loglevel = 0) {
return equilibrate(s,_equilflag(XY), tol, maxsteps, loglevel);
}
}

23
Cantera/src/equil.h Normal file
View file

@ -0,0 +1,23 @@
#ifndef CT_KERNEL_EQUIL_H
#define CT_KERNEL_EQUIL_H
//#include "ChemEquil.h"
#include "MultiPhase.h"
namespace Cantera {
//-----------------------------------------------------------
// convenience functions
//-----------------------------------------------------------
void equilibrate(thermo_t& s, const char* XY,
int solver = -1, doublereal rtol = 1.0e-9, int maxsteps = 1000,
int maxiter = 100, int loglevel = -99);
doublereal equilibrate(MultiPhase& s, const char* XY,
doublereal tol = 1.0e-9, int maxsteps = 1000, int maxiter = 100,
int loglevel = -99);
}
#endif

142
Cantera/src/equilibrate.cpp Normal file
View file

@ -0,0 +1,142 @@
/**
* @file equilibrate.cpp
*
* Driver routines for the chemical equilibrium solvers.
*
*/
#include "ChemEquil.h"
#include "MultiPhaseEquil.h"
namespace Cantera {
/**
* Set a mixture to a state of chemical equilibrium. The flag 'XY'
* determines the two properties that will be held fixed in the
* calculation.
*/
doublereal equilibrate(MultiPhase& s, const char* XY,
doublereal tol = 1.0e-9, int maxsteps = 1000, int maxiter = 100,
int loglevel = -99) {
beginLogGroup("equilibrate",loglevel);
addLogEntry("multiphase equilibrate function");
beginLogGroup("arguments");
addLogEntry("XY",XY);
addLogEntry("tol",tol);
addLogEntry("maxsteps",maxsteps);
addLogEntry("maxiter",maxiter);
addLogEntry("loglevel",loglevel);
endLogGroup("arguments");
s.init();
int ixy = _equilflag(XY);
if (ixy == TP || ixy == HP || ixy == SP || ixy == TV) {
try {
double err = s.equilibrate(ixy, tol, maxsteps, maxiter);
addLogEntry("Success. Error",err);
endLogGroup("equilibrate");
return err;
}
catch (CanteraError e) {
addLogEntry("Failure.",lastErrorMessage());
endLogGroup("equilibrate");
throw e;
}
}
else {
addLogEntry("multiphase equilibrium can be done only for TP, HP, SP, or TV");
endLogGroup("equilibrate");
throw CanteraError("equilibrate","unsupported option");
return -1.0;
}
}
/// Set a single-phase chemical solution to chemical equilibrium.
/// This is a convenience function that uses one or the other of
/// the two chemical equilibrium solvers. @param The object to
/// set to an equilibrium state @param XY An integer specifying
/// the two properties to be held constant. @param solver The
/// equilibrium solver to use. If solver = 0, the ChemEquil solver
/// will be used, and if solver = 1, the MultiPhaseEquil solver
/// will be used (slower than ChemEquil, but more stable). If
/// solver < 0 (default, then ChemEquil will be tried first, and
/// if it fails MultiPhaseEquil will be tried. @param maxsteps
/// The maximum number of steps to take to find the solution.
/// @param maxiter For the MultiPhaseEquil solver only, this is
/// the maximum number of outer temperature or pressure iterations
/// to take when T and/or P is not held fixed. @param loglevel
/// Controls amount of diagnostic output. loglevel = 0 suppresses
/// diagnostics, and increasingly-verbose messages are written as
/// loglevel increases. The messages are written to a file in HTML
/// format for viewing in a web browser.
void equilibrate(thermo_t& s, const char* XY, int solver,
doublereal rtol, int maxsteps, int maxiter, int loglevel) {
MultiPhase* m = 0;
ChemEquil* e = 0;
bool redo = true;
beginLogGroup("equilibrate", loglevel);
addLogEntry("Single-phase equilibrate function");
{
beginLogGroup("arguments");
addLogEntry("phase",s.id());
addLogEntry("XY",XY);
addLogEntry("solver",solver);
addLogEntry("rtol",rtol);
addLogEntry("maxsteps",maxsteps);
addLogEntry("maxiter",maxiter);
addLogEntry("loglevel",loglevel);
endLogGroup("arguments");
}
while (redo) {
if (solver > 0) {
m = new MultiPhase;
try {
m->addPhase(&s, 1.0);
m->init();
equilibrate(*m, XY, rtol, maxsteps, maxiter, loglevel);
redo = false;
addLogEntry("MultiPhaseEquil solver succeeded.");
delete m;
}
catch (CanteraError err) {
addLogEntry("MultiPhaseEquil solver failed.");
endLogGroup("equilibrate");
delete m;
throw err;
}
}
else { // solver <= 0
e = new ChemEquil;
try {
e->options.maxIterations = maxsteps;
e->options.relTolerance = rtol;
e->equilibrate(s,XY);
s.setElementPotentials(e->elementPotentials());
redo = false;
delete e;
addLogEntry("ChemEquil solver succeeded.");
}
catch (CanteraError err) {
delete e;
addLogEntry("ChemEquil solver failed.");
// If ChemEquil fails, try the MultiPhase solver
if (solver < 0) {
addLogEntry("Trying MultiPhaseEquil solver.");
solver = 1;
}
else {
redo = false;
endLogGroup("equilibrate");
throw err;
}
}
}
} // while (redo)
endLogGroup("equilibrate");
}
}