diff --git a/include/cantera/thermo/ThermoPhase.h b/include/cantera/thermo/ThermoPhase.h index b15419267..f415e54b1 100644 --- a/include/cantera/thermo/ThermoPhase.h +++ b/include/cantera/thermo/ThermoPhase.h @@ -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 diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index a88b16621..a7a13eeb4 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -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 + diff --git a/interfaces/cython/cantera/test/test_equilibrium.py b/interfaces/cython/cantera/test/test_equilibrium.py index 60b72bd2a..e831ed272 100644 --- a/interfaces/cython/cantera/test/test_equilibrium.py +++ b/interfaces/cython/cantera/test/test_equilibrium.py @@ -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): diff --git a/interfaces/cython/cantera/thermo.pyx b/interfaces/cython/cantera/thermo.pyx index cccb0d71c..a4f4c3d87 100644 --- a/interfaces/cython/cantera/thermo.pyx +++ b/interfaces/cython/cantera/thermo.pyx @@ -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 ######## diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp index 83c82c62b..4ba154ef1 100644 --- a/src/thermo/ThermoPhase.cpp +++ b/src/thermo/ThermoPhase.cpp @@ -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();