From 240c93ff4f386d0b0853a9a9882681fcaeeec1a4 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Fri, 20 Oct 2006 21:34:00 +0000 Subject: [PATCH] Reworking of the solver algorithm within ChemEquil: Much of the algorithm for relaxing the element potentials to equilibrium employs the Brinkley algorithm (const T, P) (p. 129 Smith and Missen). This algorithm does better at relaxing the system towards steady state from a bad initial guess. The old algorithm now does the last little bit and it does the calculation for fixed pairs other than const T, P, using the Brinkley algorithm's initial guess. Modified the damping algorithm in the old method. It now doesn't create false convergence failures when it is close to the end of the calculation. The new solver has been tested against a matrix of initial conditions for the gri30 mechanism and a 200 species jp8 mechanism. Yet to do: This routine should work but has not been checked out against: phases with electric charge nonideal phases. This routine should probably not yet work with molality based phases. More work may have to be done with the Brinkley algorithm for problems which suffer from stoichiometric constraints. A full Gauss-Jordon elimination step should probably be added to flag cases where there are effective stoichiometric constraints. --- Cantera/src/ChemEquil.cpp | 2232 +++++++++++++++++++++++++------------ Cantera/src/ChemEquil.h | 267 +++-- 2 files changed, 1654 insertions(+), 845 deletions(-) diff --git a/Cantera/src/ChemEquil.cpp b/Cantera/src/ChemEquil.cpp index 9d759c0c2..b61b8de02 100755 --- a/Cantera/src/ChemEquil.cpp +++ b/Cantera/src/ChemEquil.cpp @@ -5,6 +5,9 @@ * Chemical equilibrium. Implementation file for class * ChemEquil. * + * + * $Id$ + * * Copyright 2001 California Institute of Technology * */ @@ -27,807 +30,1580 @@ using namespace std; #include "stringUtils.h" #include "MultiPhase.h" +#ifdef DEBUG_HKM_EPEQUIL +#include "stdio.h" +int debug_prnt_lvl = 0; +#endif + namespace Cantera { - /// map property strings to integers - int _equilflag(const char* xy) { - string flag = string(xy); - if (flag == "TP") return TP; - else if (flag == "TV") return TV; - else if (flag == "HP") return HP; - else if (flag == "UV") return UV; - else if (flag == "SP") return SP; - else if (flag == "SV") return SV; - else if (flag == "UP") return UP; - else throw CanteraError("_equilflag","unknown property pair "+flag); - } + /// map property strings to integers + int _equilflag(const char* xy) { + string flag = string(xy); + if (flag == "TP") return TP; + else if (flag == "TV") return TV; + else if (flag == "HP") return HP; + else if (flag == "UV") return UV; + else if (flag == "SP") return SP; + else if (flag == "SV") return SV; + else if (flag == "UP") return UP; + else throw CanteraError("_equilflag","unknown property pair "+flag); + } - //----------------------------------------------------------- - // construction / destruction - //----------------------------------------------------------- + //----------------------------------------------------------- + // construction / destruction + //----------------------------------------------------------- - /// Default Constructor. - ChemEquil::ChemEquil() : m_skip(-1), m_p1(0), m_p2(0), m_p0(OneAtm), m_eloc(-1), - m_abscharge(Tiny) - {} + /// Default Constructor. + ChemEquil::ChemEquil() : m_skip(-1), m_p1(0), m_p2(0), m_elementTotalSum(1.0), + m_p0(OneAtm), m_eloc(-1), + m_abscharge(Tiny), + m_elemFracCutoff(1.0E-100), + m_doResPerturb(false) + {} - - /// Destructor - ChemEquil::~ChemEquil(){ - delete m_p1; - delete m_p2; + /// Destructor + ChemEquil::~ChemEquil(){ + if (m_p1) { + delete m_p1; } + if (m_p2) { + delete m_p2; + } + } + /** + * 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. + // Note: the use of two pointers is a historical artifact. + m_thermo = &s; + m_phase = &s; - /** - * 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. - // Note: the use of two pointers is a historical artifact. - m_thermo = &s; - m_phase = &s; - - m_p0 = s.refPressure(); - m_kk = m_phase->nSpecies(); - m_mm = m_phase->nElements(); - if (m_kk < m_mm) { - throw CanteraError("ChemEquil::initialize", - "number of species cannot be less than the number of elements."); - } + m_p0 = s.refPressure(); + m_kk = m_phase->nSpecies(); + m_mm = m_phase->nElements(); + if (m_kk < m_mm) { + throw CanteraError("ChemEquil::initialize", + "number of species cannot be less than the number of elements."); + } - // allocate space in internal work arrays - m_molefractions.resize(m_kk); - m_lambda.resize(m_mm, -100.0); - m_elementmolefracs.resize(m_mm); - m_comp.resize(m_mm * m_kk); - m_jwork1.resize(m_mm+2); - m_jwork2.resize(m_mm+2); - m_startSoln.resize(m_mm+1); - m_grt.resize(m_kk); - m_mu_RT.resize(m_kk); - m_component.resize(m_mm,-2); + // allocate space in internal work arrays within the ChemEquil object + m_molefractions.resize(m_kk); + m_lambda.resize(m_mm, -100.0); + m_elementmolefracs.resize(m_mm); + m_comp.resize(m_mm * m_kk); + m_jwork1.resize(m_mm+2); + m_jwork2.resize(m_mm+2); + 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; - doublereal na, ewt; - for (m = 0; m < m_mm; m++) { - for (k = 0; k < m_kk; k++) { - na = m_phase->nAtoms(k,m); + // set up elemental composition matrix + int m, k, mneg = -1; + doublereal na, ewt; + for (m = 0; m < m_mm; m++) { + 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, where the 'element' is an - // electron - if (na < 0.0) { + // 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 - // for some element other than this one, throw - // an exception - if (mneg >= 0 && mneg != m) - throw CanteraError("ChemEquil::initialize", - "negative atom numbers allowed for only one element"); - mneg = m; - ewt = m_phase->atomicWeight(m); + // if negative atom numbers have already been specified + // for some element other than this one, throw + // an exception + if (mneg >= 0 && mneg != m) + throw CanteraError("ChemEquil::initialize", + "negative atom numbers allowed for only one element"); + mneg = m; + ewt = m_phase->atomicWeight(m); - // the element should be an electron... if it isn't - // print a warning. - if (ewt > 1.0e-3) - writelog(string("WARNING: species " - +m_phase->speciesName(k) - +" has "+fp2str(m_phase->nAtoms(k,m)) - +" atoms of element " - +m_phase->elementName(m)+ - ", but this element is not an electron.\n")); - } - } - } - m_eloc = mneg; - - // set up the elemental composition matrix - for (k = 0; k < m_kk; k++) { - for (m = 0; m < m_mm; m++) { - m_comp[k*m_mm + m] = m_phase->nAtoms(k,m); - } - } + // the element should be an electron... if it isn't + // print a warning. + if (ewt > 1.0e-3) + writelog(string("WARNING: species " + +m_phase->speciesName(k) + +" has "+fp2str(m_phase->nAtoms(k,m)) + +" atoms of element " + +m_phase->elementName(m)+ + ", but this element is not an electron.\n")); + } + } } + m_eloc = mneg; - - /** - * Set mixture to an equilibrium state consistent with specified - * element potentials and temperature. - * - * @param lambda_RT vector of non-dimensional element potentials - * \f[ \lambda_m/RT \f]. - * @param t temperature in K. - * - */ - void ChemEquil::setToEquilState(thermo_t& s, - const vector_fp& lambda_RT, doublereal t) - { - // 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++) - m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m); - - // set the temperature - s.setTemperature(t); - - // call the phase-specific method to set the phase to the - // equilibrium state with the specified species chemical - // potentials. - s.setToEquilState(DATA_PTR(m_mu_RT)); - update(s); + // set up the elemental composition matrix + for (k = 0; k < m_kk; k++) { + for (m = 0; m < m_mm; m++) { + m_comp[k*m_mm + m] = m_phase->nAtoms(k,m); + } } + } - /** - * update internally stored state information. - * @todo argument not used. - */ - void ChemEquil::update(const thermo_t& s) { + /** + * Set mixture to an equilibrium state consistent with specified + * element potentials and temperature. + * + * @param lambda_RT vector of non-dimensional element potentials + * \f[ \lambda_m/RT \f]. + * @param t temperature in K. + * + */ + void ChemEquil::setToEquilState(thermo_t& s, + const vector_fp& lambda_RT, doublereal t) + { + // 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++) + m_mu_RT[k] += lambda_RT[m]*nAtoms(k,m); - // get the mole fractions, temperature, and density - m_phase->getMoleFractions(DATA_PTR(m_molefractions)); - m_temp = m_phase->temperature(); - m_dens = m_phase->density(); + // Set the temperature + s.setTemperature(t); - // compute the elemental mole fractions - doublereal sum = 0.0; - int m, k; - for (m = 0; m < m_mm; m++) { - m_elementmolefracs[m] = 0.0; - for (k = 0; k < m_kk; k++) { - m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k]; - if (m_molefractions[k] < 0.0) { - throw CanteraError("update", - "negative mole fraction for "+m_phase->speciesName(k)+ - ": "+fp2str(m_molefractions[k])); - } - } - sum += m_elementmolefracs[m]; - } + // Call the phase-specific method to set the phase to the + // equilibrium state with the specified species chemical + // potentials. + s.setToEquilState(DATA_PTR(m_mu_RT)); + update(s); + } - // normalize the element mole fractions - for (m = 0; m < m_mm; m++) m_elementmolefracs[m] /= sum; + + /** + * update internally stored state information. + * @todo argument not used. + */ + void ChemEquil::update(const thermo_t& s) { + + // get the mole fractions, temperature, and density + m_phase->getMoleFractions(DATA_PTR(m_molefractions)); + m_temp = m_phase->temperature(); + m_dens = m_phase->density(); + + // compute the elemental mole fractions + double sum = 0.0; + int m, k; + for (m = 0; m < m_mm; m++) { + m_elementmolefracs[m] = 0.0; + for (k = 0; k < m_kk; k++) { + m_elementmolefracs[m] += nAtoms(k,m) * m_molefractions[k]; + if (m_molefractions[k] < 0.0) { + throw CanteraError("update", + "negative mole fraction for "+m_phase->speciesName(k)+ + ": "+fp2str(m_molefractions[k])); + } + } + sum += m_elementmolefracs[m]; } + // Store the sum for later use + m_elementTotalSum = sum; + // normalize the element mole fractions + for (m = 0; m < m_mm; m++) m_elementmolefracs[m] /= sum; + } - /// 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(); + /// 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(); - // 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; - } - catch (CanteraError) { - delete e; - delete mp; - iok = -1; - } - endLogGroup(); - return iok; + // 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 the current values of the temp, density, and + * mole fraction, and element abundance vectors kept + * within tine ChemEquil object. + */ + update(s); + delete e; + delete mp; + iok = 0; } + catch (CanteraError) { + delete e; + delete mp; + iok = -1; + } + endLogGroup(); + return iok; + } - /** - * Generate a starting estimate for the element potentials. - */ - int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda) - { - int m, n; - 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()); + /** + * Generate a starting estimate for the element potentials. + */ + int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda) + { + int m, n; + 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()); - DenseMatrix aa(m_mm, m_mm, 0.0); - vector_fp b(m_mm, -999.0); + DenseMatrix aa(m_mm, m_mm, 0.0); + vector_fp b(m_mm, -999.0); - vector_fp mu_RT(m_kk, 0.0); + vector_fp mu_RT(m_kk, 0.0); - s.getChemPotentials(DATA_PTR(mu_RT)); - doublereal rrt = 1.0/(GasConstant*m_phase->temperature()); - scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt); - for (m = 0; m < m_mm; m++) { - for (n = 0; n < m_mm; n++) { - aa(m,n) = nAtoms(m_component[m], n); - } - b[m] = mu_RT[m_component[m]]; - } + vector_fp xMF_est(m_kk, 0.0); - int info; - try { - info = solve(aa, DATA_PTR(b)); - } - catch (CanteraError) { - addLogEntry("failed to estimate initial element potentials."); - info = -2; - } + s.getMoleFractions(DATA_PTR(xMF_est)); + for (n = 0; n < s.nSpecies(); n++) { + if (xMF_est[n] < 1.0E-20) { + xMF_est[n] = 1.0E-20; + } + } + s.setMoleFractions(DATA_PTR(xMF_est)); + - if (info == 0) { - for (m = 0; m < m_mm; m++) { - lambda[m] = b[m]; - addLogEntry(m_phase->elementName(m),b[m]); - } - } - endLogGroup(); - return info; + + s.getChemPotentials(DATA_PTR(mu_RT)); + doublereal rrt = 1.0/(GasConstant*m_phase->temperature()); + scale(mu_RT.begin(), mu_RT.end(), mu_RT.begin(), rrt); + +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + for (m = 0; m < m_mm; m++) { + int isp = m_component[m]; + string nnn = s.speciesName(isp); + printf("isp = %d, %s\n", isp, nnn.c_str()); + } + 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"); + + s.getMoleFractions(DATA_PTR(xMF_est)); + + for (n = 0; n < s.nSpecies(); n++) { + string nnn = s.speciesName(n); + printf("%10d %15s %10.5g %10.5g\n", + n, nnn.c_str(), xMF_est[n], mu_RT[n]); + } + } +#endif + + for (m = 0; m < m_mm; m++) { + for (n = 0; n < m_mm; n++) { + aa(m,n) = nAtoms(m_component[m], n); + } + b[m] = mu_RT[m_component[m]]; } - - - /** - * Equilibrate a phase, holding the elemental composition fixed - * at the initial vaollue. - */ - int ChemEquil::equilibrate(thermo_t& s, const char* XY) { - vector_fp emol(s.nElements()); - initialize(s); - update(s); - copy(m_elementmolefracs.begin(), m_elementmolefracs.end(), - emol.begin()); - return equilibrate(s, XY, emol); + int info; + try { + info = solve(aa, DATA_PTR(b)); + } + catch (CanteraError) { + addLogEntry("failed to estimate initial element potentials."); + info = -2; + } + for (m = 0; m < m_mm; m++) { + lambda[m] = b[m]; + } + if (info == 0) { + for (m = 0; m < m_mm; m++) { + addLogEntry(m_phase->elementName(m),lambda[m]); + } } +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + for (m = 0; m < m_mm; m++) { + int isp = m_component[m]; + double tmp = 0.0; + for (n = 0; n < m_mm; n++) { + tmp += nAtoms(isp, n) * lambda[n]; + } + printf("%3d %10.5g %10.5g %10.5g\n", + m, mu_RT[isp], tmp, tmp - mu_RT[isp]); - /** - * compute the equilibrium composition for 2 specified - * properties and specified element moles. + } + } +#endif + endLogGroup(); + return info; + } + + + /** + * Equilibrate a phase, holding the elemental composition fixed + * at the initial value found within the ThermoPhase object. + * + * The value of 2 specified properties are obtained by querying the + * ThermoPhase object. The properties must be already contained + * within the current thermodynamic state of the system. + */ + int ChemEquil::equilibrate(thermo_t& s, const char* XY, + bool useThermoPhaseElementPotentials) { + vector_fp emol(s.nElements()); + initialize(s); + update(s); + copy(m_elementmolefracs.begin(), m_elementmolefracs.end(), + emol.begin()); + return equilibrate(s, XY, emol, useThermoPhaseElementPotentials); + } + + + /** + * Compute the equilibrium composition for 2 specified + * properties and the specified element moles. + * + * elMoles = specified vector of element abundances. + * + * The 2 specified properties are obtained by querying the + * ThermoPhase object. The properties must be already contained + * within the current thermodynamic state of the system. + * + * Return variable: + * Successful returns are indicated by a return value of 0. + * Unsuccessful returns are indicated by a return value of -1 for + * lack of convergence or -3 for a singular jacobian. + */ + int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp& elMoles, + bool useThermoPhaseElementPotentials) + { + doublereal xval, yval, tmp; + int fail = 0; + int m, n; + + if (m_p1) delete m_p1; + if (m_p2) delete m_p2; + bool tempFixed = true; + int XY = _equilflag(XYstr); + + vector_fp state; + s.saveState(state); + +#ifdef DEBUG_HKM_EPEQUIL + const vector& eNames = s.elementNames(); +#endif + beginLogGroup("ChemEquil::equilibrate"); + initialize(s); + update(s); + switch (XY) { + case TP: case PT: + m_p1 = new TemperatureCalculator; + m_p2 = new PressureCalculator; + break; + case HP: case PH: + tempFixed = false; + m_p1 = new EnthalpyCalculator; + m_p2 = new PressureCalculator; + break; + case SP: case PS: + tempFixed = false; + m_p1 = new EntropyCalculator; + m_p2 = new PressureCalculator; + break; + case SV: case VS: + tempFixed = false; + m_p1 = new EntropyCalculator; + m_p2 = new DensityCalculator; + break; + case TV: case VT: + m_p1 = new TemperatureCalculator; + m_p2 = new DensityCalculator; + break; + case UV: case VU: + tempFixed = false; + m_p1 = new IntEnergyCalculator; + m_p2 = new DensityCalculator; + break; + default: + endLogGroup("ChemEquil::equilibrate"); + 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) { + double tfixed = s.temperature(); + if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) { + endLogGroup("ChemEquil::equilibrate"); + throw CanteraError("ChemEquil","Specified temperature (" + +fp2str(m_thermo->temperature())+" K) outside " + "valid range of "+fp2str(m_thermo->minTemp())+" K to " + +fp2str(m_thermo->maxTemp())+" K\n"); + } + } + + /* + * Before we do anything to change the ThermoPhase object, + * we calculate and store the two specified thermodynamic + * properties that we are after. */ - int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, vector_fp& elMoles) - { - doublereal xval, yval; - int fail = 0; + xval = m_p1->value(s); + yval = m_p2->value(s); - delete m_p1; - delete m_p2; - bool tempFixed = true; - int XY = _equilflag(XYstr); + int mm = m_mm; + int nvar = mm + 1; + DenseMatrix jac(nvar, nvar); // jacobian + vector_fp x(nvar, -102.0); // solution vector + vector_fp res_trial(nvar, 0.0); // residual - vector_fp state; - s.saveState(state); + /* + * Replace one of the element abundance fraction equations + * with the specified property calculation. + * + * We choose the equation of the element with the highest element + * abundance. + */ + tmp = -1.0; + for (m = 0; m < mm; m++) { + if (elMoles[m] > tmp ) { + m_skip = m; + tmp = elMoles[m]; + } + } + if (tmp <= 0.0) { + throw CanteraError("ChemEquil", + "Element Abundance Vector is zeroed"); + } - beginLogGroup("ChemEquil::equilibrate"); - initialize(s); - update(s); - switch (XY) { - case TP: case PT: - m_p1 = new TemperatureCalculator; - m_p2 = new PressureCalculator; - break; - case HP: case PH: - tempFixed = false; - m_p1 = new EnthalpyCalculator; - m_p2 = new PressureCalculator; - break; - case SP: case PS: - tempFixed = false; - m_p1 = new EntropyCalculator; - m_p2 = new PressureCalculator; - break; - case SV: case VS: - tempFixed = false; - m_p1 = new EntropyCalculator; - m_p2 = new DensityCalculator; - break; - case TV: case VT: - m_p1 = new TemperatureCalculator; - m_p2 = new DensityCalculator; - break; - case UV: case VU: - tempFixed = false; - m_p1 = new IntEnergyCalculator; - m_p2 = new DensityCalculator; - break; - default: - endLogGroup(); - throw CanteraError("equilibrate","illegal property pair."); - } + // 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) + 1.0E-32; + } + m_phase->setMoleFractions(DATA_PTR(xmm)); - 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)); + /* + * Update the internally storred values of m_temp, + * m_dens, and the element mole fractions. + */ + update(s); - // If the temperature is one of the specified variables, and - // it is outside the valid range, throw an exception. - if (tempFixed) { - double tfixed = s.temperature(); - if (tfixed > s.maxTemp() + 1.0 || tfixed < s.minTemp() - 1.0) { - endLogGroup(); - throw CanteraError("ChemEquil","Specified temperature (" - +fp2str(m_thermo->temperature())+" K) outside " - "valid range of "+fp2str(m_thermo->minTemp())+" K to " - +fp2str(m_thermo->maxTemp())+" K\n"); - } - } - xval = m_p1->value(s); - yval = m_p2->value(s); - int mm = m_mm; - - int m; - int nvar = mm + 1; - - DenseMatrix jac(nvar, nvar); // jacobian - vector_fp x(nvar, -102.0); // solution vector - - vector_fp res_trial(nvar); + // loop to estimate T + if (!tempFixed) { - for (m = 0; m < mm; m++) { - if (m_skip < 0 && elMoles[m] > 0.0 ) m_skip = m; - } + beginLogGroup("Initial T Estimate"); - // 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(DATA_PTR(xmm)); + doublereal tmax = m_thermo->maxTemp(); + doublereal tmin = m_thermo->minTemp(); + doublereal slope, phigh, plow, pval, dt; - update(s); + // 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); + phigh = m_p1->value(s); + + m_phase->setTemperature(tmin); + setInitialMoles(s); + plow = m_p1->value(s); + + // start with T at the midpoint of the range + doublereal t0 = 0.5*(tmin + tmax); + m_phase->setTemperature(t0); + + // loop up to 5 times + for (int it = 0; it < 5; it++) { + + // set the composition and get p1 + setInitialMoles(s); + pval = m_p1->value(s); - // loop to estimate T - if (!tempFixed) { + // If this value of p1 is greater than the specified + // property value, then the current temperature is too + // high. Use it as the new upper bound. Otherwise, it + // is too low, so use it as the new lower bound. + if (pval > xval) { + tmax = t0; + phigh = pval; + } + else { + tmin = t0; + plow = pval; + } - beginLogGroup("Initial T Estimate"); + // Determine the new T estimate by linearly intepolation + // between the upper and lower bounds + slope = (phigh - plow)/(tmax - tmin); + dt = (xval - plow)/slope; - doublereal tmax = m_thermo->maxTemp(); - doublereal tmin = m_thermo->minTemp(); - doublereal slope, phigh, plow, pval, dt; + // If within 100 K, terminate the search + if (fabs(dt) < 100.0) break; - // 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. + // update the T estimate + t0 = tmin + dt; + addLogEntry("new T estimate", t0); - m_phase->setTemperature(tmax); - setInitialMoles(s); - phigh = m_p1->value(s); + m_phase->setTemperature(t0); + } + endLogGroup("Initial T Estimate"); // initial T estimate + } - m_phase->setTemperature(tmin); - setInitialMoles(s); - plow = m_p1->value(s); - - // start with T at the midpoint of the range - doublereal t0 = 0.5*(tmin + tmax); - m_phase->setTemperature(t0); - - // loop up to 5 times - for (int it = 0; it < 5; it++) { - - // set the composition and get p1 - setInitialMoles(s); - pval = m_p1->value(s); - - - // If this value of p1 is greater than the specified - // property value, then the current temperature is too - // high. Use it as the new upper bound. Otherwise, it - // is too low, so use it as the new lower bound. - if (pval > xval) { - tmax = t0; - phigh = pval; - } - else { - tmin = t0; - plow = pval; - } - - // Determine the new T estimate by linearly intepolation - // between the upper and lower bounds - slope = (phigh - plow)/(tmax - tmin); - dt = (xval - plow)/slope; - - // If within 100 K, terminate the search - if (fabs(dt) < 100.0) break; - - // 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); - 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()); - - vector_fp above(nvar); - vector_fp below(nvar); - - for (m = 0; m < mm; m++) { - above[m] = 200.0; - below[m] = -2000.0; - if (elMoles[m] < Cutoff && m != m_eloc) x[m] = -1000.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 oldresid(nvar, 0.0); - doublereal f, oldf; - - int iter = 0; - int info=0; - doublereal fctr = 1.0, newval; - - goto converge; - - next: - - - // if the problem involves charged species, then the - // "electron" element equation is a charge balance. Compute - // the sum of the absolute values of the charge to use as the - // normalizing factor. - if (m_eloc >= 0) { - m_abscharge = 0.0; - int k; - for (k = 0; k < m_kk; k++) - m_abscharge += fabs(m_phase->charge(k)*m_molefractions[k]); - } - - - 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, xval, yval); - - f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); - addLogEntry("Residual norm", f); - - equilJacobian(s, x, elMoles, jac, xval, yval); - - // compute grad f = F*J - jac.leftMult(DATA_PTR(res_trial), DATA_PTR(grad)); - copy(x.begin(), x.end(), oldx.begin()); - - oldf = f; - scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0); - try { - info = solve(jac, DATA_PTR(res_trial)); - } - catch (CanteraError) { - addLogEntry("Jacobian is singular."); - endLogGroup(); // iteration - endLogGroup(); // equilibrate - s.restoreState(state); - - throw CanteraError("equilibrate", - "Jacobian is singular. \nTry adding more species, " - "changing the elemental composition slightly, \nor removing " - "unused elements."); - return -3; - } - - // find the factor by which the Newton step can be multiplied - // to keep the solution within bounds. - fctr = 1.0; - for (m = 0; m < nvar; m++) { - newval = x[m] + res_trial[m]; - if (newval > above[m]) { - fctr = fmaxx( 0.0, fminn( fctr, - 0.8*(above[m] - x[m])/(newval - x[m]))); - } - else if (newval < below[m]) { - 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 scaling factor - scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr); - - if (!dampStep(s, oldx, oldf, grad, res_trial, - x, f, elMoles , xval, yval)) - { - fail++; - if (fail > 3) { - addLogEntry("dampStep","Failed 3 times. Giving up."); - endLogGroup(); // iteration - endLogGroup(); // equilibrate - s.restoreState(state); - throw CanteraError("equilibrate", - "Cannot find an acceptable Newton damping coefficient."); - return -4; - } - } - else fail = 0; - -converge: - - // check for convergence. - - 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; - 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 ) { - writelog("Warning: Temperature (" - +fp2str(m_thermo->temperature())+" K) outside " - "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."); - return -1; - } - goto next; - } - - - 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, double xval, double yval ) - { - int nvar = x.size(); - - double slope; - double f2 = 0.0; - double oldf2 = 0.0; - double alpha = 1.e-4; - double tmpdamp = 0.0; - double rhs1; - double rhs2; - double damp = 1.0; - double damp2=0.0; - double a; - double bb; - double disc; - double minDamp = 0.0; - double xTol = 1.e-7; - - vector_fp res_new(nvar); // fix - - //slope = grad * step; - slope = dot(grad.begin(), grad.end(), step.begin()); - double temp, test = 0.0; + + setInitialMoles(s); - for (int i=0; i test) test = temp; - } - minDamp = xTol/test; - - retry: - - x = step; - scale(x, damp); - add_each(x, oldx); - - - 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) - { - return 0; // check that this is not a spurious min of f - } - else if (f <= oldf + alpha * damp * slope) - { - return 1; // good damping coefficient - } - else - { - if (damp == 1.0) // first time - { - tmpdamp = -slope/(2.0*(f - oldf - slope)); - } - else - { - rhs1 = f - oldf - damp*slope; - rhs2 = f2 - oldf2 - damp2*slope; - a = (rhs1/(damp*damp) - rhs2/(damp2*damp2))/(damp - damp2); - bb = (-damp2*rhs1/(damp*damp) + damp*rhs2/(damp2*damp2)) - /(damp - damp2); - - if (a == 0.0) - tmpdamp = -slope/(2.0*bb); - else - { - disc = bb*bb - 3.0*a*slope; - if (disc < 0.0) - tmpdamp = -slope/(2.0*bb); - else - tmpdamp = (-bb +sqrt(disc))/(3.0*a); - } - if (tmpdamp > 0.5*damp) tmpdamp = 0.5*damp; - } - - damp2 = damp; - f2 = f; - oldf2 = oldf; - damp = fmaxx(tmpdamp, 0.1*damp); - goto retry; - } - } - - - /** - * evaluates the residual vector F, of length mm + /* + * If requested, get the initial estimate for the + * chemical potentials from the ThermoPhase object + * itself. Or else, create our own estimate. */ - void ChemEquil::equilResidual(thermo_t& mix, const vector_fp& x, - const vector_fp& elmtotal, vector_fp& resid, - doublereal xval, doublereal yval) - { - beginLogGroup("ChemEquil::equilResidual"); - int n; - doublereal xx, yy; - doublereal temp = exp(x[m_mm]); - setToEquilState(mix, x, temp); + if (useThermoPhaseElementPotentials) { + bool haveEm = s.getElementPotentials(DATA_PTR(x)); + if (haveEm) { + doublereal rt = GasConstant * m_thermo->temperature(); + for (m = 0; m < m_mm; m++) { + x[m] /= rt; + } + } else { + estimateElementPotentials(s, x); + } + } else { + /* + * Calculate initial estimates of the element potentials. + * This algorithm uese the MultiPhaseEquil object's + * initialization capabilities to calculate an initial + * estimate of the mole fractions for a set of linearly + * independent component species. Then, the element + * potentials are solved for based on the chemical + * potentials of the component species. + */ + estimateElementPotentials(s, x); + } - // residuals are the total element moles - vector_fp& elm = m_elementmolefracs; - for (n=0; n < m_mm; n++) - { - // drive element potential for absent elements to -1000 - if (elmtotal[n] < Cutoff && n != m_eloc) - resid[n] = x[n] + 1000.0; - else - 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; - for (int k = 0; k < m_kk; k++) { - chrg = m_molefractions[k]*m_phase->charge(k); - 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; - string xstr = fp2str(xx)+" ("+fp2str(xval)+")"; - addLogEntry(m_p1->symbol(), xstr); - string ystr = fp2str(yy)+" ("+fp2str(yval)+")"; - addLogEntry(m_p2->symbol(), ystr); - endLogGroup(); + /* + * Do a better estimate of the element potentials. + * We have found that the current estimate may not be good + * enough to avoid drastic numerical issues associated with + * the use of a numerically generated jacobian. + * + * The Brinkley algorithm assumes a constant T, P system + * and uses a linearized analytical Jacobian that turns out + * to be very stable. + */ + int info = estimateEP_Brinkley(s, x, elMoles); + if (info != 0) { + if (info == 1) { + addLogEntry("estimateEP_Brinkley didn't converge in given max interations"); + } else if (info == -3) { + addLogEntry("estimateEP_Brinkley had a singular Jacobian. Continuing anyway"); + } + } else { + setToEquilState(s, x, m_phase->temperature()); + // Tempting -> However, nonideal is a problem. Turn on if not worried + // about nonideality and you are having problems with the main + // algorithm. + //if (XY == TP) { + // endLogGroup("ChemEquil::equilibrate"); + // return 0; + //} + } + + /* + * Install the log(temp) into the last solution unknown + * slot. + */ + x[m_mm] = log(m_phase->temperature()); + + /* + * Setting the max and min values for x[]. Also, if element + * abundance vector is zero, setting x[] to -1000. This + * effectively zeroes out all species containing that element. + */ + vector_fp above(nvar); + vector_fp below(nvar); + for (m = 0; m < mm; m++) { + above[m] = 200.0; + below[m] = -2000.0; + if (elMoles[m] < m_elemFracCutoff && m != m_eloc) x[m] = -1000.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 oldresid(nvar, 0.0); + doublereal f, oldf; + + int iter = 0; + doublereal fctr = 1.0, newval; + + goto converge; + next: + // If the problem involves charged species, then the + // "electron" element equation is a charge balance. Compute + // the sum of the absolute values of the charge to use as the + // normalizing factor. + if (m_eloc >= 0) { + m_abscharge = 0.0; + int k; + for (k = 0; k < m_kk; k++) + m_abscharge += fabs(m_phase->charge(k)*m_molefractions[k]); } - //-------------------- Jacobian evaluation --------------------------- + iter++; + if (iter > 1) endLogGroup("Iteration "+int2str(iter-1)); // iteration + beginLogGroup("Iteration "+int2str(iter)); - void ChemEquil::equilJacobian(thermo_t& mix, vector_fp& x, - const vector_fp& elmols, DenseMatrix& jac, - doublereal xval, doublereal yval) - { - beginLogGroup("equilJacobian",0); + // compute the residual and the jacobian using the current + // solution vector + equilResidual(s, x, elMoles, res_trial, xval, yval); + f = 0.5*dot(res_trial.begin(), res_trial.end(), res_trial.begin()); + addLogEntry("Residual norm", f); - int len = x.size(); - vector_fp& r0 = m_jwork1; - vector_fp& r1 = m_jwork2; - r0.resize(len); - r1.resize(len); - - int n, m; - doublereal rdx, dx, xsave; - doublereal atol = 1.e-10; + // Compute the Jacobian matrix + equilJacobian(s, x, elMoles, jac, xval, yval); - equilResidual(mix, x, elmols, r0, xval, yval); - - for (n = 0; n < len; n++) - { - // perturb x(n) - - xsave = x[n]; - dx = atol; - x[n] = xsave + dx; - dx = x[n] - xsave; - rdx = 1.0/dx; +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("Jacobian matrix %d:\n", iter); + for (m = 0; m <= m_mm; m++) { + printf(" [ "); + for (n = 0; n <= m_mm; n++) { + printf("%10.5g ", jac(m,n)); + } + printf(" ]"); + char xName[32]; + if (m < m_mm) { + string nnn = eNames[m]; + sprintf(xName, "x_%-10s", nnn.c_str()); + } else { + sprintf(xName, "x_XX"); + } + if (m_eloc == m) { + sprintf(xName, "x_ELOC"); + } + if (m == m_skip) { + sprintf(xName, "x_YY"); + } + printf("%-12s", xName); + printf(" = - (%10.5g)\n", res_trial[m]); + } + } +#endif - // calculate perturbed residual + + // compute grad f = F*J + jac.leftMult(DATA_PTR(res_trial), DATA_PTR(grad)); + copy(x.begin(), x.end(), oldx.begin()); + oldf = f; + scale(res_trial.begin(), res_trial.end(), res_trial.begin(), -1.0); + + /* + * Solve the system + */ + try { + info = solve(jac, DATA_PTR(res_trial)); + } + catch (CanteraError) { + addLogEntry("Jacobian is singular."); + endLogGroup(); // iteration + endLogGroup(); // equilibrate + s.restoreState(state); + + throw CanteraError("equilibrate", + "Jacobian is singular. \nTry adding more species, " + "changing the elemental composition slightly, \nor removing " + "unused elements."); + return -3; + } + + // find the factor by which the Newton step can be multiplied + // to keep the solution within bounds. + fctr = 1.0; + for (m = 0; m < nvar; m++) { + newval = x[m] + res_trial[m]; + if (newval > above[m]) { + fctr = fmaxx( 0.0, fminn( fctr, + 0.8*(above[m] - x[m])/(newval - x[m]))); + } + else if (newval < below[m]) { + 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 scaling factor + scale(res_trial.begin(), res_trial.end(), res_trial.begin(), fctr); + + if (!dampStep(s, oldx, oldf, grad, res_trial, + x, f, elMoles , xval, yval)) + { + fail++; + if (fail > 3) { + addLogEntry("dampStep","Failed 3 times. Giving up."); + endLogGroup(); // iteration + endLogGroup(); // equilibrate + s.restoreState(state); + throw CanteraError("equilibrate", + "Cannot find an acceptable Newton damping coefficient."); + return -4; + } + } + else fail = 0; + + converge: + + // check for convergence. + 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; + doublereal rmax = 0.0; + bool passThis = true; + for (m = 0; m < nvar; m++) { + double tval = options.relTolerance; + if (m < mm) { + tval = elMoles[m] * options.relTolerance + options.absElemTol; + } + if (fabs(res_trial[m]) > tval) { + passThis = false; + } + } + if (iter > 0 && passThis + && fabs(deltax) < options.relTolerance + && fabs(deltay) < options.relTolerance) { + options.iterations = iter; + endLogGroup("Iteration "+int2str(iter)); // iteration + 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])); + } + /* + * Save the calculated and converged element potentials + * to the original ThermoPhase object. + */ + s.setElementPotentials(m_lambda); + addLogEntry("Saving Element Potentials to ThermoPhase Object"); + endLogGroup("Element potentials"); + + if (m_thermo->temperature() > m_thermo->maxTemp() + 1.0 || + m_thermo->temperature() < m_thermo->minTemp() - 1.0 ) { + writelog("Warning: Temperature (" + +fp2str(m_thermo->temperature())+" K) outside " + "valid range of "+fp2str(m_thermo->minTemp())+" K to " + +fp2str(m_thermo->maxTemp())+" K\n"); + } + endLogGroup("Converged solution"); + endLogGroup("ChemEquil::equilibrate"); + return 0; + } + + // no convergence + + if (iter > options.maxIterations) { + addLogEntry("equilibrate","no convergence"); + endLogGroup("Iteration "+int2str(iter)); + endLogGroup("ChemEquil::equilibrate"); + s.restoreState(state); + throw CanteraError("equilibrate", + "no convergence in "+int2str(options.maxIterations) + +" iterations."); + return -1; + } + goto next; + } + + + /* + * dampStep: Come up with an acceptable step size. The original implementation + * employed a line search technique that enforced a reduction in the + * norm of the residual at every successful step. Unfortunately, + * this method created false convergence errors near the end of + * a significant number of steps, usually special conditions where + * there were stoichiometric constraints. + * + * This new method just does a delta damping approach, based on limiting + * the jump in the dimensionless element potentials. Mole fractions are + * limited to a factor of 2 jump in the values from this method. + * Near convergence, the delta damping gets out of the way. + */ + 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, double xval, double yval ) + { + int nvar = x.size(); + int m; + double damp; + + /* + * Carry out a delta damping approach on the dimensionless element potentials. + */ + damp = 1.0; + for (m = 0; m < m_mm; m++) { + if (step[m] > 0.75) { + damp = 0.75 /step[m]; + } + if (step[m] < -0.75) { + damp = -0.75 / step[m]; + } + } + + /* + * Update the solution unknown + */ + 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"); + for (m = 0; m < nvar; m++) { + printf(" %10.5g %10.5g %10.5g\n", x[m], oldx[m], step[m]); + } + } +#endif + return 1; + } + + + /** + * evaluates the residual vector F, of length mm + */ + void ChemEquil::equilResidual(thermo_t& mix, const vector_fp& x, + const vector_fp& elmFracGoal, vector_fp& resid, + doublereal xval, doublereal yval) + { + beginLogGroup("ChemEquil::equilResidual"); + int n; + doublereal xx, yy; + doublereal temp = exp(x[m_mm]); + setToEquilState(mix, x, temp); + + // residuals are the total element moles + vector_fp& elmFrac = m_elementmolefracs; + for (n = 0; n < m_mm; n++) + { + // drive element potential for absent elements to -1000 + if (elmFracGoal[n] < m_elemFracCutoff && n != m_eloc) + resid[n] = x[n] + 1000.0; + else { + /* + * Change the calculation for small element number, using + * L'Hopital's rule. + * The log formulation is unstable. + */ + if (elmFracGoal[n] < 1.0E-10 || elmFrac[n] < 1.0E-10) { + resid[n] = elmFracGoal[n] - elmFrac[n]; + } else { + resid[n] = log( (1.0 + elmFracGoal[n]) / (1.0 + elmFrac[n]) ); + } + } + addLogEntry(m_phase->elementName(n),fp2str(elmFrac[n])+" (" + +fp2str(elmFracGoal[n])+")"); + } + if (m_eloc >= 0) { + doublereal chrg, sumnet = 0.0, sumabs = 0.0; + for (int k = 0; k < m_kk; k++) { + chrg = m_molefractions[k]*m_phase->charge(k); + 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; + string xstr = fp2str(xx)+" ("+fp2str(xval)+")"; + addLogEntry(m_p1->symbol(), xstr); + string ystr = fp2str(yy)+" ("+fp2str(yval)+")"; + addLogEntry(m_p2->symbol(), ystr); + endLogGroup("ChemEquil::equilResidual"); + +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0 && !m_doResPerturb) { + printf("Residual: ElFracGoal ElFracCurrent Resid\n"); + for (n = 0; n < m_mm; n++) { + printf(" %14.9g %14.9g %10.5g\n", elmFracGoal[n], elmFrac[n], resid[n]); + } + } +#endif + } + + + //-------------------- Jacobian evaluation --------------------------- + + void ChemEquil::equilJacobian(thermo_t& s, vector_fp& x, + const vector_fp& elmols, DenseMatrix& jac, + doublereal xval, doublereal yval) + { + beginLogGroup("equilJacobian"); + int len = x.size(); + vector_fp& r0 = m_jwork1; + vector_fp& r1 = m_jwork2; + r0.resize(len); + r1.resize(len); + int n, m; + doublereal rdx, dx, xsave, dx2; + doublereal atol = 1.e-10; + + equilResidual(s, x, elmols, r0, xval, yval); + + m_doResPerturb = true; + for (n = 0; n < len; n++) { + xsave = x[n]; + dx = atol; + dx2 = fabs(xsave) * 1.0E-7; + if (dx2 > dx) dx = dx2; + x[n] = xsave + dx; + dx = x[n] - xsave; + rdx = 1.0/dx; + + // calculate perturbed residual - equilResidual(mix, x, elmols, r1, xval, yval); + equilResidual(s, x, elmols, r1, xval, yval); - // compute nth column of Jacobian + // compute nth column of Jacobian - for (m = 0; m < len; m++) { - jac(m, n) = (r1[m] - r0[m])*rdx; - } - x[n] = xsave; - } - endLogGroup(); + for (m = 0; m < len; m++) { + jac(m, n) = (r1[m] - r0[m])*rdx; + } + x[n] = xsave; + } + m_doResPerturb = false; + endLogGroup("equilJacobian"); + } + + /** + * Do a calculation of the element potentials using + * the Brinkley method, p. 129 Smith and Missen. + * + * We have found that the previous estimate may not be good + * enough to avoid drastic numerical issues associated with + * the use of a numerically generated jacobian used in the + * main algorithm. + * + * The Brinkley algorithm, here, assumes a constant T, P system + * and uses a linearized analytical Jacobian that turns out + * to be very stable even given bad initial guesses. + * + * The pressure and temperature to be used are in the + * ThermoPhase object input into the routine. + * + * The initial guess for the element potentials + * used by this routine is taken from the + * input vector, x. + * + * elMoles is the input element abundance vector to be matched. + * + * Nonideal phases are handled in principle. This is done by + * calculating the activity coefficients and adding them + * into the formula in the correct position. However, + * these are treated as a rhs contribution only. Therefore, + * convergence might be a problem. This has not been tested. + * Also molality based unit systems aren't handled. + * + * On return, int return value contains the success code: + * 0 - successful + * 1 - unsuccessful, max num iterations exceeded + * -3 - unsuccessful, singular jacobian + * + * NOTE: update for activity coefficients. + */ + int ChemEquil::estimateEP_Brinkley(thermo_t& s, vector_fp& x, + vector_fp& elMoles) { + /* + * Before we do anything, we will save the state of the solution. + * Then, if things go drastically wrong, we will restore the + * saved state. + */ + vector_fp state; + s.saveState(state); + double tmp, sum; + bool modifiedMatrix = false; + int neq = m_mm+1; + int retn = 1; + int m, n, k, info; + DenseMatrix a1(neq, neq, 0.0); + vector_fp b(neq, 0.0); + vector_fp muSS_RT(m_kk, 0.0); + vector_fp n_i(m_kk,0.0); + vector_fp n_i_calc(m_kk,0.0); + vector_fp actCoeff(m_kk, 1.0); + vector_fp muSS_RT_mod(m_kk, 0.0); + double beta = 1.0; + + s.getMoleFractions(DATA_PTR(n_i)); + + vector_fp x_old(m_mm+1, 0.0); + vector_fp resid(m_mm+1, 0.0); + vector_int lumpSum(m_mm+1, 0); + + /* + * Get the nondimensional Gibbs functions for the species + * at their standard states of solution at the current T and P + * of the solution. + */ + s.getGibbs_RT(DATA_PTR(muSS_RT)); + copy(muSS_RT.begin(), muSS_RT.end(), muSS_RT_mod.begin()); + + + vector_fp eMolesCalc(m_mm, 0.0); + vector_fp eMolesFix(m_mm, 0.0); + double elMolesTotal = 0.0; + for (m = 0; m < m_mm; m++) { + elMolesTotal += elMoles[m]; + for (k = 0; k < m_kk; k++) { + eMolesFix[m] += nAtoms(k,m) * n_i[k]; + } + } + + for (m = 0; m < m_mm; m++) { + if (x[m] > 50.0) { + x[m] = 50.; + } + if (elMoles[m] > 1.0E-70) { + if (x[m] < -100) { + x[m] = -100.; + } + } else { + if (x[m] < -1000.) { + x[m] = -1000.; + } + } } + double n_t = 0.0; + double sum2 = 0.0; + double nAtomsMax = 1.0; + s.setMoleFractions(DATA_PTR(x)); + s.getActivityCoefficients(DATA_PTR(actCoeff)); + for (k = 0; k < m_kk; k++) { + tmp = - (muSS_RT[k] + log(actCoeff[k])); + sum2 = 0.0; + for (m = 0; m < m_mm; m++) { + sum = nAtoms(k,m); + tmp += sum * x[m]; + sum2 += sum; + if (sum2 > nAtomsMax) { + nAtomsMax = sum2; + } + } + if (tmp > 100.) { + n_t += 2.8E43; + } else { + n_t += exp(tmp); + } + } + +#ifdef DEBUG_HKM_EPEQUIL + const vector& eNames = s.elementNames(); + if (debug_prnt_lvl > 0) { + printf("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"); + for (k = 0; k < m_kk; k++) { + string nnn = s.speciesName(k); + printf("%15s %13.5g %13.5g %13.5g\n", + nnn.c_str(), n_i[k], muSS_RT[k], actCoeff[k]); + } + printf("Initial n_t = %10.5g\n", n_t); + printf("Comparison of Goal Element Abundance with Initial Guess:\n"); + printf(" 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]); + } + } +#endif + for (m = 0; m < m_mm; m++) { + if (m != m_eloc) { + if (elMoles[m] <= options.absElemTol) { + x[m] = -200.; + } + } + } + /* + * ------------------------------------------------------------------- + * Main Loop. + */ + for (int iter = 0; iter < 2* options.maxIterations; iter++) { + /* + * Save the old solution + */ + for (m = 0; m < m_mm; m++) { + x_old[m] = x[m]; + } + x_old[m_mm] = n_t; + /* + * Calculate the mole numbers of species + */ +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("START ITERATION %d:\n", iter); + } +#endif + + double n_t_calc = 0.0; + s.setMoleFractions(DATA_PTR(x)); + s.getActivityCoefficients(DATA_PTR(actCoeff)); + sum2 = 0.0; + for (k = 0; k < m_kk; k++) { + tmp = - (muSS_RT[k] + log(actCoeff[k])); + for (m = 0; m < m_mm; m++) { + tmp += nAtoms(k,m) * x[m]; + } + if (tmp > 100.) tmp = 100.; + if (tmp < -300.) { + n_i_calc[k] = 0.0; + } else { + n_i_calc[k] = n_t * exp(tmp); + } + n_t_calc += n_i_calc[k]; +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + string nnn = s.speciesName(k); + printf("%15s: %10.5g (%10.5g)\n", nnn.c_str(), + n_i_calc[k], tmp); + } +#endif + } +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("%15s: %10.5g\n", "Total Molar Sum", n_t_calc); + } +#endif + for (m = 0; m < m_mm; m++) { + eMolesCalc[m] = 0.0; + for (k = 0; k < m_kk; k++) { + eMolesCalc[m] += nAtoms(k,m) * n_i_calc[k]; + } + } +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("(iter %d) element moles bal: Goal Calculated\n", iter); + 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]); + } + } +#endif + + bool normalStep = true; + /* + * Decide if we are to do a normal step or a modified step + */ + int iM = -1; + for (m = 0; m < m_mm; m++) { + if (elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + normalStep = false; + iM = m; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + normalStep = false; + iM = m; + } + } + } +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + if (!normalStep) { + printf(" NOTE: iter(%d) Doing an abnormal step due to row %d\n", iter, iM); + } + } +#endif + if (!normalStep) { + beta = 1.0; + resid[m_mm] = 0.0; + for (m = 0; m < m_mm; m++) { + resid[m] = 0.0; + if (elMoles[m] > 0.001 * elMolesTotal) { + if (eMolesCalc[m] > 1000. * elMoles[m]) { + resid[m] = -0.5; + resid[m_mm] -= 0.5; + } + if (1000 * eMolesCalc[m] < elMoles[m]) { + resid[m] = 0.5; + resid[m_mm] += 0.5; + } + } + } + if (n_t < (elMolesTotal / nAtomsMax)) { + if (resid[m_mm] < 0.0) { + resid[m_mm] = 0.1; + } + } else if (n_t > elMolesTotal) { + if (resid[m_mm] > 0.0) { + resid[m_mm] = 0.0; + } + } + goto updateSolnVector; + } + + + /* + * Determine whether the matrix should be dumbed down because + * the coefficient matrix of species (with significant concentrations) + * is rank deficient. + * + * The basic idea is that at any time during the calculation only a + * small subset of species with sufficient concentration matters. + * If the rank of the element coefficient matrix for that subset of species + * is less than the number of elements, then the matrix created by + * the Brinkley method below may become singular. + * + * The logic below looks for obvious cases where the current element + * coefficient matrix is rank deficient. + * + * The way around rank-deficiency is to lump-sum the corresponding row of the + * matrix. Note, lump-summing seems to work very well in terms of + * its stability properties, i.e., it heads in the right direction, + * albeit with lousy convergence rates. + * + * NOTE: This probably should be extended to a full blown Gauss-Jordon + * factorization scheme in the future. For Example + * the scheme below would fail for the set: HCl NH4Cl, NH3. + * Hopefully, it's caught by the equal rows logic below. + */ + for (m = 0; m < m_mm; m++) { + lumpSum[m] = 1; + } + double nCutoff = 1.0E-4 * n_t_calc; +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf(" Lump Sum Elements Calculation: \n"); + } +#endif + for (m = 0; m < m_mm; m++) { + int kMSp = -1; + int kMSp2 = -1; + int nSpeciesWithElem = 0; + for (k = 0; k < m_kk; k++) { + if (n_i_calc[k] > nCutoff) { + if (nAtoms(k,m) > 0) { + nSpeciesWithElem++; + if (kMSp != -1) { + kMSp2 = k; + double factor = nAtoms(kMSp,m) / nAtoms(kMSp2,m); + for (n = 0; n < m_mm; n++) { + if (fabs(factor * nAtoms(kMSp2,n) - nAtoms(kMSp,n)) > 1.0E-8) { + lumpSum[m] = 0; + break; + } + } + } else { + kMSp = k; + } + } + } + } +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + string nnn = eNames[m]; + printf(" %5s %3d : %5d %5d\n",nnn.c_str(), lumpSum[m], kMSp, kMSp2); + + } +#endif + } + + + + for (m = 0; m < m_mm; m++) { + for (n = 0; n < m_mm; n++) { + a1(m,n) = 0.0; + for (k = 0; k < m_kk; k++) { + a1(m,n) += nAtoms(k,m) * nAtoms(k,n) * n_i_calc[k]; + } + } + a1(m,m_mm) = eMolesCalc[m]; + a1(m_mm, m) = eMolesCalc[m]; + a1(m_mm, m_mm) = 0.0; + } + + sum = 0.0; + for (m = 0; m < m_mm; m++) { + resid[m] = elMoles[m] - eMolesCalc[m]; + tmp = resid[m] / (elMoles[m] + options.absElemTol); + sum += tmp * tmp; + } + + 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); + } +#endif + for (n = 0; n < m_mm; n++) { + a1(m,n) = 0.0; + } + a1(m,m) = 1.0; + if (resid[m] > 0.0) { + resid[m] = 1.0; + } else if (resid[m] < 0.0) { + resid[m] = -1.0; + } else { + resid[m] = 0.0; + } + } + } + + + resid[m_mm] = n_t - n_t_calc; + +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + printf(" ["); + for (n = 0; n <= m_mm; n++) { + printf(" %10.5g", a1(m,n)); + } + printf("] = %10.5g\n", resid[m]); + } + } +#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); + } +#endif + if (sum < 100. * options.relTolerance) { + retn = 0; + goto exit; + } + + + /* + * Row Sum scaling + */ + for (m = 0; m <= m_mm; m++) { + tmp = 0.0; + for (n = 0; n <= m_mm; n++) { + 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); + } +#endif + for (n = 0; n <= m_mm; n++) { + if (n != m) { + a1(m,n) = 0.0; + a1(n,m) = 0.0; + } + } + } + tmp = 1.0/tmp; + for (n = 0; n <= m_mm; n++) { + a1(m,n) *= tmp; + } + resid[m] *= tmp; + } + +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("Row Summed Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + printf(" ["); + for (n = 0; n <= m_mm; n++) { + printf(" %10.5g", a1(m,n)); + } + printf("] = %10.5g\n", resid[m]); + } + } +#endif + + /* + * Next Step: We have row-summed the equations. + * However, there are some degenerate cases where two + * rows will be multiplies of each other in terms of + * 0 < m, 0 < m part of the matrix. This occurs on a case + * by case basis, and depends upon the current state of the + * element potential values, which affect the concentrations + * of species. + * So, the way we have found to eliminate this problem is to + * lump-sum one of the rows of the matrix, except for the + * last column, and stick it all on the diagonal. + * Then, we at least have a non-singular matrix, and the + * modified equation moves the corresponding unknown in the + * correct direction. + * The previous row-sum operation has made the identification + * of identical rows much simpler. + * + * Note at least 6E-4 is necessary for the comparison. + * I'm guessing 1.0E-3. If two rows are anywhere close to being + * equivalent, the algorithm can get stuck in an oscillatory mode. + */ + modifiedMatrix = false; + for (m = 0; m < m_mm; m++) { + int sameAsRow = -1; + for (int im = 0; im < m; im++) { + bool theSame = true; + for (n = 0; n < m_mm; n++) { + if (fabs(a1(m,n) - a1(im,n)) > 1.0E-3) { + theSame = false; + break; + } + } + if (theSame) { + sameAsRow = im; + } + } + if (sameAsRow >= 0 || lumpSum[m]) { +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + if (lumpSum[m]) { + printf("Lump summing row %d, due to rank deficiency analysis\n", m); + } else if (sameAsRow >= 0) { + printf("Identified that rows %d and %d are the same\n", m, sameAsRow); + } + } +#endif + modifiedMatrix = true; + for (n = 0; n < m_mm; n++) { + if (n != m) { + a1(m,m) += fabs(a1(m,n)); + a1(m,n) = 0.0; + } + } + } + } + +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0 && modifiedMatrix) { + printf("Row Summed, MODIFIED Matrix:\n"); + for (m = 0; m <= m_mm; m++) { + printf(" ["); + for (n = 0; n <= m_mm; n++) { + printf(" %10.5g", a1(m,n)); + } + printf("] = %10.5g\n", resid[m]); + } + } +#endif + + try { + info = solve(a1, DATA_PTR(resid)); + } + catch (CanteraError) { + addLogEntry("estimateEP_Brinkley:Jacobian is singular."); +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + printf("Matrix is SINGULAR.ERROR\n"); + } +#endif + s.restoreState(state); + throw CanteraError("equilibrate:estimateEP_Brinkley()", + "Jacobian is singular. \nTry adding more species, " + "changing the elemental composition slightly, \nor removing " + "unused elements."); + return -3; + } + + /* + * Figure out the damping coefficient: Use a delta damping + * coefficient formulation: magnitude of change is capped + * to exp(1). + */ + beta = 1.0; + for (m = 0; m < m_mm; m++) { + if (resid[m] > 1.0) { + double betat = 1.0 / resid[m]; + if (betat < beta) { + beta = betat; + } + } + if (resid[m] < -1.0) { + double betat = -1.0 / resid[m]; + if (betat < beta) { + beta = betat; + } + } + } +#ifdef DEBUG_HKM_EPEQUIL + if (debug_prnt_lvl > 0) { + if (beta != 1.0) { + printf("(it %d) Beta = %g\n", iter, beta); + } + } +#endif + + /* + * Update the solution vector + */ + updateSolnVector: + for (m = 0; m < m_mm; m++) { + x[m] += beta * resid[m]; + } + 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); + 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]); + } + printf(" n_t %10.5g %10.5g %10.5g \n", x_old[m_mm], n_t, exp(resid[m_mm])); + } +#endif + } + exit: + return retn; + } } // namespace - - -// $Log: ChemEquil.cpp,v diff --git a/Cantera/src/ChemEquil.h b/Cantera/src/ChemEquil.h index b7c47e9a1..06a1a90f4 100755 --- a/Cantera/src/ChemEquil.h +++ b/Cantera/src/ChemEquil.h @@ -27,149 +27,182 @@ namespace Cantera { - int _equilflag(const char* xy); + int _equilflag(const char* xy); - /** - * Chemical equilibrium options. Used internally by class ChemEquil. - */ - class EquilOpt { - public: - EquilOpt() : relTolerance(1.e-10), maxIterations(1000), iterations(0), - maxStepSize(10.0), propertyPair(TP), contin(false) {} + /** + * Chemical equilibrium options. Used internally by class ChemEquil. + */ + class EquilOpt { + public: + EquilOpt() : relTolerance(1.e-8), absElemTol(1.0E-70),maxIterations(1000), + iterations(0), + maxStepSize(10.0), propertyPair(TP), contin(false) {} - doublereal relTolerance; ///< Relative tolerance - int maxIterations; ///< Maximum number of iterations - int iterations; ///< Iteration counter - - /** - * Maximum step size. Largest change in any element potential or - * in log(T) allowed in one Newton step. Default: 10.0 - */ - doublereal maxStepSize; - - /** - * Property pair flag. Determines which two thermodynamic properties - * are fixed. - */ - int propertyPair; - - /** - * 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. - */ - bool contin; - }; - - template - class PropertyCalculator; + doublereal relTolerance; ///< Relative tolerance + doublereal absElemTol; ///< Abs Tol in element number + int maxIterations; ///< Maximum number of iterations + int iterations; ///< Iteration counter /** - * @defgroup equil Chemical Equilibrium - * + * Maximum step size. Largest change in any element potential or + * in log(T) allowed in one Newton step. Default: 10.0 */ + doublereal maxStepSize; + + /** + * Property pair flag. Determines which two thermodynamic properties + * are fixed. + */ + int propertyPair; + + /** + * 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. + */ + bool contin; + }; + + template + class PropertyCalculator; + + /** + * @defgroup equil 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. + * @ingroup equil + */ + class ChemEquil { + + public: + ChemEquil(); + virtual ~ChemEquil(); + + int equilibrate(thermo_t& s, const char* XY, + bool useThermoPhaseElementPotentials = false); + int equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles, + bool useThermoPhaseElementPotentials = false); + const vector_fp& elementPotentials() const { return m_lambda; } /** - * 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. - * @ingroup equil - */ - class ChemEquil { - - public: - ChemEquil(); - virtual ~ChemEquil(); - - 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; } - - /** - * Options controlling how the calculation is carried out. - * @see EquilOptions - */ - EquilOpt options; + * Options controlling how the calculation is carried out. + * @see EquilOptions + */ + EquilOpt options; - protected: + protected: - thermo_t* m_phase; - thermo_t* m_thermo; + thermo_t* m_phase; + thermo_t* m_thermo; - /// number of atoms of element m in species k. - doublereal nAtoms(int k, int m) const { return m_comp[k*m_mm + m]; } + /// number of atoms of element m in species k. + doublereal nAtoms(int k, int m) const { return m_comp[k*m_mm + m]; } - void initialize(thermo_t& s); + void initialize(thermo_t& s); - void setToEquilState(thermo_t& s, - const vector_fp& x, doublereal t); + void setToEquilState(thermo_t& s, + const vector_fp& x, doublereal t); - int setInitialMoles(thermo_t& s); + int setInitialMoles(thermo_t& s); - int estimateElementPotentials(thermo_t& s, - vector_fp& lambda); + 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, double xval, double yval ); + int estimateEP_Brinkley(thermo_t&s, vector_fp& lambda, vector_fp& elMoles); - void equilResidual(thermo_t& s, const vector_fp& x, - const vector_fp& elmtotal, vector_fp& resid, - double xval, double yval); + int dampStep(thermo_t& s, vector_fp& oldx, + double oldf, vector_fp& grad, vector_fp& step, vector_fp& x, + double& f, vector_fp& elmols, double xval, double yval ); - void equilJacobian(thermo_t& s, vector_fp& x, - const vector_fp& elmols, DenseMatrix& jac, - double xval, double yval); + void equilResidual(thermo_t& s, const vector_fp& x, + const vector_fp& elmtotal, vector_fp& resid, + double xval, double yval); - void update(const thermo_t& s); + void equilJacobian(thermo_t& s, vector_fp& x, + const vector_fp& elmols, DenseMatrix& jac, + double xval, double yval); - int m_mm; - int m_kk; - int m_skip; + void update(const thermo_t& s); - PropertyCalculator *m_p1, *m_p2; + int m_mm; + int m_kk; + int m_skip; - vector_fp m_molefractions; - vector_fp m_lambda; - vector_fp m_elementmolefracs; - vector_fp m_reswork; - vector_fp m_jwork1; - vector_fp m_jwork2; - vector_fp m_comp; - doublereal m_temp, m_dens; - doublereal m_p0; - int m_eloc; - doublereal m_abscharge; + PropertyCalculator *m_p1, *m_p2; - doublereal m_startTemp, m_startDens; - vector_fp m_startSoln; + /** + * Current value of the mole fractions in the single phase. + * -> length = m_kk. + */ + vector_fp m_molefractions; + /** + * Current value of the dimensional element potentials + * -> length = m_mm + */ + vector_fp m_lambda; - vector_fp m_grt; - vector_fp m_mu_RT; - vector_int m_component; - }; + /* + * Current value of the sum of the element abundances given the + * current element potentials. + */ + doublereal m_elementTotalSum; + /* + * Current value of the element mole fractions. Note these aren't + * the goal element mole fractions. + */ + vector_fp m_elementmolefracs; + vector_fp m_reswork; + vector_fp m_jwork1; + vector_fp m_jwork2; + /* + * Storage of the element compositions + * natom(k,m) = m_comp[k*m_mm+ m]; + */ + vector_fp m_comp; + doublereal m_temp, m_dens; + doublereal m_p0; + int m_eloc; + doublereal m_abscharge; + + doublereal m_startTemp, m_startDens; + vector_fp m_startSoln; + + vector_fp m_grt; + vector_fp m_mu_RT; + vector_int m_component; + + /* + * element fractional cutoff, below which the element will be + * zeroed. + */ + double m_elemFracCutoff; + bool m_doResPerturb; + }; } #endif -