ECP-WarpX · aeriforme · Jul 27, 2024 · Jul 27, 2024 · Jul 27, 2024 · Jul 30, 2024
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -265,6 +265,10 @@ Overall simulation parameters
         In electromagnetic mode, this solver can be used to initialize the species' self fields
         (``<species_name>.initialize_self_fields=1``) provided that the field BCs are PML (``boundary.field_lo,hi = PML``).
 
+          * warpx.use_2d_slices_fft_solver (`bool`, default: 0): if 0, solve Poisson equation in full 3D geometry with the
+            Integrated Green Function solver; if 1, solve Poisson equation in a quasi 3D geometry. In the latter case,
+            the code performes many 2D Poisson solves for each slice at a given :math:`z`.
+
 * ``warpx.self_fields_required_precision`` (`float`, default: 1.e-11)
     The relative precision with which the electrostatic space-charge fields should
     be calculated. More specifically, the space-charge fields are

diff --git a/Examples/Tests/open_bc_poisson_solver/CMakeLists.txt b/Examples/Tests/open_bc_poisson_solver/CMakeLists.txt
@@ -24,3 +24,15 @@ if(WarpX_HEFFTE)
         OFF  # dependency
     )
 endif()
+
+if(WarpX_FFT)
+    add_warpx_test(
+        test_3d_open_bc_poisson_solver_sliced  # name
+        3  # dims
+        2  # nprocs
+        inputs_test_3d_open_bc_poisson_solver_sliced  # inputs
+        analysis.py  # analysis
+        diags/diag1000001  # output
+        OFF  # dependency
+    )
+endif()
diff --git a/Examples/Tests/open_bc_poisson_solver/inputs_test_3d_open_bc_poisson_solver_sliced b/Examples/Tests/open_bc_poisson_solver/inputs_test_3d_open_bc_poisson_solver_sliced
@@ -0,0 +1,3 @@
+FILE = inputs_test_3d_open_bc_poisson_solver
+
+warpx.use_2d_slices_fft_solver = 1
diff --git a/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver_sliced.json b/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver_sliced.json
@@ -0,0 +1,21 @@
+{
+  "lev=0": {
+    "Bx": 100915975.15330593,
+    "By": 157610677.31480408,
+    "Bz": 6.060146872040567e-14,
+    "Ex": 4.725066923360721e+16,
+    "Ey": 3.0253961492956904e+16,
+    "Ez": 2787755.831700608,
+    "rho": 10994013582437.197
+  },
+  "electron": {
+    "particle_momentum_x": 5.701279599510526e-19,
+    "particle_momentum_y": 3.6504531723833074e-19,
+    "particle_momentum_z": 1.145432768297242e-10,
+    "particle_position_x": 17.31408691249785,
+    "particle_position_y": 0.2583691267187801,
+    "particle_position_z": 10066.329600000008,
+    "particle_weight": 19969036501.910976
+  }
+}
+
diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H
@@ -97,7 +97,8 @@ public:
         amrex::Real required_precision,
         amrex::Real absolute_tolerance,
         int max_iters,
-        int verbosity
+        int verbosity,
+        bool is_2d_slices
     ) const;
 
     /**
@@ -158,6 +159,10 @@ public:
      *  2 : convergence progress at every MLMG iteration
      */
     int self_fields_verbosity = 2;
+
+    /** Parameters for FFT Poisson solver aka IGF */
+    // 0: full 3D, 1: many 2D z-slices (quasi-3D)
+    bool is_igf_2d_slices = false;
 };
 
 #endif // WARPX_ELECTROSTATICSOLVER_H_
diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
@@ -28,13 +28,13 @@ void ElectrostaticSolver::ReadParameters () {
 
     // Note that with the relativistic version, these parameters would be
     // input for each species.
-    utils::parser::queryWithParser(
-        pp_warpx, "self_fields_required_precision", self_fields_required_precision);
-    utils::parser::queryWithParser(
-        pp_warpx, "self_fields_absolute_tolerance", self_fields_absolute_tolerance);
-    utils::parser::queryWithParser(
-        pp_warpx, "self_fields_max_iters", self_fields_max_iters);
+    pp_warpx.query("self_fields_required_precision", self_fields_required_precision);
+    pp_warpx.query("self_fields_absolute_tolerance", self_fields_absolute_tolerance);
+    pp_warpx.query("self_fields_max_iters", self_fields_max_iters);
     pp_warpx.query("self_fields_verbosity", self_fields_verbosity);
+
+    // FFT solver flags
+    pp_warpx.query("use_2d_slices_fft_solver", is_igf_2d_slices);
 }
 
 void
