Pushed the parameters for the cropping algorithm to the XML input file.

Cantera/src/base/ct_defs.h

@@ -9,7 +9,7 @@
  * The module physConstants is defined here.
  */
 
-/* $Author$
+/*
  * $Revision$
  * $Date$
  */

Cantera/src/base/ctml.cpp

@@ -621,7 +621,9 @@ namespace ctml {
    *                 and "" , for no conversion. The default value is ""
    *                 which implies that no conversion is allowed.
    */
-  doublereal getFloat(const Cantera::XML_Node& parent, std::string name, std::string type) {
+  doublereal getFloat(const Cantera::XML_Node& parent,
+		      const std::string &name,
+		      const std::string type) {
     if (!parent.hasChild(name)) 
       throw CanteraError("getFloat (called from XML Node \"" +
 			 parent.name() + "\"): ",
@@ -681,6 +683,53 @@ namespace ctml {
   }
 
   
+  bool getOptionalFloat(const Cantera::XML_Node& parent,
+			const std::string &name,
+			doublereal &fltRtn,
+			const std::string type) {
+    if (parent.hasChild(name)) {
+      fltRtn= getFloat(parent, name, type);
+      return true;
+    }
+    return false;
+  }
+
+  //  Get an optional floating-point value from a child element. 
+  /* 
+   *  Returns a doublereal value for the child named 'name' of element 'parent'. If
+   *  'type' is supplied and matches a known unit type, unit
+   *  conversion to SI will be done if the child element has an attribute
+   *  'units'.
+   *
+   * 
+   *
+   *  Example:  
+   *
+   * Code snipet:
+   *       @verbatim
+   const XML_Node &State_XMLNode;
+   doublereal pres = OneAtm;
+   bool exists = getOptionalFloat(State_XMLNode, "pressure", pres, "toSI");
+   @endverbatim
+   *
+   *  reads the corresponding XML file:
+   *  @verbatim
+   <state>
+     <pressure units="Pa"> 101325.0 </pressure>
+   <\state>
+   @endverbatim
+   *
+   *   @param parent reference to the XML_Node object of the parent XML element
+   *   @param name   Name of the XML child element
+   *   @param fltRtn Float Return. It will be overridden if the XML 
+   *                 element exists.
+   *   @param type   String type. Currently known types are "toSI"
+   *                 and "actEnergy",
+   *                 and "" , for no conversion. The default value is "",
+   *                 which implies that no conversion is allowed.
+   *
+   * @return returns true if the child element named "name" exists
+   */
   doublereal getFloatDefaultUnits(const Cantera::XML_Node& parent, std::string name,
 				  std::string defaultUnits, std::string type) {
 

Cantera/src/base/ctml.h

@@ -20,7 +20,8 @@
 #include "Array.h"
 
 //! The ctml namespace adds functionality to the XML object, by providing
-//! standard functions that read, write, and interpret XML files and object trees.
+//! standard functions that read, write, and interpret XML files and
+//! object trees.
 namespace ctml {
 
   //! const Specifying the CTML version number
@@ -542,7 +543,7 @@ namespace ctml {
    const XML_Node &State_XMLNode;
    doublereal pres = OneAtm;
    if (state_XMLNode.hasChild("pressure")) {
-   pres = getFloat(State_XMLNode, "pressure", "toSI");
+     pres = getFloat(State_XMLNode, "pressure", "toSI");
    }
    @endverbatim
    *
@@ -559,8 +560,48 @@ namespace ctml {
    *                 and "" , for no conversion. The default value is "",
    *                 which implies that no conversion is allowed.
    */
-  doublereal getFloat(const Cantera::XML_Node& parent, std::string name,
+  doublereal getFloat(const Cantera::XML_Node& parent, const std::string &name,
-		      std::string type="");
+		      const std::string type=""); 
+
+  //!  Get an optional floating-point value from a child element. 
+  /*! 
+   *  Returns a doublereal value for the child named 'name' of element 'parent'. If
+   *  'type' is supplied and matches a known unit type, unit
+   *  conversion to SI will be done if the child element has an attribute
+   *  'units'.
+   *
+   * 
+   *
+   *  Example:  
+   *
+   * Code snipet:
+   *       @verbatim
+   const XML_Node &State_XMLNode;
+   doublereal pres = OneAtm;
+   bool exists = getOptionalFloat(State_XMLNode, "pressure", pres, "toSI");
+   @endverbatim
+   *
+   *  reads the corresponding XML file:
+   *  @verbatim
+   <state>
+     <pressure units="Pa"> 101325.0 </pressure>
+   <\state>
+   @endverbatim
+   *
+   *   @param parent reference to the XML_Node object of the parent XML element
+   *   @param name   Name of the XML child element
+   *   @param fltRtn Float Return. It will be overridden if the XML 
+   *                 element exists.
+   *   @param type   String type. Currently known types are "toSI"
+   *                 and "actEnergy",
+   *                 and "" , for no conversion. The default value is "",
+   *                 which implies that no conversion is allowed.
+   *
+   * @return returns true if the child element named "name" exists
+   */
+  bool getOptionalFloat(const Cantera::XML_Node& parent, const std::string &name,
+			doublereal &fltRtn, const std::string type="");
+
 
   //!  Get a vector of floating-point values from a child element. 
   /*! 

Cantera/src/thermo/HMWSoln.cpp

@@ -29,6 +29,7 @@
 #include "ThermoFactory.h"
 #include "WaterProps.h"
 #include "PDSS_Water.h"
+
 #include <cmath>
 #include <cstdlib>
 
@@ -75,7 +76,7 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
-    CROP_ln_gamma_o_min(-10.0),
+    CROP_ln_gamma_o_min(-6.0),
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
@@ -132,7 +133,7 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
-    CROP_ln_gamma_o_min(-10.0),
+    CROP_ln_gamma_o_min(-6.0),
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
@@ -183,7 +184,7 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
-    CROP_ln_gamma_o_min(-10.0),
+    CROP_ln_gamma_o_min(-6.0),
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
@@ -240,7 +241,7 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
-    CROP_ln_gamma_o_min(-10.0),
+    CROP_ln_gamma_o_min(-6.0),
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
@@ -400,6 +401,7 @@ namespace Cantera {
       CROP_ln_gamma_o_max   = b.CROP_ln_gamma_o_max;
       CROP_ln_gamma_k_min   = b.CROP_ln_gamma_k_min;
       CROP_ln_gamma_k_max   = b.CROP_ln_gamma_k_max;
+      CROP_speciesCropped_  = b.CROP_speciesCropped_;
       m_CounterIJ           = b.m_CounterIJ;
       m_molalitiesCropped   = b.m_molalitiesCropped;
       m_molalitiesAreCropped= b.m_molalitiesAreCropped;
@@ -473,7 +475,7 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
-    CROP_ln_gamma_o_min(-10.0),
+    CROP_ln_gamma_o_min(-6.0),
     CROP_ln_gamma_o_max(3.0),
     CROP_ln_gamma_k_min(-5.0),
     CROP_ln_gamma_k_max(15.0),
@@ -1761,6 +1763,7 @@ namespace Cantera {
     m_gamma_tmp.resize(leng, 0.0);
 
     IMS_lnActCoeffMolal_.resize(m_kk, 0.0);
+    CROP_speciesCropped_.resize(m_kk, 0);
 
     counterIJ_setup();
   }
@@ -1823,24 +1826,33 @@ namespace Cantera {
     double lnxs = log(xx);
 
     for (int k = 1; k < m_kk; k++) {
+      CROP_speciesCropped_[k] = 0;
       m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
-      if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max)) {
+      if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max- 2.5 *lnxs)) {
-	m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max;
+	CROP_speciesCropped_[k] = 2;
+	m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max - 2.5 * lnxs;
       }
       if (m_lnActCoeffMolal_Unscaled[k] < (CROP_ln_gamma_k_min - 2.5 *lnxs)) {
 	// -1.0 and -1.5 caused multiple solutions
+	CROP_speciesCropped_[k] = 2;
 	m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_min - 2.5 * lnxs;
       }
     }
-
+    CROP_speciesCropped_[0] = 0;
     m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
     if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min) {
+      CROP_speciesCropped_[0] = 2;
       m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min;
     }
     if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max) {
+      CROP_speciesCropped_[0] = 2;
 	// -0.5 caused multiple solutions
       m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max;
     }
+    if (m_lnActCoeffMolal_Unscaled[0] > CROP_ln_gamma_o_max - 0.5 * lnxs) {
+      CROP_speciesCropped_[0] = 2;
+      m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_max - 0.5 * lnxs;
+    }
 
     /*
      * Now do the pH Scaling
@@ -3398,17 +3410,17 @@ namespace Cantera {
      */
     s_updatePitzer_dlnMolalityActCoeff_dT();
 
-    double xmolSolvent = moleFraction(m_indexSolvent);
+    //double xmolSolvent = moleFraction(m_indexSolvent);
-    double xx = MAX(m_xmolSolventMIN, xmolSolvent);
+    //double xx = MAX(m_xmolSolventMIN, xmolSolvent);
-    double lnxs = log(xx);
+    //    double lnxs = log(xx);
 
     for (int k = 1; k < m_kk; k++) {
-      if (m_lnActCoeffMolal_Unscaled[k] >= (CROP_ln_gamma_k_max + lnxs)) {
+      if (CROP_speciesCropped_[k] == 2) {
 	m_dlnActCoeffMolaldT_Unscaled[k] = 0.0;
       }
     }
 
-    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+    if (CROP_speciesCropped_[0]) {
       m_dlnActCoeffMolaldT_Unscaled[0] = 0.0;
     }
 
@@ -4265,18 +4277,19 @@ namespace Cantera {
      */
     s_updatePitzer_d2lnMolalityActCoeff_dT2();
 
-    double xmolSolvent = moleFraction(m_indexSolvent);
+    //double xmolSolvent = moleFraction(m_indexSolvent);
-    double xx = MAX(m_xmolSolventMIN, xmolSolvent);
+    //double xx = MAX(m_xmolSolventMIN, xmolSolvent);
-    double lnxs = log(xx);
+    //double lnxs = log(xx);
+
     for (int k = 1; k < m_kk; k++) {
-      if (m_lnActCoeffMolal_Unscaled[k] >= (CROP_ln_gamma_k_max + lnxs)) {
+      if (CROP_speciesCropped_[k] == 2) {
 	m_d2lnActCoeffMolaldT2_Unscaled[k] = 0.0;
       }
     }
 
-    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+    if (CROP_speciesCropped_[0]) {
       m_d2lnActCoeffMolaldT2_Unscaled[0] = 0.0;
-    }
+    }    
 
     /*
      * Scale the 2nd derivatives 
@@ -5153,16 +5166,14 @@ namespace Cantera {
   
     s_updatePitzer_dlnMolalityActCoeff_dP();
 
-    double xmolSolvent = moleFraction(m_indexSolvent);
+
-    double xx = MAX(m_xmolSolventMIN, xmolSolvent);
-    double lnxs = log(xx);
     for (int k = 1; k < m_kk; k++) {
-      if (m_lnActCoeffMolal_Unscaled[k] >= (CROP_ln_gamma_k_max + lnxs)) {
+      if (CROP_speciesCropped_[k] == 2) {
 	m_dlnActCoeffMolaldP_Unscaled[k] = 0.0;
       }
     }
 
-    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+    if (CROP_speciesCropped_[0]) {
       m_dlnActCoeffMolaldP_Unscaled[0] = 0.0;
     }
 

Cantera/src/thermo/HMWSoln.h

@@ -3153,7 +3153,8 @@ namespace Cantera {
      */
     doublereal IMS_slopegCut_;
 
-    //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
+    //! Parameter in the polyExp cutoff treatment having to do with 
+    //! rate of exp decay
     doublereal IMS_dgCut_;
 
     //! Parameter in the polyExp cutoff treatment having to do with rate of exp decay
@@ -3199,6 +3200,18 @@ namespace Cantera {
     doublereal CROP_ln_gamma_k_min;
     doublereal CROP_ln_gamma_k_max;
 
+    //! This is a boolean-type vector indicating whether
+    //! a species's activity coefficient is in the cropped regime
+    /*!
+     *
+     *  0 = Not in cropped regime
+     *  1 = In a transition regime where it is altered but there
+     *      still may be a temperature or pressure dependence
+     *  2 = In a cropped regime where there is no temperature 
+     *      or pressure dependence
+     */
+    mutable std::vector<int> CROP_speciesCropped_;
+
 
     //! Local error routine
     /*!
@@ -3209,7 +3222,8 @@ namespace Cantera {
     //!  Initialize all of the species - dependent lengths in the object
     void initLengths();
 
-    //! Apply the current phScale to a set of activity Coefficients or activities
+    //! Apply the current phScale to a set of activity Coefficients or 
+    //! activities
     /*!
      *  See the Eq3/6 Manual for a thorough discussion.
      *
@@ -3446,7 +3460,14 @@ namespace Cantera {
      *                 containing the
      *                 neutral - cation - anion interaction
      */
-    void readXMLZetaCation(XML_Node &BinSalt);
+    void readXMLZetaCation(const XML_Node &BinSalt);
+
+    //! Process an XML node called "croppingCoefficients"
+    //! for the cropping coefficients values
+    /*!
+     * @param acNode Activity Coefficient XML Node
+     */
+    void readXMLCroppingCoefficients(const XML_Node &acNode);
 
     //! Precalculate the IMS Cutoff parameters for typeCutoff = 2
     void  calcIMSCutoffParams_();

Cantera/src/thermo/HMWSoln_input.cpp

@@ -902,7 +902,7 @@ namespace Cantera {
    * This node contains all of the parameters necessary to describe
    * the ternary interactions between a neutral, a cation and an anion
    */
-  void HMWSoln::readXMLZetaCation(XML_Node &BinSalt) {
+  void HMWSoln::readXMLZetaCation(const XML_Node &BinSalt) {
     string xname = BinSalt.name();
     if (xname != "zetaCation") {
       throw CanteraError("HMWSoln::readXMLZetaCation",
@@ -1002,7 +1002,35 @@ namespace Cantera {
     }
   }
    
+  // Process an XML node called "croppingCoefficients"
+  // for the cropping coefficients values
+  /*
+   * @param acNode Activity Coefficient XML Node
+   */
+  void HMWSoln::readXMLCroppingCoefficients(const XML_Node &acNode) {
 
+    if (acNode.hasChild("croppingCoefficients")) {
+      XML_Node &cropNode = acNode.child("croppingCoefficients");
+      if (cropNode.hasChild("ln_gamma_k_min")) {
+	XML_Node &gkminNode = cropNode.child("ln_gamma_k_min");
+	getOptionalFloat(gkminNode, "pureSolventValue", CROP_ln_gamma_k_min);
+      }
+      if (cropNode.hasChild("ln_gamma_k_max")) {
+	XML_Node &gkmaxNode = cropNode.child("ln_gamma_k_max");
+	getOptionalFloat(gkmaxNode, "pureSolventValue", CROP_ln_gamma_k_max);
+      }
+
+      if (cropNode.hasChild("ln_gamma_o_min")) {
+	XML_Node &gominNode = cropNode.child("ln_gamma_o_min");
+	getOptionalFloat(gominNode, "pureSolventValue", CROP_ln_gamma_o_min);
+      }
+
+      if (cropNode.hasChild("ln_gamma_o_max")) {
+	XML_Node &gomaxNode = cropNode.child("ln_gamma_o_max");
+	getOptionalFloat(gomaxNode, "pureSolventValue", CROP_ln_gamma_o_max);
+      } 
+    }
+  }
 
   /*
    *  Initialization routine for a HMWSoln phase.
@@ -1482,10 +1510,11 @@ namespace Cantera {
 	}
 	if (jmap > -1) {
 	  const XML_Node& sp = *xspecies[jmap];
-	  if (sp.hasChild("stoichIsMods")) {
+	  getOptionalFloat(sp, "stoichIsMods",  m_speciesCharge_Stoich[k]);
-	    double val = getFloat(sp, "stoichIsMods");
+	  // if (sp.hasChild("stoichIsMods")) {
-	    m_speciesCharge_Stoich[k] = val;
+	  // double val = getFloat(sp, "stoichIsMods");
-	  }
+	  //m_speciesCharge_Stoich[k] = val;
+	  //}
 	}
       }
    
@@ -1513,7 +1542,8 @@ namespace Cantera {
 	  }
 	}
       }
-
+      
+    
       /*
        * Loop through the children getting multiple instances of
        * parameters
@@ -1547,6 +1577,8 @@ namespace Cantera {
 	}
       }
 
+      // Go look up the optional Cropping parameters
+      readXMLCroppingCoefficients(acNode);
 
     }
 

Cantera/src/thermo/SurfPhase.cpp

@@ -454,8 +454,8 @@ namespace Cantera {
 
   void SurfPhase::setStateFromXML(const XML_Node& state) {
 
-    if (state.hasChild("temperature")) {
+    double t;
-      double t = getFloat(state, "temperature", "temperature");
+    if (getOptionalFloat(state, "temperature", t, "temperature")) {
       setTemperature(t);
     }