diff --git a/Python/pywarpx/fields.py b/Python/pywarpx/fields.py
index 97110ce697d..10850e97870 100644
--- a/Python/pywarpx/fields.py
+++ b/Python/pywarpx/fields.py
@@ -737,7 +737,7 @@ def GFPWrapper(level=0, include_ghosts=False):
 
 def AxFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
-        mf_name="Afield_fp_nodal[x]",
+        mf_name="Afield_fp[x]",
         level=level,
         include_ghosts=include_ghosts,
     )
@@ -745,7 +745,7 @@ def AxFPWrapper(level=0, include_ghosts=False):
 
 def AyFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
-        mf_name="Afield_fp_nodal[y]",
+        mf_name="Afield_fp[y]",
         level=level,
         include_ghosts=include_ghosts,
     )
@@ -753,7 +753,7 @@ def AyFPWrapper(level=0, include_ghosts=False):
 
 def AzFPWrapper(level=0, include_ghosts=False):
     return _MultiFABWrapper(
-        mf_name="Afield_fp_nodal[z]",
+        mf_name="Afield_fp[z]",
         level=level,
         include_ghosts=include_ghosts,
     )
diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp
index 7bd93b572c6..0346f0f226d 100644
--- a/Source/Diagnostics/FullDiagnostics.cpp
+++ b/Source/Diagnostics/FullDiagnostics.cpp
@@ -671,7 +671,7 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
         } else if ( m_varnames[comp] == "jz_displacement" ) {
                 m_all_field_functors[lev][comp] = std::make_unique<JdispFunctor>(2, lev, m_crse_ratio, true);
         } else if ( m_varnames[comp] == "Az" ){
-            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 2), lev, m_crse_ratio);
+            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 2), lev, m_crse_ratio);
         } else if ( m_varnames[comp] == "rho" ){
             // Initialize rho functor to dump total rho
             m_all_field_functors[lev][comp] = std::make_unique<RhoFunctor>(lev, m_crse_ratio, true);
@@ -720,9 +720,9 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
             } else if  (m_varnames[comp] == "jt_displacement" ){
                 m_all_field_functors[lev][comp] = std::make_unique<JdispFunctor>(1, lev, m_crse_ratio, true);
             } else if ( m_varnames[comp] == "Ar" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 0), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 0), lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "At" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 1), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 1), lev, m_crse_ratio);
             } else {
                 WARPX_ABORT_WITH_MESSAGE(m_varnames[comp] + " is not a known field output type for RZ geometry");
             }
