Cantera  2.5.1
Public Member Functions | Public Attributes | Protected Member Functions | Protected Attributes | List of all members
ChemEquil Class Reference

Class ChemEquil implements a chemical equilibrium solver for single-phase solutions. More...

#include <ChemEquil.h>

Collaboration diagram for ChemEquil:
[legend]

Public Member Functions

 ChemEquil (thermo_t &s)
 Constructor combined with the initialization function. More...
 
int equilibrate (thermo_t &s, const char *XY, int loglevel=0)
 
int equilibrate (thermo_t &s, const char *XY, vector_fp &elMoles, int loglevel=0)
 

Public Attributes

EquilOpt options
 Options controlling how the calculation is carried out. More...
 

Protected Member Functions

doublereal nAtoms (size_t k, size_t m) const
 number of atoms of element m in species k. More...
 
void initialize (thermo_t &s)
 
void setToEquilState (thermo_t &s, const vector_fp &x, doublereal t)
 
int setInitialMoles (thermo_t &s, vector_fp &elMoleGoal, int loglevel=0)
 Estimate the initial mole numbers. More...
 
int estimateElementPotentials (thermo_t &s, vector_fp &lambda, vector_fp &elMolesGoal, int loglevel=0)
 Generate a starting estimate for the element potentials. More...
 
int estimateEP_Brinkley (thermo_t &s, vector_fp &lambda, vector_fp &elMoles)
 
int dampStep (thermo_t &s, vector_fp &oldx, double oldf, vector_fp &grad, vector_fp &step, vector_fp &x, double &f, vector_fp &elmols, double xval, double yval)
 Find an acceptable step size and take it. More...
 
void equilResidual (thermo_t &s, const vector_fp &x, const vector_fp &elmtotal, vector_fp &resid, double xval, double yval, int loglevel=0)
 Evaluates the residual vector F, of length m_mm. More...
 
void equilJacobian (thermo_t &s, vector_fp &x, const vector_fp &elmols, DenseMatrix &jac, double xval, double yval, int loglevel=0)
 
void adjustEloc (thermo_t &s, vector_fp &elMolesGoal)
 
void update (const thermo_t &s)
 Update internally stored state information. More...
 
double calcEmoles (thermo_t &s, vector_fp &x, const double &n_t, const vector_fp &Xmol_i_calc, vector_fp &eMolesCalc, vector_fp &n_i_calc, double pressureConst)
 Given a vector of dimensionless element abundances, this routine calculates the moles of the elements and the moles of the species. More...
 

Protected Attributes

thermo_tm_phase
 Pointer to the ThermoPhase object used to initialize this object. More...
 
size_t m_mm
 number of elements in the phase More...
 
size_t m_kk
 number of species in the phase More...
 
size_t m_skip
 
size_t m_nComponents
 This is equal to the rank of the stoichiometric coefficient matrix when it is computed. More...
 
std::function< double(ThermoPhase &)> m_p1
 
std::function< double(ThermoPhase &)> m_p2
 
vector_fp m_molefractions
 Current value of the mole fractions in the single phase. length = m_kk. More...
 
doublereal m_elementTotalSum
 Current value of the sum of the element abundances given the current element potentials. More...
 
vector_fp m_elementmolefracs
 Current value of the element mole fractions. More...
 
vector_fp m_reswork
 
vector_fp m_jwork1
 
vector_fp m_jwork2
 
vector_fp m_comp
 Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];. More...
 
doublereal m_temp
 
doublereal m_dens
 
doublereal m_p0
 
size_t m_eloc
 Index of the element id corresponding to the electric charge of each species. More...
 
vector_fp m_startSoln
 
vector_fp m_grt
 
vector_fp m_mu_RT
 
vector_fp m_muSS_RT
 Dimensionless values of the Gibbs free energy for the standard state of each species, at the temperature and pressure of the solution (the star standard state). More...
 
std::vector< size_t > m_component
 
double m_elemFracCutoff
 element fractional cutoff, below which the element will be zeroed. More...
 
bool m_doResPerturb
 
