From 103fd934a1be745b3484066a070cf601ce76ef61 Mon Sep 17 00:00:00 2001
From: Harry Moffat <hkmoffa@sandia.gov>
Date: Mon, 22 Dec 2008 22:40:44 +0000
Subject: [PATCH] Changed the algorithm for calculation of reference state
 properties. It will now calculate the reference state at an elevated pressure
 for elevated temperature for the real water equation of state. Core dumps
 were happening around the critical point without this change.

---
 Cantera/src/thermo/PDSS_HKFT.cpp  | 10 +++++-----
 Cantera/src/thermo/PDSS_Water.cpp | 12 ++++++++++++
 Cantera/src/thermo/PDSS_Water.h   | 11 +++++++++++
 3 files changed, 28 insertions(+), 5 deletions(-)

diff --git a/Cantera/src/thermo/PDSS_HKFT.cpp b/Cantera/src/thermo/PDSS_HKFT.cpp
index 52135b9fd..28ecd9399 100644
--- a/Cantera/src/thermo/PDSS_HKFT.cpp
+++ b/Cantera/src/thermo/PDSS_HKFT.cpp
@@ -392,7 +392,7 @@ namespace Cantera {
   doublereal 
   PDSS_HKFT::gibbs_RT_ref() const {
     doublereal m_psave = m_pres;
-    m_pres = OneAtm;
+    m_pres = m_waterSS->pref_safe(m_temp);
     doublereal ee = gibbs_RT();
     m_pres = m_psave;
     return ee;
@@ -401,7 +401,7 @@ namespace Cantera {
   doublereal 
   PDSS_HKFT::enthalpy_RT_ref() const {
     doublereal m_psave = m_pres;
-    m_pres = OneAtm;
+    m_pres = m_waterSS->pref_safe(m_temp);
     doublereal hh = enthalpy_RT();
     m_pres = m_psave;
     return hh;
@@ -410,7 +410,7 @@ namespace Cantera {
   doublereal 
   PDSS_HKFT::entropy_R_ref() const {
     doublereal m_psave = m_pres;
-    m_pres = OneAtm;
+    m_pres = m_waterSS->pref_safe(m_temp);
     doublereal ee = entropy_R();
     m_pres = m_psave;
     return ee;
@@ -419,7 +419,7 @@ namespace Cantera {
   doublereal 
   PDSS_HKFT::cp_R_ref() const {
     doublereal m_psave = m_pres;
-    m_pres = OneAtm;
+    m_pres = m_waterSS->pref_safe(m_temp);
     doublereal ee = cp_R();
     m_pres = m_psave;
     return ee;
@@ -428,7 +428,7 @@ namespace Cantera {
   doublereal
   PDSS_HKFT::molarVolume_ref() const {
     doublereal m_psave = m_pres;
-    m_pres = OneAtm;
+    m_pres = m_waterSS->pref_safe(m_temp);
     doublereal ee = molarVolume();
     m_pres = m_psave;
     return ee;
diff --git a/Cantera/src/thermo/PDSS_Water.cpp b/Cantera/src/thermo/PDSS_Water.cpp
index 17015c511..e3903cb98 100644
--- a/Cantera/src/thermo/PDSS_Water.cpp
+++ b/Cantera/src/thermo/PDSS_Water.cpp
@@ -473,6 +473,18 @@ namespace Cantera {
     m_sub->setState_TR(m_temp, m_dens);
   }
 
+  doublereal PDSS_Water::pref_safe(doublereal temp) const {
+    if (temp < m_sub->Tcrit()) {
+      doublereal pp = m_sub->psat_est(temp);
+      if (pp > OneAtm) {
+	return pp;
+      }
+    } else  {
+      return m_sub->Pcrit();
+    }
+    return OneAtm;
+  }
+
   // saturation pressure
   doublereal PDSS_Water::satPressure(doublereal t){
     doublereal pp = m_sub->psat(t, WATER_LIQUID);
diff --git a/Cantera/src/thermo/PDSS_Water.h b/Cantera/src/thermo/PDSS_Water.h
index b6ac06c0e..058c5e5f5 100644
--- a/Cantera/src/thermo/PDSS_Water.h
+++ b/Cantera/src/thermo/PDSS_Water.h
@@ -234,6 +234,17 @@ namespace Cantera {
      * @{
      */
 
+    //! Returns a reference pressure value that can be safely calculated by the
+    //! underlying real equation of state for water
+    /*!
+     *  Note, this function is needed because trying to calculate a one atm
+     *  value around the critical point will cause a crash
+     *
+     * @param temp  Temperature (Kelvin)
+     */
+    doublereal PDSS_Water::pref_safe(doublereal temp) const;
+
+
     //! Return the molar gibbs free energy divided by RT at reference pressure
     /*!
      * Returns the species reference state gibbs free energy divided by RT at the