Add differential luminosity reduced diagnostic (ECP-WarpX#5161)

* adds differential luminosity red diag * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * fix bugs and add test * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * add missing dt factor * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Change formatting of energy bins * Add profiling * Replace return with continue * Add energy spread * Add automated test * Remove unused warnings * Update Warning level * Update scripts and analysis * [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci * Avoid copy data to GPU at each step * Make analysis.py executable * Change file format * Fix const-ness * Add aux. file * Do computation at every step, and dump only at certain intervals * Fix clang-tidy errors * Update checksum * Update documentation * Fix documentation * Update documentation * Update comments --------- Co-authored-by: pre-commit-ci[bot] <66853113+pre-commit-ci[bot]@users.noreply.github.com> Co-authored-by: Remi Lehe <[email protected]>
JustinRayAngus · Aug 27, 2024 · 0089ba6 · 0089ba6
1 parent 143bd4e
commit 0089ba6
Show file tree

Hide file tree

Showing 10 changed files with 591 additions and 0 deletions.
diff --git a/Docs/source/usage/parameters.rst b/Docs/source/usage/parameters.rst
@@ -3435,6 +3435,40 @@ Reduced Diagnostics
         For 1D-Z, :math:`x`-related and :math:`y`-related quantities are not outputted.
         RZ geometry is not supported yet.
 
+* ``DifferentialLuminosity``
+    This type computes the differential luminosity between two species, defined as:
+
+    .. math::
+
+        \frac{d\mathcal{L}}{d\mathcal{E}^*}(\mathcal{E}^*, t) = \int_0^t dt'\int d\boldsymbol{x}\,d\boldsymbol{p}_1 d\boldsymbol{p}_2\;
+         \sqrt{ |\boldsymbol{v}_1 - \boldsymbol{v}_2|^2 - |\boldsymbol{v}_1\times\boldsymbol{v}_2|^2/c^2} \\ f_1(\boldsymbol{x}, \boldsymbol{p}_1, t')f_2(\boldsymbol{x}, \boldsymbol{p}_2, t') \delta(\mathcal{E}^* - \mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2))
+
+    where :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2(m_1 m_2 c^4
+    \gamma_1 \gamma_2 - \boldsymbol{p}_1\cdot\boldsymbol{p}_2 c^2)}` is the energy in the center-of-mass frame,
+    and :math:`f_i` is the distribution function of species :math:`i`. Note that, if :math:`\sigma^*(\mathcal{E}^*)`
+    is the center-of-mass cross-section of a given collision process, then
+    :math:`\int d\mathcal{E}^* \frac{d\mathcal{L}}{d\mathcal{E}^*} (\mathcal{E}^*, t)\sigma^*(\mathcal{E}^*)`
+    gives the total number of collisions of that process (from the beginning of the simulation up until time :math:`t`).
+
+    The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{J}^{-1}`. For collider-relevant WarpX simulations
+    involving two crossing, high-energy beams of particles, the differential luminosity in :math:`\text{s}^{-1}.\text{m}^{-2}.\text{J}^{-1}`
+    can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams.
+
+    In practice, the above expression of the differential luminosity is evaluated over discrete bins in energy :math:`\mathcal{E}^*`,
+    and by summing over macroparticles.
+
+    * ``<reduced_diags_name>.species`` (`list of two strings`)
+        The names of the two species for which the differential luminosity is computed.
+
+    * ``<reduced_diags_name>.bin_number`` (`int` > 0)
+        The number of bins in energy :math:`\mathcal{E}^*`
+
+    * ``<reduced_diags_name>.bin_max`` (`float`, in Joules)
+        The minimum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
+
+    * ``<reduced_diags_name>.bin_min`` (`float`, in Joules)
+        The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
+
 * ``<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/diff_lumi_diag/analysis.py b/Examples/Tests/diff_lumi_diag/analysis.py
