PSATD: More Options for Time Dependency of J, Rho (#3242)

* Rename PsatdAlgorithm as PsatdAlgorithmJConstantInTime

* Add New Inputs: psatd.J_in_time, psatd.rho_in_time

* Update PSATD Classes

* Fix Bug for FDTD Build

* Fix Warning for RZ PSATD Build

* Fix Multi-J Loop w/ J Constant in Time

* Clean up

* Fix Error Message

* Fix Time of Charge Deposition for Rho Linear

* Add 3D Langmuir CI Tests w/ Multi-J Algo

* Fix Checksums of New CI Tests

* Remove Extra CI Tests (added in #3363)

* Fix CI Tests w/ Multi-J PSATD

* Add Docs for New Input Parameters

* Galilean/Comoving PSATD Not Compatible w/ J Linear

12 files changed, 81 insertions, 58 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
index eb88b1f4e..a370d4b2d 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/CMakeLists.txt
@@ -1,6 +1,6 @@
 target_sources(WarpX
   PRIVATE
-    PsatdAlgorithm.cpp
+    PsatdAlgorithmJConstantInTime.cpp
     PsatdAlgorithmJLinearInTime.cpp
     PsatdAlgorithmPml.cpp
     SpectralBaseAlgorithm.cpp
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
index 7c2bf4281..c798ffb01 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/Make.package
@@ -1,5 +1,5 @@
 CEXE_sources += SpectralBaseAlgorithm.cpp
-CEXE_sources += PsatdAlgorithm.cpp
+CEXE_sources += PsatdAlgorithmJConstantInTime.cpp
 CEXE_sources += PsatdAlgorithmJLinearInTime.cpp
 CEXE_sources += PsatdAlgorithmPml.cpp
 CEXE_sources += PsatdAlgorithmComoving.cpp
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H
index dd9c6a7fd..0d2d67434 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H
@@ -4,8 +4,8 @@
  *
  * License: BSD-3-Clause-LBNL
  */
-#ifndef WARPX_PSATD_ALGORITHM_H_
-#define WARPX_PSATD_ALGORITHM_H_
+#ifndef WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_
+#define WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_
 
 #include "FieldSolver/SpectralSolver/SpectralFieldData.H"
 #include "FieldSolver/SpectralSolver/SpectralKSpace.H"
@@ -24,12 +24,12 @@
 /* \brief Class that updates the field in spectral space
  * and stores the coefficients of the corresponding update equation.
  */
-class PsatdAlgorithm : public SpectralBaseAlgorithm
+class PsatdAlgorithmJConstantInTime : public SpectralBaseAlgorithm
 {
     public:
 
         /**
-         * \brief Constructor of the class PsatdAlgorithm
+         * \brief Constructor of the class PsatdAlgorithmJConstantInTime
          *
          * \param[in] spectral_kspace spectral space
          * \param[in] dm distribution mapping
@@ -45,7 +45,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          * \param[in] dive_cleaning Update F as part of the field update, so that errors in divE=rho propagate away at the speed of light
          * \param[in] divb_cleaning Update G as part of the field update, so that errors in divB=0 propagate away at the speed of light
          */
-        PsatdAlgorithm (
+        PsatdAlgorithmJConstantInTime (
             const SpectralKSpace& spectral_kspace,
             const amrex::DistributionMapping& dm,
             const SpectralFieldIndex& spectral_index,
@@ -142,4 +142,4 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
         bool m_is_galilean;
 };
 #endif // WARPX_USE_PSATD
-#endif // WARPX_PSATD_ALGORITHM_H_
+#endif // WARPX_PSATD_ALGORITHM_J_CONSTANT_IN_TIME_H_
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp
index 1cbc27f0b..8971061f6 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmJConstantInTime.cpp
@@ -4,7 +4,7 @@
  *
  * License: BSD-3-Clause-LBNL
  */
-#include "PsatdAlgorithm.H"
+#include "PsatdAlgorithmJConstantInTime.H"
 
 #include "Utils/TextMsg.H"
 #include "Utils/WarpXConst.H"
@@ -27,7 +27,7 @@
 
 using namespace amrex;
 
-PsatdAlgorithm::PsatdAlgorithm(
+PsatdAlgorithmJConstantInTime::PsatdAlgorithmJConstantInTime(
     const SpectralKSpace& spectral_kspace,
     const DistributionMapping& dm,
     const SpectralFieldIndex& spectral_index,
@@ -110,7 +110,7 @@ PsatdAlgorithm::PsatdAlgorithm(
 }
 
 void
-PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
+PsatdAlgorithmJConstantInTime::pushSpectralFields (SpectralFieldData& f) const
 {
     const bool update_with_rho = m_update_with_rho;
     const bool time_averaging  = m_time_averaging;
@@ -340,7 +340,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
     }
 }
 
-void PsatdAlgorithm::InitializeSpectralCoefficients (
+void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficients (
     const SpectralKSpace& spectral_kspace,
     const amrex::DistributionMapping& dm,
     const amrex::Real dt)
@@ -542,7 +542,7 @@ void PsatdAlgorithm::InitializeSpectralCoefficients (
     }
 }
 
-void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging (
+void PsatdAlgorithmJConstantInTime::InitializeSpectralCoefficientsAveraging (
     const SpectralKSpace& spectral_kspace,
     const amrex::DistributionMapping& dm,
     const amrex::Real dt)
@@ -733,10 +733,10 @@ void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging (
     }
 }
 
-void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data)
+void PsatdAlgorithmJConstantInTime::CurrentCorrection (SpectralFieldData& field_data)
 {
     // Profiling
-    BL_PROFILE("PsatdAlgorithm::CurrentCorrection");
+    BL_PROFILE("PsatdAlgorithmJConstantInTime::CurrentCorrection");
 
     const SpectralFieldIndex& Idx = m_spectral_index;
 
@@ -833,10 +833,10 @@ void PsatdAlgorithm::CurrentCorrection (SpectralFieldData& field_data)
 }
 
 void
-PsatdAlgorithm::VayDeposition (SpectralFieldData& field_data)
+PsatdAlgorithmJConstantInTime::VayDeposition (SpectralFieldData& field_data)
 {
     // Profiling
-    BL_PROFILE("PsatdAlgorithm::VayDeposition()");
+    BL_PROFILE("PsatdAlgorithmJConstantInTime::VayDeposition()");
 
     const SpectralFieldIndex& Idx = m_spectral_index;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
index 608da5fd5..097a1a9d6 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
@@ -23,7 +23,8 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
                          bool const nodal, amrex::Real const dt_step,
                          bool const update_with_rho,
                          const bool time_averaging,
-                         const bool do_multi_J,
+                         const int J_in_time,
+                         const int rho_in_time,
                          const bool dive_cleaning,
                          const bool divb_cleaning);
         // Redefine functions from base class
@@ -62,7 +63,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
         amrex::Real m_dt;
         bool m_update_with_rho;
         bool m_time_averaging;
-        bool m_do_multi_J;
+        int m_J_in_time;
         bool m_dive_cleaning;
         bool m_divb_cleaning;
         SpectralRealCoefficients C_coef, S_ck_coef, X1_coef, X2_coef, X3_coef;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
index fec0ebc8f..55b58821c 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
@@ -6,6 +6,7 @@
  */
 #include "PsatdAlgorithmRZ.H"
 #include "Utils/TextMsg.H"
+#include "Utils/WarpXAlgorithmSelection.H"
 #include "Utils/WarpXConst.H"
 #include "Utils/WarpXProfilerWrapper.H"
 #include "WarpX.H"
@@ -23,7 +24,8 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
                                     bool const nodal, amrex::Real const dt,
                                     bool const update_with_rho,
                                     const bool time_averaging,
-                                    const bool do_multi_J,
+                                    const int J_in_time,
+                                    const int rho_in_time,
                                     const bool dive_cleaning,
                                     const bool divb_cleaning)
      // Initialize members of base class
@@ -32,10 +34,11 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
        m_dt(dt),
        m_update_with_rho(update_with_rho),
        m_time_averaging(time_averaging),
-       m_do_multi_J(do_multi_J),
+       m_J_in_time(J_in_time),
        m_dive_cleaning(dive_cleaning),
        m_divb_cleaning(divb_cleaning)
 {
+    amrex::ignore_unused(rho_in_time);
 
     // Allocate the arrays of coefficients
     amrex::BoxArray const & ba = spectral_kspace.spectralspace_ba;
@@ -47,28 +50,28 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
 
     coefficients_initialized = false;
 
-    if (time_averaging && do_multi_J)
+    if (time_averaging && J_in_time == JInTime::Linear)
     {
         X5_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
         X6_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
     }
 
-    if (time_averaging && !do_multi_J)
+    if (time_averaging && J_in_time != JInTime::Linear)
     {
         amrex::Abort(Utils::TextMsg::Err(
-            "RZ PSATD: psatd.do_time_averaging = 1 implemented only with warpx.do_multi_J = 1"));
+            "RZ PSATD: psatd.do_time_averaging=1 implemented only with psatd.J_in_time=linear"));
     }
 
-    if (dive_cleaning && !do_multi_J)
+    if (dive_cleaning && J_in_time != JInTime::Linear)
     {
         amrex::Abort(Utils::TextMsg::Err(
-            "RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with warpx.do_multi_J = 1"));
+            "RZ PSATD: warpx.do_dive_cleaning=1 implemented only with psatd.J_in_time=linear"));
     }
 
-    if (divb_cleaning && !do_multi_J)
+    if (divb_cleaning && J_in_time != JInTime::Linear)
     {
         amrex::Abort(Utils::TextMsg::Err(
-            "RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with warpx.do_multi_J = 1"));
+            "RZ PSATD: warpx.do_divb_cleaning=1 implemented only with psatd.J_in_time=linear"));
     }
 }
 