@@ -747,9 +747,9 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
             } else if ( m_varnames[comp] == "jy_displacement" ){
                 m_all_field_functors[lev][comp] = std::make_unique<JdispFunctor>(1, lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "Ax" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 0), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 0), lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "Ay" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::vector_potential_fp, lev, 1), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.getFieldPointer(FieldType::Afield_fp, lev, 1), lev, m_crse_ratio);
             } else {
                 std::cout << "Error on component " << m_varnames[comp] << std::endl;
                 WARPX_ABORT_WITH_MESSAGE(m_varnames[comp] + " is not a known field output type for this geometry");
diff --git a/Source/FieldSolver/Fields.H b/Source/FieldSolver/Fields.H
index 7c9cfe5285e..d26df49f3a1 100644
--- a/Source/FieldSolver/Fields.H
+++ b/Source/FieldSolver/Fields.H
@@ -23,7 +23,7 @@ namespace warpx::fields
         F_fp,
         G_fp,
         phi_fp,
-        vector_potential_fp,
+        Afield_fp,
         Efield_cp,
         Bfield_cp,
         current_cp,
diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H
index cbde5cdd56c..8e4c593ea7d 100644
--- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H
+++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H
@@ -2,34 +2,169 @@
  *
  * This file is part of WarpX.
  *
- * Authors: S. Eric Clark (LLNL)
+ * Authors: S. Eric Clark (LLNL), Roelof Groenewald (TAE)
  *
  * License: BSD-3-Clause-LBNL
  */
 #ifndef WARPX_MAGNETOSTATICSOLVER_H_
 #define WARPX_MAGNETOSTATICSOLVER_H_
 
+#include "MagnetostaticSolver_fwd.H"
+
+#include "Utils/Parser/ParserUtils.H"
+#include "Utils/TextMsg.H"
+#include "Utils/WarpXAlgorithmSelection.H"
+
 #include <AMReX_Array.H>
 #include <AMReX_MultiFab.H>
 #include <AMReX_MLMG.H>
 #include <AMReX_REAL.H>
 
+using namespace amrex;
+
+
+/**
+ * \brief This class contains functionality involved in solving the self
+ * magnetic field of a plasma due to internal currents, i.e., the magnetostatic
+ * solver. The solver uses the plasma current to calculate the vector potential
+ * (in the Coulomb gauge) and evaluates B = ∇ x A.
+ */
+class MagnetostaticSolver
+{
+public:
+    MagnetostaticSolver (int nlevs_max); // constructor
+
+    /** Read user-defined model parameters. Called in constructor. */
+    void ReadParameters ();
+
+    /** Allocate member multifabs. Called in constructor. */
+    void AllocateMFs (int nlevs_max);
+    void AllocateLevelMFs (int lev, const amrex::BoxArray& ba, const amrex::DistributionMapping& dm,
+                           int ncomps, const amrex::IntVect& ngJ, const amrex::IntVect& ngRho,
+                           const amrex::IntVect& jx_nodal_flag, const amrex::IntVect& jy_nodal_flag,
+                           const amrex::IntVect& jz_nodal_flag, const amrex::IntVect& rho_nodal_flag);
+
+    /** Helper function to clear values from member multifabs. */
+    void ClearLevel (int lev);
 
-namespace MagnetostaticSolver {
+    void InitData ();
 
     /** Boundary Handler for the Vector Potential Poisson Solver
      * This only will handle homogeneous Dirichlet boundary conditions on
-     * embedded boundaries, and homogeneous dirichlet/Neumann or periodic boundary conditions
+     * embedded boundaries, and homogeneous dirichlet/Neumann or periodic
+     * boundary conditions on domain boundaries.
      */
-    class VectorPoissonBoundaryHandler {
-    public:
+    struct BoundaryHandler {
         amrex::Array<amrex::Array<amrex::LinOpBCType, AMREX_SPACEDIM>, 3> lobc, hibc;
         bool bcs_set = false;
         std::array<std::array<bool, AMREX_SPACEDIM * 2>, 3> dirichlet_flag;
         bool has_non_periodic = false;
 
-        void defineVectorPotentialBCs ();
+        void defineVectorPotentialBCs (
+            Vector<FieldBoundaryType> field_boundary_lo,
+            Vector<FieldBoundaryType> field_boundary_hi,
+            [[maybe_unused]] Real prob_lo
+        ) {
+            for (int adim = 0; adim < 3; adim++) {
+#ifdef WARPX_DIM_RZ
+                int dim_start = 0;
+                if (prob_lo == 0){
+                    lobc[adim][0] = LinOpBCType::Neumann;
+                    dirichlet_flag[adim][0] = false;
+                    dim_start = 1;
+
+                    // handle the r_max boundary explicitly
+                    if (field_boundary_hi[0] == FieldBoundaryType::PEC) {
+                        if (adim == 0) {
+                            hibc[adim][0] = LinOpBCType::Neumann;
+                            dirichlet_flag[adim][1] = false;
+                        } else{
+                            hibc[adim][0] = LinOpBCType::Dirichlet;
+                            dirichlet_flag[adim][1] = true;
+                        }
+                    }
+                    else if (field_boundary_hi[0] == FieldBoundaryType::Neumann) {
+                        hibc[adim][0] = LinOpBCType::Neumann;
+                        dirichlet_flag[adim][1] = false;
+                    }
+                }
+#else
+                const int dim_start = 0;
+#endif
+                bool ndotA = false;
+                for (int idim=dim_start; idim<AMREX_SPACEDIM; idim++){
+                    ndotA = (adim == idim);
+
+#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
+                    if (idim == 1) { ndotA = (adim == 2); }
+#endif
+
+                    if ( field_boundary_lo[idim] == FieldBoundaryType::Periodic
+                        && field_boundary_hi[idim] == FieldBoundaryType::Periodic ) {
+                        lobc[adim][idim] = LinOpBCType::Periodic;
+                        hibc[adim][idim] = LinOpBCType::Periodic;
+                        dirichlet_flag[adim][idim*2] = false;
+                        dirichlet_flag[adim][idim*2+1] = false;
+                    }
+                    else {
+                        has_non_periodic = true;
+                        if ( field_boundary_lo[idim] == FieldBoundaryType::PEC ) {
+                            if (ndotA) {
+                                lobc[adim][idim] = LinOpBCType::Neumann;
+                                dirichlet_flag[adim][idim*2] = false;
+                            } else {
+                                lobc[adim][idim] = LinOpBCType::Dirichlet;
+                                dirichlet_flag[adim][idim*2] = true;
+                            }
+                        }
+
+                        else if ( field_boundary_lo[idim] == FieldBoundaryType::Neumann ) {
+                            lobc[adim][idim] = LinOpBCType::Neumann;
+                            dirichlet_flag[adim][idim*2] = false;
+                        }
+                        else {
+                            WARPX_ABORT_WITH_MESSAGE(
+                                "Field boundary conditions have to be either periodic, PEC or neumann "
+                                "when using the magnetostatic solver."
+                            );
+                        }
+
+                        if ( field_boundary_hi[idim] == FieldBoundaryType::PEC ) {
+                            if (ndotA) {
+                                hibc[adim][idim] = LinOpBCType::Neumann;
+                                dirichlet_flag[adim][idim*2+1] = false;
+                            } else {
+                                hibc[adim][idim] = LinOpBCType::Dirichlet;
+                                dirichlet_flag[adim][idim*2+1] = true;
+                            }
+                        }
+                        else if ( field_boundary_hi[idim] == FieldBoundaryType::Neumann ) {
+                            hibc[adim][idim] = LinOpBCType::Neumann;
+                            dirichlet_flag[adim][idim*2+1] = false;
+                        }
+                        else {
+                            WARPX_ABORT_WITH_MESSAGE(
+                                "Field boundary conditions have to be either periodic, PEC or neumann "
+                                "when using the magnetostatic solver."
+                            );
+                        }
+                    }
+                }
+            }
+            bcs_set = true;
+        }
     };
-} // namespace MagnetostaticSolver
 
+    BoundaryHandler m_boundary_handler = BoundaryHandler();
+
+    // numerical parameters
+    amrex::Real required_precision = 1.e-11_rt;
+    amrex::Real absolute_tolerance = 0.0_rt;
+    int max_iters = 200;
+    int verbosity = 2;
+
+    // Declare multifabs specifically needed for the magnetostatic solver
+    Vector<std::array< std::unique_ptr<MultiFab>, 3 > > Afield_fp_nodal;
+    Vector<std::array< std::unique_ptr<MultiFab>, 3 > > current_fp_temp;
+};
 #endif //WARPX_MAGNETOSTATICSOLVER_H_
diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp
index 682edcf2df7..4732974d134 100644
--- a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp
+++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver.cpp
@@ -6,17 +6,15 @@
  *
  * License: BSD-3-Clause-LBNL
  */
-#include "WarpX.H"
-
 #include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H"
 #include "Parallelization/GuardCellManager.H"
 #include "Particles/MultiParticleContainer.H"
 #include "Python/callbacks.H"
 #include "Utils/WarpXConst.H"
-#include "Utils/TextMsg.H"
 #include "Utils/WarpXUtil.H"
 #include "Utils/WarpXProfilerWrapper.H"
 #include "Parallelization/WarpXComm_K.H"
+#include "WarpX.H"
 
 #include <ablastr/utils/Communication.H>
 #include <ablastr/warn_manager/WarnManager.H>
@@ -87,7 +85,7 @@ WarpX::ComputeMagnetostaticFieldLabFrame()
     // Perform current deposition at t_{n+1/2} (source of Poisson solver).
     mypc->DepositCurrent(current_fp, dt[0], -0.5_rt * dt[0]);
 
-    // Synchronize J and rho:
+    // Synchronize J:
     // filter (if used), exchange guard cells, interpolate across MR levels
     // and apply boundary conditions
     SyncCurrent(current_fp, current_cp, current_buf);
@@ -97,12 +95,6 @@ WarpX::ComputeMagnetostaticFieldLabFrame()
             lev, current_fp[lev][0].get(), current_fp[lev][1].get(),
             current_fp[lev][2].get(), PatchType::fine
         );
-        if (lev > 0) {
-            ApplyJfieldBoundary(
-                lev, current_cp[lev][0].get(), current_cp[lev][1].get(),
-                current_cp[lev][2].get(), PatchType::coarse
-            );
-        }
     }
 
     // SyncCurrent does not include a call to FillBoundary
