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

[Equil] Implement equilibrate as a method of ThermoPhase

Browse source 
The Cython interface now calls this equilibrate function.

Unlike the free function 'equilibrate', this new equilibrate method actually
does what the documentation describes regarding the 'solver' argument. This
means that there are cases where failures of the MultiphaseEquil solver were
previously being hidden by a subsequent successful solution using the ChemEquil
solver. These tests have been marked as known failures.
This commit is contained in:
Ray Speth 2014-11-08 00:53:39 +00:00
parent 4bcbe862a9
commit e16e998dd7
5 changed files with 164 additions and 16 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

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

@ -1137,6 +1137,41 @@ public:
* @{
*/
//! Equilibrate a ThermoPhase object
/*!
* Set this phase to chemical equilibrium by calling one of several
* equilibrium solvers. The XY parameter indicates what two thermodynamic
* quantities are to be held constant during the equilibration process.
*
* @param XY String representation of what two properties are being
* held constant
* @param solver Name of the solver to be used to equilibrate the phase.
* If solver = 'element_potential', the ChemEquil element potential
* solver will be used. If solver = 'vcs', the VCS solver will be used.
* If solver = 'gibbs', the MultiPhaseEquil solver will be used. If
* solver = 'auto', the solvers will be tried in order if the initial
* solver(s) fail.
* @param rtol Relative tolerance
* @param max_steps Maximum number of steps to take to find the solution
* @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 log_level loglevel Controls amount of diagnostic output.
* log_level=0 suppresses diagnostics, and increasingly-verbose
* messages are written as loglevel increases.
*
* @ingroup equilfunctions
*/
void equilibrate(const std::string& XY, const std::string& solver="auto",
double rtol=1e-9, int max_steps=50000, int max_iter=100,
int estimate_equil=0, int log_level=0);
//!This method is used by the ChemEquil equilibrium solver.
/*!
* It sets the state such that the chemical potentials satisfy

1
interfaces/cython/cantera/_cantera.pxd
View file

@ -54,6 +54,7 @@ cdef extern from "cantera/thermo/ThermoPhase.h" namespace "Cantera":
double maxTemp() except +
double refPressure() except +
cbool getElementPotentials(double*) except +
void equilibrate(string, string, double, int, int, int, int) except +
# basic thermodynamic properties
double temperature() except +

14
interfaces/cython/cantera/test/test_equilibrium.py
View file

@ -85,6 +85,20 @@ class MultiphaseEquilTest(EquilTestCases, utilities.CanteraTest):
EquilTestCases.__init__(self, 'gibbs')
unittest.TestCase.__init__(self, *args, **kwargs)
@unittest.expectedFailure
def test_equil_gri_stoichiometric(self):
gas = ct.Solution('gri30.xml')
gas.TPX = 301, 100000, 'CH4:1.0, O2:2.0'
gas.equilibrate('TP', self.solver)
self.check(gas, CH4=0, O2=0, H2O=2, CO2=1)
@unittest.expectedFailure
def test_equil_gri_lean(self):
gas = ct.Solution('gri30.xml')
gas.TPX = 301, 100000, 'CH4:1.0, O2:3.0'
gas.equilibrate('TP', self.solver)
self.check(gas, CH4=0, O2=1, H2O=2, CO2=1)
class VCS_EquilTest(EquilTestCases, utilities.CanteraTest):
def __init__(self, *args, **kwargs):

42
interfaces/cython/cantera/thermo.pyx
View file

@ -1,3 +1,5 @@
import warnings
cdef enum Thermasis:
mass_basis = 0
molar_basis = 1
@ -99,7 +101,8 @@ cdef class ThermoPhase(_SolutionBase):
return 1.0
def equilibrate(self, XY, solver='auto', double rtol=1e-9,
int maxsteps=1000, int maxiter=100, int loglevel=0):
int maxsteps=1000, int maxiter=100, int estimate_equil=0,
int loglevel=0):
"""
Set to a state of chemical equilibrium holding property pair
*XY* constant.
@ -127,26 +130,33 @@ cdef class ThermoPhase(_SolutionBase):
For the Gibbs minimization solver, this specifies the number of
'outer' iterations on T or P when some property pair other
than TP is specified.
: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 loglevel:
Set to a value > 0 to write diagnostic output.
"""
cdef int iSolver
if isinstance(solver, int):
iSolver = solver
elif solver == 'auto':
iSolver = -1
elif solver == 'element_potential':
iSolver = 0
elif solver == 'gibbs':
iSolver = 1
elif solver == 'vcs':
iSolver = 2
else:
raise ValueError('Invalid equilibrium solver specified')
warnings.warn('ThermoPhase.equilibrate: Using integer solver '
'flags is deprecated, and will be disabled after Cantera 2.2.')
if solver == -1:
solver = 'auto'
elif solver == 0:
solver = 'element_potential'
elif solver == 1:
solver = 'gibbs'
elif solver == 2:
solver = 'vcs'
else:
raise ValueError('Invalid equilibrium solver specified: '
'"{0}"'.format(solver))
XY = XY.upper()
equilibrate(deref(self.thermo), stringify(XY).c_str(),
iSolver, rtol, maxsteps, maxiter, loglevel)
self.thermo.equilibrate(stringify(XY.upper()), stringify(solver), rtol,
maxsteps, maxiter, estimate_equil, loglevel)
####### Composition, species, and elements ########

88
src/thermo/ThermoPhase.cpp
View file

@ -12,6 +12,9 @@
#include "cantera/thermo/ThermoFactory.h"
#include "cantera/thermo/SpeciesThermoInterpType.h"
#include "cantera/thermo/GeneralSpeciesThermo.h"
#include "cantera/equil/ChemEquil.h"
#include "cantera/equil/MultiPhase.h"
#include "cantera/equil/vcs_MultiPhaseEquil.h"
#include "cantera/base/ctml.h"
#include "cantera/base/vec_functions.h"
@ -749,6 +752,91 @@ void ThermoPhase::setStateFromXML(const XML_Node& state)
}
}
void ThermoPhase::equilibrate(const std::string& XY, const std::string& solver,
double rtol, int max_steps, int max_iter,
int estimate_equil, int log_level)
{
vector_fp initial_state;
saveState(initial_state);
if (solver == "auto" || solver == "element_potential") {
writelog("Trying ChemEquil solver\n", log_level);
try {
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);
if (ret < 0) {
throw CanteraError("ThermoPhase::equilibrate",
"ChemEquil solver failed. Return code: " + int2str(ret));
}
setElementPotentials(E.elementPotentials());
writelog("ChemEquil solver succeeded\n", log_level);
return;
} catch (std::exception& err) {
writelog("ChemEquil solver failed.\n", log_level);
writelog(err.what(), log_level);
restoreState(initial_state);
if (solver == "auto") {
} else {
throw;
}
}
}
int ixy = _equilflag(XY.c_str());
if (solver == "auto" || solver == "vcs") {
try {
writelog("Trying VCS equilibrium solver\n", log_level);
MultiPhase M;
M.addPhase(this, 1.0);
M.init();
VCSnonideal::vcs_MultiPhaseEquil eqsolve(&M, log_level-1);
int ret = eqsolve.equilibrate(ixy, estimate_equil, log_level-1,
rtol, max_steps);
if (ret) {
throw CanteraError("ThermoPhase::equilibrate",
"VCS solver failed. Return code: " + int2str(ret));
}
writelog("VCS solver succeeded\n");
return;
} catch (std::exception& err) {
writelog("VCS solver failed.\n", log_level);
writelog(err.what(), log_level);
restoreState(initial_state);
if (solver == "auto") {
} else {
throw;
}
}
}
if (solver == "auto" || solver == "gibbs") {
try {
writelog("Trying MultiPhaseEquil (Gibbs) equilibrium solver\n",
log_level);
MultiPhase M;
M.addPhase(this, 1.0);
M.init();
M.equilibrate(ixy, rtol, max_steps, max_iter, log_level-1);
writelog("MultiPhaseEquil solver succeeded\n");
return;
} catch (std::exception& err) {
writelog("MultiPhaseEquil solver failed.\n", log_level);
writelog(err.what(), log_level);
restoreState(initial_state);
throw;
}
}
if (solver != "auto") {
throw CanteraError("ThermoPhase::equilibrate",
"Invalid solver specified: '" + solver + "'");
}
}
void ThermoPhase::setElementPotentials(const vector_fp& lambda)
{
size_t mm = nElements();