BLAST-WarpX · RemiLehe · Aug 27, 2026 · Aug 27, 2026 · Aug 27, 2026 · Aug 27, 2026
diff --git a/.azure-pipelines.yml b/.azure-pipelines.yml
@@ -54,7 +54,7 @@ jobs:
       # Cartesian 3D
       cartesian_3d:
         CDASH_BUILD_NAME: CPU-3D
-        WARPX_CMAKE_FLAGS: -DWarpX_DIMS=3 -DWarpX_FFT=ON -DWarpX_PYTHON=ON
+        WARPX_CMAKE_FLAGS: -DWarpX_DIMS=3 -DWarpX_FFT=ON -DWarpX_PYTHON=ON -DWarpX_PETSC=ON
       # Cylindrical RZ
       cylindrical_rz:
         CDASH_BUILD_NAME: CPU-RZ

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -450,6 +450,33 @@ Overall simulation parameters
             The extended simulation box size in real space is :math:`2n_x-1, 2n_y-1, 2n_z-1` with the 3D solver, :math:`2n_x-1, 2n_y -1, n_z` with the 2D solver.
             The extended simulation box size in spectral space is :math:`n_x, 2n_y-1, 2n_z-1` with the 3D solver, :math:`n_x, 2n_y-1, n_z` with the 2D solver.
 
+    * ``petsc``: Poisson's equation is solved by PETSc's GMRES, right-preconditioned by the AMReX multigrid.
+        This is the electrostatic counterpart of the PETSc interface that the implicit electromagnetic solvers
+        offer for the curl-curl equation (``newton.linear_solver = petsc_ksp``): the linear system is
+        handed to PETSc as a matrix-free operator, whose action, as well as that of the multigrid preconditioner,
+        is computed by AMReX. It therefore discretizes Poisson's equation exactly like ``multigrid`` does, and
+        accepts the same boundary conditions; only the outer iteration differs.
+        It requires the compilation flag ``-DWarpX_PETSC=ON``, and PETSc itself must be built with
+        CUDA (or HIP) support in order to run on GPUs.
+        It is not supported in ``labframe-effective-potential`` mode.
+        Note that in 1D with ``warpx.do_electrostatic = labframe``, and with the ``poissonsolver`` Python callback,
+        Poisson's equation is solved by a dedicated solver and this option has no effect.
+
+        Further customization of the Krylov solver (e.g. ``-ksp_type``, ``-ksp_gmres_restart``) is available
+        through PETSc's own runtime options (command line or option file), which take precedence over the
+        parameters below.
+
+          * ``petsc_poisson.use_mlmg_preconditioner`` (``bool``) optional (default: 1): Whether to use the AMReX
+            multigrid V-cycles as a preconditioner, through PETSc's ``PCShell`` interface. Setting this to 0 runs
+            the Krylov solver unpreconditioned, which is much slower and mostly useful for debugging.
+
+          * ``petsc_poisson.precond_num_iters`` (``int``) optional (default: 1): The number of multigrid V-cycles
+            per application of the preconditioner.
+
+          * ``petsc_poisson.verbose`` (``int``) optional (default: value of
+            :pp:param:`warpx.self_fields_verbosity`): 0 is silent, 1 prints the exit status of every solve,
+            2 additionally prints the residual at every Krylov iteration.
+
 .. pp:param:: warpx.self_fields_required_precision
     :type: ``float``
     :default: 1.e-11

diff --git a/Examples/Tests/electrostatic_sphere/CMakeLists.txt b/Examples/Tests/electrostatic_sphere/CMakeLists.txt
@@ -21,6 +21,21 @@ add_warpx_test(
     OFF  # dependency
 )
 
+if(WarpX_PETSC)
+    # No checksum: the reference values still have to be generated on a build
+    # with PETSc. The physics analysis below already checks the solution against
+    # the analytical solution of the expanding sphere.
+    add_warpx_test(
+        test_3d_electrostatic_sphere_lab_frame_petsc  # name
+        3  # dims
+        2  # nprocs
+        inputs_test_3d_electrostatic_sphere_lab_frame_petsc  # inputs
+        "analysis_electrostatic_sphere.py diags/diag1000030"  # analysis
+        OFF  # checksum
+        OFF  # dependency
+    )
+endif()
+
 add_warpx_test(
     test_3d_electrostatic_sphere_lab_frame_mr_emass_10  # name
     3  # dims

diff --git a/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_petsc b/Examples/Tests/electrostatic_sphere/inputs_test_3d_electrostatic_sphere_lab_frame_petsc
@@ -0,0 +1,10 @@
+# base input parameters
+FILE = inputs_base_3d
+
+# test input parameters
+diag2.electron.variables = x y z ux uy uz w phi
+warpx.do_electrostatic = labframe
+
+# solve the Poisson equation with PETSc's GMRES, preconditioned by the AMReX
+# multigrid, instead of using the multigrid directly
+warpx.poisson_solver = petsc
diff --git a/Examples/Tests/electrostatic_sphere_eb/CMakeLists.txt b/Examples/Tests/electrostatic_sphere_eb/CMakeLists.txt
@@ -25,6 +25,21 @@ if(WarpX_EB)
     )
 endif()
 