@@ -112,24 +104,39 @@ WarpX::ComputeMagnetostaticFieldLabFrame()
         }
     }
 
-    // Store the boundary conditions for the field solver if they haven't been
-    // stored yet
-    if (!m_vector_poisson_boundary_handler.bcs_set) {
-        m_vector_poisson_boundary_handler.defineVectorPotentialBCs();
+    // Reference magnetostatic multifabs
+    auto& Afield_fp_nodal = m_magnetostatic_solver->Afield_fp_nodal;
+    auto& current_fp_temp = m_magnetostatic_solver->current_fp_temp;
+
+    // At this point J^{n-1/2} is stored in `current_fp_temp` and J^{n+1/2}
+    // in `current_fp`. In order to calculate B^{n+1} we need to obtain
+    // J^{n+1}. Use extrapolation for this:
+    // J^{n+1} = 1/2 * J_i^{n-1/2} + 3/2 * J_i^{n+1/2}.
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
+        for (int idim = 0; idim < 3; ++idim) {
+            MultiFab::LinComb(
+                *current_fp_temp[lev][idim],
+                0.5_rt, *current_fp_temp[lev][idim], 0,
+                1.5_rt, *current_fp[lev][idim], 0,
+                0, 1, current_fp_temp[lev][idim]->nGrowVect()
+            );
+        }
     }
