Create a multi-EOS framework

CMakeLists.txt

@@ -56,12 +56,12 @@ cmake_dependent_option(SINGULARITY_BUILD_FORTRAN_BACKEND
 
 option(SINGULARITY_USE_KOKKOS "Use Kokkos for portability" OFF)
 option(SINGULARITY_USE_EOSPAC "Enable eospac backend" OFF)
-option(SINGULARITY_EOSPAC_ENABLE_SHMEM
+cmake_dependent_option(SINGULARITY_EOSPAC_ENABLE_SHMEM
   "Support shared memory with EOSPAC backend. Requires EOSPAC version 6.5.7."
-  OFF)
-option(SINGULARITY_EOSPAC_USE_MODERN_ERROR_CODES
+  OFF "SINGULARITY_USE_EOSPAC" OFF)
+cmake_dependent_option(SINGULARITY_EOSPAC_USE_MODERN_ERROR_CODES
   "Use the current eospac implementation of error codes, which must be compared with a function. Requires EOSPAC version 6.3.0 or greater."
-  ON)
+  ON "SINGULARITY_USE_EOSPAC" OFF)
 
 # TODO This should be dependent option (or fortran option)
 option(SINGULARITY_BUILD_CLOSURE "Mixed Cell Closure" ON)

singularity-eos/CMakeLists.txt

@@ -25,6 +25,7 @@ register_headers(
 
     # Normal files
     base/fast-math/logs.hpp
+    base/generic_indexer.hpp
     base/indexable_types.hpp
     base/robust_utils.hpp
     base/root-finding-1d/root_finding.hpp
@@ -38,6 +39,7 @@ register_headers(
     base/sp5/singularity_eos_sp5.hpp
     eos/default_variant.hpp
     base/hermite.hpp
+    base/tuple_utils.hpp
     eos/eos_variant.hpp
     eos/eos_stellar_collapse.hpp
     eos/eos_ideal.hpp
@@ -72,6 +74,7 @@ register_headers(
     eos/variant_utils.hpp
     eos/eos_carnahan_starling.hpp
     eos/eos_electrons.hpp
+    eos/eos_multi_eos.hpp
 )
 
 if (SINGULARITY_BUILD_CLOSURE)

singularity-eos/base/generic_indexer.hpp

@@ -0,0 +1,48 @@
+//------------------------------------------------------------------------------
+// © 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.
+//------------------------------------------------------------------------------
+
+#ifndef SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_
+#define SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_
+
+#include <utility>
+
+namespace singularity {
+
+// Index into an arbitrary array with an arbitrary map
+template <typename arrT, typename mapT>
+struct GenericIndexer {
+  arrT arr_;
+  mapT map_;
+
+  template <typename arrT_, typename mapT_>
+  constexpr GenericIndexer(arrT_ &&arr_in, mapT_ &&map_in)
+      : arr_{std::forward<arrT_>(arr_in)}, map_{std::forward<mapT_>(map_in)} {}
+
+  template <typename idxT>
+  constexpr auto &operator[](const idxT i) {
+    return arr_[map_[i]];
+  }
+
+  template <typename idxT>
+  constexpr const auto &operator[](const idxT i) const {
+    return arr_[map_[i]];
+  }
+};
+// CTAD
+template <typename arrT_, typename mapT_>
+GenericIndexer(arrT_ &&arr_in, mapT_ &&map_in) -> GenericIndexer<arrT_, mapT_>;
+
+} // namespace singularity
+
+#endif // #ifndef SINGULARITY_EOS_BASE_GENERIC_INDEXER_HPP_

singularity-eos/base/indexable_types.hpp

@@ -22,6 +22,21 @@
 #include <ports-of-call/portability.hpp>
 #include <singularity-eos/base/variadic_utils.hpp>
 
+// SFINAE helper macro that checks if a given indexer object has the requested
+// indexable type.
+// TODO: The error is a simple template substitution failure error, but there
+// are ways short of C++20 concepts to add more readability to the error. We may
+// want to add this at some point
+#define SINGULARITY_INDEXER_HAS_INDEXABLE_TYPE(Indexer, IndexableType)                   \
+  typename = std::enable_if_t<                                                           \
+      singularity::variadic_utils::is_indexable_v<Indexer, IndexableType>>
+
+// Just a small wrapper for the above macro specifically for mass fractions
+// NOTE: it's assumed that matnum is a count, not an index (hence the minus 1)
+#define SINGULARITY_INDEXER_HAS_MASS_FRAC(Indexer, matnum)                               \
+  SINGULARITY_INDEXER_HAS_INDEXABLE_TYPE(                                                \
+      Indexer, singularity::IndexableTypes::MassFraction<matnum - 1>)
+
 namespace singularity {
 namespace IndexerUtils {
 // Convenience function for accessing an indexer by either type or
@@ -40,25 +55,22 @@ PORTABLE_FORCEINLINE_FUNCTION auto &Get(Indexer_t &&lambda, std::size_t idx = 0)
 template <typename Data_t, typename... Ts>
 class VariadicIndexerBase {
  public:
-  VariadicIndexerBase() = default;
-  PORTABLE_FORCEINLINE_FUNCTION
-  VariadicIndexerBase(const Data_t &data) : data_(data) {}
+  constexpr VariadicIndexerBase() = default;
+  constexpr VariadicIndexerBase(const Data_t &data) : data_(data) {}
   template <typename T,
             typename = std::enable_if_t<variadic_utils::contains<T, Ts...>::value>>
-  PORTABLE_FORCEINLINE_FUNCTION Real &operator[](const T &t) {
+  constexpr Real &operator[](const T &t) {
     constexpr std::size_t idx = variadic_utils::GetIndexInTL<T, Ts...>();
     return data_[idx];
   }
-  PORTABLE_FORCEINLINE_FUNCTION
-  Real &operator[](const std::size_t idx) { return data_[idx]; }
+  constexpr Real &operator[](const std::size_t idx) { return data_[idx]; }
   template <typename T,
             typename = std::enable_if_t<variadic_utils::contains<T, Ts...>::value>>
-  PORTABLE_FORCEINLINE_FUNCTION const Real &operator[](const T &t) const {
+  constexpr const Real &operator[](const T &t) const {
     constexpr std::size_t idx = variadic_utils::GetIndexInTL<T, Ts...>();
     return data_[idx];
   }
-  PORTABLE_FORCEINLINE_FUNCTION
-  const Real &operator[](const std::size_t idx) const { return data_[idx]; }
+  constexpr const Real &operator[](const std::size_t idx) const { return data_[idx]; }
   static inline constexpr std::size_t size() { return sizeof...(Ts); }
 
  private:
@@ -81,6 +93,8 @@ struct MeanAtomicNumber {};
 struct ElectronFraction {};
 struct RootStatus {};
 struct TableStatus {};
+template <int mat_idx>
+struct MassFraction {};
 } // namespace IndexableTypes
 } // namespace singularity
 #endif // SINGULARITY_EOS_BASE_INDEXABLE_TYPES_

singularity-eos/base/math_utils.hpp

@@ -15,6 +15,8 @@
 #ifndef SINGULARITY_EOS_BASE_MATH_UTILS_HPP_
 #define SINGULARITY_EOS_BASE_MATH_UTILS_HPP_
 
+#include <cmath>
+
 #include <ports-of-call/portability.hpp>
 
 namespace singularity {
@@ -58,6 +60,33 @@ PORTABLE_FORCEINLINE_FUNCTION auto pow10(const Real x) {
   return std::exp(ln10 * x);
 }
 
+/*
+ * Kahan summation, with the Neumaier correction
+ * https://onlinelibrary.wiley.com/doi/10.1002/zamm.19740540106
+ */
+struct IdentityOperator {
+  PORTABLE_FORCEINLINE_FUNCTION Real operator()(const Real x) const { return x; }
+};
+template <typename Data_t, typename Operator_t = IdentityOperator>
+PORTABLE_FORCEINLINE_FUNCTION Real
+sum_neumaier(Data_t &&data, std::size_t n, std::size_t offset = 0, std::size_t iskip = -1,
+             const Operator_t &op = IdentityOperator()) {
+  Real sum = 0;
+  Real c = 0; // correction
+  for (std::size_t i = 0; i < n; ++i) {
+    if (i == iskip) continue;
+    Real x = op(data[i + offset]);
+    Real t = sum + x;
+    if (std::abs(sum) >= std::abs(x)) {
+      c += (sum - t) + x;
+    } else {
+      c += (x - t) + sum;
+    }
+    sum = t;
+  }
+  return sum + c;
+}
+
 } // namespace math_utils
 } // namespace singularity
 

singularity-eos/base/tuple_utils.hpp

@@ -0,0 +1,47 @@
+//------------------------------------------------------------------------------
+// © 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.
+//------------------------------------------------------------------------------
+
+#ifndef SINGULARITY_EOS_BASE_TUPLE_UTILS_HPP_
+#define SINGULARITY_EOS_BASE_TUPLE_UTILS_HPP_
+
+#include <tuple>
+
+namespace singularity {
+namespace tuple_utils {
+
+namespace impl {
+
+// Helper that copies a tuple into a new tuple while also applying a
+// transformation function to each element. Wrapped in `impl` namespace due to
+// index sequence and called below with a more normal function signature
+template <typename TupleT, typename FunctionT, std::size_t... I>
+constexpr auto tuple_transform_iteration(TupleT &&tup, FunctionT &&function,
+                                         std::index_sequence<I...>) {
+  return std::make_tuple(function(std::get<I>(std::forward<TupleT>(tup)))...);
+}
+} // namespace impl
+
+// Returns a copy of the tuple with a transform function called on each element
+template <typename TupleT, typename FunctionT>
+constexpr auto tuple_transform(TupleT &&tup, FunctionT &&function) {
+  constexpr std::size_t N = std::tuple_size_v<std::remove_reference_t<TupleT>>;
+  return impl::tuple_transform_iteration(std::forward<TupleT>(tup),
+                                         std::forward<FunctionT>(function),
+                                         std::make_index_sequence<N>{});
+}
+
+} // namespace tuple_utils
+} // namespace singularity
+
+#endif // #ifndef SINGULARITY_EOS_BASE_TUPLE_UTILS_HPP_

singularity-eos/closure/mixed_cell_models.hpp

@@ -18,6 +18,7 @@
 #include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_errors.hpp>
 #include <singularity-eos/base/fast-math/logs.hpp>
+#include <singularity-eos/base/math_utils.hpp>
 #include <singularity-eos/base/robust_utils.hpp>
 #include <singularity-eos/base/variadic_utils.hpp>
 #include <singularity-eos/eos/eos.hpp>
@@ -213,33 +214,6 @@ bool solve_Ax_b_wscr(const std::size_t n, Real *a, Real *b, Real *scr) {
   return retval;
 }
 
-/*
- * Kahan summation, with the Neumaier correction
- * https://onlinelibrary.wiley.com/doi/10.1002/zamm.19740540106
- */
-struct IdentityOperator {
-  PORTABLE_FORCEINLINE_FUNCTION Real operator()(const Real x) const { return x; }
-};
-template <typename Data_t, typename Operator_t = IdentityOperator>
-PORTABLE_FORCEINLINE_FUNCTION Real
-sum_neumaier(Data_t &&data, std::size_t n, std::size_t offset = 0, std::size_t iskip = -1,
-             const Operator_t &op = IdentityOperator()) {
-  Real sum = 0;
-  Real c = 0; // correction
-  for (std::size_t i = 0; i < n; ++i) {
-    if (i == iskip) continue;
-    Real x = op(data[i + offset]);
-    Real t = sum + x;
-    if (std::abs(sum) >= std::abs(x)) {
-      c += (sum - t) + x;
-    } else {
-      c += (x - t) + sum;
-    }
-    sum = t;
-  }
-  return sum + c;
-}
-
 class CacheAccessor {
  public:
   CacheAccessor() = default;
@@ -323,12 +297,12 @@ class PTESolverBase {
       PORTABLE_REQUIRE(rho[m] > 0., "Non-positive density provided to PTE solver");
       rhobar[m] = rho[m] * vfrac[m];
     }
-    rho_total = sum_neumaier(rhobar, nmat);
+    rho_total = math_utils::sum_neumaier(rhobar, nmat);
   }
 
   PORTABLE_FORCEINLINE_FUNCTION
-  void NormalizeVfrac() const {
-    Real vfrac_sum = sum_neumaier(vfrac, nmat);
+  void NormalizeVfrac() {
+    Real vfrac_sum = math_utils::sum_neumaier(vfrac, nmat);
     for (std::size_t m = 0; m < nmat; ++m) {
       vfrac[m] *= robust::ratio(vfrac_total, vfrac_sum);
     }
@@ -559,7 +533,8 @@ class PTESolverBase {
 
   PORTABLE_INLINE_FUNCTION
   Real ResidualNorm() const {
-    return 0.5 * sum_neumaier(residual, neq, 0, -1, [](const Real x) { return x * x; });
+    return 0.5 * math_utils::sum_neumaier(residual, neq, 0, -1,
+                                          [](const Real x) { return x * x; });
   }
 
   PORTABLE_FORCEINLINE_FUNCTION
@@ -900,7 +875,7 @@ class PTESolverRhoT
 
   PORTABLE_INLINE_FUNCTION
   void Residual() const {
-    Real esum = mix_impl::sum_neumaier(u, nmat);
+    Real esum = math_utils::sum_neumaier(u, nmat);
     residual[0] = utotal_scale - esum;
     std::size_t ires = 1;
     for (std::size_t m = 0; m < nmat; ++m) {
@@ -916,8 +891,8 @@ class PTESolverRhoT
     for (std::size_t m = 0; m < nmat; ++m) {
       mean_p += vfrac[m] * press[m];
     }
-    Real error_p = std::sqrt(
-        sum_neumaier(residual, neq - 1, 1, -1, [](const Real x) { return x * x; }));
+    Real error_p = std::sqrt(math_utils::sum_neumaier(
+        residual, neq - 1, 1, -1, [](const Real x) { return x * x; }));
     Real error_u = std::abs(residual[0]);
     // Check for convergence
     bool converged_p = (error_p < params_.pte_rel_tolerance_p * std::abs(mean_p) ||
@@ -1002,7 +977,7 @@ class PTESolverRhoT
     // control how big of a step toward vfrac = 0 is allowed
     // 0th material is special as delta volume fraction for it is
     // minus the sum of the deltas for the others
-    Real dalphaskip = -mix_impl::sum_neumaier(dx, nmat - 1, 1);
+    Real dalphaskip = -math_utils::sum_neumaier(dx, nmat - 1, 1);
     std::size_t idx = 1;
     for (std::size_t m = 0; m < nmat; ++m) {
       Real mydx = (m == ms) ? dalphaskip : dx[idx++];
@@ -1058,7 +1033,7 @@ class PTESolverRhoT
       if (ms == m) continue;
       vfrac[m] = vtemp[m] + scale * dx[idx++];
     }
-    vfrac[ms] = mix_impl::sum_neumaier(vfrac, nmat, 0, ms);
+    vfrac[ms] = math_utils::sum_neumaier(vfrac, nmat, 0, ms);
     vfrac[ms] = 1 - vfrac[ms];
     for (std::size_t m = 0; m < nmat; ++m) {
       rho[m] = robust::ratio(rhobar[m], vfrac[m]);
@@ -1165,7 +1140,7 @@ class PTESolverPT
       // always approach from >0 side
       Pequil += std::abs(press[m]) * vfrac[m];
     }
-    Real vsum = mix_impl::sum_neumaier(vfrac, nmat);
+    Real vsum = math_utils::sum_neumaier(vfrac, nmat);
     Pequil /= vsum;
     Tequil = 1; // Because it's = Tnorm = initial guess
 
@@ -1201,8 +1176,8 @@ class PTESolverPT
 
   PORTABLE_INLINE_FUNCTION
   void Residual() const {
-    Real vsum = mix_impl::sum_neumaier(vfrac, nmat);
-    Real esum = mix_impl::sum_neumaier(u, nmat);
+    Real vsum = math_utils::sum_neumaier(vfrac, nmat);
+    Real esum = math_utils::sum_neumaier(u, nmat);
     residual[RV] = vfrac_total - vsum;
     residual[RSIE] = utotal_scale - esum;
   }
@@ -1433,7 +1408,7 @@ class PTESolverFixedT
 
   PORTABLE_INLINE_FUNCTION
   void Residual() const {
-    Real vsum = mix_impl::sum_neumaier(vfrac, nmat);
+    Real vsum = math_utils::sum_neumaier(vfrac, nmat);
     residual[0] = vfrac_total - vsum;
     for (std::size_t m = 0; m < nmat - 1; ++m) {
       residual[1 + m] = press[m] - press[m + 1];
@@ -1657,7 +1632,7 @@ class PTESolverFixedP
 
   PORTABLE_INLINE_FUNCTION
   void Residual() const {
-    Real vsum = mix_impl::sum_neumaier(vfrac, nmat);
+    Real vsum = math_utils::sum_neumaier(vfrac, nmat);
     for (std::size_t m = 0; m < nmat; ++m) {
       residual[m] = robust::ratio(Pequil, uscale) - press[m];
     }
@@ -1878,8 +1853,8 @@ class PTESolverRhoU
 
   PORTABLE_INLINE_FUNCTION
   void Residual() const {
-    Real vsum = mix_impl::sum_neumaier(vfrac, nmat);
-    Real esum = mix_impl::sum_neumaier(u, nmat);
+    Real vsum = math_utils::sum_neumaier(vfrac, nmat);
+    Real esum = math_utils::sum_neumaier(u, nmat);
     residual[0] = vfrac_total - vsum;
     residual[1] = utotal_scale - esum;
     for (std::size_t m = 0; m < nmat - 1; ++m) {