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

@@ -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) {
+       * Scale the total moles if necessary:
-	    //m_molNumSpecies_old[i] *= 1.0E3;
+       *  First find out the total moles
-	    m_molNumSpecies_old[i] *= 1.0;
+       */
+      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) {
+	vcs_tmoles();
-	  //m_elemAbundancesGoal[i] *= 1.0E3;
-	  m_elemAbundancesGoal[i] *= 1.0;
-	}
       }
     }
   }
@@ -166,19 +224,33 @@ namespace VCSnonideal {
       } 
       m_Faraday_dim *= tf;
     }
-    if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
+    if (m_totalMoleScale != 1.0) {
-      for (i = 0; i < m_numSpeciesTot; ++i) {
+      if (m_VCS_UnitsFormat == VCS_UNITS_MKS) {
-	if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
+#ifdef DEBUG_MODE
-	  //m_molNumSpecies_old[i] /= 1.0E3;
+	if (m_debug_print_lvl >= 2) {
-	  m_molNumSpecies_old[i] /= 1.0;
+	  plogf("  --- vcs_redim_TP() called: getting rid of mole scale of %g", m_totalMoleScale);
+	  plogendl();
 	}
-      }
+#endif
-      for (i = 0; i < m_numElemConstraints; ++i) {
+	for (i = 0; i < m_numSpeciesTot; ++i) {
-	//m_elemAbundancesGoal[i] /= 1.0E3;
+	  if (m_speciesUnknownType[i] != VCS_SPECIES_TYPE_INTERFACIALVOLTAGE) {
-	m_elemAbundancesGoal[i] /= 1.0;
+	    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

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.
-     *  vcs_report:
+   *     This prints out the current contents of the VCS_SOLVE class, V.
-     *
+   *     The "old" solution vector is printed out.
-     *     Print out a report on the state of the equilibrium problem to 
+   ***************************************************************************/
-     *     standard output.
+  int VCS_SOLVE::vcs_report(int iconv) {
-     *     This prints out the current contents of the VCS_SOLVE class, V.
+    bool printActualMoles = true;
-     *     The "old" solution vector is printed out.
-     ***************************************************************************/
-  {
     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\tTemperature  = %15.2g Kelvin\n", m_temperature);
-    plogf("\t\tPressure    = %15.5g Pa \n", m_pressurePA); 
+    plogf("\t\tPressure     = %15.5g Pa \n", m_pressurePA); 
-    plogf("\t\tVolume      = %15.5g m**3\n", m_totalVol);
+    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], 
+      plogf("%14.7E     %14.7E    %12.4E", m_molNumSpecies_old[i] * molScale, 
-	    m_molNumSpecies_new[i], m_feSpecies_old[i]);
+	    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], 
+	plogf("%14.7E     %14.7E    %12.4E", m_molNumSpecies_old[l] * molScale, 
-	      m_molNumSpecies_new[l], m_feSpecies_old[l]);
+	      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");
     }

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),

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
   /*!

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;
   }
   /*****************************************************************************/