Separate Class for Multi-J PSATD Algo (#2748)

* Separate Class for Multi-J PSATD Algo

* Cleaning

* X1,...,X6 Real, not Complex

* Cleaning

* Improve Comments, Rename Jx as Jx_old (etc.)

1 files changed, 8 insertions, 168 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 2b1e09c2d..0515aa598 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -39,7 +39,6 @@ PsatdAlgorithm::PsatdAlgorithm(
     const amrex::Real dt,
     const bool update_with_rho,
     const bool time_averaging,
-    const bool do_multi_J,
     const bool dive_cleaning,
     const bool divb_cleaning)
     // Initializer list
@@ -59,7 +58,6 @@ PsatdAlgorithm::PsatdAlgorithm(
     m_dt(dt),
     m_update_with_rho(update_with_rho),
     m_time_averaging(time_averaging),
-    m_do_multi_J(do_multi_J),
     m_dive_cleaning(dive_cleaning),
     m_divb_cleaning(divb_cleaning)
 {
@@ -84,7 +82,7 @@ PsatdAlgorithm::PsatdAlgorithm(
     InitializeSpectralCoefficients(spectral_kspace, dm, dt);
 
     // Allocate these coefficients only with time averaging
-    if (time_averaging && !do_multi_J)
+    if (time_averaging)
     {
         Psi1_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
         Psi2_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
@@ -94,14 +92,6 @@ PsatdAlgorithm::PsatdAlgorithm(
         Y4_coef = SpectralComplexCoefficients(ba, dm, 1, 0);
         InitializeSpectralCoefficientsAveraging(spectral_kspace, dm, dt);
     }
-    // Allocate these coefficients only with time averaging
-    // 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);
-        InitializeSpectralCoefficientsAvgLin(spectral_kspace, dm, dt);
-    }
 
     if (dive_cleaning && m_is_galilean)
     {
@@ -121,12 +111,11 @@ PsatdAlgorithm::PsatdAlgorithm(
 void
 PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
 {
-    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 dive_cleaning    = m_dive_cleaning;
-    const bool divb_cleaning    = m_divb_cleaning;
-    const bool is_galilean      = m_is_galilean;
+    const bool update_with_rho = m_update_with_rho;
+    const bool time_averaging  = m_time_averaging;
+    const bool dive_cleaning   = m_dive_cleaning;
+    const bool divb_cleaning   = m_divb_cleaning;
+    const bool is_galilean     = m_is_galilean;
 
     const amrex::Real dt = m_dt;
 
@@ -163,7 +152,7 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
         amrex::Array4<const Complex> Y3_arr;
         amrex::Array4<const Complex> Y4_arr;
 
-        if (time_averaging && !do_multi_J)
+        if (time_averaging)
         {
             Psi1_arr = Psi1_coef[mfi].array();
             Psi2_arr = Psi2_coef[mfi].array();
@@ -173,14 +162,6 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
             Y4_arr = Y4_coef[mfi].array();
         }
 
-        Array4<const Complex> X5_arr;
-        Array4<const Complex> X6_arr;
-        if (time_averaging && do_multi_J)
-        {
-            X5_arr = X5_coef[mfi].array();
-            X6_arr = X6_coef[mfi].array();
-        }
-
         // Extract pointers for the k vectors
         const amrex::Real* modified_kx_arr = modified_kx_vec[mfi].dataPtr();
 #if defined(WARPX_DIM_3D)
@@ -229,7 +210,6 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
 #endif
             // Physical constants and imaginary unit
             constexpr Real c2 = PhysConst::c * PhysConst::c;
-            constexpr Real ep0 = PhysConst::ep0;
             constexpr Real inv_ep0 = 1._rt / PhysConst::ep0;
             constexpr Complex I = Complex{0._rt, 1._rt};
 
@@ -320,78 +300,8 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
                 fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B;
             }
 
-            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);
-                const Complex Jz_new = fields(i,j,k,Idx.Jz_new);
-
-                fields(i,j,k,Idx.Ex) += -X1 * (Jx_new - Jx) / dt;
-                fields(i,j,k,Idx.Ey) += -X1 * (Jy_new - Jy) / dt;
-                fields(i,j,k,Idx.Ez) += -X1 * (Jz_new - Jz) / dt;
-
-                fields(i,j,k,Idx.Bx) += I * X2/c2 * (ky * (Jz_new - Jz) - kz * (Jy_new - Jy));
-                fields(i,j,k,Idx.By) += I * X2/c2 * (kz * (Jx_new - Jx) - kx * (Jz_new - Jz));
-                fields(i,j,k,Idx.Bz) += I * X2/c2 * (kx * (Jy_new - Jy) - ky * (Jx_new - Jx));
-
-                if (dive_cleaning)
-                {
-                    const Complex k_dot_dJ = kx * (Jx_new - Jx) + ky * (Jy_new - Jy) + kz * (Jz_new - Jz);
-
-                    fields(i,j,k,Idx.F) += -I * X2/c2 * k_dot_dJ;
-                }
-
-                if (time_averaging)
-                {
-                    const Complex X5 = X5_arr(i,j,k);
-                    const Complex X6 = X6_arr(i,j,k);
-
-                    // TODO: Here the code is *accumulating* the average,
-                    // because it is meant to be used with sub-cycling
-                    // maybe this should be made more generic
-
-                    fields(i,j,k,Idx.Ex_avg) += S_ck * Ex_old
-                        + I * c2 * ep0 * X1 * (ky * Bz_old - kz * By_old)
-                        + I * X5 * rho_old * kx + I * X6 * rho_new * kx + X3/c2 * Jx - X2/c2 * Jx_new;
-
-                    fields(i,j,k,Idx.Ey_avg) += S_ck * Ey_old
-                        + I * c2 * ep0 * X1 * (kz * Bx_old - kx * Bz_old)
-                        + I * X5 * rho_old * ky + I * X6 * rho_new * ky + X3/c2 * Jy - X2/c2 * Jy_new;
-
-                    fields(i,j,k,Idx.Ez_avg) += S_ck * Ez_old
-                        + I * c2 * ep0 * X1 * (kx * By_old - ky * Bx_old)
-                        + I * X5 * rho_old * kz + I * X6 * rho_new * kz + X3/c2 * Jz - X2/c2 * Jz_new;
-
-                    fields(i,j,k,Idx.Bx_avg) += S_ck * Bx_old
-                        - I * ep0 * X1 * (ky * Ez_old - kz * Ey_old)
-                        - I * X5/c2 * (ky * Jz - kz * Jy) - I * X6/c2 * (ky * Jz_new - kz * Jy_new);
-
-                    fields(i,j,k,Idx.By_avg) += S_ck * By_old
-                        - I * ep0 * X1 * (kz * Ex_old - kx * Ez_old)
-                        - I * X5/c2 * (kz * Jx - kx * Jz) - I * X6/c2 * (kz * Jx_new - kx * Jz_new);
-
-                    fields(i,j,k,Idx.Bz_avg) += S_ck * Bz_old
-                        - I * ep0 * X1 * (kx * Ey_old - ky * Ex_old)
-                        - I * X5/c2 * (kx * Jy - ky * Jx) - I * X6/c2 * (kx * Jy_new - ky * Jx_new);
-
-                    if (dive_cleaning)
-                    {
-                        fields(i,j,k,Idx.Ex_avg) += I * c2 * ep0 * X1 * F_old * kx;
-                        fields(i,j,k,Idx.Ey_avg) += I * c2 * ep0 * X1 * F_old * ky;
-                        fields(i,j,k,Idx.Ez_avg) += I * c2 * ep0 * X1 * F_old * kz;
-                    }
-
-                    if (divb_cleaning)
-                    {
-                        fields(i,j,k,Idx.Bx_avg) += I * ep0 * X1 * G_old * kx;
-                        fields(i,j,k,Idx.By_avg) += I * ep0 * X1 * G_old * ky;
-                        fields(i,j,k,Idx.Bz_avg) += I * ep0 * X1 * G_old * kz;
-                    }
-                }
-            }
-
             // Additional update equations for averaged Galilean algorithm
-            if (time_averaging && !do_multi_J)
+            if (time_averaging)
             {
                 // These coefficients are initialized in the function InitializeSpectralCoefficients below
                 const Complex Psi1 = Psi1_arr(i,j,k);
@@ -822,76 +732,6 @@ void PsatdAlgorithm::InitializeSpectralCoefficientsAveraging (
     }
 }
 
-void PsatdAlgorithm::InitializeSpectralCoefficientsAvgLin (
-    const SpectralKSpace& spectral_kspace,
-    const amrex::DistributionMapping& dm,
-    const amrex::Real dt)
-{
-    const BoxArray& ba = spectral_kspace.spectralspace_ba;
-
-    // Loop over boxes and allocate the corresponding coefficients for each box
-    for (MFIter mfi(ba, dm); mfi.isValid(); ++mfi)
-    {
-        const Box& bx = ba[mfi];
-
-        // Extract pointers for the k vectors
-        const Real* kx_s = modified_kx_vec[mfi].dataPtr();
-#if defined(WARPX_DIM_3D)
-        const Real* ky_s = modified_ky_vec[mfi].dataPtr();
-#endif
-        const Real* kz_s = modified_kz_vec[mfi].dataPtr();
-
-        Array4<Real const> C = C_coef[mfi].array();
-        Array4<Real const> S_ck = S_ck_coef[mfi].array();
-
-        Array4<Complex> X5 = X5_coef[mfi].array();
-        Array4<Complex> X6 = X6_coef[mfi].array();
-
-        // Loop over indices within one box
-        ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k) noexcept
-        {
-            // Calculate norm of k vector
-            const Real knorm_s = std::sqrt(
-                std::pow(kx_s[i], 2) +
-#if defined(WARPX_DIM_3D)
-                std::pow(ky_s[j], 2) + std::pow(kz_s[k], 2));
-#else
-                std::pow(kz_s[j], 2));
-#endif
-            // Physical constants and imaginary unit
-            constexpr Real c = PhysConst::c;
-            constexpr Real c2 = c*c;
-            constexpr Real ep0 = PhysConst::ep0;
-
-            // Auxiliary coefficients
-            const Real dt3 = dt * dt * dt;
-
-            const Real om_s  = c * knorm_s;
-            const Real om2_s = om_s * om_s;
-            const Real om4_s = om2_s * om2_s;
-
-            if (om_s != 0.)
-            {
-                X5(i,j,k) = c2 / ep0 * (S_ck(i,j,k) / om2_s - (1._rt - C(i,j,k)) / (om4_s * dt)
-                                        - 0.5_rt * dt / om2_s);
-            }
-            else
-            {
-                X5(i,j,k) = - c2 * dt3 / (8._rt * ep0);
-            }
-
-            if (om_s != 0.)
-            {
-                X6(i,j,k) = c2 / ep0 * ((1._rt - C(i,j,k)) / (om4_s * dt) - 0.5_rt * dt / om2_s);
-            }
-            else
-            {
-                X6(i,j,k) = - c2 * dt3 / (24._rt * ep0);
-            }
-        });
-    }
-}
-
 void
 PsatdAlgorithm::CurrentCorrection (
     const int lev,