Initialization on GPU

CMakeLists.txt

@@ -106,6 +106,10 @@ endif()
 set(BUILD_DOCUMENTATION OFF)
 
 ## Use Kokkos from `vendor/`
+if("${Kokkos_ENABLE_CUDA}")
+  option(Kokkos_ENABLE_CUDA_CONSTEXPR "Whether to activate experimental relaxed constexpr functions" ON)
+  option(Kokkos_ENABLE_CUDA_RELOCATABLE_DEVICE_CODE "Whether to enable relocatable device code (RDC) for CUDA" ON)
+endif()
 add_subdirectory("vendor/kokkos/" "kokkos")
 
 find_package (Eigen3 3.3 NO_MODULE)

simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp

@@ -157,30 +157,31 @@ int main(int argc, char** argv)
             std::move(masses),
             std::move(init_perturb_amplitude),
             std::move(init_perturb_mode));
-    DFieldSpVx allfequilibrium(meshSpVx);
+    device_t<DFieldSpVx> allfequilibrium_device(meshSpVx);
     BumpontailEquilibrium const init_fequilibrium(
             std::move(epsilon_bot),
             std::move(temperature_bot),
             std::move(mean_velocity_bot));
-    init_fequilibrium(allfequilibrium);
-    DFieldSpXVx allfdistribu(meshSpXVx);
+    init_fequilibrium(allfequilibrium_device);
 
     PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
     PDI_init(conf_pdi);
-
     ddc::expose_to_pdi("iter_start", iter_start);
 
+    device_t<DFieldSpXVx> allfdistribu_device(meshSpXVx);
     double time_start(0);
     if (iter_start == 0) {
         SingleModePerturbInitialization const
-                init(allfequilibrium,
+                init(allfequilibrium_device,
                      ddc::discrete_space<IDimSp>().perturb_modes(),
                      ddc::discrete_space<IDimSp>().perturb_amplitudes());
-        init(allfdistribu);
+        init(allfdistribu_device);
     } else {
         RestartInitialization const restart(iter_start, time_start);
-        restart(allfdistribu);
+        restart(allfdistribu_device);
     }
+    auto allfequilibrium = ddc::create_mirror_view_and_copy(allfequilibrium_device.span_view());
+    auto allfdistribu = ddc::create_mirror_view_and_copy(allfdistribu_device.span_view());
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -247,7 +248,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, time_start, deltat, nbiter);
+    predcorr(allfdistribu_device, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp

@@ -149,30 +149,31 @@ int main(int argc, char** argv)
             std::move(masses),
             std::move(init_perturb_amplitude),
             std::move(init_perturb_mode));
-    DFieldSpVx allfequilibrium(meshSpVx);
+    device_t<DFieldSpVx> allfequilibrium_device(meshSpVx);
     BumpontailEquilibrium const init_fequilibrium(
             std::move(epsilon_bot),
             std::move(temperature_bot),
             std::move(mean_velocity_bot));
-    init_fequilibrium(allfequilibrium);
-    DFieldSpXVx allfdistribu(meshSpXVx);
+    init_fequilibrium(allfequilibrium_device);
 
     PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
     PDI_init(conf_pdi);
-
     ddc::expose_to_pdi("iter_start", iter_start);
 
+    device_t<DFieldSpXVx> allfdistribu_device(meshSpXVx);
     double time_start(0);
     if (iter_start == 0) {
         SingleModePerturbInitialization const
-                init(allfequilibrium,
+                init(allfequilibrium_device,
                      ddc::discrete_space<IDimSp>().perturb_modes(),
                      ddc::discrete_space<IDimSp>().perturb_amplitudes());
-        init(allfdistribu);
+        init(allfdistribu_device);
     } else {
         RestartInitialization const restart(iter_start, time_start);
-        restart(allfdistribu);
+        restart(allfdistribu_device);
     }
+    auto allfequilibrium = ddc::create_mirror_view_and_copy(allfequilibrium_device.span_view());
+    auto allfdistribu = ddc::create_mirror_view_and_copy(allfdistribu_device.span_view());
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -236,7 +237,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, time_start, deltat, nbiter);
+    predcorr(allfdistribu_device, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

simulations/geometryXVx/landau/landau_fem_uniform.cpp

@@ -155,30 +155,31 @@ int main(int argc, char** argv)
             std::move(masses),
             std::move(init_perturb_amplitude),
             std::move(init_perturb_mode));
