PSATD: Implement First-Order Equations (#3466)

* Implement First-Order PSATD Equations

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

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

* Fix Unused Parameter Warning

* Fix RZ Build

* Fix Normalization of G to Match PML

* Add CI Test: 3D Uniform Plasma

* Cleaning

* Update 2D CI Checksums

* Update 3D CI Checksums

* Add F,G to CI Checksums of `uniform_plasma_multiJ`

* Allow User to Choose First-Order v. Second-Order

* Add WARPX_ALWAYS_ASSERT_WITH_MESSAGE

* Rename New Class `PsatdAlgorithmFirstOrder`

* Remove Inline Comment

* Update RZ CI Checksums

* Fix inline comment

* Use auxiliary variables to avoid divisions

* Use auxiliary variables to avoid divisions

* Make `nci_psatd_stability` dir and merge inputs

* Move all Galilean tests to `nci_psatd_stability`

* Fix CI

* Fix index for backward FFT of J

Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com>

1 files changed, 49 insertions, 0 deletions

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)