Merge pull request #304 from lanl/hls/newton-raphson

Newton-Raphson root find for Helmholtz EOS

CHANGELOG.md

@@ -22,6 +22,7 @@
 - [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables
 
 ### Added (new features/APIs/variables/...)
+- [[PR304]](https://github.com/lanl/singularity-eos/pull/304) added a Newton-Raphson root find for use with the Helmholtz EoS
 - [[PR265]](https://github.com/lanl/singularity-eos/pull/265) Add missing UnitSystem modifier combinations to variant and EOSPAC
 - [[PR279]](https://github.com/lanl/singularity-eos/pull/279) added noble-abel EoS
 - [[PR274]](https://github.com/lanl/singularity-eos/pull/274) added a stiffened gas EoS

doc/sphinx/src/models.rst

@@ -1399,14 +1399,18 @@ The constructor for the ``Helmholtz`` EOS class looks like
   Helmholtz(const std::string &filename, const bool rad = true,
             const bool gas = true, const bool coul = true,
             const bool ion = true, const bool ele = true,
-            const bool verbose = false)
+            const bool verbose = false, const bool newton_raphson = true)
 
 where ``filename`` is the file containing the tabulated model. The
 optional arguments ``rad``, ``gas``, ``coul``, ``ion``, and ``ele``
 specify whether to include the radiation, ideal gas, coulomb correction,
 ionization, and electron contributions, respectively. The default is to
 include all terms. The optional argument ``verbose`` specifies whether to print
-out additional information, e.g. when the root find fails to converge.
+out additional information, e.g. when the root find fails to converge. The
+optional argument ``newton_raphson`` specifies whether to use the Newton-Raphson
+method or the regula falsi method for root finding. The default is to use the
+Newton-Raphson method (note that the regula falsi method is used as a fallback
+in case the Newton-Raphson method does not converge).
 
 EOSPAC EOS
 ````````````

python/module.cpp

@@ -108,7 +108,9 @@ PYBIND11_MODULE(singularity_eos, m) {
   eos_class<Helmholtz>(m, "Helmholtz")
     .def(py::init())
     .def(py::init<std::string&>())
-    .def(py::init<std::string&,bool,bool,bool,bool,bool>());
+    .def(py::init<std::string&,bool,bool,bool,bool,bool>())
+    .def(py::init<std::string&,bool,bool,bool,bool,bool,bool>())
+    .def(py::init<std::string&,bool,bool,bool,bool,bool,bool,bool>());
 #endif
 #endif
 

singularity-eos/base/root-finding-1d/root_finding.hpp

@@ -28,6 +28,7 @@
 #include <math.h>
 #include <ports-of-call/portability.hpp>
 #include <stdio.h>
+#include <tuple>
 
 #define SINGULARITY_ROOT_DEBUG (0)
 #define SINGULARITY_ROOT_VERBOSE (0)
@@ -38,6 +39,7 @@ namespace RootFinding1D {
 constexpr const int SECANT_NITER_MAX{1000};
 constexpr const int BISECT_NITER_MAX{1000};
 constexpr const int BISECT_REG_MAX{1000};
+constexpr const int NEWTON_RAPHSON_NITER_MAX{100};
 enum class Status { SUCCESS = 0, FAIL = 1 };
 
 /*
@@ -235,6 +237,76 @@ PORTABLE_INLINE_FUNCTION Status regula_falsi(const T &f, const Real ytarget,
   return status;
 }
 
+// solves for f(x,params) - ytarget = 0
+// WARNING: this root finding expects a different callable f than the other
+// root finding methods. f should return a tuple of (f(x), f'(x)) where f'(x)
+// is the derivative of f with respect to x.
+template <typename T>
+PORTABLE_INLINE_FUNCTION Status newton_raphson(const T &f, const Real ytarget,
+                                               const Real guess, const Real a,
+                                               const Real b, const Real ytol, Real &xroot,
+                                               const RootCounts *counts = nullptr,
+                                               const bool &verbose = false,
+                                               const bool &fail_on_bound_root = true) {
+
+  constexpr int max_iter = NEWTON_RAPHSON_NITER_MAX;
+  Real _x = guess;
+  Real _xold = 0.0;
+  auto status = Status::SUCCESS;
+
+  Real yg;
+  Real dfunc;
+
+  int iter;
+
+  for (iter = 0; iter < max_iter; iter++) {
+    std::tie(yg, dfunc) = f(_x); // C++11 tuple unpacking
+
+    // check if we are converged already
+    if (std::abs(yg - ytarget) < (ytol * ytarget)) break;
+
+    // not converged; compute the next step
+    _xold = _x;
+    _x = _x - (yg - ytarget) / dfunc;
+
+    // check if we are out of bounds
+    // CAUTION: we do not set the root to the boundary value in this case
+    // because one might want to handle this on a case by case basis
+    // (e.g. if the boundary is a physical boundary, then one might want to
+    // set the root to the boundary value).
+    // Per default, we fail if the root is out of bounds controlled by
+    // fail_on_bound_root.
+    if ((_x <= a && _xold <= a) || (_x >= b && _xold >= b)) {
+      if (verbose) {
+        printf("newton_raphson out of bounds! %.14e %.14e %.14e %.14e\n", ytarget, guess,
+               a, b);
+      }
+      if (fail_on_bound_root) {
+        status = Status::FAIL;
+      }
+      break;
+    }
+    _x = std::max(std::min(_x, b), a);
+  }
+  if (iter >= max_iter) {
+    if (verbose) {
+      printf("root finding reached the maximum number of iterations.  likely not "
+             "converged\n");
+    }
+    status = Status::FAIL;
+  }
+
+  if (counts != nullptr) {
+    if (iter < counts->nBins()) {
+      counts->increment(iter);
+    } else {
+      counts->increment(counts->more());
+    }
+  }
+  xroot = _x;
+  return status;
+}
+
 template <typename T>
 PORTABLE_INLINE_FUNCTION Status findRoot(const T &f, const Real ytarget, Real xguess,
                                          const Real xmin, const Real xmax,

singularity-eos/eos/eos_helmholtz.hpp

@@ -32,6 +32,7 @@
 #include <fstream>
 #include <sstream>
 #include <string>
+#include <tuple>
 #include <utility>
 
 // ports of call
@@ -49,8 +50,6 @@
 // spiner
 #include <spiner/databox.hpp>
 
-#define ROOT_FINDER (RootFinding1D::regula_falsi)
-
 /*
  * The singularity-eos implementation of the Helmhotlz equation of state
  * provided by Timmes, and Swesty
@@ -388,12 +387,18 @@ class Helmholtz : public EosBase<Helmholtz> {
             const bool electron, const bool verbose)
         : ENABLE_RAD(rad), ENABLE_GAS(gas), ENABLE_COULOMB_CORRECTIONS(coul),
           GAS_IONIZED(ion), GAS_DEGENERATE(electron), VERBOSE(verbose) {}
+    Options(const bool rad, const bool gas, const bool coul, const bool ion,
+            const bool electron, const bool verbose, const bool newton_raphson)
+        : ENABLE_RAD(rad), ENABLE_GAS(gas), ENABLE_COULOMB_CORRECTIONS(coul),
+          GAS_IONIZED(ion), GAS_DEGENERATE(electron), VERBOSE(verbose),
+          USE_NEWTON_RAPHSON(newton_raphson) {}
     bool ENABLE_RAD = true;
     bool ENABLE_GAS = true;
     bool ENABLE_COULOMB_CORRECTIONS = true;
     bool GAS_IONIZED = true;
     bool GAS_DEGENERATE = true;
     bool VERBOSE = false;
+    bool USE_NEWTON_RAPHSON = true;
   };
 
   Helmholtz() = default;
@@ -406,6 +411,10 @@ class Helmholtz : public EosBase<Helmholtz> {
   Helmholtz(const std::string &filename, const bool rad, const bool gas, const bool coul,
             const bool ion, const bool ele, const bool verbose)
       : electrons_(filename), options_(rad, gas, coul, ion, ele, verbose) {}
+  Helmholtz(const std::string &filename, const bool rad, const bool gas, const bool coul,
+            const bool ion, const bool ele, const bool verbose, const bool newton_raphson)
+      : electrons_(filename),
+        options_(rad, gas, coul, ion, ele, verbose, newton_raphson) {}
 
   PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return 3; }
   static constexpr unsigned long PreferredInput() {
@@ -659,6 +668,7 @@ class Helmholtz : public EosBase<Helmholtz> {
                      const Real lDe, Real *lambda) const;
 
   static constexpr Real ROOT_THRESH = 1e-14;
+  static constexpr Real HELM_EOS_EPS = 1e-10;
   Options options_;
   HelmRad rad_;
   HelmIon ions_;
@@ -736,19 +746,51 @@ Real Helmholtz::lTFromRhoSie_(const Real rho, const Real e, const Real abar,
     }
     Real Tguess = math_utils::pow10(lTguess);
     auto &copy = *this; // stupid C++17 workaround
-    auto status = ROOT_FINDER(
-        [&](Real T) {
-          Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
-          copy.GetFromDensityLogTemperature_(rho, T, abar, zbar, ye, ytot, ywot, De, lDe,
-                                             p, e, s, etaele, nep, true);
-          return e[VAL];
-        },
-        e, Tguess, math_utils::pow10(electrons_.lTMin()),
-        math_utils::pow10(electrons_.lTMax()), ROOT_THRESH, ROOT_THRESH, T, nullptr,
-        options_.VERBOSE);
-    if (status != RootFinding1D::Status::SUCCESS) {
-      lT = lTAnalytic_(rho, e, ni, options_.GAS_IONIZED * ne);
-      T = math_utils::pow10(lT);
+    if (options_.USE_NEWTON_RAPHSON) {
+      auto status = RootFinding1D::newton_raphson(
+          [&](Real T) {
+            Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
+            copy.GetFromDensityLogTemperature_(rho, T, abar, zbar, ye, ytot, ywot, De,
+                                               lDe, p, e, s, etaele, nep, true);
+            return std::make_tuple(e[VAL], e[DDT]);
+          },
+          e, Tguess, math_utils::pow10(electrons_.lTMin()),
+          math_utils::pow10(electrons_.lTMax()), HELM_EOS_EPS, T, nullptr,
+          options_.VERBOSE, false);
+      if (status != RootFinding1D::Status::SUCCESS) {
+        if (options_.VERBOSE) {
+          printf("Newton-Raphson failed to converge, falling back to regula falsi\n");
+        }
+        status = RootFinding1D::regula_falsi(
+            [&](Real T) {
+              Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
+              copy.GetFromDensityLogTemperature_(rho, T, abar, zbar, ye, ytot, ywot, De,
+                                                 lDe, p, e, s, etaele, nep, true);
+              return e[VAL];
+            },
+            e, Tguess, math_utils::pow10(electrons_.lTMin()),
+            math_utils::pow10(electrons_.lTMax()), ROOT_THRESH, ROOT_THRESH, T, nullptr,
+            options_.VERBOSE);
+        if (status != RootFinding1D::Status::SUCCESS) {
+          lT = lTAnalytic_(rho, e, ni, options_.GAS_IONIZED * ne);
+          T = math_utils::pow10(lT);
+        }
+      }
+    } else {
+      auto status = RootFinding1D::regula_falsi(
+          [&](Real T) {
+            Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
+            copy.GetFromDensityLogTemperature_(rho, T, abar, zbar, ye, ytot, ywot, De,
+                                               lDe, p, e, s, etaele, nep, true);
+            return e[VAL];
+          },
+          e, Tguess, math_utils::pow10(electrons_.lTMin()),
+          math_utils::pow10(electrons_.lTMax()), ROOT_THRESH, ROOT_THRESH, T, nullptr,
+          options_.VERBOSE);
+      if (status != RootFinding1D::Status::SUCCESS) {
+        lT = lTAnalytic_(rho, e, ni, options_.GAS_IONIZED * ne);
+        T = math_utils::pow10(lT);
+      }
     }
   } else {
     lT = lTAnalytic_(rho, e, ni, options_.GAS_IONIZED * ne);