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
 

singularity-eos/eos/eos_base.hpp

@@ -21,6 +21,8 @@
 #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 {
@@ -836,6 +838,44 @@ 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 {
+    // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's
+    // overwrite.
+    constexpr Real MAXFAC = 1e8;
+    constexpr Real EPS = 10 * robust::EPS();
+    constexpr Real MINR_DEFAULT = 1e-4;
+    constexpr Real DEFAULT_RHO_GUESS = 4;
+    using RootFinding1D::findRoot; // more robust but slower. Better default.
+    using RootFinding1D::Status;
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    auto PofRT = [&](const Real r) {
+      return copy.PressureFromDensityTemperature(r, temp, lambda);
+    };
+    // JMM: This can't be zero, in case MinimumDensity is zero
+    Real rhomin = std::max(MINR_DEFAULT, copy.MinimumDensity());
+    Real rhomax = MAXFAC * rhomin;
+    Real rhoguess = rho; // use input density
+    if ((rhoguess < rhomin) || (rhoguess > rhomax)) {
+      if ((rhomin < DEFAULT_RHO_GUESS) && (rhomax > DEFAULT_RHO_GUESS)) {
+        rhoguess = DEFAULT_RHO_GUESS;
+      } else {
+        rhoguess = 0.5 * (rhomin + rhomax);
+      }
+    }
+    auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho);
+    // JMM: This needs to not fail and instead return something sane
+    if (status != Status::SUCCESS) {
+      PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n");
+      rho = rhoguess;
+    }
+    sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda);
+    return;
+  }
+
   // Serialization
   /*
     The methodology here is there are *three* size methods all EOS's provide:

singularity-eos/eos/eos_helmholtz.hpp

@@ -676,16 +676,37 @@ 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");
+    // TODO(JMM): Is this right???
+    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(); }

singularity-eos/eos/eos_mgusup.hpp

@@ -125,6 +125,47 @@ class MGUsup : public EosBase<MGUsup> {
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,
           const unsigned long output,
           Indexer_t &&lambda = static_cast<Real *>(nullptr)) const;
+
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumDensity() const { return 10 * robust::EPS(); }
+
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MaximumDensity() const {
+    if (_s > 1) {
+      return 0.99 * robust::ratio(_s * _rho0, _s - 1);
+    } else { // for s <= 1, no maximum, but we need to pick something.
+      return 1e3 * _rho0;
+    }
+  }
+
+  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;
+    using RootFinding1D::Status;
+    // JMM: diverges for rho -> 0
+    // and for s > 1 and rho -> rho0
+    Real rhomin = MinimumDensity();
+    Real rhomax = MaximumDensity();
+    auto PofRT = [&](const Real r) {
+      return PressureFromDensityTemperature(r, temp, lambda);
+    };
+    Real rhoguess = rho; // use input density
+    if ((rhoguess < rhomin) || (rhoguess > rhomax)) {
+      rhoguess = 0.5 * (rhomin + rhomax);
+    }
+    // JMM: regula_falsi does not respect bounds
+    auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(),
+                           robust::EPS(), rho);
+    if (status != Status::SUCCESS) {
+      PORTABLE_THROW_OR_ABORT(
+          "DensityEnergyFromPressureTemperature failed to find root\n");
+    }
+    sie = InternalEnergyFromDensityTemperature(rho, temp, lambda);
+    return;
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv,
@@ -148,11 +189,6 @@ class MGUsup : public EosBase<MGUsup> {
     _AZbar.PrintParams();
     printf("\n\n");
   }
-  // Density/Energy from P/T not unique, if used will give error
-  template <typename Indexer_t>
-  PORTABLE_INLINE_FUNCTION void
-  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                       Indexer_t &&lambda, Real &rho, Real &sie) const;
   inline void Finalize() {}
   static std::string EosType() { return std::string("MGUsup"); }
   static std::string EosPyType() { return EosType(); }
@@ -200,6 +236,8 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const
   Real eta = 1.0 - robust::ratio(_rho0, rho);
   Real f1 = 1.0 - _s * eta;
   if (f1 <= 0.0) {
+    printf("f1, eta, rho, rho0, s = %.14e %.14e %.14e %.14e %.14e\n", f1, eta, rho, _rho0,
+           _s);
     PORTABLE_ALWAYS_THROW_OR_ABORT("MGUsup model parameters s and rho0 together with rho "
                                    "give a negative argument for a logarithm.");
   }
@@ -350,13 +388,6 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityInternalEnergy(
   value = robust::ratio(_rho0, rho) * value;
   return value;
 }
-// AEM: Give error since function is not well defined
-template <typename Indexer_t>
-PORTABLE_INLINE_FUNCTION void MGUsup::DensityEnergyFromPressureTemperature(
-    const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
-  EOS_ERROR("MGUsup::DensityEnergyFromPressureTemperature: "
-            "Not implemented.\n");
-}
 // AEM: We should add entropy and Gruneissen parameters here so that it is complete
 // If we add also alpha and BT, those should also be in here.
 template <typename Indexer_t>

singularity-eos/eos/eos_powermg.hpp

@@ -157,11 +157,6 @@ class PowerMG : public EosBase<PowerMG> {
     }
     printf("\n\n");
   }
-  // Density/Energy from P/T not unique, if used will give error
-  template <typename Indexer_t>
-  PORTABLE_INLINE_FUNCTION void
-  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                       Indexer_t &&lambda, Real &rho, Real &sie) const;
   inline void Finalize() {}
   static std::string EosType() { return std::string("PowerMG"); }
   static std::string EosPyType() { return EosType(); }
@@ -431,13 +426,6 @@ PORTABLE_INLINE_FUNCTION Real PowerMG::EntropyFromDensityInternalEnergy(
   }
   return value;
 }
-// AEM: Give error since function is not well defined
-template <typename Indexer_t>
-PORTABLE_INLINE_FUNCTION void PowerMG::DensityEnergyFromPressureTemperature(
-    const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
-  EOS_ERROR("PowerMG::DensityEnergyFromPressureTemperature: "
-            "Not implemented.\n");
-}
 // AEM: We should add entropy and Gruneissen parameters here so that it is complete
 // If we add also alpha and BT, those should also be in here.
 template <typename Indexer_t>

singularity-eos/eos/eos_sap_polynomial.hpp

@@ -209,14 +209,6 @@ class SAP_Polynomial : public EosBase<SAP_Polynomial> {
     printf("      b3  = %g\n", _b3);
     _AZbar.PrintParams();
   }
-  template <typename Indexer_t>
-  PORTABLE_INLINE_FUNCTION void
-  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                       Indexer_t &&lambda, Real &rho, Real &sie) const {
-    PORTABLE_WARN("This function is a stub for an incomplete EoS.");
-    sie = 0.0;
-    rho = 0.0;
-  }
   inline void Finalize() {}
   static std::string EosType() { return std::string("SAP_Polynomial"); }
   static std::string EosPyType() { return EosType(); }

singularity-eos/eos/eos_stellar_collapse.hpp

@@ -653,7 +653,25 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::GruneisenParamFromDensityInternal
 template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperature(
     const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
-  EOS_ERROR("StellarCollapse::DensityEnergyFromPRessureTemperature is a stub");
+  using RootFinding1D::regula_falsi;
+  using RootFinding1D::Status;
+  checkLambda_(lambda);
+  Real lrguess = lRho_(rho);
+  Real lT = lT_(temp);
+  Real lP = P2lP_(press);
+  Real Ye = lambda[Lambda::Ye];
+
+  if ((lrguess < lRhoMin_) || (lrguess > lRhoMax_)) {
+    lrguess = lRho_(rhoNormal_);
+  }
+  auto lPofRT = [&](Real lR) { return lP_.interpToReal(Ye, lT, lR); };
+  auto status = regula_falsi(lPofRT, lP, lrguess, lRhoMin_, lRhoMax_, ROOT_THRESH,
+                             ROOT_THRESH, lrguess);
+
+  Real lE = lE_.interpToReal(Ye, lT, lrguess);
+  rho = rho_(lrguess);
+  sie = le2e_(lE);
+  lambda[Lambda::lT] = lT;
 }
 
 template <typename Indexer_t>

singularity-eos/eos/eos_vinet.hpp

@@ -131,6 +131,10 @@ class Vinet : public EosBase<Vinet> {
   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;
+
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; }
+
   // Generic functions provided by the base class. These contain e.g. the vector
   // overloads that use the scalar versions declared here
   SG_ADD_BASE_CLASS_USINGS(Vinet)
@@ -152,11 +156,6 @@ class Vinet : public EosBase<Vinet> {
     }
     printf("\n\n");
   }
-  // Density/Energy from P/T not unique, if used will give error
-  template <typename Indexer_t = Real *>
-  PORTABLE_INLINE_FUNCTION void
-  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
-                                       Indexer_t &&lambda, Real &rho, Real &sie) const;
   inline void Finalize() {}
   static std::string EosType() { return std::string("Vinet"); }
   static std::string EosPyType() { return EosType(); }
@@ -375,13 +374,6 @@ PORTABLE_INLINE_FUNCTION Real Vinet::BulkModulusFromDensityInternalEnergy(
   Vinet_F_DT_func(rho, temp, output);
   return output[7] * output[7] * rho;
 }
-// AEM: Give error since function is not well defined
-template <typename Indexer_t>
-PORTABLE_INLINE_FUNCTION void Vinet::DensityEnergyFromPressureTemperature(
-    const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const {
-  EOS_ERROR("Vinet::DensityEnergyFromPressureTemperature: "
-            "Not implemented.\n");
-}
 // AEM: We should add entropy and Gruneissen parameters here so that it is complete
 // If we add also alpha and BT, those should also be in here.
 template <typename Indexer_t>

singularity-eos/eos/modifiers/zsplit_eos.hpp

@@ -212,6 +212,16 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
     // dvdt, dpde unchanged
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  DensityEnergyFromPressureTemperature(const Real press, const Real temp,
+                                       Indexer_t &&lambda, Real &rho, Real &sie) const {
+    const Real scale = GetScale_(lambda);
+    const Real iscale = GetInvScale_(lambda);
+    t_.DensityEnergyFromPressureTemperature(iscale * press, temp, lambda, rho, sie);
+    sie *= scale;
+  }
+
   PORTABLE_INLINE_FUNCTION
   int nlambda() const noexcept { return NL + t_.nlambda(); }
   static constexpr unsigned long PreferredInput() { return T::PreferredInput(); }

test/eos_unit_test_helpers.hpp

@@ -34,6 +34,7 @@
 #include <memory>
 
 #include <ports-of-call/portability.hpp>
+#include <singularity-eos/base/robust_utils.hpp>
 
 inline std::string demangle(const char *name) {
 
@@ -143,6 +144,34 @@ inline void compare_two_eoss(const E1 &test_e, const E2 &ref_e) {
   return;
 }
 
+template <typename EOS, typename Indexer_t = Real *>
+PORTABLE_INLINE_FUNCTION bool
+CheckRhoSieFromPT(EOS eos, Real rho, Real T,
+                  Indexer_t &&lambda = static_cast<Real *>(nullptr)) {
+  const Real P = eos.PressureFromDensityTemperature(rho, T, lambda);
+  const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T, lambda);
+  Real rtest, etest;
+  eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest);
+  Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda);
+  bool results_good = (isClose(rho, rtest, bmod * 1e-8) || isClose(rho, rtest, 1e-8)) &&
+                      (isClose(sie, etest, bmod * 1e-8) || isClose(sie, etest, 1e-8));
+  if (!results_good) {
+    Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda);
+    Real residual = P_test - P;
+    printf("RhoSie of PT failure!\n"
+           "\trho_true = %.14e\n"
+           "\tsie_true = %.14e\n"
+           "\tP        = %.14e\n"
+           "\tT        = %.14e\n"
+           "\trho      = %.14e\n"
+           "\tsie      = %.14e\n"
+           "\tP_test   = %.14e\n"
+           "\tresidual = %.14e\n",
+           rho, sie, P, T, rtest, etest, P_test, residual);
+  }
+  return results_good;
+}
+
 // Macro that checks for an exception or is a no-op depending on
 // whether or not a non-serial backend is supplied
 #ifdef PORTABILITY_STRATEGY_NONE

test/test_eos_carnahan_starling.cpp

@@ -161,6 +161,18 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") {
                         "Energy");
         }
       }
+
+      WHEN("We check RhoSie from PT") {
+        int nwrong = 0;
+        portableReduce(
+            "Check RhoSieFromPT", 0, 1,
+            PORTABLE_LAMBDA(const int i, int &nw) {
+              nw +=
+                  !CheckRhoSieFromPT(eos, 7.8290736890381501e-03, 1.5320999999999999e+03);
+            },
+            nwrong);
+        REQUIRE(nwrong == 0);
+      }
     }
   }
 }