+if(WarpX_EB AND WarpX_PETSC)
+    # No checksum: the reference values still have to be generated on a build
+    # with PETSc. Same setup as test_3d_electrostatic_sphere_eb_mixed_bc, but
+    # solved with PETSc, so the two can be compared against each other.
+    add_warpx_test(
+        test_3d_electrostatic_sphere_eb_petsc  # name
+        3  # dims
+        2  # nprocs
+        inputs_test_3d_electrostatic_sphere_eb_petsc  # inputs
+        OFF  # analysis
+        OFF  # checksum
+        OFF  # dependency
+    )
+endif()
+
 if(WarpX_EB)
     add_warpx_test(
         test_3d_electrostatic_sphere_eb_picmi  # name

diff --git a/Examples/Tests/electrostatic_sphere_eb/inputs_test_3d_electrostatic_sphere_eb_petsc b/Examples/Tests/electrostatic_sphere_eb/inputs_test_3d_electrostatic_sphere_eb_petsc
@@ -0,0 +1,7 @@
+# base input parameters: same physics as the mixed-BC test (embedded boundary
+# with a non-zero potential, mixed PEC/Neumann domain boundaries)
+FILE = inputs_test_3d_electrostatic_sphere_eb_mixed_bc
+
+# test input parameters: solve the Poisson equation with PETSc's GMRES,
+# preconditioned by the AMReX multigrid, instead of using the multigrid directly
+warpx.poisson_solver = petsc
diff --git a/Python/pywarpx/PETScPoisson.py b/Python/pywarpx/PETScPoisson.py
@@ -0,0 +1,11 @@
+# Copyright 2026 The WarpX Community
+#
+# This file is part of WarpX.
+#
+# License: BSD-3-Clause-LBNL
+
+from .Bucket import Bucket
+
+# Options of the PETSc Krylov solver that is used for the Poisson equation
+# when warpx.poisson_solver = petsc (requires compiling with -DWarpX_PETSC=ON)
+petsc_poisson = Bucket("petsc_poisson")
diff --git a/Python/pywarpx/WarpX.py b/Python/pywarpx/WarpX.py
@@ -23,6 +23,7 @@
 from .Interpolation import interpolation
 from .Lasers import lasers, lasers_list
 from .Particles import particles, particles_list
+from .PETScPoisson import petsc_poisson
 from .PSATD import psatd
 
 
@@ -55,6 +56,7 @@ def create_argv_list(self, **kw):
         argv += interpolation.attrlist()
         argv += psatd.attrlist()
         argv += eb2.attrlist()
+        argv += petsc_poisson.attrlist()
 
         argv += particles.attrlist()
         for particle in particles_list:
@@ -223,6 +225,7 @@ def finalize(self, finalize_mpi=1):
             lasers,
             my_constants,
             particles,
+            petsc_poisson,
             psatd,
             reduced_diagnostics,
             self,

diff --git a/Python/pywarpx/__init__.py b/Python/pywarpx/__init__.py
@@ -38,6 +38,7 @@
 from .Lasers import lasers  # noqa
 from .LoadThirdParty import load_cupy  # noqa
 from .Particles import new_species, particles  # noqa
+from .PETScPoisson import petsc_poisson  # noqa
 from .PSATD import psatd  # noqa
 from .WarpX import warpx  # noqa
 

diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.H
@@ -13,6 +13,7 @@
 #include "Fluids/MultiFluidContainer.H"
 #include "Particles/MultiParticleContainer.H"
 
+#include <ablastr/fields/PETScPoissonSolver.H>
 #include <ablastr/profiler/ProfilerWrapper.H>
 #include <AMReX_Array.H>
 
@@ -159,6 +160,9 @@ public:
     /** MLGM number of smoothing sweeps */
     int self_fields_num_final_sweeps = 8;
 
+    /** Settings of the optional PETSc Krylov solver for the Poisson equation */
+    ablastr::fields::PETScPoissonOptions m_petsc_options;
+
     /** Parameters for FFT Poisson solver aka IGF */
     // 0: full 3D, 1: many 2D z-slices (quasi-3D)
     bool is_igf_2d_slices = false;

diff --git a/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp b/Source/FieldSolver/ElectrostaticSolvers/ElectrostaticSolver.cpp
@@ -47,6 +47,20 @@ void ElectrostaticSolver::ReadParameters () {
             self_fields_num_final_sweeps > 0,
             "warpx.self_fields_num_final_sweeps must be > 0");
     }
