From deeaaed03f619e936b10e908bf9835e58dc10da7 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Sat, 2 Dec 2017 16:49:34 -0500
Subject: [PATCH] [Reactor] Fix using pure substances near temperature limits

Fixes #475
---
 .../cython/cantera/test/test_reactor.py       | 26 +++++++++++++++++++
 src/zeroD/Reactor.cpp                         | 17 ++++++++----
 2 files changed, 38 insertions(+), 5 deletions(-)

diff --git a/interfaces/cython/cantera/test/test_reactor.py b/interfaces/cython/cantera/test/test_reactor.py
index fe840dd11..209180559 100644
--- a/interfaces/cython/cantera/test/test_reactor.py
+++ b/interfaces/cython/cantera/test/test_reactor.py
@@ -1468,6 +1468,32 @@ class PureFluidReactorTest(utilities.CanteraTest):
         self.assertEqual(states.X[-1], 1)
         self.assertNear(states.X[30], 0.54806, 1e-4)
 
+    def test_Reactor_2(self):
+        phase = ct.PureFluid('liquidvapor.xml', 'carbondioxide')
+        air = ct.Solution('air.xml')
+
+        phase.TP = 218, 5e6
+        r1 = ct.Reactor(phase)
+        r1.volume = 0.1
+
+        air.TP = 500, 5e6
+        r2 = ct.Reactor(air)
+        r2.volume = 10.0
+
+        w1 = ct.Wall(r1, r2, U=10000, A=1)
+        w1.expansion_rate_coeff = 1e-3
+        net = ct.ReactorNet([r1,r2])
+
+        states = ct.SolutionArray(phase, extra='t')
+        for t in np.arange(0.0, 60.0, 1):
+            net.advance(t)
+            states.append(TD=r1.thermo.TD, t=net.time)
+
+        self.assertEqual(states.X[0], 0)
+        self.assertEqual(states.X[-1], 1)
+        self.assertNear(states.X[20], 0.644865, 1e-4)
+
+
     def test_ConstPressureReactor(self):
         phase = ct.Nitrogen()
         air = ct.Solution('air.xml')
diff --git a/src/zeroD/Reactor.cpp b/src/zeroD/Reactor.cpp
index 82c87b4a0..abb1301b4 100644
--- a/src/zeroD/Reactor.cpp
+++ b/src/zeroD/Reactor.cpp
@@ -144,13 +144,20 @@ void Reactor::updateState(doublereal* y)
             TT = bmt::bracket_and_solve_root(
                 u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
         } catch (std::exception& err) {
-            // Set m_thermo back to a reasonable state if root finding fails
-            m_thermo->setState_TR(T, m_mass / m_vol);
-            throw CanteraError("Reactor::updateState", err.what());
+            // Try full-range bisection if bracketing fails (e.g. near
+            // temperature limits for the phase's equation of state)
+            try {
+                TT = bmt::bisect(u_err, m_thermo->minTemp(), m_thermo->maxTemp(),
+                    bmt::eps_tolerance<double>(48), maxiter);
+            } catch (std::exception& err2) {
+                // Set m_thermo back to a reasonable state if root finding fails
+                m_thermo->setState_TR(T, m_mass / m_vol);
+                throw CanteraError("Reactor::updateState",
+                    "{}\nat U = {}, rho = {}", err2.what(), U, m_mass / m_vol);
+            }
         }
         if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
-            throw CanteraError("Reactor::updateState",
-                "bracket_and_solve_root failed");
+            throw CanteraError("Reactor::updateState", "root finding failed");
         }
         m_thermo->setState_TR(TT.second, m_mass / m_vol);
     } else {