Resolve "[GYSELAX WEEK] 1X1V sheath : Add the simulation : periodic s…

…pline in X + Lagrange in V "

simulations/geometryXVx/sheath/sheath.cpp

@@ -28,6 +28,8 @@
 #else
 #include "femnonperiodicpoissonsolver.hpp"
 #endif
+#include "Lagrange_interpolator.hpp"
+#include "fftpoissonsolver_splinex.hpp"
 #include "geometry.hpp"
 #include "irighthandside.hpp"
 #include "kinetic_source.hpp"
@@ -202,10 +204,15 @@ int main(int argc, char** argv)
     SplineEvaluator<BSplinesVx> const spline_vx_evaluator(bv_v_min, bv_v_max);
 
     PreallocatableSplineInterpolator const spline_vx_interpolator(builder_vx, spline_vx_evaluator);
+    IVectVx static constexpr gwvx {0};
+    LagrangeInterpolator<IDimVx, BCond::DIRICHLET, BCond::DIRICHLET> const
+            lagrange_vx_evaluator(3, interpolation_domain_vx, gwvx);
+    PreallocatableLagrangeInterpolator<IDimVx, BCond::DIRICHLET, BCond::DIRICHLET> const
+            lagrange_vx_interpolator(lagrange_vx_evaluator);
 
     BslAdvectionSpatial<GeometryXVx, IDimX> const advection_x(spline_x_interpolator);
 
-    BslAdvectionVelocity<GeometryXVx, IDimVx> const advection_vx(spline_vx_interpolator);
+    BslAdvectionVelocity<GeometryXVx, IDimVx> const advection_vx(lagrange_vx_interpolator);
 
     // Creating of mesh for output saving
     IDomainX const gridx = ddc::select<IDimX>(meshSpXVx);
@@ -285,16 +292,20 @@ int main(int argc, char** argv)
         collisions_inter.emplace(meshSpXVx, PCpp_double(conf_voicexx, ".CollisionsInfo.nustar0"));
         rhs_operators.emplace_back(*collisions_inter);
     }
-
     SplitVlasovSolver const vlasov(advection_x, advection_vx);
     SplitRightHandSideSolver const boltzmann(vlasov, rhs_operators);
 
+    DFieldVx const quadrature_coeffs
+            = simpson_quadrature_coefficients_1d(allfdistribu.domain<IDimVx>());
+    Quadrature<IDimVx> const integrate_v(quadrature_coeffs);
 #ifdef PERIODIC_RDIMX
-    using FemPoissonSolverX = FemPeriodicPoissonSolver;
+    ddc::init_fourier_space<RDimX>(ddc::select<IDimX>(meshSpXVx));
+    FftPoissonSolverSplineX const poisson(builder_x, spline_x_evaluator, integrate_v);
 #else
-    using FemPoissonSolverX = FemNonPeriodicPoissonSolver;
+    FemNonPeriodicPoissonSolver const
+            poisson(builder_x, spline_x_evaluator, builder_vx, spline_vx_evaluator);
 #endif
-    FemPoissonSolverX const poisson(builder_x, spline_x_evaluator, builder_vx, spline_vx_evaluator);
+
 
     PredCorr const predcorr(boltzmann, poisson);
 

src/geometryXVx/poisson/CMakeLists.txt

@@ -4,6 +4,7 @@ foreach(GEOMETRY_VARIANT IN LISTS GEOMETRY_XVx_VARIANTS_LIST)
 
 add_library("poisson_${GEOMETRY_VARIANT}" STATIC
     chargedensitycalculator.cpp
+    chargedensitycalculator_splineless.cpp
     electricfield.cpp
     nullpoissonsolver.cpp
 )
@@ -21,6 +22,7 @@ target_link_libraries("poisson_${GEOMETRY_VARIANT}"
     PUBLIC
         DDC::DDC
         sll::splines
+        gslx::quadrature
         gslx::geometry_${GEOMETRY_VARIANT}
         gslx::speciesinfo
 )
@@ -32,6 +34,7 @@ endforeach()
 target_sources(poisson_xperiod_vx
     PRIVATE
         fftpoissonsolver.cpp
+        fftpoissonsolver_splinex.cpp
         femperiodicpoissonsolver.cpp
 )
 

src/geometryXVx/poisson/chargedensitycalculator_splineless.cpp

@@ -0,0 +1,32 @@
+// SPDX-License-Identifier: MIT
+
+#include <ddc/ddc.hpp>
+
+#include "chargedensitycalculator_splineless.hpp"
+#include "quadrature.hpp"
+#include "simpson_quadrature.hpp"
+
+
+
+ChargeDensityCalculatorSplineless::ChargeDensityCalculatorSplineless(Quadrature<IDimVx> const& quad)
+    : m_quad(quad)
+{
+}
+
+void ChargeDensityCalculatorSplineless::operator()(DSpanX const rho, DViewSpXVx const allfdistribu)
+        const
+{
+    IndexSp const last_kin_species = allfdistribu.domain<IDimSp>().back();
+    IndexSp const last_species = ddc::discrete_space<IDimSp>().charges().domain().back();
+    double chargedens_adiabspecies = 0.;
+    if (last_kin_species != last_species) {
+        chargedens_adiabspecies = double(charge(last_species));
+    }
+
+    ddc::for_each(rho.domain(), [&](IndexX const ix) {
+        rho(ix) = chargedens_adiabspecies;
+        ddc::for_each(ddc::get_domain<IDimSp>(allfdistribu), [&](IndexSp const isp) {
+            rho(ix) += charge(isp) * m_quad(allfdistribu[isp][ix]);
+        });
+    });
+}