@@ -80,7 +83,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 
     const bool update_with_rho = m_update_with_rho;
     const bool time_averaging = m_time_averaging;
-    const bool do_multi_J = m_do_multi_J;
+    const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false;
     const bool dive_cleaning = m_dive_cleaning;
     const bool divb_cleaning = m_divb_cleaning;
 
@@ -109,7 +112,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 
         amrex::Array4<const amrex::Real> X5_arr;
         amrex::Array4<const amrex::Real> X6_arr;
-        if (time_averaging && do_multi_J)
+        if (time_averaging && J_in_time_linear)
         {
             X5_arr = X5_coef[mfi].array();
             X6_arr = X6_coef[mfi].array();
@@ -235,7 +238,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
                 G_old = fields(i,j,k,G_m);
             }
 
-            if (do_multi_J)
+            if (J_in_time_linear)
             {
                 const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode;
                 const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode;
@@ -332,7 +335,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f)
 {
     const bool time_averaging = m_time_averaging;
-    const bool do_multi_J = m_do_multi_J;
+    const bool J_in_time_linear = (m_J_in_time == JInTime::Linear) ? true : false;
 
     // Fill them with the right values:
     // Loop over boxes and allocate the corresponding coefficients
@@ -353,7 +356,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
 
         amrex::Array4<amrex::Real> X5;
         amrex::Array4<amrex::Real> X6;
-        if (time_averaging && do_multi_J)
+        if (time_averaging && J_in_time_linear)
         {
             X5 = X5_coef[mfi].array();
             X6 = X6_coef[mfi].array();
@@ -392,7 +395,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
                 X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0);
             }
 
