From 6942cb3260685aa8d2b6b25452d44b02cdc22706 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Sat, 17 Jan 2009 05:42:50 +0000
Subject: [PATCH] Added pure cropping on the molar gamma_o and gamma_k. this
 seems to be the most robust and to-the-point.

---
 Cantera/src/thermo/HMWSoln.cpp       | 66 +++++++++++++++++++++++++++-
 Cantera/src/thermo/HMWSoln.h         |  3 ++
 Cantera/src/thermo/HMWSoln_input.cpp |  4 +-
 3 files changed, 69 insertions(+), 4 deletions(-)

diff --git a/Cantera/src/thermo/HMWSoln.cpp b/Cantera/src/thermo/HMWSoln.cpp
index 5e7c9c72b..d0b8cc08d 100644
--- a/Cantera/src/thermo/HMWSoln.cpp
+++ b/Cantera/src/thermo/HMWSoln.cpp
@@ -75,6 +75,8 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
+    CROP_ln_gamma_o_min(-25.0),
+    CROP_ln_gamma_k_max(23.0),
     m_debugCalc(0)
   {
     for (int i = 0; i < 17; i++) {
@@ -128,6 +130,8 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
+    CROP_ln_gamma_o_min(-25.0),
+    CROP_ln_gamma_k_max(23.0),
     m_debugCalc(0)
   {
     for (int i = 0; i < 17; i++) {
@@ -175,6 +179,8 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
+    CROP_ln_gamma_o_min(-25.0),
+   CROP_ln_gamma_k_max(23.0),
     m_debugCalc(0)
   {
     for (int i = 0; i < 17; i++) {
@@ -228,6 +234,8 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
+    CROP_ln_gamma_o_min(-25.0),
+   CROP_ln_gamma_k_max(23.0),
     m_debugCalc(0)
   {
     /*
@@ -380,7 +388,8 @@ namespace Cantera {
       MC_apCut_             = b.MC_apCut_;
       MC_bpCut_             = b.MC_bpCut_;
       MC_cpCut_             = b.MC_cpCut_;
-
+      CROP_ln_gamma_o_min   = b.CROP_ln_gamma_o_min;
+      CROP_ln_gamma_k_max   = b.CROP_ln_gamma_k_max;
       m_CounterIJ           = b.m_CounterIJ;
       m_molalitiesCropped   = b.m_molalitiesCropped;
       m_molalitiesAreCropped= b.m_molalitiesAreCropped;
@@ -454,6 +463,8 @@ namespace Cantera {
     MC_apCut_(0.0),
     MC_bpCut_(0.0),
     MC_cpCut_(0.0),
+    CROP_ln_gamma_o_min(-25.0),
+    CROP_ln_gamma_k_max(23.0),
     m_debugCalc(0)
   {
     if (testProb != 1) {
@@ -1800,13 +1811,19 @@ namespace Cantera {
     double xmolSolvent = moleFraction(m_indexSolvent);
     double xx = MAX(m_xmolSolventMIN, xmolSolvent);
     double lnActCoeffMolal0 = - log(xx) + (xx - 1.0)/xx;
+    double lnxs = log(xx);
 
     for (int k = 1; k < m_kk; k++) {
       m_lnActCoeffMolal_Unscaled[k] += IMS_lnActCoeffMolal_[k];
+      if (m_lnActCoeffMolal_Unscaled[k] > (CROP_ln_gamma_k_max + lnxs)) {
+	m_lnActCoeffMolal_Unscaled[k] = CROP_ln_gamma_k_max + lnxs;
+      }
     }
 
     m_lnActCoeffMolal_Unscaled[0] += (IMS_lnActCoeffMolal_[0] - lnActCoeffMolal0);
-
+    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+      m_lnActCoeffMolal_Unscaled[0] = CROP_ln_gamma_o_min - lnxs;
+    }
 
     /*
      * Now do the pH Scaling
@@ -3345,10 +3362,27 @@ namespace Cantera {
      *  Do the actual calculation of the unscaled temperature derivatives
      */
     s_updatePitzer_dlnMolalityActCoeff_dT();
+
+    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)) {
+	m_dlnActCoeffMolaldT_Unscaled[k] = 0.0;
+      }
+    }
+
+    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+      m_dlnActCoeffMolaldT_Unscaled[0] = 0.0;
+    }
+
     /*
      *  Do the pH scaling to the derivatives 
      */
     s_updateScaling_pHScaling_dT();
+
+
   }
 
   /*************************************************************************************/
@@ -4197,6 +4231,20 @@ namespace Cantera {
      * Calculate the unscaled 2nd derivatives
      */
     s_updatePitzer_d2lnMolalityActCoeff_dT2();
+
+    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)) {
+	m_d2lnActCoeffMolaldT2_Unscaled[k] = 0.0;
+      }
+    }
+
+    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+      m_d2lnActCoeffMolaldT2_Unscaled[0] = 0.0;
+    }
+
     /*
      * Scale the 2nd derivatives 
      */
@@ -5074,6 +5122,20 @@ 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)) {
+	m_dlnActCoeffMolaldP_Unscaled[k] = 0.0;
+      }
+    }
+
+    if (m_lnActCoeffMolal_Unscaled[0] < CROP_ln_gamma_o_min - lnxs) {
+      m_dlnActCoeffMolaldP_Unscaled[0] = 0.0;
+    }
+
+
     s_updateScaling_pHScaling_dP();
   }
 
diff --git a/Cantera/src/thermo/HMWSoln.h b/Cantera/src/thermo/HMWSoln.h
index 2ed9dc5ab..97d2015ef 100644
--- a/Cantera/src/thermo/HMWSoln.h
+++ b/Cantera/src/thermo/HMWSoln.h
@@ -3193,6 +3193,9 @@ namespace Cantera {
     //! Parameter in the Molality Exp cutoff treatment
     doublereal MC_cpCut_;
 
+    doublereal CROP_ln_gamma_o_min;
+
+    doublereal CROP_ln_gamma_k_max;
 
     //! Local error routine
     /*!
diff --git a/Cantera/src/thermo/HMWSoln_input.cpp b/Cantera/src/thermo/HMWSoln_input.cpp
index f3e6235a0..599f3c919 100644
--- a/Cantera/src/thermo/HMWSoln_input.cpp
+++ b/Cantera/src/thermo/HMWSoln_input.cpp
@@ -1685,9 +1685,9 @@ namespace Cantera {
 
   // Precalculate the MC Cutoff parameters
   void  HMWSoln::calcMCCutoffParams_() {
-    MC_X_o_min_ = 0.3;
+    MC_X_o_min_ = 0.35;
     MC_X_o_cutoff_ = 0.6;
-    MC_slopepCut_ = 0.05;
+    MC_slopepCut_ = 0.02;
     MC_cpCut_ = 0.25;
 
     // Initial starting values