Include `J` in diagnostic output when an electromagnetic solver is not used (#4116)

* include current density in diagnostic output even if an electromagnetic solver is not used

* do not redeposit current for the magnetostatic solver

* update CI benchmarks for tests that previously did not register current density

* fix remaining failing CI test

* Apply suggestions from code review

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

* do not output current density in collision tests

* update `JFunctor` constructor doc-string

* code cleanup - reduce code duplication

* specify `Ex Ey Ez Bx By Bz` fields to plot for collision CI tests

* specify `Er Et Ez Br Bt Bz` as output for rz collision test

* rename `InterpolateToDst` to `InterpolateMFForDiag`

* only deposit current density once per diagnostic output (if not already deposited)

* write deposited current for all directions into `current_fp` during diagnostic step deposition

* use `amrex::make_alias` to directly deposit current density into `warpx.current_fp` during diagnostic step deposition

* add class variable `solver_deposits_current` to `FullDiagnostics` so that there is only one place where the determination is made whether current should be deposited during diagnostic output

* correct logic to determine `m_solver_deposits_current`

---------

Co-authored-by: Edoardo Zoni <59625522+EZoni@users.noreply.github.com>

19 files changed, 239 insertions, 121 deletions

diff --git a/Examples/Tests/collision/inputs_2d b/Examples/Tests/collision/inputs_2d
index 5f81205f3..69d23dfb0 100644
--- a/Examples/Tests/collision/inputs_2d
+++ b/Examples/Tests/collision/inputs_2d
@@ -72,6 +72,7 @@ collision3.CoulombLog = 15.9
 diagnostics.diags_names = diag1 diag_parser_filter diag_uniform_filter diag_random_filter
 diag1.intervals = 10
 diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz
 
 ## diag_parser_filter is a diag used to test the particle filter function.
 diag_parser_filter.intervals = 150:150:
diff --git a/Examples/Tests/collision/inputs_3d b/Examples/Tests/collision/inputs_3d
index 5058dcb8f..3cc06061b 100644
--- a/Examples/Tests/collision/inputs_3d
+++ b/Examples/Tests/collision/inputs_3d
@@ -75,6 +75,7 @@ collision3.ndt = 10
 diagnostics.diags_names = diag1 diag_parser_filter diag_uniform_filter diag_random_filter
 diag1.intervals = 10
 diag1.diag_type = Full
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz
 
 ## diag_parser_filter is a diag used to test the particle filter function.
 diag_parser_filter.intervals = 150:150:
diff --git a/Examples/Tests/collision/inputs_rz b/Examples/Tests/collision/inputs_rz
index de6c8517c..ec1a234b1 100644
--- a/Examples/Tests/collision/inputs_rz
+++ b/Examples/Tests/collision/inputs_rz
@@ -56,3 +56,4 @@ collision.CoulombLog = 15.9
 diagnostics.diags_names = diag1
 diag1.intervals = 10
 diag1.diag_type = Full
+diag1.fields_to_plot = Er Et Ez Br Bt Bz
diff --git a/Examples/Tests/particle_boundary_process/PICMI_inputs_reflection.py b/Examples/Tests/particle_boundary_process/PICMI_inputs_reflection.py
index d8d508313..18d91ad45 100755
--- a/Examples/Tests/particle_boundary_process/PICMI_inputs_reflection.py
+++ b/Examples/Tests/particle_boundary_process/PICMI_inputs_reflection.py
@@ -74,13 +74,14 @@ electrons = picmi.Species(
 # diagnostics
 ##########################
 
-field_diag = picmi.ParticleDiagnostic(
-    species=electrons,
+field_diag = picmi.FieldDiagnostic(
+    grid=grid,
     name = 'diag1',
-    data_list=['previous_positions'],
+    data_list=['E'],
     period = 10,
     write_dir = '.',
-    warpx_file_prefix = 'Python_particle_reflection_plt'
+    warpx_file_prefix = 'Python_particle_reflection_plt',
+    warpx_write_species=False
 )
 
 ##########################
diff --git a/Regression/Checksum/benchmarks_json/Python_prev_positions.json b/Regression/Checksum/benchmarks_json/Python_prev_positions.json
index 44f959cb7..e6e39d227 100644
--- a/Regression/Checksum/benchmarks_json/Python_prev_positions.json
+++ b/Regression/Checksum/benchmarks_json/Python_prev_positions.json
@@ -1,23 +1,23 @@
 {
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 3710588.849989976,
+    "Ey": 0.0,
+    "Ez": 646727.8074440088,
+    "jx": 4745.078379617619,
+    "jy": 368.4331779923921,
+    "jz": 632.1508106460103
+  },
   "electrons": {
     "particle_momentum_x": 8.78764082600202e-23,
     "particle_momentum_y": 3.78128494181519e-24,
-    "particle_momentum_z": 5.396907550474659e-24,
+    "particle_momentum_z": 5.3969075504746586e-24,
     "particle_position_x": 0.6227948811898449,
     "particle_position_y": 0.5553787034074176,
     "particle_prev_x": 0.6152547961491118,
     "particle_prev_z": 0.5559363018627721,
     "particle_weight": 8569335937.5
-  },
-  "lev=0": {
-    "Bx": 0.0,
-    "By": 0.0,
-    "Bz": 0.0,
-    "Ex": 3710588.849989976,
-    "Ey": 0.0,
-    "Ez": 646727.8074440088,
-    "jx": 0.0,
-    "jy": 0.0,
-    "jz": 0.0
   }
-}
\ No newline at end of file
+}
diff --git a/Regression/Checksum/benchmarks_json/collisionRZ.json b/Regression/Checksum/benchmarks_json/collisionRZ.json
index d0dee71a4..7b6edb345 100644
--- a/Regression/Checksum/benchmarks_json/collisionRZ.json
+++ b/Regression/Checksum/benchmarks_json/collisionRZ.json
@@ -5,9 +5,6 @@
     "Bz": 0.0,
     "Er": 0.0,
     "Et": 0.0,
-    "Ez": 0.0,
-    "jr": 0.0,
-    "jt": 0.0,
-    "jz": 0.0
+    "Ez": 0.0
   }
