From 5a492fd730cbb40bb52b42287b0934204ac35296 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 11:43:34 -0600
Subject: [PATCH 01/30] Add ability to check if a type accepts an whole number
 index

---
 singularity-eos/base/variadic_utils.hpp | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index fc282e4c6c8..a96bd504942 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -111,6 +111,16 @@ struct is_indexable<T, Index,
 template <typename T, typename Index>
 constexpr bool is_indexable_v = is_indexable<T, Index>::value;
 
+//
+template<class T, class=void>
+struct has_whole_num_index : std::false_type {};
+template<class T>
+struct has_whole_num_index<T,
+                         std::void_t<decltype(std::declval<T>()[std::declval<std::size_t>()])>>
+    : std::true_type {};
+template <typename T>
+constexpr bool has_whole_num_index_v = has_whole_num_index<T>::value;
+
 // this flattens a typelist of typelists to a single typelist
 
 // first parameter - accumulator

From 02a7396a7143f025b79f0e41c05dce295fc80bd1 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 11:46:01 -0600
Subject: [PATCH 02/30] Add safeSet() and safeGet() helpers

---
 singularity-eos/base/indexable_types.hpp | 101 ++++++++++++++++-
 test/test_indexable_types.cpp            | 131 ++++++++++++++++++++++-
 2 files changed, 228 insertions(+), 4 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 25cc4928d95..c5ddac164b1 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -24,8 +24,85 @@
 
 namespace singularity {
 namespace IndexerUtils {
-// Convenience function for accessing an indexer by either type or
-// natural number index depending on what is available
+
+// Identifies an indexer as a type-based indexer
+template<class, class=void>
+struct is_type_indexer : std::false_type {};
+template<class Indexer_t>
+struct is_type_indexer<Indexer_t,
+  std::void_t<decltype(std::decay_t<Indexer_t>::is_type_indexable)>>
+  : std::bool_constant<std::decay_t<Indexer_t>::is_type_indexable> {};
+template<class Indexer_t>
+constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
+
+// The "safe" version of Get(). This function will ONLY return a value IF that
+// type-based index is present in the Indexer OR if the Indexer doesn't support
+// type-based indexing.
+template <typename T, typename Indexer_t>
+inline bool safeGet(Indexer_t const& lambda, std::size_t const idx, Real& out) noexcept {
+  // If null then nothing happens
+  if (variadic_utils::is_nullptr(lambda)) {
+    return false;
+  }
+
+  // Return value if type index is available
+  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
+    out = lambda[T{}];
+    return true;
+  }
+
+  // Do nothing if lambda has type indexing BUT doesn't have this type index
+  if constexpr (is_type_indexer_v<Indexer_t>) {
+    return false;
+  }
+
+  // Fall back to numerical indexing if no type indexing
+  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
+    out = lambda[idx];
+    return true;
+  }
+
+  // Something else...
+  return false;
+}
+
+// Break out "Set" functionality from "Get". The original "Get()" did both, but
+// the "safe" version needs to separate that functionality for setting the
+// values in a lambda
+template <typename T, typename Indexer_t>
+inline bool safeSet(Indexer_t& lambda, std::size_t const idx, Real const in) noexcept {
+  // If null then nothing happens
+  if (variadic_utils::is_nullptr(lambda)) {
+    return false;
+  }
+
+  // Return value if type index is available
+  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
+    lambda[T{}] = in;
+    return true;
+  }
+
+  // Do nothing if lambda has type indexing BUT doesn't have this type index
+  if constexpr (is_type_indexer_v<Indexer_t>) {
+    return false;
+  }
+
+  // Fall back to numerical indexing if no type indexing
+  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
+    lambda[idx] = in;
+    return true;
+  }
+
+  // Something else...
+  return false;
+}
+
+// NOTE: this Get is "unsafe" because it can allow you to overwrite a type-based
+// index since it automatically falls back to numeric indexing if the type
+// index isn't present.
+
+// Convenience function for accessing an indexer by either type or natural
+// number index depending on what is available.
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION auto &Get(Indexer_t &&lambda, std::size_t idx = 0) {
   if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
@@ -40,25 +117,44 @@ PORTABLE_FORCEINLINE_FUNCTION auto &Get(Indexer_t &&lambda, std::size_t idx = 0)
 template <typename Data_t, typename... Ts>
 class VariadicIndexerBase {
  public:
+  // Any class that wants to be recognized as indexable (so that we don't
+  // accidentally fall back to integer indexing when we don't want to) needs to
+  // include this.
+  constexpr static bool is_type_indexable = true;
+
+  // JHP: another option for the `is_type_indexable` flag is to take the ADL
+  // route. Essentially this would involve defining a friend function that
+  // could be defined in an appropriate namesapce so that theoretically a host
+  // code could use a TPL container with type-based indexing and allow that
+  // container to be flagged in our code as acceptable. This seems like a bit
+  // of a heavy hammer for what we need here though. We can easily change this
+  // if a TPL provides a type that is being used for this purpose.
+
   VariadicIndexerBase() = default;
+
   PORTABLE_FORCEINLINE_FUNCTION
   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 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]; }
+
   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 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]; }
+
   static inline constexpr std::size_t size() { return sizeof...(Ts); }
 
  private:
@@ -70,6 +166,7 @@ using VariadicIndexer = VariadicIndexerBase<std::array<Real, sizeof...(Ts)>, Ts.
 // uses a Real*
 template <typename... Ts>
 using VariadicPointerIndexer = VariadicIndexerBase<Real *, Ts...>;
+
 } // namespace IndexerUtils
 
 namespace IndexableTypes {
diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index d4b2231e731..40edc76c634 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -29,14 +29,22 @@ using Lambda_t = VariadicIndexer<MeanIonizationState, ElectronFraction, LogDensi
 
 class ManualLambda_t {
  public:
+  constexpr static size_t length = 3;
   ManualLambda_t() = default;
+  const Real &operator[](const std::size_t idx) const { return data_[idx]; }
   Real &operator[](const std::size_t idx) { return data_[idx]; }
   Real &operator[](const MeanIonizationState &zbar) { return data_[0]; }
   Real &operator[](const ElectronFraction &Ye) { return data_[1]; }
   Real &operator[](const LogDensity &lRho) { return data_[2]; }
 
  private:
-  std::array<Real, 3> data_;
+  std::array<Real, length> data_;
+};
+
+class NewManualLambda_t : public ManualLambda_t {
+ public:
+  // Enable recognition that this is type-indexable
+  constexpr static bool is_type_indexable = true;
 };
 
 SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer]") {
@@ -73,13 +81,73 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
         REQUIRE(lRho == static_cast<Real>(2));
       }
     }
+    WHEN("We use the safeGet functionality") {
+      constexpr Real unmodified = -1.0;
+      Real destination = unmodified;
+      WHEN("We request a type that exists") {
+        const bool modified = safeGet<MeanIonizationState>(lambda, 2, destination);
+        THEN("The destination value will be modified") {
+          CHECK(modified);
+          REQUIRE(destination == lambda[MeanIonizationState{}]);
+        }
+      }
+      WHEN("We request a type that doesn't exist") {
+        const bool modified = safeGet<TableStatus>(lambda, 2, destination);
+        THEN("The destination value will remain UNmodified") {
+          CHECK(!modified);
+          REQUIRE(destination == unmodified);
+        }
+      }
+      WHEN("A normal array-like lambda is used") {
+        std::array<Real, Lambda_t::size()> lambda_arr{1, 2, 3};
+        constexpr size_t my_index = 2;
+        const bool modified = safeGet<LogDensity>(lambda_arr, my_index, destination);
+        THEN("The destination value will reflect the index from the array") {
+          CHECK(modified);
+          REQUIRE(destination == lambda_arr[my_index]);
+        }
+      }
+    }
+    WHEN("We use the safeSet functionality") {
+      constexpr Real new_value = -1.0;
+      WHEN("We want to set a value for a type index that exists") {
+        const bool modified = safeSet<MeanIonizationState>(lambda, 2, new_value);
+        THEN("The lambda index was modified") {
+          CHECK(modified);
+          REQUIRE(lambda[MeanIonizationState{}] == new_value);
+        }
+      }
+      WHEN("We want to set a value for a type index that doesn't exist") {
+        const bool modified = safeSet<TableStatus>(lambda, 2, new_value);
+        Lambda_t old_lambda;
+        for (std::size_t i = 0; i < Lambda_t::size(); ++i) {
+          old_lambda[i] = lambda[i];
+        }
+        THEN("None of the lambda values were modified") {
+          CHECK(!modified);
+          for (std::size_t i = 0; i < Lambda_t::size(); ++i) {
+            INFO("i: " << i);
+            CHECK(lambda[i] == old_lambda[i]);
+          }
+        }
+      }
+      WHEN("A normal array-like lambda is used") {
+        std::array<Real, Lambda_t::size()> lambda_arr{4, 5, 6};
+        constexpr size_t my_index = 1;
+        const bool modified = safeSet<LogDensity>(lambda_arr, my_index, new_value);
+        THEN("The lambda value at the appropriate index has been modified") {
+          CHECK(modified);
+          REQUIRE(lambda_arr[my_index] == new_value);
+        }
+      }
+    }
   }
 }
 
 SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
   GIVEN("A manually written indexer, filled with indices 0, 1, 2") {
     ManualLambda_t lambda;
-    for (std::size_t i = 0; i < 3; ++i) {
+    for (std::size_t i = 0; i < lambda.length; ++i) {
       lambda[i] = static_cast<Real>(i);
     }
     WHEN("We use the Get functionality") {
@@ -92,5 +160,64 @@ SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
         REQUIRE(lRho == static_cast<Real>(2));
       }
     }
+    WHEN("We use the safeGet functionality") {
+      constexpr Real unmodified = -1.0;
+      Real destination = unmodified;
+      WHEN("We request a type that doesn't exist in the manual indexer") {
+        constexpr size_t my_index = 1;
+        const bool modified = safeGet<TableStatus>(lambda, my_index, destination);
+        THEN("The destination WILL be modified since the manual indexer doesn't have the "
+             " `is_type_indexable` data member") {
+          CHECK(modified);
+          REQUIRE(destination == lambda[my_index]);
+        }
+      }
+      WHEN("We define a new manual indexer that has the `is_type_indexable` "
+           "data member") {
+        NewManualLambda_t lambda_new;
+        for (std::size_t i = 0; i < lambda.length; ++i) {
+          lambda_new[i] = lambda[i];
+        }
+        WHEN("We request a type that doesn't exist in the manual indexer") {
+          destination = unmodified;
+          constexpr size_t my_index = 1;
+          const bool modified = safeGet<TableStatus>(lambda_new, my_index, destination);
+          THEN("The destination will NOT be modified") {
+            CHECK(!modified);
+            REQUIRE(destination == unmodified);
+          }
+        }
+      }
+    }
+    WHEN("We use the safeSet functionality") {
+      constexpr Real new_value = -1.0;
+      WHEN("We want to set a value for a type index that doesn't exist") {
+        constexpr size_t my_index = 1;
+        const bool modified = safeSet<TableStatus>(lambda, my_index, new_value);
+        THEN("The lambda value WILL be modified since the manual indexer doesn't have "
+             "the `is_type_indexable` data member") {
+          CHECK(modified);
+          REQUIRE(lambda[my_index] == new_value);
+        }
+      }
+      WHEN("We define a new manual indexer that has the `is_type_indexable` data "
+           "member") {
+        NewManualLambda_t lambda_new;
+        for (std::size_t i = 0; i < lambda.length; ++i) {
+          lambda_new[i] = lambda[i];
+        }
+        WHEN("We request a type that doesn't exist in the manual indexer") {
+          constexpr size_t my_index = 1;
+          const bool modified = safeSet<TableStatus>(lambda_new, my_index, new_value);
+          THEN("The lambda will NOT be modified") {
+            CHECK(!modified);
+            for (std::size_t i = 0; i < lambda.length; ++i) {
+              INFO("i: " << i);
+              CHECK(lambda_new[i] == lambda[i]);
+            }
+          }
+        }
+      }
+    }
   }
 }

From 8895d8062cad94fd9b27a37e1e3596597266bbcd Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 12:25:06 -0600
Subject: [PATCH 03/30] Add safeGet and safeSet that don't take numerical
 indices

---
 singularity-eos/base/indexable_types.hpp | 50 +++++++++++++++++++++++-
 1 file changed, 48 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index c5ddac164b1..6e25f84f6ba 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -39,7 +39,7 @@ constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
 // type-based index is present in the Indexer OR if the Indexer doesn't support
 // type-based indexing.
 template <typename T, typename Indexer_t>
-inline bool safeGet(Indexer_t const& lambda, std::size_t const idx, Real& out) noexcept {
+inline bool safeGet(Indexer_t const& lambda, std::size_t const idx, Real& out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -66,11 +66,34 @@ inline bool safeGet(Indexer_t const& lambda, std::size_t const idx, Real& out) n
   return false;
 }
 
+// Overload when no index is provided
+template <typename T, typename Indexer_t>
+inline bool safeGet(Indexer_t const& lambda, Real& out) {
+  // If null then nothing happens
+  if (variadic_utils::is_nullptr(lambda)) {
+    return false;
+  }
+
+  // Return value if type index is available
+  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
+    out = lambda[T{}];
+    return true;
+  }
+
+  // Do nothing if lambda has type indexing BUT doesn't have this type index
+  if constexpr (is_type_indexer_v<Indexer_t>) {
+    return false;
+  }
+
+  // Something else...
+  return false;
+}
+
 // Break out "Set" functionality from "Get". The original "Get()" did both, but
 // the "safe" version needs to separate that functionality for setting the
 // values in a lambda
 template <typename T, typename Indexer_t>
