Changes to exit() and additions of cstdlib

Cantera/src/thermo/DebyeHuckel.cpp

@@ -26,6 +26,7 @@
 #include "WaterProps.h"
 #include "PDSS_Water.h"
 #include <cstring>
+#include <cstdlib>
 
 using namespace std;
 
@@ -997,11 +998,11 @@ namespace Cantera {
 	if (formString == "unity") {
 	  m_formGC = 0;
 	  printf("exit standardConc = unity not done\n");
-	  exit(-1);
+	  exit(EXIT_FAILURE);
 	} else if (formString == "molar_volume") {
 	  m_formGC = 1;
 	  printf("exit standardConc = molar_volume not done\n");
-	  exit(-1);
+	  exit(EXIT_FAILURE);
 	} else if (formString == "solvent_volume") {
 	  m_formGC = 2;
 	} else {
@@ -1113,11 +1114,11 @@ namespace Cantera {
 	if (formString == "unity") {
 	  m_formGC = 0;
 	  printf("exit standardConc = unity not done\n");
-	  exit(-1);
+	  exit(EXIT_FAILURE);
 	} else if (formString == "molar_volume") {
 	  m_formGC = 1;
 	  printf("exit standardConc = molar_volume not done\n");
-	  exit(-1);
+	  exit(EXIT_FAILURE);
 	} else if (formString == "solvent_volume") {
 	  m_formGC = 2;
 	} else {
@@ -1682,7 +1683,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
     }
     return A;
   }
@@ -1716,7 +1717,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
     }
     return dAdT;
   }
@@ -1750,7 +1751,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
     }
     return d2AdT2;
   }
@@ -1784,7 +1785,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
     }
     return dAdP;
   }
