Thread execution spaces through vector API

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)
 

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
 ~~~~~~~~~~~~~~~~~~~~~~
 

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

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

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