Merge branch 'main' of github.com:lanl/singularity-eos into jhp/MinEn…

…ergyBound

CHANGELOG.md

@@ -20,9 +20,12 @@
 - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS
 - [[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
+- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Removed E0 from Davis Products EOS in favor of using the shifted EOS modifier. CHANGES API!
 
 ### Infrastructure (changes irrelevant to downstream codes)
 - [[PR329]](https://github.com/lanl/singularity-eos/pull/329) Move vinet tests into analytic test suite

doc/sphinx/src/models.rst

@@ -1254,17 +1254,44 @@ Here, there are four dimensionless parameters that are settable by the user,
 :math:`e_\mathrm{C}`, :math:`V_\mathrm{C}` and :math:`T_\mathrm{C}` are tuning
 parameters with units related to their non-subscripted counterparts.
 
+Note that the energy zero (i.e. the reference energy) for the Davis products EOS
+is arbitrary. For the isentrope to properly pass through the CJ state of a
+reacting material, the energy release of the reaction needs to be accounted for
+properly. If done external to the EOS, an energy source term is required in the
+Euler equations. However, a common convention is to specify the reactants and
+product EOS in a consistent way such that the reference energy corresponds to
+the rest state of the material *before* it reacts.
+
+The energy at the CJ state can be calculated as
+
+.. math::
+
+    e_\mathrm{CJ} = \frac{P_0 + P_\mathrm{CJ}}{2(V_0 - V_\mathrm{CJ})},
+
+relative to :math:`e = 0` at the reference state of the *reactants*. Therefore
+the energy offset of the products EOS is given by
+
+.. math::
+
+    e_0 = e_S(V_\mathrm{CJ}) - e_\mathrm{CJ}.
+
+Practically, this means :math:`e_0` should be positive for any energetic material.
+
+To provide the energy offset to the Davis Products EOS, `the energy shift
+modifier<modifiers shifted EOS>`_ should be used. Note that the convention there
+is that the shift is positive, so :math:`-e_0` should be provided to the shift
+modifier.
+
 The constructor for the Davis Products EOS is
 
 .. code-block:: cpp
 
   DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
-                const Real pc, const Real Cv, const Real E0)
+                const Real pc, const Real Cv)
 
 where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`,
 ``n`` is :math:`n`, ``vc`` is :math:`V_\mathrm{C}`, ``pc`` is
-:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`, and ``E0`` is
-:math:`e_\mathrm{C}`.
+:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`.
 
 Spiner EOS
 ````````````

doc/sphinx/src/modifiers.rst

@@ -22,6 +22,8 @@ We list below the available modifiers and their constructors.
 The Shifted EOS
 -----------------
 
+.. _modifiers shifted EOS:
+
 The shifted equation of state modifies zero point energy of an
 underlying model by some shift. So for example, it transforms
 

python/module.cpp

@@ -60,9 +60,9 @@ PYBIND11_MODULE(singularity_eos, m) {
   eos_class<DavisProducts>(m, "DavisProducts")
     .def(py::init())
     .def(
-      py::init<Real, Real, Real, Real, Real, Real, Real, Real>(),
+      py::init<Real, Real, Real, Real, Real, Real, Real>(),
       py::arg("a"), py::arg("b"), py::arg("k"), py::arg("n"), py::arg("vc"),
-      py::arg("pc"), py::arg("Cv"), py::arg("E0")
+      py::arg("pc"), py::arg("Cv")
     );
 
 #ifdef SPINER_USE_HDF

singularity-eos/closure/mixed_cell_models.hpp

@@ -182,6 +182,7 @@ class PTESolverBase {
   void Finalize() {
     for (int m = 0; m < nmat; m++) {
       temp[m] *= Tnorm;
+      PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned");
       u[m] *= uscale;
       press[m] *= uscale;
     }

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;
@@ -177,8 +191,8 @@ class DavisProducts : public EosBase<DavisProducts> {
   DavisProducts() = default;
   PORTABLE_INLINE_FUNCTION
   DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc,
-                const Real pc, const Real Cv, const Real E0)
-      : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _E0(E0) {}
+                const Real pc, const Real Cv)
+      : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {}
   DavisProducts GetOnDevice() { return *this; }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy(
@@ -286,36 +300,48 @@ class DavisProducts : public EosBase<DavisProducts> {
   static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; }
   PORTABLE_INLINE_FUNCTION void PrintParams() const {
     static constexpr char s1[]{"DavisProducts Params: "};
-    printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc,
-           _pc, _Cv, _E0);
+    printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc,
+           _Cv);
   }
   inline void Finalize() {}
   static std::string EosType() { return std::string("DavisProducts"); }
   static std::string EosPyType() { return EosType(); }
 
  private:
   static constexpr Real onethird = 1.0 / 3.0;
-  Real _a, _b, _k, _n, _vc, _pc, _Cv, _E0;
+  Real _a, _b, _k, _n, _vc, _pc, _Cv;
   // static constexpr const char _eos_type[] = "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);

singularity-eos/eos/get_sg_eos_p_t.cpp

@@ -39,6 +39,10 @@ void get_sg_eos_p_t(const char *name, int ncell, int nmat, indirection_v &offset
         // for small loops
         const int32_t token{tokens.acquire()};
         const int32_t tid{small_loop ? iloop : token};
+        // need to initialize the scratch before it's used to avoid undefined behavior
+        for (int idx = 0; idx < solver_scratch.extent(1); ++idx) {
+          solver_scratch(tid, idx) = 0.0;
+        }
         // caching mechanism
         singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0));
         double mass_sum{0.0};

singularity-eos/eos/get_sg_eos_rho_e.cpp

@@ -39,6 +39,10 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v,
         int npte{0};
         // initialize values for solver / lookup
         i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0);
+        // need to initialize the scratch before it's used to avoid undefined behavior
+        for (int idx = 0; idx < solver_scratch.extent(1); ++idx) {
+          solver_scratch(tid, idx) = 0.0;
+        }
         // get cache from offsets into scratch
         const int neq = npte + 1;
         singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) +