properly thread lambda through sie from

CHANGELOG.md

@@ -6,6 +6,7 @@
 - [[PR658]](https://github.com/lanl/singularity-eos/pull/XXX) Add `MinimumInternalEnergy`/`MaximumInternalEnergy` to the EOS introspection API, so energy bounds are reachable through modifiers and the `singularity::EOS` variant
 
 ### Fixed (Repair bugs, etc)
+- [[PR661]](https://github.com/lanl/singularity-eos/pull/661) Fix bug in InternalEnergyFromDensityPressure where the lambda wasn't threaded all the way through.
 
 ### Changed (changing behavior/API/variables/...)
 - [[PR658]](https://github.com/lanl/singularity-eos/pull/XXX) `ScaledEOS::CheckParams` now requires a strictly positive scale factor, where it previously accepted any nonzero value.

singularity-eos/eos/eos_base.hpp

@@ -719,11 +719,11 @@ class EosBase {
       return eos.PressureFromDensityInternalEnergy(rho, sie, lambda);
     };
     const Real sie_min =
-        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature());
+        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature(), lambda);
     // temp not bounded. just pick something huge.
-    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20);
+    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20, lambda);
     Real sie_guess =
-        (((sie_min < sie) && (sie < sie_max)) ? 0.5 * (sie_min + sie_max) : sie);
+        (((sie_min < sie) && (sie < sie_max)) ? sie : 0.5 * (sie_min + sie_max));
     auto status = regula_falsi(f, P, sie_guess, sie_min, sie_max, robust::EPS(),
                                robust::EPS(), sie);
     if (status == Status::FAIL) {
@@ -737,9 +737,9 @@ class EosBase {
       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());
+        eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature(), lambda);
     // temp not bounded. just pick something huge.
-    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20);
+    const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20, lambda);
     Real sie = 0.5 * (sie_min + sie_max);
     eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
     return sie;

singularity-eos/eos/eos_helmholtz.hpp

@@ -6,7 +6,7 @@
 // Original work is open-sourced under the CC-By license
 // https://creativecommons.org/licenses/by/4.0/
 //------------------------------------------------------------------------------
-// © 2023-2025. Triad National Security, LLC. All rights reserved.  This
+// © 2023-2026. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract 89233218CNA000001
 // for Los Alamos National Laboratory (LANL), which is operated by Triad
 // National Security, LLC for the U.S.  Department of Energy/National
@@ -19,6 +19,8 @@
 // publicly and display publicly, and to permit others to do so.
 //------------------------------------------------------------------------------
 
+// This file was partially modified by AI.
+
 #ifndef _SINGULARITY_EOS_EOS_HELMHOLTZ_HPP_
 #define _SINGULARITY_EOS_EOS_HELMHOLTZ_HPP_
 
@@ -458,6 +460,8 @@ class Helmholtz : public EosBase<Helmholtz> {
         options_(rad, gas, coul, ion, ele, verbose, newton_raphson) {}
 
   PORTABLE_INLINE_FUNCTION void CheckParams() const { electrons_.CheckParams(); }
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real MinimumTemperature() const { return electrons_.MinimumTemperature(); }
   constexpr static inline int nlambda() noexcept { return 3; }
   template <typename T>
   static inline constexpr bool NeedsLambda() {
@@ -955,7 +959,7 @@ PORTABLE_INLINE_FUNCTION Real Helmholtz::lTFromRhoSie_(const Real rho, const Rea
     }
     lT = electrons_.lTMax();
   }
-  IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   return lT;
 }
 
@@ -1004,6 +1008,7 @@ void HelmElectrons::GetFromDensityTemperature(Real rho, Real lT, Real Ye, Real Y
   rho = std::min(rhoMax(), std::max(rhoMin(), rho));
   De = std::min(rhoMax(), std::max(rhoMin(), De));
   lDe = std::min(lRhoMax(), std::max(lRhoMin(), lDe));
+  lT = std::min(lTMax(), std::max(lTMin(), lT));
   Real T = math_utils::pow10(lT);
 
   // Find central indexes in table

singularity-eos/eos/modifiers/zsplit_eos.hpp

@@ -178,7 +178,7 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
     const Real scale = GetScale_(lambda);
     const Real iscale = GetInvScale_(lambda);
     sie *= iscale;
-    t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda);
+    t_.InternalEnergyFromDensityPressure(rho, P * iscale, sie, lambda);
     sie *= scale;
   }
 

test/test_eos_helmholtz.cpp

