Merge branch '343-allow-the-operator-of-the-mappings-to-be-callled-on…

…-gpu' into 'main'

Resolve "Allow the operator() of the mappings to be callled on GPU."

Closes #343

See merge request gysela-developpers/gyselalibxx!677

vendor/sll/include/sll/mapping/analytical_invertible_curvilinear2d_to_cartesian.hpp

@@ -17,9 +17,10 @@ template <class X, class Y, class R, class Theta>
 class AnalyticalInvertibleCurvilinear2DToCartesian : public Curvilinear2DToCartesian<X, Y, R, Theta>
 {
 public:
-    virtual ~AnalyticalInvertibleCurvilinear2DToCartesian() {};
+    KOKKOS_FUNCTION virtual ~AnalyticalInvertibleCurvilinear2DToCartesian() {}
 
-    virtual ddc::Coordinate<X, Y> operator()(ddc::Coordinate<R, Theta> const& coord) const = 0;
+    KOKKOS_FUNCTION virtual ddc::Coordinate<X, Y> operator()(
+            ddc::Coordinate<R, Theta> const& coord) const = 0;
 
     /**
      * @brief Compute the logical coordinates from the physical coordinates.
@@ -34,5 +35,6 @@ class AnalyticalInvertibleCurvilinear2DToCartesian : public Curvilinear2DToCarte
      *
      * @see Curvilinear2DToCartesian::operator()
      */
-    virtual ddc::Coordinate<R, Theta> operator()(ddc::Coordinate<X, Y> const& coord) const = 0;
+    KOKKOS_FUNCTION virtual ddc::Coordinate<R, Theta> operator()(
+            ddc::Coordinate<X, Y> const& coord) const = 0;
 };

vendor/sll/include/sll/mapping/circular_to_cartesian.hpp

@@ -68,7 +68,7 @@ class CircularToCartesian : public AnalyticalInvertibleCurvilinear2DToCartesian<
      * @param[in] other
      * 		CircularToCartesian mapping used to instantiate the new one.
      */
-    CircularToCartesian(CircularToCartesian const& other) = default;
+    KOKKOS_FUNCTION CircularToCartesian(CircularToCartesian const& other) {}
 
     /**
      * @brief Instantiate a Curvilinear2DToCartesian from another temporary CircularToCartesian (rvalue).
@@ -100,21 +100,23 @@ class CircularToCartesian : public AnalyticalInvertibleCurvilinear2DToCartesian<
      */
     CircularToCartesian& operator=(CircularToCartesian&& x) = default;
 
-    ddc::Coordinate<X, Y> operator()(ddc::Coordinate<R, Theta> const& coord) const
+    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(
+            ddc::Coordinate<R, Theta> const& coord) const final
     {
         const double r = ddc::get<R>(coord);
         const double theta = ddc::get<Theta>(coord);
-        const double x = r * std::cos(theta);
-        const double y = r * std::sin(theta);
+        const double x = r * Kokkos::cos(theta);
+        const double y = r * Kokkos::sin(theta);
         return ddc::Coordinate<X, Y>(x, y);
     }
 
-    ddc::Coordinate<R, Theta> operator()(ddc::Coordinate<X, Y> const& coord) const
+    KOKKOS_FUNCTION ddc::Coordinate<R, Theta> operator()(
+            ddc::Coordinate<X, Y> const& coord) const final
     {
         const double x = ddc::get<X>(coord);
         const double y = ddc::get<Y>(coord);
-        const double r = std::sqrt(x * x + y * y);
-        const double theta = std::atan2(y, x);
+        const double r = Kokkos::sqrt(x * x + y * y);
+        const double theta = Kokkos::atan2(y, x);
         return ddc::Coordinate<R, Theta>(r, theta);
     }
 

vendor/sll/include/sll/mapping/curvilinear2d_to_cartesian.hpp

@@ -45,7 +45,7 @@ class Curvilinear2DToCartesian
      * @param[in] other
      * 		Curvilinear2DToCartesian mapping used to instantiate the new one.
      */
-    Curvilinear2DToCartesian(Curvilinear2DToCartesian const& other) = default;
+    KOKKOS_FUNCTION Curvilinear2DToCartesian(Curvilinear2DToCartesian const& other) {}
 
     /**
      * @brief Instantiate a Curvilinear2DToCartesian from another temporary
@@ -56,7 +56,7 @@ class Curvilinear2DToCartesian
      */
     Curvilinear2DToCartesian(Curvilinear2DToCartesian&& x) = default;
 
-    virtual ~Curvilinear2DToCartesian() = default;
+    KOKKOS_FUNCTION virtual ~Curvilinear2DToCartesian() {}
 
     /**
      * @brief Assign a Curvilinear2DToCartesian from another Curvilinear2DToCartesian (lvalue).
@@ -92,7 +92,8 @@ class Curvilinear2DToCartesian
      * @return The coordinates in the physical domain.
      *
      */
-    virtual ddc::Coordinate<X, Y> operator()(ddc::Coordinate<R, Theta> const& coord) const = 0;
+    KOKKOS_FUNCTION virtual ddc::Coordinate<X, Y> operator()(
+            ddc::Coordinate<R, Theta> const& coord) const = 0;
 
     /**
      * @brief Compute the Jacobian, the determinant of the Jacobian matrix of the mapping.

vendor/sll/include/sll/mapping/czarny_to_cartesian.hpp

@@ -90,7 +90,11 @@ class CzarnyToCartesian : public AnalyticalInvertibleCurvilinear2DToCartesian<X,
      * @param[in] other
      * 		CzarnyToCartesian mapping used to instantiate the new one.
      */
-    CzarnyToCartesian(CzarnyToCartesian const& other) = default;
+    KOKKOS_FUNCTION CzarnyToCartesian(CzarnyToCartesian const& other)
+        : m_epsilon(other.epsilon())
+        , m_e(other.e())
+    {
+    }
 
     /**
      * @brief Instantiate a CzarnyToCartesian from another temporary CzarnyToCartesian (rvalue).
@@ -129,7 +133,7 @@ class CzarnyToCartesian : public AnalyticalInvertibleCurvilinear2DToCartesian<X,
      *
      * @see CzarnyToCartesian
      */
-    double epsilon() const
+    KOKKOS_FUNCTION double epsilon() const
     {
         return m_epsilon;
     }
@@ -141,35 +145,38 @@ class CzarnyToCartesian : public AnalyticalInvertibleCurvilinear2DToCartesian<X,
      *
      * @see CzarnyToCartesian
      */
-    double e() const
+    KOKKOS_FUNCTION double e() const
     {
         return m_e;
     }
 
-    ddc::Coordinate<X, Y> operator()(ddc::Coordinate<R, Theta> const& coord) const
+    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(
+            ddc::Coordinate<R, Theta> const& coord) const final
     {
         const double r = ddc::get<R>(coord);
         const double theta = ddc::get<Theta>(coord);
-        const double tmp1 = std::sqrt(m_epsilon * (m_epsilon + 2.0 * r * std::cos(theta)) + 1.0);
+        const double tmp1
+                = Kokkos::sqrt(m_epsilon * (m_epsilon + 2.0 * r * Kokkos::cos(theta)) + 1.0);
 
         const double x = (1.0 - tmp1) / m_epsilon;
-        const double y = m_e * r * std::sin(theta)
-                         / (std::sqrt(1.0 - 0.25 * m_epsilon * m_epsilon) * (2.0 - tmp1));
+        const double y = m_e * r * Kokkos::sin(theta)
+                         / (Kokkos::sqrt(1.0 - 0.25 * m_epsilon * m_epsilon) * (2.0 - tmp1));
 
         return ddc::Coordinate<X, Y>(x, y);
     }
 
-    ddc::Coordinate<R, Theta> operator()(ddc::Coordinate<X, Y> const& coord) const
+    KOKKOS_FUNCTION ddc::Coordinate<R, Theta> operator()(
+            ddc::Coordinate<X, Y> const& coord) const final
     {
         const double x = ddc::get<X>(coord);
         const double y = ddc::get<Y>(coord);
         const double ex = 1. + m_epsilon * x;
         const double ex2 = (m_epsilon * x * x - 2. * x - m_epsilon);
         const double xi2 = 1. / (1. - m_epsilon * m_epsilon * 0.25);
-        const double xi = std::sqrt(xi2);
-        const double r = std::sqrt(y * y * ex * ex / (m_e * m_e * xi2) + ex2 * ex2 * 0.25);
+        const double xi = Kokkos::sqrt(xi2);
+        const double r = Kokkos::sqrt(y * y * ex * ex / (m_e * m_e * xi2) + ex2 * ex2 * 0.25);
         double theta
-                = std::atan2(2. * y * ex, (m_e * xi * (m_epsilon * x * x - 2. * x - m_epsilon)));
+                = Kokkos::atan2(2. * y * ex, (m_e * xi * (m_epsilon * x * x - 2. * x - m_epsilon)));
         if (theta < 0) {
             theta = 2 * M_PI + theta;
         }

vendor/sll/include/sll/mapping/discrete_mapping_builder.hpp

@@ -79,26 +79,18 @@ class DiscreteToCartesianBuilder
         , m_curvilinear_to_y_spline_alloc(builder.spline_domain())
         , m_evaluator(evaluator)
     {
-        static_assert(
-                std::is_same_v<ExecSpace, Kokkos::DefaultHostExecutionSpace>,
-                "GPU version is not yet implemented as mappings are not ported to GPU.");
-        using CurvilinearCoeff = ddc::
-                Coordinate<typename Mapping::circular_tag_r, typename Mapping::circular_tag_theta>;
-        using CartesianCoeff = ddc::Coordinate<X, Y>;
-
         SplineCoeffs curvilinear_to_x_spline = m_curvilinear_to_x_spline_alloc.span_view();
         SplineCoeffs curvilinear_to_y_spline = m_curvilinear_to_y_spline_alloc.span_view();
         InterpolationFieldMem curvilinear_to_x_vals_alloc(builder.interpolation_domain());
         InterpolationFieldMem curvilinear_to_y_vals_alloc(builder.interpolation_domain());
         InterpolationField curvilinear_to_x_vals = curvilinear_to_x_vals_alloc.span_view();
         InterpolationField curvilinear_to_y_vals = curvilinear_to_y_vals_alloc.span_view();
 
-        ddc::for_each(builder.interpolation_domain(), [&](IdxInterpolationPoints el) {
-            CurvilinearCoeff polar_coord(ddc::coordinate(el));
-            CartesianCoeff cart_coord = analytical_mapping(polar_coord);
-            curvilinear_to_x_vals(el) = ddc::select<X>(cart_coord);
-            curvilinear_to_y_vals(el) = ddc::select<Y>(cart_coord);
-        });
+        set_curvilinear_to_cartesian_values(
+                curvilinear_to_x_vals,
+                curvilinear_to_y_vals,
+                analytical_mapping,
+                builder.interpolation_domain());
 
         builder(curvilinear_to_x_spline.span_view(), curvilinear_to_x_vals.span_cview());
         builder(curvilinear_to_y_spline.span_view(), curvilinear_to_y_vals.span_cview());
@@ -119,6 +111,38 @@ class DiscreteToCartesianBuilder
                 m_evaluator,
                 m_idx_range_theta);
     }
+
+    /**
+     * @brief Fill in the curvilinear fields with interpolation 
+     * points mapped with the given analytical mapping. 
+     * Public method used for the constructor of the class. 
+     * @tparam Mapping Type of the analytical mapping. 
+     * @param[out] curvilinear_to_x_vals Field of coordinate on X. 
+     * @param[out] curvilinear_to_y_vals Field of coordinate on Y. 
+     * @param[in] analytical_mapping Analytical mapping. 
+     * @param[in] interpolation_idx_range Index range of an interpolation grid. 
+     */
+    template <class Mapping>
+    void set_curvilinear_to_cartesian_values(
+            InterpolationField const& curvilinear_to_x_vals,
+            InterpolationField const& curvilinear_to_y_vals,
+            Mapping const& analytical_mapping,
+            IdxRangeInterpolationPoints const& interpolation_idx_range)
+    {
+        using CurvilinearCoeff = ddc::
+                Coordinate<typename Mapping::circular_tag_r, typename Mapping::circular_tag_theta>;
+        using CartesianCoeff = ddc::Coordinate<X, Y>;
+
+        ddc::parallel_for_each(
+                ExecSpace(),
+                interpolation_idx_range,
+                KOKKOS_LAMBDA(IdxInterpolationPoints el) {
+                    CurvilinearCoeff polar_coord(ddc::coordinate(el));
+                    CartesianCoeff cart_coord = analytical_mapping(polar_coord);
+                    curvilinear_to_x_vals(el) = ddc::select<X>(cart_coord);
+                    curvilinear_to_y_vals(el) = ddc::select<Y>(cart_coord);
+                });
+    }
 };
 
 /**

vendor/sll/include/sll/mapping/discrete_to_cartesian.hpp

@@ -134,7 +134,7 @@ class DiscreteToCartesian
      *
      * @see SplineEvaluator2D
      */
-    ddc::Coordinate<X, Y> operator()(
+    KOKKOS_FUNCTION ddc::Coordinate<X, Y> operator()(
             ddc::Coordinate<circular_tag_r, circular_tag_theta> const& coord) const final
     {
         const double x = m_spline_evaluator(coord, m_x_spline_representation.span_cview());

vendor/sll/tests/CMakeLists.txt

@@ -117,6 +117,19 @@ target_link_libraries(refined_discrete_mapping_test
 )
 gtest_discover_tests(refined_discrete_mapping_test DISCOVERY_MODE PRE_TEST)
 
+add_executable(mapping_execution_space_access
+  main.cpp
+	mapping_execution_space_access.cpp
+)
+target_compile_features(mapping_execution_space_access PUBLIC cxx_std_17)
+target_link_libraries(mapping_execution_space_access
+    PUBLIC
+        GTest::gtest
+        sll::SLL
+)
+gtest_discover_tests(mapping_execution_space_access DISCOVERY_MODE PRE_TEST)
+
+
 foreach(CONTINUITY RANGE -1 1)
   math(EXPR MIN_DEGREE "${CONTINUITY}+1")
   if (${MIN_DEGREE} LESS 1)