-inline bool safeSet(Indexer_t& lambda, std::size_t const idx, Real const in) noexcept {
+inline bool safeSet(Indexer_t& lambda, std::size_t const idx, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -97,6 +120,29 @@ inline bool safeSet(Indexer_t& lambda, std::size_t const idx, Real const in) noe
   return false;
 }
 
+// Overload without numeric index
+template <typename T, typename Indexer_t>
+inline bool safeSet(Indexer_t& lambda, Real const in) {
+  // If null then nothing happens
+  if (variadic_utils::is_nullptr(lambda)) {
+    return false;
+  }
+
+  // Return value if type index is available
+  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
+    lambda[T{}] = in;
+    return true;
+  }
+
+  // Do nothing if lambda has type indexing BUT doesn't have this type index
+  if constexpr (is_type_indexer_v<Indexer_t>) {
+    return false;
+  }
+
+  // Something else...
+  return false;
+}
+
 // NOTE: this Get is "unsafe" because it can allow you to overwrite a type-based
 // index since it automatically falls back to numeric indexing if the type
 // index isn't present.

From 58a920bd7ce03a3375e697a7b583f502a1f2a023 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 12:25:42 -0600
Subject: [PATCH 04/30] Switch to safeGet and safeSet for indexable types

---
 singularity-eos/eos/eos_spiner_rho_sie.hpp  | 30 +++----
 singularity-eos/eos/eos_spiner_rho_temp.hpp | 88 ++++++++-------------
 2 files changed, 40 insertions(+), 78 deletions(-)

diff --git a/singularity-eos/eos/eos_spiner_rho_sie.hpp b/singularity-eos/eos/eos_spiner_rho_sie.hpp
index dd54c5b22e1..44db0c76f0d 100644
--- a/singularity-eos/eos/eos_spiner_rho_sie.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_sie.hpp
@@ -687,9 +687,7 @@ PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoSieTransformable<TransformerT>:
     }
   } else {
     lRho = to_log(rho, lRhoOffset_);
-    if (!variadic_utils::is_nullptr(lambda)) {
-      IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    }
+    IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
   }
   if (output & thermalqs::temperature) {
 
@@ -747,9 +745,7 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::interpRhoT_(const Real rho,
                                                                Indexer_t &&lambda) const {
   const Real lRho = spiner_common::to_log(rho, lRhoOffset_);
   const Real lT = spiner_common::to_log(T, lTOffset_);
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
   return db.interpToReal(lRho, lT);
 }
 
@@ -761,9 +757,7 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::interpRhoSie_(
   const Real lRho = spiner_common::to_log(rho, lRhoOffset_);
   const Real sie_transformed = transformer_.transform(sie, rho);
   const Real lE = spiner_common::to_log(sie_transformed, lEOffset_);
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
   return db.interpToReal(lRho, lE);
 }
 
@@ -786,12 +780,10 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::lRhoFromPlT_(
     }
   } else {
     Real lRhoGuess = reproducible_ ? lRhoMin_ : 0.5 * (lRhoMin_ + lRhoMax_);
-    if (!variadic_utils::is_nullptr(lambda)) {
-      Real lRho_cache =
-          IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho);
-      if ((lRhoMin_ <= lRho_cache) && (lRho_cache <= lRhoMax_)) {
-        lRhoGuess = lRho_cache;
-      }
+    Real lRho_cache;
+    IndexerUtils::safeGet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho_cache);
+    if ((lRhoMin_ <= lRho_cache) && (lRho_cache <= lRhoMax_)) {
+      lRhoGuess = lRho_cache;
     }
     const callable_interp::l_interp PFunc(dependsRhoT_.P, lT);
     status = SP_ROOT_FINDER(PFunc, P, lRhoGuess, lRhoMin_, lRhoMax_, robust::EPS(),
@@ -809,12 +801,8 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::lRhoFromPlT_(
       lRho = reproducible_ ? lRhoMin_ : lRhoGuess;
     }
   }
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t, IndexableTypes::RootStatus>) {
-      lambda[IndexableTypes::RootStatus()] = static_cast<Real>(status);
-    }
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
   return lRho;
 }
 
diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index d6bf1b9d408..947e70e574b 100644
--- a/singularity-eos/eos/eos_spiner_rho_temp.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp
@@ -880,10 +880,8 @@ SpinerEOSDependsRhoT::FillEos(Real &rho, Real &temp, Real &energy, Real &press,
   if (output & thermalqs::bulk_modulus) {
     bmod = bModFromRholRhoTlT_(rho, lRho, temp, lT, whereAmI);
   }
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
 template <typename Indexer_t>
@@ -914,10 +912,8 @@ SpinerEOSDependsRhoT::getLogsRhoT_(const Real rho, const Real temperature, Real
                                    Real &lT, Indexer_t &&lambda) const {
   lRho = lRho_(rho);
   lT = lT_(temperature);
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
 template <typename Indexer_t>
@@ -931,11 +927,10 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lRhoFromPlT_(
   const RootFinding1D::RootCounts *pcounts =
       (memoryStatus_ == DataStatus::OnDevice) ? nullptr : &counts;
 
-  if (!variadic_utils::is_nullptr(lambda)) {
-    Real lRho_cache = IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho);
-    if ((lRhoMin_ <= lRho_cache) && (lRho_cache <= lRhoMax_)) {
-      lRhoGuess = lRho_cache;
-    }
+  Real lRho_cache;
+  IndexerUtils::safeGet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho_cache);
+  if ((lRhoMin_ <= lRho_cache) && (lRho_cache <= lRhoMax_)) {
+    lRhoGuess = lRho_cache;
   }
 
   if (lT <= lTMin_) { // cold curve
@@ -972,17 +967,11 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lRhoFromPlT_(
 #endif // SPINER_EOS_VERBOSE
     lRho = reproducible_ ? lRhoMax_ : lRhoGuess;
   }
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t, IndexableTypes::RootStatus>) {
-      lambda[IndexableTypes::RootStatus()] = static_cast<Real>(status);
-    }
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t,
-                                                 IndexableTypes::TableStatus>) {
-      lambda[IndexableTypes::TableStatus()] = static_cast<Real>(whereAmI);
-    }
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  // No numerical index: only set if type-based indexing is used
+  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::safeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
   return lRho;
 }
 
@@ -1013,12 +1002,10 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoSie_(
     }
   } else {
     Real lTGuess = reproducible_ ? lTMin_ : 0.5 * (lTMin_ + lTMax_);
-    if (!variadic_utils::is_nullptr(lambda)) {
-      Real lT_cache =
-          IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT);
-      if ((lTMin_ <= lT_cache) && (lT_cache <= lTMax_)) {
-        lTGuess = lT_cache;
-      }
+    Real lT_cache;
+    IndexerUtils::safeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT_cache);
+    if ((lTMin_ <= lT_cache) && (lT_cache <= lTMax_)) {
+      lTGuess = lT_cache;
     }
     const callable_interp::r_interp sieFunc(sie_, lRho);
     status = SP_ROOT_FINDER(sieFunc, sie, lTGuess, lTMin_, lTMax_, ROOT_THRESH,
@@ -1040,17 +1027,11 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoSie_(
       lT = reproducible_ ? lTMin_ : lTGuess;
     }
   }
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t, IndexableTypes::RootStatus>) {
-      lambda[IndexableTypes::RootStatus()] = static_cast<Real>(status);
-    }
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t,
-                                                 IndexableTypes::TableStatus>) {
-      lambda[IndexableTypes::TableStatus()] = static_cast<Real>(whereAmI);
-    }
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  // No numerical index: only set if type-based indexing is used
+  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::safeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
   return lT;
 }
 
@@ -1080,11 +1061,10 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoP_(
   } else {
     whereAmI = TableStatus::OnTable;
     lTGuess = 0.5 * (lTMin_ + lTMax_);
-    if (!variadic_utils::is_nullptr(lambda)) {
-      Real lT_cache = IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, lT);
-      if ((lTMin_ <= lT_cache) && (lT_cache <= lTMax_)) {
-        lTGuess = lT_cache;
-      }
+    Real lT_cache;
+    IndexerUtils::safeGet<IndexableTypes::LogTemperature>(lambda, lT, lT_cache);
+    if ((lTMin_ <= lT_cache) && (lT_cache <= lTMax_)) {
+      lTGuess = lT_cache;
     }
     const callable_interp::r_interp PFunc(P_, lRho);
     status = SP_ROOT_FINDER(PFunc, press, lTGuess, lTMin_, lTMax_, ROOT_THRESH,
@@ -1102,17 +1082,11 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoP_(
       lT = reproducible_ ? lTMin_ : lTGuess;
     }
   }
-  if (!variadic_utils::is_nullptr(lambda)) {
-    IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, Lambda::lRho) = lRho;
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t, IndexableTypes::RootStatus>) {
-      lambda[IndexableTypes::RootStatus()] = static_cast<Real>(status);
-    }
-    if constexpr (variadic_utils::is_indexable_v<Indexer_t,
-                                                 IndexableTypes::TableStatus>) {
-      lambda[IndexableTypes::TableStatus()] = static_cast<Real>(whereAmI);
-    }
-  }
+  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  // No numerical index: only set if type-based indexing is used
+  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::safeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
   return lT;
 }
 

From fce8f7ce1f247b8bff4687c62f4ba1d2b70d26f9 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 13:02:21 -0600
Subject: [PATCH 05/30] Whoops... forgot comment

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

diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index a96bd504942..a1d7fc2d669 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -111,7 +111,7 @@ struct is_indexable<T, Index,
 template <typename T, typename Index>
 constexpr bool is_indexable_v = is_indexable<T, Index>::value;
 
-//
+// Check if a type can accept a size_t index
 template<class T, class=void>
 struct has_whole_num_index : std::false_type {};
 template<class T>

From 4959500b7ccf635bc63f29079bec0a81f26db736 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 13:02:50 -0600
Subject: [PATCH 06/30] Clang format

---
 singularity-eos/base/indexable_types.hpp | 18 +++++++++---------
 singularity-eos/base/variadic_utils.hpp  |  8 ++++----
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 6e25f84f6ba..d92bb9af5bc 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -26,20 +26,20 @@ namespace singularity {
 namespace IndexerUtils {
 
 // Identifies an indexer as a type-based indexer
-template<class, class=void>
+template <class, class = void>
 struct is_type_indexer : std::false_type {};
-template<class Indexer_t>
+template <class Indexer_t>
 struct is_type_indexer<Indexer_t,
-  std::void_t<decltype(std::decay_t<Indexer_t>::is_type_indexable)>>
-  : std::bool_constant<std::decay_t<Indexer_t>::is_type_indexable> {};
-template<class Indexer_t>
+                       std::void_t<decltype(std::decay_t<Indexer_t>::is_type_indexable)>>
+    : std::bool_constant<std::decay_t<Indexer_t>::is_type_indexable> {};
+template <class Indexer_t>
 constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
 
 // The "safe" version of Get(). This function will ONLY return a value IF that
 // type-based index is present in the Indexer OR if the Indexer doesn't support
 // type-based indexing.
 template <typename T, typename Indexer_t>
-inline bool safeGet(Indexer_t const& lambda, std::size_t const idx, Real& out) {
+inline bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -68,7 +68,7 @@ inline bool safeGet(Indexer_t const& lambda, std::size_t const idx, Real& out) {
 
 // Overload when no index is provided
 template <typename T, typename Indexer_t>
-inline bool safeGet(Indexer_t const& lambda, Real& out) {
+inline bool safeGet(Indexer_t const &lambda, Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -93,7 +93,7 @@ inline bool safeGet(Indexer_t const& lambda, Real& out) {
 // the "safe" version needs to separate that functionality for setting the
 // values in a lambda
 template <typename T, typename Indexer_t>
-inline bool safeSet(Indexer_t& lambda, std::size_t const idx, Real const in) {
+inline bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -122,7 +122,7 @@ inline bool safeSet(Indexer_t& lambda, std::size_t const idx, Real const in) {
 
 // Overload without numeric index
 template <typename T, typename Indexer_t>
-inline bool safeSet(Indexer_t& lambda, Real const in) {
+inline bool safeSet(Indexer_t &lambda, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index a1d7fc2d669..c7d7255313b 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -112,11 +112,11 @@ template <typename T, typename Index>
 constexpr bool is_indexable_v = is_indexable<T, Index>::value;
 
 // Check if a type can accept a size_t index
-template<class T, class=void>
+template <class T, class = void>
 struct has_whole_num_index : std::false_type {};
-template<class T>
-struct has_whole_num_index<T,
-                         std::void_t<decltype(std::declval<T>()[std::declval<std::size_t>()])>>
+template <class T>
+struct has_whole_num_index<
+    T, std::void_t<decltype(std::declval<T>()[std::declval<std::size_t>()])>>
     : std::true_type {};
 template <typename T>
 constexpr bool has_whole_num_index_v = has_whole_num_index<T>::value;

From c2b1dbb89c955bf36789862a00ddc1b9e85d8d82 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 13:30:51 -0600
Subject: [PATCH 07/30] Update changelog

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

diff --git a/CHANGELOG.md b/CHANGELOG.md
index edc7977fd1c..cc29993678e 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -7,6 +7,7 @@
 
 ### Fixed (Repair bugs, etc)
 - [[PR561]](https://github.com/lanl/singularity-eos/pull/561) Fix logic for kokkos-kernels in spackage so that it is only required for closure models on GPU
+- [[PR564]](https://github.com/lanl/singularity-eos/pull/564) Fix logic for numerical vs type indices by adding safeGet() and safeSet() helpers
 
 ### Changed (changing behavior/API/variables/...)
 

From 589e1dcd7b7b52846f2f103d6ef671efc91c09cf Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 13:31:06 -0600
Subject: [PATCH 08/30] Update doc

---
 doc/sphinx/src/using-eos.rst | 50 +++++++++++++++++++++++++++++++-----
 1 file changed, 44 insertions(+), 6 deletions(-)

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index f5e26b924d2..b14c4bed80b 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -699,6 +699,8 @@ of the ``[]`` operator that takes your type. For example:
 
   class MyLambda_t {
   public:
+    // Enable recognition that this is type-indexable
+    constexpr static bool is_type_indexable = true;
     MyLambda_t() = default;
     PORTABLE_FORCEINLINE_FUNCTION
     Real &operator[](const std::size_t idx) const {
@@ -724,17 +726,53 @@ which might be used as
 where ``MeanIonizationState`` is shorthand for index 2, since you
 defined that overload. Note that the ``operator[]`` must be marked
 ``const``. To more easily enable mixing and matching integer-based
-indexing with type-based indexing, the function
+indexing with type-based indexing, the functions
 
 .. code-block:: cpp
 
-  template<typename Name_t, typename Indexer_t>
+  template <typename T, typename Indexer_t>
   PORTABLE_FORCEINLINE_FUNCTION
-  Real &Get(Indexer_t &&lambda, std::size_t idx = 0);
+  inline bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out);
 
-will return a reference to the value at named index ``Name_t()`` if
-that overload is defined in ``Indexer_t`` and otherwise return a
-reference at index ``idx``.
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  inline bool safeGet(Indexer_t const &lambda, Real &out);
+
+will update the value of ``out`` with the value at either the appropriate
+type-based index, ``T``, or the numerical index, ``idx``, if the ``Indexer_t``
+doesn't accept type-based indexing. If the ``Indexer_t`` **does** accept
+type-based indexing but **doesn't** have the requested index, then the
+``out`` value is not updated. The same is true for when ``Indexer_t`` is the
+``nullptr``. The overload that doesn't take a numerical index will *only*
+return the value at a type-based index.
+
+Similarly, the functions
+
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  inline bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in);
+
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  inline bool safeSet(Indexer_t &lambda, Real const in)
+
+can modify the values in the ``Indexer_t`` and behave the same way. In this way,
+if a type-based index isn't present in the container, then another index won't
+be overwritten.
+
+.. note::
+
+  Previous versions defined a ``Get()`` function that was "unsafe" in the
+  sense that it would fall back on the numerical index even if a type-based
+  indexer was used. This could result in retrieving and overwriting incorrect
+  values in the indexer. We recommend not using this function and instead using
+  the "safe" versions.
 
 As a convenience tool, the struct
 

From fd74c257d475d8eb582799b0a26328623df395cc Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 14:29:57 -0600
Subject: [PATCH 09/30] Make functions PORTABLE and add required get/set

---
 singularity-eos/base/indexable_types.hpp | 60 ++++++++++++++++++++++--
 1 file changed, 56 insertions(+), 4 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index d92bb9af5bc..961689c8c73 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -39,7 +39,8 @@ constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
 // type-based index is present in the Indexer OR if the Indexer doesn't support
 // type-based indexing.
 template <typename T, typename Indexer_t>
-inline bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
+PORTABLE_FORCEINLINE_FUNCTION
+bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -68,7 +69,8 @@ inline bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
 
 // Overload when no index is provided
 template <typename T, typename Indexer_t>
-inline bool safeGet(Indexer_t const &lambda, Real &out) {
+PORTABLE_FORCEINLINE_FUNCTION
+bool safeGet(Indexer_t const &lambda, Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -89,11 +91,36 @@ inline bool safeGet(Indexer_t const &lambda, Real &out) {
   return false;
 }
 
+// Same as above but causes an error condition (static or dynamic) if the value
+// can't be obtained
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION
+Real safeMustGet(Indexer_t const &lambda, std::size_t const idx) {
+  // Error on null pointer
+  PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
+                          "Indexer can't be nullptr");
+
+  // Return type-based index. Static assert that type MUST exist in indexer
+  if constexpr (is_type_indexer_v<Indexer_t>) {
+    static_assert(variadic_utils::is_indexable_v<Indexer_t, T>);
+    return lambda[T{}];
+  }
+
+  // Fall back to numerical indexing if no type indexing
+  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
+    return lambda[idx];
+  }
+
+  // Something else...
+  PORTABLE_ALWAYS_THROW_OR_ABORT("Cannot obtain value from unknown indexer");
+}
+
 // Break out "Set" functionality from "Get". The original "Get()" did both, but
 // the "safe" version needs to separate that functionality for setting the
 // values in a lambda
 template <typename T, typename Indexer_t>
-inline bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+PORTABLE_FORCEINLINE_FUNCTION
+bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -122,7 +149,8 @@ inline bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
 
 // Overload without numeric index
 template <typename T, typename Indexer_t>
-inline bool safeSet(Indexer_t &lambda, Real const in) {
+PORTABLE_FORCEINLINE_FUNCTION
+bool safeSet(Indexer_t &lambda, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -143,6 +171,30 @@ inline bool safeSet(Indexer_t &lambda, Real const in) {
   return false;
 }
 
+// Same as above but causes an error condition (static or dynamic) if the value
+// can't be obtained
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION
+void safeMustSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+  // Error on null pointer
+  PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
+                          "Indexer can't be nullptr");
+
+  // Return type-based index. Static assert that type MUST exist in indexer
+  if constexpr (is_type_indexer_v<Indexer_t>) {
+    static_assert(variadic_utils::is_indexable_v<Indexer_t, T>);
+    lambda[T{}] = in;
+  }
+
+  // Fall back to numerical indexing if no type indexing
+  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
+    lambda[idx] = in;
+  }
+
+  // Something else...
+  PORTABLE_ALWAYS_THROW_OR_ABORT("Cannot obtain value from unknown indexer");
+}
+
 // NOTE: this Get is "unsafe" because it can allow you to overwrite a type-based
 // index since it automatically falls back to numeric indexing if the type
 // index isn't present.

From a3d71fcc7c331b22d4e5e6f1f4fb0b0dc0564e89 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 20:35:16 -0600
Subject: [PATCH 10/30] Make code more DRY and rename things a bit

---
 singularity-eos/base/indexable_types.hpp | 77 +++++++++---------------
 1 file changed, 27 insertions(+), 50 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 961689c8c73..f6457bf757a 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -20,6 +20,7 @@
 #include <utility>
 
 #include <ports-of-call/portability.hpp>
+#include <ports-of-call/portable_errors.hpp>
 #include <singularity-eos/base/variadic_utils.hpp>
 
 namespace singularity {
@@ -40,7 +41,7 @@ constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
 // type-based indexing.
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION
-bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
+bool SafeGet(Indexer_t const &lambda, Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -57,45 +58,33 @@ bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
     return false;
   }
 
-  // Fall back to numerical indexing if no type indexing
-  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
-    out = lambda[idx];
-    return true;
-  }
-
   // Something else...
   return false;
 }
 
-// Overload when no index is provided
+// Overload when numerical index is provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION
-bool safeGet(Indexer_t const &lambda, Real &out) {
-  // If null then nothing happens
-  if (variadic_utils::is_nullptr(lambda)) {
-    return false;
-  }
-
-  // Return value if type index is available
-  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
-    out = lambda[T{}];
-    return true;
-  }
+bool SafeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
+  const auto modified = SafeGet<T>(lambda, out);
 
-  // Do nothing if lambda has type indexing BUT doesn't have this type index
-  if constexpr (is_type_indexer_v<Indexer_t>) {
-    return false;
+  // Fall back to numerical indexing if no type indexing
+  if constexpr (variadic_utils::has_whole_num_index_v<Indexer_t> ||
+                variadic_utils::has_int_index_v<Indexer_t>) {
+    if (!modified) {
+      out = lambda[idx];
+      return true;
+    }
   }
 
-  // Something else...
-  return false;
+  return modified;
 }
 
 // Same as above but causes an error condition (static or dynamic) if the value
 // can't be obtained
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION
-Real safeMustGet(Indexer_t const &lambda, std::size_t const idx) {
+Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) {
   // Error on null pointer
   PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
                           "Indexer can't be nullptr");
@@ -120,7 +109,7 @@ Real safeMustGet(Indexer_t const &lambda, std::size_t const idx) {
 // values in a lambda
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION
-bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+bool SafeSet(Indexer_t &lambda, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -137,45 +126,33 @@ bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
     return false;
   }
 
-  // Fall back to numerical indexing if no type indexing
-  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
-    lambda[idx] = in;
-    return true;
-  }
-
   // Something else...
   return false;
 }
 
-// Overload without numeric index
+// Overload when numerical index is provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION
-bool safeSet(Indexer_t &lambda, Real const in) {
-  // If null then nothing happens
-  if (variadic_utils::is_nullptr(lambda)) {
-    return false;
-  }
+bool SafeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+  const auto modified = SafeSet<T>(lambda, in);
 
-  // Return value if type index is available
-  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
-    lambda[T{}] = in;
-    return true;
-  }
-
-  // Do nothing if lambda has type indexing BUT doesn't have this type index
-  if constexpr (is_type_indexer_v<Indexer_t>) {
-    return false;
+  // Fall back to numerical indexing if no type indexing
+  if constexpr (variadic_utils::has_whole_num_index_v<Indexer_t> ||
+                variadic_utils::has_int_index_v<Indexer_t>) {
+    if (!modified) {
+      lambda[idx] = in;
+      return true;
+    }
   }
 
-  // Something else...
-  return false;
+  return modified;
 }
 
 // Same as above but causes an error condition (static or dynamic) if the value
 // can't be obtained
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION
-void safeMustSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+void SafeMustSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
   // Error on null pointer
   PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
                           "Indexer can't be nullptr");

From 856531fec68d9382ca91e9e984b4611f321c053a Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 20:35:30 -0600
Subject: [PATCH 11/30] Add int index check

---
 singularity-eos/base/variadic_utils.hpp | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index c7d7255313b..22c6cf9f108 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -121,6 +121,16 @@ struct has_whole_num_index<
 template <typename T>
 constexpr bool has_whole_num_index_v = has_whole_num_index<T>::value;
 
+// Check if a type can accept an int index
+template <class T, class = void>
+struct has_int_index : std::false_type {};
+template <class T>
+struct has_int_index<
+    T, std::void_t<decltype(std::declval<T>()[std::declval<int>()])>>
+    : std::true_type {};
+template <typename T>
+constexpr bool has_int_index_v = has_int_index<T>::value;
+
 // this flattens a typelist of typelists to a single typelist
 
 // first parameter - accumulator

From 220df73d9696a9bdf4a66c21d96908d07271b18c Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 20:36:36 -0600
Subject: [PATCH 12/30] Rename safeGet/Set to SafeGet/Set and remove direct
 indexing or regular Get

---
 singularity-eos/eos/eos_electrons.hpp        |  4 +--
 singularity-eos/eos/eos_helmholtz.hpp        | 27 ++++++++++---------
 singularity-eos/eos/eos_spiner_rho_sie.hpp   | 12 ++++-----
 singularity-eos/eos/eos_stellar_collapse.hpp | 24 +++++++++--------
 singularity-eos/eos/modifiers/zsplit_eos.hpp | 11 +++-----
 test/test_indexable_types.cpp                | 28 ++++++++++----------
 6 files changed, 52 insertions(+), 54 deletions(-)

diff --git a/singularity-eos/eos/eos_electrons.hpp b/singularity-eos/eos/eos_electrons.hpp
index 0f5c262c6d0..d05b3c662f2 100644
--- a/singularity-eos/eos/eos_electrons.hpp
+++ b/singularity-eos/eos/eos_electrons.hpp
@@ -209,8 +209,8 @@ class IdealElectrons : public EosBase<IdealElectrons> {
  private:
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real _Cv(Indexer_t &&lambda) const {
-    const Real Z =
-        IndexerUtils::Get<IndexableTypes::MeanIonizationState>(lambda, Lambda::Zi);
+    Real Z = IndexerUtils::SafeMustGet<IndexableTypes::MeanIonizationState>(
+      lambda, Lambda::Zi);
     return _Cvbase * std::max(Z, static_cast<Real>(0.0));
   }
 
diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp
index 08157606b45..30d53facdcc 100644
--- a/singularity-eos/eos/eos_helmholtz.hpp
+++ b/singularity-eos/eos/eos_helmholtz.hpp
@@ -630,8 +630,8 @@ class Helmholtz : public EosBase<Helmholtz> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     using namespace HelmUtils;
     Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
-    Real abar = IndexerUtils::Get<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-    Real zbar = IndexerUtils::Get<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
     Real ytot, ye, ywot, De, lDe;
     GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
     Real lT = lTFromRhoSie_(rho, sie, abar, zbar, ye, ytot, ywot, De, lDe, lambda);
@@ -657,14 +657,14 @@ class Helmholtz : public EosBase<Helmholtz> {
       const Real rho, const Real T,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     using namespace HelmUtils;
-    return IndexerUtils::Get<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    return IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
       const Real rho, const Real T,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     using namespace HelmUtils;
-    return IndexerUtils::Get<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    return IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
   }
 
   template <typename Indexer_t = Real *>
@@ -758,10 +758,10 @@ class Helmholtz : public EosBase<Helmholtz> {
   GetFromDensityTemperature_(const Real rho, const Real temperature, Indexer_t &&lambda,
                              Real p[NDERIV], Real e[NDERIV], Real s[NDERIV],
                              Real etaele[NDERIV], Real nep[NDERIV]) const {
-    Real abar = IndexerUtils::Get<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-    Real zbar = IndexerUtils::Get<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
     Real lT = std::log10(temperature);
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+    IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
     Real ytot, ye, ywot, De, lDe;
     GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
     GetFromDensityLogTemperature_(rho, temperature, abar, zbar, ye, ytot, ywot, De, lDe,
@@ -773,8 +773,8 @@ class Helmholtz : public EosBase<Helmholtz> {
   GetFromDensityInternalEnergy_(const Real rho, const Real sie, Indexer_t &&lambda,
                                 Real p[NDERIV], Real e[NDERIV], Real s[NDERIV],
                                 Real etaele[NDERIV], Real nep[NDERIV]) const {
-    Real abar = IndexerUtils::Get<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-    Real zbar = IndexerUtils::Get<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
     Real ytot, ye, ywot, De, lDe;
     GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
     Real lT = lTFromRhoSie_(rho, sie, abar, zbar, ye, ytot, ywot, De, lDe, lambda);
@@ -819,8 +819,8 @@ Helmholtz::FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, R
   PORTABLE_ALWAYS_REQUIRE(
       !(need_temp && need_sie),
       "Either specific internal energy or temperature must be provided.");
-  Real abar = IndexerUtils::Get<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-  Real zbar = IndexerUtils::Get<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+  Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+  Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
   Real ytot, ye, ywot, De, lDe, lT;
   GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
   if (need_temp) {
@@ -870,7 +870,8 @@ PORTABLE_INLINE_FUNCTION Real Helmholtz::lTFromRhoSie_(const Real rho, const Rea
 
   if (options_.ENABLE_RAD || options_.GAS_DEGENERATE ||
       options_.ENABLE_COULOMB_CORRECTIONS) {
-    Real lTguess = IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT);
+    Real lTguess;
+    IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lTguess);
     if (!((electrons_.lTMin() <= lTguess) && (lTguess <= electrons_.lTMax()))) {
       lTguess = lTAnalytic_(rho, e, ni, ne);
       if (!((electrons_.lTMin() <= lTguess) && (lTguess <= electrons_.lTMax()))) {
@@ -944,7 +945,7 @@ PORTABLE_INLINE_FUNCTION Real Helmholtz::lTFromRhoSie_(const Real rho, const Rea
     }
     lT = electrons_.lTMax();
   }
-  IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+  IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   return lT;
 }
 
diff --git a/singularity-eos/eos/eos_spiner_rho_sie.hpp b/singularity-eos/eos/eos_spiner_rho_sie.hpp
index 44db0c76f0d..8c0a035b042 100644
--- a/singularity-eos/eos/eos_spiner_rho_sie.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_sie.hpp
@@ -687,7 +687,7 @@ PORTABLE_INLINE_FUNCTION void SpinerEOSDependsRhoSieTransformable<TransformerT>:
     }
   } else {
     lRho = to_log(rho, lRhoOffset_);
-    IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+    IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
   }
   if (output & thermalqs::temperature) {
 
@@ -745,7 +745,7 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::interpRhoT_(const Real rho,
                                                                Indexer_t &&lambda) const {
   const Real lRho = spiner_common::to_log(rho, lRhoOffset_);
   const Real lT = spiner_common::to_log(T, lTOffset_);
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
   return db.interpToReal(lRho, lT);
 }
 
@@ -757,7 +757,7 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::interpRhoSie_(
   const Real lRho = spiner_common::to_log(rho, lRhoOffset_);
   const Real sie_transformed = transformer_.transform(sie, rho);
   const Real lE = spiner_common::to_log(sie_transformed, lEOffset_);
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
   return db.interpToReal(lRho, lE);
 }
 
@@ -781,7 +781,7 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::lRhoFromPlT_(
   } else {
     Real lRhoGuess = reproducible_ ? lRhoMin_ : 0.5 * (lRhoMin_ + lRhoMax_);
     Real lRho_cache;
-    IndexerUtils::safeGet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho_cache);
+    IndexerUtils::SafeGet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho_cache);
     if ((lRhoMin_ <= lRho_cache) && (lRho_cache <= lRhoMax_)) {
       lRhoGuess = lRho_cache;
     }
@@ -801,8 +801,8 @@ SpinerEOSDependsRhoSieTransformable<TransformerT>::lRhoFromPlT_(
       lRho = reproducible_ ? lRhoMin_ : lRhoGuess;
     }
   }
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
-  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
   return lRho;
 }
 
diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp
index 6f98e252741..27608787cb6 100644
--- a/singularity-eos/eos/eos_stellar_collapse.hpp
+++ b/singularity-eos/eos/eos_stellar_collapse.hpp
@@ -372,8 +372,8 @@ class StellarCollapse : public EosBase<StellarCollapse> {
     checkLambda_(lambda);
     lRho = lRho_(rho);
     lT = lT_(temp);
-    Ye = IndexerUtils::Get<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+    Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+    IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   }
   template <typename Indexer_t>
   PORTABLE_INLINE_FUNCTION __attribute__((always_inline)) void
@@ -381,7 +381,7 @@ class StellarCollapse : public EosBase<StellarCollapse> {
                      Real &lT, Real &Ye) const noexcept {
     lRho = lRho_(rho);
     lT = lTFromlRhoSie_(lRho, sie, lambda);
-    Ye = IndexerUtils::Get<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+    Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
     return;
   }
 
@@ -577,7 +577,8 @@ template <typename Indexer_t>
 PORTABLE_INLINE_FUNCTION Real
 StellarCollapse::MinInternalEnergyFromDensity(const Real rho, Indexer_t &&lambda) const {
   Real lRho = lRho_(rho);
-  Real Ye = IndexerUtils::Get<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+  Real Ye;
+  Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
   return eCold_.interpToReal(Ye, lRho);
 }
 
@@ -667,7 +668,7 @@ PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperat
   Real lrguess = lRho_(rho);
   Real lT = lT_(temp);
   Real lP = P2lP_(press);
-  Real Ye = IndexerUtils::Get<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+  Real Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
 
   if ((lrguess < lRhoMin_) || (lrguess > lRhoMax_)) {
     lrguess = lRho_(rhoNormal_);
@@ -680,7 +681,7 @@ PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperat
   Real lE = lE_.interpToReal(Ye, lT, lrguess);
   rho = rho_(lrguess);
   sie = le2e_(lE);
-  IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
 template <typename Indexer_t>
@@ -762,8 +763,8 @@ StellarCollapse::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &
   dpde = dPdENormal_;
   dvdt = dVdTNormal_;
   Real lT = lT_(temp);
-  IndexerUtils::Get<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye) = YeNormal_;
-  IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+  IndexerUtils::SafeMustSet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye, YeNormal_);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
 inline void StellarCollapse::LoadFromSP5File_(const std::string &filename) {
@@ -1177,8 +1178,9 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::lTFromlRhoSie_(
   RootFinding1D::Status status = RootFinding1D::Status::SUCCESS;
   using RootFinding1D::regula_falsi;
   Real lT;
-  Real Ye = IndexerUtils::Get<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
-  Real lTGuess = IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT);
+  Real Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+  Real lTGuess;
+  IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lTGuess);
 
   const RootFinding1D::RootCounts *pcounts =
       (memoryStatus_ == DataStatus::OnDevice) ? nullptr : &counts;
@@ -1224,7 +1226,7 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::lTFromlRhoSie_(
     status_ = status;
   }
 #endif // PORTABILITY_STRATEGY_NONE
-  IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   return lT;
 }
 } // namespace singularity
diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
index 6192f17c63a..c1715cba71a 100644
--- a/singularity-eos/eos/modifiers/zsplit_eos.hpp
+++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -275,14 +275,9 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
   template <typename Indexer_t = Real *>
   PORTABLE_FORCEINLINE_FUNCTION Real GetIonizationState_(Indexer_t &&lambda) const {
     using namespace variadic_utils;
-    if (is_nullptr(lambda)) {
-      PORTABLE_THROW_OR_ABORT("ZSplitEOS: lambda must contain mean ionization state!\n");
-    }
-    if constexpr (is_indexable_v<Indexer_t, IndexableTypes::MeanIonizationState>) {
-      return std::max(0.0, lambda[IndexableTypes::MeanIonizationState()]);
-    } else {
-      return std::max(0.0, lambda[T::nlambda()]);
-    }
+    return std::max(0.0, IndexerUtils::SafeMustGet<IndexableTypes::MeanIonizationState>(
+      lambda, T::nlambda()
+    ));
   }
   // TODO(JMM): Runtime?
   template <typename Indexer_t = Real *>
diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index 40edc76c634..8f1be8d913f 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -81,18 +81,18 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
         REQUIRE(lRho == static_cast<Real>(2));
       }
     }
-    WHEN("We use the safeGet functionality") {
+    WHEN("We use the SafeGet functionality") {
       constexpr Real unmodified = -1.0;
       Real destination = unmodified;
       WHEN("We request a type that exists") {
-        const bool modified = safeGet<MeanIonizationState>(lambda, 2, destination);
+        const bool modified = SafeGet<MeanIonizationState>(lambda, 2, destination);
         THEN("The destination value will be modified") {
           CHECK(modified);
           REQUIRE(destination == lambda[MeanIonizationState{}]);
         }
       }
       WHEN("We request a type that doesn't exist") {
-        const bool modified = safeGet<TableStatus>(lambda, 2, destination);
+        const bool modified = SafeGet<TableStatus>(lambda, 2, destination);
         THEN("The destination value will remain UNmodified") {
           CHECK(!modified);
           REQUIRE(destination == unmodified);
@@ -101,24 +101,24 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
       WHEN("A normal array-like lambda is used") {
         std::array<Real, Lambda_t::size()> lambda_arr{1, 2, 3};
         constexpr size_t my_index = 2;
-        const bool modified = safeGet<LogDensity>(lambda_arr, my_index, destination);
+        const bool modified = SafeGet<LogDensity>(lambda_arr, my_index, destination);
         THEN("The destination value will reflect the index from the array") {
           CHECK(modified);
           REQUIRE(destination == lambda_arr[my_index]);
         }
       }
     }
-    WHEN("We use the safeSet functionality") {
+    WHEN("We use the SafeSet functionality") {
       constexpr Real new_value = -1.0;
       WHEN("We want to set a value for a type index that exists") {
-        const bool modified = safeSet<MeanIonizationState>(lambda, 2, new_value);
+        const bool modified = SafeSet<MeanIonizationState>(lambda, 2, new_value);
         THEN("The lambda index was modified") {
           CHECK(modified);
           REQUIRE(lambda[MeanIonizationState{}] == new_value);
         }
       }
       WHEN("We want to set a value for a type index that doesn't exist") {
-        const bool modified = safeSet<TableStatus>(lambda, 2, new_value);
+        const bool modified = SafeSet<TableStatus>(lambda, 2, new_value);
         Lambda_t old_lambda;
         for (std::size_t i = 0; i < Lambda_t::size(); ++i) {
           old_lambda[i] = lambda[i];
@@ -134,7 +134,7 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
       WHEN("A normal array-like lambda is used") {
         std::array<Real, Lambda_t::size()> lambda_arr{4, 5, 6};
         constexpr size_t my_index = 1;
-        const bool modified = safeSet<LogDensity>(lambda_arr, my_index, new_value);
+        const bool modified = SafeSet<LogDensity>(lambda_arr, my_index, new_value);
         THEN("The lambda value at the appropriate index has been modified") {
           CHECK(modified);
           REQUIRE(lambda_arr[my_index] == new_value);
@@ -160,12 +160,12 @@ SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
         REQUIRE(lRho == static_cast<Real>(2));
       }
     }
-    WHEN("We use the safeGet functionality") {
+    WHEN("We use the SafeGet functionality") {
       constexpr Real unmodified = -1.0;
       Real destination = unmodified;
       WHEN("We request a type that doesn't exist in the manual indexer") {
         constexpr size_t my_index = 1;
-        const bool modified = safeGet<TableStatus>(lambda, my_index, destination);
+        const bool modified = SafeGet<TableStatus>(lambda, my_index, destination);
         THEN("The destination WILL be modified since the manual indexer doesn't have the "
              " `is_type_indexable` data member") {
           CHECK(modified);
@@ -181,7 +181,7 @@ SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
         WHEN("We request a type that doesn't exist in the manual indexer") {
           destination = unmodified;
           constexpr size_t my_index = 1;
-          const bool modified = safeGet<TableStatus>(lambda_new, my_index, destination);
+          const bool modified = SafeGet<TableStatus>(lambda_new, my_index, destination);
           THEN("The destination will NOT be modified") {
             CHECK(!modified);
             REQUIRE(destination == unmodified);
@@ -189,11 +189,11 @@ SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
         }
       }
     }
-    WHEN("We use the safeSet functionality") {
+    WHEN("We use the SafeSet functionality") {
       constexpr Real new_value = -1.0;
       WHEN("We want to set a value for a type index that doesn't exist") {
         constexpr size_t my_index = 1;
-        const bool modified = safeSet<TableStatus>(lambda, my_index, new_value);
+        const bool modified = SafeSet<TableStatus>(lambda, my_index, new_value);
         THEN("The lambda value WILL be modified since the manual indexer doesn't have "
              "the `is_type_indexable` data member") {
           CHECK(modified);
@@ -208,7 +208,7 @@ SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
         }
         WHEN("We request a type that doesn't exist in the manual indexer") {
           constexpr size_t my_index = 1;
-          const bool modified = safeSet<TableStatus>(lambda_new, my_index, new_value);
+          const bool modified = SafeSet<TableStatus>(lambda_new, my_index, new_value);
           THEN("The lambda will NOT be modified") {
             CHECK(!modified);
             for (std::size_t i = 0; i < lambda.length; ++i) {

From 271e52e7327fe19ed3441327a9c94b6b298f0b97 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 20:37:29 -0600
Subject: [PATCH 13/30] Make indexer const correct

---
 test/test_pte_ideal.cpp | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/test_pte_ideal.cpp b/test/test_pte_ideal.cpp
index 31c91bfaa6d..f318472b251 100644
--- a/test/test_pte_ideal.cpp
+++ b/test/test_pte_ideal.cpp
@@ -51,8 +51,12 @@ using singularity::IndexableTypes::MeanIonizationState;
 struct LambdaIndexerSingle {
   PORTABLE_FORCEINLINE_FUNCTION
   Real &operator[](const int i) { return z; }
+  PORTABLE_INLINE_FUNCTION
+  const Real &operator[](const int i) const { return z; }
   PORTABLE_FORCEINLINE_FUNCTION
   Real &operator[](const MeanIonizationState &s) { return z; }
+  PORTABLE_INLINE_FUNCTION
+  const Real &operator[](const MeanIonizationState &s) const { return z; }
   Real z = 0.9;
 };
 

From 63d5a3978e560bc57f07c9ae6c8b28dd49ebd283 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 20:37:50 -0600
Subject: [PATCH 14/30] Rename safeGet/Set

---
 singularity-eos/eos/eos_spiner_rho_temp.hpp | 38 ++++++++++-----------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/singularity-eos/eos/eos_spiner_rho_temp.hpp b/singularity-eos/eos/eos_spiner_rho_temp.hpp
index 947e70e574b..3015fbd7a55 100644
--- a/singularity-eos/eos/eos_spiner_rho_temp.hpp
+++ b/singularity-eos/eos/eos_spiner_rho_temp.hpp
@@ -880,8 +880,8 @@ SpinerEOSDependsRhoT::FillEos(Real &rho, Real &temp, Real &energy, Real &press,
   if (output & thermalqs::bulk_modulus) {
     bmod = bModFromRholRhoTlT_(rho, lRho, temp, lT, whereAmI);
   }
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
-  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
 template <typename Indexer_t>
@@ -912,8 +912,8 @@ SpinerEOSDependsRhoT::getLogsRhoT_(const Real rho, const Real temperature, Real
                                    Real &lT, Indexer_t &&lambda) const {
   lRho = lRho_(rho);
   lT = lT_(temperature);
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
-  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
 template <typename Indexer_t>
@@ -928,7 +928,7 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lRhoFromPlT_(
       (memoryStatus_ == DataStatus::OnDevice) ? nullptr : &counts;
 
   Real lRho_cache;
-  IndexerUtils::safeGet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho_cache);
+  IndexerUtils::SafeGet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho_cache);
   if ((lRhoMin_ <= lRho_cache) && (lRho_cache <= lRhoMax_)) {
     lRhoGuess = lRho_cache;
   }
@@ -967,11 +967,11 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lRhoFromPlT_(
 #endif // SPINER_EOS_VERBOSE
     lRho = reproducible_ ? lRhoMax_ : lRhoGuess;
   }
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
-  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   // No numerical index: only set if type-based indexing is used
-  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
-  IndexerUtils::safeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
+  IndexerUtils::SafeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::SafeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
   return lRho;
 }
 
@@ -1003,7 +1003,7 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoSie_(
   } else {
     Real lTGuess = reproducible_ ? lTMin_ : 0.5 * (lTMin_ + lTMax_);
     Real lT_cache;
-    IndexerUtils::safeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT_cache);
+    IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT_cache);
     if ((lTMin_ <= lT_cache) && (lT_cache <= lTMax_)) {
       lTGuess = lT_cache;
     }
@@ -1027,11 +1027,11 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoSie_(
       lT = reproducible_ ? lTMin_ : lTGuess;
     }
   }
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
-  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   // No numerical index: only set if type-based indexing is used
-  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
-  IndexerUtils::safeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
+  IndexerUtils::SafeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::SafeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
   return lT;
 }
 
@@ -1062,7 +1062,7 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoP_(
     whereAmI = TableStatus::OnTable;
     lTGuess = 0.5 * (lTMin_ + lTMax_);
     Real lT_cache;
-    IndexerUtils::safeGet<IndexableTypes::LogTemperature>(lambda, lT, lT_cache);
+    IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, lT, lT_cache);
     if ((lTMin_ <= lT_cache) && (lT_cache <= lTMax_)) {
       lTGuess = lT_cache;
     }
@@ -1082,11 +1082,11 @@ PORTABLE_INLINE_FUNCTION Real SpinerEOSDependsRhoT::lTFromlRhoP_(
       lT = reproducible_ ? lTMin_ : lTGuess;
     }
   }
-  IndexerUtils::safeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
-  IndexerUtils::safeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
+  IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, Lambda::lRho, lRho);
+  IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   // No numerical index: only set if type-based indexing is used
-  IndexerUtils::safeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
-  IndexerUtils::safeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
+  IndexerUtils::SafeSet<IndexableTypes::RootStatus>(lambda, static_cast<Real>(status));
+  IndexerUtils::SafeSet<IndexableTypes::TableStatus>(lambda, static_cast<Real>(whereAmI));
   return lT;
 }
 

From 8a893a79b2d8d654768da3579d8c68ef3a925729 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 20:38:31 -0600
Subject: [PATCH 15/30] Clang format

---
 singularity-eos/base/indexable_types.hpp     | 22 +++++++-------
 singularity-eos/base/variadic_utils.hpp      |  3 +-
 singularity-eos/eos/eos_electrons.hpp        |  4 +--
 singularity-eos/eos/eos_helmholtz.hpp        | 30 +++++++++++++-------
 singularity-eos/eos/eos_stellar_collapse.hpp |  9 ++++--
 singularity-eos/eos/modifiers/zsplit_eos.hpp |  3 +-
 6 files changed, 40 insertions(+), 31 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index f6457bf757a..15e9c15c914 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -40,8 +40,7 @@ constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
 // type-based index is present in the Indexer OR if the Indexer doesn't support
 // type-based indexing.
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION
-bool SafeGet(Indexer_t const &lambda, Real &out) {
+PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -64,8 +63,8 @@ bool SafeGet(Indexer_t const &lambda, Real &out) {
 
 // Overload when numerical index is provided
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION
-bool SafeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
+PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx,
+                                           Real &out) {
   const auto modified = SafeGet<T>(lambda, out);
 
   // Fall back to numerical indexing if no type indexing
@@ -83,8 +82,8 @@ bool SafeGet(Indexer_t const &lambda, std::size_t const idx, Real &out) {
 // Same as above but causes an error condition (static or dynamic) if the value
 // can't be obtained
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION
-Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) {
+PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
+                                               std::size_t const idx) {
   // Error on null pointer
   PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
                           "Indexer can't be nullptr");
@@ -108,8 +107,7 @@ Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) {
 // the "safe" version needs to separate that functionality for setting the
 // values in a lambda
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION
-bool SafeSet(Indexer_t &lambda, Real const in) {
+PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -132,8 +130,8 @@ bool SafeSet(Indexer_t &lambda, Real const in) {
 
 // Overload when numerical index is provided
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION
-bool SafeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx,
+                                           Real const in) {
   const auto modified = SafeSet<T>(lambda, in);
 
   // Fall back to numerical indexing if no type indexing
@@ -151,8 +149,8 @@ bool SafeSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
 // Same as above but causes an error condition (static or dynamic) if the value
 // can't be obtained
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION
-void SafeMustSet(Indexer_t &lambda, std::size_t const idx, Real const in) {
+PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx,
+                                               Real const in) {
   // Error on null pointer
   PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
                           "Indexer can't be nullptr");
diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index 22c6cf9f108..ee32a04908c 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -125,8 +125,7 @@ constexpr bool has_whole_num_index_v = has_whole_num_index<T>::value;
 template <class T, class = void>
 struct has_int_index : std::false_type {};
 template <class T>
-struct has_int_index<
-    T, std::void_t<decltype(std::declval<T>()[std::declval<int>()])>>
+struct has_int_index<T, std::void_t<decltype(std::declval<T>()[std::declval<int>()])>>
     : std::true_type {};
 template <typename T>
 constexpr bool has_int_index_v = has_int_index<T>::value;
diff --git a/singularity-eos/eos/eos_electrons.hpp b/singularity-eos/eos/eos_electrons.hpp
index d05b3c662f2..7ac3f61842d 100644
--- a/singularity-eos/eos/eos_electrons.hpp
+++ b/singularity-eos/eos/eos_electrons.hpp
@@ -209,8 +209,8 @@ class IdealElectrons : public EosBase<IdealElectrons> {
  private:
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real _Cv(Indexer_t &&lambda) const {
-    Real Z = IndexerUtils::SafeMustGet<IndexableTypes::MeanIonizationState>(
-      lambda, Lambda::Zi);
+    Real Z = IndexerUtils::SafeMustGet<IndexableTypes::MeanIonizationState>(lambda,
+                                                                            Lambda::Zi);
     return _Cvbase * std::max(Z, static_cast<Real>(0.0));
   }
 
diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp
index 30d53facdcc..566de55aa7c 100644
--- a/singularity-eos/eos/eos_helmholtz.hpp
+++ b/singularity-eos/eos/eos_helmholtz.hpp
@@ -630,8 +630,10 @@ class Helmholtz : public EosBase<Helmholtz> {
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     using namespace HelmUtils;
     Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
-    Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-    Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    Real abar =
+        IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    Real zbar =
+        IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
     Real ytot, ye, ywot, De, lDe;
     GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
     Real lT = lTFromRhoSie_(rho, sie, abar, zbar, ye, ytot, ywot, De, lDe, lambda);
@@ -657,14 +659,16 @@ class Helmholtz : public EosBase<Helmholtz> {
       const Real rho, const Real T,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     using namespace HelmUtils;
-    return IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    return IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda,
+                                                                     Lambda::Abar);
   }
   template <typename Indexer_t = Real *>
   PORTABLE_INLINE_FUNCTION Real MeanAtomicNumberFromDensityTemperature(
       const Real rho, const Real T,
       Indexer_t &&lambda = static_cast<Real *>(nullptr)) const {
     using namespace HelmUtils;
-    return IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    return IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda,
+                                                                       Lambda::Zbar);
   }
 
   template <typename Indexer_t = Real *>
@@ -758,8 +762,10 @@ class Helmholtz : public EosBase<Helmholtz> {
   GetFromDensityTemperature_(const Real rho, const Real temperature, Indexer_t &&lambda,
                              Real p[NDERIV], Real e[NDERIV], Real s[NDERIV],
                              Real etaele[NDERIV], Real nep[NDERIV]) const {
-    Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-    Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    Real abar =
+        IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    Real zbar =
+        IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
     Real lT = std::log10(temperature);
     IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
     Real ytot, ye, ywot, De, lDe;
@@ -773,8 +779,10 @@ class Helmholtz : public EosBase<Helmholtz> {
   GetFromDensityInternalEnergy_(const Real rho, const Real sie, Indexer_t &&lambda,
                                 Real p[NDERIV], Real e[NDERIV], Real s[NDERIV],
                                 Real etaele[NDERIV], Real nep[NDERIV]) const {
-    Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-    Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+    Real abar =
+        IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+    Real zbar =
+        IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
     Real ytot, ye, ywot, De, lDe;
     GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
     Real lT = lTFromRhoSie_(rho, sie, abar, zbar, ye, ytot, ywot, De, lDe, lambda);
@@ -819,8 +827,10 @@ Helmholtz::FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, R
   PORTABLE_ALWAYS_REQUIRE(
       !(need_temp && need_sie),
       "Either specific internal energy or temperature must be provided.");
-  Real abar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
-  Real zbar = IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
+  Real abar =
+      IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicMass>(lambda, Lambda::Abar);
+  Real zbar =
+      IndexerUtils::SafeMustGet<IndexableTypes::MeanAtomicNumber>(lambda, Lambda::Zbar);
   Real ytot, ye, ywot, De, lDe, lT;
   GetElectronDensities_(rho, abar, zbar, ytot, ye, ywot, De, lDe);
   if (need_temp) {
diff --git a/singularity-eos/eos/eos_stellar_collapse.hpp b/singularity-eos/eos/eos_stellar_collapse.hpp
index 27608787cb6..eaa79e03b88 100644
--- a/singularity-eos/eos/eos_stellar_collapse.hpp
+++ b/singularity-eos/eos/eos_stellar_collapse.hpp
@@ -668,7 +668,8 @@ PORTABLE_INLINE_FUNCTION void StellarCollapse::DensityEnergyFromPressureTemperat
   Real lrguess = lRho_(rho);
   Real lT = lT_(temp);
   Real lP = P2lP_(press);
-  Real Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+  Real Ye =
+      IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
 
   if ((lrguess < lRhoMin_) || (lrguess > lRhoMax_)) {
     lrguess = lRho_(rhoNormal_);
@@ -763,7 +764,8 @@ StellarCollapse::ValuesAtReferenceState(Real &rho, Real &temp, Real &sie, Real &
   dpde = dPdENormal_;
   dvdt = dVdTNormal_;
   Real lT = lT_(temp);
-  IndexerUtils::SafeMustSet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye, YeNormal_);
+  IndexerUtils::SafeMustSet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye,
+                                                              YeNormal_);
   IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
 }
 
@@ -1178,7 +1180,8 @@ PORTABLE_INLINE_FUNCTION Real StellarCollapse::lTFromlRhoSie_(
   RootFinding1D::Status status = RootFinding1D::Status::SUCCESS;
   using RootFinding1D::regula_falsi;
   Real lT;
-  Real Ye = IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
+  Real Ye =
+      IndexerUtils::SafeMustGet<IndexableTypes::ElectronFraction>(lambda, Lambda::Ye);
   Real lTGuess;
   IndexerUtils::SafeGet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lTGuess);
 
diff --git a/singularity-eos/eos/modifiers/zsplit_eos.hpp b/singularity-eos/eos/modifiers/zsplit_eos.hpp
index c1715cba71a..ce850acd4e4 100644
--- a/singularity-eos/eos/modifiers/zsplit_eos.hpp
+++ b/singularity-eos/eos/modifiers/zsplit_eos.hpp
@@ -276,8 +276,7 @@ class ZSplit : public EosBase<ZSplit<ztype, T>> {
   PORTABLE_FORCEINLINE_FUNCTION Real GetIonizationState_(Indexer_t &&lambda) const {
     using namespace variadic_utils;
     return std::max(0.0, IndexerUtils::SafeMustGet<IndexableTypes::MeanIonizationState>(
-      lambda, T::nlambda()
-    ));
+                             lambda, T::nlambda()));
   }
   // TODO(JMM): Runtime?
   template <typename Indexer_t = Real *>

From d6a9f689bfb48c394ca1d457b8ea1a64d81afba8 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 21:31:44 -0600
Subject: [PATCH 16/30] Switch to template-based decision to use integer index
 or not

---
 singularity-eos/base/indexable_types.hpp | 135 +++++++++++++++--------
 singularity-eos/base/variadic_utils.hpp  |  10 --
 2 files changed, 90 insertions(+), 55 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 15e9c15c914..1c5ec9d2e99 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -36,11 +36,17 @@ struct is_type_indexer<Indexer_t,
 template <class Indexer_t>
 constexpr bool is_type_indexer_v = is_type_indexer<Indexer_t>::value;
 
+namespace impl {
+
+// Simple way to switch between pure type indexing or also allowing intergers
+enum class AllowedIndexing { Numeric, TypeOnly };
+
 // The "safe" version of Get(). This function will ONLY return a value IF that
 // type-based index is present in the Indexer OR if the Indexer doesn't support
 // type-based indexing.
-template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out) {
+template <AllowedIndexing AI, typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx,
+                                           Real &out) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -57,31 +63,21 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out) {
     return false;
   }
 
-  // Something else...
-  return false;
-}
-
-// Overload when numerical index is provided
-template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx,
-                                           Real &out) {
-  const auto modified = SafeGet<T>(lambda, out);
-
-  // Fall back to numerical indexing if no type indexing
-  if constexpr (variadic_utils::has_whole_num_index_v<Indexer_t> ||
-                variadic_utils::has_int_index_v<Indexer_t>) {
-    if (!modified) {
+  // Fall back to numeric indexing if allowed
+  if constexpr (AI == AllowedIndexing::Numeric) {
+    if constexpr (variadic_utils::has_int_index_v<Indexer_t>) {
       out = lambda[idx];
       return true;
     }
   }
 
-  return modified;
+  // Something else...
+  return false;
 }
 
 // Same as above but causes an error condition (static or dynamic) if the value
 // can't be obtained
-template <typename T, typename Indexer_t>
+template <AllowedIndexing AI, typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
                                                std::size_t const idx) {
   // Error on null pointer
@@ -94,9 +90,11 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
     return lambda[T{}];
   }
 
-  // Fall back to numerical indexing if no type indexing
-  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
-    return lambda[idx];
+  // Fall back to numerical indexing if allowed
+  if constexpr (AI == AllowedIndexing::Numeric) {
+    if constexpr (variadic_utils::has_int_index<Indexer_t>::value) {
+      return lambda[idx];
+    }
   }
 
   // Something else...
@@ -106,8 +104,9 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
 // Break out "Set" functionality from "Get". The original "Get()" did both, but
 // the "safe" version needs to separate that functionality for setting the
 // values in a lambda
-template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) {
+template <AllowedIndexing AI, typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx,
+                                           Real const in) {
   // If null then nothing happens
   if (variadic_utils::is_nullptr(lambda)) {
     return false;
@@ -124,31 +123,21 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) {
     return false;
   }
 
-  // Something else...
-  return false;
-}
-
-// Overload when numerical index is provided
-template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx,
-                                           Real const in) {
-  const auto modified = SafeSet<T>(lambda, in);
-
-  // Fall back to numerical indexing if no type indexing
-  if constexpr (variadic_utils::has_whole_num_index_v<Indexer_t> ||
-                variadic_utils::has_int_index_v<Indexer_t>) {
-    if (!modified) {
+  // Fall back to numeric indexing if allowed
+  if constexpr (AI == AllowedIndexing::Numeric) {
+    if constexpr (variadic_utils::has_int_index_v<Indexer_t>) {
       lambda[idx] = in;
       return true;
     }
   }
 
-  return modified;
+  // Something else...
+  return false;
 }
 
 // Same as above but causes an error condition (static or dynamic) if the value
 // can't be obtained
-template <typename T, typename Indexer_t>
+template <AllowedIndexing AI, typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx,
                                                Real const in) {
   // Error on null pointer
@@ -161,21 +150,77 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t co
     lambda[T{}] = in;
   }
 
-  // Fall back to numerical indexing if no type indexing
-  if constexpr (variadic_utils::has_whole_num_index<Indexer_t>::value) {
-    lambda[idx] = in;
+  // Fall back to numerical indexing if allowed
+  if constexpr (AI == AllowedIndexing::Numeric) {
+    if constexpr (variadic_utils::has_int_index<Indexer_t>::value) {
+      lambda[idx] = in;
+    }
   }
 
   // Something else...
   PORTABLE_ALWAYS_THROW_OR_ABORT("Cannot obtain value from unknown indexer");
 }
 
+} // namespace impl
+
+// Overload when numerical index is provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t const idx,
+                                           Real &out) {
+  return impl::SafeGet<impl::AllowedIndexing::Numeric, T>(lambda, idx, out);
+}
+
+// Overload when numerical index isn't provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, Real &out) {
+  std::size_t idx = 0;
+  return impl::SafeGet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx, out);
+}
+
+// Overload when numerical index is provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const idx,
+                                           Real const in) {
+  return impl::SafeSet<impl::AllowedIndexing::Numeric, T>(lambda, idx, in);
+}
+
+// Overload when numerical index isn't provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) {
+  std::size_t idx = 0;
+  return impl::SafeSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx, in);
+}
+
+// Overload when numerical index is provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) {
+  return impl::SafeMustGet<impl::AllowedIndexing::Numeric, T>(lambda, idx);
+}
+
+// Overload when numerical index isn't provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) {
+  std::size_t idx = 0;
+  return impl::SafeMustGet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx);
+}
+
+// Overload when numerical index is provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx,
+                                           Real const in) {
+  return impl::SafeMustSet<impl::AllowedIndexing::Numeric, T>(lambda, idx, in);
+}
+
+// Overload when numerical index isn't provided
+template <typename T, typename Indexer_t>
+PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, Real const in) {
+  std::size_t idx = 0;
+  return impl::SafeMustSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx, in);
+}
+
 // NOTE: this Get is "unsafe" because it can allow you to overwrite a type-based
 // index since it automatically falls back to numeric indexing if the type
 // index isn't present.
-
-// Convenience function for accessing an indexer by either type or natural
-// number index depending on what is available.
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION auto &Get(Indexer_t &&lambda, std::size_t idx = 0) {
   if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index ee32a04908c..d2aa3aba3cc 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -111,16 +111,6 @@ struct is_indexable<T, Index,
 template <typename T, typename Index>
 constexpr bool is_indexable_v = is_indexable<T, Index>::value;
 
-// Check if a type can accept a size_t index
-template <class T, class = void>
-struct has_whole_num_index : std::false_type {};
-template <class T>
-struct has_whole_num_index<
-    T, std::void_t<decltype(std::declval<T>()[std::declval<std::size_t>()])>>
-    : std::true_type {};
-template <typename T>
-constexpr bool has_whole_num_index_v = has_whole_num_index<T>::value;
-
 // Check if a type can accept an int index
 template <class T, class = void>
 struct has_int_index : std::false_type {};

From cf8bf93f2cc1c8e7e422e1f87216bb195a63b2cb Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 21:33:35 -0600
Subject: [PATCH 17/30] Whoops... forgot to return

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

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 1c5ec9d2e99..8725f46f073 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -148,12 +148,14 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t co
   if constexpr (is_type_indexer_v<Indexer_t>) {
     static_assert(variadic_utils::is_indexable_v<Indexer_t, T>);
     lambda[T{}] = in;
+    return;
   }
 
   // Fall back to numerical indexing if allowed
   if constexpr (AI == AllowedIndexing::Numeric) {
     if constexpr (variadic_utils::has_int_index<Indexer_t>::value) {
       lambda[idx] = in;
+      return;
     }
   }
 

From 2d4153a870d8e29e172c68771dc2507aedf8cb4c Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 21:34:49 -0600
Subject: [PATCH 18/30] Clang format

---
 singularity-eos/base/indexable_types.hpp | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 8725f46f073..dd96ced7999 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -195,7 +195,8 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) {
 
 // Overload when numerical index is provided
 template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx) {
+PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
+                                               std::size_t const idx) {
   return impl::SafeMustGet<impl::AllowedIndexing::Numeric, T>(lambda, idx);
 }
 
@@ -209,7 +210,7 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) {
 // Overload when numerical index is provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx,
-                                           Real const in) {
+                                               Real const in) {
   return impl::SafeMustSet<impl::AllowedIndexing::Numeric, T>(lambda, idx, in);
 }
 

From 0a4baa13c1e24e0116097cd6a2631507db1c3e63 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Fri, 12 Sep 2025 21:44:20 -0600
Subject: [PATCH 19/30] Add docs for SafeMustGet() and SafeMustSet()

---
 doc/sphinx/src/using-eos.rst | 42 ++++++++++++++++++++++++++++++------
 1 file changed, 36 insertions(+), 6 deletions(-)

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index b14c4bed80b..bcdce9792fd 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -732,13 +732,25 @@ indexing with type-based indexing, the functions
 
   template <typename T, typename Indexer_t>
   PORTABLE_FORCEINLINE_FUNCTION
-  inline bool safeGet(Indexer_t const &lambda, std::size_t const idx, Real &out);
+  bool SafeGet(Indexer_t const &lambda, std::size_t const idx, Real &out);
 
 .. code-block:: cpp
 
   template <typename T, typename Indexer_t>
   PORTABLE_FORCEINLINE_FUNCTION
-  inline bool safeGet(Indexer_t const &lambda, Real &out);
+  bool SafeGet(Indexer_t const &lambda, Real &out);
+
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real SafeMustGet(Indexer_t const &lambda, std::size_t const idx)
+
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  Real SafeMustGet(Indexer_t const &lambda)
 
 will update the value of ``out`` with the value at either the appropriate
 type-based index, ``T``, or the numerical index, ``idx``, if the ``Indexer_t``
@@ -746,7 +758,11 @@ doesn't accept type-based indexing. If the ``Indexer_t`` **does** accept
 type-based indexing but **doesn't** have the requested index, then the
 ``out`` value is not updated. The same is true for when ``Indexer_t`` is the
 ``nullptr``. The overload that doesn't take a numerical index will *only*
-return the value at a type-based index.
+return the value at a type-based index. The ``SafeMustGet()`` version
+is intended to generate errors if a value can't be retrieved. In the case of
+type-based indexing, the error will be at compile time if the type isn't located
+in the indexer. A runtime abort will occur if either the null pointer is passed
+or if integer indexing isn't allowed for some reason.
 
 Similarly, the functions
 
@@ -754,17 +770,31 @@ Similarly, the functions
 
   template <typename T, typename Indexer_t>
   PORTABLE_FORCEINLINE_FUNCTION
-  inline bool safeSet(Indexer_t &lambda, std::size_t const idx, Real const in);
+  inline bool SafeSet(Indexer_t &lambda, std::size_t const idx, Real const in);
+
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  inline bool SafeSet(Indexer_t &lambda, Real const in)
+
+.. code-block:: cpp
+
+  template <typename T, typename Indexer_t>
+  PORTABLE_FORCEINLINE_FUNCTION
+  inline bool SafeMustSet(Indexer_t &lambda, std::size_t const idx, Real const in);
 
 .. code-block:: cpp
 
   template <typename T, typename Indexer_t>
   PORTABLE_FORCEINLINE_FUNCTION
-  inline bool safeSet(Indexer_t &lambda, Real const in)
+  inline bool SafeMustSet(Indexer_t &lambda, Real const in)
 
 can modify the values in the ``Indexer_t`` and behave the same way. In this way,
 if a type-based index isn't present in the container, then another index won't
-be overwritten.
+be overwritten. Again, the ``SafeMustSet()`` version will compile-time fail
+or runtime abort if the lambda value can't be modified for the same reasons as
+``SafeMustGet()``.
 
 .. note::
 

From 4c89631785310cc70d28b068bd77d26cc16320e9 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 18:17:39 -0600
Subject: [PATCH 20/30] Get rid of Get and have wrappers use GetSet. Also
 update comments and make const correct

---
 singularity-eos/base/indexable_types.hpp | 69 +++++++-----------------
 1 file changed, 19 insertions(+), 50 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index dd96ced7999..54247da4b7b 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -75,32 +75,6 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeGet(Indexer_t const &lambda, std::size_t
   return false;
 }
 
-// Same as above but causes an error condition (static or dynamic) if the value
-// can't be obtained
-template <AllowedIndexing AI, typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
-                                               std::size_t const idx) {
-  // Error on null pointer
-  PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
-                          "Indexer can't be nullptr");
-
-  // Return type-based index. Static assert that type MUST exist in indexer
-  if constexpr (is_type_indexer_v<Indexer_t>) {
-    static_assert(variadic_utils::is_indexable_v<Indexer_t, T>);
-    return lambda[T{}];
-  }
-
-  // Fall back to numerical indexing if allowed
-  if constexpr (AI == AllowedIndexing::Numeric) {
-    if constexpr (variadic_utils::has_int_index<Indexer_t>::value) {
-      return lambda[idx];
-    }
-  }
-
-  // Something else...
-  PORTABLE_ALWAYS_THROW_OR_ABORT("Cannot obtain value from unknown indexer");
-}
-
 // Break out "Set" functionality from "Get". The original "Get()" did both, but
 // the "safe" version needs to separate that functionality for setting the
 // values in a lambda
@@ -136,10 +110,15 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, std::size_t const
 }
 
 // Same as above but causes an error condition (static or dynamic) if the value
-// can't be obtained
+// can't be obtained. Note that the `decltype(auto)` is intended to preserve the
+// value category of the square bracket operator of the `Indexer_t` type. This
+// allows references to be returned since there is also no possibility of the
+// call doing nothing (i.e. like SafeGet and SafeSet), and thus it can be used
+// for either setting or getting values. This should also allow for `const`
+// correctness downstream in the wrappers where `lambda` is `const`
 template <AllowedIndexing AI, typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx,
-                                               Real const in) {
+PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda,
+                                                            std::size_t const idx) {
   // Error on null pointer
   PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
                           "Indexer can't be nullptr");
@@ -147,15 +126,17 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t co
   // Return type-based index. Static assert that type MUST exist in indexer
   if constexpr (is_type_indexer_v<Indexer_t>) {
     static_assert(variadic_utils::is_indexable_v<Indexer_t, T>);
-    lambda[T{}] = in;
-    return;
+    // Use std::forward to maintain value category for lambda, and use
+    // parentheses to do the same for the output of the lambda[] operation
+    return (std::forward<Indexer_t>(lambda)[T{}]);
   }
 
   // Fall back to numerical indexing if allowed
   if constexpr (AI == AllowedIndexing::Numeric) {
-    if constexpr (variadic_utils::has_int_index<Indexer_t>::value) {
-      lambda[idx] = in;
-      return;
+    if constexpr (variadic_utils::has_int_index_v<Indexer_t>) {
+      // Use std::forward to maintain value category for lambda, and use
+      // parentheses to do the same for the output of the lambda[] operation
+      return (std::forward<Indexer_t>(lambda)[idx]);
     }
   }
 
@@ -197,40 +178,28 @@ PORTABLE_FORCEINLINE_FUNCTION bool SafeSet(Indexer_t &lambda, Real const in) {
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
                                                std::size_t const idx) {
-  return impl::SafeMustGet<impl::AllowedIndexing::Numeric, T>(lambda, idx);
+  return impl::SafeMustGetSet<impl::AllowedIndexing::Numeric, T>(lambda, idx);
 }
 
 // Overload when numerical index isn't provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) {
   std::size_t idx = 0;
-  return impl::SafeMustGet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx);
+  return impl::SafeMustGetSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx);
 }
 
 // Overload when numerical index is provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t const idx,
                                                Real const in) {
-  return impl::SafeMustSet<impl::AllowedIndexing::Numeric, T>(lambda, idx, in);
+  impl::SafeMustGetSet<impl::AllowedIndexing::Numeric, T>(lambda, idx) = in;
 }
 
 // Overload when numerical index isn't provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, Real const in) {
   std::size_t idx = 0;
