Diagnostics : EM fields on particles

Docs/source/usage/parameters.rst

@@ -2784,7 +2784,7 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a
 * ``<diag_name>.openpmd_encoding`` (optional, ``v`` (variable based), ``f`` (file based) or ``g`` (group based) ) only read if ``<diag_name>.format = openpmd``.
      openPMD `file output encoding <https://openpmd-api.readthedocs.io/en/0.16.1/usage/concepts.html#iteration-and-series>`__.
      File based: one file per timestep (slower), group/variable based: one file for all steps (faster)).
-     ``variable based`` is an `experimental feature with ADIOS2 <https://openpmd-api.readthedocs.io/en/0.16.1/backends/adios2.html#experimental-new-adios2-schema>`__ and not supported for back-transformed diagnostics.
+     ``variable based`` is an `experimental feature with ADIOS2 <https://openpmd-api.readthedocs.io/en/0.16.1/backends/adios2.html#experimental-new-adios2-schema>`__ and not supported for back-transformed diagnostics. This format is also not supported by OpenPMDTimeSeries at the moment, the script :download:`read_variable_based_adios2.py <../../../Tools/PostProcessing/read_variable_based_adios2.py>` provides basic functions to read variable-based particle diagnostics directly with the ``adios2`` module.
      Default: ``f`` (full diagnostics)
 
 * ``<diag_name>.adios2_operator.type`` (``zfp``, ``blosc``) optional,