+
+    // With warpx.poisson_solver = petsc, the Poisson equation is solved by one
+    // of PETSc's Krylov solvers, preconditioned by MLMG
+    m_petsc_options.use_petsc_ksp =
+        (WarpX::poisson_solver_id == PoissonSolverAlgo::PETSc);
+    if (m_petsc_options.use_petsc_ksp) {
+        ParmParse const pp_petsc("petsc_poisson");
+        pp_petsc.query("use_mlmg_preconditioner", m_petsc_options.use_mlmg_preconditioner);
+        utils::parser::queryWithParser(
+            pp_petsc, "precond_num_iters", m_petsc_options.precond_num_iters);
+        m_petsc_options.verbosity = self_fields_verbosity;
+        utils::parser::queryWithParser(pp_petsc, "verbose", m_petsc_options.verbosity);
+    }
+
     // FFT solver flags
     utils::parser::queryWithParser(
         pp_warpx, "use_2d_slices_fft_solver", is_igf_2d_slices);
@@ -221,6 +235,7 @@ ElectrostaticSolver::computePhi (
         WarpX::do_single_precision_comms,
         warpx.refRatio(),
         self_fields_num_final_sweeps,
+        m_petsc_options,
         post_phi_calculation,
         *m_poisson_boundary_handler,
         warpx.gett_new(0),

diff --git a/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.H b/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.H
@@ -14,7 +14,9 @@
 #include <AMReX_Array.H>
 #include <AMReX_Geometry.H>
 #include <AMReX_MultiFab.H>
+#include <AMReX_MLLinOp.H>
 #include <AMReX_MLMG.H>
+#include <AMReX_MLNodeLinOp.H>
 #include <AMReX_REAL.H>
 #include <AMReX_Parser.H>
 
@@ -136,6 +138,31 @@ class EBCalcEfromPhiPerLevel {
             field->mult(-1._rt);
         }
     }
+
+    /** \brief Same, but computed directly from the linear operator and `phi`
+     *
+     * `amrex::MLMG::getGradSolution` takes the potential from the internal state
+     * that `amrex::MLMG::solve` leaves behind, so it can only be used right after
+     * an MLMG solve. This overload is for the solvers that compute `phi` without
+     * going through `amrex::MLMG::solve`, such as the PETSc Poisson solver.
+     *
+     * \param[in] linop the linear operator of the Poisson solve
+     * \param[inout] phi the potential; only its ghost nodes are modified
+     * \param[in] lev the mesh refinement level
+     */
+    void operator()(amrex::MLNodeLinOp & linop, amrex::MultiFab & phi, int const lev) {
+        using namespace amrex::literals;
+
+        // Fill the ghost nodes of `phi`, in particular the ones that lie outside
+        // of a Neumann boundary, which `compGrad` reads. `amrex::MLMG::solve`
+        // does this internally on the potential that it hands to `compGrad`.
+        linop.applyBC(0, 0, phi, amrex::MLNodeLinOp::BCMode::Inhomogeneous,
+                      amrex::MLNodeLinOp::StateMode::Solution);
+        linop.compGrad(0, m_e_field[lev], phi, amrex::MLNodeLinOp::Location::FaceCenter);
+        for (auto &field: m_e_field[lev]) {
+            field->mult(-1._rt);
+        }
+    }
 };
 
 #endif // WARPX_BOUNDARYHANDLER_H_
