From 6a7e18d918e00def3b305769d831cfc6b9005312 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Thu, 30 Jul 2020 17:29:39 -0500
Subject: [PATCH 01/17] [purefluid] Fix/clarify error message for critical
 point

- Logic is updated for PQ setter in order to avoid 'no convergence' at
critical point and clarify error for supercritical conditions
- Error message is updated to properly describe failure for an attempt to
assign the critical state via TQ setters.
---
 src/tpx/Sub.cpp | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 718e02b2511..ecda63c7200 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -370,14 +370,15 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
                x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
         break;
     case PropertyPair::PX:
-        temp = Tsat(x0);
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (temp >= Tcrit()) {
+        } else if (x0 >= Pcrit()) {
             throw CanteraError("Substance::Set",
-                               "Can't set vapor fraction above the critical point");
+                               "Vapor fraction can only be set below the "
+                               "critical point");
         } else {
+            temp = Tsat(x0);
             set_T(temp);
             update_sat();
             Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
@@ -389,7 +390,8 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
                                "Invalid vapor fraction, {}", y0);
         } else if (x0 >= Tcrit()) {
             throw CanteraError("Substance::Set",
-                               "Can't set vapor fraction above the critical point");
+                               "Vapor fraction can only be set below the "
+                               "critical point");
         } else {
             set_T(x0);
             update_sat();

From d5f25563fc0e7a101a98eee0d9df8d8995d922f3 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Aug 2020 16:41:06 -0500
Subject: [PATCH 02/17] [purefluid] Fix saturation pressure calculation

- Avoid illegal temperatures in saturation pressure gradient
- Fix logic at critical point
---
 src/tpx/Sub.cpp   | 14 +++++++++++---
 src/tpx/Water.cpp |  2 +-
 2 files changed, 12 insertions(+), 4 deletions(-)

diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index ecda63c7200..6732d7a560d 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -202,8 +202,16 @@ double Substance::dPsdT()
 {
     double tsave = T;
     double ps1 = Ps();
-    T += DeltaT;
-    double dpdt = (Ps() - ps1)/DeltaT;
+    double dpdt;
+    if (T + DeltaT < Tcrit()) {
+        T += DeltaT;
+        dpdt = (Ps() - ps1)/DeltaT;
+    } else if (T - DeltaT > Tmin()) {
+        T -= DeltaT;
+        dpdt = (ps1 - Ps())/DeltaT;
+    } else {
+        throw CanteraError("Substance::dPsdT", "Illegal temperature value");
+    }
     T = tsave;
     return dpdt;
 }
@@ -445,7 +453,7 @@ double Substance::Ps()
 
 void Substance::update_sat()
 {
-    if ((T != Tslast) && (T < Tcrit())) {
+    if ((T != Tslast) && (T <= Tcrit())) {
         double Rho_save = Rho;
         double pp = Psat();
         double lps = log(pp);
diff --git a/src/tpx/Water.cpp b/src/tpx/Water.cpp
index 9079ce5f46d..8b2abf569e2 100644
--- a/src/tpx/Water.cpp
+++ b/src/tpx/Water.cpp
@@ -176,7 +176,7 @@ double water::ldens()
 {
     double sum=0;
     int i;
-    if ((T < Tmn) || (T >= Tc)) {
+    if ((T < Tmn) || (T > Tc)) {
         throw CanteraError("water::ldens",
                            "Temperature out of range. T = {}", T);
     }

From dbfa823de7c14b114b09573184cd9c7d36b0bd12 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Aug 2020 17:43:19 -0500
Subject: [PATCH 03/17] [purefluid] Prevent silent failures of substance
 updates

- do not allow for for silent failure of Substance::update_sat for
  illegal temperatures
- small formatting changes of error messages
---
 src/tpx/Sub.cpp | 31 ++++++++++++++++++-------------
 1 file changed, 18 insertions(+), 13 deletions(-)

diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 6732d7a560d..b9b8d32e803 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -210,7 +210,7 @@ double Substance::dPsdT()
         T -= DeltaT;
         dpdt = (ps1 - Ps())/DeltaT;
     } else {
-        throw CanteraError("Substance::dPsdT", "Illegal temperature value");
+        throw CanteraError("Substance::dPsdT", "Illegal temperature value: {}", T);
     }
     T = tsave;
     return dpdt;
@@ -246,7 +246,7 @@ double Substance::x()
 double Substance::Tsat(double p)
 {
     if (p <= 0.0 || p > Pcrit()) {
-        throw CanteraError("Substance::Tsat", "illegal pressure value");
+        throw CanteraError("Substance::Tsat", "Illegal pressure value: {}", p);
     }
     int LoopCount = 0;
     double tol = 1.e-6*p;
@@ -427,7 +427,7 @@ void Substance::set_T(double t0)
     if ((t0 >= Tmin()) && (t0 <= Tmax())) {
         T = t0;
     } else {
-        throw CanteraError("Substance::set_T", "illegal temperature: {}", t0);
+        throw CanteraError("Substance::set_T", "Illegal temperature: {}", t0);
     }
 }
 
@@ -437,7 +437,7 @@ void Substance::set_v(double v0)
         Rho = 1.0/v0;
     } else {
         throw CanteraError("Substance::set_v",
-                           "negative specific volume: {}", v0);
+                           "Negative specific volume: {}", v0);
     }
 }
 
@@ -445,7 +445,7 @@ double Substance::Ps()
 {
     if (T < Tmin() || T > Tcrit()) {
         throw CanteraError("Substance::Ps",
-                           "illegal temperature value {}", T);
+                           "Illegal temperature value: {}", T);
     }
     update_sat();
     return Pst;
@@ -453,7 +453,12 @@ double Substance::Ps()
 
 void Substance::update_sat()
 {
-    if ((T != Tslast) && (T <= Tcrit())) {
+    if (T < Tmin() || T > Tcrit()) {
+        // prevent silent failures
+        throw CanteraError("Substance::update_sat",
+                           "Temperature value '{}' is outside of valid range: "
+                           "{} < T < {}", T, Tmin(), Tcrit());
+    } else if (T != Tslast) {
         double Rho_save = Rho;
         double pp = Psat();
         double lps = log(pp);
@@ -506,11 +511,11 @@ void Substance::update_sat()
         }
         if (Rhf <= Rhv) {
             throw CanteraError("Substance::update_sat",
-                "wrong root found for sat. liquid or vapor at P = {}", pp);
+                "Wrong root found for sat. liquid or vapor at P = {}", pp);
         }
 
         if (i >= 20) {
-            throw CanteraError("Substance::update_sat", "no convergence");
+            throw CanteraError("Substance::update_sat", "No convergence");
         } else {
             Pst = pp;
             Tslast = T;
@@ -533,7 +538,7 @@ double Substance::vprop(propertyFlag::type ijob)
     case propertyFlag::P:
         return Pp();
     default:
-        throw CanteraError("Substance::vprop", "invalid job index");
+        throw CanteraError("Substance::vprop", "Invalid job index");
     }
 }
 
@@ -562,7 +567,7 @@ int Substance::Lever(int itp, double sat, double val, propertyFlag::type ifunc)
             return 0;
         }
     } else {
-        throw CanteraError("Substance::Lever","general error");
+        throw CanteraError("Substance::Lever", "General error");
     }
     Set(PropertyPair::TX, T, 1.0);
     double Valg = vprop(ifunc);
