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

4 files changed, 32 insertions, 31 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;