From 30b01ab8fdef40b2b08ecf75f7f8c2cacd7871d9 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Thu, 18 Apr 2013 22:07:59 +0000 Subject: [PATCH] Cleaned up Doxygen documentation for class VCS_SOLVE --- include/cantera/equil/vcs_solve.h | 1463 +++++++++++------------- src/equil/vcs_Gibbs.cpp | 25 - src/equil/vcs_TP.cpp | 65 +- src/equil/vcs_elem.cpp | 99 +- src/equil/vcs_elem_rearrange.cpp | 47 - src/equil/vcs_inest.cpp | 22 - src/equil/vcs_nondim.cpp | 51 - src/equil/vcs_phaseStability.cpp | 72 +- src/equil/vcs_prep.cpp | 74 -- src/equil/vcs_rank.cpp | 50 +- src/equil/vcs_rearrange.cpp | 8 - src/equil/vcs_report.cpp | 31 +- src/equil/vcs_rxnadj.cpp | 129 +-- src/equil/vcs_setMolesLinProg.cpp | 13 - src/equil/vcs_solve.cpp | 130 +-- src/equil/vcs_solve_TP.cpp | 612 +--------- src/equil/vcs_solve_phaseStability.cpp | 33 +- 17 files changed, 719 insertions(+), 2205 deletions(-) diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index 05eacd2f9..d122ccaa1 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -1,7 +1,7 @@ /** - * @file vcs_solve.h - * Header file for the internal object that holds the vcs equilibrium problem - * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). + * @file vcs_solve.h Header file for the internal object that holds the vcs + * equilibrium problem (see Class \link VCSnonideal::VCS_SOLVE + * VCS_SOLVE\endlink and \ref equilfunctions ). */ /* * Copyright (2005) Sandia Corporation. Under the terms of @@ -9,7 +9,6 @@ * U.S. Government retains certain rights in this software. */ - #ifndef _VCS_SOLVE_H #define _VCS_SOLVE_H @@ -43,40 +42,35 @@ class VCS_PROB; class VCS_COUNTERS; -//! This is the main structure used to hold the internal data -//! used in vcs_solve_TP(), and to solve TP systems. +//! This is the main structure used to hold the internal data +//! used in vcs_solve_TP(), and to solve TP systems. /*! - * The indices of information in this - * structure may change when the species basis changes or when - * phases pop in and out of existence. Both of these operations - * change the species ordering. - * + * The indices of information in this structure may change when the species + * basis changes or when phases pop in and out of existence. Both of these + * operations change the species ordering. */ class VCS_SOLVE { public: - //! Constructor for the VCS_SOLVE class VCS_SOLVE(); - //! Destructor ~VCS_SOLVE(); - //! Initialize the sizes within the VCS_SOLVE object /*! - * This resizes all of the internal arrays within the object. This routine - * operates in two modes. If all of the parameters are the same as it - * currently exists in the object, nothing is done by this routine; a quick - * exit is carried out and all of the data in the object persists. + * This resizes all of the internal arrays within the object. This + * routine operates in two modes. If all of the parameters are the same + * as it currently exists in the object, nothing is done by this routine; + * a quick exit is carried out and all of the data in the object + * persists. * - * IF any of the parameters are different than currently exists in the - * object, then all of the data in the object must be redone. It may not - * be zeroed, but it must be redone. + * If any of the parameters are different than currently exists in the + * object, then all of the data in the object must be redone. It may not + * be zeroed, but it must be redone. * * @param nspecies0 Number of species within the object * @param nelements Number of element constraints within the problem * @param nphase0 Number of phases defined within the problem. - * */ void vcs_initSizes(const size_t nspecies0, const size_t nelements, const size_t nphase0); @@ -127,113 +121,104 @@ public: * Any number of single-species phases and multi-species phases * can be handled by the present version. * - * Input - * ------------ - * @param print_lvl 1 -> Print results to standard output - * 0 -> don't report on anything + * @param print_lvl 1 -> Print results to standard output; + * 0 -> don't report on anything + * @param printDetails 1 -> Print intermediate results. + * @param maxit Maximum number of iterations for the algorithm * - * @param printDetails 1 -> Print intermediate results. - * - * @param maxit Maximum number of iterations for the algorithm - * - * @return 0 = Equilibrium Achieved - * 1 = Range space error encountered. The element abundance criteria are - * only partially satisfied. Specifically, the first NC= (number of - * components) conditions are satisfied. However, the full NE - * (number of elements) conditions are not satisfied. The equilibrium - * condition is returned. - * -1 = Maximum number of iterations is exceeded. Convergence was not - * found. + * @return + * * 0 = Equilibrium Achieved + * * 1 = Range space error encountered. The element abundance criteria + * are only partially satisfied. Specifically, the first NC= (number + * of components) conditions are satisfied. However, the full NE + * (number of elements) conditions are not satisfied. The + * equilibrium condition is returned. + * * -1 = Maximum number of iterations is exceeded. Convergence was + * not found. */ int vcs_solve_TP(int print_lvl, int printDetails, int maxit); - int vcs_PS(VCS_PROB* vprob, int iph, int printLvl, double& feStable); + /*! + * We make decisions on the initial mole number, and major-minor status + * here. We also fix up the total moles in a phase. + * + * irxn = id of the noncomponent species formation reaction for the + * species to be added in. + * + * The algorithm proceeds to implement these decisions in the previous + * position of the species. Then, vcs_switch_pos is called to move the + * species into the last active species slot, incrementing the number + * of active species at the same time. + * + * This routine is responsible for the global data manipulation only. + */ void vcs_reinsert_deleted(size_t kspec); - //! Choose the optimum species basis for the calculations + //! Choose the optimum species basis for the calculations /*! - * Choose the optimum component species basis for the calculations. - * This is done by choosing the species with the largest mole fraction - * not currently a linear combination of the previous components. - * Then, calculate the stoichiometric coefficient matrix for that - * basis. + * This is done by choosing the species with the largest mole fraction not + * currently a linear combination of the previous components. Then, + * calculate the stoichiometric coefficient matrix for that basis. * * Rearranges the solution data to put the component data at the * front of the species list. * - * Then, calculates m_stoichCoeffRxnMatrix[irxn][jcomp] the formation reactions - * for all noncomponent species in the mechanism. - * Also calculates DNG(I) and DNL(I), the net mole change for each - * formation reaction. - * Also, initializes IR(I) to the default state. + * Then, calculates m_stoichCoeffRxnMatrix[irxn][jcomp] the formation + * reactions for all noncomponent species in the mechanism. Also + * calculates DNG(I) and DNL(I), the net mole change for each formation + * reaction. Also, initializes IR(I) to the default state. * - * Input - * --------- - * @param doJustComponents If true, the m_stoichCoeffRxnMatrix[][] and - * m_deltaMolNumPhase[] are not calculated. + * @param[in] doJustComponents If true, the m_stoichCoeffRxnMatrix[][] and + * m_deltaMolNumPhase[] are not calculated. * - * @param aw Vector of mole fractions which will be used to construct an + * @param[in] aw Vector of mole fractions which will be used to construct an * optimal basis from. * - * @param sa Gram-Schmidt orthog work space (nc in length) sa[j] - * @param ss Gram-Schmidt orthog work space (nc in length) ss[j] - * @param sm QR matrix work space (nc*ne in length) sm[i+j*ne] - * @param test This is a small negative number dependent upon whether + * @param[in] sa Gram-Schmidt orthog work space (nc in length) sa[j] + * @param[in] ss Gram-Schmidt orthog work space (nc in length) ss[j] + * @param[in] sm QR matrix work space (nc*ne in length) sm[i+j*ne] + * @param[in] test This is a small negative number dependent upon whether * an estimate is supplied or not. - * - * Output - * --------- - * @param usedZeroedSpecies = If true, then a species with a zero concentration - * was used as a component. The problem may be - * converged. Or, the problem may have a range space - * error and may not have a proper solution. - * - * Internal Variables calculated by this routine: - * ----------------------------------------------- - * - * m_numComponents - * Number of component species - * - * component species - * This routine calculates the m_numComponent species. It switches - * their positions in the species vector so that they occupy - * the first m_numComponent spots in the species vector. - * - * m_stoichCoeffRxnMatrix[irxn][jcomp] - * Stoichiometric coefficient matrix for the reaction mechanism - * expressed in Reduced Canonical Form. - * jcomp refers to the component number, and irxn - * refers to the irxn_th non-component species. - * - * m_deltaMolNumPhase[irxn] - * Change in the number of total number of moles of species in all phases - * due to the noncomponent formation reaction, irxn. - * - * m_deltaMolNumPhase[irxn][iphase] - * Change in the number of moles in phase, iphase, due to the - * noncomponent formation reaction, irxn. - * - * m_phaseParticipation[irxn] - * This is 1 if the phase, iphase, participates in the - * formation reaction, irxn, and zero otherwise. - * + * @param[out] usedZeroedSpecies If true, then a species with a zero + * concentration was used as a component. + * The problem may be converged. Or, the + * problem may have a range space error and + * may not have a proper solution. * @return Returns VCS_SUCCESS if everything went ok. Returns * VCS_FAILED_CONVERGENCE if there is a problem. + * + * ### Internal Variables calculated by this routine: + * + * - #m_numComponents: Number of component species. This routine + * calculates the #m_numComponents species. It switches their positions + * in the species vector so that they occupy the first #m_numComponents + * spots in the species vector. + * - #m_stoichCoeffRxnMatrix[irxn][jcomp] Stoichiometric coefficient + * matrix for the reaction mechanism expressed in Reduced Canonical + * Form. jcomp refers to the component number, and irxn refers to the + * irxn_th non-component species. + * - #m_deltaMolNumPhase[irxn]: Change in the number of total number of + * moles of species in all phases due to the noncomponent formation + * reaction, irxn. + * - #m_deltaMolNumPhase[irxn][iphase]: Change in the number of moles in + * phase, iphase, due to the noncomponent formation reaction, irxn. + * - #m_phaseParticipation[irxn]: This is 1 if the phase, iphase, + * participates in the formation reaction, irxn, and zero otherwise. */ int vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], double ss[], double test, bool* const usedZeroedSpecies); //! Choose a species to test for the next component /*! - * We make the choice based on testing (molNum[i] * spSize[i]) for its maximum value. - * Preference for single species phases is also made. + * We make the choice based on testing (molNum[i] * spSize[i]) for its + * maximum value. Preference for single species phases is also made. * * @param molNum Mole number vector - * @param j index into molNum[] that indicates where the search will start from - * Previous successful components are swapped into the front of - * molNum[]. + * @param j index into molNum[] that indicates where the search + * will start from Previous successful components are + * swapped into the front of molNum[]. * @param n Length of molNum[] */ size_t vcs_basisOptMax(const double* const molNum, const size_t j, const size_t n); @@ -248,104 +233,120 @@ public: */ int vcs_species_type(const size_t kspec) const; + //! This routine evaluates the species type for all species + /*! + * - #VCS_SPECIES_MAJOR: Major species + * - #VCS_SPECIES_MINOR: Minor species + * - #VCS_SPECIES_SMALLMS: The species lies in a multicomponent phase + * that exists. Its concentration is currently very low, necessitating + * a different method of calculation. + * - #VCS_SPECIES_ZEROEDMS: The species lies in a multicomponent phase + * which currently doesn't exist. Its concentration is currently zero. + * - #VCS_SPECIES_ZEROEDSS: Species lies in a single-species phase which + * is currently zeroed out. + * - #VCS_SPECIES_DELETED: Species has such a small mole fraction it is + * deleted even though its phase may possibly exist. The species is + * believed to have such a small mole fraction that it best to throw + * the calculation of it out. It will be added back in at the end of + * the calculation. + * - #VCS_SPECIES_INTERFACIALVOLTAGE: Species refers to an electron in + * the metal The unknown is equal to the interfacial voltage drop + * across the interface on the SHE (standard hydrogen electrode) scale + * (volts). + * - #VCS_SPECIES_ZEROEDPHASE: Species lies in a multicomponent phase + * that is zeroed atm and will stay deleted due to a choice from a + * higher level. These species will formally always have zero mole + * numbers in the solution vector. + * - #VCS_SPECIES_ACTIVEBUTZERO: The species lies in a multicomponent + * phase which currently does exist. Its concentration is currently + * identically zero, though the phase exists. Note, this is a temporary + * condition that exists at the start of an equilibrium problem. The + * species is soon "birthed" or "deleted". + * - #VCS_SPECIES_STOICHZERO: The species lies in a multicomponent phase + * which currently does exist. Its concentration is currently + * identically zero, though the phase exists. This is a permanent + * condition due to stoich constraints + */ bool vcs_evaluate_speciesType(); //! We calculate the dimensionless chemical potentials of all species //! in a single phase. /*! - * We calculate the dimensionless chemical potentials of all species - * in a single phase. + * Note, for multispecies phases which are currently zeroed out, the + * chemical potential is filled out with the standard chemical potential. * - * Note, for multispecies phases which are currently zeroed out, - * the chemical potential is filled out with the standard chemical - * potential. - * - * For species in multispecies phases whose concentration is zero, - * we need to set the mole fraction to a very low value. - * Its chemical potential - * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration + * For species in multispecies phases whose concentration is zero, we need + * to set the mole fraction to a very low value. Its chemical potential is + * then calculated using the #VCS_DELETE_MINORSPECIES_CUTOFF concentration * to keep numbers positive. * - * Formula: - * --------------- + * # Formula: * - * Ideal Mixtures: + * ## Ideal Mixtures: * - * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * (This is equivalent to the adding the log of the mole fraction onto + * the standard chemical potential. ) + * + * ## Non-Ideal Mixtures: + * + * ### ActivityConvention = 0: molarity activity formulation + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * ( This is equivalent to the adding the log of the mole fraction + * multiplied by the activity coefficient onto the standard chemical + * potential. ) + * + * ### ActivityConvention = 1: molality activity formulation + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) + * - ln(Mnaught * m_units) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * + * Note: `m_SSfeSpecies(I)` is the molality based standard state. However, + * `ActCoeff[I]` is the molar based activity coefficient We have used the + * formulas: + * + * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] + * + * where `Xmol[N]` is the mole fraction of the solvent and `ActCoeff_M[I]` + * is the molality based act coeff. + * + * Note: This is equivalent to the "normal" molality formulation: + * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff_M[I] * m(I)) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase] + * + * where `m[I]` is the molality of the ith solute + * + * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) + * + * `z(I)/tPhMoles_ptr[iph] = Xmol[i]` is the mole fraction of i in the phase. * * - * ( This is equivalent to the adding the log of the - * mole fraction onto the standard chemical - * potential. ) + * NOTE: As per the discussion in vcs_dfe(), for small species where the + * mole fraction is small: * - * Non-Ideal Mixtures: - * ActivityConvention = 0: + * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * The chemical potential is calculated as: * - * ( This is equivalent to the adding the log of the - * mole fraction multiplied by the activity coefficient - * onto the standard chemical potential. ) - * - * ActivityConvention = 1: -> molality activity formulation - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * - ln(Mnaught * m_units) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * note: m_SSfeSpecies(I) is the molality based standard state. - * However, ActCoeff[I] is the molar based activity coefficient - * We have used the formulas; - * - * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] - * where Xmol[N] is the mole fraction of the solvent - * ActCoeff_M[I] is the molality based act coeff. - * - * note: This is equivalent to the "normal" molality formulation: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff_M[I] * m(I)) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase] - * where m[I] is the molality of the ith solute - * - * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) - * - * - * note: z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction - * of i in the phase. - * - * - * NOTE: - * As per the discussion in vcs_dfe(), for small species where the mole - * fraction is small: - * - * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF - * - * The chemical potential is calculated as: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) - * - * Input - * -------- - * iph : Phase to be calculated - * molNum(i) : Number of moles of species i - * (VCS species order) - * ff : standard state chemical potentials. These are the - * chemical potentials of the standard states at - * the same T and P as the solution. - * (VCS species order) - * Output - * ------- - * ac[] : Activity coefficients for species in phase - * (VCS species order) - * mu_i[] : Dimensionless chemical potentials for phase species - * (VCS species order) + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) * + * @param[in] iph Phase to be calculated + * @param[in] molNum(i) Number of moles of species i (VCS species order) + * @param[out] ac Activity coefficients for species in phase (VCS + * species order) + * @param[out] mu_i Dimensionless chemical potentials for phase + * species (VCS species order) */ void vcs_chemPotPhase(const int stateCalc, const size_t iph, const double* const molNum, double* const ac, double* const mu_i, @@ -363,129 +364,72 @@ public: * the chemical potential is filled out with the standard chemical * potential. * - * For species in multispecies phases whose concentration is zero, - * we need to set the mole fraction to a very low value. - * Its chemical potential - * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration + * For species in multispecies phases whose concentration is zero, we need + * to set the mole fraction to a very low value. Its chemical potential is + * then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration * to keep numbers positive. * - * - * Formula: - * --------------- - * - * Ideal Mixtures: - * - * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) - * + Charge[I] * Faraday_dim * phasePhi[iphase]; - * - * ( This is equivalent to the adding the log of the - * mole fraction onto the standard chemical - * potential. ) - * - * Non-Ideal Mixtures: -> molar activity formulation - * ActivityConvention = 0: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * + Charge[I] * Faraday_dim * phasePhi[iphase]; - * - * ( This is equivalent to the adding the log of the - * mole fraction multiplied by the activity coefficient - * onto the standard chemical potential. ) - * - * note: z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction - * of i in the phase. - * - * ActivityConvention = 1: -> molality activity formulation - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * - ln(Mnaught * m_units) - * + Charge[I] * Faraday_dim * phasePhi[iphase]; - * - * note: m_SSfeSpecies(I) is the molality based standard state. - * However, ActCoeff[I] is the molar based activity coefficient - * We have used the formulas; - * - * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] - * where Xmol[N] is the mole fraction of the solvent - * ActCoeff_M[I] is the molality based act coeff. - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff_M[I] * m(I)) - * + Charge[I] * Faraday_dim * phasePhi[iphase]; - * where m[I] is the molality of the ith solute - * - * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) - * + * For formulas, see vcs_chemPotPhase(). * * Handling of Small Species: * ------------------------------ - * As per the discussion above, for small species where the mole - * fraction + * As per the discussion above, for small species where the mole fraction * - * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF + * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF * - * The chemical potential is calculated as: + * The chemical potential is calculated as: * - * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) + * m_feSpecies(I)(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) * - * Species in the following categories are treated as "small species" + * Species in the following categories are treated as "small species" + * - #VCS_SPECIES_DELETED + * - #VCS_SPECIES_ACTIVEBUTZERO * - * - VCS_SPECIES_DELETED - * - VCS_SPECIES_ACTIVEBUTZERO - * . + * For species in multispecies phases which are currently not active, the + * treatment is different. These species are in the following species + * categories: + * - #VCS_SPECIES_ZEROEDMS + * - #VCS_SPECIES_ZEROEDPHASE * - * Handling of Small Species: - * ------------------------------ - * For species in multispecies phases which are currently not active, the - * treatment is different. These species are in the following species categories: + * For these species, the `ln( ActCoeff[I] X[I])` term is dropped + * altogether. The following equation is used: * - * - VCS_SPECIES_ZEROEDMS - * - VCS_SPECIES_ZEROEDPHASE - * . + * m_feSpecies(I) = m_SSfeSpecies(I) + * + Charge[I] * Faraday_dim * phasePhi[iphase]; * - * For these species, the ln( ActCoeff[I] X[I]) term is - * dropped altogether. The following equation is used. + * Handling of "Species" Representing Interfacial Voltages + * --------------------------------------------------------- * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + Charge[I] * Faraday_dim * phasePhi[iphase]; + * These species have species types of + * #VCS_SPECIES_TYPE_INTERFACIALVOLTAGE The chemical potentials for these + * "species" refer to electrons in metal electrodes. They have the + * following formula * + * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT * - * Handling of "Species" Representing Interfacial Voltages - * --------------------------------------------------------- + * - `F` is Faraday's constant. + * - `R` = gas constant + * - `T` = temperature + * - `V` = potential of the interface = phi_electrode - phi_solution * - * These species have species types of VCS_SPECIES_TYPE_INTERFACIALVOLTAGE - * The chemical potentials for these "species" refer to electrons in - * metal electrodes. They have the following formula + * For these species, the solution vector unknown, z[I], is V, the phase + * voltage, in volts. * - * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT - * - * F is Faraday's constant. - * R = gas constant - * T = temperature - * V = potential of the interface = phi_electrode - phi_solution - * - * For these species, the solution vector unknown, z[I], is V, the phase voltage, in volts. - * - * Input - * -------- * @param ll Determine which group of species gets updated - * ll = 0: Calculate for all species - * < 0: calculate for components and for major non-components - * 1: calculate for components and for minor non-components + * - `ll = 0`: Calculate for all species + * - `ll < 0`: calculate for components and for major non-components + * - `ll = 1`: calculate for components and for minor non-components * * @param lbot Restricts the calculation of the chemical potential * to the species between LBOT <= i < LTOP. Usually * LBOT and LTOP will be equal to 0 and MR, respectively. * @param ltop Top value of the loops * - * * @param stateCalc Determines whether z is old or new or tentative: - * 1: Use the tentative values for the total number of - * moles in the phases, i.e., use TG1 instead of TG etc. - * 0: Use the base values of the total number of - * moles in each system. + * - 1: Use the tentative values for the total number of + * moles in the phases, i.e., use TG1 instead of TG etc. + * - 0: Use the base values of the total number of + * moles in each system. * * Also needed: * ff : standard state chemical potentials. These are the @@ -515,6 +459,14 @@ public: //! Utility function that evaluates whether a phase can be popped //! into existence /*! + * A phase can be popped iff the stoichiometric coefficients for the + * component species, whose concentrations will be lowered during the + * process, are positive by at least a small degree. + * + * If one of the phase species is a zeroed component, then the phase can + * be popped if the component increases in mole number as the phase moles + * are increased. + * * @param iphasePop id of the phase, which is currently zeroed, * * @return Returns true if the phase can come into existence @@ -522,21 +474,15 @@ public: */ bool vcs_popPhasePossible(const size_t iphasePop) const; - //! Determine the list of problems that need to be checked to see if there are any phases pops /*! - * This routine evaluates and fills in the following quantities - * phasePopProblemLists_ + * This routine evaluates and fills in #phasePopProblemLists_. Need to + * work in species that are zeroed by element constraints. * * @return Returns the number of problems that must be checked. */ int vcs_phasePopDeterminePossibleList(); - - - - - //! Decision as to whether a phase pops back into existence /*! * @param phasePopPhaseIDs Vector containing the phase ids of the phases @@ -550,6 +496,10 @@ public: //! Calculates the deltas of the reactions due to phases popping //! into existence /*! + * Updates #m_deltaMolNumSpecies[irxn] : reaction adjustments, where irxn + * refers to the irxn'th species formation reaction. This adjustment is + * for species irxn + M, where M is the number of components. + * * @param iphasePop Phase id of the phase that will come into existence * * @return Returns an int representing the status of the step @@ -557,36 +507,37 @@ public: * - 1 : A single species phase species has been zeroed out * in this routine. The species is a noncomponent * - 2 : Same as one but, the zeroed species is a component. + * - 3 : Nothing was done because the phase couldn't be birthed + * because a needed component is zero. */ int vcs_popPhaseRxnStepSizes(const size_t iphasePop); //! Calculates formation reaction step sizes. /*! - * This is equation 6.4-16, p. 143 in Smith and Missen. + * This is equation 6.4-16, p. 143 in Smith and Missen. * * Output * ------- - * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers - * to the irxn'th species - * formation reaction. This adjustment is for species - * irxn + M, where M is the number of components. + * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers to + * the irxn'th species formation reaction. + * This adjustment is for species irxn + M, + * where M is the number of components. * * Special branching occurs sometimes. This causes the component basis * to be reevaluated * - * @param forceComponentCalc integer flagging whether a component recalculation needs - * to be carried out. - * @param kSpecial species number of phase being zeroed. + * @param forceComponentCalc integer flagging whether a component + * recalculation needs to be carried out. + * @param kSpecial species number of phase being zeroed. * * @return Returns an int representing which phase may need to be zeroed */ size_t vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial); - //! Calculates the total number of moles of species in all phases. + //! Calculates the total number of moles of species in all phases. /*! - * Calculates the total number of moles in all phases and updates - * the variable m_totalMolNum. - * Reconciles Phase existence flags with total moles in each phase. + * Also updates the variable m_totalMolNum and Reconciles Phase existence + * flags with total moles in each phase. */ double vcs_tmoles(); #ifdef DEBUG_MODE @@ -602,40 +553,37 @@ public: * coefficients for these reactions. A stoichiometric coefficient of * one is assumed for species irxn in this reaction. * - * INPUT * @param l - * L < 0 : Calculate reactions corresponding to - * major noncomponent and zeroed species only - * L = 0 : Do all noncomponent reactions, i, between - * 0 <= i < irxnl - * L > 0 : Calculate reactions corresponding to - * minor noncomponent and zeroed species only + * - `L < 0`: Calculate reactions corresponding to major noncomponent + * and zeroed species only + * - `L = 0`: Do all noncomponent reactions, i, between + * 0 <= i < irxnl + * - `L > 0`: Calculate reactions corresponding to minor noncomponent + * and zeroed species only * - * @param doDeleted Do deleted species - * @param vcsState Calculate deltaG corresponding to either old or new - * free energies - * @param alterZeroedPhases boolean indicating whether we should - * add in a special section for zeroed phases. + * @param doDeleted Do deleted species + * @param vcsState Calculate deltaG corresponding to either old or new + * free energies + * @param alterZeroedPhases boolean indicating whether we should + * add in a special section for zeroed phases. * - * Note we special case one important issue. - * If the component has zero moles, then we do not - * allow deltaG < 0.0 for formation reactions which - * would lead to the loss of more of that same component. - * This dG < 0.0 condition feeds back into the algorithm in several - * places, and leads to a infinite loop in at least one case. + * Note we special case one important issue. If the component has zero + * moles, then we do not allow deltaG < 0.0 for formation reactions which + * would lead to the loss of more of that same component. This dG < 0.0 + * condition feeds back into the algorithm in several places, and leads + * to a infinite loop in at least one case. */ void vcs_deltag(const int l, const bool doDeleted, const int vcsState, const bool alterZeroedPhases = true); void vcs_printDeltaG(const int stateCalc); - //! Calculate deltag of formation for all species in a single phase. + //! Calculate deltag of formation for all species in a single phase. /*! - * Calculate deltag of formation for all species in a single - * phase. It is assumed that the fe[] is up to date for all species. - * However, if the phase is currently zeroed out, a subproblem - * is calculated to solve for AC[i] and pseudo-X[i] for that - * phase. + * Calculate deltag of formation for all species in a single phase. It is + * assumed that the fe[] is up to date for all species. However, if the + * phase is currently zeroed out, a subproblem is calculated to solve for + * AC[i] and pseudo-X[i] for that phase. * * @param iphase phase index of the phase to be calculated * @param doDeleted boolean indicating whether to do deleted @@ -654,52 +602,66 @@ public: //! Swaps the indices for all of the global data for two species, k1 //! and k2. /*! - * - * @param ifunc: If true, switch the species data and the noncomponent reaction - * data. This must be called for a non-component species only. - * If false, switch the species data only. Typically, we use this - * option when determining the component species and at the - * end of the calculation, when we want to return unscrambled - * results. All rxn data will be out-of-date. + * @param ifunc: If true, switch the species data and the noncomponent + * reaction data. This must be called for a non-component + * species only. If false, switch the species data only. + * Typically, we use this option when determining the + * component species and at the end of the calculation, + * when we want to return unscrambled results. All rxn + * data will be out-of-date. * * @param k1 First species index - * * @param k2 Second species index */ void vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2); - - //! Birth guess returns the number of moles of a species - //! that is coming back to life. + //! Birth guess returns the number of moles of a species + //! that is coming back to life. /*! - * Birth guess returns the number of moles of a species - * that is coming back to life. - * Note, this routine is not applicable if the whole phase is coming - * back to life, not just one species in that phase. + * Birth guess returns the number of moles of a species that is coming + * back to life. Note, this routine is not applicable if the whole phase + * is coming back to life, not just one species in that phase. * - * Do a minor alt calculation. But, cap the mole numbers at - * 1.0E-15. - * For SS phases use VCS_DELETE_SPECIES_CUTOFF * 100. - * - * The routine makes sure the guess doesn't reduce the concentration - * of a component by more than 1/3. Note this may mean that - * the vlaue coming back from this routine is zero or a - * very small number. + * Do a minor alt calculation. But, cap the mole numbers at 1.0E-15. For + * SS phases use VCS_DELETE_SPECIES_CUTOFF * 100. * + * The routine makes sure the guess doesn't reduce the concentration of a + * component by more than 1/3. Note this may mean that the vlaue coming + * back from this routine is zero or a very small number. * * @param kspec Species number that is coming back to life - * - * @return Returns the number of kmol that the species should - * have. + * @return Returns the number of kmol that the species should have. */ double vcs_birthGuess(const int kspec); + //! Routine that independently determines whether a phase should be popped + //! under the current conditions. + /* + * This is the main routine that solves for equilibrium at constant T and + * P using a variant of the VCS method. Nonideal phases can be + * accommodated as well. Any number of single-species phases and multi- + * species phases can be handled by the present version. + * + * @param print_lvl 1 -> Print results to standard output; -> don't + * report on anything + * @param printDetails 1 -> Print intermediate results. + * @param maxit Maximum number of iterations for the algorithm + * + * @return + * - 0 = Equilibrium Achieved + * - 1 = Range space error encountered. The element abundance criteria are + * only partially satisfied. Specifically, the first NC= (number of + * components) conditions are satisfied. However, the full NE (number of + * elements) conditions are not satisfied. The equilibrium condition is + * returned. + * - -1 = Maximum number of iterations is exceeded. Convergence was not + * found. + */ int vcs_solve_phaseStability(const int iphase, int ifunc, double& funcval, int print_lvl); //! Main program to test whether a deleted phase should be brought //! back into existence /*! - * * @param iph Phase id of the deleted phase */ double vcs_phaseStabilityTest(const size_t iph); @@ -708,62 +670,65 @@ public: //! and pressure /*! * The actual problem statement is assumed to be in the structure - * already. This is a wrapper around the solve_TP() function. - * In this wrapper, we nondimensionalize the system - * we calculate the standard state gibbs free energies of the - * species, and we decide whether to we need to use the - * initial guess algorithm. + * already. This is a wrapper around the solve_TP() function. In this + * wrapper, we nondimensionalize the system we calculate the standard + * state gibbs free energies of the species, and we decide whether to we + * need to use the initial guess algorithm. * - * @param ipr = 1 -> Print results to standard output + * @param ipr = 1 -> Print results to standard output; * 0 -> don't report on anything - * @param ip1 = 1 -> Print intermediate results. + * @param ip1 = 1 -> Print intermediate results; * 0 -> Dont print any intermediate results * @param maxit Maximum number of iterations for the algorithm * @param T Value of the Temperature (Kelvin) - * Param pres Value of the Pressure (units given by m_VCS_UnitsFormat variable + * @param pres Value of the Pressure (units given by m_VCS_UnitsFormat variable * * @return Returns an integer representing the success of the algorithm - * 0 = Equilibrium Achieved - * 1 = Range space error encountered. The element abundance criteria are - * only partially satisfied. Specifically, the first NC= (number of - * components) conditions are satisfied. However, the full NE - * (number of elements) conditions are not satisfied. The equilibrium - * condition is returned. - * -1 = Maximum number of iterations is exceeded. Convergence was not - * found. + * * 0 = Equilibrium Achieved + * * 1 = Range space error encountered. The element abundance criteria are + * only partially satisfied. Specifically, the first NC= (number of + * components) conditions are satisfied. However, the full NE (number of + * elements) conditions are not satisfied. The equilibrium condition is + * returned. + * * -1 = Maximum number of iterations is exceeded. Convergence was not + * found. */ int vcs_TP(int ipr, int ip1, int maxit, double T, double pres); + /*! + * Evaluate the standard state free energies at the current temperature + * and pressure. Ideal gas pressure contribution is added in here. + * + * @param ipr 1 -> Print results to standard output; 0 -> don't report + * on anything + * @param ip1 1 -> Print intermediate results; 0 -> don't. + * @param Temp Temperature (Kelvin) + * @param pres Pressure (Pascal) + */ int vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres); //! Initialize the chemical potential of single species phases /*! - * For single species phases, initialize the chemical - * potential with the value of the standard state chemical - * potential. This value doesn't change during the calculation + * For single species phases, initialize the chemical potential with the + * value of the standard state chemical potential. This value doesn't + * change during the calculation */ void vcs_fePrep_TP(); //! Calculation of the total volume and the partial molar volumes /*! - * This function calculates the partial molar volume - * for all species, kspec, in the thermo problem - * at the temperature TKelvin and pressure, Pres, pres is in atm. - * And, it calculates the total volume of the combined system. + * This function calculates the partial molar volume for all species, + * kspec, in the thermo problem at the temperature TKelvin and pressure, + * Pres, pres is in atm. And, it calculates the total volume of the + * combined system. * - * Input - * --------------- - * @param tkelvin Temperature in kelvin() - * @param pres Pressure in Pascal - * @param w w[] is the vector containing the current mole numbers - * in units of kmol. - * - * Output - * ---------------- - * @param volPM[] For species in all phase, the entries are the + * @param[in] tkelvin Temperature in kelvin() + * @param[in] pres Pressure in Pascal + * @param[in] w w[] is the vector containing the current mole + * numbers in units of kmol. + * @param[out] volPM[] For species in all phase, the entries are the * partial molar volumes units of M**3 / kmol. - * - * @return The return value is the total volume of + * @return The return value is the total volume of * the entire system in units of m**3. */ double vcs_VolTotal(const double tkelvin, const double pres, @@ -779,28 +744,25 @@ public: * * The problem structure refers to: * - * the number and identity of the species. - * the formula matrix and thus the number of components. - * the number and identity of the phases. - * the equation of state - * the method and parameters for determining the standard state - * The method and parameters for determining the activity coefficients. + * - the number and identity of the species. + * - the formula matrix and thus the number of components. + * - the number and identity of the phases. + * - the equation of state + * - the method and parameters for determining the standard state + * - The method and parameters for determining the activity coefficients. * * Tasks: - * 0) Fill in the SSPhase[] array. - * 1) Check to see if any multispecies phases actually have only one + * 1. Fill in the SSPhase[] array. + * 2. Check to see if any multispecies phases actually have only one * species in that phase. If true, reassign that phase and species * to be a single-species phase. - * 2) Determine the number of components in the problem if not already + * 3. Determine the number of components in the problem if not already * done so. During this process the order of the species is changed * in the private data structure. All references to the species * properties must employ the ind[] index vector. * * @param printLvl Print level of the routine - * - * @return the return code - * VCS_SUCCESS = everything went OK - * + * @return VCS_SUCCESS = everything went OK */ int vcs_prep_oneTime(int printLvl); @@ -811,18 +773,18 @@ public: * every invocation of the vcs_solve() except for the cleanup invocation. * * Tasks: - * 1) Initialization of arrays to zero. + * 1. Initialization of arrays to zero. * - * return code - * VCS_SUCCESS = everything went OK + * @return + * VCS_SUCCESS = everything went OK; * VCS_PUB_BAD = There is an irreconcilable difference in the * public data structure from when the problem was * initially set up. */ int vcs_prep(); - //! In this routine, we check for things that will cause the algorithm - //! to fail. + //! In this routine, we check for things that will cause the algorithm + //! to fail. /*! * We check to see if the problem is well posed. If it is not, we return * false and print out error conditions. @@ -841,39 +803,33 @@ public: //! Rearrange the constraint equations represented by the Formula //! Matrix so that the operational ones are in the front /*! + * This subroutine handles the rearrangement of the constraint equations + * represented by the Formula Matrix. Rearrangement is only necessary when + * the number of components is less than the number of elements. For this + * case, some constraints can never be satisfied exactly, because the + * range space represented by the Formula Matrix of the components can't + * span the extra space. These constraints, which are out of the range + * space of the component Formula matrix entries, are migrated to the back + * of the Formula matrix. * - * This subroutine handles the rearrangement of the constraint - * equations represented by the Formula Matrix. Rearrangement is only - * necessary when the number of components is less than the number of - * elements. For this case, some constraints can never be satisfied - * exactly, because the range space represented by the Formula - * Matrix of the components can't span the extra space. These - * constraints, which are out of the range space of the component - * Formula matrix entries, are migrated to the back of the Formula - * matrix. + * A prototypical example is an extra element column in FormulaMatrix[], + * which is identically zero. For example, let's say that argon is has an + * element column in FormulaMatrix[], but no species in the mechanism + * actually contains argon. Then, nc < ne. Also, without perturbation of + * FormulaMatrix[] vcs_basopt[] would produce a zero pivot because the + * matrix would be singular (unless the argon element column was already + * the last column of FormulaMatrix[]. * - * A prototypical example is an extra element column in - * FormulaMatrix[], - * which is identically zero. For example, let's say that argon is - * has an element column in FormulaMatrix[], but no species in the - * mechanism - * actually contains argon. Then, nc < ne. Also, without perturbation - * of FormulaMatrix[] vcs_basopt[] would produce a zero pivot - * because the matrix - * would be singular (unless the argon element column was already the - * last column of FormulaMatrix[]. - * This routine borrows heavily from vcs_basopt's algorithm. It - * finds nc constraints which span the range space of the Component - * Formula matrix, and assigns them as the first nc components in the - * formula matrix. This guarantees that vcs_basopt[] has a - * nonsingular matrix to invert. + * This routine borrows heavily from vcs_basopt's algorithm. It finds nc + * constraints which span the range space of the Component Formula matrix, + * and assigns them as the first nc components in the formula matrix. This + * guarantees that vcs_basopt[] has a nonsingular matrix to invert. * * Other Variables - * @param aw aw[i[ Mole fraction work space (ne in length) - * @param sa sa[j] = Gram-Schmidt orthog work space (ne in length) - * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) - * @param ss ss[j] = Gram-Schmidt orthog work space (ne in length) - * + * @param aw Mole fraction work space (ne in length) + * @param sa Gram-Schmidt orthog work space (ne in length) + * @param sm QR matrix work space (ne*ne in length) + * @param ss Gram-Schmidt orthog work space (ne in length) */ int vcs_elem_rearrange(double* const aw, double* const sa, double* const sm, double* const ss); @@ -891,10 +847,9 @@ public: //! Calculates reaction adjustments using a full Hessian approximation /*! - * Calculates reaction adjustments. This does what equation 6.4-16, p. 143 - * in Smith and Missen is suppose to do. However, a full matrix is - * formed and then solved via a conjugate gradient algorithm. No - * preconditioning is done. + * This does what equation 6.4-16, p. 143 in Smith and Missen is supposed + * to do. However, a full matrix is formed and then solved via a conjugate + * gradient algorithm. No preconditioning is done. * * If special branching is warranted, then the program bails out. * @@ -938,20 +893,27 @@ public: */ double vcs_Hessian_actCoeff_diag(size_t irxn); + //! Recalculate all of the activity coefficients in all of the phases + //! based on input mole numbers + /* + * @param moleSpeciesVCS kmol of species to be used in the update. + * + * NOTE: This routine needs to be regulated. + */ void vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS); #ifdef DEBUG_MODE //! A line search algorithm is carried out on one reaction /*! - * In this routine we carry out a rough line search algorithm - * to make sure that the m_deltaGRxn_new doesn't switch signs prematurely. + * In this routine we carry out a rough line search algorithm to make + * sure that the m_deltaGRxn_new doesn't switch signs prematurely. * * @param irxn Reaction number * @param dx_orig Original step length * * @param ANOTE Output character string stating the conclusions of the * line search - * + * @return Returns the optimized step length found by the search */ double vcs_line_search(const size_t irxn, const double dx_orig, char* const ANOTE); @@ -959,7 +921,6 @@ public: double vcs_line_search(const size_t irxn, const double dx_orig); #endif - //! Print out a report on the state of the equilibrium problem to //! standard output. /*! @@ -1006,8 +967,6 @@ public: * or in nondimensional form. This routine switches the data from * dimensional form into nondimensional form. * - * What we do is to divide by RT. - * * @todo Add a scale factor based on the total mole numbers. * The algorithm contains hard coded numbers based on the * total mole number. If we ever were faced with a problem @@ -1023,8 +982,6 @@ public: * Essentially the internal data can either be in dimensional form * or in nondimensional form. This routine switches the data from * nondimensional form into dimensional form. - * - * What we do is to multiply by RT. */ void vcs_redim_TP(); @@ -1043,10 +1000,71 @@ public: */ void vcs_elab(); + /*! + * Checks to see if the element abundances are in compliance. If they are, + * then TRUE is returned. If not, FALSE is returned. Note the number of + * constraints checked is usually equal to the number of components in the + * problem. This routine can check satisfaction of all of the constraints + * in the problem, which is equal to ne. However, the solver can't fix + * breakage of constraints above nc, because that nc is the range space by + * definition. Satisfaction of extra constraints would have had to occur + * in the problem specification. + * + * The constraints should be broken up into 2 sections. If a constraint + * involves a formula matrix with positive and negative signs, and eaSet = + * 0.0, then you can't expect that the sum will be zero. There may be + * roundoff that inhibits this. However, if the formula matrix is all of + * one sign, then this requires that all species with nonzero entries in + * the formula matrix be identically zero. We put this into the logic + * below. + * + * @param ibound 1: Checks constraints up to the number of elements; + * 0: Checks constraints up to the number of components. + */ bool vcs_elabcheck(int ibound); - void vcs_elabPhase(size_t iphase, double* const elemAbundPhase); - int vcs_elcorr(double aa[], double x[]); + /*! + * Computes the elemental abundances vector for a single phase, + * elemAbundPhase[], and returns it through the argument list. The mole + * numbers of species are taken from the current value in + * m_molNumSpecies_old[]. + */ + void vcs_elabPhase(size_t iphase, double* const elemAbundPhase); + + /*! + * This subroutine corrects for element abundances. At the end of the + * subroutine, the total moles in all phases are recalculated again, + * because we have changed the number of moles in this routine. + * + * @param aa temporary work vector, length ne*ne + * @param x temporary work vector, length ne + * + * @returns + * - 0 = Nothing of significance happened, Element abundances were and + * still are good. + * - 1 = The solution changed significantly; The element abundances are + * now good. + * - 2 = The solution changed significantly, The element abundances are + * still bad. + * - 3 = The solution changed significantly, The element abundances are + * still bad and a component species got zeroed out. + * + * Internal data to be worked on:: + * - ga Current element abundances + * - m_elemAbundancesGoal Required elemental abundances + * - m_molNumSpecies_old Current mole number of species. + * - m_formulaMatrix[][] Formula matrix of the species + * - ne Number of elements + * - nc Number of components. + * + * NOTES: + * This routine is turning out to be very problematic. There are + * lots of special cases and problems with zeroing out species. + * + * Still need to check out when we do loops over nc vs. ne. + * + */ + int vcs_elcorr(double aa[], double x[]); //! Create an initial estimate of the solution to the thermodynamic //! equilibrium problem. @@ -1060,24 +1078,28 @@ public: //! Estimate the initial mole numbers by constrained linear programming /*! - * This is done by running - * each reaction as far forward or backward as possible, subject - * to the constraint that all mole numbers remain - * non-negative. Reactions for which \f$ \Delta \mu^0 \f$ are - * positive are run in reverse, and ones for which it is negative - * are run in the forward direction. The end result is equivalent - * to solving the linear programming problem of minimizing the - * linear Gibbs function subject to the element and - * non-negativity constraints. + * This is done by running each reaction as far forward or backward as + * possible, subject to the constraint that all mole numbers remain non- + * negative. Reactions for which \f$ \Delta \mu^0 \f$ are positive are run + * in reverse, and ones for which it is negative are run in the forward + * direction. The end result is equivalent to solving the linear + * programming problem of minimizing the linear Gibbs function subject to + * the element and non-negativity constraints. */ int vcs_setMolesLinProg(); + //! Calculate the total dimensionless Gibbs free energy + /*! + * Inert species are handled as if they had a standard free energy of + * zero. Note, for this algorithm this function should be monotonically + * decreasing. + */ double vcs_Total_Gibbs(double* w, double* fe, double* tPhMoles); //! Calculate the total dimensionless Gibbs free energy of a single phase /*! - * -> Inert species are handled as if they had a standard free - * energy of zero and if they obeyed ideal solution/gas theory + * Inert species are handled as if they had a standard free energy of + * zero and if they obeyed ideal solution/gas theory. * * @param iphase ID of the phase * @param w Species mole number vector for all species @@ -1089,7 +1111,7 @@ public: //! Transfer the results of the equilibrium calculation back to VCS_PROB /*! - * The VCS_PUB structure is returned to the user. + * The VCS_PROB structure is returned to the user. * * @param pub Pointer to VCS_PROB object that will get the results of the * equilibrium calculation transfered to it. @@ -1122,10 +1144,10 @@ private: //! Zero out the concentration of a species. /*! - * Zero out the concentration of a species. Make sure to conserve - * elements and keep track of the total moles in all phases. - * w[] - * m_tPhaseMoles_old[] + * Make sure to conserveelements and keep track of the total moles in all + * phases. + * - w[] + * - m_tPhaseMoles_old[] * * @param kspec Species index * @@ -1171,11 +1193,12 @@ private: //! Change the concentration of a species by delta moles. /*! - * Make sure to conserve elements and keep track of the total kmoles in all phases. + * Make sure to conserve elements and keep track of the total kmoles in + * all phases. * * @param kspec The species index - * @delta_ptr pointer to the delta for the species. This may change during - * the calculation + * @param delta_ptr pointer to the delta for the species. This may + * change during the calculation * * @return * 1: succeeded without change of dx @@ -1216,25 +1239,23 @@ private: //! Recheck deletion condition for multispecies phases. /*! - * We assume here that DG_i_0 has been calculated for deleted species correctly + * We assume here that DG_i_0 has been calculated for deleted species correctly * + * m_feSpecies(I) = m_SSfeSpecies(I) + * + ln(ActCoeff[I]) + * - ln(Mnaught * m_units) + * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I]) - * - ln(Mnaught * m_units) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; + * sum_u = sum_j_comp [ sigma_i_j * u_j ] + * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) * - * sum_u = sum_j_comp [ sigma_i_j * u_j ] - * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) + * DG_i_0 = m_feSpecies(I) - sum_m{ a_i_m DG_m } * - * DG_i_0 = m_feSpecies(I) - sum_m{ a_i_m DG_m } + * by first evaluating: * + * DG_i_O = u_i_O - sum_u. * - * by first evaluating: - * - * DG_i_O = u_i_O - sum_u. - * - * Then, the phase pops into existence iff + * Then, the phase pops into existence iff * * phaseDG = 1.0 - sum_i{exp(-DG_i_O)} < 0.0 * @@ -1243,59 +1264,53 @@ private: * * @param iphase Phase index number * - * @return Returns true if the phase is currently deleted - * but should be reinstated. Returns false otherwise. + * @return Returns true if the phase is currently deleted + * but should be reinstated. Returns false otherwise. * * NOTE: this routine is currently not used in the code, and * contains some basic changes that are incompatible. * - * assumptions: - * 1) Vphase Existence is up to date - * 2) Vphase->IndSpecies is up to date - * 3) m_deltaGRxn_old[irxn] is up to date + * assumptions: + * 1. Vphase Existence is up to date + * 2. Vphase->IndSpecies is up to date + * 3. m_deltaGRxn_old[irxn] is up to date */ bool recheck_deleted_phase(const int iphase); - //! Minor species alternative calculation + //! Minor species alternative calculation /*! - * This is based upon the following approximation: - * The mole fraction changes due to these reactions don't affect - * the mole numbers of the component species. Therefore the following - * approximation is valid for a small component of an ideal phase: + * This is based upon the following approximation: The mole fraction + * changes due to these reactions don't affect the mole numbers of the + * component species. Therefore the following approximation is valid for + * a small component of an ideal phase: * * 0 = m_deltaGRxn_old(I) + log(molNum_new(I)/molNum_old(I)) * - * m_deltaGRxn_old contains the contribution from + * `m_deltaGRxn_old` contains the contribution from * - * m_feSpecies_old(I) = + * m_feSpecies_old(I) = * m_SSfeSpecies(I) + * log(ActCoeff[i] * molNum_old(I) / m_tPhaseMoles_old(iph)) - * Thus, + * Thus, * - * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I)) + * molNum_new(I)= molNum_old(I) * EXP(-m_deltaGRxn_old(I)) * - * Most of this section is mainly restricting the update to reasonable - * values. - * We restrict the update a factor of 1.0E10 up and 1.0E-10 down - * because we run into trouble with the addition operator due to roundoff - * if we go larger than ~1.0E15. Roundoff will then sometimes produce - * zero mole fractions. + * Most of this section is mainly restricting the update to reasonable + * values. We restrict the update a factor of 1.0E10 up and 1.0E-10 down + * because we run into trouble with the addition operator due to roundoff + * if we go larger than ~1.0E15. Roundoff will then sometimes produce + * zero mole fractions. * - * Note: This routine was generalized to incorporate - * nonideal phases and phases on the molality basis + * Note: This routine was generalized to incorporate nonideal phases and + * phases on the molality basis * - * Input: - * ------ - * @param kspec The current species and corresponding formation - * reaction number. - * @param irxn The current species and corresponding formation - * reaction number. - * - * Output: - * --------- - * @param do_delete: BOOLEAN which if true on return, then we branch - * to the section that deletes a species from the - * current set of active species. + * @param[in] kspec The current species and corresponding formation + * reaction number. + * @param[in] irxn The current species and corresponding formation + * reaction number. + * @param[out] do_delete: BOOLEAN which if true on return, then we + * branch to the section that deletes a species + * from the current set of active species. */ double vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete #ifdef DEBUG_MODE @@ -1303,40 +1318,38 @@ private: #endif ) const; - //! This routine optimizes the minimization of the total gibbs free - //! energy by making sure the slope of the following functional stays - //! negative: + //! This routine optimizes the minimization of the total gibbs free energy + //! by making sure the slope of the following functional stays negative: /*! * The slope of the following functional is equivalent to the slope * of the total Gibbs free energy of the system: * - * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] ) + * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] ) * * along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0 0), * does this code section kick in. It finds the point on the parabola * where the slope is equal to zero. - * */ bool vcs_globStepDamp(); //! Switch rows and columns of a square matrix /*! - * Switches the row and column of a matrix. - * So that after + * Switches the row and column of a matrix. So that after * - * J[k1][j] = J_old[k2][j] and J[j][k1] = J_old[j][k2] - * J[k2][j] = J_old[k1][j] and J[j][k2] = J_old[j][k1] + * J[k1][j] = J_old[k2][j] and J[j][k1] = J_old[j][k2] + * J[k2][j] = J_old[k1][j] and J[j][k2] = J_old[j][k1] * - * @param Jac Double pointer to the Jacobian - * @param k1 first row/column value to be switched - * @param k2 second row/column value to be switched + * @param Jac Double pointer to the Jacobian + * @param k1 first row/column value to be switched + * @param k2 second row/column value to be switched */ void vcs_switch2D(double* const* const Jac, const size_t k1, const size_t k2) const; @@ -1359,7 +1372,6 @@ private: //! Estimate equilibrium compositions /*! - * Estimates equilibrium compositions. * Algorithm covered in a section of Smith and Missen's Book. * * Linear programming module is based on using dbolm. @@ -1373,33 +1385,30 @@ private: void vcs_inest(double* const aw, double* const sa, double* const sm, double* const ss, double test); - //! Calculate the status of single species phases. void vcs_SSPhase(void); //! This function recalculates the deltaG for reaction, irxn /*! - * This function recalculates the deltaG for reaction irxn, - * given the mole numbers in molNum. It uses the temporary - * space mu_i, to hold the recalculated chemical potentials. - * It only recalculates the chemical potentials for species in phases - * which participate in the irxn reaction. + * This function recalculates the deltaG for reaction irxn, given the mole + * numbers in molNum. It uses the temporary space mu_i, to hold the + * recalculated chemical potentials. It only recalculates the chemical + * potentials for species in phases which participate in the irxn + * reaction. This function is used by the vcs_line_search algorithm() and + * should not be used widely due to the unknown state it leaves the + * system. * - * Input - * ------------ - * @param irxn Reaction number - * @param molNum Current mole numbers of species to be used as - * input to the calculation (units = kmol) - * (length = totalNuMSpecies) + * @param[in] irxn Reaction number + * @param[in] molNum Current mole numbers of species to be used as input + * to the calculation (units = kmol), (length = totalNuMSpecies) * - * Output - * ------------ - * @param ac output Activity coefficients (length = totalNumSpecies) - * Note this is only partially formed. Only species in - * phases that participate in the reaction will be updated - * @param mu_i dimensionless chemical potentials (length - totalNumSpecies - * Note this is only partially formed. Only species in - * phases that participate in the reaction will be updated + * @param[out] ac Activity coefficients (length = totalNumSpecies) Note + * this is only partially formed. Only species in phases that + * participate in the reaction will be updated + * @param[out] mu_i dimensionless chemical potentials (length + * totalNumSpecies) Note this is only partially formed. Only + * species in phases that participate in the reaction will be + * updated * * @return Returns the dimensionless deltaG of the reaction */ @@ -1446,49 +1455,39 @@ private: void vcs_updateMolNumVolPhases(const int stateCalc); public: - //! Calculate the rank of a matrix and return the rows and columns that will generate an independent basis - //! for that rank + //! Calculate the rank of a matrix and return the rows and columns that + //! will generate an independent basis for that rank /* - * Choose the optimum component species basis for the calculations, finding the rank and - * set of linearly independent rows for that calculation. - * Then find the set of linearly indepedent element columns that can support that rank. - * This is done by taking the transpose of the matrix and redoing the same calculation. - * (there may be a better way to do this. I don't know.) + * Choose the optimum component species basis for the calculations, + * finding the rank and set of linearly independent rows for that + * calculation. Then find the set of linearly indepedent element columns + * that can support that rank. This is done by taking the transpose of the + * matrix and redoing the same calculation. (there may be a better way to + * do this. I don't know.) * + * @param[in] awtmp Vector of mole numbers which will be used to + * construct a ranking for how to pick the basis species. This is + * largely ignored here. + * @param[in] numSpecies Number of species. This is the number of rows in + * the matrix. + * @param[in] matrix This is the formula matrix. Nominally, the rows are + * species, while the columns are element compositions. However, + * this routine is totally general, so that the rows and columns + * can be anything. + * @param[in] numElemConstraints Number of element constraints * - * Input - * --------- + * @param[out] usedZeroedSpecies If true, then a species with a zero + * concentration was used as a component. + * @param[out] compRes Vector of rows which are linearly independent. + * (these are the components) + * @param[out] elemComp Vector of columns which are linearly independent + * (These are the actionable element constraints). * - * @param awtmp Vector of mole numbers which will be used to construct a - * ranking for how to pick the basis species. This is largely ignored - * here. - * - * @param numSpecies Number of species. This is the number of rows in the matrix. - * - * @param matrix Matrix. This is the formula matrix. Nominally, the rows are species, while - * the columns are element compositions. However, this routine - * is totally general, so that the rows and columns can be anything. - * - * @param numElemConstraints Number of element constraints - * - * Output - * --------- - * @param usedZeroedSpecies = If true, then a species with a zero concentration - * was used as a component. - * - * - * @param compRes Vector of rows which are linearly independent. (these are the components) - * - * @param elemComp Vector of columns which are linearly independent (These are the actionable element - * constraints). - * - * @return Returns number of components. This is the rank of the matrix + * @return Returns number of components. This is the rank of the matrix */ int vcs_rank(const double* awtmp, size_t numSpecies, const double* matrix, size_t numElemConstraints, std::vector &compRes, std::vector &elemComp, int* const usedZeroedSpecies) const; - -public: //! value of the number of species used to malloc data structures size_t NSPECIES0; @@ -1545,7 +1544,7 @@ public: * reaction which forms species kspec from the component species. A * stoichiometric coefficient of one is assumed for the species kspec in this mechanism. * - * NOTE: kspec = irxn + m_numComponents + * NOTE: kspec = irxn + m_numComponents * * m_stoichCoeffRxnMatrix[irxn][j] : * j refers to the component number, and irxn refers to the irxn_th non-component species. @@ -1600,42 +1599,40 @@ public: //! Setting for whether to do an initial estimate /*! - * Initial estimate: 0 Do not estimate the solution at all. Use the - * supplied mole numbers as is. - * 1 Only do an estimate if the element abundances - * aren't satisfied. - * -1 Force an estimate of the soln. Throw out the input - * mole numbers. + * Initial estimate: + * * 0 Do not estimate the solution at all. Use the + * * supplied mole numbers as is. + * * 1 Only do an estimate if the element abundances + * * aren't satisfied. + * * -1 Force an estimate of the soln. Throw out the input + * * mole numbers. */ int m_doEstimateEquil; //! Total moles of the species /*! - * Total number of moles of the kth species. - * Length = Total number of species = m + * Total number of moles of the kth species. + * Length = Total number of species = m */ std::vector m_molNumSpecies_old; //! Specifies the species unknown type /*! - * There are two types. One is the straightforward - * species, with the mole number w[k], as the - * unknown. The second is the an interfacial - * voltage where w[k] refers to the interfacial - * voltage in volts. - * These species types correspond to metallic - * electrons corresponding to electrodes. - * The voltage and other interfacial conditions - * sets up an interfacial current, which is - * set to zero in this initial treatment. - * Later we may have non-zero interfacial currents. + * There are two types. One is the straightforward species, with the mole + * number w[k], as the unknown. The second is the an interfacial voltage + * where w[k] refers to the interfacial voltage in volts. + * + * These species types correspond to metallic electrons corresponding to + * electrodes. The voltage and other interfacial conditions sets up an + * interfacial current, which is set to zero in this initial treatment. + * Later we may have non-zero interfacial currents. */ std::vector m_speciesUnknownType; - //! Change in the number of moles of phase, iphase, due to the noncomponent formation - //! reaction, irxn, for species, k: + //! Change in the number of moles of phase, iphase, due to the + //! noncomponent formation reaction, irxn, for species, k: /*! - * m_deltaMolNumPhase[irxn][iphase] = k = nc + irxn + * m_deltaMolNumPhase[irxn][iphase] = k = nc + irxn */ DoubleStarStar m_deltaMolNumPhase; @@ -1647,7 +1644,7 @@ public: std::vector m_phasePhi; //! Tentative value of the mole number vector. It's also used to store the - //! mole fraction vector. + //! mole fraction vector. //std::vector wt; std::vector m_molNumSpecies_new; @@ -1695,20 +1692,17 @@ public: //! Element abundances vector Goals /*! - * Vector of moles of each element that are the goals of the - * simulation. This is a constant in the problem. - * Note other constraint conditions are added to this vector. - * This is input from the input file and - * is considered a constant from thereon. - * units = kmoles + * Vector of moles of each element that are the goals of the simulation. + * This is a constant in the problem. Note other constraint conditions + * are added to this vector. This is input from the input file and is + * considered a constant from thereon. units = kmoles */ std::vector m_elemAbundancesGoal; //! Total number of kmoles in all phases /*! * This number includes the inerts. - * -> Don't use this except for scaling - * purposes + * -> Don't use this except for scaling purposes */ double m_totalMolNum; @@ -1746,13 +1740,15 @@ public: //! Pressure (units are determined by m_VCS_UnitsFormat /*! - * Values units - * -1: atm - * 0: atm - * 1: atm - * 2: atm - * 3: Pa - * Units being changed to Pa + * | Values | units | + * | ------ | ----- + * | -1: | atm | + * | 0: | atm | + * | 1: | atm | + * | 2: | atm | + * | 3: | Pa | + * + * Units being changed to Pa */ double m_pressurePA; @@ -1805,90 +1801,35 @@ public: //! Index vector that keeps track of the rearrangement of the elements /*! - * At the end of each run, the element vector and associated data gets put back - * in the original order. + * At the end of each run, the element vector and associated data gets + * put back in the original order. * * Example * - * e = m_elementMapIndex[eNum] - * - * eNum = current order in the vcs_solve object - * e = original order in the vcs_prob object and in the MultiPhase object + * e = m_elementMapIndex[eNum] + * eNum = current order in the vcs_solve object + * e = original order in the vcs_prob object and in the MultiPhase object */ std::vector m_elementMapIndex; //! Mapping between the species index for noncomponent species and the //! full species index. /*! - * ir[irxn] = Mapping between the reaction index for - * noncomponent formation reaction of a species - * and the full species + * ir[irxn] = Mapping between the reaction index for noncomponent + * formation reaction of a species and the full species * index. - * - Initially set to a value of K = NC + I - * This vector has length equal to number - * of noncomponent species in the mechanism. - * It starts with the first current - * noncomponent species in the mechanism. - * kspec = ir[irxn] + * + * Initially set to a value of K = NC + I This vector has length equal to + * number of noncomponent species in the mechanism. It starts with the + * first current noncomponent species in the mechanism. kspec = ir[irxn] */ std::vector m_indexRxnToSpecies; //! Major -Minor status vector for the species in the problem /*! - * The index for this is species. The reaction that this is referring - * to is - * kspec = irxn + m_numComponents - * - * kspec : 2 -> Component species VCS_SPECIES_COMPONENT - * -> deprecated, want to assign -2 to some - * component species. We can already determine - * whether the species is a component from - * its position in the species vector. - * 1 -> Major species VCS_SPECIES_MAJOR - * 0 -> Minor species VCS_SPECIES_MINOR - * -1 -> The species lies in a multicomponent phase - * that exists. Its concentration is currently - * very low, necessitating a different method - * of calculation. - * - VCS_SPECIES_ZEROEDPHASE - * -2 -> The species lies in a multicomponent phase - * which currently doesn't exist. - * Its concentration is currently zero. - * - VCS_SPECIES_ZEROEDMS - * -3 -> Species lies in a single-species phase which - * is currently zeroed out. - * - VCS_SPECIES_ZEROEDSS - * -4 -> Species has such a small mole fraction it is - * deleted even though its phase may possibly exist. - * The species is believed to have such a small - * mole fraction that it best to throw the - * calculation of it out. It will be added back in - * at the end of the calculation. - * - VCS_SPECIES_DELETED - * -5 -> Species refers to an electron in the metal - * The unknown is equal to the interfacial voltage - * drop across the interface on the SHE (standard - * hydrogen electrode) scale (volts). - * - VCS_SPECIES_INTERFACIALVOLTAGE - * -6 -> Species lies in a multicomponent phase that - * is zeroed atm and will stay deleted due to a - * choice from a higher level. - * These species will formally always have zero - * mole numbers in the solution vector. - * - VCS_SPECIES_ZEROEDPHASE - * -7 -> The species lies in a multicomponent phase which - * currently does exist. Its concentration is currently - * identically zero, though the phase exists. Note, this - * is a temporary condition that exists at the start - * of an equilibrium problem. - * The species is soon "birthed" or "deleted". - * - VCS_SPECIES_ACTIVEBUTZERO - * -8 -> The species lies in a multicomponent phase which - * currently does exist. Its concentration is currently - * identically zero, though the phase exists. This is - * a permanent condition due to stoich constraints - * - VCS_SPECIES_STOICHZERO - * + * The index for this is species. The reaction that this is referring to + * is `kspec = irxn + m_numComponents`. For possible values and their + * meanings, see vcs_evaluate_speciesType(). */ std::vector m_speciesStatus; @@ -1900,9 +1841,6 @@ public: std::vector m_SSPhase; //! Species string name for the kth species - /*! - * Species string name for the kth species - */ std::vector m_speciesName; //! Vector of strings containing the element names @@ -1913,41 +1851,35 @@ public: //! Type of the element constraint /*! - * m_elType[j] = type of the element - * 0 VCS_ELEM_TYPE_ABSPOS Normal element that is positive - * or zero in all species. - * 1 VCS_ELEM_TPYE_ELECTRONCHARGE element dof that corresponds - * to the electronic charge DOF. - * 2 VCS_ELEM_TYPE_CHARGENEUTRALITY element dof that - * corresponds to a required charge - * neutrality constraint on the phase. - * The element abundance is always exactly zero. - * 3 VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may - * mean that a species has neg 0 or pos value - * of that constraint (other than charge) + * * 0 - #VCS_ELEM_TYPE_ABSPOS Normal element that is positive or zero in + * all species. + * * 1 - #VCS_ELEM_TYPE_ELECTRONCHARGE element dof that corresponds to + * the electronic charge DOF. + * * 2 - #VCS_ELEM_TYPE_CHARGENEUTRALITY element dof that corresponds to + * a required charge neutrality constraint on the phase. The element + * abundance is always exactly zero. + * * 3 - #VCS_ELEM_TYPE_OTHERCONSTRAINT Other constraint which may mean + * that a species has neg 0 or pos value of that constraint (other than + * charge) */ std::vector m_elType; //! Specifies whether an element constraint is active /*! - * The default is true - * Length = nelements + * The default is true. Length = nelements */ std::vector m_elementActive; - //! Array of Phase Structures - /*! - * Length = number of phases - */ + //! Array of Phase Structures. Length = number of phases. std::vector m_VolPhaseList; //! String containing the title of the run std::string m_title; - //! This specifies the current state of units for the Gibbs free energy - //! properties in the program. + //! This specifies the current state of units for the Gibbs free energy + //! properties in the program. /*! - *. The default is to have this unitless + * The default is to have this unitless */ char m_unitsState; @@ -1962,41 +1894,37 @@ public: //! specifies the activity convention of the phase containing the species /*! - * SpecActConvention[kspec] - * 0 = molar based - * 1 = molality based - * length = number of species + * * 0 = molar based + * * 1 = molality based + * + * length = number of species */ std::vector m_actConventionSpecies; //! specifies the activity convention of the phase. /*! - * 0 = molar based - * 1 = molality based - * length = number of phases + * * 0 = molar based + * * 1 = molality based + * + * length = number of phases */ std::vector m_phaseActConvention; //! specifies the ln(Mnaught) used to calculate the chemical potentials /*! - * For molar based activity conventions - * this will be equal to 0.0 - * length = number of species + * For molar based activity conventions this will be equal to 0.0. + * length = number of species. */ std::vector m_lnMnaughtSpecies; - //! Molar-based Activity Coefficients for Species - /*! - * - * Length = number of species - */ + //! Molar-based Activity Coefficients for Species. + //! Length = number of species std::vector m_actCoeffSpecies_new; //! Molar-based Activity Coefficients for Species based on old mole numbers /*! - * These activity coefficients are based on the m_molNumSpecies_old values - * Molar based activity coeffients. - * Length = number of species + * These activity coefficients are based on the m_molNumSpecies_old + * values Molar based activity coeffients. Length = number of species */ std::vector m_actCoeffSpecies_old; @@ -2012,17 +1940,13 @@ public: //! Molecular weight of each species /*! - * units = kg/kmol - * length = number of species + * units = kg/kmol. length = number of species. * * note: this is a candidate for removal. I don't think we use it. */ std::vector m_wtSpecies; - //! Charge of each species - /*! - * Length = number of species - */ + //! Charge of each species. Length = number of species. std::vector m_chargeSpecies; std::vector > phasePopProblemLists_; @@ -2043,10 +1967,7 @@ public: */ int m_useActCoeffJac; - //! Total volume of all phases - /*! - * units are m^3 - */ + //! Total volume of all phases. Units are m^3 double m_totalVol; //! Partial molar volumes of the species @@ -2056,27 +1977,23 @@ public: */ std::vector m_PMVolumeSpecies; - //! dimensionless value of Faraday's constant - /*! - * F / RT (1/volt) - */ + //! dimensionless value of Faraday's constant, F / RT (1/volt) double m_Faraday_dim; //! Timing and iteration counters for the vcs object VCS_COUNTERS* m_VCount; - //! Debug printing lvl /*! * Levels correspond to the following guidlines - * - 0 No printing at all - * - 1 Serious warnings or fatal errors get one line - * - 2 one line per eacdh successful vcs package call - * - 3 one line per every successful solve_TP calculation - * - 4 one line for every successful operation -> solve_TP gets a summary report - * - 5 each iteration in solve_TP gets a report with one line per species - * - 6 Each decision in solve_TP gets a line per species in addition to 4 - * - 10 Additionally Hessian matrix is printed out + * * 0 No printing at all + * * 1 Serious warnings or fatal errors get one line + * * 2 one line per eacdh successful vcs package call + * * 3 one line per every successful solve_TP calculation + * * 4 one line for every successful operation -> solve_TP gets a summary report + * * 5 each iteration in solve_TP gets a report with one line per species + * * 6 Each decision in solve_TP gets a line per species in addition to 4 + * * 10 Additionally Hessian matrix is printed out * * Levels of printing above 4 are only accessible when DEBUG_MODE is turned on */ @@ -2084,33 +2001,25 @@ public: //! printing level of timing information /*! - * 1 allowing printing of timing - * 0 do not allow printing of timing -> everything is printed - * as a NA. + * * 1 allowing printing of timing + * * 0 do not allow printing of timing -> everything is printed as a NA. */ int m_timing_print_lvl; - //! Units for the chemical potential data: + //! Units for the chemical potential data /*! - * VCS_UnitsFormat = Units for the chemical potential data: - * -1: kcal/mol - * 0: MU/RT - * 1: kJ/mol - * 2: Kelvin - * 3: J / kmol - * and pressure data: - * -1: Pa - * 0: Pa - * 1: Pa - * 2: pa - * 3: Pa + * | Value | chemical potential units | pressure units | + * | ----- | ------------------------ | -------------- | + * | -1 | kcal/mol | Pa | + * | 0 | MU/RT | Pa | + * | 1 | kJ/mol | Pa | + * | 2 | Kelvin | Pa | + * | 3 | J / kmol | Pa | */ int m_VCS_UnitsFormat; friend class vcs_phaseStabilitySolve; - }; } #endif - diff --git a/src/equil/vcs_Gibbs.cpp b/src/equil/vcs_Gibbs.cpp index 62735d0ad..8541f096d 100644 --- a/src/equil/vcs_Gibbs.cpp +++ b/src/equil/vcs_Gibbs.cpp @@ -15,24 +15,8 @@ namespace VCSnonideal { - -/***************************************************************************/ -/***************************************************************************/ -/***************************************************************************/ - double VCS_SOLVE::vcs_Total_Gibbs(double* molesSp, double* chemPot, double* tPhMoles) - -/************************************************************************* - * - * vcs_Total_Gibbs: - * - * Calculate the total dimensionless Gibbs free energy - * -> Inert species are handled as if they had a standard free - * energy of zero. - * Note, for this algorithm this function should be MONOTONICALLY - * DECREASING. - *************************************************************************/ { double g = 0.0; @@ -56,15 +40,6 @@ double VCS_SOLVE::vcs_Total_Gibbs(double* molesSp, double* chemPot, return g; } -// Calculate the total dimensionless Gibbs free energy of a single phase -/* - * -> Inert species are handled as if they had a standard free - * energy of zero and if they obeyed ideal solution/gas theory - * - * @param iphase ID of the phase - * @param w Species mole number vector - * @param fe vector of partial molar free energies of the species. - */ double VCS_SOLVE::vcs_GibbsPhase(size_t iphase, const double* const w, const double* const fe) { diff --git a/src/equil/vcs_TP.cpp b/src/equil/vcs_TP.cpp index c7be748b1..9e1b74d19 100644 --- a/src/equil/vcs_TP.cpp +++ b/src/equil/vcs_TP.cpp @@ -9,38 +9,7 @@ namespace VCSnonideal { - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg) - -/************************************************************************** - * - * vcs_TP: - * - * Solve an equilibrium problem at a particular fixed temperature - * and pressure - * - * ipr = 1 -> Print results to standard output - * 0 -> don't report on anything - * ip1 = 1 -> Print intermediate results. - * maxit -> Maximum number of iterations for the algorithm - * T = Temperature (Kelvin) - * pres = Pressure (pascal) - * - * Return Codes - * ------------------ - * 0 = Equilibrium Achieved - * 1 = Range space error encountered. The element abundance criteria are - * only partially satisfied. Specifically, the first NC= (number of - * components) conditions are satisfied. However, the full NE - * (number of elements) conditions are not satisfied. The equilibrium - * condition is returned. - * -1 = Maximum number of iterations is exceeded. Convergence was not - * found. - ***************************************************************************/ { int retn, iconv; /* @@ -91,26 +60,8 @@ int VCS_SOLVE::vcs_TP(int ipr, int ip1, int maxit, double T_arg, double pres_arg */ return iconv; } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ -/*ARGSUSED*/ -int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) -/************************************************************************** - * - * vcs_evalSS_TP: - * - * IPR = 1 -> Print results to standard output - * 0 -> don't report on anything - * IP1 = 1 -> Print intermediate results. - * T = Temperature (Kelvin) - * Pres = Pressure (Pascal) - * - * Evaluate the standard state free energies at the current temperature - * and pressure. Ideal gas pressure contribution is added in here. - * - ***************************************************************************/ +int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) { // int i; //double R; @@ -156,18 +107,9 @@ int VCS_SOLVE::vcs_evalSS_TP(int ipr, int ip1, double Temp, double pres) } } return VCS_SUCCESS; -} /***************************************************************************/ - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ +} void VCS_SOLVE::vcs_fePrep_TP(void) - -/************************************************************************** - * - * - ***************************************************************************/ { for (size_t i = 0; i < m_numSpeciesTot; ++i) { /* @@ -180,7 +122,6 @@ void VCS_SOLVE::vcs_fePrep_TP(void) m_feSpecies_new[i] = m_SSfeSpecies[i]; } } -} /* vcs_fePrep_TP() ********************************************************/ - } +} diff --git a/src/equil/vcs_elem.cpp b/src/equil/vcs_elem.cpp index 155adadad..e2e2522be 100644 --- a/src/equil/vcs_elem.cpp +++ b/src/equil/vcs_elem.cpp @@ -11,13 +11,6 @@ namespace VCSnonideal { - - -//! Computes the current elemental abundances vector -/*! - * Computes the elemental abundances vector, m_elemAbundances[], and stores it - * back into the global structure - */ void VCS_SOLVE::vcs_elab() { for (size_t j = 0; j < m_numElemConstraints; ++j) { @@ -30,36 +23,6 @@ void VCS_SOLVE::vcs_elab() } } - -/* - * - * vcs_elabcheck: - * - * This function checks to see if the element abundances are in - * compliance. If they are, then TRUE is returned. If not, - * FALSE is returned. Note the number of constraints checked is - * usually equal to the number of components in the problem. This - * routine can check satisfaction of all of the constraints in the - * problem, which is equal to ne. However, the solver can't fix - * breakage of constraints above nc, because that nc is the - * range space by definition. Satisfaction of extra constraints would - * have had to occur in the problem specification. - * - * The constraints should be broken up into 2 sections. If - * a constraint involves a formula matrix with positive and - * negative signs, and eaSet = 0.0, then you can't expect that the - * sum will be zero. There may be roundoff that inhibits this. - * However, if the formula matrix is all of one sign, then - * this requires that all species with nonzero entries in the - * formula matrix be identically zero. We put this into - * the logic below. - * - * Input - * ------- - * ibound = 1 : Checks constraints up to the number of elements - * 0 : Checks constraints up to the number of components. - * - */ bool VCS_SOLVE::vcs_elabcheck(int ibound) { size_t top = m_numComponents; @@ -127,23 +90,9 @@ bool VCS_SOLVE::vcs_elabcheck(int ibound) } } return true; -} /* vcs_elabcheck() *********************************************************/ - -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ +} void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase) - -/************************************************************************* - * - * vcs_elabPhase: - * - * Computes the elemental abundances vector for a single phase, - * elemAbundPhase[], and returns it through the argument list. - * The mole numbers of species are taken from the current value - * in m_molNumSpecies_old[]. - *************************************************************************/ { for (size_t j = 0; j < m_numElemConstraints; ++j) { elemAbundPhase[j] = 0.0; @@ -157,52 +106,7 @@ void VCS_SOLVE::vcs_elabPhase(size_t iphase, double* const elemAbundPhase) } } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ - int VCS_SOLVE::vcs_elcorr(double aa[], double x[]) - -/************************************************************************** - * - * vcs_elcorr: - * - * This subroutine corrects for element abundances. At the end of the - * subroutine, the total moles in all phases are recalculated again, - * because we have changed the number of moles in this routine. - * - * Input - * -> temporary work vectors: - * aa[ne*ne] - * x[ne] - * - * Return Values: - * 0 = Nothing of significance happened, - * Element abundances were and still are good. - * 1 = The solution changed significantly; - * The element abundances are now good. - * 2 = The solution changed significantly, - * The element abundances are still bad. - * 3 = The solution changed significantly, - * The element abundances are still bad and a component - * species got zeroed out. - * - * Internal data to be worked on:: - * - * ga Current element abundances - * m_elemAbundancesGoal Required elemental abundances - * m_molNumSpecies_old Current mole number of species. - * m_formulaMatrix[][] Formula matrix of the species - * ne Number of elements - * nc Number of components. - * - * NOTES: - * This routine is turning out to be very problematic. There are - * lots of special cases and problems with zeroing out species. - * - * Still need to check out when we do loops over nc vs. ne. - * - *************************************************************************/ { int retn = 0, its; bool goodSpec; @@ -608,4 +512,3 @@ L_CLEANUP: } } - diff --git a/src/equil/vcs_elem_rearrange.cpp b/src/equil/vcs_elem_rearrange.cpp index 29733da09..1c5309140 100644 --- a/src/equil/vcs_elem_rearrange.cpp +++ b/src/equil/vcs_elem_rearrange.cpp @@ -20,43 +20,6 @@ namespace VCSnonideal { -// Rearrange the constraint equations represented by the Formula -// Matrix so that the operational ones are in the front -/* - * - * This subroutine handles the rearrangement of the constraint - * equations represented by the Formula Matrix. Rearrangement is only - * necessary when the number of components is less than the number of - * elements. For this case, some constraints can never be satisfied - * exactly, because the range space represented by the Formula - * Matrix of the components can't span the extra space. These - * constraints, which are out of the range space of the component - * Formula matrix entries, are migrated to the back of the Formula - * matrix. - * - * A prototypical example is an extra element column in - * FormulaMatrix[], - * which is identically zero. For example, let's say that argon is - * has an element column in FormulaMatrix[], but no species in the - * mechanism - * actually contains argon. Then, nc < ne. Also, without perturbation - * of FormulaMatrix[] vcs_basopt[] would produce a zero pivot - * because the matrix - * would be singular (unless the argon element column was already the - * last column of FormulaMatrix[]. - * This routine borrows heavily from vcs_basopt's algorithm. It - * finds nc constraints which span the range space of the Component - * Formula matrix, and assigns them as the first nc components in the - * formula matrix. This guarantees that vcs_basopt[] has a - * nonsingular matrix to invert. - * - * Other Variables - * aw[i] = Mole fraction work space (ne in length) - * sa[j] = Gramm-Schmidt orthog work space (ne in length) - * ss[j] = Gramm-Schmidt orthog work space (ne in length) - * sm[i+j*ne] = QR matrix work space (ne*ne in length) - * - */ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, double* const sm, double* const ss) { @@ -217,15 +180,6 @@ int VCS_SOLVE::vcs_elem_rearrange(double* const aw, double* const sa, return VCS_SUCCESS; } -// Swaps the indices for all of the global data for two elements, ipos -// and jpos. -/* - * This function knows all of the element information with VCS_SOLVE, and - * can therefore switch element positions - * - * @param ipos first global element index - * @param jpos second global element index - */ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) { if (ipos == jpos) { @@ -268,4 +222,3 @@ void VCS_SOLVE::vcs_switch_elem_pos(size_t ipos, size_t jpos) std::swap(m_elementName[ipos], m_elementName[jpos]); } } - diff --git a/src/equil/vcs_inest.cpp b/src/equil/vcs_inest.cpp index 8f5980274..ba8272a98 100644 --- a/src/equil/vcs_inest.cpp +++ b/src/equil/vcs_inest.cpp @@ -23,19 +23,6 @@ namespace VCSnonideal static char pprefix[20] = " --- vcs_inest: "; -// Estimate equilibrium compositions -/* - * Estimates equilibrium compositions. - * Algorithm covered in a section of Smith and Missen's Book. - * - * Linear programming module is based on using dbolm. - * - * @param aw aw[i[ Mole fraction work space (ne in length) - * @param sa sa[j] = Gramm-Schmidt orthog work space (ne in length) - * @param sm sm[i+j*ne] = QR matrix work space (ne*ne in length) - * @param ss ss[j] = Gramm-Schmidt orthog work space (ne in length) - * @param test This is a small negative number. - */ void VCS_SOLVE::vcs_inest(double* const aw, double* const sa, double* const sm, double* const ss, double test) { @@ -358,16 +345,7 @@ finished: } #endif } -/***************************************************************************/ -// Create an initial estimate of the solution to the thermodynamic -// equilibrium problem. -/* - * @return Return value indicates success: - * - 0: successful initial guess - * - -1: Unsuccessful initial guess; the elemental abundances aren't - * satisfied. - */ int VCS_SOLVE::vcs_inest_TP() { int retn = 0; diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index 72eeae77d..e529b4978 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -20,16 +20,6 @@ namespace VCSnonideal { - -// Returns the multiplier for electric charge terms -/* - * This is basically equal to F/RT - * - * @param mu_units integer representing the dimensional units system - * @param TKelvin double Temperature in Kelvin - * - * @return Returns the value of F/RT - */ double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const { double Farad; @@ -57,15 +47,6 @@ double VCS_SOLVE::vcs_nondim_Farad(int mu_units, double TKelvin) const return Farad; } -// Returns the multiplier for the nondimensionalization of the equations -/* - * This is basically equal to RT - * - * @param mu_units integer representing the dimensional units system - * @param TKelvin double Temperature in Kelvin - * - * @return Returns the value of RT - */ double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const { double rt; @@ -96,22 +77,6 @@ double VCS_SOLVE::vcs_nondimMult_TP(int mu_units, double TKelvin) const return rt; } -// Nondimensionalize the problem data -/* - * Nondimensionalize the free energies using the divisor, R * T - * - * Essentially the internal data can either be in dimensional form - * or in nondimensional form. This routine switches the data from - * dimensional form into nondimensional form. - * - * What we do is to divide by RT. - * - * @todo Add a scale factor based on the total mole numbers. - * The algorithm contains hard coded numbers based on the - * total mole number. If we ever were faced with a problem - * with significantly different total kmol numbers than one - * the algorithm would have problems. - */ void VCS_SOLVE::vcs_nondim_TP() { double tf; @@ -203,16 +168,6 @@ void VCS_SOLVE::vcs_nondim_TP() } } -// Redimensionalize the problem data -/* - * Redimensionalize the free energies using the multiplier R * T - * - * Essentially the internal data can either be in dimensional form - * or in nondimensional form. This routine switches the data from - * nondimensional form into dimensional form. - * - * What we do is to multiply by RT. - */ void VCS_SOLVE::vcs_redim_TP(void) { double tf; @@ -260,11 +215,6 @@ void VCS_SOLVE::vcs_redim_TP(void) } } -// Computes the current elemental abundances vector -/* - * Computes the elemental abundances vector, m_elemAbundances[], and stores it - * back into the global structure - */ void VCS_SOLVE::vcs_printChemPotUnits(int unitsFormat) const { switch (unitsFormat) { @@ -290,4 +240,3 @@ void VCS_SOLVE::vcs_printChemPotUnits(int unitsFormat) const } } - diff --git a/src/equil/vcs_phaseStability.cpp b/src/equil/vcs_phaseStability.cpp index d6d58cc06..2e7158885 100644 --- a/src/equil/vcs_phaseStability.cpp +++ b/src/equil/vcs_phaseStability.cpp @@ -1,7 +1,7 @@ /** * @file vcs_phaseStability.cpp * Implementation class for functions associated with determining the stability of a phase - * (see Class \link Cantera::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). + * (see Class \link VCSnonideal::VCS_SOLVE VCS_SOLVE\endlink and \ref equilfunctions ). */ #include "cantera/equil/vcs_solve.h" #include "cantera/equil/vcs_internal.h" @@ -20,25 +20,8 @@ using namespace std; namespace VCSnonideal { -//==================================================================================================================== -// Utility function that evaluates whether a phase can be popped into existence -/* - * A phase can be popped iff the stoichiometric coefficients for the - * component species, whose concentrations will be lowered during the - * process, are positive by at least a small degree. - * - * If one of the phase species is a zeroed component, then the phase can - * be popped if the component increases in mole number as the phase moles - * are increased. - * - * @param iphasePop id of the phase, which is currently zeroed, - * - * @return Returns true if the phase can come into existence - * and false otherwise. - */ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const { - vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop]; #ifdef DEBUG_MODE @@ -133,16 +116,6 @@ bool VCS_SOLVE::vcs_popPhasePossible(const size_t iphasePop) const return false; } -//==================================================================================================================== -// Determine the list of problems that need to be checked to see if there are any phases pops -/* - * This routine evaluates and fills in the following quantities - * phasePopProblemLists_ - * - * Need to work in species that are zeroed by element constraints - * - * @return Returns the number of problems that must be checked. - */ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() { int nfound = 0; @@ -270,13 +243,6 @@ int VCS_SOLVE::vcs_phasePopDeterminePossibleList() return nfound; } - -//==================================================================================================================== -// Decision as to whether a phase pops back into existence -/* - * @return returns the phase id of the phases that pops back into - * existence. Returns npos if there are no phases - */ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) { size_t iphasePop = npos; @@ -407,29 +373,7 @@ size_t VCS_SOLVE::vcs_popPhaseID(std::vector & phasePopPhaseIDs) #endif return iphasePop; } -//==================================================================================================================== -// Calculates the deltas of the reactions due to phases popping -// into existence -/* - * @param iphasePop Phase id of the phase that will come into existence - * - * Output - * ------- - * m_deltaMolNumSpecies(irxn) : reaction adjustments, where irxn refers - * to the irxn'th species - * formation reaction. This adjustment - * is for species - * irxn + M, where M is the number - * of components. - * - * @return Returns an int representing the status of the step - * - 0 : normal return - * - 1 : A single species phase species has been zeroed out - * in this routine. The species is a noncomponent - * - 2 : Same as one but, the zeroed species is a component. - * - 3 : Nothing was done because the phase couldn't be birthed - * because a needed component is zero. - */ + int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) { vcs_VolPhase* Vphase = m_VolPhaseList[iphasePop]; @@ -626,19 +570,8 @@ int VCS_SOLVE::vcs_popPhaseRxnStepSizes(const size_t iphasePop) return 0; } -//====================================================================================================================== -// Main program to test whether a deleted phase should be brought -// back into existence -/* - * - * @param iph Phase id of the deleted phase - * - * So far this algorithm seems to be stable. I haven't run across any instance where it hasn't - * converged. However, it's probably only a matter of time - */ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) { - /* * We will use the _new state calc here */ @@ -970,4 +903,3 @@ double VCS_SOLVE::vcs_phaseStabilityTest(const size_t iph) } } - diff --git a/src/equil/vcs_prep.cpp b/src/equil/vcs_prep.cpp index 7dc3e9fc1..98dea4a08 100644 --- a/src/equil/vcs_prep.cpp +++ b/src/equil/vcs_prep.cpp @@ -21,9 +21,6 @@ namespace VCSnonideal { - - -// Calculate the status of single species phases. void VCS_SOLVE::vcs_SSPhase() { size_t iph; @@ -68,41 +65,7 @@ void VCS_SOLVE::vcs_SSPhase() } } } -/*****************************************************************************/ -// This routine is mostly concerned with changing the private data -// to be consistent with what's needed for solution. It is called one -// time for each new problem structure definition. -/* - * This routine is always followed by vcs_prep(). Therefore, tasks - * that need to be done for every call to vcsc() should be placed in - * vcs_prep() and not in this routine. - * - * The problem structure refers to: - * - * the number and identity of the species. - * the formula matrix and thus the number of components. - * the number and identity of the phases. - * the equation of state - * the method and parameters for determining the standard state - * The method and parameters for determining the activity coefficients. - * - * Tasks: - * 0) Fill in the SSPhase[] array. - * 1) Check to see if any multispecies phases actually have only one - * species in that phase. If true, reassign that phase and species - * to be a single-species phase. - * 2) Determine the number of components in the problem if not already - * done so. During this process the order of the species is changed - * in the private data structure. All references to the species - * properties must employ the ind[] index vector. - * - * @param printLvl Print level of the routine - * - * @return the return code - * VCS_SUCCESS = everything went OK - * - */ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) { size_t kspec, i; @@ -114,7 +77,6 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) m_debug_print_lvl = printLvl; - /* * Calculate the Single Species status of phases * Also calculate the number of species per phase @@ -258,24 +220,6 @@ int VCS_SOLVE::vcs_prep_oneTime(int printLvl) return VCS_SUCCESS; } -/*****************************************************************************/ - -// Prepare the object for re-solution -/* - * This routine is mostly concerned with changing the private data - * to be consistent with that needed for solution. It is called for - * every invocation of the vcs_solve() except for the cleanup invocation. - * - * Tasks: - * 1) Initialization of arrays to zero. - * 2) Calculate total number of moles in all phases - * - * return code - * VCS_SUCCESS = everything went OK - * VCS_PUB_BAD = There is an irreconcilable difference in the - * public data structure from when the problem was - * initially set up. - */ int VCS_SOLVE::vcs_prep() { /* @@ -295,23 +239,6 @@ int VCS_SOLVE::vcs_prep() return VCS_SUCCESS; } -/*****************************************************************************/ - -// In this routine, we check for things that will cause the algorithm -// to fail. -/* - * We check to see if the problem is well posed. If it is not, we return - * false and print out error conditions. - * - * Current there is one condition. If all the element abundances are - * zero, the algorithm will fail - * - * @param vprob VCS_PROB pointer to the definition of the equilibrium - * problem - * - * @return If true, the problem is well-posed. If false, the problem - * is not well posed. - */ bool VCS_SOLVE::vcs_wellPosed(VCS_PROB* vprob) { double sum = 0.0; @@ -325,5 +252,4 @@ bool VCS_SOLVE::vcs_wellPosed(VCS_PROB* vprob) return true; } -/*****************************************************************************/ } diff --git a/src/equil/vcs_rank.cpp b/src/equil/vcs_rank.cpp index 8a292e9b3..aef6c2ec2 100644 --- a/src/equil/vcs_rank.cpp +++ b/src/equil/vcs_rank.cpp @@ -1,18 +1,14 @@ - /*! * @file vcs_rank.cpp * Header file for the internal class that holds the problem. */ + /* - * $Id: vcs_solve.cpp 735 2011-07-25 14:44:41Z hkmoffa $ - */ -/* - * Copywrite (2005) Sandia Corporation. Under the terms of + * Copyright (2005) Sandia Corporation. Under the terms of * Contract DE-AC04-94AL85000 with Sandia Corporation, the * U.S. Government retains certain rights in this software. */ - #include "cantera/equil/vcs_solve.h" #include "cantera/equil/vcs_internal.h" #include "cantera/equil/vcs_prob.h" @@ -30,7 +26,6 @@ using namespace std; namespace VCSnonideal { - //==================================================================================================================== static int basisOptMax1(const double * const molNum, const int n) { // int largest = 0; @@ -49,45 +44,7 @@ namespace VCSnonideal { } return n-1; } - //==================================================================================================================== - // Calculate the rank of a matrix and return the rows and columns that will generate an independent basis - // for that rank - /* - * Choose the optimum component species basis for the calculations, finding the rank and - * set of linearly independent rows for that calculation. - * Then find the set of linearly indepedent element columns that can support that rank. - * This is done by taking the transpose of the matrix and redoing the same calculation. - * (there may be a better way to do this. I don't know.) - * - * - * Input - * --------- - * - * @param awtmp Vector of mole numbers which will be used to construct a - * ranking for how to pick the basis species. This is largely ignored - * here. - * - * @param numSpecies Number of species. This is the number of rows in the matrix. - * - * @param matrix Matrix. This is the formula matrix. Nominally, the rows are species, while - * the columns are element compositions. However, this routine - * is totally general, so that the rows and columns can be anything. - * - * @param numElemConstraints Number of element constraints - * - * Output - * --------- - * @param usedZeroedSpecies = If true, then a species with a zero concentration - * was used as a component. - * - * - * @param compRes Vector of rows which are linearly independent. (these are the components) - * - * @param elemComp Vector of columns which are linearly independent (These are the actionable element - * constraints). - * - * @return Returns number of components. This is the rank of the matrix - */ + int VCS_SOLVE::vcs_rank(const double * awtmp, size_t numSpecies, const double matrix[], size_t numElemConstraints, std::vector &compRes, std::vector& elemComp, int * const usedZeroedSpecies) const { @@ -343,6 +300,5 @@ namespace VCSnonideal { } return numComponents; } - } diff --git a/src/equil/vcs_rearrange.cpp b/src/equil/vcs_rearrange.cpp index 18777d1f1..8bfe87db7 100644 --- a/src/equil/vcs_rearrange.cpp +++ b/src/equil/vcs_rearrange.cpp @@ -15,15 +15,8 @@ #include #include - namespace VCSnonideal { - - -// Switch all species data back to the original order. -/* - * This destroys the data based on reaction ordering. - */ int VCS_SOLVE::vcs_rearrange() { size_t i, l, j; @@ -53,4 +46,3 @@ int VCS_SOLVE::vcs_rearrange() } } - diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index f23d6afd9..80b38ec3d 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -1,11 +1,9 @@ - /* * Copyright (2005) Sandia Corporation. Under the terms of * Contract DE-AC04-94AL85000 with Sandia Corporation, the * U.S. Government retains certain rights in this software. */ - #include "cantera/equil/vcs_solve.h" #include "cantera/equil/vcs_internal.h" #include "cantera/equil/vcs_VolPhase.h" @@ -16,8 +14,6 @@ namespace VCSnonideal { - -/*****************************************************************************/ static void print_space(int num) { for (int j = 0; j < num; j++) { @@ -33,17 +29,6 @@ static void print_line(const std::string& schar, size_t num) plogf("\n"); } -/*****************************************************************************/ - -/************************************************************************** - * - * vcs_report: - * - * Print out a report on the state of the equilibrium problem to - * standard output. - * This prints out the current contents of the VCS_SOLVE class, V. - * The "old" solution vector is printed out. - ***************************************************************************/ int VCS_SOLVE::vcs_report(int iconv) { bool printActualMoles = true, inertYes = false; @@ -409,19 +394,9 @@ int VCS_SOLVE::vcs_report(int iconv) * Return a successful completion flag */ return VCS_SUCCESS; -} /* vcs_report() ************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ +} void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl) - -/************************************************************************** - * - * vcs_TCounters_report: - * - * Print out the total Its and time counters to standard output - ***************************************************************************/ { plogf("\nTCounters: Num_Calls Total_Its Total_Time (seconds)\n"); if (timing_print_lvl > 0) { @@ -447,8 +422,4 @@ void VCS_SOLVE::vcs_TCounters_report(int timing_print_lvl) } } -/*****************************************************************************/ -/*****************************************************************************/ -/*****************************************************************************/ } - diff --git a/src/equil/vcs_rxnadj.cpp b/src/equil/vcs_rxnadj.cpp index de659a1c3..895020482 100644 --- a/src/equil/vcs_rxnadj.cpp +++ b/src/equil/vcs_rxnadj.cpp @@ -19,28 +19,6 @@ namespace VCSnonideal { -// Calculates formation reaction step sizes. -/* - * This is equation 6.4-16, p. 143 in Smith and Missen. - * - * Output - * ------- - * m_deltaMolNumSpecies[kspec] : reaction adjustments, where irxn refers - * to the irxn'th species - * formation reaction. This adjustment - * is for species - * irxn + M, where M is the number - * of components. - * - * Special branching occurs sometimes. This causes the component basis - * to be reevaluated - * - * @return Returns an int representing the status of the step - * - 0 : normal return - * - 1 : A single species phase species has been zeroed out - * in this routine. The species is a noncomponent - * - 2 : Same as one but, the zeroed species is a component. - */ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) { size_t kspec, iph; @@ -404,32 +382,6 @@ size_t VCS_SOLVE::vcs_RxnStepSizes(int& forceComponentCalc, size_t& kSpecial) return iphDel; } -//==================================================================================================================== -// Calculates reaction adjustments using a full Hessian approximation -/* - * Calculates reaction adjustments. This does what equation 6.4-16, p. 143 - * in Smith and Missen is suppose to do. However, a full matrix is - * formed and then solved via a conjugate gradient algorithm. No - * preconditioning is done. - * - * If special branching is warranted, then the program bails out. - * - * Output - * ------- - * DS(I) : reaction adjustment, where I refers to the Ith species - * Special branching occurs sometimes. This causes the component basis - * to be reevaluated - * return = 0 : normal return - * 1 : A single species phase species has been zeroed out - * in this routine. The species is a noncomponent - * 2 : Same as one but, the zeroed species is a component. - * - * Special attention is taken to flag cases where the direction of the - * update is contrary to the steepest descent rule. This is an important - * attribute of the regular vcs algorithm. We don't want to violate this. - * - * NOTE: currently this routine is not used. - */ int VCS_SOLVE::vcs_rxn_adj_cg() { size_t irxn, j; @@ -640,18 +592,6 @@ int VCS_SOLVE::vcs_rxn_adj_cg() return soldel; } -//==================================================================================================================== -// Calculates the diagonal contribution to the Hessian due to -// the dependence of the activity coefficients on the mole numbers. -/* - * (See framemaker notes, Eqn. 20 - VCS Equations document) - * - * We allow the diagonal to be increased positively to any degree. - * We allow the diagonal to be decreased to 1/3 of the ideal solution - * value, but no more -> it must remain positive. - * - * NOTE: currently this routine is not used - */ double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal) { double diag = hessianDiag_Ideal; @@ -670,14 +610,6 @@ double VCS_SOLVE::vcs_Hessian_diag_adj(size_t irxn, double hessianDiag_Ideal) return diag; } -//==================================================================================================================== -//! Calculates the diagonal contribution to the Hessian due to -//! the dependence of the activity coefficients on the mole numbers. -/*! - * (See framemaker notes, Eqn. 20 - VCS Equations document) - * - * NOTE: currently this routine is not used - */ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) { size_t kspec, k, l, kph; @@ -716,15 +648,7 @@ double VCS_SOLVE::vcs_Hessian_actCoeff_diag(size_t irxn) } return s; } -//==================================================================================================================== -// Recalculate all of the activity coefficients in all of the phases -// based on input mole numbers -/* - * - * @param moleSpeciesVCS kmol of species to be used in the update. - * - * NOTE: This routine needs to be regulated. - */ + void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS) { /* @@ -749,41 +673,7 @@ void VCS_SOLVE::vcs_CalcLnActCoeffJac(const double* const moleSpeciesVCS) } } } -/*****************************************************************************/ -//! This function recalculates the deltaG for reaction, irxn -/*! - * This function recalculates the deltaG for reaction irxn, - * given the mole numbers in molNum. It uses the temporary - * space mu_i, to hold the recalculated chemical potentials. - * It only recalculates the chemical potentials for species in phases - * which participate in the irxn reaction. - * - * This function is used by the vcs_line_search algorithm() and - * should not be used widely due to the unknown state it leaves the - * system. - * - * Input - * ------------ - * @param irxn Reaction number - * @param molNum Current mole numbers of species to be used as - * input to the calculation (units = kmol) - * (length = totalNuMSpecies) - * - * Output - * ------------ - * @param ac output Activity coefficients (length = totalNumSpecies) - * Note this is only partially formed. Only species in - * phases that participate in the reaction will be updated - * @param mu_i dimensionless chemical potentials (length - totalNumSpecies - * Note this is only partially formed. Only species in - * phases that participate in the reaction will be updated - * - * @return Returns the dimensionless deltaG of the reaction - * - * Note, this is a dangerous routine that leaves the underlying objects in - * an unknown state. - */ double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, const size_t irxn, const double* const molNum, double* const ac, double* const mu_i) { @@ -801,22 +691,8 @@ double VCS_SOLVE::deltaG_Recalc_Rxn(const int stateCalc, const size_t irxn, cons } return deltaG; } -/*****************************************************************************/ #ifdef DEBUG_MODE -// A line search algorithm is carried out on one reaction -/* - * In this routine we carry out a rough line search algorithm - * to make sure that the m_deltaGRxn_new doesn't switch signs prematurely. - * - * @param irxn Reaction number - * @param dx_orig Original step length - * - * @param ANOTE Output character string stating the conclusions of the - * line search - * - * @return Returns the optimized step length found by the search - */ double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig, char* const ANOTE) #else double VCS_SOLVE::vcs_line_search(const size_t irxn, const double dx_orig) @@ -955,6 +831,5 @@ finalize: return dx; } -/*****************************************************************************/ -} +} diff --git a/src/equil/vcs_setMolesLinProg.cpp b/src/equil/vcs_setMolesLinProg.cpp index 193f82282..b0e6e5948 100644 --- a/src/equil/vcs_setMolesLinProg.cpp +++ b/src/equil/vcs_setMolesLinProg.cpp @@ -1,6 +1,5 @@ /*! * @file vcs_setMolesLinProg.cpp - * */ /* * Copyright (2005) Sandia Corporation. Under the terms of @@ -45,18 +44,6 @@ static void printProgress(const vector &spName, } #endif -//! Estimate the initial mole numbers. -/*! - * This is done by running - * each reaction as far forward or backward as possible, subject - * to the constraint that all mole numbers remain - * non-negative. Reactions for which \f$ \Delta \mu^0 \f$ are - * positive are run in reverse, and ones for which it is negative - * are run in the forward direction. The end result is equivalent - * to solving the linear programming problem of minimizing the - * linear Gibbs function subject to the element and - * non-negativity constraints. - */ int VCS_SOLVE::vcs_setMolesLinProg() { size_t ik, irxn; diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 6da265b49..529d16b8d 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -1,6 +1,6 @@ /*! - * @file vcs_solve.h - * Header file for the internal class that holds the problem. + * @file vcs_solve.h Header file for the internal class that holds the + * problem. */ /* * Copyright (2005) Sandia Corporation. Under the terms of @@ -8,7 +8,6 @@ * U.S. Government retains certain rights in this software. */ - #include "cantera/equil/vcs_solve.h" #include "cantera/base/ctexceptions.h" #include "cantera/equil/vcs_internal.h" @@ -62,27 +61,9 @@ VCS_SOLVE::VCS_SOLVE() : { } -// Initialize the sizes within the VCS_SOLVE object - -/* - * This resizes all of the internal arrays within the object. This routine - * operates in two modes. If all of the parameters are the same as it - * currently exists in the object, nothing is done by this routine; a quick - * exit is carried out and all of the data in the object persists. - * - * IF any of the parameters are different than currently exists in the - * object, then all of the data in the object must be redone. It may not - * be zeroed, but it must be redone. - * - * @param nspecies0 Number of species within the object - * @param nelements Number of element constraints within the problem - * @param nphase0 Number of phases defined within the problem. - * - */ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, const size_t nphase0) { - if (NSPECIES0 != 0) { if ((nspecies0 != NSPECIES0) || (nelements != m_numElemConstraints) || (nphase0 != NPHASE0)) { vcs_delete_memory(); @@ -237,19 +218,12 @@ void VCS_SOLVE::vcs_initSizes(const size_t nspecies0, const size_t nelements, return; } -/****************************************************************************/ -// Destructor VCS_SOLVE::~VCS_SOLVE() { vcs_delete_memory(); } -/*****************************************************************************/ -// Delete memory that isn't just resizeable STL containers -/* - * This gets called by the destructor or by InitSizes(). - */ void VCS_SOLVE::vcs_delete_memory() { size_t j; @@ -273,44 +247,7 @@ void VCS_SOLVE::vcs_delete_memory() m_numElemConstraints = 0; m_numComponents = 0; } -/*****************************************************************************/ -// Solve an equilibrium problem -/* - * This is the main interface routine to the equilibrium solver - * - * Input: - * @param vprob Object containing the equilibrium Problem statement - * - * @param ifunc Determines the operation to be done: Valid values: - * 0 -> Solve a new problem by initializing structures - * first. An initial estimate may or may not have - * been already determined. This is indicated in the - * VCS_PROB structure. - * 1 -> The problem has already been initialized and - * set up. We call this routine to resolve it - * using the problem statement and - * solution estimate contained in - * the VCS_PROB structure. - * 2 -> Don't solve a problem. Destroy all the private - * structures. - * - * @param ipr Printing of results - * ipr = 1 -> Print problem statement and final results to - * standard output - * 0 -> don't report on anything - * @param ip1 Printing of intermediate results - * ip1 = 1 -> Print intermediate results. - * = 0 -> No intermediate results printing - * - * @param maxit Maximum number of iterations for the algorithm - * - * Output: - * - * @return - * nonzero value: failure to solve the problem at hand. - * zero : success - */ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) { int retn = 0, iconv = 0; @@ -444,16 +381,7 @@ int VCS_SOLVE::vcs(VCS_PROB* vprob, int ifunc, int ipr, int ip1, int maxit) } return iconv; } -/*****************************************************************************/ -// Fully specify the problem to be solved using VCS_PROB -/* - * Use the contents of the VCS_PROB to specify the contents of the - * private data, VCS_SOLVE. - * - * @param pub Pointer to VCS_PROB that will be used to - * initialize the current equilibrium problem - */ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) { vcs_VolPhase* Vphase = 0; @@ -835,18 +763,7 @@ int VCS_SOLVE::vcs_prob_specifyFully(const VCS_PROB* pub) */ return VCS_SUCCESS; } -/*****************************************************************************/ -// Specify the problem to be solved using VCS_PROB, incrementally -/* - * Use the contents of the VCS_PROB to specify the contents of the - * private data, VCS_SOLVE. - * - * It's assumed we are solving the same problem. - * - * @param pub Pointer to VCS_PROdB that will be used to - * initialize the current equilibrium problem - */ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub) { size_t kspec, k, i, j, iph; @@ -971,15 +888,7 @@ int VCS_SOLVE::vcs_prob_specify(const VCS_PROB* pub) return retn; } -/*****************************************************************************/ -// Transfer the results of the equilibrium calculation back to VCS_PROB -/* - * The VCS_PUB structure is returned to the user. - * - * @param pub Pointer to VCS_PROB that will get the results of the - * equilibrium calculation transfered to it. - */ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) { size_t k1 = 0; @@ -1063,17 +972,7 @@ int VCS_SOLVE::vcs_prob_update(VCS_PROB* pub) return VCS_SUCCESS; } -/*****************************************************************************/ -// Initialize the internal counters -/* - * Initialize the internal counters containing the subroutine call - * values and times spent in the subroutines. - * - * ifunc = 0 Initialize only those counters appropriate for the top of - * vcs_solve_TP(). - * = 1 Initialize all counters. - */ void VCS_SOLVE::vcs_counters_init(int ifunc) { m_VCount->Its = 0; @@ -1091,30 +990,7 @@ void VCS_SOLVE::vcs_counters_init(int ifunc) m_VCount->T_Time_vcs = 0.0; } } -/**************************************************************************/ -// Calculation of the total volume and the partial molar volumes -/* - * This function calculates the partial molar volume - * for all species, kspec, in the thermo problem - * at the temperature TKelvin and pressure, Pres, pres is in atm. - * And, it calculates the total volume of the combined system. - * - * Input - * --------------- - * @param tkelvin Temperature in kelvin() - * @param pres Pressure in Pascal - * @param w w[] is thevector containing the current mole numbers - * in units of kmol. - * - * Output - * ---------------- - * @param volPM[] For species in all phase, the entries are the - * partial molar volumes units of M**3 / kmol. - * - * @return The return value is the total volume of - * the entire system in units of m**3. - */ double VCS_SOLVE::vcs_VolTotal(const double tkelvin, const double pres, const double w[], double volPM[]) { @@ -1128,7 +1004,5 @@ double VCS_SOLVE::vcs_VolTotal(const double tkelvin, const double pres, } return VolTot; } -/****************************************************************************/ - } diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 6126e7b1c..d44af7889 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -32,15 +32,9 @@ using namespace Cantera; namespace VCSnonideal { - -/***************************************************************************/ /************ Prototypes for static functions ******************************/ - - static void print_space(size_t num); - - #ifdef DEBUG_MODE # ifdef DEBUG_NOT static void prneav(void); @@ -49,7 +43,6 @@ static int prnfm(void); #endif /*****************************************************************************/ - #ifdef DEBUG_MODE void VCS_SOLVE::checkDelta1(double* const dsLocal, double* const delTPhMoles, int kspec) @@ -70,36 +63,7 @@ void VCS_SOLVE::checkDelta1(double* const dsLocal, } } #endif -/*****************************************************************************/ -// Main routine that solves for equilibrium at constant T and P -// using a variant of the VCS method -/* - * This is the main routine that solves for equilibrium at constant T and P - * using a variant of the VCS method. Nonideal phases can be accommodated - * as well. - * - * Any number of single-species phases and multi-species phases - * can be handled by the present version. - * - * Input - * ------------ - * @param print_lvl 1 -> Print results to standard output - * 0 -> don't report on anything - * - * @param printDetails 1 -> Print intermediate results. - * - * @param maxit Maximum number of iterations for the algorithm - * - * @return 0 = Equilibrium Achieved - * 1 = Range space error encountered. The element abundance criteria are - * only partially satisfied. Specifically, the first NC= (number of - * components) conditions are satisfied. However, the full NE - * (number of elements) conditions are not satisfied. The equilibrium - * condition is returned. - * -1 = Maximum number of iterations is exceeded. Convergence was not - * found. - */ int VCS_SOLVE::vcs_solve_TP(int print_lvl, int printDetails, int maxit) { int retn = VCS_SUCCESS, soldel, solveFail; @@ -2031,7 +1995,6 @@ L_RETURN_BLOCK_B: */ return solveFail; } -/*********************************************************************************/ double VCS_SOLVE::vcs_minor_alt_calc(size_t kspec, size_t irxn, bool* do_delete #ifdef DEBUG_MODE @@ -2151,21 +2114,7 @@ L_ZERO_SPECIES: } return dx; } -/*****************************************************************************/ -// Change the concentration of a species by delta moles. -/* - * Make sure to conserve elements and keep track of the total kmoles in all phases. - * - * - * @param kspec The species index - * @param delta_ptr pointer to the delta for the species. This may change during - * the calculation - * - * @return - * 1: succeeded without change of dx - * 0: Had to adjust dx, perhaps to zero, in order to do the delta. - */ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) { size_t irxn = kspec - m_numComponents; @@ -2229,21 +2178,7 @@ int VCS_SOLVE::delta_species(const size_t kspec, double* const delta_ptr) } return retn; } -/*****************************************************************************/ -// Zero out the concentration of a species. -/* - * Zero out the concentration of a species. Make sure to conserve - * elements and keep track of the total moles in all phases. - * w[] - * m_tPhaseMoles_old[] - * - * @param kspec Species index - * - * @return: - * 1: succeeded - * 0: failed. - */ int VCS_SOLVE::vcs_zero_species(const size_t kspec) { int retn = 1; @@ -2268,22 +2203,7 @@ int VCS_SOLVE::vcs_zero_species(const size_t kspec) } return retn; } -/**************************************************************************/ -// Change a single species from active to inactive status -/* - * Rearrange data when species is added or removed. The kspec species is - * moved to the back of the active species vector. The back of the species - * vector is indicated by the value of m_numSpeciesRdc, the current - * number of active species in the mechanism. - * - * @param kspec Species Index - * @return - * Returns 0 unless. - * The return is 1 when the current number of - * noncomponent species is equal to zero. A recheck of deleted species - * is carried out in the main code. - */ int VCS_SOLVE::vcs_delete_species(const size_t kspec) { const size_t klast = m_numSpeciesRdc - 1; @@ -2361,25 +2281,7 @@ int VCS_SOLVE::vcs_delete_species(const size_t kspec) */ return (m_numRxnRdc == 0); } -/***************************************************************************/ -/* - * - * reinsert_deleted(): - * - * irxn = id of the noncomponent species formation reaction for the - * species to be added in. - * - * We make decisions on the initial mole number, and major-minor status - * here. We also fix up the total moles in a phase. - * - * The algorithm proceeds to implement these decisions in the previous - * position of the species. Then, vcs_switch_pos is called to move the - * species into the last active species slot, incrementing the number - * of active species at the same time. - * - * This routine is responsible for the global data manipulation only. - */ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) { size_t k; @@ -2444,22 +2346,7 @@ void VCS_SOLVE::vcs_reinsert_deleted(size_t kspec) vcs_switch_pos(true, (m_numSpeciesRdc - 1), kspec); } } -/****************************************************************************/ -// This routine handles the bookkeeping involved with the -// deletion of multiphase phases from the problem. -/* - * When they are deleted, all of their species become active - * species, even though their mole numbers are set to zero. - * This routine does not make the decision to eliminate multiphases. - * - * Note, species in phases with zero mole numbers are still - * considered active. Whether the phase pops back into - * existence or not is checked as part of the main iteration - * loop. - * - * @param iph Phase to be deleted - */ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) { size_t kspec, irxn; @@ -2625,25 +2512,7 @@ bool VCS_SOLVE::vcs_delete_multiphase(const size_t iph) return successful; } -/**********************************************************************************/ -// Recheck deleted species in multispecies phases. -/* - * We are checking the equation: - * - * sum_u = sum_j_comp [ sigma_i_j * u_j ] - * = u_i_O + log((AC_i * W_i)/m_tPhaseMoles_old) - * - * by first evaluating: - * - * DG_i_O = u_i_O - sum_u. - * - * Then, if TL is zero, the phase pops into existence if DG_i_O < 0. - * Also, if the phase exists, then we check to see if the species - * can have a mole number larger than VCS_DELETE_SPECIES_CUTOFF - * (default value = 1.0E-32). - * - */ int VCS_SOLVE::vcs_recheck_deleted() { int npb; @@ -2730,7 +2599,6 @@ int VCS_SOLVE::vcs_recheck_deleted() } return npb; } -/***********************************************************************************/ bool VCS_SOLVE::recheck_deleted_phase(const int iphase) { @@ -2771,13 +2639,7 @@ bool VCS_SOLVE::recheck_deleted_phase(const int iphase) } return false; } -//==================================================================================================================== -// Provide an estimate for the deleted species in phases that are not zeroed out -/* - * Try to add back in all deleted species. An estimate of the kmol numbers - * are obtained and the species is added back into the equation system, - * into the old state vector. - */ + size_t VCS_SOLVE::vcs_add_all_deleted() { size_t iph, kspec, retn; @@ -2890,28 +2752,7 @@ size_t VCS_SOLVE::vcs_add_all_deleted() } return retn; } -/***********************************************************************************/ -// This routine optimizes the minimization of the total gibbs free -// energy by making sure the slope of the Gibbs free energy stays negative -/* - * The slope of the following functional is equivalent to the slope of the total - * Gibbs free energy of the system: - * - * d_Gibbs/ds = sum_k( m_deltaGRxn * m_deltaMolNumSpecies[k] ) - * - * along the current direction m_deltaMolNumSpecies[], by choosing a value, al: (0 0), - * does this code section kick in. It finds the point on the parabola - * where the slope is equal to zero. - * - */ bool VCS_SOLVE::vcs_globStepDamp() { double s1, s2, al; @@ -3054,7 +2895,6 @@ bool VCS_SOLVE::vcs_globStepDamp() #endif return true; } -/****************************************************************************************/ int VCS_SOLVE::vcs_basopt(const bool doJustComponents, double aw[], double sa[], double sm[], double ss[], double test, bool* const usedZeroedSpecies) @@ -3635,36 +3475,25 @@ L_CLEANUP: (m_VCount->Basis_Opts)++; return VCS_SUCCESS; } -/***************************************************************************************/ -// Choose a species to test for the next component -/* - * We make the choice based on testing (molNum[i] * spSize[i]) for its maximum value. - * Preference for single species phases is also made. - * - * The factors of 1.01 and 1.001 are placed in this routine for a purpose. - * The purpose is to ensure that roundoff errors don't influence major decisions. - * This means that the optimized and non-optimized versions of the code remain - * close to each other. - * - * ( we try to avoid the logic: a = b - * if (a > b) { do this } - * else { do something else } - * because roundoff error makes a difference in the inequality evaluation) - * - * Mole numbers are frequently equal to each other in equilibrium problems - * due to constraints. Swaps are only done if there are a 1% difference in the - * mole numbers. Of course this logic isn't foolproof. - * - * @param molNum Mole number vector - * @param j index into molNum[] that indicates where the search will start from - * Previous successful components are swapped into the front of - * molNum[]. - * @param n Length of molNum[] - */ size_t VCS_SOLVE::vcs_basisOptMax(const double* const molNum, const size_t j, const size_t n) { + /* + * The factors of 1.01 and 1.001 are placed in this routine for a purpose. + * The purpose is to ensure that roundoff errors don't influence major + * decisions. This means that the optimized and non-optimized versions of + * the code remain close to each other. + * + * (we try to avoid the logic: a = b + * if (a > b) { do this } + * else { do something else } + * because roundoff error makes a difference in the inequality evaluation) + * + * Mole numbers are frequently equal to each other in equilibrium problems + * due to constraints. Swaps are only done if there are a 1% difference in + * the mole numbers. Of course this logic isn't foolproof. + */ size_t largest = j; double big = molNum[j] * m_spSize[j] * 1.01; if (m_spSize[j] <= 0.0) { @@ -3701,16 +3530,7 @@ size_t VCS_SOLVE::vcs_basisOptMax(const double* const molNum, const size_t j, } return largest; } -/**********************************************************************************/ -// Evaluate the species category for the indicated species -/* - * All evaluations are done using the "old" version of the solution. - * - * @param kspec Species to be evaluated - * - * @return Returns the calculated species type - */ int VCS_SOLVE::vcs_species_type(const size_t kspec) const { @@ -3909,103 +3729,7 @@ int VCS_SOLVE::vcs_species_type(const size_t kspec) const } return VCS_SPECIES_MINOR; } -/*****************************************************************************/ -//! We calculate the dimensionless chemical potentials of all species -//! in a single phase. -/*! - * - * We calculate the dimensionless chemical potentials of all species - * in a single phase. - * - * Note, for multispecies phases which are currently zeroed out, - * the chemical potential is filled out with the standard chemical - * potential. - * - * For species in multispecies phases whose concentration is zero, - * we need to set the mole fraction to a very low value. - * Its chemical potential - * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration - * to keep numbers positive. - * - * Formula: - * --------------- - * - * Ideal Mixtures: - * - * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * - * ( This is equivalent to the adding the log of the - * mole fraction onto the standard chemical - * potential. ) - * - * Non-Ideal Mixtures: - * ActivityConvention = 0: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * ( This is equivalent to the adding the log of the - * mole fraction multiplied by the activity coefficient - * onto the standard chemical potential. ) - * - * ActivityConvention = 1: -> molality activity formulation - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * - ln(Mnaught * m_units) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * note: m_SSfeSpecies(I) is the molality based standard state. - * However, ActCoeff[I] is the molar based activity coefficient - * We have used the formulas; - * - * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] - * where Xmol[N] is the mole fraction of the solvent - * ActCoeff_M[I] is the molality based act coeff. - * - * note: This is equivalent to the "normal" molality formulation: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff_M[I] * m(I)) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase] - * where m[I] is the molality of the ith solute - * - * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) - * - * - * note: z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction - * of i in the phase. - * - * NOTE: - * As per the discussion in vcs_dfe(), for small species where the mole - * fraction is small: - * - * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF - * - * The chemical potential is calculated as: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) - * - * Input - * -------- - * @param iph Phase to be calculated - * @param molNum molNum[i] is the number of moles of species i - * (VCS species order) - * @param do_deleted Do species that are deleted (default = false) - * - * Output - * ----------- - * @param ac Activity coefficients for species in phase - * (VCS species order) - * @param mu_i Dimensionless chemical potentials for phase species - * (VCS species order) - * - */ void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc, const size_t iph, const double* const molNum, double* const ac, double* const mu_i, @@ -4071,136 +3795,7 @@ void VCS_SOLVE::vcs_chemPotPhase(const int stateCalc, } } } -/*********************************************************************************/ -// Calculate the dimensionless chemical potentials of all species or -// of certain groups of species, at a fixed temperature and pressure. -/* - * We calculate the dimensionless chemical potentials of all species - * or certain groups of species here, at a fixed temperature and pressure, - * for the input mole vector z[] in the parameter list. - * Nondimensionalization is achieved by division by RT. - * - * Note, for multispecies phases which are currently zeroed out, - * the chemical potential is filled out with the standard chemical - * potential. - * - * For species in multispecies phases whose concentration is zero, - * we need to set the mole fraction to a very low value. - * Its chemical potential - * is then calculated using the VCS_DELETE_MINORSPECIES_CUTOFF concentration - * to keep numbers positive. - * - * - * Formula: - * --------------- - * - * Ideal Mixtures: - * - * m_feSpecies(I) = m_SSfeSpecies(I) + ln(z(I)) - ln(m_tPhaseMoles[iph]) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * ( This is equivalent to the adding the log of the - * mole fraction onto the standard chemical - * potential. ) - * - * Non-Ideal Mixtures: -> molar activity formulation - * ActivityConvention = 0: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * ( This is equivalent to the adding the log of the - * mole fraction multiplied by the activity coefficient - * onto the standard chemical potential. ) - * - * note: z(I)/tPhMoles_ptr[iph] = Xmol[i] is the mole fraction - * of i in the phase. - * - * ActivityConvention = 1: -> molality activity formulation - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff[I] * z(I)) - ln(m_tPhaseMoles[iph]) - * - ln(Mnaught * m_units) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase]; - * - * note: m_SSfeSpecies(I) is the molality based standard state. - * However, ActCoeff[I] is the molar based activity coefficient - * We have used the formulas; - * - * ActCoeff_M[I] = ActCoeff[I] / Xmol[N] - * where Xmol[N] is the mole fraction of the solvent - * ActCoeff_M[I] is the molality based act coeff. - * - * note: This is equivalent to the "normal" molality formulation below: - * - * m_feSpecies(I) = m_SSfeSpecies(I) - * + ln(ActCoeff_M[I] * m(I)) - * + m_chargeSpecies[I] * Faraday_dim * m_phasePhi[iphase] - * where m[I] is the molality of the ith solute - * - * - * m[I] = Xmol[I] / ( Xmol[N] * Mnaught * m_units) - * - * - * Handling of Small Species: - * ------------------------------ - * As per the discussion above, for small species where the mole - * fraction - * - * z(i) < VCS_DELETE_MINORSPECIES_CUTOFF - * - * The chemical potential is calculated as: - * - * m_feSpecies(I) = m_SSfeSpecies(I) + ln(ActCoeff[i](VCS_DELETE_MINORSPECIES_CUTOFF)) - * - * Handling of "Species" Representing Interfacial Voltages - * --------------------------------------------------------- - * - * These species have species types of VCS_SPECIES_TYPE_INTERFACIALVOLTAGE - * The chemical potentials for these "species" refer to electrons in - * metal electrodes. They have the following formula - * - * m_feSpecies(I) = m_SSfeSpecies(I) - F z[I] / RT - * - * F is Faraday's constant. - * R = gas constant - * T = temperature - * V = potential of the interface = phi_electrode - phi_solution - * - * For these species, the solution vector unknown, z[I], is V, the phase voltage, - * in volts. - * - * Input - * -------- - * @param ll Determine which group of species gets updated - * ll = 0: Calculate for all species - * < 0: calculate for components and for major non-components - * 1: calculate for components and for minor non-components - * - * @param lbot Restricts the calculation of the chemical potential - * to the species between LBOT <= i < LTOP. Usually - * LBOT and LTOP will be equal to 0 and MR, respectively. - * @param ltop Top value of the loops - * - * @param molNum molNum[i] : Number of moles of species i - * -> This can either be the old solution vector - * or the new solution vector depending upon the - * stateCalc value - * - * @param stateCalc Determines whether z is old or new or tmp: - * VCS_STATECALC_NEW: Use the tentative values for the total number of - * moles in the phases, i.e., use TG1 instead of TG etc. - * VCS_STATECALC_OLD: Use the base values of the total number of - * moles in each system. - * - * Also needed: - * m_SSfeSpecies[kspec] : standard state chemical potentials. These are the - * chemical potentials of the standard states at - * the same T and P as the solution. - * tg : Total Number of moles in the phase. - */ void VCS_SOLVE::vcs_dfe(const int stateCalc, const int ll, const size_t lbot, const size_t ltop) { @@ -4474,7 +4069,6 @@ void VCS_SOLVE::vcs_dfe(const int stateCalc, } -//==================================================================================================================== void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const { @@ -4551,10 +4145,8 @@ void VCS_SOLVE::vcs_printSpeciesChemPot(const int stateCalc) const printf(" "); vcs_print_line("-", 132); } -/*****************************************************************************/ #ifdef DEBUG_MODE -//! Print out and check the elemental abundance vector void VCS_SOLVE::prneav() const { size_t j; @@ -4594,7 +4186,6 @@ void VCS_SOLVE::prneav() const plogendl(); } #endif -/*****************************************************************************/ double VCS_SOLVE::l2normdg(double dgLocal[]) const { @@ -4614,13 +4205,7 @@ double VCS_SOLVE::l2normdg(double dgLocal[]) const } return std::sqrt(tmp / m_numRxnRdc); } -/*****************************************************************************/ -// Calculates the total number of moles of species in all phases. -/* - * Calculates the total number of moles in all phases. - * Reconciles Phase existence flags with total moles in each phase. - */ double VCS_SOLVE::vcs_tmoles() { double sum; @@ -4674,7 +4259,6 @@ void VCS_SOLVE::check_tmoles() const } } #endif -/*****************************************************************************/ void VCS_SOLVE::vcs_updateVP(const int vcsState) { @@ -4700,55 +4284,6 @@ void VCS_SOLVE::vcs_updateVP(const int vcsState) } } -// This routine evaluates the species type for all species -/* - * kspec - * 1 -> Major species VCS_SPECIES_MAJOR - * 0 -> Minor species VCS_SPECIES_MINOR - * -1 -> The species lies in a multicomponent phase - * that exists. Its concentration is currently - * very low, necessitating a different method - * of calculation. - * - VCS_SPECIES_ZEROEDPHASE - * -2 -> The species lies in a multicomponent phase - * which currently doesn't exist. - * Its concentration is currently zero. - * - VCS_SPECIES_ZEROEDMS - * -3 -> Species lies in a single-species phase which - * is currently zeroed out. - * - VCS_SPECIES_ZEROEDSS - * -4 -> Species has such a small mole fraction it is - * deleted even though its phase may possibly exist. - * The species is believed to have such a small - * mole fraction that it best to throw the - * calculation of it out. It will be added back in - * at the end of the calculation. - * - VCS_SPECIES_DELETED - * -5 -> Species refers to an electron in the metal - * The unknown is equal to the interfacial voltage - * drop across the interface on the SHE (standard - * hydrogen electrode) scale (volts). - * - VCS_SPECIES_INTERFACIALVOLTAGE - * -6 -> Species lies in a multicomponent phase that - * is zeroed atm and will stay deleted due to a - * choice from a higher level. - * These species will formally always have zero - * mole numbers in the solution vector. - * - VCS_SPECIES_ZEROEDPHASE - * -7 -> The species lies in a multicomponent phase which - * currently does exist. Its concentration is currently - * identically zero, though the phase exists. Note, this - * is a temporary condition that exists at the start - * of an equilibrium problem. - * The species is soon "birthed" or "deleted". - * - VCS_SPECIES_ACTIVEBUTZERO - * -8 -> The species lies in a multicomponent phase which - * currently does exist. Its concentration is currently - * identically zero, though the phase exists. This is - * a permanent condition due to stoich constraints - * - VCS_SPECIES_STOICHZERO - * - */ bool VCS_SOLVE::vcs_evaluate_speciesType() { m_numRxnMinorZeroed = 0; @@ -4829,20 +4364,7 @@ bool VCS_SOLVE::vcs_evaluate_speciesType() return (m_numRxnMinorZeroed >= m_numRxnRdc); } -/*****************************************************************************/ -// Switch rows and columns of a square matrix -/* - * Switches the row and column of a matrix. - * So that after - * - * J[k1][j] = J_old[k2][j] and J[j][k1] = J_old[j][k2] - * J[k2][j] = J_old[k1][j] and J[j][k2] = J_old[j][k1] - * - * @param Jac Double pointer to the Jacobian - * @param k1 first row/column value to be switched - * @param k2 second row/column value to be switched - */ void VCS_SOLVE::vcs_switch2D(double* const* const Jac, const size_t k1, const size_t k2) const { @@ -4857,7 +4379,6 @@ void VCS_SOLVE::vcs_switch2D(double* const* const Jac, std::swap(Jac[i][k1], Jac[i][k2]); } } -/*****************************************************************************/ static void print_space(size_t num) { @@ -4866,37 +4387,7 @@ static void print_space(size_t num) plogf(" "); } } -/********************************************************************************/ -// This subroutine calculates reaction free energy changes for -// all noncomponent formation reactions. -/* - * Formation reactions are - * reactions which create each noncomponent species from the component - * species. m_stoichCoeffRxnMatrix[irxn][jcomp] are the stoichiometric - * coefficients for these reactions. A stoichiometric coefficient of - * one is assumed for species irxn in this reaction. - * - * INPUT - * @param l - * L < 0 : Calculate reactions corresponding to - * major noncomponent and zeroed species only - * L = 0 : Do all noncomponent reactions, i, between - * 0 <= i < irxnl - * L > 0 : Calculate reactions corresponding to - * minor noncomponent and zeroed species only - * - * @param doDeleted Do deleted species - * @param stateCalc Calculate deltaG corresponding to either old or new - * free energies - * - * Note we special case one important issue. - * If the component has zero moles, then we do not - * allow deltaG < 0.0 for formation reactions which - * would lead to the loss of more of that same component. - * This dG < 0.0 condition feeds back into the algorithm in several - * places, and leads to a infinite loop in at least one case. - */ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, const int vcsState, const bool alterZeroedPhases) { @@ -5107,7 +4598,7 @@ void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, } #endif } -//==================================================================================================================== + void VCS_SOLVE::vcs_printDeltaG(const int stateCalc) { size_t j; @@ -5241,22 +4732,7 @@ void VCS_SOLVE::vcs_printDeltaG(const int stateCalc) vcs_print_line("-", 132); } -//==================================================================================================================== -// Calculate deltag of formation for all species in a single phase. -/* - * Calculate deltag of formation for all species in a single - * phase. It is assumed that the fe[] is up to date for all species. - * However, if the phase is currently zeroed out, a subproblem - * is calculated to solve for AC[i] and pseudo-X[i] for that - * phase. - * - * @param iphase phase index of the phase to be calculated - * @param doDeleted boolean indicating whether to do deleted species or not - * @param stateCalc integer describing which set of free energies - * to use and where to stick the results. - * - * NOTE: this is currently not used used anywhere. It may be in the future? - */ + void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, const int stateCalc, const bool alterZeroedPhases) { @@ -5410,22 +4886,7 @@ void VCS_SOLVE::vcs_deltag_Phase(const size_t iphase, const bool doDeleted, } } } -/****************************************************************************/ -// Swaps the indices for all of the global data for two species, k1 -// and k2. -/* - * @param ifunc: If true, switch the species data and the noncomponent reaction - * data. This must be called for a non-component species only. - * If false, switch the species data only. Typically, we use this - * option when determining the component species and at the - * end of the calculation, when we want to return unscrambled - * results. All rxn data will be out-of-date. - * - * @param k1 First species index - * - * @param k2 Second species index - */ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k2) { size_t i1, i2, iph, kp1, kp2; @@ -5530,32 +4991,7 @@ void VCS_SOLVE::vcs_switch_pos(const bool ifunc, const size_t k1, const size_t k */ } } -/*******************************************************************************/ - -// Birth guess returns the number of moles of a species -// that is coming back to life. -/* - * Birth guess returns the number of moles of a species - * that is coming back to life. - * Note, this routine is not applicable if the whole phase is coming - * back to life, not just one species in that phase. - * - * Do a minor alt calculation. But, cap the mole numbers at - * 1.0E-15. - * For SS phases use VCS_DELETE_SPECIES_CUTOFF * 100. - * - * The routine makes sure the guess doesn't reduce the concentration - * of a component by more than 1/3. Note this may mean that - * the vlaue coming back from this routine is zero or a - * very small number. - * - * - * @param kspec Species number that is coming back to life - * - * @return Returns the number of kmol that the species should - * have. - */ double VCS_SOLVE::vcs_birthGuess(const int kspec) { size_t irxn = kspec - m_numComponents; @@ -5627,7 +5063,6 @@ double VCS_SOLVE::vcs_birthGuess(const int kspec) } return dx; } -/*******************************************************************************/ void VCS_SOLVE::vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc) { @@ -5644,7 +5079,6 @@ void VCS_SOLVE::vcs_setFlagsVolPhases(const bool upToDate, const int stateCalc) } } } -/*******************************************************************************/ void VCS_SOLVE::vcs_setFlagsVolPhase(const size_t iph, const bool upToDate, const int stateCalc) @@ -5656,16 +5090,7 @@ void VCS_SOLVE::vcs_setFlagsVolPhase(const size_t iph, const bool upToDate, Vphase->setMolesCurrent(stateCalc); } } -/*******************************************************************************/ -// Update all underlying vcs_VolPhase objects -/* - * Update the mole numbers and the phase voltages of all phases in the - * vcs problem - * - * @param stateCalc Location of the update (either VCS_STATECALC_NEW or - * VCS_STATECALC_OLD). - */ void VCS_SOLVE::vcs_updateMolNumVolPhases(const int stateCalc) { vcs_VolPhase* Vphase; @@ -5674,6 +5099,5 @@ void VCS_SOLVE::vcs_updateMolNumVolPhases(const int stateCalc) Vphase->updateFromVCS_MoleNumbers(stateCalc); } } -/*******************************************************************************/ } diff --git a/src/equil/vcs_solve_phaseStability.cpp b/src/equil/vcs_solve_phaseStability.cpp index 143a7a7a3..941f27825 100644 --- a/src/equil/vcs_solve_phaseStability.cpp +++ b/src/equil/vcs_solve_phaseStability.cpp @@ -27,7 +27,6 @@ using namespace std; namespace VCSnonideal { - int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStable) { @@ -163,38 +162,8 @@ int VCS_SOLVE::vcs_PS(VCS_PROB* vprob, int iphase, int printLvl, double& feStabl * Return the convergence success flag. */ return iStab; - - } -//==================================================================================================================== -// Routine that independently determines whether a phase should be popped -// under the current conditions. -/* - * This is the main routine that solves for equilibrium at constant T and P - * using a variant of the VCS method. Nonideal phases can be accommodated - * as well. - * - * Any number of single-species phases and multi-species phases - * can be handled by the present version. - * - * Input - * ------------ - * @param print_lvl 1 -> Print results to standard output - * 0 -> don't report on anything - * - * @param printDetails 1 -> Print intermediate results. - * - * @param maxit Maximum number of iterations for the algorithm - * - * @return 0 = Equilibrium Achieved - * 1 = Range space error encountered. The element abundance criteria are - * only partially satisfied. Specifically, the first NC= (number of - * components) conditions are satisfied. However, the full NE - * (number of elements) conditions are not satisfied. The equilibrium - * condition is returned. - * -1 = Maximum number of iterations is exceeded. Convergence was not - * found. - */ + int VCS_SOLVE::vcs_solve_phaseStability(const int iph, const int ifunc, double& funcVal, int printLvl)