+
     // set the boundary and current density potentials
-    setVectorPotentialBC(Afield_fp_nodal);
+    setVectorPotentialBC(m_magnetostatic_solver->Afield_fp_nodal);
 
     // Compute the vector potential A, by solving the Poisson equation
     WARPX_ALWAYS_ASSERT_WITH_MESSAGE( !IsPythonCallbackInstalled("poissonsolver"),
         "Python Level Poisson Solve not supported for Magnetostatic implementation.");
 
-    const amrex::Real magnetostatic_absolute_tolerance = self_fields_absolute_tolerance*PhysConst::c;
-
     computeVectorPotential(
-        current_fp, Afield_fp_nodal,
-        self_fields_required_precision, magnetostatic_absolute_tolerance,
-        self_fields_max_iters, self_fields_verbosity
+        current_fp_temp, Afield_fp_nodal,
+        m_magnetostatic_solver->required_precision,
+        m_magnetostatic_solver->absolute_tolerance,
+        m_magnetostatic_solver->max_iters,
+        m_magnetostatic_solver->verbosity
     );
 
     // Grab Interpolation Coefficients
@@ -169,7 +176,20 @@ WarpX::ComputeMagnetostaticFieldLabFrame()
         }
     }
 
+    // At this point Afield_fp contains the vector potential at t = n+1 and
+    // we are ready to obtain B^{n+1}.
     ComputeBfromVectorPotential();
