Port the geometryRTheta interpolator to gpu

simulations/geometryRTheta/diocotron/diocotron.cpp

@@ -140,7 +140,7 @@ int main(int argc, char** argv)
 
     ddc::init_discrete_space<PolarBSplinesRTheta>(discrete_mapping_host);
 
-    IdxRangeBSRTheta const idx_range_bsplinesRTheta = get_spline_idx_range(builder_host);
+    IdxRangeBSRTheta const idx_range_bsplinesRTheta = get_spline_idx_range(builder);
 
 
     // --- Time integration method --------------------------------------------------------------------
@@ -183,7 +183,7 @@ int main(int argc, char** argv)
             p_extrapolation_rule,
             p_extrapolation_rule);
 
-    PreallocatableSplineInterpolatorRTheta interpolator(builder_host, spline_evaluator_host);
+    PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator);
 
     SplinePolarFootFinder find_feet(
             time_stepper,

simulations/geometryRTheta/vortex_merger/vortex_merger.cpp

@@ -158,7 +158,7 @@ int main(int argc, char** argv)
             p_extrapolation_rule,
             p_extrapolation_rule);
 
-    PreallocatableSplineInterpolatorRTheta interpolator(builder_host, spline_evaluator_host);
+    PreallocatableSplineInterpolatorRTheta interpolator(builder, spline_evaluator);
 
     SplinePolarFootFinder find_feet(
             time_stepper,

src/geometryRTheta/advection/bsl_advection_rp.hpp

@@ -95,7 +95,7 @@ class BslAdvectionRTheta : public IAdvectionRTheta
     /**
      * @brief Allocate a Field of the advected function.
      *
-     * @param [in, out] allfdistribu
+     * @param [in, out] allfdistribu_host
      *      A Field containing the values of the function we want to advect.
      * @param [in] advection_field_xy
      *      A DConstVectorFieldRTheta containing the values of the advection field
@@ -106,33 +106,41 @@ class BslAdvectionRTheta : public IAdvectionRTheta
      * @return A Field to allfdistribu advected on the time step given.
      */
     host_t<DFieldRTheta> operator()(
-            host_t<DFieldRTheta> allfdistribu,
+            host_t<DFieldRTheta> allfdistribu_host,
             host_t<DConstVectorFieldRTheta<X, Y>> advection_field_xy,
             double dt) const override
     {
         // Pre-allocate some memory to prevent allocation later in loop
         std::unique_ptr<IInterpolatorRTheta> const interpolator_ptr = m_interpolator.preallocate();
 
         // Initialise the feet
-        host_t<FieldMemRTheta<CoordRTheta>> feet_rp(get_idx_range(advection_field_xy));
+        host_t<FieldMemRTheta<CoordRTheta>> feet_rp_host(get_idx_range(advection_field_xy));
         ddc::for_each(get_idx_range(advection_field_xy), [&](IdxRTheta const irp) {
-            feet_rp(irp) = ddc::coordinate(irp);
+            feet_rp_host(irp) = ddc::coordinate(irp);
         });
 
         // Compute the characteristic feet at tn ----------------------------------------------------
-        m_find_feet(get_field(feet_rp), advection_field_xy, dt);
+        m_find_feet(get_field(feet_rp_host), advection_field_xy, dt);
+
+        auto feet_rp = ddc::create_mirror_view_and_copy(
+                Kokkos::DefaultExecutionSpace(),
+                get_field(feet_rp_host));
+        auto allfdistribu = ddc::
+                create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), allfdistribu_host);
 
         // Interpolate the function on the characteristic feet. -------------------------------------
-        (*interpolator_ptr)(allfdistribu, get_const_field(feet_rp));
+        (*interpolator_ptr)(get_field(allfdistribu), get_const_field(feet_rp));
+
+        ddc::parallel_deepcopy(allfdistribu_host, get_const_field(allfdistribu));
 