src/geometryXVx/poisson/chargedensitycalculator_splineless.hpp

@@ -0,0 +1,30 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <ddc/ddc.hpp>
+
+#include <geometry.hpp>
+
+#include "quadrature.hpp"
+#include "simpson_quadrature.hpp"
+
+/**
+ * @brief A class which computes charges density .      
+ */
+class ChargeDensityCalculatorSplineless
+{
+    const Quadrature<IDimVx>& m_quad;
+
+public:
+    /**
+ * @brief Create a ChargeDensityCalculatorSplineless object, integration is done using composite Simpson rule.
+ */
+    ChargeDensityCalculatorSplineless(const Quadrature<IDimVx>& quad);
+    /**
+ * @brief Computes the charge density rho from the distribution function.
+ * @param[in, out] rho
+ * @param[in] allfdistribu 
+ */
+    void operator()(DSpanX rho, DViewSpXVx allfdistribu) const;
+};

src/geometryXVx/poisson/fftpoissonsolver_splinex.cpp

@@ -0,0 +1,74 @@
+// SPDX-License-Identifier: MIT
+
+#include <cassert>
+#include <cmath>
+#include <complex>
+#include <iostream>
+
+#include <ddc/ddc.hpp>
+#include <ddc/kernels/fft.hpp>
+
+#include <sll/bsplines_uniform.hpp>
+#include <sll/null_boundary_value.hpp>
+#include <sll/spline_builder.hpp>
+#include <sll/spline_evaluator.hpp>
+
+#include <geometry.hpp>
+
+#include "fftpoissonsolver_splinex.hpp"
+
+FftPoissonSolverSplineX::FftPoissonSolverSplineX(
+        SplineXBuilder const& spline_x_builder,
+        SplineEvaluator<BSplinesX> const& spline_x_evaluator,
+        Quadrature<IDimVx> const& quad)
+    : m_compute_rho(quad)
+    , m_electric_field(spline_x_builder, spline_x_evaluator)
+{
+}
+
+// 1- Inner solvers shall be passed in the constructor
+// 2- Should it take an array of distribution functions ?
+void FftPoissonSolverSplineX::operator()(
+        DSpanX const electrostatic_potential,
+        DSpanX const electric_field,
+        DViewSpXVx const allfdistribu) const
+{
+    assert(electrostatic_potential.domain() == ddc::get_domain<IDimX>(allfdistribu));
+    IDomainX const x_dom = electrostatic_potential.domain();
+
+    // Compute the RHS of the Poisson equation.
+    ddc::Chunk<double, IDomainX> rho(x_dom);
+    DFieldVx contiguous_slice_vx(allfdistribu.domain<IDimVx>());
+    m_compute_rho(rho, allfdistribu);
+
+    // Build a mesh in the fourier space, for N points
+    IDomainFx const k_mesh = ddc::FourierMesh(x_dom, false);
+
+    ddc::Chunk<Kokkos::complex<double>, IDomainFx> intermediate_chunk
+            = ddc::Chunk(k_mesh, ddc::HostAllocator<Kokkos::complex<double>>());
+
+    // Compute FFT(rho)
+    ddc::FFT_Normalization norm = ddc::FFT_Normalization::BACKWARD;
+    ddc::
+            fft(Kokkos::DefaultHostExecutionSpace(),
+                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.remove_first(IVectFx(1)), [&](IndexFx const ikx) {
+        intermediate_chunk(ikx) = intermediate_chunk(ikx) / (coordinate(ikx) * coordinate(ikx));
+    });
+
+    // Perform the inverse 1D FFT of the solution to deduce the electrostatic potential
+    ddc::
+            ifft(Kokkos::DefaultHostExecutionSpace(),
+                 electrostatic_potential.span_view(),
+                 intermediate_chunk.span_view(),
+                 ddc::kwArgs_fft {norm});
+
+    // Compute efield = -dPhi/dx where Phi is the electrostatic potential
+    m_electric_field(electric_field, electrostatic_potential);
+}

src/geometryXVx/poisson/fftpoissonsolver_splinex.hpp

@@ -0,0 +1,27 @@
+// SPDX-License-Identifier: MIT
+
+#pragma once
+
+#include <sll/spline_builder.hpp>
+#include <sll/spline_evaluator.hpp>
+
+#include "chargedensitycalculator_splineless.hpp"
+#include "electricfield.hpp"
+#include "ipoissonsolver.hpp"
+
+class FftPoissonSolverSplineX : public IPoissonSolver
+{
+    ChargeDensityCalculatorSplineless m_compute_rho;
+    ElectricField m_electric_field;
+
+public:
+    FftPoissonSolverSplineX(
+            SplineXBuilder const& spline_x_builder,
+            SplineEvaluator<BSplinesX> const& spline_x_evaluator,
+            Quadrature<IDimVx> const& quad);
+
+    ~FftPoissonSolverSplineX() override = default;
+
+    void operator()(DSpanX electrostatic_potential, DSpanX electric_field, DViewSpXVx allfdistribu)
+            const override;
+};