+
+    // copy J^{n+1/2} to current_fp_temp so that it is available at the next
+    // step as J^{n-1/2}
+    for (int lev = 0; lev <= finest_level; ++lev)
+    {
+        // copy 1 component value starting at index 0 to index 0
+        for (int idim = 0; idim < 3; ++idim) {
+            MultiFab::Copy(*current_fp_temp[lev][idim], *current_fp[lev][idim],
+                           0, 0, 1, current_fp_temp[lev][idim]->nGrowVect());
+        }
+    }
 }
 
 /* Compute the vector potential `A` by solving the Poisson equation with `J` as
@@ -228,7 +248,7 @@ WarpX::computeVectorPotential (const amrex::Vector<amrex::Array<std::unique_ptr<
         this->geom,
         this->dmap,
         this->grids,
-        this->m_vector_poisson_boundary_handler,
+        m_magnetostatic_solver->m_boundary_handler,
         WarpX::do_single_precision_comms,
         this->ref_ratio,
         gett_new(0),
@@ -249,9 +269,9 @@ void
 WarpX::setVectorPotentialBC ( amrex::Vector<amrex::Array<std::unique_ptr<amrex::MultiFab>,3>>& A ) const
 {
     // check if any dimension has non-periodic boundary conditions
-    if (!m_vector_poisson_boundary_handler.has_non_periodic) { return; }
+    if (!m_magnetostatic_solver->m_boundary_handler.has_non_periodic) { return; }
 
-    auto dirichlet_flag = m_vector_poisson_boundary_handler.dirichlet_flag;
+    auto dirichlet_flag = m_magnetostatic_solver->m_boundary_handler.dirichlet_flag;
 
     // loop over all mesh refinement levels and set the boundary values
     for (int lev=0; lev <= max_level; lev++) {
@@ -299,95 +319,84 @@ WarpX::setVectorPotentialBC ( amrex::Vector<amrex::Array<std::unique_ptr<amrex::
     } // lev
 }
 
-void MagnetostaticSolver::VectorPoissonBoundaryHandler::defineVectorPotentialBCs ( )
+
+MagnetostaticSolver::MagnetostaticSolver ( int nlevs_max )
 {
-    for (int adim = 0; adim < 3; adim++) {
-#ifdef WARPX_DIM_RZ
-        int dim_start = 0;
-        WarpX& warpx = WarpX::GetInstance();
-        auto geom = warpx.Geom(0);
-        if (geom.ProbLo(0) == 0){
-            lobc[adim][0] = LinOpBCType::Neumann;
-            dirichlet_flag[adim][0] = false;
-            dim_start = 1;
-
-            // handle the r_max boundary explicitly
-            if (WarpX::field_boundary_hi[0] == FieldBoundaryType::PEC) {
-                if (adim == 0) {
-                    hibc[adim][0] = LinOpBCType::Neumann;
-                    dirichlet_flag[adim][1] = false;
-                } else{
-                    hibc[adim][0] = LinOpBCType::Dirichlet;
-                    dirichlet_flag[adim][1] = true;
-                }
-            }
-            else if (WarpX::field_boundary_hi[0] == FieldBoundaryType::Neumann) {
-                hibc[adim][0] = LinOpBCType::Neumann;
-                dirichlet_flag[adim][1] = false;
-            }
-        }
-#else
-        const int dim_start = 0;
-#endif
-        bool ndotA = false;
-        for (int idim=dim_start; idim<AMREX_SPACEDIM; idim++){
-            ndotA = (adim == idim);
+    ReadParameters();
+    AllocateMFs(nlevs_max);
+}
 
-#if defined(WARPX_DIM_XZ) || defined(WARPX_DIM_RZ)
-            if (idim == 1) { ndotA = (adim == 2); }
-#endif
+void MagnetostaticSolver::ReadParameters ()
+{
+    const ParmParse pp_warpx("warpx");
 
-            if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::Periodic
-                && WarpX::field_boundary_hi[idim] == FieldBoundaryType::Periodic ) {
-                lobc[adim][idim] = LinOpBCType::Periodic;
-                hibc[adim][idim] = LinOpBCType::Periodic;
-                dirichlet_flag[adim][idim*2] = false;
-                dirichlet_flag[adim][idim*2+1] = false;
-            }
-            else {
-                has_non_periodic = true;
-                if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::PEC ) {
-                    if (ndotA) {
-                        lobc[adim][idim] = LinOpBCType::Neumann;
-                        dirichlet_flag[adim][idim*2] = false;
-                    } else {
-                        lobc[adim][idim] = LinOpBCType::Dirichlet;
-                        dirichlet_flag[adim][idim*2] = true;
-                    }
-                }
-
-                else if ( WarpX::field_boundary_lo[idim] == FieldBoundaryType::Neumann ) {
-                    lobc[adim][idim] = LinOpBCType::Neumann;
-                    dirichlet_flag[adim][idim*2] = false;
-                }
-                else {
-                    WARPX_ABORT_WITH_MESSAGE(
-                        "Field boundary conditions have to be either periodic, PEC or neumann "
-                        "when using the magnetostatic solver,  but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim])
-                    );
-                }
-
-                if ( WarpX::field_boundary_hi[idim] == FieldBoundaryType::PEC ) {
-                    if (ndotA) {
-                        hibc[adim][idim] = LinOpBCType::Neumann;
-                        dirichlet_flag[adim][idim*2+1] = false;
-                    } else {
-                        hibc[adim][idim] = LinOpBCType::Dirichlet;
-                        dirichlet_flag[adim][idim*2+1] = true;
-                    }
-                }
-                else if ( WarpX::field_boundary_hi[idim] == FieldBoundaryType::Neumann ) {
-                    hibc[adim][idim] = LinOpBCType::Neumann;
-                    dirichlet_flag[adim][idim*2+1] = false;
-                }
-                else {
-                    WARPX_ABORT_WITH_MESSAGE(
-                        "Field boundary conditions have to be either periodic, PEC or neumann "
-                        "when using the magnetostatic solver,  but they are " + GetFieldBCTypeString(WarpX::field_boundary_lo[idim])
-                    );
-                }
-            }
-        }
+    // Grab solver tolerances from input parameters
+    utils::parser::queryWithParser(
+        pp_warpx, "self_fields_required_precision", required_precision);
+    utils::parser::queryWithParser(
+        pp_warpx, "self_fields_absolute_tolerance", absolute_tolerance);
+    absolute_tolerance *= PhysConst::c;
+    utils::parser::queryWithParser(pp_warpx, "self_fields_max_iters", max_iters);
+
+    pp_warpx.query("self_fields_verbosity", verbosity);
+}
+
+void MagnetostaticSolver::AllocateMFs (int nlevs_max)
+{
+    Afield_fp_nodal.resize(nlevs_max);
+    current_fp_temp.resize(nlevs_max);
+}
+
+void MagnetostaticSolver::AllocateLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm,
+                                            const int ncomps, const IntVect& ngJ, const IntVect& ngRho,
+                                            const IntVect& jx_nodal_flag,
+                                            const IntVect& jy_nodal_flag,
+                                            const IntVect& jz_nodal_flag,
+                                            const IntVect& rho_nodal_flag)
+{
+    // The "Afield_fp_nodal" multifab is used to store the vector potential
+    // as obtained from the linear solver.
+    WarpX::AllocInitMultiFab(Afield_fp_nodal[lev][0], amrex::convert(ba, rho_nodal_flag),
+        dm, ncomps, ngRho, lev, "Afield_fp_nodal[x]", 0.0_rt);
+    WarpX::AllocInitMultiFab(Afield_fp_nodal[lev][1], amrex::convert(ba, rho_nodal_flag),
+        dm, ncomps, ngRho, lev, "Afield_fp_nodal[y]", 0.0_rt);
+    WarpX::AllocInitMultiFab(Afield_fp_nodal[lev][2], amrex::convert(ba, rho_nodal_flag),
+        dm, ncomps, ngRho, lev, "Afield_fp_nodal[z]", 0.0_rt);
+
+    // The "current_fp_temp" multifab is used to store the current density
+    // interpolated or extrapolated to appropriate timesteps.
+    WarpX::AllocInitMultiFab(current_fp_temp[lev][0], amrex::convert(ba, jx_nodal_flag),
+        dm, ncomps, ngJ, lev, "current_fp_temp[x]", 0.0_rt);
+    WarpX::AllocInitMultiFab(current_fp_temp[lev][1], amrex::convert(ba, jy_nodal_flag),
+        dm, ncomps, ngJ, lev, "current_fp_temp[y]", 0.0_rt);
+    WarpX::AllocInitMultiFab(current_fp_temp[lev][2], amrex::convert(ba, jz_nodal_flag),
+        dm, ncomps, ngJ, lev, "current_fp_temp[z]", 0.0_rt);
+
+    // the external current density multifab is made nodal to avoid needing to interpolate
+    // to a nodal grid as has to be done for the ion and total current density multifabs
+    // this also allows the external current multifab to not have any ghost cells
+    // WarpX::AllocInitMultiFab(current_fp_external[lev][0], amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))),
+    //     dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[x]", 0.0_rt);
+    // WarpX::AllocInitMultiFab(current_fp_external[lev][1], amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))),
+    //     dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[y]", 0.0_rt);
+    // WarpX::AllocInitMultiFab(current_fp_external[lev][2], amrex::convert(ba, IntVect(AMREX_D_DECL(1,1,1))),
+    //     dm, ncomps, IntVect(AMREX_D_DECL(0,0,0)), lev, "current_fp_external[z]", 0.0_rt);
+
+}
+
+void MagnetostaticSolver::ClearLevel (int lev)
+{
+    for (int i = 0; i < 3; ++i) {
+        Afield_fp_nodal[lev][i].reset();
+        current_fp_temp[lev][i].reset();
     }
-    bcs_set = true;
+}
+
+void MagnetostaticSolver::InitData () {
+    // Store the boundary conditions for the field solver
+    auto& warpx = WarpX::GetInstance();
+    m_boundary_handler.defineVectorPotentialBCs(
+        warpx.field_boundary_lo, warpx.field_boundary_hi, warpx.Geom(0).ProbLo(0)
+    );
+
 }
diff --git a/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver_fwd.H b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver_fwd.H
new file mode 100644
index 00000000000..b0c89088ed7
--- /dev/null
+++ b/Source/FieldSolver/MagnetostaticSolver/MagnetostaticSolver_fwd.H
@@ -0,0 +1,15 @@
+/* Copyright 2024 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Roelof Groenewald (TAE Technologies)
+ *
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_MAGNETOSTATICSOLVER_FWD_H
+#define WARPX_MAGNETOSTATICSOLVER_FWD_H
+
+class MagnetostaticSolver;
+
+#endif /* WARPX_MAGNETOSTATICSOLVER_FWD_H */
diff --git a/Source/Initialization/WarpXInitData.cpp b/Source/Initialization/WarpXInitData.cpp
index 90cf0a7296b..62aee277936 100644
--- a/Source/Initialization/WarpXInitData.cpp
+++ b/Source/Initialization/WarpXInitData.cpp
@@ -19,6 +19,7 @@
 #include "FieldSolver/Fields.H"
 #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
 #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H"
