From 3820ca3147eb9fc5770364cb1b9ca62624115e58 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Wed, 23 Aug 2017 14:34:45 -0400 Subject: [PATCH] [Equil] Eliminate "total moles" scaling in VCS solver --- include/cantera/equil/vcs_solve.h | 9 ---- src/equil/vcs_nondim.cpp | 76 ------------------------------- src/equil/vcs_report.cpp | 40 +++++++--------- src/equil/vcs_solve.cpp | 1 - src/equil/vcs_solve_TP.cpp | 3 -- 5 files changed, 16 insertions(+), 113 deletions(-) diff --git a/include/cantera/equil/vcs_solve.h b/include/cantera/equil/vcs_solve.h index d18e9d81b..ae6fc97f2 100644 --- a/include/cantera/equil/vcs_solve.h +++ b/include/cantera/equil/vcs_solve.h @@ -1426,15 +1426,6 @@ public: */ char m_unitsState; - //! Multiplier for the mole numbers within the nondimensional formulation - /*! - * All numbers within the main routine are on an absolute basis. This - * presents some problems wrt very large and very small mole numbers. We get - * around this by using a multiplier coming into and coming out of the - * equilibrium routines - */ - double m_totalMoleScale; - //! specifies the activity convention of the phase containing the species /*! * * 0 = molar based diff --git a/src/equil/vcs_nondim.cpp b/src/equil/vcs_nondim.cpp index e4863217e..fa151678e 100644 --- a/src/equil/vcs_nondim.cpp +++ b/src/equil/vcs_nondim.cpp @@ -30,60 +30,6 @@ void VCS_SOLVE::vcs_nondim_TP() } m_Faraday_dim = ElectronCharge * Avogadro / (m_temperature * GasConstant); - - // Scale the total moles if necessary: First find out the total moles - double tmole_orig = vcs_tmoles(); - - // Then add in the total moles of elements that are goals. Either one or - // the other is specified here. - double esum = 0.0; - for (size_t i = 0; i < m_nelem; ++i) { - if (m_elType[i] == VCS_ELEM_TYPE_ABSPOS) { - esum += fabs(m_elemAbundancesGoal[i]); - } - } - tmole_orig += esum; - - // Ok now test out the bounds on the total moles that this program can - // handle. These are a bit arbitrary. However, it would seem that any - // reasonable input would be between these two numbers below. - if (tmole_orig < 1.0E-200 || tmole_orig > 1.0E200) { - throw CanteraError("VCS_SOLVE::vcs_nondim_TP", - "Total input moles, {} is outside the range handled by vcs.\n", - tmole_orig); - } - - // Determine the scale of the problem - if (tmole_orig > 1.0E4) { - m_totalMoleScale = tmole_orig / 1.0E4; - } else if (tmole_orig < 1.0E-4) { - m_totalMoleScale = tmole_orig / 1.0E-4; - } else { - m_totalMoleScale = 1.0; - } - - if (m_totalMoleScale != 1.0) { - if (m_debug_print_lvl >= 2) { - plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice\n", m_totalMoleScale); - } - for (size_t i = 0; i < m_nsp; ++i) { - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale); - } - } - for (size_t i = 0; i < m_nelem; ++i) { - m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale); - } - - for (size_t iph = 0; iph < m_numPhases; iph++) { - TPhInertMoles[iph] *= (1.0 / m_totalMoleScale); - if (TPhInertMoles[iph] != 0.0) { - vcs_VolPhase* vphase = m_VolPhaseList[iph].get(); - vphase->setTotalMolesInert(TPhInertMoles[iph]); - } - } - vcs_tmoles(); - } } } @@ -103,28 +49,6 @@ void VCS_SOLVE::vcs_redim_TP() } m_Faraday_dim *= tf; } - if (m_totalMoleScale != 1.0) { - if (m_debug_print_lvl >= 2) { - plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g\n", m_totalMoleScale); - } - for (size_t i = 0; i < m_nsp; ++i) { - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - m_molNumSpecies_old[i] *= m_totalMoleScale; - } - } - for (size_t i = 0; i < m_nelem; ++i) { - m_elemAbundancesGoal[i] *= m_totalMoleScale; - } - - for (size_t iph = 0; iph < m_numPhases; iph++) { - TPhInertMoles[iph] *= m_totalMoleScale; - if (TPhInertMoles[iph] != 0.0) { - vcs_VolPhase* vphase = m_VolPhaseList[iph].get(); - vphase->setTotalMolesInert(TPhInertMoles[iph]); - } - } - vcs_tmoles(); - } } } diff --git a/src/equil/vcs_report.cpp b/src/equil/vcs_report.cpp index 1f27388bc..8cefba882 100644 --- a/src/equil/vcs_report.cpp +++ b/src/equil/vcs_report.cpp @@ -11,7 +11,7 @@ namespace Cantera { int VCS_SOLVE::vcs_report(int iconv) { - bool printActualMoles = true, inertYes = false; + bool inertYes = false; char originalUnitsState = m_unitsState; std::vector sortindex(m_nsp, 0); vector_fp xy(m_nsp, 0.0); @@ -38,10 +38,6 @@ int VCS_SOLVE::vcs_report(int iconv) if (m_unitsState == VCS_DIMENSIONAL_G) { vcs_nondim_TP(); } - double molScale = 1.0; - if (printActualMoles) { - molScale = m_totalMoleScale; - } vcs_setFlagsVolPhases(false, VCS_STATECALC_OLD); vcs_dfe(VCS_STATECALC_OLD, 0, 0, m_nsp); @@ -66,11 +62,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature); plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA); - plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol * molScale); - if (!printActualMoles) { - plogf("\t\tMole Scale = %15.5g kmol (all mole numbers and volumes are scaled by this value)\n", - molScale); - } + plogf("\t\ttotal Volume = %15.5g m**3\n", m_totalVol); // TABLE OF SPECIES IN DECREASING MOLE NUMBERS plogf("\n\n"); @@ -81,8 +73,8 @@ int VCS_SOLVE::vcs_report(int iconv) for (size_t i = 0; i < m_numComponents; ++i) { plogf(" %-12.12s", m_speciesName[i]); writeline(' ', 13, false); - plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i] * molScale, - m_molNumSpecies_new[i] * molScale, m_feSpecies_old[i]); + plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i], + m_molNumSpecies_new[i], m_feSpecies_old[i]); plogf(" %3d", m_speciesUnknownType[i]); plogf("\n"); } @@ -92,12 +84,12 @@ int VCS_SOLVE::vcs_report(int iconv) writeline(' ', 13, false); if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_MOLNUM) { - plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j] * molScale, - m_molNumSpecies_new[j] * molScale, m_feSpecies_old[j]); + plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[j], + m_molNumSpecies_new[j], m_feSpecies_old[j]); plogf(" KMolNum "); } else if (m_speciesUnknownType[j] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[j]); - plogf(" Voltage = %14.7E", m_molNumSpecies_old[j] * molScale); + plogf(" Voltage = %14.7E", m_molNumSpecies_old[j]); } else { throw CanteraError("VCS_SOLVE::vcs_report", "we have a problem"); } @@ -112,7 +104,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" Inert Species in phase %16s ", m_VolPhaseList[i]->PhaseName); } - plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i] * molScale, + plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i], TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0); } } @@ -122,8 +114,8 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" %-12.12s", m_speciesName[kspec]); // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve plogf(" %14.7E %14.7E %12.4E", - m_molNumSpecies_old[kspec]*molScale, - m_molNumSpecies_new[kspec]*molScale, m_deltaGRxn_new[kspec]); + m_molNumSpecies_old[kspec], + m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]); if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_MOLNUM) { plogf(" KMol_Num"); } else if (m_speciesUnknownType[kspec] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { @@ -151,7 +143,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" | |\n"); plogf(" NonComponent | Moles |"); for (size_t j = 0; j < m_numComponents; j++) { - plogf(" %10.3g", m_molNumSpecies_old[j] * molScale); + plogf(" %10.3g", m_molNumSpecies_old[j]); } plogf(" | DG/RT Rxn |\n"); writeline('-', m_numComponents*10 + 45); @@ -159,7 +151,7 @@ int VCS_SOLVE::vcs_report(int iconv) size_t kspec = m_indexRxnToSpecies[irxn]; plogf(" %3d ", kspec); plogf("%-10.10s", m_speciesName[kspec]); - plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale); + plogf("|%10.3g |", m_molNumSpecies_old[kspec]); for (size_t j = 0; j < m_numComponents; j++) { plogf(" %6.2f", m_stoichCoeffRxnMatrix(j,irxn)); } @@ -198,7 +190,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" %3d ", iphase); vcs_VolPhase* VPhase = m_VolPhaseList[iphase].get(); plogf("%-12.12s |",VPhase->PhaseName); - plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale); + plogf("%10.3e |", m_tPhaseMoles_old[iphase]); totalMoles += m_tPhaseMoles_old[iphase]; if (m_tPhaseMoles_old[iphase] != VPhase->totalMoles() && !vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->totalMoles())) { @@ -240,7 +232,7 @@ int VCS_SOLVE::vcs_report(int iconv) for (size_t i = 0; i < m_nelem; ++i) { writeline(' ', 26, false); plogf("%-2.2s", m_elementName[i]); - plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale); + plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]); plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]); } plogf("\n"); @@ -258,7 +250,7 @@ int VCS_SOLVE::vcs_report(int iconv) size_t j = sortindex[i]; size_t pid = m_phaseID[j]; plogf(" %-12.12s", m_speciesName[j]); - plogf(" %14.7E ", m_molNumSpecies_old[j]*molScale); + plogf(" %14.7E ", m_molNumSpecies_old[j]); plogf("%14.7E ", m_SSfeSpecies[j]); plogf("%14.7E ", log(m_actCoeffSpecies_old[j])); double tpmoles = m_tPhaseMoles_old[pid]; @@ -291,7 +283,7 @@ int VCS_SOLVE::vcs_report(int iconv) plogf(" "); } - plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j] * molScale); + plogf("| %20.9E |", m_feSpecies_old[j] * m_molNumSpecies_old[j]); plogf("\n"); } for (size_t i = 0; i < 125; i++) { diff --git a/src/equil/vcs_solve.cpp b/src/equil/vcs_solve.cpp index 8171656d7..3477081bd 100644 --- a/src/equil/vcs_solve.cpp +++ b/src/equil/vcs_solve.cpp @@ -41,7 +41,6 @@ VCS_SOLVE::VCS_SOLVE(MultiPhase* mphase, int printLvl) : m_tolmaj2(1.0E-10), m_tolmin2(1.0E-8), m_unitsState(VCS_DIMENSIONAL_G), - m_totalMoleScale(1.0), m_useActCoeffJac(0), m_totalVol(mphase->volume()), m_Faraday_dim(ElectronCharge * Avogadro), diff --git a/src/equil/vcs_solve_TP.cpp b/src/equil/vcs_solve_TP.cpp index 16641fa9d..d3b7b5318 100644 --- a/src/equil/vcs_solve_TP.cpp +++ b/src/equil/vcs_solve_TP.cpp @@ -942,9 +942,6 @@ void VCS_SOLVE::solve_tp_inner(size_t& iti, size_t& it1, } else { plogf(" (only major species):"); } - if (m_totalMoleScale != 1.0) { - plogf(" (Total Mole Scale = %g)", m_totalMoleScale); - } plogf("\n"); plogf(" --- Species Status Initial_KMoles Final_KMoles Initial_Mu/RT"); plogf(" Mu/RT Init_Del_G/RT Delta_G/RT\n");