From 4a210b571a4ede748e36f2bb233548042c9df18b Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 23 Dec 2024 12:28:13 -0700 Subject: [PATCH 01/45] add RhoE of PT to base class and implement where relevant --- singularity-eos/eos/eos_base.hpp | 32 ++++++++++++++++++++ singularity-eos/eos/eos_helmholtz.hpp | 29 +++++++++++++++--- singularity-eos/eos/eos_mgusup.hpp | 12 -------- singularity-eos/eos/eos_powermg.hpp | 12 -------- singularity-eos/eos/eos_sap_polynomial.hpp | 8 ----- singularity-eos/eos/eos_stellar_collapse.hpp | 20 +++++++++++- singularity-eos/eos/eos_vinet.hpp | 12 -------- 7 files changed, 75 insertions(+), 50 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 5ca9ed928e1..17ecc5997be 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -836,6 +836,38 @@ class EosBase { PORTABLE_ALWAYS_THROW_OR_ABORT(msg); } + // JMM: This method is often going to be overloaded for special cases. + PORTABLE_INLINE_FUNCTION + template + void DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, + Real &sie) const { + + constexpr Real MAXFAC = 1e8; + constexpr Real EPS = 1e-8; + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + CRTP copy = *(static_cast(this)); + auto PofRT = [&](const Real r) { + return copy.PressureFromDensityTemperature(r, temp, lambda); + }; + Real rhoguess = rho; // use input density + if (rhoguess < 0) { + Real tguess, sieguess, pguess, cvguess, bmodguess, dpde, dvdt; + copy.ValuesAtReferenceState(rhoguess, tguess, sieguess, pguess, cvguess, bmodguess, + dpde, dvdt); + } + Real rhomin = copyMinimumDensity(); + Real rhomax = MAXFAC * rhomin; + auto status = regula_falsi(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); + if (status != Status::SUCCESS) { + PORTABLE_FAIL_OR_ABORT( + "DensityEnergyFromPressureTemperature failed to find root\n"); + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + return; + } + // Serialization /* The methodology here is there are *three* size methods all EOS's provide: diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index b6ade5ac19e..7b0f9de4d52 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/singularity-eos/eos/eos_helmholtz.hpp @@ -676,16 +676,35 @@ class Helmholtz : public EosBase { ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const { - // JMM: I'm not sure what to put here or if it matters. Some - // reference state, maybe stellar denity, would be appropriate. - PORTABLE_ALWAYS_ABORT("Stub"); + // TODO(JMM): Is this right??? + rho = 1e7; + temp = 1.5e9; + FillEos(rho, temp, sie, press, cv, bmod, + thermalqs::specific_internal_energy | thermalqs::pressure | + thermalqs::specific_heat | thermalqs::bulk_modulus, + lambda); + dpde = GruneisenParamFromDensityTemperature(rho, temp, lambda) * rho; + dvdt = 0; } template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - // This is only used for mixed cell closures. Stubbing it out for now. - PORTABLE_ALWAYS_ABORT("Stub"); + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + PORTABLE_REQUIRE(temp > = 0, "Non-negative temperature required"); + PofRT = [&](const Real r) { return PressureFromDensityTemperature(r, temp, lambda); }; + auto status = regula_falsi(PofRT, press, 1e7, electrons_.rhoMin(), + electrons_.rhoMax(), 1.0e-8, 1.0e-8, rho); + if (status != Status::SUCCESS) { + PORTABLE_THROW_OR_ABORT("Helmholtz::DensityEnergyFromPressureTemperature: " + "Root find failed to find a solution given P, T\n"); + } + if (rho < 0.) { + PORTABLE_THROW_OR_ABORT("Helmholtz::DensityEnergyFromPressureTemperature: " + "Root find produced negative energy solution\n"); + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); } static std::string EosType() { return std::string("Helmholtz"); } static std::string EosPyType() { return EosType(); } diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index c315c1619aa..491c301c63a 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -148,11 +148,6 @@ class MGUsup : public EosBase { _AZbar.PrintParams(); printf("\n\n"); } - // Density/Energy from P/T not unique, if used will give error - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const; inline void Finalize() {} static std::string EosType() { return std::string("MGUsup"); } static std::string EosPyType() { return EosType(); } @@ -350,13 +345,6 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::BulkModulusFromDensityInternalEnergy( value = robust::ratio(_rho0, rho) * value; return value; } -// AEM: Give error since function is not well defined -template -PORTABLE_INLINE_FUNCTION void MGUsup::DensityEnergyFromPressureTemperature( - const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("MGUsup::DensityEnergyFromPressureTemperature: " - "Not implemented.\n"); -} // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. template diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index 7e400b89a7c..ed89d662ec1 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -157,11 +157,6 @@ class PowerMG : public EosBase { } printf("\n\n"); } - // Density/Energy from P/T not unique, if used will give error - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const; inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } static std::string EosPyType() { return EosType(); } @@ -431,13 +426,6 @@ PORTABLE_INLINE_FUNCTION Real PowerMG::EntropyFromDensityInternalEnergy( } return value; } -// AEM: Give error since function is not well defined -template -PORTABLE_INLINE_FUNCTION void PowerMG::DensityEnergyFromPressureTemperature( - const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("PowerMG::DensityEnergyFromPressureTemperature: " - "Not implemented.\n"); -} // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. template diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp index 7d361190c9c..081e434d0fc 100644 --- a/singularity-eos/eos/eos_sap_polynomial.hpp +++ b/singularity-eos/eos/eos_sap_polynomial.hpp @@ -209,14 +209,6 @@ class SAP_Polynomial : public EosBase { printf(" b3 = %g\n", _b3); _AZbar.PrintParams(); } - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const { - PORTABLE_WARN("This function is a stub for an incomplete EoS."); - sie = 0.0; - rho = 0.0; - } inline void Finalize() {} static std::string EosType() { return std::string("SAP_Polynomial"); } static std::string EosPyType() { return EosType(); } diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 4ff7788dda2..ee282136c75 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -653,7 +653,25 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::GruneisenParamFromDensityInternal template PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("StellarCollapse::DensityEnergyFromPRessureTemperature is a stub"); + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + CheckLambda_(lambda); + Real lrguess = lRho_(rho); + Real lT = lT_(temp); + Real lP = P2lP_(press); + Real Ye = lambda[Lambda::Ye]; + + if ((lrguess < lRhoMin_) || (lrguess > lrMax_)) { + lrguess = lRho_(rhoNormal_); + } + auto lPofRT = [&](Real lR) { return lP_.interpToReal(Ye, lT, lR); }; + auto status = regula_falsi(lPofRT, lP, lrguess, lRhoMin_, lRhoMax_, ROOT_THRESH, + ROOT_THRESH, lrguess); + + Real lE = lE_.interpToReal(Ye, lT, lrguess); + rho = rho_(lrguess); + sie = le2e_(lE); + lambda[Lambda::lT] = lT; } template diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 8ea913e4049..c942735d1bc 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -152,11 +152,6 @@ class Vinet : public EosBase { } printf("\n\n"); } - // Density/Energy from P/T not unique, if used will give error - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const; inline void Finalize() {} static std::string EosType() { return std::string("Vinet"); } static std::string EosPyType() { return EosType(); } @@ -375,13 +370,6 @@ PORTABLE_INLINE_FUNCTION Real Vinet::BulkModulusFromDensityInternalEnergy( Vinet_F_DT_func(rho, temp, output); return output[7] * output[7] * rho; } -// AEM: Give error since function is not well defined -template -PORTABLE_INLINE_FUNCTION void Vinet::DensityEnergyFromPressureTemperature( - const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - EOS_ERROR("Vinet::DensityEnergyFromPressureTemperature: " - "Not implemented.\n"); -} // AEM: We should add entropy and Gruneissen parameters here so that it is complete // If we add also alpha and BT, those should also be in here. template From 50c08e7e780a7bc2f4112dadd54cf94c3e416b15 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Mon, 23 Dec 2024 12:37:34 -0700 Subject: [PATCH 02/45] formatting --- singularity-eos/eos/eos_base.hpp | 14 +++++++------- singularity-eos/eos/eos_stellar_collapse.hpp | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 17ecc5997be..4cc12aa310b 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -21,6 +21,7 @@ #include #include #include +#include #include namespace singularity { @@ -837,11 +838,10 @@ class EosBase { } // JMM: This method is often going to be overloaded for special cases. - PORTABLE_INLINE_FUNCTION - template - void DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, - Real &sie) const { + template + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { constexpr Real MAXFAC = 1e8; constexpr Real EPS = 1e-8; @@ -857,11 +857,11 @@ class EosBase { copy.ValuesAtReferenceState(rhoguess, tguess, sieguess, pguess, cvguess, bmodguess, dpde, dvdt); } - Real rhomin = copyMinimumDensity(); + Real rhomin = copy.MinimumDensity(); Real rhomax = MAXFAC * rhomin; auto status = regula_falsi(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); if (status != Status::SUCCESS) { - PORTABLE_FAIL_OR_ABORT( + PORTABLE_THROW_OR_ABORT( "DensityEnergyFromPressureTemperature failed to find root\n"); } sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index ee282136c75..8beab384d2c 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -661,7 +661,7 @@ PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperat Real lP = P2lP_(press); Real Ye = lambda[Lambda::Ye]; - if ((lrguess < lRhoMin_) || (lrguess > lrMax_)) { + if ((lrguess < lRhoMin_) || (lrguess > lRhoMax_)) { lrguess = lRho_(rhoNormal_); } auto lPofRT = [&](Real lR) { return lP_.interpToReal(Ye, lT, lR); }; From acd9deff435b6a38828d4a90a8f5a106aa7d4344 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 15:39:53 -0700 Subject: [PATCH 03/45] Tests pass on macbook. Still need to add a few tests --- singularity-eos/eos/eos_base.hpp | 12 ++++-- singularity-eos/eos/eos_mgusup.hpp | 43 +++++++++++++++++++ singularity-eos/eos/eos_stellar_collapse.hpp | 2 +- singularity-eos/eos/modifiers/zsplit_eos.hpp | 10 +++++ test/eos_unit_test_helpers.hpp | 29 +++++++++++++ test/test_eos_carnahan_starling.cpp | 12 ++++++ test/test_eos_helmholtz.cpp | 5 ++- test/test_eos_ideal.cpp | 44 ++++++++++++++++++-- test/test_eos_mgusup.cpp | 11 +++++ test/test_eos_noble_abel.cpp | 12 ++++++ test/test_eos_powermg.cpp | 11 +++++ test/test_eos_stellar_collapse.cpp | 3 ++ test/test_eos_stiff.cpp | 10 +++++ 13 files changed, 194 insertions(+), 10 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 4cc12aa310b..5cd7dbf3fbe 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -21,6 +21,7 @@ #include #include #include +#include #include #include @@ -842,9 +843,11 @@ class EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - + // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's + // overwrite. constexpr Real MAXFAC = 1e8; - constexpr Real EPS = 1e-8; + constexpr Real EPS = 10 * robust::EPS(); + constexpr Real MINR_DEFAULT = 1e-4; using RootFinding1D::regula_falsi; using RootFinding1D::Status; CRTP copy = *(static_cast(this)); @@ -857,14 +860,15 @@ class EosBase { copy.ValuesAtReferenceState(rhoguess, tguess, sieguess, pguess, cvguess, bmodguess, dpde, dvdt); } - Real rhomin = copy.MinimumDensity(); + // JMM: This can't be zero, in case MinimumDensity is zero + Real rhomin = std::max(MINR_DEFAULT, copy.MinimumDensity()); Real rhomax = MAXFAC * rhomin; auto status = regula_falsi(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); if (status != Status::SUCCESS) { PORTABLE_THROW_OR_ABORT( "DensityEnergyFromPressureTemperature failed to find root\n"); } - sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda); return; } diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index 491c301c63a..72dfcec805f 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -125,6 +125,47 @@ class MGUsup : public EosBase { FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda = static_cast(nullptr)) const; + + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return 10 * robust::EPS(); } + + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { + if (_s > 1) { + return 0.99 * robust::ratio(_s * _rho0, _s - 1); + } else { // for s <= 1, no maximum, but we need to pick something. + return 1e3 * _rho0; + } + } + + template + PORTABLE_INLINE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { + using RootFinding1D::findRoot; + using RootFinding1D::Status; + // JMM: diverges for rho -> 0 + // and for s > 1 and rho -> rho0 + Real rhomin = MinimumDensity(); + Real rhomax = MaximumDensity(); + auto PofRT = [&](const Real r) { + return PressureFromDensityTemperature(r, temp, lambda); + }; + Real rhoguess = rho; // use input density + if ((rhoguess < rhomin) || (rhoguess > rhomax)) { + rhoguess = 0.5 * (rhomin + rhomax); + } + // JMM: regula_falsi does not respect bounds + auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), + robust::EPS(), rho); + if (status != Status::SUCCESS) { + PORTABLE_THROW_OR_ABORT( + "DensityEnergyFromPressureTemperature failed to find root\n"); + } + sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); + return; + } + template PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, @@ -195,6 +236,8 @@ PORTABLE_INLINE_FUNCTION Real MGUsup::HugTemperatureFromDensity(Real rho) const Real eta = 1.0 - robust::ratio(_rho0, rho); Real f1 = 1.0 - _s * eta; if (f1 <= 0.0) { + printf("f1, eta, rho, rho0, s = %.14e %.14e %.14e %.14e %.14e\n", f1, eta, rho, _rho0, + _s); PORTABLE_ALWAYS_THROW_OR_ABORT("MGUsup model parameters s and rho0 together with rho " "give a negative argument for a logarithm."); } diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 8beab384d2c..1ff556fadeb 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -655,7 +655,7 @@ PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperat const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { using RootFinding1D::regula_falsi; using RootFinding1D::Status; - CheckLambda_(lambda); + checkLambda_(lambda); Real lrguess = lRho_(rho); Real lT = lT_(temp); Real lP = P2lP_(press); diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index eb9d68313cf..cef562a34f8 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -212,6 +212,16 @@ class ZSplit : public EosBase> { // dvdt, dpde unchanged } + template + PORTABLE_FUNCTION void + DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const { + const Real scale = GetScale_(lambda); + const Real iscale = GetInvScale_(lambda); + t_.DensityEnergyFromPressureTemperature(iscale * press, temp, lambda, rho, sie); + sie *= scale; + } + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return NL + t_.nlambda(); } static constexpr unsigned long PreferredInput() { return T::PreferredInput(); } diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index b5bf7bde90f..2e84dd792f3 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -34,6 +34,7 @@ #include #include +#include inline std::string demangle(const char *name) { @@ -143,6 +144,34 @@ inline void compare_two_eoss(const E1 &test_e, const E2 &ref_e) { return; } +template +PORTABLE_INLINE_FUNCTION bool +CheckRhoSieFromPT(EOS eos, Real rho, Real T, + Indexer_t &&lambda = static_cast(nullptr)) { + const Real P = eos.PressureFromDensityTemperature(rho, T, lambda); + const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T, lambda); + Real rtest, etest; + eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest); + Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); + bool results_good = (isClose(rho, rtest, bmod * 1e-8) || isClose(rho, rtest, 1e-8)) && + (isClose(sie, etest, bmod * 1e-8) || isClose(sie, etest, 1e-8)); + if (!results_good) { + Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); + Real residual = P_test - P; + printf("RhoSie of PT failure!\n" + "\trho_true = %.14e\n" + "\tsie_true = %.14e\n" + "\tP = %.14e\n" + "\tT = %.14e\n" + "\trho = %.14e\n" + "\tsie = %.14e\n" + "\tP_test = %.14e\n" + "\tresidual = %.14e\n", + rho, sie, P, T, rtest, etest, P_test, residual); + } + return results_good; +} + // Macro that checks for an exception or is a no-op depending on // whether or not a non-serial backend is supplied #ifdef PORTABILITY_STRATEGY_NONE diff --git a/test/test_eos_carnahan_starling.cpp b/test/test_eos_carnahan_starling.cpp index 77735704505..da6bc18771d 100644 --- a/test/test_eos_carnahan_starling.cpp +++ b/test/test_eos_carnahan_starling.cpp @@ -161,6 +161,18 @@ SCENARIO("CarnahanStarling1", "[CarnahanStarling][CarnahanStarling1]") { "Energy"); } } + + WHEN("We check RhoSie from PT") { + int nwrong = 0; + portableReduce( + "Check RhoSieFromPT", 0, 1, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += + !CheckRhoSieFromPT(eos, 7.8290736890381501e-03, 1.5320999999999999e+03); + }, + nwrong); + REQUIRE(nwrong == 0); + } } } } diff --git a/test/test_eos_helmholtz.cpp b/test/test_eos_helmholtz.cpp index 01a2aa11038..5a200825ea4 100644 --- a/test/test_eos_helmholtz.cpp +++ b/test/test_eos_helmholtz.cpp @@ -125,13 +125,16 @@ SCENARIO("Helmholtz equation of state - Table interpolation (tgiven)", "[Helmhol eos.BulkModulusFromDensityTemperature(rho_in[i], temp_in[j], lambda); Real gruen = eos.GruneisenParamFromDensityTemperature(rho_in[i], temp_in[j], lambda); - if (!isClose(ein, ein_ref[k], 1e-10)) nwrong += 1; if (!isClose(press, press_ref[k], 1e-10)) nwrong += 1; if (!isClose(cv, cv_ref[k], 1e-6)) nwrong += 1; if (!isClose(bulkmod, bulkmod_ref[k], 1e-8)) nwrong += 1; if (!isClose(gruen, gruen_ref[k], 1e-6)) nwrong += 1; + // RhoSie of PT + nwrong += CheckRhoSieFromPT(eos, rho_in[i], temp_in[i], lambda); + + // Deserialized EOS ein = eos_2.InternalEnergyFromDensityTemperature(rho_in[i], temp_in[j], lambda); press = eos_2.PressureFromDensityTemperature(rho_in[i], temp_in[j], lambda); diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index b24a94bb3fa..58a35149298 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -106,10 +106,10 @@ SCENARIO("Ideal gas mean atomic properties", portableReduce( "Check mean atomic number", 0, N, PORTABLE_LAMBDA(const int i, int &nw) { - double rho = i; - double T = 100.0 * i; - double Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); - double Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); + Real rho = i; + Real T = 100.0 * i; + Real Ab_eval = device_eos.MeanAtomicMassFromDensityTemperature(rho, T); + Real Zb_eval = device_eos.MeanAtomicNumberFromDensityTemperature(rho, T); nw += !(isClose(Ab_eval, Abar, 1e-12)) + !(isClose(Zb_eval, Zbar, 1e-12)); }, nwrong); @@ -120,6 +120,29 @@ SCENARIO("Ideal gas mean atomic properties", } } +SCENARIO("Ideal gas density energy from prssure temperature", + "[IdealGas][DensityEnergyFromPressureTemperature]") { + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + GIVEN("An ideal gas") { + EOS host_eos = IdealGas(gm1, Cv); + auto device_eos = host_eos.GetOnDevice(); + WHEN("We compute density and energy from pressure and temperature") { + constexpr int N = 100; + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 1, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real rho = i; + Real T = 100.0 * i; + nw += !CheckRhoSieFromPT(device_eos, rho, T); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + } +} + // A non-standard pattern where we do a reduction class CheckPofRE { public: @@ -283,5 +306,18 @@ SCENARIO("Ideal electron gas", "[IdealGas][IdealEelctrons]") { REQUIRE(nwrong == 0); } } + + WHEN("We compute density and energy from pressure and temperature") { + constexpr int N = 100; + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 1, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real ll[1] = {static_cast(i)}; + nw += !CheckRhoSieFromPT(eos, rho, T, ll); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } } } diff --git a/test/test_eos_mgusup.cpp b/test/test_eos_mgusup.cpp index df35715eac5..162709cc0cb 100644 --- a/test/test_eos_mgusup.cpp +++ b/test/test_eos_mgusup.cpp @@ -408,6 +408,17 @@ SCENARIO("MGUsup EOS rho T", "[VectorEOS][MGUsupEOS]") { } } + WHEN("a [rho,sie](P, T) lookup is performed") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, v_density[i], v_temperature[i]); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + portableFor( "entropy", 0, num, PORTABLE_LAMBDA(const int i) { v_entropy[i] = diff --git a/test/test_eos_noble_abel.cpp b/test/test_eos_noble_abel.cpp index eba52ef49bb..c7f3d75b4b0 100644 --- a/test/test_eos_noble_abel.cpp +++ b/test/test_eos_noble_abel.cpp @@ -125,6 +125,18 @@ SCENARIO("NobleAbel1", "[NobleAbel][NobleAbel1]") { } } + WHEN("A [rho,e](P, T) loopup is performed") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += + !CheckRhoSieFromPT(eos, 7.2347087589269467e-03, 4.0000000000000000e+03); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + WHEN("A B_S(rho, e) lookup is performed") { eos.BulkModulusFromDensityInternalEnergy(v_density, v_energy, v_bulkmodulus, num); #ifdef PORTABILITY_STRATEGY_KOKKOS diff --git a/test/test_eos_powermg.cpp b/test/test_eos_powermg.cpp index d6565a11105..2c0e7e7eace 100644 --- a/test/test_eos_powermg.cpp +++ b/test/test_eos_powermg.cpp @@ -271,6 +271,17 @@ SCENARIO("PowerMG EOS rho sie", "[VectorEOS][PowerMGEOS]") { } } + WHEN("a [rho,sie](P, T) lookup is performed") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, 1, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, 7.285, 7.78960970661031411e+02); + }, + nwrong); + THEN("There are no errors") { REQUIRE(nwrong == 0); } + } + WHEN("A S(rho, e) lookup is performed") { portableFor( "Test entropy", 0, num, PORTABLE_LAMBDA(const int i) { diff --git a/test/test_eos_stellar_collapse.cpp b/test/test_eos_stellar_collapse.cpp index 53524019bdf..7087e4398cb 100644 --- a/test/test_eos_stellar_collapse.cpp +++ b/test/test_eos_stellar_collapse.cpp @@ -278,6 +278,9 @@ SCENARIO("Stellar Collapse EOS", "[StellarCollapse]") { if (!isClose(b1, b2)) { nwrong_d() += 1; } + if (!CheckRhoSieFromPT(sc_d, R, T, lambda)) { + nwrong_d() += 1; + } }); #ifdef PORTABILITY_STRATEGY_KOKKOS Kokkos::deep_copy(nwrong_h, nwrong_d); diff --git a/test/test_eos_stiff.cpp b/test/test_eos_stiff.cpp index a34be2a2949..3f741daf4b5 100644 --- a/test/test_eos_stiff.cpp +++ b/test/test_eos_stiff.cpp @@ -583,6 +583,16 @@ SCENARIO("Test Stiff Gas Entropy Calls", "[StiffGas][StiffGas5]") { array_compare(num, density, energy, h_entropy, entropy_true, "Density", "Energy"); } + THEN("[rho, e](P, T) also works") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, v_density[i], v_local_temp[i]); + }, + nwrong); + AND_THEN("There are no errors") { REQUIRE(nwrong == 0); } + } } } } From 608a772a9a6adbb19549fdfebbf9104425c0e4d7 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:02:59 -0700 Subject: [PATCH 04/45] swap default to findRoot because it respects bounds --- singularity-eos/eos/eos_base.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 5cd7dbf3fbe..27bd4c14f22 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -848,7 +848,7 @@ class EosBase { constexpr Real MAXFAC = 1e8; constexpr Real EPS = 10 * robust::EPS(); constexpr Real MINR_DEFAULT = 1e-4; - using RootFinding1D::regula_falsi; + using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; CRTP copy = *(static_cast(this)); auto PofRT = [&](const Real r) { @@ -863,7 +863,7 @@ class EosBase { // JMM: This can't be zero, in case MinimumDensity is zero Real rhomin = std::max(MINR_DEFAULT, copy.MinimumDensity()); Real rhomax = MAXFAC * rhomin; - auto status = regula_falsi(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); + auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); if (status != Status::SUCCESS) { PORTABLE_THROW_OR_ABORT( "DensityEnergyFromPressureTemperature failed to find root\n"); From 0bdead94d8b46e7bc3f3f2ded0308322afc8bb35 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:03:22 -0700 Subject: [PATCH 05/45] tests complete for everything except eospac and spiner --- test/test_eos_vinet.cpp | 10 ++++++++++ test/test_eos_zsplit.cpp | 5 +++++ 2 files changed, 15 insertions(+) diff --git a/test/test_eos_vinet.cpp b/test/test_eos_vinet.cpp index 0c7f15d9960..916f41aabbc 100644 --- a/test/test_eos_vinet.cpp +++ b/test/test_eos_vinet.cpp @@ -346,6 +346,16 @@ SCENARIO("Vinet EOS rho T", "[VectorEOS][VinetEOS]") { array_compare(num, density, temperature, h_pressure, pressure_true, "Density", "Temperature"); } + THEN("[rho, e](P, T) also works") { + int nwrong = 0; + portableReduce( + "Check density energy from pressure temperature", 0, num, + PORTABLE_LAMBDA(const int i, int &nw) { + nw += !CheckRhoSieFromPT(eos, v_density[i], v_temperature[i]); + }, + nwrong); + AND_THEN("There are no errors") { REQUIRE(nwrong == 0); } + } } portableFor( diff --git a/test/test_eos_zsplit.cpp b/test/test_eos_zsplit.cpp index aeaff8c18c6..1fd39e571a8 100644 --- a/test/test_eos_zsplit.cpp +++ b/test/test_eos_zsplit.cpp @@ -88,6 +88,11 @@ SCENARIO("ZSplit of Ideal Gas", "[ZSplit][IdealGas][IdealElectrons]") { printf("Bad sie Sum! %d %.14e %.14e %.14e\n", i, et, ei, ee); nw += 1; } + + if (Z > 0) { + nw += !CheckRhoSieFromPT(eos_zi, rho, temp, lambda); + nw += !CheckRhoSieFromPT(eos_ze, rho, temp, lambda); + } }, nwrong); THEN("Pe + Pi = Pt") { REQUIRE(nwrong == 0); } From 0a4b61ade0fadeb760a7209e8c8da33ab1f4319e Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:13:07 -0700 Subject: [PATCH 06/45] update zsplit --- test/test_eos_zsplit.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/test_eos_zsplit.cpp b/test/test_eos_zsplit.cpp index 1fd39e571a8..f4eecd35f64 100644 --- a/test/test_eos_zsplit.cpp +++ b/test/test_eos_zsplit.cpp @@ -124,6 +124,9 @@ SCENARIO("ZSplit of Ideal Gas", "[ZSplit][IdealGas][IdealElectrons]") { printf("Bad sie ionized! %d %.14e %.14e\n", i, e_fi, e_ie); nw += 1; } + + nw += !CheckRhoSieFromPT(eos_ze, rho, T, lambda); + nw += !CheckRhoSieFromPT(eos_ie, rho, T, lambda); }, nwrong); THEN("Pe_fi == PE_ie") { REQUIRE(nwrong == 0); } From f1d1fc4d5d609b3ca25e106c263af2812623642c Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:21:37 -0700 Subject: [PATCH 07/45] changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7e75ccc12dc..cffd00fc4f3 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -6,6 +6,7 @@ - [[PR444]](https://github.com/lanl/singularity-eos/pull/444) Add Z split modifier and electron ideal gas EOS ### Fixed (Repair bugs, etc) +- [[PR449]](https://github.com/lanl/singularity-eos/pull/449) Ensure that DensityEnergyFromPressureTemperature works for all equations of state and is properly tested - [[PR439]](https://github.com/lanl/singularity-eos/pull/439) Add mean atomic mass and number to EOS API - [[PR437]](https://github.com/lanl/singularity-eos/pull/437) Fix segfault on HIP, clean up warnings, add strict sanitizer test From c5f75a5da5e2e6a493a6cf9b9d11ade603d24555 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:24:55 -0700 Subject: [PATCH 08/45] missing auto --- singularity-eos/eos/eos_helmholtz.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index 7b0f9de4d52..d4a1af1a971 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/singularity-eos/eos/eos_helmholtz.hpp @@ -693,7 +693,9 @@ class Helmholtz : public EosBase { using RootFinding1D::regula_falsi; using RootFinding1D::Status; PORTABLE_REQUIRE(temp > = 0, "Non-negative temperature required"); - PofRT = [&](const Real r) { return PressureFromDensityTemperature(r, temp, lambda); }; + auto PofRT = [&](const Real r) { + return PressureFromDensityTemperature(r, temp, lambda); + }; auto status = regula_falsi(PofRT, press, 1e7, electrons_.rhoMin(), electrons_.rhoMax(), 1.0e-8, 1.0e-8, rho); if (status != Status::SUCCESS) { From 838b7208779a1ecfe066ec6e9f05c57d6aff1cef Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:30:30 -0700 Subject: [PATCH 09/45] it works on my machine... --- .github/workflows/tests_minimal.yml | 1 + singularity-eos/eos/eos_base.hpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/.github/workflows/tests_minimal.yml b/.github/workflows/tests_minimal.yml index 781802b3e9f..c72aaa8d6eb 100644 --- a/.github/workflows/tests_minimal.yml +++ b/.github/workflows/tests_minimal.yml @@ -37,3 +37,4 @@ jobs: make -j4 make install make test + ctest --output-on-failure --rerun-failed diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 27bd4c14f22..6d2da9c9d1f 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -846,7 +846,7 @@ class EosBase { // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's // overwrite. constexpr Real MAXFAC = 1e8; - constexpr Real EPS = 10 * robust::EPS(); + constexpr Real EPS = robust::EPS(); constexpr Real MINR_DEFAULT = 1e-4; using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; From 365dd65f88d1bea6ee085c311fd8164387d2ed94 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:34:42 -0700 Subject: [PATCH 10/45] output on failure please --- .github/workflows/tests_minimal.yml | 3 +-- .github/workflows/tests_minimal_kokkos.yml | 2 +- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/tests_minimal.yml b/.github/workflows/tests_minimal.yml index c72aaa8d6eb..52f2e70e361 100644 --- a/.github/workflows/tests_minimal.yml +++ b/.github/workflows/tests_minimal.yml @@ -36,5 +36,4 @@ jobs: .. make -j4 make install - make test - ctest --output-on-failure --rerun-failed + ctest --output-on-failure diff --git a/.github/workflows/tests_minimal_kokkos.yml b/.github/workflows/tests_minimal_kokkos.yml index b9666138452..5038f44b4b1 100644 --- a/.github/workflows/tests_minimal_kokkos.yml +++ b/.github/workflows/tests_minimal_kokkos.yml @@ -36,4 +36,4 @@ jobs: .. make -j4 make install - make test + ctest --output-on-failure From 376223f1e87c144139b1654d4c3e2372d2373bde Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:49:37 -0700 Subject: [PATCH 11/45] just change tests to output-on-failure... and also loosen EPS in case thats the issue --- 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 6d2da9c9d1f..dd268d2e68a 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -846,7 +846,7 @@ class EosBase { // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's // overwrite. constexpr Real MAXFAC = 1e8; - constexpr Real EPS = robust::EPS(); + constexpr Real EPS = 100 * robust::EPS(); constexpr Real MINR_DEFAULT = 1e-4; using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; From 9b083fa81505e3de0a5e897ddbc88e6f879d4d07 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 16:53:05 -0700 Subject: [PATCH 12/45] lets try this... --- 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 dd268d2e68a..73070681e67 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -846,7 +846,7 @@ class EosBase { // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's // overwrite. constexpr Real MAXFAC = 1e8; - constexpr Real EPS = 100 * robust::EPS(); + constexpr Real EPS = std::max(1e-8, 100 * robust::EPS()); constexpr Real MINR_DEFAULT = 1e-4; using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; From c3182ff963d09fcf502792b26dc884341d138bdf Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 17:04:05 -0700 Subject: [PATCH 13/45] make default implementation more robust maybe I hope --- singularity-eos/eos/eos_base.hpp | 15 +++++++-------- singularity-eos/eos/eos_vinet.hpp | 4 ++++ 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 73070681e67..7953b037d4d 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -854,19 +854,18 @@ class EosBase { auto PofRT = [&](const Real r) { return copy.PressureFromDensityTemperature(r, temp, lambda); }; - Real rhoguess = rho; // use input density - if (rhoguess < 0) { - Real tguess, sieguess, pguess, cvguess, bmodguess, dpde, dvdt; - copy.ValuesAtReferenceState(rhoguess, tguess, sieguess, pguess, cvguess, bmodguess, - dpde, dvdt); - } // JMM: This can't be zero, in case MinimumDensity is zero Real rhomin = std::max(MINR_DEFAULT, copy.MinimumDensity()); Real rhomax = MAXFAC * rhomin; + Real rhoguess = rho; // use input density + if ((rhoguess < rhomin) || (rhoguess > rhomax)) { + rhoguess = 0.5 * (rhomin + rhomax); + } auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); + // JMM: This needs to not fail and instead return something sane if (status != Status::SUCCESS) { - PORTABLE_THROW_OR_ABORT( - "DensityEnergyFromPressureTemperature failed to find root\n"); + PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n"); + rho = rhoguess; } sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda); return; diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index c942735d1bc..047b79f1393 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -131,6 +131,10 @@ class Vinet : public EosBase { ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const; + + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; } + // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here SG_ADD_BASE_CLASS_USINGS(Vinet) From 1034b38cc37d9bed1638d64fdf35549c8fb1d61d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 17:07:31 -0700 Subject: [PATCH 14/45] put default rho gues close to STP --- singularity-eos/eos/eos_base.hpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 7953b037d4d..11959bbde42 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -848,6 +848,7 @@ class EosBase { constexpr Real MAXFAC = 1e8; constexpr Real EPS = std::max(1e-8, 100 * robust::EPS()); constexpr Real MINR_DEFAULT = 1e-4; + constexpr Real DEFAULT_RHO_GUESS = 4; using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; CRTP copy = *(static_cast(this)); @@ -859,7 +860,11 @@ class EosBase { Real rhomax = MAXFAC * rhomin; Real rhoguess = rho; // use input density if ((rhoguess < rhomin) || (rhoguess > rhomax)) { - rhoguess = 0.5 * (rhomin + rhomax); + if ((rhomin < DEFAULT_RHO_GUESS) && (rhomax > DEFAULT_RHOG_UESS)) { + rhoguess = DEFAULT_RHO_GUESS; + } else { + rhoguess = 0.5 * (rhomin + rhomax); + } } auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); // JMM: This needs to not fail and instead return something sane From d6e80eee989077ad93373f6bf7a4ee43304b76b4 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 17:08:03 -0700 Subject: [PATCH 15/45] put default rho gues close to STP --- 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 11959bbde42..6e1ed9727ea 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -860,7 +860,7 @@ class EosBase { Real rhomax = MAXFAC * rhomin; Real rhoguess = rho; // use input density if ((rhoguess < rhomin) || (rhoguess > rhomax)) { - if ((rhomin < DEFAULT_RHO_GUESS) && (rhomax > DEFAULT_RHOG_UESS)) { + if ((rhomin < DEFAULT_RHO_GUESS) && (rhomax > DEFAULT_RHO_GUESS)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); From 5f7930785b62bc90320f7952be4310a89fa4b362 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 17:13:44 -0700 Subject: [PATCH 16/45] re-tighten bounds after initial guess problem resolved --- 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 6e1ed9727ea..91c3b1ad2dd 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -846,7 +846,7 @@ class EosBase { // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's // overwrite. constexpr Real MAXFAC = 1e8; - constexpr Real EPS = std::max(1e-8, 100 * robust::EPS()); + constexpr Real EPS = 10 * robust::EPS(); constexpr Real MINR_DEFAULT = 1e-4; constexpr Real DEFAULT_RHO_GUESS = 4; using RootFinding1D::findRoot; // more robust but slower. Better default. From 0e61f9871daf2f0f35c70b0a92036b82294c5291 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Tue, 24 Dec 2024 17:32:23 -0700 Subject: [PATCH 17/45] oops missing boolean not --- test/test_eos_helmholtz.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_eos_helmholtz.cpp b/test/test_eos_helmholtz.cpp index 5a200825ea4..037835cacaf 100644 --- a/test/test_eos_helmholtz.cpp +++ b/test/test_eos_helmholtz.cpp @@ -132,7 +132,7 @@ SCENARIO("Helmholtz equation of state - Table interpolation (tgiven)", "[Helmhol if (!isClose(gruen, gruen_ref[k], 1e-6)) nwrong += 1; // RhoSie of PT - nwrong += CheckRhoSieFromPT(eos, rho_in[i], temp_in[i], lambda); + nwrong += !CheckRhoSieFromPT(eos, rho_in[i], temp_in[i], lambda); // Deserialized EOS ein = eos_2.InternalEnergyFromDensityTemperature(rho_in[i], temp_in[j], From cc8db8fdcb4ddf47c8a48e3faf2b0527ca369496 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 10:25:08 -0700 Subject: [PATCH 18/45] add several comments based on Jeff's feedback --- singularity-eos/eos/eos_helmholtz.hpp | 3 ++- singularity-eos/eos/eos_mgusup.hpp | 4 +++- singularity-eos/eos/eos_stellar_collapse.hpp | 1 + singularity-eos/eos/eos_vinet.hpp | 2 ++ 4 files changed, 8 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index d4a1af1a971..444fea08f38 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/singularity-eos/eos/eos_helmholtz.hpp @@ -676,7 +676,8 @@ class Helmholtz : public EosBase { ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const { - // TODO(JMM): Is this right??? + // JMM: Conditions for an oxygen burning shell in a stellar + // core. Not sure if that's the best choice. rho = 1e7; temp = 1.5e9; FillEos(rho, temp, sie, press, cv, bmod, diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index 72dfcec805f..39890c29e74 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -126,6 +126,8 @@ class MGUsup : public EosBase { const unsigned long output, Indexer_t &&lambda = static_cast(nullptr)) const; + // Since reference isotherm scales as _rho0/rho, EOS diverges at + // zero density. This bounds that to some degree. PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return 10 * robust::EPS(); } @@ -135,7 +137,7 @@ class MGUsup : public EosBase { return 0.99 * robust::ratio(_s * _rho0, _s - 1); } else { // for s <= 1, no maximum, but we need to pick something. return 1e3 * _rho0; - } + } // note that s < 0 implies unphysical shock derivative. } template diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 1ff556fadeb..325d69e9ca5 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -664,6 +664,7 @@ PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperat if ((lrguess < lRhoMin_) || (lrguess > lRhoMax_)) { lrguess = lRho_(rhoNormal_); } + // Since EOS is tabulated in logrho, use that for root find. auto lPofRT = [&](Real lR) { return lP_.interpToReal(Ye, lT, lR); }; auto status = regula_falsi(lPofRT, lP, lrguess, lRhoMin_, lRhoMax_, ROOT_THRESH, ROOT_THRESH, lrguess); diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 047b79f1393..63ef97bec56 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -132,6 +132,8 @@ class Vinet : public EosBase { Real &bmod, Real &dpde, Real &dvdt, Indexer_t &&lambda = static_cast(nullptr)) const; + // EOS depends on (rho0/rho)^(1/3) so this is a reasonable bound + // before the EOS will become ill-behaved. PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; } From 58bc9c90fe1b13b45259ac7bd24664cc944b8659 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 10:36:42 -0700 Subject: [PATCH 19/45] try 1e-8 MINR_DEFAULT --- singularity-eos/eos/eos_base.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 91c3b1ad2dd..6165ffaa020 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -845,9 +845,10 @@ class EosBase { Indexer_t &&lambda, Real &rho, Real &sie) const { // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's // overwrite. - constexpr Real MAXFAC = 1e8; + constexpr Real MAXFAC = 1e12; // For MINR_DEFAULT, produces max + // density of 1e4 constexpr Real EPS = 10 * robust::EPS(); - constexpr Real MINR_DEFAULT = 1e-4; + constexpr Real MINR_DEFAULT = 1e-8; constexpr Real DEFAULT_RHO_GUESS = 4; using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; From 99aa6cdc74f551eda7ce8bc7a1399740b7214fa9 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 11:50:02 -0700 Subject: [PATCH 20/45] Make base class use MaximumDensity --- singularity-eos/eos/eos_base.hpp | 36 ++++++++++++------- singularity-eos/eos/eos_helmholtz.hpp | 6 ++++ singularity-eos/eos/eos_mgusup.hpp | 28 --------------- singularity-eos/eos/eos_spiner.hpp | 2 ++ singularity-eos/eos/eos_stellar_collapse.hpp | 1 + .../eos/modifiers/eos_unitsystem.hpp | 3 ++ singularity-eos/eos/modifiers/ramps_eos.hpp | 10 ++++++ .../eos/modifiers/relativistic_eos.hpp | 3 ++ singularity-eos/eos/modifiers/scaled_eos.hpp | 3 ++ singularity-eos/eos/modifiers/shifted_eos.hpp | 3 ++ singularity-eos/eos/modifiers/zsplit_eos.hpp | 3 ++ 11 files changed, 58 insertions(+), 40 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 6165ffaa020..0a2a3ab8d1b 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -16,6 +16,7 @@ #define _SINGULARITY_EOS_EOS_EOS_BASE_ #include +#include #include #include @@ -783,6 +784,15 @@ class EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return 0; } + // Report maximum value of density. Default is unbounded. + // JMM: Should we use actual infinity, the largest real, or just a + // big number? For comparisons, actual infinity is better. It also + // has the advantage of being projective with modifiers that modify + // the max. On the other hand, it's more fraught if someone tries to + // put it into a formula without guarding against it. + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { return 1e100; } + PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } @@ -843,31 +853,33 @@ class EosBase { PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { - // TODO(JMM): A lot hardcoded in here... Hopefully relevent EOS's - // overwrite. - constexpr Real MAXFAC = 1e12; // For MINR_DEFAULT, produces max - // density of 1e4 - constexpr Real EPS = 10 * robust::EPS(); - constexpr Real MINR_DEFAULT = 1e-8; - constexpr Real DEFAULT_RHO_GUESS = 4; using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; + + // Pressure is not monotone in density at low densities, due to a + // phase transition, which can prevent convergence. We want to + // approach tension from above, not below. Choose close to, but + // above, normal density for a metal like copper. + constexpr Real DEFAULT_RHO_GUESS = 12; + CRTP copy = *(static_cast(this)); + // P(rho) not monotone. When relevant, bound rhopmin. + Real rhomin = std::max(copy.rhoPmin(temp), copy.MinimumDensity()); + Real rhomax = copy.MaximumDensity(); auto PofRT = [&](const Real r) { return copy.PressureFromDensityTemperature(r, temp, lambda); }; - // JMM: This can't be zero, in case MinimumDensity is zero - Real rhomin = std::max(MINR_DEFAULT, copy.MinimumDensity()); - Real rhomax = MAXFAC * rhomin; Real rhoguess = rho; // use input density if ((rhoguess < rhomin) || (rhoguess > rhomax)) { - if ((rhomin < DEFAULT_RHO_GUESS) && (rhomax > DEFAULT_RHO_GUESS)) { + if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); } } - auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, EPS, EPS, rho); + printf("%.14e %.14e %.14e\n", rhomin, rhomax, rhoguess); + auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), + robust::EPS(), rho); // JMM: This needs to not fail and instead return something sane if (status != Status::SUCCESS) { PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n"); diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp index 444fea08f38..512f41ff31d 100644 --- a/singularity-eos/eos/eos_helmholtz.hpp +++ b/singularity-eos/eos/eos_helmholtz.hpp @@ -264,6 +264,12 @@ class HelmElectrons { std::size_t numTemp() const { return NTEMP; } PORTABLE_FORCEINLINE_FUNCTION std::size_t numRho() const { return NRHO; } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return rhoMin(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumTemperature() const { return TMin_; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { return rhoMax(); } private: inline void InitDataFile_(const std::string &filename); diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index 39890c29e74..d9b9ea47c57 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -140,34 +140,6 @@ class MGUsup : public EosBase { } // note that s < 0 implies unphysical shock derivative. } - template - PORTABLE_INLINE_FUNCTION void - DensityEnergyFromPressureTemperature(const Real press, const Real temp, - Indexer_t &&lambda, Real &rho, Real &sie) const { - using RootFinding1D::findRoot; - using RootFinding1D::Status; - // JMM: diverges for rho -> 0 - // and for s > 1 and rho -> rho0 - Real rhomin = MinimumDensity(); - Real rhomax = MaximumDensity(); - auto PofRT = [&](const Real r) { - return PressureFromDensityTemperature(r, temp, lambda); - }; - Real rhoguess = rho; // use input density - if ((rhoguess < rhomin) || (rhoguess > rhomax)) { - rhoguess = 0.5 * (rhomin + rhomax); - } - // JMM: regula_falsi does not respect bounds - auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), - robust::EPS(), rho); - if (status != Status::SUCCESS) { - PORTABLE_THROW_OR_ABORT( - "DensityEnergyFromPressureTemperature failed to find root\n"); - } - sie = InternalEnergyFromDensityTemperature(rho, temp, lambda); - return; - } - template PORTABLE_INLINE_FUNCTION void ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index 7a9bb2e2853..ae779f3efc4 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -204,6 +204,7 @@ class SpinerEOSDependsRhoT : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return T_(lTMin_); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } PORTABLE_INLINE_FUNCTION @@ -472,6 +473,7 @@ class SpinerEOSDependsRhoSie : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index 325d69e9ca5..d7c5ecb3190 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -237,6 +237,7 @@ class StellarCollapse : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } inline RootFinding1D::Status rootStatus() const { return status_; } diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 2e7519f7509..c2f03b9032b 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -248,6 +248,9 @@ class UnitSystem : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return inv_temp_unit_ * t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return inv_rho_unit_ * t_.MaximumDensity(); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index fb072f7d7d6..fa9485d07af 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -243,6 +243,16 @@ class BilinearRampEOS : public EosBase> { return; } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { + return t_.MinimumDensity(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { + return t_.MinimumTemperature(); + } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } + template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( const Real rho, const Real temperature, diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index f497e359246..153572ed111 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -155,6 +155,9 @@ class RelativisticEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } static constexpr unsigned long PreferredInput() { return T::PreferredInput(); } diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index 60d522026fe..11c5db28b57 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -341,6 +341,9 @@ class ScaledEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return inv_scale_ * t_.MaximumDensity(); + } PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { return inv_scale_ * t_.MeanAtomicMass(); } diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index f42003a4121..68107908a9f 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -351,6 +351,9 @@ class ShiftedEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index cef562a34f8..4cf4f1eea12 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -239,6 +239,9 @@ class ZSplit : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return t_.MinimumTemperature(); } + PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { + return t_.MaximumDensity(); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( From 1e72ff51c68ece2fc7309d09d0b73be4121a9329 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 11:50:10 -0700 Subject: [PATCH 21/45] update docs --- doc/sphinx/src/using-eos.rst | 20 +++++++++++++++++--- 1 file changed, 17 insertions(+), 3 deletions(-) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index e9bed5ddfb8..237b13583c8 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1143,16 +1143,30 @@ quantities as outputs. Methods Used for Mixed Cell Closures -------------------------------------- -Several methods were developed in support of mixed cell closures. In particular: +Several methods were developed in support of mixed cell closures. In particular the function .. cpp:function:: Real MinimumDensity() const; -and +the function .. cpp:function:: Real MinimumTemperature() const; +and the function + +.. cpp:function:: Real MaximumDensity() const; + provide bounds for valid inputs into a table, which can be used by a -root finder to meaningful bound the root search. Similarly, +root finder to meaningful bound the root search. + +.. warning:: + + For unbounded equations of state, ``MinimumDensity`` and + ``MinimumTemperature`` will return zero, while ``MaximumDensity`` + will return a very large finite number. Which number you get, + however, is not guaranteed. You may wish to apply more sensible + bounds in your own code. + +Similarly, .. cpp:function:: Real RhoPmin(const Real temp) const; From 26dd8f5f6539f0b260f85314304844552b0822fd Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 11:55:53 -0700 Subject: [PATCH 22/45] typo --- 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 0a2a3ab8d1b..74fc286085c 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -864,7 +864,7 @@ class EosBase { CRTP copy = *(static_cast(this)); // P(rho) not monotone. When relevant, bound rhopmin. - Real rhomin = std::max(copy.rhoPmin(temp), copy.MinimumDensity()); + Real rhomin = std::max(copy.RhoPmin(temp), copy.MinimumDensity()); Real rhomax = copy.MaximumDensity(); auto PofRT = [&](const Real r) { return copy.PressureFromDensityTemperature(r, temp, lambda); From be2f5b9c7c01cbf162b1e2ea5d5bf369a03bb687 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:08:54 -0700 Subject: [PATCH 23/45] ok... lets try this --- singularity-eos/eos/eos_base.hpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 74fc286085c..600cace9cf4 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -856,10 +856,10 @@ class EosBase { using RootFinding1D::findRoot; // more robust but slower. Better default. using RootFinding1D::Status; - // Pressure is not monotone in density at low densities, due to a - // phase transition, which can prevent convergence. We want to - // approach tension from above, not below. Choose close to, but - // above, normal density for a metal like copper. + // Pressure is not monotone in density at low densities, which can + // prevent convergence. We want to approach tension from above, + // not below. Choose close to, but above, normal density for a + // metal like copper. constexpr Real DEFAULT_RHO_GUESS = 12; CRTP copy = *(static_cast(this)); @@ -870,14 +870,13 @@ class EosBase { return copy.PressureFromDensityTemperature(r, temp, lambda); }; Real rhoguess = rho; // use input density - if ((rhoguess < rhomin) || (rhoguess > rhomax)) { - if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { + if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { + if ((rhomin <= DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS <= rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); } } - printf("%.14e %.14e %.14e\n", rhomin, rhomax, rhoguess); auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), robust::EPS(), rho); // JMM: This needs to not fail and instead return something sane From 44e738b3a54175fc92ac145d2baef8dd95560280 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:25:45 -0700 Subject: [PATCH 24/45] make tests more effectively check for accurate pressure --- singularity-eos/eos/eos_base.hpp | 4 ++-- test/eos_unit_test_helpers.hpp | 15 +++++++++------ 2 files changed, 11 insertions(+), 8 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 600cace9cf4..5f3a88f8fb4 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -870,8 +870,8 @@ class EosBase { return copy.PressureFromDensityTemperature(r, temp, lambda); }; Real rhoguess = rho; // use input density - if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { - if ((rhomin <= DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS <= rhomax)) { + if ((rhoguess < rhomin) || (rhoguess > rhomax)) { + if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 2e84dd792f3..7abf5ead7d1 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -153,11 +153,13 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Real rtest, etest; eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest); Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); - bool results_good = (isClose(rho, rtest, bmod * 1e-8) || isClose(rho, rtest, 1e-8)) && - (isClose(sie, etest, bmod * 1e-8) || isClose(sie, etest, 1e-8)); + Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); + Real residual = P_test - P; + Real frac_residual = residual / P; + bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) + // This is as good as it will get sometimes. + || (std::abs(frac_residual) < 10*robust::EPS()); if (!results_good) { - Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); - Real residual = P_test - P; printf("RhoSie of PT failure!\n" "\trho_true = %.14e\n" "\tsie_true = %.14e\n" @@ -166,8 +168,9 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, "\trho = %.14e\n" "\tsie = %.14e\n" "\tP_test = %.14e\n" - "\tresidual = %.14e\n", - rho, sie, P, T, rtest, etest, P_test, residual); + "\tresidual = %.14e\n" + "\fracres = %.14e\n", + rho, sie, P, T, rtest, etest, P_test, residual, frac_residual); } return results_good; } From fed072a6371ca7509f8d560b4e7117f7d4538e95 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:27:25 -0700 Subject: [PATCH 25/45] residual --- test/eos_unit_test_helpers.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 7abf5ead7d1..68bae3703ce 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -155,7 +155,7 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); Real residual = P_test - P; - Real frac_residual = residual / P; + Real frac_residual = robust::ratio(residual, P); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) // This is as good as it will get sometimes. || (std::abs(frac_residual) < 10*robust::EPS()); From 8cc2508b0a39314158c45a2cd633272492b6254c Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:27:43 -0700 Subject: [PATCH 26/45] formatting --- test/eos_unit_test_helpers.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 68bae3703ce..d54e6200016 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -157,8 +157,8 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Real residual = P_test - P; Real frac_residual = robust::ratio(residual, P); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) - // This is as good as it will get sometimes. - || (std::abs(frac_residual) < 10*robust::EPS()); + // This is as good as it will get sometimes. + || (std::abs(frac_residual) < 10 * robust::EPS()); if (!results_good) { printf("RhoSie of PT failure!\n" "\trho_true = %.14e\n" From 6dc8cae4befee42ac75da1e85ef6fda63e2784d9 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:30:35 -0700 Subject: [PATCH 27/45] namespace --- test/eos_unit_test_helpers.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index d54e6200016..de89e950e52 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -155,10 +155,10 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); Real residual = P_test - P; - Real frac_residual = robust::ratio(residual, P); + Real frac_residual = singularity::robust::ratio(residual, P); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) // This is as good as it will get sometimes. - || (std::abs(frac_residual) < 10 * robust::EPS()); + || (std::abs(frac_residual) < 10 * singularity::robust::EPS()); if (!results_good) { printf("RhoSie of PT failure!\n" "\trho_true = %.14e\n" From bc5f10eac206973cb4d53576f0d6dd690a9d9ba8 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:35:15 -0700 Subject: [PATCH 28/45] tweak how to check this residual to work when P = 0 --- test/eos_unit_test_helpers.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index de89e950e52..887c7cab8f8 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -155,7 +155,7 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); Real residual = P_test - P; - Real frac_residual = singularity::robust::ratio(residual, P); + Real frac_residual = std::min(singularity::robust::ratio(residual, P), residual); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) // This is as good as it will get sometimes. || (std::abs(frac_residual) < 10 * singularity::robust::EPS()); @@ -169,7 +169,7 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, "\tsie = %.14e\n" "\tP_test = %.14e\n" "\tresidual = %.14e\n" - "\fracres = %.14e\n", + "\tfracres = %.14e\n", rho, sie, P, T, rtest, etest, P_test, residual, frac_residual); } return results_good; From 813e7730e6e901c24f8900f7afdc0a0608e03256 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 12:35:53 -0700 Subject: [PATCH 29/45] oops need to do magnitudes --- test/eos_unit_test_helpers.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 887c7cab8f8..55f13b9b7c6 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -155,7 +155,8 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); Real residual = P_test - P; - Real frac_residual = std::min(singularity::robust::ratio(residual, P), residual); + Real frac_residual = + std::min(std::abs(singularity::robust::ratio(residual, P)), std::abs(residual)); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) // This is as good as it will get sometimes. || (std::abs(frac_residual) < 10 * singularity::robust::EPS()); From e0e587c6849e407f4f73a1c3788b78cc30c2c8b3 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 13:44:42 -0700 Subject: [PATCH 30/45] come up with sensible bounds for power series EOSs. --- singularity-eos/eos/eos_base.hpp | 9 +++++---- singularity-eos/eos/eos_powermg.hpp | 19 +++++++++++++++++++ singularity-eos/eos/eos_vinet.hpp | 9 +++++++++ test/eos_unit_test_helpers.hpp | 2 +- 4 files changed, 34 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 5f3a88f8fb4..c646f69c923 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -877,12 +877,13 @@ class EosBase { rhoguess = 0.5 * (rhomin + rhomax); } } - auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), - robust::EPS(), rho); - // JMM: This needs to not fail and instead return something sane + auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, 10 * robust::EPS(), + 10 * robust::EPS(), rho); + // JMM: This needs to not fail and instead return something sane. + // If root find failed to converge, density will at least be + // within bounds. if (status != Status::SUCCESS) { PORTABLE_WARN("DensityEnergyFromPressureTemperature failed to find root\n"); - rho = rhoguess; } sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda); return; diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index ed89d662ec1..2f8ca0a7660 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -157,6 +157,25 @@ class PowerMG : public EosBase { } printf("\n\n"); } + + // In principle, PowerMG should be valid for all densities. In + // practice, however, since the EOS is parametrized in terms of the + // compression, eta = 1 - rho0/rho, things will fail when rho is + // much less than rho0. Worse, since this is a power law in + // compression, it potentially depends on eta^M where M is the + // maximum index in the power series. In other words, the lower + // bound given by the M^th root of the smallest double times rho0, and + // the upper bound is given by the mth root of the maximum number, + // divided by rho0. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { + return std::pow(std::numeric_limits::min(), 1. / _M) * _rho0; + } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { + return robust::ratio(std::pow(std::numeric_limits::max(), 1. / _M), _rho0); + } + inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } static std::string EosPyType() { return EosType(); } diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 63ef97bec56..0450216bbd2 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -137,6 +137,15 @@ class Vinet : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; } + // In principle, density should be unbounded from above In + // practice, however, since since this is a power law in + // (rho0/rho)^(1/3), it may fail when rho is sufficiently large. + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { + Real maxind = VinetInternalParametersSize; + return robust::ratio(std::pow(std::numeric_limits::max(), 1. / maxind), _rho0); + } + // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here SG_ADD_BASE_CLASS_USINGS(Vinet) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 55f13b9b7c6..dd78417ffe9 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -159,7 +159,7 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, std::min(std::abs(singularity::robust::ratio(residual, P)), std::abs(residual)); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) // This is as good as it will get sometimes. - || (std::abs(frac_residual) < 10 * singularity::robust::EPS()); + || (std::abs(frac_residual) <= 10 * singularity::robust::EPS()); if (!results_good) { printf("RhoSie of PT failure!\n" "\trho_true = %.14e\n" From cf286de1257b56e0be74496a9fd57824efe92f80 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 13:59:46 -0700 Subject: [PATCH 31/45] I give up. 1e-4 it is --- singularity-eos/eos/eos_powermg.hpp | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index 2f8ca0a7660..1c8a80ae6ff 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -163,13 +163,11 @@ class PowerMG : public EosBase { // compression, eta = 1 - rho0/rho, things will fail when rho is // much less than rho0. Worse, since this is a power law in // compression, it potentially depends on eta^M where M is the - // maximum index in the power series. In other words, the lower - // bound given by the M^th root of the smallest double times rho0, and - // the upper bound is given by the mth root of the maximum number, - // divided by rho0. + // maximum index in the power series. PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { - return std::pow(std::numeric_limits::min(), 1. / _M) * _rho0; + // JMM: For whatever reason this seems a good heuristic. + return 1e-4 * _rho0; } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { From 2cdcc083f00ef9c30ed6dbb2aee193c4dd3bb907 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:04:22 -0700 Subject: [PATCH 32/45] comment --- singularity-eos/eos/eos_powermg.hpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index 1c8a80ae6ff..2a75e639778 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -166,7 +166,10 @@ class PowerMG : public EosBase { // maximum index in the power series. PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { - // JMM: For whatever reason this seems a good heuristic. + // JMM: I think formally this should be the mth root of machine + // epsilon times rho0, but for m = 20 or something that's of order + // 1. Things seem reasonably well behaved with this bound. They do + // NOT seem well behaved for, e.g., 10*rho0*machine epsilon. return 1e-4 * _rho0; } PORTABLE_FORCEINLINE_FUNCTION From 1d1873f81908acef9c0e43b0c0abc23793349b6f Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:05:11 -0700 Subject: [PATCH 33/45] remove unused variable --- test/eos_unit_test_helpers.hpp | 1 - 1 file changed, 1 deletion(-) diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index dd78417ffe9..39a51873b59 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -152,7 +152,6 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T, lambda); Real rtest, etest; eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest); - Real bmod = eos.BulkModulusFromDensityTemperature(rho, T, lambda); Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); Real residual = P_test - P; Real frac_residual = From ceec265e1860ca9a5c6018a04c95d937860ccb8b Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:08:13 -0700 Subject: [PATCH 34/45] tighter tolerances --- singularity-eos/eos/eos_base.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index c646f69c923..88124677e32 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -877,8 +877,9 @@ class EosBase { rhoguess = 0.5 * (rhomin + rhomax); } } - auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, 10 * robust::EPS(), - 10 * robust::EPS(), rho); + // JMM: Demand as much tolerance as we can, but don't reset rho below. + auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), + robust::EPS(), rho); // JMM: This needs to not fail and instead return something sane. // If root find failed to converge, density will at least be // within bounds. From 8424efc79edf9c63a75743d375b9cba68ccec368 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:14:34 -0700 Subject: [PATCH 35/45] relax tolerance for chekc --- singularity-eos/eos/eos_base.hpp | 1 - test/eos_unit_test_helpers.hpp | 2 +- 2 files changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 88124677e32..d9fe29d20f1 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -877,7 +877,6 @@ class EosBase { rhoguess = 0.5 * (rhomin + rhomax); } } - // JMM: Demand as much tolerance as we can, but don't reset rho below. auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), robust::EPS(), rho); // JMM: This needs to not fail and instead return something sane. diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 39a51873b59..80db84884ba 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -158,7 +158,7 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, std::min(std::abs(singularity::robust::ratio(residual, P)), std::abs(residual)); bool results_good = (isClose(rho, rtest, 1e-8) && isClose(sie, etest, 1e-8)) // This is as good as it will get sometimes. - || (std::abs(frac_residual) <= 10 * singularity::robust::EPS()); + || (std::abs(frac_residual) <= 1e-8); if (!results_good) { printf("RhoSie of PT failure!\n" "\trho_true = %.14e\n" From 7d9c90a01bed64e08e5ad2a8a13c494f788e6ed9 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:42:16 -0700 Subject: [PATCH 36/45] remove the maximum densities... we're never going to hit those bounds anyway and I think they might be interacting badly with GPUs --- singularity-eos/eos/eos_base.hpp | 4 ++-- singularity-eos/eos/eos_mgusup.hpp | 7 ++++--- singularity-eos/eos/eos_powermg.hpp | 4 ---- singularity-eos/eos/eos_vinet.hpp | 9 --------- 4 files changed, 6 insertions(+), 18 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index d9fe29d20f1..64cd955eb20 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -869,8 +869,8 @@ class EosBase { auto PofRT = [&](const Real r) { return copy.PressureFromDensityTemperature(r, temp, lambda); }; - Real rhoguess = rho; // use input density - if ((rhoguess < rhomin) || (rhoguess > rhomax)) { + Real rhoguess = rho; // use input density + if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { // avoid edge effects if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index d9b9ea47c57..c54c2ad9c20 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -126,10 +126,10 @@ class MGUsup : public EosBase { const unsigned long output, Indexer_t &&lambda = static_cast(nullptr)) const; - // Since reference isotherm scales as _rho0/rho, EOS diverges at - // zero density. This bounds that to some degree. + // Since reference isotherm scales as eta = 1 -_rho0/rho, EOS + // diverges for rho << rho0. eta^2 is used, so... PORTABLE_FORCEINLINE_FUNCTION - Real MinimumDensity() const { return 10 * robust::EPS(); } + Real MinimumDensity() const { return _RHOMINFAC * _rho0; } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { @@ -170,6 +170,7 @@ class MGUsup : public EosBase { private: static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + const Real RHOMINFAC = std::sqrt(robust::EPS()); Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; MeanAtomicProperties _AZbar; }; diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index 2a75e639778..616a31f2950 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -172,10 +172,6 @@ class PowerMG : public EosBase { // NOT seem well behaved for, e.g., 10*rho0*machine epsilon. return 1e-4 * _rho0; } - PORTABLE_FORCEINLINE_FUNCTION - Real MaximumDensity() const { - return robust::ratio(std::pow(std::numeric_limits::max(), 1. / _M), _rho0); - } inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 0450216bbd2..63ef97bec56 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -137,15 +137,6 @@ class Vinet : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; } - // In principle, density should be unbounded from above In - // practice, however, since since this is a power law in - // (rho0/rho)^(1/3), it may fail when rho is sufficiently large. - PORTABLE_FORCEINLINE_FUNCTION - Real MaximumDensity() const { - Real maxind = VinetInternalParametersSize; - return robust::ratio(std::pow(std::numeric_limits::max(), 1. / maxind), _rho0); - } - // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here SG_ADD_BASE_CLASS_USINGS(Vinet) From ee5e3748a25cfae2774e6b15a8e081134aace77d Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:43:08 -0700 Subject: [PATCH 37/45] also add sanity check for rhomax vs rhomin --- singularity-eos/eos/eos_base.hpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 64cd955eb20..950a5d1e63f 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -863,9 +863,12 @@ class EosBase { constexpr Real DEFAULT_RHO_GUESS = 12; CRTP copy = *(static_cast(this)); + // P(rho) not monotone. When relevant, bound rhopmin. Real rhomin = std::max(copy.RhoPmin(temp), copy.MinimumDensity()); Real rhomax = copy.MaximumDensity(); + PORTABLE_REQUIRE(rhomax > rhomin, "max bound > min bound"); + auto PofRT = [&](const Real r) { return copy.PressureFromDensityTemperature(r, temp, lambda); }; From 38bba460710872b55a35c3a4db1cc86f5bb547ac Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:44:48 -0700 Subject: [PATCH 38/45] missed leading underscore --- singularity-eos/eos/eos_mgusup.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index c54c2ad9c20..96cc84be520 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -170,7 +170,7 @@ class MGUsup : public EosBase { private: static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; - const Real RHOMINFAC = std::sqrt(robust::EPS()); + const Real _RHOMINFAC = std::sqrt(robust::EPS()); Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; MeanAtomicProperties _AZbar; }; From a74bd269e5aa39704b758922157a1189d1c3f833 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 14:48:33 -0700 Subject: [PATCH 39/45] stupid C++ deleting copy constructor for const correctness --- singularity-eos/eos/eos_base.hpp | 1 + singularity-eos/eos/eos_mgusup.hpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 950a5d1e63f..d662d8049b1 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -854,6 +854,7 @@ class EosBase { DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { using RootFinding1D::findRoot; // more robust but slower. Better default. + using RootFinding1D::regula_falsi; using RootFinding1D::Status; // Pressure is not monotone in density at low densities, which can diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index 96cc84be520..a44baddb48c 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -170,7 +170,7 @@ class MGUsup : public EosBase { private: static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; - const Real _RHOMINFAC = std::sqrt(robust::EPS()); + Real _RHOMINFAC = std::sqrt(robust::EPS()); Real _rho0, _T0, _Cs, _s, _G0, _Cv0, _E0, _S0; MeanAtomicProperties _AZbar; }; From 5f8e8e0e5732fab2db886733a36a2626e97ac668 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 16:43:52 -0700 Subject: [PATCH 40/45] choose worse, but more generic guess... does that help? --- singularity-eos/eos/eos_base.hpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index d662d8049b1..723cbef4610 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -854,7 +854,6 @@ class EosBase { DensityEnergyFromPressureTemperature(const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { using RootFinding1D::findRoot; // more robust but slower. Better default. - using RootFinding1D::regula_falsi; using RootFinding1D::Status; // Pressure is not monotone in density at low densities, which can @@ -875,11 +874,14 @@ class EosBase { }; Real rhoguess = rho; // use input density if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { // avoid edge effects + /* if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); } + */ + rhoguess = 0.5 * (rhomin + rhomax); } auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), robust::EPS(), rho); From 74ee37788e13d3177f73542c7774de4e18162449 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 16:55:04 -0700 Subject: [PATCH 41/45] is it because rho is uninitialized? --- singularity-eos/eos/eos_base.hpp | 3 --- test/eos_unit_test_helpers.hpp | 3 ++- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 723cbef4610..950a5d1e63f 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -874,14 +874,11 @@ class EosBase { }; Real rhoguess = rho; // use input density if ((rhoguess <= rhomin) || (rhoguess >= rhomax)) { // avoid edge effects - /* if ((rhomin < DEFAULT_RHO_GUESS) && (DEFAULT_RHO_GUESS < rhomax)) { rhoguess = DEFAULT_RHO_GUESS; } else { rhoguess = 0.5 * (rhomin + rhomax); } - */ - rhoguess = 0.5 * (rhomin + rhomax); } auto status = findRoot(PofRT, press, rhoguess, rhomin, rhomax, robust::EPS(), robust::EPS(), rho); diff --git a/test/eos_unit_test_helpers.hpp b/test/eos_unit_test_helpers.hpp index 80db84884ba..3b73e574815 100644 --- a/test/eos_unit_test_helpers.hpp +++ b/test/eos_unit_test_helpers.hpp @@ -150,7 +150,8 @@ CheckRhoSieFromPT(EOS eos, Real rho, Real T, Indexer_t &&lambda = static_cast(nullptr)) { const Real P = eos.PressureFromDensityTemperature(rho, T, lambda); const Real sie = eos.InternalEnergyFromDensityTemperature(rho, T, lambda); - Real rtest, etest; + Real rtest = 12; // set these to something + Real etest = 1; eos.DensityEnergyFromPressureTemperature(P, T, lambda, rtest, etest); Real P_test = eos.PressureFromDensityTemperature(rtest, T, lambda); Real residual = P_test - P; From 9e4f844bd4d438e88647301a2bbe895c47e60df3 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Fri, 3 Jan 2025 17:31:00 -0700 Subject: [PATCH 42/45] update docs --- doc/sphinx/src/using-eos.rst | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 237b13583c8..144148e05cd 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1172,3 +1172,35 @@ Similarly, returns the density at which pressure is minimized for a given temperature. This is again useful for root finds. + +The function + +.. code-block:: cpp + + template + void DensityEnergyFromPressureTemperature(const Real press, const Real temp, + Indexer_t &&lambda, Real &rho, Real &sie) const; + +is designed for working in Pressure-Temperature space. Given a +pressure ``press`` and temperature ``temp``, it sets a density ``rho`` +and specific internal energy ``sie``. + +.. note:: + + Note that ``lambda`` must be passed in, whether or not a given + equation of state requires one. You may pass in ``nullptr`` safely, + however. + +Typically this operation requires a root find. You may pass in an +initial guess for the density ``rho`` in-place and most EOS models +will use it. + +.. warning:: + + Pressure is not necessarily monotone in density and it may be double + valued. Thus you are not guaranteed to find the correct root and the + value of your initial guess may determine correctness. The fact that + ``rho`` may be used as an initial guess means you **must** pass in + an initialized variable, even if it is zero-initialized. Do not pass + uninitialized memory into this function! + From 1cd2f9d107bea5568a4f8732666702e87697477e Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 8 Jan 2025 14:58:21 -0700 Subject: [PATCH 43/45] add introspection for pressure --- doc/sphinx/src/using-eos.rst | 13 ++++- eospac-wrapper/eospac_wrapper.cpp | 9 ++-- eospac-wrapper/eospac_wrapper.hpp | 1 + singularity-eos/eos/eos_base.hpp | 10 ++++ singularity-eos/eos/eos_eospac.hpp | 4 ++ singularity-eos/eos/eos_gruneisen.hpp | 11 +++++ singularity-eos/eos/eos_mgusup.hpp | 8 +++ singularity-eos/eos/eos_powermg.hpp | 7 +++ singularity-eos/eos/eos_sap_polynomial.hpp | 8 +++ singularity-eos/eos/eos_spiner.hpp | 49 +++++++++++++++---- singularity-eos/eos/eos_variant.hpp | 17 +++++++ singularity-eos/eos/eos_vinet.hpp | 8 +++ .../eos/modifiers/eos_unitsystem.hpp | 7 +++ singularity-eos/eos/modifiers/ramps_eos.hpp | 7 +++ .../eos/modifiers/relativistic_eos.hpp | 7 +++ singularity-eos/eos/modifiers/scaled_eos.hpp | 7 +++ singularity-eos/eos/modifiers/shifted_eos.hpp | 7 +++ singularity-eos/eos/modifiers/zsplit_eos.hpp | 7 +++ 18 files changed, 173 insertions(+), 14 deletions(-) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 144148e05cd..cdfc1161445 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1171,7 +1171,18 @@ Similarly, .. cpp:function:: Real RhoPmin(const Real temp) const; returns the density at which pressure is minimized for a given -temperature. This is again useful for root finds. +temperature. The function + +.. cpp:function:: Real MinimumPressure() const; + +provides the minimum pressure an equation of state supports, which may +be the most negative tension state. The function + +.. cpp:function:: Real MaximumPressureFromTemperature() const; + +provides a maximum possible pressure an equation of state +supports. (Most models are unbounded in pressure.) This is again +useful for root finds. The function diff --git a/eospac-wrapper/eospac_wrapper.cpp b/eospac-wrapper/eospac_wrapper.cpp index 3a2fdb8ee84..24ec3c61c09 100644 --- a/eospac-wrapper/eospac_wrapper.cpp +++ b/eospac-wrapper/eospac_wrapper.cpp @@ -25,9 +25,9 @@ namespace EospacWrapper { void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) { - constexpr int NT = 2; + constexpr int NT = 3; EOS_INTEGER tableHandle[NT]; - EOS_INTEGER tableType[NT] = {EOS_Info, EOS_Ut_DT}; + EOS_INTEGER tableType[NT] = {EOS_Info, EOS_Ut_DT, EOS_Pt_DT}; EOS_INTEGER commentsHandle[1]; EOS_INTEGER commentsType[1] = {EOS_Comment}; @@ -48,7 +48,8 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) { } } - eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Info", "EOS_Ut_DT"}, eospacWarn); + eosSafeLoad(NT, matid, tableType, tableHandle, {"EOS_Info", "EOS_Ut_DT", "EOS_Pt_DT"}, + eospacWarn); for (int i = 0; i < numInfoTables; i++) { eosSafeTableInfo(&(tableHandle[i]), NI[i], infoItems[i].data(), infoVals[i].data(), @@ -66,6 +67,8 @@ void eosGetMetadata(int matid, SesameMetadata &metadata, Verbosity eospacWarn) { metadata.TMax = temperatureFromSesame(infoVals[1][3]); metadata.sieMin = sieFromSesame(infoVals[1][4]); metadata.sieMax = sieFromSesame(infoVals[1][5]); + metadata.PMin = pressureFromSesame(infoVals[2][4]); + metadata.PMax = pressureFromSesame(infoVals[2][5]); metadata.numRho = static_cast(infoVals[1][6]); metadata.numT = static_cast(infoVals[1][7]); metadata.rhoConversionFactor = infoVals[1][8]; diff --git a/eospac-wrapper/eospac_wrapper.hpp b/eospac-wrapper/eospac_wrapper.hpp index a9c285d6567..169c6b1e459 100644 --- a/eospac-wrapper/eospac_wrapper.hpp +++ b/eospac-wrapper/eospac_wrapper.hpp @@ -59,6 +59,7 @@ class SesameMetadata { double rhoMin, rhoMax; double TMin, TMax; double sieMin, sieMax; + double PMin, PMax; double rhoConversionFactor; double TConversionFactor; double sieConversionFactor; diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 950a5d1e63f..c6e1ae44a88 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -793,6 +793,16 @@ class EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return 1e100; } + // These are for the PT space PTE solver to bound the iterations in + // a safe range. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return 0; } + // Gruneisen EOS's often have a maximum density, which implies a maximum pressure. + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { + return 1e100; + } + PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 763cb536955..0e6ef06e0cd 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -1104,6 +1104,7 @@ class EOSPAC : public EosBase { } PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rho_min_; } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return temp_min_; } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return press_min_; } private: static constexpr const unsigned long _preferred_input = @@ -1128,6 +1129,7 @@ class EOSPAC : public EosBase { Real dvdt_ref_ = 1; Real rho_min_ = 0; Real temp_min_ = 0; + Real press_min_ = 0; // TODO(JMM): Is the fact that EOS_INTEGER isn't a size_t a // problem? Could it ever realistically overflow? EOS_INTEGER shared_size_, packed_size_; @@ -1199,6 +1201,8 @@ inline EOSPAC::EOSPAC(const int matid, bool invert_at_setup, Real insert_data, rho_ref_ = m.normalDensity; rho_min_ = m.rhoMin; temp_min_ = m.TMin; + press_min_ = m.PMin; + // use std::max to hydrogen, in case of bad table AZbar_.Abar = std::max(1.0, m.meanAtomicMass); AZbar_.Zbar = std::max(1.0, m.meanAtomicNumber); diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index a2892e9a2af..ff05389692e 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -170,6 +170,17 @@ class Gruneisen : public EosBase { s1, _C0, _s1, _s2, _s3, _G0, _b, _rho0, _T0, _P0, _Cv, _rho_max); _AZbar.PrintParams(); } + + // Report minimum values of density and temperature + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { return _rho_max; } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return PressureFromDensityInternalEnergy(0, 0); } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature(const Real T) const { + return MaxStableDensityAtTemperature(T); + } + template PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, const Real temp, diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index a44baddb48c..ed43f6271b6 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -139,6 +139,14 @@ class MGUsup : public EosBase { return 1e3 * _rho0; } // note that s < 0 implies unphysical shock derivative. } + // Hugoniot pressure ill defined at reference density. On one side, + // negative. On the other positive. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { + return 1e100; + } template PORTABLE_INLINE_FUNCTION void diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index 616a31f2950..bcfd3553bf2 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -172,6 +172,13 @@ class PowerMG : public EosBase { // NOT seem well behaved for, e.g., 10*rho0*machine epsilon. return 1e-4 * _rho0; } + // Essentially unbounded... I think. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { + return 1e100; + } inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp index 081e434d0fc..d37db9bc85a 100644 --- a/singularity-eos/eos/eos_sap_polynomial.hpp +++ b/singularity-eos/eos/eos_sap_polynomial.hpp @@ -162,6 +162,14 @@ class SAP_Polynomial : public EosBase { return rho / _rho0 - 1; } + // Essentially unbounded... I think. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { + return 1e100; + } + template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, diff --git a/singularity-eos/eos/eos_spiner.hpp b/singularity-eos/eos/eos_spiner.hpp index ae779f3efc4..6af2e644d1b 100644 --- a/singularity-eos/eos/eos_spiner.hpp +++ b/singularity-eos/eos/eos_spiner.hpp @@ -205,6 +205,9 @@ class SpinerEOSDependsRhoT : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return T_(lTMin_); } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return PMin_; } + PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } PORTABLE_INLINE_FUNCTION @@ -306,6 +309,7 @@ class SpinerEOSDependsRhoT : public EosBase { Real lRhoMin_, lRhoMax_, rhoMax_; Real lRhoMinSearch_; Real lTMin_, lTMax_, TMax_; + Real PMin_; Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoOffset_, lTOffset_; // offsets must be non-negative @@ -474,6 +478,12 @@ class SpinerEOSDependsRhoSie : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return rhoMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MinimumTemperature() const { return TMin(); } PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return rhoMax(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return PMin_; } + PORTABLE_INLINE_FUNCTION + Real RhoPmin(const Real temp) const { + return rho_at_pmin_.interpToReal(toLog_(temp, lTOffset_)); + } PORTABLE_INLINE_FUNCTION int nlambda() const noexcept { return _n_lambda; } @@ -517,22 +527,24 @@ class SpinerEOSDependsRhoSie : public EosBase { DataBox sie_; // depends on (rho,T) DataBox T_; // depends on (rho, sie) + DataBox rho_at_pmin_; SP5Tables dependsRhoT_; SP5Tables dependsRhoSie_; - int numRho_; + int numRho_, numT_; Real rhoNormal_, TNormal_, sieNormal_, PNormal_; Real CvNormal_, bModNormal_, dPdENormal_, dVdTNormal_; Real lRhoMin_, lRhoMax_, rhoMax_; DataBox PlRhoMax_, dPdRhoMax_; - + Real PMin_; Real lRhoOffset_, lTOffset_, lEOffset_; // offsets must be non-negative #define DBLIST \ - &sie_, &T_, &(dependsRhoT_.P), &(dependsRhoT_.bMod), &(dependsRhoT_.dPdRho), \ - &(dependsRhoT_.dPdE), &(dependsRhoT_.dTdRho), &(dependsRhoT_.dTdE), \ - &(dependsRhoT_.dEdRho), &(dependsRhoSie_.P), &(dependsRhoSie_.bMod), \ - &(dependsRhoSie_.dPdRho), &(dependsRhoSie_.dPdE), &(dependsRhoSie_.dTdRho), \ - &(dependsRhoSie_.dTdE), &(dependsRhoSie_.dEdRho), &PlRhoMax_, &dPdRhoMax_ + &sie_, &T_, &rho_at_pmin_, &(dependsRhoT_.P), &(dependsRhoT_.bMod), \ + &(dependsRhoT_.dPdRho), &(dependsRhoT_.dPdE), &(dependsRhoT_.dTdRho), \ + &(dependsRhoT_.dTdE), &(dependsRhoT_.dEdRho), &(dependsRhoSie_.P), \ + &(dependsRhoSie_.bMod), &(dependsRhoSie_.dPdRho), &(dependsRhoSie_.dPdE), \ + &(dependsRhoSie_.dTdRho), &(dependsRhoSie_.dTdE), &(dependsRhoSie_.dEdRho), \ + &PlRhoMax_, &dPdRhoMax_ std::vector GetDataBoxPointers_() const { return std::vector{DBLIST}; } @@ -773,11 +785,11 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, rho_at_pmin_.resize(numT_); rho_at_pmin_.setRange(0, P_.range(0)); for (int i = 0; i < numT_; i++) { - Real pmin = std::numeric_limits::max(); + PMin_ = std::numeric_limits::max(); int jmax = -1; for (int j = 0; j < numRho_; j++) { - if (P_(j, i) < pmin) { - pmin = P_(j, i); + if (P_(j, i) < PMin_) { + PMin_ = P_(j, i); jmax = j; } } @@ -1597,6 +1609,7 @@ herr_t SpinerEOSDependsRhoSie::loadDataboxes_(const std::string &matid_str, hid_ // Metadata for root finding extrapolation numRho_ = sie_.dim(2); + numT_ = sie_.dim(1); lRhoMin_ = sie_.range(1).min(); lRhoMax_ = sie_.range(1).max(); rhoMax_ = fromLog_(lRhoMax_, lRhoOffset_); @@ -1605,6 +1618,22 @@ herr_t SpinerEOSDependsRhoSie::loadDataboxes_(const std::string &matid_str, hid_ PlRhoMax_ = dependsRhoT_.P.slice(numRho_ - 1); dPdRhoMax_ = dependsRhoT_.dPdRho.slice(numRho_ - 1); + // fill in minimum pressure as a function of temperature + rho_at_pmin_.resize(numT_); + rho_at_pmin_.setRange(0, sie_.range(0)); + for (int i = 0; i < numT_; i++) { + PMin_ = std::numeric_limits::max(); + int jmax = -1; + for (int j = 0; j < numRho_; j++) { + if (dependsRhoT_.P(j, i) < PMin_) { + PMin_ = dependsRhoT_.P(j, i); + jmax = j; + } + } + if (jmax < 0) printf("Failed to find minimum pressure.\n"); + rho_at_pmin_(i) = fromLog_(dependsRhoT_.P.range(1).x(jmax), lRhoOffset_); + } + // reference state Real lRhoNormal = toLog_(rhoNormal_, lRhoOffset_); // if rho normal not on the table, set it to the middle diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 31dcb7867a3..88c57c7a770 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -333,6 +333,23 @@ class Variant { return mpark::visit([](const auto &eos) { return eos.MinimumTemperature(); }, eos_); } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumDensity() const { + return mpark::visit([](const auto &eos) { return eos.MaximumDensity(); }, eos_); + } + + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { + return mpark::visit([](const auto &eos) { return eos.MinimumPressure(); }, eos_); + } + + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature(const Real temp) const { + return mpark::visit( + [&temp](const auto &eos) { return eos.MaximumPressureFromTemperature(temp); }, + eos_); + } + // Atomic mass/atomic number functions PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index 63ef97bec56..ed14dff02cb 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -137,6 +137,14 @@ class Vinet : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumDensity() const { return std::cbrt(10 * robust::EPS()) * _rho0; } + // Essentially unbounded... I think. + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumPressure() const { return -1e100; } + PORTABLE_FORCEINLINE_FUNCTION + Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { + return 1e100; + } + // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here SG_ADD_BASE_CLASS_USINGS(Vinet) diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index c2f03b9032b..106440959dd 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -251,6 +251,13 @@ class UnitSystem : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return inv_rho_unit_ * t_.MaximumDensity(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return inv_press_unit_ * t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real + MaximumPressureFromTemperature(const Real temp) const { + return inv_press_unit_ * t_.MaximumPressureFromTemperature(temp_unit_ * temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index fa9485d07af..405991cf6eb 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -252,6 +252,13 @@ class BilinearRampEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return t_.MaximumDensity(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real + MaximumPressureFromTemperature(const Real temp) const { + return t_.MaximumPressureFromTemperature(temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index 153572ed111..e0110486016 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -158,6 +158,13 @@ class RelativisticEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return t_.MaximumDensity(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real + MaximumPressureFromTemperature(const Real temp) const { + return t_.MaximumPressureFromTemperature(temp); + } static constexpr unsigned long PreferredInput() { return T::PreferredInput(); } diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index 11c5db28b57..f4a40986e30 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -344,6 +344,13 @@ class ScaledEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return inv_scale_ * t_.MaximumDensity(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real + MaximumPressureFromTemperature(const Real temp) const { + return t_.MaximumPressureFromTemperature(temp); + } PORTABLE_INLINE_FUNCTION Real MeanAtomicMass() const { return inv_scale_ * t_.MeanAtomicMass(); } diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 68107908a9f..41c5cc5be14 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -354,6 +354,13 @@ class ShiftedEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return t_.MaximumDensity(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real + MaximumPressureFromTemperature(const Real temp) const { + return t_.MaximumPressureFromTemperature(temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index 4cf4f1eea12..96ba71badaa 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -242,6 +242,13 @@ class ZSplit : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MaximumDensity() const { return t_.MaximumDensity(); } + PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { + return t_.MinimumPressure(); + } + PORTABLE_FORCEINLINE_FUNCTION Real + MaximumPressureFromTemperature(const Real temp) const { + return t_.MaximumPressureFromTemperature(temp); + } template PORTABLE_INLINE_FUNCTION Real MeanAtomicMassFromDensityTemperature( From 6967e5516fb3c9207811087bfe229179ec813fd5 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 9 Jan 2025 08:58:30 -0700 Subject: [PATCH 44/45] add default nullptr argument --- doc/sphinx/src/using-eos.rst | 9 ++------- singularity-eos/eos/eos_base.hpp | 8 ++++++++ singularity-eos/eos/eos_variant.hpp | 10 ++++++++++ 3 files changed, 20 insertions(+), 7 deletions(-) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index cdfc1161445..eb794a7f230 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1194,13 +1194,8 @@ The function is designed for working in Pressure-Temperature space. Given a pressure ``press`` and temperature ``temp``, it sets a density ``rho`` -and specific internal energy ``sie``. - -.. note:: - - Note that ``lambda`` must be passed in, whether or not a given - equation of state requires one. You may pass in ``nullptr`` safely, - however. +and specific internal energy ``sie``. The ``lambda`` is optional and +defaults to a ``nullptr``. Typically this operation requires a root find. You may pass in an initial guess for the density ``rho`` in-place and most EOS models diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index c6e1ae44a88..5c41bd4244f 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -901,6 +901,14 @@ class EosBase { sie = copy.InternalEnergyFromDensityTemperature(rho, temp, lambda); return; } + PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, + const Real temp, + Real &rho, + Real &sie) const { + CRTP copy = *(static_cast(this)); + copy.DensityEnergyFromPressureTemperature(press, temp, static_cast(nullptr), + rho, sie); + } // Serialization /* diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 88c57c7a770..2a1c61a84fe 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -317,6 +317,16 @@ class Variant { }, eos_); } + PORTABLE_INLINE_FUNCTION void DensityEnergyFromPressureTemperature(const Real press, + const Real temp, + Real &rho, + Real &sie) const { + return mpark::visit( + [&press, &temp, &rho, &sie](const auto &eos) { + return eos.DensityEnergyFromPressureTemperature(press, temp, rho, sie); + }, + eos_); + } PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { From 144b90b2525b76c1133a0908edb94ef32578c7eb Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 9 Jan 2025 09:44:31 -0700 Subject: [PATCH 45/45] change name --- doc/sphinx/src/using-eos.rst | 8 ++++---- singularity-eos/eos/eos_base.hpp | 4 +--- singularity-eos/eos/eos_gruneisen.hpp | 2 +- singularity-eos/eos/eos_mgusup.hpp | 4 +--- singularity-eos/eos/eos_powermg.hpp | 4 +--- singularity-eos/eos/eos_sap_polynomial.hpp | 4 +--- singularity-eos/eos/eos_variant.hpp | 4 ++-- singularity-eos/eos/eos_vinet.hpp | 4 +--- singularity-eos/eos/modifiers/eos_unitsystem.hpp | 5 ++--- singularity-eos/eos/modifiers/ramps_eos.hpp | 5 ++--- singularity-eos/eos/modifiers/relativistic_eos.hpp | 5 ++--- singularity-eos/eos/modifiers/scaled_eos.hpp | 5 ++--- singularity-eos/eos/modifiers/shifted_eos.hpp | 5 ++--- singularity-eos/eos/modifiers/zsplit_eos.hpp | 5 ++--- 14 files changed, 24 insertions(+), 40 deletions(-) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index eb794a7f230..ae1db924a5a 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1178,11 +1178,11 @@ temperature. The function provides the minimum pressure an equation of state supports, which may be the most negative tension state. The function -.. cpp:function:: Real MaximumPressureFromTemperature() const; +.. cpp:function:: Real MaximumPressureAtTemperature(const Real temp) const; -provides a maximum possible pressure an equation of state -supports. (Most models are unbounded in pressure.) This is again -useful for root finds. +provides a maximum possible pressure an equation of state supports at +a given temperature. (Most models are unbounded in pressure.) This is +again useful for root finds. The function diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 5c41bd4244f..7f685dd88e2 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -799,9 +799,7 @@ class EosBase { Real MinimumPressure() const { return 0; } // Gruneisen EOS's often have a maximum density, which implies a maximum pressure. PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { - return 1e100; - } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return 0.0; } diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp index ff05389692e..d343705b847 100644 --- a/singularity-eos/eos/eos_gruneisen.hpp +++ b/singularity-eos/eos/eos_gruneisen.hpp @@ -177,7 +177,7 @@ class Gruneisen : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return PressureFromDensityInternalEnergy(0, 0); } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature(const Real T) const { + Real MaximumPressureAtTemperature(const Real T) const { return MaxStableDensityAtTemperature(T); } diff --git a/singularity-eos/eos/eos_mgusup.hpp b/singularity-eos/eos/eos_mgusup.hpp index ed43f6271b6..b85cd9175fd 100644 --- a/singularity-eos/eos/eos_mgusup.hpp +++ b/singularity-eos/eos/eos_mgusup.hpp @@ -144,9 +144,7 @@ class MGUsup : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return -1e100; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { - return 1e100; - } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } template PORTABLE_INLINE_FUNCTION void diff --git a/singularity-eos/eos/eos_powermg.hpp b/singularity-eos/eos/eos_powermg.hpp index bcfd3553bf2..76da7607d07 100644 --- a/singularity-eos/eos/eos_powermg.hpp +++ b/singularity-eos/eos/eos_powermg.hpp @@ -176,9 +176,7 @@ class PowerMG : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return -1e100; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { - return 1e100; - } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } inline void Finalize() {} static std::string EosType() { return std::string("PowerMG"); } diff --git a/singularity-eos/eos/eos_sap_polynomial.hpp b/singularity-eos/eos/eos_sap_polynomial.hpp index d37db9bc85a..c7358c48fa5 100644 --- a/singularity-eos/eos/eos_sap_polynomial.hpp +++ b/singularity-eos/eos/eos_sap_polynomial.hpp @@ -166,9 +166,7 @@ class SAP_Polynomial : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return -1e100; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { - return 1e100; - } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } template PORTABLE_INLINE_FUNCTION void diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 2a1c61a84fe..7b48a2e7e04 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -354,9 +354,9 @@ class Variant { } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature(const Real temp) const { + Real MaximumPressureAtTemperature(const Real temp) const { return mpark::visit( - [&temp](const auto &eos) { return eos.MaximumPressureFromTemperature(temp); }, + [&temp](const auto &eos) { return eos.MaximumPressureAtTemperature(temp); }, eos_); } diff --git a/singularity-eos/eos/eos_vinet.hpp b/singularity-eos/eos/eos_vinet.hpp index ed14dff02cb..67d3e722159 100644 --- a/singularity-eos/eos/eos_vinet.hpp +++ b/singularity-eos/eos/eos_vinet.hpp @@ -141,9 +141,7 @@ class Vinet : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return -1e100; } PORTABLE_FORCEINLINE_FUNCTION - Real MaximumPressureFromTemperature([[maybe_unused]] const Real T) const { - return 1e100; - } + Real MaximumPressureAtTemperature([[maybe_unused]] const Real T) const { return 1e100; } // Generic functions provided by the base class. These contain e.g. the vector // overloads that use the scalar versions declared here diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 106440959dd..853f1b2f944 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -254,9 +254,8 @@ class UnitSystem : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return inv_press_unit_ * t_.MinimumPressure(); } - PORTABLE_FORCEINLINE_FUNCTION Real - MaximumPressureFromTemperature(const Real temp) const { - return inv_press_unit_ * t_.MaximumPressureFromTemperature(temp_unit_ * temp); + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return inv_press_unit_ * t_.MaximumPressureAtTemperature(temp_unit_ * temp); } template diff --git a/singularity-eos/eos/modifiers/ramps_eos.hpp b/singularity-eos/eos/modifiers/ramps_eos.hpp index 405991cf6eb..67a5f5185aa 100644 --- a/singularity-eos/eos/modifiers/ramps_eos.hpp +++ b/singularity-eos/eos/modifiers/ramps_eos.hpp @@ -255,9 +255,8 @@ class BilinearRampEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return t_.MinimumPressure(); } - PORTABLE_FORCEINLINE_FUNCTION Real - MaximumPressureFromTemperature(const Real temp) const { - return t_.MaximumPressureFromTemperature(temp); + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); } template diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index e0110486016..d6cbe4d76c7 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -161,9 +161,8 @@ class RelativisticEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return t_.MinimumPressure(); } - PORTABLE_FORCEINLINE_FUNCTION Real - MaximumPressureFromTemperature(const Real temp) const { - return t_.MaximumPressureFromTemperature(temp); + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); } static constexpr unsigned long PreferredInput() { return T::PreferredInput(); } diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index f4a40986e30..d48e6d9504a 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -347,9 +347,8 @@ class ScaledEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return t_.MinimumPressure(); } - PORTABLE_FORCEINLINE_FUNCTION Real - MaximumPressureFromTemperature(const Real temp) const { - return t_.MaximumPressureFromTemperature(temp); + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); } PORTABLE_INLINE_FUNCTION diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index 41c5cc5be14..933f655366b 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -357,9 +357,8 @@ class ShiftedEOS : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return t_.MinimumPressure(); } - PORTABLE_FORCEINLINE_FUNCTION Real - MaximumPressureFromTemperature(const Real temp) const { - return t_.MaximumPressureFromTemperature(temp); + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); } template diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index 96ba71badaa..9b3505a3b71 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -245,9 +245,8 @@ class ZSplit : public EosBase> { PORTABLE_FORCEINLINE_FUNCTION Real MinimumPressure() const { return t_.MinimumPressure(); } - PORTABLE_FORCEINLINE_FUNCTION Real - MaximumPressureFromTemperature(const Real temp) const { - return t_.MaximumPressureFromTemperature(temp); + PORTABLE_FORCEINLINE_FUNCTION Real MaximumPressureAtTemperature(const Real temp) const { + return t_.MaximumPressureAtTemperature(temp); } template