Generalize differential luminosity for photons

Docs/source/usage/parameters.rst

@@ -3466,14 +3466,15 @@ Reduced Diagnostics
             \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}^*)`
+        where :math:`f_i` is the distribution function of species :math:`i` and
+        :math:`\mathcal{E}^*(\boldsymbol{p}_1, \boldsymbol{p}_2) = \sqrt{m_1^2c^4 + m_2^2c^4 + 2 c^2{p_1}^\mu {p_2}_\mu}`
+        is the energy in the center-of-mass frame, where :math:`p^\mu = (\sqrt{m^2 c^2 + \boldsymbol{p}^2}, \boldsymbol{p})`
+        represents the 4-momentum. 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{eV}^{-1}`. For collider-relevant WarpX simulations
+        The differential luminosity is given in units of :math:`\text{m}^{-2}.\text{eV}^{-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{eV}^{-1}`
         can be obtained by multiplying the above differential luminosity by the expected repetition rate of the beams.
 

Examples/Tests/diff_lumi_diag/CMakeLists.txt

@@ -2,10 +2,20 @@
 #
 
 add_warpx_test(
-    test_3d_diff_lumi_diag  # name
+    test_3d_diff_lumi_diag_leptons  # name
     3  # dims
     2  # nprocs
-    inputs_test_3d_diff_lumi_diag  # inputs
+    inputs_test_3d_diff_lumi_diag_leptons  # inputs
+    analysis.py  # analysis
+    diags/diag1000080  # output
+    OFF  # dependency
+)
+
+add_warpx_test(
+    test_3d_diff_lumi_diag_photons  # name
+    3  # dims
+    2  # nprocs
+    inputs_test_3d_diff_lumi_diag_photons  # inputs
     analysis.py  # analysis
     diags/diag1000080  # output
     OFF  # dependency

Examples/Tests/diff_lumi_diag/analysis.py

@@ -37,16 +37,28 @@
     * np.exp(-((E_bin - 2 * E_beam) ** 2) / (2 * sigma_E**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
+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
+
 # 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
 
 # compare checksums
 evaluate_checksum(
-    test_name=os.path.split(os.getcwd())[1],
+    test_name=test_name,
     output_file=sys.argv[1],
     rtol=1e-2,
 )

Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag → Examples/Tests/diff_lumi_diag/inputs_base_3d

@@ -6,12 +6,11 @@ my_constants.mc2_eV = m_e*clight*clight/q_e
 # BEAMS
 my_constants.beam_energy_eV = 125.e9
 my_constants.beam_gamma = beam_energy_eV/(mc2_eV)
-my_constants.beam_charge = 1.2e10*q_e
+my_constants.beam_N = 1.2e10
 my_constants.sigmax = 500e-9
 my_constants.sigmay = 10e-9
 my_constants.sigmaz = 300e-3
-my_constants.muz = -4*sigmaz
-my_constants.nmacropart = 2e5
+my_constants.muz = 4*sigmaz
 
 # BOX
 my_constants.Lx = 8*sigmax
@@ -62,17 +61,6 @@ warpx.poisson_solver = fft
 #################################
 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
@@ -82,17 +70,6 @@ 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
@@ -108,7 +85,7 @@ beam2.do_not_deposit = 1
 # FULL
 diagnostics.diags_names = diag1
 
-diag1.intervals = 1
+diag1.intervals = 80
 diag1.diag_type = Full
 diag1.write_species = 1
 diag1.fields_to_plot = rho_beam1 rho_beam2

Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_leptons

@@ -0,0 +1,31 @@
+# base input parameters
+FILE = inputs_base_3d
+
+# Test with electrons/positrons: use gaussian beam distribution
+# by providing the total charge (q_tot)
+
+my_constants.nmacropart = 2e5
+
+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_N*q_e
+beam1.z_cut = 4
+
+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_N*q_e
+beam2.z_cut = 4

Examples/Tests/diff_lumi_diag/inputs_test_3d_diff_lumi_diag_photons

@@ -0,0 +1,28 @@
+# base input parameters
+FILE = inputs_base_3d
+
+# Test with electrons/positrons: use parse_density_function
+
+beam1.species_type = electron
+beam1.injection_style = "NUniformPerCell"
+beam1.num_particles_per_cell_each_dim = 1 1 1
+beam1.profile = parse_density_function
+beam1.density_function(x,y,z) = "beam_N/(sqrt(2*pi)*2*pi*sigmax*sigmay*sigmaz)*exp(-x*x/(2*sigmax*sigmax)-y*y/(2*sigmay*sigmay)-(z+muz)*(z+muz)/(2*sigmaz*sigmaz))"
+beam1.xmin = -4*sigmax
+beam1.xmax =  4*sigmax
+beam1.ymin = -4*sigmay
+beam1.ymax =  4*sigmay
+beam1.zmin =-muz-4*sigmaz
+beam1.zmax =-muz+4*sigmaz
+
+beam2.species_type = positron
+beam2.injection_style = "NUniformPerCell"
+beam2.num_particles_per_cell_each_dim = 1 1 1
+beam2.profile = parse_density_function
+beam2.xmin = -4*sigmax
+beam2.xmax =  4*sigmax
+beam2.ymin = -4*sigmay
+beam2.ymax =  4*sigmay
+beam2.zmin = muz-4*sigmaz
+beam2.zmax = muz+4*sigmaz
+beam2.density_function(x,y,z) = "beam_N/(sqrt(2*pi)*2*pi*sigmax*sigmay*sigmaz)*exp(-x*x/(2*sigmax*sigmax)-y*y/(2*sigmay*sigmay)-(z-muz)*(z-muz)/(2*sigmaz*sigmaz))"

Regression/Checksum/benchmarks_json/test_3d_diff_lumi_diag_photons.json

@@ -0,0 +1,24 @@
+{
+  "lev=0": {
+    "rho_beam1": 656097367.2335038,
+    "rho_beam2": 656097367.2335038
+  },
+  "beam1": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 1.7512476113279403e-11,
+    "particle_position_x": 0.2621440000000001,
+    "particle_position_y": 0.005242880000000001,
+    "particle_position_z": 314572.79999473685,
+    "particle_weight": 11997744756.90957
+  },
+  "beam2": {
+    "particle_momentum_x": 0.0,
+    "particle_momentum_y": 0.0,
+    "particle_momentum_z": 1.7513431895752007e-11,
+    "particle_position_x": 0.2621440000000001,
+    "particle_position_y": 0.005242880000000001,
+    "particle_position_z": 314572.79999472946,
+    "particle_weight": 11997744756.909573
+  }
+}

Source/Diagnostics/ReducedDiags/DifferentialLuminosity.cpp

@@ -132,9 +132,8 @@ void DifferentialLuminosity::ComputeDiags (int step)
     // Since this diagnostic *accumulates* the luminosity in the
     // array d_data, we add contributions at *each timestep*, but
     // we only write the data to file at intervals specified by the user.
-
-    const Real c2_over_qe = PhysConst::c*PhysConst::c/PhysConst::q_e;
-    const Real inv_c2 = 1._rt/(PhysConst::c*PhysConst::c);
+    const Real c_sq = PhysConst::c*PhysConst::c;
+    const Real c_over_qe = PhysConst::c/PhysConst::q_e;
 
     // get a reference to WarpX instance
     auto& warpx = WarpX::GetInstance();
@@ -187,6 +186,7 @@ void DifferentialLuminosity::ComputeDiags (int step)
             amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
             amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy]; // v*gamma=p/m
             amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
+            bool const species1_is_photon = species_1.AmIA<PhysicalSpecies::photon>();
 
             const auto soa_2 = ptile_2.getParticleTileData();
             index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr();
@@ -196,6 +196,7 @@ void DifferentialLuminosity::ComputeDiags (int step)
             amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
             amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
             amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
+            bool const species2_is_photon = species_2.AmIA<PhysicalSpecies::photon>();
 
             // Extract low-level data
             auto const n_cells = static_cast<int>(bins_1.numBins());
@@ -218,34 +219,59 @@ void DifferentialLuminosity::ComputeDiags (int step)
                         index_type const j_1 = indices_1[i_1];
                         index_type const j_2 = indices_2[i_2];
 
-                        Real const u1_square =  u1x[j_1]*u1x[j_1] + u1y[j_1]*u1y[j_1] + u1z[j_1]*u1z[j_1];
-                        Real const gamma1 = std::sqrt(1._rt + u1_square*inv_c2);
-                        Real const u2_square = u2x[j_2]*u2x[j_2] + u2y[j_2]*u2y[j_2] + u2z[j_2]*u2z[j_2];
-                        Real const gamma2 = std::sqrt(1._rt + u2_square*inv_c2);
-                        Real const u1_dot_u2 = u1x[j_1]*u2x[j_2] + u1y[j_1]*u2y[j_2] + u1z[j_1]*u2z[j_2];
+                        Real p1t=0, p1x=0, p1y=0, p1z=0; // components of 4-momentum of particle 1
+                        Real const u1_sq =  u1x[j_1]*u1x[j_1] + u1y[j_1]*u1y[j_1] + u1z[j_1]*u1z[j_1];
+                        if (species1_is_photon) {
+                            // photon case (momentum is normalized by m_e in WarpX)
+                            p1t = PhysConst::m_e*std::sqrt( u1_sq );
+                            p1x = PhysConst::m_e*u1x[j_1];
+                            p1y = PhysConst::m_e*u1y[j_1];
+                            p1z = PhysConst::m_e*u1z[j_1];
+                        } else {
+                            p1t = m1*std::sqrt( c_sq + u1_sq );
+                            p1x = m1*u1x[j_1];
+                            p1y = m1*u1y[j_1];
+                            p1z = m1*u1z[j_1];
+                        }
+
+                        Real p2t=0, p2x=0, p2y=0, p2z=0; // components of 4-momentum of particle 2
+                        Real const u2_sq =  u2x[j_2]*u2x[j_2] + u2y[j_2]*u2y[j_2] + u2z[j_2]*u2z[j_2];
+                        if (species2_is_photon) {
+                            // photon case (momentum is normalized by m_e in WarpX)
+                            p2t = PhysConst::m_e*std::sqrt(u2_sq);
+                            p2x = PhysConst::m_e*u2x[j_2];
+                            p2y = PhysConst::m_e*u2y[j_2];
+                            p2z = PhysConst::m_e*u2z[j_2];
+                        } else {
+                            p2t = m2*std::sqrt( c_sq + u2_sq );
+                            p2x = m2*u2x[j_2];
+                            p2y = m2*u2y[j_2];
+                            p2z = m2*u2z[j_2];
+                        }
 
                         // center of mass energy in eV
-                        Real const E_com = c2_over_qe * std::sqrt(m1*m1 + m2*m2 + 2*m1*m2* (gamma1*gamma2 - u1_dot_u2*inv_c2));
+                        Real const E_com = c_over_qe * std::sqrt(m1*m1*c_sq + m2*m2*c_sq + 2*(p1t*p2t - p1x*p2x - p1y*p2y - p1z*p2z));
 
                         // determine particle bin
                         int const bin = int(Math::floor((E_com-bin_min)/bin_size));
 
                         if ( bin<0 || bin>=num_bins ) { continue; } // discard if out-of-range
 
-                        Real const v1_minus_v2_x = u1x[j_1]/gamma1 - u2x[j_2]/gamma2;
-                        Real const v1_minus_v2_y = u1y[j_1]/gamma1 - u2y[j_2]/gamma2;
-                        Real const v1_minus_v2_z = u1z[j_1]/gamma1 - u2z[j_2]/gamma2;
-                        Real const v1_minus_v2_square = v1_minus_v2_x*v1_minus_v2_x + v1_minus_v2_y*v1_minus_v2_y + v1_minus_v2_z*v1_minus_v2_z;
+                        Real const inv_p1t = 1.0_rt/p1t;
+                        Real const inv_p2t = 1.0_rt/p2t;
 
-                        Real const u1_cross_u2_x = u1y[j_1]*u2z[j_2] - u1z[j_1]*u2y[j_2];
-                        Real const u1_cross_u2_y = u1z[j_1]*u2x[j_2] - u1x[j_1]*u2z[j_2];
-                        Real const u1_cross_u2_z = u1x[j_1]*u2y[j_2] - u1y[j_1]*u2x[j_2];
+                        Real const beta1_sq = (p1x*p1x + p1y*p1y + p1z*p1z) * inv_p1t*inv_p1t;
+                        Real const beta2_sq = (p2x*p2x + p2y*p2y + p2z*p2z) * inv_p2t*inv_p2t;
+                        Real const beta1_dot_beta2 = (p1x*p2x + p1y*p2y + p1z*p2z) * inv_p1t*inv_p2t;
 
-                        Real const v1_cross_v2_square = (u1_cross_u2_x*u1_cross_u2_x + u1_cross_u2_y*u1_cross_u2_y + u1_cross_u2_z*u1_cross_u2_z) / (gamma1*gamma1*gamma2*gamma2);
+                        // Here we use the fact that:
+                        // (v1 - v2)^2 = v1^2 + v2^2 - 2 v1.v2
+                        // and (v1 x v2)^2 = v1^2 v2^2 - (v1.v2)^2
+                        // we also use beta=v/c instead of v
 
-                        Real const radicand = v1_minus_v2_square - v1_cross_v2_square * inv_c2;
+                        Real const radicand = beta1_sq + beta2_sq - 2*beta1_dot_beta2 - beta1_sq*beta2_sq + beta1_dot_beta2*beta1_dot_beta2;
 
-                        Real const dL_dEcom = std::sqrt( radicand ) * w1[j_1] * w2[j_2] / dV / bin_size * dt; // m^-2 eV^-1
+                        Real const dL_dEcom = PhysConst::c * std::sqrt( radicand ) * w1[j_1] * w2[j_2] / dV / bin_size * dt; // m^-2 eV^-1
 
                         amrex::HostDevice::Atomic::Add(&dptr_data[bin], dL_dEcom);