From fbf35a3999c08d8a42587b8688a5ce42c2bb0eb6 Mon Sep 17 00:00:00 2001 From: Emily Bourne Date: Fri, 7 Jun 2024 17:06:21 +0000 Subject: [PATCH] Decouple ChargeDensityCalculator and FFT PoissonSolver Part 1 of #218 See merge request gysela-developpers/gyselalibxx!476 -------------------------------------------- --- .../bump_on_tail/bumpontail_fft.cpp | 10 +- simulations/geometryXVx/landau/landau_fft.cpp | 10 +- simulations/geometryXVx/neutrals/neutrals.cpp | 8 +- simulations/geometryXVx/sheath/sheath.cpp | 9 +- .../geometryXYVxVy/landau/landau4d_fft.cpp | 13 +- src/CMakeLists.txt | 1 + src/README.md | 1 + src/geometryXVx/geometry/README.md | 1 - src/geometryXVx/geometry/geometry.hpp | 9 - src/geometryXVx/poisson/CMakeLists.txt | 5 +- src/geometryXVx/poisson/fftqnsolver.cpp | 94 ----- src/geometryXVx/poisson/qnsolver.cpp | 36 ++ .../poisson/{fftqnsolver.hpp => qnsolver.hpp} | 14 +- src/geometryXYVxVy/geometry/README.md | 1 - src/geometryXYVxVy/geometry/geometry.hpp | 21 -- src/geometryXYVxVy/poisson/CMakeLists.txt | 3 +- src/geometryXYVxVy/poisson/fftqnsolver.cpp | 119 ------- src/geometryXYVxVy/poisson/qnsolver.cpp | 46 +++ .../poisson/{fftqnsolver.hpp => qnsolver.hpp} | 16 +- src/pde_solvers/CMakeLists.txt | 21 ++ src/pde_solvers/README.md | 17 + src/pde_solvers/fft_poisson_solver.hpp | 329 ++++++++++++++++++ src/pde_solvers/ipoisson_solver.hpp | 90 +++++ src/utils/ddc_helper.hpp | 9 + src/utils/vector_field_span.hpp | 17 +- tests/CMakeLists.txt | 1 + tests/geometryXVx/CMakeLists.txt | 1 - tests/geometryXVx/fftpoissonsolver.cpp | 193 ---------- .../geometryMX/predcorr_hybrid.cpp | 9 +- tests/geometryXYVxVy/CMakeLists.txt | 3 - tests/geometryXYVxVy/fft_poisson.cpp | 155 --------- tests/pde_solvers/CMakeLists.txt | 21 ++ tests/pde_solvers/fftpoissonsolver.cpp | 244 +++++++++++++ 33 files changed, 883 insertions(+), 644 deletions(-) delete mode 100644 src/geometryXVx/poisson/fftqnsolver.cpp create mode 100644 src/geometryXVx/poisson/qnsolver.cpp rename src/geometryXVx/poisson/{fftqnsolver.hpp => qnsolver.hpp} (71%) delete mode 100644 src/geometryXYVxVy/poisson/fftqnsolver.cpp create mode 100644 src/geometryXYVxVy/poisson/qnsolver.cpp rename src/geometryXYVxVy/poisson/{fftqnsolver.hpp => qnsolver.hpp} (72%) create mode 100644 src/pde_solvers/CMakeLists.txt create mode 100644 src/pde_solvers/README.md create mode 100644 src/pde_solvers/fft_poisson_solver.hpp create mode 100644 src/pde_solvers/ipoisson_solver.hpp delete mode 100644 tests/geometryXVx/fftpoissonsolver.cpp delete mode 100644 tests/geometryXYVxVy/fft_poisson.cpp create mode 100644 tests/pde_solvers/CMakeLists.txt create mode 100644 tests/pde_solvers/fftpoissonsolver.cpp diff --git a/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp b/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp index a4462fb60..e17645ec6 100644 --- a/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp +++ b/simulations/geometryXVx/bump_on_tail/bumpontail_fft.cpp @@ -18,13 +18,14 @@ #include "bsl_advection_x.hpp" #include "bumpontailequilibrium.hpp" #include "chargedensitycalculator.hpp" -#include "fftqnsolver.hpp" +#include "fft_poisson_solver.hpp" #include "geometry.hpp" #include "neumann_spline_quadrature.hpp" #include "paraconfpp.hpp" #include "params.yaml.hpp" #include "pdi_out.yml.hpp" #include "predcorr.hpp" +#include "qnsolver.hpp" #include "restartinitialization.hpp" #include "singlemodeperturbinitialization.hpp" #include "species_info.hpp" @@ -209,16 +210,15 @@ int main(int argc, char** argv) SplitVlasovSolver const vlasov(advection_x, advection_vx); - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXVx))); - host_t const quadrature_coeffs_host = neumann_spline_quadrature_coefficients(gridvx, builder_vx_poisson); auto const quadrature_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), quadrature_coeffs_host.span_view()); + FFTPoissonSolver fft_poisson_solver( + interpolation_domain_x); ChargeDensityCalculator rhs(quadrature_coeffs); - FftQNSolver const poisson(rhs); + QNSolver const poisson(fft_poisson_solver, rhs); PredCorr const predcorr(vlasov, poisson); diff --git a/simulations/geometryXVx/landau/landau_fft.cpp b/simulations/geometryXVx/landau/landau_fft.cpp index b3bb3c950..2c060c46e 100644 --- a/simulations/geometryXVx/landau/landau_fft.cpp +++ b/simulations/geometryXVx/landau/landau_fft.cpp @@ -17,7 +17,7 @@ #include "bsl_advection_vx.hpp" #include "bsl_advection_x.hpp" #include "chargedensitycalculator.hpp" -#include "fftqnsolver.hpp" +#include "fft_poisson_solver.hpp" #include "geometry.hpp" #include "maxwellianequilibrium.hpp" #include "neumann_spline_quadrature.hpp" @@ -25,6 +25,7 @@ #include "params.yaml.hpp" #include "pdi_out.yml.hpp" #include "predcorr.hpp" +#include "qnsolver.hpp" #include "restartinitialization.hpp" #include "singlemodeperturbinitialization.hpp" #include "species_info.hpp" @@ -209,9 +210,6 @@ int main(int argc, char** argv) SplitVlasovSolver const vlasov(advection_x, advection_vx); - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXVx))); - host_t const quadrature_coeffs_host = neumann_spline_quadrature_coefficients(gridvx, builder_vx_poisson); @@ -219,7 +217,9 @@ int main(int argc, char** argv) Kokkos::DefaultExecutionSpace(), quadrature_coeffs_host.span_view()); ChargeDensityCalculator rhs(quadrature_coeffs); - FftQNSolver const poisson(rhs); + FFTPoissonSolver fft_poisson_solver( + interpolation_domain_x); + QNSolver const poisson(fft_poisson_solver, rhs); PredCorr const predcorr(vlasov, poisson); diff --git a/simulations/geometryXVx/neutrals/neutrals.cpp b/simulations/geometryXVx/neutrals/neutrals.cpp index 02eb3e176..aca03efc6 100644 --- a/simulations/geometryXVx/neutrals/neutrals.cpp +++ b/simulations/geometryXVx/neutrals/neutrals.cpp @@ -31,7 +31,7 @@ #endif #include "Lagrange_interpolator.hpp" #include "chargedensitycalculator.hpp" -#include "fftqnsolver.hpp" +#include "fft_poisson_solver.hpp" #include "geometry.hpp" #include "irighthandside.hpp" #include "kinetic_source.hpp" @@ -43,6 +43,7 @@ #include "paraconfpp.hpp" #include "pdi_out_neutrals.yml.hpp" #include "predcorr_hybrid.hpp" +#include "qnsolver.hpp" #include "restartinitialization.hpp" #include "singlemodeperturbinitialization.hpp" #include "species_info.hpp" @@ -375,9 +376,8 @@ int main(int argc, char** argv) quadrature_coeffs_host.span_view()); ChargeDensityCalculator rhs(quadrature_coeffs); #ifdef PERIODIC_RDIMX - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXVx))); - FftQNSolver const poisson(rhs); + FFTPoissonSolver fft_poisson_solver(gridx); + QNSolver const poisson(fft_poisson_solver, rhs); #else FemNonPeriodicQNSolver const poisson(builder_x_poisson, spline_x_evaluator_poisson, rhs); #endif diff --git a/simulations/geometryXVx/sheath/sheath.cpp b/simulations/geometryXVx/sheath/sheath.cpp index 3c5604379..d57568f6d 100644 --- a/simulations/geometryXVx/sheath/sheath.cpp +++ b/simulations/geometryXVx/sheath/sheath.cpp @@ -28,7 +28,7 @@ #endif #include "Lagrange_interpolator.hpp" #include "chargedensitycalculator.hpp" -#include "fftqnsolver.hpp" +#include "fft_poisson_solver.hpp" #include "geometry.hpp" #include "irighthandside.hpp" #include "kinetic_source.hpp" @@ -39,6 +39,7 @@ #include "paraconfpp.hpp" #include "pdi_out.yml.hpp" #include "predcorr.hpp" +#include "qnsolver.hpp" #include "restartinitialization.hpp" #include "sheath.yaml.hpp" #include "singlemodeperturbinitialization.hpp" @@ -309,9 +310,9 @@ int main(int argc, char** argv) quadrature_coeffs_host.span_view()); ChargeDensityCalculator rhs(quadrature_coeffs); #ifdef PERIODIC_RDIMX - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXVx))); - FftQNSolver const poisson(rhs); + FFTPoissonSolver fft_poisson_solver( + interpolation_domain_x); + QNSolver const poisson(fft_poisson_solver, rhs); #else FemNonPeriodicQNSolver const poisson(builder_x_poisson, spline_x_evaluator_poisson, rhs); #endif diff --git a/simulations/geometryXYVxVy/landau/landau4d_fft.cpp b/simulations/geometryXYVxVy/landau/landau4d_fft.cpp index c480ee429..e1fcc6917 100644 --- a/simulations/geometryXYVxVy/landau/landau4d_fft.cpp +++ b/simulations/geometryXYVxVy/landau/landau4d_fft.cpp @@ -15,7 +15,7 @@ #include "bsl_advection_vx.hpp" #include "bsl_advection_x.hpp" -#include "fftqnsolver.hpp" +#include "fft_poisson_solver.hpp" #include "geometry.hpp" #include "maxwellianequilibrium.hpp" #include "neumann_spline_quadrature.hpp" @@ -23,6 +23,7 @@ #include "params.yaml.hpp" #include "pdi_out.yml.hpp" #include "predcorr.hpp" +#include "qnsolver.hpp" #include "singlemodeperturbinitialization.hpp" //#include "species_info.hpp" #include "spline_interpolator.hpp" @@ -92,6 +93,7 @@ int main(int argc, char** argv) IDomainX interpolation_domain_x(SplineInterpPointsX::get_domain()); IDomainY interpolation_domain_y(SplineInterpPointsY::get_domain()); + IDomainXY interpolation_domain_xy(interpolation_domain_x, interpolation_domain_y); IDomainVx interpolation_domain_vx(SplineInterpPointsVx::get_domain()); IDomainVy interpolation_domain_vy(SplineInterpPointsVy::get_domain()); IDomainVxVy interpolation_domain_vxvy(interpolation_domain_vx, interpolation_domain_vy); @@ -209,11 +211,6 @@ int main(int argc, char** argv) SplitVlasovSolver const vlasov(advection_x, advection_y, advection_vx, advection_vy); - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXYVxVy))); - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXYVxVy))); - host_t const quadrature_coeffs_host = neumann_spline_quadrature_coefficients( interpolation_domain_vxvy, builder_vx_1d, @@ -221,8 +218,10 @@ int main(int argc, char** argv) auto quadrature_coeffs = ddc::create_mirror_view_and_copy( Kokkos::DefaultExecutionSpace(), quadrature_coeffs_host.span_view()); + FFTPoissonSolver fft_poisson_solver( + interpolation_domain_xy); ChargeDensityCalculator const rhs(quadrature_coeffs); - FftQNSolver const poisson(rhs); + QNSolver const poisson(fft_poisson_solver, rhs); // Create predcorr operator PredCorr const predcorr(vlasov, poisson); diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 96ebe054a..8f7ef1602 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -13,3 +13,4 @@ add_subdirectory(geometryXYVxVy) add_subdirectory(interpolation) add_subdirectory(advection) add_subdirectory(timestepper) +add_subdirectory(pde_solvers) diff --git a/src/README.md b/src/README.md index a40ce24dc..b6a8a218d 100644 --- a/src/README.md +++ b/src/README.md @@ -9,6 +9,7 @@ The `src/` folder contains all the code necessary to build a gyrokinetic semi-La - [geometry5D](./geometry5D/README.md) - Code describing methods which are specific to a simulation with 3 spatial dimensions and 2 velocity dimension. This will be the general geometry used for GyselaX++ development. - [interpolation](./interpolation/README.md) - Code describing interpolation methods. +- [pde\_solvers](./pde_solvers/README.md) - Code describing different methods for solving PDEs (e.g. Poisson's equation). - [quadrature](./quadrature/README.md) - Code describing different quadrature methods. - [speciesinfo](./speciesinfo/README.md) - Code used to describe the different species. - [timestepper](./timestepper/README.md) - Code used to describe time-stepping methods (e.g. Runge-Kutta methods). diff --git a/src/geometryXVx/geometry/README.md b/src/geometryXVx/geometry/README.md index b49d82ef4..6b5eae087 100644 --- a/src/geometryXVx/geometry/README.md +++ b/src/geometryXVx/geometry/README.md @@ -18,6 +18,5 @@ The shortcuts defined in the geometry file represent: 13. The type of a span of doubles defined on each of the domains (e.g. DSpanX). 12. The templated type of a field view defined on each of the domains (e.g. ViewX). A view is very similar to a span. The only difference is that a view is constant. In other words, the data in the underlying field cannot be modified using a view. 13. The type of a view of doubles defined on each of the domains (e.g. DViewX). -14. A dimension in real space representing the Fourier mode of the X-dimension (RDimFx). 15. Types representing coordinates, and the grid points as well as their indices, distances and domains for the Fourier mode. 16. A class GeometryXVx detailing some of the above types in a generic way which allows them to be accessed from a context where the final geometry selected is unknown. diff --git a/src/geometryXVx/geometry/geometry.hpp b/src/geometryXVx/geometry/geometry.hpp index 18c60d9d6..f802fe516 100644 --- a/src/geometryXVx/geometry/geometry.hpp +++ b/src/geometryXVx/geometry/geometry.hpp @@ -427,15 +427,6 @@ using DViewSpVx = ViewSpVx; using DViewSpXVx = ViewSpXVx; -using RDimFx = ddc::Fourier; -using CoordFx = ddc::Coordinate; -struct IDimFx : ddc::PeriodicSampling -{ -}; -using IndexFx = ddc::DiscreteElement; -using IVectFx = ddc::DiscreteVector; -using IDomainFx = ddc::DiscreteDomain; - /** * @brief A class providing aliases for useful subdomains of the geometry. It is used as template parameter for generic dimensionality-agnostic operators such as advections. */ diff --git a/src/geometryXVx/poisson/CMakeLists.txt b/src/geometryXVx/poisson/CMakeLists.txt index 9165f6d9f..2f68ec4b1 100644 --- a/src/geometryXVx/poisson/CMakeLists.txt +++ b/src/geometryXVx/poisson/CMakeLists.txt @@ -19,8 +19,9 @@ target_include_directories("poisson_${GEOMETRY_VARIANT}" target_link_libraries("poisson_${GEOMETRY_VARIANT}" PUBLIC DDC::DDC - gslx::quadrature gslx::geometry_${GEOMETRY_VARIANT} + gslx::pde_solvers + gslx::quadrature gslx::speciesinfo ) @@ -30,7 +31,7 @@ endforeach() target_sources(poisson_xperiod_vx PRIVATE - fftqnsolver.cpp + qnsolver.cpp femperiodicqnsolver.cpp ) diff --git a/src/geometryXVx/poisson/fftqnsolver.cpp b/src/geometryXVx/poisson/fftqnsolver.cpp deleted file mode 100644 index 94b3cb805..000000000 --- a/src/geometryXVx/poisson/fftqnsolver.cpp +++ /dev/null @@ -1,94 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include -#include -#include -#include - -#include -#include - -#include - -#include "fftqnsolver.hpp" - -FftQNSolver::FftQNSolver(IChargeDensityCalculator const& compute_rho) : m_compute_rho(compute_rho) -{ -} - -// 1- Inner solvers shall be passed in the constructor -// 2- Should it take an array of distribution functions ? -void FftQNSolver::operator()( - DSpanX const electrostatic_potential, - DSpanX const electric_field, - DViewSpXVx const allfdistribu) const -{ - Kokkos::Profiling::pushRegion("QNSolver"); - assert(electrostatic_potential.domain() == ddc::get_domain(allfdistribu)); - IDomainX const x_dom = electrostatic_potential.domain(); - // Compute the RHS of the Quasi-Neutrality equation. - DFieldX rho(x_dom); - - m_compute_rho(rho, allfdistribu); - - // Build a mesh in the fourier space, for N points - IDomainFx const k_mesh = ddc::FourierMesh(x_dom, false); - device_t, IDomainFx>> 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::DefaultExecutionSpace(), - intermediate_chunk, - rho.span_view(), - ddc::kwArgs_fft {norm}); - - // Solve Poisson's equation -d2Phi/dx2 = rho - // in Fourier space as kx*kx*FFT(Phi)=FFT(rho) - - // First, 0 mode of Phi is set to 0 to avoid divergency. Note: to allow writing in the GPU memory, starting a device kernel is necessary which is performed by iterating on a 0D domain (single element). - ddc::parallel_for_each( - Kokkos::DefaultExecutionSpace(), - ddc::DiscreteDomain<>(), - KOKKOS_LAMBDA(ddc::DiscreteElement<>) { - intermediate_chunk(k_mesh.front()) = Kokkos::complex(0.); - }); - - ddc::parallel_for_each( - Kokkos::DefaultExecutionSpace(), - k_mesh.remove_first(IVectFx(1)), - KOKKOS_LAMBDA(IndexFx const ikx) { - intermediate_chunk(ikx) - = intermediate_chunk(ikx) / (coordinate(ikx) * coordinate(ikx)); - }); - - // Find the electric field in Fourier space - // FFT(efield) = -kx*i*FFT(Phi) - Kokkos::Profiling::pushRegion("ElectricField"); - Kokkos::complex imaginary_unit(0.0, 1.0); - device_t, IDomainFx>> fourier_efield_alloc(k_mesh); - ddc::ChunkSpan fourier_efield = fourier_efield_alloc.span_view(); - - ddc::parallel_for_each( - Kokkos::DefaultExecutionSpace(), - k_mesh, - KOKKOS_LAMBDA(IndexFx const ikx) { - fourier_efield(ikx) = -imaginary_unit * coordinate(ikx) * intermediate_chunk(ikx); - }); - - // Perform the inverse 1D FFT of the solution to deduce the electric field - ddc:: - ifft(Kokkos::DefaultExecutionSpace(), - electric_field, - fourier_efield, - ddc::kwArgs_fft {norm}); - Kokkos::Profiling::popRegion(); - - // Perform the inverse 1D FFT of the solution to deduce the electrostatic potential - ddc:: - ifft(Kokkos::DefaultExecutionSpace(), - electrostatic_potential.span_view(), - intermediate_chunk, - ddc::kwArgs_fft {norm}); - Kokkos::Profiling::popRegion(); -} diff --git a/src/geometryXVx/poisson/qnsolver.cpp b/src/geometryXVx/poisson/qnsolver.cpp new file mode 100644 index 000000000..8335b8cc9 --- /dev/null +++ b/src/geometryXVx/poisson/qnsolver.cpp @@ -0,0 +1,36 @@ +// SPDX-License-Identifier: MIT + +#include +#include +#include +#include + +#include + +#include + +#include "qnsolver.hpp" + +QNSolver::QNSolver(PoissonSolver const& solve_poisson, IChargeDensityCalculator const& compute_rho) + : m_solve_poisson(solve_poisson) + , m_compute_rho(compute_rho) +{ +} + +void QNSolver::operator()( + DSpanX const electrostatic_potential, + DSpanX const electric_field, + DViewSpXVx const allfdistribu) const +{ + Kokkos::Profiling::pushRegion("QNSolver"); + assert(electrostatic_potential.domain() == ddc::get_domain(allfdistribu)); + IDomainX const x_dom = electrostatic_potential.domain(); + // Compute the RHS of the Quasi-Neutrality equation. + DFieldX rho(x_dom); + + m_compute_rho(rho, allfdistribu); + + m_solve_poisson(electrostatic_potential, electric_field, rho); + + Kokkos::Profiling::popRegion(); +} diff --git a/src/geometryXVx/poisson/fftqnsolver.hpp b/src/geometryXVx/poisson/qnsolver.hpp similarity index 71% rename from src/geometryXVx/poisson/fftqnsolver.hpp rename to src/geometryXVx/poisson/qnsolver.hpp index 23e9ad83e..613247d68 100644 --- a/src/geometryXVx/poisson/fftqnsolver.hpp +++ b/src/geometryXVx/poisson/qnsolver.hpp @@ -3,6 +3,7 @@ #pragma once #include "ichargedensitycalculator.hpp" +#include "ipoisson_solver.hpp" #include "iqnsolver.hpp" /** @@ -17,19 +18,26 @@ * The electric field, @f$ \frac{d \phi}{dx} @f$ is calculated using * a spline interpolation implemented in ElectricField. */ -class FftQNSolver : public IQNSolver +class QNSolver : public IQNSolver { + using PoissonSolver = IPoissonSolver< + IDomainX, + IDomainX, + std::experimental::layout_right, + typename Kokkos::DefaultExecutionSpace::memory_space>; + PoissonSolver const& m_solve_poisson; IChargeDensityCalculator const& m_compute_rho; public: /** * Construct the FftQNSolver operator. * + * @param solve_poisson The operator which solves the Poisson solver. * @param compute_rho The operator which calculates the charge density, the right hand side of the equation. */ - FftQNSolver(IChargeDensityCalculator const& compute_rho); + QNSolver(PoissonSolver const& solve_poisson, IChargeDensityCalculator const& compute_rho); - ~FftQNSolver() override = default; + ~QNSolver() override = default; /** * The operator which solves the equation using the method described by the class. diff --git a/src/geometryXYVxVy/geometry/README.md b/src/geometryXYVxVy/geometry/README.md index b627a7f89..068c990ba 100644 --- a/src/geometryXYVxVy/geometry/README.md +++ b/src/geometryXYVxVy/geometry/README.md @@ -18,6 +18,5 @@ The shortcuts defined in the geometry file represent: 13. The type of a span of doubles defined on each of the domains (e.g. `DSpanX`). 12. The templated type of a field view defined on each of the domains (e.g. `ViewX`). A view is very similar to a span. The only difference is that a view is constant. In other words, the data in the underlying field cannot be modified using a view. 13. The type of a view of doubles defined on each of the domains (e.g. `DViewX`). -14. A dimension in real space representing the Fourier mode of the spatial dimensions (`RDimFx`, `RDimFy`). 15. Types representing coordinates, and the grid points as well as their indices, distances and domains for the Fourier modes. 16. A class GeometryXVx detailing some of the above types in a generic way which allows them to be accessed from a context where the final geometry selected is unknown. diff --git a/src/geometryXYVxVy/geometry/geometry.hpp b/src/geometryXYVxVy/geometry/geometry.hpp index 9d98b8857..f05524e30 100644 --- a/src/geometryXYVxVy/geometry/geometry.hpp +++ b/src/geometryXYVxVy/geometry/geometry.hpp @@ -427,27 +427,6 @@ template using ViewSpXYVxVy = device_t>; using DViewSpXYVxVy = ViewSpXYVxVy; -// For Fourier -using RDimFx = ddc::Fourier; -using RDimFy = ddc::Fourier; -using CoordFx = ddc::Coordinate; -using CoordFy = ddc::Coordinate; -struct IDimFx : ddc::PeriodicSampling -{ -}; -struct IDimFy : ddc::PeriodicSampling -{ -}; -using IndexFx = ddc::DiscreteElement; -using IndexFy = ddc::DiscreteElement; -using IndexFxFy = ddc::DiscreteElement; -using IVectFx = ddc::DiscreteVector; -using IVectFy = ddc::DiscreteVector; -using IVectFxFy = ddc::DiscreteVector; -using IDomainFx = ddc::DiscreteDomain; -using IDomainFy = ddc::DiscreteDomain; -using IDomainFxFy = ddc::DiscreteDomain; - /** * @brief A class providing aliases for useful subdomains of the geometry. It is used as template parameter for generic dimensionality-agnostic operat ors such as advections. diff --git a/src/geometryXYVxVy/poisson/CMakeLists.txt b/src/geometryXYVxVy/poisson/CMakeLists.txt index 0f021f12d..4134f41d3 100644 --- a/src/geometryXYVxVy/poisson/CMakeLists.txt +++ b/src/geometryXYVxVy/poisson/CMakeLists.txt @@ -3,7 +3,7 @@ add_library("poisson_xy" STATIC chargedensitycalculator.cpp nullqnsolver.cpp - fftqnsolver.cpp + qnsolver.cpp ) target_compile_features("poisson_xy" @@ -20,6 +20,7 @@ target_link_libraries("poisson_xy" DDC::DDC sll::SLL gslx::geometry_xyvxvy + gslx::pde_solvers gslx::speciesinfo gslx::utils ) diff --git a/src/geometryXYVxVy/poisson/fftqnsolver.cpp b/src/geometryXYVxVy/poisson/fftqnsolver.cpp deleted file mode 100644 index 200a97c3e..000000000 --- a/src/geometryXYVxVy/poisson/fftqnsolver.cpp +++ /dev/null @@ -1,119 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include -#include -#include -#include - -#include - -#include - -#include "fftqnsolver.hpp" - -FftQNSolver::FftQNSolver(IChargeDensityCalculator const& compute_rho) : m_compute_rho(compute_rho) -{ -} - -// 1- Inner solvers shall be passed in the constructor -// 2- Should it take an array of distribution functions ? -void FftQNSolver::operator()( - DSpanXY const electrostatic_potential, - DSpanXY const electric_field_x, - DSpanXY const electric_field_y, - DViewSpXYVxVy const allfdistribu) const -{ - Kokkos::Profiling::pushRegion("QNSolver"); - assert((electrostatic_potential.domain() == ddc::get_domain(allfdistribu))); - IDomainXY const xy_dom = electrostatic_potential.domain(); - - // Compute the RHS of the Quasi-Neutrality equation. - 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); - - 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::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)) - 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 fourier_efield_alloc - = ddc::Chunk(k_mesh, ddc::DeviceAllocator>()); - ddc::ChunkSpan fourier_efield = fourier_efield_alloc.span_view(); - - // Calculate x component - 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::DefaultExecutionSpace(), - electric_field_x.span_view(), - fourier_efield.span_view(), - ddc::kwArgs_fft {norm}); - - // Calculate y component - 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::DefaultExecutionSpace(), - electric_field_y.span_view(), - fourier_efield.span_view(), - ddc::kwArgs_fft {norm}); - Kokkos::Profiling::popRegion(); - - // Perform the inverse 1D FFT of the solution to deduce the electrostatic potential - ddc:: - ifft(Kokkos::DefaultExecutionSpace(), - electrostatic_potential.span_view(), - intermediate_chunk.span_view(), - ddc::kwArgs_fft {norm}); - - Kokkos::Profiling::popRegion(); -} diff --git a/src/geometryXYVxVy/poisson/qnsolver.cpp b/src/geometryXYVxVy/poisson/qnsolver.cpp new file mode 100644 index 000000000..539c19313 --- /dev/null +++ b/src/geometryXYVxVy/poisson/qnsolver.cpp @@ -0,0 +1,46 @@ +// SPDX-License-Identifier: MIT + +#include +#include +#include +#include + +#include + +#include +#include + +#include "qnsolver.hpp" + +QNSolver::QNSolver(PoissonSolver const& solve_poisson, IChargeDensityCalculator const& compute_rho) + : m_solve_poisson(solve_poisson) + , m_compute_rho(compute_rho) +{ +} + +void QNSolver::operator()( + DSpanXY const electrostatic_potential, + DSpanXY const electric_field_x, + DSpanXY const electric_field_y, + DViewSpXYVxVy const allfdistribu) const +{ + Kokkos::Profiling::pushRegion("QNSolver"); + assert((electrostatic_potential.domain() == ddc::get_domain(allfdistribu))); + IDomainXY const xy_dom = electrostatic_potential.domain(); + + // Compute the RHS of the Quasi-Neutrality equation. + DFieldXY rho(xy_dom); + DFieldVxVy contiguous_slice_vxvy(allfdistribu.domain()); + m_compute_rho(rho, allfdistribu); + + VectorFieldSpan< + double, + IDomainXY, + NDTag, + typename DFieldXY::layout_type, + Kokkos::DefaultExecutionSpace::memory_space> + electric_field(electric_field_x, electric_field_y); + m_solve_poisson(electrostatic_potential, electric_field, rho); + + Kokkos::Profiling::popRegion(); +} diff --git a/src/geometryXYVxVy/poisson/fftqnsolver.hpp b/src/geometryXYVxVy/poisson/qnsolver.hpp similarity index 72% rename from src/geometryXYVxVy/poisson/fftqnsolver.hpp rename to src/geometryXYVxVy/poisson/qnsolver.hpp index 8d9f3c5c4..c0acec17c 100644 --- a/src/geometryXYVxVy/poisson/fftqnsolver.hpp +++ b/src/geometryXYVxVy/poisson/qnsolver.hpp @@ -3,6 +3,7 @@ #pragma once #include "chargedensitycalculator.hpp" +#include "ipoisson_solver.hpp" #include "iqnsolver.hpp" /** @@ -17,19 +18,26 @@ * The electric field, @f$ \frac{d \phi}{dx} @f$ is calculated using * a spline interpolation implemented in ElectricField. */ -class FftQNSolver : public IQNSolver +class QNSolver : public IQNSolver { + using PoissonSolver = IPoissonSolver< + IDomainXY, + IDomainXY, + std::experimental::layout_right, + typename Kokkos::DefaultExecutionSpace::memory_space>; + PoissonSolver const& m_solve_poisson; IChargeDensityCalculator const& m_compute_rho; public: /** - * Construct the FftQNSolver operator. + * Construct the QNSolver operator. * + * @param solve_poisson The operator which solves the Poisson solver. * @param compute_rho The operator which calculates the charge density, the right hand side of the equation. */ - FftQNSolver(IChargeDensityCalculator const& compute_rho); + QNSolver(PoissonSolver const& solve_poisson, IChargeDensityCalculator const& compute_rho); - ~FftQNSolver() override = default; + ~QNSolver() override = default; /** * The operator which solves the equation using the method described by the class. diff --git a/src/pde_solvers/CMakeLists.txt b/src/pde_solvers/CMakeLists.txt new file mode 100644 index 000000000..46ed6dad5 --- /dev/null +++ b/src/pde_solvers/CMakeLists.txt @@ -0,0 +1,21 @@ +# SPDX-License-Identifier: MIT + +add_library("pde_solvers" INTERFACE) + +target_compile_features("pde_solvers" + INTERFACE + cxx_std_17 +) + +target_include_directories("pde_solvers" + INTERFACE "${CMAKE_CURRENT_SOURCE_DIR}" +) + +target_link_libraries("pde_solvers" + INTERFACE + DDC::DDC + sll::SLL + gslx::utils +) + +add_library("gslx::pde_solvers" ALIAS "pde_solvers") diff --git a/src/pde_solvers/README.md b/src/pde_solvers/README.md new file mode 100644 index 000000000..a365e03d1 --- /dev/null +++ b/src/pde_solvers/README.md @@ -0,0 +1,17 @@ +# PDE Solvers + +This folder contains all generalised code describing different methods for solving PDEs (Partial Differential Equations). Solvers exist for the following PDEs: + +- Poisson's equation : $-\Delta \phi = \rho$ + +## Poisson's equation + +The following methods exist for solving Poisson's equations: +- FFTPoissonSolver + +These solvers implement 2 interfaces: +- `chunk_span_type operator()(chunk_span_type phi, chunk_span_type rho) const` +- `chunk_span_type operator()(chunk_span_type phi, vector_span_type E, chunk_span_type rho) const` + +The second interface calculates $\phi$ the solution to the equation but also $E = - \nabla \phi$. + diff --git a/src/pde_solvers/fft_poisson_solver.hpp b/src/pde_solvers/fft_poisson_solver.hpp new file mode 100644 index 000000000..9ce42f0d9 --- /dev/null +++ b/src/pde_solvers/fft_poisson_solver.hpp @@ -0,0 +1,329 @@ +#pragma once + +#include +#include + +#include +#include + +#include "ipoisson_solver.hpp" + +template < + class LaplacianDomain, + class FullDomain, + class ExecSpace, + class LayoutSpace = std::experimental::layout_right> +class FFTPoissonSolver; + +/** + * A class solve the follwing equation: + * @f$ -\Delta \phi = \rho @f$ + * using a Fourier transform. + * + * @tparam LaplacianDomain The domain on which the equation is defined. + * @tparam FullDomain The domain on which the operator() acts. This is equal to the + * LaplacianDomain plus any batched dimensions. + * @tparam ExecSpace The space (CPU/GPU) where the calculations will take place. + * @tparam LayoutSpace The layout space of the ChunkSpans passed to operator(). + */ +template +class FFTPoissonSolver, FullDomain, ExecSpace, LayoutSpace> + : public IPoissonSolver< + ddc::DiscreteDomain, + FullDomain, + LayoutSpace, + typename ExecSpace::memory_space> +{ +private: + using base_type = IPoissonSolver< + ddc::DiscreteDomain, + FullDomain, + LayoutSpace, + typename ExecSpace::memory_space>; + +public: + template + struct IDimFourier : ddc::PeriodicSampling> + { + }; + +public: + /// @brief The ChunkSpan type of the arguments to operator(). + using chunk_span_type = typename base_type::chunk_span_type; + /// @brief The const ChunkSpan type of the arguments to operator(). + using chunk_view_type = typename base_type::chunk_view_type; + + /// @brief The type of the derivative of @f$ \phi @f$. + using vector_span_type = typename base_type::vector_span_type; + + /// @brief The DiscreteDomain describing the batch dimensions. + using batch_domain_type = typename base_type::batch_domain_type; + /// @brief The DiscreteElement for indexing a batch dimension. + using batch_element_type = typename base_type::batch_element_type; + + /// @brief The type of the domain on which the equation is defined. + using laplacian_domain_type = typename base_type::laplacian_domain_type; + + /// @brief The layout space of the ChunkSpans passed to operator(). + using layout_space = typename base_type::layout_space; + /// @brief The space (CPU/GPU) where the ChunkSpans passed to operator() are saved. + using memory_space = typename base_type::memory_space; + + /// @brief The type of the Fourier space domain. + using fourier_domain_type + = ddc::DiscreteDomain...>; + /// @brief The type of an index of the Fourier space domain. + using fourier_element_type = typename fourier_domain_type::discrete_element_type; + + /// @brief The type of a Chunk storing the Fourier transform of a function. + using fourier_chunk_type = ddc::Chunk< + Kokkos::complex, + fourier_domain_type, + ddc::KokkosAllocator, memory_space>>; + /// @brief The type of a ChunkSpan storing the Fourier transform of a function. + using fourier_span_type = typename fourier_chunk_type::span_type; + +private: + /// @brief The normalisation used for the Fourier transform + static constexpr ddc::FFT_Normalization m_norm = ddc::FFT_Normalization::BACKWARD; + +private: + /** + * @brief The multiplicative factor corresponding to the Laplace operator @f$ \Delta @f$ + * expressed in Fourier space at a given mode. + * + * @param index The index of the Fourier mode. + */ + template + KOKKOS_FUNCTION static double get_laplace_operator(ddc::DiscreteElement index) + { + return (((double)ddc::coordinate(ddc::select(index)) + * (double)ddc::coordinate(ddc::select(index))) + + ...); + } + + /** + * @brief Differentiate an expression from its representation in Fourier space by multiplying + * by -i * k and then converting back to real space. + * + * @param[out] derivative The ChunkSpan where the derivative will be saved. + * @param[out] fourier_derivative The derivative of the function in Fourier space will be saved. + * @param[in] values The ChunkSpan containing the values of the function in Fourier space. + * + * @tparam Dim The dimension along which the expression is differentiated. + */ + template + void differentiate_and_invert_fourier_values( + ddc::ChunkSpan derivative, + fourier_span_type fourier_derivative, + fourier_span_type values) const + { + negative_differentiate_equation(fourier_derivative, values); + // Perform the inverse 1D FFT of the solution to deduce the electric field + ddc:: + ifft(ExecSpace(), + derivative.span_view(), + fourier_derivative.span_view(), + ddc::kwArgs_fft {m_norm}); + } + + /** + * @brief Get the gradient of a 1D expression from its representation in Fourier space. + * + * @param[out] gradient The ChunkSpan where the derivative will be saved. + * @param[out] fourier_derivative The derivative of the function in Fourier space will be saved. + * @param[in] values The ChunkSpan containing the values of the function in Fourier space. + */ + void get_gradient( + ddc::ChunkSpan gradient, + fourier_span_type fourier_derivative, + fourier_span_type values) const + { + using RDim = + typename ddc::type_seq_element_t<0, ddc::to_type_seq_t>:: + continuous_dimension_type; + using FourierDim = IDimFourier; + differentiate_and_invert_fourier_values(gradient, fourier_derivative, values); + } + + /** + * @brief Get the gradient of a multi-dimensional expression from its representation in Fourier space. + * + * @param[out] gradient The VectorFieldSpan where the gradient will be saved. + * @param[out] fourier_derivative The derivative of the function in Fourier space will be saved. + * @param[in] values The ChunkSpan containing the values of the function in Fourier space. + */ + template + void get_gradient( + VectorFieldSpan< + double, + laplacian_domain_type, + NDTag, + layout_space, + memory_space> gradient, + fourier_span_type fourier_derivative, + fourier_span_type values) const + { + ((differentiate_and_invert_fourier_values< + IDimFourier>(ddcHelper::get(gradient), fourier_derivative, values)), + ...); + } + + template + void init_fourier_space(ddc::DiscreteDomain dom) + { + using IDimF = IDimFourier; + ddc::init_discrete_space(ddc::init_fourier_space(dom)); + } + +public: + /** + * @brief A function to solve the Poisson equation in Fourier space + * This function should be private. It is not due to the inclusion of a KOKKOS_LAMBDA + * + * @param[out] intermediate_chunk The solution to the Poisson equation in Fourier space. + * @param[in] rho The right-hand side of the Poisson equation in Fourier space. + */ + template + void solve_poisson_equation( + fourier_span_type intermediate_chunk, + ddc::ChunkSpan rho) const + { + // Compute FFT(rho) + ddc::fft(ExecSpace(), intermediate_chunk, rho, ddc::kwArgs_fft {m_norm}); + + fourier_domain_type const k_mesh = intermediate_chunk.domain(); + + // Solve Poisson's equation -\Delta phi = -(\sum_j \partial_j^2) \phi = rho + // in Fourier space as -(\sum_j i*k_i * i*k_i) FFT(Phi) = FFT(rho)) + ddc::parallel_for_each( + ExecSpace(), + k_mesh, + KOKKOS_LAMBDA(fourier_element_type const ik) { + if (ik != k_mesh.front()) { + intermediate_chunk(ik) = intermediate_chunk(ik) / get_laplace_operator(ik); + } else { + intermediate_chunk(ik) = 0.; + } + }); + } + + /** + * @brief Differentiate and multiply by -1 an expression in Fourier space by multiplying by -i * k + * This function should be private. It is not due to the inclusion of a KOKKOS_LAMBDA + * + * @param derivative The ChunkSpan where the derivative will be saved. + * @param values The ChunkSpan containing the values of the function in Fourier space. + * + * @tparam Dim The dimension along which the expression is differentiated. + */ + template + void negative_differentiate_equation(fourier_span_type derivative, fourier_span_type values) + const + { + Kokkos::complex imaginary_unit(0.0, 1.0); + ddc::parallel_for_each( + ExecSpace(), + values.domain(), + KOKKOS_LAMBDA(fourier_element_type const ik) { + ddc::DiscreteElement const ikx = ddc::select(ik); + derivative(ik) = -imaginary_unit * ddc::coordinate(ikx) * values(ik); + }); + } + +public: + /** + * @brief A constructor for the FFT Poisson solver. + * This constructor calls ddc::init_discrete_space so it should only be called once per + * simulation. + * + * @param laplacian_domain The domain on which the equation should be solved. + */ + FFTPoissonSolver(laplacian_domain_type laplacian_domain) + { + ((init_fourier_space(ddc::DiscreteDomain(laplacian_domain))), ...); + } + + /** + * @brief An operator which calculates the solution @f$\phi@f$ to Poisson's equation: + * @f$ - \Delta \phi = \rho @f$ + * + * @param[out] phi The solution to Poisson's equation. + * @param[in] rho The right-hand side of Poisson's equation. + * + * @return A reference to the solution to Poisson's equation. + */ + virtual chunk_span_type operator()(chunk_span_type phi, chunk_span_type rho) const final + { + Kokkos::Profiling::pushRegion("FFTPoissonSolver"); + + laplacian_domain_type domain(phi.domain()); + batch_domain_type batch_domain(phi.domain()); + + // Build a mesh in the fourier space, for N points + fourier_domain_type const k_mesh = ddc::FourierMesh< + IDimFourier...>(domain, false); + + fourier_chunk_type intermediate_chunk_alloc(k_mesh); + fourier_span_type intermediate_chunk = intermediate_chunk_alloc.span_view(); + + ddc::for_each(batch_domain, [&](batch_element_type ib) { + solve_poisson_equation(intermediate_chunk, rho[ib]); + + // Perform the inverse 1D FFT of the solution to deduce the electrostatic potential + ddc:: + ifft(ExecSpace(), + phi[ib], + intermediate_chunk.span_view(), + ddc::kwArgs_fft {m_norm}); + }); + + Kokkos::Profiling::popRegion(); + return phi; + } + + /** + * @brief An operator which calculates the solution @f$\phi@f$ to Poisson's equation and + * its derivative: + * @f$ - \Delta \phi = \rho @f$ + * @f$ E = - \nabla \phi @f$ + * + * @param[out] phi The solution to Poisson's equation. + * @param[out] E The derivative of the solution to Poisson's equation. + * @param[in] rho The right-hand side of Poisson's equation. + * + * @return A reference to the solution to Poisson's equation. + */ + virtual chunk_span_type operator()(chunk_span_type phi, vector_span_type E, chunk_span_type rho) + const final + { + Kokkos::Profiling::pushRegion("FFTPoissonSolver"); + + laplacian_domain_type domain(phi.domain()); + batch_domain_type batch_domain(phi.domain()); + + // Build a mesh in the fourier space, for N points + fourier_domain_type const k_mesh = ddc::FourierMesh< + IDimFourier...>(domain, false); + + fourier_chunk_type intermediate_chunk_alloc(k_mesh); + fourier_chunk_type fourier_efield_alloc(k_mesh); + + fourier_span_type intermediate_chunk = intermediate_chunk_alloc.span_view(); + fourier_span_type fourier_efield = fourier_efield_alloc.span_view(); + + ddc::for_each(batch_domain, [&](batch_element_type ib) { + solve_poisson_equation(intermediate_chunk, rho[ib]); + get_gradient(E[ib], fourier_efield, intermediate_chunk); + + // Perform the inverse 1D FFT of the solution to deduce the electrostatic potential + ddc:: + ifft(ExecSpace(), + phi[ib], + intermediate_chunk.span_view(), + ddc::kwArgs_fft {m_norm}); + }); + Kokkos::Profiling::popRegion(); + return phi; + } +}; diff --git a/src/pde_solvers/ipoisson_solver.hpp b/src/pde_solvers/ipoisson_solver.hpp new file mode 100644 index 000000000..4f55bed2b --- /dev/null +++ b/src/pde_solvers/ipoisson_solver.hpp @@ -0,0 +1,90 @@ +#pragma once + +#include + +#include + +template +class IPoissonSolver; + +/** + * An abstract class from which a Poisson solver can inherit. + * Classes inheriting from this must implement a way to solve the following equation: + * @f$ -\Delta \phi = \rho @f$ + * + * @tparam LaplacianDomain The domain on which the equation is defined. + * @tparam FullDomain The domain on which the operator() acts. This is equal to the + * LaplacianDomain plus any batched dimensions. + * @tparam LayoutSpace The layout space of the ChunkSpans passed to operator(). + * @tparam MemorySpace The space (CPU/GPU) where the ChunkSpans passed to operator() + * are saved. + */ +template +class IPoissonSolver, FullDomain, LayoutSpace, MemorySpace> +{ +protected: + /// @brief The tags describing the real dimensions in the equation. + using real_laplacian_tags = ddc::detail::TypeSeq; + /// @brief The tags describing the discrete dimensions in the equation. + using laplacian_tags = ddc::detail::TypeSeq; + /// @brief The tags describing the dimensions of the domain on which the operator acts. + using space_tags = ddc::to_type_seq_t; + /// @brief The tags describing the batched dimensions. + using batch_tags = ddc::type_seq_remove_t; + +protected: + /// @brief Indicates whether the gradient is represented by a VectorFieldSpan or a ChunkSpan. + static constexpr bool using_vector_span = ddcHelper::type_seq_length_v == 1; + +public: + /// @brief The ChunkSpan type of the arguments to operator(). + using chunk_span_type = ddc::ChunkSpan; + /// @brief The const ChunkSpan type of the arguments to operator(). + using chunk_view_type = ddc::ChunkView; + + /// @brief The type of the derivative of @f$ \phi @f$. + using vector_span_type = std::conditional_t< + ddcHelper::type_seq_length_v == 1, + chunk_span_type, + VectorFieldSpan>; + + /// @brief The DiscreteDomain describing the batch dimensions. + using batch_domain_type = typename ddc::detail::convert_type_seq_to_discrete_domain; + /// @brief The DiscreteElement for indexing a batch dimension. + using batch_element_type = typename batch_domain_type::discrete_element_type; + + /// @brief The type of the domain on which the equation is defined. + using laplacian_domain_type = ddc::DiscreteDomain; + + /// @brief The layout space of the ChunkSpans passed to operator(). + using layout_space = LayoutSpace; + /// @brief The space (CPU/GPU) where the ChunkSpans passed to operator() are saved. + using memory_space = MemorySpace; + +public: + /** + * @brief An operator which calculates the solution @f$\phi@f$ to Poisson's equation: + * @f$ -\Delta \phi = \rho @f$ + * + * @param[out] phi The solution to Poisson's equation. + * @param[in] rho The right-hand side of Poisson's equation. + * + * @return A reference to the solution to Poisson's equation. + */ + virtual chunk_span_type operator()(chunk_span_type phi, chunk_span_type rho) const = 0; + + /** + * @brief An operator which calculates the solution @f$\phi@f$ to Poisson's equation and + * its derivative: + * @f$ - \Delta \phi = \rho @f$ + * @f$ E = - \nabla \phi @f$ + * + * @param[out] phi The solution to Poisson's equation. + * @param[out] E The derivative of the solution to Poisson's equation. + * @param[in] rho The right-hand side of Poisson's equation. + * + * @return A reference to the solution to Poisson's equation. + */ + virtual chunk_span_type operator()(chunk_span_type phi, vector_span_type E, chunk_span_type rho) + const = 0; +}; diff --git a/src/utils/ddc_helper.hpp b/src/utils/ddc_helper.hpp index 62972ed87..6874a4bde 100644 --- a/src/utils/ddc_helper.hpp +++ b/src/utils/ddc_helper.hpp @@ -155,3 +155,12 @@ using host_t = on_memory_space_t; */ template using device_t = on_memory_space_t; + +namespace ddcHelper { +/// A helper to determine the number of tags in a type sequence. +template +constexpr std::size_t type_seq_length_v = std::numeric_limits::max(); + +template +constexpr std::size_t type_seq_length_v> = sizeof...(Tags); +} // namespace ddcHelper diff --git a/src/utils/vector_field_span.hpp b/src/utils/vector_field_span.hpp index 1f85e4f54..d25336158 100644 --- a/src/utils/vector_field_span.hpp +++ b/src/utils/vector_field_span.hpp @@ -117,13 +117,6 @@ class VectorFieldSpan { } - template < - class... Chunks, - class = std::enable_if_t...>>> - KOKKOS_FUNCTION constexpr VectorFieldSpan(Chunks&&... chunks) : base_type(std::move(chunks)...) - { - } - /** Constructs a new VectorFieldSpan by copy of a chunk, yields a new view to the same data * @param other the VectorFieldSpan to move */ @@ -236,6 +229,16 @@ class VectorFieldSpan { } + /** Constructs a new VectorFieldSpan containing references to ChunkSpans. + * @param chunks The ChunkSpans. + */ + template < + class... ChunkSpan, + class = std::enable_if_t...>>> + KOKKOS_FUNCTION constexpr VectorFieldSpan(ChunkSpan... chunks) : base_type(std::move(chunks)...) + { + } + /** Copy-assigns a new value to this VectorFieldSpan, yields a new view to the same data * @param other the VectorFieldSpan to copy * @return *this diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9b49a7a8a..09f9aecd6 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -18,5 +18,6 @@ gtest_discover_tests(unit_tests_common DISCOVERY_MODE PRE_TEST) add_subdirectory(geometryXVx) add_subdirectory(geometryXYVxVy) add_subdirectory(geometryRTheta) +add_subdirectory(pde_solvers) add_subdirectory(timestepper) add_subdirectory(utils) diff --git a/tests/geometryXVx/CMakeLists.txt b/tests/geometryXVx/CMakeLists.txt index e1bb4f27b..8a2fa1356 100644 --- a/tests/geometryXVx/CMakeLists.txt +++ b/tests/geometryXVx/CMakeLists.txt @@ -85,7 +85,6 @@ gtest_discover_tests(unit_tests_advection target_sources(unit_tests_xperiod_vx PRIVATE femperiodicpoissonsolver.cpp - fftpoissonsolver.cpp ) target_sources(unit_tests_xnonperiod_vx PRIVATE femnonperiodicpoissonsolver.cpp) diff --git a/tests/geometryXVx/fftpoissonsolver.cpp b/tests/geometryXVx/fftpoissonsolver.cpp deleted file mode 100644 index 782995e2a..000000000 --- a/tests/geometryXVx/fftpoissonsolver.cpp +++ /dev/null @@ -1,193 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include - -#include -#include - -#include -#include - -#include "chargedensitycalculator.hpp" -#include "fftqnsolver.hpp" -#include "geometry.hpp" -#include "neumann_spline_quadrature.hpp" -#include "species_info.hpp" - -TEST(FftQNSolver, CosineSource) -{ - CoordX const x_min(0.0); - CoordX const x_max(2.0 * M_PI); - IVectX const x_size(100); - - CoordVx const vx_min(-0.5); - CoordVx const vx_max(0.5); - IVectVx const vx_size(11); - - IVectSp const nb_species(2); - IDomainSp const dom_sp(IndexSp(0), nb_species); - IndexSp const my_ielec = dom_sp.front(); - IndexSp const my_iion = dom_sp.back(); - - // Creating mesh & supports - ddc::init_discrete_space(x_min, x_max, x_size); - - ddc::init_discrete_space(vx_min, vx_max, vx_size); - - ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); - ddc::init_discrete_space(SplineInterpPointsVx::get_sampling()); - ddc::DiscreteDomain interpolation_domain_x(SplineInterpPointsX::get_domain()); - ddc::DiscreteDomain interpolation_domain_vx(SplineInterpPointsVx::get_domain()); - - SplineVxBuilder_1d const builder_vx(interpolation_domain_vx); - - IDomainX const gridx = interpolation_domain_x; - IDomainVx const gridvx = interpolation_domain_vx; - IDomainSp const gridsp = IDomainSp(my_iion, IVectSp(1)); - - IDomainSpXVx const mesh(gridsp, gridx, gridvx); - - ddc::init_discrete_space(ddc::init_fourier_space(ddc::select(mesh))); - - // Creating operators - - host_t> charges(dom_sp); - charges(my_ielec) = -1; - charges(my_iion) = 1; - host_t masses(dom_sp); - ddc::parallel_fill(masses, 1); - - // Initialization of the distribution function - ddc::init_discrete_space(std::move(charges), std::move(masses)); - - host_t const quadrature_coeffs_host - = neumann_spline_quadrature_coefficients(gridvx, builder_vx); - auto const quadrature_coeffs = ddc::create_mirror_view_and_copy( - Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); - ChargeDensityCalculator rhs(quadrature_coeffs); - FftQNSolver poisson(rhs); - - host_t electrostatic_potential_host(gridx); - host_t electric_field_host(gridx); - host_t allfdistribu_host(mesh); - - // Initialization of the distribution function --> fill values - for (IndexSp const isp : gridsp) { - for (IndexX const ix : gridx) { - double fdistribu_val = cos(ddc::coordinate(ix)); - for (IndexVx const iv : gridvx) { - allfdistribu_host(isp, ix, iv) = fdistribu_val; - } - } - } - DFieldX electrostatic_potential(gridx); - DFieldX electric_field(gridx); - DFieldSpXVx allfdistribu(mesh); - - ddc::parallel_deepcopy(allfdistribu, allfdistribu_host); - poisson(electrostatic_potential, electric_field, allfdistribu); - ddc::parallel_deepcopy(electric_field_host, electric_field); - ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); - - double error_pot = 0.0; - double error_field = 0.0; - - for (IndexX const ix : gridx) { - double const exact_pot = cos(ddc::coordinate(ix)); - error_pot = fmax(fabs(electrostatic_potential_host(ix) - exact_pot), error_pot); - double const exact_field = sin(ddc::coordinate(ix)); - error_field = fmax(fabs(electric_field_host(ix) - exact_field), error_field); - } - EXPECT_LE(error_pot, 1e-8); - EXPECT_LE(error_field, 1e-6); -} - -TEST(FftQNSolver, CosineSourceParallel) -{ - CoordX const x_min(0.0); - CoordX const x_max(2.0 * M_PI); - IVectX const x_size(100); - - CoordVx const vx_min(-0.5); - CoordVx const vx_max(0.5); - IVectVx const vx_size(11); - - IVectSp const nb_species(2); - IDomainSp const dom_sp(IndexSp(0), nb_species); - IndexSp const my_ielec = dom_sp.front(); - IndexSp const my_iion = dom_sp.back(); - - // Creating mesh & supports - ddc::init_discrete_space(x_min, x_max, x_size); - - ddc::init_discrete_space(vx_min, vx_max, vx_size); - - ddc::init_discrete_space(SplineInterpPointsX::get_sampling()); - ddc::init_discrete_space(SplineInterpPointsVx::get_sampling()); - ddc::DiscreteDomain interpolation_domain_x(SplineInterpPointsX::get_domain()); - ddc::DiscreteDomain interpolation_domain_vx(SplineInterpPointsVx::get_domain()); - - SplineVxBuilder_1d const builder_vx(interpolation_domain_vx); - - IDomainX const gridx = interpolation_domain_x; - IDomainVx const gridvx = interpolation_domain_vx; - IDomainSp const gridsp = IDomainSp(my_iion, IVectSp(1)); - - IDomainSpXVx const mesh(gridsp, gridx, gridvx); - - ddc::init_discrete_space(ddc::init_fourier_space(ddc::select(mesh))); - - // Creating operators - - host_t> charges(dom_sp); - charges(my_ielec) = -1; - charges(my_iion) = 1; - host_t masses(dom_sp); - ddc::parallel_fill(masses, 1); - - // Initialization of the distribution function - ddc::init_discrete_space(std::move(charges), std::move(masses)); - - host_t const quadrature_coeffs_host - = neumann_spline_quadrature_coefficients(gridvx, builder_vx); - DFieldVx quadrature_coeffs(quadrature_coeffs_host.domain()); - - ddc::parallel_deepcopy(quadrature_coeffs, quadrature_coeffs_host); - ChargeDensityCalculator rhs(quadrature_coeffs); - FftQNSolver poisson(rhs); - - host_t electrostatic_potential_host(gridx); - host_t electric_field_host(gridx); - host_t allfdistribu_host(mesh); - - // Initialization of the distribution function --> fill values - for (IndexSp const isp : gridsp) { - for (IndexX const ix : gridx) { - double fdistribu_val = cos(ddc::coordinate(ix)); - for (IndexVx const iv : gridvx) { - allfdistribu_host(isp, ix, iv) = fdistribu_val; - } - } - } - DFieldX electrostatic_potential(gridx); - DFieldX electric_field(gridx); - DFieldSpXVx allfdistribu(mesh); - - ddc::parallel_deepcopy(allfdistribu, allfdistribu_host); - poisson(electrostatic_potential, electric_field, allfdistribu); - ddc::parallel_deepcopy(electric_field_host, electric_field); - ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); - - double error_pot = 0.0; - double error_field = 0.0; - - for (IndexX const ix : gridx) { - double const exact_pot = cos(ddc::coordinate(ix)); - error_pot = fmax(fabs(electrostatic_potential_host(ix) - exact_pot), error_pot); - double const exact_field = sin(ddc::coordinate(ix)); - error_field = fmax(fabs(electric_field_host(ix) - exact_field), error_field); - } - EXPECT_LE(error_pot, 1e-8); - EXPECT_LE(error_field, 1e-6); -} diff --git a/tests/geometryXVx/geometryMX/predcorr_hybrid.cpp b/tests/geometryXVx/geometryMX/predcorr_hybrid.cpp index 79e8e549a..8c0cc86b1 100644 --- a/tests/geometryXVx/geometryMX/predcorr_hybrid.cpp +++ b/tests/geometryXVx/geometryMX/predcorr_hybrid.cpp @@ -20,7 +20,7 @@ #include "femnonperiodicqnsolver.hpp" #endif #include "Lagrange_interpolator.hpp" -#include "fftqnsolver.hpp" +#include "fft_poisson_solver.hpp" #include "geometry.hpp" #include "irighthandside.hpp" #include "maxwellianequilibrium.hpp" @@ -28,6 +28,7 @@ #include "nullfluidsolver.hpp" #include "predcorr.hpp" #include "predcorr_hybrid.hpp" +#include "qnsolver.hpp" #include "quadrature.hpp" #include "singlemodeperturbinitialization.hpp" #include "species_info.hpp" @@ -221,10 +222,8 @@ TEST(GeometryXM, PredCorrHybrid) quadrature_coeffs_host.span_view()); ChargeDensityCalculator rhs(quadrature_coeffs); #ifdef PERIODIC_RDIMX - IDomainSpXVx const meshSpXVx(dom_kinsp, meshXVx); - ddc::init_discrete_space( - ddc::init_fourier_space(ddc::select(meshSpXVx))); - FftQNSolver const poisson(rhs); + FFTPoissonSolver fft_poisson_solver(meshX); + QNSolver const poisson(fft_poisson_solver, rhs); #else FemNonPeriodicQNSolver const poisson(builder_x_poisson, spline_x_evaluator_poisson, rhs); #endif diff --git a/tests/geometryXYVxVy/CMakeLists.txt b/tests/geometryXYVxVy/CMakeLists.txt index 08943d829..3ffa66bf5 100644 --- a/tests/geometryXYVxVy/CMakeLists.txt +++ b/tests/geometryXYVxVy/CMakeLists.txt @@ -5,7 +5,6 @@ cmake_minimum_required(VERSION 3.15) add_executable(unit_tests_xy_vxvy ../main.cpp quadrature.cpp - fft_poisson.cpp ) target_compile_features(unit_tests_xy_vxvy PUBLIC cxx_std_17) target_link_libraries(unit_tests_xy_vxvy @@ -16,8 +15,6 @@ target_link_libraries(unit_tests_xy_vxvy paraconf::paraconf gslx::geometry_xyvxvy sll::SLL - gslx::advection - gslx::poisson_xy gslx::quadrature ) diff --git a/tests/geometryXYVxVy/fft_poisson.cpp b/tests/geometryXYVxVy/fft_poisson.cpp deleted file mode 100644 index 1a528ac8c..000000000 --- a/tests/geometryXYVxVy/fft_poisson.cpp +++ /dev/null @@ -1,155 +0,0 @@ -// SPDX-License-Identifier: MIT - -#include -#include - -#include "chargedensitycalculator.hpp" -#include "fftqnsolver.hpp" -#include "geometry.hpp" -#include "neumann_spline_quadrature.hpp" -#include "quadrature.hpp" -#include "species_info.hpp" - -namespace { - -static void TestFftQNSolverCosineSource() -{ - CoordX const x_min(0.0); - CoordX const x_max(2.0 * M_PI); - IVectX const x_size(10); - - CoordY const y_min(0.0); - CoordY const y_max(2.0 * M_PI); - IVectY const y_size(10); - - CoordVx const vx_min(-0.5); - CoordVx const vx_max(0.5); - IVectVx const vx_size(11); - - CoordVy const vy_min(-0.5); - CoordVy const vy_max(0.5); - IVectVy const vy_size(11); - - IVectSp const nb_species(2); - IDomainSp const dom_sp(IndexSp(0), nb_species); - IndexSp const my_ielec = dom_sp.front(); - IndexSp const my_iion = dom_sp.back(); - - // 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(SplineInterpPointsVx::get_sampling()); - ddc::init_discrete_space(SplineInterpPointsVy::get_sampling()); - - ddc::init_discrete_space(IDimX::init(x_min, x_max, x_size + 1)); - ddc::init_discrete_space(IDimY::init(y_min, y_max, y_size + 1)); - - IDomainSp const gridsp = IDomainSp(my_iion, IVectSp(1)); - - IDomainX gridx = IDomainX(IndexX(0), x_size); - IDomainY gridy = IDomainY(IndexY(0), y_size); - - IDomainXY gridxy(gridx, gridy); - - IDomainVx gridvx = SplineInterpPointsVx::get_domain(); - IDomainVy gridvy = SplineInterpPointsVy::get_domain(); - - IDomainVxVy gridvxvy(gridvx, gridvy); - - IDomainXYVxVy const meshxyvxvy(gridxy, gridvxvy); - - SplineVxBuilder_1d const builder_vx_1d(gridvx); - SplineVyBuilder_1d const builder_vy_1d(gridvy); - - IDomainSpXYVxVy const mesh(gridsp, meshxyvxvy); - - ddc::init_discrete_space(ddc::init_fourier_space(gridx)); - ddc::init_discrete_space(ddc::init_fourier_space(gridy)); - - // Initialise infomation about species - host_t> charges(dom_sp); - charges(my_ielec) = -1; - charges(my_iion) = 1; - host_t masses(dom_sp); - ddc::parallel_fill(masses, 1); - - ddc::init_discrete_space(std::move(charges), std::move(masses)); - - host_t const quadrature_coeffs_host - = neumann_spline_quadrature_coefficients(gridvxvy, builder_vx_1d, builder_vy_1d); - auto quadrature_coeffs = ddc::create_mirror_view_and_copy( - Kokkos::DefaultExecutionSpace(), - quadrature_coeffs_host.span_view()); - ChargeDensityCalculator rhs(quadrature_coeffs); - FftQNSolver poisson(rhs); - - 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 - 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 = 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(FftQNSolver, CosineSource) -{ - TestFftQNSolverCosineSource(); -} diff --git a/tests/pde_solvers/CMakeLists.txt b/tests/pde_solvers/CMakeLists.txt new file mode 100644 index 000000000..b93fd1d8d --- /dev/null +++ b/tests/pde_solvers/CMakeLists.txt @@ -0,0 +1,21 @@ +# SPDX-License-Identifier: MIT + +cmake_minimum_required(VERSION 3.15) + +include(GoogleTest) + +add_executable(unit_tests_pde_solvers + fftpoissonsolver.cpp + ../main.cpp +) +target_compile_features(unit_tests_pde_solvers PUBLIC cxx_std_17) +target_link_libraries(unit_tests_pde_solvers + PUBLIC + DDC::DDC + GTest::gtest + GTest::gmock + gslx::utils + gslx::pde_solvers +) + +gtest_discover_tests(unit_tests_pde_solvers) diff --git a/tests/pde_solvers/fftpoissonsolver.cpp b/tests/pde_solvers/fftpoissonsolver.cpp new file mode 100644 index 000000000..2c938da0d --- /dev/null +++ b/tests/pde_solvers/fftpoissonsolver.cpp @@ -0,0 +1,244 @@ +// SPDX-License-Identifier: MIT + +#include + +#include + +#include +#include + +#include "fft_poisson_solver.hpp" + +namespace FFTPoissonSolverTest { + +struct RDimX +{ + /// @brief A boolean indicating if the dimension is periodic. + static bool constexpr PERIODIC = true; +}; + +struct RDimY +{ + /// @brief A boolean indicating if the dimension is periodic. + static bool constexpr PERIODIC = true; +}; + +using CoordX = ddc::Coordinate; +using CoordY = ddc::Coordinate; + +struct IDimX : ddc::UniformPointSampling +{ +}; +struct IDimY : ddc::UniformPointSampling +{ +}; + +using IndexX = ddc::DiscreteElement; +using IndexY = ddc::DiscreteElement; +using IndexXY = ddc::DiscreteElement; + +using IVectX = ddc::DiscreteVector; +using IVectY = ddc::DiscreteVector; + +using IDomainX = ddc::DiscreteDomain; +using IDomainY = ddc::DiscreteDomain; +using IDomainXY = ddc::DiscreteDomain; + +using DFieldX = device_t>; +using DFieldY = device_t>; +using DFieldXY = device_t>; + +using DSpanX = device_t>; +using DSpanY = device_t>; +using DSpanXY = device_t>; + + +TEST(FftPoissonSolver, CosineSource) +{ + CoordX const x_min(0.0); + CoordX const x_max(2.0 * M_PI); + IVectX const x_size(100); + + // Creating mesh & supports + ddc::init_discrete_space(IDimX::init(x_min, x_max, x_size + 1)); + IDomainX gridx = IDomainX(IndexX(0), x_size); + + // Creating operators + FFTPoissonSolver poisson(gridx); + + host_t electrostatic_potential_host(gridx); + host_t electric_field_host(gridx); + host_t rhs_host(gridx); + + // Initialization of the distribution function --> fill values + for (IndexX const ix : gridx) { + rhs_host(ix) = cos(ddc::coordinate(ix)); + } + DFieldX electrostatic_potential(gridx); + DFieldX electric_field(gridx); + DFieldX rhs(gridx); + + ddc::parallel_deepcopy(rhs, rhs_host); + poisson(electrostatic_potential, electric_field, rhs); + ddc::parallel_deepcopy(electric_field_host, electric_field); + ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); + + double error_pot = 0.0; + double error_field = 0.0; + + for (IndexX const ix : gridx) { + double const exact_pot = cos(ddc::coordinate(ix)); + error_pot = fmax(fabs(electrostatic_potential_host(ix) - exact_pot), error_pot); + double const exact_field = sin(ddc::coordinate(ix)); + error_field = fmax(fabs(electric_field_host(ix) - exact_field), error_field); + } + EXPECT_LE(error_pot, 1e-8); + EXPECT_LE(error_field, 1e-6); +} + +TEST(FftPoissonSolver, BatchedCosineSource) +{ + CoordX const x_min(0.0); + CoordX const x_max(2.0 * M_PI); + IVectX const x_size(10); + + CoordY const y_min(0.0); + CoordY const y_max(2.0 * M_PI); + IVectY const y_size(100); + + // Creating mesh & supports + ddc::init_discrete_space(IDimX::init(x_min, x_max, x_size + 1)); + ddc::init_discrete_space(IDimY::init(y_min, y_max, y_size + 1)); + IDomainX gridx = IDomainX(IndexX(0), x_size); + IDomainY gridy = IDomainY(IndexY(0), y_size); + IDomainXY gridxy(gridx, gridy); + + // Creating operators + FFTPoissonSolver poisson(gridy); + + host_t electrostatic_potential_host(gridxy); + host_t electric_field_host(gridxy); + host_t rhs_host(gridxy); + + // Initialization of the distribution function --> fill values + for (IndexX const ix : gridx) { + double const c = ix.uid() + 1; + for (IndexY const iy : gridy) { + rhs_host(ix, iy) = ipow(c, 2) * cos(c * ddc::coordinate(iy)); + } + } + DFieldXY electrostatic_potential(gridxy); + DFieldXY electric_field(gridxy); + DFieldXY rhs(gridxy); + + ddc::parallel_deepcopy(rhs, rhs_host); + poisson(electrostatic_potential, electric_field, rhs); + ddc::parallel_deepcopy(electric_field_host, electric_field); + ddc::parallel_deepcopy(electrostatic_potential_host, electrostatic_potential); + + double error_pot = 0.0; + double error_field = 0.0; + + for (IndexX const ix : gridx) { + double const c = ix.uid() + 1; + for (IndexY const iy : gridy) { + double const exact_pot = cos(c * ddc::coordinate(iy)); + error_pot = fmax(fabs(electrostatic_potential_host(ix, iy) - exact_pot), error_pot); + double const exact_field = c * sin(c * ddc::coordinate(iy)); + error_field = fmax(fabs(electric_field_host(ix, iy) - exact_field), error_field); + } + } + EXPECT_LE(error_pot, 1e-8); + EXPECT_LE(error_field, 1e-6); +} + +static void TestFftPoissonSolver2DCosineSource() +{ + CoordX const x_min(0.0); + CoordX const x_max(2.0 * M_PI); + IVectX const x_size(10); + + CoordY const y_min(0.0); + CoordY const y_max(2.0 * M_PI); + IVectY const y_size(10); + + // Creating mesh & supports + ddc::init_discrete_space(IDimX::init(x_min, x_max, x_size + 1)); + ddc::init_discrete_space(IDimY::init(y_min, y_max, y_size + 1)); + + IDomainX gridx = IDomainX(IndexX(0), x_size); + IDomainY gridy = IDomainY(IndexY(0), y_size); + + IDomainXY gridxy(gridx, gridy); + + FFTPoissonSolver poisson(gridxy); + + DFieldXY electrostatic_potential_alloc(gridxy); + VectorField, ddc::DeviceAllocator> + electric_field_alloc(gridxy); + DFieldXY rhs_alloc(gridxy); + + DSpanXY electrostatic_potential = electrostatic_potential_alloc.span_view(); + VectorFieldSpan electric_field = electric_field_alloc.span_view(); + DSpanXY rhs = rhs_alloc.span_view(); + + // Initialization of the distribution function --> fill values + ddc::parallel_for_each( + Kokkos::DefaultExecutionSpace(), + gridxy, + KOKKOS_LAMBDA(IndexXY const ixy) { + IndexX ix = ddc::select(ixy); + IndexY iy = ddc::select(ixy); + double x = ddc::coordinate(ix); + double y = ddc::coordinate(iy); + rhs(ixy) = Kokkos::cos(x) + Kokkos::cos(y); + }); + + poisson(electrostatic_potential_alloc, electric_field_alloc, rhs_alloc); + + 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); + }); + ddc::ChunkSpan electric_field_x = electric_field.get(); + 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); + }); + ddc::ChunkSpan electric_field_y = electric_field.get(); + 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); +} + +TEST(FftPoissonSolver2D, CosineSource) +{ + TestFftPoissonSolver2DCosineSource(); +} + +} // namespace FFTPoissonSolverTest