Make sure that DensityEnergyFromPressureTemperature works for all equations of state

.github/workflows/tests_minimal.yml

@@ -36,4 +36,4 @@ jobs:
                   ..
             make -j4
             make install
-            make test
+            ctest --output-on-failure

.github/workflows/tests_minimal_kokkos.yml

@@ -36,4 +36,4 @@ jobs:
                   ..
             make -j4
             make install
-            make test
+            ctest --output-on-failure

CHANGELOG.md

@@ -6,6 +6,7 @@
 - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS
 
 ### Fixed (Repair bugs, etc)
+- [[PR449]](https://github.com/lanl/singularity-eos/pull/449) Ensure that DensityEnergyFromPressureTemperature works for all equations of state and is properly tested
 - [[PR439]](https://github.com/lanl/singularity-eos/pull/439) Add mean atomic mass and number to EOS API
 - [[PR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test
 

doc/sphinx/src/using-eos.rst

@@ -1143,18 +1143,70 @@ quantities as outputs.
 Methods Used for Mixed Cell Closures
 --------------------------------------
 
-Several methods were developed in support of mixed cell closures. In particular:
+Several methods were developed in support of mixed cell closures. In particular the function
 
 .. cpp:function:: Real MinimumDensity() const;
 
-and 
+the function
 
 .. cpp:function:: Real MinimumTemperature() const;
 
+and the function
+
+.. cpp:function:: Real MaximumDensity() const;
+
 provide bounds for valid inputs into a table, which can be used by a
-root finder to meaningful bound the root search. Similarly,
+root finder to meaningful bound the root search.
+
+.. warning::
+
+  For unbounded equations of state, ``MinimumDensity`` and
+  ``MinimumTemperature`` will return zero, while ``MaximumDensity``
+  will return a very large finite number. Which number you get,
+  however, is not guaranteed. You may wish to apply more sensible
+  bounds in your own code.
+
+Similarly,
 
 .. cpp:function:: Real RhoPmin(const Real temp) const;
 
 returns the density at which pressure is minimized for a given
-temperature. This is again useful for root finds.
+temperature. The function
+
+.. cpp:function:: Real MinimumPressure() const;
+
+provides the minimum pressure an equation of state supports, which may
+be the most negative tension state. The function
+
+.. cpp:function:: Real MaximumPressureAtTemperature(const Real temp) const;
+
+provides a maximum possible pressure an equation of state supports at
+a given temperature. (Most models are unbounded in pressure.) This is
+again useful for root finds.
+
+The function
+
+.. code-block:: cpp
+
+  template <typename Indexer_t = Real*>
+  void DensityEnergyFromPressureTemperature(const Real press, const Real temp,
+                                            Indexer_t &&lambda, Real &rho, Real &sie) const;
+
+is designed for working in Pressure-Temperature space. Given a
+pressure ``press`` and temperature ``temp``, it sets a density ``rho``
+and specific internal energy ``sie``. The ``lambda`` is optional and
+defaults to a ``nullptr``.
+
+Typically this operation requires a root find. You may pass in an
+initial guess for the density ``rho`` in-place and most EOS models
+will use it.
+
+.. warning::
+
+  Pressure is not necessarily monotone in density and it may be double
+  valued. Thus you are not guaranteed to find the correct root and the
+  value of your initial guess may determine correctness. The fact that
+  ``rho`` may be used as an initial guess means you **must** pass in
+  an initialized variable, even if it is zero-initialized. Do not pass
+  uninitialized memory into this function!
+

eospac-wrapper/eospac_wrapper.cpp

@@ -25,9 +25,9 @@
 namespace EospacWrapper {
 
 void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) {
-  constexpr int NT = 2;
+  constexpr int NT = 3;
   EOS_INTEGER tableHandle[NT];
-  EOS_INTEGER tableType[NT] = {EOS_Info, EOS_Ut_DT};
+  EOS_INTEGER tableType[NT] = {EOS_Info, EOS_Ut_DT, EOS_Pt_DT};
 
   EOS_INTEGER commentsHandle[1];
   EOS_INTEGER commentsType[1] = {EOS_Comment};
@@ -48,7 +48,8 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) {
     }
   }
 
-  eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Info", "EOS_Ut_DT"}, eospacWarn);
+  eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Info", "EOS_Ut_DT", "EOS_Pt_DT"},
+              eospacWarn);
 
   for (int i = 0; i < numInfoTables; i++) {
     eosSafeTableInfo(&(tableHandle[i]), NI[i], infoItems[i].data(), infoVals[i].data(),
@@ -66,6 +67,8 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) {
   metadata.TMax = temperatureFromSesame(infoVals[1][3]);
   metadata.sieMin = sieFromSesame(infoVals[1][4]);
   metadata.sieMax = sieFromSesame(infoVals[1][5]);
+  metadata.PMin = pressureFromSesame(infoVals[2][4]);
+  metadata.PMax = pressureFromSesame(infoVals[2][5]);
   metadata.numRho = static_cast<int>(infoVals[1][6]);
   metadata.numT = static_cast<int>(infoVals[1][7]);
   metadata.rhoConversionFactor = infoVals[1][8];

eospac-wrapper/eospac_wrapper.hpp

@@ -59,6 +59,7 @@ class SesameMetadata {
   double rhoMin, rhoMax;
   double TMin, TMax;
   double sieMin, sieMax;
+  double PMin, PMax;
   double rhoConversionFactor;
   double TConversionFactor;
   double sieConversionFactor;

singularity-eos/eos/eos_base.hpp

@@ -16,11 +16,14 @@
 #define _SINGULARITY_EOS_EOS_EOS_BASE_
 
 #include <cstring>
+#include <limits>
 #include <string>
 
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_errors.hpp>
 #include <singularity-eos/base/constants.hpp>
+#include <singularity-eos/base/robust_utils.hpp>
+#include <singularity-eos/base/root-finding-1d/root_finding.hpp>
 #include <singularity-eos/base/variadic_utils.hpp>
 
 namespace singularity {
@@ -781,6 +784,23 @@ class EosBase {
   PORTABLE_FORCEINLINE_FUNCTION
   Real MinimumTemperature() const { return 0; }
 
+  // Report maximum value of density. Default is unbounded.
+  // JMM: Should we use actual infinity, the largest real, or just a
+  // big number?  For comparisons, actual infinity is better. It also
+  // has the advantage of being projective with modifiers that modify
+  // the max. On the other hand, it's more fraught if someone tries to
+  // put it into a formula without guarding against it.
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumDensity() const { return 1e100; }
+
+  // These are for the PT space PTE solver to bound the iterations in
+  // a safe range.
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumPressure() const { return 0; }
+  // Gruneisen EOS's often have a maximum density, which implies a maximum pressure.
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; }
+
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const { return 0.0; }
 
@@ -836,6 +856,58 @@ class EosBase {
     PORTABLE_ALWAYS_THROW_OR_ABORT(msg);
   }
 
+  // JMM: This method is often going to be overloaded for special cases.
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
+                                       Indexer_t &&lambda, Real &rho, Real &sie) const {
+    using RootFinding1D::findRoot; // more robust but slower. Better default.
+    using RootFinding1D::Status;
+
+    // Pressure is not monotone in density at low densities, which can
+    // prevent convergence. We want to approach tension from above,
+    // not below. Choose close to, but above, normal density for a
+    // metal like copper.
+    constexpr Real DEFAULT_RHO_GUESS = 12;
+
+    CRTP copy = *(static_cast<CRTP const *>(this));
+
+    // P(rho) not monotone. When relevant, bound rhopmin.
+    Real rhomin = std::max(copy.RhoPmin(temp), copy.MinimumDensity());
+    Real rhomax = copy.MaximumDensity();
+    PORTABLE_REQUIRE(rhomax > rhomin, "max bound > min bound");
+
+    auto PofRT = [&](const Real r) {
+      return copy.PressureFromDensityTemperature(r, temp, lambda);
+    };
+    Real rhoguess = rho;                                // use input density
+    if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { // avoid edge effects
+      if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) {
+        rhoguess = DEFAULT_RHO_GUESS;
+      } else {
+        rhoguess = 0.5 * (rhomin + rhomax);
+      }
+    }
+    auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(),
+                           robust::EPS(), rho);
+    // JMM: This needs to not fail and instead return something sane.
+    // If root find failed to converge, density will at least be
+    // within bounds.
+    if (status != Status::SUCCESS) {
+      PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n");
+    }
+    sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda);
+    return;
+  }
+  PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press,
+                                                                     const Real temp,
+                                                                     Real &rho,
+                                                                     Real &sie) const {
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    copy.DensityEnergyFromPressureTemperature(press, temp, static_cast<Real *>(nullptr),
+                                              rho, sie);
+  }
+
   // Serialization
   /*
     The methodology here is there are *three* size methods all EOS's provide:

singularity-eos/eos/eos_eospac.hpp

@@ -1104,6 +1104,7 @@ class EOSPAC : public EosBase<EOSPAC> {
   }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; }
   PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; }
+  PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return press_min_; }
 
  private:
   static constexpr const unsigned long _preferred_input =
@@ -1128,6 +1129,7 @@ class EOSPAC : public EosBase<EOSPAC> {
   Real dvdt_ref_ = 1;
   Real rho_min_ = 0;
   Real temp_min_ = 0;
+  Real press_min_ = 0;
   // TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a
   // problem? Could it ever realistically overflow?
   EOS_INTEGER shared_size_, packed_size_;
@@ -1199,6 +1201,8 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data,
   rho_ref_ = m.normalDensity;
   rho_min_ = m.rhoMin;
   temp_min_ = m.TMin;
+  press_min_ = m.PMin;
+
   // use std::max to hydrogen, in case of bad table
   AZbar_.Abar = std::max(1.0, m.meanAtomicMass);
   AZbar_.Zbar = std::max(1.0, m.meanAtomicNumber);

singularity-eos/eos/eos_gruneisen.hpp

@@ -170,6 +170,17 @@ class Gruneisen : public EosBase<Gruneisen> {
            s1, _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max);
     _AZbar.PrintParams();
   }
+
+  // Report minimum values of density and temperature
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumDensity() const { return _rho_max; }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumPressure() const { return PressureFromDensityInternalEnergy(0, 0); }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumPressureAtTemperature(const Real T) const {
+    return MaxStableDensityAtTemperature(T);
+  }
+
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,

singularity-eos/eos/eos_helmholtz.hpp

@@ -264,6 +264,12 @@ class HelmElectrons {
   std::size_t numTemp() const { return NTEMP; }
   PORTABLE_FORCEINLINE_FUNCTION
   std::size_t numRho() const { return NRHO; }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumDensity() const { return rhoMin(); }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumTemperature() const { return TMin_; }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumDensity() const { return rhoMax(); }
 
  private:
   inline void InitDataFile_(const std::string &filename);
@@ -676,16 +682,38 @@ class Helmholtz : public EosBase<Helmholtz> {
   ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
                          Real &bmod, Real &dpde, Real &dvdt,
                          Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
-    // JMM: I'm not sure what to put here or if it matters. Some
-    // reference state, maybe stellar denity, would be appropriate.
-    PORTABLE_ALWAYS_ABORT("Stub");
+    // JMM: Conditions for an oxygen burning shell in a stellar
+    // core. Not sure if that's the best choice.
+    rho = 1e7;
+    temp = 1.5e9;
+    FillEos(rho, temp, sie, press, cv, bmod,
+            thermalqs::specific_internal_energy | thermalqs::pressure |
+                thermalqs::specific_heat | thermalqs::bulk_modulus,
+            lambda);
+    dpde = GruneisenParamFromDensityTemperature(rho, temp, lambda) * rho;
+    dvdt = 0;
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   DensityEnergyFromPressureTemperature(const Real press, const Real temp,
                                        Indexer_t &&lambda, Real &rho, Real &sie) const {
-    // This is only used for mixed cell closures. Stubbing it out for now.
-    PORTABLE_ALWAYS_ABORT("Stub");
+    using RootFinding1D::regula_falsi;
+    using RootFinding1D::Status;
+    PORTABLE_REQUIRE(temp > = 0, "Non-negative temperature required");
+    auto PofRT = [&](const Real r) {
+      return PressureFromDensityTemperature(r, temp, lambda);
+    };
+    auto status = regula_falsi(PofRT, press, 1e7, electrons_.rhoMin(),
+                               electrons_.rhoMax(), 1.0e-8, 1.0e-8, rho);
+    if (status != Status::SUCCESS) {
+      PORTABLE_THROW_OR_ABORT("Helmholtz::DensityEnergyFromPressureTemperature: "
+                              "Root find failed to find a solution given P, T\n");
+    }
+    if (rho < 0.) {
+      PORTABLE_THROW_OR_ABORT("Helmholtz::DensityEnergyFromPressureTemperature: "
+                              "Root find produced negative energy solution\n");
+    }
+    sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
   }
   static std::string EosType() { return std::string("Helmholtz"); }
   static std::string EosPyType() { return EosType(); }