diff --git a/Cantera/src/ChemEquil.cpp b/Cantera/src/ChemEquil.cpp index ab5f67888..78929d5ed 100755 --- a/Cantera/src/ChemEquil.cpp +++ b/Cantera/src/ChemEquil.cpp @@ -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; - m_p2 = new PressureCalculator; break; + m_p2 = new PressureCalculator; + break; case HP: case PH: tempFixed = false; m_p1 = new EnthalpyCalculator; - m_p2 = new PressureCalculator; break; + m_p2 = new PressureCalculator; + break; case SP: case PS: tempFixed = false; m_p1 = new EntropyCalculator; - m_p2 = new PressureCalculator; break; + m_p2 = new PressureCalculator; + break; case SV: case VS: tempFixed = false; m_p1 = new EntropyCalculator; - m_p2 = new DensityCalculator; break; + m_p2 = new DensityCalculator; + break; case TV: case VT: m_p1 = new TemperatureCalculator; - m_p2 = new DensityCalculator; break; + m_p2 = new DensityCalculator; + break; case UV: case VU: tempFixed = false; m_p1 = new IntEnergyCalculator; - m_p2 = new DensityCalculator; break; + m_p2 = new DensityCalculator; + 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 diff --git a/Cantera/src/ChemEquil.h b/Cantera/src/ChemEquil.h index bcff25c07..b59a5809c 100755 --- a/Cantera/src/ChemEquil.h +++ b/Cantera/src/ChemEquil.h @@ -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); - } - } diff --git a/Cantera/src/MultiPhase.cpp b/Cantera/src/MultiPhase.cpp index bd747ffb9..373e03c30 100644 --- a/Cantera/src/MultiPhase.cpp +++ b/Cantera/src/MultiPhase.cpp @@ -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; } diff --git a/Cantera/src/MultiPhase.h b/Cantera/src/MultiPhase.h index e2caf96b9..fae29569d 100644 --- a/Cantera/src/MultiPhase.h +++ b/Cantera/src/MultiPhase.h @@ -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]. diff --git a/Cantera/src/MultiPhaseEquil.cpp b/Cantera/src/MultiPhaseEquil.cpp index dd3cf029f..e99dbd581 100644 --- a/Cantera/src/MultiPhaseEquil.cpp +++ b/Cantera/src/MultiPhaseEquil.cpp @@ -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; } diff --git a/Cantera/src/MultiPhaseEquil.h b/Cantera/src/MultiPhaseEquil.h index acfb3523b..5335c600d 100644 --- a/Cantera/src/MultiPhaseEquil.h +++ b/Cantera/src/MultiPhaseEquil.h @@ -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); - } } diff --git a/Cantera/src/equil.h b/Cantera/src/equil.h new file mode 100644 index 000000000..43867adb3 --- /dev/null +++ b/Cantera/src/equil.h @@ -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 diff --git a/Cantera/src/equilibrate.cpp b/Cantera/src/equilibrate.cpp new file mode 100644 index 000000000..5cc2efaac --- /dev/null +++ b/Cantera/src/equilibrate.cpp @@ -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"); + } +}