@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2026. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract 89233218CNA000001
 // for Los Alamos National Laboratory (LANL), which is operated by Triad
 // National Security, LLC for the U.S.  Department of Energy/National
@@ -12,6 +12,8 @@
 // publicly and display publicly, and to permit others to do so.
 //------------------------------------------------------------------------------
 
+// This file was partially modified by AI.
+
 #ifdef SINGULARITY_TEST_HELMHOLTZ
 
 #include <cmath>
@@ -36,6 +38,51 @@
 
 using singularity::Helmholtz;
 const std::string filename = "../data/helmholtz/helm_table.dat";
+
+SCENARIO("Helmholtz internal energy from density and pressure",
+         "[HelmholtzEOS][SieFromRhoP]") {
+  Helmholtz host_eos(filename, true, true, false, true, true);
+  auto eos = host_eos.GetOnDevice();
+
+  int nwrong = 0;
+  portableReduce(
+      "Helmholtz internal energy from density and pressure", 0, 3,
+      PORTABLE_LAMBDA(const int k, int &nw) {
+        constexpr Real abar[3] = {1.0, 4.0, 12.0};
+        constexpr Real zbar[3] = {1.0, 2.0, 6.0};
+        constexpr Real densities[4] = {1e-3, 1e1, 1e5, 1e9};
+        constexpr Real temperatures[4] = {1e4, 1e6, 1e8, 1e10};
+        for (int i = 0; i < 4; ++i) {
+          for (int j = 0; j < 4; ++j) {
+            Real lambda[3] = {abar[k], zbar[k], -1.0};
+            const Real rho = densities[i];
+            const Real temp = temperatures[j];
+            const Real expected =
+                eos.InternalEnergyFromDensityTemperature(rho, temp, lambda);
+            const Real pressure = eos.PressureFromDensityTemperature(rho, temp, lambda);
+
+            Real sie = 0.0;
+            eos.InternalEnergyFromDensityPressure(rho, pressure, sie, lambda);
+            const Real sie_returned =
+                eos.InternalEnergyFromDensityPressure(rho, pressure, lambda);
+            if (!isClose(sie, expected, 1e-6) || !isClose(sie_returned, expected, 1e-6)) {
+              printf("Helmholtz sie mismatch: Abar=%g Zbar=%g rho=%g T=%g "
+                     "expected=%.14e output=%.14e returned=%.14e\n",
+                     abar[k], zbar[k], rho, temp, expected, sie, sie_returned);
+              nw += 1;
+            }
+          }
+        }
+      },
+      nwrong);
+  THEN("Both scalar overloads recover the energy with a composition lambda") {
+    REQUIRE(nwrong == 0);
+  }
+
+  eos.Finalize();
+  host_eos.Finalize();
+}
+
 SCENARIO("Helmholtz equation of state - Table interpolation (tgiven)", "[HelmholtzEOS]") {
   GIVEN("A Helmholtz EOS") {
     /* We only test the EOS without Coulomb corrections since those are

test/test_eos_stellar_collapse.cpp

@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2021-2024. Triad National Security, LLC. All rights reserved.  This
+// © 2021-2026. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract 89233218CNA000001
 // for Los Alamos National Laboratory (LANL), which is operated by Triad
 // National Security, LLC for the U.S.  Department of Energy/National
@@ -12,6 +12,8 @@
 // publicly and display publicly, and to permit others to do so.
 //------------------------------------------------------------------------------
 
+// This file was partially modified by AI.
+
 #include <array>
 #include <cmath>
 #include <cstdio>
@@ -40,6 +42,43 @@
 
 #include <singularity-eos/base/spiner_table_utils.hpp>
 
+SCENARIO("StellarCollapse internal energy from density and pressure",
+         "[StellarCollapse][SieFromRhoP]") {
+  singularity::StellarCollapse sc("./goldfiles/stellar_collapse_ideal.h5", false, false);
+  auto eos = sc.GetOnDevice();
+  const Real yemin = sc.YeMin();
+  const Real yemax = sc.YeMax();
+  const Real lrhomin = std::log10(sc.rhoMin());
+  const Real lrhomax = std::log10(sc.rhoMax());
+  const Real ltmin = std::log10(sc.TMin());
+  const Real ltmax = std::log10(sc.TMax());
+
+  int nwrong = 0;
+  portableReduce(
+      "StellarCollapse internal energy from density and pressure", 0, 3,
+      PORTABLE_LAMBDA(const int i, int &nw) {
+        const Real fraction = (i + 1) / 4.0;
+        const Real rho = std::pow(10., lrhomin + fraction * (lrhomax - lrhomin));
+        const Real temp = std::pow(10., ltmin + fraction * (ltmax - ltmin));
+        Real lambda[2] = {yemin + fraction * (yemax - yemin), 0.0};
+        const Real expected = eos.InternalEnergyFromDensityTemperature(rho, temp, lambda);
+        const Real pressure = eos.PressureFromDensityTemperature(rho, temp, lambda);
+
+        Real sie = 0.0;
+        eos.InternalEnergyFromDensityPressure(rho, pressure, sie, lambda);
+        nw += !isClose(sie, expected, 1e-6);
+        nw += !isClose(eos.InternalEnergyFromDensityPressure(rho, pressure, lambda),
+                       expected, 1e-6);
+      },
+      nwrong);
+  THEN("Both scalar overloads recover the energy with an electron-fraction lambda") {
+    REQUIRE(nwrong == 0);
+  }
+
+  eos.Finalize();
+  sc.Finalize();
+}
+
 template <typename EOS_t>
 void CompareStellarCollapse(EOS_t sc, EOS_t sc2) {
   Real yemin = sc.YeMin();

test/test_eos_zsplit.cpp

@@ -1,5 +1,5 @@
 //------------------------------------------------------------------------------
-// © 2024-2025. Triad National Security, LLC. All rights reserved.  This
+// © 2024-2026. Triad National Security, LLC. All rights reserved.  This
 // program was produced under U.S. Government contract 89233218CNA000001
 // for Los Alamos National Laboratory (LANL), which is operated by Triad
 // National Security, LLC for the U.S.  Department of Energy/National
@@ -12,6 +12,8 @@
 // publicly and display publicly, and to permit others to do so.
 //------------------------------------------------------------------------------
 
+// This file was partially modified by AI.
+
 #include <cmath>
 #include <cstdio>
 #include <cstdlib>
@@ -42,6 +44,50 @@ using Lambda_t = singularity::IndexerUtils::VariadicIndexer<MeanIonizationState>
 using EOS =
     singularity::Variant<IdealGas, IdealElectrons, ZSplitI<IdealGas>, ZSplitE<IdealGas>>;
 
+SCENARIO("ZSplit internal energy from density and pressure", "[ZSplit][SieFromRhoP]") {
+  constexpr Real gm1 = 2. / 3.;
+  constexpr Real Cv = 1e8;
+  auto eos_i = ZSplitI<IdealGas>(IdealGas(gm1, Cv));
+  auto eos_e = ZSplitE<IdealGas>(IdealGas(gm1, Cv));
+  auto ions = eos_i.GetOnDevice();
+  auto electrons = eos_e.GetOnDevice();
+
+  int nwrong = 0;
+  portableReduce(
+      "ZSplit internal energy from density and pressure", 0, 3,
+      PORTABLE_LAMBDA(const int i, int &nw) {
+        const Real rho = 1.0 + i;
+        const Real temp = 1e3 * (i + 1);
+        const Real Z = 0.5 + i;
+        Lambda_t lambda;
+        lambda[MeanIonizationState()] = Z;
+        const Real ei = Cv * temp / (Z + 1);
+        const Real ee = Z * ei;
+        const Real Pi = gm1 * rho * ei;
+        const Real Pe = gm1 * rho * ee;
+
+        Real sie_i = 0.0;
+        Real sie_e = 0.0;
+        ions.InternalEnergyFromDensityPressure(rho, Pi, sie_i, lambda);
+        electrons.InternalEnergyFromDensityPressure(rho, Pe, sie_e, lambda);
+        nw += !isClose(sie_i, ei, 1e-12);
+        nw += !isClose(sie_e, ee, 1e-12);
+        nw +=
+            !isClose(ions.InternalEnergyFromDensityPressure(rho, Pi, lambda), ei, 1e-12);
+        nw += !isClose(electrons.InternalEnergyFromDensityPressure(rho, Pe, lambda), ee,
+                       1e-12);
+      },
+      nwrong);
+  THEN("Both scalar overloads recover the ion and electron energies") {
+    REQUIRE(nwrong == 0);
+  }
+
+  ions.Finalize();
+  electrons.Finalize();
+  eos_i.Finalize();
+  eos_e.Finalize();
+}
+
 SCENARIO("ZSplit of Ideal Gas", "[ZSplit][IdealGas][IdealElectrons]") {
   GIVEN("An ideal gas EOS") {
     constexpr Real gm1 = (5. / 3.) - 1.;