Add reduced diagnostic: 2d differential luminosity (#5545)

Adds a luminosity diagnostic differentiated in the energies of two
colliding species, called `DifferentialLuminosity2D`.
It is defined as follows:
```math
\begin{align*}
\frac{d^2\mathcal{L}}{dE_1 dE_2}(E_1, E_2, t)  = \int_0^t dt'\int d\boldsymbol{x}\, & \int d\boldsymbol{p}_1 \int 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(E_1 - E_1(\boldsymbol{p}_1) \delta( E_2 - E_2(\boldsymbol{p}_2))
\end{align*}
```
where:
*  $\boldsymbol{p}_i$ is the momentum of a particle of species $i$
* $E_i$ is the energy of a particle of species $i$, $E_i
(\boldsymbol{p}_i) = \sqrt{m_1^2c^4 + c^2 |\boldsymbol{p}_i|^2}$
* $f_i$ is the distribution function of species $i$, normalized such
that $\int \int f(\boldsymbol{x} \boldsymbol{p}, t )d\boldsymbol{x}
d\boldsymbol{p} = N$, the number of particles in species $i$ at time $t$

The 2D differential luminosity is given in units of $\text{m}^{-2} \
\text{eV}^{-2}$.

The user must specify the minimum, maximum, and number of bins to
discretize the $E_1$ and $E_2$ axes.
The computation of this diagnostic is similar to that of
`ParticleHistogram2D`.
 
The output is a folder containing a set of openPMD files.
The values of the diagnostic are stored in a record labeled
`d2L_dE1_dE2`, with axes `E1` and `E2`.

---------

Co-authored-by: Remi Lehe <remi.lehe@normalesup.org>

Docs/source/usage/parameters.rst

@@ -3626,6 +3626,52 @@ This shifts analysis from post-processing to runtime calculation of reduction op
         * ``<reduced_diags_name>.bin_min`` (`float`, in eV)
             The maximum value of :math:`\mathcal{E}^*` for which the differential luminosity is computed.
 
+    * ``DifferentialLuminosity2D``
+        This type computes the two-dimensional differential luminosity between two species, defined as:
+
+        .. math::
+
+            \frac{d^2\mathcal{L}}{dE_1 dE_2}(E_1, E_2, t) = \int_0^t dt'\int d\boldsymbol{x}\, \int d\boldsymbol{p}_1 \int 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(E_1 - E_1(\boldsymbol{p}_1)) \delta(E_2 - E_2(\boldsymbol{p}_2))
+
+        where :math:`f_i` is the distribution function of species :math:`i`
+        (normalized such that :math:`\int \int f(\boldsymbol{x} \boldsymbol{p}, t )d\boldsymbol{x} d\boldsymbol{p} = N`
+        is the number of particles in species :math:`i` at time :math:`t`),
+        :math:`\boldsymbol{p}_i` and :math:`E_i (\boldsymbol{p}_i) = \sqrt{m_1^2c^4 + c^2 |\boldsymbol{p}_i|^2}`
+        are, respectively, the momentum and the energy of a particle of the :math:`i`-th species.
+        The 2D differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-2}`.
+
+        * ``<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_1`` (`int` > 0)
+            The number of bins in energy :math:`E_1`
+
+        * ``<reduced_diags_name>.bin_max_1`` (`float`, in eV)
+            The minimum value of :math:`E_1` for which the 2D differential luminosity is computed.
+
+        * ``<reduced_diags_name>.bin_min_1`` (`float`, in eV)
+            The maximum value of :math:`E_2` for which the 2D differential luminosity is compute
+
+        * ``<reduced_diags_name>.bin_number_2`` (`int` > 0)
+            The number of bins in energy :math:`E_2`
+
+        * ``<reduced_diags_name>.bin_max_2`` (`float`, in eV)
+            The minimum value of :math:`E_2` for which the 2D differential luminosity is computed.
+
+        * ``<reduced_diags_name>.bin_min_2`` (`float`, in eV)
+            The minimum value of :math:`E_2` for which the 2D differential luminosity is computed.
+
+        * ``<reduced_diags_name>.file_min_digits`` (`int`) optional (default `6`)
+            The minimum number of digits used for the iteration number appended to the diagnostic file names.
+
+        The output is a ``<reduced_diags_name>`` folder containing a set of openPMD files.
+        The values of the diagnostic are stored in a record labeled `d2L_dE1_dE2`.
+        An example input file and a loading python script of
+        using the DifferentialLuminosity2D reduced diagnostics
+        are given in ``Examples/Tests/diff_lumi_diag/``.
+
     * ``Timestep``
         This type outputs the simulation's physical timestep (in seconds) at each mesh refinement level.
 

Examples/Tests/diff_lumi_diag/CMakeLists.txt

@@ -1,6 +1,6 @@
 # Add tests (alphabetical order) ##############################################
 #
-
+if(WarpX_FFT)
 add_warpx_test(
     test_3d_diff_lumi_diag_leptons  # name
     3  # dims
@@ -10,7 +10,9 @@ add_warpx_test(
     "analysis_default_regression.py --path diags/diag1000080 --rtol 1e-2"  # checksum
     OFF  # dependency
 )
+endif()
 
+if(WarpX_FFT)
 add_warpx_test(
     test_3d_diff_lumi_diag_photons  # name
     3  # dims
@@ -20,3 +22,4 @@ add_warpx_test(
     "analysis_default_regression.py --path diags/diag1000080 --rtol 1e-2"  # checksum
     OFF  # dependency
 )
+endif()

Examples/Tests/diff_lumi_diag/analysis.py

@@ -5,15 +5,20 @@
 # In that case, the differential luminosity can be calculated analytically.
 
 import os
+import re
 
 import numpy as np
-from read_raw_data import read_reduced_diags_histogram
+from openpmd_viewer import OpenPMDTimeSeries
 
-# 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
+# Extract the 1D differential luminosity from the file
+filename = "./diags/reducedfiles/DifferentialLuminosity_beam1_beam2.txt"
+with open(filename) as f:
+    # First line: header, contains the energies
+    line = f.readline()
+    E_bin = np.array(list(map(float, re.findall("=(.*?)\(", line))))
+data = np.loadtxt(filename)
+dE_bin = E_bin[1] - E_bin[0]
+dL_dE_sim = data[-1, 2:]  # Differential luminosity at the end of the simulation
 
 # Beam parameters
 N = 1.2e10
@@ -33,21 +38,47 @@
     * np.exp(-((E_bin - 2 * E_beam) ** 2) / (2 * sigma_E**2))
 )
 
+# Extract the 2D differential luminosity from the file
+series = OpenPMDTimeSeries("./diags/reducedfiles/DifferentialLuminosity2d_beam1_beam2/")
+d2L_dE1_dE2_sim, info = series.get_field("d2L_dE1_dE2", iteration=80)
+
+# Compute the analytical 2D differential luminosity for 2 Gaussian beams
+assert info.axes[0] == "E2"
+assert info.axes[1] == "E1"
+E2, E1 = np.meshgrid(info.E2, info.E1, indexing="ij")
+d2L_dE1_dE2_th = (
+    N**2
+    / (2 * (2 * np.pi) ** 2 * sigma_x * sigma_y * sigma_E1 * sigma_E2)
+    * np.exp(
+        -((E1 - E_beam) ** 2) / (2 * sigma_E1**2)
+        - (E2 - E_beam) ** 2 / (2 * sigma_E2**2)
+    )
+)
+
 # Extract test name from path
 test_name = os.path.split(os.getcwd())[1]
 print("test_name", test_name)
 
 # Pick tolerance
 if "leptons" in test_name:
-    tol = 1e-2
+    tol1 = 0.02
+    tol2 = 0.04
 elif "photons" in test_name:
     # In the photons case, the particles are
     # initialized from a density distribution ;
     # tolerance is larger due to lower particle statistics
-    tol = 6e-2
+    tol1 = 0.021
+    tol2 = 0.06
+
+# Check that the 1D diagnostic and analytical result match
+error1 = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max()
+print("Relative error: ", error1)
+print("Tolerance: ", tol1)
+
+# Check that the 2D and 1D diagnostics match
+error2 = abs(d2L_dE1_dE2_sim - d2L_dE1_dE2_th).max() / abs(d2L_dE1_dE2_th).max()
+print("Relative error: ", error2)
+print("Tolerance: ", tol2)
 
-# Check that the simulation result and analytical result match
-error = abs(dL_dE_sim - dL_dE_th).max() / abs(dL_dE_th).max()
-print("Relative error: ", error)
-print("Tolerance: ", tol)
-assert error < tol
+assert error1 < tol1
+assert error2 < tol2

Examples/Tests/diff_lumi_diag/inputs_base_3d

@@ -28,6 +28,7 @@ my_constants.dt = sigmaz/clight/10.
 #################################
 ####### GENERAL PARAMETERS ######
 #################################
+
 stop_time = T
 amr.n_cell = nx ny nz
 amr.max_grid_size = 128
@@ -93,11 +94,21 @@ diag1.dump_last_timestep = 1
 diag1.species = beam1 beam2
 
 # REDUCED
-warpx.reduced_diags_names = DifferentialLuminosity_beam1_beam2
+warpx.reduced_diags_names = DifferentialLuminosity_beam1_beam2 DifferentialLuminosity2d_beam1_beam2
 
 DifferentialLuminosity_beam1_beam2.type = DifferentialLuminosity
-DifferentialLuminosity_beam1_beam2.intervals = 5
+DifferentialLuminosity_beam1_beam2.intervals = 80
 DifferentialLuminosity_beam1_beam2.species = beam1 beam2
 DifferentialLuminosity_beam1_beam2.bin_number = 128
 DifferentialLuminosity_beam1_beam2.bin_max = 2.1*beam_energy_eV
-DifferentialLuminosity_beam1_beam2.bin_min = 1.9*beam_energy_eV
+DifferentialLuminosity_beam1_beam2.bin_min = 0
+
+DifferentialLuminosity2d_beam1_beam2.type = DifferentialLuminosity2D
+DifferentialLuminosity2d_beam1_beam2.intervals = 80
+DifferentialLuminosity2d_beam1_beam2.species = beam1 beam2
+DifferentialLuminosity2d_beam1_beam2.bin_number_1 = 128
+DifferentialLuminosity2d_beam1_beam2.bin_max_1 = 1.45*beam_energy_eV
+DifferentialLuminosity2d_beam1_beam2.bin_min_1 = 0
+DifferentialLuminosity2d_beam1_beam2.bin_number_2 = 128
+DifferentialLuminosity2d_beam1_beam2.bin_max_2 = 1.45*beam_energy_eV
+DifferentialLuminosity2d_beam1_beam2.bin_min_2 = 0

Source/Diagnostics/ReducedDiags/CMakeLists.txt

@@ -6,6 +6,7 @@ foreach(D IN LISTS WarpX_DIMS)
         ChargeOnEB.cpp
         ColliderRelevant.cpp
         DifferentialLuminosity.cpp
+        DifferentialLuminosity2D.cpp
         FieldEnergy.cpp
         FieldMaximum.cpp
         FieldMomentum.cpp

Source/Diagnostics/ReducedDiags/DifferentialLuminosity2D.H

@@ -0,0 +1,70 @@
+/* Copyright 2023 The WarpX Community
+ *
+ * This file is part of WarpX.
+ *
+ * Authors: Arianna Formenti, Remi Lehe
+ * License: BSD-3-Clause-LBNL
+ */
+
+#ifndef WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_
+#define WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_
+
+#include "ReducedDiags.H"
+#include <AMReX_GpuContainers.H>
+#include <AMReX_TableData.H>
+
+#include <map>
+#include <string>
+#include <vector>
+
+/**
+ *  This class contains the differential luminosity diagnostics.
+ */
+class DifferentialLuminosity2D : public ReducedDiags
+{
+public:
+
+    /**
+     * constructor
+     * @param[in] rd_name reduced diags names
+     */
+    DifferentialLuminosity2D(const std::string& rd_name);
+
+    /// File type
+    std::string m_openpmd_backend {"default"};
+
+    /// minimum number of digits for file suffix (file-based only supported for now) */
+    int m_file_min_digits = 6;
+
+    /// name of the two colliding species
+    std::vector<std::string> m_beam_name;
+
+    /// number of bins for the c.o.m. energy of the 2 species
+    int m_bin_num_1;
+    int m_bin_num_2;
+
+    /// max and min bin values
+    amrex::Real m_bin_max_1;
+    amrex::Real m_bin_min_1;
+    amrex::Real m_bin_max_2;
+    amrex::Real m_bin_min_2;
+
+    /// bin size
+    amrex::Real m_bin_size_1;
+    amrex::Real m_bin_size_2;
+
+    /// output data
+    amrex::TableData<amrex::Real,2> m_h_data_2D;
+
+    void ComputeDiags(int step) final;
+
+    void WriteToFile (int step) const final;
+
+private:
+
+    /// output table in which to accumulate the luminosity across timesteps
+    amrex::TableData<amrex::Real,2> m_d_data_2D;
+
+};
+
+#endif  // WARPX_DIAGNOSTICS_REDUCEDDIAGS_DIFFERENTIALLUMINOSITY2D_H_