@@ -2916,9 +2916,9 @@ In-situ capabilities can be used by turning on Sensei or Ascent (provided they a
 * ``<diag_name>.<species_name>.variables`` (list of `strings` separated by spaces, optional)
     List of particle quantities to write to output.
     Choices are ``x``, ``y``, ``z`` for the particle positions (3D and RZ), ``x`` & ``z`` in 2D, ``z`` in 1D,
-    ``w`` for the particle weight and ``ux``, ``uy``, ``uz`` for the particle momenta.
+    ``w`` for the particle weight, ``ux``, ``uy``, ``uz`` for the particle momenta, and ``EM`` for the electromagnetic fields.
     When using the lab-frame electrostatic solver, ``phi`` (electrostatic potential, on the macroparticles) is also available.
-    By default, all particle quantities (except ``phi``) are written.
+    By default, all particle quantities (except ``phi`` and ``EM``) are written.
     If ``<diag_name>.<species_name>.variables = none``, no particle data are written.
 
 * ``<diag_name>.<species_name>.random_fraction`` (`float`) optional

Examples/Tests/CMakeLists.txt

@@ -19,6 +19,7 @@ add_subdirectory(embedded_boundary_em_particle_absorption)
 add_subdirectory(embedded_boundary_python_api)
 add_subdirectory(embedded_boundary_rotated_cube)
 add_subdirectory(embedded_circle)
+add_subdirectory(EM_fields_on_particles)
 add_subdirectory(energy_conserving_thermal_plasma)
 add_subdirectory(field_probe)
 add_subdirectory(flux_injection)

Examples/Tests/EM_fields_on_particles/CMakeLists.txt

@@ -0,0 +1,11 @@
+# add tests
+
+add_warpx_test(
+    test_2d_particle_EM_diagnostics # name
+    2 # dims
+    2 # nprocs
+    inputs_test_2d_particle_EM_diagnostics # inputs
+    "analysis_2d.py" # analysis
+    "analysis_default_regression.py --path diags/diag_checksum/" # checksum
+    OFF # dependency
+)

Examples/Tests/EM_fields_on_particles/analysis_2d.py

@@ -0,0 +1,158 @@
+#!/usr/bin/env python3
+
+import numpy as np
+from openpmd_viewer import OpenPMDTimeSeries
+from scipy.interpolate import interp1d
+from scipy.optimize import minimize
+from scipy.signal import hilbert
+
+c = 3e8
+lambda_laser = 800e-9
+T_laser = lambda_laser / c
+
+E_max = 3.26e10
+T_peak = 10e-15
+tau = 5e-15
+delta_t = 1.6e-6 / c  # time for the laser to reach the particle
+
+ts_particle = OpenPMDTimeSeries("diags/diag1/")
+
+ex_t, ey_t, ez_t, bx_t, by_t, bz_t = ts_particle.iterate(
+    ts_particle.get_particle, ["ex", "ey", "ez", "bx", "by", "bz"], species="electron"
+)
+
+ex_t = np.squeeze(ex_t)
+ey_t = np.squeeze(ey_t)
+ez_t = np.squeeze(ez_t)
+bx_t = np.squeeze(bx_t)
+by_t = np.squeeze(by_t)
+bz_t = np.squeeze(bz_t)
+
+n_diags = 401
+
+DT = 6.324524234e-17  # @ CFL = 0.99
+iterations = np.linspace(0, n_diags - 1, n_diags)
+T = (iterations * DT) - DT / 2
+
+T_peak_part = T_peak + delta_t
+Ey_max = E_max * 1 / np.sqrt(5)  # polarisation vector (0 1 2)
+Ez_max = E_max * np.sqrt(1 - 1 / 5)
+By_max = Ez_max / c
+Bz_max = Ey_max / c
+
+
+def get_laser_th(E_0, T_p, T, tau, T_laser, phi=0):
+    alpha = np.exp(-((T - T_peak_part) ** 2) / (tau**2))
+    return E_0 * alpha * np.cos(2 * np.pi * (T - T_p) / T_laser + phi)
+
+
+Ey_th = get_laser_th(Ey_max, T_peak_part, T, tau, T_laser)
+Ez_th = get_laser_th(Ez_max, T_peak_part, T, tau, T_laser)
+By_th = get_laser_th(By_max, T_peak_part, T, tau, T_laser, phi=np.pi)
+Bz_th = get_laser_th(Bz_max, T_peak_part, T, tau, T_laser)
+
+Ey_th = np.where(T >= delta_t, Ey_th, 0)
+Ez_th = np.where(T >= delta_t, Ez_th, 0)
+By_th = np.where(T >= delta_t, By_th, 0)
+Bz_th = np.where(T >= delta_t, Bz_th, 0)
+
+
+# due to injection method the simulated fields have a slight dephasing
+# wrt the theoretical field
+# here we compute the dephasing and correct it on the theoretical field
+def align_cost(shift, sig1, sig2, time):
+    # Create an interpolation function for the signal to be shifted
+    interp_func = interp1d(time, sig1, kind="linear", fill_value="extrapolate")
+    shifted_sig1 = interp_func(time + shift)
+    # MSE
+    cost = np.mean((shifted_sig1 - sig2) ** 2)
+    return cost
+
+
+init_shift = 0.0
+res = minimize(
+    align_cost,
+    init_shift,
+    args=(Ey_th, ey_t, T),
+    method="Nelder-Mead",
+    options={"xatol": 1e-20, "fatol": 1e-20, "maxiter": 10000},
+)
+opt_shift = res.x[0]
+
+interp_func_ey = interp1d(T, Ey_th, kind="cubic", fill_value="extrapolate")
+interp_func_ez = interp1d(T, Ez_th, kind="cubic", fill_value="extrapolate")
+interp_func_by = interp1d(T, By_th, kind="cubic", fill_value="extrapolate")
+interp_func_bz = interp1d(T, Bz_th, kind="cubic", fill_value="extrapolate")
+
+Ey_aligned = interp_func_ey(T + opt_shift)
+Ez_aligned = interp_func_ez(T + opt_shift)
+By_aligned = interp_func_by(T + opt_shift)
+Bz_aligned = interp_func_bz(T + opt_shift)
+
+# verif
+
+# E.1
+
+assert np.allclose(np.zeros(ex_t.shape), ex_t, rtol=0, atol=5e-4)
+
+# E.2
+
+ey_p = ey_t[T >= delta_t]
+Ey_p = Ey_aligned[T >= delta_t]
+ey_p += 1e11
+Ey_p += 1e11
+
+ez_p = ez_t[T >= delta_t]
+Ez_p = Ez_aligned[T >= delta_t]
+ez_p += 1e11
+Ez_p += 1e11
+
+assert np.allclose(ey_p, Ey_p, rtol=8e-4, atol=0)
+assert np.allclose(ez_p, Ez_p, rtol=2e-3, atol=0)
+
+# E.3
+
+ey_m = ey_t[T < delta_t - DT * 4]
+ez_m = ez_t[T < delta_t - DT * 4]
+assert np.allclose(ey_m, np.zeros(ey_m.shape), atol=1e-10, rtol=0)
+assert np.allclose(ez_m, np.zeros(ey_m.shape), atol=1e-10, rtol=0)
+
+# E.4
+
+env_ez_t = np.abs(hilbert(ez_t))
+env_ey_t = np.abs(hilbert(ey_t))
+assert np.isclose(np.max(env_ez_t) / np.max(env_ey_t), 2, rtol=2e-3, atol=0)
+assert np.allclose(2 * env_ey_t[85:390], env_ez_t[85:390], rtol=8e-3, atol=0)
+
+# B.1
+
+assert np.allclose(np.zeros(bx_t.shape), bx_t, rtol=0, atol=5e-4)
+
+# B.2
+
+by_p = by_t[T >= delta_t]
+By_p = By_aligned[T >= delta_t]
+by_p += 1e3
+By_p += 1e3
+
+bz_p = bz_t[T >= delta_t]
+Bz_p = Bz_aligned[T >= delta_t]
+bz_p += 1e3
+Bz_p += 1e3
+
+assert np.allclose(by_p, By_p, rtol=8e-4, atol=0)
+assert np.allclose(bz_p, Bz_p, rtol=4e-4, atol=0)
+
+# B.3
+
+by_m = by_t[T < delta_t - DT * 4]
+bz_m = bz_t[T < delta_t - DT * 4]
+assert np.allclose(by_m, np.zeros(by_m.shape), atol=1e-10, rtol=0)
+assert np.allclose(bz_m, np.zeros(by_m.shape), atol=1e-10, rtol=0)
+
+# B.4
+
+env_bz_t = np.abs(hilbert(bz_t))
+env_by_t = np.abs(hilbert(by_t))
+assert np.isclose(np.max(env_by_t) / np.max(env_bz_t), 2, rtol=2e-3, atol=0)
+assert np.allclose(0.5 * env_by_t[85:390], env_bz_t[85:390], rtol=8e-3, atol=0)

Examples/Tests/EM_fields_on_particles/analysis_default_regression.py

@@ -0,0 +1 @@
+../../analysis_default_regression.py

Examples/Tests/EM_fields_on_particles/inputs_test_2d_particle_EM_diagnostics

@@ -0,0 +1,84 @@
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+
+# Parameters
+max_step = 400
+amr.n_cell = 256 256
+warpx.numprocs = 2 1
+amr.max_level = 0
+geometry.dims = 2
+geometry.prob_lo = 0 0   # physical domain
+geometry.prob_hi = 5e-6 2.5e-5
+
+# BC
+
+boundary.field_lo = pml pml
+boundary.field_hi = pml pml
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.serialize_initial_conditions = 1
+warpx.verbose = 1
+warpx.cfl = 0.99
+warpx.use_filter = 0
+
+# Order of particle shape factors
+algo.particle_shape = 2
+
+#################################
+########### "PLASMA" ############
+#################################
+
+particles.species_names = electron
+
+electron.species_type = electron
+electron.injection_style = "SingleParticle"
+electron.single_particle_pos = 2.5e-6 0 1.25e-5
+electron.single_particle_u = 0 0 0
+electron.single_particle_weight = 0
+
+electron.do_not_push = 1
+electron.do_not_deposit = 1
+
+#################################
+######### LASER #################
+#################################
+
+lasers.names = las1
+
+las1.position = 0.9e-6 326 1.25e-5
+las1.polarization = 0 1 2
+las1.direction = 1 0 0
+las1.e_max = 3.26e10
+las1.wavelength = 800e-9
+las1.profile = "Gaussian"
+las1.profile_t_peak = 10e-15
+las1.profile_duration = 5e-15
+las1.profile_waist = 5e-6
+las1.profile_focal_distance = 1.6e-6
+
+#################################
+############# DIAGS #############
+#################################
+
+diagnostics.diags_names = diag1 diag_checksum
+
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.format = openpmd
+diag1.openpmd_backend = h5
+diag1.fields_to_plot = none
+diag1.write_species = 1
+diag1.species = electron
+diag1.electron.variables = EM # will emit a warning because no position
+
+diag_checksum.intervals = 400:400
+diag_checksum.diag_type = Full
+diag_checksum.format = openpmd
+diag_checksum.openpmd_backend = h5
+diag_checksum.fields_to_plot = Ex Ey Ez Bx By Bz
+diag_checksum.write_species = 1
+diag_checksum.species = electron
+diag_checksum.electron.variables = x z ux uy uz w EM # will emit a warning because no position

Regression/Checksum/benchmarks_json/test_2d_particle_EM_diagnostics.json

@@ -0,0 +1,25 @@
+{
+  "electron": {
+    "particle_bx": 3.890980310565517e-07,
+    "particle_by": 0.34760630481899957,
+    "particle_bz": 0.1737192637399916,
+    "particle_ex": 0.00013839455641573295,
+    "particle_ey": 50608854.55959022,
+    "particle_ez": 111088826.0524848,
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 0.0,
+    "particle_position_x": 2.5e-06,
+    "particle_position_y": 0.0,
+    "particle_position_z": 1.25e-05,
+    "particle_weight": 0.0
+  },
+  "lev=0": {
+    "Bx": 3864.690332686649,
+    "By": 239757.36697091258,
+    "Bz": 119878.42982987352,
+    "Ex": 2396689176285.191,
+    "Ey": 35934018061852.125,
+    "Ez": 71884648325008.73
+  }
+}

Source/Diagnostics/ParticleDiag/ParticleDiag.H

@@ -24,6 +24,7 @@ public:
     [[nodiscard]] std::string getSpeciesName() const { return m_name; }
     amrex::Vector<int> m_plot_flags;
     bool m_plot_phi = false; // Whether to output the potential phi on the particles
+    bool m_plot_EM = false; // Whether to output the E and B fields on the particles
 
     bool m_do_random_filter  = false;
     bool m_do_uniform_filter = false;

Source/Diagnostics/ParticleDiag/ParticleDiag.cpp

@@ -42,10 +42,9 @@ ParticleDiag::ParticleDiag (
                 if (var == "y") { var = "theta"; }
 #endif
                 if (var == "phi") {
-                    // User requests phi on particle. This is *not* part of the variables that
-                    // the particle container carries, and is only added to particles during output.
-                    // Therefore, this case needs to be treated specifically.
                     m_plot_phi = true;
+                } else if (var=="EM") {
+                    m_plot_EM = true;
                 } else {
                     WARPX_ALWAYS_ASSERT_WITH_MESSAGE(pc->HasRealComp(var),
                         "variables argument '" + var