From da8726006279fa4e02872dd41ce08731b5e63993 Mon Sep 17 00:00:00 2001
From: Ray Speth <speth@mit.edu>
Date: Thu, 3 Nov 2016 22:47:18 -0400
Subject: [PATCH] [Thermo] Replace RootFind with Boost in RedlichKwongMFTP

As previously implemented, calling densSpinodalLiquid or densSpinodalGas would
crash since the pointer fdpdv_ was uninitialized, and furthermore could not even
be initialized with an instance of the spinodalFunc class because this class was
abstract.

Also fix densSpinodalLiquid and densSpinodalGas to actually calculate the roots
of the equation of state which are used as bounds for the search interval.
---
 include/cantera/thermo/MixtureFugacityTP.h | 10 ----
 src/thermo/MixtureFugacityTP.cpp           | 14 -----
 src/thermo/RedlichKwongMFTP.cpp            | 60 +++++++++++-----------
 3 files changed, 31 insertions(+), 53 deletions(-)

diff --git a/include/cantera/thermo/MixtureFugacityTP.h b/include/cantera/thermo/MixtureFugacityTP.h
index e2b9fc1ed..de9aa1d87 100644
--- a/include/cantera/thermo/MixtureFugacityTP.h
+++ b/include/cantera/thermo/MixtureFugacityTP.h
@@ -551,14 +551,6 @@ protected:
 
     //@}
 
-    class spinodalFunc : public ResidEval
-    {
-    public:
-        spinodalFunc(MixtureFugacityTP* tp);
-        virtual int evalSS(const doublereal t, const doublereal* const y, doublereal* const r);
-        MixtureFugacityTP* m_tp;
-    };
-
 protected:
     virtual void invalidateCache();
 
@@ -604,8 +596,6 @@ protected:
 
     //! Temporary storage for dimensionless reference state entropies
     mutable vector_fp m_s0_R;
-
-    spinodalFunc* fdpdv_;
 };
 }
 
diff --git a/src/thermo/MixtureFugacityTP.cpp b/src/thermo/MixtureFugacityTP.cpp
index 91b70ce93..96eafe94c 100644
--- a/src/thermo/MixtureFugacityTP.cpp
+++ b/src/thermo/MixtureFugacityTP.cpp
@@ -546,20 +546,6 @@ void MixtureFugacityTP::updateMixingExpressions()
 {
 }
 
-MixtureFugacityTP::spinodalFunc::spinodalFunc(MixtureFugacityTP* tp) :
-    m_tp(tp)
-{
-}
-
-int MixtureFugacityTP::spinodalFunc::evalSS(const doublereal t, const doublereal* const y,
-        doublereal* const r)
-{
-    doublereal molarVol = y[0];
-    doublereal pp;
-    r[0] = m_tp->dpdVCalc(m_tp->temperature(), molarVol, pp);
-    return 0;
-}
-
 int MixtureFugacityTP::corr0(doublereal TKelvin, doublereal pres, doublereal& densLiqGuess,
                              doublereal& densGasGuess, doublereal& liqGRT, doublereal& gasGRT)
 {
diff --git a/src/thermo/RedlichKwongMFTP.cpp b/src/thermo/RedlichKwongMFTP.cpp
index 67d850778..73a739542 100644
--- a/src/thermo/RedlichKwongMFTP.cpp
+++ b/src/thermo/RedlichKwongMFTP.cpp
@@ -6,11 +6,13 @@
 #include "cantera/thermo/RedlichKwongMFTP.h"
 #include "cantera/thermo/mix_defs.h"
 #include "cantera/thermo/ThermoFactory.h"
-#include "cantera/numerics/RootFind.h"
 #include "cantera/base/stringUtils.h"
 #include "cantera/base/ctml.h"
 
+#include <boost/math/tools/roots.hpp>
+
 using namespace std;
+namespace bmt = boost::math::tools;
 
 namespace Cantera
 {
@@ -899,46 +901,46 @@ doublereal RedlichKwongMFTP::densityCalc(doublereal TKelvin, doublereal presPa,
 
 doublereal RedlichKwongMFTP::densSpinodalLiquid() const
 {
-    if (NSolns_ != 3) {
+    double Vroot[3];
+    double T = temperature();
+    int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
+    if (nsol != 3) {
         return critDensity();
     }
-    double vmax = Vroot_[1];
-    double vmin = Vroot_[0];
-    RootFind rf(fdpdv_);
-    rf.setPrintLvl(10);
-    rf.setTol(1.0E-5, 1.0E-10);
-    rf.setFuncIsGenerallyDecreasing(true);
 
-    double vbest = 0.5 * (Vroot_[0]+Vroot_[1]);
-    double funcNeeded = 0.0;
-    int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
-    if (status != ROOTFIND_SUCCESS) {
-        throw CanteraError("  RedlichKwongMFTP::densSpinodalLiquid() ", "didn't converge");
-    }
+    auto resid = [this, T](double v) {
+        double pp;
+        return dpdVCalc(T, v, pp);
+    };
+
+    boost::uintmax_t maxiter = 100;
+    std::pair<double, double> vv = bmt::toms748_solve(
+        resid, Vroot[0], Vroot[1], bmt::eps_tolerance<double>(48), maxiter);
+
     doublereal mmw = meanMolecularWeight();
-    return mmw / vbest;
+    return mmw / (0.5 * (vv.first + vv.second));
 }
 
 doublereal RedlichKwongMFTP::densSpinodalGas() const
 {
-    if (NSolns_ != 3) {
+    double Vroot[3];
+    double T = temperature();
+    int nsol = NicholsSolve(T, pressure(), m_a_current, m_b_current, Vroot);
+    if (nsol != 3) {
         return critDensity();
     }
-    double vmax = Vroot_[2];
-    double vmin = Vroot_[1];
-    RootFind rf(fdpdv_);
-    rf.setPrintLvl(10);
-    rf.setTol(1.0E-5, 1.0E-10);
-    rf.setFuncIsGenerallyIncreasing(true);
 
-    double vbest = 0.5 * (Vroot_[1]+Vroot_[2]);
-    double funcNeeded = 0.0;
-    int status = rf.solve(vmin, vmax, 100, funcNeeded, &vbest);
-    if (status != ROOTFIND_SUCCESS) {
-        throw CanteraError(" RedlichKwongMFTP::densSpinodalGas() ", "didn't converge");
-    }
+    auto resid = [this, T](double v) {
+        double pp;
+        return dpdVCalc(T, v, pp);
+    };
+
+    boost::uintmax_t maxiter = 100;
+    std::pair<double, double> vv = bmt::toms748_solve(
+        resid, Vroot[1], Vroot[2], bmt::eps_tolerance<double>(48), maxiter);
+
     doublereal mmw = meanMolecularWeight();
-    return mmw / vbest;
+    return mmw / (0.5 * (vv.first + vv.second));
 }
 
 doublereal RedlichKwongMFTP::pressureCalc(doublereal TKelvin, doublereal molarVol) const