-}
\ No newline at end of file
+}
diff --git a/Regression/Checksum/benchmarks_json/collisionXYZ.json b/Regression/Checksum/benchmarks_json/collisionXYZ.json
index 148bb80c1..927848745 100644
--- a/Regression/Checksum/benchmarks_json/collisionXYZ.json
+++ b/Regression/Checksum/benchmarks_json/collisionXYZ.json
@@ -5,9 +5,6 @@
     "Bz": 0.0,
     "Ex": 0.0,
     "Ey": 0.0,
-    "Ez": 0.0,
-    "jx": 0.0,
-    "jy": 0.0,
-    "jz": 0.0
+    "Ez": 0.0
   }
 }
diff --git a/Regression/Checksum/benchmarks_json/collisionXZ.json b/Regression/Checksum/benchmarks_json/collisionXZ.json
index d24938471..927848745 100644
--- a/Regression/Checksum/benchmarks_json/collisionXZ.json
+++ b/Regression/Checksum/benchmarks_json/collisionXZ.json
@@ -5,9 +5,6 @@
     "Bz": 0.0,
     "Ex": 0.0,
     "Ey": 0.0,
-    "Ez": 0.0,
-    "jx": 0.0,
-    "jy": 0.0,
-    "jz": 0.0
+    "Ez": 0.0
   }
-}
\ No newline at end of file
+}
diff --git a/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles.json b/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles.json
index 59894c50e..80a2dc879 100644
--- a/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles.json
+++ b/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles.json
@@ -6,9 +6,9 @@
     "Ex": 9.882421146615367e-06,
     "Ey": 1.0440261046714249e-05,
     "Ez": 1.003739697324731e-05,
