[Reactor] Use Boost's bracket_and_solve_root function

This algorithm is more robust than the simple Newton's method that it
replaces. Among its advantages is that it works with PureFluidPhase models.

Resolves #257

src/zeroD/Reactor.cpp

@@ -10,9 +10,10 @@
 #include "cantera/zeroD/ReactorNet.h"
 #include "cantera/zeroD/ReactorSurface.h"
 
-#include <cfloat>
+#include <boost/math/tools/roots.hpp>
 
 using namespace std;
+namespace bmt = boost::math::tools;
 
 namespace Cantera
 {
@@ -139,38 +140,29 @@ void Reactor::updateState(doublereal* y)
     m_thermo->setMassFractions_NoNorm(y+3);
 
     if (m_energy) {
-        // Use a damped Newton's method to determine the mixture temperature.
+        double U = y[2];
-        // Tight tolerances are required both for Jacobian evaluation and for
+        // Residual function: error in internal energy as a function of T
-        // sensitivity analysis to work correctly.
+        auto u_err = [this, U](double T) {
-        doublereal U = y[2];
-        doublereal T = temperature();
-        double dT = 100;
-        double dUprev = 1e10;
-        double dU = 1e10;
-        int i = 0;
-        double damp = 1.0;
-        while (abs(dT / T) > 10 * DBL_EPSILON) {
-            dUprev = dU;
             m_thermo->setState_TR(T, m_mass / m_vol);
-            double dUdT = m_thermo->cv_mass() * m_mass;
+            return m_thermo->intEnergy_mass() * m_mass - U;
-            dU = m_thermo->intEnergy_mass() * m_mass - U;
+        };
-            dT = dU / dUdT;
+
-            // Reduce the damping coefficient if the magnitude of the error
+        double T = m_thermo->temperature();
-            // isn't decreasing
+        boost::uintmax_t maxiter = 100;
-            if (std::abs(dU) < std::abs(dUprev)) {
+        std::pair<double, double> TT;
-                damp = 1.0;
+        try {
-            } else {
+            TT = bmt::bracket_and_solve_root(
-                damp *= 0.8;
+                u_err, T, 1.2, true, bmt::eps_tolerance<double>(48), maxiter);
-            }
+        } catch (std::exception& err) {
-            dT = std::min(dT, 0.5 * T) * damp;
+            // Set m_thermo back to a reasonable state if root finding fails
-            T -= dT;
+            m_thermo->setState_TR(T, m_mass / m_vol);
-            i++;
+            throw CanteraError("Reactor::updateState", err.what());
-            if (i > 100) {
-                throw CanteraError("Reactor::updateState",
-                    "no convergence\nU/m = {}\nT = {}\nrho = {}\n",
-                    U / m_mass, T, m_mass / m_vol);
-            }
         }
+        if (fabs(TT.first - TT.second) > 1e-7*TT.first) {
+            throw CanteraError("Reactor::updateState",
+                "bracket_and_solve_root failed");
+        }
+        m_thermo->setState_TR(TT.second, m_mass / m_vol);
     } else {
         m_thermo->setDensity(m_mass/m_vol);
     }