+#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H"
 #include "Filter/BilinearFilter.H"
 #include "Filter/NCIGodfreyFilter.H"
 #include "Initialization/ExternalField.H"
@@ -557,6 +558,10 @@ WarpX::InitData ()
         m_hybrid_pic_model->InitData();
     }
 
+    if (do_magnetostatic_solve) {
+        m_magnetostatic_solver->InitData();
+    }
+
     if (ParallelDescriptor::IOProcessor()) {
         std::cout << "\nGrids Summary:\n";
         printGridSummary(std::cout, 0, finestLevel());
diff --git a/Source/Python/WarpX.cpp b/Source/Python/WarpX.cpp
index bdc1bc0bf1f..e537a08bec7 100644
--- a/Source/Python/WarpX.cpp
+++ b/Source/Python/WarpX.cpp
@@ -15,6 +15,7 @@
 #include <FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H>
 #include <FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H>
 #include <FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H>
+#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H"
 #ifdef WARPX_USE_FFT
 #   include <FieldSolver/SpectralSolver/SpectralKSpace.H>
 #   ifdef WARPX_DIM_RZ
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 5091a6fe6ea..cc46be2aed7 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -16,10 +16,15 @@
 #include "Diagnostics/MultiDiagnostics_fwd.H"
 #include "Diagnostics/ReducedDiags/MultiReducedDiags_fwd.H"
 #include "EmbeddedBoundary/WarpXFaceInfoBox_fwd.H"
+#include "FieldSolver/ElectrostaticSolver.H"
 #include "FieldSolver/ElectrostaticSolver/SemiImplicitSolver_fwd.H"
+#include "FieldSolver/Fields.H"
 #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver_fwd.H"
 #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties_fwd.H"
 #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel_fwd.H"
+#include "FieldSolver/ImplicitSolvers/ImplicitSolver.H"
+#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
+#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver_fwd.H"
 #include "Filter/NCIGodfreyFilter_fwd.H"
 #include "Initialization/ExternalField_fwd.H"
 #include "Particles/ParticleBoundaryBuffer_fwd.H"
@@ -39,11 +44,6 @@
 #include "AcceleratorLattice/AcceleratorLattice.H"
 #include "Evolve/WarpXDtType.H"
 #include "Evolve/WarpXPushType.H"
-#include "FieldSolver/Fields.H"
-#include "FieldSolver/ElectrostaticSolver.H"
-#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H"
-#include "FieldSolver/ImplicitSolvers/ImplicitSolver.H"
-#include "FieldSolver/ImplicitSolvers/WarpXSolverVec.H"
 #include "Filter/BilinearFilter.H"
 #include "Parallelization/GuardCellManager.H"
 #include "Utils/Parser/IntervalsParser.H"
@@ -1012,7 +1012,6 @@ public:
                                       amrex::Vector<std::unique_ptr<amrex::MultiFab> >& phi) const;
 
     // Magnetostatic Solver Interface
