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>

2 files changed, 14 insertions, 4 deletions

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