From 388dacb859fc50cbbc103695670924b6aa2ca5ee Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 16:47:09 -0500 Subject: [PATCH 01/11] add InternalEnergyFromDensityPressure call to base class and variant --- singularity-eos/eos/eos_base.hpp | 57 +++++++++++++++++++++++++++ singularity-eos/eos/eos_variant.hpp | 60 +++++++++++++++++++++++++++++ 2 files changed, 117 insertions(+) diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 739dbc62764..3668f5d5cd3 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -82,6 +82,7 @@ char *StrCat(char *destination, const char *source) { using EosBase<__VA_ARGS__>::BulkModulusFromDensityInternalEnergy; \ using EosBase<__VA_ARGS__>::GruneisenParamFromDensityTemperature; \ using EosBase<__VA_ARGS__>::GruneisenParamFromDensityInternalEnergy; \ + using EosBase<__VA_ARGS__>::InternalEnergyFromDensityPressure; \ using EosBase<__VA_ARGS__>::FillEos; \ using EosBase<__VA_ARGS__>::EntropyFromDensityTemperature; \ using EosBase<__VA_ARGS__>::EntropyFromDensityInternalEnergy; \ @@ -928,6 +929,62 @@ class EosBase { rho, sie); } + // JMM: Another set of calls that are often overloaded for special cases + // TODO(JMM): Do we also want TemperatureFromDensityPressure? That's + // the more likely fundamental call, but the less likely "useful" + // call... + // We pass in sie by value to allow for initial guesses + template + PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure( + const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + using RootFinding1D::regula_falsi; + using RootFinding1D::Status; + const CRTP &eos = *(static_cast(this)); + auto f = [&](const Real sie) { + return eos.PressureFromDensityInternalEnergy(rho, sie, lambda); + }; + const Real sie_min = + eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature()); + // temp not bounded. just pick something huge. + const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20); + Real sie_guess = + (((sie_min < sie) && (sie < sie_max)) ? 0.5 * (sie_min + sie_max) : sie_guess); + auto status = regula_falsi(f, P, sie_guess, sie_min, sie_max, robust::EPS(), + robust::EPS(), sie); + if (status == Status::FAIL) { + PORTABLE_WARN("InternalEnergyFromDensityPressure: failed to find root\n"); + } + } + // Classes like EOSPAC probably wish to overload/shadow this version + template + inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, + ConstRealIndexer &&Ps, RealIndexer &&sies, + const int num, + LambdaIndexer &&lambdas) const { + static auto const name = SG_MEMBER_FUNC_NAME(); + static auto const cname = name.c_str(); + CRTP copy = *(static_cast(this)); + portableFor( + cname, 0, num, PORTABLE_LAMBDA(const int i) { + copy.InternalEnergyFromDensityPressure(rhos[i], Ps[i], sies[i], lambdas[i]); + }); + } + template + inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, + ConstRealIndexer &&Ps, RealIndexer &&sies, + Real * /*scratch */, const int num, + LambdaIndexer &&lambdas) const { + return InternalEnergyFromDensityPressure(rhos, Ps, sies, num, lambdas); + } + template + inline void InternalEnergyFromDensityPressure(const Real *rhos, const Real *Ps, + Real *sies, Real * /*scratch */, + const int num, LambdaIndexer &&lambdas, + Transform && = Transform()) const { + return InternalEnergyFromDensityPressure(rhos, Ps, sies, num, lambdas); + } + // Serialization /* The methodology here is there are *three* size methods all EOS's provide: diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index a73ed8e5c44..6bb5ce9f4cc 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -334,6 +334,17 @@ class Variant { eos_); } + template + PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure( + const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + return mpark::visit( + [&rho, &P, &sie, &lambda](const auto &eos) { + return eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda); + }, + eos_); + } + PORTABLE_INLINE_FUNCTION Real RhoPmin(const Real temp) const { return mpark::visit([&temp](const auto &eos) { return eos.RhoPmin(temp); }, eos_); @@ -1188,6 +1199,55 @@ class Variant { eos_); } + template + inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, + ConstRealIndexer &&Ps, RealIndexer &&sies, + const int num) const { + NullIndexer lambdas{}; // Returns null pointer for every index + return InternalEnergyFromDensityPressure( + std::forward(rhos), std::forward(Ps), + std::forward(sies), num, lambdas); + } + + template + inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, + ConstRealIndexer &&Ps, RealIndexer &&sies, + const int num, + LambdaIndexer &&lambdas) const { + return mpark::visit( + [&](const auto &eos) { + return eos.InternalEnergyFromDensityPressure( + std::forward(rhos), std::forward(Ps), + std::forward(sies), num, std::forward(lambdas)); + }, + eos_); + } + + template + inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, + ConstRealIndexer &&Ps, RealIndexer &&sies, + Real *scratch, const int num) const { + NullIndexer lambdas{}; // Returns null pointer for every index + return InternalEnergyFromDensityPressure( + std::forward(rhos), std::forward(Ps), + std::forward(sies), scratch, num, lambdas); + } + + template + inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, + ConstRealIndexer &&Ps, RealIndexer &&sies, + Real *scratch, const int num, + LambdaIndexer &&lambdas) const { + return mpark::visit( + [&rhos, &Ps, &sies, &scratch, &num, &lambdas](const auto &eos) { + return eos.InternalEnergyFromDensityPressure( + std::forward(rhos), std::forward(Ps), + std::forward(sies), scratch, num, + std::forward(lambdas)); + }, + eos_); + } + template inline void FillEos(RealIndexer &&rhos, RealIndexer &&temps, RealIndexer &&energies, RealIndexer &&presses, RealIndexer &&cvs, RealIndexer &&bmods, From e8996710f6d38b87167d4d837924fce778aff202 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 16:49:37 -0500 Subject: [PATCH 02/11] add special case for ideal gas. use root find for gruneisen. tell eospac to use the base class FOR NOW --- singularity-eos/eos/eos_eospac.hpp | 2 ++ singularity-eos/eos/eos_ideal.hpp | 8 ++++++++ 2 files changed, 10 insertions(+) diff --git a/singularity-eos/eos/eos_eospac.hpp b/singularity-eos/eos/eos_eospac.hpp index 032a5334b13..246144083ed 100644 --- a/singularity-eos/eos/eos_eospac.hpp +++ b/singularity-eos/eos/eos_eospac.hpp @@ -494,6 +494,8 @@ class EOSPAC : public EosBase { // TODO(JMM): Add performant Gibbs Free Energy using EosBase::FillEos; + // TODO(JMM): Add a performant internal energy call + using EosBase::InternalEnergyFromDensityPressure; SG_ADD_DEFAULT_MEAN_ATOMIC_FUNCTIONS(AZbar_) diff --git a/singularity-eos/eos/eos_ideal.hpp b/singularity-eos/eos/eos_ideal.hpp index 364bca4aec5..73c6bea53a3 100644 --- a/singularity-eos/eos/eos_ideal.hpp +++ b/singularity-eos/eos/eos_ideal.hpp @@ -186,6 +186,14 @@ class IdealGas : public EosBase { sie = MYMAX(0.0, _Cv * temp); rho = MYMAX(0.0, press / (_gm1 * sie)); } + template + PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure( + const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = static_cast(nullptr)) const { + // P = gm1 rho sie + sie = P / (_gm1 * rho); + } + inline void Finalize() {} static std::string EosType() { return std::string("IdealGas"); } static std::string EosPyType() { return EosType(); } From 58218cc10fdaa77e6309e19e14634874744e929a Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 16:49:57 -0500 Subject: [PATCH 03/11] add tests to gruneisen and ideal gas. tests pass. --- test/test_eos_gruneisen.cpp | 61 +++++++++++++++++++++++++++++++++++++ test/test_eos_ideal.cpp | 49 +++++++++++++++++++++++++++++ 2 files changed, 110 insertions(+) diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 376760c64b1..4bdb9af5248 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -208,6 +208,67 @@ SCENARIO("Gruneisen EOS", "[VectorEOS][GruneisenEOS]") { } } +SCENARIO("Gruneisen energy from pressure and temperature", + "[Gruneisen][VectorEOS][SieFromRhoP]") { + GIVEN("An ideal gas object") { + // Unit conversions + constexpr Real cm = 1.; + constexpr Real us = 1e-06; + constexpr Real Mbcc_per_g = 1e12; + // Gruneisen parameters for copper + constexpr Real C0 = 0.394 * cm / us; + constexpr Real S1 = 1.489; + constexpr Real S2 = 0.; + constexpr Real S3 = 0.; + constexpr Real Gamma0 = 2.02; + constexpr Real b = 0.47; + constexpr Real rho0 = 8.93; + constexpr Real T0 = 298.; + constexpr Real P0 = 0.; + constexpr Real Cv = 0.383e-05 * Mbcc_per_g; + // Create the EOS + EOS eos_host = Gruneisen(C0, S1, S2, S3, Gamma0, b, rho0, T0, P0, Cv); + EOS eos = eos_host.GetOnDevice(); + + WHEN("We compute pressure from energy for several energies") { + constexpr int N = 20; + Real *rhos = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + Real *sies = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + Real *Ps = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + portableFor( + "Set rho, sie, and P", 0, N, PORTABLE_LAMBDA(const int i) { + rhos[i] = 0.5 * (eos.MinimumDensity() + eos.MaximumDensity()); + sies[i] = Cv * (298 + i); + Ps[i] = eos.PressureFromDensityInternalEnergy(rhos[i], sies[i]); + }); + + THEN("We can compute the internal energy from the pressure") { + Real *sies_new = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + eos_host.InternalEnergyFromDensityPressure(rhos, Ps, sies_new, N); + + AND_THEN("They agree with the previous pressures") { + std::size_t nwrong = 0; + portableReduce( + "Compute sie diff", 0, N, + PORTABLE_LAMBDA(const int i, std::size_t &nw) { + nw += !isClose(sies[i], sies_new[i], 1e-12); + }, + nwrong); + REQUIRE(nwrong == 0); + } + + PORTABLE_FREE(sies_new); + } + + PORTABLE_FREE(Ps); + PORTABLE_FREE(sies); + PORTABLE_FREE(rhos); + } + + eos.Finalize(); + } +} + SCENARIO("Aluminum Gruneisen EOS", "[GruneisenEOS]") { GIVEN("Parameters for an aluminum Gruneisen EOS") { // Unit conversions diff --git a/test/test_eos_ideal.cpp b/test/test_eos_ideal.cpp index 6685370744c..4b90f524894 100644 --- a/test/test_eos_ideal.cpp +++ b/test/test_eos_ideal.cpp @@ -89,6 +89,55 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy][GibbsFreeEnergy]") { } } +SCENARIO("Ideal gas energy from pressure and temperature", + "[IdealGas][VectorEOS][SieFromRhoP]") { + GIVEN("An ideal gas object") { + constexpr Real Cv = 5.0; + constexpr Real gm1 = 0.4; + EOS eos_host = IdealGas(gm1, Cv); + EOS eos = eos_host.GetOnDevice(); + + WHEN("We compute pressure from energy for several energies") { + constexpr int N = 20; + constexpr Real sie_min = 0; + constexpr Real dsie = 0.5; + Real *rhos = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + Real *sies = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + Real *Ps = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + portableFor( + "Set rho, sie, and P", 0, N, PORTABLE_LAMBDA(const int i) { + rhos[i] = 1; + sies[i] = sie_min + dsie * i; + Ps[i] = eos.PressureFromDensityInternalEnergy(rhos[i], sies[i]); + }); + + THEN("We can compute the internal energy from the pressure") { + Real *sies_new = (Real *)PORTABLE_MALLOC(sizeof(Real) * N); + eos_host.InternalEnergyFromDensityPressure(rhos, Ps, sies_new, N); + + AND_THEN("They agree with the previous pressures") { + std::size_t nwrong = 0; + portableReduce( + "Compute sie diff", 0, N, + PORTABLE_LAMBDA(const int i, std::size_t &nw) { + nw += !isClose(sies[i], sies_new[i], 1e-12); + }, + nwrong); + REQUIRE(nwrong == 0); + } + + PORTABLE_FREE(sies_new); + } + + PORTABLE_FREE(Ps); + PORTABLE_FREE(sies); + PORTABLE_FREE(rhos); + } + + eos.Finalize(); + } +} + SCENARIO("Ideal gas mean atomic properties", "[IdealGas][MeanAtomicMass][MeanAtomicNumber]") { constexpr Real Cv = 5.0; From fdd5faff9653d605d98715e960822954c0428ba6 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 16:53:36 -0500 Subject: [PATCH 04/11] docs --- doc/sphinx/src/using-eos.rst | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index 326c4173367..f987a615d0e 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1181,6 +1181,18 @@ given density in :math:`g/cm^3` and temperature in Kelvin. returns the unitless Gruneisen parameter given density in :math:`g/cm^3` and specific internal energy in :math:`erg/g`. +.. code-block:: cpp + + template + void InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const; + +returns the specific internal energy in :math:`erg/g` given a density +in :math:`g/cm^3` and a pressure in Barye. This call is often a root +find, and thus an initial guess for specific internal energy may be +passed in by reference. It will be set to the new value by the +function. + .. code-block:: cpp template From bf29d4e0e510151daa85b640210b27eee4007044 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 16:56:55 -0500 Subject: [PATCH 05/11] changelog + CC --- CHANGELOG.md | 1 + test/test_eos_gruneisen.cpp | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 15571cf144e..663d444628a 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,6 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) +- [[PR589]](https://github.com/lanl/singularity-eos/pull/589) PressureFromDensityInternalEnergy - [[PR556]](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant - [[PR564]](https://github.com/lanl/singularity-eos/pull/564) Removed Get() function from IndexableTypes since it could have unexpected consequences when a type wasn't present - [[PR583]](https://github.com/lanl/singularity-eos/pull/583) Added GenericIndexer class to provide more complex array indirection diff --git a/test/test_eos_gruneisen.cpp b/test/test_eos_gruneisen.cpp index 4bdb9af5248..ac1df8f9e6f 100644 --- a/test/test_eos_gruneisen.cpp +++ b/test/test_eos_gruneisen.cpp @@ -1,5 +1,5 @@ //------------------------------------------------------------------------------ -// © 2021-2024. Triad National Security, LLC. All rights reserved. This +// © 2021-2025. Triad National Security, LLC. All rights reserved. This // program was produced under U.S. Government contract 89233218CNA000001 // for Los Alamos National Laboratory (LANL), which is operated by Triad // National Security, LLC for the U.S. Department of Energy/National From fb251abf82a902d368d36a1936dc09927a11351a Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 17:38:53 -0500 Subject: [PATCH 06/11] try to add a scalar function and python vector functionality --- doc/sphinx/src/using-eos.rst | 15 +++++++++---- python/module.hpp | 21 ++++++++++++++++++- singularity-eos/eos/eos_base.hpp | 14 +++++++++++++ .../eos/modifiers/eos_unitsystem.hpp | 8 +++++++ .../eos/modifiers/floored_energy.hpp | 6 ++++++ .../eos/modifiers/relativistic_eos.hpp | 6 ++++++ singularity-eos/eos/modifiers/scaled_eos.hpp | 6 ++++++ singularity-eos/eos/modifiers/shifted_eos.hpp | 8 +++++++ singularity-eos/eos/modifiers/zsplit_eos.hpp | 11 ++++++++++ 9 files changed, 90 insertions(+), 5 deletions(-) diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index f987a615d0e..ed42b204eb1 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/doc/sphinx/src/using-eos.rst @@ -1184,14 +1184,21 @@ returns the unitless Gruneisen parameter given density in .. code-block:: cpp template - void InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Real InternalEnergyFromDensityPressure(const Real rho, const Real P, Indexer_t &&lambda = nullptr) const; returns the specific internal energy in :math:`erg/g` given a density in :math:`g/cm^3` and a pressure in Barye. This call is often a root -find, and thus an initial guess for specific internal energy may be -passed in by reference. It will be set to the new value by the -function. +find, and thus an alternative signature allows an initial guess to be +passed in by reference: + +.. code-block:: cpp + + template + void InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const; + +In this case, the guess will be set to the new value by the function. .. code-block:: cpp diff --git a/python/module.hpp b/python/module.hpp index 5520ddf25cd..0867c783f1a 100644 --- a/python/module.hpp +++ b/python/module.hpp @@ -43,6 +43,16 @@ Real two_params_no_lambda(const T& self, const Real a, const Real b) { return (self.*Func)(a, b, np()); } +template +Real three_params(const T& self, const Real a, const Real b, Real &c, py::array_t lambda) { + return (self.*Func)(a, b, c, lambda.mutable_data()); +} + +template +Real three_params_no_lambda(const T& self, const Real a, const Real b, Real &c) { + return (self.*Func)(a, b, c, np()); +} + class LambdaHelper { py::array_t & lambdas; public: @@ -169,6 +179,7 @@ EOS_VEC_FUNC_TMPL(BulkModulusFromDensityTemperature, rhos, temperatures, bmods) EOS_VEC_FUNC_TMPL(BulkModulusFromDensityInternalEnergy, rhos, sies, bmods) EOS_VEC_FUNC_TMPL(GruneisenParamFromDensityTemperature, rhos, temperatures, gm1s) EOS_VEC_FUNC_TMPL(GruneisenParamFromDensityInternalEnergy, rhos, sies, gm1s) +EOS_VEC_FUNC_TMPL(InternalEnergyFromDensityPressure, rhos, Ps, sies) struct EOSState { Real density; @@ -229,6 +240,7 @@ struct VectorFunctions { .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergy, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("lmbdas")) .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperature, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("lmbdas")) .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergy, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("lmbdas")) + .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressure, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("lmbdas")) .def("TemperatureFromDensityInternalEnergy", &TemperatureFromDensityInternalEnergyNoLambda, py::arg("rhos"), py::arg("sies"), py::arg("temperatures")) .def("InternalEnergyFromDensityTemperature", &InternalEnergyFromDensityTemperatureNoLambda, py::arg("rhos"), py::arg("temperatures"), py::arg("sies")) @@ -240,7 +252,7 @@ struct VectorFunctions { .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyNoLambda, py::arg("rhos"), py::arg("sies"), py::arg("bmods")) .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureNoLambda, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s")) .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyNoLambda, py::arg("rhos"), py::arg("sies"), py::arg("gm1s")) - + .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureNoLambda, py::arg("rhos"), py::arg("Ps"), py::arg("sies")) .def("FillEos", [](const T & self, py::array_t rhos, py::array_t temperatures, py::array_t sies, @@ -300,6 +312,7 @@ struct VectorFunctions { .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyWithScratch, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("scratch"), py::arg("lmbdas")) .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureWithScratch, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("scratch"), py::arg("lmbdas")) .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyWithScratch, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("scratch"), py::arg("lmbdas")) + .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureWithScratch, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("scratch"), py::arg("lmbdas")) .def("TemperatureFromDensityInternalEnergy", &TemperatureFromDensityInternalEnergyNoLambdaWithScratch, py::arg("rhos"), py::arg("sies"), py::arg("temperatures"), py::arg("scratch")) .def("InternalEnergyFromDensityTemperature", &InternalEnergyFromDensityTemperatureNoLambdaWithScratch, py::arg("rhos"), py::arg("temperatures"), py::arg("sies"), py::arg("scratch")) @@ -311,6 +324,7 @@ struct VectorFunctions { .def("BulkModulusFromDensityInternalEnergy", &BulkModulusFromDensityInternalEnergyNoLambdaWithScratch, py::arg("rhos"), py::arg("sies"), py::arg("bmods"), py::arg("scratch")) .def("GruneisenParamFromDensityTemperature", &GruneisenParamFromDensityTemperatureNoLambdaWithScratch, py::arg("rhos"), py::arg("temperatures"), py::arg("gm1s"), py::arg("scratch")) .def("GruneisenParamFromDensityInternalEnergy", &GruneisenParamFromDensityInternalEnergyNoLambdaWithScratch, py::arg("rhos"), py::arg("sies"), py::arg("gm1s"), py::arg("scratch")) + .def("InternalEnergyFromDensityPressure", &InternalEnergyFromDensityPressureNoLambdaWithScratch, py::arg("rhos"), py::arg("Ps"), py::arg("sies"), py::arg("scratch")) .def("FillEos", [](const T & self, py::array_t rhos, py::array_t temperatures, py::array_t sies, @@ -365,6 +379,9 @@ struct VectorFunctions { } }; +// TODO(JMM): I do not understand why my definitions below fail, but +// template substituion fails. Perhaps there are too many overloads of +// this method to disambiguate. template py::class_ eos_class(py::module_ & m, std::string name) { py::class_ cls(m, name.c_str()); @@ -378,6 +395,7 @@ py::class_ eos_class(py::module_ & m, std::string name) { .def("BulkModulusFromDensityInternalEnergy", &two_params, py::arg("rho"), py::arg("sie"), py::arg("lmbda")) .def("GruneisenParamFromDensityTemperature", &two_params, py::arg("rho"), py::arg("temperature"), py::arg("lmbda")) .def("GruneisenParamFromDensityInternalEnergy", &two_params, py::arg("rho"), py::arg("sie"), py::arg("lmbda")) + //.def("InternalEnergyFromDensityPressure", &three_params, py::arg("rho"), py::arg("P"), py::arg("sie"), py::arg("lmbda")) .def("TemperatureFromDensityInternalEnergy", &two_params_no_lambda, py::arg("rho"), py::arg("sie")) .def("InternalEnergyFromDensityTemperature", &two_params_no_lambda, py::arg("rho"), py::arg("temperature")) @@ -389,6 +407,7 @@ py::class_ eos_class(py::module_ & m, std::string name) { .def("BulkModulusFromDensityInternalEnergy", &two_params_no_lambda, py::arg("rho"), py::arg("sie")) .def("GruneisenParamFromDensityTemperature", &two_params_no_lambda, py::arg("rho"), py::arg("temperature")) .def("GruneisenParamFromDensityInternalEnergy", &two_params_no_lambda, py::arg("rho"), py::arg("sie")) + //.def("InternalEnergyFromDensityPressure", &three_params_no_lambda, py::arg("rho"), py::arg("P"), py::arg("sie")) .def("FillEos", [](const T & self, const py::kwargs& kwargs) { unsigned long output = thermalqs::none; diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index 3668f5d5cd3..7a46b65b291 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -956,6 +956,20 @@ class EosBase { PORTABLE_WARN("InternalEnergyFromDensityPressure: failed to find root\n"); } } + // This version foregoes the initial guess. You get what you get. + template + PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityPressure( + const Real rho, const Real P, + Indexer_t &&lambda = static_cast(nullptr)) const { + const CRTP &eos = *(static_cast(this)); + const Real sie_min = + eos.InternalEnergyFromDensityTemperature(rho, eos.MinimumTemperature()); + // temp not bounded. just pick something huge. + const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20); + Real sie = 0.5 * (sie_min + sie_max); + eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda); + return sie; + } // Classes like EOSPAC probably wish to overload/shadow this version template inline void InternalEnergyFromDensityPressure(ConstRealIndexer &&rhos, diff --git a/singularity-eos/eos/modifiers/eos_unitsystem.hpp b/singularity-eos/eos/modifiers/eos_unitsystem.hpp index 98bf594ef84..8acbf19025b 100644 --- a/singularity-eos/eos/modifiers/eos_unitsystem.hpp +++ b/singularity-eos/eos/modifiers/eos_unitsystem.hpp @@ -193,6 +193,14 @@ class UnitSystem : public EosBase> { temp * temp_unit_, lambda); return gm1; } + template + PORTABLE_FUNCTION void + InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const { + sie *= sie_unit_; + t_.InternalEnergyFromDensityPressure(rho_unit_ * rho, press_unit_ * P, sie, lambda); + sie *= inv_sie_unit_; + } template PORTABLE_FUNCTION void diff --git a/singularity-eos/eos/modifiers/floored_energy.hpp b/singularity-eos/eos/modifiers/floored_energy.hpp index b3b926df57c..314cd33a432 100644 --- a/singularity-eos/eos/modifiers/floored_energy.hpp +++ b/singularity-eos/eos/modifiers/floored_energy.hpp @@ -121,6 +121,12 @@ class FlooredEnergy : public EosBase> { return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda); } template + PORTABLE_FUNCTION void + InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const { + t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda); + } + template PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda = nullptr) const { diff --git a/singularity-eos/eos/modifiers/relativistic_eos.hpp b/singularity-eos/eos/modifiers/relativistic_eos.hpp index fafa2f95a62..8863a7295a9 100644 --- a/singularity-eos/eos/modifiers/relativistic_eos.hpp +++ b/singularity-eos/eos/modifiers/relativistic_eos.hpp @@ -140,6 +140,12 @@ class RelativisticEOS : public EosBase> { return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda); } template + PORTABLE_FUNCTION void + InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const { + t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda); + } + template PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda = nullptr) const { diff --git a/singularity-eos/eos/modifiers/scaled_eos.hpp b/singularity-eos/eos/modifiers/scaled_eos.hpp index c7c2292a3e9..63851ded5fb 100644 --- a/singularity-eos/eos/modifiers/scaled_eos.hpp +++ b/singularity-eos/eos/modifiers/scaled_eos.hpp @@ -136,6 +136,12 @@ class ScaledEOS : public EosBase> { return t_.GruneisenParamFromDensityTemperature(scale_ * rho, temperature, lambda); } template + PORTABLE_FUNCTION void + InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const { + t_.InternalEnergyFromDensityPressure(scale_ * rho, P, sie, lambda); + } + template PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda = nullptr) const { diff --git a/singularity-eos/eos/modifiers/shifted_eos.hpp b/singularity-eos/eos/modifiers/shifted_eos.hpp index ff52b2eac26..cd0bad539b6 100644 --- a/singularity-eos/eos/modifiers/shifted_eos.hpp +++ b/singularity-eos/eos/modifiers/shifted_eos.hpp @@ -127,6 +127,14 @@ class ShiftedEOS : public EosBase> { return t_.GruneisenParamFromDensityTemperature(rho, temperature, lambda); } template + PORTABLE_FUNCTION void + InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const { + sie -= shift_; + t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda); + sie += shift_; + } + template PORTABLE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda = nullptr) const { diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp index ce850acd4e4..ee9f260e400 100644 --- a/singularity-eos/eos/modifiers/zsplit_eos.hpp +++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp @@ -169,6 +169,17 @@ class ZSplit : public EosBase> { return t_.GruneisenParamFromDensityInternalEnergy(rho, iscale * sie, lambda); } + template + PORTABLE_FUNCTION void + InternalEnergyFromDensityPressure(const Real rho, const Real P, Real &sie, + Indexer_t &&lambda = nullptr) const { + const Real scale = GetScale_(lambda); + const Real iscale = GetInvScale_(lambda); + sie *= iscale; + t_.InternalEnergyFromDensityPressure(rho, P, sie, lambda); + sie *= scale; + } + template PORTABLE_INLINE_FUNCTION void FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, Real &bmod, From ab16f91fbc2bb9d1e97a763886205cfffc0b887b Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Wed, 26 Nov 2025 17:42:48 -0500 Subject: [PATCH 07/11] oops fix a warning --- 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 7a46b65b291..2e8581454b4 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -949,7 +949,7 @@ class EosBase { // temp not bounded. just pick something huge. const Real sie_max = eos.InternalEnergyFromDensityTemperature(rho, 1e20); Real sie_guess = - (((sie_min < sie) && (sie < sie_max)) ? 0.5 * (sie_min + sie_max) : sie_guess); + (((sie_min < sie) && (sie < sie_max)) ? 0.5 * (sie_min + sie_max) : sie); auto status = regula_falsi(f, P, sie_guess, sie_min, sie_max, robust::EPS(), robust::EPS(), sie); if (status == Status::FAIL) { From a102531aca7a303023932563e3edf59d8040c1c6 Mon Sep 17 00:00:00 2001 From: Richard Berger Date: Thu, 4 Dec 2025 13:49:38 -0700 Subject: [PATCH 08/11] fix python bindings for InternalEnergyFromDensityPressure --- python/module.hpp | 20 ++++++++------------ 1 file changed, 8 insertions(+), 12 deletions(-) diff --git a/python/module.hpp b/python/module.hpp index 0867c783f1a..88f6597b39b 100644 --- a/python/module.hpp +++ b/python/module.hpp @@ -43,16 +43,6 @@ Real two_params_no_lambda(const T& self, const Real a, const Real b) { return (self.*Func)(a, b, np()); } -template -Real three_params(const T& self, const Real a, const Real b, Real &c, py::array_t lambda) { - return (self.*Func)(a, b, c, lambda.mutable_data()); -} - -template -Real three_params_no_lambda(const T& self, const Real a, const Real b, Real &c) { - return (self.*Func)(a, b, c, np()); -} - class LambdaHelper { py::array_t & lambdas; public: @@ -395,7 +385,6 @@ py::class_ eos_class(py::module_ & m, std::string name) { .def("BulkModulusFromDensityInternalEnergy", &two_params, py::arg("rho"), py::arg("sie"), py::arg("lmbda")) .def("GruneisenParamFromDensityTemperature", &two_params, py::arg("rho"), py::arg("temperature"), py::arg("lmbda")) .def("GruneisenParamFromDensityInternalEnergy", &two_params, py::arg("rho"), py::arg("sie"), py::arg("lmbda")) - //.def("InternalEnergyFromDensityPressure", &three_params, py::arg("rho"), py::arg("P"), py::arg("sie"), py::arg("lmbda")) .def("TemperatureFromDensityInternalEnergy", &two_params_no_lambda, py::arg("rho"), py::arg("sie")) .def("InternalEnergyFromDensityTemperature", &two_params_no_lambda, py::arg("rho"), py::arg("temperature")) @@ -407,7 +396,6 @@ py::class_ eos_class(py::module_ & m, std::string name) { .def("BulkModulusFromDensityInternalEnergy", &two_params_no_lambda, py::arg("rho"), py::arg("sie")) .def("GruneisenParamFromDensityTemperature", &two_params_no_lambda, py::arg("rho"), py::arg("temperature")) .def("GruneisenParamFromDensityInternalEnergy", &two_params_no_lambda, py::arg("rho"), py::arg("sie")) - //.def("InternalEnergyFromDensityPressure", &three_params_no_lambda, py::arg("rho"), py::arg("P"), py::arg("sie")) .def("FillEos", [](const T & self, const py::kwargs& kwargs) { unsigned long output = thermalqs::none; @@ -471,6 +459,14 @@ py::class_ eos_class(py::module_ & m, std::string name) { self.DensityEnergyFromPressureTemperature(press, temp, np(), rho, sie); return std::pair(rho, sie); }, py::arg("press"), py::arg("temp")) + .def("InternalEnergyFromDensityPressure", [](const T & self, const Real rho, const Real P, Real sie) { + self.InternalEnergyFromDensityPressure(rho, P, sie, np()); + return sie; + }, py::arg("rho"), py::arg("P"), py::arg("sie")) + .def("InternalEnergyFromDensityPressure", [](const T & self, const Real rho, const Real P, Real sie, py::array_t lambda) { + self.InternalEnergyFromDensityPressure(rho, P, sie, lambda.mutable_data()); + return sie; + }, py::arg("rho"), py::arg("P"), py::arg("sie"), py::arg("lmbda")) .def("Finalize", &T::Finalize) .def_property_readonly_static("EosType", [](py::object) { return T::EosType(); }) .def_static("scratch_size", &T::scratch_size, py::arg("method_name"), py::arg("nelements")) From 0c4cf0ca16ff7949d4521fe2a570468b9c35feab Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 4 Dec 2025 17:32:06 -0500 Subject: [PATCH 09/11] Update CHANGELOG.md --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 663d444628a..538c0ffe653 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -3,7 +3,7 @@ ## Current develop ### Added (new features/APIs/variables/...) -- [[PR589]](https://github.com/lanl/singularity-eos/pull/589) PressureFromDensityInternalEnergy +- [[PR589]](https://github.com/lanl/singularity-eos/pull/589) InternalEnergyFromDensityPressure - [[PR556]](https://github.com/lanl/singularity-eos/pull/556) Add introspection into types available in the variant - [[PR564]](https://github.com/lanl/singularity-eos/pull/564) Removed Get() function from IndexableTypes since it could have unexpected consequences when a type wasn't present - [[PR583]](https://github.com/lanl/singularity-eos/pull/583) Added GenericIndexer class to provide more complex array indirection From ea9d4628e38f1e6bc50b24e24bad63c187c52475 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 18 Dec 2025 17:25:16 -0500 Subject: [PATCH 10/11] Update module.hpp --- python/module.hpp | 3 --- 1 file changed, 3 deletions(-) diff --git a/python/module.hpp b/python/module.hpp index 88f6597b39b..278cff68277 100644 --- a/python/module.hpp +++ b/python/module.hpp @@ -369,9 +369,6 @@ struct VectorFunctions { } }; -// TODO(JMM): I do not understand why my definitions below fail, but -// template substituion fails. Perhaps there are too many overloads of -// this method to disambiguate. template py::class_ eos_class(py::module_ & m, std::string name) { py::class_ cls(m, name.c_str()); From e83ed8737a821e7e9fb5b9c705cd06b62bd02aa0 Mon Sep 17 00:00:00 2001 From: Jonah Miller Date: Thu, 18 Dec 2025 18:07:34 -0500 Subject: [PATCH 11/11] mpark -> portsofcall in new code --- singularity-eos/eos/eos_variant.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 063a64707be..d04b7c3674f 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -340,7 +340,7 @@ class Variant { PORTABLE_INLINE_FUNCTION void InternalEnergyFromDensityPressure( const Real rho, const Real P, Real &sie, Indexer_t &&lambda = static_cast(nullptr)) const { - return mpark::visit( + return PortsOfCall::visit( [&rho, &P, &sie, &lambda](const auto &eos) { return eos.InternalEnergyFromDensityPressure(rho, P, sie, lambda); }, @@ -1220,7 +1220,7 @@ class Variant { ConstRealIndexer &&Ps, RealIndexer &&sies, const int num, LambdaIndexer &&lambdas) const { - return mpark::visit( + return PortsOfCall::visit( [&](const auto &eos) { return eos.InternalEnergyFromDensityPressure( std::forward(rhos), std::forward(Ps), @@ -1244,7 +1244,7 @@ class Variant { ConstRealIndexer &&Ps, RealIndexer &&sies, Real *scratch, const int num, LambdaIndexer &&lambdas) const { - return mpark::visit( + return PortsOfCall::visit( [&rhos, &Ps, &sies, &scratch, &num, &lambdas](const auto &eos) { return eos.InternalEnergyFromDensityPressure( std::forward(rhos), std::forward(Ps),