diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index c1df197075b..b272aa27d23 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -241,6 +241,8 @@ Overall simulation parameters
         non-zero value is specified by the user via
         ``warpx.self_fields_absolute_tolerance``).
 
+    * ``gmres``: Poisson's equation is solved using an MLMG-preconditioned GMRES solver.
+
     * ``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``.
diff --git a/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.cpp b/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.cpp
index efcdac3404c..4f50f043848 100644
--- a/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.cpp
+++ b/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;
diff --git a/Source/Initialization/DivCleaner/ProjectionDivCleaner.cpp b/Source/Initialization/DivCleaner/ProjectionDivCleaner.cpp
index 660dbb753cb..6cffb1e72b4 100644
--- a/Source/Initialization/DivCleaner/ProjectionDivCleaner.cpp
+++ b/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
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 900a359b8b4..3f275cd381c 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/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";
diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H
index dea40d4b2df..a1e8f73b595 100644
--- a/Source/Utils/WarpXAlgorithmSelection.H
+++ b/Source/Utils/WarpXAlgorithmSelection.H
@@ -67,6 +67,7 @@ AMREX_ENUM(ElectrostaticSolverAlgo,
 
 AMREX_ENUM(PoissonSolverAlgo,
            Multigrid,
+           GMRES,
            IntegratedGreenFunction,
            fft = IntegratedGreenFunction,
            Default = Multigrid);
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index 0fa3656c9c0..b6b69612b89 100644
--- a/Source/WarpX.cpp
+++ b/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;
         }
diff --git a/Source/ablastr/fields/PoissonSolver.H b/Source/ablastr/fields/PoissonSolver.H
index 72034b667c0..1cff057e654 100755
--- a/Source/ablastr/fields/PoissonSolver.H
+++ b/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)