-    "jx": 0.0,
-    "jy": 0.0,
-    "jz": 0.0
+    "jx": 2.9148662809570633e-10,
+    "jy": 8.46582291468749e-19,
+    "jz": 3.823492756863969e-08
   },
   "electron": {
     "particle_momentum_x": 2.0819392991319055e-25,
diff --git a/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles_moving.json b/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles_moving.json
index f00caa736..813eead73 100644
--- a/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles_moving.json
+++ b/Regression/Checksum/benchmarks_json/hard_edged_quadrupoles_moving.json
@@ -1,22 +1,22 @@
 {
-  "electron": {
-    "particle_momentum_x": 2.0819392998019253e-25,
-    "particle_momentum_y": 2.3316091155621773e-33,
-    "particle_momentum_z": 2.730924530757351e-23,
-    "particle_position_x": 0.034923287741880055,
-    "particle_position_y": 9.300626840674331e-11,
-    "particle_position_z": 1.491521766472402,
-    "particle_weight": 1.0
-  },
   "lev=0": {
     "Bx": 0.0,
     "By": 0.0,
     "Bz": 0.0,
-    "Ex": 6.028256519009052e-05,
-    "Ey": 6.384796595673982e-05,
-    "Ez": 7.88045921306518e-05,
-    "jx": 0.0,
-    "jy": 0.0,
-    "jz": 0.0
+    "Ex": 6.0282565190090465e-05,
+    "Ey": 6.38479659567398e-05,
+    "Ez": 7.880459213065183e-05,
+    "jx": 1.1659465170956616e-09,
+    "jy": 2.6115239381823616e-17,
+    "jz": 1.5293971084510288e-07
+  },
+  "electron": {
+    "particle_momentum_x": 2.0819392998019267e-25,
+    "particle_momentum_y": 2.3316091155621786e-33,
+    "particle_momentum_z": 2.730924530757353e-23,
+    "particle_position_x": 0.03492328774188008,
+    "particle_position_y": 9.300626840674328e-11,
+    "particle_position_z": 1.4915217664724023,
+    "particle_weight": 1.0
   }
 }
diff --git a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt
index 01f780a57..2ef6af16b 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt
+++ b/Source/Diagnostics/ComputeDiagFunctors/CMakeLists.txt
@@ -5,6 +5,7 @@ foreach(D IN LISTS WarpX_DIMS)
         CellCenterFunctor.cpp
         DivBFunctor.cpp
         DivEFunctor.cpp
+        JFunctor.cpp
         RhoFunctor.cpp
         PartPerCellFunctor.cpp
         PartPerGridFunctor.cpp
diff --git a/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp
index f43714a9c..fe5037131 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp
+++ b/Source/Diagnostics/ComputeDiagFunctors/CellCenterFunctor.cpp
@@ -1,11 +1,6 @@
 #include "CellCenterFunctor.H"
 
-#include "Utils/TextMsg.H"
-#ifdef WARPX_DIM_RZ
-#   include "WarpX.H"
-#endif
-
-#include <ablastr/coarsen/sample.H>
+#include "WarpX.H"
 
 #include <AMReX.H>
 #include <AMReX_IntVect.H>
@@ -21,29 +16,7 @@ CellCenterFunctor::CellCenterFunctor(amrex::MultiFab const * mf_src, int lev,
 void
 CellCenterFunctor::operator()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer*/) const
 {
-#ifdef WARPX_DIM_RZ
-    if (m_convertRZmodes2cartesian) {
-        // In cylindrical geometry, sum real part of all modes of m_mf_src in
-        // temporary multifab mf_dst_stag, and cell-center it to mf_dst.
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            nComp()==1,
-            "The RZ averaging over modes must write into 1 single component");
-        auto& warpx = WarpX::GetInstance();
-        amrex::MultiFab mf_dst_stag(m_mf_src->boxArray(), warpx.DistributionMap(m_lev), 1, m_mf_src->nGrowVect());
-        // Mode 0
-        amrex::MultiFab::Copy(mf_dst_stag, *m_mf_src, 0, 0, 1, m_mf_src->nGrowVect());
-        for (int ic=1 ; ic < m_mf_src->nComp() ; ic += 2) {
-            // All modes > 0
-            amrex::MultiFab::Add(mf_dst_stag, *m_mf_src, ic, 0, 1, m_mf_src->nGrowVect());
-        }
-        ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0,  m_crse_ratio);
-    } else {
-        ablastr::coarsen::sample::Coarsen( mf_dst, *m_mf_src, dcomp, 0, nComp(), 0, m_crse_ratio);
-    }
-#else
-    // In cartesian geometry, coarsen and interpolate from simulation MultiFab, m_mf_src,
-    // to output diagnostic MultiFab, mf_dst.
-    ablastr::coarsen::sample::Coarsen(mf_dst, *m_mf_src, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio);
-    amrex::ignore_unused(m_lev, m_convertRZmodes2cartesian);
-#endif
+    auto& warpx = WarpX::GetInstance();
+    InterpolateMFForDiag(mf_dst, *m_mf_src, dcomp, warpx.DistributionMap(m_lev),
+                         m_convertRZmodes2cartesian);
 }