std::vector< size_t > m_orderVectorElements
 
std::vector< size_t > m_orderVectorSpecies
 
int m_loglevel
 Verbosity of printed output. More...
 

Detailed Description

Class ChemEquil implements a chemical equilibrium solver for single-phase solutions.

It is a "non-stoichiometric" solver in the terminology of Smith and Missen, meaning that every intermediate state is a valid chemical equilibrium state, but does not necessarily satisfy the element constraints. In contrast, the solver implemented in class MultiPhaseEquil uses a "stoichiometric" algorithm, in which each intermediate state satisfies the element constraints but is not a state of chemical equilibrium. Non- stoichiometric methods are faster when they converge, but stoichiometric ones tend to be more robust and can be used also for problems with multiple condensed phases. As expected, the ChemEquil solver is faster than MultiPhaseEquil for many single-phase equilibrium problems (particularly if there are only a few elements but very many species), but can be less stable. Problem situations include low temperatures where only a few species have non-zero mole fractions, precisely stoichiometric compositions (e.g. 2 H2 + O2). In general, if speed is important, this solver should be tried first, and if it fails then use MultiPhaseEquil.

Definition at line 82 of file ChemEquil.h.

Constructor & Destructor Documentation

◆ ChemEquil()

ChemEquil ( thermo_t s)

Constructor combined with the initialization function.

This constructor initializes the ChemEquil object with everything it needs to start solving equilibrium problems.

Parameters
sThermoPhase object that will be used in the equilibrium calls.

Definition at line 50 of file ChemEquil.cpp.

References ChemEquil::initialize().

Member Function Documentation

◆ equilibrate() [1/2]

int equilibrate ( thermo_t s,
const char *  XY,
int  loglevel = 0 
)

Equilibrate a phase, holding the elemental composition fixed at the initial value found within the ThermoPhase object s.

The value of two specified properties are obtained by querying the ThermoPhase object. The properties must be already contained within the current thermodynamic state of the system.

Definition at line 308 of file ChemEquil.cpp.

References ChemEquil::initialize(), ChemEquil::m_elementmolefracs, and ChemEquil::update().

◆ equilibrate() [2/2]

int equilibrate ( thermo_t s,
const char *  XY,
vector_fp elMoles,
int  loglevel = 0 
)

Compute the equilibrium composition for two specified properties and the specified element moles.

The two specified properties are obtained by querying the ThermoPhase object. The properties must be already contained within the current thermodynamic state of the system.

Parameters
sphase object to be equilibrated
XYproperty pair to hold constant
elMolesspecified vector of element abundances.
loglevelSpecify amount of debug logging (0 to disable)
Returns
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.

Definition at line 316 of file ChemEquil.cpp.

References Cantera::_equilflag(), EquilOpt::absElemTol, Cantera::clip(), ChemEquil::dampStep(), Array2D::data(), Phase::density(), Cantera::dot(), Phase::elementName(), ThermoPhase::enthalpy_mass(), ThermoPhase::entropy_mass(), ChemEquil::equilResidual(), ChemEquil::estimateElementPotentials(), ChemEquil::estimateEP_Brinkley(), CanteraError::getMessage(), ChemEquil::initialize(), ThermoPhase::intEnergy_mass(), EquilOpt::iterations, ChemEquil::m_elemFracCutoff, ChemEquil::m_eloc, ChemEquil::m_kk, ChemEquil::m_loglevel, ChemEquil::m_mm, ChemEquil::m_nComponents, EquilOpt::maxIterations, ThermoPhase::maxTemp(), ThermoPhase::minTemp(), Phase::moleFraction(), Phase::nElements(), Cantera::npos, Phase::nSpecies(), ChemEquil::options, Phase::pressure(), EquilOpt::relTolerance, Phase::restoreState(), Phase::saveState(), Cantera::scale(), ChemEquil::setInitialMoles(), Phase::setMoleFractions(), Phase::setTemperature(), ChemEquil::setToEquilState(), Cantera::solve(), Phase::temperature(), ChemEquil::update(), Cantera::warn_user(), Cantera::writelog(), and Cantera::writelogf().