@@ -116,7 +116,8 @@ ElectrostaticSolver::computePhi (const amrex::Vector<std::unique_ptr<amrex::Mult
                    Real const required_precision,
                    Real absolute_tolerance,
                    int const max_iters,
-                   int const verbosity) const {
+                   int const verbosity,
+                   bool is_igf_2d_slices) const {
     // create a vector to our fields, sorted by level
     amrex::Vector<amrex::MultiFab *> sorted_rho;
     amrex::Vector<amrex::MultiFab *> sorted_phi;
@@ -195,6 +196,7 @@ ElectrostaticSolver::computePhi (const amrex::Vector<std::unique_ptr<amrex::Mult
         WarpX::grid_type,
         *m_poisson_boundary_handler,
         is_solver_igf_on_lev0,
+        is_igf_2d_slices,
         EB::enabled(),
         WarpX::do_single_precision_comms,
         warpx.refRatio(),

diff --git a/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.cpp
@@ -74,7 +74,7 @@ void LabFrameExplicitES::ComputeSpaceChargeField (
         // Use the AMREX MLMG or the FFT (IGF) solver otherwise
         computePhi(rho_fp, phi_fp, beta, self_fields_required_precision,
                    self_fields_absolute_tolerance, self_fields_max_iters,
-                   self_fields_verbosity);
+                   self_fields_verbosity, is_igf_2d_slices);
 #endif
 
     }

diff --git a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp
@@ -124,7 +124,7 @@ void RelativisticExplicitES::AddSpaceChargeField (
     // Compute the potential phi, by solving the Poisson equation
     computePhi( rho, phi, beta, pc.self_fields_required_precision,
                 pc.self_fields_absolute_tolerance, pc.self_fields_max_iters,
-                pc.self_fields_verbosity );
+                pc.self_fields_verbosity, is_igf_2d_slices );
 
     // Compute the corresponding electric and magnetic field, from the potential phi
     computeE( Efield_fp, phi, beta );
@@ -161,7 +161,7 @@ void RelativisticExplicitES::AddBoundaryField (amrex::Vector<std::array< std::un
     // Compute the potential phi, by solving the Poisson equation
     computePhi( rho, phi, beta, self_fields_required_precision,
                 self_fields_absolute_tolerance, self_fields_max_iters,
-                self_fields_verbosity );
+                self_fields_verbosity, is_igf_2d_slices );
 
     // Compute the corresponding electric field, from the potential phi.
     computeE( Efield_fp, phi, beta );

diff --git a/Source/ablastr/fields/IntegratedGreenFunctionSolver.H b/Source/ablastr/fields/IntegratedGreenFunctionSolver.H
@@ -15,13 +15,15 @@
 #include <AMReX_REAL.H>
 #include <AMReX_Vector.H>
 
+#include <ablastr/constant.H>
 
 #include <array>
 #include <cmath>
 
 
 namespace ablastr::fields
 {
+    using namespace amrex::literals;
 
     /** @brief Implements equation 2 in https://doi.org/10.1103/PhysRevSTAB.10.129901
      *         with some modification to symmetrize the function.
@@ -30,14 +32,12 @@ namespace ablastr::fields
      * @param[in] y y-coordinate of given location
      * @param[in] z z-coordinate of given location
      *
-     * @return the integrated Green function G
+     * @return the integrated Green function G in 3D
      */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     amrex::Real
-    IntegratedPotential (amrex::Real x, amrex::Real y, amrex::Real z)
+    IntegratedPotential3D (amrex::Real x, amrex::Real y, amrex::Real z)
     {
-        using namespace amrex::literals;
-
         amrex::Real const r = std::sqrt( x*x + y*y + z*z );
         amrex::Real const G =
             - 0.5_rt * z*z * std::atan( x*y/(z*r) )
@@ -49,35 +49,79 @@ namespace ablastr::fields
         return G;
     }
 
+
     /** @brief add
      *
      * @param[in] x x-coordinate of given location
      * @param[in] y y-coordinate of given location
      * @param[in] z z-coordinate of given location
+     * @param[in] dx cell size along x
+     * @param[in] dy cell size along y
+     * @param[in] dz cell size along z
+     *
+     * @return the sum of integrated Green function G in 3D
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real
+    SumOfIntegratedPotential3D (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real dx, amrex::Real dy, amrex::Real dz)
+    {
+        amrex::Real const G_value = 1._rt/(4._rt*ablastr::constant::math::pi*ablastr::constant::SI::ep0) * (
+            IntegratedPotential3D( x+0.5_rt*dx, y+0.5_rt*dy, z+0.5_rt*dz )
+          - IntegratedPotential3D( x-0.5_rt*dx, y+0.5_rt*dy, z+0.5_rt*dz )
+          - IntegratedPotential3D( x+0.5_rt*dx, y-0.5_rt*dy, z+0.5_rt*dz )
+          + IntegratedPotential3D( x-0.5_rt*dx, y-0.5_rt*dy, z+0.5_rt*dz )
+          - IntegratedPotential3D( x+0.5_rt*dx, y+0.5_rt*dy, z-0.5_rt*dz )
+          + IntegratedPotential3D( x-0.5_rt*dx, y+0.5_rt*dy, z-0.5_rt*dz )
+          + IntegratedPotential3D( x+0.5_rt*dx, y-0.5_rt*dy, z-0.5_rt*dz )
+          - IntegratedPotential3D( x-0.5_rt*dx, y-0.5_rt*dy, z-0.5_rt*dz )
+        );
+        return G_value;
+    }
+
+
+    /** @brief Implements equation 2 in https://doi.org/10.1103/PhysRevSTAB.10.129901
+     *         with some modification to symmetrize the function.
      *
-     * @return the sum of integrated Green function G
+     * @param[in] x x-coordinate of given location
+     * @param[in] y y-coordinate of given location
+     *
+     * @return the integrated Green function G in 2D
      */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     amrex::Real
-    SumOfIntegratedPotential (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real dx, amrex::Real dy, amrex::Real dz)
+    IntegratedPotential2D (amrex::Real x, amrex::Real y)
     {
-        using namespace amrex::literals;
+        amrex::Real const G = 3_rt*x*y
+                              - x*x * std::atan(y/x)
+                              - y*y * std::atan(x/y)
+                              - x*y * std::log(x*x + y*y);
+        return G;
+    }
 
 
+    /** @brief add
+     *
+     * @param[in] x x-coordinate of given location
+     * @param[in] y y-coordinate of given location
+     * @param[in] dx cell size along x
+     * @param[in] dy cell size along y
+     *
+     * @return the sum of integrated Green function G in 2D
+     */
+    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    amrex::Real
+    SumOfIntegratedPotential2D (amrex::Real x, amrex::Real y, amrex::Real dx, amrex::Real dy)
+    {
         amrex::Real const G_value = 1._rt/(4._rt*ablastr::constant::math::pi*ablastr::constant::SI::ep0) * (
-            IntegratedPotential( x+0.5_rt*dx, y+0.5_rt*dy, z+0.5_rt*dz )
-          - IntegratedPotential( x-0.5_rt*dx, y+0.5_rt*dy, z+0.5_rt*dz )
-          - IntegratedPotential( x+0.5_rt*dx, y-0.5_rt*dy, z+0.5_rt*dz )
-          + IntegratedPotential( x-0.5_rt*dx, y-0.5_rt*dy, z+0.5_rt*dz )
-          - IntegratedPotential( x+0.5_rt*dx, y+0.5_rt*dy, z-0.5_rt*dz )
-          + IntegratedPotential( x-0.5_rt*dx, y+0.5_rt*dy, z-0.5_rt*dz )
-          + IntegratedPotential( x+0.5_rt*dx, y-0.5_rt*dy, z-0.5_rt*dz )
-          - IntegratedPotential( x-0.5_rt*dx, y-0.5_rt*dy, z-0.5_rt*dz )
+              IntegratedPotential2D( x+0.5_rt*dx, y+0.5_rt*dy )
+            - IntegratedPotential2D( x+0.5_rt*dx, y-0.5_rt*dy )
+            - IntegratedPotential2D( x-0.5_rt*dx, y+0.5_rt*dy )
+            + IntegratedPotential2D( x-0.5_rt*dx, y-0.5_rt*dy )
         );
-
         return G_value;
     }
 
+
     /** @brief Compute the electrostatic potential using the Integrated Green Function method
      *         as in http://dx.doi.org/10.1103/PhysRevSTAB.9.044204
      *
@@ -90,7 +134,8 @@ namespace ablastr::fields
     computePhiIGF (amrex::MultiFab const & rho,
                    amrex::MultiFab & phi,
                    std::array<amrex::Real, 3> const & cell_size,
-                   amrex::BoxArray const & ba);
+                   amrex::BoxArray const & ba,
+                   bool is_2d_slices);
 
 } // namespace ablastr::fields