lanl · Yurlungur · Dec 19, 2025 · Nov 26, 2025 · Nov 26, 2025 · Nov 26, 2025
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -3,6 +3,7 @@
 ## Current develop
 
 ### Added (new features/APIs/variables/...)
+- [[PR589]](https://github.com/lanl/singularity-eos/pull/589) InternalEnergyFromDensityPressure
 
 ### Fixed (Repair bugs, etc)
 

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
@@ -1181,6 +1181,25 @@ given density in :math:`g/cm^3` and temperature in Kelvin.
 returns the unitless Gruneisen parameter given density in
 :math:`g/cm^3` and specific internal energy in :math:`erg/g`.
 
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   Real InternalEnergyFromDensityPressure(const Real rho, const Real P, 
+                                          Indexer_t &&lambda = nullptr) const;
+
+returns the specific internal energy in :math:`erg/g` given a density
+in :math:`g/cm^3` and a pressure in Barye. This call is often a root
+find, and thus an alternative signature allows an initial guess to be
+passed in by reference:
+
+.. code-block:: cpp
+
+   template <typename Indexer_t = Real*>
+   void InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                          Indexer_t &&lambda = nullptr) const;
+
+In this case, the guess will be set to the new value by the function.
+
 .. code-block:: cpp
 
    template <typename Indexer_t = Real*>

diff --git a/python/module.hpp b/python/module.hpp
@@ -169,6 +169,7 @@ EOS_VEC_FUNC_TMPL(BulkModulusFromDensityTemperature, rhos, temperatures, bmods)
 EOS_VEC_FUNC_TMPL(BulkModulusFromDensityInternalEnergy, rhos, sies, bmods)
 EOS_VEC_FUNC_TMPL(GruneisenParamFromDensityTemperature, rhos, temperatures, gm1s)
 EOS_VEC_FUNC_TMPL(GruneisenParamFromDensityInternalEnergy, rhos, sies, gm1s)
+EOS_VEC_FUNC_TMPL(InternalEnergyFromDensityPressure, rhos, Ps, sies)
 
 struct EOSState {
   Real density;
@@ -229,6 +230,7 @@ struct VectorFunctions {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergy<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperature<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergy<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("lmbdas"))
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressure<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("lmbdas"))
 
     .def("TemperatureFromDensityInternalEnergy", &TemperatureFromDensityInternalEnergyNoLambda<T>, py::arg("rhos"), py::arg("sies"), py::arg("temperatures"))
     .def("InternalEnergyFromDensityTemperature", &InternalEnergyFromDensityTemperatureNoLambda<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("sies"))
