From 26a1fe644672a849138a12eb9e540fd938ab4dce Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 09:44:54 -0500
Subject: [PATCH 01/16] test PTE should have a fixed seed for all solvers
 called, but that varies per trial

---
 test/test_pte.cpp | 7 ++++---
 1 file changed, 4 insertions(+), 3 deletions(-)

diff --git a/test/test_pte.cpp b/test/test_pte.cpp
index bee30f1adbb..92f757cbfa4 100644
--- a/test/test_pte.cpp
+++ b/test/test_pte.cpp
@@ -50,7 +50,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,6 +121,10 @@ 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
+  srand(123456789); 
   for (int n = 0; n < NTRIAL; n++) {
     Indexer2D<decltype(rho_hm)> r(n, rho_hm);
     Indexer2D<decltype(vfrac_hm)> vf(n, vfrac_hm);
@@ -218,8 +221,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;

From 3c9305a76c7a65b9e485ef45774ad110c874526a Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 09:45:21 -0500
Subject: [PATCH 02/16] Guard what volume fractions can be passed into the PTE
 solver

---
 test/pte_test_3mat_analytic.hpp | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/pte_test_3mat_analytic.hpp b/test/pte_test_3mat_analytic.hpp
index 741680a30e8..6d55e4a0933 100644
--- a/test/pte_test_3mat_analytic.hpp
+++ b/test/pte_test_3mat_analytic.hpp
@@ -48,7 +48,6 @@ 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) {
-
   rho[0] = 8.93;
   rho[1] = 1.89;
   rho[2] = 2.5;
@@ -57,7 +56,7 @@ inline void set_state(RealIndexer &&rho, RealIndexer &&vfrac, RealIndexer &&sie,
   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);
+    vfrac[i] = std::min(std::max(rand() / (1.0 * RAND_MAX), 1e-15), 1.0);
     vsum += vfrac[i];
   }
 

From 5799e857af52e9a5a86f683111bac4baf51253ec Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 09:46:20 -0500
Subject: [PATCH 03/16] address warnings from -Wall

---
 singularity-eos/base/spiner_table_utils.hpp  | 3 ++-
 singularity-eos/eos/eos_stellar_collapse.hpp | 2 --
 test/test_spiner_transform.cpp               | 1 -
 3 files changed, 2 insertions(+), 4 deletions(-)

diff --git a/singularity-eos/base/spiner_table_utils.hpp b/singularity-eos/base/spiner_table_utils.hpp
index af88b4579b1..4ca1c51f1ae 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 <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;
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/test_spiner_transform.cpp b/test/test_spiner_transform.cpp
index c904ae14b25..add2f78811f 100644
--- a/test/test_spiner_transform.cpp
+++ b/test/test_spiner_transform.cpp
@@ -150,7 +150,6 @@ SCENARIO("NullTransform behave correctly", "[TransformTest]") {
 
   GIVEN("A density with a known linear cold-curve in log-space") {
     Real rho = 10.0;
-    Real lRho = to_log(rho, data.lRhoOffset);
     Real e_actual = 42.0;
 
     THEN("NullTransform is identity throughout when run on device") {

From a037692386aaab3a8d6ad5b9427a97b0a8c8a35f Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 09:47:04 -0500
Subject: [PATCH 04/16] fix bug for bulk modulus from density internal energy

---
 singularity-eos/eos/eos_spiner_rho_temp.hpp | 7 ++-----
 1 file changed, 2 insertions(+), 5 deletions(-)

diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index 79332d90e5a..4fa9612c441 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<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;
 }

From c685ec5e3281a33a68e15d0959c51cfe1b5811c7 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 10:01:18 -0500
Subject: [PATCH 05/16] move rho_max into PORTABLE_REQUIRE to get rid of unused
 variable warning

---
 singularity-eos/closure/mixed_cell_models.hpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp
index d93ff92fe47..61a432b1ff3 100644
--- a/singularity-eos/closure/mixed_cell_models.hpp
+++ b/singularity-eos/closure/mixed_cell_models.hpp
@@ -495,9 +495,9 @@ 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();
         PORTABLE_REQUIRE(rho_min < rho_max, "Valid density range must exist!");
-        PORTABLE_REQUIRE(rho[m] < rho_max, "Density must be less than rho_max");
+        PORTABLE_REQUIRE(rho[m] < eos[m].MaximumDensity(),
+                         "Density must be less than rho_max");
         if (rho[m] < rho_min) {
           bad_vfrac_guess = true;
           Tguess *= Tfactor;

From f28c4885e53f7111ec573c974a4a217a067f38b1 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 14:23:04 -0500
Subject: [PATCH 06/16] Add protection for 0 density in Davis EOS. Also make
 davis EOS report that 0 is not a safe value for root finds.

---
 singularity-eos/eos/eos_davis.hpp | 15 +++++++++++++++
 1 file changed, 15 insertions(+)

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

From 1f0cd168557819ab75e9ec963392f26b34bd76aa Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 14:23:30 -0500
Subject: [PATCH 07/16] add test for davis EOS

---
 test/CMakeLists.txt     |   1 +
 test/test_eos_davis.cpp | 150 ++++++++++++++++++++++++++++++++++++++++
 2 files changed, 151 insertions(+)
 create mode 100644 test/test_eos_davis.cpp

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/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 <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();
+  }
+}

From 61badf52ca10c0796e2b6fde3a70de2b78e718e6 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 14:24:26 -0500
Subject: [PATCH 08/16] Switch from using rand to using a modern C++-based
 random number generator

---
 test/pte_test_3mat_analytic.hpp | 21 +++++++++++++++++----
 test/test_pte.cpp               |  5 +++--
 2 files changed, 20 insertions(+), 6 deletions(-)

diff --git a/test/pte_test_3mat_analytic.hpp b/test/pte_test_3mat_analytic.hpp
index 6d55e4a0933..f1437370f7f 100644
--- a/test/pte_test_3mat_analytic.hpp
+++ b/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,21 +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] = std::min(std::max(rand() / (1.0 * RAND_MAX), 1e-15), 1.0);
+    // 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_pte.cpp b/test/test_pte.cpp
index 92f757cbfa4..fc8f1406bb6 100644
--- a/test/test_pte.cpp
+++ b/test/test_pte.cpp
@@ -15,6 +15,7 @@
 #include <chrono>
 #include <iostream>
 #include <memory>
