Filter diagnostic rho (#1295)

* First cut of adding filtering of rho diagnostic - not yet working

* Fixed the handling of parallel transfers in RhoFunctor when filtering is turned on

* Fixed end of line space

* In RhoFunctor, only do kspace filter with RZ and PSATD

* PICMI analytic applied fields (#1306)

* PICMI Added applied fields

* PICMI fixed typo in Constants

* Change picmistandard version requirement to 0.0.9

* Increment required picmi version

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

* Updated checksums with filter of diagnostic rho

* Fix comment in rho diagnostic

* Updated LaserIonAcc2d benchmark after filtering rho

* Update galilean_rz_psatd benchmark with filtering of rho

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

8 files changed, 60 insertions, 18 deletions

diff --git a/Regression/Checksum/benchmarks_json/LaserAcceleration.json b/Regression/Checksum/benchmarks_json/LaserAcceleration.json
index 5af011038..7a5bc9475 100644
--- a/Regression/Checksum/benchmarks_json/LaserAcceleration.json
+++ b/Regression/Checksum/benchmarks_json/LaserAcceleration.json
@@ -20,6 +20,6 @@
     "jx": 555905302485608.94,
     "jy": 1596092839924503.8,
     "jz": 1045523884130288.4,
-    "rho": 2213786555.2727532
+    "rho": 2211743042.5513821
   }
-}
\ No newline at end of file
+}
diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json b/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json
index 71d1f37cf..02281dd21 100644
--- a/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json
+++ b/Regression/Checksum/benchmarks_json/LaserAccelerationMR.json
@@ -29,7 +29,7 @@
     "jx": 1131265582036586.2,
     "jy": 2.325136123562568e+18,
     "jz": 6769514298165733.0,
-    "rho": 153059939.2876033
+    "rho": 152898720.2638083
   },
   "lev=1": {
     "Bx": 53303255.15572796,
@@ -41,6 +41,6 @@
     "jx": 48050235573831.96,
     "jy": 9.274838877361125e+18,
     "jz": 1.92627519812193e+16,
-    "rho": 67241940.12364514
+    "rho": 67304284.98422323
   }
-}
\ No newline at end of file
+}
diff --git a/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json b/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json
index c4ebe719d..567e0d908 100644
--- a/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json
+++ b/Regression/Checksum/benchmarks_json/LaserAccelerationRZ.json
@@ -31,6 +31,6 @@
     "jx": 913724361200.7623,
     "jy": 2.0581594676890458e+17,
     "jz": 1772971850765479.5,
-    "rho": 39161887.06938790
+    "rho": 38578930.57390296
   }
 }
diff --git a/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json b/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json
index 7a2aca2b6..9e7693862 100644
--- a/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json
+++ b/Regression/Checksum/benchmarks_json/LaserIonAcc2d.json
@@ -29,8 +29,8 @@
     "jx": 1.1795616964119449e+20,
     "jy": 0.0,
     "jz": 7.413959422336308e+19,
-    "rho": 893318471772.7588,
+    "rho": 344505246673.0273,
     "rho_electrons": 74759628935311.0,
     "rho_hydrogen": 74759628935311.02
   }
-}
\ No newline at end of file
+}
diff --git a/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json b/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json
index e49a6a66b..3b6b66679 100644
--- a/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json
+++ b/Regression/Checksum/benchmarks_json/galilean_rz_psatd.json
@@ -29,6 +29,6 @@
     "divE": 1349004.6304044977,
     "jx": 163.91848307748262,
     "jz": 50271.11378286561,
-    "rho": 0.0001686919398337135
+    "rho": 0.0000120119242959871
   }
-}
\ No newline at end of file
+}
diff --git a/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp
index c8ccb2ded..e48703e91 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp
+++ b/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp
@@ -2,6 +2,10 @@
 #include "RhoFunctor.H"
 #include "Utils/CoarsenIO.H"
 