@@ -771,7 +776,7 @@ void Substance::set_TPp(double Temp, double Pressure)
                     Pmax = P_here;
                 }
                 if (Vmin >= Vmax) {
-                    throw CanteraError("Substance::set_TPp","Vmin >= Vmax");
+                    throw CanteraError("Substance::set_TPp", "Vmin >= Vmax");
                 } else if ((Vmin > 0.0) && (Vmax < Big)) {
                     kbr = 1;
                 }
@@ -803,14 +808,14 @@ void Substance::set_TPp(double Temp, double Pressure)
         }
         v_here += dv;
         if (dv == 0.0) {
-            throw CanteraError("Substance::set_TPp","dv = 0 and no convergence");
+            throw CanteraError("Substance::set_TPp", "dv = 0 and no convergence");
         }
         Set(PropertyPair::TV, Temp, v_here);
         LoopCount++;
         if (LoopCount > 100) {
             Set(PropertyPair::TV, Temp, v_save);
             throw CanteraError("Substance::set_TPp",
-                               "no convergence for P* = {}, V* = {}",
+                               "No convergence for P* = {}, V* = {}",
                                Pressure/Pcrit(), v_save/Vcrit());
         }
     }

From 9ebee40977bf9f86410248e045119997454bef98 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Aug 2020 17:59:10 -0500
Subject: [PATCH 04/17] [purefluid] Simplify exception handling

---
 src/tpx/Sub.cpp | 31 +++++++++++++------------------
 1 file changed, 13 insertions(+), 18 deletions(-)

diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index b9b8d32e803..a3b871325e2 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -206,11 +206,10 @@ double Substance::dPsdT()
     if (T + DeltaT < Tcrit()) {
         T += DeltaT;
         dpdt = (Ps() - ps1)/DeltaT;
-    } else if (T - DeltaT > Tmin()) {
+    } else {
         T -= DeltaT;
+        // Ps() will fail for illegal temperature
         dpdt = (ps1 - Ps())/DeltaT;
-    } else {
-        throw CanteraError("Substance::dPsdT", "Illegal temperature value: {}", T);
     }
     T = tsave;
     return dpdt;
