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. diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 17c63e224c0..70832d732c8 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/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(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; 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 diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index 7114f9bf2c4..9781c1bde81 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -178,7 +178,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_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 diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index a0fc54123e6..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 @@ -40,6 +42,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..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 @@ -42,6 +44,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.;