Cantera  2.1.2
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 VCSnonideal::VCS_SOLVE
4  * VCS_SOLVE\endlink and \ref equilfunctions ).
5  */
6 /*
7  * Copyright (2005) Sandia Corporation. Under the terms of
8  * Contract DE-AC04-94AL85000 with Sandia Corporation, the
9  * U.S. Government retains certain rights in this software.
10  */
11
12 #ifndef _VCS_SOLVE_H
13 #define _VCS_SOLVE_H
14
15 /*
16 * Index of Symbols
17 * -------------------
18 * irxn -> refers to the species or rxn between the species and
19 * the components in the problem
20 * k -> refers to the species
21 * j -> refers to the element or component
22 *
23 * ### -> to be eliminated
24 */
25
26 #include "cantera/base/ct_defs.h"
27 #include "cantera/equil/vcs_defs.h"
31
32 namespace VCSnonideal
33 {
34 /*
35  * Forward references
36  */
37 class vcs_VolPhase;
38 class VCS_SPECIES_THERMO;
39 class VCS_PROB;
40 class VCS_COUNTERS;
41
42
43 //! This is the main structure used to hold the internal data
44 //! used in vcs_solve_TP(), and to solve TP systems.
45 /*!
46  * The indices of information in this structure may change when the species
47  * basis changes or when phases pop in and out of existence. Both of these
48  * operations change the species ordering.
49  */
50 class VCS_SOLVE
51 {
52 public:
53  VCS_SOLVE();
54
55  ~VCS_SOLVE();
56
57  //! Initialize the sizes within the VCS_SOLVE object
58  /*!
59  * This resizes all of the internal arrays within the object. This
60  * routine operates in two modes. If all of the parameters are the same
61  * as it currently exists in the object, nothing is done by this routine;
62  * a quick exit is carried out and all of the data in the object
63  * persists.
64  *
65  * If any of the parameters are different than currently exists in the
66  * object, then all of the data in the object must be redone. It may not
67  * be zeroed, but it must be redone.
68  *
69  * @param nspecies0 Number of species within the object
70  * @param nelements Number of element constraints within the problem
71  * @param nphase0 Number of phases defined within the problem.
72  */
73  void vcs_initSizes(const size_t nspecies0, const size_t nelements, const size_t nphase0);
74
75  //! Solve an equilibrium problem
76  /*!
77  * This is the main interface routine to the equilibrium solver
78  *
79  * Input:
80  * @param vprob Object containing the equilibrium Problem statement
81  *
82  * @param ifunc Determines the operation to be done: Valid values:
83  * 0 -> Solve a new problem by initializing structures
84  * first. An initial estimate may or may not have
85  * been already determined. This is indicated in the
86  * VCS_PROB structure.
87  * 1 -> The problem has already been initialized and
88  * set up. We call this routine to resolve it
89  * using the problem statement and
90  * solution estimate contained in
91  * the VCS_PROB structure.
92  * 2 -> Don't solve a problem. Destroy all the private
93  * structures.
94  *
95  * @param ipr Printing of results
96  * ipr = 1 -> Print problem statement and final results to
97  * standard output
98  * 0 -> don't report on anything
99  * @param ip1 Printing of intermediate results
100  * IP1 = 1 -> Print intermediate results.
101  *
102  * @param maxit Maximum number of iterations for the algorithm
103  *
104  * Output:
105  *
106  * @return
107  * nonzero value: failure to solve the problem at hand.
108  * zero : success
109  */
110  int vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit);
111
112  //! Main routine that solves for equilibrium at constant T and P
113  //! using a variant of the VCS method
114  /*!
115  * This is the main routine that solves for equilibrium at constant T and P
116  * using a variant of the VCS method. Nonideal phases can be accommodated
117  * as well.
118  *
119  * Any number of single-species phases and multi-species phases
120  * can be handled by the present version.
121  *
122  * @param print_lvl 1 -> Print results to standard output;
123  * 0 -> don't report on anything
124  * @param printDetails 1 -> Print intermediate results.
125  * @param maxit Maximum number of iterations for the algorithm
126  *
127  * @return
128  * * 0 = Equilibrium Achieved
129  * * 1 = Range space error encountered. The element abundance criteria
130  * are only partially satisfied. Specifically, the first NC= (number
131  * of components) conditions are satisfied. However, the full NE
132  * (number of elements) conditions are not satisfied. The
133  * equilibrium condition is returned.
134  * * -1 = Maximum number of iterations is exceeded. Convergence was
136  */
137  int vcs_solve_TP(int print_lvl, int printDetails, int maxit);
138
139  int vcs_PS(VCS_PROB* vprob, int iph, int printLvl, double& feStable);
140
141  /*!
142  * We make decisions on the initial mole number, and major-minor status
143  * here. We also fix up the total moles in a phase.
144  *
145  * irxn = id of the noncomponent species formation reaction for the
146  * species to be added in.
147  *
148  * The algorithm proceeds to implement these decisions in the previous
149  * position of the species. Then, vcs_switch_pos is called to move the
150  * species into the last active species slot, incrementing the number
151  * of active species at the same time.
152  *
153  * This routine is responsible for the global data manipulation only.
154  */
155  void vcs_reinsert_deleted(size_t kspec);
156
157  //! Choose the optimum species basis for the calculations
158  /*!
159  * This is done by choosing the species with the largest mole fraction not
160  * currently a linear combination of the previous components. Then,
161  * calculate the stoichiometric coefficient matrix for that basis.
162  *
163  * Rearranges the solution data to put the component data at the
164  * front of the species list.
165  *
166  * Then, calculates m_stoichCoeffRxnMatrix[irxn][jcomp] the formation
167  * reactions for all noncomponent species in the mechanism. Also
168  * calculates DNG(I) and DNL(I), the net mole change for each formation
169  * reaction. Also, initializes IR(I) to the default state.
170  *
171  * @param[in] doJustComponents If true, the m_stoichCoeffRxnMatrix[][] and
172  * m_deltaMolNumPhase[] are not calculated.
173  *
174  * @param[in] aw Vector of mole fractions which will be used to construct an
175  * optimal basis from.
176  *
177  * @param[in] sa Gram-Schmidt orthog work space (nc in length) sa[j]
178  * @param[in] ss Gram-Schmidt orthog work space (nc in length) ss[j]
179  * @param[in] sm QR matrix work space (nc*ne in length) sm[i+j*ne]
180  * @param[in] test This is a small negative number dependent upon whether
181  * an estimate is supplied or not.
182  * @param[out] usedZeroedSpecies If true, then a species with a zero
183  * concentration was used as a component.
184  * The problem may be converged. Or, the
185  * problem may have a range space error and
186  * may not have a proper solution.
187  * @return Returns VCS_SUCCESS if everything went ok. Returns
188  * VCS_FAILED_CONVERGENCE if there is a problem.
189  *
190  * ### Internal Variables calculated by this routine:
191  *
192  * - #m_numComponents: Number of component species. This routine
193  * calculates the #m_numComponents species. It switches their positions
194  * in the species vector so that they occupy the first #m_numComponents
195  * spots in the species vector.
196  * - #m_stoichCoeffRxnMatrix[irxn][jcomp] Stoichiometric coefficient
197  * matrix for the reaction mechanism expressed in Reduced Canonical
198  * Form. jcomp refers to the component number, and irxn refers to the
199  * irxn_th non-component species.
200  * - #m_deltaMolNumPhase[irxn]: Change in the number of total number of
201  * moles of species in all phases due to the noncomponent formation
202  * reaction, irxn.
203  * - #m_deltaMolNumPhase[irxn][iphase]: Change in the number of moles in
204  * phase, iphase, due to the noncomponent formation reaction, irxn.
205  * - #m_phaseParticipation[irxn]: This is 1 if the phase, iphase,
206  * participates in the formation reaction, irxn, and zero otherwise.
207  */
208  int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[],
209  double ss[], double test, bool* const usedZeroedSpecies);
210
211  //! Choose a species to test for the next component
212  /*!
213  * We make the choice based on testing (molNum[i] * spSize[i]) for its
214  * maximum value. Preference for single species phases is also made.
215  *
216  * @param molNum Mole number vector
217  * @param j index into molNum[] that indicates where the search
218  * will start from Previous successful components are
219  * swapped into the front of molNum[].
220  * @param n Length of molNum[]
221  */
222  size_t vcs_basisOptMax(const double* const molNum, const size_t j, const size_t n);
223
224  //! Evaluate the species category for the indicated species
225  /*!
226  * All evaluations are done using the "old" version of the solution.
227  *
228  * @param kspec Species to be evaluated
229  *
230  * @return Returns the calculated species type
231  */
232  int vcs_species_type(const size_t kspec) const;
233
234  //! This routine evaluates the species type for all species
235  /*!
236  * - #VCS_SPECIES_MAJOR: Major species
237  * - #VCS_SPECIES_MINOR: Minor species
238  * - #VCS_SPECIES_SMALLMS: The species lies in a multicomponent phase
239  * that exists. Its concentration is currently very low, necessitating
240  * a different method of calculation.
241  * - #VCS_SPECIES_ZEROEDMS: The species lies in a multicomponent phase
242  * which currently doesn't exist. Its concentration is currently zero.
243  * - #VCS_SPECIES_ZEROEDSS: Species lies in a single-species phase which
244  * is currently zeroed out.
245  * - #VCS_SPECIES_DELETED: Species has such a small mole fraction it is
246  * deleted even though its phase may possibly exist. The species is
247  * believed to have such a small mole fraction that it best to throw
248  * the calculation of it out. It will be added back in at the end of
249  * the calculation.
250  * - #VCS_SPECIES_INTERFACIALVOLTAGE: Species refers to an electron in
251  * the metal The unknown is equal to the interfacial voltage drop
252  * across the interface on the SHE (standard hydrogen electrode) scale
253  * (volts).
254  * - #VCS_SPECIES_ZEROEDPHASE: Species lies in a multicomponent phase
255  * that is zeroed atm and will stay deleted due to a choice from a
256  * higher level. These species will formally always have zero mole
257  * numbers in the solution vector.
258  * - #VCS_SPECIES_ACTIVEBUTZERO: The species lies in a multicomponent
259  * phase which currently does exist. Its concentration is currently
260  * identically zero, though the phase exists. Note, this is a temporary
261  * condition that exists at the start of an equilibrium problem. The
262  * species is soon "birthed" or "deleted".
263  * - #VCS_SPECIES_STOICHZERO: The species lies in a multicomponent phase
264  * which currently does exist. Its concentration is currently
265  * identically zero, though the phase exists. This is a permanent
266  * condition due to stoich constraints
267  */
269
270  //! We calculate the dimensionless chemical potentials of all species
271  //! in a single phase.
272  /*!
273  * Note, for multispecies phases which are currently zeroed out, the
274  * chemical potential is filled out with the standard chemical potential.
275  *
276  * For species in multispecies phases whose concentration is zero, we need
277  * to set the mole fraction to a very low value. Its chemical potential is
278  * then calculated using the #VCS_DELETE_MINORSPECIES_CUTOFF concentration
279  * to keep numbers positive.
280  *
281  * # Formula:
282  *
283  * ## Ideal Mixtures:
284  *
285  * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph])
286  * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
287  *
288  * (This is equivalent to the adding the log of the mole fraction onto
289  * the standard chemical potential. )
290  *
291  * ## Non-Ideal Mixtures:
292  *
293  * ### ActivityConvention = 0: molarity activity formulation
294  *
295  * m_feSpecies(I) = m_SSfeSpecies(I)
296  * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
297  * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
298  *
299  * ( This is equivalent to the adding the log of the mole fraction
300  * multiplied by the activity coefficient onto the standard chemical
301  * potential. )
302  *
303  * ### ActivityConvention = 1: molality activity formulation
304  *
305  * m_feSpecies(I) = m_SSfeSpecies(I)
306  * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph])
307  * - ln(Mnaught * m_units)
308  * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
309  *
310  * Note: m_SSfeSpecies(I) is the molality based standard state. However,
311  * ActCoeff[I] is the molar based activity coefficient We have used the
312  * formulas:
313  *
314  * ActCoeff_M[I] = ActCoeff[I] / Xmol[N]
315  *
316  * where Xmol[N] is the mole fraction of the solvent and ActCoeff_M[I]
317  * is the molality based act coeff.
318  *
319  * Note: This is equivalent to the "normal" molality formulation:
320  *
321  * m_feSpecies(I) = m_SSfeSpecies(I)
322  * + ln(ActCoeff_M[I] * m(I))
323  * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]
324  *
325  * where m[I] is the molality of the ith solute
326  *
327  * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units)
328  *
329  * z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction of i in the phase.
330  *
331  *
332  * NOTE: As per the discussion in vcs_dfe(), for small species where the
333  * mole fraction is small:
334  *
335  * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF
336  *
337  * The chemical potential is calculated as:
338  *
339  * m_feSpecies(I) = m_SSfeSpecies(I)
340  * + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
341  *
342  * @param[in] iph Phase to be calculated
343  * @param[in] molNum(i) Number of moles of species i (VCS species order)
344  * @param[out] ac Activity coefficients for species in phase (VCS
345  * species order)
346  * @param[out] mu_i Dimensionless chemical potentials for phase
347  * species (VCS species order)
348  */
349  void vcs_chemPotPhase(const int stateCalc, const size_t iph, const double* const molNum,
350  double* const ac, double* const mu_i,
351  const bool do_deleted = false);
352
353  //! Calculate the dimensionless chemical potentials of all species or
354  //! of certain groups of species, at a fixed temperature and pressure.
355  /*!
356  * We calculate the dimensionless chemical potentials of all species
357  * or certain groups of species here, at a fixed temperature and pressure,
358  * for the input mole vector z[] in the parameter list.
359  * Nondimensionalization is achieved by division by RT.
360  *
361  * Note, for multispecies phases which are currently zeroed out,
362  * the chemical potential is filled out with the standard chemical
363  * potential.
364  *
365  * For species in multispecies phases whose concentration is zero, we need
366  * to set the mole fraction to a very low value. Its chemical potential is
367  * then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration
368  * to keep numbers positive.
369  *
370  * For formulas, see vcs_chemPotPhase().
371  *
372  * Handling of Small Species:
373  * ------------------------------
374  * As per the discussion above, for small species where the mole fraction
375  *
376  * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF
377  *
378  * The chemical potential is calculated as:
379  *
380  * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF))
381  *
382  * Species in the following categories are treated as "small species"
383  * - #VCS_SPECIES_DELETED
384  * - #VCS_SPECIES_ACTIVEBUTZERO
385  *
386  * For species in multispecies phases which are currently not active, the
387  * treatment is different. These species are in the following species
388  * categories:
389  * - #VCS_SPECIES_ZEROEDMS
390  * - #VCS_SPECIES_ZEROEDPHASE
391  *
392  * For these species, the ln( ActCoeff[I] X[I]) term is dropped
393  * altogether. The following equation is used:
394  *
395  * m_feSpecies(I) = m_SSfeSpecies(I)
396  * + Charge[I] * Faraday_dim * phasePhi[iphase];
397  *
398  * Handling of "Species" Representing Interfacial Voltages
399  * ---------------------------------------------------------
400  *
401  * These species have species types of
402  * #VCS_SPECIES_TYPE_INTERFACIALVOLTAGE The chemical potentials for these
403  * "species" refer to electrons in metal electrodes. They have the
404  * following formula
405  *
406  * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT
407  *
408  * - F is Faraday's constant.
409  * - R = gas constant
410  * - T = temperature
411  * - V = potential of the interface = phi_electrode - phi_solution
412  *
413  * For these species, the solution vector unknown, z[I], is V, the phase
414  * voltage, in volts.
415  *
416  * @param ll Determine which group of species gets updated
417  * - ll = 0: Calculate for all species
418  * - ll < 0: calculate for components and for major non-components
419  * - ll = 1: calculate for components and for minor non-components
420  *
421  * @param lbot Restricts the calculation of the chemical potential
422  * to the species between LBOT <= i < LTOP. Usually
423  * LBOT and LTOP will be equal to 0 and MR, respectively.
424  * @param ltop Top value of the loops
425  *
426  * @param stateCalc Determines whether z is old or new or tentative:
427  * - 1: Use the tentative values for the total number of
428  * moles in the phases, i.e., use TG1 instead of TG etc.
429  * - 0: Use the base values of the total number of
430  * moles in each system.
431  *
432  * Also needed:
433  * ff : standard state chemical potentials. These are the
434  * chemical potentials of the standard states at
435  * the same T and P as the solution.
436  * tg : Total Number of moles in the phase.
437  */
438  void vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop);
439
440  //! Print out a table of chemical potentials
441  /*!
442  * @param stateCalc Determines where to get the mole numbers from.
443  * - VCS_STATECALC_OLD -> from m_molNumSpecies_old
444  * - VCS_STATECALC_NEW -> from m_molNumSpecies_new
445  */
446  void vcs_printSpeciesChemPot(const int stateCalc) const;
447
448  //! This routine uploads the state of the system into all of the
449  //! vcs_VolumePhase objects in the current problem.
450  /*!
451  * @param stateCalc Determines where to get the mole numbers from.
452  * - VCS_STATECALC_OLD -> from m_molNumSpecies_old
453  * - VCS_STATECALC_NEW -> from m_molNumSpecies_new
454  */
455  void vcs_updateVP(const int stateCalc);
456
457  //! Utility function that evaluates whether a phase can be popped
458  //! into existence
459  /*!
460  * A phase can be popped iff the stoichiometric coefficients for the
461  * component species, whose concentrations will be lowered during the
462  * process, are positive by at least a small degree.
463  *
464  * If one of the phase species is a zeroed component, then the phase can
465  * be popped if the component increases in mole number as the phase moles
466  * are increased.
467  *
468  * @param iphasePop id of the phase, which is currently zeroed,
469  *
470  * @return Returns true if the phase can come into existence
471  * and false otherwise.
472  */
473  bool vcs_popPhasePossible(const size_t iphasePop) const;
474
475  //! Determine the list of problems that need to be checked to see if there are any phases pops
476  /*!
477  * This routine evaluates and fills in #phasePopProblemLists_. Need to
478  * work in species that are zeroed by element constraints.
479  *
480  * @return Returns the number of problems that must be checked.
481  */
483
484  //! Decision as to whether a phase pops back into existence
485  /*!
486  * @param phasePopPhaseIDs Vector containing the phase ids of the phases
487  * that will be popped this step.
488  *
489  * @return returns the phase id of the phase that pops back into
490  * existence. Returns -1 if there are no phases
491  */
492  size_t vcs_popPhaseID(std::vector<size_t> &phasePopPhaseIDs);
493
494  //! Calculates the deltas of the reactions due to phases popping
495  //! into existence
496  /*!
498  * refers to the irxn'th species formation reaction. This adjustment is
499  * for species irxn + M, where M is the number of components.
500  *
501  * @param iphasePop Phase id of the phase that will come into existence
502  *
503  * @return Returns an int representing the status of the step
504  * - 0 : normal return
505  * - 1 : A single species phase species has been zeroed out
506  * in this routine. The species is a noncomponent
507  * - 2 : Same as one but, the zeroed species is a component.
508  * - 3 : Nothing was done because the phase couldn't be birthed
509  * because a needed component is zero.
510  */
511  int vcs_popPhaseRxnStepSizes(const size_t iphasePop);
512
513  //! Calculates formation reaction step sizes.
514  /*!
515  * This is equation 6.4-16, p. 143 in Smith and Missen.
516  *
517  * Output
518  * -------
519  * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers to
520  * the irxn'th species formation reaction.
521  * This adjustment is for species irxn + M,
522  * where M is the number of components.
523  *
524  * Special branching occurs sometimes. This causes the component basis
525  * to be reevaluated
526  *
527  * @param forceComponentCalc integer flagging whether a component
528  * recalculation needs to be carried out.
529  * @param kSpecial species number of phase being zeroed.
530  *
531  * @return Returns an int representing which phase may need to be zeroed
532  */
533  size_t vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial);
534
535  //! Calculates the total number of moles of species in all phases.
536  /*!
537  * Also updates the variable m_totalMolNum and Reconciles Phase existence
538  * flags with total moles in each phase.
539  */
540  double vcs_tmoles();
541 #ifdef DEBUG_MODE
542  void check_tmoles() const;
543 #endif
544
545  //! This subroutine calculates reaction free energy changes for
546  //! all noncomponent formation reactions.
547  /*!
548  * Formation reactions are
549  * reactions which create each noncomponent species from the component
550  * species. m_stoichCoeffRxnMatrix[irxn][jcomp] are the stoichiometric
551  * coefficients for these reactions. A stoichiometric coefficient of
552  * one is assumed for species irxn in this reaction.
553  *
554  * @param l
555  * - L < 0: Calculate reactions corresponding to major noncomponent
556  * and zeroed species only
557  * - L = 0: Do all noncomponent reactions, i, between
558  * 0 <= i < irxnl
559  * - L > 0: Calculate reactions corresponding to minor noncomponent
560  * and zeroed species only
561  *
562  * @param doDeleted Do deleted species
563  * @param vcsState Calculate deltaG corresponding to either old or new
564  * free energies
565  * @param alterZeroedPhases boolean indicating whether we should
566  * add in a special section for zeroed phases.
567  *
568  * Note we special case one important issue. If the component has zero
569  * moles, then we do not allow deltaG < 0.0 for formation reactions which
570  * would lead to the loss of more of that same component. This dG < 0.0
571  * condition feeds back into the algorithm in several places, and leads
572  * to a infinite loop in at least one case.
573  */
574  void vcs_deltag(const int l, const bool doDeleted, const int vcsState,
575  const bool alterZeroedPhases = true);
576
577  void vcs_printDeltaG(const int stateCalc);
578
579  //! Calculate deltag of formation for all species in a single phase.
580  /*!
581  * Calculate deltag of formation for all species in a single phase. It is
582  * assumed that the fe[] is up to date for all species. However, if the
583  * phase is currently zeroed out, a subproblem is calculated to solve for
584  * AC[i] and pseudo-X[i] for that phase.
585  *
586  * @param iphase phase index of the phase to be calculated
587  * @param doDeleted boolean indicating whether to do deleted
588  * species or not
589  * @param stateCalc integer describing which set of free energies
590  * to use and where to stick the results.
591  * @param alterZeroedPhases boolean indicating whether we should
592  * add in a special section for zeroed phases.
593  *
594  * NOTE: this is currently not used used anywhere.
595  * It may be in the future?
596  */
597  void vcs_deltag_Phase(const size_t iphase, const bool doDeleted,
598  const int stateCalc, const bool alterZeroedPhases = true);
599
600  //! Swaps the indices for all of the global data for two species, k1
601  //! and k2.
602  /*!
603  * @param ifunc: If true, switch the species data and the noncomponent
604  * reaction data. This must be called for a non-component
605  * species only. If false, switch the species data only.
606  * Typically, we use this option when determining the
607  * component species and at the end of the calculation,
608  * when we want to return unscrambled results. All rxn
609  * data will be out-of-date.
610  *
611  * @param k1 First species index
612  * @param k2 Second species index
613  */
614  void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2);
615
616  //! Birth guess returns the number of moles of a species
617  //! that is coming back to life.
618  /*!
619  * Birth guess returns the number of moles of a species that is coming
620  * back to life. Note, this routine is not applicable if the whole phase
621  * is coming back to life, not just one species in that phase.
622  *
623  * Do a minor alt calculation. But, cap the mole numbers at 1.0E-15. For
624  * SS phases use VCS_DELETE_SPECIES_CUTOFF * 100.
625  *
626  * The routine makes sure the guess doesn't reduce the concentration of a
627  * component by more than 1/3. Note this may mean that the vlaue coming
628  * back from this routine is zero or a very small number.
629  *
630  * @param kspec Species number that is coming back to life
631  * @return Returns the number of kmol that the species should have.
632  */
633  double vcs_birthGuess(const int kspec);
634
635  //! Routine that independently determines whether a phase should be popped
636  //! under the current conditions.
637  /*
638  * This is the main routine that solves for equilibrium at constant T and
639  * P using a variant of the VCS method. Nonideal phases can be
640  * accommodated as well. Any number of single-species phases and multi-
641  * species phases can be handled by the present version.
642  *
643  * @param print_lvl 1 -> Print results to standard output; -> don't
644  * report on anything
645  * @param printDetails 1 -> Print intermediate results.
646  * @param maxit Maximum number of iterations for the algorithm
647  *
648  * @return
649  * - 0 = Equilibrium Achieved
650  * - 1 = Range space error encountered. The element abundance criteria are
651  * only partially satisfied. Specifically, the first NC= (number of
652  * components) conditions are satisfied. However, the full NE (number of
653  * elements) conditions are not satisfied. The equilibrium condition is
654  * returned.
655  * - -1 = Maximum number of iterations is exceeded. Convergence was not
656  * found.
657  */
658  int vcs_solve_phaseStability(const int iphase, int ifunc, double& funcval, int print_lvl);
659
660  //! Main program to test whether a deleted phase should be brought
661  //! back into existence
662  /*!
663  * @param iph Phase id of the deleted phase
664  */
665  double vcs_phaseStabilityTest(const size_t iph);
666
667  //! Solve an equilibrium problem at a particular fixed temperature
668  //! and pressure
669  /*!
670  * The actual problem statement is assumed to be in the structure
671  * already. This is a wrapper around the solve_TP() function. In this
672  * wrapper, we nondimensionalize the system we calculate the standard
673  * state gibbs free energies of the species, and we decide whether to we
674  * need to use the initial guess algorithm.
675  *
676  * @param ipr = 1 -> Print results to standard output;
677  * 0 -> don't report on anything
678  * @param ip1 = 1 -> Print intermediate results;
679  * 0 -> Dont print any intermediate results
680  * @param maxit Maximum number of iterations for the algorithm
681  * @param T Value of the Temperature (Kelvin)
682  * @param pres Value of the Pressure (units given by m_VCS_UnitsFormat variable
683  *
684  * @return Returns an integer representing the success of the algorithm
685  * * 0 = Equilibrium Achieved
686  * * 1 = Range space error encountered. The element abundance criteria are
687  * only partially satisfied. Specifically, the first NC= (number of
688  * components) conditions are satisfied. However, the full NE (number of
689  * elements) conditions are not satisfied. The equilibrium condition is
690  * returned.
691  * * -1 = Maximum number of iterations is exceeded. Convergence was not
692  * found.
693  */
694  int vcs_TP(int ipr, int ip1, int maxit, double T, double pres);
695
696  /*!
697  * Evaluate the standard state free energies at the current temperature
698  * and pressure. Ideal gas pressure contribution is added in here.
699  *
700  * @param ipr 1 -> Print results to standard output; 0 -> don't report
701  * on anything
702  * @param ip1 1 -> Print intermediate results; 0 -> don't.
703  * @param Temp Temperature (Kelvin)
704  * @param pres Pressure (Pascal)
705  */
706  int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres);
707
708  //! Initialize the chemical potential of single species phases
709  /*!
710  * For single species phases, initialize the chemical potential with the
711  * value of the standard state chemical potential. This value doesn't
712  * change during the calculation
713  */
714  void vcs_fePrep_TP();
715
716  //! Calculation of the total volume and the partial molar volumes
717  /*!
718  * This function calculates the partial molar volume for all species,
719  * kspec, in the thermo problem at the temperature TKelvin and pressure,
720  * Pres, pres is in atm. And, it calculates the total volume of the
721  * combined system.
722  *
723  * @param[in] tkelvin Temperature in kelvin()
724  * @param[in] pres Pressure in Pascal
725  * @param[in] w w[] is the vector containing the current mole
726  * numbers in units of kmol.
727  * @param[out] volPM[] For species in all phase, the entries are the
728  * partial molar volumes units of M**3 / kmol.
729  * @return The return value is the total volume of
730  * the entire system in units of m**3.
731  */
732  double vcs_VolTotal(const double tkelvin, const double pres,
733  const double w[], double volPM[]);
734
735  //! This routine is mostly concerned with changing the private data
736  //! to be consistent with what's needed for solution. It is called one
737  //! time for each new problem structure definition.
738  /*!
739  * This routine is always followed by vcs_prep(). Therefore, tasks
740  * that need to be done for every call to vcsc() should be placed in
741  * vcs_prep() and not in this routine.
742  *
743  * The problem structure refers to:
744  *
745  * - the number and identity of the species.
746  * - the formula matrix and thus the number of components.
747  * - the number and identity of the phases.
748  * - the equation of state
749  * - the method and parameters for determining the standard state
750  * - The method and parameters for determining the activity coefficients.
751  *
753  * 1. Fill in the SSPhase[] array.
754  * 2. Check to see if any multispecies phases actually have only one
755  * species in that phase. If true, reassign that phase and species
756  * to be a single-species phase.
757  * 3. Determine the number of components in the problem if not already
758  * done so. During this process the order of the species is changed
759  * in the private data structure. All references to the species
760  * properties must employ the ind[] index vector.
761  *
762  * @param printLvl Print level of the routine
763  * @return VCS_SUCCESS = everything went OK
764  */
765  int vcs_prep_oneTime(int printLvl);
766
767  //! Prepare the object for solution
768  /*!
769  * This routine is mostly concerned with changing the private data
770  * to be consistent with that needed for solution. It is called for
771  * every invocation of the vcs_solve() except for the cleanup invocation.
772  *
774  * 1. Initialization of arrays to zero.
775  *
776  * @return
777  * VCS_SUCCESS = everything went OK;
778  * VCS_PUB_BAD = There is an irreconcilable difference in the
779  * public data structure from when the problem was
780  * initially set up.
781  */
782  int vcs_prep();
783
784  //! In this routine, we check for things that will cause the algorithm
785  //! to fail.
786  /*!
787  * We check to see if the problem is well posed. If it is not, we return
788  * false and print out error conditions.
789  *
790  * Current there is one condition. If all the element abundances are
791  * zero, the algorithm will fail.
792  *
793  * @param vprob VCS_PROB pointer to the definition of the equilibrium
794  * problem
795  *
796  * @return If true, the problem is well-posed. If false, the problem
797  * is not well posed.
798  */
799  bool vcs_wellPosed(VCS_PROB* vprob);
800
801  //! Rearrange the constraint equations represented by the Formula
802  //! Matrix so that the operational ones are in the front
803  /*!
804  * This subroutine handles the rearrangement of the constraint equations
805  * represented by the Formula Matrix. Rearrangement is only necessary when
806  * the number of components is less than the number of elements. For this
807  * case, some constraints can never be satisfied exactly, because the
808  * range space represented by the Formula Matrix of the components can't
809  * span the extra space. These constraints, which are out of the range
810  * space of the component Formula matrix entries, are migrated to the back
811  * of the Formula matrix.
812  *
813  * A prototypical example is an extra element column in FormulaMatrix[],
814  * which is identically zero. For example, let's say that argon is has an
815  * element column in FormulaMatrix[], but no species in the mechanism
816  * actually contains argon. Then, nc < ne. Also, without perturbation of
817  * FormulaMatrix[] vcs_basopt[] would produce a zero pivot because the
818  * matrix would be singular (unless the argon element column was already
819  * the last column of FormulaMatrix[].
820  *
821  * This routine borrows heavily from vcs_basopt's algorithm. It finds nc
822  * constraints which span the range space of the Component Formula matrix,
823  * and assigns them as the first nc components in the formula matrix. This
824  * guarantees that vcs_basopt[] has a nonsingular matrix to invert.
825  *
826  * Other Variables
827  * @param aw Mole fraction work space (ne in length)
828  * @param sa Gram-Schmidt orthog work space (ne in length)
829  * @param sm QR matrix work space (ne*ne in length)
830  * @param ss Gram-Schmidt orthog work space (ne in length)
831  */
832  int vcs_elem_rearrange(double* const aw, double* const sa,
833  double* const sm, double* const ss);
834
835  //! Swaps the indices for all of the global data for two elements, ipos
836  //! and jpos.
837  /*!
838  * This function knows all of the element information with VCS_SOLVE, and
839  * can therefore switch element positions
840  *
841  * @param ipos first global element index
842  * @param jpos second global element index
843  */
844  void vcs_switch_elem_pos(size_t ipos, size_t jpos);
845
846  //! Calculates reaction adjustments using a full Hessian approximation
847  /*!
848  * This does what equation 6.4-16, p. 143 in Smith and Missen is supposed
849  * to do. However, a full matrix is formed and then solved via a conjugate
850  * gradient algorithm. No preconditioning is done.
851  *
852  * If special branching is warranted, then the program bails out.
853  *
854  * Output
855  * -------
856  * DS(I) : reaction adjustment, where I refers to the Ith species
857  * Special branching occurs sometimes. This causes the component basis
858  * to be reevaluated
859  * return = 0 : normal return
860  * 1 : A single species phase species has been zeroed out
861  * in this routine. The species is a noncomponent
862  * 2 : Same as one but, the zeroed species is a component.
863  *
864  * Special attention is taken to flag cases where the direction of the
865  * update is contrary to the steepest descent rule. This is an important
866  * attribute of the regular vcs algorithm. We don't want to violate this.
867  *
868  * NOTE: currently this routine is not used.
869  */
871
872  //! Calculates the diagonal contribution to the Hessian due to
873  //! the dependence of the activity coefficients on the mole numbers.
874  /*!
875  * (See framemaker notes, Eqn. 20 - VCS Equations document)
876  *
877  * We allow the diagonal to be increased positively to any degree.
878  * We allow the diagonal to be decreased to 1/3 of the ideal solution
879  * value, but no more -> it must remain positive.
880  *
881  * NOTE: currently this routine is not used
882  */
883  double vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal);
884
885  //! Calculates the diagonal contribution to the Hessian due to
886  //! the dependence of the activity coefficients on the mole numbers.
887  /*!
888  * (See framemaker notes, Eqn. 20 - VCS Equations document)
889  *
890  * NOTE: currently this routine is not used
891  */
892  double vcs_Hessian_actCoeff_diag(size_t irxn);
893
894  //! Recalculate all of the activity coefficients in all of the phases
895  //! based on input mole numbers
896  /*
897  * @param moleSpeciesVCS kmol of species to be used in the update.
898  *
899  * NOTE: This routine needs to be regulated.
900  */
901  void vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS);
902
903 #ifdef DEBUG_MODE
904  //! A line search algorithm is carried out on one reaction
905  /*!
906  * In this routine we carry out a rough line search algorithm to make
907  * sure that the m_deltaGRxn_new doesn't switch signs prematurely.
908  *
909  * @param irxn Reaction number
910  * @param dx_orig Original step length
911  *
912  * @param ANOTE Output character string stating the conclusions of the
913  * line search
914  * @return Returns the optimized step length found by the search
915  */
916  double vcs_line_search(const size_t irxn, const double dx_orig,
917  char* const ANOTE);
918 #else
919  double vcs_line_search(const size_t irxn, const double dx_orig);
920 #endif
921
922  //! Print out a report on the state of the equilibrium problem to
923  //! standard output.
924  /*!
925  * @param iconv Indicator of convergence, to be printed out in the report:
926  * - 0 converged
927  * - 1 range space error
928  * - -1 not converged
929  */
930  int vcs_report(int iconv);
931
932  //! Switch all species data back to the original order.
933  /*!
934  * This destroys the data based on reaction ordering.
935  */
936  int vcs_rearrange();
937
938  //! Returns the multiplier for electric charge terms
939  /*
940  * This is basically equal to F/RT
941  *
942  * @param mu_units integer representing the dimensional units system
943  * @param TKelvin double Temperature in Kelvin
944  *
945  * @return Returns the value of F/RT
946  */
947  double vcs_nondim_Farad(int mu_units, double TKelvin) const;
948
949  //! Returns the multiplier for the nondimensionalization of the equations
950  /*!
951  * This is basically equal to RT
952  *
953  * @param mu_units integer representing the dimensional units system
954  * @param TKelvin double Temperature in Kelvin
955  *
956  * @return Returns the value of RT
957  */
958  double vcs_nondimMult_TP(int mu_units, double TKelvin) const;
959
960  //! Nondimensionalize the problem data
961  /*!
962  * Nondimensionalize the free energies using the divisor, R * T
963  *
964  * Essentially the internal data can either be in dimensional form
965  * or in nondimensional form. This routine switches the data from
966  * dimensional form into nondimensional form.
967  *
968  * @todo Add a scale factor based on the total mole numbers.
969  * The algorithm contains hard coded numbers based on the
970  * total mole number. If we ever were faced with a problem
971  * with significantly different total kmol numbers than one
972  * the algorithm would have problems.
973  */
974  void vcs_nondim_TP();
975
976  //! Redimensionalize the problem data
977  /*!
978  * Reddimensionalize the free energies using the multiplier R * T
979  *
980  * Essentially the internal data can either be in dimensional form
981  * or in nondimensional form. This routine switches the data from
982  * nondimensional form into dimensional form.
983  */
984  void vcs_redim_TP();
985
986  //! Print the string representing the Chemical potential units
987  /*!
988  * This gets printed using plogf()
989  *
990  * @param unitsFormat Integer representing the units system
991  */
992  void vcs_printChemPotUnits(int unitsFormat) const;
993
994  //! Computes the current elemental abundances vector
995  /*!
996  * Computes the elemental abundances vector, m_elemAbundances[], and stores it
997  * back into the global structure
998  */
999  void vcs_elab();
1000
1001  /*!
1002  * Checks to see if the element abundances are in compliance. If they are,
1003  * then TRUE is returned. If not, FALSE is returned. Note the number of
1004  * constraints checked is usually equal to the number of components in the
1005  * problem. This routine can check satisfaction of all of the constraints
1006  * in the problem, which is equal to ne. However, the solver can't fix
1007  * breakage of constraints above nc, because that nc is the range space by
1008  * definition. Satisfaction of extra constraints would have had to occur
1009  * in the problem specification.
1010  *
1011  * The constraints should be broken up into 2 sections. If a constraint
1012  * involves a formula matrix with positive and negative signs, and eaSet =
1013  * 0.0, then you can't expect that the sum will be zero. There may be
1014  * roundoff that inhibits this. However, if the formula matrix is all of
1015  * one sign, then this requires that all species with nonzero entries in
1016  * the formula matrix be identically zero. We put this into the logic
1017  * below.
1018  *
1019  * @param ibound 1: Checks constraints up to the number of elements;
1020  * 0: Checks constraints up to the number of components.
1021  */
1022  bool vcs_elabcheck(int ibound);
1023
1024  /*!
1025  * Computes the elemental abundances vector for a single phase,
1026  * elemAbundPhase[], and returns it through the argument list. The mole
1027  * numbers of species are taken from the current value in
1028  * m_molNumSpecies_old[].
1029  */
1030  void vcs_elabPhase(size_t iphase, double* const elemAbundPhase);
1031
1032  /*!
1033  * This subroutine corrects for element abundances. At the end of the
1034  * subroutine, the total moles in all phases are recalculated again,
1035  * because we have changed the number of moles in this routine.
1036  *
1037  * @param aa temporary work vector, length ne*ne
1038  * @param x temporary work vector, length ne
1039  *
1040  * @returns
1041  * - 0 = Nothing of significance happened, Element abundances were and
1042  * still are good.
1043  * - 1 = The solution changed significantly; The element abundances are
1044  * now good.
1045  * - 2 = The solution changed significantly, The element abundances are
1047  * - 3 = The solution changed significantly, The element abundances are
1048  * still bad and a component species got zeroed out.
1049  *
1050  * Internal data to be worked on::
1051  * - ga Current element abundances
1052  * - m_elemAbundancesGoal Required elemental abundances
1053  * - m_molNumSpecies_old Current mole number of species.
1054  * - m_formulaMatrix[][] Formula matrix of the species
1055  * - ne Number of elements
1056  * - nc Number of components.
1057  *
1058  * NOTES:
1059  * This routine is turning out to be very problematic. There are
1060  * lots of special cases and problems with zeroing out species.
1061  *
1062  * Still need to check out when we do loops over nc vs. ne.
1063  *
1064  */
1065  int vcs_elcorr(double aa[], double x[]);
1066
1067  //! Create an initial estimate of the solution to the thermodynamic
1068  //! equilibrium problem.
1069  /*!
1070  * @return Return value indicates success:
1071  * - 0: successful initial guess
1072  * - -1: Unsuccessful initial guess; the elemental abundances aren't
1073  * satisfied.
1074  */
1075  int vcs_inest_TP();
1076
1077  //! Estimate the initial mole numbers by constrained linear programming
1078  /*!
1079  * This is done by running each reaction as far forward or backward as
1080  * possible, subject to the constraint that all mole numbers remain non-
1081  * negative. Reactions for which \f$\Delta \mu^0 \f$ are positive are run
1082  * in reverse, and ones for which it is negative are run in the forward
1083  * direction. The end result is equivalent to solving the linear
1084  * programming problem of minimizing the linear Gibbs function subject to
1085  * the element and non-negativity constraints.
1086  */
1087  int vcs_setMolesLinProg();
1088
1089  //! Calculate the total dimensionless Gibbs free energy
1090  /*!
1091  * Inert species are handled as if they had a standard free energy of
1092  * zero. Note, for this algorithm this function should be monotonically
1093  * decreasing.
1094  */
1095  double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles);
1096
1097  //! Calculate the total dimensionless Gibbs free energy of a single phase
1098  /*!
1099  * Inert species are handled as if they had a standard free energy of
1100  * zero and if they obeyed ideal solution/gas theory.
1101  *
1102  * @param iphase ID of the phase
1103  * @param w Species mole number vector for all species
1104  * @param fe vector of partial molar free energies of all of the
1105  * species
1106  */
1107  double vcs_GibbsPhase(size_t iphase, const double* const w,
1108  const double* const fe);
1109
1110  //! Transfer the results of the equilibrium calculation back to VCS_PROB
1111  /*!
1112  * The VCS_PROB structure is returned to the user.
1113  *
1114  * @param pub Pointer to VCS_PROB object that will get the results of the
1115  * equilibrium calculation transfered to it.
1116  */
1117  int vcs_prob_update(VCS_PROB* pub);
1118
1119  //! Fully specify the problem to be solved using VCS_PROB
1120  /*!
1121  * Use the contents of the VCS_PROB to specify the contents of the
1122  * private data, VCS_SOLVE.
1123  *
1124  * @param pub Pointer to VCS_PROB that will be used to
1125  * initialize the current equilibrium problem
1126  */
1127  int vcs_prob_specifyFully(const VCS_PROB* pub);
1128
1129  //! Specify the problem to be solved using VCS_PROB, incrementally
1130  /*!
1131  * Use the contents of the VCS_PROB to specify the contents of the
1132  * private data, VCS_SOLVE.
1133  *
1134  * It's assumed we are solving the same problem.
1135  *
1136  * @param pub Pointer to VCS_PROB that will be used to
1137  * initialize the current equilibrium problem
1138  */
1139  int vcs_prob_specify(const VCS_PROB* pub);
1140
1141 private:
1142
1143  //! Zero out the concentration of a species.
1144  /*!
1145  * Make sure to conserveelements and keep track of the total moles in all
1146  * phases.
1147  * - w[]
1148  * - m_tPhaseMoles_old[]
1149  *
1150  * @param kspec Species index
1151  *
1152  * @return:
1153  * 1: succeeded
1154  * 0: failed.
1155  */
1156  int vcs_zero_species(const size_t kspec);
1157
1158  //! Change a single species from active to inactive status
1159  /*!
1160  * Rearrange data when species is added or removed. The Lth species is
1161  * moved to the back of the species vector. The back of the species
1162  * vector is indicated by the value of MR, the current number of
1163  * active species in the mechanism.
1164  *
1165  * @param kspec Species Index
1166  * @return
1167  * Returns 0 unless.
1168  * The return is 1 when the current number of
1169  * noncomponent species is equal to zero. A recheck of deleted species
1170  * is carried out in the main code.
1171  */
1172  int vcs_delete_species(const size_t kspec);
1173
1174  //! This routine handles the bookkeeping involved with the
1175  //! deletion of multiphase phases from the problem.
1176  /*!
1177  * When they are deleted, all of their species become active
1178  * species, even though their mole numbers are set to zero.
1179  * The routine does not make the decision to eliminate multiphases.
1180  *
1181  * Note, species in phases with zero mole numbers are still
1182  * considered active. Whether the phase pops back into
1183  * existence or not is checked as part of the main iteration
1184  * loop.
1185  *
1186  * @param iph Phase to be deleted
1187  *
1188  * @return Returns whether the operation was successful or not
1189  */
1190  bool vcs_delete_multiphase(const size_t iph);
1191
1192  //! Change the concentration of a species by delta moles.
1193  /*!
1194  * Make sure to conserve elements and keep track of the total kmoles in
1195  * all phases.
1196  *
1197  * @param kspec The species index
1198  * @param delta_ptr pointer to the delta for the species. This may
1199  * change during the calculation
1200  *
1201  * @return
1202  * 1: succeeded without change of dx
1203  * 0: Had to adjust dx, perhaps to zero, in order to do the delta.
1204  */
1205  int delta_species(const size_t kspec, double* const delta_ptr);
1206
1207  //! Provide an estimate for the deleted species in phases that
1208  //! are not zeroed out
1209  /*!
1210  * Try to add back in all deleted species. An estimate of the kmol numbers
1211  * are obtained and the species is added back into the equation system,
1212  * into the old state vector.
1213  *
1214  * This routine is called at the end of the calculation, just before
1215  * returning to the user.
1216  */
1218
1219  //! Recheck deleted species in multispecies phases.
1220  /*!
1221  * We are checking the equation:
1222  *
1223  * sum_u = sum_j_comp [ sigma_i_j * u_j ]
1224  * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
1225  *
1226  * by first evaluating:
1227  *
1228  * DG_i_O = u_i_O - sum_u.
1229  *
1230  * Then, if TL is zero, the phase pops into existence if DG_i_O < 0.
1231  * Also, if the phase exists, then we check to see if the species
1232  * can have a mole number larger than VCS_DELETE_SPECIES_CUTOFF
1233  * (default value = 1.0E-32).
1234  *
1235  */
1236  int vcs_recheck_deleted();
1237
1238  //! Recheck deletion condition for multispecies phases.
1239  /*!
1240  * We assume here that DG_i_0 has been calculated for deleted species correctly
1241  *
1242  * m_feSpecies(I) = m_SSfeSpecies(I)
1243  * + ln(ActCoeff[I])
1244  * - ln(Mnaught * m_units)
1245  * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase];
1246  *
1247  * sum_u = sum_j_comp [ sigma_i_j * u_j ]
1248  * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old)
1249  *
1250  * DG_i_0 = m_feSpecies(I) - sum_m{ a_i_m DG_m }
1251  *
1252  * by first evaluating:
1253  *
1254  * DG_i_O = u_i_O - sum_u.
1255  *
1256  * Then, the phase pops into existence iff
1257  *
1258  * phaseDG = 1.0 - sum_i{exp(-DG_i_O)} < 0.0
1259  *
1260  * This formula works for both single species phases and for multispecies
1261  * phases. It's an overkill for single species phases.
1262  *
1263  * @param iphase Phase index number
1264  *
1265  * @return Returns true if the phase is currently deleted
1266  * but should be reinstated. Returns false otherwise.
1267  *
1268  * NOTE: this routine is currently not used in the code, and
1269  * contains some basic changes that are incompatible.
1270  *
1271  * assumptions:
1272  * 1. Vphase Existence is up to date
1273  * 2. Vphase->IndSpecies is up to date
1274  * 3. m_deltaGRxn_old[irxn] is up to date
1275  */
1276  bool recheck_deleted_phase(const int iphase);
1277
1278  //! Minor species alternative calculation
1279  /*!
1280  * This is based upon the following approximation: The mole fraction
1281  * changes due to these reactions don't affect the mole numbers of the
1282  * component species. Therefore the following approximation is valid for
1283  * a small component of an ideal phase:
1284  *
1285  * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I))
1286  *
1287  * m_deltaGRxn_old contains the contribution from
1288  *
1289  * m_feSpecies_old(I) =
1290  * m_SSfeSpecies(I) +
1291  * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph))
1292  * Thus,
1293  *
1294  * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I))
1295  *
1296  * Most of this section is mainly restricting the update to reasonable
1297  * values. We restrict the update a factor of 1.0E10 up and 1.0E-10 down
1298  * because we run into trouble with the addition operator due to roundoff
1299  * if we go larger than ~1.0E15. Roundoff will then sometimes produce
1300  * zero mole fractions.
1301  *
1302  * Note: This routine was generalized to incorporate nonideal phases and
1303  * phases on the molality basis
1304  *
1305  * @param[in] kspec The current species and corresponding formation
1306  * reaction number.
1307  * @param[in] irxn The current species and corresponding formation
1308  * reaction number.
1309  * @param[out] do_delete: BOOLEAN which if true on return, then we
1310  * branch to the section that deletes a species
1311  * from the current set of active species.
1312  */
1313  double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete
1314 #ifdef DEBUG_MODE
1315  , char* ANOTE
1316 #endif
1317  ) const;
1318
1319  //! This routine optimizes the minimization of the total gibbs free energy
1320  //! by making sure the slope of the following functional stays negative:
1321  /*!
1322  * The slope of the following functional is equivalent to the slope
1323  * of the total Gibbs free energy of the system:
1324  *
1325  * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] )
1326  *
1327  * along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0<al<1)
1328  * such that the a parabola approximation to Gibbs(al) fit to the
1329  * end points al = 0 and al = 1 is minimized.
1330  * - s1 = slope of Gibbs function at al = 0, which is the previous
1331  * solution = d(Gibbs)/d(al).
1332  * - s2 = slope of Gibbs function at al = 1, which is the current
1333  * solution = d(Gibbs)/d(al).
1334  *
1335  * Only if there has been an inflection point (i.e., s1 < 0 and s2 > 0),
1336  * does this code section kick in. It finds the point on the parabola
1337  * where the slope is equal to zero.
1338  */
1339  bool vcs_globStepDamp();
1340
1341  //! Switch rows and columns of a square matrix
1342  /*!
1343  * Switches the row and column of a matrix. So that after
1344  *
1345  * J[k1][j] = J_old[k2][j] and J[j][k1] = J_old[j][k2]
1346  * J[k2][j] = J_old[k1][j] and J[j][k2] = J_old[j][k1]
1347  *
1348  * @param Jac Double pointer to the Jacobian
1349  * @param k1 first row/column value to be switched
1350  * @param k2 second row/column value to be switched
1351  */
1352  void vcs_switch2D(double* const* const Jac,
1353  const size_t k1, const size_t k2) const;
1354
1355  //! Calculate the norm of a deltaGibbs free energy vector
1356  /*!
1357  * Positive DG for species which don't exist are ignored.
1358  *
1359  * @param dg Vector of local delta G's.
1360  */
1361  double l2normdg(double dg[]) const;
1362
1363 #ifdef DEBUG_MODE
1364
1365  //! Print out and check the elemental abundance vector
1366  void prneav() const;
1367
1368  void checkDelta1(double* const ds, double* const delTPhMoles, int kspec);
1369 #endif
1370
1371  //! Estimate equilibrium compositions
1372  /*!
1373  * Algorithm covered in a section of Smith and Missen's Book.
1374  *
1375  * Linear programming module is based on using dbolm.
1376  *
1377  * @param aw aw[i[ Mole fraction work space (ne in length)
1378  * @param sa sa[j] = Gram-Schmidt orthog work space (ne in length)
1379  * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length)
1380  * @param ss ss[j] = Gram-Schmidt orthog work space (ne in length)
1381  * @param test This is a small negative number.
1382  */
1383  void vcs_inest(double* const aw, double* const sa, double* const sm,
1384  double* const ss, double test);
1385
1386  //! Calculate the status of single species phases.
1387  void vcs_SSPhase(void);
1388
1389  //! This function recalculates the deltaG for reaction, irxn
1390  /*!
1391  * This function recalculates the deltaG for reaction irxn, given the mole
1392  * numbers in molNum. It uses the temporary space mu_i, to hold the
1393  * recalculated chemical potentials. It only recalculates the chemical
1394  * potentials for species in phases which participate in the irxn
1395  * reaction. This function is used by the vcs_line_search algorithm() and
1396  * should not be used widely due to the unknown state it leaves the
1397  * system.
1398  *
1399  * @param[in] irxn Reaction number
1400  * @param[in] molNum Current mole numbers of species to be used as input
1401  * to the calculation (units = kmol), (length = totalNuMSpecies)
1402  *
1403  * @param[out] ac Activity coefficients (length = totalNumSpecies) Note
1404  * this is only partially formed. Only species in phases that
1405  * participate in the reaction will be updated
1406  * @param[out] mu_i dimensionless chemical potentials (length
1407  * totalNumSpecies) Note this is only partially formed. Only
1408  * species in phases that participate in the reaction will be
1409  * updated
1410  *
1411  * @return Returns the dimensionless deltaG of the reaction
1412  */
1413  double deltaG_Recalc_Rxn(const int stateCalc,
1414  const size_t irxn, const double* const molNum,
1415  double* const ac, double* const mu_i);
1416
1417  //! Delete memory that isn't just resizable STL containers
1418  /*!
1419  * This gets called by the destructor or by InitSizes().
1420  */
1421  void vcs_delete_memory();
1422
1423  //! Initialize the internal counters
1424  /*!
1425  * Initialize the internal counters containing the subroutine call
1426  * values and times spent in the subroutines.
1427  *
1428  * ifunc = 0 Initialize only those counters appropriate for the top of
1429  * vcs_solve_TP().
1430  * = 1 Initialize all counters.
1431  */
1432  void vcs_counters_init(int ifunc);
1433
1434  //! Create a report on the plog file containing timing and its information
1435  /*!
1436  * @param timing_print_lvl If 0, just report the iteration count.
1437  * If larger than zero, report the timing information
1438  */
1439  void vcs_TCounters_report(int timing_print_lvl = 1);
1440
1441  void vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc);
1442
1443  void vcs_setFlagsVolPhase(const size_t iph, const bool upToDate, const int stateCalc);
1444
1445  //! Update all underlying vcs_VolPhase objects
1446  /*!
1447  * Update the mole numbers and the phase voltages of all phases in the
1448  * vcs problem
1449  *
1450  * @param stateCalc Location of the update (either VCS_STATECALC_NEW or
1451  * VCS_STATECALC_OLD).
1452  */
1453  void vcs_updateMolNumVolPhases(const int stateCalc);
1454
1455 public:
1456  //! Calculate the rank of a matrix and return the rows and columns that
1457  //! will generate an independent basis for that rank
1458  /*
1459  * Choose the optimum component species basis for the calculations,
1460  * finding the rank and set of linearly independent rows for that
1461  * calculation. Then find the set of linearly indepedent element columns
1462  * that can support that rank. This is done by taking the transpose of the
1463  * matrix and redoing the same calculation. (there may be a better way to
1464  * do this. I don't know.)
1465  *
1466  * @param[in] awtmp Vector of mole numbers which will be used to
1467  * construct a ranking for how to pick the basis species. This is
1468  * largely ignored here.
1469  * @param[in] numSpecies Number of species. This is the number of rows in
1470  * the matrix.
1471  * @param[in] matrix This is the formula matrix. Nominally, the rows are
1472  * species, while the columns are element compositions. However,
1473  * this routine is totally general, so that the rows and columns
1474  * can be anything.
1475  * @param[in] numElemConstraints Number of element constraints
1476  *
1477  * @param[out] usedZeroedSpecies If true, then a species with a zero
1478  * concentration was used as a component.
1479  * @param[out] compRes Vector of rows which are linearly independent.
1480  * (these are the components)
1481  * @param[out] elemComp Vector of columns which are linearly independent
1482  * (These are the actionable element constraints).
1483  *
1484  * @return Returns number of components. This is the rank of the matrix
1485  * @deprecated
1486  */
1487  int vcs_rank(const double* awtmp, size_t numSpecies, const double* matrix, size_t numElemConstraints,
1488  std::vector<size_t> &compRes, std::vector<size_t> &elemComp, int* const usedZeroedSpecies) const;
1489
1490  //! value of the number of species used to malloc data structures
1491  size_t NSPECIES0;
1492
1493  //! value of the number of phases used to malloc data structures
1494  size_t NPHASE0;
1495
1496  //! Total number of species in the problems
1498
1499  //! Number of element constraints in the problem
1500  /*!
1501  * This is typically equal to the number of elements in the problem
1502  */
1504
1505  //! Number of components calculated for the problem
1507
1508  //! Total number of non-component species in the problem
1509  size_t m_numRxnTot;
1510
1511  //! Current number of species in the problems
1512  /*!
1513  * Species can be deleted if they aren't
1514  * stable under the current conditions
1515  */
1517
1518  //! Current number of non-component species in the problem
1519  /*!
1520  * Species can be deleted if they aren't
1521  * stable under the current conditions
1522  */
1523  size_t m_numRxnRdc;
1524
1525  //! Number of active species which are currently either treated as
1526  //! minor species
1528
1529  //! Number of Phases in the problem
1530  size_t m_numPhases;
1531
1532  //! Formula matrix for the problem
1533  /*!
1534  * FormulaMatrix[j][kspec] = Number of elements, j, in the kspec species
1535  *
1536  * Both element and species indices are swapped.
1537  */
1539
1540  //! Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form.
1541  /*!
1542  * This is the stoichiometric coefficient matrix for the
1543  * reaction which forms species kspec from the component species. A
1544  * stoichiometric coefficient of one is assumed for the species kspec in this mechanism.
1545  *
1546  * NOTE: kspec = irxn + m_numComponents
1547  *
1548  * m_stoichCoeffRxnMatrix[irxn][j] :
1549  * j refers to the component number, and irxn refers to the irxn_th non-component species.
1550  * The stoichiometric coefficients multilplied by the Formula coefficients of the
1551  * component species add up to the negative value of the number of elements in
1552  * the species kspec.
1553  *
1554  * length = [nspecies0][nelements0]
1555  */
1557
1558  //! Absolute size of the stoichiometric coefficients
1559  /*!
1560  * scSize[irxn] = abs(Size) of the stoichiometric
1561  * coefficients. These are used to determine
1562  * whether a given species should be
1563  * handled by the alt_min treatment or
1564  * should be handled as a major species.
1565  */
1566  std::vector<double> m_scSize;
1567
1568  //! total size of the species
1569  /*!
1570  * This is used as a multiplier to the mole number in figuring out which
1571  * species should be components.
1572  */
1573  std::vector<double> m_spSize;
1574
1575  //! Standard state chemical potentials for species K at the current
1576  //! temperature and pressure.
1577  /*!
1578  * The first NC entries are for components. The following NR entries are
1579  * for the current non-component species in the mechanism.
1580  */
1581  std::vector<double> m_SSfeSpecies;
1582
1583  //! Free energy vector from the start of the current iteration
1584  /*!
1585  * The free energies are saved at the start of the current iteration.
1586  * Length = number of species
1587  */
1588  std::vector<double> m_feSpecies_old;
1589
1590  //! Dimensionless new free energy for all the species in the mechanism
1591  //! at the new tentatite T, P, and mole numbers.
1592  /*!
1593  * The first NC entries are for components. The following
1594  * NR entries are for the current non-component species in the mechanism.
1595  * Length = number of species
1596  */
1597  std::vector<double> m_feSpecies_new;
1598
1599  //! Setting for whether to do an initial estimate
1600  /*!
1601  * Initial estimate:
1602  * * 0 Do not estimate the solution at all. Use the
1603  * * supplied mole numbers as is.
1604  * * 1 Only do an estimate if the element abundances
1605  * * aren't satisfied.
1606  * * -1 Force an estimate of the soln. Throw out the input
1607  * * mole numbers.
1608  */
1610
1611  //! Total moles of the species
1612  /*!
1613  * Total number of moles of the kth species.
1614  * Length = Total number of species = m
1615  */
1616  std::vector<double> m_molNumSpecies_old;
1617
1618  //! Specifies the species unknown type
1619  /*!
1620  * There are two types. One is the straightforward species, with the mole
1621  * number w[k], as the unknown. The second is the an interfacial voltage
1622  * where w[k] refers to the interfacial voltage in volts.
1623  *
1624  * These species types correspond to metallic electrons corresponding to
1625  * electrodes. The voltage and other interfacial conditions sets up an
1626  * interfacial current, which is set to zero in this initial treatment.
1627  * Later we may have non-zero interfacial currents.
1628  */
1629  std::vector<int> m_speciesUnknownType;
1630
1631  //! Change in the number of moles of phase, iphase, due to the
1632  //! noncomponent formation reaction, irxn, for species, k:
1633  /*!
1634  * m_deltaMolNumPhase[irxn][iphase] = k = nc + irxn
1635  */
1637
1638  //! This is 1 if the phase, iphase, participates in the formation reaction
1639  //! irxn, and zero otherwise. PhaseParticipation[irxn][iphase]
1641
1642  //! electric potential of the iph phase
1643  std::vector<double> m_phasePhi;
1644
1645  //! Tentative value of the mole number vector. It's also used to store the
1646  //! mole fraction vector.
1647  //std::vector<double> wt;
1648  std::vector<double> m_molNumSpecies_new;
1649
1650  //! Delta G(irxn) for the noncomponent species in the mechanism.
1651  /*!
1652  * Computed by the subroutine deltaG. m_deltaGRxn is the free
1653  * energy change for the reaction which forms species K from the
1654  * component species. This vector has length equal to the number
1655  * of noncomponent species in the mechanism. It starts with
1656  * the first current noncomponent species in the mechanism.
1657  */
1658  std::vector<double> m_deltaGRxn_new;
1659
1660  //! Last deltag[irxn] from the previous step
1661  std::vector<double> m_deltaGRxn_old;
1662
1663  //! Last deltag[irxn] from the previous step with additions for
1664  //! possible births of zeroed phases.
1665  std::vector<double> m_deltaGRxn_Deficient;
1666
1667  //! Temporary vector of Rxn DeltaG's
1668  /*!
1669  * This is used from time to time, for printing purposes
1670  */
1671  std::vector<double> m_deltaGRxn_tmp;
1672
1673  //! Reaction Adjustments for each species during the current step
1674  /*!
1675  * delta Moles for each species during the current step.
1676  * Length = number of species
1677  */
1678  std::vector<double> m_deltaMolNumSpecies;
1679
1680  //! Element abundances vector
1681  /*!
1682  * Vector of moles of each element actually in the solution
1683  * vector. Except for certain parts of the algorithm,
1684  * this is a constant.
1685  * Note other constraint conditions are added to this vector.
1686  * This is input from the input file and
1687  * is considered a constant from thereon.
1688  * units = kmoles
1689  */
1690  std::vector<double> m_elemAbundances;
1691
1692  //! Element abundances vector Goals
1693  /*!
1694  * Vector of moles of each element that are the goals of the simulation.
1695  * This is a constant in the problem. Note other constraint conditions
1696  * are added to this vector. This is input from the input file and is
1697  * considered a constant from thereon. units = kmoles
1698  */
1699  std::vector<double> m_elemAbundancesGoal;
1700
1701  //! Total number of kmoles in all phases
1702  /*!
1703  * This number includes the inerts.
1704  * -> Don't use this except for scaling purposes
1705  */
1707
1708  //! Total kmols of species in each phase
1709  /*!
1710  * This contains the total number of moles of species in each phase
1711  *
1712  * Length = number of phases
1713  */
1714  std::vector<double> m_tPhaseMoles_old;
1715
1716  //! total kmols of species in each phase in the tentative soln vector
1717  /*!
1718  * This contains the total number of moles of species in each phase
1719  * in the tentative solution vector
1720  *
1721  * Length = number of phases
1722  */
1723  std::vector<double> m_tPhaseMoles_new;
1724
1725  //! Temporary vector of length NPhase
1726  mutable std::vector<double> m_TmpPhase;
1727
1728  //! Temporary vector of length NPhase
1729  mutable std::vector<double> m_TmpPhase2;
1730
1731  //! Change in the total moles in each phase
1732  /*!
1733  * Length number of phases.
1734  */
1735  std::vector<double> m_deltaPhaseMoles;
1736
1737  //! Temperature (Kelvin)
1739
1740  //! Pressure (units are determined by m_VCS_UnitsFormat
1741  /*!
1742  * | Values | units |
1743  * | ------ | -----
1744  * | -1: | atm |
1745  * | 0: | atm |
1746  * | 1: | atm |
1747  * | 2: | atm |
1748  * | 3: | Pa |
1749  *
1750  * Units being changed to Pa
1751  */
1753
1754  //! Total kmoles of inert to add to each phase
1755  /*!
1756  * TPhInertMoles[iph] = Total kmoles of inert to add to each phase
1757  * length = number of phases
1758  */
1759  std::vector<double> TPhInertMoles;
1760
1761  //! Tolerance requirement for major species
1762  double m_tolmaj;
1763
1764  //! Tolerance requirements for minor species
1765  double m_tolmin;
1766
1767  //! Below this, major species aren't refined any more
1768  double m_tolmaj2;
1769
1770  //! Below this, minor species aren't refined any more
1771  double m_tolmin2;
1772
1773  //! Index vector that keeps track of the species vector rearrangement
1774  /*!
1775  * At the end of each run, the species vector and associated data gets put back
1776  * in the original order.
1777  *
1778  * Example
1779  *
1780  * k = m_speciesMapIndex[kspec]
1781  *
1782  * kspec = current order in the vcs_solve object
1783  * k = original order in the vcs_prob object and in the MultiPhase object
1784  */
1785  std::vector<size_t> m_speciesMapIndex;
1786
1787  //! Index that keeps track of the index of the species within the local
1788  //! phase
1789  /*!
1790  * This returns the local index of the species within the phase. Its argument
1791  * is the global species index within the VCS problem.
1792  *
1793  * k = m_speciesLocalPhaseIndex[kspec]
1794  *
1795  * k varies between 0 and the nSpecies in the phase
1796  *
1797  * Length = number of species
1798  */
1799  std::vector<size_t> m_speciesLocalPhaseIndex;
1800
1801  //! Index vector that keeps track of the rearrangement of the elements
1802  /*!
1803  * At the end of each run, the element vector and associated data gets
1804  * put back in the original order.
1805  *
1806  * Example
1807  *
1808  * e = m_elementMapIndex[eNum]
1809  * eNum = current order in the vcs_solve object
1810  * e = original order in the vcs_prob object and in the MultiPhase object
1811  */
1812  std::vector<size_t> m_elementMapIndex;
1813
1814  //! Mapping between the species index for noncomponent species and the
1815  //! full species index.
1816  /*!
1817  * ir[irxn] = Mapping between the reaction index for noncomponent
1818  * formation reaction of a species and the full species
1819  * index.
1820  *
1821  * Initially set to a value of K = NC + I This vector has length equal to
1822  * number of noncomponent species in the mechanism. It starts with the
1823  * first current noncomponent species in the mechanism. kspec = ir[irxn]
1824  */
1825  std::vector<size_t> m_indexRxnToSpecies;
1826
1827  //! Major -Minor status vector for the species in the problem
1828  /*!
1829  * The index for this is species. The reaction that this is referring to
1830  * is kspec = irxn + m_numComponents. For possible values and their
1831  * meanings, see vcs_evaluate_speciesType().
1832  */
1833  std::vector<int> m_speciesStatus;
1834
1835  //! Mapping from the species number to the phase number
1836  std::vector<size_t> m_phaseID;
1837
1838  //! Boolean indicating whether a species belongs to a single-species phase
1839  // vector<bool> can't be used here because it doesn't work with std::swap
1840  std::vector<char> m_SSPhase;
1841
1842  //! Species string name for the kth species
1843  std::vector<std::string> m_speciesName;
1844
1845  //! Vector of strings containing the element names
1846  /*!
1847  * ElName[j] = String containing element names
1848  */
1849  std::vector<std::string> m_elementName;
1850
1851  //! Type of the element constraint
1852  /*!
1853  * * 0 - #VCS_ELEM_TYPE_ABSPOS Normal element that is positive or zero in
1854  * all species.
1855  * * 1 - #VCS_ELEM_TYPE_ELECTRONCHARGE element dof that corresponds to
1856  * the electronic charge DOF.
1857  * * 2 - #VCS_ELEM_TYPE_CHARGENEUTRALITY element dof that corresponds to
1858  * a required charge neutrality constraint on the phase. The element
1859  * abundance is always exactly zero.
1860  * * 3 - #VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may mean
1861  * that a species has neg 0 or pos value of that constraint (other than
1862  * charge)
1863  */
1864  std::vector<int> m_elType;
1865
1866  //! Specifies whether an element constraint is active
1867  /*!
1868  * The default is true. Length = nelements
1869  */
1870  std::vector<int> m_elementActive;
1871
1872  //! Array of Phase Structures. Length = number of phases.
1873  std::vector<vcs_VolPhase*> m_VolPhaseList;
1874
1875  //! String containing the title of the run
1876  std::string m_title;
1877
1878  //! This specifies the current state of units for the Gibbs free energy
1879  //! properties in the program.
1880  /*!
1881  * The default is to have this unitless
1882  */
1884
1885  //! Multiplier for the mole numbers within the nondimensionless formulation
1886  /*!
1887  * All numbers within the main routine are on an absolute basis. This
1888  * presents some problems wrt very large and very small mole numbers.
1889  * We get around this by using a multiplier coming into and coming
1890  * out of the equilibrium routines
1891  */
1893
1894  //! specifies the activity convention of the phase containing the species
1895  /*!
1896  * * 0 = molar based
1897  * * 1 = molality based
1898  *
1899  * length = number of species
1900  */
1901  std::vector<int> m_actConventionSpecies;
1902
1903  //! specifies the activity convention of the phase.
1904  /*!
1905  * * 0 = molar based
1906  * * 1 = molality based
1907  *
1908  * length = number of phases
1909  */
1910  std::vector<int> m_phaseActConvention;
1911
1912  //! specifies the ln(Mnaught) used to calculate the chemical potentials
1913  /*!
1914  * For molar based activity conventions this will be equal to 0.0.
1915  * length = number of species.
1916  */
1917  std::vector<double> m_lnMnaughtSpecies;
1918
1919  //! Molar-based Activity Coefficients for Species.
1920  //! Length = number of species
1921  std::vector<double> m_actCoeffSpecies_new;
1922
1923  //! Molar-based Activity Coefficients for Species based on old mole numbers
1924  /*!
1925  * These activity coefficients are based on the m_molNumSpecies_old
1926  * values Molar based activity coeffients. Length = number of species
1927  */
1928  std::vector<double> m_actCoeffSpecies_old;
1929
1930  //! Change in the log of the activity coefficient with respect to the mole number
1931  //! multiplied by the phase mole number
1932  /*!
1933  * length = [nspecies][nspecies]
1934  *
1935  * This is a temporary array that gets regenerated every time it's
1936  * needed. It is not swapped wrt species.
1937  */
1939
1940  //! Molecular weight of each species
1941  /*!
1942  * units = kg/kmol. length = number of species.
1943  *
1944  * note: this is a candidate for removal. I don't think we use it.
1945  */
1946  std::vector<double> m_wtSpecies;
1947
1948  //! Charge of each species. Length = number of species.
1949  std::vector<double> m_chargeSpecies;
1950
1951  std::vector<std::vector<size_t> > phasePopProblemLists_;
1952
1953  //! Vector of pointers to thermostructures which identify the model
1954  //! and parameters for evaluating the thermodynamic functions for that
1955  //! particular species.
1956  /*!
1957  * SpeciesThermo[k] pointer to the thermo information for the kth species
1958  */
1959  std::vector<VCS_SPECIES_THERMO*> m_speciesThermoList;
1960
1961  //! Choice of Hessians
1962  /*!
1963  * If this is true, then we will use a better approximation to the
1964  * Hessian based on Jacobian of the ln(ActCoeff) with respect to mole
1965  * numbers
1966  */
1968
1969  //! Total volume of all phases. Units are m^3
1970  double m_totalVol;
1971
1972  //! Partial molar volumes of the species
1973  /*!
1974  * units = mks (m^3/kmol) -determined by m_VCS_UnitsFormat
1975  * Length = number of species
1976  */
1977  std::vector<double> m_PMVolumeSpecies;
1978
1979  //! dimensionless value of Faraday's constant, F / RT (1/volt)
1981
1982  //! Timing and iteration counters for the vcs object
1984
1985  //! Debug printing lvl
1986  /*!
1987  * Levels correspond to the following guidlines
1988  * * 0 No printing at all
1989  * * 1 Serious warnings or fatal errors get one line
1990  * * 2 one line per eacdh successful vcs package call
1991  * * 3 one line per every successful solve_TP calculation
1992  * * 4 one line for every successful operation -> solve_TP gets a summary report
1993  * * 5 each iteration in solve_TP gets a report with one line per species
1994  * * 6 Each decision in solve_TP gets a line per species in addition to 4
1995  * * 10 Additionally Hessian matrix is printed out
1996  *
1997  * Levels of printing above 4 are only accessible when DEBUG_MODE is turned on
1998  */
2000
2001  //! printing level of timing information
2002  /*!
2003  * * 1 allowing printing of timing
2004  * * 0 do not allow printing of timing -> everything is printed as a NA.
2005  */
2007
2008  //! Units for the chemical potential data
2009  /*!
2010  * | Value | chemical potential units | pressure units |
2011  * | ----- | ------------------------ | -------------- |
2012  * | -1 | kcal/mol | Pa |
2013  * | 0 | MU/RT | Pa |
2014  * | 1 | kJ/mol | Pa |
2015  * | 2 | Kelvin | Pa |
2016  * | 3 | J / kmol | Pa |
2017  */
2019
2020  friend class vcs_phaseStabilitySolve;
2021 };
2022
2023 }
2024 #endif
std::vector< double > m_SSfeSpecies
Standard state chemical potentials for species K at the current temperature and pressure.
Definition: vcs_solve.h:1581
void vcs_elab()
Computes the current elemental abundances vector.
Definition: vcs_elem.cpp:14
std::vector< double > m_deltaGRxn_tmp
Temporary vector of Rxn DeltaG's.
Definition: vcs_solve.h:1671
char m_unitsState
This specifies the current state of units for the Gibbs free energy properties in the program...
Definition: vcs_solve.h:1883
size_t m_numRxnTot
Total number of non-component species in the problem.
Definition: vcs_solve.h:1509
std::vector< int > m_phaseActConvention
specifies the activity convention of the phase.
Definition: vcs_solve.h:1910
int vcs_phasePopDeterminePossibleList()
Determine the list of problems that need to be checked to see if there are any phases pops...
size_t m_numComponents
Number of components calculated for the problem.
Definition: vcs_solve.h:1506
std::vector< double > m_phasePhi
electric potential of the iph phase
Definition: vcs_solve.h:1643
std::vector< double > m_actCoeffSpecies_old
Molar-based Activity Coefficients for Species based on old mole numbers.
Definition: vcs_solve.h:1928
int vcs_prob_update(VCS_PROB *pub)
Transfer the results of the equilibrium calculation back to VCS_PROB.
Definition: vcs_solve.cpp:888
void vcs_fePrep_TP()
Initialize the chemical potential of single species phases.
Definition: vcs_TP.cpp:109
IntStarStar m_phaseParticipation
This is 1 if the phase, iphase, participates in the formation reaction irxn, and zero otherwise...
Definition: vcs_solve.h:1640
int vcs_inest_TP()
Create an initial estimate of the solution to the thermodynamic equilibrium problem.
Definition: vcs_inest.cpp:345
int vcs_report(int iconv)
Print out a report on the state of the equilibrium problem to standard output.
Definition: vcs_report.cpp:29
std::vector< VCS_SPECIES_THERMO * > m_speciesThermoList
Vector of pointers to thermostructures which identify the model and parameters for evaluating the the...
Definition: vcs_solve.h:1959
double m_tolmin
Tolerance requirements for minor species.
Definition: vcs_solve.h:1765
double l2normdg(double dg[]) const
Calculate the norm of a deltaGibbs free energy vector.
void vcs_SSPhase(void)
Calculate the status of single species phases.
Definition: vcs_prep.cpp:20
std::vector< size_t > m_indexRxnToSpecies
Mapping between the species index for noncomponent species and the full species index.
Definition: vcs_solve.h:1825
void vcs_reinsert_deleted(size_t kspec)
size_t m_numRxnMinorZeroed
Number of active species which are currently either treated as minor species.
Definition: vcs_solve.h:1527
std::vector< double > m_tPhaseMoles_new
total kmols of species in each phase in the tentative soln vector
Definition: vcs_solve.h:1723
void vcs_deltag_Phase(const size_t iphase, const bool doDeleted, const int stateCalc, const bool alterZeroedPhases=true)
Calculate deltag of formation for all species in a single phase.
void vcs_switch2D(double *const *const Jac, const size_t k1, const size_t k2) const
Switch rows and columns of a square matrix.
std::vector< int > m_actConventionSpecies
specifies the activity convention of the phase containing the species
Definition: vcs_solve.h:1901
void vcs_updateMolNumVolPhases(const int stateCalc)
Update all underlying vcs_VolPhase objects.
bool vcs_evaluate_speciesType()
This routine evaluates the species type for all species.
int m_debug_print_lvl
Debug printing lvl.
Definition: vcs_solve.h:1999
This file contains definitions of terms that are used in internal routines and are unlikely to need m...
double m_temperature
Temperature (Kelvin)
Definition: vcs_solve.h:1738
std::vector< double > m_lnMnaughtSpecies
specifies the ln(Mnaught) used to calculate the chemical potentials
Definition: vcs_solve.h:1917
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...
void vcs_initSizes(const size_t nspecies0, const size_t nelements, const size_t nphase0)
Initialize the sizes within the VCS_SOLVE object.
Definition: vcs_solve.cpp:60
int vcs_prob_specify(const VCS_PROB *pub)
Specify the problem to be solved using VCS_PROB, incrementally.
Definition: vcs_solve.cpp:763
int m_timing_print_lvl
printing level of timing information
Definition: vcs_solve.h:2006
std::vector< double > m_feSpecies_old
Free energy vector from the start of the current iteration.
Definition: vcs_solve.h:1588
int m_VCS_UnitsFormat
Units for the chemical potential data.
Definition: vcs_solve.h:2018
void vcs_inest(double *const aw, double *const sa, double *const sm, double *const ss, double test)
Estimate equilibrium compositions.
Definition: vcs_inest.cpp:22
int vcs_zero_species(const size_t kspec)
Zero out the concentration of a species.
double m_tolmin2
Below this, minor species aren't refined any more.
Definition: vcs_solve.h:1771
Provide an estimate for the deleted species in phases that are not zeroed out.
bool vcs_wellPosed(VCS_PROB *vprob)
In this routine, we check for things that will cause the algorithm to fail.
Definition: vcs_prep.cpp:238
Class to keep track of time and iterations.
Definition: vcs_internal.h:55
double deltaG_Recalc_Rxn(const int stateCalc, const size_t irxn, const double *const molNum, double *const ac, double *const mu_i)
This function recalculates the deltaG for reaction, irxn.
bool vcs_delete_multiphase(const size_t iph)
This routine handles the bookkeeping involved with the deletion of multiphase phases from the problem...
int vcs_prep()
Prepare the object for solution.
Definition: vcs_prep.cpp:219
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.
void vcs_CalcLnActCoeffJac(const double *const moleSpeciesVCS)
Recalculate all of the activity coefficients in all of the phases based on input mole numbers...
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...
std::vector< double > m_deltaGRxn_new
Delta G(irxn) for the noncomponent species in the mechanism.
Definition: vcs_solve.h:1658
void vcs_counters_init(int ifunc)
Initialize the internal counters.
Definition: vcs_solve.cpp:972
double m_totalMolNum
Total number of kmoles in all phases.
Definition: vcs_solve.h:1706
std::vector< double > m_elemAbundancesGoal
Element abundances vector Goals.
Definition: vcs_solve.h:1699
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.
int vcs_popPhaseRxnStepSizes(const size_t iphasePop)
Calculates the deltas of the reactions due to phases popping into existence.
std::string m_title
String containing the title of the run.
Definition: vcs_solve.h:1876
int m_useActCoeffJac
Choice of Hessians.
Definition: vcs_solve.h:1967
std::vector< double > m_elemAbundances
Element abundances vector.
Definition: vcs_solve.h:1690
std::vector< size_t > m_speciesMapIndex
Index vector that keeps track of the species vector rearrangement.
Definition: vcs_solve.h:1785
Defines and definitions within the vcs package.
int vcs_elcorr(double aa[], double x[])
Definition: vcs_elem.cpp:109
std::vector< double > m_wtSpecies
Molecular weight of each species.
Definition: vcs_solve.h:1946
int m_doEstimateEquil
Setting for whether to do an initial estimate.
Definition: vcs_solve.h:1609
std::vector< double > m_molNumSpecies_new
Tentative value of the mole number vector.
Definition: vcs_solve.h:1648
void vcs_printChemPotUnits(int unitsFormat) const
Print the string representing the Chemical potential units.
Definition: vcs_nondim.cpp:214
std::vector< double > m_deltaGRxn_Deficient
Last deltag[irxn] from the previous step with additions for possible births of zeroed phases...
Definition: vcs_solve.h:1665
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.
Internal declarations for the VCSnonideal package.
double vcs_birthGuess(const int kspec)
Birth guess returns the number of moles of a species that is coming back to life. ...
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:9
double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool *do_delete, char *ANOTE) const
Minor species alternative calculation.
DoubleStarStar m_formulaMatrix
Formula matrix for the problem.
Definition: vcs_solve.h:1538
std::vector< double > m_molNumSpecies_old
Total moles of the species.
Definition: vcs_solve.h:1616
std::vector< double > m_actCoeffSpecies_new
Molar-based Activity Coefficients for Species.
Definition: vcs_solve.h:1921
size_t m_numSpeciesTot
Total number of species in the problems.
Definition: vcs_solve.h:1497
DoubleStarStar m_deltaMolNumPhase
Change in the number of moles of phase, iphase, due to the noncomponent formation reaction...
Definition: vcs_solve.h:1636
dimensionless value of Faraday's constant, F / RT (1/volt)
Definition: vcs_solve.h:1980
std::vector< std::string > m_elementName
Vector of strings containing the element names.
Definition: vcs_solve.h:1849
void vcs_delete_memory()
Delete memory that isn't just resizable STL containers.
Definition: vcs_solve.cpp:223
std::vector< int > m_elementActive
Specifies whether an element constraint is active.
Definition: vcs_solve.h:1870
std::vector< double > m_scSize
Absolute size of the stoichiometric coefficients.
Definition: vcs_solve.h:1566
double m_tolmaj2
Below this, major species aren't refined any more.
Definition: vcs_solve.h:1768
int vcs_species_type(const size_t kspec) const
Evaluate the species category for the indicated species.
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:396
void prneav() const
Print out and check the elemental abundance vector.
void vcs_chemPotPhase(const int stateCalc, const size_t iph, const double *const molNum, double *const ac, double *const mu_i, const bool do_deleted=false)
We calculate the dimensionless chemical potentials of all species in a single phase.
std::vector< double > m_deltaGRxn_old
Last deltag[irxn] from the previous step.
Definition: vcs_solve.h:1661
std::vector< double > m_deltaMolNumSpecies
Reaction Adjustments for each species during the current step.
Definition: vcs_solve.h:1678
size_t m_numElemConstraints
Number of element constraints in the problem.
Definition: vcs_solve.h:1503
void vcs_printSpeciesChemPot(const int stateCalc) const
Print out a table of chemical potentials.
int vcs_solve_phaseStability(const int iphase, int ifunc, double &funcval, int print_lvl)
Routine that independently determines whether a phase should be popped under the current conditions...
std::vector< int > m_speciesStatus
Major -Minor status vector for the species in the problem.
Definition: vcs_solve.h:1833
std::vector< vcs_VolPhase * > m_VolPhaseList
Array of Phase Structures. Length = number of phases.
Definition: vcs_solve.h:1873
size_t vcs_popPhaseID(std::vector< size_t > &phasePopPhaseIDs)
Decision as to whether a phase pops back into existence.
double vcs_Total_Gibbs(double *w, double *fe, double *tPhMoles)
Calculate the total dimensionless Gibbs free energy.
Definition: vcs_Gibbs.cpp:18
DoubleStarStar 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:1938
double m_totalMoleScale
Multiplier for the mole numbers within the nondimensionless formulation.
Definition: vcs_solve.h:1892
A class for 2D int arrays stored in column-major (Fortran-compatible) form.
size_t m_numPhases
Number of Phases in the problem.
Definition: vcs_solve.h:1530
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:50
int vcs_prob_specifyFully(const VCS_PROB *pub)
Fully specify the problem to be solved using VCS_PROB.
Definition: vcs_solve.cpp:381
double vcs_nondim_Farad(int mu_units, double TKelvin) const
Returns the multiplier for electric charge terms.
Definition: vcs_nondim.cpp:19
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
bool vcs_globStepDamp()
This routine optimizes the minimization of the total gibbs free energy by making sure the slope of th...
std::vector< size_t > m_phaseID
Mapping from the species number to the phase number.
Definition: vcs_solve.h:1836
void vcs_nondim_TP()
Nondimensionalize the problem data.
Definition: vcs_nondim.cpp:76
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< double > m_feSpecies_new
Dimensionless new free energy for all the species in the mechanism at the new tentatite T...
Definition: vcs_solve.h:1597
bool recheck_deleted_phase(const int iphase)
Recheck deletion condition for multispecies phases.
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:990
int vcs_rearrange()
Switch all species data back to the original order.
std::vector< double > m_tPhaseMoles_old
Total kmols of species in each phase.
Definition: vcs_solve.h:1714
double vcs_nondimMult_TP(int mu_units, double TKelvin) const
Returns the multiplier for the nondimensionalization of the equations.
Definition: vcs_nondim.cpp:46
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...
double vcs_tmoles()
Calculates the total number of moles of species in all phases.
void vcs_elabPhase(size_t iphase, double *const elemAbundPhase)
Definition: vcs_elem.cpp:95
DoubleStarStar m_stoichCoeffRxnMatrix
Stoichiometric coefficient matrix for the reaction mechanism expressed in Reduced Canonical Form...
Definition: vcs_solve.h:1556
void vcs_redim_TP()
Redimensionalize the problem data.
Definition: vcs_nondim.cpp:167
double vcs_phaseStabilityTest(const size_t iph)
Main program to test whether a deleted phase should be brought back into existence.
std::vector< double > TPhInertMoles
Total kmoles of inert to add to each phase.
Definition: vcs_solve.h:1759
std::vector< char > m_SSPhase
Boolean indicating whether a species belongs to a single-species phase.
Definition: vcs_solve.h:1840
double vcs_line_search(const size_t irxn, const double dx_orig, char *const ANOTE)
A line search algorithm is carried out on one reaction.
double vcs_Hessian_actCoeff_diag(size_t irxn)
Calculates the diagonal contribution to the Hessian due to the dependence of the activity coefficient...
bool vcs_elabcheck(int ibound)
Definition: vcs_elem.cpp:26
int vcs_recheck_deleted()
Recheck deleted species in multispecies phases.
std::vector< double > m_PMVolumeSpecies
Partial molar volumes of the species.
Definition: vcs_solve.h:1977
Interface class for the vcs thermo equilibrium solver package, which generally describes the problem ...
Definition: vcs_prob.h:27
std::vector< size_t > m_elementMapIndex
Index vector that keeps track of the rearrangement of the elements.
Definition: vcs_solve.h:1812
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:43
int vcs_rank(const double *awtmp, size_t numSpecies, const double *matrix, size_t numElemConstraints, std::vector< size_t > &compRes, std::vector< size_t > &elemComp, int *const usedZeroedSpecies) const
Calculate the rank of a matrix and return the rows and columns that will generate an independent basi...
Definition: vcs_rank.cpp:46
std::vector< int > m_elType
Type of the element constraint.
Definition: vcs_solve.h:1864
double m_pressurePA
Pressure (units are determined by m_VCS_UnitsFormat.
Definition: vcs_solve.h:1752
std::vector< double > m_TmpPhase2
Temporary vector of length NPhase.
Definition: vcs_solve.h:1729
size_t vcs_RxnStepSizes(int &forceComponentCalc, size_t &kSpecial)
Calculates formation reaction step sizes.
size_t NSPECIES0
value of the number of species used to malloc data structures
Definition: vcs_solve.h:1491
std::vector< int > m_speciesUnknownType
Specifies the species unknown type.
Definition: vcs_solve.h:1629
std::vector< double > m_chargeSpecies
Charge of each species. Length = number of species.
Definition: vcs_solve.h:1949
std::vector< size_t > m_speciesLocalPhaseIndex
Index that keeps track of the index of the species within the local phase.
Definition: vcs_solve.h:1799
int vcs_prep_oneTime(int printLvl)
This routine is mostly concerned with changing the private data to be consistent with what's needed f...
Definition: vcs_prep.cpp:65
VCS_COUNTERS * m_VCount
Timing and iteration counters for the vcs object.
Definition: vcs_solve.h:1983
std::vector< double > m_spSize
total size of the species
Definition: vcs_solve.h:1573
Calculates reaction adjustments using a full Hessian approximation.
bool vcs_popPhasePossible(const size_t iphasePop) const
Utility function that evaluates whether a phase can be popped into existence.
std::vector< double > m_TmpPhase
Temporary vector of length NPhase.
Definition: vcs_solve.h:1726
size_t m_numRxnRdc
Current number of non-component species in the problem.
Definition: vcs_solve.h:1523
double m_totalVol
Total volume of all phases. Units are m^3.
Definition: vcs_solve.h:1970
int vcs_delete_species(const size_t kspec)
Change a single species from active to inactive status.
int delta_species(const size_t kspec, double *const delta_ptr)
Change the concentration of a species by delta moles.
double m_tolmaj
Tolerance requirement for major species.
Definition: vcs_solve.h:1762
int vcs_setMolesLinProg()
Estimate the initial mole numbers by constrained linear programming.
std::vector< double > m_deltaPhaseMoles
Change in the total moles in each phase.
Definition: vcs_solve.h:1735
int vcs(VCS_PROB *vprob, int ifunc, int ipr, int ip1, int maxit)
Solve an equilibrium problem.
Definition: vcs_solve.cpp:247
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...
A class for 2D double arrays stored in column-major (Fortran-compatible) form.
size_t NPHASE0
value of the number of phases used to malloc data structures
Definition: vcs_solve.h:1494
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.
std::vector< std::string > m_speciesName
Species string name for the kth species.
Definition: vcs_solve.h:1843
size_t m_numSpeciesRdc
Current number of species in the problems.
Definition: vcs_solve.h:1516
int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres)
Definition: vcs_TP.cpp:61