From 74167cc3eb81dc63b306c045e2a7b1e86fb5b08b Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 19 May 2018 22:27:32 -0400 Subject: [PATCH] [Equil] Don't use saved element potentials as initial guess Saved element potentials are only valid at equilibrium, and may not be a good guess for calls to equilibrate() after the state has changed. By always using the estimation method for the element potentials at the start of the ChemEquil algorithm, the results of equilibrate() are repeatable, and do not depend on the results of previous calls to equilibrate(). Resolves #524 --- include/cantera/equil/ChemEquil.h | 8 ++---- include/cantera/thermo/ThermoPhase.h | 12 ++++----- src/equil/ChemEquil.cpp | 39 +++++++--------------------- src/thermo/ThermoPhase.cpp | 3 +-- 4 files changed, 18 insertions(+), 44 deletions(-) diff --git a/include/cantera/equil/ChemEquil.h b/include/cantera/equil/ChemEquil.h index fb211e841..c5e35e966 100644 --- a/include/cantera/equil/ChemEquil.h +++ b/include/cantera/equil/ChemEquil.h @@ -103,8 +103,7 @@ public: * ThermoPhase object. The properties must be already contained within the * current thermodynamic state of the system. */ - int equilibrate(thermo_t& s, const char* XY, - bool useThermoPhaseElementPotentials = false, int loglevel = 0); + int equilibrate(thermo_t& s, const char* XY, int loglevel = 0); /*! * Compute the equilibrium composition for two specified properties and the @@ -117,16 +116,13 @@ public: * @param s phase object to be equilibrated * @param XY property pair to hold constant * @param elMoles specified vector of element abundances. - * @param useThermoPhaseElementPotentials get the initial estimate for the - * chemical potentials from the ThermoPhase object (true) or create - * our own estimate (false) * @param loglevel Specify amount of debug logging (0 to disable) * @return 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 equilibrate(thermo_t& s, const char* XY, vector_fp& elMoles, - bool useThermoPhaseElementPotentials = false, int loglevel = 0); + int loglevel = 0); const vector_fp& elementPotentials() const { return m_lambda; } diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index 861fd36fd..e4810f64f 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -1196,12 +1196,12 @@ public: * @param max_iter For the 'gibbs' and 'vcs' solvers, this is the maximum * number of outer temperature or pressure iterations to take when T * and/or P is not held fixed. - * @param estimate_equil integer indicating whether the solver should - * estimate its own initial condition. If 0, the initial mole fraction - * vector in the ThermoPhase object is used as the initial condition. - * If 1, the initial mole fraction vector is used if the element - * abundances are satisfied. If -1, the initial mole fraction vector is - * thrown out, and an estimate is formulated. + * @param estimate_equil For MultiPhaseEquil solver, an integer indicating + * whether the solver should estimate its own initial condition. If 0, + * the initial mole fraction vector in the ThermoPhase object is used + * as the initial condition. If 1, the initial mole fraction vector is + * used if the element abundances are satisfied. If -1, the initial + * mole fraction vector is thrown out, and an estimate is formulated. * @param log_level loglevel Controls amount of diagnostic output. * log_level=0 suppresses diagnostics, and increasingly-verbose * messages are written as loglevel increases. diff --git a/src/equil/ChemEquil.cpp b/src/equil/ChemEquil.cpp index 9d9fdb687..83cb8de42 100644 --- a/src/equil/ChemEquil.cpp +++ b/src/equil/ChemEquil.cpp @@ -305,20 +305,16 @@ int ChemEquil::estimateElementPotentials(thermo_t& s, vector_fp& lambda_RT, return info; } -int ChemEquil::equilibrate(thermo_t& s, const char* XY, - bool useThermoPhaseElementPotentials, int loglevel) +int ChemEquil::equilibrate(thermo_t& s, const char* XY, int loglevel) { initialize(s); update(s); vector_fp elMolesGoal = m_elementmolefracs; - return equilibrate(s, XY, elMolesGoal, useThermoPhaseElementPotentials, - loglevel-1); + return equilibrate(s, XY, elMolesGoal, loglevel-1); } int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, - vector_fp& elMolesGoal, - bool useThermoPhaseElementPotentials, - int loglevel) + vector_fp& elMolesGoal, int loglevel) { int fail = 0; bool tempFixed = true; @@ -501,29 +497,12 @@ int ChemEquil::equilibrate(thermo_t& s, const char* XYstr, setInitialMoles(s, elMolesGoal,loglevel); - // If requested, get the initial estimate for the chemical potentials from - // the ThermoPhase object itself. Or else, create our own estimate. - if (useThermoPhaseElementPotentials) { - bool haveEm = s.getElementPotentials(x.data()); - if (haveEm) { - if (s.temperature() < 100.) { - writelog("we are here {:g}\n", s.temperature()); - } - for (size_t m = 0; m < m_mm; m++) { - x[m] *= 1.0 / s.RT(); - } - } else { - estimateElementPotentials(s, x, elMolesGoal); - } - } else { - // Calculate initial estimates of the element potentials. This algorithm - // uses 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, elMolesGoal); - } + // Calculate initial estimates of the element potentials. This algorithm + // uses 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, elMolesGoal); // 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 diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp index f5e74e3fd..a5e677408 100644 --- a/src/thermo/ThermoPhase.cpp +++ b/src/thermo/ThermoPhase.cpp @@ -697,8 +697,7 @@ void ThermoPhase::equilibrate(const std::string& XY, const std::string& solver, ChemEquil E; E.options.maxIterations = max_steps; E.options.relTolerance = rtol; - bool use_element_potentials = (estimate_equil == 0); - int ret = E.equilibrate(*this, XY.c_str(), use_element_potentials, log_level-1); + int ret = E.equilibrate(*this, XY.c_str(), log_level-1); if (ret < 0) { throw CanteraError("ThermoPhase::equilibrate", "ChemEquil solver failed. Return code: {}", ret);