8 files changed, 106 insertions, 60 deletions
diff --git a/Examples/Tests/galilean/analysis.py b/Examples/Tests/nci_psatd_stability/analysis_galilean.py
index 9fe6bab72..9fe6bab72 100755
--- a/Examples/Tests/galilean/analysis.py
+++ b/Examples/Tests/nci_psatd_stability/analysis_galilean.py
diff --git a/Examples/Tests/nci_psatd_stability/analysis_multiJ.py b/Examples/Tests/nci_psatd_stability/analysis_multiJ.py
new file mode 100755
index 000000000..1c68b114c
--- /dev/null
+++ b/Examples/Tests/nci_psatd_stability/analysis_multiJ.py
@@ -0,0 +1,49 @@
+#!/usr/bin/env python3
+"""
+This script is used to test the results of the multi-J PSATD
+first-order equations, with one J deposition. It compares the
+energy of the electric field with a given reference energy.
+
+The reference energy is computed by running the same test with J constant
+in time, rho linear in time, and without divergence cleaning. The reference
+energy corresponds to unstable results due to NCI (suppressed by the use of
+both J and rho constant in time, and with divergence cleaning).
+"""
+import os
+import sys
+
+import numpy as np
+import scipy.constants as scc
+
+import yt ; yt.funcs.mylog.setLevel(0)
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+filename = sys.argv[1]
+
+ds = yt.load(filename)
+
+# yt 4.0+ has rounding issues with our domain data:
+# RuntimeError: yt attempted to read outside the boundaries
+# of a non-periodic domain along dimension 0.
+if 'force_periodicity' in dir(ds): ds.force_periodicity()
+
+all_data = ds.covering_grid(level=0, left_edge=ds.domain_left_edge, dims=ds.domain_dimensions)
+Ex = all_data['boxlib', 'Ex'].squeeze().v
+Ey = all_data['boxlib', 'Ey'].squeeze().v
+Ez = all_data['boxlib', 'Ez'].squeeze().v
+
+# Set reference energy values, and tolerances for numerical stability and charge conservation
+tol_energy = 1e-8
+energy_ref = 66e6
+
+# Check numerical stability by comparing electric field energy to reference energy
+energy = np.sum(scc.epsilon_0/2*(Ex**2+Ey**2+Ez**2))
+err_energy = energy / energy_ref
+print('\nCheck numerical stability:')
+print(f'err_energy = {err_energy}')
+print(f'tol_energy = {tol_energy}')
+assert(err_energy < tol_energy)
+
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)
diff --git a/Examples/Tests/galilean/inputs_2d b/Examples/Tests/nci_psatd_stability/inputs_2d
index f19bfdde5..f19bfdde5 100644
--- a/Examples/Tests/galilean/inputs_2d
+++ b/Examples/Tests/nci_psatd_stability/inputs_2d
diff --git a/Examples/Tests/galilean/inputs_2d_hybrid b/Examples/Tests/nci_psatd_stability/inputs_2d_hybrid
index a9d1ecd14..a9d1ecd14 100644
--- a/Examples/Tests/galilean/inputs_2d_hybrid
+++ b/Examples/Tests/nci_psatd_stability/inputs_2d_hybrid
diff --git a/Examples/Tests/galilean/inputs_3d b/Examples/Tests/nci_psatd_stability/inputs_3d
index 7a0d80748..53305af12 100644
--- a/Examples/Tests/galilean/inputs_3d
+++ b/Examples/Tests/nci_psatd_stability/inputs_3d
@@ -1,84 +1,81 @@
-#################################
-####### GENERAL PARAMETERS ######
-#################################
-max_step = 300
-
-amr.n_cell = 32 32 32
-warpx.numprocs = 1 1 2
-amr.max_level = 0
-psatd.v_galilean = 0. 0. 0.99498743710662
-
-geometry.dims    = 3
-geometry.prob_lo = -9.67 -9.67 -9.67
-geometry.prob_hi = 9.67 9.67 9.67
-
-#################################
-###### Boundary condition #######
-#################################
-boundary.field_lo = periodic periodic periodic
-boundary.field_hi = periodic periodic periodic
-
-#################################
-############ NUMERICS ###########
-#################################
-warpx.verbose = 1
-
+# algo
 algo.current_deposition = direct
 algo.maxwell_solver = psatd
 algo.particle_pusher = vay
-
-warpx.cfl = 1.
 algo.particle_shape = 3
 