◆ nAtoms()

doublereal nAtoms ( size_t  k,
size_t  m 
) const
inlineprotected

number of atoms of element m in species k.

Definition at line 143 of file ChemEquil.h.

References ChemEquil::m_comp, and ChemEquil::m_mm.

Referenced by ChemEquil::calcEmoles(), ChemEquil::estimateElementPotentials(), ChemEquil::setToEquilState(), and ChemEquil::update().

◆ initialize()

void initialize ( thermo_t s)
protected

◆ setToEquilState()

void setToEquilState ( thermo_t s,
const vector_fp x,
doublereal  t 
)
protected

Set mixture to an equilibrium state consistent with specified element potentials and temperature.

Parameters
smixture to be updated
xvector of non-dimensional element potentials

\[ \lambda_m/RT \]

.
ttemperature in K.

Definition at line 132 of file ChemEquil.cpp.

References ChemEquil::m_kk, ChemEquil::m_mm, ChemEquil::nAtoms(), Phase::setTemperature(), ThermoPhase::setToEquilState(), and ChemEquil::update().

Referenced by ChemEquil::equilibrate(), and ChemEquil::equilResidual().

◆ setInitialMoles()

int setInitialMoles ( thermo_t s,
vector_fp elMoleGoal,
int  loglevel = 0 
)
protected

◆ estimateElementPotentials()

int estimateElementPotentials ( thermo_t s,
vector_fp lambda,
vector_fp elMolesGoal,
int  loglevel = 0 
)
protected

◆ estimateEP_Brinkley()

int estimateEP_Brinkley ( thermo_t s,
vector_fp lambda,
vector_fp elMoles 
)
protected

Do a calculation of the element potentials using the Brinkley method, p. 129 Smith and Missen.

We have found that the previous estimate may not be good enough to avoid drastic numerical issues associated with the use of a numerically generated Jacobian used in the main algorithm.

The Brinkley algorithm, here, assumes a constant T, P system and uses a linearized analytical Jacobian that turns out to be very stable even given bad initial guesses.

The pressure and temperature to be used are in the ThermoPhase object input into the routine.

The initial guess for the element potentials used by this routine is taken from the input vector, x.

elMoles is the input element abundance vector to be matched.

Nonideal phases are handled in principle. This is done by calculating the activity coefficients and adding them into the formula in the correct position. However, these are treated as a RHS contribution only. Therefore, convergence might be a problem. This has not been tested. Also molality based unit systems aren't handled.

On return, int return value contains the success code:

  • 0 - successful
  • 1 - unsuccessful, max num iterations exceeded
  • -3 - unsuccessful, singular Jacobian

NOTE: update for activity coefficients.

Definition at line 843 of file ChemEquil.cpp.

References ChemEquil::m_kk, ChemEquil::m_mm, and Phase::saveState().

Referenced by ChemEquil::equilibrate().

◆ dampStep()

int dampStep ( thermo_t s,
vector_fp oldx,
double  oldf,
vector_fp grad,
vector_fp step,
vector_fp x,
double &  f,
vector_fp elmols,
double  xval,
double  yval 
)
protected

Find an acceptable step size and take it.

The original implementation employed a line search technique that enforced a reduction in the norm of the residual at every successful step. Unfortunately, this method created false convergence errors near the end of a significant number of steps, usually special conditions where there were stoichiometric constraints.

This new method just does a delta damping approach, based on limiting the jump in the dimensionless element potentials. Mole fractions are limited to a factor of 2 jump in the values from this method. Near convergence, the delta damping gets out of the way.

Definition at line 689 of file ChemEquil.cpp.

References ChemEquil::m_eloc, ChemEquil::m_loglevel, ChemEquil::m_mm, Cantera::writelog(), and Cantera::writelogf().

Referenced by ChemEquil::equilibrate().

◆ equilResidual()

void equilResidual ( thermo_t s,
const vector_fp x,
const vector_fp elmtotal,
vector_fp resid,
double  xval,
double  yval,
int  loglevel = 0 
)
protected

