Various robustness additions to the Davis EOS

CHANGELOG.md

@@ -21,6 +21,7 @@
 - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver
 - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface
 - [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions
+- [[PR374]](https://github.com/lanl/singularity-eos/pull/374) Make the Davis EOS more numerically robust
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls

singularity-eos/eos/eos_davis.hpp

@@ -18,8 +18,10 @@
 #include <cmath>
 #include <cstdio>
 
+#include <ports-of-call/portable_errors.hpp>
 #include <singularity-eos/base/constants.hpp>
 #include <singularity-eos/base/math_utils.hpp>
+#include <singularity-eos/base/robust_utils.hpp>
 #include <singularity-eos/base/root-finding-1d/root_finding.hpp>
 #include <singularity-eos/eos/eos_base.hpp>
 
@@ -41,17 +43,21 @@ class DavisReactants : public EosBase<DavisReactants> {
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    const Real es = Es(rho);
-    const Real tmp = std::pow((1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0,
-                              1.0 / (1.0 + _alpha));
-    if (tmp > 0) return Ts(rho) * tmp;
-    return Ts(rho) + (sie - es) / _Cv0; // This branch is a negative temperature
+    const Real power_base = DimlessEdiff(rho, sie);
+    if (power_base <= 0) {
+      // This case would result in an imaginary temperature (i.e. negative), but we won't
+      // allow that so return zero
+      return 0.;
+    }
+    const Real tmp = std::pow(power_base, 1.0 / (1.0 + _alpha));
+    return Ts(rho) * tmp;
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature(
       const Real rho, const Real temp,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     const Real t_s = Ts(rho);
+    PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided");
     return Es(rho) +
            _Cv0 * t_s / (1.0 + _alpha) * (std::pow(temp / t_s, 1.0 + _alpha) - 1.0);
   }
@@ -72,8 +78,11 @@ class DavisReactants : public EosBase<DavisReactants> {
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real MinInternalEnergyFromDensity(
       const Real rho, Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    MinInternalEnergyIsNotEnabled("DavisReactants");
-    return 0.0;
+    // Minimum enegy is when the returned temperature is zero. This only happens
+    // when the base to the exponent is zero (see T(rho, e) equation)
+    const Real es = Es(rho);
+    const Real ts = Ts(rho);
+    return es - (_Cv0 * ts) / (1 + _alpha);
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real
@@ -100,8 +109,12 @@ class DavisReactants : public EosBase<DavisReactants> {
   PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy(
       const Real rho, const Real sie,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    return _Cv0 / std::pow((1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1,
-                           -_alpha / (1 + _alpha));
+    const Real power_base = DimlessEdiff(rho, sie);
+    if (power_base <= 0) {
+      // Return zero heat capacity instead of an imaginary value
+      return 0.;
+    }
+    return _Cv0 / std::pow(power_base, -_alpha / (1 + _alpha));
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature(
@@ -166,6 +179,7 @@ class DavisReactants : public EosBase<DavisReactants> {
   // static constexpr const char _eos_type[] = "DavisReactants";
   static constexpr unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
+  PORTABLE_FORCEINLINE_FUNCTION Real DimlessEdiff(const Real rho, const Real sie) const;
   PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const;
   PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const;
   PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const;
@@ -300,22 +314,34 @@ class DavisProducts : public EosBase<DavisProducts> {
   static constexpr const unsigned long _preferred_input =
       thermalqs::density | thermalqs::specific_internal_energy;
   PORTABLE_INLINE_FUNCTION Real F(const Real rho) const {
+    if (rho <= 0) {
+      return 0.;
+    }
     const Real vvc = 1.0 / (rho * _vc);
     return 2.0 * _a / (std::pow(vvc, 2 * _n) + 1.0);
   }
   PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const {
+    if (rho <= 0) {
+      return 0.;
+    }
     const Real vvc = 1 / (rho * _vc);
     return _pc * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) /
            std::pow(vvc, _k + _a) * (_k - 1.0 + F(rho)) / (_k - 1.0 + _a);
   }
   PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const {
+    if (rho <= 0) {
+      return 0.;
+    }
     const Real vvc = 1 / (rho * _vc);
     const Real ec = _pc * _vc / (_k - 1.0 + _a);
     // const Real de = ecj-(Es(rho0)-_E0);
     return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) /
            std::pow(vvc, _k - 1.0 + _a);
   }
   PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const {
+    if (rho <= 0) {
+      return 0.;
+    }
     const Real vvc = 1 / (rho * _vc);
     return std::pow(2.0, -_a * _b / _n) * _pc * _vc / (_Cv * (_k - 1 + _a)) *
            std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n * (1 - _b)) /
@@ -326,9 +352,14 @@ class DavisProducts : public EosBase<DavisProducts> {
   }
 };
 
+PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho,
+                                                                const Real sie) const {
+  return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0;
+}
+
 PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const {
   using namespace math_utils;
-  const Real y = 1.0 - _rho0 / rho;
+  const Real y = 1.0 - robust::ratio(_rho0, std::max(rho, 0.));
   const Real phat = 0.25 * _A * _A / _B * _rho0;
   const Real b4y = 4.0 * _B * y;
 
@@ -342,7 +373,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const {
   }
 }
 PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const {
-  const Real y = 1 - _rho0 / rho;
+  const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.));
   const Real phat = 0.25 * _A * _A / _B * _rho0;
   const Real b4y = 4 * _B * y;
   Real e_s;
@@ -354,19 +385,17 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const {
   } else {
     e_s = -y - (1.0 - std::exp(b4y)) / (4.0 * _B);
   }
-  return _e0 + _P0 * (1.0 / _rho0 - 1.0 / rho) + phat / _rho0 * e_s;
+  return _e0 + _P0 * (1.0 / _rho0 - robust::ratio(1.0, std::max(rho, 0.))) +
+         phat / _rho0 * e_s;
 }
 PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const {
-  if (rho >= _rho0) {
-    const Real y = 1 - _rho0 / rho;
-    return _T0 * std::exp(-_Z * y) * std::pow(_rho0 / rho, -_G0 - _Z);
-  } else {
-    return _T0 * std::pow(_rho0 / rho, -_G0);
-  }
+  const Real rho0overrho = robust::ratio(_rho0, std::max(rho, 0.));
+  const Real y = 1 - rho0overrho;
+  return _T0 * std::exp(-_Z * y) * std::pow(rho0overrho, -_G0 - _Z);
 }
 PORTABLE_INLINE_FUNCTION Real DavisReactants::Gamma(const Real rho) const {
   if (rho >= _rho0) {
-    const Real y = 1 - _rho0 / rho;
+    const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.));
     return _G0 + _Z * y;
   } else {
     return _G0;
@@ -377,30 +406,33 @@ template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEnergy(
     const Real rho, const Real sie, Indexer_t &&lambda) const {
   using namespace math_utils;
-  const Real y = 1 - _rho0 / rho;
+  const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.));
   const Real phat = 0.25 * _A * _A / _B * _rho0;
   const Real b4y = 4 * _B * y;
-  const Real gamma = Gamma(rho);
-  const Real esv = -Ps(rho);
+  const Real gamma = Gamma(std::max(rho, 0.));
+  const Real esv = -Ps(std::max(rho, 0.));
   const Real psv =
       (rho >= _rho0)
           ? -phat * _rho0 *
                 (4 * _B * (1 + b4y + 0.5 * (pow<2>(b4y) + _C / 3 * pow<3>(b4y))) +
                  3 * y / pow<4>(1 - y) + 4 * pow<2>(y) / pow<5>(1 - y))
           : -phat * 4 * _B * _rho0 * std::exp(b4y);
   const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0;
-  return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * rho) - gamma * rho * esv) /
-         rho;
+  const Real numerator =
+      -(psv + (sie - Es(rho)) * std::max(rho, 0.) * (gammav - gamma * std::max(rho, 0.)) -
+        gamma * std::max(rho, 0.) * esv);
+  return robust::ratio(numerator, std::max(rho, 0.));
 }
 
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperature(
     const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
   // First, solve P=P(rho,T) for rho.  Note P(rho,e) has an sie-es term, which is only a
   // function of T
+  PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided");
   auto PofRatT = [&](const Real r) {
     return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) *
-                        (std::pow(temp / Ts(r), 1 + _alpha) - 1.0));
+                        (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0));
   };
   using RootFinding1D::regula_falsi;
   using RootFinding1D::Status;
