From 3bd5fdbd36abdb0c99b1b6dc279d92396745a2b3 Mon Sep 17 00:00:00 2001 From: Thomas Padioleau Date: Thu, 16 Nov 2023 17:01:59 +0000 Subject: [PATCH] Initialization on GPU --- CMakeLists.txt | 4 ++ .../bump_on_tail/bumpontail_fem_uniform.cpp | 17 +++---- .../bump_on_tail/bumpontail_fft.cpp | 17 +++---- .../geometryXVx/landau/landau_fem_uniform.cpp | 17 +++---- simulations/geometryXVx/landau/landau_fft.cpp | 17 +++---- simulations/geometryXVx/sheath/sheath.cpp | 17 +++---- src/geometryXVx/geometry/CMakeLists.txt | 1 + src/geometryXVx/geometry/geometry.hpp | 3 +- .../initialization/bumpontailequilibrium.cpp | 41 +++++++++------- .../initialization/bumpontailequilibrium.hpp | 7 ++- .../initialization/iequilibrium.hpp | 2 +- .../initialization/iinitialization.hpp | 2 +- .../initialization/maxwellianequilibrium.cpp | 31 +++++++----- .../initialization/maxwellianequilibrium.hpp | 4 +- .../initialization/restartinitialization.cpp | 7 ++- .../initialization/restartinitialization.hpp | 2 +- .../singlemodeperturbinitialization.cpp | 47 +++++++++++-------- .../singlemodeperturbinitialization.hpp | 14 +++--- src/geometryXVx/rhs/collisions_inter.cpp | 6 ++- src/geometryXVx/rhs/krook_source_adaptive.cpp | 5 +- src/geometryXVx/rhs/krook_source_constant.cpp | 5 +- .../time_integration/itimesolver.hpp | 4 +- src/geometryXVx/time_integration/predcorr.cpp | 10 ++-- src/geometryXVx/time_integration/predcorr.hpp | 7 ++- src/utils/ddc_helper.hpp | 30 ++++++++++++ tests/geometryXVx/CMakeLists.txt | 3 +- tests/geometryXVx/collisions_inter.cpp | 5 +- tests/geometryXVx/collisions_intra_gridvx.cpp | 1 - .../collisions_intra_maxwellian.cpp | 10 ++-- tests/geometryXVx/fluid_moments.cpp | 6 ++- tests/geometryXVx/krooksource.cpp | 32 +++++++------ 31 files changed, 233 insertions(+), 141 deletions(-) diff --git a/CMakeLists.txt b/CMakeLists.txt index e460b11a9..82adfcd68 100644 --- a/CMakeLists.txt +++ b/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) diff --git a/simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp b/simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp index ea2dcc650..5d7759e2c 100644 --- a/simulations/geometryXVx/bump_on_tail/bumpontail_fem_uniform.cpp +++ b/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 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 allfdistribu_device(meshSpXVx); double time_start(0); if (iter_start == 0) { SingleModePerturbInitialization const - init(allfequilibrium, + init(allfequilibrium_device, ddc::discrete_space().perturb_modes(), ddc::discrete_space().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(); diff --git a/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp b/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp index a14935d08..faa74de6d 100644 --- a/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp +++ b/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 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 allfdistribu_device(meshSpXVx); double time_start(0); if (iter_start == 0) { SingleModePerturbInitialization const - init(allfequilibrium, + init(allfequilibrium_device, ddc::discrete_space().perturb_modes(), ddc::discrete_space().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(); diff --git a/simulations/geometryXVx/landau/landau_fem_uniform.cpp b/simulations/geometryXVx/landau/landau_fem_uniform.cpp index 9c7c8ce99..e2dff3232 100644 --- a/simulations/geometryXVx/landau/landau_fem_uniform.cpp +++ b/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 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 allfdistribu_device(meshSpXVx); double time_start(0); if (iter_start == 0) { SingleModePerturbInitialization const - init(allfequilibrium, + init(allfequilibrium_device, ddc::discrete_space().perturb_modes(), ddc::discrete_space().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(); diff --git a/simulations/geometryXVx/landau/landau_fft.cpp b/simulations/geometryXVx/landau/landau_fft.cpp index 7b288e30a..0573df227 100644 --- a/simulations/geometryXVx/landau/landau_fft.cpp +++ b/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 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 allfdistribu_device(meshSpXVx); double time_start(0); if (iter_start == 0) { SingleModePerturbInitialization const - init(allfequilibrium, + init(allfequilibrium_device, ddc::discrete_space().perturb_modes(), ddc::discrete_space().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(); diff --git a/simulations/geometryXVx/sheath/sheath.cpp b/simulations/geometryXVx/sheath/sheath.cpp index fe15d6255..147b5b4a7 100644 --- a/simulations/geometryXVx/sheath/sheath.cpp +++ b/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 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 allfdistribu_device(meshSpXVx); double time_start(0); if (iter_start == 0) { SingleModePerturbInitialization const - init(allfequilibrium, + init(allfequilibrium_device, ddc::discrete_space().perturb_modes(), ddc::discrete_space().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(); diff --git a/src/geometryXVx/geometry/CMakeLists.txt b/src/geometryXVx/geometry/CMakeLists.txt index ab8df4e9a..b3d38e861 100644 --- a/src/geometryXVx/geometry/CMakeLists.txt +++ b/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}") diff --git a/src/geometryXVx/geometry/geometry.hpp b/src/geometryXVx/geometry/geometry.hpp index f38929a97..eed1a8e72 100644 --- a/src/geometryXVx/geometry/geometry.hpp +++ b/src/geometryXVx/geometry/geometry.hpp @@ -3,7 +3,7 @@ #pragma once #include -#include // Strange but needed to get access to Fourier tag, to be clarified +#include #include #include @@ -11,6 +11,7 @@ #include #include +#include #include /** diff --git a/src/geometryXVx/initialization/bumpontailequilibrium.cpp b/src/geometryXVx/initialization/bumpontailequilibrium.cpp index 919316424..5161fa196 100644 --- a/src/geometryXVx/initialization/bumpontailequilibrium.cpp +++ b/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,13 +14,15 @@ BumpontailEquilibrium::BumpontailEquilibrium( { } -DSpanSpVx BumpontailEquilibrium::operator()(DSpanSpVx const allfequilibrium) const +device_t BumpontailEquilibrium::operator()( + device_t const allfequilibrium) const { IDomainVx const gridvx = allfequilibrium.domain(); IDomainSp const gridsp = allfequilibrium.domain(); // Initialization of the maxwellian - DFieldVx maxwellian(gridvx); + device_t maxwellian_alloc(gridvx); + ddc::ChunkSpan maxwellian = maxwellian_alloc.span_view(); ddc::for_each(gridsp, [&](IndexSp const isp) { compute_twomaxwellian( maxwellian, @@ -30,7 +30,10 @@ DSpanSpVx BumpontailEquilibrium::operator()(DSpanSpVx const allfequilibrium) con 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,7 +47,7 @@ with f2(v) = epsilon/sqrt(2*PI*T0)[exp(-(v-v0)**2/2*T0)] */ void BumpontailEquilibrium::compute_twomaxwellian( - DSpanVx const fMaxwellian, + device_t const fMaxwellian, double const epsilon_bot, double const temperature_bot, double const mean_velocity_bot) const @@ -52,15 +55,19 @@ void BumpontailEquilibrium::compute_twomaxwellian( 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; + }); } diff --git a/src/geometryXVx/initialization/bumpontailequilibrium.hpp b/src/geometryXVx/initialization/bumpontailequilibrium.hpp index 8e46ddd23..2595addeb 100644 --- a/src/geometryXVx/initialization/bumpontailequilibrium.hpp +++ b/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 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 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 operator()(device_t allfequilibrium) const override; ViewSp epsilon_bot() const { diff --git a/src/geometryXVx/initialization/iequilibrium.hpp b/src/geometryXVx/initialization/iequilibrium.hpp index df5e921f1..57b0d5e2c 100644 --- a/src/geometryXVx/initialization/iequilibrium.hpp +++ b/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 operator()(device_t allfequilibrium) const = 0; }; diff --git a/src/geometryXVx/initialization/iinitialization.hpp b/src/geometryXVx/initialization/iinitialization.hpp index dcba17a32..f2a5acd13 100644 --- a/src/geometryXVx/initialization/iinitialization.hpp +++ b/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 operator()(device_t allfdistribu) const = 0; }; diff --git a/src/geometryXVx/initialization/maxwellianequilibrium.cpp b/src/geometryXVx/initialization/maxwellianequilibrium.cpp index fbc0ef973..c57659e1d 100644 --- a/src/geometryXVx/initialization/maxwellianequilibrium.cpp +++ b/src/geometryXVx/initialization/maxwellianequilibrium.cpp @@ -14,13 +14,15 @@ MaxwellianEquilibrium::MaxwellianEquilibrium( { } -DSpanSpVx MaxwellianEquilibrium::operator()(DSpanSpVx const allfequilibrium) const +device_t MaxwellianEquilibrium::operator()( + device_t const allfequilibrium) const { IDomainVx const gridvx = allfequilibrium.domain(); IDomainSp const gridsp = allfequilibrium.domain(); // Initialization of the maxwellian - DFieldVx maxwellian(gridvx); + device_t maxwellian_alloc(gridvx); + ddc::ChunkSpan maxwellian = maxwellian_alloc.span_view(); ddc::for_each(gridsp, [&](IndexSp const isp) { compute_maxwellian( maxwellian, @@ -28,7 +30,10 @@ DSpanSpVx MaxwellianEquilibrium::operator()(DSpanSpVx const allfequilibrium) con m_temperature_eq(isp), m_mean_velocity_eq(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; } @@ -40,17 +45,21 @@ DSpanSpVx MaxwellianEquilibrium::operator()(DSpanSpVx const allfequilibrium) con where u is the mean velocity */ void MaxwellianEquilibrium::compute_maxwellian( - DSpanVx const fMaxwellian, + device_t const fMaxwellian, double const density, double const temperature, double const mean_velocity) { - double const inv_sqrt_2piT = 1. / std::sqrt(2. * M_PI * temperature); + double const inv_sqrt_2piT = 1. / Kokkos::sqrt(2. * M_PI * temperature); IDomainVx const gridvx = fMaxwellian.domain(); - for (IndexVx const iv : gridvx) { - CoordVx const v = ddc::coordinate(iv); - fMaxwellian(iv) - = density * inv_sqrt_2piT - * std::exp(-(v - mean_velocity) * (v - mean_velocity) / (2. * temperature)); - } + ddc::for_each( + ddc::policies::parallel_device, + gridvx, + DDC_LAMBDA(IndexVx const ivx) { + CoordVx const vx = ddc::coordinate(ivx); + fMaxwellian(ivx) = density * inv_sqrt_2piT + * Kokkos::exp( + -(vx - mean_velocity) * (vx - mean_velocity) + / (2. * temperature)); + }); } diff --git a/src/geometryXVx/initialization/maxwellianequilibrium.hpp b/src/geometryXVx/initialization/maxwellianequilibrium.hpp index b5f8661b7..9d20cda07 100644 --- a/src/geometryXVx/initialization/maxwellianequilibrium.hpp +++ b/src/geometryXVx/initialization/maxwellianequilibrium.hpp @@ -24,7 +24,7 @@ class MaxwellianEquilibrium : public IEquilibrium ~MaxwellianEquilibrium() override = default; - DSpanSpVx operator()(DSpanSpVx allfequilibrium) const override; + device_t operator()(device_t allfequilibrium) const override; /** * Computing the non-centered Maxwellian function as @@ -33,7 +33,7 @@ class MaxwellianEquilibrium : public IEquilibrium * where u is the mean velocity */ static void compute_maxwellian( - DSpanVx const fMaxwellian, + device_t const fMaxwellian, double const density, double const temperature, double const mean_velocity); diff --git a/src/geometryXVx/initialization/restartinitialization.cpp b/src/geometryXVx/initialization/restartinitialization.cpp index bc9b2f734..e556acc82 100644 --- a/src/geometryXVx/initialization/restartinitialization.cpp +++ b/src/geometryXVx/initialization/restartinitialization.cpp @@ -10,8 +10,11 @@ RestartInitialization::RestartInitialization(int iter_start, double& time_start) { } -DSpanSpXVx RestartInitialization::operator()(DSpanSpXVx const allfdistribu) const +device_t RestartInitialization::operator()( + device_t const allfdistribu_device) const { + auto allfdistribu = ddc::create_mirror_view_and_copy(allfdistribu_device.span_view()); ddc::PdiEvent("restart").with("time_saved", m_time_start).with("fdistribu", allfdistribu); - return allfdistribu; + ddc::deepcopy(allfdistribu_device, allfdistribu); + return allfdistribu_device; } \ No newline at end of file diff --git a/src/geometryXVx/initialization/restartinitialization.hpp b/src/geometryXVx/initialization/restartinitialization.hpp index 875cac2c0..1d98e5b96 100644 --- a/src/geometryXVx/initialization/restartinitialization.hpp +++ b/src/geometryXVx/initialization/restartinitialization.hpp @@ -41,5 +41,5 @@ class RestartInitialization : public IInitialization * read from an external file. * @return The initialized distribution function. */ - DSpanSpXVx operator()(DSpanSpXVx allfdistribu) const override; + device_t operator()(device_t allfdistribu_device) const override; }; \ No newline at end of file diff --git a/src/geometryXVx/initialization/singlemodeperturbinitialization.cpp b/src/geometryXVx/initialization/singlemodeperturbinitialization.cpp index d4ac0c5e5..49b577561 100644 --- a/src/geometryXVx/initialization/singlemodeperturbinitialization.cpp +++ b/src/geometryXVx/initialization/singlemodeperturbinitialization.cpp @@ -6,7 +6,7 @@ #include "singlemodeperturbinitialization.hpp" SingleModePerturbInitialization::SingleModePerturbInitialization( - DViewSpVx fequilibrium, + device_t fequilibrium, ViewSp const init_perturb_mode, DViewSp const init_perturb_amplitude) : m_fequilibrium(fequilibrium) @@ -15,36 +15,40 @@ SingleModePerturbInitialization::SingleModePerturbInitialization( { } -DSpanSpXVx SingleModePerturbInitialization::operator()(DSpanSpXVx const allfdistribu) const +device_t SingleModePerturbInitialization::operator()( + device_t const allfdistribu) const { - IDomainX const gridx = allfdistribu.domain(); - IDomainVx const gridvx = allfdistribu.domain(); IDomainSp const gridsp = allfdistribu.domain(); + IDomainX const gridx = allfdistribu.domain(); + IDomainXVx const gridxvx = allfdistribu.domain(); // Initialization of the perturbation - DFieldX perturbation(gridx); + device_t perturbation_alloc(gridx); + ddc::ChunkSpan perturbation = perturbation_alloc.span_view(); ddc::for_each(gridsp, [&](IndexSp const isp) { perturbation_initialization( perturbation, m_init_perturb_mode(isp), m_init_perturb_amplitude(isp)); - // Initialization of the distribution function --> fill values - ddc::for_each(gridx, [&](IndexX const ix) { - ddc::for_each(gridvx, [&](IndexVx const iv) { - double fdistribu_val = m_fequilibrium(isp, iv) * (1. + perturbation(ix)); - if (fdistribu_val < 1.e-60) { - fdistribu_val = 1.e-60; - } - allfdistribu(isp, ix, iv) = fdistribu_val; - }); - }); + ddc::for_each( + ddc::policies::parallel_device, + gridxvx, + KOKKOS_CLASS_LAMBDA(IndexXVx const ixvx) { + IndexX const ix = ddc::select(ixvx); + IndexVx const ivx = ddc::select(ixvx); + double fdistribu_val = m_fequilibrium(isp, ivx) * (1. + perturbation(ix)); + if (fdistribu_val < 1.e-60) { + fdistribu_val = 1.e-60; + } + allfdistribu(isp, ix, ivx) = fdistribu_val; + }); }); return allfdistribu; } void SingleModePerturbInitialization::perturbation_initialization( - DSpanX const perturbation, + device_t const perturbation, int const mode, double const perturb_amplitude) const { @@ -52,8 +56,11 @@ void SingleModePerturbInitialization::perturbation_initialization( double const Lx = ddcHelper::total_interval_length(gridx); double const kx = mode * 2. * M_PI / Lx; - for (IndexX const ix : gridx) { - CoordX const x = ddc::coordinate(ix); - perturbation(ix) = perturb_amplitude * cos(kx * x); - } + ddc::for_each( + ddc::policies::parallel_device, + gridx, + DDC_LAMBDA(IndexX const ix) { + CoordX const x = ddc::coordinate(ix); + perturbation(ix) = perturb_amplitude * Kokkos::cos(kx * x); + }); } diff --git a/src/geometryXVx/initialization/singlemodeperturbinitialization.hpp b/src/geometryXVx/initialization/singlemodeperturbinitialization.hpp index 6f63ed642..387331e8b 100644 --- a/src/geometryXVx/initialization/singlemodeperturbinitialization.hpp +++ b/src/geometryXVx/initialization/singlemodeperturbinitialization.hpp @@ -10,26 +10,28 @@ /// Initialization operator with a sinusoidal perturbation of a Maxwellian. This initializes all species. class SingleModePerturbInitialization : public IInitialization { - DViewSpVx m_fequilibrium; + device_t m_fequilibrium; ViewSp m_init_perturb_mode; DViewSp m_init_perturb_amplitude; -private: +public: /* Initialization of the perturbation */ - void perturbation_initialization(DSpanX perturbation, int mode, double perturb_amplitude) const; + void perturbation_initialization( + device_t perturbation, + int mode, + double perturb_amplitude) const; -public: // Warning: all variables shall remain valid until the last call to `operator()` SingleModePerturbInitialization( - DViewSpVx fequilibrium, + device_t fequilibrium, ViewSp init_perturb_mode, DViewSp init_perturb_amplitude); ~SingleModePerturbInitialization() override = default; - DSpanSpXVx operator()(DSpanSpXVx allfdistribu) const override; + device_t operator()(device_t allfdistribu) const override; }; diff --git a/src/geometryXVx/rhs/collisions_inter.cpp b/src/geometryXVx/rhs/collisions_inter.cpp index 7cb64926f..67d5c68a8 100644 --- a/src/geometryXVx/rhs/collisions_inter.cpp +++ b/src/geometryXVx/rhs/collisions_inter.cpp @@ -82,12 +82,14 @@ void CollisionsInter::compute_rhs(DSpanSpXVx const rhs, DViewSpXVx const allfdis temperature); ddc::for_each(dom_sp, [&](IndexSp const isp) { - DFieldVx fmaxwellian(gridvx); + device_t fmaxwellian_device(gridvx); MaxwellianEquilibrium::compute_maxwellian( - fmaxwellian.span_view(), + fmaxwellian_device.span_view(), density(isp), temperature(isp), fluid_velocity(isp)); + auto fmaxwellian = ddc::create_mirror_view_and_copy(fmaxwellian_device.span_view()); + ddc::for_each(gridvx, [&](IndexVx const ivx) { double const coordv = ddc::coordinate(ivx); diff --git a/src/geometryXVx/rhs/krook_source_adaptive.cpp b/src/geometryXVx/rhs/krook_source_adaptive.cpp index 9a2a89256..03cb59014 100644 --- a/src/geometryXVx/rhs/krook_source_adaptive.cpp +++ b/src/geometryXVx/rhs/krook_source_adaptive.cpp @@ -58,7 +58,10 @@ KrookSourceAdaptive::KrookSourceAdaptive( } // target distribution function - MaxwellianEquilibrium::compute_maxwellian(m_ftarget, m_density, m_temperature, 0.); + device_t ftarget_device(gridvx); + MaxwellianEquilibrium::compute_maxwellian(ftarget_device, m_density, m_temperature, 0.); + auto ftarget = ddc::create_mirror_view_and_copy(ftarget_device.span_view()); + ddc::deepcopy(m_ftarget, ftarget); switch (m_type) { case RhsType::Source: diff --git a/src/geometryXVx/rhs/krook_source_constant.cpp b/src/geometryXVx/rhs/krook_source_constant.cpp index 162a3f23f..3cb63e016 100644 --- a/src/geometryXVx/rhs/krook_source_constant.cpp +++ b/src/geometryXVx/rhs/krook_source_constant.cpp @@ -51,7 +51,10 @@ KrookSourceConstant::KrookSourceConstant( break; } // target distribution function - MaxwellianEquilibrium::compute_maxwellian(m_ftarget, m_density, m_temperature, 0.); + device_t ftarget_device(gridvx); + MaxwellianEquilibrium::compute_maxwellian(ftarget_device, m_density, m_temperature, 0.); + auto ftarget = ddc::create_mirror_view_and_copy(ftarget_device.span_view()); + ddc::deepcopy(m_ftarget, ftarget); switch (m_type) { case RhsType::Source: diff --git a/src/geometryXVx/time_integration/itimesolver.hpp b/src/geometryXVx/time_integration/itimesolver.hpp index f9d9fd2c3..2ddee166a 100644 --- a/src/geometryXVx/time_integration/itimesolver.hpp +++ b/src/geometryXVx/time_integration/itimesolver.hpp @@ -22,8 +22,8 @@ class ITimeSolver * @param[in] steps The number of iterations to be performed by the solver. * @return The distribution function after solving the system. */ - virtual DSpanSpXVx operator()( - DSpanSpXVx allfdistribu, + virtual device_t operator()( + device_t allfdistribu, double time_start, double dt, int steps = 1) const = 0; diff --git a/src/geometryXVx/time_integration/predcorr.cpp b/src/geometryXVx/time_integration/predcorr.cpp index 09440ae7b..bf8f755f1 100644 --- a/src/geometryXVx/time_integration/predcorr.cpp +++ b/src/geometryXVx/time_integration/predcorr.cpp @@ -16,12 +16,15 @@ PredCorr::PredCorr(IBoltzmannSolver const& boltzmann_solver, IPoissonSolver cons { } -DSpanSpXVx PredCorr::operator()( - DSpanSpXVx const allfdistribu, +device_t PredCorr::operator()( + device_t const allfdistribu_device, double const time_start, double const dt, int const steps) const { + auto allfdistribu_alloc = ddc::create_mirror_view_and_copy(allfdistribu_device); + ddc::ChunkSpan allfdistribu = allfdistribu_alloc.span_view(); + // electrostatic potential and electric field (depending only on x) DFieldX electrostatic_potential(allfdistribu.domain()); DFieldX electric_field(allfdistribu.domain()); @@ -66,5 +69,6 @@ DSpanSpXVx PredCorr::operator()( .and_with("fdistribu", allfdistribu) .and_with("electrostatic_potential", electrostatic_potential); - return allfdistribu; + ddc::deepcopy(allfdistribu_device, allfdistribu); + return allfdistribu_device; } diff --git a/src/geometryXVx/time_integration/predcorr.hpp b/src/geometryXVx/time_integration/predcorr.hpp index 057ceda8b..4ec4536a1 100644 --- a/src/geometryXVx/time_integration/predcorr.hpp +++ b/src/geometryXVx/time_integration/predcorr.hpp @@ -46,6 +46,9 @@ class PredCorr : public ITimeSolver * @param[in] steps The number of iterations to be performed by the predictor-corrector. * @return The distribution function after solving the system. */ - DSpanSpXVx operator()(DSpanSpXVx allfdistribu, double time_start, double dt, int steps = 1) - const override; + device_t operator()( + device_t allfdistribu, + double time_start, + double dt, + int steps = 1) const override; }; diff --git a/src/utils/ddc_helper.hpp b/src/utils/ddc_helper.hpp index c4ef7097f..0b9313e77 100644 --- a/src/utils/ddc_helper.hpp +++ b/src/utils/ddc_helper.hpp @@ -39,3 +39,33 @@ constexpr std::enable_if_t total_interval_length( return std::fabs(ddc::rlength(dom) + ddc::step>()); } }; // namespace ddcHelper + +// Template type giving the "device" version of a Chunk/Chunkspan/View +template +struct Device +{ +}; + +template +struct Device> +{ + using type = typename ddc::Chunk< + ElementType, + SupportType, + ddc::KokkosAllocator>; +}; + +template +struct Device> +{ + using type = typename ddc::ChunkSpan< + ElementType, + SupportType, + Layout, + Kokkos::DefaultExecutionSpace::memory_space>; +}; + +template +using device_t = typename Device::type; + +// namespace ddcHelper diff --git a/tests/geometryXVx/CMakeLists.txt b/tests/geometryXVx/CMakeLists.txt index df5bd59e6..f16c119e2 100644 --- a/tests/geometryXVx/CMakeLists.txt +++ b/tests/geometryXVx/CMakeLists.txt @@ -13,11 +13,12 @@ add_executable(unit_tests_${GEOMETRY_VARIANT} fluid_moments.cpp kineticsource.cpp krooksource.cpp - ../main.cpp masks.cpp quadrature.cpp splitvlasovsolver.cpp + ../main.cpp ) + target_compile_features(unit_tests_${GEOMETRY_VARIANT} PUBLIC cxx_std_17) target_link_libraries(unit_tests_${GEOMETRY_VARIANT} PUBLIC diff --git a/tests/geometryXVx/collisions_inter.cpp b/tests/geometryXVx/collisions_inter.cpp index 309685540..1a7f6d146 100644 --- a/tests/geometryXVx/collisions_inter.cpp +++ b/tests/geometryXVx/collisions_inter.cpp @@ -86,12 +86,13 @@ TEST(CollisionsInter, CollisionsInter) temperature_init(my_ielec) = 1.2; double const fluid_velocity_init(0.); ddc::for_each(ddc::get_domain(allfdistribu), [&](IndexSpX const ispx) { - DFieldVx finit(gridvx); + device_t finit_device(gridvx); MaxwellianEquilibrium::compute_maxwellian( - finit.span_view(), + finit_device.span_view(), density_init, temperature_init(ddc::select(ispx)), fluid_velocity_init); + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); ddc::deepcopy(allfdistribu[ispx], finit); }); diff --git a/tests/geometryXVx/collisions_intra_gridvx.cpp b/tests/geometryXVx/collisions_intra_gridvx.cpp index 93d4b3159..8e61717a1 100644 --- a/tests/geometryXVx/collisions_intra_gridvx.cpp +++ b/tests/geometryXVx/collisions_intra_gridvx.cpp @@ -9,7 +9,6 @@ #include #include #include -#include #include #include #include diff --git a/tests/geometryXVx/collisions_intra_maxwellian.cpp b/tests/geometryXVx/collisions_intra_maxwellian.cpp index 337aae87a..3242fa007 100644 --- a/tests/geometryXVx/collisions_intra_maxwellian.cpp +++ b/tests/geometryXVx/collisions_intra_maxwellian.cpp @@ -104,12 +104,13 @@ TEST(CollisionsIntraMaxwellian, CollisionsIntraMaxwellian) temperature_init(ispx) = temperature + temperature_ampl * std::sin(2 * M_PI * coordx / ddc::coordinate(gridx.back())); - DFieldVx finit(gridvx); + device_t finit_device(gridvx); MaxwellianEquilibrium::compute_maxwellian( - finit.span_view(), + finit_device.span_view(), density_init(ispx), temperature_init(ispx), mean_velocity_init(ispx)); + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); ddc::deepcopy(allfdistribu[ispx], finit); }); @@ -240,12 +241,13 @@ TEST(CollisionsIntraMaxwellian, CollisionsIntraMaxwellian) double const temperature_loc = temperature + temperature_ampl * std::sin(2 * M_PI * coordx / ddc::coordinate(gridx.back())); - DFieldVx finit(gridvx); + device_t finit_device(gridvx); MaxwellianEquilibrium::compute_maxwellian( - finit.span_view(), + finit_device.span_view(), density_loc, temperature_loc, mean_velocity_loc); + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); ddc::deepcopy(allfdistribu[ispx], finit); }); diff --git a/tests/geometryXVx/fluid_moments.cpp b/tests/geometryXVx/fluid_moments.cpp index 1c46b93a9..d096b1b38 100644 --- a/tests/geometryXVx/fluid_moments.cpp +++ b/tests/geometryXVx/fluid_moments.cpp @@ -93,12 +93,14 @@ TEST(Physics, FluidMoments) temperature_init(ispx) = temperature + temperature_ampl * std::sin(2 * M_PI * coordx / ddc::coordinate(gridx.back())); - DFieldVx finit(gridvx); + device_t finit_device(gridvx); MaxwellianEquilibrium::compute_maxwellian( - finit.span_view(), + finit_device.span_view(), density_init(ispx), temperature_init(ispx), mean_velocity_init(ispx)); + + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); ddc::deepcopy(allfdistribu[ispx], finit); }); diff --git a/tests/geometryXVx/krooksource.cpp b/tests/geometryXVx/krooksource.cpp index 18e8ff647..a6f37950f 100644 --- a/tests/geometryXVx/krooksource.cpp +++ b/tests/geometryXVx/krooksource.cpp @@ -108,16 +108,16 @@ TEST(KrookSource, Adaptive) double const temperature_init = 1.; DFieldSpXVx allfdistribu(mesh); ddc::for_each(ddc::get_domain(allfdistribu), [&](IndexSpX const ispx) { - DFieldVx finit(gridvx); + device_t finit_device(gridvx); if (charge(ddc::select(ispx)) >= 0.) { MaxwellianEquilibrium:: - compute_maxwellian(finit, density_init_ion, temperature_init, 0.); - ddc::deepcopy(allfdistribu[ispx], finit); + compute_maxwellian(finit_device, density_init_ion, temperature_init, 0.); } else { MaxwellianEquilibrium:: - compute_maxwellian(finit, density_init_elec, temperature_init, 0.); - ddc::deepcopy(allfdistribu[ispx], finit); + compute_maxwellian(finit_device, density_init_elec, temperature_init, 0.); } + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); + ddc::deepcopy(allfdistribu[ispx], finit); }); // error with a given deltat @@ -139,16 +139,16 @@ TEST(KrookSource, Adaptive) // reinitialization of the distribution function ddc::for_each(ddc::get_domain(allfdistribu), [&](IndexSpX const ispx) { - DFieldVx finit(gridvx); + device_t finit_device(gridvx); if (charge(ddc::select(ispx)) >= 0.) { MaxwellianEquilibrium:: - compute_maxwellian(finit, density_init_ion, temperature_init, 0.); - ddc::deepcopy(allfdistribu[ispx], finit); + compute_maxwellian(finit_device, density_init_ion, temperature_init, 0.); } else { MaxwellianEquilibrium:: - compute_maxwellian(finit, density_init_elec, temperature_init, 0.); - ddc::deepcopy(allfdistribu[ispx], finit); + compute_maxwellian(finit_device, density_init_elec, temperature_init, 0.); } + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); + ddc::deepcopy(allfdistribu[ispx], finit); }); // error with a deltat 10 times smaller @@ -252,8 +252,9 @@ TEST(KrookSource, Constant) // Initialization of the distribution function : maxwellian double const density_init = 1.; double const temperature_init = 1.; - DFieldVx finit(gridvx); - MaxwellianEquilibrium::compute_maxwellian(finit, density_init, temperature_init, 0.); + device_t finit_device(gridvx); + MaxwellianEquilibrium::compute_maxwellian(finit_device, density_init, temperature_init, 0.); + auto finit = ddc::create_mirror_view_and_copy(finit_device.span_view()); DFieldSpXVx allfdistribu(mesh); ddc::for_each(ddc::get_domain(allfdistribu), [&](IndexSpX const ispx) { ddc::deepcopy(allfdistribu[ispx], finit); @@ -265,8 +266,11 @@ TEST(KrookSource, Constant) }; // tests if distribution function matches theoretical prediction - DFieldVx ftarget(gridvx); - MaxwellianEquilibrium::compute_maxwellian(ftarget, density_target, temperature_target, 0.); + device_t ftarget_device(gridvx); + MaxwellianEquilibrium:: + compute_maxwellian(ftarget_device, density_target, temperature_target, 0.); + auto ftarget = ddc::create_mirror_view_and_copy(ftarget_device.span_view()); + ddc::for_each(allfdistribu.domain(), [&](IndexSpXVx const ispxvx) { // predicted distribution function value double const allfdistribu_pred