Add timers in routines that depend on cell-related work (#1692)

* Add timers

* eol

* AtomicAdd

* lev argument for getCosts

* style

* style

* wip

* eol

* .ipynb

* passing down lev

* eol

* passing lev

* eol

* Update Source/Particles/Collision/PairWiseCoulombCollision.cpp

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

* Add for QED and ionization routines

* eol

* remove unneeded

* mfi-->pti

* move cost

* eol

Co-authored-by: Axel Huebl <axel.huebl@plasma.ninja>

26 files changed, 348 insertions, 186 deletions

diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
index 338cc31f9..2f66796f3 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveB.cpp
@@ -5,6 +5,7 @@
  * License: BSD-3-Clause-LBNL
  */
 
+#include "WarpX.H"
 #include "Utils/WarpXAlgorithmSelection.H"
 #include "FiniteDifferenceSolver.H"
 #ifdef WARPX_DIM_RZ
@@ -24,27 +25,27 @@ using namespace amrex;
 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 ) {
+    int lev, 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 );
+        EvolveBCylindrical <CylindricalYeeAlgorithm> ( Bfield, Efield, lev, dt );
 
 #else
     if (m_do_nodal) {
 
-        EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, dt );
+        EvolveBCartesian <CartesianNodalAlgorithm> ( Bfield, Efield, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, dt );
+        EvolveBCartesian <CartesianYeeAlgorithm> ( Bfield, Efield, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, dt );
+        EvolveBCartesian <CartesianCKCAlgorithm> ( Bfield, Efield, lev, dt );
 
 #endif
     } else {
@@ -60,13 +61,20 @@ template<typename T_Algo>
 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 ) {
+    int lev, amrex::Real const dt ) {
+
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
     for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
 
         // Extract field data for this grid/tile
         Array4<Real> const& Bx = Bfield[0]->array(mfi);
@@ -109,8 +117,13 @@ void FiniteDifferenceSolver::EvolveBCartesian (
 
         );
 
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     }
-
 }
 
 #else // corresponds to ifndef WARPX_DIM_RZ
@@ -119,13 +132,20 @@ template<typename T_Algo>
 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 ) {
+    int lev, amrex::Real const dt ) {
+
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
     for ( MFIter mfi(*Bfield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
 
         // Extract field data for this grid/tile
         Array4<Real> const& Br = Bfield[0]->array(mfi);
@@ -214,8 +234,13 @@ void FiniteDifferenceSolver::EvolveBCylindrical (
 
         );
 
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     }
-
 }
 
 #endif // corresponds to ifndef WARPX_DIM_RZ
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
index 3762d3f71..cd1e3bdf9 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
+++ b/Source/FieldSolver/FiniteDifferenceSolver/EvolveE.cpp
@@ -5,6 +5,7 @@
  * License: BSD-3-Clause-LBNL
  */
 
+#include "WarpX.H"
 #include "Utils/WarpXAlgorithmSelection.H"
 #include "FiniteDifferenceSolver.H"
 #ifdef WARPX_DIM_RZ
@@ -28,27 +29,27 @@ void FiniteDifferenceSolver::EvolveE (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
-    amrex::Real const dt ) {
+    int lev, 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){
 
-        EvolveECylindrical <CylindricalYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, dt );
+        EvolveECylindrical <CylindricalYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
 
 #else
     if (m_do_nodal) {
 
-        EvolveECartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, Ffield, dt );
+        EvolveECartesian <CartesianNodalAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::Yee) {
 
-        EvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, dt );
+        EvolveECartesian <CartesianYeeAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
 
     } else if (m_fdtd_algo == MaxwellSolverAlgo::CKC) {
 
-        EvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, Ffield, dt );
+        EvolveECartesian <CartesianCKCAlgorithm> ( Efield, Bfield, Jfield, Ffield, lev, dt );
 
 #endif
     } else {
@@ -66,8 +67,9 @@ void FiniteDifferenceSolver::EvolveECartesian (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
-    amrex::Real const dt ) {
+    int lev, amrex::Real const dt ) {
 
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
     Real constexpr c2 = PhysConst::c * PhysConst::c;
 
     // Loop through the grids, and over the tiles within each grid
@@ -75,6 +77,11 @@ void FiniteDifferenceSolver::EvolveECartesian (
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
     for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
 
         // Extract field data for this grid/tile
         Array4<Real> const& Ex = Efield[0]->array(mfi);
@@ -150,6 +157,12 @@ void FiniteDifferenceSolver::EvolveECartesian (
 
         }
 
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     }
 
 }
@@ -162,13 +175,20 @@ void FiniteDifferenceSolver::EvolveECylindrical (
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
     std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
     std::unique_ptr<amrex::MultiFab> const& Ffield,
-    amrex::Real const dt ) {
+    int lev, amrex::Real const dt ) {
+
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
 
     // Loop through the grids, and over the tiles within each grid
 #ifdef AMREX_USE_OMP
 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
 #endif
     for ( MFIter mfi(*Efield[0], TilingIfNotGPU()); mfi.isValid(); ++mfi ) {
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
 
         // Extract field data for this grid/tile
         Array4<Real> const& Er = Efield[0]->array(mfi);
@@ -339,6 +359,12 @@ void FiniteDifferenceSolver::EvolveECylindrical (
 
         } // end of if condition for F
 
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     } // end of loop over grid/tiles
 
 }
diff --git a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
index bb96dee45..72fabb54e 100644
--- a/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
+++ b/Source/FieldSolver/FiniteDifferenceSolver/FiniteDifferenceSolver.H
@@ -38,13 +38,13 @@ class FiniteDifferenceSolver
 
         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 );
+                       int lev, amrex::Real const dt );
 
         void EvolveE ( std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Efield,
                        std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
                        std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
                        std::unique_ptr<amrex::MultiFab> const& Ffield,
-                       amrex::Real const dt );
+                       int lev, amrex::Real const dt );
 
         void EvolveF ( std::unique_ptr<amrex::MultiFab>& Ffield,
                        std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
@@ -114,6 +114,7 @@ class FiniteDifferenceSolver
         void EvolveBCylindrical (
             std::array< std::unique_ptr<amrex::MultiFab>, 3 >& Bfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Efield,
+            const int lev,
             amrex::Real const dt );
 
         template< typename T_Algo >
@@ -122,6 +123,7 @@ class FiniteDifferenceSolver
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
             std::unique_ptr<amrex::MultiFab> const& Ffield,
+            const int lev,
             amrex::Real const dt );
 
         template< typename T_Algo >
@@ -142,7 +144,7 @@ class FiniteDifferenceSolver
         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 );
