From 9301bb663b5dd154cd9bd18af1f070d533b5afb9 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 16 Jul 2026 12:53:48 -0400 Subject: [PATCH 1/4] intermediate commit with JWL --- test/CMakeLists.txt | 1 + test/test_eos_jwl.cpp | 163 ++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 164 insertions(+) create mode 100644 test/test_eos_jwl.cpp diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index 4289a9c0893..9b40ff1972f 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable( eos_unit_test_helpers.hpp test_eos_ideal.cpp test_eos_gruneisen.cpp + test_eos_jwl.cpp test_eos_davis.cpp test_eos_sap_polynomial.cpp test_eos_noble_abel.cpp diff --git a/test/test_eos_jwl.cpp b/test/test_eos_jwl.cpp new file mode 100644 index 00000000000..c44c97b41e3 --- /dev/null +++ b/test/test_eos_jwl.cpp @@ -0,0 +1,163 @@ +//------------------------------------------------------------------------------ +// © 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 +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +// This file was made in part with generative AI. + +#include + +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include + +using singularity::JWL; +using EOS = singularity::Variant; + +namespace { + +constexpr Real A = 6.3207e14; +constexpr Real B = -4.472e10; +constexpr Real R1 = 11.3; +constexpr Real R2 = 1.13; +constexpr Real w = 0.8938; +constexpr Real rho0 = 1.895; +constexpr Real Cv = 1.31240105540897101e7; + +EOS MakeJWL() { return JWL(A, B, R1, R2, w, rho0, Cv); } + +} // namespace + +SCENARIO("JWL EOS basic calls", "[JWL]") { + GIVEN("A JWL EOS for detonation products") { + EOS eos_host = MakeJWL(); + auto eos = eos_host.GetOnDevice(); + + WHEN("Thermodynamic properties are evaluated at representative states") { + constexpr int num = 4; + int nwrong = 0; + portableReduce( + "Check JWL thermodynamic properties", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + Real rho, temp, sie_true, press_true, bmod_true; + if (i == 0) { + rho = 1.42125000000000012e0; + temp = 5.00000000000000000e2; + sie_true = 1.94155962893205547e9; + press_true = -1.39539537633559513e9; + bmod_true = 3.57742905528225708e9; + } else if (i == 1) { + rho = 1.89500000000000002e0; + temp = 1.00000000000000000e3; + sie_true = 6.74299306354014778e9; + press_true = 1.56033130109453297e10; + bmod_true = 1.14145142327138275e11; + } else if (i == 2) { + rho = 2.84250000000000025e0; + temp = 2.00000000000000000e3; + sie_true = 3.22063102789073410e10; + press_true = 3.83759661376449402e11; + bmod_true = 2.65765485356934912e12; + } else { + rho = 3.79000000000000004e0; + temp = 3.00000000000000000e3; + sie_true = 1.31330196851926880e11; + press_true = 2.33127259815571191e12; + bmod_true = 1.27999609124014121e13; + } + + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, temp); + nw += !isClose(sie, sie_true, 1e-12); + nw += + !isClose(eos.TemperatureFromDensityInternalEnergy(rho, sie), temp, 1e-12); + nw += !isClose(eos.PressureFromDensityTemperature(rho, temp), press_true, + 1e-12); + nw += !isClose(eos.PressureFromDensityInternalEnergy(rho, sie), press_true, + 1e-12); + nw += !isClose(eos.SpecificHeatFromDensityTemperature(rho, temp), Cv, 1e-12); + nw += + !isClose(eos.SpecificHeatFromDensityInternalEnergy(rho, sie), Cv, 1e-12); + nw += !isClose(eos.BulkModulusFromDensityTemperature(rho, temp), bmod_true, + 1e-12); + nw += !isClose(eos.BulkModulusFromDensityInternalEnergy(rho, sie), bmod_true, + 1e-12); + nw += !isClose(eos.GruneisenParamFromDensityTemperature(rho, temp), w, 1e-12); + nw += + !isClose(eos.GruneisenParamFromDensityInternalEnergy(rho, sie), w, 1e-12); + }, + nwrong); + + THEN("The values agree with the reference results") { REQUIRE(nwrong == 0); } + } + + WHEN("Density and energy are recovered from pressure and temperature") { + constexpr int num = 16; + int nwrong = 0; + portableReduce( + "Check JWL pressure-temperature inversion", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + const Real rho = rho0 * (1.0 + static_cast(i) / num); + const Real temp = 500.0 + 100.0 * i; + nw += !CheckRhoSieFromPT(eos, rho, temp); + }, + nwrong); + + THEN("The original density and energy are recovered") { REQUIRE(nwrong == 0); } + } + + eos.Finalize(); + eos_host.Finalize(); + } +} + +SCENARIO("JWL isentropic bulk modulus agrees with finite differences", + "[JWL][BulkModulus]") { + GIVEN("A JWL EOS and a range of physical states") { + EOS eos_host = MakeJWL(); + auto eos = eos_host.GetOnDevice(); + + constexpr int num = 32; + int nwrong = 0; + portableReduce( + "Compare JWL bulk modulus to finite differences", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + const Real rho = rho0 * (0.8 + 0.04 * i); + const Real temp = 500.0 + 75.0 * i; + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, temp); + const Real press = eos.PressureFromDensityInternalEnergy(rho, sie); + const Real drho = 1e-6 * rho; + const Real dsie = 1e-6 * std::abs(sie); + + const Real dpdr_e = (eos.PressureFromDensityInternalEnergy(rho + drho, sie) - + eos.PressureFromDensityInternalEnergy(rho - drho, sie)) / + (2.0 * drho); + const Real dpde_r = (eos.PressureFromDensityInternalEnergy(rho, sie + dsie) - + eos.PressureFromDensityInternalEnergy(rho, sie - dsie)) / + (2.0 * dsie); + const Real bmod_fd = rho * dpdr_e + press / rho * dpde_r; + const Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, sie); + nw += !isClose(bmod, bmod_fd, 1e-8); + }, + nwrong); + + THEN("The analytic and finite-difference results agree") { REQUIRE(nwrong == 0); } + + eos.Finalize(); + eos_host.Finalize(); + } +} From 7cd896f5faa2b438c4b7bed3dc7436d4de9d2f42 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 16 Jul 2026 15:47:38 -0400 Subject: [PATCH 2/4] JWL test and fix for issue with root finder --- singularity-eos/eos/eos_jwl.hpp | 18 +++++++++++++++ test/test_eos_jwl.cpp | 39 +++++++++++++++++++++++++++++++++ 2 files changed, 57 insertions(+) diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index cb5106fc26c..8bb333f3a55 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -73,6 +73,18 @@ class JWL : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const; template + PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure( + const Real rho, const Real press, Real &sie, + Indexer_t &&lambda = static_cast(nullptr)) const; + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( + const Real rho, const Real press, + Indexer_t &&lambda = static_cast(nullptr)) const { + Real sie; + INternalEnergyFromDensityPressure(rho, press, sie, lambda); + return sie; + } + template PORTABLE_INLINE_FUNCTION Real EntropyFromDensityTemperature(const Real rho, const Real temperature, Indexer_t &&lambda = static_cast(nullptr)) const; @@ -163,6 +175,12 @@ PORTABLE_INLINE_FUNCTION Real JWL::PressureFromDensityInternalEnergy( return ReferencePressure(rho) + _w * rho * (sie - ReferenceEnergy(rho)); } template +PORTABLE_INLINE_FUNCTION void +JWL::InternalEnergyFromDensityPressure(const Real rho, const Real press, Real &sie, + Indexer_t &&lambda) const { + sie = ReferenceEnergy(rho) + robust::ratio(press - ReferencePressure(rho), _w * rho); +} +template PORTABLE_INLINE_FUNCTION Real JWL::EntropyFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda) const { EntropyIsNotEnabled("JWL"); diff --git a/test/test_eos_jwl.cpp b/test/test_eos_jwl.cpp index c44c97b41e3..9018d3ab725 100644 --- a/test/test_eos_jwl.cpp +++ b/test/test_eos_jwl.cpp @@ -161,3 +161,42 @@ SCENARIO("JWL isentropic bulk modulus agrees with finite differences", eos_host.Finalize(); } } + +SCENARIO("JWL internal energy from density and pressure", "[JWL][SieFromRhoP]") { + GIVEN("A state from Lee et al. (2013). 10.1016/j.jcp.2013.03.046") { + // The source uses g, cm, and microseconds, so pressures in Mbar and specific + // energies in (cm/us)^2 both acquire a factor of 1e12 when converted to cgs. + constexpr Real paper_units_to_cgs = 1e12; + constexpr Real A = 632.1 * paper_units_to_cgs; + constexpr Real B = -0.04472 * paper_units_to_cgs; + constexpr Real R1 = 11.3; + constexpr Real R2 = 1.13; + constexpr Real w = 0.8938; + constexpr Real rho0 = 1.905; + constexpr Real Cv = 0.00003888 * paper_units_to_cgs; + constexpr Real rho = 3.810; + constexpr Real press = 2.0 * paper_units_to_cgs; + constexpr Real sie_true = 0.0333355 * paper_units_to_cgs; + + EOS eos_host = JWL(A, B, R1, R2, w, rho0, Cv); + auto eos = eos_host.GetOnDevice(); + + WHEN("Internal energy is computed from density and pressure") { + Real sie = 0.0; + portableReduce( + "Compute JWL internal energy from density and pressure", 0, 1, + PORTABLE_LAMBDA(const int, Real &result) { + eos.InternalEnergyFromDensityPressure(rho, press, result); + }, + sie); + + THEN("It agrees with the published state") { + INFO("Computed sie: " << sie << ", expected sie: " << sie_true); + REQUIRE(isClose(sie, sie_true, 1e-6)); + } + } + + eos.Finalize(); + eos_host.Finalize(); + } +} From a15c28e722a3813a4f1ea2ef6e4a391ea5799397 Mon Sep 17 00:00:00 2001 From: Jonah Maxwell Miller Date: Thu, 16 Jul 2026 15:54:49 -0400 Subject: [PATCH 3/4] changelog, CC, AI disclaimer --- CHANGELOG.md | 1 + singularity-eos/eos/eos_jwl.hpp | 4 +++- 2 files changed, 4 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 71f298d1e89..7e332a331cc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -11,6 +11,7 @@ - [[PR638]](https://github.com/lanl/singularity-eos/pull/638) Fix Bugs in Davis Reactants and Products. Bulk modulus is now isentropic. ### Fixed (Repair bugs, etc) +- [[PR645]](https://github.com/lanl/singularity-eos/pull/645) Fix internal energy from density and pressure for JWL - [[PR63]](https://github.com/lanl/singularity-eos/pull/639) Fixed UnitSystem temperature bug. Temperature is now treated in the same as the time, mass, and length unit factors. ### Changed (changing behavior/API/variables/...) diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 8bb333f3a55..50b6bed9287 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2025. 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 @@ -30,6 +30,8 @@ #include #include +// This file was made in part with generative AI. + namespace singularity { using namespace eos_base; From b51fb2d1fc51fb40165cbe08f89afbf719a3f876 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 16 Jul 2026 16:30:31 -0400 Subject: [PATCH 4/4] Update singularity-eos/eos/eos_jwl.hpp Co-authored-by: Adam M. Dempsey --- singularity-eos/eos/eos_jwl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_jwl.hpp b/singularity-eos/eos/eos_jwl.hpp index 50b6bed9287..135762f911f 100644 --- a/singularity-eos/eos/eos_jwl.hpp +++ b/singularity-eos/eos/eos_jwl.hpp @@ -83,7 +83,7 @@ class JWL : public EosBase { const Real rho, const Real press, Indexer_t &&lambda = static_cast(nullptr)) const { Real sie; - INternalEnergyFromDensityPressure(rho, press, sie, lambda); + InternalEnergyFromDensityPressure(rho, press, sie, lambda); return sie; } template