From 8c934ab6785b6cffc73f2676019a0bc6b070c725 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Sat, 8 Nov 2014 00:53:42 +0000 Subject: [PATCH] [Equil] Implement equilibrate as a method of MultiPhase This differs from the current equilibrate method in that it can use either of the multiphase equilibrium solvers. --- include/cantera/equil/MultiPhase.h | 33 +++++++++++++ interfaces/cython/cantera/_cantera.pxd | 8 +-- interfaces/cython/cantera/mixture.pyx | 44 +++++++++-------- src/equil/MultiPhase.cpp | 67 ++++++++++++++++++++++++++ src/thermo/ThermoPhase.cpp | 56 ++++----------------- 5 files changed, 135 insertions(+), 73 deletions(-) diff --git a/include/cantera/equil/MultiPhase.h b/include/cantera/equil/MultiPhase.h index da63a8424..61ca73858 100644 --- a/include/cantera/equil/MultiPhase.h +++ b/include/cantera/equil/MultiPhase.h @@ -350,6 +350,39 @@ public: doublereal equilibrate(int XY, doublereal err = 1.0e-9, int maxsteps = 1000, int maxiter = 200, int loglevel = -99); + //! Equilibrate a MultiPhase object + /*! + * Set this mixture to chemical equilibrium by calling one of Cantera's + * 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 = 'vcs', the vcs_MultiPhaseEquil solver will be used. If + * solver = 'gibbs', the MultiPhaseEquil solver will be used. If solver + * = 'auto', the 'vcs' solver will be tried first, followed by the + * 'gibbs' solver if the first one fails. + * @param rtol Relative tolerance + * @param max_steps Maximum number of steps to take to find the solution + * @param max_iter 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); + /// Set the temperature [K]. /*! * @param T value of the temperature (Kelvin) diff --git a/interfaces/cython/cantera/_cantera.pxd b/interfaces/cython/cantera/_cantera.pxd index a7a13eeb4..3a34ab125 100644 --- a/interfaces/cython/cantera/_cantera.pxd +++ b/interfaces/cython/cantera/_cantera.pxd @@ -211,6 +211,8 @@ cdef extern from "cantera/equil/MultiPhase.h" namespace "Cantera": void addPhase(CxxThermoPhase*, double) except + void init() except + + void equilibrate(string, string, double, int, int, int, int) except + + size_t nSpecies() size_t nElements() size_t nPhases() @@ -243,12 +245,6 @@ cdef extern from "cantera/equil/MultiPhase.h" namespace "Cantera": double cp() except + double volume() except + -cdef extern from "cantera/equil/equil.h" namespace "Cantera": - int equilibrate(CxxThermoPhase&, char*, int, double, int, int, int) except + - -cdef extern from "cantera/equil/vcs_MultiPhaseEquil.h" namespace "Cantera": - int vcs_equilibrate(CxxMultiPhase&, char*, int, int, int, double, int, int, int) except + - cdef extern from "cantera/zeroD/ReactorBase.h" namespace "Cantera": cdef cppclass CxxWall "Cantera::Wall" diff --git a/interfaces/cython/cantera/mixture.pyx b/interfaces/cython/cantera/mixture.pyx index 86392e5ba..36637b7cf 100644 --- a/interfaces/cython/cantera/mixture.pyx +++ b/interfaces/cython/cantera/mixture.pyx @@ -1,3 +1,5 @@ +import warnings + cdef class Mixture: """ @@ -271,8 +273,8 @@ cdef class Mixture: self.mix.getChemPotentials(&data[0]) return data - def equilibrate(self, XY, solver='vcs', rtol=1e-9, max_steps=1000, - max_iter=100, estimate_equil=0, print_level=0, log_level=0): + def equilibrate(self, XY, solver='auto', rtol=1e-9, max_steps=1000, + max_iter=100, estimate_equil=0, log_level=0): """ Set to a state of chemical equilibrium holding property pair *XY* constant. This method uses a version of the VCS algorithm to find the @@ -284,11 +286,11 @@ cdef class Mixture: A two-letter string, which must be one of the set:: ['TP', 'HP', 'SP'] - :param solver: - Set to either 'vcs' or 'gibbs' to choose implementation - of the solver to use. 'vcs' uses the solver implemented in the - C++ class 'VCSnonideal', and 'gibbs' uses the one implemented - in class 'MultiPhaseEquil'. + :param solver: Set to either 'auto', 'vcs', or 'gibbs' to choose + implementation of the solver to use. 'vcs' uses the solver + implemented in the C++ class 'VCSnonideal', 'gibbs' uses the one + implemented in class 'MultiPhaseEquil'. 'auto' will try the 'vcs' + solver first and then the 'gibbs' solver if that fails. :param rtol: Error tolerance. Iteration will continue until (Delta mu)/RT is less than this value for each reaction. Note that this default is @@ -308,21 +310,23 @@ cdef class Mixture: 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 print_level: + :param log_level: Determines the amount of output displayed during the solution process. 0 indicates no output, while larger numbers produce successively more verbose information. - :param log_level: - Controls the amount of diagnostic output written. """ - if solver == 'vcs': - iSolver = 2 - elif solver == 'gibbs': - iSolver = 1 - else: - raise ValueError('Unrecognized equilibrium solver ' - 'specified: "{0}"'.format(solver)) + if isinstance(solver, int): + warnings.warn('Mixture.equilibrate: Using integer solver flags is ' + 'deprecated, and will be disabled after Cantera 2.2.') + if solver == -1: + solver = 'auto' + elif solver == 1: + solver = 'gibbs' + elif solver == 2: + solver = 'vcs' + else: + raise ValueError('Unrecognized equilibrium solver ' + 'specified: "{0}"'.format(solver)) - vcs_equilibrate(deref(self.mix), stringify(XY).c_str(), estimate_equil, - print_level, iSolver, rtol, max_steps, max_iter, - log_level) + self.mix.equilibrate(stringify(XY.upper()), stringify(solver), rtol, + max_steps, max_iter, estimate_equil, log_level) diff --git a/src/equil/MultiPhase.cpp b/src/equil/MultiPhase.cpp index 9b6cdd0c4..6d7e5d988 100644 --- a/src/equil/MultiPhase.cpp +++ b/src/equil/MultiPhase.cpp @@ -3,8 +3,11 @@ * Definitions for the \link Cantera::MultiPhase MultiPhase\endlink * object that is used to set up multiphase equilibrium problems (see \ref equilfunctions). */ + +#include "cantera/equil/ChemEquil.h" #include "cantera/equil/MultiPhase.h" #include "cantera/equil/MultiPhaseEquil.h" +#include "cantera/equil/vcs_MultiPhaseEquil.h" #include "cantera/base/stringUtils.h" using namespace std; @@ -744,6 +747,70 @@ doublereal MultiPhase::equilibrate(int XY, doublereal err, return -1.0; } +void MultiPhase::equilibrate(const std::string& XY, const std::string& solver, + double rtol, int max_steps, int max_iter, + int estimate_equil, int log_level) +{ + // Save the initial state so that it can be restored in case one of the + // solvers fails + vector_fp initial_moleFractions = m_moleFractions; + vector_fp initial_moles = m_moles; + double initial_T = m_temp; + double initial_P = m_press; + + int ixy = _equilflag(XY.c_str()); + if (solver == "auto" || solver == "vcs") { + try { + writelog("Trying VCS equilibrium solver\n", log_level); + VCSnonideal::vcs_MultiPhaseEquil eqsolve(this, log_level-1); + int ret = eqsolve.equilibrate(ixy, estimate_equil, log_level-1, + rtol, max_steps); + if (ret) { + throw CanteraError("MultiPhase::equilibrate", + "VCS solver failed. Return code: " + int2str(ret)); + } + writelog("VCS solver succeeded\n", log_level); + return; + } catch (std::exception& err) { + writelog("VCS solver failed.\n", log_level); + writelog(err.what(), log_level); + m_moleFractions = initial_moleFractions; + m_moles = initial_moles; + m_temp = initial_T; + m_press = initial_P; + updatePhases(); + if (solver == "auto") { + } else { + throw; + } + } + } + + if (solver == "auto" || solver == "gibbs") { + try { + writelog("Trying MultiPhaseEquil (Gibbs) equilibrium solver\n", + log_level); + equilibrate(ixy, rtol, max_steps, max_iter, log_level-1); + writelog("MultiPhaseEquil solver succeeded\n", log_level); + return; + } catch (std::exception& err) { + writelog("MultiPhaseEquil solver failed.\n", log_level); + writelog(err.what(), log_level); + m_moleFractions = initial_moleFractions; + m_moles = initial_moles; + m_temp = initial_T; + m_press = initial_P; + updatePhases(); + throw; + } + } + + if (solver != "auto") { + throw CanteraError("MultiPhase::equilibrate", + "Invalid solver specified: '" + solver + "'"); + } +} + #ifdef MULTIPHASE_DEVEL void importFromXML(string infile, string id) { diff --git a/src/thermo/ThermoPhase.cpp b/src/thermo/ThermoPhase.cpp index 4ba154ef1..dfe3ff56b 100644 --- a/src/thermo/ThermoPhase.cpp +++ b/src/thermo/ThermoPhase.cpp @@ -14,7 +14,6 @@ #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" @@ -756,10 +755,10 @@ 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") { + vector_fp initial_state; + saveState(initial_state); writelog("Trying ChemEquil solver\n", log_level); try { ChemEquil E; @@ -785,50 +784,13 @@ void ThermoPhase::equilibrate(const std::string& XY, const std::string& solver, } } - 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" || solver == "vcs" || solver == "gibbs") { + MultiPhase M; + M.addPhase(this, 1.0); + M.init(); + M.equilibrate(XY, solver, rtol, max_steps, max_iter, + estimate_equil, log_level); + return; } if (solver != "auto") {