Cantera  2.5.1
MultiPhase.h
Go to the documentation of this file.
1 /**
2  * @file MultiPhase.h
3  * Headers for the \link Cantera::MultiPhase MultiPhase\endlink
4  * object that is used to set up multiphase equilibrium problems (see \ref equilfunctions).
5  */
6 
7 // This file is part of Cantera. See License.txt in the top-level directory or
8 // at https://cantera.org/license.txt for license and copyright information.
9 
10 #ifndef CT_MULTIPHASE_H
11 #define CT_MULTIPHASE_H
12 
15 
16 namespace Cantera
17 {
18 
19 //! @defgroup equilfunctions
20 
21 //! A class for multiphase mixtures. The mixture can contain any
22 //! number of phases of any type.
23 /*!
24  * This object is the basic tool used by Cantera for use in Multiphase
25  * equilibrium calculations.
26  *
27  * It is a container for a set of phases. Each phase has a given number of
28  * kmoles. Therefore, MultiPhase may be considered an "extrinsic"
29  * thermodynamic object, in contrast to the ThermoPhase object, which is an
30  * "intrinsic" thermodynamic object.
31  *
32  * MultiPhase may be considered to be "upstream" of the ThermoPhase objects in
33  * the sense that setting a property within MultiPhase, such as temperature,
34  * pressure, or species mole number, affects the underlying ThermoPhase
35  * object, but not the other way around.
36  *
37  * All phases have the same temperature and pressure, and a specified number
38  * of moles for each phase. The phases do not need to have the same elements.
39  * For example, a mixture might consist of a gaseous phase with elements (H,
40  * C, O, N), a solid carbon phase containing only element C, etc. A master
41  * element set will be constructed for the mixture that is the intersection of
42  * the elements of each phase.
43  *
44  * Below, reference is made to global species and global elements. These refer
45  * to the collective species and elements encompassing all of the phases
46  * tracked by the object.
47  *
48  * The global element list kept by this object is an intersection of the
49  * element lists of all the phases that comprise the MultiPhase.
50  *
51  * The global species list kept by this object is a concatenated list of all
52  * of the species in all the phases that comprise the MultiPhase. The ordering
53  * of species is contiguous with respect to the phase id.
54  *
55  * @ingroup equilfunctions
56  */
58 {
59 public:
60  //! Constructor.
61  /*!
62  * The constructor takes no arguments, since phases are added using
63  * method addPhase().
64  */
65  MultiPhase();
66 
67  //! Destructor. Does nothing. Class MultiPhase does not take "ownership"
68  //! (i.e. responsibility for destroying) the phase objects.
69  virtual ~MultiPhase() {}
70 
71  //! Add a vector of phases to the mixture
72  /*!
73  * See the single addPhases command. This just does a bunch of phases
74  * at one time
75  * @param phases Vector of pointers to phases
76  * @param phaseMoles Vector of mole numbers in each phase (kmol)
77  */
78  void addPhases(std::vector<ThermoPhase*>& phases, const vector_fp& phaseMoles);
79 
80  //! Add all phases present in 'mix' to this mixture.
81  /*!
82  * @param mix Add all of the phases in another MultiPhase
83  * object to the current object.
84  */
85  void addPhases(MultiPhase& mix);
86 
87  //! Add a phase to the mixture.
88  /*!
89  * This function must be called before the init() function is called,
90  * which serves to freeze the MultiPhase.
91  *
92  * @param p pointer to the phase object
93  * @param moles total number of moles of all species in this phase
94  */
95  void addPhase(ThermoPhase* p, doublereal moles);
96 
97  //! Number of elements.
98  size_t nElements() const {
99  return m_nel;
100  }
101 
102  //! Check that the specified element index is in range.
103  //! Throws an exception if m is greater than nElements()-1
104  void checkElementIndex(size_t m) const;
105 
106  //! Check that an array size is at least nElements().
107  //! Throws an exception if mm is less than nElements(). Used before calls
108  //! which take an array pointer.
109  void checkElementArraySize(size_t mm) const;
110 
111  //! Returns the name of the global element *m*.
112  /*!
113  * @param m index of the global element
114  */
115  std::string elementName(size_t m) const;
116 
117  //! Returns the index of the element with name \a name.
118  /*!
119  * @param name String name of the global element
120  */
121  size_t elementIndex(const std::string& name) const;
122 
123  //! Number of species, summed over all phases.
124  size_t nSpecies() const {
125  return m_nsp;
126  }
127 
128  //! Check that the specified species index is in range.
129  //! Throws an exception if k is greater than nSpecies()-1
130  void checkSpeciesIndex(size_t k) const;
131 
132  //! Check that an array size is at least nSpecies().
133  //! Throws an exception if kk is less than nSpecies(). Used before calls
134  //! which take an array pointer.
135  void checkSpeciesArraySize(size_t kk) const;
136 
137  //! Name of species with global index \a kGlob
138  /*!
139  * @param kGlob global species index
140  */
141  std::string speciesName(const size_t kGlob) const;
142 
143  //! Returns the Number of atoms of global element \a mGlob in
144  //! global species \a kGlob.
145  /*!
146  * @param kGlob global species index
147  * @param mGlob global element index
148  * @returns the number of atoms.
149  */
150  doublereal nAtoms(const size_t kGlob, const size_t mGlob) const;
151 
152  /// Returns the global Species mole fractions.
153  /*!
154  * Write the array of species mole
155  * fractions into array \c x. The mole fractions are
156  * normalized to sum to one in each phase.
157  *
158  * @param x vector of mole fractions. Length = number of global species.
159  */
160  void getMoleFractions(doublereal* const x) const;
161 
162  //! Process phases and build atomic composition array.
163  /*!
164  * This method must be called after all phases are added, before doing
165  * anything else with the mixture. After init() has been called, no more
166  * phases may be added.
167  */
168  void init();
169 
170  //! Returns the name of the n'th phase
171  /*!
172  * @param iph phase Index
173  */
174  std::string phaseName(const size_t iph) const;
175 
176  //! Returns the index, given the phase name
177  /*!
178  * @param pName Name of the phase
179  * @returns the index. A value of -1 means the phase isn't in the object.
180  */
181  int phaseIndex(const std::string& pName) const;
182 
183  //! Return the number of moles in phase n.
184  /*!
185  * @param n Index of the phase.
186  */
187  doublereal phaseMoles(const size_t n) const;
188 
189  //! Set the number of moles of phase with index n.
190  /*!
191  * @param n Index of the phase
192  * @param moles Number of moles in the phase (kmol)
193  */
194  void setPhaseMoles(const size_t n, const doublereal moles);
195 
196  /// Return a reference to phase n.
197  /*!
198  * The state of phase n is also updated to match the state stored locally
199  * in the mixture object.
200  *
201  * @param n Phase Index
202  * @return Reference to the ThermoPhase object for the phase
203  */
204  thermo_t& phase(size_t n);
205 
206  //! Check that the specified phase index is in range
207  //! Throws an exception if m is greater than nPhases()
208  void checkPhaseIndex(size_t m) const;
209 
210  //! Check that an array size is at least nPhases()
211  //! Throws an exception if mm is less than nPhases(). Used before calls
212  //! which take an array pointer.
213  void checkPhaseArraySize(size_t mm) const;
214 
215  //! Returns the moles of global species \c k. units = kmol
216  /*!
217  * @param kGlob Global species index k
218  */
219  doublereal speciesMoles(size_t kGlob) const;
220 
221  //! Return the global index of the species belonging to phase number \c p
222  //! with local index \c k within the phase.
223  /*!
224  * @param k local index of the species within the phase
225  * @param p index of the phase
226  */
227  size_t speciesIndex(size_t k, size_t p) const {
228  return m_spstart[p] + k;
229  }
230 
231  //! Return the global index of the species belonging to phase name \c phaseName
232  //! with species name \c speciesName
233  /*!
234  * @param speciesName Species Name
235  * @param phaseName Phase Name
236  *
237  * @returns the global index
238  *
239  * If the species or phase name is not recognized, this routine throws a
240  * CanteraError.
241  */
242  size_t speciesIndex(const std::string& speciesName, const std::string& phaseName);
243 
244  /// Minimum temperature for which all solution phases have valid thermo
245  /// data. Stoichiometric phases are not considered, since they may have
246  /// thermo data only valid for conditions for which they are stable.
247  doublereal minTemp() const {
248  return m_Tmin;
249  }
250 
251  /// Maximum temperature for which all solution phases have valid thermo
252  /// data. Stoichiometric phases are not considered, since they may have
253  /// thermo data only valid for conditions for which they are stable.
254  doublereal maxTemp() const {
255  return m_Tmax;
256  }
257 
258  //! Total charge summed over all phases (Coulombs).
259  doublereal charge() const;
260 
261  /// Charge (Coulombs) of phase with index \a p.
262  /*!
263  * The net charge is computed as \f[ Q_p = N_p \sum_k F z_k X_k \f]
264  * where the sum runs only over species in phase \a p.
265  * @param p index of the phase for which the charge is desired.
266  */
267  doublereal phaseCharge(size_t p) const;
268 
269  //! Total moles of global element \a m, summed over all phases.
270  /*!
271  * @param m Index of the global element
272  */
273  doublereal elementMoles(size_t m) const;
274 
275  //! Returns a vector of Chemical potentials.
276  /*!
277  * Write into array \a mu the chemical potentials of all species
278  * [J/kmol]. The chemical potentials are related to the activities by
279  *
280  * \f$
281  * \mu_k = \mu_k^0(T, P) + RT \ln a_k.
282  * \f$.
283  *
284  * @param mu Chemical potential vector. Length = num global species. Units
285  * = J/kmol.
286  */
287  void getChemPotentials(doublereal* mu) const;
288 
289  /// Returns a vector of Valid chemical potentials.
290  /*!
291  * Write into array \a mu the chemical potentials of all species with
292  * thermo data valid for the current temperature [J/kmol]. For other
293  * species, set the chemical potential to the value \a not_mu. If \a
294  * standard is set to true, then the values returned are standard chemical
295  * potentials.
296  *
297  * This method is designed for use in computing chemical equilibrium by
298  * Gibbs minimization. For solution phases (more than one species), this
299  * does the same thing as getChemPotentials. But for stoichiometric
300  * phases, this writes into array \a mu the user-specified value \a not_mu
301  * instead of the chemical potential if the temperature is outside the
302  * range for which the thermo data for the one species in the phase are
303  * valid. The need for this arises since many condensed phases have thermo
304  * data fit only for the temperature range for which they are stable. For
305  * example, in the NASA database, the fits for H2O(s) are only done up to
306  * 0 C, the fits for H2O(L) are only done from 0 C to 100 C, etc. Using
307  * the polynomial fits outside the range for which the fits were done can
308  * result in spurious chemical potentials, and can lead to condensed
309  * phases appearing when in fact they should be absent.
310  *
311  * By setting \a not_mu to a large positive value, it is possible to force
312  * routines which seek to minimize the Gibbs free energy of the mixture to
313  * zero out any phases outside the temperature range for which their
314  * thermo data are valid.
315  *
316  * @param not_mu Value of the chemical potential to set species in phases,
317  * for which the thermo data is not valid
318  * @param mu Vector of chemical potentials. length = Global species,
319  * units = J kmol-1
320  * @param standard If this method is called with \a standard set to true,
321  * then the composition-independent standard chemical
322  * potentials are returned instead of the composition-
323  * dependent chemical potentials.
324  */
325  void getValidChemPotentials(doublereal not_mu, doublereal* mu,
326  bool standard = false) const;
327 
328  //! Temperature [K].
329  doublereal temperature() const {
330  return m_temp;
331  }
332 
333  //! Equilibrate a MultiPhase object
334  /*!
335  * Set this mixture to chemical equilibrium by calling one of Cantera's
336  * equilibrium solvers. The XY parameter indicates what two thermodynamic
337  * quantities are to be held constant during the equilibration process.
338  *
339  * @param XY String representation of what two properties are being
340  * held constant
341  * @param solver Name of the solver to be used to equilibrate the phase.
342  * If solver = 'vcs', the vcs_MultiPhaseEquil solver will be used. If
343  * solver = 'gibbs', the MultiPhaseEquil solver will be used. If solver
344  * = 'auto', the 'vcs' solver will be tried first, followed by the
345  * 'gibbs' solver if the first one fails.
346  * @param rtol Relative tolerance
347  * @param max_steps Maximum number of steps to take to find the solution
348  * @param max_iter The maximum number of outer temperature or pressure
349  * iterations to take when T and/or P is not held fixed.
350  * @param estimate_equil integer indicating whether the solver should
351  * estimate its own initial condition. If 0, the initial mole fraction
352  * vector in the ThermoPhase object is used as the initial condition.
353  * If 1, the initial mole fraction vector is used if the element
354  * abundances are satisfied. If -1, the initial mole fraction vector is
355  * thrown out, and an estimate is formulated.
356  * @param log_level loglevel Controls amount of diagnostic output.
357  * log_level=0 suppresses diagnostics, and increasingly-verbose
358  * messages are written as loglevel increases.
359  *
360  * @ingroup equilfunctions
361  */
362  void equilibrate(const std::string& XY, const std::string& solver="auto",
363  double rtol=1e-9, int max_steps=50000, int max_iter=100,
364  int estimate_equil=0, int log_level=0);
365 
366  /// Set the temperature [K].
367  /*!
368  * @param T value of the temperature (Kelvin)
369  */
370  void setTemperature(const doublereal T);
371 
372  //! Set the state of the underlying ThermoPhase objects in one call
373  /*!
374  * @param T Temperature of the system (kelvin)
375  * @param Pres pressure of the system (pascal)
376  */
377  void setState_TP(const doublereal T, const doublereal Pres);
378 
379  //! Set the state of the underlying ThermoPhase objects in one call
380  /*!
381  * @param T Temperature of the system (kelvin)
382  * @param Pres pressure of the system (pascal)
383  * @param Moles Vector of mole numbers of all the species in all the phases
384  * (kmol)
385  */
386  void setState_TPMoles(const doublereal T, const doublereal Pres, const doublereal* Moles);
387 
388  /// Pressure [Pa].
389  doublereal pressure() const {
390  return m_press;
391  }
392 
393  /// The total mixture volume [m^3].
394  /*!
395  * Returns the cumulative sum of the volumes of all the phases in the
396  * mixture.
397  */
398  doublereal volume() const;
399 
400  //! Set the pressure [Pa].
401  /*!
402  * @param P Set the pressure in the MultiPhase object (Pa)
403  */
404  void setPressure(doublereal P) {
405  m_press = P;
406  updatePhases();
407  }
408 
409  //! The enthalpy of the mixture [J].
410  doublereal enthalpy() const;
411 
412  //! The internal energy of the mixture [J].
413  doublereal IntEnergy() const;
414 
415  //! The entropy of the mixture [J/K].
416  doublereal entropy() const;
417 
418  //! The Gibbs function of the mixture [J].
419  doublereal gibbs() const;
420 
421  //! Heat capacity at constant pressure [J/K]. Note that this does not
422  //! account for changes in composition of the mixture with temperature.
423  doublereal cp() const;
424 
425  //! Number of phases.
426  size_t nPhases() const {
427  return m_phase.size();
428  }
429 
430  //! Return true is species \a kGlob is a species in a multicomponent
431  //! solution phase.
432  /*!
433  * @param kGlob index of the global species
434  */
435  bool solutionSpecies(size_t kGlob) const;
436 
437  //! Returns the phase index of the Kth "global" species
438  /*!
439  * @param kGlob Global species index.
440  * @returns the index of the owning phase.
441  */
442  size_t speciesPhaseIndex(const size_t kGlob) const;
443 
444  //! Returns the mole fraction of global species k
445  /*!
446  * @param kGlob Index of the global species.
447  */
448  doublereal moleFraction(const size_t kGlob) const;
449 
450  //! Set the Mole fractions of the nth phase
451  /*!
452  * This function sets the mole fractions of the nth phase. Note, the mole
453  * number of the phase stays constant
454  *
455  * @param n index of the phase
456  * @param x Vector of input mole fractions.
457  */
458  void setPhaseMoleFractions(const size_t n, const doublereal* const x);
459 
460  //! Set the number of moles of species in the mixture
461  /*!
462  * @param xMap CompositionMap of the species with nonzero mole numbers.
463  * Mole numbers that are less than or equal to zero will be
464  * set to zero. units = kmol.
465  */
466  void setMolesByName(const compositionMap& xMap);
467 
468  //! Set the moles via a string containing their names.
469  /*!
470  * The string x is in the form of a composition map. Species which are not
471  * listed are set to zero.
472  *
473  * @param x string x in the form of a composition map
474  * where values are the moles of the species.
475  */
476  void setMolesByName(const std::string& x);
477 
478  //! Get the mole numbers of all species in the multiphase object
479  /*!
480  * @param[out] molNum Vector of doubles of length nSpecies containing the
481  * global mole numbers (kmol).
482  */
483  void getMoles(doublereal* molNum) const;
484 
485  //! Sets all of the global species mole numbers
486  /*!
487  * The state of each phase object is also updated to have the specified
488  * composition and the mixture temperature and pressure.
489  *
490  * @param n Vector of doubles of length nSpecies containing the global
491  * mole numbers (kmol).
492  */
493  void setMoles(const doublereal* n);
494 
495  //! Adds moles of a certain species to the mixture
496  /*!
497  * @param indexS Index of the species in the MultiPhase object
498  * @param addedMoles Value of the moles that are added to the species.
499  */
500  void addSpeciesMoles(const int indexS, const doublereal addedMoles);
501 
502  //! Retrieves a vector of element abundances
503  /*!
504  * @param elemAbundances Vector of element abundances
505  * Length = number of elements in the MultiPhase object.
506  * Index is the global element index. Units is in kmol.
507  */
508  void getElemAbundances(doublereal* elemAbundances) const;
509 
510  //! Return true if the phase \a p has valid thermo data for the current
511  //! temperature.
512  /*!
513  * @param p Index of the phase.
514  */
515  bool tempOK(size_t p) const;
516 
517  // These methods are meant for internal use.
518 
519  //! Update the locally-stored composition within this object to match the
520  //! current compositions of the phase objects.
521  /*!
522  * Query the underlying ThermoPhase objects for their mole fractions and
523  * fill in the mole fraction vector of this current object. Adjust element
524  * compositions within this object to match.
525  *
526  * This is an upload operation in the sense that we are taking downstream
527  * information (ThermoPhase object info) and applying it to an upstream
528  * object (MultiPhase object).
529  */
531 
532  //! Set the states of the phase objects to the locally-stored
533  //! state within this MultiPhase object.
534  /*!
535  * This method sets each phase to the mixture temperature and pressure,
536  * and sets the phase mole fractions based on the mixture mole numbers.
537  *
538  * This is an download operation in the sense that we are taking upstream
539  * object information (MultiPhase object) and applying it to downstream
540  * objects (ThermoPhase object information)
541  *
542  * Therefore, the term, "update", is appropriate for a downstream operation.
543  */
544  void updatePhases() const;
545 
546 private:
547  //! Calculate the element abundance vector
548  void calcElemAbundances() const;
549 
550  //! Set the mixture to a state of chemical equilibrium using the
551  //! MultiPhaseEquil solver.
552  /*!
553  * @param XY Integer flag specifying properties to hold fixed.
554  * @param err Error tolerance for \f$\Delta \mu/RT \f$ for all reactions.
555  * Also used as the relative error tolerance for the outer loop.
556  * @param maxsteps Maximum number of steps to take in solving the fixed TP
557  * problem.
558  * @param maxiter Maximum number of "outer" iterations for problems holding
559  * fixed something other than (T,P).
560  * @param loglevel Level of diagnostic output
561  */
562  double equilibrate_MultiPhaseEquil(int XY, doublereal err, int maxsteps,
563  int maxiter, int loglevel);
564 
565  //! Vector of the number of moles in each phase.
566  /*!
567  * Length = m_np, number of phases.
568  */
570 
571  //! Vector of the ThermoPhase pointers.
572  std::vector<ThermoPhase*> m_phase;
573 
574  //! Global Stoichiometric Coefficient array
575  /*!
576  * This is a two dimensional array m_atoms(m, k). The first index is the
577  * global element index. The second index, k, is the global species index.
578  * The value is the number of atoms of type m in species k.
579  */
581 
582  //! Locally stored vector of mole fractions of all species comprising the
583  //! MultiPhase object.
585 
586  //! Mapping between the global species number and the phase ID
587  /*!
588  * m_spphase[kGlobal] = iPhase
589  * Length = number of global species
590  */
591  std::vector<size_t> m_spphase;
592 
593  //! Vector of ints containing of first species index in the global list of
594  //! species for each phase
595  /*!
596  * kfirst = m_spstart[ip], kfirst is the index of the first species in
597  * the ip'th phase.
598  */
599  std::vector<size_t> m_spstart;
600 
601  //! String names of the global elements. This has a length equal to the
602  //! number of global elements.
603  std::vector<std::string> m_enames;
604 
605  //! Atomic number of each global element.
607 
608  //! Vector of species names in the problem. Vector is over all species
609  //! defined in the object, the global species index.
610  std::vector<std::string> m_snames;
611 
612  //! Returns the global element index, given the element string name
613  /*!
614  * -> used in the construction. However, wonder if it needs to be global.
615  */
616  std::map<std::string, size_t> m_enamemap;
617 
618  //! Current value of the temperature (kelvin)
619  doublereal m_temp;
620 
621  //! Current value of the pressure (Pa)
622  doublereal m_press;
623 
624  //! Number of distinct elements in all of the phases
625  size_t m_nel;
626 
627  //! Number of distinct species in all of the phases
628  size_t m_nsp;
629 
630  //! True if the init() routine has been called, and the MultiPhase frozen
631  bool m_init;
632 
633  //! Global ID of the element corresponding to the electronic charge. If
634  //! there is none, then this is equal to -1
635  size_t m_eloc;
636 
637  //! Vector of bools indicating whether temperatures are ok for phases.
638  /*!
639  * If the current temperature is outside the range of valid temperatures
640  * for the phase thermodynamics, the phase flag is set to false.
641  */
642  mutable std::vector<bool> m_temp_OK;
643 
644  //! Minimum temperature for which thermo parameterizations are valid.
645  //! Stoichiometric phases are ignored in this determination. units Kelvin
646  doublereal m_Tmin;
647 
648  //! Minimum temperature for which thermo parameterizations are valid.
649  //! Stoichiometric phases are ignored in this determination. units Kelvin
650  doublereal m_Tmax;
651 
652  //! Vector of element abundances
653  /*!
654  * m_elemAbundances[mGlobal] = kmol of element mGlobal summed over all
655  * species in all phases.
656  */
658 };
659 
660 //! Function to output a MultiPhase description to a stream
661 /*!
662  * Writes out a description of the contents of each phase of the
663  * MultiPhase using the report function.
664  *
665  * @param s ostream
666  * @param x Reference to a MultiPhase
667  * @returns a reference to the ostream
668  */
669 inline std::ostream& operator<<(std::ostream& s, MultiPhase& x)
670 {
671  x.updatePhases();
672  for (size_t ip = 0; ip < x.nPhases(); ip++) {
673  if (x.phase(ip).name() != "") {
674  s << "*************** " << x.phase(ip).name() << " *****************" << std::endl;
675  } else {
676  s << "*************** Phase " << ip << " *****************" << std::endl;
677  }
678  s << "Moles: " << x.phaseMoles(ip) << std::endl;
679 
680  s << x.phase(ip).report() << std::endl;
681  }
682  return s;
683 }
684 
685 //! Choose the optimum basis of species for the equilibrium calculations.
686 /*!
687  * This is done by choosing the species with the largest mole fraction not
688  * currently a linear combination of the previous components. Then, calculate
689  * the stoichiometric coefficient matrix for that basis.
690  *
691  * Calculates the identity of the component species in the mechanism. Rearranges
692  * the solution data to put the component data at the front of the species list.
693  *
694  * Then, calculates SC(J,I) the formation reactions for all noncomponent
695  * species in the mechanism.
696  *
697  * @param[in] mphase Pointer to the multiphase object. Contains the species
698  * mole fractions, which are used to pick the current optimal species
699  * component basis.
700  * @param[in] orderVectorElements Order vector for the elements. The element
701  * rows in the formula matrix are rearranged according to this vector.
702  * @param[in] orderVectorSpecies Order vector for the species. The species are
703  * rearranged according to this formula. The first nCompoments of this
704  * vector contain the calculated species components on exit.
705  * @param[in] doFormRxn If true, the routine calculates the formation
706  * reaction matrix based on the calculated component species. If
707  * false, this step is skipped.
708  * @param[out] usedZeroedSpecies = If true, then a species with a zero
709  * concentration was used as a component. The problem may be converged.
710  * @param[out] formRxnMatrix
711  * @return The number of components.
712  *
713  * @ingroup equilfunctions
714  */
715 size_t BasisOptimize(int* usedZeroedSpecies, bool doFormRxn,
716  MultiPhase* mphase, std::vector<size_t>& orderVectorSpecies,
717  std::vector<size_t>& orderVectorElements,
718  vector_fp& formRxnMatrix);
719 
720 //! Handles the potential rearrangement of the constraint equations
721 //! represented by the Formula Matrix.
722 /*!
723  * Rearrangement is only necessary when the number of components is less
724  * than the number of elements. For this case, some constraints can never
725  * be satisfied exactly, because the range space represented by the Formula
726  * Matrix of the components can't span the extra space. These constraints,
727  * which are out of the range space of the component Formula matrix
728  * entries, are migrated to the back of the Formula matrix.
729  *
730  * A prototypical example is an extra element column in FormulaMatrix[], which
731  * is identically zero. For example, let's say that argon is has an element
732  * column in FormulaMatrix[], but no species in the mechanism actually
733  * contains argon. Then, nc < ne. Unless the entry for desired element
734  * abundance vector for Ar is zero, then this element abundance constraint can
735  * never be satisfied. The constraint vector is not in the range space of the
736  * formula matrix.
737  *
738  * Also, without perturbation of FormulaMatrix[], BasisOptimize[] would
739  * produce a zero pivot because the matrix would be singular (unless the argon
740  * element column was already the last column of FormulaMatrix[].
741  *
742  * This routine borrows heavily from BasisOptimize algorithm. It finds nc
743  * constraints which span the range space of the Component Formula matrix, and
744  * assigns them as the first nc components in the formula matrix. This
745  * guarantees that BasisOptimize has a nonsingular matrix to invert.
746  *
747  * @param[in] nComponents Number of components calculated previously.
748  * @param[in] elementAbundances Current value of the element abundances
749  * @param[in] mphase Input pointer to a MultiPhase object
750  * @param[in] orderVectorSpecies input vector containing the ordering of the
751  * global species in mphase. This is used to extract the component
752  * basis of the mphase object.
753  * @param[out] orderVectorElements Output vector containing the order of the
754  * elements that is necessary for calculation of the formula matrix.
755  *
756  * @ingroup equilfunctions
757  */
758 void ElemRearrange(size_t nComponents, const vector_fp& elementAbundances,
759  MultiPhase* mphase,
760  std::vector<size_t>& orderVectorSpecies,
761  std::vector<size_t>& orderVectorElements);
762 
763 //! External int that is used to turn on debug printing for the
764 //! BasisOptimze program.
765 /*!
766  * Set this to 1 if you want debug printing from BasisOptimize.
767  */
768 extern int BasisOptimize_print_lvl;
769 }
770 
771 #endif
Headers for the DenseMatrix object, which deals with dense rectangular matrices and description of th...
Header file for class ThermoPhase, the base class for phases with thermodynamic properties,...
A class for full (non-sparse) matrices with Fortran-compatible data storage, which adds matrix operat...
Definition: DenseMatrix.h:51
A class for multiphase mixtures.
Definition: MultiPhase.h:58
doublereal maxTemp() const
Maximum temperature for which all solution phases have valid thermo data.
Definition: MultiPhase.h:254
void init()
Process phases and build atomic composition array.
Definition: MultiPhase.cpp:116
void getMoleFractions(doublereal *const x) const
Returns the global Species mole fractions.
Definition: MultiPhase.cpp:769
void setState_TPMoles(const doublereal T, const doublereal Pres, const doublereal *Moles)
Set the state of the underlying ThermoPhase objects in one call.
Definition: MultiPhase.cpp:434
void checkPhaseIndex(size_t m) const
Check that the specified phase index is in range Throws an exception if m is greater than nPhases()
Definition: MultiPhase.cpp:170
size_t speciesIndex(size_t k, size_t p) const
Return the global index of the species belonging to phase number p with local index k within the phas...
Definition: MultiPhase.h:227
bool solutionSpecies(size_t kGlob) const
Return true is species kGlob is a species in a multicomponent solution phase.
Definition: MultiPhase.cpp:271
doublereal m_press
Current value of the pressure (Pa)
Definition: MultiPhase.h:622
void setMoles(const doublereal *n)
Sets all of the global species mole numbers.
Definition: MultiPhase.cpp:385
DenseMatrix m_atoms
Global Stoichiometric Coefficient array.
Definition: MultiPhase.h:580
void getMoles(doublereal *molNum) const
Get the mole numbers of all species in the multiphase object.
Definition: MultiPhase.cpp:370
void checkSpeciesIndex(size_t k) const
Check that the specified species index is in range.
Definition: MultiPhase.cpp:745
size_t m_nel
Number of distinct elements in all of the phases.
Definition: MultiPhase.h:625
vector_fp m_moleFractions
Locally stored vector of mole fractions of all species comprising the MultiPhase object.
Definition: MultiPhase.h:584
doublereal minTemp() const
Minimum temperature for which all solution phases have valid thermo data.
Definition: MultiPhase.h:247
double equilibrate_MultiPhaseEquil(int XY, doublereal err, int maxsteps, int maxiter, int loglevel)
Set the mixture to a state of chemical equilibrium using the MultiPhaseEquil solver.
Definition: MultiPhase.cpp:482
void calcElemAbundances() const
Calculate the element abundance vector.
Definition: MultiPhase.cpp:450
size_t nSpecies() const
Number of species, summed over all phases.
Definition: MultiPhase.h:124
doublereal gibbs() const
The Gibbs function of the mixture [J].
Definition: MultiPhase.cpp:280
doublereal cp() const
Heat capacity at constant pressure [J/K].
Definition: MultiPhase.cpp:328
doublereal pressure() const
Pressure [Pa].
Definition: MultiPhase.h:389
doublereal enthalpy() const
The enthalpy of the mixture [J].
Definition: MultiPhase.cpp:292
std::map< std::string, size_t > m_enamemap
Returns the global element index, given the element string name.
Definition: MultiPhase.h:616
doublereal speciesMoles(size_t kGlob) const
Returns the moles of global species k. units = kmol.
Definition: MultiPhase.cpp:184
void addSpeciesMoles(const int indexS, const doublereal addedMoles)
Adds moles of a certain species to the mixture.
Definition: MultiPhase.cpp:415
doublereal m_temp
Current value of the temperature (kelvin)
Definition: MultiPhase.h:619
doublereal IntEnergy() const
The internal energy of the mixture [J].
Definition: MultiPhase.cpp:304
void checkPhaseArraySize(size_t mm) const
Check that an array size is at least nPhases() Throws an exception if mm is less than nPhases().
Definition: MultiPhase.cpp:177
vector_fp m_elemAbundances
Vector of element abundances.
Definition: MultiPhase.h:657
std::vector< size_t > m_spstart
Vector of ints containing of first species index in the global list of species for each phase.
Definition: MultiPhase.h:599
doublereal m_Tmax
Minimum temperature for which thermo parameterizations are valid.
Definition: MultiPhase.h:650
doublereal m_Tmin
Minimum temperature for which thermo parameterizations are valid.
Definition: MultiPhase.h:646
vector_fp m_moles
Vector of the number of moles in each phase.
Definition: MultiPhase.h:569
size_t nPhases() const
Number of phases.
Definition: MultiPhase.h:426
size_t m_eloc
Global ID of the element corresponding to the electronic charge.
Definition: MultiPhase.h:635
doublereal moleFraction(const size_t kGlob) const
Returns the mole fraction of global species k.
Definition: MultiPhase.cpp:805
doublereal phaseMoles(const size_t n) const
Return the number of moles in phase n.
Definition: MultiPhase.cpp:790
std::vector< size_t > m_spphase
Mapping between the global species number and the phase ID.
Definition: MultiPhase.h:591
std::vector< std::string > m_snames
Vector of species names in the problem.
Definition: MultiPhase.h:610
std::string speciesName(const size_t kGlob) const
Name of species with global index kGlob.
Definition: MultiPhase.cpp:759
vector_int m_atomicNumber
Atomic number of each global element.
Definition: MultiPhase.h:606
size_t elementIndex(const std::string &name) const
Returns the index of the element with name name.
Definition: MultiPhase.cpp:735
bool tempOK(size_t p) const
Return true if the phase p has valid thermo data for the current temperature.
Definition: MultiPhase.cpp:810
std::string elementName(size_t m) const
Returns the name of the global element m.
Definition: MultiPhase.cpp:730
doublereal nAtoms(const size_t kGlob, const size_t mGlob) const
Returns the Number of atoms of global element mGlob in global species kGlob.
Definition: MultiPhase.cpp:764
void checkElementArraySize(size_t mm) const
Check that an array size is at least nElements().
Definition: MultiPhase.cpp:723
size_t speciesPhaseIndex(const size_t kGlob) const
Returns the phase index of the Kth "global" species.
Definition: MultiPhase.cpp:800
void addPhases(std::vector< ThermoPhase * > &phases, const vector_fp &phaseMoles)
Add a vector of phases to the mixture.
Definition: MultiPhase.cpp:40
std::vector< std::string > m_enames
String names of the global elements.
Definition: MultiPhase.h:603
void setPressure(doublereal P)
Set the pressure [Pa].
Definition: MultiPhase.h:404
doublereal elementMoles(size_t m) const
Total moles of global element m, summed over all phases.
Definition: MultiPhase.cpp:190
void setPhaseMoleFractions(const size_t n, const doublereal *const x)
Set the Mole fractions of the nth phase.
Definition: MultiPhase.cpp:340
doublereal charge() const
Total charge summed over all phases (Coulombs).
Definition: MultiPhase.cpp:205
std::string phaseName(const size_t iph) const
Returns the name of the n'th phase.
Definition: MultiPhase.cpp:774
void getValidChemPotentials(doublereal not_mu, doublereal *mu, bool standard=false) const
Returns a vector of Valid chemical potentials.
Definition: MultiPhase.cpp:251
void setState_TP(const doublereal T, const doublereal Pres)
Set the state of the underlying ThermoPhase objects in one call.
Definition: MultiPhase.cpp:424
doublereal entropy() const
The entropy of the mixture [J/K].
Definition: MultiPhase.cpp:316
void setMolesByName(const compositionMap &xMap)
Set the number of moles of species in the mixture.
Definition: MultiPhase.cpp:353
size_t m_nsp
Number of distinct species in all of the phases.
Definition: MultiPhase.h:628
void uploadMoleFractionsFromPhases()
Update the locally-stored composition within this object to match the current compositions of the pha...
Definition: MultiPhase.cpp:815
bool m_init
True if the init() routine has been called, and the MultiPhase frozen.
Definition: MultiPhase.h:631
doublereal volume() const
The total mixture volume [m^3].
Definition: MultiPhase.cpp:472
void updatePhases() const
Set the states of the phase objects to the locally-stored state within this MultiPhase object.
Definition: MultiPhase.cpp:826
void setPhaseMoles(const size_t n, const doublereal moles)
Set the number of moles of phase with index n.
Definition: MultiPhase.cpp:795
void checkSpeciesArraySize(size_t kk) const
Check that an array size is at least nSpecies().
Definition: MultiPhase.cpp:752
void getElemAbundances(doublereal *elemAbundances) const
Retrieves a vector of element abundances.
Definition: MultiPhase.cpp:442
doublereal temperature() const
Temperature [K].
Definition: MultiPhase.h:329
size_t nElements() const
Number of elements.
Definition: MultiPhase.h:98
virtual ~MultiPhase()
Destructor.
Definition: MultiPhase.h:69
std::vector< bool > m_temp_OK
Vector of bools indicating whether temperatures are ok for phases.
Definition: MultiPhase.h:642
int phaseIndex(const std::string &pName) const
Returns the index, given the phase name.
Definition: MultiPhase.cpp:780
std::vector< ThermoPhase * > m_phase
Vector of the ThermoPhase pointers.
Definition: MultiPhase.h:572
MultiPhase()
Constructor.
Definition: MultiPhase.cpp:21
void addPhase(ThermoPhase *p, doublereal moles)
Add a phase to the mixture.
Definition: MultiPhase.cpp:49
void getChemPotentials(doublereal *mu) const
Returns a vector of Chemical potentials.
Definition: MultiPhase.cpp:241
void checkElementIndex(size_t m) const
Check that the specified element index is in range.
Definition: MultiPhase.cpp:716
doublereal phaseCharge(size_t p) const
Charge (Coulombs) of phase with index p.
Definition: MultiPhase.cpp:230
thermo_t & phase(size_t n)
Return a reference to phase n.
Definition: MultiPhase.cpp:159
void setTemperature(const doublereal T)
Set the temperature [K].
Definition: MultiPhase.cpp:707
std::string name() const
Return the name of the phase.
Definition: Phase.cpp:84
Base class for a phase with thermodynamic properties.
Definition: ThermoPhase.h:102
virtual std::string report(bool show_thermo=true, doublereal threshold=-1e-14) const
returns a summary of the state of the phase as a string
size_t BasisOptimize(int *usedZeroedSpecies, bool doFormRxn, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements, vector_fp &formRxnMatrix)
Choose the optimum basis of species for the equilibrium calculations.
void ElemRearrange(size_t nComponents, const vector_fp &elementAbundances, MultiPhase *mphase, std::vector< size_t > &orderVectorSpecies, std::vector< size_t > &orderVectorElements)
Handles the potential rearrangement of the constraint equations represented by the Formula Matrix.
void equilibrate(const std::string &XY, const std::string &solver="auto", double rtol=1e-9, int max_steps=50000, int max_iter=100, int estimate_equil=0, int log_level=0)
Equilibrate a MultiPhase object.
Definition: MultiPhase.cpp:643
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
std::map< std::string, double > compositionMap
Map connecting a string name with a double.
Definition: ct_defs.h:172
Namespace for the Cantera kernel.
Definition: AnyMap.cpp:264
int BasisOptimize_print_lvl
External int that is used to turn on debug printing for the BasisOptimze program.