+#include <random>
 #include <stdlib.h>
 #include <time.h>
 #include <utility>
@@ -124,13 +125,13 @@ void TestPTE(const std::string name, const std::size_t nscratch_vars,
 
   // we want a repeatable random seed every time this function is
   // called
-  srand(123456789); 
+  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;

From 80f19e38326875a23d98175032dc24884e60768c Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 14:24:58 -0500
Subject: [PATCH 09/16] formatting

---
 singularity-eos/base/spiner_table_utils.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/base/spiner_table_utils.hpp b/singularity-eos/base/spiner_table_utils.hpp
index 4ca1c51f1ae..52f5a2c0209 100644
--- a/singularity-eos/base/spiner_table_utils.hpp
+++ b/singularity-eos/base/spiner_table_utils.hpp
@@ -293,7 +293,7 @@ class Bounds {
 // spiner-based EOS's.
 template <typename EOS>
 class SpinerTricks { // class instead of struct for -Wmismatched-tags
-public:
+ public:
   static auto GetOnDevice(EOS *peos_h) {
     // trivially copy all but dynamic memory
     EOS eos_d = *peos_h;

From 52a567616883164abced798436ba91fe72b6193b Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 14:25:43 -0500
Subject: [PATCH 10/16] formatting

---
 singularity-eos/eos/eos_gruneisen.hpp       | 2 +-
 singularity-eos/eos/eos_spiner_rho_temp.hpp | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp
index 8b79e79978a..b4eb41745f8 100644
--- a/singularity-eos/eos/eos_gruneisen.hpp
+++ b/singularity-eos/eos/eos_gruneisen.hpp
@@ -292,7 +292,7 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::ComputeRhoMax(const Real s1, const Real
           maxbound = std::min(min_extremum, maxbound);
         }
       } // s3 > 0 ; else
-    }   // discriminant >= 0
+    } // discriminant >= 0
     if (poly(minbound) * poly(maxbound) < 0.) {
       // Root is appropriately bounded
       using RootFinding1D::regula_falsi;
diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index 4fa9612c441..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<SpinerEOSDependsRhoT> {
                              reproducibility_mode, pmin_vapor_dome) {}
   PORTABLE_INLINE_FUNCTION
   SpinerEOSDependsRhoT()
-    : split_(TableSplit::Total), memoryStatus_(DataStatus::Deallocated) {}
+      : split_(TableSplit::Total), memoryStatus_(DataStatus::Deallocated) {}
 
   inline SpinerEOSDependsRhoT GetOnDevice();
 

From 1bc7475fe268194c2b121ef37d85b1fc075c5274 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 16:02:44 -0500
Subject: [PATCH 11/16] undo formatting change. apparently our version of clang
 format likes this one

---
 singularity-eos/eos/eos_gruneisen.hpp | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/singularity-eos/eos/eos_gruneisen.hpp b/singularity-eos/eos/eos_gruneisen.hpp
index b4eb41745f8..8b79e79978a 100644
--- a/singularity-eos/eos/eos_gruneisen.hpp
+++ b/singularity-eos/eos/eos_gruneisen.hpp
@@ -292,7 +292,7 @@ PORTABLE_INLINE_FUNCTION Real Gruneisen::ComputeRhoMax(const Real s1, const Real
           maxbound = std::min(min_extremum, maxbound);
         }
       } // s3 > 0 ; else
-    } // discriminant >= 0
+    }   // discriminant >= 0
     if (poly(minbound) * poly(maxbound) < 0.) {
       // Root is appropriately bounded
       using RootFinding1D::regula_falsi;

From a9d9f58612a0e0338c9842cfa4ecdb2279b2d5c9 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 16:03:37 -0500
Subject: [PATCH 12/16] turn on warnings in clang build

---
 .github/workflows/tests_clang.yml | 1 +
 1 file changed, 1 insertion(+)

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

From ac5c1b4ceac7737164e506f32d2d1d1d5a9afc51 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 16:05:01 -0500
Subject: [PATCH 13/16] changelog

---
 CHANGELOG.md | 1 +
 1 file changed, 1 insertion(+)

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

From d79b3c2a9e0bdb0bb813a47651fd16ae374ba56f Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 16:09:22 -0500
Subject: [PATCH 14/16] rho max present, but inside NDEBUG

---
 singularity-eos/closure/mixed_cell_models.hpp | 6 ++++--
 1 file changed, 4 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp
index 61a432b1ff3..042c81ca1a8 100644
--- a/singularity-eos/closure/mixed_cell_models.hpp
+++ b/singularity-eos/closure/mixed_cell_models.hpp
@@ -495,9 +495,11 @@ class PTESolverBase {
       bad_vfrac_guess = false;
       for (std::size_t m = 0; m < nmat; ++m) {
         const Real rho_min = eos[m].RhoPmin(Tguess);
+#ifndef NDEBUG
+        const Real rho_max = eos[m].MaximumDensity();
         PORTABLE_REQUIRE(rho_min < rho_max, "Valid density range must exist!");
-        PORTABLE_REQUIRE(rho[m] < eos[m].MaximumDensity(),
-                         "Density must be less than rho_max");
+        PORTABLE_REQUIRE(rho[m] < rho_max, "Density must be less than rho_max");
+#endif // NDEBUG
         if (rho[m] < rho_min) {
           bad_vfrac_guess = true;
           Tguess *= Tfactor;

From 652ece8ef0663f68db0cac20b8c9f0229c21e78d Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 16:13:38 -0500
Subject: [PATCH 15/16] explanatory comment

---
 singularity-eos/closure/mixed_cell_models.hpp | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp
index 042c81ca1a8..3bc7ab40b76 100644
--- a/singularity-eos/closure/mixed_cell_models.hpp
+++ b/singularity-eos/closure/mixed_cell_models.hpp
@@ -495,6 +495,8 @@ class PTESolverBase {
       bad_vfrac_guess = false;
       for (std::size_t m = 0; m < nmat; ++m) {
         const Real rho_min = eos[m].RhoPmin(Tguess);
+        // rho max only active inside the requires, which are only active
+        // when asserts are
 #ifndef NDEBUG
         const Real rho_max = eos[m].MaximumDensity();
         PORTABLE_REQUIRE(rho_min < rho_max, "Valid density range must exist!");

From 2a3576f2e081779ea40f934e13c45e09e023cc76 Mon Sep 17 00:00:00 2001
From: Jonah Miller <jonahm@lanl.gov>
Date: Fri, 12 Dec 2025 16:14:44 -0500
Subject: [PATCH 16/16] rho_max now marked maybe_unused

---
 singularity-eos/closure/mixed_cell_models.hpp | 6 +-----
 1 file changed, 1 insertion(+), 5 deletions(-)

diff --git a/singularity-eos/closure/mixed_cell_models.hpp b/singularity-eos/closure/mixed_cell_models.hpp
index 3bc7ab40b76..b4d6b6d0c83 100644
--- a/singularity-eos/closure/mixed_cell_models.hpp
+++ b/singularity-eos/closure/mixed_cell_models.hpp
@@ -495,13 +495,9 @@ class PTESolverBase {
       bad_vfrac_guess = false;
       for (std::size_t m = 0; m < nmat; ++m) {
         const Real rho_min = eos[m].RhoPmin(Tguess);
-        // rho max only active inside the requires, which are only active
-        // when asserts are
-#ifndef NDEBUG
-        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");
-#endif // NDEBUG
         if (rho[m] < rho_min) {
           bad_vfrac_guess = true;
           Tguess *= Tfactor;