WarpX
ProtonBoronFusionInitializeMomentum.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_PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
9 #define WARPX_PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
10 
11 #include "TwoProductFusionUtil.H"
13 #include "Utils/ParticleUtils.H"
14 #include "Utils/WarpXConst.H"
15 
16 #include <AMReX_DenseBins.H>
17 #include <AMReX_Random.H>
18 #include <AMReX_REAL.H>
19 
20 #include <cmath>
21 #include <limits>
22 
23 namespace {
24  // Define shortcuts for frequently-used type names
30  using index_type = typename ParticleBins::index_type;
31 
60  void ProtonBoronFusionInitializeMomentum (
61  const SoaData_type& soa_1, const SoaData_type& soa_2,
62  SoaData_type& soa_alpha,
63  const index_type& idx_1, const index_type& idx_2,
64  const index_type& idx_alpha_start,
65  const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
66  const amrex::RandomEngine& engine)
67  {
68  // General notations in this function:
69  // x_sq denotes the square of x
70  // x_Bestar denotes the value of x in the Beryllium rest frame
71 
72  using namespace amrex::literals;
73 
74  constexpr amrex::ParticleReal mev_to_joule = PhysConst::q_e*1.e6_prt;
75  // Energy produced in the fusion reaction proton + boron11 -> Beryllium8 + alpha
76  // cf. Janis book of proton-induced cross-sections (2019)
77  constexpr amrex::ParticleReal E_fusion = 8.59009_prt*mev_to_joule;
78  // Energy produced when Beryllium8 decays into two alphas
79  // cf. JEFF-3.3 radioactive decay data library (2017)
80  constexpr amrex::ParticleReal E_decay = 0.0918984_prt*mev_to_joule;
81 
82  // The constexprs ma_sq and mBe_sq underflow in single precision because we use SI units,
83  // which can cause compilation to fail or generate a warning, so we're explicitly setting
84  // them as double. Note that nuclear fusion module does not currently work with single
85  // precision anyways.
86  constexpr double m_alpha = PhysConst::m_u * 4.00260325413_prt;
87  constexpr double m_beryllium = PhysConst::m_p * 7.94748_prt;
88  constexpr double mBe_sq = m_beryllium*m_beryllium;
89  constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c;
90 
91  // Compute the resulting momenta of the alpha and beryllium particles
92  // produced in the fusion reaction
93  amrex::ParticleReal ux_alpha1 = 0.0, uy_alpha1 = 0.0, uz_alpha1 = 0.0;
94  amrex::ParticleReal ux_Be = 0.0, uy_Be = 0.0, uz_Be = 0.0;
95  TwoProductFusionComputeProductMomenta (
96  soa_1.m_rdata[PIdx::ux][idx_1],
97  soa_1.m_rdata[PIdx::uy][idx_1],
98  soa_1.m_rdata[PIdx::uz][idx_1], m1,
99  soa_2.m_rdata[PIdx::ux][idx_2],
100  soa_2.m_rdata[PIdx::uy][idx_2],
101  soa_2.m_rdata[PIdx::uz][idx_2], m2,
102  ux_alpha1, uy_alpha1, uz_alpha1, static_cast<amrex::ParticleReal>(m_alpha),
103  ux_Be, uy_Be, uz_Be, static_cast<amrex::ParticleReal>(m_beryllium),
104  E_fusion, engine);
105 
106  // Compute momentum of beryllium in lab frame
107  const amrex::ParticleReal px_Be = static_cast<amrex::ParticleReal>(m_beryllium) * ux_Be;
108  const amrex::ParticleReal py_Be = static_cast<amrex::ParticleReal>(m_beryllium) * uy_Be;
109  const amrex::ParticleReal pz_Be = static_cast<amrex::ParticleReal>(m_beryllium) * uz_Be;
110 
111  // Compute momentum norm of second and third alphas in Beryllium rest frame
112  // Factor 0.5 is here because each alpha only gets half of the decay energy
113  constexpr amrex::ParticleReal gamma_Bestar = (1._prt + 0.5_prt*E_decay/(m_alpha*c_sq));
114  constexpr amrex::ParticleReal gamma_Bestar_sq_minus_one = gamma_Bestar*gamma_Bestar - 1._prt;
115  const amrex::ParticleReal p_Bestar = static_cast<amrex::ParticleReal>(m_alpha)*PhysConst::c*std::sqrt(gamma_Bestar_sq_minus_one);
116 
117  // Compute momentum of second alpha in Beryllium rest frame, assuming isotropic distribution
118  amrex::ParticleReal px_Bestar, py_Bestar, pz_Bestar;
119  ParticleUtils::RandomizeVelocity(px_Bestar, py_Bestar, pz_Bestar, p_Bestar, engine);
120 
121  // Next step is to convert momentum of second alpha to lab frame
122  amrex::ParticleReal px_alpha2, py_alpha2, pz_alpha2;
123  // Preliminary calculation: compute Beryllium velocity v_Be
124  const amrex::ParticleReal p_Be_sq = px_Be*px_Be + py_Be*py_Be + pz_Be*pz_Be;
125  const amrex::ParticleReal g_Be = std::sqrt(1._prt + p_Be_sq / (static_cast<amrex::ParticleReal>(mBe_sq)*c_sq));
126  const amrex::ParticleReal mg_Be = static_cast<amrex::ParticleReal>(m_beryllium)*g_Be;
127  const amrex::ParticleReal v_Bex = px_Be / mg_Be;
128  const amrex::ParticleReal v_Bey = py_Be / mg_Be;
129  const amrex::ParticleReal v_Bez = pz_Be / mg_Be;
130  const amrex::ParticleReal v_Be_sq = v_Bex*v_Bex + v_Bey*v_Bey + v_Bez*v_Bez;
131 
132  // Convert momentum of second alpha to lab frame, using equation (13) of F. Perez et al.,
133  // Phys.Plasmas.19.083104 (2012)
134  if ( v_Be_sq > std::numeric_limits<amrex::ParticleReal>::min() )
135  {
136  const amrex::ParticleReal vcDps = v_Bex*px_Bestar + v_Bey*py_Bestar + v_Bez*pz_Bestar;
137  const amrex::ParticleReal factor0 = (g_Be-1._prt)/v_Be_sq;
138  const amrex::ParticleReal factor = factor0*vcDps + static_cast<amrex::ParticleReal>(m_alpha)*gamma_Bestar*g_Be;
139  px_alpha2 = px_Bestar + v_Bex * factor;
140  py_alpha2 = py_Bestar + v_Bey * factor;
141  pz_alpha2 = pz_Bestar + v_Bez * factor;
142  }
143  else // If Beryllium velocity is zero, we are already in the lab frame
144  {
145  px_alpha2 = px_Bestar;
146  py_alpha2 = py_Bestar;
147  pz_alpha2 = pz_Bestar;
148  }
149 
150  // Compute momentum of third alpha in lab frame, using total momentum conservation
151  const amrex::ParticleReal px_alpha3 = px_Be - px_alpha2;
152  const amrex::ParticleReal py_alpha3 = py_Be - py_alpha2;
153  const amrex::ParticleReal pz_alpha3 = pz_Be - pz_alpha2;
154 
155  // Fill alpha species momentum data with the computed momentum (note that we actually
156  // create 6 alphas, 3 at the position of the proton and 3 at the position of the boron, so
157  // each computed momentum is used twice)
158  soa_alpha.m_rdata[PIdx::ux][idx_alpha_start] = ux_alpha1;
159  soa_alpha.m_rdata[PIdx::uy][idx_alpha_start] = uy_alpha1;
160  soa_alpha.m_rdata[PIdx::uz][idx_alpha_start] = uz_alpha1;
161  soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 1] = ux_alpha1;
162  soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 1] = uy_alpha1;
163  soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 1] = uz_alpha1;
164  soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 2] = px_alpha2/static_cast<amrex::ParticleReal>(m_alpha);
165  soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 2] = py_alpha2/static_cast<amrex::ParticleReal>(m_alpha);
166  soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 2] = pz_alpha2/static_cast<amrex::ParticleReal>(m_alpha);
167  soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 3] = px_alpha2/static_cast<amrex::ParticleReal>(m_alpha);
168  soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 3] = py_alpha2/static_cast<amrex::ParticleReal>(m_alpha);
169  soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 3] = pz_alpha2/static_cast<amrex::ParticleReal>(m_alpha);
170  soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 4] = px_alpha3/static_cast<amrex::ParticleReal>(m_alpha);
171  soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 4] = py_alpha3/static_cast<amrex::ParticleReal>(m_alpha);
172  soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 4] = pz_alpha3/static_cast<amrex::ParticleReal>(m_alpha);
173  soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 5] = px_alpha3/static_cast<amrex::ParticleReal>(m_alpha);
174  soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 5] = py_alpha3/static_cast<amrex::ParticleReal>(m_alpha);
175  soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 5] = pz_alpha3/static_cast<amrex::ParticleReal>(m_alpha);
176  }
177 
178 }
179 
180 #endif // WARPX_PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > ParticleTileType
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
typename WarpXParticleContainer::ParticleTileType ParticleTileType
Definition: ParticleUtils.cpp:32
typename ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition: ParticleUtils.cpp:33
typename ParticleBins::index_type index_type
Definition: ParticleUtils.cpp:35
DenseBins< ParticleTileDataType > ParticleBins
Definition: ParticleUtils.cpp:34
typename WarpXParticleContainer::ParticleType ParticleType
Definition: ParticleUtils.cpp:31
static constexpr auto c
vacuum speed of light [m/s]
Definition: constant.H:44
static constexpr auto m_p
proton mass [kg]
Definition: constant.H:54
static constexpr auto m_u
dalton: unified atomic mass unit [kg]
Definition: constant.H:56
static constexpr auto q_e
elementary charge [C]
Definition: constant.H:50
@ uz
Definition: NamedComponentParticleContainer.H:34
@ uy
Definition: NamedComponentParticleContainer.H:34
@ ux
Definition: NamedComponentParticleContainer.H:34
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType