WarpX
BinaryCollisionUtils.H
Go to the documentation of this file.
1 /* Copyright 2021 Neil Zaim
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_BINARY_COLLISION_UTILS_H_
9 #define WARPX_BINARY_COLLISION_UTILS_H_
10 
11 #include <string>
12 
14 
15 #include <AMReX_Math.H>
16 
22  DSMC,
24  Undefined };
25 
26 enum struct NuclearFusionType {
32  Undefined };
33 
34 namespace BinaryCollisionUtils{
35 
36  NuclearFusionType get_nuclear_fusion_type (const std::string& collision_name,
37  MultiParticleContainer const * mypc);
38 
39  CollisionType get_collision_type (const std::string& collision_name,
40  MultiParticleContainer const * mypc);
41 
43 
51  const amrex::ParticleReal& u1x, const amrex::ParticleReal& u1y,
52  const amrex::ParticleReal& u1z, const amrex::ParticleReal& u2x,
53  const amrex::ParticleReal& u2y, const amrex::ParticleReal& u2z,
54  const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
55  amrex::ParticleReal& E_kin_COM, amrex::ParticleReal& v_rel_COM,
56  amrex::ParticleReal& lab_to_COM_lorentz_factor )
57  {
58  // General notations in this function:
59  // x_sq denotes the square of x
60  // x_star denotes the value of x in the center of mass frame
61 
62  using namespace amrex::literals;
63  using namespace amrex::Math;
64 
65  constexpr auto one_pr = amrex::ParticleReal(1.);
66  constexpr auto inv_four_pr = amrex::ParticleReal(1./4.);
67  constexpr double c_sq = PhysConst::c * PhysConst::c;
68  constexpr double inv_csq = 1.0 / c_sq;
69 
70  const amrex::ParticleReal m1_sq = m1*m1;
71  const amrex::ParticleReal m2_sq = m2*m2;
72 
73  // Compute Lorentz factor gamma in the lab frame
74  const double g1 = std::sqrt( 1.0 + static_cast<double>(u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq );
75  const double g2 = std::sqrt( 1.0 + static_cast<double>(u2x*u2x+u2y*u2y+u2z*u2z)*inv_csq );
76 
77  // Compute momenta
78  const amrex::ParticleReal p1x = u1x * m1;
79  const amrex::ParticleReal p1y = u1y * m1;
80  const amrex::ParticleReal p1z = u1z * m1;
81  const amrex::ParticleReal p2x = u2x * m2;
82  const amrex::ParticleReal p2y = u2y * m2;
83  const amrex::ParticleReal p2z = u2z * m2;
84  // Square norm of the total (sum between the two particles) momenta in the lab frame
85  const auto p_total_sq = static_cast<double>(
86  powi<2>(p1x + p2x) + powi<2>(p1y + p2y) + powi<2>(p1z + p2z)
87  );
88 
89  // Total energy in the lab frame
90  // Note the use of `double` for energy since this calculation is
91  // prone to error with single precision.
92  const auto m1_dbl = static_cast<double>(m1);
93  const auto m2_dbl = static_cast<double>(m2);
94  const double E_lab = (m1_dbl * g1 + m2_dbl * g2) * c_sq;
95  // Total energy squared in the center of mass frame, calculated using the Lorentz invariance
96  // of the four-momentum norm
97  const double E_star_sq = E_lab*E_lab - c_sq*p_total_sq;
98 
99  // Kinetic energy in the center of mass frame
100  const double E_star = std::sqrt(E_star_sq);
101  E_kin_COM = static_cast<amrex::ParticleReal>(E_star - (m1_dbl + m2_dbl)*c_sq);
102 
103  // Square of the norm of the momentum of one of the particles in the center of mass frame
104  // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle
105  // The expression below is specifically written in a form that avoids returning
106  // small negative numbers due to machine precision errors, for low-energy particles
107  const auto E_ratio = static_cast<amrex::ParticleReal>(E_star/((m1 + m2)*c_sq));
108  const auto p_star_sq = static_cast<amrex::ParticleReal>(
109  m1*m2*c_sq * ( powi<2>(E_ratio) - one_pr )
110  + powi<2>(m1 - m2)*c_sq*inv_four_pr * powi<2>( E_ratio - 1._prt/E_ratio)
111  );
112 
113  // Lorentz factors in the center of mass frame
114  const auto g1_star = std::sqrt(one_pr + p_star_sq / static_cast<amrex::ParticleReal>(m1_sq*c_sq));
115  const auto g2_star = std::sqrt(one_pr + p_star_sq / static_cast<amrex::ParticleReal>(m2_sq*c_sq));
116 
117  // relative velocity in the center of mass frame
118  v_rel_COM = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) + one_pr/(m2*g2_star));
119 
120  // Cross sections and relative velocity are computed in the center of mass frame.
121  // On the other hand, the particle densities (weight over volume) in the lab frame are used.
122  // To take this disrepancy into account, it is needed to multiply the
123  // collision probability by the ratio between the Lorentz factors in the
124  // COM frame and the Lorentz factors in the lab frame (see
125  // Perez et al., Phys.Plasmas.19.083104 (2012)). The correction factor
126  // is calculated here.
127  lab_to_COM_lorentz_factor = g1_star*g2_star/static_cast<amrex::ParticleReal>(g1*g2);
128  }
129 }
130 
131 #endif // WARPX_BINARY_COLLISION_UTILS_H_
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
CollisionType
Definition: BinaryCollisionUtils.H:17
@ ProtonBoronToAlphasFusion
@ DeuteriumDeuteriumToProtonTritiumFusion
@ DeuteriumDeuteriumToNeutronHeliumFusion
@ DeuteriumTritiumToNeutronHeliumFusion
@ DeuteriumHeliumToProtonHeliumFusion
NuclearFusionType
Definition: BinaryCollisionUtils.H:26
@ DeuteriumDeuteriumToProtonTritium
@ DeuteriumDeuteriumToNeutronHelium
Definition: MultiParticleContainer.H:66
Definition: BinaryCollisionUtils.cpp:18
CollisionType get_collision_type(const std::string &collision_name, MultiParticleContainer const *const mypc)
Definition: BinaryCollisionUtils.cpp:20
CollisionType nuclear_fusion_type_to_collision_type(const NuclearFusionType fusion_type)
Definition: BinaryCollisionUtils.cpp:125
NuclearFusionType get_nuclear_fusion_type(const std::string &collision_name, MultiParticleContainer const *const mypc)
Definition: BinaryCollisionUtils.cpp:40
AMREX_GPU_HOST_DEVICE AMREX_INLINE void get_collision_parameters(const amrex::ParticleReal &u1x, const amrex::ParticleReal &u1y, const amrex::ParticleReal &u1z, const amrex::ParticleReal &u2x, const amrex::ParticleReal &u2y, const amrex::ParticleReal &u2z, const amrex::ParticleReal &m1, const amrex::ParticleReal &m2, amrex::ParticleReal &E_kin_COM, amrex::ParticleReal &v_rel_COM, amrex::ParticleReal &lab_to_COM_lorentz_factor)
Return (relativistic) collision energy, collision speed and Lorentz factor for transforming between t...
Definition: BinaryCollisionUtils.H:50
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44