From 2ddf9a409f2ee808cbd37be3375c6ac43fe0fa9e Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Mon, 2 Jun 2008 23:24:32 +0000 Subject: [PATCH] Worked on making sure that vcs_deltag() adhered to the STATECALC concept. --- Cantera/src/equil/vcs_defs.h | 6 +- Cantera/src/equil/vcs_inest.cpp | 46 +++++----- Cantera/src/equil/vcs_internal.h | 32 +++---- Cantera/src/equil/vcs_report.cpp | 3 +- Cantera/src/equil/vcs_solve.h | 32 ++++++- Cantera/src/equil/vcs_solve_TP.cpp | 134 ++++++++++++++++------------- 6 files changed, 152 insertions(+), 101 deletions(-) diff --git a/Cantera/src/equil/vcs_defs.h b/Cantera/src/equil/vcs_defs.h index 0e8d70308..dcdf290dd 100644 --- a/Cantera/src/equil/vcs_defs.h +++ b/Cantera/src/equil/vcs_defs.h @@ -82,11 +82,11 @@ namespace VCSnonideal { * @{ */ - //! Cutoff relative mole number value, + //! Cutoff relative mole fraction value, //! below which species are deleted from the equilibrium problem. -#ifndef VCS_DELETE_SPECIES_CUTOFF -#define VCS_DELETE_SPECIES_CUTOFF 1.0e-32 +#ifndef VCS_RELDELETE_SPECIES_CUTOFF +#define VCS_RELDELETE_SPECIES_CUTOFF 1.0e-64 #endif //! Cutoff relative mole number value, diff --git a/Cantera/src/equil/vcs_inest.cpp b/Cantera/src/equil/vcs_inest.cpp index 7a2a8a8fd..5bafe6642 100644 --- a/Cantera/src/equil/vcs_inest.cpp +++ b/Cantera/src/equil/vcs_inest.cpp @@ -51,7 +51,7 @@ namespace VCSnonideal { int nrxn = m_numRxnTot; vcs_VolPhase *Vphase = 0; - double *molNum = VCS_DATA_PTR(m_molNumSpecies_old); + // double *molNum = VCS_DATA_PTR(m_molNumSpecies_old); double TMolesMultiphase; double *xtphMax = VCS_DATA_PTR(m_TmpPhase); double *xtphMin = VCS_DATA_PTR(m_TmpPhase2); @@ -104,7 +104,7 @@ namespace VCSnonideal { plogf("%s SPECIES MOLE_NUMBER -SS_ChemPotential\n", pprefix); for (kspec = 0; kspec < nspecies; ++kspec) { plogf("%s ", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str()); - plogf(" %15.5g %12.3g\n", molNum[kspec], -m_SSfeSpecies[kspec]); + plogf(" %15.5g %12.3g\n", m_molNumSpecies_old[kspec], -m_SSfeSpecies[kspec]); } plogf("%s Element Abundance Agreement returned from linear " "programming (vcs_inest initial guess):", @@ -116,7 +116,7 @@ namespace VCSnonideal { if (m_elementActive[j]) { double tmp = 0.0; for (kspec = 0; kspec < nspecies; ++kspec) { - tmp += m_formulaMatrix[j][kspec] * molNum[kspec]; + tmp += m_formulaMatrix[j][kspec] * m_molNumSpecies_old[kspec]; } plogf("%s ", pprefix); plogf(" %-9.9s", (m_elementName[j]).c_str()); plogf(" %12.3g %12.3g\n", m_elemAbundancesGoal[j], tmp); @@ -137,18 +137,18 @@ namespace VCSnonideal { iph = m_phaseID[kspec]; Vphase = m_VolPhaseList[iph]; if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (molNum[kspec] <= 0.0) { + if (m_molNumSpecies_old[kspec] <= 0.0) { /* * HKM Should eventually include logic here for non SS phases */ if (!m_SSPhase[kspec]) { - molNum[kspec] = 1.0e-30; + m_molNumSpecies_old[kspec] = 1.0e-30; } } } else { - molNum[kspec] = 0.0; + m_molNumSpecies_old[kspec] = 0.0; } - if (molNum[kspec] > 0.0) { + if (m_molNumSpecies_old[kspec] > 0.0) { if (Vphase->Existence == 0) { Vphase->Existence = 1; } @@ -179,7 +179,7 @@ namespace VCSnonideal { } for (kspec = 0; kspec < m_numComponents; ++kspec) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { - m_tPhaseMoles_new[m_phaseID[kspec]] += molNum[kspec]; + m_tPhaseMoles_new[m_phaseID[kspec]] += m_molNumSpecies_old[kspec]; } } TMolesMultiphase = 0.0; @@ -188,36 +188,36 @@ namespace VCSnonideal { TMolesMultiphase += m_tPhaseMoles_new[iph]; } } - vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), molNum, nspecies); + vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_new), VCS_DATA_PTR(m_molNumSpecies_old), nspecies); for (kspec = 0; kspec < m_numComponents; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_MOLNUM) { m_molNumSpecies_new[kspec] = 0.0; } } - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(m_SSfeSpecies), + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_SSfeSpecies), nspecies); for (kspec = 0; kspec < m_numComponents; ++kspec) { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { if (! m_SSPhase[kspec]) { iph = m_phaseID[kspec]; - m_feSpecies_curr[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]); + m_feSpecies_new[kspec] += log(m_molNumSpecies_new[kspec] / m_tPhaseMoles_old[iph]); } } else { m_molNumSpecies_new[kspec] = 0.0; } } - vcs_deltag(0, true); + vcs_deltag(0, true, VCS_STATECALC_NEW); #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { for (kspec = 0; kspec < nspecies; ++kspec) { plogf("%s", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str()); if (kspec < m_numComponents) - plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_curr[kspec], + plogf("fe* = %15.5g ff = %15.5g\n", m_feSpecies_new[kspec], m_SSfeSpecies[kspec]); else plogf("fe* = %15.5g ff = %15.5g dg* = %15.5g\n", - m_feSpecies_curr[kspec], m_SSfeSpecies[kspec], m_deltaGRxn_new[kspec-m_numComponents]); + m_feSpecies_new[kspec], m_SSfeSpecies[kspec], m_deltaGRxn_new[kspec-m_numComponents]); } } #endif @@ -269,8 +269,8 @@ namespace VCSnonideal { plogf("%sdirection (", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str()); plogf(") = %g", m_deltaMolNumSpecies[kspec]); if (m_SSPhase[kspec]) { - if (molNum[kspec] > 0.0) { - plogf(" (ssPhase exists at w = %g moles)", molNum[kspec]); + if (m_molNumSpecies_old[kspec] > 0.0) { + plogf(" (ssPhase exists at w = %g moles)", m_molNumSpecies_old[kspec]); } else { plogf(" (ssPhase doesn't exist -> stability not checked)"); } @@ -304,14 +304,16 @@ namespace VCSnonideal { do { for (kspec = 0; kspec < m_numComponents; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - molNum[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec]; + m_molNumSpecies_old[kspec] = m_molNumSpecies_new[kspec] + par * m_deltaMolNumSpecies[kspec]; } else { m_deltaMolNumSpecies[kspec] = 0.0; } } for (kspec = m_numComponents; kspec < nspecies; ++kspec) { if (m_speciesUnknownType[kspec] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - if (m_deltaMolNumSpecies[kspec] != 0.0) molNum[kspec] = m_deltaMolNumSpecies[kspec] * par; + if (m_deltaMolNumSpecies[kspec] != 0.0) { + m_molNumSpecies_old[kspec] = m_deltaMolNumSpecies[kspec] * par; + } } } /* @@ -323,9 +325,9 @@ namespace VCSnonideal { /* ******************************************* */ /* **** CONVERGENCE FORCING SECTION ********** */ /* ******************************************* */ - vcs_dfe(molNum, VCS_STATECALC_OLD, 0, 0, nspecies); + vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, nspecies); for (kspec = 0, s = 0.0; kspec < nspecies; ++kspec) { - s += m_deltaMolNumSpecies[kspec] * m_feSpecies_curr[kspec]; + s += m_deltaMolNumSpecies[kspec] * m_feSpecies_old[kspec]; } if (s == 0.0) { finished = TRUE; continue; @@ -377,7 +379,7 @@ namespace VCSnonideal { plogf("%s SPECIES MOLE_NUMBER\n", pprefix); for (kspec = 0; kspec < nspecies; ++kspec) { plogf("%s ", pprefix); plogf("%-12.12s", m_speciesName[kspec].c_str()); - plogf(" %g", molNum[kspec]); + plogf(" %g", m_molNumSpecies_old[kspec]); plogendl(); } } @@ -520,7 +522,7 @@ namespace VCSnonideal { #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { plogf("%sTotal Dimensionless Gibbs Free Energy = %15.7E", pprefix, - vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), + vcs_Total_Gibbs(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_tPhaseMoles_old))); plogendl(); } diff --git a/Cantera/src/equil/vcs_internal.h b/Cantera/src/equil/vcs_internal.h index bb7743873..f00cee2d0 100644 --- a/Cantera/src/equil/vcs_internal.h +++ b/Cantera/src/equil/vcs_internal.h @@ -325,7 +325,7 @@ namespace VCSnonideal { * @param vec_to vector of doubles * @param length length of the vector to zero. */ - inline void vcs_dzero(double *vec_to, int length) { + inline void vcs_dzero(double * const vec_to, const int length) { (void) memset((void *) vec_to, 0, length * sizeof(double)); } @@ -334,7 +334,7 @@ namespace VCSnonideal { * @param vec_to vector of ints * @param length length of the vector to zero. */ - inline void vcs_izero(int *vec_to, int length) { + inline void vcs_izero(int * const vec_to, const int length) { (void) memset((void *) vec_to, 0, length * sizeof(int)); } @@ -345,7 +345,8 @@ namespace VCSnonideal { * @param vec_from Vector to copy from * @param length Number of doubles to copy. */ - inline void vcs_dcopy(double *vec_to, const double *vec_from, int length) { + inline void vcs_dcopy(double * const vec_to, + const double * const vec_from, const int length) { (void) memcpy((void *) vec_to, (const void *) vec_from, (length) * sizeof(double)); } @@ -358,7 +359,8 @@ namespace VCSnonideal { * @param vec_from Vector to copy from * @param length Number of int to copy. */ - inline void vcs_icopy(int *vec_to, const int *vec_from, int length) { + inline void vcs_icopy(int * const vec_to, + const int * const vec_from, const int length) { (void) memcpy((void *) vec_to, (const void *) vec_from, (length) * sizeof(int)); } @@ -368,7 +370,7 @@ namespace VCSnonideal { * @param vec_to vector of doubles * @param length length of the vector to zero. */ - inline void vcs_vdzero(std::vector &vec_to, int length) { + inline void vcs_vdzero(std::vector &vec_to, const int length) { (void) memset((void *)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(double)); } @@ -377,7 +379,7 @@ namespace VCSnonideal { * @param vec_to vector of ints * @param length length of the vector to zero. */ - inline void vcs_vizero(std::vector &vec_to, int length) { + inline void vcs_vizero(std::vector &vec_to, const int length) { (void) memset((void *)VCS_DATA_PTR(vec_to), 0, (length) * sizeof(int)); } @@ -409,21 +411,21 @@ namespace VCSnonideal { * @param length Number of integers to copy. */ inline void vcs_vicopy(std::vector & vec_to, - const std::vector & vec_from, int length) { + const std::vector & vec_from, const int length) { (void) memcpy((void *)&(vec_to[0]), (const void *) &(vec_from[0]), (length) * sizeof(int)); } #else - extern void vcs_dzero(double *, int); - extern void vcs_izero(int *, int); - extern void vcs_dcopy(double *, double *, int); - extern void vcs_icopy(int *, int *, int); - extern void vcs_vdzero(std::vector &vvv, int len = -1); - extern void vcs_vizero(std::vector &vvv, int len = -1); + extern void vcs_dzero(double * const, const int); + extern void vcs_izero(int * const , const int); + extern void vcs_dcopy(double * const, const double * const, const int); + extern void vcs_icopy(int * const, const int * const, const int); + extern void vcs_vdzero(std::vector &vvv, const int len = -1); + extern void vcs_vizero(std::vector &vvv, const int len = -1); void vcs_vdcopy(std::vector &vec_to, - const std::vector vec_from, int len = -1); + const std::vector vec_from, const int len = -1); void vcs_vicopy(std::vector &vec_to, - const std::vector vec_from, int len = -1); + const std::vector vec_from, const int len = -1); #endif //! Finds the location of the maximum component in a double vector diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 889c892ac..abc66a7fc 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -329,7 +329,8 @@ namespace VCSnonideal { lx = 0.0; } else { if (tpmoles > 0.0 && m_molNumSpecies_old[l] > 0.0) { - lx = log(m_molNumSpecies_old[l]) - log(tpmoles); + double tmp = MAX(VCS_DELETE_MINORSPECIES_CUTOFF, m_molNumSpecies_old[l]); + lx = log(tmp) - log(tpmoles); } else { lx = m_feSpecies_old[l] - m_SSfeSpecies[l] - log(m_actCoeffSpecies_old[l]) + m_lnMnaughtSpecies[l]; diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index b3f03cd79..b23762d41 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -504,7 +504,37 @@ public: */ void vcs_tmoles(); - void vcs_deltag(int l, bool doDeleted); + + //! 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_deltag(const int l, const bool doDeleted, const int stateCalc); //! Swaps the indecises for all of the global data for two species, k1 //! and k2. diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index 0202623bc..97c40e101 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -267,7 +267,12 @@ namespace VCSnonideal { "negative MF, setting it small", m_speciesName[i].c_str()); plogendl(); - m_molNumSpecies_old[i] = VCS_DELETE_SPECIES_CUTOFF; + iph = m_phaseID[i]; + double tmp = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10; + if ( VCS_DELETE_MINORSPECIES_CUTOFF*10. > tmp) { + tmp = VCS_DELETE_MINORSPECIES_CUTOFF*10.; + } + m_molNumSpecies_old[i] = tmp; } } @@ -390,7 +395,7 @@ namespace VCSnonideal { #endif // Update the phase objects with the contents of the soln vector vcs_updateVP(0); - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_OLD); iti = 0; goto L_MAINLOOP_ALL_SPECIES; @@ -416,7 +421,7 @@ namespace VCSnonideal { */ if (uptodate_minors == FALSE) { vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); - vcs_deltag(1, false); + vcs_deltag(1, false, VCS_STATECALC_NEW); } uptodate_minors = TRUE; } else { @@ -435,10 +440,13 @@ namespace VCSnonideal { } } - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); - vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); - vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_old), VCS_DATA_PTR(m_actCoeffSpecies_new), m_numSpeciesRdc); - vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc); + /* + * Copy the old solution into the new solution as an initial guess + */ + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_curr), VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_feSpecies_new), VCS_DATA_PTR(m_feSpecies_old), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_new), VCS_DATA_PTR(m_actCoeffSpecies_old), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_new), VCS_DATA_PTR(m_deltaGRxn_old), m_numRxnRdc); /* Go find a new reaction adjustment -> * i.e., change in extent of reaction for each reaction. @@ -1094,7 +1102,7 @@ namespace VCSnonideal { * Evaluate DeltaG for all components if ITI=0, and for * major components only if ITI NE 0 */ - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_NEW); /* *************************************************************** */ /* **** CONVERGENCE FORCER SECTION ******************************* */ @@ -1231,6 +1239,8 @@ namespace VCSnonideal { vcs_dcopy(VCS_DATA_PTR(m_tPhaseMoles_old), VCS_DATA_PTR(m_tPhaseMoles_new), m_numPhases); vcs_dcopy(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_molNumSpecies_new), m_numSpeciesRdc); + vcs_dcopy(VCS_DATA_PTR(m_actCoeffSpecies_old), + VCS_DATA_PTR(m_actCoeffSpecies_new), m_numSpeciesRdc); vcs_dcopy(VCS_DATA_PTR(m_deltaGRxn_old), VCS_DATA_PTR(m_deltaGRxn_new), m_numRxnRdc); vcs_dcopy(VCS_DATA_PTR(m_feSpecies_old), VCS_DATA_PTR(m_feSpecies_curr), m_numSpeciesRdc); @@ -1268,7 +1278,7 @@ namespace VCSnonideal { if (m_phaseID[kspec] == iph && m_molNumSpecies_old[kspec] > 0.0) { irxn = kspec - m_numComponents; if (kspec < m_numComponents) { - if (m_molNumSpecies_old[kspec] > VCS_DELETE_SPECIES_CUTOFF) { + if (m_molNumSpecies_old[kspec] > VCS_RELDELETE_SPECIES_CUTOFF) { soldel = 0; break; } @@ -1311,7 +1321,7 @@ namespace VCSnonideal { &usedZeroedSpecies); if (retn != VCS_SUCCESS) return retn; vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); - vcs_deltag(0, true); + vcs_deltag(0, true, VCS_STATECALC_OLD); uptodate_minors = TRUE; if (conv) { /* @@ -1344,7 +1354,7 @@ namespace VCSnonideal { #endif vcs_elcorr(VCS_DATA_PTR(sm), VCS_DATA_PTR(wx)); vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); - vcs_deltag(0, true); + vcs_deltag(0, true, VCS_STATECALC_OLD); uptodate_minors = TRUE; } #ifdef DEBUG_MODE @@ -1579,7 +1589,7 @@ namespace VCSnonideal { */ if (iti != 0) { vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, kspec, kspec+1); - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_OLD); } } } @@ -1660,7 +1670,7 @@ namespace VCSnonideal { */ if (iti != 0) { vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); - vcs_deltag(1, false); + vcs_deltag(1, false, VCS_STATECALC_OLD); uptodate_minors = TRUE; } #ifdef DEBUG_MODE @@ -1770,7 +1780,7 @@ namespace VCSnonideal { * Go back to evaluate the total moles of gas and liquid. */ vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesRdc); - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_OLD); /* * */ @@ -1829,7 +1839,7 @@ namespace VCSnonideal { } } // Calculate delta g's - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_OLD); // Go back to equilibrium check as a prep to eventually checking out goto L_EQUILIB_CHECK; @@ -1860,7 +1870,7 @@ namespace VCSnonideal { */ MajorSpeciesHaveConverged = true; vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_OLD); iti = 0; goto L_MAINLOOP_ALL_SPECIES; /*************************************************************************/ @@ -1879,7 +1889,7 @@ namespace VCSnonideal { */ MajorSpeciesHaveConverged = true; vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 1, 0, m_numSpeciesRdc); - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_OLD); iti = 0; goto L_MAINLOOP_ALL_SPECIES; } @@ -2336,7 +2346,7 @@ namespace VCSnonideal { */ void VCS_SOLVE::vcs_reinsert_deleted(int kspec) { int i, k, irxn = kspec - m_numComponents; - int *phaseID = VCS_DATA_PTR(m_phaseID); + int iph = m_phaseID[kspec]; double dx; #ifdef DEBUG_MODE if (m_debug_print_lvl >= 2) { @@ -2348,7 +2358,7 @@ namespace VCSnonideal { * this adjusts m_molNumSpecies_old[] and m_tPhaseMoles_old[] * HKM -> make this a relative mole number! */ - dx = VCS_DELETE_SPECIES_CUTOFF * 10.; + dx = m_tPhaseMoles_old[iph] * VCS_RELDELETE_SPECIES_CUTOFF * 10.; delta_species(kspec, &dx); m_rxnStatus[irxn] = VCS_SPECIES_MINOR; @@ -2356,7 +2366,7 @@ namespace VCSnonideal { m_rxnStatus[irxn] = VCS_SPECIES_MAJOR; --(m_numRxnMinorZeroed); } - int iph = m_phaseID[kspec]; + vcs_VolPhase *Vphase = m_VolPhaseList[iph]; Vphase->setMolesFromVCSCheck(VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_tPhaseMoles_old)); @@ -2373,7 +2383,7 @@ namespace VCSnonideal { if (Vphase->Existence == 0) { Vphase->Existence = 1; for (k = 0; k < m_numSpeciesTot; k++) { - if (phaseID[k] == iph) { + if (m_phaseID[k] == iph) { i = k - m_numComponents; if (m_rxnStatus[i] == VCS_SPECIES_ZEROEDPHASE) m_rxnStatus[i] = VCS_SPECIES_ZEROEDMS; @@ -2547,7 +2557,7 @@ namespace VCSnonideal { */ for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { iph = m_phaseID[kspec]; - m_feSpecies_curr[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) + m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); } @@ -2556,11 +2566,11 @@ namespace VCSnonideal { * Recalculate the DeltaG's of the formation reactions for the * deleted species in the mechanism */ - vcs_deltag(0, true); + vcs_deltag(0, true, VCS_STATECALC_NEW); for (iph = 0; iph < m_numPhases; iph++) { if (m_tPhaseMoles_old[iph] > 0.0) - xtcutoff[iph] = log (m_tPhaseMoles_old[iph] / VCS_DELETE_SPECIES_CUTOFF); + xtcutoff[iph] = log (1.0 / VCS_RELDELETE_SPECIES_CUTOFF); else xtcutoff[iph] = 0.0; } @@ -2706,7 +2716,7 @@ namespace VCSnonideal { */ for (kspec = m_numSpeciesRdc; kspec < m_numSpeciesTot; ++kspec) { iph = m_phaseID[kspec]; - m_feSpecies_curr[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) + m_feSpecies_new[kspec] = (m_SSfeSpecies[kspec] + log(m_actCoeffSpecies_old[kspec]) - m_lnMnaughtSpecies[kspec] + m_chargeSpecies[kspec] * m_Faraday_dim * m_phasePhi[iph]); } @@ -2714,13 +2724,13 @@ namespace VCSnonideal { * Recalculate the DeltaG's of the formation reactions for the * deleted species in the mechanism */ - vcs_deltag(0, true); + vcs_deltag(0, true, VCS_STATECALC_NEW); for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { kspec = m_indexRxnToSpecies[irxn]; iph = m_phaseID[kspec]; if (m_tPhaseMoles_old[iph] > 0.0) { - double maxDG = MIN(m_deltaGRxn_new[irxn], 300.0); + double maxDG = MIN(m_deltaGRxn_new[irxn], 690.0); double dx = m_tPhaseMoles_old[iph] * exp(- maxDG); retn = delta_species(kspec, &dx); @@ -2762,7 +2772,7 @@ namespace VCSnonideal { } vcs_dfe(VCS_DATA_PTR(m_molNumSpecies_old), VCS_STATECALC_OLD, 0, 0, m_numSpeciesTot); - vcs_deltag(0, true); + vcs_deltag(0, true, VCS_STATECALC_OLD); retn = 0; for (int irxn = m_numRxnRdc; irxn < m_numRxnTot; ++irxn) { @@ -2770,15 +2780,18 @@ namespace VCSnonideal { iph = m_phaseID[kspec]; if (m_tPhaseMoles_old[iph] > 0.0) { if (fabs(m_deltaGRxn_old[irxn]) > m_tolmin) { - retn++; + if (((m_molNumSpecies_old[kspec] * exp(-m_deltaGRxn_old[irxn])) > VCS_DELETE_MINORSPECIES_CUTOFF) || + (m_molNumSpecies_old[kspec] > VCS_DELETE_MINORSPECIES_CUTOFF)) { + retn++; #ifdef DEBUG_MODE - if (m_debug_print_lvl >= 2) { - plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", - m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], - m_deltaGRxn_old[irxn]); - plogendl(); - } + if (m_debug_print_lvl >= 2) { + plogf(" --- add_deleted(): species %s with mol number %g not converged: DG = %g", + m_speciesName[kspec].c_str(), m_molNumSpecies_old[kspec], + m_deltaGRxn_old[irxn]); + plogendl(); + } #endif + } } } } @@ -2926,7 +2939,7 @@ namespace VCSnonideal { * Evaluate DeltaG for all components if ITI=0, and for * major components only if ITI NE 0 */ - vcs_deltag(0, false); + vcs_deltag(0, false, VCS_STATECALC_NEW); dptr = VCS_DATA_PTR(m_deltaGRxn_new); s2 = 0.0; @@ -3274,34 +3287,38 @@ namespace VCSnonideal { } /*****************************************************************************/ + + // 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. * - * vcs_deltag: - * - * 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(J,I) are the stoichiometric coefficients for these - * reactions. A stoichiometric coefficient of one is assumed for - * species I in this reaction. - * - * INPUT + * 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 - * irxnl : used with L = 0 to indicate upper limit. * - * 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 the component. - * This dG < 0.0 feeds back into the algorithm in several - * places, and leads to a infinite loop in at least one case. + * @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(int l, bool doDeleted) { + void VCS_SOLVE::vcs_deltag(const int l, const bool doDeleted, + const int vcsState) { int iph; int lneed, irxn, kspec; double *dtmp_ptr; @@ -3310,7 +3327,6 @@ namespace VCSnonideal { if (doDeleted) { irxnl = m_numRxnTot; } - int vcsState = VCS_STATECALC_NEW; double * deltaGRxn; double * feSpecies; @@ -3318,10 +3334,10 @@ namespace VCSnonideal { double * actCoeffSpecies; if (vcsState == VCS_STATECALC_NEW) { deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_new); - feSpecies = VCS_DATA_PTR(m_feSpecies_curr); + feSpecies = VCS_DATA_PTR(m_feSpecies_new); molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_new); actCoeffSpecies = VCS_DATA_PTR(m_actCoeffSpecies_new); - } else if (vcsState == VCS_STATECALC_NEW) { + } else if (vcsState == VCS_STATECALC_OLD) { deltaGRxn = VCS_DATA_PTR(m_deltaGRxn_old); feSpecies = VCS_DATA_PTR(m_feSpecies_old); molNumSpecies = VCS_DATA_PTR(m_molNumSpecies_old); @@ -5292,7 +5308,7 @@ namespace VCSnonideal { if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { return dx; } - double w_kspec = VCS_DELETE_SPECIES_CUTOFF; + double w_kspec = VCS_DELETE_MINORSPECIES_CUTOFF; #ifdef DEBUG_MODE // Check to make sure that species is zero in the solution vector // If it isn't, we don't know what's happening @@ -5324,7 +5340,7 @@ namespace VCSnonideal { * There is no real way to estimate the moles. So * we set it to a small number. */ - dx = VCS_DELETE_SPECIES_CUTOFF * 100.; + dx = 1.0E-30; } /*