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