Decouple ChargeDensityCalculator and FFT PoissonSolver

Part 1 of #218

See merge request gysela-developpers/gyselalibxx!476

--------------------------------------------

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<IDimFx>(
-            ddc::init_fourier_space<IDimFx>(ddc::select<IDimX>(meshSpXVx)));
-
     host_t<DFieldVx> 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<IDomainX, IDomainX, Kokkos::DefaultExecutionSpace> 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);
 

simulations/geometryXVx/landau/landau_fft.cpp

@@ -17,14 +17,15 @@
 #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"
 #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,17 +210,16 @@ int main(int argc, char** argv)
 
     SplitVlasovSolver const vlasov(advection_x, advection_vx);
 
-    ddc::init_discrete_space<IDimFx>(
-            ddc::init_fourier_space<IDimFx>(ddc::select<IDimX>(meshSpXVx)));
-
     host_t<DFieldVx> const quadrature_coeffs_host
             = neumann_spline_quadrature_coefficients(gridvx, builder_vx_poisson);
 
     auto quadrature_coeffs = ddc::create_mirror_view_and_copy(
             Kokkos::DefaultExecutionSpace(),
             quadrature_coeffs_host.span_view());
     ChargeDensityCalculator rhs(quadrature_coeffs);
-    FftQNSolver const poisson(rhs);
+    FFTPoissonSolver<IDomainX, IDomainX, Kokkos::DefaultExecutionSpace> fft_poisson_solver(
+            interpolation_domain_x);
+    QNSolver const poisson(fft_poisson_solver, rhs);
 
     PredCorr const predcorr(vlasov, poisson);
 

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<IDimFx>(
-            ddc::init_fourier_space<IDimFx>(ddc::select<IDimX>(meshSpXVx)));
-    FftQNSolver const poisson(rhs);
+    FFTPoissonSolver<IDomainX, IDomainX, Kokkos::DefaultExecutionSpace> fft_poisson_solver(gridx);
+    QNSolver const poisson(fft_poisson_solver, rhs);
 #else
     FemNonPeriodicQNSolver const poisson(builder_x_poisson, spline_x_evaluator_poisson, rhs);
 #endif

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<IDimFx>(
-            ddc::init_fourier_space<IDimFx>(ddc::select<IDimX>(meshSpXVx)));
-    FftQNSolver const poisson(rhs);
+    FFTPoissonSolver<IDomainX, IDomainX, Kokkos::DefaultExecutionSpace> 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

simulations/geometryXYVxVy/landau/landau4d_fft.cpp

@@ -15,14 +15,15 @@
 
 #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"
 #include "paraconfpp.hpp"
 #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<IDimX>());
     IDomainY interpolation_domain_y(SplineInterpPointsY::get_domain<IDimY>());
+    IDomainXY interpolation_domain_xy(interpolation_domain_x, interpolation_domain_y);
     IDomainVx interpolation_domain_vx(SplineInterpPointsVx::get_domain<IDimVx>());
     IDomainVy interpolation_domain_vy(SplineInterpPointsVy::get_domain<IDimVy>());
     IDomainVxVy interpolation_domain_vxvy(interpolation_domain_vx, interpolation_domain_vy);
@@ -209,20 +211,17 @@ int main(int argc, char** argv)
 
     SplitVlasovSolver const vlasov(advection_x, advection_y, advection_vx, advection_vy);
 
-    ddc::init_discrete_space<IDimFx>(
-            ddc::init_fourier_space<IDimFx>(ddc::select<IDimX>(meshSpXYVxVy)));
-    ddc::init_discrete_space<IDimFy>(
-            ddc::init_fourier_space<IDimFy>(ddc::select<IDimY>(meshSpXYVxVy)));
-
     host_t<DFieldVxVy> const quadrature_coeffs_host = neumann_spline_quadrature_coefficients(
             interpolation_domain_vxvy,
             builder_vx_1d,
             builder_vy_1d);
     auto quadrature_coeffs = ddc::create_mirror_view_and_copy(
             Kokkos::DefaultExecutionSpace(),
             quadrature_coeffs_host.span_view());
+    FFTPoissonSolver<IDomainXY, IDomainXY, Kokkos::DefaultExecutionSpace> 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);

src/CMakeLists.txt

@@ -13,3 +13,4 @@ add_subdirectory(geometryXYVxVy)
 add_subdirectory(interpolation)
 add_subdirectory(advection)
 add_subdirectory(timestepper)
+add_subdirectory(pde_solvers)

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.
 <!-- - [paraconfpp](./paraconfpp/README.md) - Paraconf utility functions. -->
+- [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).

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<ElementType>). 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.

src/geometryXVx/geometry/geometry.hpp

@@ -427,15 +427,6 @@ using DViewSpVx = ViewSpVx<double>;
 using DViewSpXVx = ViewSpXVx<double>;
 
 
-using RDimFx = ddc::Fourier<RDimX>;
-using CoordFx = ddc::Coordinate<RDimFx>;
-struct IDimFx : ddc::PeriodicSampling<RDimFx>
-{
-};
-using IndexFx = ddc::DiscreteElement<IDimFx>;
-using IVectFx = ddc::DiscreteVector<IDimFx>;
-using IDomainFx = ddc::DiscreteDomain<IDimFx>;
-
 /**
  * @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.
  */

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
 )
 

src/geometryXVx/poisson/qnsolver.cpp

@@ -0,0 +1,36 @@
+// SPDX-License-Identifier: MIT
+
+#include <cassert>
+#include <cmath>
+#include <complex>
+#include <iostream>
+
+#include <ddc/ddc.hpp>
+
+#include <geometry.hpp>
+
+#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<IDimX>(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();
+}