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

[Equil] Implement equilibrate as a method of MultiPhase

Browse source 
This differs from the current equilibrate method in that it can use either of
the multiphase equilibrium solvers.
This commit is contained in:
Ray Speth 2014-11-08 00:53:42 +00:00
parent e16e998dd7
commit 8c934ab678
5 changed files with 135 additions and 73 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

33
include/cantera/equil/MultiPhase.h
View file

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

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

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

44
interfaces/cython/cantera/mixture.pyx
View file

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

67
src/equil/MultiPhase.cpp
View file

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

56
src/thermo/ThermoPhase.cpp
View file

@ -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") {