diff --git a/Source/Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H
index 1f71732d4..53a9aff8f 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H
+++ b/Source/Diagnostics/ComputeDiagFunctors/ComputeDiagFunctor.H
@@ -2,6 +2,9 @@
 #define WARPX_COMPUTEDIAGFUNCTOR_H_
 
 #include "ComputeDiagFunctor_fwd.H"
+#include "Utils/TextMsg.H"
+
+#include <ablastr/coarsen/sample.H>
 
 #include <AMReX.H>
 #include <AMReX_MultiFab.H>
@@ -57,9 +60,41 @@ public:
                                           k_index_zlab, max_box_size, snapshot_full);
                                       }
     virtual void InitData() {}
+
+    void InterpolateMFForDiag (
+        amrex::MultiFab& mf_dst, const amrex::MultiFab& mf_src, int dcomp,
+        const amrex::DistributionMapping& dm, bool convertRZmodes2cartesian ) const
+    {
+#ifdef WARPX_DIM_RZ
+        if (convertRZmodes2cartesian) {
+            // In cylindrical geometry, sum real part of all modes of mf_src in
+            // temporary MultiFab mf_dst_stag, and cell-center it to mf_dst
+            WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+                nComp()==1,
+                "The RZ averaging over modes must write into one single component");
+            amrex::MultiFab mf_dst_stag( mf_src.boxArray(), dm, 1, mf_src.nGrowVect() );
+            // Mode 0
+            amrex::MultiFab::Copy( mf_dst_stag, mf_src, 0, 0, 1, mf_src.nGrowVect() );
+            for (int ic=1 ; ic < mf_src.nComp() ; ic += 2) {
+                // Real part of all modes > 0
+                amrex::MultiFab::Add( mf_dst_stag, mf_src, ic, 0, 1, mf_src.nGrowVect() );
+            }
+            ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio );
+        } else {
+            ablastr::coarsen::sample::Coarsen( mf_dst, mf_src, dcomp, 0, nComp(), 0, m_crse_ratio );
+        }
+#else
+        // In Cartesian geometry, coarsen and interpolate from temporary MultiFab mf_src
+        // to output diagnostic MultiFab mf_dst
+        ablastr::coarsen::sample::Coarsen(mf_dst, mf_src, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio );
+        amrex::ignore_unused(convertRZmodes2cartesian, dm);
+#endif
+    }
+
 private:
     /** Number of components of mf_dst that this functor updates. */
     int m_ncomp;
+
 protected:
     /** Coarsening ratio used to interpolate fields from simulation MultiFabs to output MultiFab. */
     amrex::IntVect m_crse_ratio;