@@ -410,6 +442,10 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu
     EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: "
               "Root find failed to find a solution given P, T\n");
   }
+  if (rho < 0.) {
+    EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: "
+              "Root find resulted in a negative density\n");
+  }
   sie = InternalEnergyFromDensityTemperature(rho, temp);
 }
 
@@ -418,6 +454,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::FillEos(Real &rho, Real &temp, Rea
                                                       Real &press, Real &cv, Real &bmod,
                                                       const unsigned long output,
                                                       Indexer_t &&lambda) const {
+  // BROKEN: This can only handle density-energy inputs! MUST BE CHANGED
   if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie);
   if (output & thermalqs::temperature)
     temp = TemperatureFromDensityInternalEnergy(rho, sie);
@@ -450,6 +487,9 @@ template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnergy(
     const Real rho, const Real sie, Indexer_t &&lambda) const {
   using namespace math_utils;
+  if (rho <= 0) {
+    return 0.;
+  }
   const Real vvc = 1 / (rho * _vc);
   const Real Fx = -4 * _a * std::pow(vvc, 2 * _n - 1) / pow<2>(1 + std::pow(vvc, 2 * _n));
   const Real tmp = std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) /
@@ -470,6 +510,7 @@ PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnerg
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperature(
     const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
+  PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided");
   auto PofRatT = [&](const Real r) {
     return (Ps(r) + Gamma(r) * r * _Cv * (temp - Ts(r)));
   };
@@ -482,12 +523,17 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur
     EOS_ERROR("DavisProducts::DensityEnergyFromPressureTemperature: "
               "Root find failed to find a solution given P, T\n");
   }
+  if (rho < 0.) {
+    EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: "
+              "Root find resulted in a negative density\n");
+  }
   sie = InternalEnergyFromDensityTemperature(rho, temp);
 }
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void
 DavisProducts::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
                        Real &bmod, const unsigned long output, Indexer_t &&lambda) const {
+  // BROKEN: This can only handle density-energy inputs! MUST BE CHANGED
   if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie);
   if (output & thermalqs::temperature)
     temp = TemperatureFromDensityInternalEnergy(rho, sie);