additional robustness cleanups

.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 \
                   ..

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

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 <typename EOS>
-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;

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) {

singularity-eos/eos/eos_davis.hpp

@@ -25,6 +25,14 @@
 #include <singularity-eos/base/root-finding-1d/root_finding.hpp>
 #include <singularity-eos/eos/eos_base.hpp>
 
+/* 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<DavisReactants> {
   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<DavisProducts> {
   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);

singularity-eos/eos/eos_spiner_rho_temp.hpp

@@ -102,7 +102,7 @@ class SpinerEOSDependsRhoT : public EosBase<SpinerEOSDependsRhoT> {
                              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;
 }

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);

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

test/pte_test_3mat_analytic.hpp

@@ -17,15 +17,18 @@
 
 #include <stdlib.h>
 
+#include <random>
+
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_arrays.hpp>
+#include <singularity-eos/base/robust_utils.hpp>
 #include <singularity-eos/eos/eos_models.hpp>
 #include <singularity-eos/eos/eos_variant.hpp>
 
 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 <typename RealIndexer, typename EOSIndexer>
 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<Real> 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;
 }

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 <algorithm>
+#include <cmath>
+#include <cstdio>
+#include <cstdlib>
+#include <limits>
+
+#include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_arrays.hpp>
+#include <ports-of-call/portable_errors.hpp>
+#include <singularity-eos/base/indexable_types.hpp>
+#include <singularity-eos/eos/eos.hpp>
+
+#ifndef CATCH_CONFIG_FAST_COMPILE
+#define CATCH_CONFIG_FAST_COMPILE
+#include <catch2/catch_test_macros.hpp>
+#endif
+
+#include <test/eos_unit_test_helpers.hpp>
+
+using singularity::DavisProducts;
+using singularity::DavisReactants;
+using EOS = singularity::Variant<DavisReactants, DavisProducts>;
+
+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();
+  }
+}

test/test_pte.cpp

@@ -15,6 +15,7 @@
 #include <chrono>
 #include <iostream>
 #include <memory>
+#include <random>
 #include <stdlib.h>
 #include <time.h>
 #include <utility>
@@ -50,7 +51,6 @@ template <template <typename... Types> class Method_t>
 void TestPTE(const std::string name, const std::size_t nscratch_vars,
              std::size_t &nsuccess, std::vector<Real> &host_vals) {
   constexpr Real EPS = 1e-5;
-  Real time;
   nsuccess = 0;
 
   // EOS
@@ -122,12 +122,16 @@ void TestPTE(const std::string name, const std::size_t nscratch_vars,
 #endif
 
   // setup state
+
+  // we want a repeatable random seed every time this function is
+  // called
+  std::mt19937 gen(123456789);
   for (int n = 0; n < NTRIAL; n++) {
     Indexer2D<decltype(rho_hm)> r(n, rho_hm);
     Indexer2D<decltype(vfrac_hm)> vf(n, vfrac_hm);
     Indexer2D<decltype(sie_hm)> e(n, sie_hm);
     Indexer2D<decltype(temp_hm)> t(n, temp_hm);
-    set_state(r, vf, e, t, eos_h);
+    set_state(r, vf, e, t, eos_h, gen);
   }
   for (int i = 0; i < HIST_SIZE; ++i) {
     hist_vh[i] = 0;
@@ -218,8 +222,6 @@ void TestPTE(const std::string name, const std::size_t nscratch_vars,
 #endif
 
   Real milliseconds = sum_time.count() / 1e3;
-  time = milliseconds;
-
   std::cout << "Finished " << NTRIAL << " solves in " << milliseconds << " milliseconds"
             << std::endl;
   std::cout << "Solves/ms = " << NTRIAL / milliseconds << std::endl;