diff --git a/Source/Diagnostics/ComputeDiagFunctors/JFunctor.H b/Source/Diagnostics/ComputeDiagFunctors/JFunctor.H
new file mode 100644
index 000000000..c3be88bdb
--- /dev/null
+++ b/Source/Diagnostics/ComputeDiagFunctors/JFunctor.H
@@ -0,0 +1,50 @@
+#ifndef WARPX_JFUNCTOR_H_
+#define WARPX_JFUNCTOR_H_
+
+#include "ComputeDiagFunctor.H"
+
+#include <AMReX_BaseFwd.H>
+
+/**
+ * \brief Functor to cell-center MF for current density and store result in mf_out.
+ */
+class
+JFunctor final : public ComputeDiagFunctor
+{
+public:
+    /** Constructor.
+     *
+     * \param[in] dir direction of vector field to operate on
+     * \param[in] lev level of multifab. Used for averaging in rz.
+     * \param[in] crse_ratio for interpolating field values from the simulation MultiFab, src_mf,
+                  to the output diagnostic MultiFab, mf_dst.
+     * \param[in] convertRZmodes2cartesian (in cylindrical) whether to
+     *            sum all modes in mf_src before cell-centering into dst multifab.
+     * \param[in] ncomp Number of component of mf_src to cell-center in dst multifab.
+     */
+    JFunctor (const int dir, const int lev,
+             const amrex::IntVect crse_ratio,
+             bool convertRZmodes2cartesian=true,
+             bool deposit_current=false, int ncomp=1);
+    /** \brief Cell-center m_mf_src and write the result in mf_dst.
+     *
+     * In cylindrical geometry, by default this functor average all components
+     * of a MultiFab and writes into one single component.
+     *
+     * \param[out] mf_dst output MultiFab where the result is written
+     * \param[in] dcomp first component of mf_dst in which cell-centered
+     *            data is stored
+     */
+    virtual void operator()(amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer=0*/) const override;
+private:
+    /** direction of the current density to save */
+    const int m_dir;
+    /** level on which mf_src is defined */
+    int m_lev;
+    /** (for cylindrical) whether to average all modes into 1 comp */
+    bool m_convertRZmodes2cartesian;
+    /** whether to deposit current density before saving */
+    bool m_deposit_current;
+};
+
+#endif // WARPX_JFUNCTOR_H_
diff --git a/Source/Diagnostics/ComputeDiagFunctors/JFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/JFunctor.cpp
new file mode 100644
index 000000000..4663a5e1a
--- /dev/null
+++ b/Source/Diagnostics/ComputeDiagFunctors/JFunctor.cpp
@@ -0,0 +1,59 @@
+/* This file is part of WarpX.
+ *
+ * Authors: Roelof Groenewald
+ * License: BSD-3-Clause-LBNL
+ */
+
+#include "JFunctor.H"
+
+#include "Particles/MultiParticleContainer.H"
+#include "WarpX.H"
+
+#include <AMReX.H>
+#include <AMReX_IntVect.H>
+#include <AMReX_MultiFab.H>
+
+JFunctor::JFunctor (const int dir, int lev,
+                   amrex::IntVect crse_ratio,
+                   bool convertRZmodes2cartesian,
+                   bool deposit_current, int ncomp)
+    : ComputeDiagFunctor(ncomp, crse_ratio), m_dir(dir), m_lev(lev),
+      m_convertRZmodes2cartesian(convertRZmodes2cartesian),
+      m_deposit_current(deposit_current)
+{ }
+
+void
+JFunctor::operator() (amrex::MultiFab& mf_dst, int dcomp, const int /*i_buffer*/) const
+{
+    auto& warpx = WarpX::GetInstance();
+    /** pointer to source multifab (can be multi-component) */
+    amrex::MultiFab* m_mf_src = warpx.get_pointer_current_fp(m_lev, m_dir);
+
+    // Deposit current if no solver or the electrostatic solver is being used
+    if (m_deposit_current)
+    {
+        // allocate temporary multifab to deposit current density into
+        amrex::Vector<std::array< std::unique_ptr<amrex::MultiFab>, 3 > > current_fp_temp;
+        current_fp_temp.resize(1);
+
+        auto& current_fp_x = warpx.getcurrent_fp(m_lev,0);
+        current_fp_temp[0][0] = std::make_unique<amrex::MultiFab>(
+            current_fp_x, amrex::make_alias, 0, current_fp_x.nComp()
+        );
+
+        auto& current_fp_y = warpx.getcurrent_fp(m_lev,1);
+        current_fp_temp[0][1] = std::make_unique<amrex::MultiFab>(
+            current_fp_y, amrex::make_alias, 0, current_fp_y.nComp()
+        );
+        auto& current_fp_z = warpx.getcurrent_fp(m_lev,2);
+        current_fp_temp[0][2] = std::make_unique<amrex::MultiFab>(
+            current_fp_z, amrex::make_alias, 0, current_fp_z.nComp()
+        );
+
+        auto& mypc = warpx.GetPartContainer();
+        mypc.DepositCurrent(current_fp_temp, warpx.getdt(m_lev), 0.0);
+    }
+
+    InterpolateMFForDiag(mf_dst, *m_mf_src, dcomp, warpx.DistributionMap(m_lev),
+                         m_convertRZmodes2cartesian);
+}
diff --git a/Source/Diagnostics/ComputeDiagFunctors/Make.package b/Source/Diagnostics/ComputeDiagFunctors/Make.package
index 611ccb60d..9c329263b 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/Make.package
+++ b/Source/Diagnostics/ComputeDiagFunctors/Make.package
@@ -3,6 +3,7 @@ CEXE_sources += PartPerCellFunctor.cpp
 CEXE_sources += PartPerGridFunctor.cpp
 CEXE_sources += DivBFunctor.cpp
 CEXE_sources += DivEFunctor.cpp
