From 5eda3522a939cbcc72b7d08fcd68d29fa7662427 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 14:04:58 -0600 Subject: [PATCH 01/16] Fix DensityEnergyFromPressureTemperature and make portable lambda --- singularity-eos/eos/eos_gruneisen.hpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 1e8254614a..23dc59d30c 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -400,19 +400,20 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityTemperature( template PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - sie = _Cv * (temp - _T0); + // Energy is not a function of density + sie = InternalEnergyFromDensityTemperature(rho0, temp) // We have a branch at rho0, so we need to decide, based on our pressure, whether we // should be above or below rho0 Real Pref = PressureFromDensityInternalEnergy(_rho0, sie); if (press < Pref) { rho = (press - _P0 + _C0 * _C0 * _rho0) / (_C0 * _C0 + _G0 * sie); } else { // We are in compression; iterate - auto PofRatE = [&](const Real r) { - return PressureFromDensityInternalEnergy(r, sie); - }; + auto P_residual = PORTABLE_LAMBDA(const Real r) { + return press - PressureFromDensityInternalEnergy(r, sie); + } using RootFinding1D::regula_falsi; using RootFinding1D::Status; - auto status = regula_falsi(PofRatE, press, _rho0, 1.0e-5, 1.0e3, 1.0e-8, 1.0e-8, rho); + auto status = regula_falsi(P_residual, press, _rho0, 1.0e-5, 1.0e3, 1.0e-8, 1.0e-8, rho); if (status != Status::SUCCESS) { // Root finder failed even though the solution was bracketed... this is an error EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " From e2e167bf6c23b270cbf418e7050679e16ee497d5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 14:05:44 -0600 Subject: [PATCH 02/16] Make tol depend on Real type and add test for DensityEnergyFromPressureTemperature --- test/test_eos_gruneisen.cpp | 32 ++++++++++++++++++++++++-------- 1 file changed, 24 insertions(+), 8 deletions(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 995fafef4c..57e416d4db 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -33,6 +33,8 @@ PORTABLE_INLINE_FUNCTION Real QuadFormulaMinus(Real a, Real b, Real c) { return (-b - std::sqrt(b * b - 4 * a * c)) / (2 * a); } +constexpr Real REAL_TOL = std::numeric_limits::epsilon() * 1e4; // 1e-12 for double + // Only run exception checking when we aren't offloading to the GPUs #ifdef PORTABILITY_STRATEGY_NONE SCENARIO("Gruneisen EOS entropy is disabled", "[GruneisenEOS][Entropy]") { @@ -206,8 +208,8 @@ SCENARIO("Gruneisen EOS", "[VectorEOS][GruneisenEOS]") { } } -SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[GruneisenEOS]") { - GIVEN("Parameters for a Gruneisen EOS") { +SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { + GIVEN("Parameters for an aluminum Gruneisen EOS") { // Unit conversions // constexpr Real mm = 10.; constexpr Real cm = 1.; @@ -239,7 +241,7 @@ SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[Gruneis pres = pres / Mbar; INFO("Density: " << density << " Energy: " << energy << " Pressure: " << pres << " Mbar True pressure: " << true_pres << " Mbar"); - REQUIRE(isClose(pres, true_pres, 1e-12)); + REQUIRE(isClose(pres, true_pres, REAL_TOL)); } } WHEN("A B_S(rho, e) lookup is performed") { @@ -251,7 +253,21 @@ SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[Gruneis << " Sound speed: " << sound_speed << " cm/us True sound speed: " << true_sound_speed << " cm/us"); - REQUIRE(isClose(sound_speed, true_sound_speed, 1e-12)); + REQUIRE(isClose(sound_speed, true_sound_speed, REAL_TOL)); + } + } + WHEN("A pressure-temperature lookup is performed") { + const Real temperature = eos.TemperatureFromDensityInternalEnergy(density, energy); + Real test_density; + Real test_energy; + Real* lambda; + eos.DensityEnergyFromPressureTemperature(true_pres, temperature, lambda, test_density, test_energy); + THEN("The correct energy and density should be returned") { + INFO("Pressure: " << true_pres << " Temperature: " << temperature); + INFO("Density: " << density << " Energy: " << energy); + INFO("Calculated Density: " << test_density << " Calculated Energy:" << test_energy); + CHECK(isClose(density, test_density, REAL_TOL)); + CHECK(isClose(energy, test_energy, REAL_TOL)); } } WHEN("A finite difference approximation is used for the bulk modulus") { @@ -294,7 +310,7 @@ SCENARIO("Aluminum Gruneisen EOS Sound Speed and Pressure Comparison", "[Gruneis INFO("Density: " << density << " Energy: " << energy << " Pressure: " << pres / Mbar << " Mbar True pressure: " << true_pres / Mbar << " Mbar"); - REQUIRE(isClose(pres, true_pres, 1e-12)); + REQUIRE(isClose(pres, true_pres, REAL_TOL)); } } } @@ -399,7 +415,7 @@ SCENARIO("Gruneisen EOS density limit") { THEN("The generated rho_max parameter should be properly set") { const Real rho_max = eos.ComputeRhoMax(S1, S2, S3, rho0); INFO("True rho_max: " << rho_max_true << ", Calculated rho_max:" << rho_max); - REQUIRE(isClose(rho_max, rho_max_true, 1e-12)); + REQUIRE(isClose(rho_max, rho_max_true, REAL_TOL)); } WHEN("Lookups are performed beyond the maximum density") { // Note: there is a small safety factor to prevent us from hitting the @@ -495,7 +511,7 @@ SCENARIO("Gruneisen EOS density limit") { const Real rho_max_true = rho0 / (1 - eta_max); const Real rho_max = eos.ComputeRhoMax(S1, S2, S3, rho0); INFO("True rho_max: " << rho_max_true << ", Calculated rho_max:" << rho_max); - REQUIRE(isClose(rho_max, rho_max_true, 1e-12)); + REQUIRE(isClose(rho_max, rho_max_true, REAL_TOL)); } } WHEN("The quadratic oot multiplicity is 2") { @@ -511,7 +527,7 @@ SCENARIO("Gruneisen EOS density limit") { const Real rho_max_true = rho0 / (1 - eta_max); const Real rho_max = eos.ComputeRhoMax(S1, S2, S3, rho0); INFO("True rho_max: " << rho_max_true << ", Calculated rho_max:" << rho_max); - REQUIRE(isClose(rho_max, rho_max_true, 1e-12)); + REQUIRE(isClose(rho_max, rho_max_true, REAL_TOL)); } } WHEN("No root exists") { From 119137dabc73b27f7f384630f053b30eb9288821 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 14:06:38 -0600 Subject: [PATCH 03/16] Clang format --- singularity-eos/eos/eos_gruneisen.hpp | 9 +++++---- test/test_eos_gruneisen.cpp | 14 +++++++++----- 2 files changed, 14 insertions(+), 9 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 23dc59d30c..7759ee1ff5 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -402,9 +402,9 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { // Energy is not a function of density sie = InternalEnergyFromDensityTemperature(rho0, temp) - // We have a branch at rho0, so we need to decide, based on our pressure, whether we - // should be above or below rho0 - Real Pref = PressureFromDensityInternalEnergy(_rho0, sie); + // We have a branch at rho0, so we need to decide, based on our pressure, whether we + // should be above or below rho0 + Real Pref = PressureFromDensityInternalEnergy(_rho0, sie); if (press < Pref) { rho = (press - _P0 + _C0 * _C0 * _rho0) / (_C0 * _C0 + _G0 * sie); } else { // We are in compression; iterate @@ -413,7 +413,8 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( } using RootFinding1D::regula_falsi; using RootFinding1D::Status; - auto status = regula_falsi(P_residual, press, _rho0, 1.0e-5, 1.0e3, 1.0e-8, 1.0e-8, rho); + auto status = + regula_falsi(P_residual, press, _rho0, 1.0e-5, 1.0e3, 1.0e-8, 1.0e-8, rho); if (status != Status::SUCCESS) { // Root finder failed even though the solution was bracketed... this is an error EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 57e416d4db..d6401abd57 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -257,15 +257,19 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { } } WHEN("A pressure-temperature lookup is performed") { - const Real temperature = eos.TemperatureFromDensityInternalEnergy(density, energy); + const Real temperature = + eos.TemperatureFromDensityInternalEnergy(density, energy); Real test_density; Real test_energy; - Real* lambda; - eos.DensityEnergyFromPressureTemperature(true_pres, temperature, lambda, test_density, test_energy); + Real *lambda; + eos.DensityEnergyFromPressureTemperature(true_pres, temperature, lambda, + test_density, test_energy); THEN("The correct energy and density should be returned") { - INFO("Pressure: " << true_pres << " Temperature: " << temperature); + INFO("Pressure: " << true_pres + << " Temperature: " << temperature); INFO("Density: " << density << " Energy: " << energy); - INFO("Calculated Density: " << test_density << " Calculated Energy:" << test_energy); + INFO("Calculated Density: " << test_density + << " Calculated Energy:" << test_energy); CHECK(isClose(density, test_density, REAL_TOL)); CHECK(isClose(energy, test_energy, REAL_TOL)); } From 0cb709ae622c4461a678fcb80ec651b6e6aa306c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 14:12:06 -0600 Subject: [PATCH 04/16] Fix some bugs --- singularity-eos/eos/eos_gruneisen.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 7759ee1ff5..de7cc4ee44 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -401,7 +401,7 @@ template PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { // Energy is not a function of density - sie = InternalEnergyFromDensityTemperature(rho0, temp) + sie = InternalEnergyFromDensityTemperature(_rho0, temp); // We have a branch at rho0, so we need to decide, based on our pressure, whether we // should be above or below rho0 Real Pref = PressureFromDensityInternalEnergy(_rho0, sie); @@ -410,7 +410,7 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( } else { // We are in compression; iterate auto P_residual = PORTABLE_LAMBDA(const Real r) { return press - PressureFromDensityInternalEnergy(r, sie); - } + }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; auto status = From 8bb528fed4116871a4bcde7cc40d806d8cc0d412 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 16:39:01 -0600 Subject: [PATCH 05/16] Change bounds on root find to be actually physical --- singularity-eos/eos/eos_gruneisen.hpp | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index de7cc4ee44..ae8b30533f 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -400,26 +400,26 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityTemperature( template PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - // Energy is not a function of density - sie = InternalEnergyFromDensityTemperature(_rho0, temp); - // We have a branch at rho0, so we need to decide, based on our pressure, whether we - // should be above or below rho0 - Real Pref = PressureFromDensityInternalEnergy(_rho0, sie); + // We have a branch at rho0, so we need to decide, based on our pressure, whether we + // should be above or below rho0 + Real Pref = PressureFromDensityTemperature(_rho0, temp); if (press < Pref) { + sie = InternalEnergyFromDensityTemperature(_rho0, temp); rho = (press - _P0 + _C0 * _C0 * _rho0) / (_C0 * _C0 + _G0 * sie); } else { // We are in compression; iterate - auto P_residual = PORTABLE_LAMBDA(const Real r) { - return press - PressureFromDensityInternalEnergy(r, sie); + auto PofRatT = PORTABLE_LAMBDA(const Real r) { + return PressureFromDensityTemperature(r, temp); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; - auto status = - regula_falsi(P_residual, press, _rho0, 1.0e-5, 1.0e3, 1.0e-8, 1.0e-8, rho); + auto status = regula_falsi(PofRatT, press, _rho0, 0.9 * _rho0, + std::min(_rho_max, 1.0e4), 1.0e-8, 1.0e-8, rho); if (status != Status::SUCCESS) { // Root finder failed even though the solution was bracketed... this is an error EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " "Root find failed to find a solution given P, T\n"); } + sie = InternalEnergyFromDensityTemperature(_rho0, temp); } } template From 185e58883ae2f7b1edf784e1422042bbf8605424 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 16:39:37 -0600 Subject: [PATCH 06/16] Change structure of test --- test/test_eos_gruneisen.cpp | 32 ++++++++++++++++++-------------- 1 file changed, 18 insertions(+), 14 deletions(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index d6401abd57..7af4b879d2 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -256,22 +256,26 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { REQUIRE(isClose(sound_speed, true_sound_speed, REAL_TOL)); } } - WHEN("A pressure-temperature lookup is performed") { + WHEN("A the pressure and temperature are determined from the density and energy") { const Real temperature = eos.TemperatureFromDensityInternalEnergy(density, energy); - Real test_density; - Real test_energy; - Real *lambda; - eos.DensityEnergyFromPressureTemperature(true_pres, temperature, lambda, - test_density, test_energy); - THEN("The correct energy and density should be returned") { - INFO("Pressure: " << true_pres - << " Temperature: " << temperature); - INFO("Density: " << density << " Energy: " << energy); - INFO("Calculated Density: " << test_density - << " Calculated Energy:" << test_energy); - CHECK(isClose(density, test_density, REAL_TOL)); - CHECK(isClose(energy, test_energy, REAL_TOL)); + const Real pressure = eos.PressureFromDensityInternalEnergy(density, energy); + AND_WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") { + Real test_density; + Real test_energy; + Real *lambda; + eos.DensityEnergyFromPressureTemperature(pressure, temperature, lambda, + test_density, test_energy); + THEN("The correct energy and density should be returned") { + INFO("Pressure: " << pressure << " microbar" + << " Temperature: " << temperature << " K "); + INFO("Density: " << density << " g/cm^3 " + << " Energy: " << energy << " erg/g "); + INFO("Calc Density: " << test_density << " g/cm^3 " + << " Calc Energy: " << test_energy << " erg/g "); + CHECK(isClose(density, test_density, REAL_TOL)); + CHECK(isClose(energy, test_energy, REAL_TOL)); + } } } WHEN("A finite difference approximation is used for the bulk modulus") { From 390a7ed8ada3dc44b76203e0c1200c90f672f9dd Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 17:36:58 -0600 Subject: [PATCH 07/16] Make REfromPT iterate in expansion too --- singularity-eos/eos/eos_gruneisen.hpp | 37 +++++++++++++++------------ 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index ae8b30533f..b96a32c073 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -403,24 +403,29 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( // We have a branch at rho0, so we need to decide, based on our pressure, whether we // should be above or below rho0 Real Pref = PressureFromDensityTemperature(_rho0, temp); + Real rho_lower; + Real rho_upper; + // Pick bounds appropriate depending on whether in compression or expansion if (press < Pref) { - sie = InternalEnergyFromDensityTemperature(_rho0, temp); - rho = (press - _P0 + _C0 * _C0 * _rho0) / (_C0 * _C0 + _G0 * sie); - } else { // We are in compression; iterate - auto PofRatT = PORTABLE_LAMBDA(const Real r) { - return PressureFromDensityTemperature(r, temp); - }; - using RootFinding1D::regula_falsi; - using RootFinding1D::Status; - auto status = regula_falsi(PofRatT, press, _rho0, 0.9 * _rho0, - std::min(_rho_max, 1.0e4), 1.0e-8, 1.0e-8, rho); - if (status != Status::SUCCESS) { - // Root finder failed even though the solution was bracketed... this is an error - EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " - "Root find failed to find a solution given P, T\n"); - } - sie = InternalEnergyFromDensityTemperature(_rho0, temp); + rho_lower = 0.; + rho_upper = 1.1 * _rho0; + } else { + rho_lower = 0.9 * _rho0; + rho_upper = std::min(_rho_max, 1.0e4); + } + auto PofRatT = PORTABLE_LAMBDA(const Real r) { + return PressureFromDensityTemperature(r, temp); + }; + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + auto status = regula_falsi(PofRatT, press, _rho0, 0.9 * _rho0, + std::min(_rho_max, 1.0e4), 1.0e-8, 1.0e-8, rho); + if (status != Status::SUCCESS) { + // Root finder failed even though the solution was bracketed... this is an error + EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " + "Root find failed to find a solution given P, T\n"); } + sie = InternalEnergyFromDensityTemperature(rho, temp); } template PORTABLE_INLINE_FUNCTION void From 3889658fd7ecbae90d79fbe482956373d7c0752c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 17:47:30 -0600 Subject: [PATCH 08/16] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6c4142d09c..3f03b8d214 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ ### Fixed (Repair bugs, etc) - [[PR401]](https://github.com/lanl/singularity-eos/pull/401) Fix for internal energy scaling in PTE closure +- [[PR403]](https://github.com/lanl/singularity-eos/pull/403) Fix Gruneisen EOS DensityEnergyFromPressureTemperature function ### Changed (changing behavior/API/variables/...) From 148bc6f09b1ed26807d303f1e9dd5cd651334d08 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 8 Aug 2024 17:47:41 -0600 Subject: [PATCH 09/16] Update copyright --- test/test_eos_gruneisen.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 7af4b879d2..915c10777c 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2023. Triad National Security, LLC. All rights reserved. This +// © 2021-2024. 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 From 3b1c6fff0723dd8840fbe15fffaaf98e6701d57d Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 9 Aug 2024 16:21:07 -0600 Subject: [PATCH 10/16] Add ability to compute maximum stable density --- singularity-eos/eos/eos_gruneisen.hpp | 60 +++++++++++++++++++++++++-- 1 file changed, 57 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index b96a32c073..564d1b596d 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -63,6 +63,8 @@ class Gruneisen : public EosBase { _Cv(Cv), _rho_max(RHOMAX_SAFETY * ComputeRhoMax(s1, s2, s3, rho0)) {} static PORTABLE_INLINE_FUNCTION Real ComputeRhoMax(const Real s1, const Real s2, const Real s3, const Real rho0); + PORTABLE_INLINE_FUNCTION Real + MaxStableDensityAtTemperature(const Real temperature) const; Gruneisen GetOnDevice() { return *this; } template PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( @@ -397,6 +399,48 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::BulkModulusFromDensityTemperature( return BulkModulusFromDensityInternalEnergy( rho, InternalEnergyFromDensityTemperature(rho, temp)); } +PORTABLE_INLINE_FUNCTION Real +Gruneisen::MaxStableDensityAtTemperature(const Real temperature) const { + // Because of the constant cold curve assumption, this EOS is thermodynamically + // inconsistent, and leads to states that are effectively off the EOS surface even + // though the temperature is positive. Mathematically, this means that the \Gamma\rho(e + // - e_H) term becomes highly negative and dominates the positive P_H term at some + // density. This results in a local maximum in the pressure as a function of density for + // a given temperature. Beyond this point, the pressure decreases with increasing + // density, which is thermodynamcially unstable. In a thermodynamically consistent EOS, + // the cold curve energy would increase with density, leading an appropriate bounding at + // T=0. + + // Since E_H and P_H are monotonic up to the singularity given by _rho_max, if the + // derivative of the pressure is negative at _rho_max, a maximum exists and this should + // be the highest density for the isotherm. If this maximum doesn't exist, then the EOS + // is thermodynamically consistent up to the maximum density for this isotherm. + Real slope_at_max_density = + dPres_drho_e(_rho_max, InternalEnergyFromDensityTemperature(_rho_max, temperature)); + if (slope_at_max_density >= 0) { + // No maximum pressure before _rho_max + return _rho_max; + } + + // Maximum pressure should exist... do a root find to locate where the derivative is + // zero + auto dPdrho_T = PORTABLE_LAMBDA(const Real r) { + return dPres_drho_e(r, InternalEnergyFromDensityTemperature(r, temperature)); + }; + Real rho_lower = 0.9 * _rho0; + Real rho_upper = std::min(_rho_max, 1.0e4); + Real rho_at_max_P; + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + auto status = regula_falsi(dPdrho_T, 0., _rho0, rho_lower, rho_upper, 1.0e-8, 1.0e-8, + rho_at_max_P); + if (status != Status::SUCCESS) { + // Root finder failed even though the solution should be bracketed + EOS_ERROR("Gruneisen::MaxStableDensityAtTemperature: " + "Root find failed to find maximum P at T despite aparent bracket\n"); + } + return rho_at_max_P; +} template PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { @@ -411,15 +455,25 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( rho_upper = 1.1 * _rho0; } else { rho_lower = 0.9 * _rho0; - rho_upper = std::min(_rho_max, 1.0e4); + // Find maximum thermodynamically _stable_ density at this temperature. Use this to + // check if we're actually on the EOS surface or not + rho_upper = MaxStableDensityAtTemperature(temp); + auto pres_max = PressureFromDensityTemperature(rho_upper, temp); + if (press > pres_max) { + // We're off the EOS surface + using PortsOfCall::printf; + printf("Requested pressure, %.15g, exceeds maximum, %.15g, for temperature, %.15g", + press, pres_max, temp); + PORTABLE_ALWAYS_THROW_OR_ABORT("Input pressure is off EOS surface"); + } } auto PofRatT = PORTABLE_LAMBDA(const Real r) { return PressureFromDensityTemperature(r, temp); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; - auto status = regula_falsi(PofRatT, press, _rho0, 0.9 * _rho0, - std::min(_rho_max, 1.0e4), 1.0e-8, 1.0e-8, rho); + auto status = + regula_falsi(PofRatT, press, _rho0, rho_lower, rho_upper, 1.0e-8, 1.0e-8, rho); if (status != Status::SUCCESS) { // Root finder failed even though the solution was bracketed... this is an error EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " From 5663445f1a504d2981fae3458f4f73dd44c03839 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 9 Aug 2024 16:23:32 -0600 Subject: [PATCH 11/16] Update aluminum parameters to trigger failure and adjust tests --- test/test_eos_gruneisen.cpp | 43 ++++++++++++++++++++++++++----------- 1 file changed, 30 insertions(+), 13 deletions(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 915c10777c..b2791adf35 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -216,25 +216,42 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { constexpr Real us = 1.e-06; constexpr Real Mbar = 1.e12; constexpr Real Mbcc_per_g = 1e12; - // Gruneisen parameters for copper - constexpr Real C0 = 0.535 * cm / us; - constexpr Real S1 = 1.34; + // Gruneisen parameters for aluminum + constexpr Real C0 = 0.524 * cm / us; + constexpr Real S1 = 1.4; constexpr Real S2 = 0.; constexpr Real S3 = 0.; constexpr Real Gamma0 = 1.97; - constexpr Real b = 0.; - constexpr Real rho0 = 2.714000; + constexpr Real b = 0.48; + constexpr Real rho0 = 2.703; constexpr Real T0 = 298.; - constexpr Real P0 = 1e-06 * Mbar; + constexpr Real P0 = 0. * Mbar; constexpr Real Cv = 0.383e-05 * Mbcc_per_g; // Create the EOS EOS host_eos = Gruneisen(C0, S1, S2, S3, Gamma0, b, rho0, T0, P0, Cv); EOS eos = host_eos.GetOnDevice(); + GIVEN("Ambient pressure and temperature") { + constexpr Real P_ambient = 1.0e-06 * Mbar; + constexpr Real T_ambient = 298.; + WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") { + Real test_density; + Real test_energy; + Real *lambda; + THEN("The root find should converge") { + eos.DensityEnergyFromPressureTemperature(P_ambient, T_ambient, lambda, + test_density, test_energy); + } + } + } GIVEN("Density and energy") { - constexpr Real density = 5.92418956756592; // g/cm^3 - constexpr Real energy = 792486007.804619; // erg/g - constexpr Real true_pres = 2.620656373250729; // Mbar - constexpr Real true_sound_speed = 1.5247992468363685; // cm/us + constexpr Real density = 5.92418956756592; // g/cm^3 + constexpr Real energy = 792486007.804619; // erg/g + constexpr Real true_pres = 2.19200396047868; // Mbar + constexpr Real true_sound_speed = 1.29592658233752; // cm/us + THEN("The density should not exceeed the computed rho_max") { + const Real density_max = Gruneisen::ComputeRhoMax(S1, S2, S3, rho0); + REQUIRE(density < density_max); + } WHEN("A P(rho, e) lookup is performed") { Real pres = eos.PressureFromDensityInternalEnergy(density, energy); THEN("The correct pressure should be returned") { @@ -256,17 +273,17 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { REQUIRE(isClose(sound_speed, true_sound_speed, REAL_TOL)); } } - WHEN("A the pressure and temperature are determined from the density and energy") { + WHEN("The pressure and temperature are determined from the density and energy") { const Real temperature = eos.TemperatureFromDensityInternalEnergy(density, energy); - const Real pressure = eos.PressureFromDensityInternalEnergy(density, energy); + const Real pressure = eos.PressureFromDensityInternalEnergy(density, energy); AND_WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") { Real test_density; Real test_energy; Real *lambda; eos.DensityEnergyFromPressureTemperature(pressure, temperature, lambda, test_density, test_energy); - THEN("The correct energy and density should be returned") { + THEN("The consistent energy and density should be returned") { INFO("Pressure: " << pressure << " microbar" << " Temperature: " << temperature << " K "); INFO("Density: " << density << " g/cm^3 " From e4e1c67e8939d923080550a60fd1d693ae3fc249 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 9 Aug 2024 16:31:02 -0600 Subject: [PATCH 12/16] Clang format --- test/test_eos_gruneisen.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index b2791adf35..6499043f6b 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -276,7 +276,7 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { WHEN("The pressure and temperature are determined from the density and energy") { const Real temperature = eos.TemperatureFromDensityInternalEnergy(density, energy); - const Real pressure = eos.PressureFromDensityInternalEnergy(density, energy); + const Real pressure = eos.PressureFromDensityInternalEnergy(density, energy); AND_WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") { Real test_density; Real test_energy; From 1032d9ce0d6226ce0d16ae7099b48036b6d2a586 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 12 Aug 2024 18:24:17 -0600 Subject: [PATCH 13/16] Make better initial guesses --- singularity-eos/eos/eos_gruneisen.hpp | 15 ++++++++++----- 1 file changed, 10 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 564d1b596d..590efd55ac 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -427,13 +427,14 @@ Gruneisen::MaxStableDensityAtTemperature(const Real temperature) const { auto dPdrho_T = PORTABLE_LAMBDA(const Real r) { return dPres_drho_e(r, InternalEnergyFromDensityTemperature(r, temperature)); }; - Real rho_lower = 0.9 * _rho0; - Real rho_upper = std::min(_rho_max, 1.0e4); + const Real rho_lower = _rho0; + const Real rho_upper = std::min(_rho_max, 1.0e4); + const Real rho_guess = (rho_lower + rho_upper) / 2.; Real rho_at_max_P; using RootFinding1D::regula_falsi; using RootFinding1D::Status; - auto status = regula_falsi(dPdrho_T, 0., _rho0, rho_lower, rho_upper, 1.0e-8, 1.0e-8, - rho_at_max_P); + auto status = regula_falsi(dPdrho_T, 0., rho_guess, rho_lower, rho_upper, 1.0e-8, + 1.0e-8, rho_at_max_P); if (status != Status::SUCCESS) { // Root finder failed even though the solution should be bracketed EOS_ERROR("Gruneisen::MaxStableDensityAtTemperature: " @@ -449,16 +450,20 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( Real Pref = PressureFromDensityTemperature(_rho0, temp); Real rho_lower; Real rho_upper; + Real rho_guess; // Pick bounds appropriate depending on whether in compression or expansion if (press < Pref) { rho_lower = 0.; rho_upper = 1.1 * _rho0; + rho_guess = _rho0; } else { rho_lower = 0.9 * _rho0; // Find maximum thermodynamically _stable_ density at this temperature. Use this to // check if we're actually on the EOS surface or not rho_upper = MaxStableDensityAtTemperature(temp); auto pres_max = PressureFromDensityTemperature(rho_upper, temp); + const Real slope = (rho_upper - _rho0) / (pres_max - Pref); + rho_guess = _rho0 + slope * (press - Pref); if (press > pres_max) { // We're off the EOS surface using PortsOfCall::printf; @@ -473,7 +478,7 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( using RootFinding1D::regula_falsi; using RootFinding1D::Status; auto status = - regula_falsi(PofRatT, press, _rho0, rho_lower, rho_upper, 1.0e-8, 1.0e-8, rho); + regula_falsi(PofRatT, press, rho_guess, rho_lower, rho_upper, 1.0e-8, 1.0e-8, rho); if (status != Status::SUCCESS) { // Root finder failed even though the solution was bracketed... this is an error EOS_ERROR("Gruneisen::DensityEnergyFromPressureTemperature: " From a721e045b52eaf1515330d54876f01b797f8fd78 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 12 Aug 2024 18:24:44 -0600 Subject: [PATCH 14/16] Adjust initial conditions --- test/test_eos_gruneisen.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 6499043f6b..376760c64b 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -216,6 +216,7 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { constexpr Real us = 1.e-06; constexpr Real Mbar = 1.e12; constexpr Real Mbcc_per_g = 1e12; + constexpr Real K_per_eV = 1.602176565E-19 / 1.380648800E-23; // Gruneisen parameters for aluminum constexpr Real C0 = 0.524 * cm / us; constexpr Real S1 = 1.4; @@ -231,8 +232,8 @@ SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { EOS host_eos = Gruneisen(C0, S1, S2, S3, Gamma0, b, rho0, T0, P0, Cv); EOS eos = host_eos.GetOnDevice(); GIVEN("Ambient pressure and temperature") { - constexpr Real P_ambient = 1.0e-06 * Mbar; - constexpr Real T_ambient = 298.; + constexpr Real P_ambient = 1.0e6; // 1 bar + constexpr Real T_ambient = 0.025 * K_per_eV; // approximately 290 K WHEN("A DensityEnergyFromPressureTemperature() lookup is performed") { Real test_density; Real test_energy; From 69470b64295d5c5a374150a0c352104074c70ea5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 14 Aug 2024 12:41:03 -0600 Subject: [PATCH 15/16] Error out when all pressures are unstable above the reference density --- singularity-eos/eos/eos_gruneisen.hpp | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 590efd55ac..6d295d5eeb 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -417,10 +417,21 @@ Gruneisen::MaxStableDensityAtTemperature(const Real temperature) const { // is thermodynamically consistent up to the maximum density for this isotherm. Real slope_at_max_density = dPres_drho_e(_rho_max, InternalEnergyFromDensityTemperature(_rho_max, temperature)); + Real slope_at_ref_density = + dPres_drho_e(_rho0, InternalEnergyFromDensityTemperature(_rho0, temperature)); if (slope_at_max_density >= 0) { // No maximum pressure before _rho_max return _rho_max; } + if (slope_at_ref_density < 0) { + // Something is very wrong in the construction of this EOS... just error out + using PortsOfCall::printf; + printf("ERROR: The pressure is decreasing as density increases at the reference\n" + " density, %.15g, for temperature, %.15g. Check that the reference\n" + " temperature is set correctly. This is an unstable state." _rho0, + temp); + PORTABLE_ALWAYS_THROW_OR_ABORT("Input pressure is off EOS surface"); + } // Maximum pressure should exist... do a root find to locate where the derivative is // zero @@ -462,15 +473,17 @@ PORTABLE_INLINE_FUNCTION void Gruneisen::DensityEnergyFromPressureTemperature( // check if we're actually on the EOS surface or not rho_upper = MaxStableDensityAtTemperature(temp); auto pres_max = PressureFromDensityTemperature(rho_upper, temp); - const Real slope = (rho_upper - _rho0) / (pres_max - Pref); - rho_guess = _rho0 + slope * (press - Pref); if (press > pres_max) { // We're off the EOS surface using PortsOfCall::printf; - printf("Requested pressure, %.15g, exceeds maximum, %.15g, for temperature, %.15g", + printf("ERROR: Requested pressure, %.15g, exceeds maximum, %.15g, for \n" + " temperature, %.15g", press, pres_max, temp); PORTABLE_ALWAYS_THROW_OR_ABORT("Input pressure is off EOS surface"); } + // Construct a reasonable guess for the density + const Real slope = (rho_upper - _rho0) / (pres_max - Pref); + rho_guess = _rho0 + slope * (press - Pref); } auto PofRatT = PORTABLE_LAMBDA(const Real r) { return PressureFromDensityTemperature(r, temp); From 8e0977b146261c20ee98a2264dabc2714d14247c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Wed, 14 Aug 2024 12:59:12 -0600 Subject: [PATCH 16/16] Add missing comma --- singularity-eos/eos/eos_gruneisen.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index 6d295d5eeb..f5ff57c958 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -428,8 +428,8 @@ Gruneisen::MaxStableDensityAtTemperature(const Real temperature) const { using PortsOfCall::printf; printf("ERROR: The pressure is decreasing as density increases at the reference\n" " density, %.15g, for temperature, %.15g. Check that the reference\n" - " temperature is set correctly. This is an unstable state." _rho0, - temp); + " temperature is set correctly. This is an unstable state.", + _rho0, temperature); PORTABLE_ALWAYS_THROW_OR_ABORT("Input pressure is off EOS surface"); }