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