GPU Properties POC

CMakeLists_files.cmake

@@ -1372,6 +1372,7 @@ list(APPEND PUBLIC_HEADER_FILES
   opm/material/fluidmatrixinteractions/EclMaterialLawHystParams.hpp
   opm/material/fluidmatrixinteractions/EclMaterialLawInitParams.hpp
   opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp
+  opm/material/fluidmatrixinteractions/GpuEclMaterialLawManager.hpp
   opm/material/fluidmatrixinteractions/EclMaterialLawReadEffectiveParams.hpp
   opm/material/fluidmatrixinteractions/EclMaterialLawTwoPhaseTypes.hpp
   opm/material/fluidmatrixinteractions/EclMultiplexerMaterial.hpp
@@ -1487,6 +1488,7 @@ list(APPEND PUBLIC_HEADER_FILES
   opm/material/thermal/EclThermalConductionLawMultiplexer.hpp
   opm/material/thermal/EclThermalConductionLawMultiplexerParams.hpp
   opm/material/thermal/EclThermalLawManager.hpp
+  opm/material/thermal/GpuEclThermalLawManager.hpp
   opm/material/thermal/EnergyModuleType.hpp
   opm/material/thermal/FluidThermalConductionLaw.hpp
   opm/material/thermal/FluidThermalConductionLawParams.hpp

opm/material/common/Tabulated1DFunction.cpp

@@ -28,8 +28,8 @@
 
 namespace Opm {
 
-template<class Scalar>
-void Tabulated1DFunction<Scalar>::
+template<class Scalar, template<class> class Storage>
+void Tabulated1DFunction<Scalar, Storage>::
 printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const
 {
     Scalar x0 = std::min(xi0, xi1);
@@ -63,10 +63,10 @@ printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const
 }
 
 template void
-Tabulated1DFunction<double>::printCSV(double,double,
-                                      unsigned,std::ostream&) const;
+Tabulated1DFunction<double, VectorWithDefaultAllocator>::printCSV(
+    double, double, unsigned, std::ostream&) const;
 template void
-Tabulated1DFunction<float>::printCSV(float,float,
-                                     unsigned,std::ostream&) const;
+Tabulated1DFunction<float, VectorWithDefaultAllocator>::printCSV(
+    float, float, unsigned, std::ostream&) const;
 
 } // namespace Opm

opm/material/common/Tabulated1DFunction.hpp

@@ -27,14 +27,19 @@
 #ifndef OPM_TABULATED_1D_FUNCTION_HPP
 #define OPM_TABULATED_1D_FUNCTION_HPP
 
+#include <opm/common/ErrorMacros.hpp>
 #include <opm/common/OpmLog/OpmLog.hpp>
+#include <opm/common/utility/VectorWithDefaultAllocator.hpp>
+#include <opm/common/utility/gpuDecorators.hpp>
+#include <opm/common/utility/gpuistl_if_available.hpp>
 #include <opm/material/densead/Math.hpp>
 
 #include <algorithm>
 #include <cassert>
 #include <cstddef>
 #include <iosfwd>
 #include <stdexcept>
+#include <utility>
 #include <vector>
 
 namespace Opm {
@@ -43,20 +48,41 @@ struct SegmentIndex {
     std::size_t value;
 };
 
+template <class Scalar,
+          template <class> class Storage = VectorWithDefaultAllocator>
+class Tabulated1DFunction;
+
+#if HAVE_CUDA
+namespace gpuistl {
+
+template <class Scalar, template <class> class ContainerT>
+Tabulated1DFunction<Scalar, GpuView>
+make_view(Tabulated1DFunction<Scalar, ContainerT>& gpuBuffers);
+
+} // namespace gpuistl
+#endif // HAVE_CUDA
+
 /*!
  * \brief Implements a linearly interpolated scalar function that depends on one
  *        variable.
  */
-template <class Scalar>
+template <class Scalar, template <class> class Storage>
 class Tabulated1DFunction
 {
 public:
+    using ValueVector = Storage<Scalar>;
+
     /*!
      * \brief Default constructor for a piecewise linear function.
      *
      * To specfiy the acutal curve, use one of the set() methods.
      */
-    Tabulated1DFunction()
+    OPM_HOST_DEVICE Tabulated1DFunction() = default;
+
+    OPM_HOST_DEVICE Tabulated1DFunction(ValueVector xValues,
+                                         ValueVector yValues)
+        : xValues_(std::move(xValues))
+        , yValues_(std::move(yValues))
     {}
 
     /*!
@@ -213,7 +239,7 @@ class Tabulated1DFunction
     /*!
      * \brief Returns the number of sampling points.
      */
-    std::size_t numSamples() const
+    OPM_HOST_DEVICE std::size_t numSamples() const
     { return xValues_.size(); }
 
     /*!
@@ -234,10 +260,10 @@ class Tabulated1DFunction
     Scalar xAt(std::size_t i) const
     { return xValues_[i]; }
 
-    const std::vector<Scalar>& xValues() const
+    const ValueVector& xValues() const
     { return xValues_; }
 
-    const std::vector<Scalar>& yValues() const
+    const ValueVector& yValues() const
     { return yValues_; }
 
     /*!
@@ -250,7 +276,7 @@ class Tabulated1DFunction
      * \brief Return true iff the given x is in range [x1, xn].
      */
     template <class Evaluation>
-    bool applies(const Evaluation& x) const
+    OPM_HOST_DEVICE bool applies(const Evaluation& x) const
     { return xValues_[0] <= x && x <= xValues_[numSamples() - 1]; }
 
     /*!
@@ -263,22 +289,21 @@ class Tabulated1DFunction
      *                    failed assertation.
      */
     template <class Evaluation>
-    Evaluation eval(const Evaluation& x, bool extrapolate = false) const
+    OPM_HOST_DEVICE Evaluation eval(const Evaluation& x, bool extrapolate = false) const
     {
         SegmentIndex segIdx = findSegmentIndex(x, extrapolate);
         return eval(x, segIdx);
     }
 
     template <class Evaluation>
-    Evaluation eval(const Evaluation& x, SegmentIndex segIdxIn) const
+    OPM_HOST_DEVICE Evaluation eval(const Evaluation& x, SegmentIndex segIdxIn) const
     {
         std::size_t segIdx = segIdxIn.value;
         Scalar x0 = xValues_[segIdx];
         Scalar x1 = xValues_[segIdx + 1];
 
         Scalar y0 = yValues_[segIdx];
         Scalar y1 = yValues_[segIdx + 1];
-
         return y0 + (y1 - y0)*(x - x0)/(x1 - x0);
     }
 
@@ -430,31 +455,24 @@ class Tabulated1DFunction
     */
     void printCSV(Scalar xi0, Scalar xi1, unsigned k, std::ostream& os) const;
 
-    bool operator==(const Tabulated1DFunction<Scalar>& data) const {
+    bool operator==(const Tabulated1DFunction<Scalar, Storage>& data) const {
         return xValues_ == data.xValues_ &&
                yValues_ == data.yValues_;
     }
 
     template <class Evaluation>
-    SegmentIndex findSegmentIndex(const Evaluation& x, bool extrapolate = false) const
+    OPM_HOST_DEVICE SegmentIndex findSegmentIndex(const Evaluation& x, bool extrapolate = false) const
     {
-        if (!isfinite(x)) {
-            throw std::runtime_error("We can not search for extrapolation/interpolation "
-                                     "segment in an 1D table for non-finite value " +
-                                     std::to_string(getValue(x)) + " .");
-        }
 
-        if (!extrapolate && !applies(x))
-            throw std::logic_error("Trying to evaluate a tabulated function outside of its range");
+#if OPM_IS_INSIDE_HOST_FUNCTION
+        // relies on std::isfinite() which is not supported in device code
+        OPM_ERROR_IF(!isfinite(x), "Trying to evaluate a tabulated function at a non-finite value");
+#endif
 
-        // we need at least two sampling points!
-        if (numSamples() < 2) {
-            throw std::logic_error("We need at least two sampling points to "
-                                   "do interpolation/extrapolation, "
-                                   "and the table only contains " +
-                                   std::to_string(numSamples()) +
-                                   " sampling points");
-        }
+        OPM_ERROR_IF(numSamples() < 2, "Trying to evaluate a tabulated function with less than two samples");
+
+        OPM_ERROR_IF(!extrapolate && !applies(x),
+                     "Trying to evaluate a tabulated function outside of its range");
 
         if (x <= xValues_[1])
             return SegmentIndex{0};
@@ -473,6 +491,7 @@ class Tabulated1DFunction
             }
 
             if (xValues_[lowerIdx] > x || x > xValues_[lowerIdx + 1]) {
+#if OPM_IS_INSIDE_HOST_FUNCTION
                 std::string msg = "Problematic interpolation/extrapolation "
                                   "segment is found for the input value " +
                                   std::to_string(Opm::getValue(x)) +
@@ -499,6 +518,9 @@ class Tabulated1DFunction
                 msg += "\n";
                 OpmLog::debug(msg);
                 throw std::runtime_error(msg);
+#else
+                OPM_THROW(std::runtime_error, "Problematic interpolation/extrapolation segment found");
+#endif
             }
             return SegmentIndex{lowerIdx};
         }
@@ -555,14 +577,14 @@ class Tabulated1DFunction
      */
     struct ComparatorX_
     {
-        explicit ComparatorX_(const std::vector<Scalar>& x)
+        explicit ComparatorX_(const ValueVector& x)
             : x_(x)
         {}
 
         bool operator ()(std::size_t idxA, std::size_t idxB) const
         { return x_.at(idxA) < x_.at(idxB); }
 
-        const std::vector<Scalar>& x_;
+        const ValueVector& x_;
     };
 
     /*!
@@ -583,7 +605,7 @@ class Tabulated1DFunction
         std::ranges::sort(idxVector, cmp);
 
         // reorder the sample points
-        std::vector<Scalar> tmpX(n), tmpY(n);
+        ValueVector tmpX(n), tmpY(n);
         for (std::size_t i = 0; i < idxVector.size(); ++ i) {
             tmpX[i] = xValues_[idxVector[i]];
             tmpY[i] = yValues_[idxVector[i]];
@@ -615,10 +637,40 @@ class Tabulated1DFunction
         yValues_.resize(nSamples);
     }
 
-    std::vector<Scalar> xValues_;
-    std::vector<Scalar> yValues_;
+#if HAVE_CUDA
+    template <class ScalarT, template <class> class ContainerT>
+    friend Tabulated1DFunction<ScalarT, gpuistl::GpuView>
+    gpuistl::make_view(Tabulated1DFunction<ScalarT, ContainerT>& gpuBuffers);
+#endif
+
+    ValueVector xValues_;
+    ValueVector yValues_;
 };
 
 } // namespace Opm
 
+#if HAVE_CUDA
+namespace Opm::gpuistl {
+
+template <class Scalar>
+Tabulated1DFunction<Scalar, GpuBuffer>
+copy_to_gpu(const Tabulated1DFunction<Scalar>& cpu)
+{
+    return Tabulated1DFunction<Scalar, GpuBuffer>(
+        GpuBuffer<Scalar>(cpu.xValues()),
+        GpuBuffer<Scalar>(cpu.yValues()));
+}
+
+template <class Scalar, template <class> class ContainerT>
+Tabulated1DFunction<Scalar, GpuView>
+make_view(Tabulated1DFunction<Scalar, ContainerT>& gpuBuffers)
+{
+    return Tabulated1DFunction<Scalar, GpuView>(
+        make_view(gpuBuffers.xValues_),
+        make_view(gpuBuffers.yValues_));
+}
+
+} // namespace Opm::gpuistl
+#endif // HAVE_CUDA
+
 #endif