@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+# This test checks the differential luminosity of the beam-beam interaction
+# in the case of two Gaussian beams crossing rigidly at the interaction point.
+# The beams have a Gaussian distribution both in energy and in transverse positions.
+# In that case, the differential luminosity can be calculated analytically.
+
+import os
+import sys
+
+import numpy as np
+from read_raw_data import read_reduced_diags_histogram
+from scipy.constants import eV
+
+sys.path.insert(1, "../../../../warpx/Regression/Checksum/")
+import checksumAPI
+
+# Extract the differential luminosity from the file
+_, _, E_bin, bin_data = read_reduced_diags_histogram(
+    "./diags/reducedfiles/DifferentialLuminosity_beam1_beam2.txt"
+)
+dL_dE_sim = bin_data[-1]  # Differential luminosity at the end of the simulation
+
+# Beam parameters
+N = 1.2e10
+E_beam = 125e9 * eV
+sigma_x = 500e-9
+sigma_y = 10e-9
+
+# Compute the analytical differential luminosity for 2 Gaussian beams
+sigma_E1 = 0.02 * E_beam
+sigma_E2 = 0.03 * E_beam
+sigma_E = np.sqrt(
+    sigma_E1**2 + sigma_E2**2
+)  # Expected width of the differential luminosity
+dL_dE_th = (
+    N**2
+    / (2 * (2 * np.pi) ** 1.5 * sigma_x * sigma_y * sigma_E)
+    * np.exp(-((E_bin - 2 * E_beam) ** 2) / (2 * sigma_E**2))
+)
+
+# Check that the simulation result and analytical result match
+error = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max()
+tol = 1e-2
+print("Relative error: ", error)
+print("Tolerance: ", tol)
+assert error < tol
+
+# Get name of the test
+fn = sys.argv[1]
+test_name = os.path.split(os.getcwd())[1]
+
+# Run checksum regression test
+checksumAPI.evaluate_checksum(test_name, fn, rtol=1e-2)
diff --git a/Examples/Tests/diff_lumi_diag/inputs b/Examples/Tests/diff_lumi_diag/inputs
@@ -0,0 +1,127 @@
+#################################
+########## MY CONSTANTS #########
+#################################
+my_constants.mc2 = m_e*clight*clight
+my_constants.GeV = q_e*1.e9
+
+# BEAMS
+my_constants.beam_energy = 125.*GeV
+my_constants.beam_gamma = beam_energy/(mc2)
+my_constants.beam_charge = 1.2e10*q_e
+my_constants.sigmax = 500e-9
+my_constants.sigmay = 10e-9
+my_constants.sigmaz = 300e-3
+my_constants.muz = -4*sigmaz
+my_constants.nmacropart = 2e5
+
+# BOX
+my_constants.Lx = 8*sigmax
+my_constants.Ly = 8*sigmay
+my_constants.Lz = 16*sigmaz
+
+my_constants.nx = 64
+my_constants.ny = 64
+my_constants.nz = 128
+
+# TIME
+my_constants.T = 0.5*Lz/clight
+my_constants.dt = sigmaz/clight/10.
+
+#################################
+####### GENERAL PARAMETERS ######
+#################################
+stop_time = T
+amr.n_cell = nx ny nz
+amr.max_grid_size = 128
+amr.blocking_factor = 2
+amr.max_level = 0
+geometry.dims = 3
+geometry.prob_lo = -0.5*Lx -0.5*Ly -0.5*Lz
+geometry.prob_hi =  0.5*Lx  0.5*Ly  0.5*Lz
+
+#################################
+######## BOUNDARY CONDITION #####
+#################################
+boundary.field_lo = open open open
+boundary.field_hi = open open open
+boundary.particle_lo = Absorbing Absorbing Absorbing
+boundary.particle_hi = Absorbing Absorbing Absorbing
+
+#################################
+############ NUMERICS ###########
+#################################
+warpx.do_electrostatic = relativistic
+warpx.const_dt = dt
+warpx.grid_type = collocated
+algo.particle_shape = 3
+algo.load_balance_intervals=100
+algo.particle_pusher = vay
+warpx.poisson_solver = fft
+
+#################################
+########### PARTICLES ###########
+#################################
+particles.species_names = beam1 beam2
+
+beam1.species_type = electron
+beam1.injection_style = gaussian_beam
+beam1.x_rms = sigmax
+beam1.y_rms = sigmay
+beam1.z_rms = sigmaz
+beam1.x_m = 0
+beam1.y_m = 0
+beam1.z_m = muz
+beam1.npart = nmacropart
+beam1.q_tot = -beam_charge
+beam1.z_cut = 4
+beam1.momentum_distribution_type = gaussian
+beam1.uz_m = beam_gamma
+beam1.uy_m = 0.0
+beam1.ux_m = 0.0
+beam1.ux_th = 0
+beam1.uy_th = 0
+beam1.uz_th = 0.02*beam_gamma
+beam1.do_not_deposit = 1
+
+beam2.species_type = positron
+beam2.injection_style = gaussian_beam
+beam2.x_rms = sigmax
+beam2.y_rms = sigmay
+beam2.z_rms = sigmaz
+beam2.x_m = 0
+beam2.y_m = 0
+beam2.z_m = -muz
+beam2.npart = nmacropart
+beam2.q_tot = beam_charge
+beam2.z_cut = 4
+beam2.momentum_distribution_type = gaussian
+beam2.uz_m = -beam_gamma
+beam2.uy_m = 0.0
+beam2.ux_m = 0.0
+beam2.ux_th = 0
+beam2.uy_th = 0
+beam2.uz_th = 0.03*beam_gamma
+beam2.do_not_deposit = 1
+
+#################################
+######### DIAGNOSTICS ###########
+#################################
+# FULL
+diagnostics.diags_names = diag1
+
+diag1.intervals = 1
+diag1.diag_type = Full
+diag1.write_species = 1
+diag1.fields_to_plot = rho_beam1 rho_beam2
+diag1.dump_last_timestep = 1
+diag1.species = beam1 beam2
+
+# REDUCED
+warpx.reduced_diags_names = DifferentialLuminosity_beam1_beam2
+
+DifferentialLuminosity_beam1_beam2.type = DifferentialLuminosity
+DifferentialLuminosity_beam1_beam2.intervals = 5
+DifferentialLuminosity_beam1_beam2.species = beam1 beam2
+DifferentialLuminosity_beam1_beam2.bin_number = 128
+DifferentialLuminosity_beam1_beam2.bin_max = 2.1*beam_energy
+DifferentialLuminosity_beam1_beam2.bin_min = 1.9*beam_energy
diff --git a/Regression/Checksum/benchmarks_json/diff_lumi_diag.json b/Regression/Checksum/benchmarks_json/diff_lumi_diag.json
@@ -0,0 +1,24 @@
+{
+  "lev=0": {
+    "rho_beam1": 656097730.398061,
+    "rho_beam2": 656073566.2600528
+  },
+  "beam2": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 1.3357156607472946e-11,
+    "particle_position_x": 0.07980004798883583,
+    "particle_position_y": 0.0015964156825818534,
+    "particle_position_z": 240184.81987149065,
+    "particle_weight": 11997240000.0
+  },
+  "beam1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 1.3357730733347833e-11,
+    "particle_position_x": 0.0796200211671279,
+    "particle_position_y": 0.001592794429510446,
+    "particle_position_z": 239913.37896451252,
+    "particle_weight": 11997780000.0
+  }
+}
diff --git a/Regression/WarpX-tests.ini b/Regression/WarpX-tests.ini
@@ -3866,3 +3866,18 @@ numprocs = 2
 useOMP = 1
 numthreads = 1
 analysisRoutine = Examples/Tests/openbc_poisson_solver/analysis.py