-    DFieldSpVx allfequilibrium(meshSpVx);
+    device_t<DFieldSpVx> allfequilibrium_device(meshSpVx);
     MaxwellianEquilibrium const init_fequilibrium(
             std::move(density_eq),
             std::move(temperature_eq),
             std::move(mean_velocity_eq));
-    init_fequilibrium(allfequilibrium);
-    DFieldSpXVx allfdistribu(meshSpXVx);
+    init_fequilibrium(allfequilibrium_device);
 
     PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
     PDI_init(conf_pdi);
-
     ddc::expose_to_pdi("iter_start", iter_start);
 
+    device_t<DFieldSpXVx> allfdistribu_device(meshSpXVx);
     double time_start(0);
     if (iter_start == 0) {
         SingleModePerturbInitialization const
-                init(allfequilibrium,
+                init(allfequilibrium_device,
                      ddc::discrete_space<IDimSp>().perturb_modes(),
                      ddc::discrete_space<IDimSp>().perturb_amplitudes());
-        init(allfdistribu);
+        init(allfdistribu_device);
     } else {
         RestartInitialization const restart(iter_start, time_start);
-        restart(allfdistribu);
+        restart(allfdistribu_device);
     }
+    auto allfequilibrium = ddc::create_mirror_view_and_copy(allfequilibrium_device.span_view());
+    auto allfdistribu = ddc::create_mirror_view_and_copy(allfdistribu_device.span_view());
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -246,7 +247,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, time_start, deltat, nbiter);
+    predcorr(allfdistribu_device, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

simulations/geometryXVx/landau/landau_fft.cpp

@@ -152,30 +152,31 @@ int main(int argc, char** argv)
             std::move(masses),
             std::move(init_perturb_amplitude),
             std::move(init_perturb_mode));
-    DFieldSpVx allfequilibrium(meshSpVx);
+    device_t<DFieldSpVx> allfequilibrium_device(meshSpVx);
     MaxwellianEquilibrium const init_fequilibrium(
             std::move(density_eq),
             std::move(temperature_eq),
             std::move(mean_velocity_eq));
-    init_fequilibrium(allfequilibrium);
-    DFieldSpXVx allfdistribu(meshSpXVx);
+    init_fequilibrium(allfequilibrium_device);
 
     PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
     PDI_init(conf_pdi);
-
     ddc::expose_to_pdi("iter_start", iter_start);
 
+    device_t<DFieldSpXVx> allfdistribu_device(meshSpXVx);
     double time_start(0);
     if (iter_start == 0) {
         SingleModePerturbInitialization const
-                init(allfequilibrium,
+                init(allfequilibrium_device,
                      ddc::discrete_space<IDimSp>().perturb_modes(),
                      ddc::discrete_space<IDimSp>().perturb_amplitudes());
-        init(allfdistribu);
+        init(allfdistribu_device);
     } else {
         RestartInitialization const restart(iter_start, time_start);
-        restart(allfdistribu);
+        restart(allfdistribu_device);
     }
+    auto allfequilibrium = ddc::create_mirror_view_and_copy(allfequilibrium_device.span_view());
+    auto allfdistribu = ddc::create_mirror_view_and_copy(allfdistribu_device.span_view());
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -239,7 +240,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, time_start, deltat, nbiter);
+    predcorr(allfdistribu_device, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

simulations/geometryXVx/sheath/sheath.cpp

@@ -157,30 +157,31 @@ int main(int argc, char** argv)
             std::move(masses),
             std::move(init_perturb_amplitude),
             std::move(init_perturb_mode));
-    DFieldSpVx allfequilibrium(meshSpVx);
+    device_t<DFieldSpVx> allfequilibrium_device(meshSpVx);
     MaxwellianEquilibrium const init_fequilibrium(
             std::move(density_eq),
             std::move(temperature_eq),
             std::move(mean_velocity_eq));
-    init_fequilibrium(allfequilibrium);
-    DFieldSpXVx allfdistribu(meshSpXVx);
+    init_fequilibrium(allfequilibrium_device);
 
     PC_tree_t conf_pdi = PC_parse_string(PDI_CFG);
     PDI_init(conf_pdi);
