WarpX
TwoProductFusionUtil.H
Go to the documentation of this file.
1 /* Copyright 2022 Remi Lehe
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 
8 #ifndef WARPX_TWO_PRODUCT_FUSION_UTIL_H
9 #define WARPX_TWO_PRODUCT_FUSION_UTIL_H
10 
11 #include "Utils/ParticleUtils.H"
12 #include "Utils/WarpXConst.H"
13 
14 #include <AMReX_Math.H>
15 #include <AMReX_Random.H>
16 #include <AMReX_REAL.H>
17 
18 #include <cmath>
19 #include <limits>
20 
21 namespace {
48  void TwoProductFusionComputeProductMomenta (
49  const amrex::ParticleReal& u1x_in,
50  const amrex::ParticleReal& u1y_in,
51  const amrex::ParticleReal& u1z_in,
52  const amrex::ParticleReal& m1_in,
53  const amrex::ParticleReal& u2x_in,
54  const amrex::ParticleReal& u2y_in,
55  const amrex::ParticleReal& u2z_in,
56  const amrex::ParticleReal& m2_in,
57  amrex::ParticleReal& u1x_out,
58  amrex::ParticleReal& u1y_out,
59  amrex::ParticleReal& u1z_out,
60  const amrex::ParticleReal& m1_out,
61  amrex::ParticleReal& u2x_out,
62  amrex::ParticleReal& u2y_out,
63  amrex::ParticleReal& u2z_out,
64  const amrex::ParticleReal& m2_out,
65  const amrex::ParticleReal& E_fusion,
66  const amrex::RandomEngine& engine )
67  {
68  using namespace amrex::literals;
69  using namespace amrex::Math;
70 
71  constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c;
72  constexpr amrex::ParticleReal inv_csq = 1._prt / ( c_sq );
73  // Rest energy of incident particles
74  const amrex::ParticleReal E_rest_in = (m1_in + m2_in)*c_sq;
75  // Rest energy of products
76  const amrex::ParticleReal E_rest_out = (m1_out + m2_out)*c_sq;
77 
78  // Compute Lorentz factor gamma in the lab frame
79  const amrex::ParticleReal g1_in = std::sqrt( 1._prt
80  + (u1x_in*u1x_in + u1y_in*u1y_in + u1z_in*u1z_in)*inv_csq );
81  const amrex::ParticleReal g2_in = std::sqrt( 1._prt
82  + (u2x_in*u2x_in + u2y_in*u2y_in + u2z_in*u2z_in)*inv_csq );
83 
84  // Compute momenta
85  const amrex::ParticleReal p1x_in = u1x_in * m1_in;
86  const amrex::ParticleReal p1y_in = u1y_in * m1_in;
87  const amrex::ParticleReal p1z_in = u1z_in * m1_in;
88  const amrex::ParticleReal p2x_in = u2x_in * m2_in;
89  const amrex::ParticleReal p2y_in = u2y_in * m2_in;
90  const amrex::ParticleReal p2z_in = u2z_in * m2_in;
91  // Square norm of the total (sum between the two particles) momenta in the lab frame
92  const amrex::ParticleReal p_total_sq = powi<2>(p1x_in+p2x_in) +
93  powi<2>(p1y_in+p2y_in) +
94  powi<2>(p1z_in+p2z_in);
95 
96  // Total energy of incident macroparticles in the lab frame
97  const amrex::ParticleReal E_lab = (m1_in * g1_in + m2_in * g2_in) * c_sq;
98  // Total energy squared of proton+boron in the center of mass frame, calculated using the
99  // Lorentz invariance of the four-momentum norm
100  const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq;
101  // Total energy squared of the products in the center of mass frame
102  // In principle, the term - E_rest_in + E_rest_out + E_fusion is not needed and equal to
103  // zero (i.e. the energy liberated during fusion is equal to the mass difference). However,
104  // due to possible inconsistencies in how the mass is defined in the code, it is
105  // probably more robust to subtract the rest masses and to add the fusion energy to the
106  // total kinetic energy.
107  const amrex::ParticleReal E_star_f_sq = powi<2>(std::sqrt(E_star_sq)
108  - E_rest_in + E_rest_out + E_fusion);
109 
110  // Square of the norm of the momentum of the products in the center of mass frame
111  // Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle
112  // The expression below is specifically written in a form that avoids returning
113  // small negative numbers due to machine precision errors, for low-energy particles
114  const amrex::ParticleReal E_ratio = std::sqrt(E_star_f_sq)/((m1_out + m2_out)*c_sq);
115  const amrex::ParticleReal p_star_f_sq = m1_out*m2_out*c_sq * ( powi<2>(E_ratio) - 1._prt )
116  + powi<2>(m1_out - m2_out)*c_sq*0.25_prt * powi<2>( E_ratio - 1._prt/E_ratio );
117 
118  // Compute momentum of first product in the center of mass frame, assuming isotropic
119  // distribution
120  amrex::ParticleReal px_star, py_star, pz_star;
121  ParticleUtils::RandomizeVelocity(px_star, py_star, pz_star, std::sqrt(p_star_f_sq),
122  engine);
123 
124  // Next step is to convert momenta to lab frame
125  amrex::ParticleReal p1x_out, p1y_out, p1z_out;
126  // Preliminary calculation: compute center of mass velocity vc
127  const amrex::ParticleReal mass_g = m1_in * g1_in + m2_in * g2_in;
128  const amrex::ParticleReal vcx = (p1x_in+p2x_in) / mass_g;
129  const amrex::ParticleReal vcy = (p1y_in+p2y_in) / mass_g;
130  const amrex::ParticleReal vcz = (p1z_in+p2z_in) / mass_g;
131  const amrex::ParticleReal vc_sq = vcx*vcx + vcy*vcy + vcz*vcz;
132 
133  // Convert momentum of first product to lab frame, using equation (13) of F. Perez et al.,
134  // Phys.Plasmas.19.083104 (2012)
135  if ( vc_sq > std::numeric_limits<amrex::ParticleReal>::min() )
136  {
137  const amrex::ParticleReal gc = 1._prt / std::sqrt( 1._prt - vc_sq*inv_csq );
138  const amrex::ParticleReal g_star = std::sqrt(1._prt + p_star_f_sq / (m1_out*m1_out*c_sq));
139  const amrex::ParticleReal vcDps = vcx*px_star + vcy*py_star + vcz*pz_star;
140  const amrex::ParticleReal factor0 = (gc-1._prt)/vc_sq;
141  const amrex::ParticleReal factor = factor0*vcDps + m1_out*g_star*gc;
142  p1x_out = px_star + vcx * factor;
143  p1y_out = py_star + vcy * factor;
144  p1z_out = pz_star + vcz * factor;
145  }
146  else // If center of mass velocity is zero, we are already in the lab frame
147  {
148  p1x_out = px_star;
149  p1y_out = py_star;
150  p1z_out = pz_star;
151  }
152 
153  // Compute momentum of beryllium in lab frame, using total momentum conservation
154  const amrex::ParticleReal p2x_out = p1x_in + p2x_in - p1x_out;
155  const amrex::ParticleReal p2y_out = p1y_in + p2y_in - p1y_out;
156  const amrex::ParticleReal p2z_out = p1z_in + p2z_in - p1z_out;
157 
158  // Compute the momentum of the product macroparticles
159  u1x_out = p1x_out/m1_out;
160  u1y_out = p1y_out/m1_out;
161  u1z_out = p1z_out/m1_out;
162  u2x_out = p2x_out/m2_out;
163  u2y_out = p2y_out/m2_out;
164  u2z_out = p2z_out/m2_out;
165  }
166 }
167 
168 #endif // WARPX_TWO_PRODUCT_FUSION_UTIL_H
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
AMREX_GPU_HOST_DEVICE AMREX_INLINE void RandomizeVelocity(amrex::ParticleReal &ux, amrex::ParticleReal &uy, amrex::ParticleReal &uz, const amrex::ParticleReal vp, amrex::RandomEngine const &engine)
Function to perform scattering of a particle that results in a random velocity vector with given magn...
Definition: ParticleUtils.H:162
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44