Cantera 2.6.0
ChemEquil.h
Go to the documentation of this file.
1/**
2 * @file ChemEquil.h Chemical equilibrium.
3 */
4
5// This file is part of Cantera. See License.txt in the top-level directory or
6// at https://cantera.org/license.txt for license and copyright information.
7
8#ifndef CT_CHEM_EQUIL_H
9#define CT_CHEM_EQUIL_H
10
12#include <functional>
13
14namespace Cantera
15{
16
17class DenseMatrix;
18class ThermoPhase;
19/// map property strings to integers
20int _equilflag(const char* xy);
21
22/**
23 * Chemical equilibrium options. Used internally by class ChemEquil.
24 */
26{
27public:
28 EquilOpt() : relTolerance(1.e-8), absElemTol(1.0E-70),maxIterations(1000),
29 iterations(0),
30 maxStepSize(10.0), propertyPair(TP), contin(false) {}
31
32 doublereal relTolerance; ///< Relative tolerance
33 doublereal absElemTol; ///< Abs Tol in element number
34 int maxIterations; ///< Maximum number of iterations
35 int iterations; ///< Iteration counter
36
37 /**
38 * Maximum step size. Largest change in any element potential or
39 * in log(T) allowed in one Newton step. Default: 10.0
40 */
41 doublereal maxStepSize;
42
43 /**
44 * Property pair flag. Determines which two thermodynamic properties
45 * are fixed.
46 */
48
49 /**
50 * Continuation flag. Set true if the calculation should be initialized from
51 * the last calculation. Otherwise, the calculation will be started from
52 * scratch and the initial composition and element potentials estimated.
53 */
54 bool contin;
55};
56
57/**
58 * @defgroup equil Chemical Equilibrium
59 */
60
61/**
62 * Class ChemEquil implements a chemical equilibrium solver for single-phase
63 * solutions. It is a "non-stoichiometric" solver in the terminology of Smith
64 * and Missen, meaning that every intermediate state is a valid chemical
65 * equilibrium state, but does not necessarily satisfy the element constraints.
66 * In contrast, the solver implemented in class MultiPhaseEquil uses a
67 * "stoichiometric" algorithm, in which each intermediate state satisfies the
68 * element constraints but is not a state of chemical equilibrium. Non-
69 * stoichiometric methods are faster when they converge, but stoichiometric ones
70 * tend to be more robust and can be used also for problems with multiple
71 * condensed phases. As expected, the ChemEquil solver is faster than
72 * MultiPhaseEquil for many single-phase equilibrium problems (particularly if
73 * there are only a few elements but very many species), but can be less stable.
74 * Problem situations include low temperatures where only a few species have
75 * non-zero mole fractions, precisely stoichiometric compositions (for example,
76 * 2 H2 + O2). In general, if speed is important, this solver should be tried first,
77 * and if it fails then use MultiPhaseEquil.
78 * @ingroup equil
79 */
81{
82public:
83 ChemEquil();
84
85 //! Constructor combined with the initialization function
86 /*!
87 * This constructor initializes the ChemEquil object with everything it
88 * needs to start solving equilibrium problems.
89 *
90 * @param s ThermoPhase object that will be used in the equilibrium calls.
91 */
93
94 virtual ~ChemEquil();
95
96 /*!
97 * Equilibrate a phase, holding the elemental composition fixed at the
98 * initial value found within the ThermoPhase object *s*.
99 *
100 * The value of two specified properties are obtained by querying the
101 * ThermoPhase object. The properties must be already contained within the
102 * current thermodynamic state of the system.
103 */
104 int equilibrate(ThermoPhase& s, const char* XY, int loglevel = 0);
105
106 /*!
107 * Compute the equilibrium composition for two specified properties and the
108 * specified element moles.
109 *
110 * The two specified properties are obtained by querying the ThermoPhase
111 * object. The properties must be already contained within the current
112 * thermodynamic state of the system.
113 *
114 * @param s phase object to be equilibrated
115 * @param XY property pair to hold constant
116 * @param elMoles specified vector of element abundances.
117 * @param loglevel Specify amount of debug logging (0 to disable)
118 * @return Successful returns are indicated by a return value of 0.
119 * Unsuccessful returns are indicated by a return value of -1 for lack
120 * of convergence or -3 for a singular Jacobian.
121 */
122 int equilibrate(ThermoPhase& s, const char* XY, vector_fp& elMoles,
123 int loglevel = 0);
124
125 /**
126 * Options controlling how the calculation is carried out.
127 * @see EquilOptions
128 */
130
131protected:
132 //! Pointer to the ThermoPhase object used to initialize this object.
133 /*!
134 * This ThermoPhase object must be compatible with the ThermoPhase objects
135 * input from the equilibrate function. Currently, this means that the 2
136 * ThermoPhases have to have consist of the same species and elements.
137 */
139
140 //! number of atoms of element m in species k.
141 doublereal nAtoms(size_t k, size_t m) const {
142 return m_comp[k*m_mm + m];
143 }
144
145 /*!
146 * Prepare for equilibrium calculations.
147 * @param s object representing the solution phase.
148 */
149 void initialize(ThermoPhase& s);
150
151 /*!
152 * Set mixture to an equilibrium state consistent with specified element
153 * potentials and temperature.
154 *
155 * @param s mixture to be updated
156 * @param x vector of non-dimensional element potentials
157 * \f[ \lambda_m/RT \f].
158 * @param t temperature in K.
159 */
161 const vector_fp& x, doublereal t);
162
163 //! Estimate the initial mole numbers. This version borrows from the
164 //! MultiPhaseEquil solver.
165 int setInitialMoles(ThermoPhase& s, vector_fp& elMoleGoal, int loglevel = 0);
166
167 //! Generate a starting estimate for the element potentials.
169 vector_fp& elMolesGoal, int loglevel = 0);
170
171 /*!
172 * Do a calculation of the element potentials using the Brinkley method,
173 * p. 129 Smith and Missen.
174 *
175 * We have found that the previous estimate may not be good enough to
176 * avoid drastic numerical issues associated with the use of a numerically
177 * generated Jacobian used in the main algorithm.
178 *
179 * The Brinkley algorithm, here, assumes a constant T, P system and uses a
180 * linearized analytical Jacobian that turns out to be very stable even
181 * given bad initial guesses.
182 *
183 * The pressure and temperature to be used are in the ThermoPhase object
184 * input into the routine.
185 *
186 * The initial guess for the element potentials used by this routine is
187 * taken from the input vector, x.
188 *
189 * elMoles is the input element abundance vector to be matched.
190 *
191 * Nonideal phases are handled in principle. This is done by calculating
192 * the activity coefficients and adding them into the formula in the
193 * correct position. However, these are treated as a RHS contribution
194 * only. Therefore, convergence might be a problem. This has not been
195 * tested. Also molality based unit systems aren't handled.
196 *
197 * On return, int return value contains the success code:
198 * - 0 - successful
199 * - 1 - unsuccessful, max num iterations exceeded
200 * - -3 - unsuccessful, singular Jacobian
201 *
202 * NOTE: update for activity coefficients.
203 */
204 int estimateEP_Brinkley(ThermoPhase& s, vector_fp& lambda, vector_fp& elMoles);
205
206 //! Find an acceptable step size and take it.
207 /*!
208 * The original implementation employed a line search technique that
209 * enforced a reduction in the norm of the residual at every successful
210 * step. Unfortunately, this method created false convergence errors near
211 * the end of a significant number of steps, usually special conditions
212 * where there were stoichiometric constraints.
213 *
214 * This new method just does a delta damping approach, based on limiting
215 * the jump in the dimensionless element potentials. Mole fractions are
216 * limited to a factor of 2 jump in the values from this method. Near
217 * convergence, the delta damping gets out of the way.
218 */
219 int dampStep(ThermoPhase& s, vector_fp& oldx,
220 double oldf, vector_fp& grad, vector_fp& step, vector_fp& x,
221 double& f, vector_fp& elmols, double xval, double yval);
222
223 /**
224 * Evaluates the residual vector F, of length #m_mm
225 */
226 void equilResidual(ThermoPhase& s, const vector_fp& x,
227 const vector_fp& elmtotal, vector_fp& resid,
228 double xval, double yval, int loglevel = 0);
229
230 void equilJacobian(ThermoPhase& s, vector_fp& x,
231 const vector_fp& elmols, DenseMatrix& jac,
232 double xval, double yval, int loglevel = 0);
233
234 void adjustEloc(ThermoPhase& s, vector_fp& elMolesGoal);
235
236 //! Update internally stored state information.
237 void update(const ThermoPhase& s);
238
239 /**
240 * Given a vector of dimensionless element abundances, this routine
241 * calculates the moles of the elements and the moles of the species.
242 *
243 * @param[in] x = current dimensionless element potentials..
244 */
245 double calcEmoles(ThermoPhase& s, vector_fp& x,
246 const double& n_t, const vector_fp& Xmol_i_calc,
247 vector_fp& eMolesCalc, vector_fp& n_i_calc,
248 double pressureConst);
249
250 size_t m_mm; //!< number of elements in the phase
251 size_t m_kk; //!< number of species in the phase
252 size_t m_skip;
253
254 //! This is equal to the rank of the stoichiometric coefficient matrix when
255 //! it is computed. It's initialized to #m_mm.
257
258 std::function<double(ThermoPhase&)> m_p1, m_p2;
259
260 //! Current value of the mole fractions in the single phase. length = #m_kk.
262
263 //! Current value of the sum of the element abundances given the current
264 //! element potentials.
266
267 //! Current value of the element mole fractions. Note these aren't the goal
268 //! element mole fractions.
270 vector_fp m_reswork;
271 vector_fp m_jwork1;
272 vector_fp m_jwork2;
273
274 //! Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];
276 doublereal m_temp, m_dens;
277 doublereal m_p0;
278
279 //! Index of the element id corresponding to the electric charge of each
280 //! species. Equal to -1 if there is no such element id.
281 size_t m_eloc;
282
283 vector_fp m_startSoln;
284
285 vector_fp m_grt;
286 vector_fp m_mu_RT;
287
288 //! Dimensionless values of the Gibbs free energy for the standard state of
289 //! each species, at the temperature and pressure of the solution (the star
290 //! standard state).
292 std::vector<size_t> m_component;
293
294 //! element fractional cutoff, below which the element will be zeroed.
296 bool m_doResPerturb;
297
298 std::vector<size_t> m_orderVectorElements;
299 std::vector<size_t> m_orderVectorSpecies;
300
301 //! Verbosity of printed output. No messages when m_loglevel == 0. More
302 //! output as level increases.
304};
305
306}
307
308#endif
Class ChemEquil implements a chemical equilibrium solver for single-phase solutions.
Definition: ChemEquil.h:81
double calcEmoles(ThermoPhase &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...
Definition: ChemEquil.cpp:808
vector_fp m_muSS_RT
Dimensionless values of the Gibbs free energy for the standard state of each species,...
Definition: ChemEquil.h:291
int estimateElementPotentials(ThermoPhase &s, vector_fp &lambda, vector_fp &elMolesGoal, int loglevel=0)
Generate a starting estimate for the element potentials.
Definition: ChemEquil.cpp:219
void setToEquilState(ThermoPhase &s, const vector_fp &x, doublereal t)
Definition: ChemEquil.cpp:133
int equilibrate(ThermoPhase &s, const char *XY, int loglevel=0)
Definition: ChemEquil.cpp:309
size_t m_kk
number of species in the phase
Definition: ChemEquil.h:251
int m_loglevel
Verbosity of printed output.
Definition: ChemEquil.h:303
size_t m_nComponents
This is equal to the rank of the stoichiometric coefficient matrix when it is computed.
Definition: ChemEquil.h:256
ThermoPhase * m_phase
Pointer to the ThermoPhase object used to initialize this object.
Definition: ChemEquil.h:138
void update(const ThermoPhase &s)
Update internally stored state information.
Definition: ChemEquil.cpp:154
int estimateEP_Brinkley(ThermoPhase &s, vector_fp &lambda, vector_fp &elMoles)
Definition: ChemEquil.cpp:844
size_t m_eloc
Index of the element id corresponding to the electric charge of each species.
Definition: ChemEquil.h:281
int setInitialMoles(ThermoPhase &s, vector_fp &elMoleGoal, int loglevel=0)
Estimate the initial mole numbers.
Definition: ChemEquil.cpp:183
int dampStep(ThermoPhase &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.
Definition: ChemEquil.cpp:690
void initialize(ThermoPhase &s)
Definition: ChemEquil.cpp:65
EquilOpt options
Options controlling how the calculation is carried out.
Definition: ChemEquil.h:129
vector_fp m_molefractions
Current value of the mole fractions in the single phase. length = m_kk.
Definition: ChemEquil.h:261
vector_fp m_comp
Storage of the element compositions. natom(k,m) = m_comp[k*m_mm+ m];.
Definition: ChemEquil.h:275
vector_fp m_elementmolefracs
Current value of the element mole fractions.
Definition: ChemEquil.h:269
doublereal nAtoms(size_t k, size_t m) const
number of atoms of element m in species k.
Definition: ChemEquil.h:141
double m_elemFracCutoff
element fractional cutoff, below which the element will be zeroed.
Definition: ChemEquil.h:295
doublereal m_elementTotalSum
Current value of the sum of the element abundances given the current element potentials.
Definition: ChemEquil.h:265
size_t m_mm
number of elements in the phase
Definition: ChemEquil.h:250
void equilResidual(ThermoPhase &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.
Definition: ChemEquil.cpp:729
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:51
Chemical equilibrium options.
Definition: ChemEquil.h:26
int iterations
Iteration counter.
Definition: ChemEquil.h:35
doublereal maxStepSize
Maximum step size.
Definition: ChemEquil.h:41
doublereal relTolerance
Relative tolerance.
Definition: ChemEquil.h:32
doublereal absElemTol
Abs Tol in element number.
Definition: ChemEquil.h:33
int maxIterations
Maximum number of iterations.
Definition: ChemEquil.h:34
int propertyPair
Property pair flag.
Definition: ChemEquil.h:47
bool contin
Continuation flag.
Definition: ChemEquil.h:54
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:102
This file contains definitions of constants, types and terms that are used in internal routines and a...
Namespace for the Cantera kernel.
Definition: AnyMap.h:29
std::vector< double > vector_fp
Turn on the use of stl vectors for the basic array type within cantera Vector of doubles.
Definition: ct_defs.h:184
int _equilflag(const char *xy)
map property strings to integers
Definition: ChemEquil.cpp:22