GMRES for Poisson Equation

Docs/source/usage/parameters.rst

@@ -241,6 +241,11 @@ Overall simulation parameters
         non-zero value is specified by the user via
         ``warpx.self_fields_absolute_tolerance``).
 
+    * ``gmres``: Poisson's equation is solved using a standard GMRES solver.
+        The GMRES algorithm is the classic right-preconditioned variation presented in the original GMRES
+        paper in `Y. Saad et al. (SIAM J. Sci. Stat. Comput., 7, 3, 1986) <https://doi.org/10.1137/0907058>`__.
+        The preconditioner is 1 cycle of the AMReX-based MLMG solvers.
+
     * ``fft``: Poisson's equation is solved using an Integrated Green Function method (which requires FFT calculations).
         See these references for more details :cite:t:`param-QiangPhysRevSTAB2006`, :cite:t:`param-QiangPhysRevSTAB2006err`.
         It only works in 3D and it requires the compilation flag ``-DWarpX_FFT=ON``.

Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.cpp

@@ -81,7 +81,7 @@ void PoissonBoundaryHandler::DefinePhiBCs (const amrex::Geometry& geom)
     amrex::ignore_unused(geom);
 #endif
     for (int idim=dim_start; idim<AMREX_SPACEDIM; idim++){
-        if (WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid){
+        if (WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid || WarpX::poisson_solver_id == PoissonSolverAlgo::GMRES){
             if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::Periodic
                     && WarpX::field_boundary_hi[idim] == FieldBoundaryType::Periodic ) {
                 lobc[idim] = LinOpBCType::Periodic;

Source/Initialization/DivCleaner/ProjectionDivCleaner.cpp

@@ -401,7 +401,7 @@ WarpX::ProjectionCleanDivB() {
             ||  WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC
             ||  ( (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame
                 || WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic)
-                && WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid))
+                && (WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid || WarpX::poisson_solver_id == PoissonSolverAlgo::GMRES)))
 #if defined(WARPX_DIM_RZ)
                 && WarpX::grid_type == GridType::Staggered
 #endif

Source/Initialization/WarpXInitData.cpp

@@ -584,6 +584,9 @@ WarpX::PrintMainPICparameters ()
         else if(poisson_solver_id == PoissonSolverAlgo::Multigrid){
             amrex::Print() << "Poisson solver:       | multigrid" << "\n";
         }
+        else if(poisson_solver_id == PoissonSolverAlgo::GMRES){
+            amrex::Print() << "Poisson solver:       | gmres" << "\n";
+        }
     }
 
     amrex::Print() << "-------------------------------------------------------------------------------\n";

Source/Utils/WarpXAlgorithmSelection.H

@@ -67,6 +67,7 @@ AMREX_ENUM(ElectrostaticSolverAlgo,
 
 AMREX_ENUM(PoissonSolverAlgo,
            Multigrid,
+           GMRES,
            IntegratedGreenFunction,
            fft = IntegratedGreenFunction,
            Default = Multigrid);

Source/WarpX.cpp

@@ -1120,7 +1120,7 @@ WarpX::ReadParameters ()
             ||  WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC
             ||  ( (WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrame
                 || WarpX::electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic)
-                && WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid)))
+                && (WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid || WarpX::poisson_solver_id == PoissonSolverAlgo::GMRES))))
         {
             m_do_initial_div_cleaning = true;
         }

Source/ablastr/fields/PoissonSolver.H

@@ -51,6 +51,8 @@
 #ifdef AMREX_USE_EB
 #   include <AMReX_EBFabFactory.H>
 #endif
+#include <AMReX_GMRES_MLMG.H>
+#include "WarpX.H"
 
 #include <algorithm>
 #include <array>
@@ -414,8 +416,14 @@ computePhi (
         }
 
         // Solve Poisson equation at lev
-        mlmg.solve( {phi[lev]}, {rho[lev]},
-                     relative_tolerance, absolute_tolerance );
+        if (WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid) {
+            mlmg.solve( {phi[lev]}, {rho[lev]},
+                        relative_tolerance, absolute_tolerance );
+        } else {
+            amrex::GMRESMLMG gmsolve(mlmg);
+            gmsolve.solve( *phi[lev], *rho[lev], relative_tolerance, absolute_tolerance );
+            linop->postSolve({phi[lev]});
+        }
 
         const amrex::IntVect& refratio = rel_ref_ratio.value()[lev];
         const int ncomp = linop->getNComp();
@@ -435,12 +443,47 @@ computePhi (
                                         ng);
         }
 
-        // Run additional operations, such as calculation of the E field for embedded boundaries
-        if constexpr (!std::is_same_v<T_PostPhiCalculationFunctor, std::nullopt_t>) {
-            if (post_phi_calculation.has_value()) {
-                post_phi_calculation.value()(mlmg, lev);
+        if (WarpX::poisson_solver_id == PoissonSolverAlgo::Multigrid) {
+            // Run additional operations, such as calculation of the E field for embedded boundaries
+            if constexpr (!std::is_same_v<T_PostPhiCalculationFunctor, std::nullopt_t>) {
+                if (post_phi_calculation.has_value()) {
+                    post_phi_calculation.value()(mlmg, lev);
+                }
             }
         }
+
+        if (WarpX::poisson_solver_id == PoissonSolverAlgo::GMRES && eb_enabled) {
+
+            // create an alieas to the WarpX Efield
+            auto & warpx = WarpX::GetInstance();
+
+            amrex::Vector< amrex::Array<amrex::MultiFab *, AMREX_SPACEDIM> > e_field;
+            e_field.push_back(
+#if defined(WARPX_DIM_1D_Z)
+                amrex::Array<amrex::MultiFab*, 1>{
+                    warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{2}, lev)
+                }
+#elif defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+                amrex::Array<amrex::MultiFab*, 2>{
+                    warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{0}, lev),
+                    warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{2}, lev)
+                        }
+#elif defined(WARPX_DIM_3D)
+                amrex::Array<amrex::MultiFab *, 3>{
+                    warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{0}, lev),
+                    warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{1}, lev),
+                    warpx.m_fields.get(warpx::fields::FieldType::Efield_fp, Direction{2}, lev)
+                        }
+#endif
+                );
+
+            // compute E = grad(phi)
+            linop->compGrad(lev, e_field[lev], *phi[lev], amrex::MLLinOp::Location::CellCenter);
+            for(int dir=0; dir<AMREX_SPACEDIM; ++dir) {
+                e_field[lev][dir]->mult(-1.);
+            }
+        }
+
         rho[lev]->mult(-ablastr::constant::SI::epsilon_0);  // Multiply rho by epsilon again
 
     } // loop over lev(els)