+            int lev, amrex::Real const dt );
 
         template< typename T_Algo >
         void EvolveECartesian (
@@ -150,7 +152,7 @@ class FiniteDifferenceSolver
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Bfield,
             std::array< std::unique_ptr<amrex::MultiFab>, 3 > const& Jfield,
             std::unique_ptr<amrex::MultiFab> const& Ffield,
-            amrex::Real const dt );
+            int lev, amrex::Real const dt );
 
         template< typename T_Algo >
         void EvolveFCartesian (
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
index 2623a3ea6..22b04aeba 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.H
@@ -52,7 +52,8 @@ class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm
          *                           the three components of the current density
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection (SpectralFieldData& field_data,
+        virtual void CurrentCorrection (const int lev,
+                                        SpectralFieldData& field_data,
                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                         const std::unique_ptr<amrex::MultiFab>& rho) override final;
 
@@ -66,7 +67,8 @@ class ComovingPsatdAlgorithm : public SpectralBaseAlgorithm
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldData& field_data,
+        virtual void VayDeposition (const int lev,
+                                    SpectralFieldData& field_data,
                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
     private:
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
index ca82d5a8c..96b993dcf 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/ComovingPsatdAlgorithm.cpp
@@ -391,7 +391,8 @@ void ComovingPsatdAlgorithm::InitializeSpectralCoefficients (const SpectralKSpac
 }
 
 void
-ComovingPsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
+ComovingPsatdAlgorithm::CurrentCorrection (const int lev,
+                                           SpectralFieldData& field_data,
                                            std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                            const std::unique_ptr<amrex::MultiFab>& rho)
 {
@@ -401,11 +402,11 @@ ComovingPsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
     using Idx = SpectralFieldIndex;
 
     // Forward Fourier transform of J and rho
-    field_data.ForwardTransform(*current[0], Idx::Jx, 0);
-    field_data.ForwardTransform(*current[1], Idx::Jy, 0);
-    field_data.ForwardTransform(*current[2], Idx::Jz, 0);
-    field_data.ForwardTransform(*rho, Idx::rho_old, 0);
-    field_data.ForwardTransform(*rho, Idx::rho_new, 1);
+    field_data.ForwardTransform(lev, *current[0], Idx::Jx, 0);
+    field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0);
+    field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0);
+    field_data.ForwardTransform(lev, *rho, Idx::rho_old, 0);
+    field_data.ForwardTransform(lev, *rho, Idx::rho_new, 1);
 
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
@@ -487,13 +488,14 @@ ComovingPsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data,
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform(*current[0], Idx::Jx, 0);
-    field_data.BackwardTransform(*current[1], Idx::Jy, 0);
-    field_data.BackwardTransform(*current[2], Idx::Jz, 0);
+    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0);
+    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0);
+    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0);
 }
 
 void
-ComovingPsatdAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/,
+ComovingPsatdAlgorithm::VayDeposition (const int /*lev*/,
+                                       SpectralFieldData& /*field_data*/,
                                        std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
 {
     amrex::Abort("Vay deposition not implemented for comoving PSATD");
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H
index 18a7f2dcd..c5680a4ce 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.H
@@ -42,7 +42,8 @@ class GalileanPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
          *                           storing the three components of the current density
          * \param[in]     rho        unique pointer to MultiFab storing the charge density
          */
-        virtual void CurrentCorrection ( SpectralFieldDataRZ& field_data,
+        virtual void CurrentCorrection ( const int lev,
+                                         SpectralFieldDataRZ& field_data,
                                          std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                          const std::unique_ptr<amrex::MultiFab>& rho ) override final;
 
@@ -56,7 +57,7 @@ class GalileanPsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldDataRZ& field_data,
+        virtual void VayDeposition (const int lev, SpectralFieldDataRZ& field_data,
                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
     private:
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp
index 85de8ffc1..dd608fa18 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/GalileanPsatdAlgorithmRZ.cpp
@@ -286,7 +286,8 @@ void GalileanPsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldData
 }
 
 void
-GalileanPsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
+GalileanPsatdAlgorithmRZ::CurrentCorrection (const int lev,
+                                             SpectralFieldDataRZ& field_data,
                                              std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                              const std::unique_ptr<amrex::MultiFab>& rho )
 {
@@ -296,11 +297,11 @@ GalileanPsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
     using Idx = SpectralFieldIndex;
 
     // Forward Fourier transform of J and rho
-    field_data.ForwardTransform( *current[0], Idx::Jx,
+    field_data.ForwardTransform( lev, *current[0], Idx::Jx,
                                  *current[1], Idx::Jy);
-    field_data.ForwardTransform( *current[2], Idx::Jz, 0);
-    field_data.ForwardTransform( *rho, Idx::rho_old, 0 );
-    field_data.ForwardTransform( *rho, Idx::rho_new, 1 );
+    field_data.ForwardTransform( lev, *current[2], Idx::Jz, 0);
+    field_data.ForwardTransform( lev, *rho, Idx::rho_old, 0 );
+    field_data.ForwardTransform( lev, *rho, Idx::rho_new, 1 );
 
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
@@ -364,14 +365,16 @@ GalileanPsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform( *current[0], Idx::Jx,
+    field_data.BackwardTransform( lev,
+                                  *current[0], Idx::Jx,
                                   *current[1], Idx::Jy);
-    field_data.BackwardTransform( *current[2], Idx::Jz, 0 );
+    field_data.BackwardTransform( lev, *current[2], Idx::Jz, 0 );
 
 }
 
 void
-GalileanPsatdAlgorithmRZ::VayDeposition (SpectralFieldDataRZ& /*field_data*/,
+GalileanPsatdAlgorithmRZ::VayDeposition (const int /*lev*/,
+                                         SpectralFieldDataRZ& /*field_data*/,
                                          std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
 {
     amrex::Abort("Vay deposition not implemented in RZ geometry");
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
index a4cccbc55..368ea27ad 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.H
@@ -46,7 +46,8 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
          *                           the three components of the current density
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection (SpectralFieldData& field_data,
+        virtual void CurrentCorrection (const int lev,
+                                        SpectralFieldData& field_data,
                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                         const std::unique_ptr<amrex::MultiFab>& rho) override final;
 
@@ -61,7 +62,8 @@ class PMLPsatdAlgorithm : public SpectralBaseAlgorithm
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldData& field_data,
+        virtual void VayDeposition (const int lev,
+                                    SpectralFieldData& field_data,
                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
     private:
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
index 6b9397807..aafa3ba8d 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PMLPsatdAlgorithm.cpp
@@ -226,7 +226,8 @@ void PMLPsatdAlgorithm::InitializeSpectralCoefficients (
 }
 
 void
-PMLPsatdAlgorithm::CurrentCorrection (SpectralFieldData& /*field_data*/,
+PMLPsatdAlgorithm::CurrentCorrection (const int /*lev*/,
+                                      SpectralFieldData& /*field_data*/,
                                       std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/,
                                       const std::unique_ptr<amrex::MultiFab>& /*rho*/)
 {
@@ -234,7 +235,8 @@ PMLPsatdAlgorithm::CurrentCorrection (SpectralFieldData& /*field_data*/,
 }
 
 void
-PMLPsatdAlgorithm::VayDeposition (SpectralFieldData& /*field_data*/,
+PMLPsatdAlgorithm::VayDeposition (const int /*lev*/,
+                                  SpectralFieldData& /*field_data*/,
                                   std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
 {
     amrex::Abort("Vay deposition not implemented for PML PSATD");
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index 1ca4aaf31..d5efe5b74 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -62,6 +62,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
         virtual void CurrentCorrection (
+            const int lev,
             SpectralFieldData& field_data,
             std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
             const std::unique_ptr<amrex::MultiFab>& rho) override final;
@@ -78,6 +79,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          *                           the three components of the current density
          */
         virtual void VayDeposition (
+            const int lev,
             SpectralFieldData& field_data,
             std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 0a278a54a..379bb6369 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -881,6 +881,7 @@ void PsatdAlgorithm::InitializeSpectralCoefficients (
 
 void
 PsatdAlgorithm::CurrentCorrection (
+    const int lev,
     SpectralFieldData& field_data,
     std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
     const std::unique_ptr<amrex::MultiFab>& rho)
@@ -891,11 +892,11 @@ PsatdAlgorithm::CurrentCorrection (
     using Idx = SpectralFieldIndex;
 
     // Forward Fourier transform of J and rho
-    field_data.ForwardTransform(*current[0], Idx::Jx, 0);
-    field_data.ForwardTransform(*current[1], Idx::Jy, 0);
-    field_data.ForwardTransform(*current[2], Idx::Jz, 0);
-    field_data.ForwardTransform(*rho, Idx::rho_old, 0);
-    field_data.ForwardTransform(*rho, Idx::rho_new, 1);
+    field_data.ForwardTransform(lev, *current[0], Idx::Jx, 0);
+    field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0);
+    field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0);
+    field_data.ForwardTransform(lev, *rho, Idx::rho_old, 0);
+    field_data.ForwardTransform(lev, *rho, Idx::rho_new, 1);
 
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
@@ -989,13 +990,14 @@ PsatdAlgorithm::CurrentCorrection (
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform(*current[0], Idx::Jx, 0);
-    field_data.BackwardTransform(*current[1], Idx::Jy, 0);
-    field_data.BackwardTransform(*current[2], Idx::Jz, 0);
+    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0);
+    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0);
+    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0);
 }
 
 void
 PsatdAlgorithm::VayDeposition (
+    const int lev,
     SpectralFieldData& field_data,
     std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
 {
@@ -1007,9 +1009,9 @@ PsatdAlgorithm::VayDeposition (
     // Forward Fourier transform of D (temporarily stored in current):
     // D is nodal and does not match the staggering of J, therefore we pass the
     // actual staggering of D (IntVect(1)) to the ForwardTransform function
-    field_data.ForwardTransform(*current[0], Idx::Jx, 0, IntVect(1));
-    field_data.ForwardTransform(*current[1], Idx::Jy, 0, IntVect(1));
-    field_data.ForwardTransform(*current[2], Idx::Jz, 0, IntVect(1));
+    field_data.ForwardTransform(lev, *current[0], Idx::Jx, 0, IntVect(1));
+    field_data.ForwardTransform(lev, *current[1], Idx::Jy, 0, IntVect(1));
+    field_data.ForwardTransform(lev, *current[2], Idx::Jz, 0, IntVect(1));
 
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi)
@@ -1065,9 +1067,9 @@ PsatdAlgorithm::VayDeposition (
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform(*current[0], Idx::Jx, 0);
-    field_data.BackwardTransform(*current[1], Idx::Jy, 0);
-    field_data.BackwardTransform(*current[2], Idx::Jz, 0);
+    field_data.BackwardTransform(lev, *current[0], Idx::Jx, 0);
+    field_data.BackwardTransform(lev, *current[1], Idx::Jy, 0);
+    field_data.BackwardTransform(lev, *current[2], Idx::Jz, 0);
 }
 
 #endif // WARPX_USE_PSATD
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
index 74bf71c29..9d6ed9cfa 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
@@ -41,7 +41,8 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
          *                           the three components of the current density
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection (SpectralFieldDataRZ& field_data,
+        virtual void CurrentCorrection (const int lev,
+                                        SpectralFieldDataRZ& field_data,
                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                         const std::unique_ptr<amrex::MultiFab>& rho) override final;
 
@@ -56,7 +57,8 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldDataRZ& field_data,
+        virtual void VayDeposition (const int lev,
+                                    SpectralFieldDataRZ& field_data,
                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current) override final;
 
     private:
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
index 002f9e55f..b96d5fa93 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
@@ -215,7 +215,8 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
 }
 
 void
-PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
+PsatdAlgorithmRZ::CurrentCorrection (const int lev,
+                                     SpectralFieldDataRZ& field_data,
                                      std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                      const std::unique_ptr<amrex::MultiFab>& rho)
 {
@@ -225,11 +226,12 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
     using Idx = SpectralFieldIndex;
 
     // Forward Fourier transform of J and rho
-    field_data.ForwardTransform( *current[0], Idx::Jx,
+    field_data.ForwardTransform( lev,
+                                 *current[0], Idx::Jx,
                                  *current[1], Idx::Jy);
-    field_data.ForwardTransform( *current[2], Idx::Jz, 0);
-    field_data.ForwardTransform( *rho, Idx::rho_old, 0 );
-    field_data.ForwardTransform( *rho, Idx::rho_new, 1 );
+    field_data.ForwardTransform( lev, *current[2], Idx::Jz, 0);
+    field_data.ForwardTransform( lev, *rho, Idx::rho_old, 0 );
+    field_data.ForwardTransform( lev, *rho, Idx::rho_new, 1 );
 
     // Loop over boxes
     for (amrex::MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
@@ -289,13 +291,16 @@ PsatdAlgorithmRZ::CurrentCorrection (SpectralFieldDataRZ& field_data,
     }
 
     // Backward Fourier transform of J
-    field_data.BackwardTransform( *current[0], Idx::Jx,
+    field_data.BackwardTransform( lev,
+                                  *current[0], Idx::Jx,
                                   *current[1], Idx::Jy);
-    field_data.BackwardTransform( *current[2], Idx::Jz, 0 );
+    field_data.BackwardTransform( lev,
+                                  *current[2], Idx::Jz, 0 );
 }
 
 void
-PsatdAlgorithmRZ::VayDeposition (SpectralFieldDataRZ& /*field_data*/,
+PsatdAlgorithmRZ::VayDeposition (const int lev /**/,
+                                 SpectralFieldDataRZ& /*field_data*/,
                                  std::array<std::unique_ptr<amrex::MultiFab>,3>& /*current*/)
 {
     amrex::Abort("Vay deposition not implemented in RZ geometry");
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
index 3e909145d..598046042 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.H
@@ -42,7 +42,8 @@ class SpectralBaseAlgorithm
          *                           the three components of the current density
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection (SpectralFieldData& field_data,
+        virtual void CurrentCorrection (const int lev,
+                                        SpectralFieldData& field_data,
                                         std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                         const std::unique_ptr<amrex::MultiFab>& rho) = 0;
 
@@ -55,13 +56,15 @@ class SpectralBaseAlgorithm
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldData& field_data,
+        virtual void VayDeposition (const int lev,
+                                    SpectralFieldData& field_data,
                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0;
 
         /**
          * \brief Compute spectral divergence of E
          */
-        void ComputeSpectralDivE ( SpectralFieldData& field_data,
+        void ComputeSpectralDivE ( const int lev,
+                                   SpectralFieldData& field_data,
                                    const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
                                    amrex::MultiFab& divE );
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
index 46e9b50f7..46fe83900 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithm.cpp
@@ -14,6 +14,7 @@ using namespace amrex;
  */
 void
 SpectralBaseAlgorithm::ComputeSpectralDivE (
+    const int lev,
     SpectralFieldData& field_data,
     const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
     amrex::MultiFab& divE )
@@ -21,9 +22,9 @@ SpectralBaseAlgorithm::ComputeSpectralDivE (
     using Idx = SpectralFieldIndex;
 
     // Forward Fourier transform of E
-    field_data.ForwardTransform( *Efield[0], Idx::Ex, 0 );
-    field_data.ForwardTransform( *Efield[1], Idx::Ey, 0 );
-    field_data.ForwardTransform( *Efield[2], Idx::Ez, 0 );
+    field_data.ForwardTransform(lev, *Efield[0], Idx::Ex, 0 );
+    field_data.ForwardTransform(lev, *Efield[1], Idx::Ey, 0 );
+    field_data.ForwardTransform(lev, *Efield[2], Idx::Ez, 0 );
 
     // Loop over boxes
     for (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
@@ -64,5 +65,5 @@ SpectralBaseAlgorithm::ComputeSpectralDivE (
     }
 
     // Backward Fourier transform
-    field_data.BackwardTransform( divE, Idx::divE, 0 );
+    field_data.BackwardTransform(lev, divE, Idx::divE, 0 );
 }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
index 29a8d8e7f..c19d2fc5e 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.H
@@ -38,14 +38,15 @@ class SpectralBaseAlgorithmRZ
          *                           the three components of the current density
          * \param[in]     rho        Unique pointer to \c MultiFab storing the charge density
          */
-        virtual void CurrentCorrection ( SpectralFieldDataRZ& field_data,
+        virtual void CurrentCorrection ( const int lev, SpectralFieldDataRZ& field_data,
                                          std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                          const std::unique_ptr<amrex::MultiFab>& rho ) = 0;
 
         /**
          * \brief Compute spectral divergence of E
          */
-        void ComputeSpectralDivE ( SpectralFieldDataRZ& field_data,
+        void ComputeSpectralDivE ( const int lev,
+                                   SpectralFieldDataRZ& field_data,
                                    const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
                                    amrex::MultiFab& divE );
 
@@ -58,7 +59,8 @@ class SpectralBaseAlgorithmRZ
          * \param[in,out] current    Array of unique pointers to \c MultiFab storing
          *                           the three components of the current density
          */
-        virtual void VayDeposition (SpectralFieldDataRZ& field_data,
+        virtual void VayDeposition (const int lev,
+                                    SpectralFieldDataRZ& field_data,
                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current) = 0;
 
     protected: // Meant to be used in the subclasses
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp
index 4609b7177..f80f54208 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/SpectralBaseAlgorithmRZ.cpp
@@ -15,6 +15,7 @@ using namespace amrex;
  */
 void
 SpectralBaseAlgorithmRZ::ComputeSpectralDivE (
+    const int lev,
     SpectralFieldDataRZ& field_data,
     const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
     amrex::MultiFab& divE )
@@ -23,9 +24,11 @@ SpectralBaseAlgorithmRZ::ComputeSpectralDivE (
     using Idx = SpectralFieldIndex;
 
     // Forward Fourier transform of E
-    field_data.ForwardTransform( *Efield[0], Idx::Ex,
+    field_data.ForwardTransform( lev,
+                                 *Efield[0], Idx::Ex,
                                  *Efield[1], Idx::Ey );
-    field_data.ForwardTransform( *Efield[2], Idx::Ez, 0 );
+    field_data.ForwardTransform( lev,
+                                 *Efield[2], Idx::Ez, 0 );
 
     // Loop over boxes
     for (MFIter mfi(field_data.fields); mfi.isValid(); ++mfi){
@@ -69,5 +72,5 @@ SpectralBaseAlgorithmRZ::ComputeSpectralDivE (
     }
 
     // Backward Fourier transform
-    field_data.BackwardTransform( divE, Idx::divE, 0 );
+    field_data.BackwardTransform( lev, divE, Idx::divE, 0 );
 }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index 4990f9926..764ecc8af 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -47,24 +47,27 @@ class SpectralFieldData
 {
 
     public:
-        SpectralFieldData( const amrex::BoxArray& realspace_ba,
+        SpectralFieldData( const int lev,
+                           const amrex::BoxArray& realspace_ba,
                            const SpectralKSpace& k_space,
                            const amrex::DistributionMapping& dm,
                            const int n_field_required,
-                           const bool periodic_single_box );
+                           const bool periodic_single_box);
         SpectralFieldData() = default; // Default constructor
         SpectralFieldData& operator=(SpectralFieldData&& field_data) = default;
         ~SpectralFieldData();
 
-        void ForwardTransform (const amrex::MultiFab& mf, const int field_index,
+        void ForwardTransform (const int lev,
+                               const amrex::MultiFab& mf, const int field_index,
                                const int i_comp, const amrex::IntVect& stag);
         AMREX_FORCE_INLINE
-        void ForwardTransform (const amrex::MultiFab& mf, const int field_index, const int i_comp)
+        void ForwardTransform (const int lev,
+                               const amrex::MultiFab& mf, const int field_index, const int i_comp)
         {
-            ForwardTransform(mf, field_index, i_comp, mf.ixType().toIntVect());
+            ForwardTransform(lev, mf, field_index, i_comp, mf.ixType().toIntVect());
         }
 
-        void BackwardTransform (amrex::MultiFab& mf, const int field_index, const int i_comp);
+        void BackwardTransform (const int lev, amrex::MultiFab& mf, const int field_index, const int i_comp);
 
         // `fields` stores fields in spectral space, as multicomponent FabArray
         SpectralField fields;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index e99e7d57b..9c4e786de 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -6,6 +6,7 @@
  * License: BSD-3-Clause-LBNL
  */
 #include "SpectralFieldData.H"
+#include "WarpX.H"
 
 #include <map>
 
@@ -14,12 +15,15 @@
 using namespace amrex;
 
 /* \brief Initialize fields in spectral space, and FFT plans */
-SpectralFieldData::SpectralFieldData( const amrex::BoxArray& realspace_ba,
+SpectralFieldData::SpectralFieldData( const int lev,
+                                      const amrex::BoxArray& realspace_ba,
                                       const SpectralKSpace& k_space,
                                       const amrex::DistributionMapping& dm,
                                       const int n_field_required,
-                                      const bool periodic_single_box )
+                                      const bool periodic_single_box)
 {
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
     m_periodic_single_box = periodic_single_box;
 
     const BoxArray& spectralspace_ba = k_space.spectralspace_ba;
@@ -62,6 +66,12 @@ SpectralFieldData::SpectralFieldData( const amrex::BoxArray& realspace_ba,
     // Loop over boxes and allocate the corresponding plan
     // for each box owned by the local MPI proc
     for ( MFIter mfi(spectralspace_ba, dm); mfi.isValid(); ++mfi ){
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
+
         // Note: the size of the real-space box and spectral-space box
         // differ when using real-to-complex FFT. When initializing
         // the FFT plan, the valid dimensions are those of the real-space box.
@@ -76,6 +86,13 @@ SpectralFieldData::SpectralFieldData( const amrex::BoxArray& realspace_ba,
             fft_size, tmpRealField[mfi].dataPtr(),
             reinterpret_cast<AnyFFT::Complex*>( tmpSpectralField[mfi].dataPtr()),
             AnyFFT::direction::C2R, AMREX_SPACEDIM);
+
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     }
 }
 
@@ -94,9 +111,12 @@ SpectralFieldData::~SpectralFieldData()
  *  to spectral space, and store the corresponding result internally
  *  (in the spectral field specified by `field_index`) */
 void
-SpectralFieldData::ForwardTransform (const MultiFab& mf, const int field_index,
+SpectralFieldData::ForwardTransform (const int lev,
+                     const MultiFab& mf, const int field_index,
                                      const int i_comp, const IntVect& stag)
 {
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
     // Check field index type, in order to apply proper shift in spectral space
     const bool is_nodal_x = (stag[0] == amrex::IndexType::NODE) ? true : false;
 #if (AMREX_SPACEDIM == 3)
@@ -108,6 +128,11 @@ SpectralFieldData::ForwardTransform (const MultiFab& mf, const int field_index,
 
     // Loop over boxes
     for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
 
         // Copy the real-space field `mf` to the temporary field `tmpRealField`
         // This ensures that all fields have the same number of points
@@ -164,6 +189,13 @@ SpectralFieldData::ForwardTransform (const MultiFab& mf, const int field_index,
                 fields_arr(i,j,k,field_index) = spectral_field_value;
             });
         }
+
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     }
 }
 
@@ -171,10 +203,13 @@ SpectralFieldData::ForwardTransform (const MultiFab& mf, const int field_index,
 /* \brief Transform spectral field specified by `field_index` back to
  * real space, and store it in the component `i_comp` of `mf` */
 void
-SpectralFieldData::BackwardTransform( MultiFab& mf,
+SpectralFieldData::BackwardTransform( const int lev,
+                                      MultiFab& mf,
                                       const int field_index,
                                       const int i_comp )
 {
+    amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
+
     // Check field index type, in order to apply proper shift in spectral space
     const bool is_nodal_x = mf.is_nodal(0);
 #if (AMREX_SPACEDIM == 3)
@@ -186,6 +221,11 @@ SpectralFieldData::BackwardTransform( MultiFab& mf,
 
     // Loop over boxes
     for ( MFIter mfi(mf); mfi.isValid(); ++mfi ){
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+        }
+        Real wt = amrex::second();
 
         // Copy the spectral-space field `tmpSpectralField` to the appropriate
         // field (specified by the input argument field_index)
@@ -262,6 +302,13 @@ SpectralFieldData::BackwardTransform( MultiFab& mf,
             }
 
         }
+
+        if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
+        {
+            amrex::Gpu::synchronize();
+            wt = amrex::second() - wt;
+            amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
+        }
     }
 }
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H
index 0dff6da7d..8484d7174 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.H
@@ -36,23 +36,23 @@ class SpectralFieldDataRZ
 
         using BinomialFilter = amrex::LayoutData<SpectralBinomialFilter>;
 
-        SpectralFieldDataRZ (const amrex::BoxArray& realspace_ba,
+        SpectralFieldDataRZ (const int lev,
+                             const amrex::BoxArray& realspace_ba,
                              const SpectralKSpaceRZ& k_space,
                              const amrex::DistributionMapping& dm,
                              const int n_field_required,
-                             const int n_modes,
-                             const int lev);
+                             const int n_modes);
         SpectralFieldDataRZ () = default; // Default constructor
         SpectralFieldDataRZ& operator=(SpectralFieldDataRZ&& field_data) = default;
         ~SpectralFieldDataRZ ();
 
-        void ForwardTransform (const amrex::MultiFab& mf, const int field_index,
+        void ForwardTransform (const int lev, const amrex::MultiFab& mf, const int field_index,
                                const int i_comp=0);
-        void ForwardTransform (const amrex::MultiFab& mf_r, const int field_index_r,
+        void ForwardTransform (const int lev, const amrex::MultiFab& mf_r, const int field_index_r,
                                const amrex::MultiFab& mf_t, const int field_index_t);
-        void BackwardTransform (amrex::MultiFab& mf, const int field_index,
+        void BackwardTransform (const int lev, amrex::MultiFab& mf, const int field_index,
                                 const int i_comp=0);
-        void BackwardTransform (amrex::MultiFab& mf_r, const int field_index_r,
+        void BackwardTransform (const int lev, amrex::MultiFab& mf_r, const int field_index_r,
                                 amrex::MultiFab& mf_t, const int field_index_t);
 
         void FABZForwardTransform (amrex::MFIter const & mfi, amrex::Box const & realspace_bx,
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp
index c2d8dfe94..645d605c8 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldDataRZ.cpp
@@ -19,12 +19,12 @@ using amrex::operator""_rt;
  * \param n_field_required Specifies the number of fields that will be transformed
  * \param n_modes Number of cylindrical modes
  * */
-SpectralFieldDataRZ::SpectralFieldDataRZ (amrex::BoxArray const & realspace_ba,
+SpectralFieldDataRZ::SpectralFieldDataRZ (const int lev,
+                                          amrex::BoxArray const & realspace_ba,
                                           SpectralKSpaceRZ const & k_space,
                                           amrex::DistributionMapping const & dm,
                                           int const n_field_required,
-                                          int const n_modes,
-                                          int const lev)
+                                          int const n_modes)
     : n_rz_azimuthal_modes(n_modes)
 {
     amrex::BoxArray const & spectralspace_ba = k_space.spectralspace_ba;
@@ -394,7 +394,8 @@ SpectralFieldDataRZ::FABZBackwardTransform (amrex::MFIter const & mfi, amrex::Bo
  *  to spectral space, and store the corresponding result internally
  *  (in the spectral field specified by `field_index`) */
 void
-SpectralFieldDataRZ::ForwardTransform (amrex::MultiFab const & field_mf, int const field_index,
+SpectralFieldDataRZ::ForwardTransform (const int lev,
+                                       amrex::MultiFab const & field_mf, int const field_index,
                                        int const i_comp)
 {
     // Check field index type, in order to apply proper shift in spectral space.
@@ -440,7 +441,8 @@ SpectralFieldDataRZ::ForwardTransform (amrex::MultiFab const & field_mf, int con
  *  to spectral space, and store the corresponding result internally
  *  (in the spectral fields specified by `field_index_r` and `field_index_t`) */
 void
-SpectralFieldDataRZ::ForwardTransform (amrex::MultiFab const & field_mf_r, int const field_index_r,
+SpectralFieldDataRZ::ForwardTransform (const int lev,
+                                       amrex::MultiFab const & field_mf_r, int const field_index_r,
                                        amrex::MultiFab const & field_mf_t, int const field_index_t)
 {
     // Check field index type, in order to apply proper shift in spectral space.
@@ -488,7 +490,8 @@ SpectralFieldDataRZ::ForwardTransform (amrex::MultiFab const & field_mf_r, int c
 /* \brief Transform spectral field specified by `field_index` back to
  * real space, and store it in the component `i_comp` of `field_mf` */
 void
-SpectralFieldDataRZ::BackwardTransform (amrex::MultiFab& field_mf, int const field_index,
+SpectralFieldDataRZ::BackwardTransform (const int lev,
+                                        amrex::MultiFab& field_mf, int const field_index,
                                         int const i_comp)
 {
     // Check field index type, in order to apply proper shift in spectral space.
@@ -522,7 +525,8 @@ SpectralFieldDataRZ::BackwardTransform (amrex::MultiFab& field_mf, int const fie
 /* \brief Transform spectral fields specified by `field_index_r` and
  * `field_index_t` back to real space, and store them in `field_mf_r` and `field_mf_t` */
 void
-SpectralFieldDataRZ::BackwardTransform (amrex::MultiFab& field_mf_r, int const field_index_r,
+SpectralFieldDataRZ::BackwardTransform (const int lev,
+                                        amrex::MultiFab& field_mf_r, int const field_index_r,
                                         amrex::MultiFab& field_mf_t, int const field_index_t)
 {
     // Check field index type, in order to apply proper shift in spectral space.
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 2626c6982..c9a7bbcc4 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -28,13 +28,15 @@ class SpectralSolver
         // underlying classes `SpectralFieldData` and `PsatdAlgorithm`
 
         // Constructor
-        SpectralSolver( const amrex::BoxArray& realspace_ba,
+        SpectralSolver( const int lev,
+                        const amrex::BoxArray& realspace_ba,
                         const amrex::DistributionMapping& dm,
                         const int norder_x, const int norder_y,
                         const int norder_z, const bool nodal,
                         const amrex::Array<amrex::Real,3>& v_galilean,
                         const amrex::Array<amrex::Real,3>& v_comoving,
-                        const amrex::RealVect dx, const amrex::Real dt,
+                        const amrex::RealVect dx,
+                        const amrex::Real dt,
                         const bool pml=false,
                         const bool periodic_single_box=false,
                         const bool update_with_rho=false,
@@ -44,7 +46,8 @@ class SpectralSolver
          * \brief Transform the component `i_comp` of MultiFab `mf`
          *  to spectral space, and store the corresponding result internally
          *  (in the spectral field specified by `field_index`) */
-        void ForwardTransform( const amrex::MultiFab& mf,
+        void ForwardTransform( const int lev,
+                               const amrex::MultiFab& mf,
                                const int field_index,
                                const int i_comp=0 );
 
@@ -52,7 +55,8 @@ class SpectralSolver
          * \brief Transform spectral field specified by `field_index` back to
          * real space, and store it in the component `i_comp` of `mf`
          */
-        void BackwardTransform( amrex::MultiFab& mf,
+        void BackwardTransform( const int lev,
+                                amrex::MultiFab& mf,
                                 const int field_index,
                                 const int i_comp=0 );
 
@@ -65,9 +69,10 @@ class SpectralSolver
           * \brief Public interface to call the member function ComputeSpectralDivE
           * of the base class SpectralBaseAlgorithm from objects of class SpectralSolver
           */
-        void ComputeSpectralDivE ( const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
+        void ComputeSpectralDivE ( const int lev,
+                                   const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
                                    amrex::MultiFab& divE ) {
-            algorithm->ComputeSpectralDivE( field_data, Efield, divE );
+            algorithm->ComputeSpectralDivE( lev, field_data, Efield, divE );
         }
 
         /**
@@ -80,9 +85,10 @@ class SpectralSolver
          *                        storing the three components of the current density
          * \param[in]     rho     unique pointer to MultiFab storing the charge density
          */
-        void CurrentCorrection ( std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+        void CurrentCorrection ( const int lev,
+                                 std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                  const std::unique_ptr<amrex::MultiFab>& rho ) {
-             algorithm->CurrentCorrection( field_data, current, rho );
+            algorithm->CurrentCorrection( lev, field_data, current, rho );
         }
 
         /**
@@ -94,9 +100,9 @@ class SpectralSolver
          * \param[in,out] current Array of unique pointers to \c MultiFab storing
          *                        the three components of the current density
          */
-        void VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+        void VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
         {
-            algorithm->VayDeposition(field_data, current);
+            algorithm->VayDeposition(lev, field_data, current);
         }
 
     private:
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index 725cfcf92..ec1219f93 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -33,6 +33,7 @@
  * \param periodic_single_box Whether the full simulation domain consists of a single periodic box (i.e. the global domain is not MPI parallelized)
  */
 SpectralSolver::SpectralSolver(
+                const int lev,
                 const amrex::BoxArray& realspace_ba,
                 const amrex::DistributionMapping& dm,
                 const int norder_x, const int norder_y,
@@ -72,27 +73,29 @@ SpectralSolver::SpectralSolver(
     }
 
     // - Initialize arrays for fields in spectral space + FFT plans
-    field_data = SpectralFieldData( realspace_ba, k_space, dm,
-            algorithm->getRequiredNumberOfFields(), periodic_single_box );
+    field_data = SpectralFieldData( lev, realspace_ba, k_space, dm,
+                    algorithm->getRequiredNumberOfFields(), periodic_single_box);
 
 }
 
 void
-SpectralSolver::ForwardTransform( const amrex::MultiFab& mf,
+SpectralSolver::ForwardTransform( const int lev,
+                                  const amrex::MultiFab& mf,
                                   const int field_index,
                                   const int i_comp )
 {
     WARPX_PROFILE("SpectralSolver::ForwardTransform");
-    field_data.ForwardTransform( mf, field_index, i_comp );
+    field_data.ForwardTransform( lev, mf, field_index, i_comp );
 }
 
 void
-SpectralSolver::BackwardTransform( amrex::MultiFab& mf,
+SpectralSolver::BackwardTransform( const int lev,
+                                   amrex::MultiFab& mf,
                                    const int field_index,
                                    const int i_comp )
 {
     WARPX_PROFILE("SpectralSolver::BackwardTransform");
-    field_data.BackwardTransform( mf, field_index, i_comp );
+    field_data.BackwardTransform( lev, mf, field_index, i_comp );
 }
 
 void
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
index 8b6ee51d9..7fe4bd00c 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
@@ -25,35 +25,35 @@ class SpectralSolverRZ
         // underlying classes `SpectralFieldData` and `PsatdAlgorithm`
 
         // Constructor
-        SpectralSolverRZ (amrex::BoxArray const & realspace_ba,
+        SpectralSolverRZ (const int lev,
+                          amrex::BoxArray const & realspace_ba,
                           amrex::DistributionMapping const & dm,
                           int const n_rz_azimuthal_modes,
                           int const norder_z, bool const nodal,
                           const amrex::Array<amrex::Real,3>& v_galilean,
                           amrex::RealVect const dx, amrex::Real const dt,
-                          int const lev,
                           bool const update_with_rho);
 
         /* \brief Transform the component `i_comp` of MultiFab `field_mf`
          *  to spectral space, and store the corresponding result internally
          *  (in the spectral field specified by `field_index`) */
-        void ForwardTransform (amrex::MultiFab const & field_mf, int const field_index,
+        void ForwardTransform (const int lev, amrex::MultiFab const & field_mf, int const field_index,
                                int const i_comp=0);
 
         /* \brief Transform the two MultiFabs `field_mf1` and `field_mf2`
          *  to spectral space, and store the corresponding results internally
          *  (in the spectral field specified by `field_index1` and `field_index2`) */
-        void ForwardTransform (amrex::MultiFab const & field_mf1, int const field_index1,
+        void ForwardTransform (const int lev, amrex::MultiFab const & field_mf1, int const field_index1,
                                amrex::MultiFab const & field_mf2, int const field_index2);
 
         /* \brief Transform spectral field specified by `field_index` back to
          * real space, and store it in the component `i_comp` of `field_mf` */
-        void BackwardTransform (amrex::MultiFab& field_mf, int const field_index,
+        void BackwardTransform (const int lev, amrex::MultiFab& field_mf, int const field_index,
                                 int const i_comp=0);
 
         /* \brief Transform spectral fields specified by `field_index1` and `field_index2`
          * back to real space, and store it in `field_mf1` and `field_mf2`*/
-        void BackwardTransform (amrex::MultiFab& field_mf1, int const field_index1,
+        void BackwardTransform (const int lev, amrex::MultiFab& field_mf1, int const field_index1,
                                 amrex::MultiFab& field_mf2, int const field_index2);
 
         /* \brief Update the fields in spectral space, over one timestep */
@@ -82,7 +82,7 @@ class SpectralSolverRZ
           * \brief Public interface to call the member function ComputeSpectralDivE
           * of the base class SpectralBaseAlgorithmRZ from objects of class SpectralSolverRZ
           */
-        void ComputeSpectralDivE (const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
+        void ComputeSpectralDivE (const int lev, const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
                                   amrex::MultiFab& divE);
 
         /**
@@ -95,7 +95,7 @@ class SpectralSolverRZ
          *                        storing the three components of the current density
          * \param[in]     rho     unique pointer to MultiFab storing the charge density
          */
-        void CurrentCorrection (std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+        void CurrentCorrection (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
                                 const std::unique_ptr<amrex::MultiFab>& rho);
 
         /**
@@ -107,7 +107,7 @@ class SpectralSolverRZ
          * \param[in,out] current Array of unique pointers to \c MultiFab storing
          * the three components of the current density
          */
-        void VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current);
+        void VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current);
 
     private:
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
index 820db5a12..f59b99752 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
@@ -25,13 +25,13 @@
  * \param pml      Whether the boxes in which the solver is applied are PML boxes
  *                 PML is not supported.
  */
-SpectralSolverRZ::SpectralSolverRZ (amrex::BoxArray const & realspace_ba,
+SpectralSolverRZ::SpectralSolverRZ (const int lev,
+                                    amrex::BoxArray const & realspace_ba,
                                     amrex::DistributionMapping const & dm,
                                     int const n_rz_azimuthal_modes,
                                     int const norder_z, bool const nodal,
                                     const amrex::Array<amrex::Real,3>& v_galilean,
                                     amrex::RealVect const dx, amrex::Real const dt,
-                                    int const lev,
                                     bool const update_with_rho)
     : k_space(realspace_ba, dm, dx)
 {
@@ -55,48 +55,54 @@ SpectralSolverRZ::SpectralSolverRZ (amrex::BoxArray const & realspace_ba,
     }
 
     // - Initialize arrays for fields in spectral space + FFT plans
-    field_data = SpectralFieldDataRZ(realspace_ba, k_space, dm,
+    field_data = SpectralFieldDataRZ(lev, realspace_ba, k_space, dm,
                                      algorithm->getRequiredNumberOfFields(),
-                                     n_rz_azimuthal_modes, lev);
+                                     n_rz_azimuthal_modes);
 }
 
 /* \brief Transform the component `i_comp` of MultiFab `field_mf`
  *  to spectral space, and store the corresponding result internally
  *  (in the spectral field specified by `field_index`) */
 void
-SpectralSolverRZ::ForwardTransform (amrex::MultiFab const & field_mf, int const field_index,
+SpectralSolverRZ::ForwardTransform (const int lev,
+                                    amrex::MultiFab const & field_mf, int const field_index,
                                     int const i_comp) {
     WARPX_PROFILE("SpectralSolverRZ::ForwardTransform");
-    field_data.ForwardTransform(field_mf, field_index, i_comp);
+    field_data.ForwardTransform(lev, field_mf, field_index, i_comp);
 }
 
 /* \brief Transform the two MultiFabs `field_mf1` and `field_mf2`
  *  to spectral space, and store the corresponding results internally
  *  (in the spectral field specified by `field_index1` and `field_index2`) */
 void
-SpectralSolverRZ::ForwardTransform (amrex::MultiFab const & field_mf1, int const field_index1,
+SpectralSolverRZ::ForwardTransform (const int lev,
+                                    amrex::MultiFab const & field_mf1, int const field_index1,
                                     amrex::MultiFab const & field_mf2, int const field_index2) {
     WARPX_PROFILE("SpectralSolverRZ::ForwardTransform");
-    field_data.ForwardTransform(field_mf1, field_index1,
+    field_data.ForwardTransform(lev,
+                                field_mf1, field_index1,
                                 field_mf2, field_index2);
 }
 
 /* \brief Transform spectral field specified by `field_index` back to
  * real space, and store it in the component `i_comp` of `field_mf` */
 void
-SpectralSolverRZ::BackwardTransform (amrex::MultiFab& field_mf, int const field_index,
+SpectralSolverRZ::BackwardTransform (const int lev,
+                                     amrex::MultiFab& field_mf, int const field_index,
                                      int const i_comp) {
     WARPX_PROFILE("SpectralSolverRZ::BackwardTransform");
-    field_data.BackwardTransform(field_mf, field_index, i_comp);
+    field_data.BackwardTransform(lev, field_mf, field_index, i_comp);
 }
 
 /* \brief Transform spectral fields specified by `field_index1` and `field_index2`
  * back to real space, and store it in `field_mf1` and `field_mf2`*/
 void
-SpectralSolverRZ::BackwardTransform (amrex::MultiFab& field_mf1, int const field_index1,
+SpectralSolverRZ::BackwardTransform (const int lev,
+                                     amrex::MultiFab& field_mf1, int const field_index1,
                                      amrex::MultiFab& field_mf2, int const field_index2) {
     WARPX_PROFILE("SpectralSolverRZ::BackwardTransform");
-    field_data.BackwardTransform(field_mf1, field_index1,
+    field_data.BackwardTransform(lev,
+                                 field_mf1, field_index1,
                                  field_mf2, field_index2);
 }
 
@@ -115,9 +121,10 @@ SpectralSolverRZ::pushSpectralFields () {
   * of the base class SpectralBaseAlgorithmRZ from objects of class SpectralSolverRZ
   */
 void
-SpectralSolverRZ::ComputeSpectralDivE (const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
+SpectralSolverRZ::ComputeSpectralDivE (const int lev,
+                                       const std::array<std::unique_ptr<amrex::MultiFab>,3>& Efield,
                                        amrex::MultiFab& divE) {
-    algorithm->ComputeSpectralDivE(field_data, Efield, divE);
+    algorithm->ComputeSpectralDivE(lev, field_data, Efield, divE);
 }
 
 /**
@@ -131,13 +138,14 @@ SpectralSolverRZ::ComputeSpectralDivE (const std::array<std::unique_ptr<amrex::M
  * \param[in]     rho     unique pointer to MultiFab storing the charge density
  */
 void
-SpectralSolverRZ::CurrentCorrection (std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
-                                      const std::unique_ptr<amrex::MultiFab>& rho) {
-     algorithm->CurrentCorrection(field_data, current, rho);
+SpectralSolverRZ::CurrentCorrection (const int lev,
+                                     std::array<std::unique_ptr<amrex::MultiFab>,3>& current,
+                                     const std::unique_ptr<amrex::MultiFab>& rho) {
+    algorithm->CurrentCorrection(lev, field_data, current, rho);
 }
 
 void
-SpectralSolverRZ::VayDeposition (std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
+SpectralSolverRZ::VayDeposition (const int lev, std::array<std::unique_ptr<amrex::MultiFab>,3>& current)
 {
-  algorithm->VayDeposition(field_data, current);
+    algorithm->VayDeposition(lev, field_data, current);
 }
diff --git a/Source/FieldSolver/WarpXPushFieldsEM.cpp b/Source/FieldSolver/WarpXPushFieldsEM.cpp
index 71521ffbe..1e138b4c1 100644
--- a/Source/FieldSolver/WarpXPushFieldsEM.cpp
+++ b/Source/FieldSolver/WarpXPushFieldsEM.cpp
@@ -27,6 +27,7 @@ using namespace amrex;
 namespace {
     void
     PushPSATDSinglePatch (
+        const int lev,
 #ifdef WARPX_DIM_RZ
         SpectralSolverRZ& solver,
 #else
@@ -47,33 +48,36 @@ namespace {
 
         // Perform forward Fourier transform
 #ifdef WARPX_DIM_RZ
-        solver.ForwardTransform(*Efield[0], Idx::Ex,
+        solver.ForwardTransform(lev,
+                                *Efield[0], Idx::Ex,
                                 *Efield[1], Idx::Ey);
 #else
-        solver.ForwardTransform(*Efield[0], Idx::Ex);
-        solver.ForwardTransform(*Efield[1], Idx::Ey);
+        solver.ForwardTransform(lev, *Efield[0], Idx::Ex);
+        solver.ForwardTransform(lev, *Efield[1], Idx::Ey);
 #endif
-        solver.ForwardTransform(*Efield[2], Idx::Ez);
+        solver.ForwardTransform(lev, *Efield[2], Idx::Ez);
 #ifdef WARPX_DIM_RZ
-        solver.ForwardTransform(*Bfield[0], Idx::Bx,
+        solver.ForwardTransform(lev,
+                                *Bfield[0], Idx::Bx,
                                 *Bfield[1], Idx::By);
 #else
-        solver.ForwardTransform(*Bfield[0], Idx::Bx);
-        solver.ForwardTransform(*Bfield[1], Idx::By);
+        solver.ForwardTransform(lev, *Bfield[0], Idx::Bx);
+        solver.ForwardTransform(lev, *Bfield[1], Idx::By);
 #endif
-        solver.ForwardTransform(*Bfield[2], Idx::Bz);
+        solver.ForwardTransform(lev, *Bfield[2], Idx::Bz);
 #ifdef WARPX_DIM_RZ
-        solver.ForwardTransform(*current[0], Idx::Jx,
+        solver.ForwardTransform(lev,
+                                *current[0], Idx::Jx,
                                 *current[1], Idx::Jy);
 #else
-        solver.ForwardTransform(*current[0], Idx::Jx);
-        solver.ForwardTransform(*current[1], Idx::Jy);
+        solver.ForwardTransform(lev, *current[0], Idx::Jx);
+        solver.ForwardTransform(lev, *current[1], Idx::Jy);
 #endif
-        solver.ForwardTransform(*current[2], Idx::Jz);
+        solver.ForwardTransform(lev, *current[2], Idx::Jz);
 
         if (rho) {
-            solver.ForwardTransform(*rho, Idx::rho_old, 0);
-            solver.ForwardTransform(*rho, Idx::rho_new, 1);
+            solver.ForwardTransform(lev, *rho, Idx::rho_old, 0);
+            solver.ForwardTransform(lev, *rho, Idx::rho_new, 1);
         }
 #ifdef WARPX_DIM_RZ
         if (WarpX::use_kspace_filter) {
@@ -86,31 +90,33 @@ namespace {
         solver.pushSpectralFields();
         // Perform backward Fourier Transform
 #ifdef WARPX_DIM_RZ
-        solver.BackwardTransform(*Efield[0], Idx::Ex,
+        solver.BackwardTransform(lev,
+                                 *Efield[0], Idx::Ex,
                                  *Efield[1], Idx::Ey);
 #else
-        solver.BackwardTransform(*Efield[0], Idx::Ex);
-        solver.BackwardTransform(*Efield[1], Idx::Ey);
+        solver.BackwardTransform(lev, *Efield[0], Idx::Ex);
+        solver.BackwardTransform(lev, *Efield[1], Idx::Ey);
 #endif
-        solver.BackwardTransform(*Efield[2], Idx::Ez);
+        solver.BackwardTransform(lev, *Efield[2], Idx::Ez);
 #ifdef WARPX_DIM_RZ
-        solver.BackwardTransform(*Bfield[0], Idx::Bx,
+        solver.BackwardTransform(lev,
+                                 *Bfield[0], Idx::Bx,
                                  *Bfield[1], Idx::By);
 #else
-        solver.BackwardTransform(*Bfield[0], Idx::Bx);
-        solver.BackwardTransform(*Bfield[1], Idx::By);
+        solver.BackwardTransform(lev, *Bfield[0], Idx::Bx);
+        solver.BackwardTransform(lev, *Bfield[1], Idx::By);
 #endif
-        solver.BackwardTransform(*Bfield[2], Idx::Bz);
+        solver.BackwardTransform(lev, *Bfield[2], Idx::Bz);
 
 #ifndef WARPX_DIM_RZ
         if (WarpX::fft_do_time_averaging){
-            solver.BackwardTransform(*Efield_avg[0], Idx::Ex_avg);
-            solver.BackwardTransform(*Efield_avg[1], Idx::Ey_avg);
-            solver.BackwardTransform(*Efield_avg[2], Idx::Ez_avg);
+            solver.BackwardTransform(lev, *Efield_avg[0], Idx::Ex_avg);
+            solver.BackwardTransform(lev, *Efield_avg[1], Idx::Ey_avg);
+            solver.BackwardTransform(lev, *Efield_avg[2], Idx::Ez_avg);
 
-            solver.BackwardTransform(*Bfield_avg[0], Idx::Bx_avg);
-            solver.BackwardTransform(*Bfield_avg[1], Idx::By_avg);
-            solver.BackwardTransform(*Bfield_avg[2], Idx::Bz_avg);
+            solver.BackwardTransform(lev, *Bfield_avg[0], Idx::Bx_avg);
+            solver.BackwardTransform(lev, *Bfield_avg[1], Idx::By_avg);
+            solver.BackwardTransform(lev, *Bfield_avg[2], Idx::Bz_avg);
         }
 #endif
     }
@@ -131,7 +137,7 @@ WarpX::PushPSATD (amrex::Real a_dt)
 
         // Evolve the fields in the PML boxes
         if (do_pml && pml[lev]->ok()) {
-            pml[lev]->PushPSATD();
+            pml[lev]->PushPSATD(lev);
         }
     }
 #endif
@@ -149,10 +155,10 @@ WarpX::PushPSATD (int lev, amrex::Real /* dt */) {
                                          "WarpX::PushPSATD: only supported for PSATD solver.");
     }
     // Update the fields on the fine and coarse patch
-    PushPSATDSinglePatch( *spectral_solver_fp[lev],
+    PushPSATDSinglePatch( lev, *spectral_solver_fp[lev],
         Efield_fp[lev], Bfield_fp[lev], Efield_avg_fp[lev], Bfield_avg_fp[lev], current_fp[lev], rho_fp[lev] );
     if (spectral_solver_cp[lev]) {
-        PushPSATDSinglePatch( *spectral_solver_cp[lev],
+        PushPSATDSinglePatch( lev, *spectral_solver_cp[lev],
              Efield_cp[lev], Bfield_cp[lev], Efield_avg_cp[lev], Bfield_avg_cp[lev], current_cp[lev], rho_cp[lev] );
     }
     if (use_damp_fields_in_z_guard) {
@@ -186,9 +192,9 @@ WarpX::EvolveB (int lev, PatchType patch_type, amrex::Real a_dt)
 
     // Evolve B field in regular cells
     if (patch_type == PatchType::fine) {
-        m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], a_dt );
+        m_fdtd_solver_fp[lev]->EvolveB( Bfield_fp[lev], Efield_fp[lev], lev, a_dt );
     } else {
-        m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], a_dt );
+        m_fdtd_solver_cp[lev]->EvolveB( Bfield_cp[lev], Efield_cp[lev], lev, a_dt );
     }
 
     // Evolve B field in PML cells
@@ -230,10 +236,10 @@ WarpX::EvolveE (int lev, PatchType patch_type, amrex::Real a_dt)
     // Evolve E field in regular cells
     if (patch_type == PatchType::fine) {
         m_fdtd_solver_fp[lev]->EvolveE( Efield_fp[lev], Bfield_fp[lev],
-                                      current_fp[lev], F_fp[lev], a_dt );
+                    current_fp[lev], F_fp[lev], lev, a_dt );
     } else {
         m_fdtd_solver_cp[lev]->EvolveE( Efield_cp[lev], Bfield_cp[lev],
-                                      current_cp[lev], F_cp[lev], a_dt );
+                    current_cp[lev], F_cp[lev], lev, a_dt );
     }
 
     // Evolve E field in PML cells