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

[Equil] Don't use saved element potentials as initial guess

Browse source 
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
This commit is contained in:
Ray Speth 2018-05-19 22:27:32 -04:00
parent dab739013e
commit 74167cc3eb
4 changed files with 18 additions and 44 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

8
include/cantera/equil/ChemEquil.h
View file

@ -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;
}

12
include/cantera/thermo/ThermoPhase.h
View file

@ -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.

39
src/equil/ChemEquil.cpp
View file

@ -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

3
src/thermo/ThermoPhase.cpp
View file

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