@@ -240,7 +242,7 @@ struct VectorFunctions {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyNoLambda<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureNoLambda<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyNoLambda<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"))
-
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureNoLambda<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"))
 
     .def("FillEos", [](const T & self, py::array_t<Real> rhos,
                        py::array_t<Real> temperatures, py::array_t<Real> sies,
@@ -300,6 +302,7 @@ struct VectorFunctions<T,true> {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("scratch"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureWithScratch<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("scratch"), py::arg("lmbdas"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("scratch"), py::arg("lmbdas"))
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureWithScratch<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("scratch"), py::arg("lmbdas"))
 
     .def("TemperatureFromDensityInternalEnergy", &TemperatureFromDensityInternalEnergyNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("temperatures"), py::arg("scratch"))
     .def("InternalEnergyFromDensityTemperature", &InternalEnergyFromDensityTemperatureNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("sies"), py::arg("scratch"))
@@ -311,6 +314,7 @@ struct VectorFunctions<T,true> {
     .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("scratch"))
     .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("scratch"))
     .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("scratch"))
+    .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureNoLambdaWithScratch<T>, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("scratch"))
 
     .def("FillEos", [](const T & self, py::array_t<Real> rhos,
                        py::array_t<Real> temperatures, py::array_t<Real> sies,
@@ -452,6 +456,14 @@ py::class_<T> eos_class(py::module_ & m, std::string name) {
       self.DensityEnergyFromPressureTemperature(press, temp, np<Real>(), rho, sie);
       return std::pair<Real, Real>(rho, sie);
     }, py::arg("press"), py::arg("temp"))
+    .def("InternalEnergyFromDensityPressure", [](const T & self, const Real rho, const Real P, Real sie) {
+      self.InternalEnergyFromDensityPressure(rho, P, sie, np<Real>());
+      return sie;
+    }, py::arg("rho"), py::arg("P"), py::arg("sie"))
+    .def("InternalEnergyFromDensityPressure", [](const T & self, const Real rho, const Real P, Real sie, py::array_t<Real> lambda) {
+      self.InternalEnergyFromDensityPressure(rho, P, sie, lambda.mutable_data());
+      return sie;
+    }, py::arg("rho"), py::arg("P"), py::arg("sie"), py::arg("lmbda"))
     .def("Finalize", &T::Finalize)
     .def_property_readonly_static("EosType", [](py::object) { return T::EosType(); })
     .def_static("scratch_size", &T::scratch_size, py::arg("method_name"), py::arg("nelements"))

diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp
@@ -82,6 +82,7 @@ char *StrCat(char *destination, const char *source) {
   using EosBase<__VA_ARGS__>::BulkModulusFromDensityInternalEnergy;                      \
   using EosBase<__VA_ARGS__>::GruneisenParamFromDensityTemperature;                      \
   using EosBase<__VA_ARGS__>::GruneisenParamFromDensityInternalEnergy;                   \
+  using EosBase<__VA_ARGS__>::InternalEnergyFromDensityPressure;                         \
   using EosBase<__VA_ARGS__>::FillEos;                                                   \
   using EosBase<__VA_ARGS__>::EntropyFromDensityTemperature;                             \
   using EosBase<__VA_ARGS__>::EntropyFromDensityInternalEnergy;                          \
@@ -928,6 +929,76 @@ class EosBase {
                                               rho, sie);
   }
 
+  // JMM: Another set of calls that are often overloaded for special cases
+  // TODO(JMM): Do we also want TemperatureFromDensityPressure? That's
+  // the more likely fundamental call, but the less likely "useful"
+  // call...
+  // We pass in sie by value to allow for initial guesses
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure(
+      const Real rho, const Real P, Real &sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    using RootFinding1D::regula_falsi;
+    using RootFinding1D::Status;
+    const CRTP &eos = *(static_cast<CRTP const *>(this));
+    auto f = [&](const Real sie) {
+      return eos.PressureFromDensityInternalEnergy(rho, sie, lambda);
+    };
+    const Real sie_min =
+        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature());
+    // temp not bounded. just pick something huge.
+    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20);
+    Real sie_guess =
+        (((sie_min < sie) && (sie < sie_max)) ? 0.5 * (sie_min + sie_max) : sie);
+    auto status = regula_falsi(f, P, sie_guess, sie_min, sie_max, robust::EPS(),
+                               robust::EPS(), sie);
+    if (status == Status::FAIL) {
+      PORTABLE_WARN("InternalEnergyFromDensityPressure: failed to find root\n");
+    }
+  }
+  // This version foregoes the initial guess. You get what you get.
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure(
+      const Real rho, const Real P,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    const CRTP &eos = *(static_cast<CRTP const *>(this));
+    const Real sie_min =
+        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature());
+    // temp not bounded. just pick something huge.
+    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20);
+    Real sie = 0.5 * (sie_min + sie_max);
+    eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    return sie;
+  }
+  // Classes like EOSPAC probably wish to overload/shadow this version
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                const int num,
+                                                LambdaIndexer &&lambdas) const {
+    static auto const name = SG_MEMBER_FUNC_NAME();
+    static auto const cname = name.c_str();
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    portableFor(
+        cname, 0, num, PORTABLE_LAMBDA(const int i) {
+          copy.InternalEnergyFromDensityPressure(rhos[i], Ps[i], sies[i], lambdas[i]);
+        });
+  }
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                Real * /*scratch */, const int num,
+                                                LambdaIndexer &&lambdas) const {
+    return InternalEnergyFromDensityPressure(rhos, Ps, sies, num, lambdas);
+  }
+  template <typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(const Real *rhos, const Real *Ps,
+                                                Real *sies, Real * /*scratch */,
+                                                const int num, LambdaIndexer &&lambdas,
+                                                Transform && = Transform()) const {
+    return InternalEnergyFromDensityPressure(rhos, Ps, sies, num, lambdas);
+  }
+
   // Serialization
   /*
     The methodology here is there are *three* size methods all EOS's provide:

diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp
@@ -494,6 +494,8 @@ class EOSPAC : public EosBase<EOSPAC> {
 
   // TODO(JMM): Add performant Gibbs Free Energy
   using EosBase<EOSPAC>::FillEos;
+  // TODO(JMM): Add a performant internal energy call
+  using EosBase<EOSPAC>::InternalEnergyFromDensityPressure;
 
   SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_)
 

diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp
@@ -186,6 +186,14 @@ class IdealGas : public EosBase<IdealGas> {
     sie = MYMAX(0.0, _Cv * temp);
     rho = MYMAX(0.0, press / (_gm1 * sie));
   }
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure(
+      const Real rho, const Real P, Real &sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    // P = gm1 rho sie
+    sie = P / (_gm1 * rho);
+  }
+
   inline void Finalize() {}
   static std::string EosType() { return std::string("IdealGas"); }
   static std::string EosPyType() { return EosType(); }

diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp
@@ -336,6 +336,17 @@ class Variant {
         eos_);
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure(
+      const Real rho, const Real P, Real &sie,
+      Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
+    return PortsOfCall::visit(
+        [&rho, &P, &sie, &lambda](const auto &eos) {
+          return eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+        },
+        eos_);
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real RhoPmin(const Real temp) const {
     return PortsOfCall::visit([&temp](const auto &eos) { return eos.RhoPmin(temp); },
@@ -1194,6 +1205,55 @@ class Variant {
         eos_);
   }
 
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return InternalEnergyFromDensityPressure(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+        std::forward<RealIndexer>(sies), num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                const int num,
+                                                LambdaIndexer &&lambdas) const {
+    return PortsOfCall::visit(
+        [&](const auto &eos) {
+          return eos.InternalEnergyFromDensityPressure(
+              std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+              std::forward<RealIndexer>(sies), num, std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                Real *scratch, const int num) const {
+    NullIndexer lambdas{}; // Returns null pointer for every index
+    return InternalEnergyFromDensityPressure(
+        std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+        std::forward<RealIndexer>(sies), scratch, num, lambdas);
+  }
+
+  template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
+  inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos,
+                                                ConstRealIndexer &&Ps, RealIndexer &&sies,
+                                                Real *scratch, const int num,
+                                                LambdaIndexer &&lambdas) const {
+    return PortsOfCall::visit(
+        [&rhos, &Ps, &sies, &scratch, &num, &lambdas](const auto &eos) {
+          return eos.InternalEnergyFromDensityPressure(
+              std::forward<ConstRealIndexer>(rhos), std::forward<ConstRealIndexer>(Ps),
+              std::forward<RealIndexer>(sies), scratch, num,
+              std::forward<LambdaIndexer>(lambdas));
+        },
+        eos_);
+  }
+
   template <typename RealIndexer>
   inline void FillEos(RealIndexer &&rhos, RealIndexer &&temps, RealIndexer &&energies,
                       RealIndexer &&presses, RealIndexer &&cvs, RealIndexer &&bmods,

diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp
@@ -193,6 +193,14 @@ class UnitSystem : public EosBase<UnitSystem<T>> {
                                                              temp * temp_unit_, lambda);
     return gm1;
   }
+  template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    sie *= sie_unit_;
+    t_.InternalEnergyFromDensityPressure(rho_unit_ * rho, press_unit_ * P, sie, lambda);
+    sie *= inv_sie_unit_;
+  }
 
   template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void

diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp
@@ -121,6 +121,12 @@ class FlooredEnergy : public EosBase<FlooredEnergy<T>> {
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+  }
+  template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
                                  Indexer_t &&lambda = nullptr) const {

diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp
@@ -140,6 +140,12 @@ class RelativisticEOS : public EosBase<RelativisticEOS<T>> {
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+  }
+  template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
                                  Indexer_t &&lambda = nullptr) const {

diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp
@@ -136,6 +136,12 @@ class ScaledEOS : public EosBase<ScaledEOS<T>> {
     return t_.GruneisenParamFromDensityTemperature(scale_ * rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    t_.InternalEnergyFromDensityPressure(scale_ * rho, P, sie, lambda);
+  }
+  template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
                                  Indexer_t &&lambda = nullptr) const {

diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp
@@ -127,6 +127,14 @@ class ShiftedEOS : public EosBase<ShiftedEOS<T>> {
     return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda);
   }
   template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    sie -= shift_;
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    sie += shift_;
+  }
+  template <typename Indexer_t = Real *>
   PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press,
                                  Real &cv, Real &bmod, const unsigned long output,
                                  Indexer_t &&lambda = nullptr) const {

diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -169,6 +169,17 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
     return t_.GruneisenParamFromDensityInternalEnergy(rho, iscale * sie, lambda);
   }
 
+  template <typename Indexer_t = Real *>
+  PORTABLE_FUNCTION void
+  InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie,
+                                    Indexer_t &&lambda = nullptr) const {
+    const Real scale = GetScale_(lambda);
+    const Real iscale = GetInvScale_(lambda);
+    sie *= iscale;
+    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    sie *= scale;
+  }
+
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION void
   FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod,