lanl · Yurlungur · Apr 10, 2026 · Apr 7, 2026 · Apr 7, 2026 · Apr 7, 2026
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -5,6 +5,7 @@
 ### Added (new features/APIs/variables/...)
 - [[PR623]](https://github.com/lanl/singularity-eos/pull/623) Expanded the sesame2spiner syntax to support multiple material definitions in one input file.
 - [[PR618]](https://github.com/lanl/singularity-eos/pull/618) Add PTDerivativesFromPreferred for computing derivatives of a mixture  in a cell
+- [[PR630]](https://github.com/lanl/singularity-eos/pull/630) Thread execution spaces through vector API.
 
 ### Fixed (Repair bugs, etc)
 

diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst
@@ -475,6 +475,23 @@ The vector API is templated to accept accessors. We do note, however,
 that vectorization may suffer if your underlying data structure is not
 contiguous in memory.
 
+Execution spaces for the vector API
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+``singularity-eos`` borrows from the `Kokkos`_ notion of an *execution
+space*. Each vector API call accepts an execution space as an optional
+first argument, telling ``singularity-eos`` where to launch work,
+either on CPU or GPU (for example). These are implemented via
+`ports-of-call`_. By default, the two execution spaces available are
+``PortsOfCall::Exec::Host`` and ``PortsOfCall::Exec::Device``. If you
+build with the serial backend, these are equivalent. If you build with
+``Kokkos`` enabled, you may also pass in any valid ``Kokkos``
+execution space.
+
+.. _Kokkos: https://kokkos.org/kokkos-core-wiki/
+
+.. _ports-of-call: https://lanl.github.io/ports-of-call/main/index.html
+
 .. _eospac_vector:
 
 EOSPAC Vector Functions
@@ -506,6 +523,12 @@ available member function.
    eos.TemperatureFromDensityInternalEnergy(density.data(), energy.data(), temperature.data(),
                                             scratch.data(), density.size());
 
+.. warning::
+
+  EOSPAC does **not** support selection of execution spaces. You will
+  get whatever execution space (likely host) that the EOSPAC backend
+  has been configured for.
+
 The Evaluate Methods
 ~~~~~~~~~~~~~~~~~~~~~~
 

diff --git a/plan_histories/MR630-2026-04-09-modifier-vector-macro-proposal.md b/plan_histories/MR630-2026-04-09-modifier-vector-macro-proposal.md
@@ -0,0 +1,186 @@
+# MR630 Modifier Vector Macro Proposal
+
+## Goal
+
+Reduce the line count in `singularity-eos/eos/modifiers/*.hpp` for the raw-pointer
+vector API, especially the overloads that only forward to
+`PortsOfCall::Exec::Device()`.
+
+## Current Situation
+
+The modified vector API in the modifier headers now has two overload families for each
+raw-pointer vector method:
+
+1. A `Space` overload:
+
+```cpp
+template <typename Space, typename LambdaIndexer,
+          typename EnableIfSpace =
+              std::enable_if_t<!variadic_utils::has_int_index_v<Space>>>
+inline void NAME(const Space &s, ...);
+```
+
+2. A default-device wrapper:
+
+```cpp
+template <typename LambdaIndexer>
+inline void NAME(...) const {
+  NAME(PortsOfCall::Exec::Device(), ...);
+}
+```
+
+This repeats across:
+
+- `scaled_eos.hpp`
+- `shifted_eos.hpp`
+- `eos_unitsystem.hpp`
+- `ramps_eos.hpp`
+- `floored_energy.hpp`
+
+## Important Constraint
+
+The default-device wrappers cannot simply be deleted while preserving the current API.
+The inherited wrappers in `eos/eos_base.hpp` do not automatically route into the
+modifier-specific raw-pointer `Space` overloads. If we want to keep the existing
+call surface, the realistic near-term win is to generate these wrappers, not remove the
+behavior.
+
+## Recommended Macro Scheme
+
+Introduce a small modifier-only macro header, for example:
+
+- `singularity-eos/eos/modifiers/modifier_vector_macros.hpp`
+
+That header would define the default-device wrappers once.
+
+### 2-in / 1-out wrapper
+
+```cpp
+#define SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(NAME)                                           \
+  template <typename LambdaIndexer>                                                      \
+  inline void NAME(const Real *in1, const Real *in2, Real *out, Real *scratch,           \
+                   const int num, LambdaIndexer &&lambdas,                               \
+                   Transform &&transform = Transform()) const {                          \
+    NAME(PortsOfCall::Exec::Device(), in1, in2, out, scratch, num,                       \
+         std::forward<LambdaIndexer>(lambdas), std::forward<Transform>(transform));      \
+  }
+```
+
+### 1-in / 1-out wrapper
+
+```cpp
+#define SG_MODIFIER_DEVICE_WRAP_1IN_1OUT(NAME)                                           \
+  template <typename LambdaIndexer>                                                      \
+  inline void NAME(const Real *in1, Real *out, Real *scratch, const int num,             \
+                   LambdaIndexer &&lambdas,                                              \
+                   Transform &&transform = Transform()) const {                          \
+    NAME(PortsOfCall::Exec::Device(), in1, out, scratch, num,                            \
+         std::forward<LambdaIndexer>(lambdas), std::forward<Transform>(transform));      \
+  }
+```
+
+### Aggregate wrapper list
+
+```cpp
+#define SG_MODIFIER_DEVICE_WRAP_ALL()                                                    \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(TemperatureFromDensityInternalEnergy)                  \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(PressureFromDensityTemperature)                        \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(PressureFromDensityInternalEnergy)                     \
+  SG_MODIFIER_DEVICE_WRAP_1IN_1OUT(MinInternalEnergyFromDensity)                          \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(SpecificHeatFromDensityTemperature)                    \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(SpecificHeatFromDensityInternalEnergy)                 \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(BulkModulusFromDensityTemperature)                     \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(BulkModulusFromDensityInternalEnergy)                  \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(GruneisenParamFromDensityTemperature)                  \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(GruneisenParamFromDensityInternalEnergy)               \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(InternalEnergyFromDensityTemperature)                  \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(EntropyFromDensityTemperature)                         \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(EntropyFromDensityInternalEnergy)
+```
+
+## How This Would Be Used
+
+Each modifier would keep only its `Space` overload implementations. At the end of the
+vector API block, it would add:
+
+```cpp
+SG_MODIFIER_DEVICE_WRAP_ALL()
+```
+
+That removes most of the duplicate default-device wrapper code while preserving the
+current behavior.
+
+## Second Stage for Simple Forwarders
+
+For `scaled_eos.hpp` and `eos_unitsystem.hpp`, the `Space` overload bodies are also
+highly repetitive:
+
+- update `transform`
+- call `t_.NAME(s, ...)`
+
+Those two files could use a second macro family such as:
+
+```cpp
+#define SG_MODIFIER_FORWARD_2IN_1OUT(NAME, PREPARE)                                      \
+  template <typename Space, typename LambdaIndexer,                                      \
+            typename EnableIfSpace =                                                     \
+                std::enable_if_t<!variadic_utils::has_int_index_v<Space>>>               \
+  inline void NAME(const Space &s, const Real *in1, const Real *in2, Real *out,          \
+                   Real *scratch, const int num, LambdaIndexer &&lambdas,                \
+                   Transform &&transform = Transform()) const {                          \
+    PREPARE                                                                              \
+    t_.NAME(s, in1, in2, out, scratch, num,                                              \
+            std::forward<LambdaIndexer>(lambdas), std::forward<Transform>(transform));   \
+  }                                                                                      \
+  SG_MODIFIER_DEVICE_WRAP_2IN_1OUT(NAME)
+```
+
+Example usage:
+
+```cpp
+SG_MODIFIER_FORWARD_2IN_1OUT(
+    TemperatureFromDensityInternalEnergy,
+    transform.x.apply(scale_);
+    transform.y.apply(inv_scale_);)
+```
+
+## Where Not to Overuse Macros
+
+I would not force the body macros into:
+
+- `shifted_eos.hpp`
+- `floored_energy.hpp`
+- `ramps_eos.hpp`
+
+Those files contain real per-method logic:
+
+- scratch-buffer preprocessing and postprocessing
+- local `portableFor` loops
+- ramp-pressure overrides
+
+For those files, the better tradeoff is:
+
+- keep explicit `Space` overload bodies
+- macro-generate only the default-device wrappers
+
+## Alternative If True Elimination Is Desired
+
+If the actual goal is to eliminate per-file default-device overloads entirely, use a CRTP
+mixin instead of macros. Each modifier would define only the `Space` overloads, and the
+mixin would provide the no-space wrappers once.
+
+That is conceptually cleaner than making the macro scheme more aggressive, but it is a
+larger refactor than the macro-only approach above.
+
+## Recommendation
+
+Preferred implementation order:
+
+1. Add a shared modifier macro header that generates the default-device wrappers.
+2. Convert all five modified modifier headers to use that wrapper macro set.
+3. Optionally add the forwarding-body macros for `ScaledEOS` and `EOSUnitSystem`.
+4. Leave `ShiftedEOS`, `FlooredEnergy`, and `BilinearRampEOS` with explicit `Space`
+   overload bodies.
+
+This gives most of the line-count reduction without hiding the custom vector logic that
+still needs to be readable.
diff --git a/python/module.hpp b/python/module.hpp
@@ -22,6 +22,7 @@
 #include <pybind11/pybind11.h>
 #include <pybind11/numpy.h>
 
+#include <ports-of-call/portability.hpp>
 #include <ports-of-call/portable_errors.hpp>
 
 #include <singularity-eos/base/variadic_utils.hpp>
@@ -104,10 +105,11 @@ void func(const T & self, py::array_t<Real> a, py::array_t<Real> b,
   auto av = a.unchecked<1>();                                                     \
   auto bv = b.unchecked<1>();                                                     \
   auto outv = PyArrayHelper(out.mutable_unchecked<1>());                          \
+  PortsOfCall::Exec::Host host;                                                   \
   if(lambdas_info.shape[1] > 0) {                                                 \
-    self.func(av, bv, outv, a.size(), LambdaHelper(lambdas));                     \
+    self.func(host, av, bv, outv, a.size(), LambdaHelper(lambdas));               \
   } else {                                                                        \
-    self.func(av, bv, outv, a.size(), NoLambdaHelper());                          \
+    self.func(host, av, bv, outv, a.size(), NoLambdaHelper());                    \
   }                                                                               \
 }                                                                                 \
                                                                                   \
@@ -122,15 +124,16 @@ void func##WithScratch(const T & self, py::array_t<Real> a, py::array_t<Real> b,
   auto bv = PyArrayHelper(b.mutable_unchecked<1>());                              \
   auto outv = PyArrayHelper(out.mutable_unchecked<1>());                          \
   auto scrv = PyArrayHelper(scratch.mutable_unchecked<1>());                      \
+  PortsOfCall::Exec::Host host;                                                   \
   PORTABLE_REQUIRE(av.IsContiguous(), "arrays w/ scratch must be contiguous");    \
   PORTABLE_REQUIRE(bv.IsContiguous(), "arrays w/ scratch must be contiguous");    \
   PORTABLE_REQUIRE(outv.IsContiguous(), "arrays w/ scratch must be contiguous");  \
   PORTABLE_REQUIRE(scrv.IsContiguous(), "arrays w/ scratch must be contiguous");  \
   if(lambdas_info.shape[1] > 0) {                                                 \
-    self.func(av.data(), bv.data(), outv.data(), scrv.data(), a.size(),           \
+    self.func(host, av.data(), bv.data(), outv.data(), scrv.data(), a.size(),     \
               LambdaHelper(lambdas));                                             \
   } else {                                                                        \
-    self.func(av.data(), bv.data(), outv.data(), scrv.data(), a.size(),           \
+    self.func(host, av.data(), bv.data(), outv.data(), scrv.data(), a.size(),     \
               NoLambdaHelper());                                                  \
   }                                                                               \
 }                                                                                 \
@@ -141,7 +144,8 @@ void func##NoLambda(const T & self, py::array_t<Real> a, py::array_t<Real> b,
   auto av = a.unchecked<1>();                                                     \
   auto bv = b.unchecked<1>();                                                     \
   auto outv = PyArrayHelper(out.mutable_unchecked<1>());                          \
-  self.func(av, bv, outv, a.size(), NoLambdaHelper());                            \
+  PortsOfCall::Exec::Host host;                                                   \
+  self.func(host, av, bv, outv, a.size(), NoLambdaHelper());                      \
 }                                                                                 \
                                                                                   \
 template<typename T>                                                              \
@@ -151,11 +155,12 @@ void func##NoLambdaWithScratch(const T & self, py::array_t<Real> a, py::array_t<
   auto bv = PyArrayHelper(b.mutable_unchecked<1>());                              \
   auto outv = PyArrayHelper(out.mutable_unchecked<1>());                          \
   auto scrv = PyArrayHelper(scratch.mutable_unchecked<1>());                      \
+  PortsOfCall::Exec::Host host;                                                   \
   PORTABLE_REQUIRE(av.IsContiguous(), "arrays w/ scratch must be contiguous");    \
   PORTABLE_REQUIRE(bv.IsContiguous(), "arrays w/ scratch must be contiguous");    \
   PORTABLE_REQUIRE(outv.IsContiguous(), "arrays w/ scratch must be contiguous");  \
   PORTABLE_REQUIRE(scrv.IsContiguous(), "arrays w/ scratch must be contiguous");  \
-  self.func(av.data(), bv.data(), outv.data(), scrv.data(), a.size(),             \
+  self.func(host, av.data(), bv.data(), outv.data(), scrv.data(), a.size(),       \
             NoLambdaHelper());                                                    \
 }
 
@@ -260,11 +265,13 @@ struct VectorFunctions {
       auto bmodsv = PyArrayHelper(bmods.mutable_unchecked<1>());
 
       if(lambdas_info.shape[1] > 0) {
-        self.FillEos(rhosv, temperaturesv,
+        self.FillEos(PortsOfCall::Exec::Host(),
+                     rhosv, temperaturesv,
                      siesv, pressuresv, cvsv,
                      bmodsv, rhos.size(), output, LambdaHelper(lambdas));
       } else {
-        self.FillEos(rhosv, temperaturesv,
+        self.FillEos(PortsOfCall::Exec::Host(),
+                     rhosv, temperaturesv,
                      siesv, pressuresv, cvsv,
                      bmodsv, rhos.size(), output, NoLambdaHelper());
       }
@@ -282,7 +289,8 @@ struct VectorFunctions {
       auto pressuresv = PyArrayHelper(pressures.mutable_unchecked<1>());
       auto cvsv = PyArrayHelper(cvs.mutable_unchecked<1>());
       auto bmodsv = PyArrayHelper(bmods.mutable_unchecked<1>());
-      self.FillEos(rhosv, temperaturesv,
+      self.FillEos(PortsOfCall::Exec::Host(),
+                   rhosv, temperaturesv,
                     siesv, pressuresv, cvsv,
                    bmodsv, rhos.size(), output, NoLambdaHelper());
     }, py::arg("rhos"), py::arg("temperatures"), py::arg("sies"), py::arg("pressures"), py::arg("cvs"), py::arg("bmods"), py::arg("output"));
@@ -334,12 +342,14 @@ struct VectorFunctions<T,true> {
       auto scrv = PyArrayHelper(scratch.mutable_unchecked<1>());
 
       if(lambdas_info.shape[1] > 0) {
-        self.FillEos(rhosv.data(), temperaturesv.data(),
+        self.FillEos(PortsOfCall::Exec::Host(),
+                     rhosv.data(), temperaturesv.data(),
                      siesv.data(), pressuresv.data(), cvsv.data(),
                      bmodsv.data(), scrv.data(), rhos.size(), output,
                      LambdaHelper(lambdas));
       } else {
-        self.FillEos(rhosv.data(), temperaturesv.data(),
+        self.FillEos(PortsOfCall::Exec::Host(),
+                     rhosv.data(), temperaturesv.data(),
                      siesv.data(), pressuresv.data(), cvsv.data(),
                      bmodsv.data(), scrv.data(), rhos.size(), output,
                      NoLambdaHelper());
@@ -361,7 +371,8 @@ struct VectorFunctions<T,true> {
       auto bmodsv = PyArrayHelper(bmods.mutable_unchecked<1>());
       auto scrv = PyArrayHelper(scratch.mutable_unchecked<1>());
 
-      self.FillEos(rhosv.data(), temperaturesv.data(),
+      self.FillEos(PortsOfCall::Exec::Host(),
+                   rhosv.data(), temperaturesv.data(),
                    siesv.data(), pressuresv.data(), cvsv.data(),
                    bmodsv.data(), scrv.data(), rhos.size(),
                    output, NoLambdaHelper());

diff --git a/sesame2spiner/io_eospac.cpp b/sesame2spiner/io_eospac.cpp
@@ -333,7 +333,6 @@ bool eosMassFraction(int matid, const Bounds &lRhoBounds, const Bounds &lTBounds
                      std::string &phase_names, Verbosity eospacWarn) {
   using namespace EospacWrapper;
 
-  EOS_INTEGER errorCode = EOS_OK;
   constexpr int NT = 2;
   EOS_INTEGER tableHandle[NT];
   EOS_INTEGER tableType[NT] = {EOS_M_DT, EOS_Comment};

diff --git a/singularity-eos/CMakeLists.txt b/singularity-eos/CMakeLists.txt
@@ -57,6 +57,7 @@ register_headers(
     eos/eos_helmholtz.hpp
     eos/eos_sap_polynomial.hpp
     eos/eos_type_lists.hpp
+    eos/modifiers/modifier_vector_macros.hpp
     eos/modifiers/relativistic_eos.hpp
     eos/modifiers/scaled_eos.hpp
     eos/modifiers/ramps_eos.hpp

diff --git a/singularity-eos/base/indexable_types.hpp b/singularity-eos/base/indexable_types.hpp
@@ -302,8 +302,8 @@ class VariadicIndexerBase {
   PORTABLE_FORCEINLINE_FUNCTION
   VariadicIndexerBase(const Data_t &data) : data_(data) {}
 
-  template <typename T,
-            typename = std::enable_if_t<variadic_utils::contains<T, Ts...>::value>>
+  template <typename T, typename EnableIfContained =
+                            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];
@@ -312,8 +312,8 @@ class VariadicIndexerBase {
   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>>
+  template <typename T, typename EnableIfContained =
+                            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];