From a2aa7e5d1f688dfd54eef13923f5b6b90477f79c Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Tue, 19 Mar 2024 11:13:24 +0000 Subject: [PATCH] Landau4D GPU porting Landau4D is rebased on ddc splines and entirely ported on GPU. See merge request gysela-developpers/gyselalibxx!374 -------------------------------------------- Co-authored-by: Emily Bourne --- .../geometryXYVxVy/landau/landau4d_fft.cpp | 176 ++++++-------- src/geometryXYVxVy/geometry/CMakeLists.txt | 2 +- src/geometryXYVxVy/geometry/geometry.hpp | 219 ++++++++++++------ .../initialization/maxwellianequilibrium.cpp | 49 ++-- .../initialization/maxwellianequilibrium.hpp | 53 ++++- .../singlemodeperturbinitialization.cpp | 52 +++-- .../singlemodeperturbinitialization.hpp | 34 ++- src/geometryXYVxVy/poisson/CMakeLists.txt | 1 - .../poisson/chargedensitycalculator.cpp | 11 +- src/geometryXYVxVy/poisson/electricfield.cpp | 54 ----- src/geometryXYVxVy/poisson/electricfield.hpp | 29 --- .../poisson/fftpoissonsolver.cpp | 70 +++--- .../poisson/fftpoissonsolver.hpp | 4 - .../time_integration/predcorr.cpp | 33 ++- .../vlasov/splitvlasovsolver.cpp | 19 +- tests/geometryXYVxVy/fft_poisson.cpp | 114 +++++---- tests/geometryXYVxVy/quadrature.cpp | 26 ++- 17 files changed, 499 insertions(+), 447 deletions(-) delete mode 100644 src/geometryXYVxVy/poisson/electricfield.cpp delete mode 100644 src/geometryXYVxVy/poisson/electricfield.hpp diff --git a/simulations/geometryXYVxVy/landau/landau4d_fft.cpp b/simulations/geometryXYVxVy/landau/landau4d_fft.cpp index 707151638..dc063df83 100644 --- a/simulations/geometryXYVxVy/landau/landau4d_fft.cpp +++ b/simulations/geometryXYVxVy/landau/landau4d_fft.cpp @@ -10,17 +10,13 @@ #include -#include -#include -#include - -#include #include #include -#include "bsl_advection_vx.hpp" -#include "bsl_advection_x.hpp" +#include "bsl_advection_vx_batched.hpp" +#include "bsl_advection_x_batched.hpp" #include "fftpoissonsolver.hpp" +#include "geometry.hpp" #include "maxwellianequilibrium.hpp" #include "neumann_spline_quadrature.hpp" #include "paraconfpp.hpp" @@ -29,7 +25,7 @@ #include "predcorr.hpp" #include "singlemodeperturbinitialization.hpp" //#include "species_info.hpp" -#include "spline_interpolator.hpp" +#include "spline_interpolator_batched.hpp" #include "splitvlasovsolver.hpp" using std::cerr; @@ -37,25 +33,6 @@ using std::endl; using std::chrono::steady_clock; namespace fs = std::filesystem; -using PreallocatableSplineInterpolatorX - = PreallocatableSplineInterpolator; -using PreallocatableSplineInterpolatorY - = PreallocatableSplineInterpolator; -using PreallocatableSplineInterpolatorVx = PreallocatableSplineInterpolator< - IDimVx, - BSplinesVx, - BoundCond::HERMITE, - BoundCond::HERMITE>; -using PreallocatableSplineInterpolatorVy = PreallocatableSplineInterpolator< - IDimVy, - BSplinesVy, - BoundCond::HERMITE, - BoundCond::HERMITE>; -using BslAdvectionX = BslAdvectionSpatial; -using BslAdvectionY = BslAdvectionSpatial; -using BslAdvectionVx = BslAdvectionVelocity; -using BslAdvectionVy = BslAdvectionVelocity; - int main(int argc, char** argv) { // Environments variables for profiling @@ -102,58 +79,46 @@ int main(int argc, char** argv) // Creating mesh & supports ddc::init_discrete_space(x_min, x_max, x_ncells); - ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); - ddc::DiscreteDomain interpolation_domain_x(SplineInterpPointsX::get_domain()); - SplineXBuilder const builder_x(interpolation_domain_x); - ddc::init_discrete_space(y_min, y_max, y_ncells); - ddc::init_discrete_space(SplineInterpPointsY::get_sampling()); - ddc::DiscreteDomain interpolation_domain_y(SplineInterpPointsY::get_domain()); - SplineYBuilder const builder_y(interpolation_domain_y); - - ddc::DiscreteDomain - interpolation_domain_xy(interpolation_domain_x, interpolation_domain_y); - SplineXYBuilder const builder_xy(interpolation_domain_xy); - ddc::init_discrete_space(vx_min, vx_max, vx_ncells); - ddc::init_discrete_space(vx_min, vx_max, vx_ncells); - ddc::init_discrete_space(SplineInterpPointsVx::get_sampling()); - ddc::DiscreteDomain interpolation_domain_vx(SplineInterpPointsVx::get_domain()); - SplineVxBuilder const builder_vx(interpolation_domain_vx); - ddc::init_discrete_space(vy_min, vy_max, vy_ncells); - ddc::init_discrete_space(vy_min, vy_max, vy_ncells); + ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); + ddc::init_discrete_space(SplineInterpPointsY::get_sampling()); + ddc::init_discrete_space(SplineInterpPointsVx::get_sampling()); ddc::init_discrete_space(SplineInterpPointsVy::get_sampling()); - ddc::DiscreteDomain interpolation_domain_vy(SplineInterpPointsVy::get_domain()); - SplineVyBuilder const builder_vy(interpolation_domain_vy); + IVectSp const nb_kinspecies(PCpp_len(conf_voicexx, ".SpeciesInfo")); + IDomainSp const dom_kinsp(IndexSp(0), nb_kinspecies); + + IDomainX interpolation_domain_x(SplineInterpPointsX::get_domain()); + IDomainY interpolation_domain_y(SplineInterpPointsY::get_domain()); + IDomainVx interpolation_domain_vx(SplineInterpPointsVx::get_domain()); + IDomainVy interpolation_domain_vy(SplineInterpPointsVy::get_domain()); + IDomainVxVy interpolation_domain_vxvy(interpolation_domain_vx, interpolation_domain_vy); + + IDomainXYVxVy meshXYVxVy( + interpolation_domain_x, + interpolation_domain_y, + interpolation_domain_vx, + interpolation_domain_vy); + IDomainSpVxVy const meshSpVxVy(dom_kinsp, interpolation_domain_vx, interpolation_domain_vy); + IDomainSpXYVxVy const meshSpXYVxVy(dom_kinsp, meshXYVxVy); + + SplineXBuilder const builder_x(meshXYVxVy); + SplineYBuilder const builder_y(meshXYVxVy); + SplineVxBuilder const builder_vx(meshXYVxVy); + SplineVyBuilder const builder_vy(meshXYVxVy); SplineVxBuilder_1d const builder_vx_1d(interpolation_domain_vx); SplineVyBuilder_1d const builder_vy_1d(interpolation_domain_vy); - IVectSp const nb_kinspecies(PCpp_len(conf_voicexx, ".SpeciesInfo")); - IDomainSp const dom_kinsp(IndexSp(0), nb_kinspecies); + host_t> kinetic_charges(dom_kinsp); + host_t masses(dom_kinsp); + host_t density_eq(dom_kinsp); + host_t temperature_eq(dom_kinsp); + host_t mean_velocity_eq(dom_kinsp); + host_t init_perturb_amplitude(dom_kinsp); + host_t> init_perturb_mode(dom_kinsp); - IDomainSpXYVxVy const meshSpXYVxVy( - dom_kinsp, - builder_x.interpolation_domain(), - builder_y.interpolation_domain(), - builder_vx.interpolation_domain(), - builder_vy.interpolation_domain()); - - IDomainSpVxVy const meshSpVxVy( - dom_kinsp, - builder_vx.interpolation_domain(), - builder_vy.interpolation_domain()); - - IDomainVxVy const interpolation_domain_vxvy(interpolation_domain_vx, interpolation_domain_vy); - - FieldSp kinetic_charges(dom_kinsp); - DFieldSp masses(dom_kinsp); - DFieldSp density_eq(dom_kinsp); - DFieldSp temperature_eq(dom_kinsp); - DFieldSp mean_velocity_eq(dom_kinsp); - DFieldSp init_perturb_amplitude(dom_kinsp); - FieldSp init_perturb_mode(dom_kinsp); int nb_elec_adiabspecies = 1; int nb_ion_adiabspecies = 1; @@ -179,7 +144,7 @@ int main(int argc, char** argv) // Create the domain of all species including kinetic species + adiabatic species (if existing) IDomainSp const dom_allsp(IndexSp(0), nb_kinspecies + nb_elec_adiabspecies + nb_ion_adiabspecies); - FieldSp charges(dom_allsp); + host_t> charges(dom_allsp); for (IndexSp isp : dom_kinsp) { charges(isp) = kinetic_charges(isp); } @@ -201,6 +166,7 @@ int main(int argc, char** argv) init_perturb_mode.span_cview(), init_perturb_amplitude.span_cview()); init(allfdistribu); + auto allfequilibrium_host = ddc::create_mirror_view_and_copy(allfequilibrium.span_view()); // --> Algorithm info double const deltat = PCpp_double(conf_voicexx, ".Algorithm.deltat"); @@ -211,39 +177,39 @@ int main(int argc, char** argv) int const nbstep_diag = int(time_diag / deltat); // Create spline evaluator - ConstantExtrapolationBoundaryValue bv_x_min(x_min); - ConstantExtrapolationBoundaryValue bv_x_max(x_max); - SplineEvaluator const spline_x_evaluator(bv_x_min, bv_x_max); - PreallocatableSplineInterpolatorX const spline_x_interpolator(builder_x, spline_x_evaluator); - - ConstantExtrapolationBoundaryValue bv_y_min(y_min); - ConstantExtrapolationBoundaryValue bv_y_max(y_max); - SplineEvaluator const spline_y_evaluator(bv_y_min, bv_y_max); - PreallocatableSplineInterpolatorY const spline_y_interpolator(builder_y, spline_y_evaluator); - - ConstantExtrapolationBoundaryValue bv_vx_min(vx_min); - ConstantExtrapolationBoundaryValue bv_vx_max(vx_max); - SplineEvaluator const spline_vx_evaluator(bv_vx_min, bv_vx_max); - PreallocatableSplineInterpolatorVx const + ddc::PeriodicExtrapolationRule bv_x_min; + ddc::PeriodicExtrapolationRule bv_x_max; + SplineXEvaluator const spline_x_evaluator(builder_x.spline_domain(), bv_x_min, bv_x_max); + + PreallocatableSplineInterpolatorBatched const + spline_x_interpolator(builder_x, spline_x_evaluator); + + ddc::PeriodicExtrapolationRule bv_y_min; + ddc::PeriodicExtrapolationRule bv_y_max; + SplineYEvaluator const spline_y_evaluator(builder_y.spline_domain(), bv_y_min, bv_y_max); + + PreallocatableSplineInterpolatorBatched const + spline_y_interpolator(builder_y, spline_y_evaluator); + + ddc::ConstantExtrapolationRule bv_vx_min(vx_min); + ddc::ConstantExtrapolationRule bv_vx_max(vx_max); + SplineVxEvaluator const spline_vx_evaluator(builder_vx.spline_domain(), bv_vx_min, bv_vx_max); + + PreallocatableSplineInterpolatorBatched const spline_vx_interpolator(builder_vx, spline_vx_evaluator); - ConstantExtrapolationBoundaryValue bv_vy_min(vy_min); - ConstantExtrapolationBoundaryValue bv_vy_max(vy_max); - SplineEvaluator const spline_vy_evaluator(bv_vy_min, bv_vy_max); - PreallocatableSplineInterpolatorVy const - spline_vy_interpolator(builder_vy, spline_vy_evaluator); + ddc::ConstantExtrapolationRule bv_vy_min(vy_min); + ddc::ConstantExtrapolationRule bv_vy_max(vy_max); + SplineVyEvaluator const spline_vy_evaluator(builder_vy.spline_domain(), bv_vy_min, bv_vy_max); - SplineXYEvaluator const spline_xy_evaluator( - g_null_boundary_2d, - g_null_boundary_2d, - g_null_boundary_2d, - g_null_boundary_2d); + PreallocatableSplineInterpolatorBatched const + spline_vy_interpolator(builder_vy, spline_vy_evaluator); // Create advection operator - BslAdvectionX const advection_x(spline_x_interpolator); - BslAdvectionY const advection_y(spline_y_interpolator); - BslAdvectionVx const advection_vx(spline_vx_interpolator); - BslAdvectionVy const advection_vy(spline_vy_interpolator); + BslAdvectionSpatialBatched const advection_x(spline_x_interpolator); + BslAdvectionSpatialBatched const advection_y(spline_y_interpolator); + BslAdvectionVelocityBatched const advection_vx(spline_vx_interpolator); + BslAdvectionVelocityBatched const advection_vy(spline_vy_interpolator); SplitVlasovSolver const vlasov(advection_x, advection_y, advection_vx, advection_vy); @@ -263,22 +229,22 @@ int main(int argc, char** argv) PredCorr const predcorr(vlasov, poisson); // Creating of mesh for output saving - FieldX meshX_coord(interpolation_domain_x); + host_t> meshX_coord(interpolation_domain_x); ddc::for_each(interpolation_domain_x, [&](IndexX const ix) { meshX_coord(ix) = ddc::coordinate(ix); }); - FieldY meshY_coord(interpolation_domain_y); + host_t> meshY_coord(interpolation_domain_y); ddc::for_each(interpolation_domain_y, [&](IndexY const iy) { meshY_coord(iy) = ddc::coordinate(iy); }); - FieldVx meshVx_coord(interpolation_domain_vx); + host_t> meshVx_coord(interpolation_domain_vx); for (IndexVx const ivx : interpolation_domain_vx) { meshVx_coord(ivx) = ddc::coordinate(ivx); } - FieldVy meshVy_coord(interpolation_domain_vy); + host_t> meshVy_coord(interpolation_domain_vy); for (IndexVy const ivy : interpolation_domain_vy) { meshVy_coord(ivy) = ddc::coordinate(ivy); } @@ -296,11 +262,11 @@ int main(int argc, char** argv) ddc::expose_to_pdi("Nkinspecies", nb_kinspecies.value()); ddc::expose_to_pdi("fdistribu_charges", ddc::discrete_space().charges()[dom_kinsp]); ddc::expose_to_pdi("fdistribu_masses", ddc::discrete_space().masses()[dom_kinsp]); - ddc::PdiEvent("initial_state").with("fdistribu_eq", allfequilibrium); + ddc::PdiEvent("initial_state").with("fdistribu_eq", allfequilibrium_host); steady_clock::time_point const start = steady_clock::now(); - predcorr(allfdistribu, deltat, nbiter); + predcorr(allfdistribu.span_view(), deltat, nbiter); steady_clock::time_point const end = steady_clock::now(); diff --git a/src/geometryXYVxVy/geometry/CMakeLists.txt b/src/geometryXYVxVy/geometry/CMakeLists.txt index eca1f655e..307c13f1d 100644 --- a/src/geometryXYVxVy/geometry/CMakeLists.txt +++ b/src/geometryXYVxVy/geometry/CMakeLists.txt @@ -9,7 +9,7 @@ target_include_directories("geometry_xyvxvy" ) target_link_libraries("geometry_xyvxvy" INTERFACE DDC::DDC - sll::splines gslx::speciesinfo + gslx::utils ) add_library("gslx::geometry_xyvxvy" ALIAS "geometry_xyvxvy") diff --git a/src/geometryXYVxVy/geometry/geometry.hpp b/src/geometryXYVxVy/geometry/geometry.hpp index fd512f335..9c3335638 100644 --- a/src/geometryXYVxVy/geometry/geometry.hpp +++ b/src/geometryXYVxVy/geometry/geometry.hpp @@ -6,14 +6,7 @@ #include #include -#include -#include -#include -#include -#include -#include -#include - +#include #include /** @@ -81,38 +74,26 @@ bool constexpr BsplineOnUniformCellsVy = true; using BSplinesX = std::conditional_t< BsplineOnUniformCellsX, - UniformBSplines, - NonUniformBSplines>; + ddc::UniformBSplines, + ddc::NonUniformBSplines>; using BSplinesY = std::conditional_t< BsplineOnUniformCellsY, - UniformBSplines, - NonUniformBSplines>; + ddc::UniformBSplines, + ddc::NonUniformBSplines>; using BSplinesVx = std::conditional_t< - BsplineOnUniformCellsVx, - UniformBSplines, - NonUniformBSplines>; -using BSplinesVy = std::conditional_t< - BsplineOnUniformCellsVy, - UniformBSplines, - NonUniformBSplines>; -using DDCBSplinesVx = std::conditional_t< BsplineOnUniformCellsVx, ddc::UniformBSplines, ddc::NonUniformBSplines>; -using DDCBSplinesVy = std::conditional_t< +using BSplinesVy = std::conditional_t< BsplineOnUniformCellsVy, ddc::UniformBSplines, ddc::NonUniformBSplines>; -BoundCond constexpr SplineXBoundary = BoundCond::PERIODIC; -BoundCond constexpr SplineYBoundary = BoundCond::PERIODIC; -BoundCond constexpr SplineVxBoundary = BoundCond::HERMITE; -BoundCond constexpr SplineVyBoundary = BoundCond::HERMITE; -ddc::BoundCond constexpr DDCSplineXBoundary = ddc::BoundCond::PERIODIC; -ddc::BoundCond constexpr DDCSplineYBoundary = ddc::BoundCond::PERIODIC; -ddc::BoundCond constexpr DDCSplineVxBoundary = ddc::BoundCond::HERMITE; -ddc::BoundCond constexpr DDCSplineVyBoundary = ddc::BoundCond::HERMITE; +ddc::BoundCond constexpr SplineXBoundary = ddc::BoundCond::PERIODIC; +ddc::BoundCond constexpr SplineYBoundary = ddc::BoundCond::PERIODIC; +ddc::BoundCond constexpr SplineVxBoundary = ddc::BoundCond::HERMITE; +ddc::BoundCond constexpr SplineVyBoundary = ddc::BoundCond::HERMITE; bool constexpr UniformMeshX = is_spline_interpolation_mesh_uniform( BsplineOnUniformCellsX, @@ -156,40 +137,124 @@ using IDimVy = std::conditional_t< // IDim initialisers using SplineInterpPointsX - = GrevilleInterpolationPoints; + = ddc::GrevilleInterpolationPoints; using SplineInterpPointsY - = GrevilleInterpolationPoints; + = ddc::GrevilleInterpolationPoints; using SplineInterpPointsVx - = GrevilleInterpolationPoints; + = ddc::GrevilleInterpolationPoints; using SplineInterpPointsVy - = GrevilleInterpolationPoints; + = ddc::GrevilleInterpolationPoints; // SplineBuilder and SplineEvaluator definition -using SplineXBuilder = SplineBuilder; -using SplineYBuilder = SplineBuilder; -using SplineXYBuilder = SplineBuilder2D; -using SplineVxBuilder = SplineBuilder; -using SplineVyBuilder = SplineBuilder; -using SplineVxVyBuilder = SplineBuilder2D; -using SplineXYEvaluator = SplineEvaluator2D; -using SplineVxVyEvaluator = SplineEvaluator2D; +using SplineXBuilder = ddc::SplineBuilder< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesX, + IDimX, + SplineXBoundary, + SplineXBoundary, + ddc::SplineSolver::GINKGO, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineXEvaluator = ddc::SplineEvaluator< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesX, + IDimX, + ddc::PeriodicExtrapolationRule, + ddc::PeriodicExtrapolationRule, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineYBuilder = ddc::SplineBuilder< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesY, + IDimY, + SplineYBoundary, + SplineYBoundary, + ddc::SplineSolver::GINKGO, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineYEvaluator = ddc::SplineEvaluator< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesY, + IDimY, + ddc::PeriodicExtrapolationRule, + ddc::PeriodicExtrapolationRule, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineVxBuilder = ddc::SplineBuilder< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesVx, + IDimVx, + SplineVxBoundary, + SplineVxBoundary, + ddc::SplineSolver::GINKGO, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineVxEvaluator = ddc::SplineEvaluator< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesVx, + IDimVx, + ddc::ConstantExtrapolationRule, + ddc::ConstantExtrapolationRule, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineVyBuilder = ddc::SplineBuilder< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesVy, + IDimVy, + SplineVyBoundary, + SplineVyBoundary, + ddc::SplineSolver::GINKGO, + IDimX, + IDimY, + IDimVx, + IDimVy>; +using SplineVyEvaluator = ddc::SplineEvaluator< + Kokkos::DefaultExecutionSpace, + Kokkos::DefaultExecutionSpace::memory_space, + BSplinesVy, + IDimVy, + ddc::ConstantExtrapolationRule, + ddc::ConstantExtrapolationRule, + IDimX, + IDimY, + IDimVx, + IDimVy>; using SplineVxBuilder_1d = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, - DDCBSplinesVx, + BSplinesVx, IDimVx, - DDCSplineVxBoundary, - DDCSplineVxBoundary, + SplineVxBoundary, + SplineVxBoundary, ddc::SplineSolver::GINKGO, IDimVx>; using SplineVyBuilder_1d = ddc::SplineBuilder< Kokkos::DefaultHostExecutionSpace, Kokkos::DefaultHostExecutionSpace::memory_space, - DDCBSplinesVy, + BSplinesVy, IDimVy, - DDCSplineVyBoundary, - DDCSplineVyBoundary, + SplineVyBoundary, + SplineVyBoundary, ddc::SplineSolver::GINKGO, IDimVy>; @@ -201,7 +266,7 @@ using BSDomainVy = ddc::DiscreteDomain; using BSDomainVxVy = ddc::DiscreteDomain; template -using BSViewXY = ddc::ChunkSpan; +using BSViewXY = device_t>; using DBSViewXY = BSViewXY; // Index @@ -211,7 +276,9 @@ using IndexY = ddc::DiscreteElement; using IndexXY = ddc::DiscreteElement; using IndexVx = ddc::DiscreteElement; using IndexVy = ddc::DiscreteElement; +using IndexVxVy = ddc::DiscreteElement; using IndexXYVxVy = ddc::DiscreteElement; +using IndexSpXYVxVy = ddc::DiscreteElement; // IVect definition using IVectSp = ddc::DiscreteVector; @@ -234,103 +301,107 @@ using IDomainSpXYVxVy = ddc::DiscreteDomain -using FieldSp = ddc::Chunk; +using FieldSp = device_t>; using DFieldSp = FieldSp; template -using FieldX = ddc::Chunk; +using FieldX = device_t>; using DFieldX = FieldX; template -using FieldY = ddc::Chunk; +using FieldY = device_t>; using DFieldY = FieldY; template -using FieldXY = ddc::Chunk; +using FieldXY = device_t>; using DFieldXY = FieldXY; template -using FieldVx = ddc::Chunk; +using FieldVx = device_t>; template -using FieldVy = ddc::Chunk; +using FieldVy = device_t>; template -using FieldVxVy = ddc::Chunk; +using FieldVxVy = device_t>; using DFieldVxVy = FieldVxVy; template -using FieldSpVxVy = ddc::Chunk; +using FieldXYVxVy = device_t>; +using DFieldXYVxVy = FieldXYVxVy; + +template +using FieldSpVxVy = device_t>; using DFieldSpVxVy = FieldSpVxVy; template -using FieldSpXYVxVy = ddc::Chunk; +using FieldSpXYVxVy = device_t>; using DFieldSpXYVxVy = FieldSpXYVxVy; // Span definition template -using SpanX = ddc::ChunkSpan; +using SpanX = device_t>; using DSpanX = SpanX; template -using SpanY = ddc::ChunkSpan; +using SpanY = device_t>; using DSpanY = SpanY; template -using SpanXY = ddc::ChunkSpan; +using SpanXY = device_t>; using DSpanXY = SpanXY; template -using SpanVx = ddc::ChunkSpan; +using SpanVx = device_t>; using DSpanVx = SpanVx; template -using SpanVy = ddc::ChunkSpan; +using SpanVy = device_t>; using DSpanVy = SpanVy; template -using SpanVxVy = ddc::ChunkSpan; +using SpanVxVy = device_t>; using DSpanVxVy = SpanVxVy; template -using SpanSpVxVy = ddc::ChunkSpan; +using SpanSpVxVy = device_t>; using DSpanSpVxVy = SpanSpVxVy; template -using SpanSpXYVxVy = ddc::ChunkSpan; +using SpanSpXYVxVy = device_t>; using DSpanSpXYVxVy = SpanSpXYVxVy; // View definition template -using ViewSp = ddc::ChunkSpan; +using ViewSp = device_t>; using DViewSp = ViewSp; template -using ViewX = ddc::ChunkSpan; +using ViewX = device_t>; template -using ViewY = ddc::ChunkSpan; +using ViewY = device_t>; template -using ViewXY = ddc::ChunkSpan; +using ViewXY = device_t>; using DViewXY = ViewXY; template -using ViewVx = ddc::ChunkSpan; +using ViewVx = device_t>; template -using ViewVy = ddc::ChunkSpan; +using ViewVy = device_t>; template -using ViewVxVy = ddc::ChunkSpan; +using ViewVxVy = device_t>; using DViewVxVy = ViewVxVy; template -using ViewSpVxVy = ddc::ChunkSpan; +using ViewSpVxVy = device_t>; using DViewSpVxVy = ViewSpVxVy; template -using ViewSpXYVxVy = ddc::ChunkSpan; +using ViewSpXYVxVy = device_t>; using DViewSpXYVxVy = ViewSpXYVxVy; // For Fourier diff --git a/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp b/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp index 94dd868fe..ad60c11a6 100644 --- a/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp +++ b/src/geometryXYVxVy/initialization/maxwellianequilibrium.cpp @@ -5,9 +5,9 @@ #include "maxwellianequilibrium.hpp" MaxwellianEquilibrium::MaxwellianEquilibrium( - DFieldSp density_eq, - DFieldSp temperature_eq, - DFieldSp mean_velocity_eq) + host_t density_eq, + host_t temperature_eq, + host_t mean_velocity_eq) : m_density_eq(std::move(density_eq)) , m_temperature_eq(std::move(temperature_eq)) , m_mean_velocity_eq(std::move(mean_velocity_eq)) @@ -22,7 +22,8 @@ DSpanSpVxVy MaxwellianEquilibrium::operator()(DSpanSpVxVy const allfequilibrium) IDomainVxVy const gridvxvy = allfequilibrium.domain(); // Initialization of the maxwellian - DFieldVxVy maxwellian(gridvxvy); + DFieldVxVy maxwellian_alloc(gridvxvy); + DSpanVxVy maxwellian = maxwellian_alloc.span_view(); ddc::for_each(gridsp, [&](IndexSp const isp) { compute_maxwellian( maxwellian, @@ -30,11 +31,12 @@ DSpanSpVxVy MaxwellianEquilibrium::operator()(DSpanSpVxVy const allfequilibrium) m_temperature_eq(isp), m_mean_velocity_eq(isp)); - ddc::for_each(gridvy, [&](IndexVy const ivy) { - ddc::for_each(gridvx, [&](IndexVx const ivx) { - allfequilibrium(isp, ivx, ivy) = maxwellian(ivx, ivy); - }); - }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + gridvxvy, + KOKKOS_LAMBDA(IndexVxVy const ivxvy) { + allfequilibrium(isp, ivxvy) = maxwellian(ivxvy); + }); }); return allfequilibrium; } @@ -51,19 +53,20 @@ void MaxwellianEquilibrium::compute_maxwellian( double const mean_velocity) { double const inv_2pi = 1. / (2. * M_PI * temperature); - IDomainVx const gridvx = fMaxwellian.domain(); - IDomainVy const gridvy = fMaxwellian.domain(); + IDomainVxVy const gridvxvy = fMaxwellian.domain(); - for (IndexVx const ivx : gridvx) { - double const vx = ddc::coordinate(ivx); - for (IndexVy const ivy : gridvy) { - double const vy = ddc::coordinate(ivy); - fMaxwellian(ivx, ivy) - = density * inv_2pi - * std::exp( - -((vx - std::sqrt(mean_velocity)) * (vx - std::sqrt(mean_velocity)) - + (vy - std::sqrt(mean_velocity)) * (vy - std::sqrt(mean_velocity))) - / (2. * temperature)); - } - } + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + gridvxvy, + KOKKOS_LAMBDA(IndexVxVy const ivxvy) { + double const vx = ddc::coordinate(ddc::select(ivxvy)); + double const vy = ddc::coordinate(ddc::select(ivxvy)); + fMaxwellian(ivxvy) = density * inv_2pi + * Kokkos::exp( + -((vx - Kokkos::sqrt(mean_velocity)) + * (vx - Kokkos::sqrt(mean_velocity)) + + (vy - Kokkos::sqrt(mean_velocity)) + * (vy - Kokkos::sqrt(mean_velocity))) + / (2. * temperature)); + }); } diff --git a/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp b/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp index 59882880e..41656504c 100644 --- a/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp +++ b/src/geometryXYVxVy/initialization/maxwellianequilibrium.hpp @@ -11,27 +11,46 @@ class MaxwellianEquilibrium : public IEquilibrium { // equilibrium density of all kinetic species - FieldSp m_density_eq; + host_t> m_density_eq; // equilibrium temperature of all kinetic species - FieldSp m_temperature_eq; + host_t> m_temperature_eq; // equilibrium mean velocity of all kinetic species - FieldSp m_mean_velocity_eq; + host_t> m_mean_velocity_eq; public: - MaxwellianEquilibrium(DFieldSp density_eq, DFieldSp temperature_eq, DFieldSp mean_velocity_eq); + /** + * @brief The constructor for the MaxwellianEquilibrium class. + * @param[in] density_eq The density of the Maxwellian + * @param[in] temperature_eq The temperature of the Maxwellian + * @param[in] mean_velocity_eq The mean velocity of the Maxwellian + */ + MaxwellianEquilibrium( + host_t density_eq, + host_t temperature_eq, + host_t mean_velocity_eq); ~MaxwellianEquilibrium() override = default; + /** + * @brief Initializes allfequilibrium as a Maxwellian. + * @param[out] allfequilibrium A Span containing a Maxwellian distribution function. + * @return A Span containing a Maxwellian distribution function. + */ DSpanSpVxVy operator()(DSpanSpVxVy allfequilibrium) const override; /** - * Computing the non-centered Maxwellian function as - * fM(v) = n/(2*PI*T)*exp(-(v-u)**2/(2*T)) - * with n the density and T the temperature and - * where u is the mean velocity + * @brief Compute a Maxwellian distribution function. + * The Maxwellian distribution function is defined as + * $f_M(v) = n/(sqrt(2*PI*T))*exp(-(v-u)**2/(2*T))$ + * with $n$ the density, $T$ the temperature and + * $u$ is the mean velocity. + * @param[out] fMaxwellian A Maxwellian distribution function. + * @param[in] density A parameter that represents the density of Maxwellian. + * @param[in] temperature A parameter that represents the temperature of Maxwellian. + * @param[in] mean_velocity A parameter that represents the mean velocity of Maxwellian. */ static void compute_maxwellian( DSpanVxVy const fMaxwellian, @@ -39,17 +58,29 @@ class MaxwellianEquilibrium : public IEquilibrium double const temperature, double const mean_velocity); - ViewSp density_eq() const + /** + * @brief A method for accessing the m_density_eq member variable of the class. + * @return A view containing the m_density_eq value. + */ + host_t> density_eq() const { return m_density_eq; } - ViewSp temperature_eq() const + /** + * @brief A method for accessing the m_temperature_eq member variable of the class. + * @return A view containing the m_temperature_eq value. + */ + host_t> temperature_eq() const { return m_temperature_eq; } - ViewSp mean_velocity_eq() const + /** + * @brief A method for accessing the m_mean_velocity_eq member variable of the class. + * @return A view containing the m_velocity_eq value. + */ + host_t> mean_velocity_eq() const { return m_mean_velocity_eq; } diff --git a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp index ba14c1db1..84048c5e0 100644 --- a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp +++ b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.cpp @@ -9,8 +9,8 @@ SingleModePerturbInitialization::SingleModePerturbInitialization( DViewSpVxVy fequilibrium, - ViewSp const init_perturb_mode, - DViewSp const init_perturb_amplitude) + host_t> const init_perturb_mode, + host_t const init_perturb_amplitude) : m_fequilibrium(fequilibrium) , m_init_perturb_mode(init_perturb_mode) , m_init_perturb_amplitude(init_perturb_amplitude) @@ -24,7 +24,8 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al IDomainXYVxVy const gridxyvxvy = allfdistribu.domain(); // Initialization of the perturbation - DFieldXY perturbation(gridxy); + DFieldXY perturbation_alloc(gridxy); + DSpanXY perturbation = perturbation_alloc.span_view(); ddc::for_each(gridsp, [&](IndexSp const isp) { perturbation_initialization( perturbation, @@ -32,17 +33,22 @@ DSpanSpXYVxVy SingleModePerturbInitialization::operator()(DSpanSpXYVxVy const al m_init_perturb_amplitude(isp)); // Initialization of the distribution function --> fill values - ddc::for_each(gridxyvxvy, [&](IndexXYVxVy const ixyvxvy) { - IndexX const ix = ddc::select(ixyvxvy); - IndexY const iy = ddc::select(ixyvxvy); - IndexVx const ivx = ddc::select(ixyvxvy); - IndexVy const ivy = ddc::select(ixyvxvy); - double fdistribu_val = m_fequilibrium(isp, ivx, ivy) * (1. + perturbation(ix, iy)); - if (fdistribu_val < 1.e-60) { - fdistribu_val = 1.e-60; - } - allfdistribu(isp, ix, iy, ivx, ivy) = fdistribu_val; - }); + DViewSpVxVy fequilibrium_proxy = m_fequilibrium; + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + gridxyvxvy, + KOKKOS_LAMBDA(IndexXYVxVy const ixyvxvy) { + IndexX const ix = ddc::select(ixyvxvy); + IndexY const iy = ddc::select(ixyvxvy); + IndexVx const ivx = ddc::select(ixyvxvy); + IndexVy const ivy = ddc::select(ixyvxvy); + double fdistribu_val + = fequilibrium_proxy(isp, ivx, ivy) * (1. + perturbation(ix, iy)); + if (fdistribu_val < 1.e-60) { + fdistribu_val = 1.e-60; + } + allfdistribu(isp, ix, iy, ivx, ivy) = fdistribu_val; + }); }); return allfdistribu; } @@ -58,11 +64,15 @@ void SingleModePerturbInitialization::perturbation_initialization( double const ky = mode * 2. * M_PI / ddcHelper::total_interval_length(ddc::select(gridxy)); - ddc::for_each(gridxy, [&](IndexXY const ixy) { - IndexX const ix = ddc::select(ixy); - CoordX const x = ddc::coordinate(ix); - IndexY const iy = ddc::select(ixy); - CoordY const y = ddc::coordinate(iy); - perturbation(ix, iy) = perturb_amplitude * (cos(kx * x) + cos(ky * y)); - }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + gridxy, + KOKKOS_LAMBDA(IndexXY const ixy) { + IndexX const ix = ddc::select(ixy); + CoordX const x = ddc::coordinate(ix); + IndexY const iy = ddc::select(ixy); + CoordY const y = ddc::coordinate(iy); + perturbation(ix, iy) + = perturb_amplitude * (Kokkos::cos(kx * x) + Kokkos::cos(ky * y)); + }); } diff --git a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp index 9651132e3..dd99f00fd 100644 --- a/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp +++ b/src/geometryXYVxVy/initialization/singlemodeperturbinitialization.hpp @@ -12,25 +12,39 @@ class SingleModePerturbInitialization : public IInitialization { DViewSpVxVy m_fequilibrium; - ViewSp m_init_perturb_mode; + host_t> m_init_perturb_mode; - DViewSp m_init_perturb_amplitude; + host_t m_init_perturb_amplitude; -private: - /* - Initialization of the perturbation - */ +public: + /** + * @brief Initialization of the perturbation. + * @param[in, out] perturbation On input: an uninitialized array + * On output: an array containing a values that has a + * sinusoidal variation with given amplitude and mode. + * @param[in] mode The mode of the perturbation. + * @param[in] perturb_amplitude The amplitude of the perturbation. + */ void perturbation_initialization(DSpanXY perturbation, int mode, double perturb_amplitude) const; -public: - // Warning: all variables shall remain valid until the last call to `operator()` + /** + * @brief Creates an instance of the SingleModePerturbInitialization class. + * @param[in] fequilibrium A Maxwellian. + * @param[in] init_perturb_mode The perturbation mode. + * @param[in] init_perturb_amplitude The perturbation amplitude. + */ SingleModePerturbInitialization( DViewSpVxVy fequilibrium, - ViewSp init_perturb_mode, - DViewSp init_perturb_amplitude); + host_t> init_perturb_mode, + host_t init_perturb_amplitude); ~SingleModePerturbInitialization() override = default; + /** + * @brief Initializes the distribution function as as a perturbed Maxwellian. + * @param[in, out] allfdistribu The initialized distribution function. + * @return The initialized distribution function. + */ DSpanSpXYVxVy operator()(DSpanSpXYVxVy allfdistribu) const override; }; diff --git a/src/geometryXYVxVy/poisson/CMakeLists.txt b/src/geometryXYVxVy/poisson/CMakeLists.txt index dab207568..d4b25c852 100644 --- a/src/geometryXYVxVy/poisson/CMakeLists.txt +++ b/src/geometryXYVxVy/poisson/CMakeLists.txt @@ -2,7 +2,6 @@ add_library("poisson_xy" STATIC chargedensitycalculator.cpp - electricfield.cpp nullpoissonsolver.cpp fftpoissonsolver.cpp ) diff --git a/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp b/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp index db33f2f3e..1c6ebcd33 100644 --- a/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp +++ b/src/geometryXYVxVy/poisson/chargedensitycalculator.cpp @@ -12,11 +12,8 @@ ChargeDensityCalculator::ChargeDensityCalculator(const ChunkViewType& coeffs) { } -void ChargeDensityCalculator::operator()(DSpanXY rho_host, DViewSpXYVxVy allfdistribu_host) const +void ChargeDensityCalculator::operator()(DSpanXY rho, DViewSpXYVxVy allfdistribu) const { - auto rho = create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), rho_host); - auto allfdistribu - = create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), allfdistribu_host); Kokkos::Profiling::pushRegion("ChargeDensityCalculator"); IndexSp const last_kin_species = allfdistribu.domain().back(); IndexSp const last_species = ddc::discrete_space().charges().domain().back(); @@ -30,8 +27,8 @@ void ChargeDensityCalculator::operator()(DSpanXY rho_host, DViewSpXYVxVy allfdis = allfdistribu.allocation_kokkos_view(); Kokkos::View const rho_view = rho.allocation_kokkos_view(); - ViewSp const charges_host = ddc::host_discrete_space().charges(); - ViewSp const kinetic_charges_host = charges_host[allfdistribu.domain()]; + host_t> const charges_host = ddc::host_discrete_space().charges(); + host_t> const kinetic_charges_host = charges_host[allfdistribu.domain()]; auto charges_alloc = create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), kinetic_charges_host); @@ -83,6 +80,4 @@ void ChargeDensityCalculator::operator()(DSpanXY rho_host, DViewSpXYVxVy allfdis rho_view(ix, iy) = chargedens_adiabspecies + teamSum; }); Kokkos::Profiling::popRegion(); - - ddc::parallel_deepcopy(rho_host, rho); } diff --git a/src/geometryXYVxVy/poisson/electricfield.cpp b/src/geometryXYVxVy/poisson/electricfield.cpp deleted file mode 100644 index 8686d4581..000000000 --- a/src/geometryXYVxVy/poisson/electricfield.cpp +++ /dev/null @@ -1,54 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include - -#include - -#include "electricfield.hpp" - -ElectricField::ElectricField( - SplineXYBuilder const& spline_xy_builder, - SplineXYEvaluator const& spline_xy_evaluator) - : m_spline_xy_builder(spline_xy_builder) - , m_spline_xy_evaluator(spline_xy_evaluator) -{ -} - -//=========================================================================== -// Compute efield = -dPhi/dx where Phi is the electrostatic potential -// input : Phi values -//=========================================================================== -void ElectricField::operator()( - DSpanXY const electric_field_x, - DSpanXY const electric_field_y, - DBSViewXY const electrostatic_potential) const -{ - IDomainXY const& xy_dom = electric_field_x.domain(); - ddc::for_each(xy_dom, [&](IndexXY const ixy) { - IndexX const ix = ddc::select(ixy); - IndexY const iy = ddc::select(ixy); - - auto coordx = ddc::coordinate(ix); - auto coordy = ddc::coordinate(iy); - - electric_field_x(ix, iy) - = -m_spline_xy_evaluator - .deriv_dim_1(CoordXY(coordx, coordy), electrostatic_potential); - electric_field_y(ix, iy) - = -m_spline_xy_evaluator - .deriv_dim_2(CoordXY(coordx, coordy), electrostatic_potential); - }); -} - - -void ElectricField::operator()( - DSpanXY const electric_field_x, - DSpanXY const electric_field_y, - DViewXY const electrostatic_potential) const -{ - Kokkos::Profiling::pushRegion("PoissonSolver"); - ddc::Chunk elecpot_spline_coef((m_spline_xy_builder.spline_domain())); - m_spline_xy_builder(elecpot_spline_coef, electrostatic_potential); - (*this)(electric_field_x, electric_field_y, elecpot_spline_coef); - Kokkos::Profiling::popRegion(); -} diff --git a/src/geometryXYVxVy/poisson/electricfield.hpp b/src/geometryXYVxVy/poisson/electricfield.hpp deleted file mode 100644 index e4561908d..000000000 --- a/src/geometryXYVxVy/poisson/electricfield.hpp +++ /dev/null @@ -1,29 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include -#include - -//=========================================================================== -// Compute efield = -dPhi/dx where Phi is the electrostatic potential -// input : Phi values -//=========================================================================== -class ElectricField -{ - SplineXYBuilder const& m_spline_xy_builder; - - SplineXYEvaluator m_spline_xy_evaluator; - -public: - ElectricField( - SplineXYBuilder const& spline_xy_builder, - SplineXYEvaluator const& spline_xy_evaluator); - - void operator()( - DSpanXY electric_field_x, - DSpanXY electric_field_y, - DBSViewXY electrostatic_potential) const; - void operator()( - DSpanXY electric_field_x, - DSpanXY electric_field_y, - DViewXY electrostatic_potential) const; -}; diff --git a/src/geometryXYVxVy/poisson/fftpoissonsolver.cpp b/src/geometryXYVxVy/poisson/fftpoissonsolver.cpp index b190402f7..89517f698 100644 --- a/src/geometryXYVxVy/poisson/fftpoissonsolver.cpp +++ b/src/geometryXYVxVy/poisson/fftpoissonsolver.cpp @@ -34,67 +34,81 @@ void FftPoissonSolver::operator()( IDomainXY const xy_dom = electrostatic_potential.domain(); // Compute the RHS of the Poisson equation. - ddc::Chunk rho(xy_dom); + DFieldXY rho(xy_dom); DFieldVxVy contiguous_slice_vxvy(allfdistribu.domain()); m_compute_rho(rho, allfdistribu); // Build a mesh in the fourier space, for N points IDomainFxFy const k_mesh = ddc::FourierMesh(xy_dom, false); - ddc::Chunk, IDomainFxFy> intermediate_chunk - = ddc::Chunk(k_mesh, ddc::HostAllocator>()); + device_t, IDomainFxFy>> intermediate_chunk_alloc(k_mesh); + ddc::ChunkSpan intermediate_chunk = intermediate_chunk_alloc.span_view(); // Compute FFT(rho) ddc::FFT_Normalization norm = ddc::FFT_Normalization::BACKWARD; ddc:: - fft(Kokkos::DefaultHostExecutionSpace(), + fft(Kokkos::DefaultExecutionSpace(), intermediate_chunk.span_view(), rho.span_view(), ddc::kwArgs_fft {norm}); // Solve Poisson's equation -d2Phi/dx2 = rho // in Fourier space as -kx*kx*FFT(Phi)=FFT(rho)) - intermediate_chunk(k_mesh.front()) = 0.; - ddc::for_each(k_mesh, [&](IndexFxFy const ikxky) { - IndexFx const ikx = ddc::select(ikxky); - IndexFy const iky = ddc::select(ikxky); - - if (ikxky != k_mesh.front()) { - intermediate_chunk(ikxky) = intermediate_chunk(ikxky) - / ((double)coordinate(ikx) * (double)coordinate(ikx) - + (double)coordinate(iky) * (double)coordinate(iky)); - } - }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + k_mesh, + KOKKOS_LAMBDA(IndexFxFy const ikxky) { + IndexFx const ikx = ddc::select(ikxky); + IndexFy const iky = ddc::select(ikxky); + + if (ikxky != k_mesh.front()) { + intermediate_chunk(ikxky) + = intermediate_chunk(ikxky) + / ((double)coordinate(ikx) * (double)coordinate(ikx) + + (double)coordinate(iky) * (double)coordinate(iky)); + } else { + intermediate_chunk(ikxky) = 0.; + } + }); // Find the electric field in Fourier space // FFT(efield) = [-kx*i*FFT(Phi), -ky*i*FFT(Phi)] Kokkos::Profiling::pushRegion("ElectricField"); Kokkos::complex imaginary_unit(0.0, 1.0); - ddc::Chunk, IDomainFxFy> fourier_efield - = ddc::Chunk(k_mesh, ddc::HostAllocator>()); + ddc::Chunk fourier_efield_alloc + = ddc::Chunk(k_mesh, ddc::DeviceAllocator>()); + ddc::ChunkSpan fourier_efield = fourier_efield_alloc.span_view(); // Calculate x component - ddc::for_each(k_mesh, [&](IndexFxFy const ikxky) { - IndexFx const ikx = ddc::select(ikxky); - fourier_efield(ikxky) = -imaginary_unit * coordinate(ikx) * intermediate_chunk(ikxky); - }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + k_mesh, + KOKKOS_LAMBDA(IndexFxFy const ikxky) { + IndexFx const ikx = ddc::select(ikxky); + fourier_efield(ikxky) + = -imaginary_unit * coordinate(ikx) * intermediate_chunk(ikxky); + }); // Perform the inverse 1D FFT of the solution to deduce the electric field ddc:: - ifft(Kokkos::DefaultHostExecutionSpace(), + ifft(Kokkos::DefaultExecutionSpace(), electric_field_x.span_view(), fourier_efield.span_view(), ddc::kwArgs_fft {norm}); // Calculate y component - ddc::for_each(k_mesh, [&](IndexFxFy const ikxky) { - IndexFy const iky = ddc::select(ikxky); - fourier_efield(ikxky) = -imaginary_unit * coordinate(iky) * intermediate_chunk(ikxky); - }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + k_mesh, + KOKKOS_LAMBDA(IndexFxFy const ikxky) { + IndexFy const iky = ddc::select(ikxky); + fourier_efield(ikxky) + = -imaginary_unit * coordinate(iky) * intermediate_chunk(ikxky); + }); // Perform the inverse 1D FFT of the solution to deduce the electric field ddc:: - ifft(Kokkos::DefaultHostExecutionSpace(), + ifft(Kokkos::DefaultExecutionSpace(), electric_field_y.span_view(), fourier_efield.span_view(), ddc::kwArgs_fft {norm}); @@ -102,7 +116,7 @@ void FftPoissonSolver::operator()( // Perform the inverse 1D FFT of the solution to deduce the electrostatic potential ddc:: - ifft(Kokkos::DefaultHostExecutionSpace(), + ifft(Kokkos::DefaultExecutionSpace(), electrostatic_potential.span_view(), intermediate_chunk.span_view(), ddc::kwArgs_fft {norm}); diff --git a/src/geometryXYVxVy/poisson/fftpoissonsolver.hpp b/src/geometryXYVxVy/poisson/fftpoissonsolver.hpp index 187abbb72..d75ed3b9a 100644 --- a/src/geometryXYVxVy/poisson/fftpoissonsolver.hpp +++ b/src/geometryXYVxVy/poisson/fftpoissonsolver.hpp @@ -2,11 +2,7 @@ #pragma once -#include -#include - #include "chargedensitycalculator.hpp" -#include "electricfield.hpp" #include "ipoissonsolver.hpp" /** diff --git a/src/geometryXYVxVy/time_integration/predcorr.cpp b/src/geometryXYVxVy/time_integration/predcorr.cpp index e4286b0ca..4cbf8334f 100644 --- a/src/geometryXYVxVy/time_integration/predcorr.cpp +++ b/src/geometryXYVxVy/time_integration/predcorr.cpp @@ -21,15 +21,24 @@ DSpanSpXYVxVy PredCorr::operator()( double const dt, int const steps) const { + auto allfdistribu_host_alloc = ddc::create_mirror_view_and_copy(allfdistribu); + ddc::ChunkSpan allfdistribu_host = allfdistribu_host_alloc.span_view(); + // electrostatic potential and electric field (depending only on x) DFieldXY electrostatic_potential(allfdistribu.domain()); DFieldXY electric_field_x(allfdistribu.domain()); DFieldXY electric_field_y(allfdistribu.domain()); + host_t electrostatic_potential_host(allfdistribu.domain()); + // a 2D chunck of the same size as fdistribu DFieldSpXYVxVy allfdistribu_half_t(allfdistribu.domain()); - m_poisson_solver(electrostatic_potential, electric_field_x, electric_field_y, allfdistribu); + m_poisson_solver( + electrostatic_potential, + electric_field_x, + electric_field_y, + allfdistribu.span_cview()); int iter = 0; for (; iter < steps; ++iter) { @@ -37,13 +46,19 @@ DSpanSpXYVxVy PredCorr::operator()( // computation of the electrostatic potential at time tn and // the associated electric field - m_poisson_solver(electrostatic_potential, electric_field_x, electric_field_y, allfdistribu); - + m_poisson_solver( + electrostatic_potential, + electric_field_x, + electric_field_y, + allfdistribu.span_cview()); + // copies necessary to PDI + ddc::parallel_deepcopy(allfdistribu_host, allfdistribu); + ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); ddc::PdiEvent("iteration") .with("iter", iter) .and_with("time_saved", iter_time) - .and_with("fdistribu", allfdistribu) - .and_with("electrostatic_potential", electrostatic_potential); + .and_with("fdistribu", allfdistribu_host) + .and_with("electrostatic_potential", electrostatic_potential_host); // copy fdistribu ddc::parallel_deepcopy(allfdistribu_half_t, allfdistribu); @@ -65,11 +80,15 @@ DSpanSpXYVxVy PredCorr::operator()( double const final_time = iter * dt; m_poisson_solver(electrostatic_potential, electric_field_x, electric_field_y, allfdistribu); + + //copies necessary to PDI + ddc::parallel_deepcopy(allfdistribu_host, allfdistribu); + ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); ddc::PdiEvent("last_iteration") .with("iter", iter) .and_with("time_saved", final_time) - .and_with("fdistribu", allfdistribu) - .and_with("electrostatic_potential", electrostatic_potential); + .and_with("fdistribu", allfdistribu_host) + .and_with("electrostatic_potential", electrostatic_potential_host); return allfdistribu; } diff --git a/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp b/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp index fd0edbbe6..119bcf53f 100644 --- a/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp +++ b/src/geometryXYVxVy/vlasov/splitvlasovsolver.cpp @@ -17,21 +17,11 @@ SplitVlasovSolver::SplitVlasovSolver( } DSpanSpXYVxVy SplitVlasovSolver::operator()( - DSpanSpXYVxVy const allfdistribu_host, - DViewXY const electric_field_x_host, - DViewXY const electric_field_y_host, + DSpanSpXYVxVy const allfdistribu, + DViewXY const electric_field_x, + DViewXY const electric_field_y, double const dt) const { - auto allfdistribu_alloc - = ddc::create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), allfdistribu_host); - ddc::ChunkSpan allfdistribu = allfdistribu_alloc.span_view(); - auto electric_field_x_alloc = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), electric_field_x_host); - ddc::ChunkSpan electric_field_x = electric_field_x_alloc.span_view(); - auto electric_field_y_alloc = ddc:: - create_mirror_view_and_copy(Kokkos::DefaultExecutionSpace(), electric_field_y_host); - ddc::ChunkSpan electric_field_y = electric_field_y_alloc.span_view(); - m_advec_x(allfdistribu, dt / 2); m_advec_y(allfdistribu, dt / 2); m_advec_vx(allfdistribu, electric_field_x, dt / 2); @@ -40,6 +30,5 @@ DSpanSpXYVxVy SplitVlasovSolver::operator()( m_advec_y(allfdistribu, dt / 2); m_advec_x(allfdistribu, dt / 2); - ddc::parallel_deepcopy(allfdistribu_host, allfdistribu); - return allfdistribu_host; + return allfdistribu; } diff --git a/tests/geometryXYVxVy/fft_poisson.cpp b/tests/geometryXYVxVy/fft_poisson.cpp index 42f7efc7f..9013fcef3 100644 --- a/tests/geometryXYVxVy/fft_poisson.cpp +++ b/tests/geometryXYVxVy/fft_poisson.cpp @@ -12,7 +12,9 @@ #include "quadrature.hpp" #include "species_info.hpp" -TEST(FftPoissonSolver, CosineSource) +namespace { + +static void TestFftPoissonSolverCosineSource() { CoordX const x_min(0.0); CoordX const x_max(2.0 * M_PI); @@ -38,8 +40,6 @@ TEST(FftPoissonSolver, CosineSource) // Creating mesh & supports ddc::init_discrete_space(vx_min, vx_max, vx_size); ddc::init_discrete_space(vy_min, vy_max, vy_size); - ddc::init_discrete_space(vx_min, vx_max, vx_size); - ddc::init_discrete_space(vy_min, vy_max, vy_size); ddc::init_discrete_space(SplineInterpPointsVx::get_sampling()); ddc::init_discrete_space(SplineInterpPointsVy::get_sampling()); @@ -59,26 +59,21 @@ TEST(FftPoissonSolver, CosineSource) IDomainVxVy gridvxvy(gridvx, gridvy); - ddc::init_fourier_space(gridx); - ddc::init_fourier_space(gridy); - - SplineVxVyBuilder const builder_vx_vy(gridvxvy); - SplineVxVyEvaluator const evaluator_vx_vy( - g_null_boundary_2d, - g_null_boundary_2d, - g_null_boundary_2d, - g_null_boundary_2d); + IDomainXYVxVy const meshxyvxvy(gridxy, gridvxvy); SplineVxBuilder_1d const builder_vx_1d(gridvx); SplineVyBuilder_1d const builder_vy_1d(gridvy); - IDomainSpXYVxVy const mesh(gridsp, gridxy, gridvxvy); + IDomainSpXYVxVy const mesh(gridsp, meshxyvxvy); + + ddc::init_fourier_space(gridx); + ddc::init_fourier_space(gridy); // Initialise infomation about species - FieldSp charges(dom_sp); + host_t> charges(dom_sp); charges(my_ielec) = -1; charges(my_iion) = 1; - DFieldSp masses(dom_sp); + host_t masses(dom_sp); ddc::parallel_fill(masses, 1); ddc::init_discrete_space(std::move(charges), std::move(masses)); @@ -91,41 +86,72 @@ TEST(FftPoissonSolver, CosineSource) ChargeDensityCalculator rhs(quadrature_coeffs); FftPoissonSolver poisson(rhs); - DFieldXY electrostatic_potential(gridxy); - DFieldXY electric_field_x(gridxy); - DFieldXY electric_field_y(gridxy); - DFieldSpXYVxVy allfdistribu(mesh); + DFieldXY electrostatic_potential_alloc(gridxy); + DFieldXY electric_field_x_alloc(gridxy); + DFieldXY electric_field_y_alloc(gridxy); + DFieldSpXYVxVy allfdistribu_alloc(mesh); + + DSpanXY electrostatic_potential = electrostatic_potential_alloc.span_view(); + DSpanXY electric_field_x = electric_field_x_alloc.span_view(); + DSpanXY electric_field_y = electric_field_y_alloc.span_view(); + DSpanSpXYVxVy allfdistribu = allfdistribu_alloc.span_view(); // Initialization of the distribution function --> fill values - auto c_dom = ddc::remove_dims_of(mesh, gridxy); - ddc::for_each(gridxy, [&](IndexXY const ixy) { - IndexX ix = ddc::select(ixy); - IndexY iy = ddc::select(ixy); - double x = ddc::coordinate(ix); - double y = ddc::coordinate(iy); - double fdistribu_val = cos(x) + cos(y); - ddc::for_each(c_dom, [&](auto const ispvxvy) { - allfdistribu(ixy, ispvxvy) = fdistribu_val; - }); - }); + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + mesh, + KOKKOS_LAMBDA(IndexSpXYVxVy const ispxyvxvy) { + IndexX ix = ddc::select(ispxyvxvy); + IndexY iy = ddc::select(ispxyvxvy); + double x = ddc::coordinate(ix); + double y = ddc::coordinate(iy); + double fdistribu_val = Kokkos::cos(x) + Kokkos::cos(y); + allfdistribu(ispxyvxvy) = fdistribu_val; + }); poisson(electrostatic_potential, electric_field_x, electric_field_y, allfdistribu); - double error_pot = 0.0; - double error_field_x = 0.0; - double error_field_y = 0.0; - - ddc::for_each(gridxy, [&](IndexXY const ixy) { - IndexX ix = ddc::select(ixy); - IndexY iy = ddc::select(ixy); - double const exact_pot = cos(ddc::coordinate(ix)) + cos(ddc::coordinate(iy)); - error_pot = fmax(fabs(electrostatic_potential(ixy) - exact_pot), error_pot); - double const exact_field_x = sin(ddc::coordinate(ix)); - double const exact_field_y = sin(ddc::coordinate(iy)); - error_field_x = fmax(fabs(electric_field_x(ixy) - exact_field_x), error_field_x); - error_field_y = fmax(fabs(electric_field_y(ixy) - exact_field_y), error_field_y); - }); + double error_pot = ddc::parallel_transform_reduce( + Kokkos::DefaultExecutionSpace(), + gridxy, + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(IndexXY const ixy) { + IndexX ix = ddc::select(ixy); + IndexY iy = ddc::select(ixy); + double const exact_pot + = Kokkos::cos(ddc::coordinate(ix)) + Kokkos::cos(ddc::coordinate(iy)); + return Kokkos::abs(electrostatic_potential(ixy) - exact_pot); + }); + double error_field_x = ddc::parallel_transform_reduce( + Kokkos::DefaultExecutionSpace(), + gridxy, + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(IndexXY const ixy) { + IndexX ix = ddc::select(ixy); + double const exact_field_x = Kokkos::sin(ddc::coordinate(ix)); + return Kokkos::abs(electric_field_x(ixy) - exact_field_x); + }); + double error_field_y = ddc::parallel_transform_reduce( + Kokkos::DefaultExecutionSpace(), + gridxy, + 0., + ddc::reducer::max(), + KOKKOS_LAMBDA(IndexXY const ixy) { + IndexY iy = ddc::select(ixy); + double const exact_field_y = Kokkos::sin(ddc::coordinate(iy)); + return Kokkos::abs(electric_field_y(ixy) - exact_field_y); + }); + EXPECT_LE(error_pot, 1e-12); EXPECT_LE(error_field_x, 1e-10); EXPECT_LE(error_field_y, 1e-10); } + +} // namespace + +TEST(FftPoissonSolver, CosineSource) +{ + TestFftPoissonSolverCosineSource(); +} diff --git a/tests/geometryXYVxVy/quadrature.cpp b/tests/geometryXYVxVy/quadrature.cpp index 77ec62202..9c44c05b0 100644 --- a/tests/geometryXYVxVy/quadrature.cpp +++ b/tests/geometryXYVxVy/quadrature.cpp @@ -1,7 +1,5 @@ // SPDX-License-Identifier: MIT -#include - #include #include @@ -31,10 +29,10 @@ TEST(QuadratureTest, ExactForConstantFunc) IDomainXY const gridxy(gridx, gridy); - DFieldXY const quadrature_coeffs = trapezoid_quadrature_coefficients(gridxy); + host_t const quadrature_coeffs = trapezoid_quadrature_coefficients(gridxy); Quadrature const integrate(quadrature_coeffs); - DFieldXY values(gridxy); + host_t values(gridxy); ddc::for_each(gridxy, [&](ddc::DiscreteElement const idx) { values(idx) = 1.0; }); double integral = integrate(values); @@ -59,12 +57,16 @@ double compute_error(int n_elems) { using DimX = X; using DimY = Y; - using BSplinesX = UniformBSplines; - using BSplinesY = UniformBSplines; - using GrevillePointsX - = GrevilleInterpolationPoints; - using GrevillePointsY - = GrevilleInterpolationPoints; + using BSplinesX = ddc::UniformBSplines; + using BSplinesY = ddc::UniformBSplines; + using GrevillePointsX = ddc::GrevilleInterpolationPoints< + BSplinesX, + ddc::BoundCond::GREVILLE, + ddc::BoundCond::GREVILLE>; + using GrevillePointsY = ddc::GrevilleInterpolationPoints< + BSplinesY, + ddc::BoundCond::GREVILLE, + ddc::BoundCond::GREVILLE>; using IDimX = typename GrevillePointsX::interpolation_mesh_type; using IDimY = typename GrevillePointsY::interpolation_mesh_type; using IDomainX = ddc::DiscreteDomain; @@ -87,10 +89,10 @@ double compute_error(int n_elems) IDomainY const gridy(GrevillePointsY::get_domain()); IDomainXY const gridxy(gridx, gridy); - DFieldXY const quadrature_coeffs = trapezoid_quadrature_coefficients(gridxy); + host_t const quadrature_coeffs = trapezoid_quadrature_coefficients(gridxy); Quadrature const integrate(quadrature_coeffs); - DFieldXY values(gridxy); + host_t values(gridxy); ddc::for_each(gridxy, [&](ddc::DiscreteElement const idx) { double const y_cos = cos(ddc::get(ddc::coordinate(idx)));