From 02a1e2980f894be8774b1521757ff43e76a730b8 Mon Sep 17 00:00:00 2001 From: Jeffrey H Peterson Date: Fri, 10 May 2024 17:19:53 -0600 Subject: [PATCH 1/7] 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 2/7] 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 3/7] 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 4/7] 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 5/7] 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 6/7] 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 7/7] 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/...)