+
+[diff_lumi_diag]
+buildDir = .
+inputFile = Examples/Tests/diff_lumi_diag/inputs
+runtime_params = warpx.abort_on_warning_threshold = medium
+dim = 3
+addToCompileString = USE_OPENPMD=TRUE USE_FFT=TRUE
+cmakeSetupOpts = -DWarpX_DIMS=3 -DWarpX_FFT=ON -DWarpX_OPENPMD=ON
+restartTest = 0
+useMPI = 1
+numprocs = 2
+useOMP = 1
+numthreads = 1
+analysisRoutine = Examples/Tests/diff_lumi_diag/analysis.py
+aux1File = Tools/PostProcessing/read_raw_data.py
diff --git a/Source/Diagnostics/ReducedDiags/CMakeLists.txt b/Source/Diagnostics/ReducedDiags/CMakeLists.txt
@@ -4,6 +4,7 @@ foreach(D IN LISTS WarpX_DIMS)
       PRIVATE
         BeamRelevant.cpp
         ColliderRelevant.cpp
+        DifferentialLuminosity.cpp
         FieldEnergy.cpp
         FieldProbe.cpp
         FieldProbeParticleContainer.cpp

diff --git a/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.H b/Source/Diagnostics/ReducedDiags/DifferentialLuminosity.H
@@ -0,0 +1,62 @@
+/* Copyright 2023 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Arianna Formenti
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY_H_
+
+#include "ReducedDiags.H"
+#include <AMReX_GpuContainers.H>
+
+#include <map>
+#include <string>
+#include <vector>
+
+/**
+ *  This class contains the differential luminosity diagnostics.
+ */
+class DifferentialLuminosity : public ReducedDiags
+{
+public:
+
+    /**
+     * constructor
+     * @param[in] rd_name reduced diags names
+     */
+    DifferentialLuminosity(const std::string& rd_name);
+
+    /// name of the two colliding species
+    std::vector<std::string> m_beam_name;
+
+    /// number of bins
+    int m_bin_num;
+
+    /// max and min bin values
+    amrex::Real m_bin_max;
+    amrex::Real m_bin_min;
+
+    /// bin size
+    amrex::Real m_bin_size;
+
+    void ComputeDiags(int step) final;
+
+private:
+    /// auxiliary structure to store headers and indices of the reduced diagnostics
+    struct aux_header_index
+    {
+        std::string header;
+        int idx;
+    };
+
+    /// map to store header texts and indices of the reduced diagnostics
+    std::map<std::string, aux_header_index> m_headers_indices;
+
+    // Array in which to accumulate the luminosity across timesteps
+    amrex::Gpu::DeviceVector< amrex::Real > d_data;
+};
+
+#endif  // WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY_H_