-        return allfdistribu;
+        return allfdistribu_host;
     }
 
 
     /**
      * @brief Allocate a Field to the advected function.
      *
-     * @param [in, out] allfdistribu
+     * @param [in, out] allfdistribu_host
      *      A Field containing the values of the function we want to advect.
      * @param [in] advection_field_rp
      *      A DConstVectorFieldRTheta containing the values of the advection field
@@ -146,13 +154,13 @@ class BslAdvectionRTheta : public IAdvectionRTheta
      * @return A Field to allfdistribu advected on the time step given.
      */
     host_t<DFieldRTheta> operator()(
-            host_t<DFieldRTheta> allfdistribu,
+            host_t<DFieldRTheta> allfdistribu_host,
             host_t<DConstVectorFieldRTheta<R, Theta>> advection_field_rp,
             CoordXY const& advection_field_xy_center,
             double dt) const override
     {
         Kokkos::Profiling::pushRegion("PolarAdvection");
-        IdxRangeRTheta grid(get_idx_range<GridR, GridTheta>(allfdistribu));
+        IdxRangeRTheta grid(get_idx_range<GridR, GridTheta>(allfdistribu_host));
 
         const int npoints_p = IdxRangeTheta(grid).size();
         IdxRangeRTheta const grid_without_Opoint(grid.remove_first(IdxStepRTheta(1, 0)));
@@ -191,16 +199,24 @@ class BslAdvectionRTheta : public IAdvectionRTheta
         std::unique_ptr<IInterpolatorRTheta> const interpolator_ptr = m_interpolator.preallocate();
 
         // Initialise the feet
-        host_t<FieldMemRTheta<CoordRTheta>> feet_rp(grid);
-        ddc::for_each(grid, [&](IdxRTheta const irp) { feet_rp(irp) = ddc::coordinate(irp); });
+        host_t<FieldMemRTheta<CoordRTheta>> feet_rp_host(grid);
+        ddc::for_each(grid, [&](IdxRTheta const irp) { feet_rp_host(irp) = ddc::coordinate(irp); });
 
         // Compute the characteristic feet at tn ----------------------------------------------------
-        m_find_feet(get_field(feet_rp), get_const_field(advection_field_xy), dt);
+        m_find_feet(get_field(feet_rp_host), get_const_field(advection_field_xy), dt);
+
+        auto feet_rp = ddc::create_mirror_view_and_copy(
+                Kokkos::DefaultExecutionSpace(),
+                get_field(feet_rp_host));
+        auto allfdistribu = ddc::
+                create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), allfdistribu_host);
 
         // Interpolate the function on the characteristic feet. -------------------------------------
-        (*interpolator_ptr)(allfdistribu, get_const_field(feet_rp));
+        (*interpolator_ptr)(get_field(allfdistribu), get_const_field(feet_rp));
+
+        ddc::parallel_deepcopy(allfdistribu_host, get_const_field(allfdistribu));
         Kokkos::Profiling::popRegion();
 
-        return allfdistribu;
+        return allfdistribu_host;
     }
 };

src/geometryRTheta/interpolation/i_interpolator_2d_rp.hpp

@@ -31,9 +31,9 @@ class IInterpolatorRTheta
      *
      * @return A reference to the inout_data array containing the value of the function at the coordinates.
      */
-    virtual host_t<DFieldRTheta> operator()(
-            host_t<DFieldRTheta> inout_data,
-            host_t<ConstFieldRTheta<CoordRTheta>> coordinates) const = 0;
+    virtual DFieldRTheta operator()(
+            DFieldRTheta inout_data,
+            ConstFieldRTheta<CoordRTheta> coordinates) const = 0;
 };
 
 
@@ -68,9 +68,9 @@ class IPreallocatableInterpolatorRTheta : public IInterpolatorRTheta
      */
     virtual std::unique_ptr<IInterpolatorRTheta> preallocate() const = 0;
 
-    host_t<DFieldRTheta> operator()(
-            host_t<DFieldRTheta> const inout_data,
-            host_t<ConstFieldRTheta<CoordRTheta>> const coordinates) const override
+    DFieldRTheta operator()(
+            DFieldRTheta const inout_data,
+            ConstFieldRTheta<CoordRTheta> const coordinates) const override
     {
         return (*preallocate())(inout_data, coordinates);
     }

src/geometryRTheta/interpolation/spline_interpolator_2d_rp.hpp

@@ -17,8 +17,8 @@ class SplineInterpolatorRTheta : public IInterpolatorRTheta
 public:
     /// The type of the 2D Spline Evaluator used by this class
     using evaluator_type = ddc::SplineEvaluator2D<
-            Kokkos::DefaultHostExecutionSpace,
-            Kokkos::HostSpace,
+            Kokkos::DefaultExecutionSpace,
+            Kokkos::DefaultExecutionSpace::memory_space,
             BSplinesR,
             BSplinesTheta,
             GridR,
@@ -31,26 +31,23 @@ class SplineInterpolatorRTheta : public IInterpolatorRTheta
             GridTheta>;
 
 private:
-    SplineRThetaBuilder_host const& m_builder;
+    SplineRThetaBuilder const& m_builder;
 
     evaluator_type const& m_evaluator;
 
-    mutable host_t<DFieldMem<IdxRangeBSRTheta>> m_coefs;
+    mutable DFieldMem<IdxRangeBSRTheta> m_coefs;
 
