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