+CEXE_sources += JFunctor.cpp
 CEXE_sources += RhoFunctor.cpp
 CEXE_sources += BackTransformFunctor.cpp
 CEXE_sources += BackTransformParticleFunctor.cpp
diff --git a/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp b/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp
index 6f4fcda4d..73feef4ed 100644
--- a/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp
+++ b/Source/Diagnostics/ComputeDiagFunctors/RhoFunctor.cpp
@@ -8,11 +8,8 @@
 #endif
 #include "Particles/MultiParticleContainer.H"
 #include "Particles/WarpXParticleContainer.H"
-#include "Utils/TextMsg.H"
 #include "WarpX.H"
 
-#include <ablastr/coarsen/sample.H>
-
 #include <AMReX.H>
 #include <AMReX_IntVect.H>
 #include <AMReX_MultiFab.H>
@@ -69,29 +66,6 @@ RhoFunctor::operator() ( amrex::MultiFab& mf_dst, const int dcomp, const int /*i
     }
 #endif
 
-
-#ifdef WARPX_DIM_RZ
-    if (m_convertRZmodes2cartesian) {
-        // In cylindrical geometry, sum real part of all modes of rho in
-        // temporary MultiFab mf_dst_stag, and cell-center it to mf_dst
-        WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
-            nComp()==1,
-            "The RZ averaging over modes must write into one single component");
-        amrex::MultiFab mf_dst_stag( rho->boxArray(), warpx.DistributionMap(m_lev), 1, rho->nGrowVect() );
-        // Mode 0
-        amrex::MultiFab::Copy( mf_dst_stag, *rho, 0, 0, 1, rho->nGrowVect() );
-        for (int ic=1 ; ic < rho->nComp() ; ic += 2) {
-            // Real part of all modes > 0
-            amrex::MultiFab::Add( mf_dst_stag, *rho, ic, 0, 1, rho->nGrowVect() );
-        }
-        ablastr::coarsen::sample::Coarsen( mf_dst, mf_dst_stag, dcomp, 0, nComp(), 0, m_crse_ratio );
-    } else {
-        ablastr::coarsen::sample::Coarsen( mf_dst, *rho, dcomp, 0, nComp(), 0, m_crse_ratio );
-    }
-#else
-    // In Cartesian geometry, coarsen and interpolate from temporary MultiFab rho
-    // to output diagnostic MultiFab mf_dst
-    ablastr::coarsen::sample::Coarsen(mf_dst, *rho, dcomp, 0, nComp(), mf_dst.nGrowVect(), m_crse_ratio );
-    amrex::ignore_unused(m_convertRZmodes2cartesian);
-#endif
+    InterpolateMFForDiag(mf_dst, *rho, dcomp, warpx.DistributionMap(m_lev),
+                         m_convertRZmodes2cartesian);
 }
diff --git a/Source/Diagnostics/FullDiagnostics.H b/Source/Diagnostics/FullDiagnostics.H
index a20b44397..4bb3e2789 100644
--- a/Source/Diagnostics/FullDiagnostics.H
+++ b/Source/Diagnostics/FullDiagnostics.H
@@ -22,6 +22,10 @@ private:
     bool m_plot_raw_fields_guards = false;
     /** Whether to dump the RZ modes */
     bool m_dump_rz_modes = false;
+    /** Whether the field solver deposits current density, if not it will be done
+     * before writing the diagnostic.
+     */
+    bool m_solver_deposits_current = true;
     /** Flush m_mf_output and particles to file for the i^th buffer */
     void Flush (int i_buffer) override;
     /** Flush raw data */
diff --git a/Source/Diagnostics/FullDiagnostics.cpp b/Source/Diagnostics/FullDiagnostics.cpp
index 13dff5432..7abe923b3 100644
--- a/Source/Diagnostics/FullDiagnostics.cpp
+++ b/Source/Diagnostics/FullDiagnostics.cpp
@@ -3,6 +3,7 @@
 #include "ComputeDiagFunctors/CellCenterFunctor.H"
 #include "ComputeDiagFunctors/DivBFunctor.H"
 #include "ComputeDiagFunctors/DivEFunctor.H"
+#include "ComputeDiagFunctors/JFunctor.H"
 #include "ComputeDiagFunctors/PartPerCellFunctor.H"
 #include "ComputeDiagFunctors/PartPerGridFunctor.H"
 #include "ComputeDiagFunctors/ParticleReductionFunctor.H"
@@ -45,6 +46,11 @@ FullDiagnostics::FullDiagnostics (int i, std::string name)
 {
     ReadParameters();
     BackwardCompatibility();
+
+    m_solver_deposits_current = !(
+        WarpX::electromagnetic_solver_id == ElectromagneticSolverAlgo::None &&
+        WarpX::electrostatic_solver_id != ElectrostaticSolverAlgo::LabFrameElectroMagnetostatic
+    );
 }
 
 void
@@ -192,6 +198,10 @@ FullDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev)
     m_all_field_functors[lev].clear();
     m_all_field_functors[lev].resize(m_varnames_fields.size());
 
+    // Boolean flag for whether the current density should be deposited before
+    // diagnostic output
+    bool deposit_current = !m_solver_deposits_current;
+
     // Fill vector of functors for all components except individual cylindrical modes.
     for (int comp=0, n=m_varnames_fields.size(); comp<n; comp++){
         if        ( m_varnames_fields[comp] == "Er" ){
@@ -231,20 +241,23 @@ FullDiagnostics::InitializeFieldFunctorsRZopenPMD (int lev)
                 AddRZModesToOutputNames(std::string("Bz"), ncomp);
             }
         } else if ( m_varnames_fields[comp] == "jr" ){
-            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 0), lev, m_crse_ratio,
-                                                        false, ncomp);
+            m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(0, lev, m_crse_ratio,
+                                                        false, deposit_current, ncomp);
+            deposit_current = false;
             if (update_varnames) {
                 AddRZModesToOutputNames(std::string("jr"), ncomp);
             }
         } else if ( m_varnames_fields[comp] == "jt" ){
-            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 1), lev, m_crse_ratio,
-                                                        false, ncomp);
+            m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(1, lev, m_crse_ratio,
+                                                        false, deposit_current, ncomp);
+            deposit_current = false;
             if (update_varnames) {
                 AddRZModesToOutputNames(std::string("jt"), ncomp);
             }
         } else if ( m_varnames_fields[comp] == "jz" ){
-            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 2), lev, m_crse_ratio,
-                                                        false, ncomp);
+            m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(2, lev, m_crse_ratio,
+                                                        false, deposit_current, ncomp);
+            deposit_current = false;
             if (update_varnames) {
                 AddRZModesToOutputNames(std::string("jz"), ncomp);
             }