◆ update()

void update ( const thermo_t s)
protected

◆ calcEmoles()

double calcEmoles ( thermo_t s,
vector_fp x,
const double &  n_t,
const vector_fp Xmol_i_calc,
vector_fp eMolesCalc,
vector_fp n_i_calc,
double  pressureConst 
)
protected

Given a vector of dimensionless element abundances, this routine calculates the moles of the elements and the moles of the species.

Parameters
[in]x= current dimensionless element potentials..

Definition at line 807 of file ChemEquil.cpp.

References ThermoPhase::getActivityCoefficients(), ChemEquil::m_kk, ChemEquil::m_mm, ChemEquil::m_muSS_RT, ChemEquil::nAtoms(), Phase::setMoleFractions(), and Phase::setPressure().

Member Data Documentation

◆ options

EquilOpt options

Options controlling how the calculation is carried out.

See also
EquilOptions

Definition at line 131 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate().

◆ m_phase

thermo_t* m_phase
protected

Pointer to the ThermoPhase object used to initialize this object.

This ThermoPhase object must be compatible with the ThermoPhase objects input from the equilibrate function. Currently, this means that the 2 ThermoPhases have to have consist of the same species and elements.

Definition at line 140 of file ChemEquil.h.

Referenced by ChemEquil::initialize().

◆ m_mm

size_t m_mm
protected

◆ m_kk

size_t m_kk
protected

◆ m_nComponents

size_t m_nComponents
protected

This is equal to the rank of the stoichiometric coefficient matrix when it is computed.

It's initialized to m_mm.

Definition at line 258 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate(), ChemEquil::equilResidual(), ChemEquil::estimateElementPotentials(), ChemEquil::initialize(), and ChemEquil::setInitialMoles().

◆ m_molefractions

vector_fp m_molefractions
protected

Current value of the mole fractions in the single phase. length = m_kk.

Definition at line 263 of file ChemEquil.h.

Referenced by ChemEquil::initialize(), and ChemEquil::update().

◆ m_elementTotalSum

doublereal m_elementTotalSum
protected

Current value of the sum of the element abundances given the current element potentials.

Definition at line 267 of file ChemEquil.h.

Referenced by ChemEquil::update().

◆ m_elementmolefracs

vector_fp m_elementmolefracs
protected

Current value of the element mole fractions.

Note these aren't the goal element mole fractions.

Definition at line 271 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate(), ChemEquil::equilResidual(), ChemEquil::initialize(), ChemEquil::setInitialMoles(), and ChemEquil::update().

◆ m_comp

vector_fp m_comp
protected

Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.

Definition at line 277 of file ChemEquil.h.

Referenced by ChemEquil::initialize(), and ChemEquil::nAtoms().

◆ m_eloc

size_t m_eloc
protected

Index of the element id corresponding to the electric charge of each species.

Equal to -1 if there is no such element id.

Definition at line 283 of file ChemEquil.h.

Referenced by ChemEquil::dampStep(), ChemEquil::equilibrate(), ChemEquil::equilResidual(), and ChemEquil::initialize().

◆ m_muSS_RT

vector_fp m_muSS_RT
protected

Dimensionless values of the Gibbs free energy for the standard state of each species, at the temperature and pressure of the solution (the star standard state).

Definition at line 293 of file ChemEquil.h.

Referenced by ChemEquil::calcEmoles(), and ChemEquil::initialize().

◆ m_elemFracCutoff

double m_elemFracCutoff
protected

element fractional cutoff, below which the element will be zeroed.

Definition at line 297 of file ChemEquil.h.

Referenced by ChemEquil::equilibrate(), and ChemEquil::equilResidual().

◆ m_loglevel

int m_loglevel
protected

Verbosity of printed output.

No messages when m_loglevel == 0. More output as level increases.

Definition at line 305 of file ChemEquil.h.

Referenced by ChemEquil::dampStep(), ChemEquil::equilibrate(), and ChemEquil::setInitialMoles().


The documentation for this class was generated from the following files: