diff --git a/CHANGELOG.md b/CHANGELOG.md index d555491177..8aef7c8d8d 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/doc/sphinx/src/using-eos.rst b/doc/sphinx/src/using-eos.rst index bbb3e413a7..402f340e96 100644 --- a/doc/sphinx/src/using-eos.rst +++ b/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 + 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 + // 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 -------------------------------- diff --git a/singularity-eos/eos/eos_base.hpp b/singularity-eos/eos/eos_base.hpp index f206a01a58..9588a57fa9 100644 --- a/singularity-eos/eos/eos_base.hpp +++ b/singularity-eos/eos/eos_base.hpp @@ -128,6 +128,13 @@ class EosBase { struct is_raw_pointer : std::is_same>, R *> {}; + // Generic evaluator + template + constexpr void Evaluate(Functor_t &f) const { + CRTP copy = *(static_cast(this)); + f(copy); + } + // Vector member functions template inline void diff --git a/singularity-eos/eos/eos_variant.hpp b/singularity-eos/eos/eos_variant.hpp index 04dc51a9bf..3e72d87662 100644 --- a/singularity-eos/eos/eos_variant.hpp +++ b/singularity-eos/eos/eos_variant.hpp @@ -72,6 +72,11 @@ class Variant { } // Place member functions here + template + 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 { diff --git a/test/test_eos_unit.cpp b/test/test_eos_unit.cpp index b595d4e16f..04f5ff009c 100644 --- a/test/test_eos_unit.cpp +++ b/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 + 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") {