add collider-relevant reduced diags (ECP-WarpX#4024)

* added collider reduced diags * added test * fixed xy_ave and xy_std * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fixed usage of IntVect and removed print * Set up automated CI test To-do: - Fix bug (erroneous arithmetic operation) - Remove unused parameters from input file * Fix warning: set mass only if WARPX_QED * Use amrex::ParticleReal for mass * Update Make.package to fix GNU Make build * reduced_data.value() only once * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Merge `development` into `coll_rel_red_diags` * updated 3d beam test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * updated 3d test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * `m_write_header` instead of `m_IsNotRestart` * added angles * updated 3d beam test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * added 3 particle test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * added warning * updated diags names * updated analysis multiple particles * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * added positrons to test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Update CI test * Reuse input file parser from Tools * Apply suggestions from code review * Apply suggestions from code review * changed order of diags and added docs * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Apply WarpX style conventions * Fix CodeQL alerts * assert in RZ and removed fine ref levs * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix compilation in 2D * Fix compilation in RZ * Update Docs/source/usage/parameters.rst Co-authored-by: Axel Huebl <[email protected]> * Update Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp Co-authored-by: Axel Huebl <[email protected]> * Update Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp Co-authored-by: Luca Fedeli <[email protected]> * Update Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp Co-authored-by: Luca Fedeli <[email protected]> * Update Source/Diagnostics/ReducedDiags/ColliderRelevant.H Co-authored-by: Axel Huebl <[email protected]> * Update Source/Diagnostics/ReducedDiags/ColliderRelevant.H Co-authored-by: Axel Huebl <[email protected]> * minimized number of kernel launches and parallel ops * updated constructor of ReduceDiags * fixed based on comments * updated docs: if QED not enable * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fixed ReducedDiags constructor * two passes instead of one for robusteness * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fixed bugs * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Fix 1D build (typo `wtot` instead of `w_tot`) * refixed compile in RZ * check io process in chi_ave computation * updated test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Temporary fix: remove IOProcessor * Update Source/Diagnostics/ReducedDiags/ColliderRelevant.cpp Co-authored-by: Axel Huebl <[email protected]> * Apply suggestions from code review * Remove `Debug` compilation mode from CI test * CI test analysis: check all particles have same charge * Use `m_beam_name` instead of `species_names` --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Edoardo Zoni <[email protected]> Co-authored-by: Edoardo Zoni <[email protected]> Co-authored-by: Axel Huebl <[email protected]> Co-authored-by: Luca Fedeli <[email protected]>
eebasso · Sep 11, 2023 · ac52a76 · ac52a76
1 parent 2e4d6fd
commit ac52a76
Show file tree

Hide file tree

Showing 14 changed files with 999 additions and 2 deletions.
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -2988,6 +2988,57 @@ Reduced Diagnostics
         1 or 0, it is possible to compute the charge on only some part of the
         embedded boundary.
 
+    * ``ColliderRelevant``
+        This diagnostics computes properties of two colliding beams that are relevant for particle colliders.
+        Two species must be specified. Photon species are not supported yet.
+        It is assumed that the two species propagate and collide along the ``z`` direction.
+        The output columns (for 3D-XYZ) are the following, where the minimum, average and maximum
+        are done over the whole species:
+
+        [0]: simulation step (iteration).
+
+        [1]: time (s).
+
+        [2]: time derivative of the luminosity (:math:`m^{-2}s^{-1}`) defined as:
+
+        .. math::
+
+            \frac{dL}{dt} = 2 c \iiint  n_1(x,y,z) n_2(x,y,z) dx dy dz
+
+        where :math:`n_1`, :math:`n_2` are the number densities of the two colliding species.
+
+        [3], [4], [5]: If, QED is enabled, the minimum, average and maximum values of the quantum parameter :math:`\chi` of species 1:
+        :math:`\chi_{min}`,
+        :math:`\langle \chi \rangle`,
+        :math:`\chi_{max}`.
+        If QED is not enabled, these numbers are not computed.
+
+        [6], [7]: The average and standard deviation of the values of the transverse coordinate :math:`x` (m) of species 1:
+        :math:`\langle x \rangle`,
+        :math:`\sqrt{\langle x- \langle x \rangle \rangle^2}`.
+
+        [8], [9]: The average and standard deviation of the values of the transverse coordinate :math:`y` (m) of species 1:
+        :math:`\langle y \rangle`,
+        :math:`\sqrt{\langle y- \langle y \rangle \rangle^2}`.
+
+        [10], [11], [12], [13]: The minimum, average, maximum and standard deviation of the angle :math:`\theta_x = \angle (u_x, u_z)` (rad) of species 1:
+        :math:`{\theta_x}_{min}`,
+        :math:`\langle \theta_x \rangle`,
+        :math:`{\theta_x}_{max}`,
+        :math:`\sqrt{\langle \theta_x- \langle \theta_x \rangle \rangle^2}`.
+
+        [14], [15], [16], [17]:  The minimum, average, maximum and standard deviation of the angle :math:`\theta_y = \angle (u_y, u_z)` (rad) of species 1:
+        :math:`{\theta_y}_{min}`,
+        :math:`\langle \theta_y \rangle`,
+        :math:`{\theta_y}_{max}`,
+        :math:`\sqrt{\langle \theta_y- \langle \theta_y \rangle \rangle^2}`.
+
+        [18], ..., [32]: Analogous quantities for species 2.
+
+        For 2D-XZ, :math:`y`-related quantities are not outputted.
+        For 1D-Z, :math:`x`-related and :math:`y`-related quantities are not outputted.
+        RZ geometry is not supported yet.
+
 * ``<reduced_diags_name>.intervals`` (`string`)
     Using the `Intervals Parser`_ syntax, this string defines the timesteps at which reduced
     diagnostics are written to file.

diff --git a/Examples/Tests/collider_relevant_diags/analysis_multiple_particles.py b/Examples/Tests/collider_relevant_diags/analysis_multiple_particles.py
@@ -0,0 +1,150 @@
+#!/usr/bin/env python3
+
+import os
+import sys
+
+import numpy as np
+import openpmd_api as io
+import pandas as pd
+from scipy.constants import c, e, hbar, m_e
+
+sys.path.append('../../../../warpx/Regression/Checksum/')
+import checksumAPI
+
+sys.path.append('../../../../warpx/Tools/Parser/')
+from input_file_parser import parse_input_file
+
+E_crit = m_e**2*c**3/(e*hbar)
+B_crit = m_e**2*c**2/(e*hbar)
+
+def chi(ux, uy, uz, Ex, Ey, Ez, Bx, By, Bz):
+    gamma = np.sqrt(1.+ux**2+uy**2+uz**2)
+    vx = ux / gamma * c
+    vy = uy / gamma * c
+    vz = uz / gamma * c
+    tmp1x = Ex + vy*Bz - vz*By
+    tmp1y = Ey - vx*Bz + vz*Bx
+    tmp1z = Ez + vx*By - vy*Bx
+    tmp2 = (Ex*vx + Ey*vy + Ez*vz)/c
+    chi = gamma/E_crit*np.sqrt(tmp1x**2+tmp1y**2+tmp1z**2 - tmp2**2)
+    return chi
+
+def dL_dt():
+    series = io.Series("diags/diag2/openpmd_%T.h5",io.Access.read_only)
+    iterations = np.asarray(series.iterations)
+    lumi = []
+    for n,ts in enumerate(iterations):
+        it = series.iterations[ts]
+        rho1 = it.meshes["rho_beam_e"]
+        dV = np.prod(rho1.grid_spacing)
+        rho1 = it.meshes["rho_beam_e"][io.Mesh_Record_Component.SCALAR].load_chunk()
+        rho2 = it.meshes["rho_beam_p"][io.Mesh_Record_Component.SCALAR].load_chunk()
+        beam_e_charge = it.particles["beam_e"]["charge"][io.Mesh_Record_Component.SCALAR].load_chunk()
+        beam_p_charge = it.particles["beam_p"]["charge"][io.Mesh_Record_Component.SCALAR].load_chunk()
+        q1 = beam_e_charge[0]
+        if not np.all(beam_e_charge == q1):
+            sys.exit('beam_e particles do not have the same charge')
+        q2 = beam_p_charge[0]
+        if not np.all(beam_p_charge == q2):
+            sys.exit('beam_p particles do not have the same charge')
+        series.flush()
+        n1 = rho1/q1
+        n2 = rho2/q2
+        l = 2*np.sum(n1*n2)*dV*c
+        lumi.append(l)
+    return lumi
+
+input_dict = parse_input_file('inputs_3d_multiple_particles')
+Ex, Ey, Ez = [float(w) for w in input_dict['particles.E_external_particle']]
+Bx, By, Bz = [float(w) for w in input_dict['particles.B_external_particle']]
+
+CollDiagFname='diags/reducedfiles/ColliderRelevant_beam_e_beam_p.txt'
+df = pd.read_csv(CollDiagFname, sep=" ", header=0)
+
+for species in ['beam_p', 'beam_e']:
+
+    ux1, ux2, ux3 = [float(w) for w in input_dict[f'{species}.multiple_particles_ux']]
+    uy1, uy2, uy3 = [float(w) for w in input_dict[f'{species}.multiple_particles_uy']]
+    uz1, uz2, uz3 = [float(w) for w in input_dict[f'{species}.multiple_particles_uz']]
+
+    x = np.array([float(w) for w in input_dict[f'{species}.multiple_particles_pos_x']])
+    y = np.array([float(w) for w in input_dict[f'{species}.multiple_particles_pos_y']])
+
+    w = np.array([float(w) for w in input_dict[f'{species}.multiple_particles_weight']])
+
+    CHI_ANALYTICAL = np.array([chi(ux1, uy1, uz1, Ex, Ey, Ez, Bx, By, Bz),
+                               chi(ux2, uy2, uz2, Ex, Ey, Ez, Bx, By, Bz),
+                               chi(ux3, uy3, uz3, Ex, Ey, Ez, Bx, By, Bz)])
+    THETAX = np.array([np.arctan2(ux1, uz1), np.arctan2(ux2, uz2), np.arctan2(ux3, uz3)])
+    THETAY = np.array([np.arctan2(uy1, uz1), np.arctan2(uy2, uz2), np.arctan2(uy3, uz3)])
+
+    # CHI MAX
+    fname=f'diags/reducedfiles/ParticleExtrema_{species}.txt'
+    chimax_pe = np.loadtxt(fname)[:,19]
+    chimax_cr = df[[col for col in df.columns if f'chi_max_{species}' in col]].to_numpy()
+    assert np.allclose(np.max(CHI_ANALYTICAL), chimax_cr, rtol=1e-8)
+    assert np.allclose(chimax_pe, chimax_cr, rtol=1e-8)
+
+    # CHI MIN
+    fname=f'diags/reducedfiles/ParticleExtrema_{species}.txt'
+    chimin_pe = np.loadtxt(fname)[:,18]
+    chimin_cr = df[[col for col in df.columns if f'chi_min_{species}' in col]].to_numpy()
+    assert np.allclose(np.min(CHI_ANALYTICAL), chimin_cr, rtol=1e-8)
+    assert np.allclose(chimin_pe, chimin_cr, rtol=1e-8)
+
+    # CHI AVERAGE
+    chiave_cr = df[[col for col in df.columns if f'chi_ave_{species}' in col]].to_numpy()
+    assert np.allclose(np.average(CHI_ANALYTICAL, weights=w), chiave_cr, rtol=1e-8)
+
+    # X AVE STD
+    x_ave_cr = df[[col for col in df.columns if f']x_ave_{species}' in col]].to_numpy()
+    x_std_cr = df[[col for col in df.columns if f']x_std_{species}' in col]].to_numpy()
+    x_ave = np.average(x, weights=w)
+    x_std = np.sqrt(np.average((x-x_ave)**2, weights=w))
+    assert np.allclose(x_ave, x_ave_cr, rtol=1e-8)
+    assert np.allclose(x_std, x_std_cr, rtol=1e-8)
+
+    # Y AVE STD
+    y_ave_cr = df[[col for col in df.columns if f']y_ave_{species}' in col]].to_numpy()
+    y_std_cr = df[[col for col in df.columns if f']y_std_{species}' in col]].to_numpy()
+    y_ave = np.average(y, weights=w)
+    y_std = np.sqrt(np.average((y-y_ave)**2, weights=w))
+    assert np.allclose(y_ave, y_ave_cr, rtol=1e-8)
+    assert np.allclose(y_std, y_std_cr, rtol=1e-8)
+
+    # THETA X MIN AVE MAX STD
+    thetax_min_cr = df[[col for col in df.columns if f'theta_x_min_{species}' in col]].to_numpy()
+    thetax_ave_cr = df[[col for col in df.columns if f'theta_x_ave_{species}' in col]].to_numpy()
+    thetax_max_cr = df[[col for col in df.columns if f'theta_x_max_{species}' in col]].to_numpy()
+    thetax_std_cr = df[[col for col in df.columns if f'theta_x_std_{species}' in col]].to_numpy()
+    thetax_min = np.min(THETAX)
+    thetax_ave = np.average(THETAX, weights=w)
+    thetax_max = np.max(THETAX)
+    thetax_std = np.sqrt(np.average((THETAX-thetax_ave)**2, weights=w))
+    assert np.allclose(thetax_min, thetax_min_cr, rtol=1e-8)
+    assert np.allclose(thetax_ave, thetax_ave_cr, rtol=1e-8)
+    assert np.allclose(thetax_max, thetax_max_cr, rtol=1e-8)
+    assert np.allclose(thetax_std, thetax_std_cr, rtol=1e-8)
+
+    # THETA Y MIN AVE MAX STD
+    thetay_min_cr = df[[col for col in df.columns if f'theta_y_min_{species}' in col]].to_numpy()
+    thetay_ave_cr = df[[col for col in df.columns if f'theta_y_ave_{species}' in col]].to_numpy()
+    thetay_max_cr = df[[col for col in df.columns if f'theta_y_max_{species}' in col]].to_numpy()
+    thetay_std_cr = df[[col for col in df.columns if f'theta_y_std_{species}' in col]].to_numpy()
+    thetay_min = np.min(THETAY)
+    thetay_ave = np.average(THETAY, weights=w)
+    thetay_max = np.max(THETAY)
+    thetay_std = np.sqrt(np.average((THETAY-thetay_ave)**2, weights=w))
+    assert np.allclose(thetay_min, thetay_min_cr, rtol=1e-8)
+    assert np.allclose(thetay_ave, thetay_ave_cr, rtol=1e-8)
+    assert np.allclose(thetay_max, thetay_max_cr, rtol=1e-8)
+    assert np.allclose(thetay_std, thetay_std_cr, rtol=1e-8)
+
+    # dL/dt
+    dL_dt_cr = df[[col for col in df.columns if 'dL_dt' in col]].to_numpy()
+    assert np.allclose(dL_dt_cr, dL_dt(), rtol=1e-8)
+
+# Checksum analysis
+plotfile = sys.argv[1]
+test_name = os.path.split(os.getcwd())[1]
+checksumAPI.evaluate_checksum(test_name, plotfile)
diff --git a/Examples/Tests/collider_relevant_diags/inputs_3d_multiple_particles b/Examples/Tests/collider_relevant_diags/inputs_3d_multiple_particles
@@ -0,0 +1,130 @@
+#################################
+########## MY CONSTANTS #########
+#################################
+my_constants.nx = 8
+my_constants.ny = 8
+my_constants.nz = 8
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+max_step = 1
+amr.n_cell = nx ny nz
+amr.max_grid_size = 4
+amr.blocking_factor = 4
+amr.max_level = 0
+geometry.dims = 3
+geometry.prob_lo = 0 0 0
+geometry.prob_hi =  8 8 8
+particles.do_tiling = 0
+warpx.use_filter = 0
+
+#################################
+######## BOUNDARY CONDITION #####
+#################################
+boundary.field_lo = periodic periodic periodic
+boundary.field_hi = periodic periodic periodic
+boundary.particle_lo = periodic periodic periodic
+boundary.particle_hi = periodic periodic periodic
+
+#################################
+############ NUMERICS ###########
+#################################
+algo.maxwell_solver = ckc
+warpx.cfl = 0.99
+algo.particle_shape = 1
+
+#################################
+############ FIELDS #############
+#################################
+particles.E_ext_particle_init_style = constant
+particles.B_ext_particle_init_style = constant
+particles.E_external_particle = 10000. 0. 0.
+particles.B_external_particle = 0. 5000. 0.
+
+#################################
+########### PARTICLES ###########
+#################################
+particles.species_names = pho beam_p beam_e
+particles.photon_species = pho
+
+beam_e.species_type = electron
+beam_e.injection_style = MultipleParticles
+beam_e.multiple_particles_pos_x = 4.5  3.5  0.5
+beam_e.multiple_particles_pos_y = 4.5  2.5  1.5
+beam_e.multiple_particles_pos_z = 4.5  1.5  1.5
+beam_e.multiple_particles_ux = 0.3 0.2 0.1
+beam_e.multiple_particles_uy = 0.4 -0.3 -0.1
+beam_e.multiple_particles_uz = 0.3 0.1 -10.
+beam_e.multiple_particles_weight = 1. 2 3
+beam_e.initialize_self_fields = 0
+beam_e.self_fields_required_precision = 5e-10
+beam_e.do_qed_quantum_sync = 1
+beam_e.qed_quantum_sync_phot_product_species = pho
+beam_e.do_not_push = 1
+beam_e.do_not_deposit = 1
+
+beam_p.species_type = positron
+beam_p.injection_style = MultipleParticles
+beam_p.multiple_particles_pos_x = 4.5 3.5 0.5
+beam_p.multiple_particles_pos_y = 4.5 2.5 1.5
+beam_p.multiple_particles_pos_z = 4.5 1.5 1.5
+beam_p.multiple_particles_ux = 0.3 0.2 0.1
+beam_p.multiple_particles_uy = 0.4 -0.3 -0.1
+beam_p.multiple_particles_uz = 0.3 0.1 -10.
+beam_p.multiple_particles_weight = 1. 2 3
+beam_p.initialize_self_fields = 0
+beam_p.self_fields_required_precision = 5e-10
+beam_p.do_qed_quantum_sync = 1
+beam_p.qed_quantum_sync_phot_product_species = pho
+beam_p.do_not_push = 1
+beam_p.do_not_deposit = 1
+
+pho.species_type = photon
+pho.injection_style = none
+
+#################################
+############# QED ###############
+#################################
+qed_qs.photon_creation_energy_threshold = 0.
+qed_qs.lookup_table_mode = builtin
+qed_qs.chi_min = 1.e-3
+warpx.do_qed_schwinger = 0
+
+#################################
+######### DIAGNOSTICS ###########
+#################################
+# FULL
+diagnostics.diags_names = diag1 diag2
+
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.write_species = 1
+diag1.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho_beam_e rho_beam_p rho
+diag1.species = pho beam_e beam_p
+diag1.format = plotfile
+#diag1.dump_last_timestep = 1
+
+diag2.intervals = 1
+diag2.diag_type = Full
+diag2.write_species = 1
+diag2.fields_to_plot = Ex Ey Ez Bx By Bz jx jy jz rho_beam_e rho_beam_p rho
+diag2.species = pho beam_e beam_p
+diag2.format = openpmd
+diag2.openpmd_backend = h5
+#diag2.dump_last_timestep = 1
+
+# REDUCED
+warpx.reduced_diags_names = ParticleExtrema_beam_e ParticleExtrema_beam_p ColliderRelevant_beam_e_beam_p
+
+ColliderRelevant_beam_e_beam_p.type = ColliderRelevant
+ColliderRelevant_beam_e_beam_p.intervals = 1
+ColliderRelevant_beam_e_beam_p.species =beam_e beam_p
+
+ParticleExtrema_beam_e.type = ParticleExtrema
+ParticleExtrema_beam_e.intervals = 1
+ParticleExtrema_beam_e.species = beam_e
+
+ParticleExtrema_beam_p.type = ParticleExtrema
+ParticleExtrema_beam_p.intervals = 1
+ParticleExtrema_beam_p.species = beam_p
diff --git a/Regression/Checksum/benchmarks_json/collider_diagnostics.json b/Regression/Checksum/benchmarks_json/collider_diagnostics.json
@@ -0,0 +1,36 @@
+{
+  "lev=0": {
+    "Bx": 0.0,
+    "By": 0.0,
+    "Bz": 0.0,
+    "Ex": 0.0,
+    "Ey": 0.0,
+    "Ez": 0.0,
+    "jx": 0.0,
+    "jy": 0.0,
+    "jz": 0.0,
+    "rho": 0.0,
+    "rho_beam_e": 9.613059803999999e-19,
+    "rho_beam_p": 9.613059803999999e-19
+  },
+  "beam_e": {
+    "particle_momentum_x": 1.638554718442694e-22,
+    "particle_momentum_y": 2.184739624590259e-22,
+    "particle_momentum_z": 2.8401615119673365e-21,
+    "particle_opticalDepthQSR": 9.767741201839083e+00,
+    "particle_position_x": 8.5,
+    "particle_position_y": 8.5,
+    "particle_position_z": 7.5,
+    "particle_weight": 6.0
+  },
+  "beam_p": {
+    "particle_momentum_x": 1.638554718442694e-22,
+    "particle_momentum_y": 2.184739624590259e-22,
+    "particle_momentum_z": 2.8401615119673365e-21,
+    "particle_opticalDepthQSR": 8.773870620305825e+00,
+    "particle_position_x": 8.5,
+    "particle_position_y": 8.5,
+    "particle_position_z": 7.5,
+    "particle_weight": 6.0
+  }
+}
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
@@ -209,6 +209,22 @@ compileTest = 0
 doVis = 0
 analysisRoutine = Examples/Tests/btd_rz/analysis_BTD_laser_antenna.py
 
+[collider_diagnostics]
+buildDir = .
+inputFile = Examples/Tests/collider_relevant_diags/inputs_3d_multiple_particles
+runtime_params = warpx.abort_on_warning_threshold=high
+dim = 3
+addToCompileString =
+cmakeSetupOpts = -DWarpX_DIMS=3
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+compileTest = 0
+doVis = 0
+analysisRoutine = Examples/Tests/collider_relevant_diags/analysis_multiple_particles.py
+
 [collisionISO]
 buildDir = .
 inputFile = Examples/Tests/collision/inputs_3d_isotropization