div(E,B) Cleaning Options for PSATD (#2403)

* Implement div(E)/div(B) Cleaning with Standard PSATD

* Cleaning

* Update Benchmark

* Add Nodal Synchronization of F,G

* OneStep_multiJ: Nodal Syncs, Damp PML

* OneStep_multiJ: Push PSATD Fields in PML

* div Cleaning Defaults (Domain v. PML)

* Include Fix of #2429 until Merged

* Reset Benchmark of Langmuir_multi_psatd_div_cleaning

* Multi-J: Remove PML Support

* Include Fix of #2474 Until Merged

* Exchange All Guard Cells for F,G

* Fix Defaults

* Update Test, Reset Benchmark

* Fix Defaults

* Cleaning

* Default update_with_rho=1 if do_dive_cleaning=1

* Update CI Test pml_psatd_dive_divb_cleaning

* Replace Warning with Abort

2 files changed, 33 insertions, 33 deletions

diff --git a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
index 395f164cc..1f3c67629 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralAlgorithms/PsatdAlgorithm.cpp
@@ -103,14 +103,14 @@ PsatdAlgorithm::PsatdAlgorithm(
         InitializeSpectralCoefficientsAvgLin(spectral_kspace, dm, dt);
     }
 
-    if (dive_cleaning && !J_linear_in_time)
+    if (dive_cleaning && m_is_galilean)
     {
-        amrex::Abort("PSATD: warpx.do_dive_cleaning = 1 implemented only with psatd.J_linear_in_time = 1");
+        amrex::Abort("warpx.do_dive_cleaning = 1 not implemented for Galilean PSATD algorithms");
     }
 
-    if (divb_cleaning && !J_linear_in_time)
+    if (divb_cleaning && m_is_galilean)
     {
-        amrex::Abort("PSATD: warpx.do_divb_cleaning = 1 implemented only with psatd.J_linear_in_time = 1");
+        amrex::Abort("warpx.do_divb_cleaning = 1 not implemented for Galilean PSATD algorithms");
     }
     if (time_averaging && !update_with_rho)
     {
@@ -296,6 +296,30 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
                                    - I * T2 * S_ck * (kx * Ey_old - ky * Ex_old)
                                    + I * X1 * (kx * Jy - ky * Jx);
 
+            if (dive_cleaning)
+            {
+                const Complex k_dot_J  = kx * Jx + ky * Jy + kz * Jz;
+                const Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old;
+
+                fields(i,j,k,Idx.Ex) += I * c2 * S_ck * F_old * kx;
+                fields(i,j,k,Idx.Ey) += I * c2 * S_ck * F_old * ky;
+                fields(i,j,k,Idx.Ez) += I * c2 * S_ck * F_old * kz;
+
+                fields(i,j,k,Idx.F) = C * F_old + S_ck * (I * k_dot_E - rho_old * inv_ep0)
+                    - X1 * ((rho_new - rho_old) / dt + I * k_dot_J);
+            }
+
+            if (divb_cleaning)
+            {
+                const Complex k_dot_B = kx * Bx_old + ky * By_old + kz * Bz_old;
+
+                fields(i,j,k,Idx.Bx) += I * S_ck * G_old * kx;
+                fields(i,j,k,Idx.By) += I * S_ck * G_old * ky;
+                fields(i,j,k,Idx.Bz) += I * S_ck * G_old * kz;
+
+                fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B;
+            }
+
             if (J_linear_in_time)
             {
                 const Complex Jx_new = fields(i,j,k,Idx.Jx_new);
@@ -312,27 +336,9 @@ PsatdAlgorithm::pushSpectralFields (SpectralFieldData& f) const
 
                 if (dive_cleaning)
                 {
-                    const Complex k_dot_J  = kx * Jx + ky * Jy + kz * Jz;
                     const Complex k_dot_dJ = kx * (Jx_new - Jx) + ky * (Jy_new - Jy) + kz * (Jz_new - Jz);
-                    const Complex k_dot_E = kx * Ex_old + ky * Ey_old + kz * Ez_old;
-
-                    fields(i,j,k,Idx.Ex) += I * c2 * S_ck * F_old * kx;
-                    fields(i,j,k,Idx.Ey) += I * c2 * S_ck * F_old * ky;
-                    fields(i,j,k,Idx.Ez) += I * c2 * S_ck * F_old * kz;
-
-                    fields(i,j,k,Idx.F) = C * F_old + S_ck * (I * k_dot_E - rho_old * inv_ep0)
-                        - X1 * ((rho_new - rho_old) / dt + I * k_dot_J) - I * X2/c2 * k_dot_dJ;
-                }
-
-                if (divb_cleaning)
-                {
-                    const Complex k_dot_B = kx * Bx_old + ky * By_old + kz * Bz_old;
-
-                    fields(i,j,k,Idx.Bx) += I * S_ck * G_old * kx;
-                    fields(i,j,k,Idx.By) += I * S_ck * G_old * ky;
-                    fields(i,j,k,Idx.Bz) += I * S_ck * G_old * kz;
 
-                    fields(i,j,k,Idx.G) = C * G_old + I * c2 * S_ck * k_dot_B;
+                    fields(i,j,k,Idx.F) += -I * X2/c2 * k_dot_dJ;
                 }
 
                 if (time_averaging)
diff --git a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
index ebb720854..f16fc3aa2 100644
--- a/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
+++ b/Source/FieldSolver/SpectralSolver/SpectralFieldData.cpp
@@ -62,19 +62,13 @@ SpectralFieldIndex::SpectralFieldIndex (const bool update_with_rho,
             Bx_avg = c++; By_avg = c++; Bz_avg = c++;
         }
 
+        if (dive_cleaning) F = c++;
+
+        if (divb_cleaning) G = c++;
+
         if (J_linear_in_time)
         {
             Jx_new = c++; Jy_new = c++; Jz_new = c++;
-
-            if (dive_cleaning)
-            {
-                F = c++;
-            }
-
-            if (divb_cleaning)
-            {
-                G = c++;
-            }
         }
     }
     else // PML