+#ifdef WARPX_DIM_RZ
+#include "FieldSolver/SpectralSolver/SpectralFieldData.H"
+#endif
+
 #include <AMReX.H>
 
 RhoFunctor::RhoFunctor (const int lev,
@@ -21,17 +25,36 @@ RhoFunctor::operator() ( amrex::MultiFab& mf_dst, const int dcomp, const int /*i
     auto& warpx = WarpX::GetInstance();
     std::unique_ptr<amrex::MultiFab> rho;
 
+    // Deposit charge density
+    // Call this with local=true since the parallel transfers will be handled
+    // by ApplyFilterandSumBoundaryRho
+
     // Dump total rho
     if (m_species_index == -1) {
         auto& mypc = warpx.GetPartContainer();
-        rho = mypc.GetChargeDensity(m_lev);
+        rho = mypc.GetChargeDensity(m_lev, true);
     }
     // Dump rho per species
     else {
         auto& mypc = warpx.GetPartContainer().GetParticleContainer(m_species_index);
-        rho = mypc.GetChargeDensity(m_lev);
+        rho = mypc.GetChargeDensity(m_lev, true);
     }
 
+    // Handle the parallel transfers of guard cells and
+    // apply the filtering if requested.
+    warpx.ApplyFilterandSumBoundaryRho(m_lev, m_lev, *rho, 0, rho->nComp());
+
+#if (defined WARPX_DIM_RZ) && (defined WARPX_USE_PSATD)
+    using Idx = SpectralAvgFieldIndex;
+    if (WarpX::use_kspace_filter) {
+        auto & solver = warpx.get_spectral_solver_fp(m_lev);
+        solver.ForwardTransform(*rho, Idx::rho_new);
+        solver.ApplyFilter(Idx::rho_new);
+        solver.BackwardTransform(*rho, Idx::rho_new);
+    }
+#endif
+
+
 #ifdef WARPX_DIM_RZ
     if (m_convertRZmodes2cartesian) {
         // In cylindrical geometry, sum real part of all modes of rho in
diff --git a/Source/Parallelization/WarpXComm.cpp b/Source/Parallelization/WarpXComm.cpp
index 667be0696..daee24054 100644
--- a/Source/Parallelization/WarpXComm.cpp
+++ b/Source/Parallelization/WarpXComm.cpp
@@ -873,17 +873,23 @@ void
 WarpX::ApplyFilterandSumBoundaryRho (int lev, PatchType patch_type, int icomp, int ncomp)
 {
     const int glev = (patch_type == PatchType::fine) ? lev : lev-1;
-    const auto& period = Geom(glev).periodicity();
     auto& r = (patch_type == PatchType::fine) ? rho_fp[lev] : rho_cp[lev];
     if (r == nullptr) return;
+    ApplyFilterandSumBoundaryRho(lev, glev, *r, icomp, ncomp);
+}
+
+void
+WarpX::ApplyFilterandSumBoundaryRho (int lev, int glev, amrex::MultiFab& rho, int icomp, int ncomp)
+{
+    const auto& period = Geom(glev).periodicity();
     if (use_filter) {
-        IntVect ng = r->nGrowVect();
+        IntVect ng = rho.nGrowVect();
         ng += bilinear_filter.stencil_length_each_dir-1;
-        MultiFab rf(r->boxArray(), r->DistributionMap(), ncomp, ng);
-        bilinear_filter.ApplyStencil(rf, *r, icomp, 0, ncomp);
-        WarpXSumGuardCells(*r, rf, period, icomp, ncomp );
+        MultiFab rf(rho.boxArray(), rho.DistributionMap(), ncomp, ng);
+        bilinear_filter.ApplyStencil(rf, rho, icomp, 0, ncomp);
+        WarpXSumGuardCells(rho, rf, period, icomp, ncomp );
     } else {
-        WarpXSumGuardCells(*r, period, icomp, ncomp);
+        WarpXSumGuardCells(rho, period, icomp, ncomp);
     }
 }
 
diff --git a/Source/WarpX.H b/Source/WarpX.H
index 68dd45949..928e40b42 100644
--- a/Source/WarpX.H
+++ b/Source/WarpX.H
@@ -543,6 +543,8 @@ public:
      */
     void ComputeCostsHeuristic (amrex::Vector<std::unique_ptr<amrex::LayoutData<amrex::Real> > >& costs);
 
+    void ApplyFilterandSumBoundaryRho (int lev, int glev, amrex::MultiFab& rho, int icomp, int ncomp);
+
 protected:
 
     /**
@@ -877,7 +879,18 @@ private:
         amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_fp;
         amrex::Vector<std::unique_ptr<SpectralSolver>> spectral_solver_cp;
 #   endif
+
+public:
+
+#   ifdef WARPX_DIM_RZ
+        SpectralSolverRZ&
+#   else
+        SpectralSolver&
+#   endif
+            get_spectral_solver_fp (int lev) {return *spectral_solver_fp[lev];}
 #endif
+
+private:
     amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> m_fdtd_solver_fp;
     amrex::Vector<std::unique_ptr<FiniteDifferenceSolver>> m_fdtd_solver_cp;
 };