Multi-J: J Always Linear in Time (#2679)

* Remove WarpX::J_linear_in_time

* Fix Bug

* Implement Deposition of Rho at Half Time

* Implement New Equations for E

* Fix New Equations for E

* Cleaning

* Fix Limits of X2, X3, X7

* Implement New Equations for F

* Implement New Equations for F

* Revert E Equations to Avoid Roundoff Changes

* Remove Quadratic Equations

* Remove Deposition of Rho at Half Time Step

10 files changed, 61 insertions, 58 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
index 41053fc9e..5549eda07 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.H
@@ -43,8 +43,9 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
          * \param[in] dt time step of the simulation
          * \param[in] update_with_rho whether the update equation for E uses rho or not
          * \param[in] time_averaging whether to use time averaging for large time steps
-         * \param[in] J_linear_in_time whether to use two currents computed at the beginning and the end
-         *            of the time interval (instead of using one current computed at half time)
+         * \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] 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
          */
@@ -61,7 +62,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
             const amrex::Real dt,
             const bool update_with_rho,
             const bool time_averaging,
-            const bool J_linear_in_time,
+            const bool do_multi_J,
             const bool dive_cleaning,
             const bool divb_cleaning);
 
@@ -171,7 +172,7 @@ class PsatdAlgorithm : public SpectralBaseAlgorithm
         amrex::Real m_dt;
         bool m_update_with_rho;
         bool m_time_averaging;
-        bool m_J_linear_in_time;
+        bool m_do_multi_J;
         bool m_dive_cleaning;
         bool m_divb_cleaning;
         bool m_is_galilean;
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index fe9562dba..2b1e09c2d 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -39,7 +39,7 @@ PsatdAlgorithm::PsatdAlgorithm(
     const amrex::Real dt,
     const bool update_with_rho,
     const bool time_averaging,
-    const bool J_linear_in_time,
+    const bool do_multi_J,
     const bool dive_cleaning,
     const bool divb_cleaning)
     // Initializer list