@@ -245,12 +244,17 @@ double Substance::x()
 
 double Substance::Tsat(double p)
 {
-    if (p <= 0.0 || p > Pcrit()) {
-        throw CanteraError("Substance::Tsat", "Illegal pressure value: {}", p);
+    double Tsave = T;
+    T = Tmin();
+    if (p < Ps() || p > Pcrit()) {
+        throw CanteraError("Substance::Tsat",
+                           "Illegal pressure value: {}", p);
     }
+    T = Tsave;
+    Tsave = T;
+
     int LoopCount = 0;
     double tol = 1.e-6*p;
-    double Tsave = T;
     if (T < Tmin()) {
         T = 0.5*(Tcrit() - Tmin());
     }
@@ -453,14 +457,9 @@ double Substance::Ps()
 
 void Substance::update_sat()
 {
-    if (T < Tmin() || T > Tcrit()) {
-        // prevent silent failures
-        throw CanteraError("Substance::update_sat",
-                           "Temperature value '{}' is outside of valid range: "
-                           "{} < T < {}", T, Tmin(), Tcrit());
-    } else if (T != Tslast) {
+    if (T != Tslast) {
         double Rho_save = Rho;
-        double pp = Psat();
+        double pp = Psat(); // illegal temperatures are caught by Psat()
         double lps = log(pp);
         // trial value = Psat from correlation
         int i;
@@ -489,7 +488,7 @@ void Substance::update_sat()
                 dg = - dg;
             }
 
-            if (fabs(dg) < 0.001 && Rhf > Rhv) {
+            if (fabs(dg) < 0.001) {
                 break;
             }
             double dp = dg/(1.0/Rhv - 1.0/Rhf);
@@ -509,10 +508,6 @@ void Substance::update_sat()
                 lps = log(pp);
             }
         }
-        if (Rhf <= Rhv) {
-            throw CanteraError("Substance::update_sat",
-                "Wrong root found for sat. liquid or vapor at P = {}", pp);
-        }
 
         if (i >= 20) {
             throw CanteraError("Substance::update_sat", "No convergence");

From e73bd022e52505ebad934753e9f2709044b10eb5 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 24 Aug 2020 21:24:56 -0500
Subject: [PATCH 05/17] [purefluid] Allow for TQ/PQ setters in critical point

---
 src/tpx/Sub.cpp | 23 ++++++++---------------
 1 file changed, 8 insertions(+), 15 deletions(-)

diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index a3b871325e2..9dd0475f8f0 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -255,20 +255,13 @@ double Substance::Tsat(double p)
 
     int LoopCount = 0;
     double tol = 1.e-6*p;
-    if (T < Tmin()) {
-        T = 0.5*(Tcrit() - Tmin());
-    }
-    if (T >= Tcrit()) {
+    if (T < Tmin() || T > Tcrit()) {
         T = 0.5*(Tcrit() - Tmin());
     }
     double dp = 10*tol;
     while (fabs(dp) > tol) {
-        if (T > Tcrit()) {
-            T = Tcrit() - 0.001;
-        }
-        if (T < Tmin()) {
-            T = Tmin() + 0.001;
-        }
+        T = std::min(T, Tcrit());
+        T = std::max(T, Tmin());
         dp = p - Ps();
         double dt = dp/dPsdT();
         double dta = fabs(dt);
@@ -280,7 +273,7 @@ double Substance::Tsat(double p)
         LoopCount++;
         if (LoopCount > 100) {
             T = Tsave;
-            throw CanteraError("Substance::Tsat", "No convergence");
+            throw CanteraError("Substance::Tsat", "No convergence: p = {}", p);
         }
     }
     double tsat = T;
@@ -385,9 +378,9 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (x0 >= Pcrit()) {
+        } else if (x0 > Pcrit()) {
             throw CanteraError("Substance::Set",
-                               "Vapor fraction can only be set below the "
+                               "Vapor fraction cannot be set above the "
                                "critical point");
         } else {
             temp = Tsat(x0);
@@ -400,9 +393,9 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (x0 >= Tcrit()) {
+        } else if (x0 > Tcrit()) {
             throw CanteraError("Substance::Set",
-                               "Vapor fraction can only be set below the "
+                               "Vapor fraction cannot be set above the "
                                "critical point");
         } else {
             set_T(x0);

From 476068a2285ec65899861a35ba8aea30d1742585 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Aug 2020 10:55:13 -0500
Subject: [PATCH 06/17] [unittests] Add tests for critical point and saturation
 pressure limits

---
 .../cython/cantera/test/test_purefluid.py     | 50 +++++++++++++++++++
 1 file changed, 50 insertions(+)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 73501d2b27b..5ae5cbe69dc 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -168,6 +168,56 @@ def test_fd_properties_twophase(self):
         self.assertEqual(self.water.isothermal_compressibility, np.inf)
         self.assertEqual(self.water.thermal_expansion_coeff, np.inf)
 
+    def test_quality_exceptions(self):
+        # Critical point
+        self.water.TP = 300, ct.one_atm
+        self.water.TQ = self.water.critical_temperature, .5
+        self.assertNear(self.water.P, self.water.critical_pressure)
+        self.water.TP = 300, ct.one_atm
+        self.water.PQ = self.water.critical_pressure, .5
+        self.assertNear(self.water.T, self.water.critical_temperature)
+
+        # Supercritical
+        with self.assertRaisesRegex(ct.CanteraError, 'above the critical point'):
+            self.water.TQ = 1.001 * self.water.critical_temperature, 0.
+        with self.assertRaisesRegex(ct.CanteraError, 'above the critical point'):
+            self.water.PQ = 1.001 * self.water.critical_pressure, 0.
+
+        # Q negative
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.TQ = 373.15, -.001
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.PQ = ct.one_atm, -.001
+
+        # Q larger than one
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.TQ = 373.15, 1.001
+        with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
+            self.water.PQ = ct.one_atm, 1.001
+
+    def test_saturation_near_limits(self):
+        # low temperature limit (triple point)
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = self.water.min_temp, ct.one_atm
+        psat = self.water.P_sat
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = 300, psat
+        self.assertNear(self.water.T_sat, self.water.min_temp)
+
+        # high temperature limit (critical point) - saturation temperature
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = self.water.critical_temperature, self.water.critical_pressure
+        self.assertNear(self.water.T_sat, self.water.critical_temperature)
+
+        # high temperature limit (critical point) - saturation pressure
+        self.water.TP = 300, ct.one_atm
+        self.water.P_sat # ensure that solver buffers sufficiently different values
+        self.water.TP = self.water.critical_temperature, self.water.critical_pressure
+        self.assertNear(self.water.P_sat, self.water.critical_pressure)
+
     def test_TPQ(self):
         self.water.TQ = 400, 0.8
         T, P, Q = self.water.TPQ

From 5fe0cf5d4d31fbe204599dceaf299a8fa5832175 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Thu, 30 Jul 2020 17:56:42 -0500
Subject: [PATCH 07/17] [example] create Python example illustrating vapor dome

---
 .../cantera/examples/thermo/vapordome.py      | 94 +++++++++++++++++++
 1 file changed, 94 insertions(+)
 create mode 100644 interfaces/cython/cantera/examples/thermo/vapordome.py

diff --git a/interfaces/cython/cantera/examples/thermo/vapordome.py b/interfaces/cython/cantera/examples/thermo/vapordome.py
new file mode 100644
index 00000000000..e759905b3c3
--- /dev/null
+++ b/interfaces/cython/cantera/examples/thermo/vapordome.py
@@ -0,0 +1,94 @@
+"""
+This example generates a saturated steam table and plots the vapor dome. The
+steam table corresponds to data typically found in thermodynamic text books
+and uses the same customary units.
+
+Requires: Cantera >= 2.5.0, matplotlib >= 2.0, pandas >= 1.1.0, numpy >= 1.12
+"""
+
+import cantera as ct
+import pandas as pd
+import numpy as np
+from matplotlib import pyplot as plt
+
+w = ct.Water()
+
+# create colums
+columns = ['T', 'P',
+           'vf', 'vfg', 'vg',
+           'uf', 'ufg', 'ug',
+           'hf', 'hfg', 'hg',
+           'sf', 'sfg', 'sg']
+
+# temperatures correspond to Engineering Thermodynamics, Moran et al. (9th ed),
+# Table A-2; additional data points are added close to the critical point;
+# w.min_temp is equal to the triple point temperature
+degc = np.hstack([np.array([w.min_temp - 273.15, 4, 5, 6, 8]),
+                  np.arange(10, 37), np.array([38]),
+                  np.arange(40, 100, 5), np.arange(100, 300, 10),
+                  np.arange(300, 380, 20), np.arange(370, 374),
+                  np.array([w.critical_temperature - 273.15])])
+
+df = pd.DataFrame(0, index=np.arange(len(degc)), columns=columns)
+df.T = degc
+
+arr = ct.SolutionArray(w, len(degc))
+
+# saturated vapor data
+arr.TQ = degc + 273.15, 1
+df.P = arr.P_sat / 1.e5
+df.vg = arr.v
+df.ug = arr.int_energy_mass / 1.e3
+df.hg = arr.enthalpy_mass / 1.e3
+df.sg = arr.entropy_mass / 1.e3
+
+# saturated liquid data
+arr.TQ = degc + 273.15, 0
+df.vf = arr.v
+df.uf = arr.int_energy_mass / 1.e3
+df.hf = arr.enthalpy_mass / 1.e3
+df.sf = arr.entropy_mass / 1.e3
+
+# delta values
+df.vfg = df.vg - df.vf
+df.ufg = df.ug - df.uf
+df.hfg = df.hg - df.hf
+df.sfg = df.sg - df.sf
+
+# reference state (triple point; liquid state)
+w.TQ = w.min_temp, 0
+uf0 = w.int_energy_mass / 1.e3
+hf0 = w.enthalpy_mass / 1.e3
+sf0 = w.entropy_mass / 1.e3
+pv0 = w.P * w.v / 1.e3
+
+# change reference state
+df.ug -= uf0
+df.uf -= uf0
+df.hg -= hf0 - pv0
+df.hf -= hf0 - pv0
+df.sg -= sf0
+df.sf -= sf0
+
+# print and write saturated steam table to csv file
+print(df)
+df.to_csv('saturated_steam_T.csv', index=False)
+
+# illustrate the vapor dome in a P-v diagram
+plt.semilogx(df.vf.values, df.P.values, label='Saturated liquid')
+plt.semilogx(df.vg.values, df.P.values, label='Saturated vapor')
+plt.semilogx(df.vg.values[-1], df.P.values[-1], 'o', label='Critical point')
+plt.xlabel(r'Specific volume - $v$ ($\mathrm{m^3/kg}$)')
+plt.ylabel(r'Presssure - $P$ (bar)')
+plt.legend()
+
+# illustrate the vapor dome in a T-s diagram
+plt.figure()
+plt.plot(df.sf.values, df['T'].values, label='Saturated liquid')
+plt.plot(df.sg.values, df['T'].values, label='Saturated vapor')
+plt.plot(df.sg.values[-1], df['T'].values[-1], 'o', label='Critical point')
+plt.xlabel(r'Specific entropy - $s$ ($\mathrm{kJ/kg-K}$)')
+plt.ylabel(r'Temperature - $T$ (${}^\circ C$)')
+plt.legend()
+
+plt.show()

From 9f50f2d2ce13212c00bb5882c5c8f1792321d141 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sat, 22 Aug 2020 20:57:20 -0500
Subject: [PATCH 08/17] [unittests] Add tests for tpx substance exceptions

---
 .../cython/cantera/test/test_purefluid.py     | 24 ++++++++++++++++---
 1 file changed, 21 insertions(+), 3 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 5ae5cbe69dc..77de0cdf38c 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -34,6 +34,8 @@ def test_substance_set(self):
         self.water.TV = 400, 1.45
         self.assertNear(self.water.T, 400)
         self.assertNear(self.water.v, 1.45)
+        with self.assertRaisesRegex(ct.CanteraError, 'Negative specific volume'):
+            self.water.TV = 300, -1.
 
         self.water.PV = 101325, 1.45
         self.assertNear(self.water.P, 101325)
@@ -196,7 +198,7 @@ def test_quality_exceptions(self):
             self.water.PQ = ct.one_atm, 1.001
 
     def test_saturation_near_limits(self):
-        # low temperature limit (triple point)
+        # Low temperature limit (triple point)
         self.water.TP = 300, ct.one_atm
         self.water.P_sat # ensure that solver buffers sufficiently different values
         self.water.TP = self.water.min_temp, ct.one_atm
@@ -206,18 +208,34 @@ def test_saturation_near_limits(self):
         self.water.TP = 300, psat
         self.assertNear(self.water.T_sat, self.water.min_temp)
 
-        # high temperature limit (critical point) - saturation temperature
+        # High temperature limit (critical point) - saturation temperature
         self.water.TP = 300, ct.one_atm
         self.water.P_sat # ensure that solver buffers sufficiently different values
         self.water.TP = self.water.critical_temperature, self.water.critical_pressure
         self.assertNear(self.water.T_sat, self.water.critical_temperature)
 
-        # high temperature limit (critical point) - saturation pressure
+        # High temperature limit (critical point) - saturation pressure
         self.water.TP = 300, ct.one_atm
         self.water.P_sat # ensure that solver buffers sufficiently different values
         self.water.TP = self.water.critical_temperature, self.water.critical_pressure
         self.assertNear(self.water.P_sat, self.water.critical_pressure)
 
+        # Supercricital
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal temperature value'):
+            self.water.TP = 1.001 * self.water.critical_temperature, self.water.critical_pressure
+            self.water.P_sat
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
+            self.water.TP = self.water.critical_temperature, 1.001 * self.water.critical_pressure
+            self.water.T_sat
+
+        # Below triple point
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal temperature'):
+            self.water.TP = .999 * self.water.min_temp, ct.one_atm
+            self.water.P_sat
+        with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
+            self.water.TP = 300, .999 * psat
+            self.water.T_sat
+
     def test_TPQ(self):
         self.water.TQ = 400, 0.8
         T, P, Q = self.water.TPQ

From db0fd17e72003f0ad07f39facea03b25e31548c0 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Tue, 25 Aug 2020 00:30:57 -0500
Subject: [PATCH 09/17] [purefluid] Disable TP setter for saturated mixture

---
 interfaces/cython/cantera/test/test_purefluid.py |  5 +++++
 src/tpx/Sub.cpp                                  | 10 ++++++----
 2 files changed, 11 insertions(+), 4 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 77de0cdf38c..ad516c5a317 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -197,6 +197,11 @@ def test_quality_exceptions(self):
         with self.assertRaisesRegex(ct.CanteraError, 'Invalid vapor fraction'):
             self.water.PQ = ct.one_atm, 1.001
 
+    def test_saturated_mixture(self):
+        self.water.TP = 300, ct.one_atm
+        with self.assertRaisesRegex(ct.CanteraError, 'Saturated mixture detected'):
+            self.water.TP = 300, self.water.P_sat
+
     def test_saturation_near_limits(self):
         # Low temperature limit (triple point)
         self.water.TP = 300, ct.one_atm
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 9dd0475f8f0..f324cd91838 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -327,15 +327,17 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
                x0, y0, TolAbsP, TolAbsV, TolRel, TolRel);
         break;
     case PropertyPair::TP:
+        set_T(x0);
         if (x0 < Tcrit()) {
-            set_T(x0);
-            if (y0 < Ps()) {
+            if (fabs(y0 - Ps()) / y0 < TolRel) {
+                throw CanteraError("Substance::Set",
+                                   "Saturated mixture detected: use vapor "
+                                   "fraction to specify state instead");
+            } else if (y0 < Ps()) {
                 Set(PropertyPair::TX, x0, 1.0);
             } else {
                 Set(PropertyPair::TX, x0, 0.0);
             }
-        } else {
-            set_T(x0);
         }
         set_xy(propertyFlag::T, propertyFlag::P,
                x0, y0, TolAbsT, TolAbsP, TolRel, TolRel);

From ead1d55a14da13e1fd3d532e92b3f4b4e3af5006 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Tue, 25 Aug 2020 09:11:01 -0500
Subject: [PATCH 10/17] [purefluid] Ensure phase of matter of saturated vapor
 is 'gas'

---
 src/thermo/PureFluidPhase.cpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp
index 7d4a9ea79c9..ca3f9d6a776 100644
--- a/src/thermo/PureFluidPhase.cpp
+++ b/src/thermo/PureFluidPhase.cpp
@@ -90,7 +90,7 @@ std::string PureFluidPhase::phaseOfMatter() const
         return "supercritical";
     } else if (m_sub->TwoPhase() == 1) {
         return "liquid-gas-mix";
-    } else if (pressure() < m_sub->Ps()) {
+    } else if (m_sub->x() == 1) {
         return "gas";
     } else {
         return "liquid";

From 69125d14988becfc63b70d244db2fa723ac5db24 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Tue, 25 Aug 2020 08:03:14 -0500
Subject: [PATCH 11/17] [purefluid] Fix calculation of fd properties at
 saturation lines

- TP setters cannot be used at saturated conditions
- Update logic for saturated properties that are based on finite differences
---
 .../cython/cantera/test/test_purefluid.py     |  25 ++-
 src/tpx/Sub.cpp                               | 198 +++++++++++-------
 2 files changed, 142 insertions(+), 81 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index ad516c5a317..22683214621 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -180,9 +180,9 @@ def test_quality_exceptions(self):
         self.assertNear(self.water.T, self.water.critical_temperature)
 
         # Supercritical
-        with self.assertRaisesRegex(ct.CanteraError, 'above the critical point'):
+        with self.assertRaisesRegex(ct.CanteraError, 'supercritical'):
             self.water.TQ = 1.001 * self.water.critical_temperature, 0.
-        with self.assertRaisesRegex(ct.CanteraError, 'above the critical point'):
+        with self.assertRaisesRegex(ct.CanteraError, 'supercritical'):
             self.water.PQ = 1.001 * self.water.critical_pressure, 0.
 
         # Q negative
@@ -202,6 +202,27 @@ def test_saturated_mixture(self):
         with self.assertRaisesRegex(ct.CanteraError, 'Saturated mixture detected'):
             self.water.TP = 300, self.water.P_sat
 
+        # Saturated vapor
+        self.water.TQ = 373.15, 1.
+        self.assertEqual(self.water.phase_of_matter, 'gas')
+        self.assertFalse(np.isinf(self.water.cp_mass))
+        self.assertFalse(np.isinf(self.water.thermal_expansion_coeff))
+        self.assertFalse(np.isinf(self.water.isothermal_compressibility))
+
+        # Saturated mixture
+        self.water.TQ = 373.15, .5
+        self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
+        self.assertTrue(np.isinf(self.water.cp_mass))
+        self.assertTrue(np.isinf(self.water.thermal_expansion_coeff))
+        self.assertTrue(np.isinf(self.water.isothermal_compressibility))
+
+        # Saturated liquid
+        self.water.TQ = 373.15, 0.
+        self.assertEqual(self.water.phase_of_matter, 'liquid')
+        self.assertFalse(np.isinf(self.water.cp_mass))
+        self.assertFalse(np.isinf(self.water.thermal_expansion_coeff))
+        self.assertFalse(np.isinf(self.water.isothermal_compressibility))
+
     def test_saturation_near_limits(self):
         # Low temperature limit (triple point)
         self.water.TP = 300, ct.one_atm
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index f324cd91838..07a9880cb75 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -85,117 +85,155 @@ double Substance::cv()
 
 double Substance::cp()
 {
-    double Tsave = T, dt = 1.e-4*T;
-    double RhoSave = Rho;
-    double T1 = std::max(Tmin(), Tsave - dt);
-    double T2 = std::min(Tmax(), Tsave + dt);
-    double p0 = P();
-    double x0 = x();
     if (TwoPhase()) {
         // In the two-phase region, cp is infinite
         return std::numeric_limits<double>::infinity();
     }
 
-    Set(PropertyPair::TP, T1, p0);
-    double x1 = x();
-    if (x1 != x0) {
-        // If the initial state was pure liquid or pure vapor, and the state at
-        // T-dT is not, just take a one-sided difference
-        T1 = Tsave;
-        Set(PropertyPair::TP, T1, p0);
-    }
-    double s1 = s();
+    double Tsave = T, dt = 1.e-4 * T;
+    double RhoSave = Rho;
+    double T1, T2, s1, s2;
+    double p0 = P();
 
-    Set(PropertyPair::TP, T2, p0);
-    double x2 = x();
-    if (x2 != x0) {
-        // If the initial state was pure liquid or pure vapor, and the state at
-        // T+dT is not, just take a one-sided difference
-        T2 = Tsave;
+    if (Rho >= Rhf) {
+        // initial state is pure liquid
+        T1 = std::max(Tmin(), Tsave - dt);
+        Set(PropertyPair::TP, T1, p0);
+        s1 = s();
+        try {
+            T2 = std::min(Tsat(p0), Tsave + dt);
+        } catch (CanteraError&) {
+            // Tsat does not exist beyond two-phase region
+            T2 = Tsave + dt;
+        }
+        if (T2 < Tsave + dt) {
+            Set(PropertyPair::TX, T2, 0.);
+        } else {
+            Set(PropertyPair::TP, T2, p0);
+        }
+        s2 = s();
+    } else {
+        // initial state is pure vapor
+        try {
+            T1 = std::max(Tsat(p0), Tsave - dt);
+        } catch (CanteraError&) {
+            // Tsat does not exist beyond two-phase region
+            T1 = Tsave - dt;
+        }
+        if (T1 > Tsave - dt) {
+            Set(PropertyPair::TX, T1, 1.);
+        } else {
+            Set(PropertyPair::TP, T1, p0);
+        }
+        s1 = s();
+        T2 = std::min(Tmax(), Tsave + dt);
         Set(PropertyPair::TP, T2, p0);
+        s2 = s();
     }
-    double s2 = s();
 
     Set(PropertyPair::TV, Tsave, 1.0 / RhoSave);
-    return T*(s2 - s1)/(T2-T1);
+    return T * (s2 - s1) / (T2 - T1);
 }
 
 double Substance::thermalExpansionCoeff()
 {
-    double Tsave = T, dt = 1.e-4*T;
-    double RhoSave = Rho;
-    double T1 = std::max(Tmin(), Tsave - dt);
-    double T2 = std::min(Tmax(), Tsave + dt);
-    double p0 = P();
-    double x0 = x();
-
     if (TwoPhase()) {
         // In the two-phase region, the thermal expansion coefficient is
         // infinite
         return std::numeric_limits<double>::infinity();
     }
 
-    Set(PropertyPair::TP, T1, p0);
-    double x1 = x();
-    if (x1 != x0) {
-        // If the initial state was pure liquid or pure vapor, and the state at
-        // T-dT is not, just take a one-sided difference
-        T1 = Tsave;
-        Set(PropertyPair::TP, T1, p0);
-    }
-    double v1 = v();
+    double Tsave = T, dt = 1.e-4 * T;
+    double RhoSave = Rho;
+    double T1, T2, v1, v2;
+    double p0 = P();
 
-    Set(PropertyPair::TP, T2, p0);
-    double x2 = x();
-    if (x2 != x0) {
-        // If the initial state was pure liquid or pure vapor, and the state at
-        // T+dT is not, just take a one-sided difference
-        T2 = Tsave;
+    if (Rho >= Rhf) {
+        // initial state is pure liquid
+        T1 = std::max(Tmin(), Tsave - dt);
+        Set(PropertyPair::TP, T1, p0);
+        v1 = v();
+        try {
+            T2 = std::min(Tsat(p0), Tsave + dt);
+        } catch (CanteraError&) {
+            // Tsat does not exist beyond two-phase region
+            T2 = Tsave + dt;
+        }
+        if (T2 < Tsave + dt) {
+            Set(PropertyPair::TX, T2, 0.);
+        } else {
+            Set(PropertyPair::TP, T2, p0);
+        }
+        v2 = v();
+    } else {
+        // initial state is pure vapor
+        try {
+            T1 = std::max(Tsat(p0), Tsave - dt);
+        } catch (CanteraError&) {
+            // Tsat does not exist beyond two-phase region
+            T1 = Tsave - dt;
+        }
+        if (T1 > Tsave - dt) {
+            Set(PropertyPair::TX, T1, 1.);
+        } else {
+            Set(PropertyPair::TP, T1, p0);
+        }
+        v1 = v();
+        T2 = std::min(Tmax(), Tsave + dt);
         Set(PropertyPair::TP, T2, p0);
+        v2 = v();
     }
-    double v2 = v();
 
     Set(PropertyPair::TV, Tsave, 1.0 / RhoSave);
-    return 2.0*(v2 - v1)/((v2 + v1)*(T2-T1));
+    return 2.0 * (v2 - v1) / ((v2 + v1) * (T2 - T1));
 }
 
 double Substance::isothermalCompressibility()
 {
-    double Psave = P(), dp = 1.e-4*Psave;
-    double RhoSave = Rho;
-    double x0 = x();
-
     if (TwoPhase()) {
         // In the two-phase region, the isothermal compressibility is infinite
         return std::numeric_limits<double>::infinity();
     }
 
+    double Psave = P(), dp = 1.e-4 * Psave;
+    double RhoSave = Rho;
+    double P1, P2, v1, v2;
     double v0 = v();
-    double P1 = Psave - dp;
-    double P2 = Psave + dp;
 
-    Set(PropertyPair::TP, T, P1);
-    double x1 = x();
-    if (x1 != x0) {
-        // If the initial state was pure liquid or pure vapor, and the state at
-        // P-dP is not, just take a one-sided difference
-        P1 = Psave;
-        Set(PropertyPair::TP, T, P1);
-    }
-    double v1 = v();
-
-    Set(PropertyPair::TP, T, P2);
-    double x2 = x();
-    if (x2 != x0) {
-        // If the initial state was pure liquid or pure vapor, and the state at
-        // P+dP is not, just take a one-sided difference
-        P2 = Psave;
+    if (Rho >= Rhf) {
+        // initial state is pure liquid
+        P1 = Psave - dp;
+        if (T > Tmin() && T <= Tcrit()) {
+            P1 = std::max(Ps(), P1);
+        }
+        if (P1 > Psave - dp) {
+            Set(PropertyPair::PX, P1, 0.);
+        } else {
+            Set(PropertyPair::TP, T, P1);
+        }
+        v1 = v();
+        P2 = Psave + dp;
         Set(PropertyPair::TP, T, P2);
+        v2 = v();
+    } else {
+        // initial state is pure vapor
+        P1 = std::max(1.e-7, Psave - dp);
+        Set(PropertyPair::TP, T, P1);
+        v1 = v();
+        P2 = Psave + dp;
+        if (T > Tmin() && T <= Tcrit()) {
+            P2 = std::min(Ps(), P2);
+        }
+        if (P2 < Psave + dp) {
+            Set(PropertyPair::PX, P2, 1.);
+        } else {
+            Set(PropertyPair::TP, T, P2);
+        }
+        v2 = v();
     }
-    double v2 = v();
 
     Set(PropertyPair::TV, T, 1.0 / RhoSave);
-    return -(v2 - v1)/(v0*(P2-P1));
+    return -(v2 - v1) / (v0 * (P2 - P1));
 }
 
 double Substance::dPsdT()
@@ -377,13 +415,14 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
                x0, y0, TolAbsS, TolAbsH, TolRel, TolRel);
         break;
     case PropertyPair::PX:
+        temp = Tmin();
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (x0 > Pcrit()) {
+        } else if (x0 < Ps() || x0 > Pcrit()) {
             throw CanteraError("Substance::Set",
-                               "Vapor fraction cannot be set above the "
-                               "critical point");
+                               "Illegal pressure value: {} (supercritical "
+                               "or below triple line)", x0);
         } else {
             temp = Tsat(x0);
             set_T(temp);
@@ -395,10 +434,10 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (x0 > Tcrit()) {
+        } else if (x0 < Tmin() || x0 > Tcrit()) {
             throw CanteraError("Substance::Set",
-                               "Vapor fraction cannot be set above the "
-                               "critical point");
+                               "Illegal temperature value: {} "
+                               "(supercritical of below triple line)", x0);
         } else {
             set_T(x0);
             update_sat();
@@ -802,10 +841,11 @@ void Substance::set_TPp(double Temp, double Pressure)
         }
         Set(PropertyPair::TV, Temp, v_here);
         LoopCount++;
-        if (LoopCount > 100) {
+        if (LoopCount > 200) {
             Set(PropertyPair::TV, Temp, v_save);
             throw CanteraError("Substance::set_TPp",
-                               "No convergence for P* = {}, V* = {}",
+                               "No convergence for P = {}, T = {}\n"
+                               "(P* = {}, V* = {})", Pressure, Temp,
                                Pressure/Pcrit(), v_save/Vcrit());
         }
     }

From b6b6b662bd76e276850005816d04beaffcb5e3f7 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Wed, 26 Aug 2020 18:11:11 -0500
Subject: [PATCH 12/17] [purefluid] Disable calculation of cv for two-phase
 region

---
 src/tpx/Sub.cpp | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 07a9880cb75..818edbff5bd 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -54,6 +54,13 @@ const double DeltaT = 0.000001;
 
 double Substance::cv()
 {
+    if (TwoPhase()) {
+        // While cv can be calculated for the two-phase region (the state can
+        // always be continuously varied along an isochor on a T-v diagram),
+        // this calculation is currently not implemented
+        return std::numeric_limits<double>::quiet_NaN();
+    }
+
     double Tsave = T, dt = 1.e-4*T;
     double x0 = x();
     double T1 = std::max(Tmin(), Tsave - dt);

From ae570b274bde80a6500f0edc777aa490651c96a2 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Wed, 26 Aug 2020 17:40:55 -0500
Subject: [PATCH 13/17] [unittests] Update tests for fd properties at saturated
 conditions

---
 .../cython/cantera/test/test_purefluid.py     | 31 ++++++++++---------
 1 file changed, 16 insertions(+), 15 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 22683214621..082d4749f93 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -164,12 +164,6 @@ def test_thermal_expansion_coeff_TD(self):
             self.water.TD = T, 0.1
             self.assertNear(T * self.water.thermal_expansion_coeff, 1.0, 1e-2)
 
-    def test_fd_properties_twophase(self):
-        self.water.TQ = 400, 0.1
-        self.assertEqual(self.water.cp, np.inf)
-        self.assertEqual(self.water.isothermal_compressibility, np.inf)
-        self.assertEqual(self.water.thermal_expansion_coeff, np.inf)
-
     def test_quality_exceptions(self):
         # Critical point
         self.water.TP = 300, ct.one_atm
@@ -202,26 +196,33 @@ def test_saturated_mixture(self):
         with self.assertRaisesRegex(ct.CanteraError, 'Saturated mixture detected'):
             self.water.TP = 300, self.water.P_sat
 
+        w = ct.Water()
+
         # Saturated vapor
         self.water.TQ = 373.15, 1.
         self.assertEqual(self.water.phase_of_matter, 'gas')
-        self.assertFalse(np.isinf(self.water.cp_mass))
-        self.assertFalse(np.isinf(self.water.thermal_expansion_coeff))
-        self.assertFalse(np.isinf(self.water.isothermal_compressibility))
+        w.TP = self.water.T, .999 * self.water.P_sat
+        self.assertNear(self.water.cp, w.cp, 1.e-3)
+        self.assertNear(self.water.cv, w.cv, 1.e-3)
+        self.assertNear(self.water.thermal_expansion_coeff, w.thermal_expansion_coeff, 1.e-3)
+        self.assertNear(self.water.isothermal_compressibility, w.isothermal_compressibility, 1.e-3)
 
         # Saturated mixture
         self.water.TQ = 373.15, .5
         self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
-        self.assertTrue(np.isinf(self.water.cp_mass))
-        self.assertTrue(np.isinf(self.water.thermal_expansion_coeff))
-        self.assertTrue(np.isinf(self.water.isothermal_compressibility))
+        self.assertEqual(self.water.cp, np.inf)
+        self.assertTrue(np.isnan(self.water.cv))
+        self.assertEqual(self.water.isothermal_compressibility, np.inf)
+        self.assertEqual(self.water.thermal_expansion_coeff, np.inf)
 
         # Saturated liquid
         self.water.TQ = 373.15, 0.
         self.assertEqual(self.water.phase_of_matter, 'liquid')
-        self.assertFalse(np.isinf(self.water.cp_mass))
-        self.assertFalse(np.isinf(self.water.thermal_expansion_coeff))
-        self.assertFalse(np.isinf(self.water.isothermal_compressibility))
+        w.TP = self.water.T, 1.001 * self.water.P_sat
+        self.assertNear(self.water.cp, w.cp, 1.e-3)
+        self.assertNear(self.water.cv, w.cv, 1.e-3)
+        self.assertNear(self.water.thermal_expansion_coeff, w.thermal_expansion_coeff, 1.e-3)
+        self.assertNear(self.water.isothermal_compressibility, w.isothermal_compressibility, 1.e-3)
 
     def test_saturation_near_limits(self):
         # Low temperature limit (triple point)

From 24ec89be3db3189e95da1c7f1879d5571b2fe2be Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Thu, 27 Aug 2020 09:54:34 -0500
Subject: [PATCH 14/17] [purefluid] Disable low pressure check for
 Substance::Tsat

The calculation of Substance::Tsat should check for pressures below the
triple point and raise a corresponding exception. A bug causing low
temperature TP setting of heptane (GitHub issue #605) to fail prevents an
implementation of this check at the moment.
---
 interfaces/cython/cantera/test/test_purefluid.py | 7 ++++---
 src/tpx/Sub.cpp                                  | 9 +++++++--
 2 files changed, 11 insertions(+), 5 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 082d4749f93..2dbcf7ea5e2 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -259,9 +259,10 @@ def test_saturation_near_limits(self):
         with self.assertRaisesRegex(ct.CanteraError, 'Illegal temperature'):
             self.water.TP = .999 * self.water.min_temp, ct.one_atm
             self.water.P_sat
-        with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
-            self.water.TP = 300, .999 * psat
-            self.water.T_sat
+        # Test disabled pending fix of GitHub issue #605
+        # with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
+        #     self.water.TP = 300, .999 * psat
+        #     self.water.T_sat
 
     def test_TPQ(self):
         self.water.TQ = 400, 0.8
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 818edbff5bd..e3b5489d04b 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -290,8 +290,13 @@ double Substance::x()
 double Substance::Tsat(double p)
 {
     double Tsave = T;
-    T = Tmin();
-    if (p < Ps() || p > Pcrit()) {
+    if (p <= 0. || p > Pcrit()) {
+        // TODO: this check should also fail below the triple point pressure
+        // (calculated as Ps() for Tmin() as it is not available as a numerical
+        // parameter). A bug causing low-temperature TP setters to fail for
+        // Heptane (GitHub issue #605) prevents this at the moment.
+        // Without this check, a calculation attempt fails with the less
+        // meaningful error "No convergence" below.
         throw CanteraError("Substance::Tsat",
                            "Illegal pressure value: {}", p);
     }

From 2a70b7abe073b2a59c7d09a9d9c4846956708379 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 31 Aug 2020 01:48:42 -0500
Subject: [PATCH 15/17] [purefluid] Make pure saturated liquid/vapor part of
 TwoPhase region

Saturated liquid and saturated vapor form the boundary of the vapor dome and
are included in the phase-of-matter labeled 'liquid-vapor-mix'. At the same
time, properties that are based on finite difference calculations of a pure
liquid or pure vapor can be calculated. This approach corresponds to the
implementation of the CoolProp package (http://www.coolprop.org).
---
 include/cantera/tpx/Sub.h                       |  6 ++++--
 .../cython/cantera/test/test_purefluid.py       |  4 ++--
 src/thermo/PureFluidPhase.cpp                   |  2 +-
 src/tpx/Sub.cpp                                 | 17 ++++++++++-------
 4 files changed, 17 insertions(+), 12 deletions(-)

diff --git a/include/cantera/tpx/Sub.h b/include/cantera/tpx/Sub.h
index e5e3841b3ef..c445a56e8f0 100644
--- a/include/cantera/tpx/Sub.h
+++ b/include/cantera/tpx/Sub.h
@@ -146,8 +146,10 @@ class Substance
     //! is returned if v > Vcrit.
     double x();
 
-    //! Returns 1 if the current state is a liquid/vapor mixture, 0 otherwise
-    int TwoPhase();
+    //! Returns 1 if the current state is a liquid/vapor mixture, 0 otherwise.
+    //! By default, saturated vapor and saturated liquid are included; setting
+    //! the flag *strict* to true will exclude the boundaries.
+    int TwoPhase(bool strict=false);
     //! @}
 
     virtual double Pp()=0;
diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 2dbcf7ea5e2..db7e8ac59f2 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -200,7 +200,7 @@ def test_saturated_mixture(self):
 
         # Saturated vapor
         self.water.TQ = 373.15, 1.
-        self.assertEqual(self.water.phase_of_matter, 'gas')
+        self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
         w.TP = self.water.T, .999 * self.water.P_sat
         self.assertNear(self.water.cp, w.cp, 1.e-3)
         self.assertNear(self.water.cv, w.cv, 1.e-3)
@@ -217,7 +217,7 @@ def test_saturated_mixture(self):
 
         # Saturated liquid
         self.water.TQ = 373.15, 0.
-        self.assertEqual(self.water.phase_of_matter, 'liquid')
+        self.assertEqual(self.water.phase_of_matter, 'liquid-gas-mix')
         w.TP = self.water.T, 1.001 * self.water.P_sat
         self.assertNear(self.water.cp, w.cp, 1.e-3)
         self.assertNear(self.water.cv, w.cv, 1.e-3)
diff --git a/src/thermo/PureFluidPhase.cpp b/src/thermo/PureFluidPhase.cpp
index ca3f9d6a776..7d4a9ea79c9 100644
--- a/src/thermo/PureFluidPhase.cpp
+++ b/src/thermo/PureFluidPhase.cpp
@@ -90,7 +90,7 @@ std::string PureFluidPhase::phaseOfMatter() const
         return "supercritical";
     } else if (m_sub->TwoPhase() == 1) {
         return "liquid-gas-mix";
-    } else if (m_sub->x() == 1) {
+    } else if (pressure() < m_sub->Ps()) {
         return "gas";
     } else {
         return "liquid";
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index e3b5489d04b..6326178ffc6 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -54,14 +54,14 @@ const double DeltaT = 0.000001;
 
 double Substance::cv()
 {
-    if (TwoPhase()) {
+    if (TwoPhase(true)) {
         // While cv can be calculated for the two-phase region (the state can
         // always be continuously varied along an isochor on a T-v diagram),
         // this calculation is currently not implemented
         return std::numeric_limits<double>::quiet_NaN();
     }
 
-    double Tsave = T, dt = 1.e-4*T;
+    double Tsave = T, dt = 1.e-4 * T;
     double x0 = x();
     double T1 = std::max(Tmin(), Tsave - dt);
     double T2 = std::min(Tmax(), Tsave + dt);
@@ -92,7 +92,7 @@ double Substance::cv()
 
 double Substance::cp()
 {
-    if (TwoPhase()) {
+    if (TwoPhase(true)) {
         // In the two-phase region, cp is infinite
         return std::numeric_limits<double>::infinity();
     }
@@ -144,7 +144,7 @@ double Substance::cp()
 
 double Substance::thermalExpansionCoeff()
 {
-    if (TwoPhase()) {
+    if (TwoPhase(true)) {
         // In the two-phase region, the thermal expansion coefficient is
         // infinite
         return std::numeric_limits<double>::infinity();
@@ -197,7 +197,7 @@ double Substance::thermalExpansionCoeff()
 
 double Substance::isothermalCompressibility()
 {
-    if (TwoPhase()) {
+    if (TwoPhase(true)) {
         // In the two-phase region, the isothermal compressibility is infinite
         return std::numeric_limits<double>::infinity();
     }
@@ -260,13 +260,16 @@ double Substance::dPsdT()
     return dpdt;
 }
 
-int Substance::TwoPhase()
+int Substance::TwoPhase(bool strict)
 {
     if (T >= Tcrit()) {
         return 0;
     }
     update_sat();
-    return ((Rho < Rhf) && (Rho > Rhv) ? 1 : 0);
+    if (strict) {
+        return ((Rho > Rhv) && (Rho < Rhf) ? 1 : 0);
+    }
+    return ((Rho >= Rhv) && (Rho <= Rhf) ? 1 : 0);
 }
 
 double Substance::x()

From 04ec6dcccd2cd08bbd0823f0d950fb053c037958 Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Sun, 6 Dec 2020 11:17:22 -0600
Subject: [PATCH 16/17] [purefluid] Update PQ triple point checks

---
 .../cython/cantera/test/test_purefluid.py     | 23 +++++++++++++---
 src/tpx/Sub.cpp                               | 27 ++++++++++---------
 2 files changed, 33 insertions(+), 17 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index db7e8ac59f2..e5bd38042de 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -118,10 +118,10 @@ def check_fd_properties(self, T1, P1, T2, P2, tol):
         self.assertNear(k1, k2, tol)
         self.assertNear(alpha1, alpha2, tol)
 
-        # calculating these finite difference properties should not perturbe the
-        # state of the object
-        self.assertEqual(h1a, h1b)
-        self.assertEqual(h2a, h2b)
+        # calculating these finite difference properties should not perturb the
+        # state of the object (except for checks on edge cases)
+        self.assertNear(h1a, h1b, 1e-9)
+        self.assertNear(h2a, h2b, 1e-9)
 
     def test_properties_near_min(self):
         self.check_fd_properties(self.water.min_temp*(1+1e-5), 101325,
@@ -164,6 +164,21 @@ def test_thermal_expansion_coeff_TD(self):
             self.water.TD = T, 0.1
             self.assertNear(T * self.water.thermal_expansion_coeff, 1.0, 1e-2)
 
+    def test_pq_setter_triple_check(self):
+        self.water.PQ = 101325, .2
+        T = self.water.T
+        # change T such that it would result in a Psat larger than P
+        self.water.TP = 400, 101325
+        # ensure that correct triple point pressure is recalculated
+        # (necessary as this value is not stored by the C++ base class)
+        self.water.PQ = 101325, .2
+        self.assertNear(T, self.water.T, 1e-9)
+        with self.assertRaisesRegex(ct.CanteraError, 'below triple point'):
+            # min_temp is triple point temperature
+            self.water.TP = self.water.min_temp, 101325
+            P = self.water.P_sat # triple-point pressure
+            self.water.PQ = .999*P, .2
+
     def test_quality_exceptions(self):
         # Critical point
         self.water.TP = 300, ct.one_atm
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 6326178ffc6..430722ef39d 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -431,33 +431,34 @@ void Substance::Set(PropertyPair::type XY, double x0, double y0)
         break;
     case PropertyPair::PX:
         temp = Tmin();
+        set_T(temp);
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (x0 < Ps() || x0 > Pcrit()) {
+        }
+        if (x0 < Ps() || x0 > Pcrit()) {
             throw CanteraError("Substance::Set",
                                "Illegal pressure value: {} (supercritical "
-                               "or below triple line)", x0);
-        } else {
-            temp = Tsat(x0);
-            set_T(temp);
-            update_sat();
-            Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
+                               "or below triple point)", x0);
         }
+        temp = Tsat(x0);
+        set_T(temp);
+        update_sat();
+        Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
         break;
     case PropertyPair::TX:
         if (y0 > 1.0 || y0 < 0.0) {
             throw CanteraError("Substance::Set",
                                "Invalid vapor fraction, {}", y0);
-        } else if (x0 < Tmin() || x0 > Tcrit()) {
+        }
+        if (x0 < Tmin() || x0 > Tcrit()) {
             throw CanteraError("Substance::Set",
                                "Illegal temperature value: {} "
-                               "(supercritical of below triple line)", x0);
-        } else {
-            set_T(x0);
-            update_sat();
-            Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
+                               "(supercritical or below triple point)", x0);
         }
+        set_T(x0);
+        update_sat();
+        Rho = 1.0/((1.0 - y0)/Rhf + y0/Rhv);
         break;
     default:
         throw CanteraError("Substance::Set", "Invalid input.");

From 864caa04bba3c5b798557dacac5a22cdcb1f370d Mon Sep 17 00:00:00 2001
From: Ingmar Schoegl <ischoegl@lsu.edu>
Date: Mon, 18 Jan 2021 18:15:55 -0600
Subject: [PATCH 17/17] Mark @TODO items

---
 interfaces/cython/cantera/test/test_purefluid.py | 2 +-
 src/tpx/Sub.cpp                                  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index e5bd38042de..dc910f12ceb 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -274,7 +274,7 @@ def test_saturation_near_limits(self):
         with self.assertRaisesRegex(ct.CanteraError, 'Illegal temperature'):
             self.water.TP = .999 * self.water.min_temp, ct.one_atm
             self.water.P_sat
-        # Test disabled pending fix of GitHub issue #605
+        # @TODO: test disabled pending fix of GitHub issue #605
         # with self.assertRaisesRegex(ct.CanteraError, 'Illegal pressure value'):
         #     self.water.TP = 300, .999 * psat
         #     self.water.T_sat
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 430722ef39d..b5d88373f0a 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -294,7 +294,7 @@ double Substance::Tsat(double p)
 {
     double Tsave = T;
     if (p <= 0. || p > Pcrit()) {
-        // TODO: this check should also fail below the triple point pressure
+        // @TODO: this check should also fail below the triple point pressure
         // (calculated as Ps() for Tmin() as it is not available as a numerical
         // parameter). A bug causing low-temperature TP setters to fail for
         // Heptane (GitHub issue #605) prevents this at the moment.