From 7f028e78cbf71f3d57f625684ec6ca9d6feaf7d1 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 25 Sep 2026 12:02:55 -0400 Subject: [PATCH 1/8] fix lambda call in sie-from-rho-p in base --- singularity-eos/eos/eos_base.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 60823b3f0e1..c55c72c5a1c 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -694,9 +694,9 @@ 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); auto status = regula_falsi(f, P, sie_guess, sie_min, sie_max, robust::EPS(), @@ -712,9 +712,9 @@ class EosBase { Indexer_t &&lambda = static_cast(nullptr)) const { const CRTP &eos = *(static_cast(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; From b6a65db830e70089842abd62b40543dd8a9c0c27 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 25 Sep 2026 13:37:14 -0400 Subject: [PATCH 2/8] fix bug and add unit test --- singularity-eos/eos/modifiers/zsplit_eos.hpp | 2 +- test/test_eos_stellar_collapse.cpp | 37 ++++++++++++++++ test/test_eos_zsplit.cpp | 44 ++++++++++++++++++++ 3 files changed, 82 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index 827f3d9145f..f0295fea235 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -176,7 +176,7 @@ class ZSplit : public EosBase> { 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; } diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index a0fc54123e6..141fc76a61d 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/test/test_eos_stellar_collapse.cpp @@ -40,6 +40,43 @@ #include +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 void CompareStellarCollapse(EOS_t sc, EOS_t sc2) { Real yemin = sc.YeMin(); diff --git a/test/test_eos_zsplit.cpp b/test/test_eos_zsplit.cpp index 81b50ae832a..f80f001f07c 100644 --- a/test/test_eos_zsplit.cpp +++ b/test/test_eos_zsplit.cpp @@ -42,6 +42,50 @@ using Lambda_t = singularity::IndexerUtils::VariadicIndexer using EOS = singularity::Variant, ZSplitE>; +SCENARIO("ZSplit internal energy from density and pressure", "[ZSplit][SieFromRhoP]") { + constexpr Real gm1 = 2. / 3.; + constexpr Real Cv = 1e8; + auto eos_i = ZSplitI(IdealGas(gm1, Cv)); + auto eos_e = ZSplitE(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.; From a6b581503e646cfa615f29cf154bcc91bb722500 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 25 Sep 2026 13:39:02 -0400 Subject: [PATCH 3/8] AI disclaimer and CC --- singularity-eos/eos/modifiers/zsplit_eos.hpp | 2 ++ test/test_eos_stellar_collapse.cpp | 4 +++- test/test_eos_zsplit.cpp | 4 +++- 3 files changed, 8 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index f0295fea235..e59210e9dfb 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -12,6 +12,8 @@ // publicly and display publicly, and to permit others to do so. //------------------------------------------------------------------------------ +// This file was partially modified by AI. + #ifndef _SINGULARITY_EOS_EOS_ZSPLIT_EOS_ #define _SINGULARITY_EOS_EOS_ZSPLIT_EOS_ diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index 141fc76a61d..a9e60fbb9be 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/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 #include #include diff --git a/test/test_eos_zsplit.cpp b/test/test_eos_zsplit.cpp index f80f001f07c..2f3eb55fa12 100644 --- a/test/test_eos_zsplit.cpp +++ b/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 #include #include From 23eca2cca3155a2d8fa2c93a9a85460b99e1d9a7 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Fri, 25 Sep 2026 13:42:31 -0400 Subject: [PATCH 4/8] CHANGELOG --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index f10e8526e04..68081afc280 100644 --- a/CHANGELOG.md +++ b/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. From 0fd5705241925b1d86867d75aea58e383a69c294 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Tue, 29 Sep 2026 11:17:10 -0400 Subject: [PATCH 5/8] fix eos base guess --- singularity-eos/eos/eos_base.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 414f975fe3e..a4b463749d6 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -723,7 +723,7 @@ class EosBase { // temp not bounded. just pick something huge. 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) { From 623a480409f174f00a03bb1516df3920f278df93 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 1 Oct 2026 16:08:33 -0400 Subject: [PATCH 6/8] fix helmholtz minor bugs --- singularity-eos/eos/eos_helmholtz.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index 69d472410f6..3f938692028 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/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 { 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 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(lambda, Lambda::lT, lT); + IndexerUtils::SafeSet(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 From 46bc2e47e459e3c824874d3a602ce92aa265d1da Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 1 Oct 2026 16:16:52 -0400 Subject: [PATCH 7/8] add a test for Helmholtz EOS for the sie(rho, P) call, which is the only other failure mode --- test/test_eos_helmholtz.cpp | 49 ++++++++++++++++++++++++++++++++++++- 1 file changed, 48 insertions(+), 1 deletion(-) diff --git a/test/test_eos_helmholtz.cpp b/test/test_eos_helmholtz.cpp index d3735c2ea24..77ebeb79ece 100644 --- a/test/test_eos_helmholtz.cpp +++ b/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 @@ -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 From d02b01d964f4f4d9a7f36afc8fcb5a10564f34c7 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 1 Oct 2026 16:49:32 -0400 Subject: [PATCH 8/8] formatting --- singularity-eos/eos/eos_base.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index a4b463749d6..70832d732c8 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -723,7 +723,7 @@ class EosBase { // temp not bounded. just pick something huge. const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20, lambda); Real sie_guess = - (((sie_min < sie) && (sie < sie_max)) ? sie : 0.5 * (sie_min + sie_max)); + (((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) {