From 3f3b96587b22ecd6b484dfda4310a1c3c7a356b1 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 20 Oct 2016 15:50:08 -0400
Subject: [PATCH] [Thermo] Fix finite-difference properties in two-phase region

These properties are actually infinite in the two-phase region, but attempting
to compute them by finite difference would incorrectly give a finite result, so
they need to be treated as a special case.
---
 .../cython/cantera/test/test_purefluid.py     |  7 +++++
 src/tpx/Sub.cpp                               | 28 +++++++++++++++----
 2 files changed, 29 insertions(+), 6 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_purefluid.py b/interfaces/cython/cantera/test/test_purefluid.py
index 6660cae2e..b9da60a7c 100644
--- a/interfaces/cython/cantera/test/test_purefluid.py
+++ b/interfaces/cython/cantera/test/test_purefluid.py
@@ -132,6 +132,13 @@ class TestPureFluid(utilities.CanteraTest):
         self.water.TP = 450, 12
         self.assertNear(ref.thermal_expansion_coeff,
                         self.water.thermal_expansion_coeff, 1e-5)
+
+    def test_fd_properties_twophase(self):
+        self.water.TX = 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_TPX(self):
         self.water.TX = 400, 0.8
         T,P,X = self.water.TPX
diff --git a/src/tpx/Sub.cpp b/src/tpx/Sub.cpp
index 6e6f448e2..dd66db232 100644
--- a/src/tpx/Sub.cpp
+++ b/src/tpx/Sub.cpp
@@ -90,10 +90,14 @@ double Substance::cp()
     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 ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
+    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;
@@ -103,7 +107,7 @@ double Substance::cp()
 
     Set(PropertyPair::TP, T2, p0);
     double x2 = x();
-    if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
+    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;
@@ -123,9 +127,15 @@ double Substance::thermalExpansionCoeff()
     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 ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
+    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;
@@ -135,7 +145,7 @@ double Substance::thermalExpansionCoeff()
 
     Set(PropertyPair::TP, T2, p0);
     double x2 = x();
-    if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
+    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;
@@ -151,13 +161,19 @@ double Substance::isothermalCompressibility()
 {
     double Psave = P(), dp = 1.e-4*Psave;
     double x0 = x();
+
+    if (TwoPhase()) {
+        // In the two-phase region, the isothermal compressibility is infinite
+        return std::numeric_limits<double>::infinity();
+    }
+
     double v0 = v();
     double P1 = Psave - dp;
     double P2 = Psave + dp;
 
     Set(PropertyPair::TP, T, P1);
     double x1 = x();
-    if ((x0 == 1.0 || x0 == 0.0) && x1 != x0) {
+    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;
@@ -167,7 +183,7 @@ double Substance::isothermalCompressibility()
 
     Set(PropertyPair::TP, T, P2);
     double x2 = x();
-    if ((x0 == 1.0 || x0 == 0.0) && x2 != x0) {
+    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;