-#################################
-############ PLASMA #############
-#################################
-particles.species_names = electrons ions
+# amr
+amr.max_level = 0
+amr.n_cell = 32 32 32
 
-warpx.do_nodal = 1
-warpx.use_filter = 1
+# boundary
+boundary.field_hi = periodic periodic periodic
+boundary.field_lo = periodic periodic periodic
 
-psatd.nox = 8
-psatd.noy = 8
-psatd.noz = 8
+# diag1
+diag1.diag_type = Full
+diag1.intervals = 300
 
+# diagnostics
+diagnostics.diags_names = diag1
+
+# electrons
 electrons.charge = -q_e
-electrons.mass = m_e
+electrons.density = 282197938148984.7
 electrons.injection_style = "NUniformPerCell"
+electrons.mass = m_e
+electrons.momentum_distribution_type = "gaussian"
 electrons.num_particles_per_cell_each_dim = 1 1 1
 electrons.profile = constant
-electrons.density = 282197938148984.7
-electrons.momentum_distribution_type = "gaussian"
-electrons.uz_m = 9.9498743710661994
-electrons.xmin = -9.67
-electrons.xmax = 9.67
-electrons.ymin = -9.67
-electrons.ymax = 9.67
-electrons.zmin =    -9.67
-electrons.zmax =    9.67
 electrons.ux_th = 0.0001
 electrons.uy_th = 0.0001
+electrons.uz_m = 9.9498743710661994
 electrons.uz_th = 0.0001
+electrons.xmax =  9.67
+electrons.xmin = -9.67
+electrons.ymax =  9.67
+electrons.ymin = -9.67
+electrons.zmax =  9.67
+electrons.zmin = -9.67
+
+# geometry
+geometry.dims = 3
+geometry.prob_hi =  9.67  9.67  9.67
+geometry.prob_lo = -9.67 -9.67 -9.67
 
+# ions
 ions.charge = q_e
-ions.mass = m_p
+ions.density = 282197938148984.7
 ions.injection_style = "NUniformPerCell"
+ions.mass = m_p
+ions.momentum_distribution_type = "gaussian"
 ions.num_particles_per_cell_each_dim = 1 1 1
 ions.profile = constant
-ions.density =  282197938148984.7
-ions.momentum_distribution_type = "gaussian"
-ions.uz_m = 9.9498743710661994
-ions.xmin = -9.67
-ions.xmax = 9.67
-ions.ymin = -9.67
-ions.ymax = 9.67
-ions.zmin =  -9.67
-ions.zmax = 9.67
 ions.ux_th = 0.0001
 ions.uy_th = 0.0001
+ions.uz_m = 9.9498743710661994
 ions.uz_th = 0.0001
+ions.xmax =  9.67
+ions.xmin = -9.67
+ions.ymax =  9.67
+ions.ymin = -9.67
+ions.zmax =  9.67
+ions.zmin = -9.67
 
-# Diagnostics
-diagnostics.diags_names = diag1
-diag1.intervals = 100
-diag1.diag_type = Full
+# max_step
+max_step = 300
+
+# particles
+particles.species_names = electrons ions
+
+# psatd
+psatd.nox = 8
+psatd.noy = 8
+psatd.noz = 8
+
+# warpx
+warpx.cfl = 1.
+warpx.do_nodal = 1
+warpx.numprocs = 1 1 2
+warpx.use_filter = 1
+warpx.verbose = 1
diff --git a/Examples/Tests/averaged_galilean/inputs_avg_2d b/Examples/Tests/nci_psatd_stability/inputs_avg_2d
index 2a2611906..2a2611906 100644
--- a/Examples/Tests/averaged_galilean/inputs_avg_2d
+++ b/Examples/Tests/nci_psatd_stability/inputs_avg_2d
diff --git a/Examples/Tests/averaged_galilean/inputs_avg_3d b/Examples/Tests/nci_psatd_stability/inputs_avg_3d
index bb76e4c47..bb76e4c47 100644
--- a/Examples/Tests/averaged_galilean/inputs_avg_3d
+++ b/Examples/Tests/nci_psatd_stability/inputs_avg_3d
diff --git a/Examples/Tests/galilean/inputs_rz b/Examples/Tests/nci_psatd_stability/inputs_rz
index 772711ddc..772711ddc 100644
--- a/Examples/Tests/galilean/inputs_rz
+++ b/Examples/Tests/nci_psatd_stability/inputs_rz