@@ -353,6 +366,10 @@ FullDiagnostics::AddRZModesToDiags (int lev)
     // If rho is requested, all components will be written out
     const bool rho_requested = utils::algorithms::is_in( m_varnames, "rho" );
 
+    // Boolean flag for whether the current density should be deposited before
+    // diagnostic output
+    bool deposit_current = !m_solver_deposits_current;
+
     // First index of m_all_field_functors[lev] where RZ modes are stored
     int icomp = m_all_field_functors[0].size();
     const std::array<std::string, 3> coord {"r", "theta", "z"};
@@ -391,8 +408,8 @@ FullDiagnostics::AddRZModesToDiags (int lev)
     for (int dim=0; dim<3; dim++){
         // 3 components, r theta z
         m_all_field_functors[lev][icomp] =
-            std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, dim), lev,
-                              m_crse_ratio, false, ncomp_multimodefab);
+            std::make_unique<JFunctor>(dim, lev, m_crse_ratio, false, deposit_current, ncomp_multimodefab);
+        deposit_current = false;
         icomp += 1;
         AddRZModesToOutputNames(std::string("J") + coord[dim],
                                 warpx.get_pointer_current_fp(0, 0)->nComp());
@@ -583,6 +600,10 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
     const auto nspec = static_cast<int>(m_pfield_species.size());
     const auto ntot = static_cast<int>(nvar + m_pfield_varnames.size() * nspec);
 
+    // Boolean flag for whether the current density should be deposited before
+    // diagnostic output
+    bool deposit_current = !m_solver_deposits_current;
+
     m_all_field_functors[lev].resize(ntot);
     // Fill vector of functors for all components except individual cylindrical modes.
     for (int comp=0; comp<nvar; comp++){
@@ -591,7 +612,8 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
         } else if ( m_varnames[comp] == "Bz" ){
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 2), lev, m_crse_ratio);
         } else if ( m_varnames[comp] == "jz" ){
-            m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 2), lev, m_crse_ratio);
+            m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(2, lev, m_crse_ratio, true, deposit_current);
+            deposit_current = false;
         } else if ( m_varnames[comp] == "Az" ){
             m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_vector_potential_fp(lev, 2), lev, m_crse_ratio);
         } else if ( m_varnames[comp] == "rho" ){
@@ -628,9 +650,11 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
             } else if ( m_varnames[comp] == "Bt" ){
                 m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 1), lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "jr" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 0), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(0, lev, m_crse_ratio, true, deposit_current);
+                deposit_current = false;
             } else if ( m_varnames[comp] == "jt" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 1), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(1, lev, m_crse_ratio, true, deposit_current);
+                deposit_current = false;
             } else if ( m_varnames[comp] == "Ar" ){
                 m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_vector_potential_fp(lev, 0), lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "At" ){
@@ -649,9 +673,11 @@ FullDiagnostics::InitializeFieldFunctors (int lev)
             } else if ( m_varnames[comp] == "By" ){
                 m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_Bfield_aux(lev, 1), lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "jx" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 0), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(0, lev, m_crse_ratio, true, deposit_current);
+                deposit_current = false;
             } else if ( m_varnames[comp] == "jy" ){
-                m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_current_fp(lev, 1), lev, m_crse_ratio);
+                m_all_field_functors[lev][comp] = std::make_unique<JFunctor>(1, lev, m_crse_ratio, true, deposit_current);
+                deposit_current = false;
             } else if ( m_varnames[comp] == "Ax" ){
                 m_all_field_functors[lev][comp] = std::make_unique<CellCenterFunctor>(warpx.get_pointer_vector_potential_fp(lev, 0), lev, m_crse_ratio);
             } else if ( m_varnames[comp] == "Ay" ){