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