-    using r_deriv_type = host_t<DConstField<SplineRThetaBuilder_host::batched_derivs_domain_type1>>;
-    using p_deriv_type = host_t<DConstField<SplineRThetaBuilder_host::batched_derivs_domain_type2>>;
-    using mixed_deriv_type
-            = host_t<DConstField<SplineRThetaBuilder_host::batched_derivs_domain_type>>;
+    using r_deriv_type = DConstField<SplineRThetaBuilder::batched_derivs_domain_type1>;
+    using p_deriv_type = DConstField<SplineRThetaBuilder::batched_derivs_domain_type2>;
+    using mixed_deriv_type = DConstField<SplineRThetaBuilder::batched_derivs_domain_type>;
 
 public:
     /**
      * @brief Create a spline interpolator object.
      * @param[in] builder An operator which builds spline coefficients from the values of a function at known interpolation points.
      * @param[in] evaluator An operator which evaluates the value of a spline at requested coordinates.
      */
-    SplineInterpolatorRTheta(
-            SplineRThetaBuilder_host const& builder,
-            evaluator_type const& evaluator)
+    SplineInterpolatorRTheta(SplineRThetaBuilder const& builder, evaluator_type const& evaluator)
         : m_builder(builder)
         , m_evaluator(evaluator)
         , m_coefs(get_spline_idx_range(builder))
@@ -71,30 +68,29 @@ class SplineInterpolatorRTheta : public IInterpolatorRTheta
      *
      * @return A reference to the inout_data array containing the value of the function at the coordinates.
      */
-    host_t<DFieldRTheta> operator()(
-            host_t<DFieldRTheta> const inout_data,
-            host_t<Field<CoordRTheta const, IdxRangeRTheta>> const coordinates) const override
+    DFieldRTheta operator()(
+            DFieldRTheta const inout_data,
+            ConstField<CoordRTheta, IdxRangeRTheta> const coordinates) const override
     {
 #ifndef NDEBUG
         // To ensure that the interpolator is C0, we ensure that
         // the value at (r=0,theta) is the same for all theta.
         IdxRangeR r_idx_range = get_idx_range<GridR>(inout_data);
         IdxRangeTheta theta_idx_range = get_idx_range<GridTheta>(inout_data);
         if (ddc::coordinate(r_idx_range.front()) == 0) {
-            ddc::for_each(theta_idx_range, [&](IdxTheta const ip) {
-                bool const unicity_center_point
-                        = inout_data(r_idx_range.front(), ip)
-                          == inout_data(r_idx_range.front(), theta_idx_range.front());
-                if (!unicity_center_point) {
-                    std::printf("Unicity of the value at the center point is not verified.");
-                    assert(unicity_center_point);
-                }
-            });
+            ddc::parallel_for_each(
+                    theta_idx_range,
+                    KOKKOS_LAMBDA(IdxTheta const ip) {
+                        KOKKOS_ASSERT(
+                                inout_data(r_idx_range.front(), ip)
+                                == inout_data(r_idx_range.front(), theta_idx_range.front()));
+                    });
         }
 #endif
 
         m_builder(get_field(m_coefs), get_const_field(inout_data));
         m_evaluator(get_field(inout_data), coordinates, get_const_field(m_coefs));
+
         return inout_data;
     }
 };
@@ -114,8 +110,8 @@ class PreallocatableSplineInterpolatorRTheta : public IPreallocatableInterpolato
 public:
     /// The type of the 2D Spline Evaluator used by this class
     using evaluator_type = ddc::SplineEvaluator2D<
-            Kokkos::DefaultHostExecutionSpace,
-            Kokkos::HostSpace,
+            Kokkos::DefaultExecutionSpace,
+            Kokkos::DefaultExecutionSpace::memory_space,
             BSplinesR,
             BSplinesTheta,
             GridR,
@@ -128,7 +124,7 @@ class PreallocatableSplineInterpolatorRTheta : public IPreallocatableInterpolato
             GridTheta>;
 
 private:
-    SplineRThetaBuilder_host const& m_builder;
+    SplineRThetaBuilder const& m_builder;
 
     evaluator_type const& m_evaluator;
 
@@ -139,7 +135,7 @@ class PreallocatableSplineInterpolatorRTheta : public IPreallocatableInterpolato
      * @param[in] evaluator An operator which evaluates the value of a spline at requested coordinates.
      */
     PreallocatableSplineInterpolatorRTheta(
-            SplineRThetaBuilder_host const& builder,
+            SplineRThetaBuilder const& builder,
             evaluator_type const& evaluator)
         : m_builder(builder)
         , m_evaluator(evaluator)

src/math_tools/math_tools.hpp

@@ -55,7 +55,7 @@ inline T modulo(T x, T y)
     return x - y * std::floor(double(x) / y);
 }
 
-constexpr inline double ipow(double a, std::size_t i)
+KOKKOS_INLINE_FUNCTION constexpr double ipow(double a, std::size_t i)
 {
     double r(1.0);
     for (std::size_t j(0); j < i; ++j) {
@@ -64,7 +64,7 @@ constexpr inline double ipow(double a, std::size_t i)
     return r;
 }
 
-inline double ipow(double a, int i)
+KOKKOS_INLINE_FUNCTION double ipow(double a, int i)
 {
     double r(1.0);
     if (i > 0) {