-  return impl::SafeMustSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx, in);
-}
-
-// NOTE: this Get is "unsafe" because it can allow you to overwrite a type-based
-// index since it automatically falls back to numeric indexing if the type
-// index isn't present.
-template <typename T, typename Indexer_t>
-PORTABLE_FORCEINLINE_FUNCTION auto &Get(Indexer_t &&lambda, std::size_t idx = 0) {
-  if constexpr (variadic_utils::is_indexable_v<Indexer_t, T>) {
-    return lambda[T()];
-  } else {
-    return lambda[idx];
-  }
+  impl::SafeMustGetSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx) = in;
 }
 
 // This is a convenience struct to easily build a small indexer with

From 9d5296068778ffc6b3edd95ed153d4fab5f942af Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 18:18:46 -0600
Subject: [PATCH 21/30] Switch Get for Safe versions and expand tests

---
 test/test_indexable_types.cpp | 74 ++++++++++++++++++++++++++++++++---
 1 file changed, 68 insertions(+), 6 deletions(-)

diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index 8f1be8d913f..ad3277eb3ae 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -21,6 +21,9 @@
 #define CATCH_CONFIG_FAST_COMPILE
 #include <catch2/catch_test_macros.hpp>
 #endif
+#include <catch2/matchers/catch_matchers_floating_point.hpp>
+
+#include <test/eos_unit_test_helpers.hpp>
 
 using namespace singularity::IndexerUtils;
 using namespace singularity::IndexableTypes;
@@ -71,14 +74,41 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
         REQUIRE(lRho == static_cast<Real>(2));
       }
     }
-    WHEN("We use the Get functionality") {
-      // Request a type that exists, but an incorrect index
-      Real Zbar = Get<MeanIonizationState>(lambda, 2);
-      // Request a type that doesn't exist but an index that does
-      Real lRho = Get<LogTemperature>(lambda, 2);
+    WHEN("We use the SafeMustGet functionality") {
+      // Request a type that exists, but no integer index
+      const Real Zbar = SafeMustGet<MeanIonizationState>(lambda);
+      // Request a type that exists, but the wrong integer index
+      Real Zbar_1 = SafeMustGet<MeanIonizationState>(lambda, 2);
       THEN("We get the correct values") {
         REQUIRE(Zbar == static_cast<Real>(0));
-        REQUIRE(lRho == static_cast<Real>(2));
+        REQUIRE_THAT(Zbar, Catch::Matchers::WithinRel(Zbar_1, 1.0e-14));
+      }
+      // This is probably a dumb check since Zbar_1 wasn't a Real* or a Real&
+      AND_THEN("We don't modify the lambda values by changing the local values") {
+        Zbar_1 = 5.3;
+        REQUIRE_THAT(
+            Zbar,
+            Catch::Matchers::WithinRel(
+                SafeMustGet<MeanIonizationState>(lambda),
+                1.0e-14));
+      }
+    }
+    WHEN("We use the SafeMustSet functionality") {
+      // Request a type that exists, but no index
+      Real Zbar = 5.567;
+      SafeMustSet<MeanIonizationState>(lambda, Zbar);
+      // Request a type that exists, but the wrong integer index
+      Real lRho = 1.102;
+      SafeMustSet<LogDensity>(lambda, 0, lRho);
+      THEN("We get the correct values") {
+        REQUIRE_THAT(
+          Zbar,
+          Catch::Matchers::WithinRel(SafeMustGet<MeanIonizationState>(lambda),
+                                     1.0e-14));
+        REQUIRE_THAT(
+          lRho,
+          Catch::Matchers::WithinRel(SafeMustGet<LogDensity>(lambda),
+                                     1.0e-14));
       }
     }
     WHEN("We use the SafeGet functionality") {
@@ -142,6 +172,38 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
       }
     }
   }
+  GIVEN("A normal array that does not support IndexableTypes") {
+    constexpr size_t num_lambda = 4;
+    std::array<Real, num_lambda> lambda{};
+    for (size_t i = 0; i < num_lambda; i++) {
+      lambda[i] = static_cast<Real>(i);
+    }
+    WHEN("We use the SafeMustGet functionality") {
+      THEN("The type-based index is ignored and only the integer index is used") {
+        for (size_t i = 0; i < num_lambda; i++) {
+          INFO("i: " << i);
+          const Real val = SafeMustGet<LogTemperature>(lambda, i);
+          CHECK_THAT(lambda[i], Catch::Matchers::WithinRel(val, 1.0e-14));
+        }
+      }
+      THEN("If an integer index isn't provided, an exception is thrown") {
+        REQUIRE_MAYBE_THROWS(SafeMustGet<void>(lambda));
+      }
+    }
+    WHEN("We use the SafeMustSet functionality") {
+      THEN("The type-based index is ignored and only the integer index is used") {
+        for (size_t i = 0; i < num_lambda; i++) {
+          INFO("i: " << i);
+          const Real val = i * i;
+          SafeMustSet<LogDensity>(lambda, i, val);
+          CHECK_THAT(lambda[i], Catch::Matchers::WithinRel(val, 1.0e-14));
+        }
+      }
+      THEN("If an integer index isn't provided, an exception is thrown") {
+        REQUIRE_MAYBE_THROWS(SafeMustSet<void>(lambda, 1.0));
+      }
+    }
+  }
 }
 
 SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {

From 0feb8e437809a030873b07223e332238dcfb911a Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 18:23:33 -0600
Subject: [PATCH 22/30] Remove last Get in favor of SafeSet

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

diff --git a/singularity-eos/eos/eos_helmholtz.hpp b/singularity-eos/eos/eos_helmholtz.hpp
index 566de55aa7c..665bf5925ed 100644
--- a/singularity-eos/eos/eos_helmholtz.hpp
+++ b/singularity-eos/eos/eos_helmholtz.hpp
@@ -838,7 +838,7 @@ Helmholtz::FillEos(Real &rho, Real &temp, Real &energy, Real &press, Real &cv, R
     temp = math_utils::pow10(lT);
   } else {
     lT = std::log10(temp);
-    IndexerUtils::Get<IndexableTypes::LogTemperature>(lambda, Lambda::lT) = lT;
+    IndexerUtils::SafeSet<IndexableTypes::LogTemperature>(lambda, Lambda::lT, lT);
   }
   Real p[NDERIV], e[NDERIV], s[NDERIV], etaele[NDERIV], nep[NDERIV];
   GetFromDensityLogTemperature_(rho, temp, abar, zbar, ye, ytot, ywot, De, lDe, p, e, s,

From c245728679055e0095251738a332abfb0bdd13a7 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 19:50:01 -0600
Subject: [PATCH 23/30] Remove a few more instances of Get in favor of the Safe
 variety for IndexableTypes

---
 test/test_indexable_types.cpp  | 10 ----------
 test/test_spiner_transform.cpp |  4 +---
 2 files changed, 1 insertion(+), 13 deletions(-)

diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index ad3277eb3ae..d1949f8a3f6 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -212,16 +212,6 @@ SCENARIO("IndexableTypes and ManualLambda", "[IndexableTypes]") {
     for (std::size_t i = 0; i < lambda.length; ++i) {
       lambda[i] = static_cast<Real>(i);
     }
-    WHEN("We use the Get functionality") {
-      // Request a type that exists, but an incorrect index
-      Real Zbar = Get<MeanIonizationState>(lambda, 2);
-      // Request a type that doesn't exist but an index that does
-      Real lRho = Get<LogTemperature>(lambda, 2);
-      THEN("We get the correct values") {
-        REQUIRE(Zbar == static_cast<Real>(0));
-        REQUIRE(lRho == static_cast<Real>(2));
-      }
-    }
     WHEN("We use the SafeGet functionality") {
       constexpr Real unmodified = -1.0;
       Real destination = unmodified;
diff --git a/test/test_spiner_transform.cpp b/test/test_spiner_transform.cpp
index 960f608114d..0baeedaddf6 100644
--- a/test/test_spiner_transform.cpp
+++ b/test/test_spiner_transform.cpp
@@ -102,9 +102,7 @@ struct TestDataContainer {
     const Real lRho = spiner_common::to_log(rho, lRhoOffset);
     const Real lE = spiner_common::to_log(sie, lEOffset);
     // If we wanted to mock the lambda write-back (optional)
-    if (!variadic_utils::is_nullptr(lambda)) {
-      IndexerUtils::Get<IndexableTypes::LogDensity>(lambda, 0) = lRho;
-    }
+    IndexerUtils::SafeSet<IndexableTypes::LogDensity>(lambda, 0, lRho);
     return db.interpToReal(lRho, lE);
   }
 

From 734d55b12e40b10dce4b63f9f8c36c0ffb3a9cac Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 19:50:29 -0600
Subject: [PATCH 24/30] Clang format

---
 test/test_indexable_types.cpp | 17 +++++------------
 1 file changed, 5 insertions(+), 12 deletions(-)

diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index d1949f8a3f6..591a98b98eb 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -86,11 +86,8 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
       // This is probably a dumb check since Zbar_1 wasn't a Real* or a Real&
       AND_THEN("We don't modify the lambda values by changing the local values") {
         Zbar_1 = 5.3;
-        REQUIRE_THAT(
-            Zbar,
-            Catch::Matchers::WithinRel(
-                SafeMustGet<MeanIonizationState>(lambda),
-                1.0e-14));
+        REQUIRE_THAT(Zbar, Catch::Matchers::WithinRel(
+                               SafeMustGet<MeanIonizationState>(lambda), 1.0e-14));
       }
     }
     WHEN("We use the SafeMustSet functionality") {
@@ -101,14 +98,10 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
       Real lRho = 1.102;
       SafeMustSet<LogDensity>(lambda, 0, lRho);
       THEN("We get the correct values") {
+        REQUIRE_THAT(Zbar, Catch::Matchers::WithinRel(
+                               SafeMustGet<MeanIonizationState>(lambda), 1.0e-14));
         REQUIRE_THAT(
-          Zbar,
-          Catch::Matchers::WithinRel(SafeMustGet<MeanIonizationState>(lambda),
-                                     1.0e-14));
-        REQUIRE_THAT(
-          lRho,
-          Catch::Matchers::WithinRel(SafeMustGet<LogDensity>(lambda),
-                                     1.0e-14));
+            lRho, Catch::Matchers::WithinRel(SafeMustGet<LogDensity>(lambda), 1.0e-14));
       }
     }
     WHEN("We use the SafeGet functionality") {

From efd9a9f088bb7421fdea3b7c40e1f241562e443d Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 20:05:32 -0600
Subject: [PATCH 25/30] Whoops... void was a bad choice for a type index

---
 test/test_indexable_types.cpp | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index 591a98b98eb..2b8a87ef5cc 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -180,7 +180,7 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
         }
       }
       THEN("If an integer index isn't provided, an exception is thrown") {
-        REQUIRE_MAYBE_THROWS(SafeMustGet<void>(lambda));
+        REQUIRE_MAYBE_THROWS(SafeMustGet<LogTemperature>(lambda));
       }
     }
     WHEN("We use the SafeMustSet functionality") {
@@ -193,7 +193,7 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
         }
       }
       THEN("If an integer index isn't provided, an exception is thrown") {
-        REQUIRE_MAYBE_THROWS(SafeMustSet<void>(lambda, 1.0));
+        REQUIRE_MAYBE_THROWS(SafeMustSet<LogTemperature>(lambda, 1.0));
       }
     }
   }