diff --git a/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.cpp b/Source/FieldSolver/ElectrostaticSolvers/PoissonBoundaryHandler.cpp
@@ -81,7 +81,9 @@ 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){
+        // The grid-based solvers (multigrid, and PETSc preconditioned by it)
+        // share the same boundary conditions
+        if (WarpX::poisson_solver_id != PoissonSolverAlgo::IntegratedGreenFunction){
             if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::Periodic
                     && WarpX::field_boundary_hi[idim] == FieldBoundaryType::Periodic ) {
                 lobc[idim] = LinOpBCType::Periodic;
@@ -102,7 +104,7 @@ void PoissonBoundaryHandler::DefinePhiBCs (const amrex::Geometry& geom)
                 else {
                     WARPX_ABORT_WITH_MESSAGE(
                         "Field boundary conditions have to be either periodic, PEC or neumann "
-                        "when using the electrostatic Multigrid solver,  but they are " + amrex::getEnumNameString(WarpX::field_boundary_lo[idim])
+                        "when using the electrostatic Multigrid or PETSc solver,  but they are " + amrex::getEnumNameString(WarpX::field_boundary_lo[idim])
                     );
                 }
 
@@ -117,7 +119,7 @@ void PoissonBoundaryHandler::DefinePhiBCs (const amrex::Geometry& geom)
                 else {
                     WARPX_ABORT_WITH_MESSAGE(
                         "Field boundary conditions have to be either periodic, PEC or neumann "
-                        "when using the electrostatic Multigrid solver,  but they are " + amrex::getEnumNameString(WarpX::field_boundary_hi[idim])
+                        "when using the electrostatic Multigrid or PETSc solver,  but they are " + amrex::getEnumNameString(WarpX::field_boundary_hi[idim])
                     );
                 }
             }

diff --git a/Source/Initialization/DivCleaner/ProjectionDivCleaner.cpp b/Source/Initialization/DivCleaner/ProjectionDivCleaner.cpp
@@ -419,7 +419,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::IntegratedGreenFunction))
 #if defined(WARPX_DIM_RZ)
                 && WarpX::grid_type == GridType::Staggered
 #endif

diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
@@ -591,6 +591,9 @@ WarpX::PrintMainPICparameters ()
         else if(poisson_solver_id == PoissonSolverAlgo::Multigrid){
             amrex::Print() << "Poisson solver:       | multigrid" << "\n";
         }
+        else if(poisson_solver_id == PoissonSolverAlgo::PETSc){
+            amrex::Print() << "Poisson solver:       | PETSc (multigrid preconditioner)" << "\n";
+        }
     }
 
     amrex::Print() << "-------------------------------------------------------------------------------\n";

diff --git a/Source/Utils/WarpXAlgorithmSelection.H b/Source/Utils/WarpXAlgorithmSelection.H
@@ -70,6 +70,7 @@ AMREX_ENUM(PoissonSolverAlgo,
            Multigrid,
            IntegratedGreenFunction,
            fft = IntegratedGreenFunction,
+           PETSc,
            Default = Multigrid);
 
 AMREX_ENUM(ParticlePusherAlgo,

diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
@@ -802,6 +802,24 @@ WarpX::ReadParameters ()
         poisson_solver_id!=PoissonSolverAlgo::IntegratedGreenFunction,
         "To use the FFT Poisson solver, compile with -DWarpX_FFT=ON.");
 #endif
+
+#ifndef AMREX_USE_PETSC
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        poisson_solver_id!=PoissonSolverAlgo::PETSc,
+        "To use the PETSc Poisson solver, compile with -DWarpX_PETSC=ON.");
+#endif
+
+        // The effective-potential solver does not go through the same Poisson
+        // solve as the other electrostatic solvers, so it would silently ignore
+        // the PETSc solver.
+        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        (
+            electrostatic_solver_id!=ElectrostaticSolverAlgo::LabFrameEffectivePotential ||
+            poisson_solver_id!=PoissonSolverAlgo::PETSc
+        ),
+        "The PETSc Poisson solver is not implemented in labframe-effective-potential mode yet."
+        );
+
         // Read magnetostatic solver parameters
         // First use self_fields_* as defaults for backward compatibility,
         // then allow explicit magnetostatic_solver_* parameters to override
@@ -1192,7 +1210,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::IntegratedGreenFunction)))
         {
             m_do_initial_div_cleaning = true;
         }

diff --git a/Source/ablastr/fields/CMakeLists.txt b/Source/ablastr/fields/CMakeLists.txt
@@ -4,6 +4,7 @@ foreach(D IN LISTS WarpX_DIMS)
     target_sources(ablastr_${SD}
       PRIVATE
         MultiFabRegister.cpp
+        PETScPoissonSolver.cpp
     )
 
     if(ABLASTR_FFT AND D EQUAL 3)

diff --git a/Source/ablastr/fields/Make.package b/Source/ablastr/fields/Make.package
@@ -1,5 +1,6 @@
 
 CEXE_sources += MultiFabRegister.cpp
+CEXE_sources += PETScPoissonSolver.cpp
 
 ifeq ($(USE_FFT),TRUE)
     ifeq ($(DIM),3)