diff --git a/.github/workflows/tests_clang.yml b/.github/workflows/tests_clang.yml index 1b6c3ef47e6..f0d56d917b8 100644 --- a/.github/workflows/tests_clang.yml +++ b/.github/workflows/tests_clang.yml @@ -46,6 +46,7 @@ jobs: -DSINGULARITY_USE_V_AND_V_EOS=OFF \ -DSINGULARITY_USE_KOKKOS=OFF \ -DSINGULARITY_PLUGINS=$(pwd)/../example/plugin \ + -DCMAKE_CXX_FLAGS="-Wall -Werr" \ -DCMAKE_LINKER=ld.gold \ -DCMAKE_BUILD_TYPE=Release \ .. diff --git a/CHANGELOG.md b/CHANGELOG.md index caaf1b3de4b..3a56b47d14f 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -20,6 +20,7 @@ - [[PR575]](https://github.com/lanl/singularity-eos/pull/575) Pin variant submodule to the same commit as the spackage ### Infrastructure (changes irrelevant to downstream codes) +- [[PR595]](https://github.com/lanl/singularity-eos/pull/595) A number of robustness fixes/cleanups. Add a warnings build to the CI. - [[PR594]](https://github.com/lanl/singularity-eos/pull/594) clean up tests in preparation for ports-of-call variant transition - [[PR588]](https://github.com/lanl/singularity-eos/pull/588) Add DensityEnergyFromPressureTemperature in unit system test - [[PR576]](https://github.com/lanl/singularity-eos/pull/576) Clean up some header includes to use the C++ versions diff --git a/singularity-eos/base/spiner_table_utils.hpp b/singularity-eos/base/spiner_table_utils.hpp index af88b4579b1..52f5a2c0209 100644 --- a/singularity-eos/base/spiner_table_utils.hpp +++ b/singularity-eos/base/spiner_table_utils.hpp @@ -292,7 +292,8 @@ class Bounds { // static methods of `SpinerTricks` are called from class methods of // spiner-based EOS's. template -struct SpinerTricks { +class SpinerTricks { // class instead of struct for -Wmismatched-tags + public: static auto GetOnDevice(EOS *peos_h) { // trivially copy all but dynamic memory EOS eos_d = *peos_h; diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp index d93ff92fe47..b4d6b6d0c83 100644 --- a/singularity-eos/closure/mixed_cell_models.hpp +++ b/singularity-eos/closure/mixed_cell_models.hpp @@ -495,7 +495,7 @@ class PTESolverBase { bad_vfrac_guess = false; for (std::size_t m = 0; m < nmat; ++m) { const Real rho_min = eos[m].RhoPmin(Tguess); - const Real rho_max = eos[m].MaximumDensity(); + [[maybe_unused]] const Real rho_max = eos[m].MaximumDensity(); PORTABLE_REQUIRE(rho_min < rho_max, "Valid density range must exist!"); PORTABLE_REQUIRE(rho[m] < rho_max, "Density must be less than rho_max"); if (rho[m] < rho_min) { diff --git a/singularity-eos/eos/eos_davis.hpp b/singularity-eos/eos/eos_davis.hpp index 5ad23b71dc1..84973fdd26e 100644 --- a/singularity-eos/eos/eos_davis.hpp +++ b/singularity-eos/eos/eos_davis.hpp @@ -25,6 +25,14 @@ #include #include +/* Implements the equation of state from + Wescott, Stewart, and Davis, 2005 + https://doi.org/10.1063/1.2035310 + + The expansion branch of the EOS is described by Aslam, 2018 + https://doi.org/10.1063/1.5020172 + */ + namespace singularity { using namespace eos_base; @@ -183,6 +191,8 @@ class DavisReactants : public EosBase { void inline Finalize() {} static std::string EosType() { return std::string("DavisReactants"); } static std::string EosPyType() { return EosType(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return robust::EPS(); } private: static constexpr Real onethird = 1.0 / 3.0; @@ -330,6 +340,8 @@ class DavisProducts : public EosBase { inline void Finalize() {} static std::string EosType() { return std::string("DavisProducts"); } static std::string EosPyType() { return EosType(); } + PORTABLE_FORCEINLINE_FUNCTION + Real MinimumDensity() const { return robust::EPS(); } private: static constexpr Real onethird = 1.0 / 3.0; @@ -389,6 +401,7 @@ PORTABLE_FORCEINLINE_FUNCTION Real DavisReactants::DimlessEdiff(const Real rho, PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { using namespace math_utils; + if (rho <= 0) return 0; const Real y = 1.0 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A2oB * _rho0; const Real b4y = 4.0 * _B * y; @@ -403,6 +416,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Ps(const Real rho) const { } } PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { + if (rho <= 0) return 0; const Real y = 1 - robust::ratio(_rho0, std::max(rho, 0.)); const Real phat = 0.25 * _A2oB * _rho0; const Real b4y = 4 * _B * y; @@ -419,6 +433,7 @@ PORTABLE_INLINE_FUNCTION Real DavisReactants::Es(const Real rho) const { phat * _spvol0 * e_s; } PORTABLE_INLINE_FUNCTION Real DavisReactants::Ts(const Real rho) const { + if (rho <= 0) return _T0; 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); diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp index 79332d90e5a..883b768739d 100644 --- a/singularity-eos/eos/eos_spiner_rho_temp.hpp +++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp @@ -102,7 +102,7 @@ class SpinerEOSDependsRhoT : public EosBase { reproducibility_mode, pmin_vapor_dome) {} PORTABLE_INLINE_FUNCTION SpinerEOSDependsRhoT() - : memoryStatus_(DataStatus::Deallocated), split_(TableSplit::Total) {} + : split_(TableSplit::Total), memoryStatus_(DataStatus::Deallocated) {} inline SpinerEOSDependsRhoT GetOnDevice(); @@ -492,8 +492,6 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, lTMax_ = P_.range(0).max(); TMax_ = from_log(lTMax_, lTOffset_); - Real rhoMin = from_log(lRhoMin_, lRhoOffset_); - // bulk modulus can be wrong in the tables. Use FLAG's approach to // fix the table. fixBulkModulus_(); @@ -522,7 +520,6 @@ inline herr_t SpinerEOSDependsRhoT::loadDataboxes_(const std::string &matid_str, for (int j = 0; j < numRho_; j++) { Real lRho = bModCold_.range(0).x(j); Real lT = lTColdCrit_(j); - Real rho = rho_(lRho); bModCold_(j) = bMod_.interpToReal(lRho, lT); dPdECold_(j) = dPdE_.interpToReal(lRho, lT); dTdRhoCold_(j) = dTdRho_.interpToReal(lRho, lT); @@ -792,7 +789,7 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::BulkModulusFromDensityIntern const Real gm1 = gm1Max_.interpToReal(lRho); bMod = (gm1 + 1) * gm1 * rho * sie; } else { // on table - bMod = bMod_.interpToReal(rho, sie); + bMod = bMod_.interpToReal(lRho, lT); } return bMod; } diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp index eaa79e03b88..311e565dea5 100644 --- a/singularity-eos/eos/eos_stellar_collapse.hpp +++ b/singularity-eos/eos/eos_stellar_collapse.hpp @@ -1107,9 +1107,7 @@ inline Real StellarCollapse::findMedian_(Real buffer[], int size) const { inline void StellarCollapse::computeBulkModulus_() { lBMod_.copyMetadata(lP_); for (int iY = 0; iY < numYe_; ++iY) { - Real Ye = lBMod_.range(2).x(iY); for (int iT = 0; iT < numT_; ++iT) { - Real lT = lBMod_.range(1).x(iT); for (int irho = 0; irho < numRho_; ++irho) { Real lRho = lBMod_.range(0).x(irho); Real rho = std::pow(10., lRho); // rho_(lRho); diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index c7f1a7da17e..30d7c303266 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -17,6 +17,7 @@ add_executable( eos_unit_test_helpers.hpp test_eos_ideal.cpp test_eos_gruneisen.cpp + test_eos_davis.cpp test_eos_sap_polynomial.cpp test_eos_noble_abel.cpp test_eos_stiff.cpp diff --git a/test/pte_test_3mat_analytic.hpp b/test/pte_test_3mat_analytic.hpp index 741680a30e8..f1437370f7f 100644 --- a/test/pte_test_3mat_analytic.hpp +++ b/test/pte_test_3mat_analytic.hpp @@ -17,15 +17,18 @@ #include +#include + #include #include +#include #include #include constexpr int NMAT = 3; constexpr int NTRIAL = 100; constexpr int NPTS = NTRIAL * NMAT; -constexpr std::size_t HIST_SIZE = 10; +constexpr std::size_t HIST_SIZE = 13; using singularity::DavisProducts; using singularity::DavisReactants; @@ -47,22 +50,31 @@ inline void set_eos(T *eos) { template inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie, - RealIndexer &&temp, EOSIndexer &&eos) { + RealIndexer &&temp, EOSIndexer &&eos, std::mt19937 &gen) { + constexpr Real alpha_min = singularity::robust::EPS(); rho[0] = 8.93; rho[1] = 1.89; rho[2] = 2.5; + // sample volume fractions from the Dirichlet(1,1,1) distribution, + // which more evenly samples edge case volume fractions + std::exponential_distribution dist(1.0); + Real vsum = 0.; for (int i = 0; i < NMAT; i++) { temp[i] = 600.0; sie[i] = eos[i].InternalEnergyFromDensityTemperature(rho[i], temp[i]); - vfrac[i] = rand() / (1.0 * RAND_MAX); + // may be >> 1 + vfrac[i] = dist(gen) + alpha_min; vsum += vfrac[i]; } - for (int i = 0; i < NMAT; i++) + for (int i = 0; i < NMAT; i++) { + // (0, 1), sums to 1 vfrac[i] *= 1.0 / vsum; + vfrac[i] = std::max(alpha_min, std::min(vfrac[i], 1.0)); + } return; } diff --git a/test/test_eos_davis.cpp b/test/test_eos_davis.cpp new file mode 100644 index 00000000000..f16705aec87 --- /dev/null +++ b/test/test_eos_davis.cpp @@ -0,0 +1,150 @@ +//------------------------------------------------------------------------------ +// © 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 +// Nuclear Security Administration. All rights in the program are +// reserved by Triad National Security, LLC, and the U.S. Department of +// Energy/National Nuclear Security Administration. The Government is +// granted for itself and others acting on its behalf a nonexclusive, +// paid-up, irrevocable worldwide license in this material to reproduce, +// prepare derivative works, distribute copies to the public, perform +// publicly and display publicly, and to permit others to do so. +//------------------------------------------------------------------------------ + +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#ifndef CATCH_CONFIG_FAST_COMPILE +#define CATCH_CONFIG_FAST_COMPILE +#include +#endif + +#include + +using singularity::DavisProducts; +using singularity::DavisReactants; +using EOS = singularity::Variant; + +constexpr int N = 100; +SCENARIO("Davis reactants EOS basic calls", "[DavisReactants]") { + GIVEN("A davis reactants equation of state") { + constexpr Real rho0 = 1.890; + constexpr Real e0 = 4.115e10; + constexpr Real P0 = 1.0e6; + constexpr Real T0 = 297.0; + constexpr Real A = 1.8e5; + constexpr Real B = 4.6; + constexpr Real C = 0.34; + constexpr Real G0 = 0.56; + constexpr Real Z = 0.0; + constexpr Real alpha = 0.4265; + constexpr Real Cv0 = 0.001074e10; + EOS eos_host = DavisReactants(rho0, e0, P0, T0, A, B, C, G0, Z, alpha, Cv0); + auto eos = eos_host.GetOnDevice(); + + WHEN("The EOS is in mild expansion") { + THEN("The Gruneisen coefficient is G0") { + int nwrong = 0; + portableReduce( + "Check gruneisen coefficient in expansion", 0, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real rho = 0.5 * rho0; + Real T = 10 * i; + Real gamma = eos.GruneisenParamFromDensityTemperature(rho, T); + nw += !(isClose(G0, gamma, 1e-12)); + }, + nwrong); + REQUIRE(nwrong == 0); + } + } + + WHEN("The EOS is in compression") { + THEN("The Gruneisen coefficient approaches G0 + 1") { + int nwrong = 0; + portableReduce( + "Check gruneisen coefficient in compression", 0, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real rho = 10 * rho0; + Real T = 10 * i; + Real gamma = eos.GruneisenParamFromDensityTemperature(rho, T); + const Real y = 1 - (rho0 / rho); + nw += !(isClose(G0 + Z * y, gamma, 1e-12)); + }, + nwrong); + REQUIRE(nwrong == 0); + } + } + + WHEN("The EOS is in extreme expansion") { + THEN("It still produces a sensible answer") { + int nwrong = 0; + portableReduce( + "Check properties in extreme expansion", 0, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real T = 10 * i; + Real rho = 0; + Real P = eos.PressureFromDensityTemperature(rho, T); + Real E = eos.InternalEnergyFromDensityTemperature(rho, T); + Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, T); + nw += std::isnan(P); + nw += std::isnan(E); + nw += std::isnan(bmod); + nw += (bmod < 0); + }, + nwrong); + REQUIRE(nwrong == 0); + } + } + + eos_host.Finalize(); + eos.Finalize(); + } +} + +SCENARIO("Davis products EOS basic calls", "[DavisProducts]") { + GIVEN("A davis reactants equation of state") { + constexpr Real a = 0.798311; + constexpr Real b = 0.58; + constexpr Real k = 1.35; + constexpr Real n = 2.66182; + constexpr Real vc = 0.75419; + constexpr Real pc = 3.2e10; + constexpr Real Cv = 0.001072e10; + EOS eos_host = DavisProducts(a, b, k, n, vc, pc, Cv); + auto eos = eos_host.GetOnDevice(); + + WHEN("The EOS is in extreme expansion") { + THEN("It still produces a sensible answer") { + int nwrong = 0; + portableReduce( + "Check properties in extreme expansion", 0, N, + PORTABLE_LAMBDA(const int i, int &nw) { + Real T = 10 * i; + Real rho = 0; + Real P = eos.PressureFromDensityTemperature(rho, T); + Real E = eos.InternalEnergyFromDensityTemperature(rho, T); + Real bmod = eos.BulkModulusFromDensityInternalEnergy(rho, T); + nw += std::isnan(P); + nw += std::isnan(E); + nw += std::isnan(bmod); + nw += (bmod < 0); + }, + nwrong); + REQUIRE(nwrong == 0); + } + } + + eos_host.Finalize(); + eos.Finalize(); + } +} diff --git a/test/test_pte.cpp b/test/test_pte.cpp index bee30f1adbb..fc8f1406bb6 100644 --- a/test/test_pte.cpp +++ b/test/test_pte.cpp @@ -15,6 +15,7 @@ #include #include #include +#include #include #include #include @@ -50,7 +51,6 @@ template