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);