Merge pull request #306 from lanl/jdolence/generic_evaluate

Add a generic evaluator function

CHANGELOG.md

@@ -22,6 +22,7 @@
 - [[PR232]](https://github.com/lanl/singularity-eos/pull/228) Fixed uninitialized cmake path variables
 
 ### Added (new features/APIs/variables/...)
+- [[PR306]](https://github.com/lanl/singularity-eos/pull/306) Added generic Evaluate method
 - [[PR304]](https://github.com/lanl/singularity-eos/pull/304) added a Newton-Raphson root find for use with the Helmholtz EoS
 - [[PR265]](https://github.com/lanl/singularity-eos/pull/265) Add missing UnitSystem modifier combinations to variant and EOSPAC
 - [[PR279]](https://github.com/lanl/singularity-eos/pull/279) added noble-abel EoS

doc/sphinx/src/using-eos.rst

@@ -35,6 +35,10 @@ conditions. The type of parallelism used depends on how
 ``singularity-eos`` is compiled. If the ``Kokkos`` backend is used,
 any parallel dispatch supported by ``Kokkos`` is supported.
 
+A more generic version of the vector calls exists in the ``Evaluate``
+method, which allows the user to specify arbitrary parallel dispatch
+models by writing their own loops. See the relevant section below.
+
 .. _variant section:
 
 Variants
@@ -216,6 +220,121 @@ available member function.
    eos.TemperatureFromDensityInternalEnergy(density.data(), energy.data(), temperature.data(),
                                             scratch.data(), density.size());
 
+The Evaluate Method
+~~~~~~~~~~~~~~~~~~~
+
+A special call related to the vector calls is the ``Evaluate``
+method. The ``Evaluate`` method requests the EOS object to evaluate
+almost arbitrary code, but in a way where the type of the underlying
+EOS object is resolved *before* this arbitrary code is evaluated. This
+means the code required to resolve the type of the variant is only
+executed *once* per ``Evaluate`` call. This can enable composite EOS
+calls, non-standard vector calls, and vector calls with non-standard
+loop structure.
+
+The ``Evaluate`` call has the signature
+
+.. code-block:: cpp
+
+  template<typename Functor_t>
+  PORTABLE_INLINE_FUNCTION
+  void Evaluate(Functor_t f);
+
+where a ``Functor_t`` is a class that *must* provide a ``void
+operator() const`` method templated on EOS type. ``Evaluate`` is
+decorated so that it may be evaluated on either host or device,
+depending on desired use-case. Alternatively, you may use an anonymous
+function with an `auto` argument as the input, e.g.,
+
+.. code-block::
+
+   // equivalent to [=], but with device markings
+   eos.Evaluate(PORTABLE_LAMBDA(auto eos) { /* my code snippet */ });
+
+.. warning::
+
+  It can be dangerous to use functors with side-effects. Especially
+  with GPUs it can produce very unintuitive behaviour. We recommend
+  you only make the ``operator()`` non-const if you really know what
+  you're doing. And in the anonymous function case, we recommend you
+  capture by value, not reference.
+
+To see the utlity of the ``Evaluate`` function, it's probably just
+easiest to provide an example. The following code evaluates the EOS on
+device and compares against a tabulated pressure. The total difference
+is summed using the ``Kokkos::parallel_reduce`` functionality in the
+``Kokkos`` performance portability library.
+
+.. code-block:: cpp
+
+  // The functor we use is defined here.
+  // This class definition needs to be of appropriately global scope.
+  class CheckPofRE {
+   public:
+    CheckPofRE(Real *P, Real *rho, Real *sie, int N) : P_(P), rho_(rho), sie_(sie), N_(N) {}
+    template <typename T>
+    // this is a host-only call, but if you wanted to write
+    // a function that you wanted to evaluate on device
+    // you could add the
+    // PORTABLE_INLINE_FUNCTION
+    // decorator here.
+    void operator()(const T &eos) const {
+      // Capturing member functions of a class in a lambda typically causes problems
+      // when launching a GPU kernel.
+      // Better to pull out new variables to capture before launching a kernel.
+      Real *P = P_;
+      Real *rho = rho_;
+      Real *sie = sie_;
+      // reduction target
+      Real tot_diff;
+      // reduction op
+      Kokkos::parallel_reduce(
+          "MyCheckPofRE", N_,
+          KOKKOS_LAMBDA(const int i, Real &diff) {
+            diff += std::abs(P[i] - eos.PressureFromDensityInternalEnergy(rho[i], sie[i]));
+          },
+          tot_diff);
+      std::cout << "Total difference = " << tot_diff << std::endl;
+    }
+  
+   private:
+    int N_;
+    Real *P_;
+    Real *rho_;
+    Real *sie_;
+  };
+
+  // Here we construct our functor
+  // it is assumed the pointers to device memory P, rho, sie, are defined elsewhere.
+  CheckPofRE my_op(P, rho, sie, N);
+
+  // Here we call the evaluate function
+  eos.Evaluate(my_op);
+
+  // The above two lines could have been called "in-one" with:
+  // eos.Evaluate(CheckPofRE(P, rho, sie, N));
+
+Alternatively, you could eliminate the functor and use an anonymous
+function with:
+
+.. code-block:: cpp
+
+  eos.Evaluate([=](auto eos) {
+    Real tot_diff;
+    Kokkos::parallel_reduce(
+        "MyCheckPofRE", N_,
+        KOKKOS_LAMBDA(const int i, Real &diff) {
+          diff += std::abs(P[i] - eos.PressureFromDensityInternalEnergy(rho[i], sie[i]));
+        },
+        tot_diff);
+    std::cout << "Total difference = " << tot_diff << std::endl;
+  });
+
+This is not functionality that would be available with the standard
+vector calls provided by ``singularity-eos``, at least not without
+chaining multiple parallel dispatch calls. Here we can do it in a
+single call.
+
 Lambdas and Optional Parameters
 --------------------------------
 

singularity-eos/eos/eos_base.hpp

@@ -128,6 +128,13 @@ class EosBase {
   struct is_raw_pointer
       : std::is_same<std::remove_reference_t<std::remove_cv_t<T>>, R *> {};
 
+  // Generic evaluator
+  template <typename Functor_t>
+  constexpr void Evaluate(Functor_t &f) const {
+    CRTP copy = *(static_cast<CRTP const *>(this));
+    f(copy);
+  }
+
   // Vector member functions
   template <typename RealIndexer, typename ConstRealIndexer, typename LambdaIndexer>
   inline void

singularity-eos/eos/eos_variant.hpp

@@ -72,6 +72,11 @@ class Variant {
   }
 
   // Place member functions here
+  template <typename Functor_t>
+  constexpr void Evaluate(Functor_t &f) const {
+    return mpark::visit([&f](const auto &eos) { return eos.Evaluate(f); }, eos_);
+  }
+
   PORTABLE_INLINE_FUNCTION
   Real TemperatureFromDensityInternalEnergy(const Real rho, const Real sie,
                                             Real *lambda = nullptr) const {

test/test_eos_unit.cpp

@@ -509,6 +509,74 @@ SCENARIO("Ideal gas entropy", "[IdealGas][Entropy]") {
   }
 }
 
+// A non-standard pattern where we do a reduction
+class CheckPofRE {
+ public:
+  CheckPofRE(Real *P, Real *rho, Real *sie, int N) : P_(P), rho_(rho), sie_(sie), N_(N) {}
+  template <typename T>
+  void operator()(const T &eos) {
+    Real *P = P_;
+    Real *rho = rho_;
+    Real *sie = sie_;
+    portableReduce(
+        "MyCheckPofRE", 0, N_,
+        PORTABLE_LAMBDA(const int i, int &nw) {
+          nw += !(isClose(P[i],
+                          eos.PressureFromDensityInternalEnergy(rho[i], sie[i], nullptr),
+                          1e-15));
+        },
+        nwrong);
+  }
+
+  // must be initialized to zero because portableReduce is a simple for loop on host
+  int nwrong = 0;
+
+ private:
+  int N_;
+  Real *P_;
+  Real *rho_;
+  Real *sie_;
+};
+SCENARIO("Ideal gas vector Evaluate call", "[IdealGas][Evaluate]") {
+  GIVEN("An ideal gas, and some device memory") {
+    constexpr Real Cv = 2.0;
+    constexpr Real gm1 = 0.5;
+    constexpr int N = 100;
+    constexpr Real rhomin = 1;
+    constexpr Real rhomax = 11;
+    constexpr Real drho = (rhomax - rhomin) / (N - 1);
+    constexpr Real siemin = 1;
+    constexpr Real siemax = 11;
+    constexpr Real dsie = (siemax - siemin) / (N - 1);
+
+    EOS eos_host = IdealGas(gm1, Cv);
+    auto eos_device = eos_host.GetOnDevice();
+
+    Real *P = (Real *)PORTABLE_MALLOC(N * sizeof(Real));
+    Real *rho = (Real *)PORTABLE_MALLOC(N * sizeof(Real));
+    Real *sie = (Real *)PORTABLE_MALLOC(N * sizeof(Real));
+
+    WHEN("We set P, rho, sie via the scalar API") {
+      portableFor(
+          "Set rho and sie", 0, N, PORTABLE_LAMBDA(const int i) {
+            rho[i] = rhomin + drho * i; // just some diagonal in the 2d rho/sie space
+            sie[i] = siemin + dsie * i;
+            P[i] = eos_device.PressureFromDensityInternalEnergy(rho[i], sie[i]);
+          });
+      THEN("The vector Evaluate API can be used to compare") {
+        CheckPofRE my_op(P, rho, sie, N);
+        eos_device.Evaluate(my_op);
+        REQUIRE(my_op.nwrong == 0);
+      }
+    }
+
+    eos_device.Finalize();
+    PORTABLE_FREE(P);
+    PORTABLE_FREE(rho);
+    PORTABLE_FREE(sie);
+  }
+}
+
 SCENARIO("Vector EOS", "[VectorEOS][IdealGas]") {
 
   GIVEN("Parameters for an ideal gas") {