-    MagnetostaticSolver::VectorPoissonBoundaryHandler m_vector_poisson_boundary_handler;
     void ComputeMagnetostaticField ();
     void ComputeMagnetostaticFieldLabFrame ();
     void computeVectorPotential (const amrex::Vector<amrex::Array<std::unique_ptr<amrex::MultiFab>, 3> >& curr,
@@ -1521,7 +1520,6 @@ private:
     // Memory buffers for computing magnetostatic fields
     // Vector Potential A and previous step.  Time buffer needed for computing dA/dt to first order
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Afield_fp;
-    amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Afield_fp_nodal;
 
     // Same as Bfield_fp/Efield_fp for reading external field data
     amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > Efield_fp_external;
@@ -1647,12 +1645,15 @@ private:
     // Macroscopic properties
     std::unique_ptr<MacroscopicProperties> m_macroscopic_properties;
 
-    // Hybrid PIC algorithm parameters
+    // Hybrid PIC algorithm class
     std::unique_ptr<HybridPICModel> m_hybrid_pic_model;
 
     // Semi-implicit electrostatic algorithm
     std::unique_ptr<SemiImplicitSolver> m_semi_implicit_solver;
 
+    // Magnetostatic algorithm class
+    std::unique_ptr<MagnetostaticSolver> m_magnetostatic_solver;
+
     // Load balancing
     /** Load balancing intervals that reads the "load_balance_intervals" string int the input file
      * for getting steps at which load balancing is performed */
diff --git a/Source/WarpX.cpp b/Source/WarpX.cpp
index fba01440538..095a680261b 100644
--- a/Source/WarpX.cpp
+++ b/Source/WarpX.cpp
@@ -20,6 +20,7 @@
 #include "FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H"
 #include "FieldSolver/FiniteDifferenceSolver/MacroscopicProperties/MacroscopicProperties.H"
 #include "FieldSolver/FiniteDifferenceSolver/HybridPICModel/HybridPICModel.H"
+#include "FieldSolver/MagnetostaticSolver/MagnetostaticSolver.H"
 #ifdef WARPX_USE_FFT
 #   include "FieldSolver/SpectralSolver/SpectralKSpace.H"
 #   ifdef WARPX_DIM_RZ
@@ -329,18 +330,18 @@ WarpX::WarpX ()
     Efield_fp.resize(nlevs_max);
     Bfield_fp.resize(nlevs_max);
 
-    // Only allocate vector potential arrays when using the Magnetostatic Solver
-    if (do_magnetostatic_solve)
-    {
-        Afield_fp_nodal.resize(nlevs_max);
-        Afield_fp.resize(nlevs_max);
-    }
-
     // Initialize the semi-implicit electrostatic solver if required
     if (electrostatic_solver_id == ElectrostaticSolverAlgo::LabFrameSemiImplicit)
     {
         m_semi_implicit_solver = std::make_unique<SemiImplicitSolver>(nlevs_max);
     }
+    // Initialize the magnetostatic solver if required
+    // Only allocate vector potential arrays when using the Magnetostatic Solver
+    if (do_magnetostatic_solve)
+    {
+        m_magnetostatic_solver = std::make_unique<MagnetostaticSolver>(nlevs_max);
+        Afield_fp.resize(nlevs_max);
+    }
 
     if (fft_do_time_averaging)
     {
@@ -2118,12 +2119,6 @@ WarpX::ClearLevel (int lev)
             current_fp_vay[lev][i].reset();
         }
 
-        if (do_magnetostatic_solve)
-        {
-            Afield_fp_nodal[lev][i].reset();
-            Afield_fp[lev][i].reset();
-        }
-
         current_cp[lev][i].reset();
         Efield_cp [lev][i].reset();
         Bfield_cp [lev][i].reset();
@@ -2137,6 +2132,13 @@ WarpX::ClearLevel (int lev)
     {
         m_semi_implicit_solver->ClearLevel(lev);
     }
+    if (do_magnetostatic_solve)
+    {
+        m_magnetostatic_solver->ClearLevel(lev);
+        for (int i = 0; i < 3; ++i) {
+            Afield_fp[lev][i].reset();
+        }
+    }
 
     if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC)
     {
@@ -2376,9 +2378,6 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
 
     if (do_magnetostatic_solve)
     {
-        AllocInitMultiFab(Afield_fp_nodal[lev][0], amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, lev, "Afield_fp_nodal[x]", 0.0_rt);
-        AllocInitMultiFab(Afield_fp_nodal[lev][1], amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, lev, "Afield_fp_nodal[y]", 0.0_rt);
-        AllocInitMultiFab(Afield_fp_nodal[lev][2], amrex::convert(ba, rho_nodal_flag), dm, ncomps, ngRho, lev, "Afield_fp_nodal[z]", 0.0_rt);
         AllocInitMultiFab(Afield_fp[lev][0], amrex::convert(ba, Ex_nodal_flag), dm, ncomps, ngEB, lev, "Afield_fp[x]", 0.0_rt);
         AllocInitMultiFab(Afield_fp[lev][1], amrex::convert(ba, Ey_nodal_flag), dm, ncomps, ngEB, lev, "Afield_fp[y]", 0.0_rt);
         AllocInitMultiFab(Afield_fp[lev][2], amrex::convert(ba, Ez_nodal_flag), dm, ncomps, ngEB, lev, "Afield_fp[z]", 0.0_rt);
@@ -2392,6 +2391,15 @@ WarpX::AllocLevelMFs (int lev, const BoxArray& ba, const DistributionMapping& dm
         );
     }
 
+    // Allocate extra multifabs needed by the magnetostatic algorithm.
+    if (do_magnetostatic_solve)
+    {
+        m_magnetostatic_solver->AllocateLevelMFs(
+            lev, ba, dm, ncomps, ngJ, ngRho, jx_nodal_flag, jy_nodal_flag,
+            jz_nodal_flag, rho_nodal_flag
+        );
+    }
+
     // Allocate extra multifabs needed by the kinetic-fluid hybrid algorithm.
     if (WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::HybridPIC)
     {
@@ -3535,8 +3543,8 @@ WarpX::getFieldPointerUnchecked (const FieldType field_type, const int lev, cons
         case FieldType::phi_fp :
             field_pointer = phi_fp[lev].get();
             break;
-       case FieldType::vector_potential_fp :
-            field_pointer = Afield_fp_nodal[lev][direction].get();
+       case FieldType::Afield_fp :
+            field_pointer = Afield_fp[lev][direction].get();
             break;
         case FieldType::Efield_cp :
             field_pointer = Efield_cp[lev][direction].get();