From a441bbf272e6d06d2b93ae0d0a3121e736766fa1 Mon Sep 17 00:00:00 2001 From: Harry Moffat Date: Tue, 29 Jul 2008 16:23:29 +0000 Subject: [PATCH] Put in a scale factor for mole number within the nondimensionalization routines. Now, the code can solve equilibrium problems for virtually any range of total mole numbers. --- Cantera/src/equil/vcs_nondim.cpp | 112 +++++++++++++++++++++++------ Cantera/src/equil/vcs_report.cpp | 71 +++++++++--------- Cantera/src/equil/vcs_solve.cpp | 1 + Cantera/src/equil/vcs_solve.h | 11 ++- Cantera/src/equil/vcs_solve_TP.cpp | 11 ++- 5 files changed, 150 insertions(+), 56 deletions(-) diff --git a/Cantera/src/equil/vcs_nondim.cpp b/Cantera/src/equil/vcs_nondim.cpp index 0f136088b..2b6d6e99e 100644 --- a/Cantera/src/equil/vcs_nondim.cpp +++ b/Cantera/src/equil/vcs_nondim.cpp @@ -16,6 +16,8 @@ #include "vcs_solve.h" #include "vcs_internal.h" +#include "vcs_VolPhase.h" +#include "stringUtils.h" namespace VCSnonideal { @@ -82,7 +84,7 @@ namespace VCSnonideal { default: plogf("vcs_nondimMult_TP error: unknown units: %d\n", mu_units); plogendl(); - exit(-1); + std::exit(-1); } return rt; } @@ -122,17 +124,73 @@ namespace VCSnonideal { } m_Faraday_dim = vcs_nondim_Farad(m_VCS_UnitsFormat, m_temperature); - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - for (i = 0; i < m_numSpeciesTot; ++i) { - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - //m_molNumSpecies_old[i] *= 1.0E3; - m_molNumSpecies_old[i] *= 1.0; + + /* + * 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 (i = 0; i < m_numElemConstraints; ++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) { + plogf(" VCS_SOLVE::vcs_nondim_TP ERROR: Total input moles , %g, is outside the range handled by vcs. exit", + tmole_orig); + plogendl(); + throw vcsError("VCS_SOLVE::vcs_nondim_TP", " Total input moles ," + Cantera::fp2str(tmole_orig) + + "is outside the range handled by vcs.\n"); + } + + // 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_VCS_UnitsFormat == VCS_UNITS_MKS) { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_nondim_TP() called: USING A MOLE SCALE OF %g until further notice", m_totalMoleScale); + plogendl(); + } +#endif + for (i = 0; i < m_numSpeciesTot; ++i) { + if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + m_molNumSpecies_old[i] *= (1.0 / m_totalMoleScale); + } + } + for (i = 0; i < m_numElemConstraints; ++i) { + m_elemAbundancesGoal[i] *= (1.0 / m_totalMoleScale); + } + + for (int iph = 0; iph < m_numPhases; iph++) { + TPhInertMoles[iph] *= (1.0 / m_totalMoleScale); + if (TPhInertMoles[iph] != 0.0) { + vcs_VolPhase *vphase = m_VolPhaseList[iph]; + vphase->setTotalMolesInert(TPhInertMoles[iph]); + } } } - for (i = 0; i < m_numElemConstraints; ++i) { - //m_elemAbundancesGoal[i] *= 1.0E3; - m_elemAbundancesGoal[i] *= 1.0; - } + vcs_tmoles(); } } } @@ -166,19 +224,33 @@ namespace VCSnonideal { } m_Faraday_dim *= tf; } - if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { - for (i = 0; i < m_numSpeciesTot; ++i) { - if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { - //m_molNumSpecies_old[i] /= 1.0E3; - m_molNumSpecies_old[i] /= 1.0; + if (m_totalMoleScale != 1.0) { + if (m_VCS_UnitsFormat == VCS_UNITS_MKS) { +#ifdef DEBUG_MODE + if (m_debug_print_lvl >= 2) { + plogf(" --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale); + plogendl(); } - } - for (i = 0; i < m_numElemConstraints; ++i) { - //m_elemAbundancesGoal[i] /= 1.0E3; - m_elemAbundancesGoal[i] /= 1.0; +#endif + for (i = 0; i < m_numSpeciesTot; ++i) { + if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { + m_molNumSpecies_old[i] *= m_totalMoleScale; + } + } + for (i = 0; i < m_numElemConstraints; ++i) { + m_elemAbundancesGoal[i] *= m_totalMoleScale; + } + + for (int iph = 0; iph < m_numPhases; iph++) { + TPhInertMoles[iph] *= m_totalMoleScale; + if (TPhInertMoles[iph] != 0.0) { + vcs_VolPhase *vphase = m_VolPhaseList[iph]; + vphase->setTotalMolesInert(TPhInertMoles[iph]); + } + } + vcs_tmoles(); } } - } // Computes the current elemental abundances vector diff --git a/Cantera/src/equil/vcs_report.cpp b/Cantera/src/equil/vcs_report.cpp index 6c8ff892f..053ade6cc 100644 --- a/Cantera/src/equil/vcs_report.cpp +++ b/Cantera/src/equil/vcs_report.cpp @@ -32,21 +32,18 @@ namespace VCSnonideal { } /*****************************************************************************/ - /*****************************************************************************/ - /*****************************************************************************/ - - int VCS_SOLVE::vcs_report(int iconv) - - /************************************************************************** - * - * 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. - ***************************************************************************/ - { + + /************************************************************************** + * + * 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; int i, j, l, k, inertYes = FALSE, kspec; int nspecies = m_numSpeciesTot; double g; @@ -86,6 +83,11 @@ namespace VCSnonideal { 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_numSpeciesTot); /* ******************************************************** */ @@ -110,9 +112,13 @@ namespace VCSnonideal { m_totalVol = vcs_VolTotal(m_temperature, m_pressurePA, VCS_DATA_PTR(m_molNumSpecies_old), VCS_DATA_PTR(m_PMVolumeSpecies)); - plogf("\t\tTemperature = %15.2g Kelvin\n", m_temperature); - plogf("\t\tPressure = %15.5g Pa \n", m_pressurePA); - plogf("\t\tVolume = %15.5g m**3\n", m_totalVol); + 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); + } /* * -------- TABLE OF SPECIES IN DECREASING MOLE NUMBERS -------------- @@ -125,8 +131,8 @@ namespace VCSnonideal { for (i = 0; i < m_numComponents; ++i) { plogf(" %-12.12s", m_speciesName[i].c_str()); print_space(13); - plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i], - m_molNumSpecies_new[i], m_feSpecies_old[i]); + plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[i] * molScale, + m_molNumSpecies_new[i] * molScale, m_feSpecies_old[i]); plogf(" %3d", m_speciesUnknownType[i]); plogf("\n"); } @@ -136,15 +142,15 @@ namespace VCSnonideal { print_space(13); if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_MOLNUM) { - plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[l], - m_molNumSpecies_new[l], m_feSpecies_old[l]); + plogf("%14.7E %14.7E %12.4E", m_molNumSpecies_old[l] * molScale, + m_molNumSpecies_new[l] * molScale, m_feSpecies_old[l]); plogf(" KMolNum "); } else if (m_speciesUnknownType[l] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { plogf(" NA %14.7E %12.4E", 1.0, m_feSpecies_old[l]); - plogf(" Voltage = %14.7E", m_molNumSpecies_old[l]); + plogf(" Voltage = %14.7E", m_molNumSpecies_old[l] * molScale); } else { plogf("we have a problem\n"); - exit(-1); + std::exit(-1); } plogf("\n"); } @@ -157,7 +163,7 @@ namespace VCSnonideal { plogf(" Inert Species in phase %16s ", (m_VolPhaseList[i])->PhaseName.c_str()); } - plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i], + plogf("%14.7E %14.7E %12.4E\n", TPhInertMoles[i] * molScale, TPhInertMoles[i] / m_tPhaseMoles_old[i], 0.0); } } @@ -167,7 +173,8 @@ namespace VCSnonideal { plogf(" %-12.12s", m_speciesName[kspec].c_str()); // Note m_deltaGRxn_new[] stores in kspec slot not irxn slot, after solve plogf(" %14.7E %14.7E %12.4E", - m_molNumSpecies_old[kspec], m_molNumSpecies_new[kspec], m_deltaGRxn_new[kspec]); + m_molNumSpecies_old[kspec]*molScale, + m_molNumSpecies_new[kspec]*molScale, m_deltaGRxn_new[kspec]); if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_MOLNUM) { plogf(" KMol_Num"); } else if (m_speciesUnknownType[i] == VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) { @@ -199,7 +206,7 @@ namespace VCSnonideal { plogf(" | |\n"); plogf(" NonComponent | Moles |"); for (j = 0; j < m_numComponents; j++) { - plogf(" %10.3g", m_molNumSpecies_old[j]); + plogf(" %10.3g", m_molNumSpecies_old[j] * molScale); } plogf(" | DG/RT Rxn |\n"); print_line("-", m_numComponents*10 + 45); @@ -207,7 +214,7 @@ namespace VCSnonideal { int kspec = m_indexRxnToSpecies[irxn]; plogf(" %3d ", kspec); plogf("%-10.10s", m_speciesName[kspec].c_str()); - plogf("|%10.3g |", m_molNumSpecies_old[kspec]); + plogf("|%10.3g |", m_molNumSpecies_old[kspec]*molScale); for (j = 0; j < m_numComponents; j++) { plogf(" %6.2f", m_stoichCoeffRxnMatrix[irxn][j]); } @@ -248,7 +255,7 @@ namespace VCSnonideal { plogf(" %3d ", iphase); vcs_VolPhase *VPhase = m_VolPhaseList[iphase]; plogf("%-12.12s |",VPhase->PhaseName.c_str()); - plogf("%10.3e |", m_tPhaseMoles_old[iphase]); + plogf("%10.3e |", m_tPhaseMoles_old[iphase]*molScale); totalMoles += m_tPhaseMoles_old[iphase]; if (m_tPhaseMoles_old[iphase] != VPhase->TotalMoles()) { if (! vcs_doubleEqual(m_tPhaseMoles_old[iphase], VPhase->TotalMoles())) { @@ -296,7 +303,7 @@ namespace VCSnonideal { plogf(" Actual Target Type ElActive\n"); for (i = 0; i < m_numElemConstraints; ++i) { print_space(26); plogf("%-2.2s", (m_elementName[i]).c_str()); - plogf("%20.12E %20.12E", m_elemAbundances[i], m_elemAbundancesGoal[i]); + plogf("%20.12E %20.12E", m_elemAbundances[i]*molScale, m_elemAbundancesGoal[i]*molScale); plogf(" %3d %3d\n", m_elType[i], m_elementActive[i]); } plogf("\n"); @@ -322,7 +329,7 @@ namespace VCSnonideal { l = sortindex[i]; int pid = m_phaseID[l]; plogf(" %-12.12s", m_speciesName[l].c_str()); - plogf(" %14.7E ", m_molNumSpecies_old[l]); + plogf(" %14.7E ", m_molNumSpecies_old[l]*molScale); plogf("%14.7E ", m_SSfeSpecies[l]); plogf("%14.7E ", log(m_actCoeffSpecies_old[l])); double tpmoles = m_tPhaseMoles_old[pid]; @@ -356,7 +363,7 @@ namespace VCSnonideal { } #ifdef DEBUG_MODE - plogf("| %20.13E |", m_feSpecies_old[l] * m_molNumSpecies_old[l]); + plogf("| %20.13E |", m_feSpecies_old[l] * m_molNumSpecies_old[l] * molScale); #endif plogf("\n"); } diff --git a/Cantera/src/equil/vcs_solve.cpp b/Cantera/src/equil/vcs_solve.cpp index 45f17535c..3ce0a2ee6 100644 --- a/Cantera/src/equil/vcs_solve.cpp +++ b/Cantera/src/equil/vcs_solve.cpp @@ -50,6 +50,7 @@ namespace VCSnonideal { m_tolmaj2(0.0), m_tolmin2(0.0), m_unitsState(VCS_DIMENSIONAL_G), + m_totalMoleScale(1.0), m_useActCoeffJac(0), m_totalVol(0.0), m_Faraday_dim(1.602e-19 * 6.022136736e26), diff --git a/Cantera/src/equil/vcs_solve.h b/Cantera/src/equil/vcs_solve.h index 58e9d897d..f010fcf42 100644 --- a/Cantera/src/equil/vcs_solve.h +++ b/Cantera/src/equil/vcs_solve.h @@ -506,7 +506,7 @@ public: * the variable m_totalMolNum. * Reconciles Phase existence flags with total moles in each phase. */ - void vcs_tmoles(); + double vcs_tmoles(); //! This subroutine calculates reaction free energy changes for @@ -1758,6 +1758,15 @@ public: *. The default is to have this unitless */ char m_unitsState; + + //! Multiplier for the mole numbers within the nondimensionless 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 /*! diff --git a/Cantera/src/equil/vcs_solve_TP.cpp b/Cantera/src/equil/vcs_solve_TP.cpp index e518846c4..449e8fd0a 100644 --- a/Cantera/src/equil/vcs_solve_TP.cpp +++ b/Cantera/src/equil/vcs_solve_TP.cpp @@ -1157,10 +1157,14 @@ namespace VCSnonideal { plogf(" "); vcs_print_line("-", 103); plogf(" --- Summary of the Update "); if (iti == 0) { - plogf(" (all species):\n"); + plogf(" (all species):"); } else { - plogf(" (only major species):\n"); + 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"); for (i = 0; i < m_numComponents; ++i) { @@ -4936,7 +4940,7 @@ namespace VCSnonideal { * Calculates the total number of moles in all phases. * Reconciles Phase existence flags with total moles in each phase. */ - void VCS_SOLVE::vcs_tmoles() { + double VCS_SOLVE::vcs_tmoles() { int i; double sum; vcs_VolPhase *Vphase; @@ -4961,6 +4965,7 @@ namespace VCSnonideal { } } m_totalMolNum = sum; + return m_totalMolNum; } /*****************************************************************************/