Adding current density to Field Reduction reduced diagnostic (#3980)

* Adding j in FieldReduction

* Modiyfying input file with FieldReduced diag

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

* Update Source/Diagnostics/ReducedDiags/FieldReduction.H

* updating Field Reduction Docs to include current density

* Update Source/Diagnostics/ReducedDiags/FieldReduction.H

Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com>

* adding cuurent density BackwardCompatibility to Field Reduction diagnostic

* including E dot j density using the modifyed Field Reduction

* fixing picmi interface

* adding Edotj to automated python test

* [pre-commit.ci] auto fixes from pre-commit.com hooks

for more information, see https://pre-commit.ci

---------

Co-authored-by: Lizzette Corrales <lzcorrales@kong.dhcp.lbl.gov>
Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>
Co-authored-by: Neïl Zaim <49716072+NeilZaim@users.noreply.github.com>

6 files changed, 64 insertions, 12 deletions

diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
index 91239e24f..f7eba5583 100644
--- a/Docs/source/usage/parameters.rst
+++ b/Docs/source/usage/parameters.rst
@@ -2675,9 +2675,9 @@ Reduced Diagnostics
         are computed.
 
     * ``FieldReduction``
-        This type computes an arbitrary reduction of the positions and the electromagnetic fields.
+        This type computes an arbitrary reduction of the positions, the current density, and the electromagnetic fields.
 
-        * ``<reduced_diags_name>.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz)`` (`string`)
+        * ``<reduced_diags_name>.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz)`` (`string`)
             An analytic function to be reduced must be provided, using the math parser.
 
         * ``<reduced_diags_name>.reduction_type`` (`string`)
diff --git a/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py b/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py
index 76732aeec..93352d29f 100755
--- a/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py
+++ b/Examples/Tests/reduced_diags/analysis_reduced_diags_impl.py
@@ -166,6 +166,9 @@ def do_analysis(single_precision = False):
     Bx = ad[('mesh','Bx')].to_ndarray()
     By = ad[('mesh','By')].to_ndarray()
     Bz = ad[('mesh','Bz')].to_ndarray()
+    jx = ad[('mesh','jx')].to_ndarray()
+    jy = ad[('mesh','jy')].to_ndarray()
+    jz = ad[('mesh','jz')].to_ndarray()
     rho = ad[('boxlib','rho')].to_ndarray()
     rho_electrons = ad[('boxlib','rho_electrons')].to_ndarray()
     rho_protons = ad[('boxlib','rho_protons')].to_ndarray()
@@ -203,6 +206,7 @@ def do_analysis(single_precision = False):
     values_yt['protons: maximum of |rho|'] = np.amax(np.abs(rho_protons))
     values_yt['maximum of |B| from generic field reduction'] = np.amax(np.sqrt(Bx**2 + By**2 + Bz**2))
     values_yt['minimum of x*Ey*Bz'] = np.amin(x*Ey*Bz)
+    values_yt['maximum of Edotj'] = np.amax(Ex*jx + Ey*jy + Ez*jz)
 
     #--------------------------------------------------------------------------------------------------
     # Part 2: get results from reduced diagnostics (label '_rd')
@@ -222,6 +226,7 @@ def do_analysis(single_precision = False):
     FR_Maxdata = np.genfromtxt('./diags/reducedfiles/FR_Max.txt')  # Field Reduction using maximum
     FR_Mindata = np.genfromtxt('./diags/reducedfiles/FR_Min.txt')  # Field Reduction using minimum
     FR_Integraldata = np.genfromtxt('./diags/reducedfiles/FR_Integral.txt')  # Field Reduction using integral
+    Edotjdata = np.genfromtxt('./diags/reducedfiles/Edotj.txt') # E dot j maximum
 
     # First index "1" points to the values written at the last time step
     values_rd['field energy'] = EFdata[1][2]
@@ -283,6 +288,7 @@ def do_analysis(single_precision = False):
     values_rd['photons: sum of weights'] = NPdata[1][9]
     values_rd['maximum of |B| from generic field reduction'] = FR_Maxdata[1][2]
     values_rd['minimum of x*Ey*Bz'] = FR_Mindata[1][2]
+    values_rd['maximum of Edotj'] = Edotjdata[1][2]
 
     #--------------------------------------------------------------------------------------------------
     # Part 3: compare values from plotfiles and reduced diagnostics and print output
diff --git a/Examples/Tests/reduced_diags/inputs b/Examples/Tests/reduced_diags/inputs
index 9b8e45d3e..dc0c57264 100644
--- a/Examples/Tests/reduced_diags/inputs
+++ b/Examples/Tests/reduced_diags/inputs
@@ -68,7 +68,7 @@ photons.uz_th = 0.2
 #################################
 ###### REDUCED DIAGS ############
 #################################
-warpx.reduced_diags_names = EP NP EF PP PF MF MR FP FP_integrate FP_line FP_plane FR_Max FR_Min FR_Integral
+warpx.reduced_diags_names = EP NP EF PP PF MF MR FP FP_integrate FP_line FP_plane FR_Max FR_Min FR_Integral Edotj
 EP.type = ParticleEnergy
 EP.intervals = 200
 EF.type = FieldEnergy
@@ -122,17 +122,21 @@ NP.type = ParticleNumber
 NP.intervals = 200
 FR_Max.type = FieldReduction
 FR_Max.intervals = 200
-FR_Max.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz) = sqrt(Bx**2 + By**2 + Bz**2)
+FR_Max.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz) = sqrt(Bx**2 + By**2 + Bz**2)
 FR_Max.reduction_type = Maximum
 FR_Min.type = FieldReduction
 FR_Min.intervals = 200
-FR_Min.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz) = x*Ey*Bz
+FR_Min.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz) = x*Ey*Bz
 FR_Min.reduction_type = Minimum
 FR_Integral.type = FieldReduction
 FR_Integral.intervals = 200
-FR_Integral.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz) = "if(y > 0 and z < 0,
+FR_Integral.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz) = "if(y > 0 and z < 0,
                             0.5*((Ex**2 + Ey**2 + Ez**2)*epsilon0+(Bx**2 + By**2 + Bz**2)/mu0), 0)"
 FR_Integral.reduction_type = Integral
+Edotj.type = FieldReduction
+Edotj.intervals = 200
+Edotj.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz) = Ex*jx + Ey*jy + Ez*jz
+Edotj.reduction_type = Maximum
 
 # Diagnostics
 diagnostics.diags_names = diag1
diff --git a/Python/pywarpx/picmi.py b/Python/pywarpx/picmi.py
index 0f0426ca6..efa14badd 100644
--- a/Python/pywarpx/picmi.py
+++ b/Python/pywarpx/picmi.py
@@ -2478,7 +2478,7 @@ class ReducedDiagnostic(picmistandard.base._ClassWithInit, WarpXDiagnosticBase):
         self.reduction_type = kw.pop("reduction_type")
         reduced_function = kw.pop("reduced_function")
 
-        self.__setattr__("reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz)", reduced_function)
+        self.__setattr__("reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz)", reduced_function)
 
         # Check the reduced function expression for constants
         for k in list(kw.keys()):
diff --git a/Source/Diagnostics/ReducedDiags/FieldReduction.H b/Source/Diagnostics/ReducedDiags/FieldReduction.H
index 46fcfc295..c4c81c016 100644
--- a/Source/Diagnostics/ReducedDiags/FieldReduction.H
+++ b/Source/Diagnostics/ReducedDiags/FieldReduction.H
@@ -60,10 +60,16 @@ public:
      */
     virtual void ComputeDiags(int step) override final;
 
+    /**
+     * This function queries deprecated input parameters and aborts
+     * the run if one of them is specified.
+     */
+    void BackwardCompatibility();
+
 private:
     /// Parser to read expression to be reduced from the input file.
-    /// 9 elements are x, y, z, Ex, Ey, Ez, Bx, By, Bz
-    static constexpr int m_nvars = 9;
+    /// 12 elements are x, y, z, Ex, Ey, Ez, Bx, By, Bz, jx, jy, jz
+    static constexpr int m_nvars = 12;
     std::unique_ptr<amrex::Parser> m_parser;
 
     // Type of reduction (e.g. Maximum, Minimum or Sum)
@@ -99,6 +105,10 @@ public:
         const amrex::MultiFab & Bx = warpx.getBfield(lev,0);
         const amrex::MultiFab & By = warpx.getBfield(lev,1);
         const amrex::MultiFab & Bz = warpx.getBfield(lev,2);
+        const amrex::MultiFab & jx = warpx.getcurrent_fp(lev,0);
+        const amrex::MultiFab & jy = warpx.getcurrent_fp(lev,1);
+        const amrex::MultiFab & jz = warpx.getcurrent_fp(lev,2);
+
 
         // General preparation of interpolation and reduction operations
         const amrex::GpuArray<int,3> cellCenteredtype{0,0,0};
@@ -118,6 +128,9 @@ public:
         auto Bxtype = amrex::GpuArray<int,3>{0,0,0};
         auto Bytype = amrex::GpuArray<int,3>{0,0,0};
         auto Bztype = amrex::GpuArray<int,3>{0,0,0};