From d8dbcccee5e8836b91b31a42f7f537f8c6650aea Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 20:26:22 -0600
Subject: [PATCH 26/30] Let's try an unreachable return to make decltype(auto)
 happy

---
 singularity-eos/base/indexable_types.hpp | 5 +++++
 1 file changed, 5 insertions(+)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 54247da4b7b..9239c527af2 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -142,6 +142,11 @@ PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda,
 
   // Something else...
   PORTABLE_ALWAYS_THROW_OR_ABORT("Cannot obtain value from unknown indexer");
+
+  // This code is unreachable, but we need a concrete return type so that
+  // everything downstream is happy
+  Real ret = 0;
+  return ret;
 }
 
 } // namespace impl

From bfacaacbd0fed10aa4c4903d1ca3800ae4e223d8 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 20:56:25 -0600
Subject: [PATCH 27/30] Add dependent_false_v for if constexpr static asserts

---
 singularity-eos/base/variadic_utils.hpp | 10 ++++++++++
 1 file changed, 10 insertions(+)

diff --git a/singularity-eos/base/variadic_utils.hpp b/singularity-eos/base/variadic_utils.hpp
index d2aa3aba3cc..8b30c0cb5b2 100644
--- a/singularity-eos/base/variadic_utils.hpp
+++ b/singularity-eos/base/variadic_utils.hpp
@@ -24,6 +24,16 @@ namespace variadic_utils {
 // Some generic variatic utilities
 // ======================================================================
 
+// Template parameter dependent boolean suitable for causing static_assert
+// errors within `if constexpr` branches. Essentially the issue is that if the
+// static_assert _always_ evaluates to false, then it will _always_ cause a
+// compile time error even if that branch of the code will never be reached.
+// Making the evaluation (superficially) dependent on the template deduction
+// causes it to be evaluated after the `if constexpr` branching has already been
+// determined. See https://en.cppreference.com/w/cpp/language/if.html#Constexpr_if
+template <class>
+inline constexpr bool dependent_false_v = false;
+
 // Useful for generating nullptr of a specific pointer type
 template <typename T>
 inline constexpr T *np() {

From 36d42868b0db1a5fd80a1ba21882b1d033ede257 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 21:00:55 -0600
Subject: [PATCH 28/30] Move throw into SafeGet/SafeSet wrappers and provide
 helpful compile time error if neither indexable types or integer indexing are
 allowed

---
 singularity-eos/base/indexable_types.hpp | 43 ++++++++++++++----------
 1 file changed, 25 insertions(+), 18 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 9239c527af2..6c56ac1f865 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -123,30 +123,25 @@ PORTABLE_FORCEINLINE_FUNCTION decltype(auto) SafeMustGetSet(Indexer_t &&lambda,
   PORTABLE_ALWAYS_REQUIRE(!variadic_utils::is_nullptr(lambda),
                           "Indexer can't be nullptr");
 
-  // Return type-based index. Static assert that type MUST exist in indexer
   if constexpr (is_type_indexer_v<Indexer_t>) {
+    // Return type-based index. Static assert that type MUST exist in indexer
     static_assert(variadic_utils::is_indexable_v<Indexer_t, T>);
     // Use std::forward to maintain value category for lambda, and use
     // parentheses to do the same for the output of the lambda[] operation
     return (std::forward<Indexer_t>(lambda)[T{}]);
+  } else if constexpr (AI == AllowedIndexing::Numeric) {
+    // Fall back to numerical indexing if allowed
+    static_assert(variadic_utils::has_int_index_v<Indexer_t>);
+    // Use std::forward to maintain value category for lambda, and use
+    // parentheses to do the same for the output of the lambda[] operation
+    return (std::forward<Indexer_t>(lambda)[idx]);
+  } else {
+    // Something else that can't be compiled...
+    static_assert(variadic_utils::dependent_false_v<Indexer_t>,
+                  "Indexer must either be designated as type-based through a "
+                  "`is_type_indexable` boolean data member or SafeGet/SafeSet function "
+                  "must be called with a numerical index");
   }
-
-  // Fall back to numerical indexing if allowed
-  if constexpr (AI == AllowedIndexing::Numeric) {
-    if constexpr (variadic_utils::has_int_index_v<Indexer_t>) {
-      // Use std::forward to maintain value category for lambda, and use
-      // parentheses to do the same for the output of the lambda[] operation
-      return (std::forward<Indexer_t>(lambda)[idx]);
-    }
-  }
-
-  // Something else...
-  PORTABLE_ALWAYS_THROW_OR_ABORT("Cannot obtain value from unknown indexer");
-
-  // This code is unreachable, but we need a concrete return type so that
-  // everything downstream is happy
-  Real ret = 0;
-  return ret;
 }
 
 } // namespace impl
@@ -189,6 +184,12 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
 // Overload when numerical index isn't provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) {
+  if constexpr (!is_type_indexer_v<Indexer_t>) {
+    // Ill-formed since no index provided
+    PORTABLE_ALWAYS_THROW_OR_ABORT(
+        "SafeMustGet: Type-based indexing is required, but lambda does does not contain "
+        "is_type_indexable boolean data member");
+  }
   std::size_t idx = 0;
   return impl::SafeMustGetSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx);
 }