@@ -2176,7 +2177,7 @@ namespace Cantera {
 
     default:
       printf("ERROR\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
     }
     /*
      * Above, we calculated the ln(activitySolvent). Translate that
@@ -2321,7 +2322,7 @@ namespace Cantera {
 
     default:
       printf("ERROR\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
       break;
     }
 
@@ -2458,7 +2459,7 @@ namespace Cantera {
 
     default:
       printf("ERROR\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
       break;
     }
   }
@@ -2596,7 +2597,7 @@ namespace Cantera {
 
     default:
       printf("ERROR\n");
-      exit(-1);
+      exit(EXIT_FAILURE);
       break;
     }
   }

Cantera/src/thermo/HMWSoln.cpp

@@ -29,7 +29,8 @@
 #include "ThermoFactory.h"
 #include "WaterProps.h"
 #include "PDSS_Water.h"
-#include <math.h>
+#include <cmath>
+#include <cstdlib>
 
 namespace Cantera {
 
@@ -389,7 +390,7 @@ namespace Cantera {
   {
     if (testProb != 1) {
       printf("unknown test problem\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 
     constructPhaseFile("HMW_NaCl.xml", "");
@@ -1324,7 +1325,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
     return A;
   }
@@ -1359,7 +1360,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
     return dAdT;
   }
@@ -1393,7 +1394,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
     return dAdP;
   }
@@ -1500,7 +1501,7 @@ namespace Cantera {
       break;
     default:
       printf("shouldn't be here\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
     return d2AdT2;
   }
@@ -2175,7 +2176,7 @@ namespace Cantera {
      */
     if (m_indexSolvent != 0) {
       printf("Wrong index solvent value!\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 
 #ifdef DEBUG_MODE
@@ -2979,7 +2980,7 @@ namespace Cantera {
 	  if (j == (m_kk-1)) {
 	    // we should never reach this step
 	    printf("logic error 1 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] > 0.0) {
 	    /*
@@ -3009,7 +3010,7 @@ namespace Cantera {
 	  if (j == m_kk-1) {
 	    // we should never reach this step
 	    printf("logic error 2 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] < 0) {
 	    /*
@@ -3149,7 +3150,7 @@ namespace Cantera {
 #endif
     if (m_indexSolvent != 0) {
       printf("Wrong index solvent value!\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 
     double d_wateract_dT;
@@ -3787,7 +3788,7 @@ namespace Cantera {
 	  if (j == (m_kk-1)) {
 	    // we should never reach this step
 	    printf("logic error 1 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] > 0.0) {
 	    /*
@@ -3817,7 +3818,7 @@ namespace Cantera {
 	  if (j == m_kk-1) {
 	    // we should never reach this step
 	    printf("logic error 2 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] < 0) {
 	    /*
@@ -3938,7 +3939,7 @@ namespace Cantera {
 #endif
     if (m_indexSolvent != 0) {
       printf("Wrong index solvent value!\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 
     std::string sni, snj, snk; 
@@ -4583,7 +4584,7 @@ namespace Cantera {
 	  if (j == (m_kk-1)) {
 	    // we should never reach this step
 	    printf("logic error 1 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] > 0.0) {
 	    /*
@@ -4613,7 +4614,7 @@ namespace Cantera {
 	  if (j == m_kk-1) {
 	    // we should never reach this step
 	    printf("logic error 2 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] < 0) {
 	    /*
@@ -4754,7 +4755,7 @@ namespace Cantera {
 #endif
     if (m_indexSolvent != 0) {
       printf("Wrong index solvent value!\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 
     double d_wateract_dP;
@@ -5397,7 +5398,7 @@ namespace Cantera {
 	  if (j == (m_kk-1)) {
 	    // we should never reach this step
 	    printf("logic error 1 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] > 0.0) {
 	    /*
@@ -5428,7 +5429,7 @@ namespace Cantera {
 	  if (j == m_kk-1) {
 	    // we should never reach this step
 	    printf("logic error 2 in Step 9 of hmw_act");
-	    std::exit(1);
+	    exit(EXIT_FAILURE);
 	  }
 	  if (charge[k] < 0) {
 	    /*
@@ -5614,13 +5615,13 @@ namespace Cantera {
 #ifdef DEBUG_MODE
     if (i > 4 || j > 4) {
       printf("we shouldn't be here\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 #endif
 
     if ((i == 0) || (j == 0)) {
       printf("ERROR calc_thetas called with one species being neutral\n");
-      std::exit(-1);
+      exit(EXIT_FAILURE);
     }
 
     /*

Cantera/src/thermo/HMWSoln_input.cpp

@@ -21,10 +21,10 @@
 #include "WaterProps.h"
 #include "PDSS_Water.h"
 #include <cstring>
+#include <cstdlib>
 
 using namespace std;
 
-
 namespace Cantera {
 
 
@@ -1011,11 +1011,11 @@ namespace Cantera {
 	if (formString == "unity") {
 	  m_formGC = 0;
 	  printf("exit standardConc = unity not done\n");
-	  exit(-1);
+	  exit(EXIT_FAILURE);
 	} else if (formString == "molar_volume") {
 	  m_formGC = 1;
 	  printf("exit standardConc = molar_volume not done\n");
-	  exit(-1);
+	  exit(EXIT_FAILURE);
 	} else if (formString == "solvent_volume") {
 	  m_formGC = 2;
 	} else {

Cantera/src/thermo/MolalityVPSSTP.cpp

@@ -41,6 +41,8 @@ namespace Cantera {
   MolalityVPSSTP::MolalityVPSSTP() :
     VPStandardStateTP(),
     m_indexSolvent(0),
+    m_pHScalingType(PHSCALE_PITZER),
+    m_indexCLM(-1),
     m_weightSolvent(18.01528),
     m_xmolSolventMIN(0.01),
     m_Mnaught(18.01528E-3)
@@ -62,6 +64,8 @@ namespace Cantera {
   MolalityVPSSTP::MolalityVPSSTP(const MolalityVPSSTP &b) :
     VPStandardStateTP(),
     m_indexSolvent(b.m_indexSolvent),
+    m_pHScalingType(b.m_pHScalingType),
+    m_indexCLM(b.m_indexCLM),
     m_xmolSolventMIN(b.m_xmolSolventMIN),
     m_Mnaught(b.m_Mnaught),
     m_molalities(b.m_molalities)
@@ -80,6 +84,8 @@ namespace Cantera {
     if (&b != this) {
       VPStandardStateTP::operator=(b);
       m_indexSolvent     = b.m_indexSolvent;
+      m_pHScalingType    = b.m_pHScalingType;
+      m_indexCLM         = b.m_indexCLM;
       m_weightSolvent    = b.m_weightSolvent;
       m_xmolSolventMIN   = b.m_xmolSolventMIN;
       m_Mnaught          = b.m_Mnaught;
@@ -391,6 +397,24 @@ namespace Cantera {
     return cAC_CONVENTION_MOLALITY;
   }
 
+  void MolalityVPSSTP::getActivityConcentrations(doublereal* c) const {
+      err("getActivityConcentrations");
+  }
+
+  doublereal MolalityVPSSTP::standardConcentration(int k) const {
+    err("standardConcentration");
+    return -1.0;
+  }
+
+  doublereal MolalityVPSSTP::logStandardConc(int k) const {
+    err("logStandardConc");
+    return -1.0;
+  }
+
+  void MolalityVPSSTP::getActivities(doublereal* ac) const {
+    err("getActivities");
+  }
+
   /*
    * Get the array of non-dimensional activity coefficients at
    * the current solution temperature, pressure, and
@@ -416,7 +440,12 @@ namespace Cantera {
       }
     }
   }
-    
+
+
+  void MolalityVPSSTP::getMolalityActivityCoefficients(doublereal *acMolality) const {
+    err("getMolalityActivityCoefficients");
+  }
+
   /*
    * osmotic coefficient:
    * 
@@ -452,6 +481,15 @@ namespace Cantera {
     return oc;
   }
 
+
+  void MolalityVPSSTP::getElectrochemPotentials(doublereal* mu) const {
+    getChemPotentials(mu);
+    double ve = Faraday * electricPotential();
+    for (int k = 0; k < m_kk; k++) {
+      mu[k] += ve*charge(k);
+    }
+  }
+  
   /*
    * ------------ Partial Molar Properties of the Solution ------------
    */
@@ -495,7 +533,12 @@ namespace Cantera {
       if (i == 5) uA[5] = 0.0;
     }
   }
-    
+
+  void MolalityVPSSTP::setToEquilState(const doublereal* lambda_RT) {
+    updateStandardStateThermo();
+    err("setToEquilState");
+  }
+
   /*
    * Set the thermodynamic state.
    */
@@ -561,6 +604,77 @@ namespace Cantera {
     setSolvent(0);
   }
 
+  // Returns the index of the Cl- species.
+  /*
+   *  The Cl- species is special in the sense that it's single ion
+   *  molalality-based activity coefficient is used in the specification
+   *  of the pH scale for single ions. Therefore, we need to know
+   *  what species index is Cl-. If the species isn't in the species
+   *  list then this routine returns -1, and we can't use the NBS
+   *  pH scale. 
+   *   
+   *  Right now we use a restrictive interpretation. The species
+   *  must be named "Cl-". It must consist of exactly one Cl and one E 
+   *  atom. 
+   */
+  int MolalityVPSSTP::findCLMIndex() const {
+    int indexCLM = -1;
+    int eCl = -1;
+    int eE = -1;
+    int ne= nElements();
+    string sn;
+    for (int e = 0; e < ne; e++) {
+      sn = elementName(e);
+      if (sn == "Cl" || sn == "CL") {
+	eCl = e;
+	break;
+      }
+    }
+    // We have failed if we can't find the Cl element index
+    if (eCl == -1) {
+      return -1;
+    }
+    for (int e = 0; e < ne; e++) {
+      sn = elementName(e);
+      if (sn == "E" || sn == "e") {
+	eE = e;
+	break;
+      }
+    }
+    // We have failed if we can't find the E element index
+    if (eE == -1) {
+      return -1;
+    }
+    for (int k = 1; k < m_kk; k++) {
+      doublereal nCl = nAtoms(k, eCl);
+      if (nCl != 1.0) {
+	continue;
+      }
+      doublereal nE = nAtoms(k, eE);
+      if (nE != -1.0) {
+	continue;
+      }
+      for (int e = 0; e < ne; e++) {
+	if (e != eE && e != eCl) {
+	  doublereal nA = nAtoms(k, e);
+	  if (nA != 0.0) {
+	    continue;
+	  }
+	}
+      }
+      sn = speciesName(k);
+      if (sn != "Cl-" && sn != "CL-") {
+	continue;
+      }
+
+      indexCLM = k;
+      break;
+    }
+    return indexCLM;
+  }
+
+  //   Initialize lengths of local variables after all species have
+  //   been identified.
   void  MolalityVPSSTP::initLengths() {
     m_kk = nSpecies();
     m_molalities.resize(m_kk);

Cantera/src/thermo/MolalityVPSSTP.h

@@ -138,8 +138,8 @@ namespace Cantera {
    *  return the molality-based quantities. Also all functions which return 
    *  activities return the molality-based activities. The reason for this convention
    *  has been discussed in supporting memos. However, it's important because the
-   *  term in the equation above is non-trivial. For example it's equal to 2.38 kcal gmol<SUP>-1</SUP>
+   *  term in the equation above is non-trivial. For example it's equal
-   *  for water at 298 K. 
+   *  to 2.38 kcal gmol<SUP>-1</SUP> for water at 298 K. 
    *
    *
    * In order to prevent a singularity, this class includes the concept of a minimum
@@ -222,6 +222,14 @@ namespace Cantera {
     virtual int eosType() const { return 0; }
 
 
+    //! Set the pH scale, which determines the scale for single-ion activity 
+    //! coefficients.
+    /*!
+     *  Single ion activity coefficients are not unique in terms of the
+     *  representing actual measureable quantities. 
+     */
+    void setpHScale(const int pHscaleType);
+
     /**
      * @} 
      * @name  Molar Thermodynamic Properties 
@@ -417,9 +425,7 @@ namespace Cantera {
      *           units depend upon the implementation of the
      *           reaction rate expressions within the phase.
      */
-    virtual void getActivityConcentrations(doublereal* c) const {
+    virtual void getActivityConcentrations(doublereal* c) const;
-      err("getActivityConcentrations");
-    }
 
     /**
      * The standard concentration \f$ C^0_k \f$ used to normalize
@@ -434,10 +440,7 @@ namespace Cantera {
      *
      * @param k species index. Defaults to zero.
      */
-    virtual doublereal standardConcentration(int k=0) const {
+    virtual doublereal standardConcentration(int k=0) const;
-      err("standardConcentration");
-      return -1.0;
-    }
 
     /**
      * Returns the natural logarithm of the standard 
@@ -445,10 +448,7 @@ namespace Cantera {
      *
      * @param k  species index
      */
-    virtual doublereal logStandardConc(int k=0) const {
+    virtual doublereal logStandardConc(int k=0) const;
-      err("logStandardConc");
-      return -1.0;
-    }
 
     /**
      * Returns the units of the standard and generalized
@@ -492,9 +492,7 @@ namespace Cantera {
      *
      * @param ac     Output vector of molality-based activities. Length: m_kk.
      */
-    virtual void getActivities(doublereal* ac) const {
+    virtual void getActivities(doublereal* ac) const;
-      err("getActivities");
-    }
 
     //! Get the array of non-dimensional activity coefficients at
     //! the current solution temperature, pressure, and solution concentration.
@@ -540,9 +538,7 @@ namespace Cantera {
      * @param acMolality Output vector containing the molality based activity coefficients.
      *                   length: m_kk.
      */
-    virtual void getMolalityActivityCoefficients(doublereal *acMolality) const {
+    virtual void getMolalityActivityCoefficients(doublereal *acMolality) const;
-      err("getMolalityActivityCoefficients");
-    }
    
     //! Calculate the osmotic coefficient
     /*!
@@ -576,13 +572,7 @@ namespace Cantera {
      * @param mu     output vector containing the species electrochemical potentials.
      *               Length: m_kk.
      */
-    void getElectrochemPotentials(doublereal* mu) const {
+    void getElectrochemPotentials(doublereal* mu) const;
-      getChemPotentials(mu);
-      double ve = Faraday * electricPotential();
-      for (int k = 0; k < m_kk; k++) {
-	mu[k] += ve*charge(k);
-      }
-    }
 
  
     //@}
@@ -640,10 +630,7 @@ namespace Cantera {
      * @param lambda_RT Input vector containing the dimensionless 
      *                  element potentials.
      */
-    virtual void setToEquilState(const doublereal* lambda_RT) {
+    virtual void setToEquilState(const doublereal* lambda_RT);
-      updateStandardStateThermo();
-      err("setToEquilState");
-    }
 
 
     //@}
@@ -753,6 +740,23 @@ namespace Cantera {
     virtual std::string report(bool show_thermo = true) const;
 
   private:
+    //! Returns the index of the Cl- species.
+    /*!
+     *  The Cl- species is special in the sense that it's single ion
+     *  molalality-based activity coefficient is used in the specification
+     *  of the pH scale for single ions. Therefore, we need to know
+     *  what species index is Cl-. If the species isn't in the species
+     *  list then this routine returns -1, and we can't use the NBS
+     *  pH scale. 
+     *    
+     *  Right now we use a restrictive interpretation. The species
+     *  must be named "Cl-". It must consist of exactly one Cl and one E 
+     *  atom. 
+     */
+    virtual int findCLMIndex() const;
+
+    //! Initialize lengths of local variables after all species have
+    //! been identified.
     void initLengths();
             
   protected:
@@ -763,6 +767,24 @@ namespace Cantera {
      */
     int        m_indexSolvent;
 
+    //! Scaling to be used for output of single-ion species activity
+    //! coefficients.
+    /*!
+     *   Index of the species to be used in the single-ion scaling
+     *   law. This is the indentity of the Cl- species for the PHSCALE_NBS
+     *   scaling.
+     *   Either PHSCALE_PITZER or PHSCALE_NBS
+     */
+    int        m_pHScalingType;
+
+    //! Index of the phScale species
+    /*!
+     *   Index of the species to be used in the single-ion scaling
+     *   law. This is the indentity of the Cl- species for the PHSCALE_NBS
+     *   scaling
+     */
+    int        m_indexCLM;
+
     //! Molecular weight of the Solvent
     doublereal m_weightSolvent;
 
@@ -790,11 +812,68 @@ namespace Cantera {
      * units are (kg/kmol)
      */
     mutable vector_fp  m_molalities;
+
   private:
+    //! Error function
+    /*!
+     *  Print an error string and exit
+     *
+     * @param msg  Message to be printed
+     */
     doublereal err(std::string msg) const;
 
   };
 
+
+
+  //! Scale to be used for the output of single-ion activity coefficients
+  //! is that used by Pitzer.
+  /*!
+   *  This is the internal scale used within the code. One property is that
+   *  the activity coefficients for the cation and anion of a single salt
+   *  will be equal. This scale is the one presumed by the formulation of the
+   *  single-ion activity coefficients described in this report.
+   *
+   *  Activity coefficients for species k may be altered between scales s1 to s2
+   *  using the following formula
+   *
+   *   \f[
+   *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) 
+   *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
+   *   \f]
+   *
+   *  where j is any one species.
+   *
+   *
+   */
+  const int PHSCALE_PITZER = 0;
+
+  //! Scale to be used for evaluation of single-ion activity coefficients
+  //! is that used by the NBS standard for evaluation of the pH variable.
+  /*!
+   *  This is not the internal scale used within the code.
+   *
+   *  Activity coefficients for species k may be altered between scales s1 to s2
+   *  using the following formula
+   *
+   *   \f[
+   *       ln(\gamma_k^{s2}) = ln(\gamma_k^{s1}) 
+   *          + \frac{z_k}{z_j} \left(  ln(\gamma_j^{s2}) - ln(\gamma_j^{s1}) \right)
+   *   \f]
+   *
+   *  where j is any one species. For the NBS scale, j is equal to the Cl- species
+   *  and 
+   *
+   *  \f[
+   *       ln(\gamma_{Cl-}^{s2}) = \frac{-A_{\phi} \sqrt{I}}{1.0 + 1.5 \sqrt{I}}
+   *  \f]
+   *
+   *  This is the NBS pH scale, which is used in all conventional pH 
+   *  measurements. and is based on the Bates-Guggenheim quations.
+   *
+   */
+  const int PHSCALE_NBS    = 1;
+
 }
         
 #endif

Cantera/src/thermo/WaterPropsIAPWS.cpp

@@ -16,9 +16,9 @@
 #include "WaterPropsIAPWS.h"
 #include "ctexceptions.h"
 #include "stringUtils.h"
-#include <math.h>
+#include <cmath>
-#include <stdio.h>
+#include <cstdio>
-#include <stdlib.h>
+#include <cstdlib>
 /*
  * Critical Point values of water in mks units
  */

Cantera/src/thermo/WaterPropsIAPWSphi.cpp

@@ -14,8 +14,8 @@
 
 #include "WaterPropsIAPWSphi.h"
 
-#include <math.h>
+#include <cmath>
-#include <stdio.h>
+#include <cstdio>
 
 /*
  * Critical Point values in mks units: Note, these aren't used in this 

Cantera/src/thermo/phasereport.cpp

@@ -19,7 +19,7 @@
 
 #include "ThermoPhase.h"
 #include "PureFluidPhase.h"
-#include <stdio.h>
+#include <cstdio>
 #include "mix_defs.h"
 
 using namespace std;