+        auto jxtype = amrex::GpuArray<int,3>{0,0,0};
+        auto jytype = amrex::GpuArray<int,3>{0,0,0};
+        auto jztype = amrex::GpuArray<int,3>{0,0,0};
         for (int i = 0; i < AMREX_SPACEDIM; ++i){
             Extype[i] = Ex.ixType()[i];
             Eytype[i] = Ey.ixType()[i];
@@ -125,6 +138,10 @@ public:
             Bxtype[i] = Bx.ixType()[i];
             Bytype[i] = By.ixType()[i];
             Bztype[i] = Bz.ixType()[i];
+            jxtype[i] = jx.ixType()[i];
+            jytype[i] = jy.ixType()[i];
+            jztype[i] = jz.ixType()[i];
+
         }
 
         // get parser
@@ -145,6 +162,9 @@ public:
             const auto& arrBx = Bx[mfi].array();
             const auto& arrBy = By[mfi].array();
             const auto& arrBz = Bz[mfi].array();
+            const auto& arrjx = jx[mfi].array();
+            const auto& arrjy = jy[mfi].array();
+            const auto& arrjz = jz[mfi].array();
 
             reduce_op.eval(box, reduce_data,
             [=] AMREX_GPU_DEVICE (int i, int j, int k) -> ReduceTuple
@@ -176,8 +196,16 @@ public:
                                                                                reduction_coarsening_ratio, i, j, k, reduction_comp);
                 const amrex::Real Bz_interp = ablastr::coarsen::sample::Interp(arrBz, Bztype, cellCenteredtype,
                                                                                reduction_coarsening_ratio, i, j, k, reduction_comp);
+                const amrex::Real jx_interp = ablastr::coarsen::sample::Interp(arrjx, jxtype, cellCenteredtype,
+                                                                               reduction_coarsening_ratio, i, j, k, reduction_comp);
+                const amrex::Real jy_interp = ablastr::coarsen::sample::Interp(arrjy, jytype, cellCenteredtype,
+                                                                               reduction_coarsening_ratio, i, j, k, reduction_comp);
+                const amrex::Real jz_interp = ablastr::coarsen::sample::Interp(arrjz, jztype, cellCenteredtype,
+                                                                               reduction_coarsening_ratio, i, j, k, reduction_comp);
+
                 return reduction_function_parser(x, y, z, Ex_interp, Ey_interp, Ez_interp,
-                                                          Bx_interp, By_interp, Bz_interp);
+                        Bx_interp, By_interp, Bz_interp,
+                        jx_interp, jy_interp, jz_interp);
             });
         }
 
diff --git a/Source/Diagnostics/ReducedDiags/FieldReduction.cpp b/Source/Diagnostics/ReducedDiags/FieldReduction.cpp
index f4b521c10..32c8c9653 100644
--- a/Source/Diagnostics/ReducedDiags/FieldReduction.cpp
+++ b/Source/Diagnostics/ReducedDiags/FieldReduction.cpp
@@ -44,14 +44,16 @@ FieldReduction::FieldReduction (std::string rd_name)
     // resize data array
     m_data.resize(noutputs, 0.0_rt);
 
+    BackwardCompatibility();
+
     const amrex::ParmParse pp_rd_name(rd_name);
 
     // read reduced function with parser
     std::string parser_string = "";
-    utils::parser::Store_parserString(pp_rd_name,"reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz)",
+    utils::parser::Store_parserString(pp_rd_name,"reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz)",
                        parser_string);
     m_parser = std::make_unique<amrex::Parser>(
-        utils::parser::makeParser(parser_string,{"x","y","z","Ex","Ey","Ez","Bx","By","Bz"}));
+        utils::parser::makeParser(parser_string,{"x","y","z","Ex","Ey","Ez","Bx","By","Bz","jx","jy","jz"}));
 
     // Replace all newlines and possible following whitespaces with a single whitespace. This
     // should avoid weird formatting when the string is written in the header of the output file.
@@ -85,6 +87,18 @@ FieldReduction::FieldReduction (std::string rd_name)
 }
 // end constructor
 
+void FieldReduction::BackwardCompatibility ()
+{
+    amrex::ParmParse pp_rd_name(m_rd_name);
+    std::vector<std::string> backward_strings;
+    WARPX_ALWAYS_ASSERT_WITH_MESSAGE(
+        !pp_rd_name.queryarr("reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz)", backward_strings),
+        "<reduced_diag_name>.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz) is no longer a valid option. "
+        "Please use the renamed option <reduced_diag_name>.reduced_function(x,y,z,Ex,Ey,Ez,Bx,By,Bz,jx,jy,jz) instead."
+    );
+}
+
+
 // function that does an arbitrary reduction of the electromagnetic fields
 void FieldReduction::ComputeDiags (int step)
 {