@@ -59,7 +59,7 @@ PsatdAlgorithm::PsatdAlgorithm(
     m_dt(dt),
     m_update_with_rho(update_with_rho),
     m_time_averaging(time_averaging),
-    m_J_linear_in_time(J_linear_in_time),
+    m_do_multi_J(do_multi_J),
     m_dive_cleaning(dive_cleaning),
     m_divb_cleaning(divb_cleaning)
 {
@@ -84,7 +84,7 @@ PsatdAlgorithm::PsatdAlgorithm(
     InitializeSpectralCoefficients(spectral_kspace, dm, dt);
 
     // Allocate these coefficients only with time averaging
-    if (time_averaging && !J_linear_in_time)
+    if (time_averaging && !do_multi_J)
     {
         Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
         Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
@@ -95,8 +95,8 @@ PsatdAlgorithm::PsatdAlgorithm(
         InitializeSpectralCoefficientsAveraging(spectral_kspace, dm, dt);
     }
     // Allocate these coefficients only with time averaging
-    // and with the assumption that J is linear in time
-    else if (time_averaging && J_linear_in_time)
+    // and with the assumption that J is linear in time (always with multi-J algorithm)
+    else if (time_averaging && do_multi_J)
     {
         X5_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
         X6_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
@@ -123,7 +123,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
 {
     const bool update_with_rho  = m_update_with_rho;
     const bool time_averaging   = m_time_averaging;
-    const bool J_linear_in_time = m_J_linear_in_time;
+    const bool do_multi_J = m_do_multi_J;
     const bool dive_cleaning    = m_dive_cleaning;
     const bool divb_cleaning    = m_divb_cleaning;
     const bool is_galilean      = m_is_galilean;
@@ -163,7 +163,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
         amrex::Array4<const Complex> Y3_arr;
         amrex::Array4<const Complex> Y4_arr;
 
-        if (time_averaging && !J_linear_in_time)
+        if (time_averaging && !do_multi_J)
         {
             Psi1_arr = Psi1_coef[mfi].array();
             Psi2_arr = Psi2_coef[mfi].array();
@@ -175,7 +175,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
 
         Array4<const Complex> X5_arr;
         Array4<const Complex> X6_arr;
-        if (time_averaging && J_linear_in_time)
+        if (time_averaging && do_multi_J)
         {
             X5_arr = X5_coef[mfi].array();
             X6_arr = X6_coef[mfi].array();
@@ -320,7 +320,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
                 fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B;
             }
 
-            if (J_linear_in_time)
+            if (do_multi_J)
             {
                 const Complex Jx_new = fields(i,j,k,Idx.Jx_new);
                 const Complex Jy_new = fields(i,j,k,Idx.Jy_new);
@@ -391,7 +391,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
             }
 
             // Additional update equations for averaged Galilean algorithm
-            if (time_averaging && !J_linear_in_time)
+            if (time_averaging && !do_multi_J)
             {
                 // These coefficients are initialized in the function InitializeSpectralCoefficients below
                 const Complex Psi1 = Psi1_arr(i,j,k);
diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
index a23afbf4c..fba8f0881 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.H
@@ -23,7 +23,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
                          bool const nodal, amrex::Real const dt_step,
                          bool const update_with_rho,
                          const bool time_averaging,
-                         const bool J_linear_in_time,
+                         const bool do_multi_J,
                          const bool dive_cleaning,
                          const bool divb_cleaning);
         // Redefine functions from base class
@@ -74,7 +74,7 @@ class PsatdAlgorithmRZ : public SpectralBaseAlgorithmRZ
         amrex::Real m_dt;
         bool m_update_with_rho;
         bool m_time_averaging;
-        bool m_J_linear_in_time;
+        bool m_do_multi_J;
         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 60a30e6f6..efc9337af 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithmRZ.cpp
@@ -22,7 +22,7 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
                                     bool const nodal, amrex::Real const dt,
                                     bool const update_with_rho,
                                     const bool time_averaging,
-                                    const bool J_linear_in_time,
+                                    const bool do_multi_J,
                                     const bool dive_cleaning,
                                     const bool divb_cleaning)
      // Initialize members of base class
@@ -31,7 +31,7 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
        m_dt(dt),
        m_update_with_rho(update_with_rho),
        m_time_averaging(time_averaging),
-       m_J_linear_in_time(J_linear_in_time),
+       m_do_multi_J(do_multi_J),
        m_dive_cleaning(dive_cleaning),
        m_divb_cleaning(divb_cleaning)
 {
@@ -46,25 +46,25 @@ PsatdAlgorithmRZ::PsatdAlgorithmRZ (SpectralKSpaceRZ const & spectral_kspace,
 
     coefficients_initialized = false;
 
-    if (time_averaging && J_linear_in_time)
+    if (time_averaging && do_multi_J)
     {
         X5_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
         X6_coef = SpectralRealCoefficients(ba, dm, n_rz_azimuthal_modes, 0);
     }
 
-    if (time_averaging && !J_linear_in_time)
+    if (time_averaging && !do_multi_J)
     {
-        amrex::Abort("RZ PSATD: psatd.do_time_averaging = 1 implemented only with psatd.J_linear_in_time = 1");
+        amrex::Abort("RZ PSATD: psatd.do_time_averaging = 1 implemented only with warpx.do_multi_J = 1");
     }
 
-    if (dive_cleaning && !J_linear_in_time)
+    if (dive_cleaning && !do_multi_J)
     {
-        amrex::Abort("RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with psatd.J_linear_in_time = 1");
+        amrex::Abort("RZ PSATD: warpx.do_dive_cleaning = 1 implemented only with warpx.do_multi_J = 1");
     }
 
-    if (divb_cleaning && !J_linear_in_time)
+    if (divb_cleaning && !do_multi_J)
     {
-        amrex::Abort("RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with psatd.J_linear_in_time = 1");
+        amrex::Abort("RZ PSATD: warpx.do_divb_cleaning = 1 implemented only with warpx.do_multi_J = 1");
     }
 }
 
@@ -76,7 +76,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 
     const bool update_with_rho = m_update_with_rho;
     const bool time_averaging = m_time_averaging;
-    const bool J_linear_in_time = m_J_linear_in_time;
+    const bool do_multi_J = m_do_multi_J;
     const bool dive_cleaning = m_dive_cleaning;
     const bool divb_cleaning = m_divb_cleaning;
 
@@ -105,7 +105,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 
         amrex::Array4<const amrex::Real> X5_arr;
         amrex::Array4<const amrex::Real> X6_arr;
-        if (time_averaging && J_linear_in_time)
+        if (time_averaging && do_multi_J)
         {
             X5_arr = X5_coef[mfi].array();
             X6_arr = X6_coef[mfi].array();
@@ -231,7 +231,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
                 G_old = fields(i,j,k,G_m);
             }
 
-            if (J_linear_in_time)
+            if (do_multi_J)
             {
                 const int Jp_m_new = Idx.Jx_new + Idx.n_fields*mode;
                 const int Jm_m_new = Idx.Jy_new + Idx.n_fields*mode;
@@ -328,7 +328,7 @@ PsatdAlgorithmRZ::pushSpectralFields(SpectralFieldDataRZ & f)
 void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const & f)
 {
     const bool time_averaging = m_time_averaging;
-    const bool J_linear_in_time = m_J_linear_in_time;
+    const bool do_multi_J = m_do_multi_J;
 
     // Fill them with the right values:
     // Loop over boxes and allocate the corresponding coefficients
@@ -349,7 +349,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
 
         amrex::Array4<amrex::Real> X5;
         amrex::Array4<amrex::Real> X6;
-        if (time_averaging && J_linear_in_time)
+        if (time_averaging && do_multi_J)
         {
             X5 = X5_coef[mfi].array();
             X6 = X6_coef[mfi].array();
@@ -388,7 +388,7 @@ void PsatdAlgorithmRZ::InitializeSpectralCoefficients (SpectralFieldDataRZ const
                 X3(i,j,k,mode) = - c*c * dt*dt / (3._rt*ep0);
             }
 
-            if (time_averaging && J_linear_in_time)
+            if (time_averaging && do_multi_J)
             {
                 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 16f8e179c..9b748a048 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.H
@@ -39,21 +39,21 @@ class SpectralFieldIndex
          * Set integer indices to access data in spectral space
          * and total number of fields to be stored.
          *
-         * \param[in] update_with_rho   whether rho is used in the field update equations
-         * \param[in] time_averaging    whether the time averaging algorithm is used
-         * \param[in] J_linear_in_time  whether to use two currents computed at the beginning and
-         *                              the end of the time interval (instead of using one current
-         *                              computed at half time)
-         * \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
-         *                              div(B) = 0 law (new field G in the update equations)
-         * \param[in] pml               whether the indices are used to access spectral data
-         *                              for the PML spectral solver
+         * \param[in] update_with_rho whether rho is used in the field update equations
+         * \param[in] 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] 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
+         *                            div(B) = 0 law (new field G in the update equations)
+         * \param[in] pml             whether the indices are used to access spectral data
+         *                            for the PML spectral solver
          */
         SpectralFieldIndex (const bool update_with_rho,
                             const bool time_averaging,
-                            const bool J_linear_in_time,
+                            const bool do_multi_J,
                             const bool dive_cleaning,
                             const bool divb_cleaning,
                             const bool pml);
@@ -86,7 +86,7 @@ class SpectralFieldIndex
     int Ex_avg = -1, Ey_avg = -1, Ez_avg = -1;
     int Bx_avg = -1, By_avg = -1, Bz_avg = -1;
 
-    // J linear in time
+    // Multi-J, div(E) and div(B) cleaning
     int Jx_new = -1, Jy_new = -1, Jz_new = -1;
     int F = -1, G = -1;
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index 9d8563b6b..86fe8f7bd 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -34,7 +34,7 @@ using namespace amrex;
 
 SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
                                         const bool time_averaging,
-                                        const bool J_linear_in_time,
+                                        const bool do_multi_J,
                                         const bool dive_cleaning,
                                         const bool divb_cleaning,
                                         const bool pml)
@@ -66,9 +66,11 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
 
         if (divb_cleaning) G = c++;
 
-        if (J_linear_in_time)
+        if (do_multi_J)
         {
-            Jx_new = c++; Jy_new = c++; Jz_new = c++;
+            Jx_new = c++;
+            Jy_new = c++;
+            Jz_new = c++;
         }
     }
     else // PML
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.H b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
index 5c3036838..f5cfb814b 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.H
@@ -59,9 +59,9 @@ 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] J_linear_in_time whether to use two currents computed at the beginning and
-         *                             the end of the time interval (instead of using one current
-         *                             computed at half time)
+         * \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] 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
@@ -81,7 +81,7 @@ class SpectralSolver
                         const bool periodic_single_box,
                         const bool update_with_rho,
                         const bool fft_do_time_averaging,
-                        const bool J_linear_in_time,
+                        const bool do_multi_J,
                         const bool dive_cleaning,
                         const bool divb_cleaning);
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
index a6b2f0aa4..70800f732 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolver.cpp
@@ -30,7 +30,7 @@ SpectralSolver::SpectralSolver(
                 const bool pml, const bool periodic_single_box,
                 const bool update_with_rho,
                 const bool fft_do_time_averaging,
-                const bool J_linear_in_time,
+                const bool do_multi_J,
                 const bool dive_cleaning,
                 const bool divb_cleaning)
 {
@@ -42,7 +42,7 @@ SpectralSolver::SpectralSolver(
     const SpectralKSpace k_space= SpectralKSpace(realspace_ba, dm, dx);
 
     m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging,
-                                          J_linear_in_time, dive_cleaning, divb_cleaning, pml);
+                                          do_multi_J, dive_cleaning, divb_cleaning, pml);
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
@@ -63,7 +63,7 @@ SpectralSolver::SpectralSolver(
         else {
             algorithm = std::make_unique<PsatdAlgorithm>(
                 k_space, dm, m_spectral_index, norder_x, norder_y, norder_z, nodal, fill_guards,
-                v_galilean, dt, update_with_rho, fft_do_time_averaging, J_linear_in_time,
+                v_galilean, dt, update_with_rho, fft_do_time_averaging, do_multi_J,
                 dive_cleaning, divb_cleaning);
         }
     }
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
index b97487401..827d226c2 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.H
@@ -36,7 +36,7 @@ class SpectralSolverRZ
                           amrex::RealVect const dx, amrex::Real const dt,
                           bool const update_with_rho,
                           const bool fft_do_time_averaging,
-                          const bool J_linear_in_time,
+                          const bool do_multi_J,
                           const bool dive_cleaning,
                           const bool divb_cleaning);
 
diff --git a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
index b47dfa4ad..e187b9d27 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralSolverRZ.cpp
@@ -34,7 +34,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
                                     amrex::RealVect const dx, amrex::Real const dt,
                                     bool const update_with_rho,
                                     const bool fft_do_time_averaging,
-                                    const bool J_linear_in_time,
+                                    const bool do_multi_J,
                                     const bool dive_cleaning,
                                     const bool divb_cleaning)
     : k_space(realspace_ba, dm, dx)
@@ -47,7 +47,7 @@ SpectralSolverRZ::SpectralSolverRZ (const int lev,
 
     const bool pml = false;
     m_spectral_index = SpectralFieldIndex(update_with_rho, fft_do_time_averaging,
-                                          J_linear_in_time, dive_cleaning, divb_cleaning, pml);
+                                          do_multi_J, dive_cleaning, divb_cleaning, pml);
 
     // - Select the algorithm depending on the input parameters
     //   Initialize the corresponding coefficients over k space
@@ -56,7 +56,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, J_linear_in_time, dive_cleaning, divb_cleaning);
+            update_with_rho, fft_do_time_averaging, do_multi_J, dive_cleaning, divb_cleaning);
     } else {
         // Otherwise: use the Galilean algorithm
         algorithm = std::make_unique<GalileanPsatdAlgorithmRZ>(