Cantera  2.5.1
MultiPhaseEquil.h
Go to the documentation of this file.
1 //! @file MultiPhaseEquil.h
2 
3 // This file is part of Cantera. See License.txt in the top-level directory or
4 // at https://cantera.org/license.txt for license and copyright information.
5 
6 #ifndef CT_MULTIPHASE_EQUIL
7 #define CT_MULTIPHASE_EQUIL
8 
9 #include "MultiPhase.h"
10 
11 namespace Cantera
12 {
13 
14 /*!
15  * Multiphase chemical equilibrium solver. Class MultiPhaseEquil is designed
16  * to be used to set a mixture containing one or more phases to a state of
17  * chemical equilibrium. It implements the VCS algorithm, described in Smith
18  * and Missen, "Chemical Reaction Equilibrium."
19  *
20  * This class only handles chemical equilibrium at a specified temperature and
21  * pressure. To compute equilibrium holding other properties fixed, it is
22  * necessary to iterate on T and P in an "outer" loop, until the specified
23  * properties have the desired values. This is done, for example, in method
24  * equilibrate of class MultiPhase.
25  *
26  * This class is primarily meant to be used internally by the equilibrate
27  * method of class MultiPhase, although there is no reason it cannot be used
28  * directly in application programs if desired.
29  *
30  * @ingroup equil
31  */
33 {
34 public:
35  //! Construct a multiphase equilibrium manager for a multiphase mixture.
36  //! @param mix Pointer to a multiphase mixture object.
37  //! @param start If true, the initial composition will be determined by a
38  //! linear Gibbs minimization, otherwise the initial mixture
39  //! composition will be used.
40  //! @param loglevel Desired level of debug printing. loglevel = 0 suppresses
41  //! printing. Higher values request more verbose logging output.
42  MultiPhaseEquil(MultiPhase* mix, bool start=true, int loglevel = 0);
43 
44  virtual ~MultiPhaseEquil() {}
45 
46  size_t constituent(size_t m) {
47  if (m < m_nel) {
48  return m_order[m];
49  } else {
50  return npos;
51  }
52  }
53 
54  void getStoichVector(size_t rxn, vector_fp& nu) {
55  nu.resize(m_nsp, 0.0);
56  if (rxn > nFree()) {
57  return;
58  }
59  for (size_t k = 0; k < m_nsp; k++) {
60  nu[m_order[k]] = m_N(k, rxn);
61  }
62  }
63 
64  int iterations() {
65  return m_iter;
66  }
67 
68  doublereal equilibrate(int XY, doublereal err = 1.0e-9,
69  int maxsteps = 1000, int loglevel=-99);
70  doublereal error();
71 
72  std::string reactionString(size_t j) {
73  return "";
74  }
75  void setInitialMixMoles(int loglevel = 0) {
76  setInitialMoles(loglevel);
77  finish();
78  }
79 
80  size_t componentIndex(size_t n) {
81  return m_species[m_order[n]];
82  }
83 
84  void reportCSV(const std::string& reportFile);
85 
86  double phaseMoles(size_t iph) const;
87 
88 protected:
89  //! This method finds a set of component species and a complete set of
90  //! formation reactions for the non-components in terms of the components.
91  //! In most cases, many different component sets are possible, and
92  //! therefore neither the components returned by this method nor the
93  //! formation reactions are unique. The algorithm used here is described
94  //! in Smith and Missen, Chemical Reaction Equilibrium Analysis.
95  //!
96  //! The component species are taken to be the first M species in array
97  //! 'species' that have linearly-independent compositions.
98  //!
99  //! @param order On entry, vector \a order should contain species index
100  //! numbers in the order of decreasing desirability as a component.
101  //! For example, if it is desired to choose the components from among
102  //! the major species, this array might list species index numbers in
103  //! decreasing order of mole fraction. If array 'species' does not
104  //! have length = nSpecies(), then the species will be considered as
105  //! candidates to be components in declaration order, beginning with
106  //! the first phase added.
107  void getComponents(const std::vector<size_t>& order);
108 
109  //! Estimate the initial mole numbers. This is done by running each
110  //! reaction as far forward or backward as possible, subject to the
111  //! constraint that all mole numbers remain non-negative. Reactions for
112  //! which \f$ \Delta \mu^0 \f$ are positive are run in reverse, and ones
113  //! for which it is negative are run in the forward direction. The end
114  //! result is equivalent to solving the linear programming problem of
115  //! minimizing the linear Gibbs function subject to the element and non-
116  //! negativity constraints.
117  int setInitialMoles(int loglevel = 0);
118 
119  void computeN();
120 
121  //! Take one step in composition, given the gradient of G at the starting
122  //! point, and a vector of reaction steps dxi.
123  doublereal stepComposition(int loglevel = 0);
124 
125  //! Re-arrange a vector of species properties in sorted form
126  //! (components first) into unsorted, sequential form.
127  void unsort(vector_fp& x);
128 
129  void step(doublereal omega, vector_fp& deltaN, int loglevel = 0);
130 
131  //! Compute the change in extent of reaction for each reaction.
132  doublereal computeReactionSteps(vector_fp& dxi);
133 
134  void updateMixMoles();
135 
136  //! Clean up the composition. The solution algorithm can leave some
137  //! species in stoichiometric condensed phases with very small negative
138  //! mole numbers. This method simply sets these to zero.
139  void finish();
140 
141  // moles of the species with sorted index ns
142  double moles(size_t ns) const {
143  return m_moles[m_order[ns]];
144  }
145  double& moles(size_t ns) {
146  return m_moles[m_order[ns]];
147  }
148  int solutionSpecies(size_t n) const {
149  return m_dsoln[m_order[n]];
150  }
151  bool isStoichPhase(size_t n) const {
152  return (m_dsoln[m_order[n]] == 0);
153  }
154  doublereal mu(size_t n) const {
155  return m_mu[m_species[m_order[n]]];
156  }
157  std::string speciesName(size_t n) const {
158  return
159  m_mix->speciesName(m_species[m_order[n]]);
160  }
161 
162  //! Number of degrees of freedom
163  size_t nFree() const {
164  return (m_nsp > m_nel) ? m_nsp - m_nel : 0;
165  }
166 
167  size_t m_nel_mix, m_nsp_mix;
168  size_t m_nel, m_nsp;
169  size_t m_eloc;
170  int m_iter;
171  MultiPhase* m_mix;
172  doublereal m_press, m_temp;
173  std::vector<size_t> m_order;
174  DenseMatrix m_N, m_A;
175  vector_fp m_work, m_work2, m_work3;
176  vector_fp m_moles, m_lastmoles, m_dxi;
177  vector_fp m_deltaG_RT, m_mu;
178  std::vector<bool> m_majorsp;
179  std::vector<size_t> m_sortindex;
180  vector_int m_lastsort;
181  vector_int m_dsoln;
182  vector_int m_incl_element, m_incl_species;
183 
184  // Vector of indices for species that are included in the calculation.
185  // This is used to exclude pure-phase species with invalid thermo data
186  std::vector<size_t> m_species;
187  std::vector<size_t> m_element;
188  std::vector<bool> m_solnrxn;
189  bool m_force;
190 };
191 
192 }
193 
194 #endif
Headers for the MultiPhase object that is used to set up multiphase equilibrium problems (see Equilfu...
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:51
void getComponents(const std::vector< size_t > &order)
This method finds a set of component species and a complete set of formation reactions for the non-co...
doublereal computeReactionSteps(vector_fp &dxi)
Compute the change in extent of reaction for each reaction.
void unsort(vector_fp &x)
Re-arrange a vector of species properties in sorted form (components first) into unsorted,...
void finish()
Clean up the composition.
size_t nFree() const
Number of degrees of freedom.
doublereal stepComposition(int loglevel=0)
Take one step in composition, given the gradient of G at the starting point, and a vector of reaction...
MultiPhaseEquil(MultiPhase *mix, bool start=true, int loglevel=0)
Construct a multiphase equilibrium manager for a multiphase mixture.
int setInitialMoles(int loglevel=0)
Estimate the initial mole numbers.
A class for multiphase mixtures.
Definition: MultiPhase.h:58
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:759
const size_t npos
index returned by functions to indicate "no position"
Definition: ct_defs.h:188
std::vector< int > vector_int
Vector of ints.
Definition: ct_defs.h:182
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:180
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264