Cantera  2.5.1
vcs_solve.h
Go to the documentation of this file.
1 /**
2  * @file vcs_solve.h Header file for the internal object that holds the vcs
3  * equilibrium problem (see Class \link Cantera::VCS_SOLVE
4  * VCS_SOLVE\endlink and \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 _VCS_SOLVE_H
11 #define _VCS_SOLVE_H
12 
13 /*
14 * Index of Symbols
15 * -------------------
16 * irxn -> refers to the species or rxn between the species and
17 * the components in the problem
18 * k -> refers to the species
19 * j -> refers to the element or component
20 *
21 * ### -> to be eliminated
22 */
23 
24 #include "cantera/base/ct_defs.h"
25 #include "cantera/equil/vcs_defs.h"
27 #include "cantera/base/Array.h"
28 
29 namespace Cantera
30 {
31 
32 class vcs_VolPhase;
33 class VCS_SPECIES_THERMO;
34 class VCS_COUNTERS;
35 class MultiPhase;
36 
37 //! This is the main structure used to hold the internal data
38 //! used in vcs_solve_TP(), and to solve TP systems.
39 /*!
40  * The indices of information in this structure may change when the species
41  * basis changes or when phases pop in and out of existence. Both of these
42  * operations change the species ordering.
43  */
44 class VCS_SOLVE
45 {
46 public:
47  //! Initialize the sizes within the VCS_SOLVE object
48  /*!
49  * This resizes all of the internal arrays within the object.
50  */
51  VCS_SOLVE(MultiPhase* mphase, int printLvl=0);
52 
53  ~VCS_SOLVE();
54 
55  //! Solve an equilibrium problem
56  /*!
57  * This is the main interface routine to the equilibrium solver
58  *
59  * @param ipr Printing of results
60  * ipr = 1 -> Print problem statement and final results to
61  * standard output
62  * 0 -> don't report on anything
63  * @param ip1 Printing of intermediate results
64  * IP1 = 1 -> Print intermediate results.
65  * @param maxit Maximum number of iterations for the algorithm
66  * @return nonzero value: failure to solve the problem at hand. zero :
67  * success
68  */
69  int vcs(int ipr, int ip1, int maxit);
70 
71  //! Main routine that solves for equilibrium at constant T and P using a
72  //! variant of the VCS method
73  /*!
74  * This is the main routine that solves for equilibrium at constant T and P
75  * using a variant of the VCS method. Nonideal phases can be accommodated as
76  * well.
77  *
78  * Any number of single-species phases and multi-species phases can be
79  * handled by the present version.
80  *
81  * @param print_lvl 1 -> Print results to standard output;
82  * 0 -> don't report on anything
83  * @param printDetails 1 -> Print intermediate results.
84  * @param maxit Maximum number of iterations for the algorithm
85  * @return
86  * * 0 = Equilibrium Achieved
87  * * 1 = Range space error encountered. The element abundance criteria
88  * are only partially satisfied. Specifically, the first NC= (number
89  * of components) conditions are satisfied. However, the full NE
90  * (number of elements) conditions are not satisfied. The
91  * equilibrium condition is returned.
92  * * -1 = Maximum number of iterations is exceeded. Convergence was
93  * not found.
94  */
95  int vcs_solve_TP(int print_lvl, int printDetails, int maxit);
96 
97  /*!
98  * We make decisions on the initial mole number, and major-minor status
99  * here. We also fix up the total moles in a phase.
100  *
101  * irxn = id of the noncomponent species formation reaction for the
102  * species to be added in.
103  *
104  * The algorithm proceeds to implement these decisions in the previous
105  * position of the species. Then, vcs_switch_pos is called to move the
106  * species into the last active species slot, incrementing the number of
107  * active species at the same time.
108  *
109  * This routine is responsible for the global data manipulation only.
110  */
111  void vcs_reinsert_deleted(size_t kspec);
112 
113  //! Choose the optimum species basis for the calculations
114  /*!
115  * This is done by choosing the species with the largest mole fraction not
116  * currently a linear combination of the previous components. Then,
117  * calculate the stoichiometric coefficient matrix for that basis.
118  *
119  * Rearranges the solution data to put the component data at the
120  * front of the species list.
121  *
122  * Then, calculates `m_stoichCoeffRxnMatrix(jcomp,irxn)` the formation
123  * reactions for all noncomponent species in the mechanism. Also
124  * calculates DNG(I) and DNL(I), the net mole change for each formation
125  * reaction. Also, initializes IR(I) to the default state.
126  *
127  * @param[in] doJustComponents If true, the #m_stoichCoeffRxnMatrix and
128  * #m_deltaMolNumPhase are not calculated.
129  * @param[in] aw Vector of mole fractions which will be used to
130  * construct an optimal basis from.
131  * @param[in] sa Gram-Schmidt orthog work space (nc in length) sa[j]
132  * @param[in] ss Gram-Schmidt orthog work space (nc in length) ss[j]
133  * @param[in] sm QR matrix work space (nc*ne in length) sm[i+j*ne]
134  * @param[in] test This is a small negative number dependent upon whether
135  * an estimate is supplied or not.
136  * @param[out] usedZeroedSpecies If true, then a species with a zero
137  * concentration was used as a component.
138  * The problem may be converged. Or, the
139  * problem may have a range space error and
140  * may not have a proper solution.
141  * @returns VCS_SUCCESS if everything went ok. Returns
142  * VCS_FAILED_CONVERGENCE if there is a problem.
143  *
144  * ### Internal Variables calculated by this routine:
145  *
146  * - #m_numComponents: Number of component species. This routine
147  * calculates the #m_numComponents species. It switches their positions
148  * in the species vector so that they occupy the first #m_numComponents
149  * spots in the species vector.
150  * - `m_stoichCoeffRxnMatrix(jcomp,irxn)` Stoichiometric coefficient
151  * matrix for the reaction mechanism expressed in Reduced Canonical
152  * Form. jcomp refers to the component number, and irxn refers to the
153  * irxn_th non-component species.
154  * - `m_deltaMolNumPhase(iphase,irxn)`: Change in the number of moles in
155  * phase, iphase, due to the noncomponent formation reaction, irxn.
156  * - `m_phaseParticipation(iphase,irxn)`: This is 1 if the phase, iphase,
157  * participates in the formation reaction, irxn, and zero otherwise.
158  */
159  int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[],
160  double ss[], double test, bool* const usedZeroedSpecies);
161 
162  //! Choose a species to test for the next component
163  /*!
164  * We make the choice based on testing (molNum[i] * spSize[i]) for its
165  * maximum value. Preference for single species phases is also made.
166  *
167  * @param molNum Mole number vector
168  * @param j index into molNum[] that indicates where the search will
169  * start from Previous successful components are swapped into the front
170  * of molNum[].
171  * @param n Length of molNum[]
172  */
173  size_t vcs_basisOptMax(const double* const molNum, const size_t j, const size_t n);
174 
175  //! Evaluate the species category for the indicated species
176  /*!
177  * All evaluations are done using the "old" version of the solution.
178  *
179  * @param kspec Species to be evaluated
180  * @returns the calculated species type
181  */
182  int vcs_species_type(const size_t kspec) const;
183 
184  //! This routine evaluates the species type for all species
185  /*!
186  * - #VCS_SPECIES_MAJOR: Major species
187  * - #VCS_SPECIES_MINOR: Minor species
188  * - #VCS_SPECIES_SMALLMS: The species lies in a multicomponent phase
189  * that exists. Its concentration is currently very low, necessitating
190  * a different method of calculation.
191  * - #VCS_SPECIES_ZEROEDMS: The species lies in a multicomponent phase
192  * which currently doesn't exist. Its concentration is currently zero.
193  * - #VCS_SPECIES_ZEROEDSS: Species lies in a single-species phase which
194  * is currently zeroed out.
195  * - #VCS_SPECIES_DELETED: Species has such a small mole fraction it is
196  * deleted even though its phase may possibly exist. The species is
197  * believed to have such a small mole fraction that it best to throw
198  * the calculation of it out. It will be added back in at the end of
199  * the calculation.
200  * - #VCS_SPECIES_INTERFACIALVOLTAGE: Species refers to an electron in
201  * the metal The unknown is equal to the interfacial voltage drop
202  * across the interface on the SHE (standard hydrogen electrode) scale
203  * (volts).
204  * - #VCS_SPECIES_ZEROEDPHASE: Species lies in a multicomponent phase
205  * that is zeroed atm and will stay deleted due to a choice from a
206  * higher level. These species will formally always have zero mole
207  * numbers in the solution vector.
208  * - #VCS_SPECIES_ACTIVEBUTZERO: The species lies in a multicomponent
209  * phase which currently does exist. Its concentration is currently
210  * identically zero, though the phase exists. Note, this is a temporary
211  * condition that exists at the start of an equilibrium problem. The
212  * species is soon "birthed" or "deleted".
213  * - #VCS_SPECIES_STOICHZERO: The species lies in a multicomponent phase
214  * which currently does exist. Its concentration is currently
215  * identically zero, though the phase exists. This is a permanent
216  * condition due to stoich constraints
217  */
219 
220  //! Calculate the dimensionless chemical potentials of all species or
221  //! of certain groups of species, at a fixed temperature and pressure.
222  /*!
223  * We calculate the dimensionless chemical potentials of all species or
224  * certain groups of species here, at a fixed temperature and pressure, for
225  * the input mole vector z[] in the parameter list. Nondimensionalization is
226  * achieved by division by RT.
227  *
228  * Note, for multispecies phases which are currently zeroed out, the
229  * chemical potential is filled out with the standard chemical potential.
230  *
231  * For species in multispecies phases whose concentration is zero, we need
232  * to set the mole fraction to a very low value. Its chemical potential is
233  * then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration
234  * to keep numbers positive.
235  *
236  * For formulas, see vcs_chemPotPhase().
237  *
238  * Handling of Small Species:
239  * ------------------------------
240  * As per the discussion above, for small species where the mole fraction
241  *
242  * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF
243  *
244  * The chemical potential is calculated as:
245  *
246  * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
247  *
248  * Species in the following categories are treated as "small species"
249  * - #VCS_SPECIES_DELETED
250  * - #VCS_SPECIES_ACTIVEBUTZERO
251  *
252  * For species in multispecies phases which are currently not active, the
253  * treatment is different. These species are in the following species
254  * categories:
255  * - #VCS_SPECIES_ZEROEDMS
256  * - #VCS_SPECIES_ZEROEDPHASE
257  *
258  * For these species, the `ln( ActCoeff[I] X[I])` term is dropped
259  * altogether. The following equation is used:
260  *
261  * m_feSpecies(I) = m_SSfeSpecies(I)
262  * + Charge[I] * Faraday_dim * phasePhi[iphase];
263  *
264  * Handling of "Species" Representing Interfacial Voltages
265  * ---------------------------------------------------------
266  *
267  * These species have species types of
268  * #VCS_SPECIES_TYPE_INTERFACIALVOLTAGE The chemical potentials for these
269  * "species" refer to electrons in metal electrodes. They have the
270  * following formula
271  *
272  * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT
273  *
274  * - `F` is Faraday's constant.
275  * - `R` = gas constant
276  * - `T` = temperature
277  * - `V` = potential of the interface = phi_electrode - phi_solution
278  *
279  * For these species, the solution vector unknown, z[I], is V, the phase
280  * voltage, in volts.
281  *
282  * @param ll Determine which group of species gets updated
283  * - `ll = 0`: Calculate for all species
284  * - `ll < 0`: calculate for components and for major non-components
285  * - `ll = 1`: calculate for components and for minor non-components
286  *
287  * @param lbot Restricts the calculation of the chemical potential
288  * to the species between LBOT <= i < LTOP. Usually
289  * LBOT and LTOP will be equal to 0 and MR, respectively.
290  * @param ltop Top value of the loops
291  *
292  * @param stateCalc Determines whether z is old or new or tentative:
293  * - 1: Use the tentative values for the total number of
294  * moles in the phases, i.e., use TG1 instead of TG etc.
295  * - 0: Use the base values of the total number of
296  * moles in each system.
297  *
298  * Also needed:
299  * ff : standard state chemical potentials. These are the
300  * chemical potentials of the standard states at
301  * the same T and P as the solution.
302  * tg : Total Number of moles in the phase.
303  */
304  void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop);
305 
306  //! This routine uploads the state of the system into all of the
307  //! vcs_VolumePhase objects in the current problem.
308  /*!
309  * @param stateCalc Determines where to get the mole numbers from.
310  * - VCS_STATECALC_OLD -> from m_molNumSpecies_old
311  * - VCS_STATECALC_NEW -> from m_molNumSpecies_new
312  */
313  void vcs_updateVP(const int stateCalc);
314 
315  //! Utility function that evaluates whether a phase can be popped
316  //! into existence
317  /*!
318  * A phase can be popped iff the stoichiometric coefficients for the
319  * component species, whose concentrations will be lowered during the
320  * process, are positive by at least a small degree.
321  *
322  * If one of the phase species is a zeroed component, then the phase can
323  * be popped if the component increases in mole number as the phase moles
324  * are increased.
325  *
326  * @param iphasePop id of the phase, which is currently zeroed,
327  * @returns true if the phase can come into existence and false otherwise.
328  */
329  bool vcs_popPhasePossible(const size_t iphasePop) const;
330 
331  //! Decision as to whether a phase pops back into existence
332  /*!
333  * @param phasePopPhaseIDs Vector containing the phase ids of the phases
334  * that will be popped this step.
335  * @returns the phase id of the phase that pops back into existence. Returns
336  * -1 if there are no phases
337  */
338  size_t vcs_popPhaseID(std::vector<size_t> &phasePopPhaseIDs);
339 
340  //! Calculates the deltas of the reactions due to phases popping
341  //! into existence
342  /*!
343  * Updates #m_deltaMolNumSpecies[irxn] : reaction adjustments, where irxn
344  * refers to the irxn'th species formation reaction. This adjustment is
345  * for species irxn + M, where M is the number of components.
346  *
347  * @param iphasePop Phase id of the phase that will come into existence
348  * @returns an int representing the status of the step
349  * - 0 : normal return
350  * - 1 : A single species phase species has been zeroed out
351  * in this routine. The species is a noncomponent
352  * - 2 : Same as one but, the zeroed species is a component.
353  * - 3 : Nothing was done because the phase couldn't be birthed
354  * because a needed component is zero.
355  */
356  int vcs_popPhaseRxnStepSizes(const size_t iphasePop);
357 
358  //! Calculates formation reaction step sizes.
359  /*!
360  * This is equation 6.4-16, p. 143 in Smith and Missen.
361  *
362  * Output
363  * -------
364  * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers to
365  * the irxn'th species formation reaction.
366  * This adjustment is for species irxn + M,
367  * where M is the number of components.
368  *
369  * Special branching occurs sometimes. This causes the component basis
370  * to be reevaluated
371  *
372  * @param forceComponentCalc integer flagging whether a component
373  * recalculation needs to be carried out.
374  * @param kSpecial species number of phase being zeroed.
375  * @returns an int representing which phase may need to be zeroed
376  */
377  size_t vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial);
378 
379  //! Calculates the total number of moles of species in all phases.
380  /*!
381  * Also updates the variable m_totalMolNum and Reconciles Phase existence
382  * flags with total moles in each phase.
383  */
384  double vcs_tmoles();
385 
386  void check_tmoles() const;
387 
388  //! This subroutine calculates reaction free energy changes for
389  //! all noncomponent formation reactions.
390  /*!
391  * Formation reactions are
392  * reactions which create each noncomponent species from the component
393  * species. `m_stoichCoeffRxnMatrix(jcomp,irxn)` are the stoichiometric
394  * coefficients for these reactions. A stoichiometric coefficient of
395  * one is assumed for species irxn in this reaction.
396  *
397  * @param L
398  * - `L < 0`: Calculate reactions corresponding to major noncomponent
399  * and zeroed species only
400  * - `L = 0`: Do all noncomponent reactions, i, between
401  * 0 <= i < irxnl
402  * - `L > 0`: Calculate reactions corresponding to minor noncomponent
403  * and zeroed species only
404  *
405  * @param doDeleted Do deleted species
406  * @param vcsState Calculate deltaG corresponding to either old or new
407  * free energies
408  * @param alterZeroedPhases boolean indicating whether we should
409  * add in a special section for zeroed phases.
410  *
411  * Note we special case one important issue. If the component has zero
412  * moles, then we do not allow deltaG < 0.0 for formation reactions which
413  * would lead to the loss of more of that same component. This dG < 0.0
414  * condition feeds back into the algorithm in several places, and leads
415  * to a infinite loop in at least one case.
416  */
417  void vcs_deltag(const int L, const bool doDeleted, const int vcsState,
418  const bool alterZeroedPhases = true);
419 
420  void vcs_printDeltaG(const int stateCalc);
421 
422  //! Swaps the indices for all of the global data for two species, k1
423  //! and k2.
424  /*!
425  * @param ifunc If true, switch the species data and the noncomponent
426  * reaction data. This must be called for a non-component species only.
427  * If false, switch the species data only. Typically, we use this option
428  * when determining the component species and at the end of the
429  * calculation, when we want to return unscrambled results. All rxn data
430  * will be out-of-date.
431  * @param k1 First species index
432  * @param k2 Second species index
433  */
434  void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2);
435 
436  //! Main program to test whether a deleted phase should be brought
437  //! back into existence
438  /*!
439  * @param iph Phase id of the deleted phase
440  */
441  double vcs_phaseStabilityTest(const size_t iph);
442 
443  //! Solve an equilibrium problem at a particular fixed temperature
444  //! and pressure
445  /*!
446  * The actual problem statement is assumed to be in the structure already.
447  * This is a wrapper around the solve_TP() function. In this wrapper, we
448  * calculate the standard state Gibbs free energies of the species
449  * and we decide whether to we need to use the initial guess algorithm.
450  *
451  * @param ipr = 1 -> Print results to standard output;
452  * 0 -> don't report on anything
453  * @param ip1 = 1 -> Print intermediate results;
454  * 0 -> Dont print any intermediate results
455  * @param maxit Maximum number of iterations for the algorithm
456  * @param T Value of the Temperature (Kelvin)
457  * @param pres Value of the Pressure
458  * @returns an integer representing the success of the algorithm
459  * * 0 = Equilibrium Achieved
460  * * 1 = Range space error encountered. The element abundance criteria are
461  * only partially satisfied. Specifically, the first NC= (number of
462  * components) conditions are satisfied. However, the full NE (number of
463  * elements) conditions are not satisfied. The equilibrium condition is
464  * returned.
465  * * -1 = Maximum number of iterations is exceeded. Convergence was not
466  * found.
467  */
468  int vcs_TP(int ipr, int ip1, int maxit, double T, double pres);
469 
470  /*!
471  * Evaluate the standard state free energies at the current temperature
472  * and pressure. Ideal gas pressure contribution is added in here.
473  *
474  * @param ipr 1 -> Print results to standard output; 0 -> don't report
475  * on anything
476  * @param ip1 1 -> Print intermediate results; 0 -> don't.
477  * @param Temp Temperature (Kelvin)
478  * @param pres Pressure (Pascal)
479  */
480  int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres);
481 
482  //! Initialize the chemical potential of single species phases
483  /*!
484  * For single species phases, initialize the chemical potential with the
485  * value of the standard state chemical potential. This value doesn't
486  * change during the calculation
487  */
488  void vcs_fePrep_TP();
489 
490  //! Calculation of the total volume and the partial molar volumes
491  /*!
492  * This function calculates the partial molar volume for all species, kspec,
493  * in the thermo problem at the temperature TKelvin and pressure, Pres, pres
494  * is in atm. And, it calculates the total volume of the combined system.
495  *
496  * @param[in] tkelvin Temperature in kelvin()
497  * @param[in] pres Pressure in Pascal
498  * @param[in] w w[] is the vector containing the current mole
499  * numbers in units of kmol.
500  * @param[out] volPM[] For species in all phase, the entries are the
501  * partial molar volumes units of M**3 / kmol.
502  * @returns the total volume of the entire system in units of m**3.
503  */
504  double vcs_VolTotal(const double tkelvin, const double pres,
505  const double w[], double volPM[]);
506 
507  //! This routine is mostly concerned with changing the private data to be
508  //! consistent with what's needed for solution. It is called one time for
509  //! each new problem structure definition.
510  /*!
511  * The problem structure refers to:
512  *
513  * - the number and identity of the species.
514  * - the formula matrix and thus the number of components.
515  * - the number and identity of the phases.
516  * - the equation of state
517  * - the method and parameters for determining the standard state
518  * - The method and parameters for determining the activity coefficients.
519  *
520  * Tasks:
521  * 1. Fill in the SSPhase[] array.
522  * 2. Check to see if any multispecies phases actually have only one
523  * species in that phase. If true, reassign that phase and species
524  * to be a single-species phase.
525  * 3. Determine the number of components in the problem if not already
526  * done so. During this process the order of the species is changed
527  * in the private data structure. All references to the species
528  * properties must employ the ind[] index vector.
529  * 4. Initialization of arrays to zero.
530  * 5. Check to see if the problem is well posed (If all the element
531  * abundances are zero, the algorithm will fail)
532  *
533  * @param printLvl Print level of the routine
534  * @return
535  * VCS_SUCCESS = everything went OK;
536  * VCS_PUB_BAD = There is an irreconcilable difference in the
537  * public data structure from when the problem was
538  * initially set up.
539  */
540  int vcs_prep(int printLvl);
541 
542  //! Rearrange the constraint equations represented by the Formula
543  //! Matrix so that the operational ones are in the front
544  /*!
545  * This subroutine handles the rearrangement of the constraint equations
546  * represented by the Formula Matrix. Rearrangement is only necessary when
547  * the number of components is less than the number of elements. For this
548  * case, some constraints can never be satisfied exactly, because the
549  * range space represented by the Formula Matrix of the components can't
550  * span the extra space. These constraints, which are out of the range
551  * space of the component Formula matrix entries, are migrated to the back
552  * of the Formula matrix.
553  *
554  * A prototypical example is an extra element column in FormulaMatrix[],
555  * which is identically zero. For example, let's say that argon is has an
556  * element column in FormulaMatrix[], but no species in the mechanism
557  * actually contains argon. Then, nc < ne. Also, without perturbation of
558  * FormulaMatrix[] vcs_basopt[] would produce a zero pivot because the
559  * matrix would be singular (unless the argon element column was already
560  * the last column of FormulaMatrix[].
561  *
562  * This routine borrows heavily from vcs_basopt's algorithm. It finds nc
563  * constraints which span the range space of the Component Formula matrix,
564  * and assigns them as the first nc components in the formula matrix. This
565  * guarantees that vcs_basopt[] has a nonsingular matrix to invert.
566  *
567  * Other Variables
568  * @param aw Mole fraction work space (ne in length)
569  * @param sa Gram-Schmidt orthog work space (ne in length)
570  * @param sm QR matrix work space (ne*ne in length)
571  * @param ss Gram-Schmidt orthog work space (ne in length)
572  */
573  int vcs_elem_rearrange(double* const aw, double* const sa,
574  double* const sm, double* const ss);
575 
576  //! Swaps the indices for all of the global data for two elements, ipos
577  //! and jpos.
578  /*!
579  * This function knows all of the element information with VCS_SOLVE, and
580  * can therefore switch element positions
581  *
582  * @param ipos first global element index
583  * @param jpos second global element index
584  */
585  void vcs_switch_elem_pos(size_t ipos, size_t jpos);
586 
587  //! Calculates the diagonal contribution to the Hessian due to
588  //! the dependence of the activity coefficients on the mole numbers.
589  /*!
590  * (See framemaker notes, Eqn. 20 - VCS Equations document)
591  *
592  * We allow the diagonal to be increased positively to any degree.
593  * We allow the diagonal to be decreased to 1/3 of the ideal solution
594  * value, but no more -> it must remain positive.
595  */
596  double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal);
597 
598  //! Calculates the diagonal contribution to the Hessian due to
599  //! the dependence of the activity coefficients on the mole numbers.
600  /*!
601  * (See framemaker notes, Eqn. 20 - VCS Equations document)
602  */
603  double vcs_Hessian_actCoeff_diag(size_t irxn);
604 
605  //! Recalculate all of the activity coefficients in all of the phases
606  //! based on input mole numbers
607  /*
608  * @param moleSpeciesVCS kmol of species to be used in the update.
609  *
610  * NOTE: This routine needs to be regulated.
611  */
612  void vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS);
613 
614  //! Print out a report on the state of the equilibrium problem to
615  //! standard output.
616  /*!
617  * @param iconv Indicator of convergence, to be printed out in the report:
618  * - 0 converged
619  * - 1 range space error
620  * - -1 not converged
621  */
622  int vcs_report(int iconv);
623 
624  //! Computes the current elemental abundances vector
625  /*!
626  * Computes the elemental abundances vector, m_elemAbundances[], and stores
627  * it back into the global structure
628  */
629  void vcs_elab();
630 
631  /*!
632  * Checks to see if the element abundances are in compliance. If they are,
633  * then TRUE is returned. If not, FALSE is returned. Note the number of
634  * constraints checked is usually equal to the number of components in the
635  * problem. This routine can check satisfaction of all of the constraints
636  * in the problem, which is equal to ne. However, the solver can't fix
637  * breakage of constraints above nc, because that nc is the range space by
638  * definition. Satisfaction of extra constraints would have had to occur
639  * in the problem specification.
640  *
641  * The constraints should be broken up into 2 sections. If a constraint
642  * involves a formula matrix with positive and negative signs, and eaSet =
643  * 0.0, then you can't expect that the sum will be zero. There may be
644  * roundoff that inhibits this. However, if the formula matrix is all of
645  * one sign, then this requires that all species with nonzero entries in
646  * the formula matrix be identically zero. We put this into the logic
647  * below.
648  *
649  * @param ibound 1: Checks constraints up to the number of elements;
650  * 0: Checks constraints up to the number of components.
651  */
652  bool vcs_elabcheck(int ibound);
653 
654  /*!
655  * Computes the elemental abundances vector for a single phase,
656  * elemAbundPhase[], and returns it through the argument list. The mole
657  * numbers of species are taken from the current value in
658  * m_molNumSpecies_old[].
659  */
660  void vcs_elabPhase(size_t iphase, double* const elemAbundPhase);
661 
662  /*!
663  * This subroutine corrects for element abundances. At the end of the
664  * subroutine, the total moles in all phases are recalculated again,
665  * because we have changed the number of moles in this routine.
666  *
667  * @param aa temporary work vector, length ne*ne
668  * @param x temporary work vector, length ne
669  *
670  * @returns
671  * - 0 = Nothing of significance happened, Element abundances were and
672  * still are good.
673  * - 1 = The solution changed significantly; The element abundances are
674  * now good.
675  * - 2 = The solution changed significantly, The element abundances are
676  * still bad.
677  * - 3 = The solution changed significantly, The element abundances are
678  * still bad and a component species got zeroed out.
679  *
680  * Internal data to be worked on::
681  * - ga Current element abundances
682  * - m_elemAbundancesGoal Required elemental abundances
683  * - m_molNumSpecies_old Current mole number of species.
684  * - m_formulaMatrix[][] Formula matrix of the species
685  * - ne Number of elements
686  * - nc Number of components.
687  *
688  * NOTES:
689  * This routine is turning out to be very problematic. There are
690  * lots of special cases and problems with zeroing out species.
691  *
692  * Still need to check out when we do loops over nc vs. ne.
693  */
694  int vcs_elcorr(double aa[], double x[]);
695 
696  //! Create an initial estimate of the solution to the thermodynamic
697  //! equilibrium problem.
698  /*!
699  * @return value indicates success:
700  * - 0: successful initial guess
701  * - -1: Unsuccessful initial guess; the elemental abundances aren't
702  * satisfied.
703  */
704  int vcs_inest_TP();
705 
706  //! Estimate the initial mole numbers by constrained linear programming
707  /*!
708  * This is done by running each reaction as far forward or backward as
709  * possible, subject to the constraint that all mole numbers remain non-
710  * negative. Reactions for which \f$ \Delta \mu^0 \f$ are positive are run
711  * in reverse, and ones for which it is negative are run in the forward
712  * direction. The end result is equivalent to solving the linear
713  * programming problem of minimizing the linear Gibbs function subject to
714  * the element and non-negativity constraints.
715  */
716  int vcs_setMolesLinProg();
717 
718  //! Calculate the total dimensionless Gibbs free energy
719  /*!
720  * Inert species are handled as if they had a standard free energy of
721  * zero. Note, for this algorithm this function should be monotonically
722  * decreasing.
723  */
724  double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles);
725 
726  //! Calculate the total dimensionless Gibbs free energy of a single phase
727  /*!
728  * Inert species are handled as if they had a standard free energy of
729  * zero and if they obeyed ideal solution/gas theory.
730  *
731  * @param iphase ID of the phase
732  * @param w Species mole number vector for all species
733  * @param fe vector of partial molar free energies of all of the
734  * species
735  */
736  double vcs_GibbsPhase(size_t iphase, const double* const w,
737  const double* const fe);
738 
739  //! Transfer the results of the equilibrium calculation back from VCS_SOLVE
740  void vcs_prob_update();
741 
742  //! Fully specify the problem to be solved
743  void vcs_prob_specifyFully();
744 
745 private:
746  //! Zero out the concentration of a species.
747  /*!
748  * Make sure to conserve elements and keep track of the total moles in all
749  * phases.
750  * - w[]
751  * - m_tPhaseMoles_old[]
752  *
753  * @param kspec Species index
754  * @return:
755  * 1: succeeded
756  * 0: failed.
757  */
758  int vcs_zero_species(const size_t kspec);
759 
760  //! Change a single species from active to inactive status
761  /*!
762  * Rearrange data when species is added or removed. The Lth species is
763  * moved to the back of the species vector. The back of the species
764  * vector is indicated by the value of MR, the current number of
765  * active species in the mechanism.
766  *
767  * @param kspec Species Index
768  * @return
769  * Returns 0 unless.
770  * The return is 1 when the current number of
771  * noncomponent species is equal to zero. A recheck of deleted species
772  * is carried out in the main code.
773  */
774  int vcs_delete_species(const size_t kspec);
775 
776  //! This routine handles the bookkeeping involved with the deletion of
777  //! multiphase phases from the problem.
778  /*!
779  * When they are deleted, all of their species become active species, even
780  * though their mole numbers are set to zero. The routine does not make the
781  * decision to eliminate multiphases.
782  *
783  * Note, species in phases with zero mole numbers are still considered
784  * active. Whether the phase pops back into existence or not is checked as
785  * part of the main iteration loop.
786  *
787  * @param iph Phase to be deleted
788  * @returns whether the operation was successful or not
789  */
790  bool vcs_delete_multiphase(const size_t iph);
791 
792  //! Change the concentration of a species by delta moles.
793  /*!
794  * Make sure to conserve elements and keep track of the total kmoles in
795  * all phases.
796  *
797  * @param kspec The species index
798  * @param delta_ptr pointer to the delta for the species. This may
799  * change during the calculation
800  * @return
801  * 1: succeeded without change of dx
802  * 0: Had to adjust dx, perhaps to zero, in order to do the delta.
803  */
804  int delta_species(const size_t kspec, double* const delta_ptr);
805 
806  //! Provide an estimate for the deleted species in phases that are not
807  //! zeroed out
808  /*!
809  * Try to add back in all deleted species. An estimate of the kmol numbers
810  * are obtained and the species is added back into the equation system, into
811  * the old state vector.
812  *
813  * This routine is called at the end of the calculation, just before
814  * returning to the user.
815  */
816  size_t vcs_add_all_deleted();
817 
818  //! Recheck deleted species in multispecies phases.
819  /*!
820  * We are checking the equation:
821  *
822  * sum_u = sum_j_comp [ sigma_i_j * u_j ]
823  * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
824  *
825  * by first evaluating:
826  *
827  * DG_i_O = u_i_O - sum_u.
828  *
829  * Then, if TL is zero, the phase pops into existence if DG_i_O < 0. Also,
830  * if the phase exists, then we check to see if the species can have a mole
831  * number larger than VCS_DELETE_SPECIES_CUTOFF (default value = 1.0E-32).
832  */
833  int vcs_recheck_deleted();
834 
835  //! Minor species alternative calculation
836  /*!
837  * This is based upon the following approximation: The mole fraction
838  * changes due to these reactions don't affect the mole numbers of the
839  * component species. Therefore the following approximation is valid for
840  * a small component of an ideal phase:
841  *
842  * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I))
843  *
844  * `m_deltaGRxn_old` contains the contribution from
845  *
846  * m_feSpecies_old(I) =
847  * m_SSfeSpecies(I) +
848  * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph))
849  * Thus,
850  *
851  * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I))
852  *
853  * Most of this section is mainly restricting the update to reasonable
854  * values. We restrict the update a factor of 1.0E10 up and 1.0E-10 down
855  * because we run into trouble with the addition operator due to roundoff if
856  * we go larger than ~1.0E15. Roundoff will then sometimes produce zero mole
857  * fractions.
858  *
859  * Note: This routine was generalized to incorporate nonideal phases and
860  * phases on the molality basis
861  *
862  * @param[in] kspec The current species and corresponding formation
863  * reaction number.
864  * @param[in] irxn The current species and corresponding formation
865  * reaction number.
866  * @param[out] do_delete: BOOLEAN which if true on return, then we
867  * branch to the section that deletes a species
868  * from the current set of active species.
869  */
870  double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete,
871  char* ANOTE=0) const;
872 
873  //! This routine optimizes the minimization of the total Gibbs free energy
874  //! by making sure the slope of the following functional stays negative:
875  /*!
876  * The slope of the following functional is equivalent to the slope of the
877  * total Gibbs free energy of the system:
878  *
879  * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] )
880  *
881  * along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0<al<1)
882  * such that the a parabola approximation to Gibbs(al) fit to the
883  * end points al = 0 and al = 1 is minimized.
884  * - s1 = slope of Gibbs function at al = 0, which is the previous
885  * solution = d(Gibbs)/d(al).
886  * - s2 = slope of Gibbs function at al = 1, which is the current
887  * solution = d(Gibbs)/d(al).
888  *
889  * Only if there has been an inflection point (i.e., s1 < 0 and s2 > 0),
890  * does this code section kick in. It finds the point on the parabola
891  * where the slope is equal to zero.
892  */
893  bool vcs_globStepDamp();
894 
895  //! Calculate the norm of a deltaGibbs free energy vector
896  /*!
897  * Positive DG for species which don't exist are ignored.
898  *
899  * @param dg Vector of local delta G's.
900  */
901  double l2normdg(double dg[]) const;
902 
903  void checkDelta1(double* const ds, double* const delTPhMoles, size_t kspec);
904 
905  //! Estimate equilibrium compositions
906  /*!
907  * Algorithm covered in a section of Smith and Missen's Book.
908  *
909  * Linear programming module is based on using dbolm.
910  *
911  * @param aw aw[i[ Mole fraction work space (ne in length)
912  * @param sa sa[j] = Gram-Schmidt orthog work space (ne in length)
913  * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length)
914  * @param ss ss[j] = Gram-Schmidt orthog work space (ne in length)
915  * @param test This is a small negative number.
916  */
917  void vcs_inest(double* const aw, double* const sa, double* const sm,
918  double* const ss, double test);
919 
920  //! Calculate the status of single species phases.
921  void vcs_SSPhase();
922 
923  //! Delete memory that isn't just resizable STL containers
924  /*!
925  * This gets called by the destructor or by InitSizes().
926  */
927  void vcs_delete_memory();
928 
929  //! Initialize the internal counters
930  /*!
931  * Initialize the internal counters containing the subroutine call
932  * values and times spent in the subroutines.
933  *
934  * ifunc = 0 Initialize only those counters appropriate for the top of
935  * vcs_solve_TP().
936  * = 1 Initialize all counters.
937  */
938  void vcs_counters_init(int ifunc);
939 
940  //! Create a report on the plog file containing timing and its information
941  /*!
942  * @param timing_print_lvl If 0, just report the iteration count. If larger
943  * than zero, report the timing information
944  */
945  void vcs_TCounters_report(int timing_print_lvl = 1);
946 
947  void vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc);
948 
949  void vcs_setFlagsVolPhase(const size_t iph, const bool upToDate, const int stateCalc);
950 
951  //! Update all underlying vcs_VolPhase objects
952  /*!
953  * Update the mole numbers and the phase voltages of all phases in the vcs
954  * problem
955  *
956  * @param stateCalc Location of the update (either VCS_STATECALC_NEW or
957  * VCS_STATECALC_OLD).
958  */
959  void vcs_updateMolNumVolPhases(const int stateCalc);
960 
961  // Helper functions used internally by vcs_solve_TP
962  int solve_tp_component_calc(bool& allMinorZeroedSpecies);
963  void solve_tp_inner(size_t& iti, size_t& it1, bool& uptodate_minors,
964  bool& allMinorZeroedSpecies, int& forceComponentCalc,
965  int& stage, bool printDetails, char* ANOTE);
966  void solve_tp_equilib_check(bool& allMinorZeroedSpecies, bool& uptodate_minors,
967  bool& giveUpOnElemAbund, int& solveFail,
968  size_t& iti, size_t& it1, int maxit,
969  int& stage, bool& lec);
970  void solve_tp_elem_abund_check(size_t& iti, int& stage, bool& lec,
971  bool& giveUpOnElemAbund,
972  int& finalElemAbundAttempts,
973  int& rangeErrorFound);
974 
975  // data used by vcs_solve_TP and it's helper functions
976  vector_fp m_sm;
977  vector_fp m_ss;
978  vector_fp m_sa;
979  vector_fp m_aw;
980  vector_fp m_wx;
981 
982 public:
983  //! Print level for print routines
985 
986  MultiPhase* m_mix;
987 
988  //! Print out the problem specification in all generality as it currently
989  //! exists in the VCS_SOLVE object
990  /*!
991  * @param print_lvl Parameter lvl for printing
992  * * 0 - no printing
993  * * 1 - all printing
994  */
995  void prob_report(int print_lvl);
996 
997  //! Add elements to the local element list
998  /*!
999  * This routine sorts through the elements defined in the vcs_VolPhase
1000  * object. It then adds the new elements to the VCS_SOLVE object, and creates
1001  * a global map, which is stored in the vcs_VolPhase object. Id and matching
1002  * of elements is done strictly via the element name, with case not
1003  * mattering.
1004  *
1005  * The routine also fills in the position of the element in the vcs_VolPhase
1006  * object's ElGlobalIndex field.
1007  *
1008  * @param volPhase Object containing the phase to be added. The elements in
1009  * this phase are parsed for addition to the global element list
1010  */
1011  void addPhaseElements(vcs_VolPhase* volPhase);
1012 
1013  //! This routines adds entries for the formula matrix for one species
1014  /*!
1015  * This routines adds entries for the formula matrix for this object for one
1016  * species
1017  *
1018  * This object also fills in the index filed, IndSpecies, within the
1019  * volPhase object.
1020  *
1021  * @param volPhase object containing the species
1022  * @param k Species number within the volPhase k
1023  * @param kT global Species number within this object
1024  *
1025  */
1026  size_t addOnePhaseSpecies(vcs_VolPhase* volPhase, size_t k, size_t kT);
1027 
1028  //! This routine resizes the number of elements in the VCS_SOLVE object by
1029  //! adding a new element to the end of the element list
1030  /*!
1031  * The element name is added. Formula vector entries ang element abundances
1032  * for the new element are set to zero.
1033  *
1034  * @param elNameNew New name of the element
1035  * @param elType Type of the element
1036  * @param elactive boolean indicating whether the element is active
1037  * @returns the index number of the new element
1038  */
1039  size_t addElement(const char* elNameNew, int elType, int elactive);
1040 
1041  void reportCSV(const std::string& reportFile);
1042 
1043  //! Total number of species in the problems
1044  size_t m_nsp;
1045 
1046  //! Number of element constraints in the problem
1047  size_t m_nelem;
1048 
1049  //! Number of components calculated for the problem
1051 
1052  //! Total number of non-component species in the problem
1053  size_t m_numRxnTot;
1054 
1055  //! Current number of species in the problems. Species can be deleted if
1056  //! they aren't stable under the current conditions
1058 
1059  //! Current number of non-component species in the problem. Species can be
1060  //! deleted if they aren't stable under the current conditions
1061  size_t m_numRxnRdc;
1062 
1063  //! Number of active species which are currently either treated as
1064  //! minor species
1066 
1067  //! Number of Phases in the problem
1068  size_t m_numPhases;
1069 
1070  //! Formula matrix for the problem
1071  /*!
1072  * FormulaMatrix(kspec,j) = Number of elements, j, in the kspec species
1073  *
1074  * Both element and species indices are swapped.
1075  */
1077 
1078  //! Stoichiometric coefficient matrix for the reaction mechanism expressed
1079  //! in Reduced Canonical Form.
1080  /*!
1081  * This is the stoichiometric coefficient matrix for the reaction which
1082  * forms species kspec from the component species. A stoichiometric
1083  * coefficient of one is assumed for the species kspec in this mechanism.
1084  *
1085  * NOTE: kspec = irxn + m_numComponents
1086  *
1087  * `m_stoichCoeffRxnMatrix(j,irxn)` : j refers to the component number, and
1088  * irxn refers to the irxn_th non-component species. The stoichiometric
1089  * coefficients multiplied by the Formula coefficients of the component
1090  * species add up to the negative value of the number of elements in the
1091  * species kspec.
1092  *
1093  * size = nelements0 x nspecies0
1094  */
1096 
1097  //! Absolute size of the stoichiometric coefficients
1098  /*!
1099  * scSize[irxn] = abs(Size) of the stoichiometric coefficients. These are
1100  * used to determine whether a given species should be
1101  * handled by the alt_min treatment or should be handled as a
1102  * major species.
1103  */
1105 
1106  //! total size of the species
1107  /*!
1108  * This is used as a multiplier to the mole number in figuring out which
1109  * species should be components.
1110  */
1112 
1113  //! Standard state chemical potentials for species K at the current
1114  //! temperature and pressure.
1115  /*!
1116  * The first NC entries are for components. The following NR entries are
1117  * for the current non-component species in the mechanism.
1118  */
1120 
1121  //! Free energy vector from the start of the current iteration
1122  /*!
1123  * The free energies are saved at the start of the current iteration.
1124  * Length = number of species
1125  */
1127 
1128  //! Dimensionless new free energy for all the species in the mechanism
1129  //! at the new tentative T, P, and mole numbers.
1130  /*!
1131  * The first NC entries are for components. The following
1132  * NR entries are for the current non-component species in the mechanism.
1133  * Length = number of species
1134  */
1136 
1137  //! Setting for whether to do an initial estimate
1138  /*!
1139  * Initial estimate:
1140  * 0 Do not estimate the solution at all. Use the
1141  * supplied mole numbers as is.
1142  * 1 Only do an estimate if the element abundances aren't satisfied.
1143  * -1 Force an estimate of the soln. Throw out the input mole numbers.
1144  */
1146 
1147  //! Total moles of the species
1148  /*!
1149  * Total number of moles of the kth species.
1150  * Length = Total number of species = m
1151  */
1153 
1154  //! Specifies the species unknown type
1155  /*!
1156  * There are two types. One is the straightforward species, with the mole
1157  * number w[k], as the unknown. The second is the an interfacial voltage
1158  * where w[k] refers to the interfacial voltage in volts.
1159  *
1160  * These species types correspond to metallic electrons corresponding to
1161  * electrodes. The voltage and other interfacial conditions sets up an
1162  * interfacial current, which is set to zero in this initial treatment.
1163  * Later we may have non-zero interfacial currents.
1164  */
1166 
1167  //! Change in the number of moles of phase, iphase, due to the
1168  //! noncomponent formation reaction, irxn, for species, k:
1169  /*!
1170  * m_deltaMolNumPhase(iphase,irxn) = k = nc + irxn
1171  */
1173 
1174  //! This is 1 if the phase, iphase, participates in the formation reaction
1175  //! irxn, and zero otherwise. PhaseParticipation(iphase,irxn)
1177 
1178  //! electric potential of the iph phase
1180 
1181  //! Tentative value of the mole number vector. It's also used to store the
1182  //! mole fraction vector.
1184 
1185  //! Delta G(irxn) for the noncomponent species in the mechanism.
1186  /*!
1187  * Computed by the subroutine deltaG. m_deltaGRxn is the free energy change
1188  * for the reaction which forms species K from the component species. This
1189  * vector has length equal to the number of noncomponent species in the
1190  * mechanism. It starts with the first current noncomponent species in the
1191  * mechanism.
1192  */
1194 
1195  //! Last deltag[irxn] from the previous step
1197 
1198  //! Last deltag[irxn] from the previous step with additions for
1199  //! possible births of zeroed phases.
1201 
1202  //! Temporary vector of Rxn DeltaG's
1203  /*!
1204  * This is used from time to time, for printing purposes
1205  */
1207 
1208  //! Reaction Adjustments for each species during the current step
1209  /*!
1210  * delta Moles for each species during the current step.
1211  * Length = number of species
1212  */
1214 
1215  //! Element abundances vector
1216  /*!
1217  * Vector of moles of each element actually in the solution vector. Except
1218  * for certain parts of the algorithm, this is a constant. Note other
1219  * constraint conditions are added to this vector. This is input from the
1220  * input file and is considered a constant from thereon. units = kmoles
1221  */
1223 
1224  //! Element abundances vector Goals
1225  /*!
1226  * Vector of moles of each element that are the goals of the simulation.
1227  * This is a constant in the problem. Note other constraint conditions are
1228  * added to this vector. This is input from the input file and is considered
1229  * a constant from thereon. units = kmoles
1230  */
1232 
1233  //! Total number of kmoles in all phases
1234  /*!
1235  * This number includes the inerts.
1236  * -> Don't use this except for scaling purposes
1237  */
1239 
1240  //! Total kmols of species in each phase
1241  /*!
1242  * This contains the total number of moles of species in each phase
1243  *
1244  * Length = number of phases
1245  */
1247 
1248  //! total kmols of species in each phase in the tentative soln vector
1249  /*!
1250  * This contains the total number of moles of species in each phase
1251  * in the tentative solution vector
1252  *
1253  * Length = number of phases
1254  */
1256 
1257  //! Temporary vector of length NPhase
1259 
1260  //! Temporary vector of length NPhase
1262 
1263  //! Change in the total moles in each phase
1264  /*!
1265  * Length number of phases.
1266  */
1268 
1269  //! Temperature (Kelvin)
1271 
1272  //! Pressure
1274 
1275  //! Total kmoles of inert to add to each phase
1276  /*!
1277  * TPhInertMoles[iph] = Total kmoles of inert to add to each phase
1278  * length = number of phases
1279  */
1281 
1282  //! Tolerance requirement for major species
1283  double m_tolmaj;
1284 
1285  //! Tolerance requirements for minor species
1286  double m_tolmin;
1287 
1288  //! Below this, major species aren't refined any more
1289  double m_tolmaj2;
1290 
1291  //! Below this, minor species aren't refined any more
1292  double m_tolmin2;
1293 
1294  //! Index vector that keeps track of the species vector rearrangement
1295  /*!
1296  * At the end of each run, the species vector and associated data gets put
1297  * back in the original order.
1298  *
1299  * Example
1300  *
1301  * k = m_speciesMapIndex[kspec]
1302  *
1303  * kspec = current order in the vcs_solve object
1304  * k = original order in the MultiPhase object
1305  */
1306  std::vector<size_t> m_speciesMapIndex;
1307 
1308  //! Index that keeps track of the index of the species within the local
1309  //! phase
1310  /*!
1311  * This returns the local index of the species within the phase. Its
1312  * argument is the global species index within the VCS problem.
1313  *
1314  * k = m_speciesLocalPhaseIndex[kspec]
1315  *
1316  * k varies between 0 and the nSpecies in the phase
1317  *
1318  * Length = number of species
1319  */
1320  std::vector<size_t> m_speciesLocalPhaseIndex;
1321 
1322  //! Index vector that keeps track of the rearrangement of the elements
1323  /*!
1324  * At the end of each run, the element vector and associated data gets
1325  * put back in the original order.
1326  *
1327  * Example
1328  *
1329  * e = m_elementMapIndex[eNum]
1330  * eNum = current order in the vcs_solve object
1331  * e = original order in the MultiPhase object
1332  */
1333  std::vector<size_t> m_elementMapIndex;
1334 
1335  //! Mapping between the species index for noncomponent species and the
1336  //! full species index.
1337  /*!
1338  * ir[irxn] = Mapping between the reaction index for noncomponent
1339  * formation reaction of a species and the full species
1340  * index.
1341  *
1342  * Initially set to a value of K = NC + I This vector has length equal to
1343  * number of noncomponent species in the mechanism. It starts with the
1344  * first current noncomponent species in the mechanism. kspec = ir[irxn]
1345  */
1346  std::vector<size_t> m_indexRxnToSpecies;
1347 
1348  //! Major -Minor status vector for the species in the problem
1349  /*!
1350  * The index for this is species. The reaction that this is referring to
1351  * is `kspec = irxn + m_numComponents`. For possible values and their
1352  * meanings, see vcs_evaluate_speciesType().
1353  */
1355 
1356  //! Mapping from the species number to the phase number
1357  std::vector<size_t> m_phaseID;
1358 
1359  //! Boolean indicating whether a species belongs to a single-species phase
1360  // vector<bool> can't be used here because it doesn't work with std::swap
1361  std::vector<char> m_SSPhase;
1362 
1363  //! Species string name for the kth species
1364  std::vector<std::string> m_speciesName;
1365 
1366  //! Vector of strings containing the element names
1367  /*!
1368  * ElName[j] = String containing element names
1369  */
1370  std::vector<std::string> m_elementName;
1371 
1372  //! Type of the element constraint
1373  /*!
1374  * * 0 - #VCS_ELEM_TYPE_ABSPOS Normal element that is positive or zero in
1375  * all species.
1376  * * 1 - #VCS_ELEM_TYPE_ELECTRONCHARGE element dof that corresponds to
1377  * the electronic charge DOF.
1378  * * 2 - #VCS_ELEM_TYPE_CHARGENEUTRALITY element dof that corresponds to
1379  * a required charge neutrality constraint on the phase. The element
1380  * abundance is always exactly zero.
1381  * * 3 - #VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may mean
1382  * that a species has neg 0 or pos value of that constraint (other than
1383  * charge)
1384  */
1386 
1387  //! Specifies whether an element constraint is active
1388  /*!
1389  * The default is true. Length = nelements
1390  */
1392 
1393  //! Array of Phase Structures. Length = number of phases.
1394  std::vector<std::unique_ptr<vcs_VolPhase>> m_VolPhaseList;
1395 
1396  //! specifies the activity convention of the phase containing the species
1397  /*!
1398  * * 0 = molar based
1399  * * 1 = molality based
1400  *
1401  * length = number of species
1402  */
1404 
1405  //! specifies the activity convention of the phase.
1406  /*!
1407  * * 0 = molar based
1408  * * 1 = molality based
1409  *
1410  * length = number of phases
1411  */
1413 
1414  //! specifies the ln(Mnaught) used to calculate the chemical potentials
1415  /*!
1416  * For molar based activity conventions this will be equal to 0.0.
1417  * length = number of species.
1418  */
1420 
1421  //! Molar-based Activity Coefficients for Species.
1422  //! Length = number of species
1424 
1425  //! Molar-based Activity Coefficients for Species based on old mole numbers
1426  /*!
1427  * These activity coefficients are based on the m_molNumSpecies_old
1428  * values Molar based activity coefficients. Length = number of species
1429  */
1431 
1432  //! Change in the log of the activity coefficient with respect to the mole
1433  //! number multiplied by the phase mole number
1434  /*!
1435  * size = nspecies x nspecies
1436  *
1437  * This is a temporary array that gets regenerated every time it's needed.
1438  * It is not swapped wrt species.
1439  */
1441 
1442  //! Molecular weight of each species
1443  /*!
1444  * units = kg/kmol. length = number of species.
1445  *
1446  * note: this is a candidate for removal. I don't think we use it.
1447  */
1449 
1450  //! Charge of each species. Length = number of species.
1452 
1453  std::vector<std::vector<size_t> > phasePopProblemLists_;
1454 
1455  //! Vector of pointers to thermo structures which identify the model
1456  //! and parameters for evaluating the thermodynamic functions for that
1457  //! particular species.
1458  /*!
1459  * SpeciesThermo[k] pointer to the thermo information for the kth species
1460  */
1461  std::vector<std::unique_ptr<VCS_SPECIES_THERMO>> m_speciesThermoList;
1462 
1463  //! Choice of Hessians
1464  /*!
1465  * If this is true, then we will use a better approximation to the Hessian
1466  * based on Jacobian of the ln(ActCoeff) with respect to mole numbers
1467  */
1469 
1470  //! Total volume of all phases. Units are m^3
1471  double m_totalVol;
1472 
1473  //! Partial molar volumes of the species
1474  /*!
1475  * units = mks (m^3/kmol)
1476  * Length = number of species
1477  */
1479 
1480  //! dimensionless value of Faraday's constant, F / RT (1/volt)
1482 
1483  //! Timing and iteration counters for the vcs object
1485 
1486  //! Debug printing lvl
1487  /*!
1488  * Levels correspond to the following guidelines
1489  * * 0 No printing at all
1490  * * 1 Serious warnings or fatal errors get one line
1491  * * 2 one line per each successful vcs package call
1492  * * 3 one line per every successful solve_TP calculation
1493  * * 4 one line for every successful operation -> solve_TP gets a summary report
1494  * * 5 each iteration in solve_TP gets a report with one line per species
1495  * * 6 Each decision in solve_TP gets a line per species in addition to 4
1496  * * 10 Additionally Hessian matrix is printed out
1497  */
1499 
1500  //! printing level of timing information
1501  /*!
1502  * * 1 allowing printing of timing
1503  * * 0 do not allow printing of timing -> everything is printed as a NA.
1504  */
1506 
1507  //! Disable printing of timing information. Used to generate consistent
1508  //! output for tests.
1509  static void disableTiming();
1510 
1511  friend class vcs_phaseStabilitySolve;
1512 };
1513 
1514 }
1515 #endif
Header file for class Cantera::Array2D.
A class for 2D arrays stored in column-major (Fortran-compatible) form.
Definition: Array.h:32
A class for multiphase mixtures.
Definition: MultiPhase.h:58
Class to keep track of time and iterations.
Definition: vcs_internal.h:36
This is the main structure used to hold the internal data used in vcs_solve_TP(), and to solve TP sys...
Definition: vcs_solve.h:45
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1068
int m_doEstimateEquil
Setting for whether to do an initial estimate.
Definition: vcs_solve.h:1145
VCS_SOLVE(MultiPhase *mphase, int printLvl=0)
Initialize the sizes within the VCS_SOLVE object.
Definition: vcs_solve.cpp:24
int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
Definition: vcs_TP.cpp:43
int vcs_inest_TP()
Create an initial estimate of the solution to the thermodynamic equilibrium problem.
Definition: vcs_inest.cpp:278
vector_fp m_phasePhi
electric potential of the iph phase
Definition: vcs_solve.h:1179
vector_fp m_scSize
Absolute size of the stoichiometric coefficients.
Definition: vcs_solve.h:1104
void vcs_deltag(const int L, const bool doDeleted, const int vcsState, const bool alterZeroedPhases=true)
This subroutine calculates reaction free energy changes for all noncomponent formation reactions.
double m_tolmaj2
Below this, major species aren't refined any more.
Definition: vcs_solve.h:1289
double vcs_VolTotal(const double tkelvin, const double pres, const double w[], double volPM[])
Calculation of the total volume and the partial molar volumes.
Definition: vcs_solve.cpp:659
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1364
int vcs_elem_rearrange(double *const aw, double *const sa, double *const sm, double *const ss)
Rearrange the constraint equations represented by the Formula Matrix so that the operational ones are...
Array2D m_np_dLnActCoeffdMolNum
Change in the log of the activity coefficient with respect to the mole number multiplied by the phase...
Definition: vcs_solve.h:1440
vector_fp m_TmpPhase2
Temporary vector of length NPhase.
Definition: vcs_solve.h:1261
void vcs_updateMolNumVolPhases(const int stateCalc)
Update all underlying vcs_VolPhase objects.
vector_int m_elType
Type of the element constraint.
Definition: vcs_solve.h:1385
size_t m_nelem
Number of element constraints in the problem.
Definition: vcs_solve.h:1047
int m_useActCoeffJac
Choice of Hessians.
Definition: vcs_solve.h:1468
void vcs_updateVP(const int stateCalc)
This routine uploads the state of the system into all of the vcs_VolumePhase objects in the current p...
double m_totalVol
Total volume of all phases. Units are m^3.
Definition: vcs_solve.h:1471
size_t vcs_add_all_deleted()
Provide an estimate for the deleted species in phases that are not zeroed out.
bool vcs_globStepDamp()
This routine optimizes the minimization of the total Gibbs free energy by making sure the slope of th...
vector_int m_speciesStatus
Major -Minor status vector for the species in the problem.
Definition: vcs_solve.h:1354
int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], double ss[], double test, bool *const usedZeroedSpecies)
Choose the optimum species basis for the calculations.
int vcs_prep(int printLvl)
This routine is mostly concerned with changing the private data to be consistent with what's needed f...
Definition: vcs_prep.cpp:50
int vcs_elcorr(double aa[], double x[])
Definition: vcs_elem.cpp:97
vector_int m_actConventionSpecies
specifies the activity convention of the phase containing the species
Definition: vcs_solve.h:1403
size_t vcs_basisOptMax(const double *const molNum, const size_t j, const size_t n)
Choose a species to test for the next component.
vector_fp m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1152
void vcs_TCounters_report(int timing_print_lvl=1)
Create a report on the plog file containing timing and its information.
Definition: vcs_report.cpp:298
Array2D m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise.
Definition: vcs_solve.h:1176
static void disableTiming()
Disable printing of timing information.
Definition: vcs_solve.cpp:673
double vcs_phaseStabilityTest(const size_t iph)
Main program to test whether a deleted phase should be brought back into existence.
vector_fp m_molNumSpecies_new
Tentative value of the mole number vector.
Definition: vcs_solve.h:1183
void vcs_inest(double *const aw, double *const sa, double *const sm, double *const ss, double test)
Estimate equilibrium compositions.
Definition: vcs_inest.cpp:19
bool vcs_popPhasePossible(const size_t iphasePop) const
Utility function that evaluates whether a phase can be popped into existence.
vector_fp m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1222
double l2normdg(double dg[]) const
Calculate the norm of a deltaGibbs free energy vector.
int delta_species(const size_t kspec, double *const delta_ptr)
Change the concentration of a species by delta moles.
vector_int m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1165
vector_fp m_deltaGRxn_tmp
Temporary vector of Rxn DeltaG's.
Definition: vcs_solve.h:1206
vector_fp m_deltaGRxn_old
Last deltag[irxn] from the previous step.
Definition: vcs_solve.h:1196
int vcs_TP(int ipr, int ip1, int maxit, double T, double pres)
Solve an equilibrium problem at a particular fixed temperature and pressure.
Definition: vcs_TP.cpp:11
double m_tolmin2
Below this, minor species aren't refined any more.
Definition: vcs_solve.h:1292
void vcs_elab()
Computes the current elemental abundances vector.
Definition: vcs_elem.cpp:17
void prob_report(int print_lvl)
Print out the problem specification in all generality as it currently exists in the VCS_SOLVE object.
Definition: vcs_prob.cpp:25
int m_timing_print_lvl
printing level of timing information
Definition: vcs_solve.h:1505
size_t vcs_popPhaseID(std::vector< size_t > &phasePopPhaseIDs)
Decision as to whether a phase pops back into existence.
int vcs_setMolesLinProg()
Estimate the initial mole numbers by constrained linear programming.
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles)
Calculate the total dimensionless Gibbs free energy.
Definition: vcs_Gibbs.cpp:14
void vcs_CalcLnActCoeffJac(const double *const moleSpeciesVCS)
Recalculate all of the activity coefficients in all of the phases based on input mole numbers.
Definition: vcs_rxnadj.cpp:345
size_t m_numRxnMinorZeroed
Number of active species which are currently either treated as minor species.
Definition: vcs_solve.h:1065
vector_fp m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
Definition: vcs_solve.h:1255
vector_fp m_feSpecies_new
Dimensionless new free energy for all the species in the mechanism at the new tentative T,...
Definition: vcs_solve.h:1135
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Definition: vcs_elem.cpp:85
vector_fp m_tPhaseMoles_old
Total kmols of species in each phase.
Definition: vcs_solve.h:1246
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1057
int vcs_delete_species(const size_t kspec)
Change a single species from active to inactive status.
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1498
void vcs_counters_init(int ifunc)
Initialize the internal counters.
Definition: vcs_solve.cpp:641
vector_fp m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1231
void vcs_reinsert_deleted(size_t kspec)
vector_fp TPhInertMoles
Total kmoles of inert to add to each phase.
Definition: vcs_solve.h:1280
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1357
std::vector< size_t > m_indexRxnToSpecies
Mapping between the species index for noncomponent species and the full species index.
Definition: vcs_solve.h:1346
Array2D m_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1076
vector_fp m_wtSpecies
Molecular weight of each species.
Definition: vcs_solve.h:1448
vector_fp m_deltaPhaseMoles
Change in the total moles in each phase.
Definition: vcs_solve.h:1267
void addPhaseElements(vcs_VolPhase *volPhase)
Add elements to the local element list.
Definition: vcs_prob.cpp:122
int vcs_solve_TP(int print_lvl, int printDetails, int maxit)
Main routine that solves for equilibrium at constant T and P using a variant of the VCS method.
std::vector< size_t > m_speciesLocalPhaseIndex
Index that keeps track of the index of the species within the local phase.
Definition: vcs_solve.h:1320
size_t addOnePhaseSpecies(vcs_VolPhase *volPhase, size_t k, size_t kT)
This routines adds entries for the formula matrix for one species.
Definition: vcs_prob.cpp:168
Array2D m_stoichCoeffRxnMatrix
Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form.
Definition: vcs_solve.h:1095
size_t m_numRxnRdc
Current number of non-component species in the problem.
Definition: vcs_solve.h:1061
VCS_COUNTERS * m_VCount
Timing and iteration counters for the vcs object.
Definition: vcs_solve.h:1484
int vcs(int ipr, int ip1, int maxit)
Solve an equilibrium problem.
Definition: vcs_solve.cpp:487
std::vector< std::unique_ptr< VCS_SPECIES_THERMO > > m_speciesThermoList
Vector of pointers to thermo structures which identify the model and parameters for evaluating the th...
Definition: vcs_solve.h:1461
double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
Definition: vcs_rxnadj.cpp:297
double vcs_GibbsPhase(size_t iphase, const double *const w, const double *const fe)
Calculate the total dimensionless Gibbs free energy of a single phase.
Definition: vcs_Gibbs.cpp:39
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1050
void vcs_prob_specifyFully()
Fully specify the problem to be solved.
Definition: vcs_solve.cpp:539
int vcs_species_type(const size_t kspec) const
Evaluate the species category for the indicated species.
double m_pressurePA
Pressure.
Definition: vcs_solve.h:1273
void vcs_SSPhase()
Calculate the status of single species phases.
Definition: vcs_prep.cpp:15
vector_fp m_chargeSpecies
Charge of each species. Length = number of species.
Definition: vcs_solve.h:1451
vector_fp m_feSpecies_old
Free energy vector from the start of the current iteration.
Definition: vcs_solve.h:1126
int vcs_report(int iconv)
Print out a report on the state of the equilibrium problem to standard output.
Definition: vcs_report.cpp:12
vector_int m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1391
double m_Faraday_dim
dimensionless value of Faraday's constant, F / RT (1/volt)
Definition: vcs_solve.h:1481
int vcs_recheck_deleted()
Recheck deleted species in multispecies phases.
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1361
double m_tolmin
Tolerance requirements for minor species.
Definition: vcs_solve.h:1286
size_t m_nsp
Total number of species in the problems.
Definition: vcs_solve.h:1044
std::vector< std::unique_ptr< vcs_VolPhase > > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1394
vector_fp m_spSize
total size of the species
Definition: vcs_solve.h:1111
bool vcs_elabcheck(int ibound)
Definition: vcs_elem.cpp:29
std::vector< size_t > m_elementMapIndex
Index vector that keeps track of the rearrangement of the elements.
Definition: vcs_solve.h:1333
bool vcs_delete_multiphase(const size_t iph)
This routine handles the bookkeeping involved with the deletion of multiphase phases from the problem...
void vcs_delete_memory()
Delete memory that isn't just resizable STL containers.
Definition: vcs_solve.cpp:476
double vcs_Hessian_actCoeff_diag(size_t irxn)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
Definition: vcs_rxnadj.cpp:315
int vcs_zero_species(const size_t kspec)
Zero out the concentration of a species.
void vcs_prob_update()
Transfer the results of the equilibrium calculation back from VCS_SOLVE.
Definition: vcs_solve.cpp:628
size_t vcs_RxnStepSizes(int &forceComponentCalc, size_t &kSpecial)
Calculates formation reaction step sizes.
Definition: vcs_rxnadj.cpp:18
vector_fp m_deltaMolNumSpecies
Reaction Adjustments for each species during the current step.
Definition: vcs_solve.h:1213
vector_fp m_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
Definition: vcs_solve.h:1119
Array2D m_deltaMolNumPhase
Change in the number of moles of phase, iphase, due to the noncomponent formation reaction,...
Definition: vcs_solve.h:1172
std::vector< size_t > m_speciesMapIndex
Index vector that keeps track of the species vector rearrangement.
Definition: vcs_solve.h:1306
void vcs_switch_elem_pos(size_t ipos, size_t jpos)
Swaps the indices for all of the global data for two elements, ipos and jpos.
double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool *do_delete, char *ANOTE=0) const
Minor species alternative calculation.
vector_fp m_actCoeffSpecies_old
Molar-based Activity Coefficients for Species based on old mole numbers.
Definition: vcs_solve.h:1430
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2)
Swaps the indices for all of the global data for two species, k1 and k2.
vector_fp m_PMVolumeSpecies
Partial molar volumes of the species.
Definition: vcs_solve.h:1478
vector_fp m_deltaGRxn_new
Delta G(irxn) for the noncomponent species in the mechanism.
Definition: vcs_solve.h:1193
int m_printLvl
Print level for print routines.
Definition: vcs_solve.h:984
void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop)
Calculate the dimensionless chemical potentials of all species or of certain groups of species,...
vector_fp m_lnMnaughtSpecies
specifies the ln(Mnaught) used to calculate the chemical potentials
Definition: vcs_solve.h:1419
vector_int m_phaseActConvention
specifies the activity convention of the phase.
Definition: vcs_solve.h:1412
int vcs_popPhaseRxnStepSizes(const size_t iphasePop)
Calculates the deltas of the reactions due to phases popping into existence.
size_t m_numRxnTot
Total number of non-component species in the problem.
Definition: vcs_solve.h:1053
void vcs_fePrep_TP()
Initialize the chemical potential of single species phases.
Definition: vcs_TP.cpp:57
vector_fp m_deltaGRxn_Deficient
Last deltag[irxn] from the previous step with additions for possible births of zeroed phases.
Definition: vcs_solve.h:1200
double m_temperature
Temperature (Kelvin)
Definition: vcs_solve.h:1270
vector_fp m_TmpPhase
Temporary vector of length NPhase.
Definition: vcs_solve.h:1258
double m_tolmaj
Tolerance requirement for major species.
Definition: vcs_solve.h:1283
vector_fp m_actCoeffSpecies_new
Molar-based Activity Coefficients for Species.
Definition: vcs_solve.h:1423
bool vcs_evaluate_speciesType()
This routine evaluates the species type for all species.
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1370
size_t addElement(const char *elNameNew, int elType, int elactive)
This routine resizes the number of elements in the VCS_SOLVE object by adding a new element to the en...
Definition: vcs_prob.cpp:148
double m_totalMolNum
Total number of kmoles in all phases.
Definition: vcs_solve.h:1238
Phase information and Phase calculations for vcs.
Definition: vcs_VolPhase.h:82
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
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
Defines and definitions within the vcs package.
Internal declarations for the VCSnonideal package.