CI: Test New v. Legacy BTD in RigidInjection_BTD (#3327)

* Improve Testing of New v. Legacy BTD

* openPMD: Flush Series before Accessing Arrays

Co-authored-by: Ryan Sandberg <[email protected]>

* Add Checksums for RigidInjection_BTD (BTD data)

* Compare Also Particle Momenta (BTD data)

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

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

* Apply suggestions from code review

Co-authored-by: Ryan Sandberg <[email protected]>

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

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

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

Examples/Modules/RigidInjection/analysis_rigid_injection_BoostedFrame.py

@@ -14,40 +14,89 @@
 A Gaussian electron beam starts from -5 microns, propagates rigidly up to
 20 microns after which it expands due to emittance only (the focal position is
 20 microns). The beam width is measured after ~50 microns, and compared with
-the theory (with a 5% error allowed).
+the theory (with a 1% relative error allowed).
 
 The simulation runs in a boosted frame, and the analysis is done in the lab
 frame, i.e., on the back-transformed diagnostics.
 '''
 
+import os
+import sys
+
 import numpy as np
+import openpmd_api as io
 import read_raw_data
+from scipy.constants import m_e
 import yt
 
 yt.funcs.mylog.setLevel(0)
 
-# Read data from back-transformed diagnostics
+sys.path.insert(1, '../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+filename = sys.argv[1]
+
+# Tolerances to check consistency between legacy BTD and new BTD
+rtol = 1e-16
+atol = 1e-16
+
+# Read data from legacy back-transformed diagnostics
 snapshot = './lab_frame_data/snapshots/snapshot00001'
-header   = './lab_frame_data/snapshots/Header'
-allrd, info = read_raw_data.read_lab_snapshot(snapshot, header)
-z = np.mean( read_raw_data.get_particle_field(snapshot, 'beam', 'z') )
-w = np.std ( read_raw_data.get_particle_field(snapshot, 'beam', 'x') )
+x_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'x')
+z_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'z')
+ux_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'ux')
+uy_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'uy')
+uz_legacy = read_raw_data.get_particle_field(snapshot, 'beam', 'uz')
+
+# Read data from new back-transformed diagnostics (plotfile)
+ds_plotfile = yt.load(filename)
+x_plotfile = ds_plotfile.all_data()['beam', 'particle_position_x'].v
+z_plotfile = ds_plotfile.all_data()['beam', 'particle_position_y'].v
+ux_plotfile = ds_plotfile.all_data()['beam', 'particle_momentum_x'].v
+uy_plotfile = ds_plotfile.all_data()['beam', 'particle_momentum_y'].v
+uz_plotfile = ds_plotfile.all_data()['beam', 'particle_momentum_z'].v
 
-# initial parameters
+# Read data from new back-transformed diagnostics (openPMD)
+series = io.Series("./diags/diag2/openpmd_%T.h5", io.Access.read_only)
+ds_openpmd = series.iterations[1]
+x_openpmd = ds_openpmd.particles['beam']['position']['x'][:]
+z_openpmd = ds_openpmd.particles['beam']['position']['z'][:]
+ux_openpmd = ds_openpmd.particles['beam']['momentum']['x'][:]
+uy_openpmd = ds_openpmd.particles['beam']['momentum']['y'][:]
+uz_openpmd = ds_openpmd.particles['beam']['momentum']['z'][:]
+series.flush()
+
+# Sort and compare arrays to check consistency between legacy BTD and new BTD (plotfile)
+assert(np.allclose(np.sort(x_legacy), np.sort(x_plotfile), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(z_legacy), np.sort(z_plotfile), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(ux_legacy*m_e), np.sort(ux_plotfile), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(uy_legacy*m_e), np.sort(uy_plotfile), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(uz_legacy*m_e), np.sort(uz_plotfile), rtol=rtol, atol=atol))
+
+# Sort and compare arrays to check consistency between legacy BTD and new BTD (openPMD)
+assert(np.allclose(np.sort(x_legacy), np.sort(x_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(z_legacy), np.sort(z_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(ux_legacy*m_e), np.sort(ux_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(uy_legacy*m_e), np.sort(uy_openpmd), rtol=rtol, atol=atol))
+assert(np.allclose(np.sort(uz_legacy*m_e), np.sort(uz_openpmd), rtol=rtol, atol=atol))
+
+# Initial parameters
 z0 = 20.e-6
-w0 = 1.e-6
+x0 = 1.e-6
 theta0 = np.arcsin(0.1)
 
-# Theoretical beam width after propagation if rigid ON
-wth = np.sqrt( w0**2 + (z-z0)**2*theta0**2 )
-error_rel = np.abs((w-wth)/wth)
-tolerance_rel = 0.03
-
-# Print error and assert small error
-print("Beam position: " + str(z))
-print("Beam width   : " + str(w))
+# Theoretical beam width after propagation with rigid injection
+z = np.mean(z_legacy)
+x = np.std(x_legacy)
+print(f'Beam position = {z}')
+print(f'Beam width    = {x}')
 
-print("error_rel    : " + str(error_rel))
-print("tolerance_rel: " + str(tolerance_rel))
+xth = np.sqrt(x0**2 + (z-z0)**2 * theta0**2)
+err = np.abs((x-xth) / xth)
+tol = 1e-2
+print(f'error = {err}')
+print(f'tolerance = {tol}')
+assert(err < tol)
 
-assert( error_rel < tolerance_rel )
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, filename)

Examples/Modules/RigidInjection/inputs_2d_BoostedFrame

@@ -49,27 +49,25 @@ beam.zinject_plane = 20.e-6
 beam.rigid_advance = true
 
 # Diagnostics
-diagnostics.diags_names = diag1 btd_pltfile btd_openpmd
-diag1.intervals = 10000
-diag1.diag_type = Full
+diagnostics.diags_names = diag1 diag2
 
-btd_openpmd.diag_type = BackTransformed
-btd_openpmd.do_back_transformed_fields = 1
-btd_openpmd.num_snapshots_lab = 2
-btd_openpmd.dt_snapshots_lab = 1.8679589331096515e-13
-btd_openpmd.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
-btd_openpmd.format = openpmd
-btd_openpmd.openpmd_backend = h5
-btd_openpmd.buffer_size = 32
+diag1.diag_type = BackTransformed
+diag1.do_back_transformed_fields = 1
+diag1.num_snapshots_lab = 2
+diag1.dt_snapshots_lab = 1.8679589331096515e-13
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
+diag1.format = plotfile
+diag1.buffer_size = 32
+diag1.write_species = 1
 
-btd_pltfile.diag_type = BackTransformed
-btd_pltfile.do_back_transformed_fields = 1
-btd_pltfile.num_snapshots_lab = 2
-btd_pltfile.dt_snapshots_lab = 1.8679589331096515e-13
-btd_pltfile.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
-btd_pltfile.format = plotfile
-btd_pltfile.buffer_size = 32
-btd_pltfile.write_species = 1
+diag2.diag_type = BackTransformed
+diag2.do_back_transformed_fields = 1
+diag2.num_snapshots_lab = 2
+diag2.dt_snapshots_lab = 1.8679589331096515e-13
+diag2.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho
+diag2.format = openpmd
+diag2.openpmd_backend = h5
+diag2.buffer_size = 32
 
 # old BTD diagnostics
 warpx.do_back_transformed_diagnostics = 1

Regression/Checksum/benchmarks_json/RigidInjection_BTD.json

@@ -0,0 +1,22 @@
+{
+  "beam": {
+    "particle_momentum_x": 2.2080215038948936e-16,
+    "particle_momentum_y": 2.18711072170811e-16,
+    "particle_momentum_z": 2.730924530737497e-15,
+    "particle_position_x": 0.0260823588888081,
+    "particle_position_y": 0.5049438607316916,
+    "particle_weight": 62415.090744607645
+  },
+  "lev=0": {
+    "Bx": 3.721807007218884e-05,
+    "By": 0.004860056238272468,
+    "Bz": 5.5335765596325185e-06,
+    "Ex": 1466447.517373168,
+    "Ey": 11214.10223280318,
+    "Ez": 283216.0961218869,
+    "jx": 16437877.898892513,
+    "jy": 2492340.3149980404,
+    "jz": 215102423.57036877,
+    "rho": 0.7246235591902177
+  }
+}