diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst index 3787acbd639..310a7986abf 100644 --- a/Docs/source/usage/parameters.rst +++ b/Docs/source/usage/parameters.rst @@ -290,6 +290,14 @@ Overall simulation parameters In electromagnetic mode, this solver can be used to initialize the species' self fields (``.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): Select the type of Integrated Green Function solver. + If 0, solve Poisson equation in full 3D geometry. + If 1, solve Poisson equation in a quasi 3D geometry, neglecting the :math:`z` derivatives in the Laplacian of the Poisson equation. + In practice, in this case, the code performes many 2D Poisson solves on all :math:`(x,y)` slices, each slice at a given :math:`z`. + This is often a good approximation for ultra-relativistic beams propagating along the :math:`z` direction, with the relativistic solver. + As a consequence, this solver does not need to do an FFT along the :math:`z` direction, + and instead uses only transverse FFTs (along :math:`x` and :math:`y`) at each :math:`z` position (or :math:`z` "slice"). + * ``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 index c5ec4583da1..95a8d23687e 100644 --- a/Examples/Tests/open_bc_poisson_solver/CMakeLists.txt +++ b/Examples/Tests/open_bc_poisson_solver/CMakeLists.txt @@ -12,3 +12,15 @@ if(WarpX_FFT) 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/analysis.py b/Examples/Tests/open_bc_poisson_solver/analysis.py index 8ffd9ef52e2..25b55503cff 100755 --- a/Examples/Tests/open_bc_poisson_solver/analysis.py +++ b/Examples/Tests/open_bc_poisson_solver/analysis.py @@ -61,7 +61,6 @@ def evaluate_E(x, y, z): assert np.allclose(Ex_warpx, Ex_theory, rtol=0.032, atol=0) assert np.allclose(Ey_warpx, Ey_theory, rtol=0.029, atol=0) - # compare checksums evaluate_checksum( test_name=os.path.split(os.getcwd())[1], 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 new file mode 100644 index 00000000000..e2639c59e74 --- /dev/null +++ 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 new file mode 100644 index 00000000000..fd4a9afbc29 --- /dev/null +++ b/Regression/Checksum/benchmarks_json/test_3d_open_bc_poisson_solver_sliced.json @@ -0,0 +1,21 @@ +{ + "lev=0": { + "Bx": 100915933.44551668, + "By": 157610622.18551716, + "Bz": 2.598515299403035e-15, + "Ex": 4.725065270620093e+16, + "Ey": 3.0253948989229424e+16, + "Ez": 2787743.3330717986, + "rho": 10994013582437.193 + }, + "electron": { + "particle_momentum_x": 5.701277606056779e-19, + "particle_momentum_y": 3.650451663675671e-19, + "particle_momentum_z": 1.145432768297242e-10, + "particle_position_x": 17.314086912497864, + "particle_position_y": 0.25836912671877954, + "particle_position_z": 10066.329600000008, + "particle_weight": 19969036501.910976 + } +} + diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H old mode 100644 new mode 100755 index e58af394a7a..f57cfff6080 --- a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H +++ b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H @@ -92,7 +92,8 @@ public: amrex::Real required_precision, amrex::Real absolute_tolerance, int max_iters, - int verbosity + int verbosity, + bool is_igf_2d_slices ) const; /** @@ -153,6 +154,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 old mode 100644 new mode 100755 index 0b1dca675be..429e007b4d0 --- a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp @@ -38,7 +38,12 @@ void ElectrostaticSolver::ReadParameters () { 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_verbosity", self_fields_verbosity); + utils::parser::queryWithParser( + pp_warpx, "self_fields_verbosity", self_fields_verbosity); + + // FFT solver flags + utils::parser::queryWithParser( + pp_warpx, "use_2d_slices_fft_solver", is_igf_2d_slices); } void @@ -121,7 +126,9 @@ ElectrostaticSolver::computePhi ( Real const required_precision, Real absolute_tolerance, int const max_iters, - int const verbosity) const + int const verbosity, + bool const is_igf_2d +) const { using ablastr::fields::Direction; @@ -202,6 +209,7 @@ ElectrostaticSolver::computePhi ( warpx.boxArray(), WarpX::grid_type, is_solver_igf_on_lev0, + is_igf_2d, EB::enabled(), WarpX::do_single_precision_comms, warpx.refRatio(), diff --git a/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.H b/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.H old mode 100644 new mode 100755 diff --git a/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.cpp old mode 100644 new mode 100755 index e973ae66975..643efefb2f3 --- a/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/LabFrameExplicitES.cpp @@ -75,7 +75,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.H b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.H old mode 100644 new mode 100755 diff --git a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp old mode 100644 new mode 100755 index 69647da1702..0b1bcecd1e5 --- a/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp +++ b/Source/FieldSolver/ElectrostaticSolvers/RelativisticExplicitES.cpp @@ -130,7 +130,7 @@ void RelativisticExplicitES::AddSpaceChargeField ( computePhi( amrex::GetVecOfPtrs(rho), amrex::GetVecOfPtrs(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, amrex::GetVecOfPtrs(phi), beta ); @@ -168,7 +168,7 @@ void RelativisticExplicitES::AddBoundaryField (ablastr::fields::MultiLevelVector computePhi( amrex::GetVecOfPtrs(rho), amrex::GetVecOfPtrs(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, amrex::GetVecOfPtrs(phi), beta ); diff --git a/Source/ablastr/fields/IntegratedGreenFunctionSolver.H b/Source/ablastr/fields/IntegratedGreenFunctionSolver.H old mode 100644 new mode 100755 index 28885e167a3..9492cff885e --- a/Source/ablastr/fields/IntegratedGreenFunctionSolver.H +++ b/Source/ablastr/fields/IntegratedGreenFunctionSolver.H @@ -15,13 +15,14 @@ #include #include - #include #include 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,54 +31,92 @@ 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) ) - - 0.5_rt * y*y * std::atan( x*z/(y*r) ) - - 0.5_rt * x*x * std::atan( y*z/(x*r) ) - + y*z*std::asinh( x/std::sqrt(y*y + z*z) ) - + x*z*std::asinh( y/std::sqrt(x*x + z*z) ) - + x*y*std::asinh( z/std::sqrt(x*x + y*y) ); + amrex::Real const G = - 0.5_rt * z*z * std::atan( x*y/(z*r) ) + - 0.5_rt * y*y * std::atan( x*z/(y*r) ) + - 0.5_rt * x*x * std::atan( y*z/(x*r) ) + + y*z*std::asinh( x/std::sqrt(y*y + z*z) ) + + x*z*std::asinh( y/std::sqrt(x*x + z*z) ) + + x*y*std::asinh( z/std::sqrt(x*x + y*y) ); + return G; + } + + + /** @brief Implements equation 58 in https://doi.org/10.1016/j.jcp.2004.01.008 + * + * @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 + IntegratedPotential2D (amrex::Real x, amrex::Real y) + { + 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] 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 + * @return the sum of integrated Green function G in 3D */ 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) + SumOfIntegratedPotential3D (amrex::Real x, amrex::Real y, amrex::Real z, amrex::Real dx, amrex::Real dy, amrex::Real dz) { - using namespace amrex::literals; - + return 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 ) + ); + } - 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 ) - ); - return G_value; + /** @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) + { + return 1._rt/(4._rt*ablastr::constant::math::pi*ablastr::constant::SI::ep0) * ( + 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 ) + ); } + /** @brief Compute the electrostatic potential using the Integrated Green Function method * as in http://dx.doi.org/10.1103/PhysRevSTAB.9.044204 * @@ -90,7 +129,8 @@ namespace ablastr::fields computePhiIGF (amrex::MultiFab const & rho, amrex::MultiFab & phi, std::array const & cell_size, - amrex::BoxArray const & ba); + amrex::BoxArray const & ba, + bool is_igf_2d_slices); } // namespace ablastr::fields diff --git a/Source/ablastr/fields/IntegratedGreenFunctionSolver.cpp b/Source/ablastr/fields/IntegratedGreenFunctionSolver.cpp old mode 100644 new mode 100755 index b142978c8be..998bb179f5b --- a/Source/ablastr/fields/IntegratedGreenFunctionSolver.cpp +++ b/Source/ablastr/fields/IntegratedGreenFunctionSolver.cpp @@ -34,7 +34,8 @@ void computePhiIGF ( amrex::MultiFab const & rho, amrex::MultiFab & phi, std::array const & cell_size, - amrex::BoxArray const & ba) + amrex::BoxArray const & ba, + bool const is_igf_2d_slices) { using namespace amrex::literals; @@ -45,9 +46,6 @@ computePhiIGF ( amrex::MultiFab const & rho, domain.surroundingNodes(); // get nodal points, since `phi` and `rho` are nodal domain.grow( phi.nGrowVect() ); // include guard cells - // Do we grow the domain in the z-direction in the 2D mode? - bool const do_2d_fft = false; - int nprocs = amrex::ParallelDescriptor::NProcs(); { amrex::ParmParse pp("ablastr"); @@ -61,7 +59,7 @@ computePhiIGF ( amrex::MultiFab const & rho, } if (!obc_solver || obc_solver->Domain() != domain) { amrex::FFT::Info info{}; - if (do_2d_fft) { info.setBatchMode(true); } + if (is_igf_2d_slices) { info.setBatchMode(true); } // do 2D FFTs info.setNumProcs(nprocs); obc_solver = std::make_unique>(domain, info); } @@ -71,7 +69,9 @@ computePhiIGF ( amrex::MultiFab const & rho, amrex::Real const dy = cell_size[1]; amrex::Real const dz = cell_size[2]; - obc_solver->setGreensFunction( + if (!is_igf_2d_slices){ + // 2D sliced solver + obc_solver->setGreensFunction( [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::Real { int const i0 = i - lo[0]; @@ -80,9 +80,26 @@ computePhiIGF ( amrex::MultiFab const & rho, amrex::Real const x = i0*dx; amrex::Real const y = j0*dy; amrex::Real const z = k0*dz; - return SumOfIntegratedPotential(x, y, z, dx, dy, dz); + + return SumOfIntegratedPotential3D(x, y, z, dx, dy, dz); + }); + }else{ + // fully 3D solver + obc_solver->setGreensFunction( + [=] AMREX_GPU_DEVICE (int i, int j, int k) -> amrex::Real + { + int const i0 = i - lo[0]; + int const j0 = j - lo[1]; + amrex::Real const x = i0*dx; + amrex::Real const y = j0*dy; + amrex::ignore_unused(k); + + return SumOfIntegratedPotential2D(x, y, dx, dy); }); + } + obc_solver->solve(phi, rho); -} +} // computePhiIGF + } // namespace ablastr::fields diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H old mode 100644 new mode 100755 index aa9288fe950..1cc7d39e9b0 --- a/Source/ablastr/fields/PoissonSolver.H +++ b/Source/ablastr/fields/PoissonSolver.H @@ -206,6 +206,7 @@ computePhi ( amrex::Vector const& grids, utils::enums::GridType grid_type, bool is_solver_igf_on_lev0, + [[maybe_unused]] bool const is_igf_2d, bool eb_enabled = false, bool do_single_precision_comms = false, std::optional > rel_ref_ratio = std::nullopt, @@ -233,13 +234,13 @@ computePhi ( #endif #if !defined(ABLASTR_USE_FFT) - ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, - "Must compile with FFT support to use the IGF solver!"); + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, + "Must compile with FFT support to use the IGF solver!"); #endif #if !defined(WARPX_DIM_3D) - ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, - "The FFT Poisson solver is currently only implemented for 3D!"); + ABLASTR_ALWAYS_ASSERT_WITH_MESSAGE( !is_solver_igf_on_lev0, + "The FFT Poisson solver is currently only implemented for 3D!"); #endif // Set the value of beta @@ -270,7 +271,7 @@ computePhi ( if ( max_norm_b == 0 ) { phi[lev]->setVal(0); } else { - computePhiIGF( *rho[lev], *phi[lev], dx_scaled, grids[lev] ); + computePhiIGF( *rho[lev], *phi[lev], dx_scaled, grids[lev], is_igf_2d); } continue; }