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