-
     ddc::expose_to_pdi("iter_start", iter_start);
 
+    device_t<DFieldSpXVx> allfdistribu_device(meshSpXVx);
     double time_start(0);
     if (iter_start == 0) {
         SingleModePerturbInitialization const
-                init(allfequilibrium,
+                init(allfequilibrium_device,
                      ddc::discrete_space<IDimSp>().perturb_modes(),
                      ddc::discrete_space<IDimSp>().perturb_amplitudes());
-        init(allfdistribu);
+        init(allfdistribu_device);
     } else {
         RestartInitialization const restart(iter_start, time_start);
-        restart(allfdistribu);
+        restart(allfdistribu_device);
     }
+    auto allfequilibrium = ddc::create_mirror_view_and_copy(allfequilibrium_device.span_view());
+    auto allfdistribu = ddc::create_mirror_view_and_copy(allfdistribu_device.span_view());
 
     // --> Algorithm info
     double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat");
@@ -322,7 +323,7 @@ int main(int argc, char** argv)
 
     steady_clock::time_point const start = steady_clock::now();
 
-    predcorr(allfdistribu, time_start, deltat, nbiter);
+    predcorr(allfdistribu_device, time_start, deltat, nbiter);
 
     steady_clock::time_point const end = steady_clock::now();
 

src/geometryXVx/geometry/CMakeLists.txt

@@ -13,6 +13,7 @@ target_link_libraries("geometry_${GEOMETRY_VARIANT}" INTERFACE
     DDC::DDC
     sll::splines
     gslx::speciesinfo
+    gslx::utils
 )
 add_library("gslx::geometry_${GEOMETRY_VARIANT}" ALIAS "geometry_${GEOMETRY_VARIANT}")
 

src/geometryXVx/geometry/geometry.hpp

@@ -3,14 +3,15 @@
 #pragma once
 
 #include <ddc/ddc.hpp>
-#include <ddc/kernels/fft.hpp> // Strange but needed to get access to Fourier tag, to be clarified
+#include <ddc/kernels/fft.hpp>
 
 #include <sll/bsplines_non_uniform.hpp>
 #include <sll/bsplines_uniform.hpp>
 #include <sll/greville_interpolation_points.hpp>
 #include <sll/spline_boundary_conditions.hpp>
 #include <sll/spline_builder.hpp>
 
+#include <ddc_helper.hpp>
 #include <species_info.hpp>
 
 /**

src/geometryXVx/initialization/bumpontailequilibrium.cpp

@@ -4,8 +4,6 @@
 
 #include "bumpontailequilibrium.hpp"
 
-using std::sqrt, std::exp;
-
 BumpontailEquilibrium::BumpontailEquilibrium(
         DViewSp const epsilon_bot,
         DViewSp const temperature_bot,
@@ -16,21 +14,26 @@ BumpontailEquilibrium::BumpontailEquilibrium(
 {
 }
 
-DSpanSpVx BumpontailEquilibrium::operator()(DSpanSpVx const allfequilibrium) const
+device_t<DSpanSpVx> BumpontailEquilibrium::operator()(
+        device_t<DSpanSpVx> const allfequilibrium) const
 {
     IDomainVx const gridvx = allfequilibrium.domain<IDimVx>();
     IDomainSp const gridsp = allfequilibrium.domain<IDimSp>();
 
     // Initialization of the maxwellian
-    DFieldVx maxwellian(gridvx);
+    device_t<DFieldVx> maxwellian_alloc(gridvx);
+    ddc::ChunkSpan maxwellian = maxwellian_alloc.span_view();
     ddc::for_each(gridsp, [&](IndexSp const isp) {
         compute_twomaxwellian(
                 maxwellian,
                 m_epsilon_bot(isp),
                 m_temperature_bot(isp),
                 m_mean_velocity_bot(isp));
 
-        ddc::for_each(gridvx, [&](IndexVx const iv) { allfequilibrium(isp, iv) = maxwellian(iv); });
+        ddc::for_each(
+                ddc::policies::parallel_device,
+                gridvx,
+                DDC_LAMBDA(IndexVx const ivx) { allfequilibrium(isp, ivx) = maxwellian(ivx); });
     });
     return allfequilibrium;
 }
@@ -44,23 +47,27 @@ with
   f2(v) = epsilon/sqrt(2*PI*T0)[exp(-(v-v0)**2/2*T0)]
 */
 void BumpontailEquilibrium::compute_twomaxwellian(
-        DSpanVx const fMaxwellian,
+        device_t<DSpanVx> const fMaxwellian,
         double const epsilon_bot,
         double const temperature_bot,
         double const mean_velocity_bot) const
 {
     double const inv_sqrt_2pi = 1. / sqrt(2. * M_PI);
     double const norm_f2 = inv_sqrt_2pi / sqrt(temperature_bot);
     IDomainVx const gridvx = fMaxwellian.domain();
-    for (IndexVx const iv : gridvx) {
-        CoordVx const v = ddc::coordinate(iv);
-        // bulk plasma particles
-        double const f1_v = (1. - epsilon_bot) * inv_sqrt_2pi * exp(-0.5 * v * v);
-        // beam
-        double const f2_v = epsilon_bot * norm_f2
-                            * (exp(-(v - mean_velocity_bot) * (v - mean_velocity_bot)
-                                   / (2. * temperature_bot)));
-        // fM(v) = f1(v) + f2(v)
-        fMaxwellian(iv) = f1_v + f2_v;
-    }
+    ddc::for_each(
+            ddc::policies::parallel_device,
+            gridvx,
+            DDC_LAMBDA(IndexVx const ivx) {
+                CoordVx const vx = ddc::coordinate(ivx);
+                // bulk plasma particles
+                double const f1_v = (1. - epsilon_bot) * inv_sqrt_2pi * Kokkos::exp(-0.5 * vx * vx);
+                // beam
+                double const f2_v = epsilon_bot * norm_f2
+                                    * (Kokkos::exp(
+                                            -(vx - mean_velocity_bot) * (vx - mean_velocity_bot)
+                                            / (2. * temperature_bot)));
+                // fM(v) = f1(v) + f2(v)
+                fMaxwellian(ivx) = f1_v + f2_v;
+            });
 }

src/geometryXVx/initialization/bumpontailequilibrium.hpp

@@ -19,7 +19,7 @@ class BumpontailEquilibrium : public IEquilibrium
     // mean velocity of the bump-on-tail for all kinetic species
     FieldSp<double> m_mean_velocity_bot;
 
-private:
+public:
     /*
       Computation of the Maxwellian fM which is the equilibrium part 
       of the distribution function : 
@@ -30,17 +30,16 @@ class BumpontailEquilibrium : public IEquilibrium
                                            +exp(-(v+v0)**2/2*T0)]
      */
     void compute_twomaxwellian(
-            DSpanVx fMaxwellian,
+            device_t<DSpanVx> fMaxwellian,
             double epsilon_bot,
             double temperature_bot,
             double mean_velocity_bot) const;
 
-public:
     BumpontailEquilibrium(DViewSp epsilon_bot, DViewSp temperature_bot, DViewSp mean_velocity_bot);
 
     ~BumpontailEquilibrium() override = default;
 
-    DSpanSpVx operator()(DSpanSpVx allfequilibrium) const override;
+    device_t<DSpanSpVx> operator()(device_t<DSpanSpVx> allfequilibrium) const override;
 
     ViewSp<double> epsilon_bot() const
     {

src/geometryXVx/initialization/iequilibrium.hpp

@@ -9,5 +9,5 @@ class IEquilibrium
 public:
     virtual ~IEquilibrium() = default;
 
-    virtual DSpanSpVx operator()(DSpanSpVx allfequilibrium) const = 0;
+    virtual device_t<DSpanSpVx> operator()(device_t<DSpanSpVx> allfequilibrium) const = 0;
 };

src/geometryXVx/initialization/iinitialization.hpp

@@ -9,5 +9,5 @@ class IInitialization
 public:
     virtual ~IInitialization() = default;
 
-    virtual DSpanSpXVx operator()(DSpanSpXVx allfdistribu) const = 0;
+    virtual device_t<DSpanSpXVx> operator()(device_t<DSpanSpXVx> allfdistribu) const = 0;
 };