@@ -203,6 +204,12 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t co
 // Overload when numerical index isn't provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, Real const in) {
+  if constexpr (!is_type_indexer_v<Indexer_t>) {
+    // Ill-formed since no index provided
+    PORTABLE_ALWAYS_THROW_OR_ABORT(
+        "SafeMustSet: Type-based indexing is required, but lambda does does not contain "
+        "is_type_indexable boolean data member");
+  }
   std::size_t idx = 0;
   impl::SafeMustGetSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx) = in;
 }

From c783623c70cae491b424fec48ef6fecc0a4d01e1 Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 21:12:20 -0600
Subject: [PATCH 29/30] Can't test for a runtime throw when compile-time error
 will be used

---
 singularity-eos/base/indexable_types.hpp | 12 ------------
 test/test_indexable_types.cpp            |  6 ------
 2 files changed, 18 deletions(-)

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
index 6c56ac1f865..9f6d4737f0e 100644
--- a/singularity-eos/base/indexable_types.hpp
+++ b/singularity-eos/base/indexable_types.hpp
@@ -184,12 +184,6 @@ PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda,
 // Overload when numerical index isn't provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION Real SafeMustGet(Indexer_t const &lambda) {
-  if constexpr (!is_type_indexer_v<Indexer_t>) {
-    // Ill-formed since no index provided
-    PORTABLE_ALWAYS_THROW_OR_ABORT(
-        "SafeMustGet: Type-based indexing is required, but lambda does does not contain "
-        "is_type_indexable boolean data member");
-  }
   std::size_t idx = 0;
   return impl::SafeMustGetSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx);
 }
@@ -204,12 +198,6 @@ PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, std::size_t co
 // Overload when numerical index isn't provided
 template <typename T, typename Indexer_t>
 PORTABLE_FORCEINLINE_FUNCTION void SafeMustSet(Indexer_t &lambda, Real const in) {
-  if constexpr (!is_type_indexer_v<Indexer_t>) {
-    // Ill-formed since no index provided
-    PORTABLE_ALWAYS_THROW_OR_ABORT(
-        "SafeMustSet: Type-based indexing is required, but lambda does does not contain "
-        "is_type_indexable boolean data member");
-  }
   std::size_t idx = 0;
   impl::SafeMustGetSet<impl::AllowedIndexing::TypeOnly, T>(lambda, idx) = in;
 }
diff --git a/test/test_indexable_types.cpp b/test/test_indexable_types.cpp
index 2b8a87ef5cc..9b83e251688 100644
--- a/test/test_indexable_types.cpp
+++ b/test/test_indexable_types.cpp
@@ -179,9 +179,6 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
           CHECK_THAT(lambda[i], Catch::Matchers::WithinRel(val, 1.0e-14));
         }
       }
-      THEN("If an integer index isn't provided, an exception is thrown") {
-        REQUIRE_MAYBE_THROWS(SafeMustGet<LogTemperature>(lambda));
-      }
     }
     WHEN("We use the SafeMustSet functionality") {
       THEN("The type-based index is ignored and only the integer index is used") {
@@ -192,9 +189,6 @@ SCENARIO("IndexableTypes and VariadicIndexer", "[IndexableTypes][VariadicIndexer
           CHECK_THAT(lambda[i], Catch::Matchers::WithinRel(val, 1.0e-14));
         }
       }
-      THEN("If an integer index isn't provided, an exception is thrown") {
-        REQUIRE_MAYBE_THROWS(SafeMustSet<LogTemperature>(lambda, 1.0));
-      }
     }
   }
 }

From 232ecca8fae12dedaa746a7471ba632c164ddb6a Mon Sep 17 00:00:00 2001
From: Jeffrey H Peterson <jhp@lanl.gov>
Date: Tue, 16 Sep 2025 22:36:28 -0600
Subject: [PATCH 30/30] Small doc tweaks

---
 doc/sphinx/src/using-eos.rst | 21 +++++++++++++++------
 1 file changed, 15 insertions(+), 6 deletions(-)

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
index bcdce9792fd..326c4173367 100644
--- a/doc/sphinx/src/using-eos.rst
+++ b/doc/sphinx/src/using-eos.rst
@@ -758,11 +758,20 @@ doesn't accept type-based indexing. If the ``Indexer_t`` **does** accept
 type-based indexing but **doesn't** have the requested index, then the
 ``out`` value is not updated. The same is true for when ``Indexer_t`` is the
 ``nullptr``. The overload that doesn't take a numerical index will *only*
-return the value at a type-based index. The ``SafeMustGet()`` version
-is intended to generate errors if a value can't be retrieved. In the case of
-type-based indexing, the error will be at compile time if the type isn't located
-in the indexer. A runtime abort will occur if either the null pointer is passed
-or if integer indexing isn't allowed for some reason.
+return the value at a type-based index.
+
+The ``SafeMustGet()`` version is intended to generate errors if a value can't be
+retrieved with four types of errors that can occur:
+
+1. If a null pointer is passed as the indexer, a runtime abort or exception will
+   occur
+2. If a type-based indexer is passed (i.e. one with the ``constexpr static bool``
+   member ``is_type_indexable = true``), but the type doesn't exist then a **static**
+   assertion will fail
+3. If one of the overloads is used where a ``std::size_t`` index is provided, but
+   the indexer can't accept integer indexing, then a **static** assertion will fail
+4. If the indexer can't use type-based indexing but an ``std::size_t`` index
+   wasn't provided, then a **static** assertion will fail
 
 Similarly, the functions
 
@@ -791,7 +800,7 @@ Similarly, the functions
   inline bool SafeMustSet(Indexer_t &lambda, Real const in)
 
 can modify the values in the ``Indexer_t`` and behave the same way. In this way,
-if a type-based index isn't present in the container, then another index won't
+if a type-based index isn't present in the container, then a different index won't
 be overwritten. Again, the ``SafeMustSet()`` version will compile-time fail
 or runtime abort if the lambda value can't be modified for the same reasons as
 ``SafeMustGet()``.