From fa9d603aae87044fe97776abc2b035653f954b65 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 9 May 2024 19:14:04 -0600 Subject: [PATCH 01/53] Allow E0 to adjust isentrope energy --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87e2ef06f6..d9f2e4fa50 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -313,7 +313,7 @@ class DavisProducts : public EosBase { const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / - std::pow(vvc, _k - 1.0 + _a); + std::pow(vvc, _k - 1.0 + _a) - _E0; } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { const Real vvc = 1 / (rho * _vc); From 96d71fcfc6c9bfe01ff19d8f80d2e50218bffcec Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 9 May 2024 19:16:51 -0600 Subject: [PATCH 02/53] Add place-holder for CJ state calculation --- singularity-eos/eos/eos_davis.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d9f2e4fa50..0ab82cbd53 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -292,6 +292,7 @@ class DavisProducts : public EosBase { inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } static std::string EosPyType() { return EosType(); } + // TODO (JHP): Create helper function to find the CJ state given the reference state private: static constexpr Real onethird = 1.0 / 3.0; From 4defc1fb438b3ebaa40b3ea525e2d68a6ad3b5b7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 9 May 2024 19:46:51 -0600 Subject: [PATCH 03/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 0ab82cbd53..86eedd224b 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -314,7 +314,8 @@ class DavisProducts : public EosBase { const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / - std::pow(vvc, _k - 1.0 + _a) - _E0; + std::pow(vvc, _k - 1.0 + _a) - + _E0; } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { const Real vvc = 1 / (rho * _vc); From fe4a2edf645eef2caa1fa5a0163edd3674391e36 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 9 May 2024 19:47:01 -0600 Subject: [PATCH 04/53] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6800600af4..7a49bee9d5 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface +- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Fix missing utilization of E0 in Davis Products EOS ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls From 4c26d70c0152dcf41502947ac4f287ceca4c10d8 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Thu, 9 May 2024 20:16:55 -0600 Subject: [PATCH 05/53] Update documentation of E0 for Davis Products --- doc/sphinx/src/models.rst | 36 +++++++++++++++++++++++++++++++----- 1 file changed, 31 insertions(+), 5 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 12f555caac..19424ce45c 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -11,6 +11,8 @@ .. _DavisReactants: https://doi.org/10.1016/S0010-2180(99)00112-1 +.. _DavisProducts: https://doi.org/10.1063/1.2035310 + .. _ProbingOffHugoniotStates: https://doi.org/10.1063/1.4939675 .. _WillsThermo: https://www.osti.gov/biblio/1561015 @@ -1194,10 +1196,11 @@ Davis Products EOS .. warning:: Entropy is not yet available for this EOS -The Davis products EOS is created from the reference isentrope passing through -the CJ state of the high explosive along with a constant heat capacity. The -constant heat capacity leads to the energy being a simple funciton of the -temperature deviation from the reference isentrope such that +The `Davis products EOS `_ is created from the reference +isentrope passing through the CJ state of the high explosive along with a +constant heat capacity. The constant heat capacity leads to the energy being a +simple funciton of the temperature deviation from the reference isentrope such +that .. math:: @@ -1228,7 +1231,7 @@ Finally, the pressure, energy, and temperature along the isentrope are given by .. math:: - e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} + e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} - e_0 .. math:: @@ -1254,6 +1257,29 @@ Here, there are four dimensionless parameters that are settable by the user, :math:`e_\mathrm{C}`, :math:`V_\mathrm{C}` and :math:`T_\mathrm{C}` are tuning parameters with units related to their non-subscripted counterparts. +Note that the energy zero (i.e. the reference energy) for the Davis products EOS +is arbitrary. For the isentrope to properly pass through the CJ state of a +reacting material, the energy release of the reaction needs to be accounted for +properly. If done external to the EOS, an energy source term is required in the +Euler equations. However, a common convention is to specify the reactants and +product EOS in a consistent way such that the reference energy corresponds to +the rest state of the material *before* it reacts. + +The energy at the CJ state can be calculated as + +.. math:: + + e_\mathrm{CJ} = \frac{P_0 + P_\mathrm{CJ}}{2(V_0 - V_\mathrm{CJ})}, + +relative to :math:`e = 0` at the reference state of the *reactants*. Therefore +the :math:`e_0` energy offset of the products EOS is given by + +.. math:: + + e_0 = e_S(V_\mathrm{CJ}) - e_\mathrm{CJ}. + +Practically, this means :math:`e_0` should be positive for any energetic material. + The constructor for the Davis Products EOS is .. code-block:: cpp From 02a1e2980f894be8774b1521757ff43e76a730b8 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:19:53 -0600 Subject: [PATCH 06/53] Add initialization loops to the PTE/EOS cache --- singularity-eos/eos/get_sg_eos_p_t.cpp | 4 ++++ singularity-eos/eos/get_sg_eos_rho_e.cpp | 4 ++++ singularity-eos/eos/get_sg_eos_rho_p.cpp | 4 ++++ singularity-eos/eos/get_sg_eos_rho_t.cpp | 4 ++++ 4 files changed, 16 insertions(+) diff --git a/singularity-eos/eos/get_sg_eos_p_t.cpp b/singularity-eos/eos/get_sg_eos_p_t.cpp index fb1db64d2d..5845a7f634 100644 --- a/singularity-eos/eos/get_sg_eos_p_t.cpp +++ b/singularity-eos/eos/get_sg_eos_p_t.cpp @@ -39,6 +39,10 @@ void get_sg_eos_p_t(const char *name, int ncell, int nmat, indirection_v &offset // for small loops const int32_t token{tokens.acquire()}; const int32_t tid{small_loop ? iloop : token}; + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } // caching mechanism singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0)); double mass_sum{0.0}; diff --git a/singularity-eos/eos/get_sg_eos_rho_e.cpp b/singularity-eos/eos/get_sg_eos_rho_e.cpp index 43049e31a7..95cb97b2ed 100644 --- a/singularity-eos/eos/get_sg_eos_rho_e.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_e.cpp @@ -39,6 +39,10 @@ void get_sg_eos_rho_e(const char *name, int ncell, indirection_v &offsets_v, int npte{0}; // initialize values for solver / lookup i_func(i, tid, mass_sum, npte, 0.0, 1.0, 0.0); + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } // get cache from offsets into scratch const int neq = npte + 1; singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) + diff --git a/singularity-eos/eos/get_sg_eos_rho_p.cpp b/singularity-eos/eos/get_sg_eos_rho_p.cpp index dc6da80ecc..4a788b1897 100644 --- a/singularity-eos/eos/get_sg_eos_rho_p.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_p.cpp @@ -41,6 +41,10 @@ void get_sg_eos_rho_p(const char *name, int ncell, indirection_v &offsets_v, // initialize values for solver / lookup i_func(i, tid, mass_sum, npte, 0.0, 0.0, 1.0); Real sie_tot_true{0.0}; + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } const int neq = npte + 1; singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) + neq * (neq + 4) + 2 * npte); diff --git a/singularity-eos/eos/get_sg_eos_rho_t.cpp b/singularity-eos/eos/get_sg_eos_rho_t.cpp index f1b91eed1f..7d79bc1b2c 100644 --- a/singularity-eos/eos/get_sg_eos_rho_t.cpp +++ b/singularity-eos/eos/get_sg_eos_rho_t.cpp @@ -42,6 +42,10 @@ void get_sg_eos_rho_t(const char *name, int ncell, indirection_v &offsets_v, i_func(i, tid, mass_sum, npte, 1.0, 0.0, 0.0); // calculate pte condition (lookup for 1 mat cell) Real sie_tot_true{0.0}; + // need to initialize the scratch before it's used to avoid undefined behavior + for (int idx = 0; idx < solver_scratch.extent(1); ++idx) { + solver_scratch(tid, idx) = 0.0; + } const int neq = npte; singularity::mix_impl::CacheAccessor cache(&solver_scratch(tid, 0) + neq * (neq + 4) + 2 * npte); From b6ad414fd5401f058fc42f44760e935f2b31df77 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:21:17 -0600 Subject: [PATCH 07/53] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a49bee9d5..6ef62fd2dc 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -21,6 +21,7 @@ - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface - [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Fix missing utilization of E0 in Davis Products EOS +- [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls From ed042bd5e4612f622d6d2d421a76c0b8fc29fc2c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:23:48 -0600 Subject: [PATCH 08/53] Add positive temperature check for debug builds --- singularity-eos/closure/mixed_cell_models.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index cd3c094f8a..cbf8b0a58f 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -182,6 +182,7 @@ class PTESolverBase { void Finalize() { for (int m = 0; m < nmat; m++) { temp[m] *= Tnorm; + PORTABLE_REQUIRE(temp[m] > 0., "Non-positive temperature returned"); u[m] *= uscale; press[m] *= uscale; } From 8c0efef9bc8837ebc748b47b72712a93a9550853 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:31:50 -0600 Subject: [PATCH 09/53] Revert Davis EOS changes --- singularity-eos/eos/eos_davis.hpp | 4 +--- 1 file changed, 1 insertion(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 86eedd224b..87e2ef06f6 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -292,7 +292,6 @@ class DavisProducts : public EosBase { inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } static std::string EosPyType() { return EosType(); } - // TODO (JHP): Create helper function to find the CJ state given the reference state private: static constexpr Real onethird = 1.0 / 3.0; @@ -314,8 +313,7 @@ class DavisProducts : public EosBase { const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / - std::pow(vvc, _k - 1.0 + _a) - - _E0; + std::pow(vvc, _k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { const Real vvc = 1 / (rho * _vc); From dd749abf7d2057dff9ec25a6f736e41d1e9f58d0 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:32:42 -0600 Subject: [PATCH 10/53] Revert Davis doc changes --- doc/sphinx/src/models.rst | 36 +++++------------------------------- 1 file changed, 5 insertions(+), 31 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 19424ce45c..12f555caac 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -11,8 +11,6 @@ .. _DavisReactants: https://doi.org/10.1016/S0010-2180(99)00112-1 -.. _DavisProducts: https://doi.org/10.1063/1.2035310 - .. _ProbingOffHugoniotStates: https://doi.org/10.1063/1.4939675 .. _WillsThermo: https://www.osti.gov/biblio/1561015 @@ -1196,11 +1194,10 @@ Davis Products EOS .. warning:: Entropy is not yet available for this EOS -The `Davis products EOS `_ is created from the reference -isentrope passing through the CJ state of the high explosive along with a -constant heat capacity. The constant heat capacity leads to the energy being a -simple funciton of the temperature deviation from the reference isentrope such -that +The Davis products EOS is created from the reference isentrope passing through +the CJ state of the high explosive along with a constant heat capacity. The +constant heat capacity leads to the energy being a simple funciton of the +temperature deviation from the reference isentrope such that .. math:: @@ -1231,7 +1228,7 @@ Finally, the pressure, energy, and temperature along the isentrope are given by .. math:: - e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} - e_0 + e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} .. math:: @@ -1257,29 +1254,6 @@ Here, there are four dimensionless parameters that are settable by the user, :math:`e_\mathrm{C}`, :math:`V_\mathrm{C}` and :math:`T_\mathrm{C}` are tuning parameters with units related to their non-subscripted counterparts. -Note that the energy zero (i.e. the reference energy) for the Davis products EOS -is arbitrary. For the isentrope to properly pass through the CJ state of a -reacting material, the energy release of the reaction needs to be accounted for -properly. If done external to the EOS, an energy source term is required in the -Euler equations. However, a common convention is to specify the reactants and -product EOS in a consistent way such that the reference energy corresponds to -the rest state of the material *before* it reacts. - -The energy at the CJ state can be calculated as - -.. math:: - - e_\mathrm{CJ} = \frac{P_0 + P_\mathrm{CJ}}{2(V_0 - V_\mathrm{CJ})}, - -relative to :math:`e = 0` at the reference state of the *reactants*. Therefore -the :math:`e_0` energy offset of the products EOS is given by - -.. math:: - - e_0 = e_S(V_\mathrm{CJ}) - e_\mathrm{CJ}. - -Practically, this means :math:`e_0` should be positive for any energetic material. - The constructor for the Davis Products EOS is .. code-block:: cpp From 4a8e82491d60a10c6794512670d254cc0ef94760 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:43:55 -0600 Subject: [PATCH 11/53] Change to non-negative temperature check rather than strictly positive --- singularity-eos/closure/mixed_cell_models.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index cbf8b0a58f..d1ca6f89db 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -182,7 +182,7 @@ class PTESolverBase { void Finalize() { for (int m = 0; m < nmat; m++) { temp[m] *= Tnorm; - PORTABLE_REQUIRE(temp[m] > 0., "Non-positive temperature returned"); + PORTABLE_REQUIRE(temp[m] >= 0., "Non-positive temperature returned"); u[m] *= uscale; press[m] *= uscale; } From 733d12f1dff8f7b6d083229b26a24d7b5718d75f Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:45:42 -0600 Subject: [PATCH 12/53] Remove Davis EOS change --- CHANGELOG.md | 1 - 1 file changed, 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6ef62fd2dc..0220ebcb58 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,7 +20,6 @@ - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface -- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Fix missing utilization of E0 in Davis Products EOS - [[PR373]](https://github.com/lanl/singularity-eos/pull/373) Initialize cache in `get_sg_eos*` functions ### Changed (changing behavior/API/variables/...) From 83a658b43f6d97d5e321ecc88e269fbd466af1ca Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 18:54:05 -0600 Subject: [PATCH 13/53] Add robustness to calculations --- singularity-eos/eos/eos_davis.hpp | 72 ++++++++++++++++++++++--------- 1 file changed, 52 insertions(+), 20 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87e2ef06f6..1e6c7b6efb 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -18,6 +18,7 @@ #include #include +#include #include #include #include @@ -42,16 +43,21 @@ class DavisReactants : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { const Real es = Es(rho); - const Real tmp = std::pow((1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0, - 1.0 / (1.0 + _alpha)); - if (tmp > 0) return Ts(rho) * tmp; - return Ts(rho) + (sie - es) / _Cv0; // This branch is a negative temperature + const Real power_base = (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0; + if (power_base <= 0) { + // This case would result in an imaginary temperature (i.e. negative), but we won't + // allow that so return zero + return 0.; + } + const Real tmp = std::pow(power_base, 1.0 / (1.0 + _alpha)); + return Ts(rho) * tmp; } template PORTABLE_INLINE_FUNCTION Real InternalEnergyFromDensityTemperature( const Real rho, const Real temp, Indexer_t &&lambda = static_cast(nullptr)) const { const Real t_s = Ts(rho); + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); return Es(rho) + _Cv0 * t_s / (1.0 + _alpha) * (std::pow(temp / t_s, 1.0 + _alpha) - 1.0); } @@ -72,8 +78,11 @@ class DavisReactants : public EosBase { template PORTABLE_INLINE_FUNCTION Real MinInternalEnergyFromDensity( const Real rho, Indexer_t &&lambda = static_cast(nullptr)) const { - MinInternalEnergyIsNotEnabled("DavisReactants"); - return 0.0; + // Minimum enegy is when the returned temperature is zero. This only happens + // when the base to the exponent is zero (see T(rho, e) equation) + const Real es = Es(rho); + const Real ts = Ts(rho); + return es - (1 + _alpha) / (_Cv0 * ts); } template PORTABLE_INLINE_FUNCTION Real @@ -100,8 +109,12 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - return _Cv0 / std::pow((1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1, - -_alpha / (1 + _alpha)); + const Real power_base = (1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1; + if (power_base <= 0) { + // Return zero heat capacity instead of an imaginary value + return 0.; + } + return _Cv0 / std::pow(power_base, -_alpha / (1 + _alpha)); } template PORTABLE_INLINE_FUNCTION Real BulkModulusFromDensityTemperature( @@ -149,7 +162,7 @@ class DavisReactants : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisReactants Params: "}; printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e " @@ -283,7 +296,7 @@ class DavisProducts : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisProducts Params: "}; printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc, @@ -300,15 +313,24 @@ class DavisProducts : public EosBase { static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; PORTABLE_INLINE_FUNCTION Real F(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1.0 / (rho * _vc); return 2.0 * _a / (std::pow(vvc, 2 * _n) + 1.0); } PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); return _pc * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / std::pow(vvc, _k + _a) * (_k - 1.0 + F(rho)) / (_k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); @@ -316,6 +338,9 @@ class DavisProducts : public EosBase { std::pow(vvc, _k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); return std::pow(2.0, -_a * _b / _n) * _pc * _vc / (_Cv * (_k - 1 + _a)) * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n * (1 - _b)) / @@ -328,7 +353,7 @@ class DavisProducts : public EosBase { PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; - const Real y = 1.0 - _rho0 / rho; + const Real y = 1.0 - robust::ratio(_rho0, rho); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4.0 * _B * y; @@ -342,7 +367,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { } } PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { - const Real y = 1 - _rho0 / rho; + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4 * _B * y; Real e_s; @@ -354,19 +379,19 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { } else { e_s = -y - (1.0 - std::exp(b4y)) / (4.0 * _B); } - return _e0 + _P0 * (1.0 / _rho0 - 1.0 / rho) + phat / _rho0 * e_s; + return _e0 + _P0 * (1.0 / _rho0 - robust::ratio(1.0, std::max(rho, 0.))) + phat / _rho0 * e_s; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const { if (rho >= _rho0) { - const Real y = 1 - _rho0 / rho; - return _T0 * std::exp(-_Z * y) * std::pow(_rho0 / rho, -_G0 - _Z); + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); + return _T0 * std::exp(-_Z * y) * std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); } else { - return _T0 * std::pow(_rho0 / rho, -_G0); + return _T0 * std::pow(robust::ratio(_rho0, rho), -_G0); } } PORTABLE_INLINE_FUNCTION Real DavisReactants::Gamma(const Real rho) const { if (rho >= _rho0) { - const Real y = 1 - _rho0 / rho; + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); return _G0 + _Z * y; } else { return _G0; @@ -377,7 +402,7 @@ template PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda) const { using namespace math_utils; - const Real y = 1 - _rho0 / rho; + const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4 * _B * y; const Real gamma = Gamma(rho); @@ -389,8 +414,8 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner 3 * y / pow<4>(1 - y) + 4 * pow<2>(y) / pow<5>(1 - y)) : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; - return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * rho) - gamma * rho * esv) / - rho; + return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - robust::ratio(gamma * std::max(rho, 0.) * esv), + rho); } template @@ -398,6 +423,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { // First, solve P=P(rho,T) for rho. Note P(rho,e) has an sie-es term, which is only a // function of T + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided") auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) * (std::pow(temp / Ts(r), 1 + _alpha) - 1.0)); @@ -418,6 +444,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::FillEos(Real &rho, Real &temp, Rea Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda) const { + // BROKEN: This can only handle density-energy inputs! MUST BE CHANGED if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); if (output & thermalqs::temperature) temp = TemperatureFromDensityInternalEnergy(rho, sie); @@ -450,6 +477,9 @@ template PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda) const { using namespace math_utils; + if (rho <= 0) { + return 0.; + } const Real vvc = 1 / (rho * _vc); const Real Fx = -4 * _a * std::pow(vvc, 2 * _n - 1) / pow<2>(1 + std::pow(vvc, 2 * _n)); const Real tmp = std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / @@ -470,6 +500,7 @@ PORTABLE_INLINE_FUNCTION Real DavisProducts::BulkModulusFromDensityInternalEnerg template PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperature( const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * r * _Cv * (temp - Ts(r))); }; @@ -488,6 +519,7 @@ template PORTABLE_INLINE_FUNCTION void DavisProducts::FillEos(Real &rho, Real &temp, Real &sie, Real &press, Real &cv, Real &bmod, const unsigned long output, Indexer_t &&lambda) const { + // BROKEN: This can only handle density-energy inputs! MUST BE CHANGED if (output & thermalqs::pressure) press = PressureFromDensityInternalEnergy(rho, sie); if (output & thermalqs::temperature) temp = TemperatureFromDensityInternalEnergy(rho, sie); From 3ce28ef4ea2f8698c6389b1e42d577fe630bbaa2 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 18:54:41 -0600 Subject: [PATCH 14/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 1e6c7b6efb..c767f186b1 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -379,12 +379,14 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { } else { e_s = -y - (1.0 - std::exp(b4y)) / (4.0 * _B); } - return _e0 + _P0 * (1.0 / _rho0 - robust::ratio(1.0, std::max(rho, 0.))) + phat / _rho0 * e_s; + return _e0 + _P0 * (1.0 / _rho0 - robust::ratio(1.0, std::max(rho, 0.))) + + phat / _rho0 * e_s; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const { if (rho >= _rho0) { const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); - return _T0 * std::exp(-_Z * y) * std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); + return _T0 * std::exp(-_Z * y) * + std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); } else { return _T0 * std::pow(robust::ratio(_rho0, rho), -_G0); } @@ -414,8 +416,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner 3 * y / pow<4>(1 - y) + 4 * pow<2>(y) / pow<5>(1 - y)) : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; - return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - robust::ratio(gamma * std::max(rho, 0.) * esv), - rho); + return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - + robust::ratio(gamma * std::max(rho, 0.) * esv), + rho); } template From 02b089ee0cd84dac894cc2190817d294250e3df4 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 18:59:22 -0600 Subject: [PATCH 15/53] Add robust_utils.hpp header --- singularity-eos/eos/eos_davis.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index c767f186b1..093125e322 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -22,6 +22,7 @@ #include #include #include +#include #include namespace singularity { From b50971ebe41972387ff2e5d2bf56b715149c08ce Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 19:05:54 -0600 Subject: [PATCH 16/53] Update changelog --- CHANGELOG.md | 1 + 1 file changed, 1 insertion(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 6800600af4..12fbe42a13 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface +- [[PR374]](https://github.com/lanl/singularity-eos/pull/374) Make the Davis EOS more numerically robust ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls From 1feb551c46fdd3781a7aca4e6203dc4f598ea6da Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 19:06:38 -0600 Subject: [PATCH 17/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 093125e322..80482da087 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -21,8 +21,8 @@ #include #include #include -#include #include +#include #include namespace singularity { From d1714e421c55d945bd7e818df6e800574bf7885c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:45:54 -0600 Subject: [PATCH 18/53] Fix some silly mistakes --- singularity-eos/eos/eos_davis.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 80482da087..d3cb8e0eb0 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -163,7 +163,7 @@ class DavisReactants : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisReactants Params: "}; printf("%srho0:%e e0:%e P0:%e\nT0:%e A:%e B:%e\nC:%e G0:%e Z:%e\nalpha:%e " @@ -417,9 +417,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner 3 * y / pow<4>(1 - y) + 4 * pow<2>(y) / pow<5>(1 - y)) : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; - return -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - - robust::ratio(gamma * std::max(rho, 0.) * esv), - rho); + const Real numerator = -(psv + (sie - Es(rho)) * rho * (gammav - gamma * + std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); + return robust::ratio(numerator, rho); } template From 62fc6824a4a862c6613787913dd35671d2a0559e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:46:23 -0600 Subject: [PATCH 19/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d3cb8e0eb0..d6f382b1c4 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -417,8 +417,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner 3 * y / pow<4>(1 - y) + 4 * pow<2>(y) / pow<5>(1 - y)) : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; - const Real numerator = -(psv + (sie - Es(rho)) * rho * (gammav - gamma * - std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); + const Real numerator = + -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - + gamma * std::max(rho, 0.) * esv); return robust::ratio(numerator, rho); } From 4d78e474d5cf1dde7f36b6932e3baa75feeb7a47 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:48:24 -0600 Subject: [PATCH 20/53] Another silly find/replace --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d6f382b1c4..b66d1dfd13 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -297,7 +297,7 @@ class DavisProducts : public EosBase { static inline unsigned long scratch_size(std::string method, unsigned int nelements) { return 0; } - static inline unsigned long std::max_scratch_size(unsigned int nelements) { return 0; } + static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisProducts Params: "}; printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc, From ac2a8bcf21b2256ad19c9900bb661cb191f76d34 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 20:53:28 -0600 Subject: [PATCH 21/53] Last silly mistake hopefully... --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index b66d1dfd13..027138c8da 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -428,7 +428,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { // First, solve P=P(rho,T) for rho. Note P(rho,e) has an sie-es term, which is only a // function of T - PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided") + PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) * (std::pow(temp / Ts(r), 1 + _alpha) - 1.0)); From 99cfb4b325558e0b4c1c1ea53daf19c1328dbf58 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:31:56 -0600 Subject: [PATCH 22/53] Reference shift modifier for Davis Products EOS --- doc/sphinx/src/models.rst | 9 +++++++-- doc/sphinx/src/modifiers.rst | 2 ++ 2 files changed, 9 insertions(+), 2 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 19424ce45c..163b2d9721 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1231,7 +1231,7 @@ Finally, the pressure, energy, and temperature along the isentrope are given by .. math:: - e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} - e_0 + e_S(\rho) = e_{\mathrm{C}} G(\rho) \frac{1}{\rho V_{\mathrm{C}}} .. math:: @@ -1272,7 +1272,7 @@ The energy at the CJ state can be calculated as e_\mathrm{CJ} = \frac{P_0 + P_\mathrm{CJ}}{2(V_0 - V_\mathrm{CJ})}, relative to :math:`e = 0` at the reference state of the *reactants*. Therefore -the :math:`e_0` energy offset of the products EOS is given by +the energy offset of the products EOS is given by .. math:: @@ -1280,6 +1280,11 @@ the :math:`e_0` energy offset of the products EOS is given by Practically, this means :math:`e_0` should be positive for any energetic material. +To provide the energy offset to the Davis Products EOS, `the energy shift +modifier`_ should be used. Note that the convention there +is that the shift is positive, so :math:`-e_0` should be provided to the shift +modifier. + The constructor for the Davis Products EOS is .. code-block:: cpp diff --git a/doc/sphinx/src/modifiers.rst b/doc/sphinx/src/modifiers.rst index 552b16940b..75e11310af 100644 --- a/doc/sphinx/src/modifiers.rst +++ b/doc/sphinx/src/modifiers.rst @@ -22,6 +22,8 @@ We list below the available modifiers and their constructors. The Shifted EOS ----------------- +.. _modifiers shifted EOS: + The shifted equation of state modifies zero point energy of an underlying model by some shift. So for example, it transforms From 83700480eb8769959f92bfd5d5883785ab1f3f96 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:32:58 -0600 Subject: [PATCH 23/53] Remove E0 from input --- singularity-eos/eos/eos_davis.hpp | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 86eedd224b..993d8d5b5b 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -177,8 +177,8 @@ class DavisProducts : public EosBase { DavisProducts() = default; PORTABLE_INLINE_FUNCTION DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc, - const Real pc, const Real Cv, const Real E0) - : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), _E0(E0) {} + const Real pc, const Real Cv) + : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), {} DavisProducts GetOnDevice() { return *this; } template PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( @@ -287,7 +287,7 @@ class DavisProducts : public EosBase { PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisProducts Params: "}; printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc, - _pc, _Cv, _E0); + _pc, _Cv); } inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } @@ -296,7 +296,7 @@ class DavisProducts : public EosBase { private: static constexpr Real onethird = 1.0 / 3.0; - Real _a, _b, _k, _n, _vc, _pc, _Cv, _E0; + Real _a, _b, _k, _n, _vc, _pc, _Cv; // static constexpr const char _eos_type[] = "DavisProducts"; static constexpr const unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; @@ -314,8 +314,7 @@ class DavisProducts : public EosBase { const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / - std::pow(vvc, _k - 1.0 + _a) - - _E0; + std::pow(vvc, _k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { const Real vvc = 1 / (rho * _vc); From 08da5dde4d53c742269f1c9f70162f30c778964b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:33:53 -0600 Subject: [PATCH 24/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 993d8d5b5b..6d3d413211 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -314,7 +314,7 @@ class DavisProducts : public EosBase { const Real ec = _pc * _vc / (_k - 1.0 + _a); // const Real de = ecj-(Es(rho0)-_E0); return ec * std::pow(0.5 * (std::pow(vvc, _n) + std::pow(vvc, -_n)), _a / _n) / - std::pow(vvc, _k - 1.0 + _a); + std::pow(vvc, _k - 1.0 + _a); } PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const { const Real vvc = 1 / (rho * _vc); From 257e33fcd85f9486d46fe15868f1fdbe710615fe Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:35:54 -0600 Subject: [PATCH 25/53] Typo --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 6d3d413211..5691e10f39 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -178,7 +178,7 @@ class DavisProducts : public EosBase { PORTABLE_INLINE_FUNCTION DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc, const Real pc, const Real Cv) - : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv), {} + : _a(a), _b(b), _k(k), _n(n), _vc(vc), _pc(pc), _Cv(Cv) {} DavisProducts GetOnDevice() { return *this; } template PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( From 0032e6793d4c602f9032daa514d80e1926819be1 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:40:48 -0600 Subject: [PATCH 26/53] Remove E0 from Davis Products EOS --- singularity-eos/eos/singularity_eos.cpp | 10 +++++----- singularity-eos/eos/singularity_eos.f90 | 4 ++-- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index e000bfe268..c7c3769bb7 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -154,23 +154,23 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0, + const double pc, const double Cv, int const *const enabled, double *const vals) { assert(matindex >= 0); - EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); + EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv)); if (enabled[3] == 1) { singularity::pAlpha2BilinearRampParams(eosi, vals[2], vals[3], vals[4], vals[2], vals[3], vals[4], vals[5]); } - EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv, E0)); + EOS eos_ = SGAPPLYMOD(DavisProducts(a, b, k, n, vc, pc, Cv)); eos[matindex] = eos_.GetOnDevice(); return 0; } int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0) { - return init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, E0, def_en, def_v); + const double pc, const double Cv) { + return init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, def_en, def_v); } int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index 2c0ef227ec..ae131206b1 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -112,13 +112,13 @@ end function init_sg_JWL interface integer(kind=c_int) function & - init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, E0, & + init_sg_DavisProducts(matindex, eos, a, b, k, n, vc, pc, Cv, & sg_mods_enabled, sg_mods_values) & bind(C, name='init_sg_DavisProducts') import integer(c_int), value, intent(in) :: matindex type(c_ptr), value, intent(in) :: eos - real(kind=c_double), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0 + real(kind=c_double), value, intent(in) :: a, b, k, n, vc, pc, Cv type(c_ptr), value, intent(in) :: sg_mods_enabled, sg_mods_values end function init_sg_DavisProducts end interface From 01627d90fc5b2d9b8bc4d1e3c83863b727c18362 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:41:50 -0600 Subject: [PATCH 27/53] Remove E0 from Davis Products signature --- doc/sphinx/src/models.rst | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/doc/sphinx/src/models.rst b/doc/sphinx/src/models.rst index 163b2d9721..fe37f8e332 100644 --- a/doc/sphinx/src/models.rst +++ b/doc/sphinx/src/models.rst @@ -1290,12 +1290,11 @@ The constructor for the Davis Products EOS is .. code-block:: cpp DavisProducts(const Real a, const Real b, const Real k, const Real n, const Real vc, - const Real pc, const Real Cv, const Real E0) + const Real pc, const Real Cv) where ``a`` is :math:`a`, ``b`` is :math:`b`, ``k`` is :math:`k`, ``n`` is :math:`n`, ``vc`` is :math:`V_\mathrm{C}`, ``pc`` is -:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`, and ``E0`` is -:math:`e_\mathrm{C}`. +:math:`P_\mathrm{C}`, ``Cv`` is :math:`C_{V,0}`. Spiner EOS ```````````` From 234512fa199e2fff9f900055e2a279a3aca490f1 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:43:14 -0600 Subject: [PATCH 28/53] Clang format --- singularity-eos/eos/singularity_eos.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.cpp b/singularity-eos/eos/singularity_eos.cpp index c7c3769bb7..c80ab710b5 100644 --- a/singularity-eos/eos/singularity_eos.cpp +++ b/singularity-eos/eos/singularity_eos.cpp @@ -154,8 +154,8 @@ int init_sg_JWL(const int matindex, EOS *eos, const double A, const double B, int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, - int const *const enabled, double *const vals) { + const double pc, const double Cv, int const *const enabled, + double *const vals) { assert(matindex >= 0); EOS eosi = SGAPPLYMODSIMPLE(DavisProducts(a, b, k, n, vc, pc, Cv)); if (enabled[3] == 1) { From 571fe88b88e7c911071c92a8c8d509fd9e82f29f Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 21:45:30 -0600 Subject: [PATCH 29/53] Remove final Davis Products EO --- singularity-eos/eos/singularity_eos.f90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.f90 b/singularity-eos/eos/singularity_eos.f90 index ae131206b1..e325601758 100644 --- a/singularity-eos/eos/singularity_eos.f90 +++ b/singularity-eos/eos/singularity_eos.f90 @@ -460,12 +460,12 @@ integer function init_sg_JWL_f(matindex, eos, A, B, R1, R2, w, rho0, Cv, & end function init_sg_JWL_f integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & - Cv, E0, sg_mods_enabled, & + Cv, sg_mods_enabled, & sg_mods_values) & result(err) integer(c_int), value, intent(in) :: matindex type(sg_eos_ary_t), intent(in) :: eos - real(kind=8), value, intent(in) :: a, b, k, n, vc, pc, Cv, E0 + real(kind=8), value, intent(in) :: a, b, k, n, vc, pc, Cv integer(kind=c_int), dimension(:), target, optional, intent(inout) :: sg_mods_enabled real(kind=8), dimension(:), target, optional, intent(inout) :: sg_mods_values ! local vars @@ -478,7 +478,7 @@ integer function init_sg_DavisProducts_f(matindex, eos, a, b, k, n, vc, pc, & if(present(sg_mods_values)) sg_mods_values_use = sg_mods_values err = init_sg_DavisProducts(matindex-1, eos%ptr, a, b, k, n, vc, pc, Cv, & - E0, c_loc(sg_mods_enabled_use), & + c_loc(sg_mods_enabled_use), & c_loc(sg_mods_values_use)) end function init_sg_DavisProducts_f From c519d29e0f7d05ae9b739dd1460ca437074c770b Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 22:19:36 -0600 Subject: [PATCH 30/53] Remove E0 from print parameters function --- singularity-eos/eos/eos_davis.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 5691e10f39..a00a9bfd4e 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -286,8 +286,8 @@ class DavisProducts : public EosBase { static inline unsigned long max_scratch_size(unsigned int nelements) { return 0; } PORTABLE_INLINE_FUNCTION void PrintParams() const { static constexpr char s1[]{"DavisProducts Params: "}; - printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e E0:%e\n", s1, _a, _b, _k, _n, _vc, - _pc, _Cv); + printf("%sa:%e b:%e k:%e\nn:%e vc:%e pc:%e\nCv:%e \n", s1, _a, _b, _k, _n, _vc, _pc, + _Cv); } inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } From 06289f4145182a07ec33b857086a0ad7c075e818 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 22:32:52 -0600 Subject: [PATCH 31/53] Remove E0 From Davis Products in python module --- python/module.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/module.cpp b/python/module.cpp index a6b3ad1c79..711e0990b3 100644 --- a/python/module.cpp +++ b/python/module.cpp @@ -62,7 +62,7 @@ PYBIND11_MODULE(singularity_eos, m) { .def( py::init(), py::arg("a"), py::arg("b"), py::arg("k"), py::arg("n"), py::arg("vc"), - py::arg("pc"), py::arg("Cv"), py::arg("E0") + py::arg("pc"), py::arg("Cv") ); #ifdef SPINER_USE_HDF From cbe4c76f60a764cc793acbcb1431cfc10681000a Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 11:24:45 -0600 Subject: [PATCH 32/53] Remove E0 argument from Davis Products --- test/pte_test_utils.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/pte_test_utils.hpp b/test/pte_test_utils.hpp index d61b09747e..af2e68ad94 100644 --- a/test/pte_test_utils.hpp +++ b/test/pte_test_utils.hpp @@ -59,7 +59,7 @@ inline void set_eos(T *eos) { singularity::EOS dr = singularity::DavisReactants( 1.890, 4.115e10, 1.0e6, 297.0, 1.8e5, 4.6, 0.34, 0.56, 0.0, 0.4265, 0.001074e10); singularity::EOS dp = singularity::DavisProducts(0.798311, 0.58, 1.35, 2.66182, 0.75419, - 3.2e10, 0.001072e10, 0.0); + 3.2e10, 0.001072e10); eos[0] = gr.GetOnDevice(); eos[1] = dr.GetOnDevice(); eos[2] = dp.GetOnDevice(); From bbf23e1e843fd79ce46f75befe710bb79a60d369 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 11:26:31 -0600 Subject: [PATCH 33/53] Update changelong to reflect that Davis Products E0 is now gone instead of 'fixed' --- CHANGELOG.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7a49bee9d5..6bbc4c9b0c 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,10 +20,10 @@ - [[PR335]](https://github.com/lanl/singularity-eos/pull/335) Fix missing hermite.hpp in CMake install required for Helmholtz EOS - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Guard against FPEs in the PTE solver - [[PR356]](https://github.com/lanl/singularity-eos/pull/356) Update CMake for proper Kokkos linking in Fortran interface -- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Fix missing utilization of E0 in Davis Products EOS ### Changed (changing behavior/API/variables/...) - [[PR363]](https://github.com/lanl/singularity-eos/pull/363) Template lambda values for scalar calls +- [[PR372]](https://github.com/lanl/singularity-eos/pull/372) Removed E0 from Davis Products EOS in favor of using the shifted EOS modifier. CHANGES API! ### Infrastructure (changes irrelevant to downstream codes) - [[PR329]](https://github.com/lanl/singularity-eos/pull/329) Move vinet tests into analytic test suite From 4eac409d506655219d480993ac2d52f353f9531e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 16:17:46 -0600 Subject: [PATCH 34/53] Remove E0 argument from Davis Products --- test/test_f_iface.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_f_iface.f90 b/test/test_f_iface.f90 index 52a92e95cc..7b845a2119 100644 --- a/test/test_f_iface.f90 +++ b/test/test_f_iface.f90 @@ -43,7 +43,7 @@ program test_sg_fortran_interface 4.6d0, 0.34d0, 0.56d0,0.d0, 0.4265d0, 0.001074d10) mat = mat + 1 res = init_sg_DavisProducts_f(mat, eos, 0.798311d0, 0.58d0, 1.35d0, 2.66182d0,& - 0.75419d0, 3.2d10, 0.001072d10, 0.d0) + 0.75419d0, 3.2d10, 0.001072d10) ! cleanup res = finalize_sg_eos_f(nmat, eos) From c265ca4b0c46da4834af0df2693016868b914424 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 16:56:19 -0600 Subject: [PATCH 35/53] Remove E0 from Davis Products EOS --- singularity-eos/eos/singularity_eos.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 9a7ddd43f7..5cec1491f4 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -39,8 +39,8 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0, - int const *const enabled, double *const vals); + const double pc, const double Cv, int const *const enabled, + double *const vals); int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double e0, const double P0, const double T0, From a6f09ba8d93a58da633a07210fed151b977596cc Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 16:58:15 -0600 Subject: [PATCH 36/53] Remove Davis Products EOS E0 parameter --- singularity-eos/eos/singularity_eos.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/singularity_eos.hpp b/singularity-eos/eos/singularity_eos.hpp index 5cec1491f4..9a31d10f54 100644 --- a/singularity-eos/eos/singularity_eos.hpp +++ b/singularity-eos/eos/singularity_eos.hpp @@ -134,7 +134,7 @@ int init_sg_Gruneisen(const int matindex, EOS *eos, const double C0, const doubl int init_sg_DavisProducts(const int matindex, EOS *eos, const double a, const double b, const double k, const double n, const double vc, - const double pc, const double Cv, const double E0); + const double pc, const double Cv); int init_sg_DavisReactants(const int matindex, EOS *eos, const double rho0, const double e0, const double P0, const double T0, From a0b5a9ffd7de6839c6179d90a12047bec682ae2d Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 18:22:19 -0600 Subject: [PATCH 37/53] Fix number of arguments to Davis Products EOS --- test/python_bindings.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/python_bindings.py b/test/python_bindings.py index a0ff479ad7..951e90523e 100644 --- a/test/python_bindings.py +++ b/test/python_bindings.py @@ -87,7 +87,7 @@ def testDavisReactants(self): eos = singularity_eos.DavisReactants(1,1,1,1,1,1,1,1,1,1,1) def testDavisProducts(self): - eos = singularity_eos.DavisProducts(1,1,1,1,1,1,1,1) + eos = singularity_eos.DavisProducts(1,1,1,1,1,1,1) class Modifiers(unittest.TestCase, EOSTestBase): def setUp(self): From 768f629c95178b08189fc999183b8b9439471574 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:18:18 -0600 Subject: [PATCH 38/53] Finally remove E0 from DavisProducts EOS pybinding init --- python/module.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/module.cpp b/python/module.cpp index 711e0990b3..2b89795552 100644 --- a/python/module.cpp +++ b/python/module.cpp @@ -60,7 +60,7 @@ PYBIND11_MODULE(singularity_eos, m) { eos_class(m, "DavisProducts") .def(py::init()) .def( - py::init(), + py::init(), py::arg("a"), py::arg("b"), py::arg("k"), py::arg("n"), py::arg("vc"), py::arg("pc"), py::arg("Cv") ); From 61edd15006d3b4a6242ee715d96e188b0e57e3d7 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:32:33 -0600 Subject: [PATCH 39/53] Fix missing maxes and remove last negative density branch --- singularity-eos/eos/eos_davis.hpp | 14 +++++--------- 1 file changed, 5 insertions(+), 9 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 027138c8da..8dc257272d 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -83,7 +83,7 @@ class DavisReactants : public EosBase { // when the base to the exponent is zero (see T(rho, e) equation) const Real es = Es(rho); const Real ts = Ts(rho); - return es - (1 + _alpha) / (_Cv0 * ts); + return es - (_Cv0 * ts) / (1 + _alpha); } template PORTABLE_INLINE_FUNCTION Real @@ -354,7 +354,7 @@ class DavisProducts : public EosBase { PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; - const Real y = 1.0 - robust::ratio(_rho0, rho); + const Real y = 1.0 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4.0 * _B * y; @@ -384,13 +384,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { phat / _rho0 * e_s; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const { - if (rho >= _rho0) { - const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); - return _T0 * std::exp(-_Z * y) * - std::pow(robust::ratio(_rho0, std::max(rho, 0.)), -_G0 - _Z); - } else { - return _T0 * std::pow(robust::ratio(_rho0, rho), -_G0); - } + const Real rho0overrho = robust::ratio(_rho0, std::max(rho, 0.)); + const Real y = 1 - rho0overrho; + return _T0 * std::exp(-_Z * y) * std::pow(rho0overrho, -_G0 - _Z); } PORTABLE_INLINE_FUNCTION Real DavisReactants::Gamma(const Real rho) const { if (rho >= _rho0) { From 4ec06fec56b4fad04abe495bf5c149a3588cb3e3 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:34:08 -0600 Subject: [PATCH 40/53] Add missing maxes and do a clang format --- singularity-eos/eos/eos_davis.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 8dc257272d..e80aea9839 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -404,8 +404,8 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A * _A / _B * _rho0; const Real b4y = 4 * _B * y; - const Real gamma = Gamma(rho); - const Real esv = -Ps(rho); + const Real gamma = Gamma(std::max(rho, 0.)); + const Real esv = -Ps(std::max(rho, 0.)); const Real psv = (rho >= _rho0) ? -phat * _rho0 * @@ -414,9 +414,9 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; const Real numerator = - -(psv + (sie - Es(rho)) * rho * (gammav - gamma * std::max(rho, 0.)) - + -(psv + (sie - Es(std::max(rho, 0.))) * rho * (gammav - gamma * std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); - return robust::ratio(numerator, rho); + return robust::ratio(numerator, std::max(rho, 0.)); } template From 4e2c4b469440133c5e023fbcc65e71362f2c2c76 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:44:34 -0600 Subject: [PATCH 41/53] Break out dimensionless energy into its own function --- singularity-eos/eos/eos_davis.hpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index e80aea9839..ca82e25824 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -44,7 +44,7 @@ class DavisReactants : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { const Real es = Es(rho); - const Real power_base = (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0; + const Real power_base = DimlessEdiff(sie); if (power_base <= 0) { // This case would result in an imaginary temperature (i.e. negative), but we won't // allow that so return zero @@ -110,7 +110,7 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real power_base = (1 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1; + const Real power_base = DimlessEdiff(sie); if (power_base <= 0) { // Return zero heat capacity instead of an imaginary value return 0.; @@ -180,6 +180,7 @@ class DavisReactants : public EosBase { // static constexpr const char _eos_type[] = "DavisReactants"; static constexpr unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; + PORTABLE_FORCEINLINE_FUNCTION Real DimlessEdiff(const Real sie) const; PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const; PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const; PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const; @@ -352,6 +353,10 @@ class DavisProducts : public EosBase { } }; +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const real sie) const { + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0 +} + PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; const Real y = 1.0 - robust::ratio(_rho0, std::max(rho, 0.)); From ad3cb40176631ea3a3f939810db62a9a1858972e Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:44:55 -0600 Subject: [PATCH 42/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index ca82e25824..19d5cb82c4 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,8 +354,7 @@ class DavisProducts : public EosBase { }; PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const real sie) const { - return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0 -} + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; From 88f7b861428281a6c4e59e75aa7f60880665fba0 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 19:56:17 -0600 Subject: [PATCH 43/53] Add max to another rho --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 19d5cb82c4..c7f70dd974 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -418,7 +418,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::BulkModulusFromDensityInternalEner : -phat * 4 * _B * _rho0 * std::exp(b4y); const Real gammav = (rho >= _rho0) ? _Z * _rho0 : 0.0; const Real numerator = - -(psv + (sie - Es(std::max(rho, 0.))) * rho * (gammav - gamma * std::max(rho, 0.)) - + -(psv + (sie - Es(rho)) * std::max(rho, 0.) * (gammav - gamma * std::max(rho, 0.)) - gamma * std::max(rho, 0.) * esv); return robust::ratio(numerator, std::max(rho, 0.)); } From c68391ba011145ac3025cbcba46be5045dd13fda Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:03:14 -0600 Subject: [PATCH 44/53] A few more robustness changes --- singularity-eos/eos/eos_davis.hpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index c7f70dd974..ace6dfb885 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -430,8 +430,8 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu // function of T PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { - return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) * - (std::pow(temp / Ts(r), 1 + _alpha) - 1.0)); + return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv0 * Ts(r) / (1 + _alpha) * + (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; @@ -507,7 +507,7 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { - return (Ps(r) + Gamma(r) * r * _Cv * (temp - Ts(r))); + return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv * (temp - Ts(r))); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; From 5a84cab4fcaae7e2a862b88660d2191a04eca218 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:03:38 -0600 Subject: [PATCH 45/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index ace6dfb885..b23db2358b 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -431,7 +431,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv0 * Ts(r) / (1 + _alpha) * - (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); + (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; From 43fb3ad7c8959ba9e183c13dadb0e8865e831b6a Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:07:55 -0600 Subject: [PATCH 46/53] Correct real to Real and switch to error instead of maxing density --- singularity-eos/eos/eos_davis.hpp | 14 +++++++++++--- 1 file changed, 11 insertions(+), 3 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index b23db2358b..87f843a649 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -353,7 +353,7 @@ class DavisProducts : public EosBase { } }; -PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const real sie) const { +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real sie) const { return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { @@ -430,7 +430,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu // function of T PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { - return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv0 * Ts(r) / (1 + _alpha) * + return (Ps(r) + Gamma(r) * r * _Cv0 * Ts(r) / (1 + _alpha) * (std::pow(robust::ratio(temp, Ts(r)), 1 + _alpha) - 1.0)); }; using RootFinding1D::regula_falsi; @@ -441,6 +441,10 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " "Root find failed to find a solution given P, T\n"); } + if (rho < 0.) { + EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " + "Root find resulted in a negative density") + } sie = InternalEnergyFromDensityTemperature(rho, temp); } @@ -507,7 +511,7 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur const Real press, const Real temp, Indexer_t &&lambda, Real &rho, Real &sie) const { PORTABLE_REQUIRE(temp >= 0, "Negative temperature provided"); auto PofRatT = [&](const Real r) { - return (Ps(r) + Gamma(r) * std::max(r, 0.) * _Cv * (temp - Ts(r))); + return (Ps(r) + Gamma(r) * r * _Cv * (temp - Ts(r))); }; using RootFinding1D::regula_falsi; using RootFinding1D::Status; @@ -518,6 +522,10 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur EOS_ERROR("DavisProducts::DensityEnergyFromPressureTemperature: " "Root find failed to find a solution given P, T\n"); } + if (rho < 0.) { + EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " + "Root find resulted in a negative density") + } sie = InternalEnergyFromDensityTemperature(rho, temp); } template From e30a855cc18f18bd3214861dc49c9b93a7ad16dc Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:12:47 -0600 Subject: [PATCH 47/53] Add missing density to DimlessEDiff function --- singularity-eos/eos/eos_davis.hpp | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 87f843a649..0d3281ee51 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -44,7 +44,7 @@ class DavisReactants : public EosBase { const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { const Real es = Es(rho); - const Real power_base = DimlessEdiff(sie); + const Real power_base = DimlessEdiff(rho, sie); if (power_base <= 0) { // This case would result in an imaginary temperature (i.e. negative), but we won't // allow that so return zero @@ -110,7 +110,7 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real SpecificHeatFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real power_base = DimlessEdiff(sie); + const Real power_base = DimlessEdiff(rho, sie); if (power_base <= 0) { // Return zero heat capacity instead of an imaginary value return 0.; @@ -180,7 +180,7 @@ class DavisReactants : public EosBase { // static constexpr const char _eos_type[] = "DavisReactants"; static constexpr unsigned long _preferred_input = thermalqs::density | thermalqs::specific_internal_energy; - PORTABLE_FORCEINLINE_FUNCTION Real DimlessEdiff(const Real sie) const; + PORTABLE_FORCEINLINE_FUNCTION Real DimlessEdiff(const Real rho, const Real sie) const; PORTABLE_INLINE_FUNCTION Real Ps(const Real rho) const; PORTABLE_INLINE_FUNCTION Real Es(const Real rho) const; PORTABLE_INLINE_FUNCTION Real Ts(const Real rho) const; @@ -353,7 +353,7 @@ class DavisProducts : public EosBase { } }; -PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real sie) const { +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { From 98541a73e22b5c9aedaac82575e29dcd03b03c79 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:13:22 -0600 Subject: [PATCH 48/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 0d3281ee51..d0813317fc 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -353,7 +353,8 @@ class DavisProducts : public EosBase { } }; -PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { +PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, + const Real sie) const { return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { From 5377650864e40a8b8add4eec7bb29180b0d5eae5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:15:08 -0600 Subject: [PATCH 49/53] Correct es to Es(rho) --- singularity-eos/eos/eos_davis.hpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d0813317fc..bd4cdc1742 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -43,7 +43,6 @@ class DavisReactants : public EosBase { PORTABLE_INLINE_FUNCTION Real TemperatureFromDensityInternalEnergy( const Real rho, const Real sie, Indexer_t &&lambda = static_cast(nullptr)) const { - const Real es = Es(rho); const Real power_base = DimlessEdiff(rho, sie); if (power_base <= 0) { // This case would result in an imaginary temperature (i.e. negative), but we won't @@ -355,7 +354,7 @@ class DavisProducts : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { - return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - es) + 1.0} + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; From 4f380404272075749d969e01bc9661de71ee09c4 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:20:12 -0600 Subject: [PATCH 50/53] Add missing semicolon --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index bab74467df..137e74542a 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,7 +354,7 @@ class DavisProducts : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { - return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0} + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0}; PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; From 6708ddada5538e6431e4711c4f6c56dc31e3a58c Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:22:16 -0600 Subject: [PATCH 51/53] Add missing bracket --- singularity-eos/eos/eos_davis.hpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 137e74542a..d53b4b8818 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,7 +354,8 @@ class DavisProducts : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { - return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0}; + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0; +} PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; From c422370b535963ad3ee239420a43f41242b96352 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:23:06 -0600 Subject: [PATCH 52/53] Clang format --- singularity-eos/eos/eos_davis.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index d53b4b8818..2c2987ee6e 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -354,7 +354,7 @@ class DavisProducts : public EosBase { PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, const Real sie) const { - return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0; + return (1.0 + _alpha) / (Ts(rho) * _Cv0) * (sie - Es(rho)) + 1.0; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { From ff3901dd0d301607739b9a52ff79361191f4fdf5 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Mon, 13 May 2024 20:25:28 -0600 Subject: [PATCH 53/53] Add missing semicolon --- singularity-eos/eos/eos_davis.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 2c2987ee6e..fb34c2a229 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -444,7 +444,7 @@ PORTABLE_INLINE_FUNCTION void DavisReactants::DensityEnergyFromPressureTemperatu } if (rho < 0.) { EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " - "Root find resulted in a negative density") + "Root find resulted in a negative density\n"); } sie = InternalEnergyFromDensityTemperature(rho, temp); } @@ -525,7 +525,7 @@ PORTABLE_INLINE_FUNCTION void DavisProducts::DensityEnergyFromPressureTemperatur } if (rho < 0.) { EOS_ERROR("DavisReactants::DensityEnergyFromPressureTemperature: " - "Root find resulted in a negative density") + "Root find resulted in a negative density\n"); } sie = InternalEnergyFromDensityTemperature(rho, temp); }