-            if (time_averaging && do_multi_J)
+            if (time_averaging && J_in_time_linear)
             {
                 constexpr amrex::Real c2 = PhysConst::c;
                 const amrex::Real dt3 = dt * dt * dt;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
index 811fdd4d7..4ab88f1a3 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -55,7 +55,8 @@ class SpectralFieldIndex
          */
         SpectralFieldIndex (const bool update_with_rho,
                             const bool time_averaging,
-                            const bool do_multi_J,
+                            const int J_in_time,
+                            const int rho_in_time,
                             const bool dive_cleaning,
                             const bool divb_cleaning,
                             const bool pml,
@@ -89,8 +90,13 @@ class SpectralFieldIndex
     int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1;
     int Bx_avg = -1, By_avg = -1, Bz_avg = -1;
 
-    // Multi-J, div(E) and div(B) cleaning
+    // J linear in time
     int Jx_new = -1, Jy_new = -1, Jz_new = -1;
+
+    // rho quadratic in time
+    int rho_mid = -1;
+
+    // div(E) and div(B) cleaning
     int F = -1, G = -1;
 
     // PML
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index 00346a2c7..56189e0ff 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -35,7 +35,8 @@ using namespace amrex;
 
 SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
                                         const bool time_averaging,
-                                        const bool do_multi_J,
+                                        const int J_in_time,
+                                        const int rho_in_time,
                                         const bool dive_cleaning,
                                         const bool divb_cleaning,
                                         const bool pml,
@@ -68,13 +69,18 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
 
         if (divb_cleaning) G = c++;
 
-        if (do_multi_J)
+        if (J_in_time == JInTime::Linear)
         {
             Jx_new = c++;
             Jy_new = c++;
             Jz_new = c++;
         }
 
+        if (rho_in_time == RhoInTime::Quadratic)
+        {
+            rho_mid = c++;
+        }
+
         if (pml_rz)
         {
             Er_pml = c++; Et_pml = c++;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 684cf9586..38b242010 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -58,9 +58,10 @@ class SpectralSolver
          *                                (no domain decomposition)
          * \param[in] update_with_rho whether rho is used in the field update equations
          * \param[in] fft_do_time_averaging whether the time averaging algorithm is used
-         * \param[in] do_multi_J whether the multi-J algorithm is used (hence two currents
-         *                       computed at the beginning and the end of the time interval
-         *                       instead of one current computed at half time)
+         * \param[in] J_in_time integer that corresponds to the time dependency of J
+         *                      (constant, linear) for the PSATD algorithm
+         * \param[in] rho_in_time integer that corresponds to the time dependency of rho
+         *                        (linear, quadratic) for the PSATD algorithm
          * \param[in] dive_cleaning whether to use div(E) cleaning to account for errors in
          *                          Gauss law (new field F in the update equations)
          * \param[in] divb_cleaning whether to use div(B) cleaning to account for errors in
@@ -79,7 +80,8 @@ class SpectralSolver
                         const bool periodic_single_box,
                         const bool update_with_rho,
                         const bool fft_do_time_averaging,
-                        const bool do_multi_J,
+                        const int J_in_time,
+                        const int rho_in_time,
                         const bool dive_cleaning,
                         const bool divb_cleaning);
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index 75c82319c..c0d7b8941 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -8,7 +8,7 @@
 #include "FieldSolver/SpectralSolver/SpectralFieldData.H"
 #include "SpectralAlgorithms/PsatdAlgorithmComoving.H"
 #include "SpectralAlgorithms/PsatdAlgorithmPml.H"
-#include "SpectralAlgorithms/PsatdAlgorithm.H"
+#include "SpectralAlgorithms/PsatdAlgorithmJConstantInTime.H"
 #include "SpectralAlgorithms/PsatdAlgorithmJLinearInTime.H"
 #include "SpectralKSpace.H"
 #include "SpectralSolver.H"
@@ -30,7 +30,8 @@ SpectralSolver::SpectralSolver(
                 const bool pml, const bool periodic_single_box,
                 const bool update_with_rho,
                 const bool fft_do_time_averaging,
-                const bool do_multi_J,
+                const int J_in_time,
+                const int rho_in_time,
                 const bool dive_cleaning,
                 const bool divb_cleaning)
 {
@@ -41,8 +42,9 @@ SpectralSolver::SpectralSolver(
     // as well as the value of the corresponding k coordinates)
     const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx);
 
-    m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging,
-                                          do_multi_J, dive_cleaning, divb_cleaning, pml);
+    m_spectral_index = SpectralFieldIndex(
+        update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time,
+        dive_cleaning, divb_cleaning, pml);
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
@@ -64,18 +66,18 @@ SpectralSolver::SpectralSolver(
         }
         else // PSATD algorithms: standard, Galilean, averaged Galilean, multi-J
         {
-            if (do_multi_J)
+            if (J_in_time == JInTime::Constant)
             {
-                algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>(
+                algorithm = std::make_unique<PsatdAlgorithmJConstantInTime>(
                     k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
-                    dt, fft_do_time_averaging, dive_cleaning, divb_cleaning);
+                    v_galilean, dt, update_with_rho, fft_do_time_averaging,
+                    dive_cleaning, divb_cleaning);
             }
-            else // standard, Galilean, averaged Galilean
+            else // J linear in time
             {
-                algorithm = std::make_unique<PsatdAlgorithm>(
+                algorithm = std::make_unique<PsatdAlgorithmJLinearInTime>(
                     k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal,
-                    v_galilean, dt, update_with_rho, fft_do_time_averaging,
-                    dive_cleaning, divb_cleaning);
+                    dt, fft_do_time_averaging, dive_cleaning, divb_cleaning);
             }
         }
     }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
index 30c922512..45b55d9d4 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
@@ -37,7 +37,8 @@ class SpectralSolverRZ
                           bool const with_pml,
                           bool const update_with_rho,
                           const bool fft_do_time_averaging,
-                          const bool do_multi_J,
+                          const int J_in_time,
+                          const int rho_in_time,
                           const bool dive_cleaning,
                           const bool divb_cleaning);
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
index 3f4bddd8b..23d93fe04 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
@@ -35,7 +35,8 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
                                     bool const with_pml,
                                     bool const update_with_rho,
                                     const bool fft_do_time_averaging,
-                                    const bool do_multi_J,
+                                    const int J_in_time,
+                                    const int rho_in_time,
                                     const bool dive_cleaning,
                                     const bool divb_cleaning)
     : k_space(realspace_ba, dm, dx)
@@ -47,8 +48,9 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
     //   as well as the value of the corresponding k coordinates.
 
     const bool is_pml = false;
-    m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging,
-                                          do_multi_J, dive_cleaning, divb_cleaning, is_pml, with_pml);
+    m_spectral_index = SpectralFieldIndex(
+        update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time,
+        dive_cleaning, divb_cleaning, is_pml, with_pml);
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
@@ -60,7 +62,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
          // v_galilean is 0: use standard PSATD algorithm
         algorithm = std::make_unique<PsatdAlgorithmRZ>(
             k_space, dm, m_spectral_index, n_rz_azimuthal_modes, norder_z, nodal, dt,
-            update_with_rho, fft_do_time_averaging, do_multi_J, dive_cleaning, divb_cleaning);
+            update_with_rho, fft_do_time_averaging, J_in_time, rho_in_time, dive_cleaning, divb_cleaning);
     } else {
         // Otherwise: use the Galilean algorithm
         algorithm = std::make_unique<PsatdAlgorithmGalileanRZ>(