[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

@@ -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_;
 };
 }
 

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)
 {

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