Add comments

7 files changed, 151 insertions, 65 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index 7a3f5aa9b..a84cb0a53 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -11,22 +11,34 @@
 
 using namespace amrex;
 
-void FiniteDifferenceSolver::EvolveB ( VectorField& Bfield,
-                                       VectorField const& Efield,
-                                       amrex::Real const dt ) {
+/**
+ * \brief Update the B field, over one timestep
+ */
+void FiniteDifferenceSolver::EvolveB (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    amrex::Real const dt ) {
 
    // Select algorithm (The choice of algorithm is a runtime option,
    // but we compile code for each algorithm, using templates)
 #ifdef WARPX_DIM_RZ
     if (m_fdtd_algo == MaxwellSolverAlgo::Yee){
+
         EvolveBCylindrical <CylindricalYeeAlgorithm> ( Bfield, Efield, dt );
+
 #else
     if (m_do_nodal) {
+
         EvolveBCartesian <NodalAlgorithm> ( Bfield, Efield, dt );
+
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
+
         EvolveBCartesian <YeeAlgorithm> ( Bfield, Efield, dt );
+
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
+
         EvolveBCartesian <CKCAlgorithm> ( Bfield, Efield, dt );
+
 #endif
     } else {
         amrex::Abort("Unknown algorithm");
@@ -34,12 +46,14 @@ void FiniteDifferenceSolver::EvolveB ( VectorField& Bfield,
 
 }
 
+
 #ifndef WARPX_DIM_RZ
 
 template<typename T_Algo>
-void FiniteDifferenceSolver::EvolveBCartesian ( VectorField& Bfield,
-                                               VectorField const& Efield,
-                                               amrex::Real const dt ) {
+void FiniteDifferenceSolver::EvolveBCartesian (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    amrex::Real const dt ) {
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef _OPENMP
@@ -95,9 +109,10 @@ void FiniteDifferenceSolver::EvolveBCartesian ( VectorField& Bfield,
 #else // corresponds to ifndef WARPX_DIM_RZ
 
 template<typename T_Algo>
-void FiniteDifferenceSolver::EvolveBCylindrical ( VectorField& Bfield,
-                                               VectorField const& Efield,
-                                               amrex::Real const dt ) {
+void FiniteDifferenceSolver::EvolveBCylindrical (
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+    std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+    amrex::Real const dt ) {
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef _OPENMP
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H
index 5b0d5e718..771c44a22 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CKCAlgorithm.H
@@ -23,7 +23,7 @@ struct CKCAlgorithm {
         Real const inv_dx = 1./cell_size[0];
         Real const inv_dy = 1./cell_size[1];
         Real const inv_dz = 1./cell_size[2];
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             const Real delta = std::max( { inv_dx,inv_dy,inv_dz } );
             const Real rx = (inv_dx/delta)*(inv_dx/delta);
             const Real ry = (inv_dy/delta)*(inv_dy/delta);
@@ -53,7 +53,7 @@ struct CKCAlgorithm {
             gammax *= inv_dx;
             gammay *= inv_dy;
             gammaz *= inv_dz;
-        #elif defined WARPX_DIM_XZ
+#        elif defined WARPX_DIM_XZ
             const Real delta = std::max(inv_dx,inv_dz);
             const Real rx = (inv_dx/delta)*(inv_dx/delta);
             const Real rz = (inv_dz/delta)*(inv_dz/delta);
@@ -65,7 +65,7 @@ struct CKCAlgorithm {
             betazx *= inv_dz;
             alphax *= inv_dx;
             alphaz *= inv_dz;
-        #endif
+#        endif
 
         // Store the coefficients in array `stencil_coefs`, in prescribed order
         stencil_coefs_x.resize(6);
@@ -88,6 +88,8 @@ struct CKCAlgorithm {
         stencil_coefs_z[4] = gammaz;
     }
 
+    /**
+    /* Perform derivative along x on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDx(
         amrex::Array4<amrex::Real> const& F,
@@ -98,7 +100,7 @@ struct CKCAlgorithm {
         amrex::Real betaxy = coefs_x[2];
         amrex::Real betaxz = coefs_x[3];
         amrex::Real gammax = coefs_x[4];
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             return alphax * (F(i+1,j  ,k  ) - F(i,  j,  k  ))
                  + betaxy * (F(i+1,j+1,k  ) - F(i  ,j+1,k  )
                           +  F(i+1,j-1,k  ) - F(i  ,j-1,k  ))
@@ -108,15 +110,17 @@ struct CKCAlgorithm {
                           +  F(i+1,j-1,k+1) - F(i  ,j-1,k+1)
                           +  F(i+1,j+1,k-1) - F(i  ,j+1,k-1)
                           +  F(i+1,j-1,k-1) - F(i  ,j-1,k-1));
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return alphax * (F(i+1,j  ,k  ) - F(i,  j,  k  ))
                  + betaxz * (F(i+1,j+1,k  ) - F(i  ,j+1,k  )
                           +  F(i+1,j-1,k  ) - F(i  ,j-1,k  ));
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along x on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDx(
+    static amrex::Real Downwardx(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_x, int const n_coefs_x,
         int const i, int const j, int const k ) {
@@ -125,13 +129,15 @@ struct CKCAlgorithm {
         return inv_dx*( F(i,j,k) - F(i-1,j,k) );
     };
 
+    /**
+    /* Perform derivative along y on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDy(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_y, int const n_coefs_y,
         int const i, int const j, int const k ) {
 
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             amrex::Real alphay = coefs_y[1];
             amrex::Real betayz = coefs_y[2];
             amrex::Real betayx = coefs_y[3];
@@ -145,25 +151,29 @@ struct CKCAlgorithm {
                           +  F(i-1,j+1,k+1) - F(i-1,j  ,k+1)
                           +  F(i+1,j+1,k-1) - F(i+1,j  ,k-1)
                           +  F(i-1,j+1,k-1) - F(i-1,j  ,k-1));
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return 0; // 2D Cartesian: derivative along y is 0
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along y on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDy(
+    static amrex::Real Downwardy(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_y, int const n_coefs_y,
         int const i, int const j, int const k ) {
 
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             amrex::Real inv_dy = coefs_y[0];
             return inv_dy*( F(i,j,k) - F(i,j-1,k) );
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return 0; // 2D Cartesian: derivative along y is 0
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along z on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDz(
         amrex::Array4<amrex::Real> const& F,
@@ -174,7 +184,7 @@ struct CKCAlgorithm {
         amrex::Real betazx = coefs_z[2];
         amrex::Real betazy = coefs_z[3];
         amrex::Real gammaz = coefs_z[4];
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             return alphaz * (F(i  ,j  ,k+1) - F(i  ,j  ,k  ))
                  + betazx * (F(i+1,j  ,k+1) - F(i+1,j  ,k  )
                           +  F(i-1,j  ,k+1) - F(i-1,j  ,k  ))
@@ -184,25 +194,27 @@ struct CKCAlgorithm {
                           +  F(i-1,j+1,k+1) - F(i-1,j+1,k  )
                           +  F(i+1,j-1,k+1) - F(i+1,j-1,k  )
                           +  F(i-1,j-1,k+1) - F(i-1,j-1,k  ));
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return alphaz * (F(i  ,j+1,k  ) - F(i  ,j  ,k  ))
                  + betazx * (F(i+1,j+1,k  ) - F(i+1,j  ,k  )
                           +  F(i-1,j+1,k  ) - F(i-1,j  ,k  ));
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along z on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDz(
+    static amrex::Real Downwardz(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_z, int const n_coefs_z,
         int const i, int const j, int const k ) {
 
         amrex::Real inv_dz = coefs_z[0];
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             return inv_dz*( F(i,j,k) - F(i,j,k-1) );
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return inv_dz*( F(i,j,k) - F(i,j-1,k) );
-        #endif
+#        endif
     };
 
 };
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
index 2cb88fb79..70dfe23b5 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H
@@ -22,7 +22,7 @@ struct CylindricalYeeAlgorithm {
     /** Applies the differential operator `1/r * d(rF)/dr`,
      * where `F` is on a *nodal* grid in `r`
      * and the differential operator is evaluated on *cell-centered* grid.
-     * The input parameter `r` is given at the cell-centered position*/
+     * The input parameter `r` is given at the cell-centered position */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDrr_over_r(
         amrex::Array4<amrex::Real> const& F,
@@ -34,6 +34,8 @@ struct CylindricalYeeAlgorithm {
         return 1./r * inv_dr*( (r+0.5*dr)*F(i+1,j,k,comp) - (r-0.5*dr)*F(i,j,k,comp) );
     };
 
+    /**
+    /* Perform derivative along r on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDr(
         amrex::Array4<amrex::Real> const& F,
@@ -44,6 +46,8 @@ struct CylindricalYeeAlgorithm {
         return inv_dr*( F(i+1,j,k,comp) - F(i,j,k,comp) );
     };
 
+    /**
+    /* Perform derivative along r on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real DownwardDr(
         amrex::Array4<amrex::Real> const& F,
@@ -54,6 +58,8 @@ struct CylindricalYeeAlgorithm {
         return inv_dr*( F(i,j,k,comp) - F(i-1,j,k,comp) );
     };
 
+    /**
+    /* Perform derivative along z on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDz(
         amrex::Array4<amrex::Real> const& F,
@@ -64,6 +70,8 @@ struct CylindricalYeeAlgorithm {
         return inv_dz*( F(i,j+1,k,comp) - F(i,j,k,comp) );
     };
 
+    /**
+    /* Perform derivative along z on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real DownwardDz(
         amrex::Array4<amrex::Real> const& F,
@@ -74,6 +82,7 @@ struct CylindricalYeeAlgorithm {
         return inv_dz*( F(i,j,k,comp) - F(i,j-1,k,comp) );
     };
 
+    /** Divide by the radius `r` and avoid potential singularities */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real DivideByR(
         amrex::Array4<amrex::Real> const& F,
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/NodalAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/NodalAlgorithm.H
index 8ffd23b79..db07ca2b3 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/NodalAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/NodalAlgorithm.H
@@ -22,6 +22,10 @@ struct NodalAlgorithm {
         stencil_coefs_z[0] = 1./cell_size[2];
     }
 
+    /**
+    /* Perform derivative along x
+    /* (For a solver on a staggered grid, `UpwardDx` and `DownwardDx` take into
+    /* account the staggering; but for `NodalAlgorithm`, they are equivalent) */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDx(
         amrex::Array4<amrex::Real> const& F,
@@ -32,8 +36,12 @@ struct NodalAlgorithm {
         return 0.5*inv_dx*( F(i+1,j,k) - F(i-1,j,k) );
     };
 
+    /**
+    /* Perform derivative along x
+    /* (For a solver on a staggered grid, `UpwardDx` and `DownwardDx` take into
+    /* account the staggering; but for `NodalAlgorithm`, they are equivalent) */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDx(
+    static amrex::Real Downwardx(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_x, int const n_coefs_x,
         int const i, int const j, int const k ) {
@@ -42,34 +50,46 @@ struct NodalAlgorithm {
         return 0.5*inv_dx*( F(i+1,j,k) - F(i-1,j,k) );
     };
 
+    /**
+    /* Perform derivative along y
+    /* (For a solver on a staggered grid, `UpwardDy` and `DownwardDy` take into
+    /* account the staggering; but for `NodalAlgorithm`, they are equivalent) */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDy(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_y, int const n_coefs_y,
         int const i, int const j, int const k ) {
 
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             amrex::Real inv_dy = coefs_y[0];
             return 0.5*inv_dy*( F(i,j+1,k) - F(i,j-1,k) );
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return 0; // 2D Cartesian: derivative along y is 0
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along y
+    /* (For a solver on a staggered grid, `UpwardDy` and `DownwardDy` take into
+    /* account the staggering; but for `NodalAlgorithm`, they are equivalent) */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDy(
+    static amrex::Real Downwardy(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_y, int const n_coefs_y,
         int const i, int const j, int const k ) {
 
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             amrex::Real inv_dy = coefs_y[0];
             return 0.5*inv_dy*( F(i,j+1,k) - F(i,j-1,k) );
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return 0; // 2D Cartesian: derivative along y is 0
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along z
+    /* (For a solver on a staggered grid, `UpwardDz` and `DownwardDz` take into
+    /* account the staggering; but for `NodalAlgorithm`, they are equivalent) */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDz(
         amrex::Array4<amrex::Real> const& F,
@@ -77,25 +97,29 @@ struct NodalAlgorithm {
         int const i, int const j, int const k ) {
 
         amrex::Real inv_dz = coefs_z[0];
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             return 0.5*inv_dz*( F(i,j,k+1) - F(i,j,k-1) );
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return 0.5*inv_dz*( F(i,j+1,k) - F(i,j-1,k) );
-        #endif
+#        endif
     };
 
+    /**
+    /* Perform derivative along z
+    /* (For a solver on a staggered grid, `UpwardDz` and `DownwardDz` take into
+    /* account the staggering; but for `NodalAlgorithm`, they are equivalent) */
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDz(
+    static amrex::Real Downwardz(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_z, int const n_coefs_z,
         int const i, int const j, int const k ) {
 
         amrex::Real inv_dz = coefs_z[0];
-        #if defined WARPX_DIM_3D
+#        if defined WARPX_DIM_3D
             return 0.5*inv_dz*( F(i,j,k+1) - F(i,j,k-1) );
-        #elif (defined WARPX_DIM_XZ)
+#        elif (defined WARPX_DIM_XZ)
             return 0.5*inv_dz*( F(i,j+1,k) - F(i,j-1,k) );
-        #endif
+#        endif
     };
 
 };
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H
index 54057091d..d967662f1 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceAlgorithms/YeeAlgorithm.H
@@ -22,6 +22,8 @@ struct YeeAlgorithm {
         stencil_coefs_z[0] = 1./cell_size[2];
     }
 
+    /**
+    /* Perform derivative along x on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDx(
         amrex::Array4<amrex::Real> const& F,
@@ -32,8 +34,10 @@ struct YeeAlgorithm {
         return inv_dx*( F(i+1,j,k) - F(i,j,k) );
     };
 
+    /**
+    /* Perform derivative along x on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDx(
+    static amrex::Real Downwardx(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_x, int const n_coefs_x,
         int const i, int const j, int const k ) {
@@ -42,6 +46,8 @@ struct YeeAlgorithm {
         return inv_dx*( F(i,j,k) - F(i-1,j,k) );
     };
 
+    /**
+    /* Perform derivative along y on a cell-centered grid, from a nodal field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDy(
         amrex::Array4<amrex::Real> const& F,
@@ -56,8 +62,10 @@ struct YeeAlgorithm {
         #endif
     };
 
+    /**
+    /* Perform derivative along y on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDy(
+    static amrex::Real Downwardy(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_y, int const n_coefs_y,
         int const i, int const j, int const k ) {
@@ -70,7 +78,9 @@ struct YeeAlgorithm {
         #endif
     };
 
-    AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
+    /**
+    /* Perform derivative along z on a cell-centered grid, from a nodal field `F`*/
+   AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
     static amrex::Real UpwardDz(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_z, int const n_coefs_z,
@@ -84,8 +94,10 @@ struct YeeAlgorithm {
         #endif
     };
 
+    /**
+    /* Perform derivative along z on a nodal grid, from a cell-centered field `F`*/
     AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE
-    static amrex::Real DownwardDz(
+    static amrex::Real Downwardz(
         amrex::Array4<amrex::Real> const& F,
         amrex::Real const* coefs_z, int const n_coefs_z,
         int const i, int const j, int const k ) {
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index 770aa4810..3041bebbe 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -13,16 +13,14 @@ class FiniteDifferenceSolver
 {
     public:
 
-        using VectorField = std::array< std::unique_ptr<amrex::MultiFab>, 3 >;
-
         // Constructor
         FiniteDifferenceSolver (
             int const fdtd_algo,
             std::array<amrex::Real,3> cell_size,
             bool const do_nodal );
 
-        void EvolveB ( VectorField& Bfield,
-                       VectorField const& Efield,
+        void EvolveB ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+                       std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
                        amrex::Real const dt );
     private:
 
@@ -46,14 +44,16 @@ class FiniteDifferenceSolver
 
 #ifdef WARPX_DIM_RZ
         template< typename T_Algo >
-        void EvolveBCylindrical ( VectorField& Bfield,
-                               VectorField const& Efield,
-                               amrex::Real const dt );
+        void EvolveBCylindrical (
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+            amrex::Real const dt );
 #else
         template< typename T_Algo >
-        void EvolveBCartesian ( VectorField& Bfield,
-                               VectorField const& Efield,
-                               amrex::Real const dt );
+        void EvolveBCartesian (
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
+            std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+            amrex::Real const dt );
 #endif
 
 };
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
index c333c6d17..06f5960a6 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.cpp
@@ -1,15 +1,22 @@
 #include "WarpXAlgorithmSelection.H"
 #ifdef WARPX_DIM_RZ
-    #include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
+#    include "FiniteDifferenceAlgorithms/CylindricalYeeAlgorithm.H"
 #else
-    #include "FiniteDifferenceAlgorithms/YeeAlgorithm.H"
-    #include "FiniteDifferenceAlgorithms/CKCAlgorithm.H"
-    #include "FiniteDifferenceAlgorithms/NodalAlgorithm.H"
+#    include "FiniteDifferenceAlgorithms/YeeAlgorithm.H"
+#    include "FiniteDifferenceAlgorithms/CKCAlgorithm.H"
+#    include "FiniteDifferenceAlgorithms/NodalAlgorithm.H"
 #endif
 #include "FiniteDifferenceSolver.H"
 #include "WarpX.H"
 
-// Constructor
+/* \brief Initialize the finite-difference Maxwell solver (for a given refinement level)
+ *
+ * This function initializes the stencil coefficients for the chosen finite-difference algorithm
+ *
+ * \param fdtd_algo Identifies the chosen algorithm, as defined in WarpXAlgorithmSelection.H
+ * \param cell_size Cell size along each dimension, for the chosen refinement level
+ * \param do_nodal  Whether the solver is applied to a nodal or staggered grid
+ */
 FiniteDifferenceSolver::FiniteDifferenceSolver (
     int const fdtd_algo,
     std::array<amrex::Real,3> cell_size,
@@ -25,18 +32,25 @@ FiniteDifferenceSolver::FiniteDifferenceSolver (
     m_nmodes = WarpX::GetInstance().n_rz_azimuthal_modes;
     m_rmin = WarpX::GetInstance().Geom(0).ProbLo(0);
     if (fdtd_algo == MaxwellSolverAlgo::Yee) {
+
         CylindricalYeeAlgorithm::InitializeStencilCoefficients( cell_size,
             stencil_coefs_r, stencil_coefs_z );
 #else
     if (do_nodal) {
+
         NodalAlgorithm::InitializeStencilCoefficients( cell_size,
             stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+
     } else if (fdtd_algo == MaxwellSolverAlgo::Yee) {
+
         YeeAlgorithm::InitializeStencilCoefficients( cell_size,
             stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+
     } else if (fdtd_algo == MaxwellSolverAlgo::CKC) {
+
         CKCAlgorithm::InitializeStencilCoefficients( cell_size,
             stencil_coefs_x, stencil_coefs_y, stencil_coefs_z );
+
 #endif
     } else {
         amrex::Abort("Unknown algorithm");