1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
|
/* Copyright 2021 Neil Zaim
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
#define PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
#include "Particles/WarpXParticleContainer.H"
#include "Utils/ParticleUtils.H"
#include "Utils/WarpXConst.H"
#include <AMReX_DenseBins.H>
#include <AMReX_Random.H>
#include <AMReX_REAL.H>
#include <cmath>
#include <limits>
namespace {
// Define shortcuts for frequently-used type names
using SoaData_type = WarpXParticleContainer::ParticleTileType::ParticleTileDataType;
using ParticleType = WarpXParticleContainer::ParticleType;
using ParticleBins = amrex::DenseBins<ParticleType>;
using index_type = ParticleBins::index_type;
/**
* \brief This function initializes the momentum of the alpha particles produced from
* proton-boron fusion. The momentum is initialized by assuming that the fusion of a proton
* with a boron nucleus into 3 alphas takes place in two steps. In the first step, the proton
* and the boron fuse into a beryllium nucleus and an alpha particle. In the second step, the
* beryllium decays into two alpha particles. The first step produces 8.59009 MeV of kinetic
* energy while the second step produces 91.8984 keV of kinetic energy. This two-step process
* is considered to be the dominant process of proton+boron fusion into alphas (see
* Becker et al., Zeitschrift für Physik A Atomic Nuclei, 327(3), 341-355 (1987)).
* For each step, we assume in this function that the particles are emitted isotropically in
* the corresponding center of mass frame (center of mass frame of proton + boron for the
* creation of first alpha+beryllium and rest frame of beryllium for the creation of second and
* third alphas). This isotropic assumption is exact for the second step but is only an
* approximation for the first step.
*
* @param[in] soa_1 struct of array data of the first colliding species (can be either proton
* or boron)
* @param[in] soa_2 struct of array data of the second colliding species (can be either proton
* or boron)
* @param[out] soa_alpha struct of array data of the alpha species
* @param[in] idx_1 index of first colliding macroparticle
* @param[in] idx_2 index of second colliding macroparticle
* @param[in] idx_alpha_start index of first produced alpha macroparticle
* @param[in] m1 mass of first colliding species
* @param[in] m2 mass of second colliding species
* @param[in] engine the random engine
*/
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void ProtonBoronFusionInitializeMomentum (
const SoaData_type& soa_1, const SoaData_type& soa_2,
SoaData_type& soa_alpha,
const index_type& idx_1, const index_type& idx_2,
const index_type& idx_alpha_start,
const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
const amrex::RandomEngine& engine)
{
// General notations in this function:
// x_sq denotes the square of x
// x_star denotes the value of x in the proton+boron center of mass frame
// x_Bestar denotes the value of x in the Beryllium rest frame
using namespace amrex::literals;
constexpr amrex::ParticleReal mev_to_joule = PhysConst::q_e*1.e6_prt;
// Energy produced in the fusion reaction proton + boron11 -> Beryllium8 + alpha
// cf. Janis book of proton-induced cross-sections (2019)
constexpr amrex::ParticleReal E_fusion = 8.59009_prt*mev_to_joule;
// Energy produced when Beryllium8 decays into two alphas
// cf. JEFF-3.3 radioactive decay data library (2017)
constexpr amrex::ParticleReal E_decay = 0.0918984_prt*mev_to_joule;
constexpr amrex::ParticleReal m_alpha = PhysConst::m_p * 3.97369_prt;
constexpr amrex::ParticleReal m_beryllium = PhysConst::m_p * 7.94748_prt;
constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c;
constexpr amrex::ParticleReal inv_csq = 1._prt / ( c_sq );
// Rest energy of proton+boron
const amrex::ParticleReal E_rest_pb = (m1 + m2)*c_sq;
// Rest energy of alpha+beryllium
constexpr amrex::ParticleReal E_rest_abe = (m_alpha + m_beryllium)*c_sq;
// These constexprs underflow in single precision because we use SI units and sometimes the
// code won't compile. Nuclear fusion module does not currently work with single precision
// anyways so these #ifdef are just here to make sure that the code always compiles with
// single precision.
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
const amrex::ParticleReal ma_sq = m_alpha*m_alpha;
const amrex::ParticleReal mBe_sq = m_beryllium*m_beryllium;
#else
constexpr amrex::ParticleReal ma_sq = m_alpha*m_alpha;
constexpr amrex::ParticleReal mBe_sq = m_beryllium*m_beryllium;
#endif
// Normalized momentum of colliding particles (proton and boron)
const amrex::ParticleReal u1x = soa_1.m_rdata[PIdx::ux][idx_1];
const amrex::ParticleReal u1y = soa_1.m_rdata[PIdx::uy][idx_1];
const amrex::ParticleReal u1z = soa_1.m_rdata[PIdx::uz][idx_1];
const amrex::ParticleReal u2x = soa_2.m_rdata[PIdx::ux][idx_2];
const amrex::ParticleReal u2y = soa_2.m_rdata[PIdx::uy][idx_2];
const amrex::ParticleReal u2z = soa_2.m_rdata[PIdx::uz][idx_2];
// Compute Lorentz factor gamma in the lab frame
const amrex::ParticleReal g1 = std::sqrt( 1._prt + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq );
const amrex::ParticleReal g2 = std::sqrt( 1._prt + (u2x*u2x+u2y*u2y+u2z*u2z)*inv_csq );
// Compute momenta
const amrex::ParticleReal p1x = u1x * m1;
const amrex::ParticleReal p1y = u1y * m1;
const amrex::ParticleReal p1z = u1z * m1;
const amrex::ParticleReal p2x = u2x * m2;
const amrex::ParticleReal p2y = u2y * m2;
const amrex::ParticleReal p2z = u2z * m2;
// Square norm of the total (sum between the two particles) momenta in the lab frame
const amrex::ParticleReal p_total_sq = std::pow(p1x+p2x, 2) +
std::pow(p1y+p2y, 2) +
std::pow(p1z+p2z, 2);
// Total energy of proton+boron in the lab frame
const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq;
// Total energy squared of proton+boron in the center of mass frame, calculated using the
// Lorentz invariance of the four-momentum norm
const amrex::ParticleReal E_star_sq = E_lab*E_lab - c_sq*p_total_sq;
// Total energy squared of beryllium+alpha in the center of mass frame
// In principle, the term - E_rest_pb + E_rest_abe + E_fusion is not needed and equal to
// zero (i.e. the energy liberated during fusion is equal to the mass difference). However,
// due to possible inconsistencies in how the mass is defined in the code (e.g. currently,
// the mass of hydrogen is the mass of the proton, not including the electron, while the
// mass of the other elements is the atomic mass, that includes the electrons mass), it is
// probably more robust to subtract the rest masses and to add the fusion energy to the
// total kinetic energy.
const amrex::ParticleReal E_star_f_sq = std::pow(std::sqrt(E_star_sq)
- E_rest_pb + E_rest_abe + E_fusion, 2);
// Square of the norm of the momentum of Beryllium or alpha in the center of mass frame
// Formula obtained by inverting E^2 = p^2*c^2 + m^2*c^4 in the COM frame for each particle
const amrex::ParticleReal p_star_f_sq =
E_star_f_sq*0.25_prt*inv_csq - (ma_sq + mBe_sq)*c_sq*0.5_prt +
std::pow(c_sq, 3)*0.25_prt*std::pow(ma_sq - mBe_sq, 2)/E_star_f_sq;
// Compute momentum of first alpha in the center of mass frame, assuming isotropic
// distribution
amrex::ParticleReal px_star, py_star, pz_star;
ParticleUtils::RandomizeVelocity(px_star, py_star, pz_star, std::sqrt(p_star_f_sq),
engine);
// Next step is to convert momentum of first alpha to lab frame
amrex::ParticleReal px_alpha1, py_alpha1, pz_alpha1;
// Preliminary calculation: compute center of mass velocity vc
const amrex::ParticleReal mass_g = m1 * g1 + m2 * g2;
const amrex::ParticleReal vcx = (p1x+p2x) / mass_g;
const amrex::ParticleReal vcy = (p1y+p2y) / mass_g;
const amrex::ParticleReal vcz = (p1z+p2z) / mass_g;
const amrex::ParticleReal vc_sq = vcx*vcx + vcy*vcy + vcz*vcz;
// Convert momentum of first alpha to lab frame, using equation (13) of F. Perez et al.,
// Phys.Plasmas.19.083104 (2012)
if ( vc_sq > std::numeric_limits<amrex::ParticleReal>::min() )
{
const amrex::ParticleReal gc = 1._prt / std::sqrt( 1._prt - vc_sq*inv_csq );
const amrex::ParticleReal g_star = std::sqrt(1._prt + p_star_f_sq / (ma_sq*c_sq));
const amrex::ParticleReal vcDps = vcx*px_star + vcy*py_star + vcz*pz_star;
const amrex::ParticleReal factor0 = (gc-1._prt)/vc_sq;
const amrex::ParticleReal factor = factor0*vcDps + m_alpha*g_star*gc;
px_alpha1 = px_star + vcx * factor;
py_alpha1 = py_star + vcy * factor;
pz_alpha1 = pz_star + vcz * factor;
}
else // If center of mass velocity is zero, we are already in the lab frame
{
px_alpha1 = px_star;
py_alpha1 = py_star;
pz_alpha1 = pz_star;
}
// Compute momentum of beryllium in lab frame, using total momentum conservation
const amrex::ParticleReal px_Be = p1x + p2x - px_alpha1;
const amrex::ParticleReal py_Be = p1y + p2y - py_alpha1;
const amrex::ParticleReal pz_Be = p1z + p2z - pz_alpha1;
// Compute momentum norm of second and third alphas in Beryllium rest frame
// Factor 0.5 is here because each alpha only gets half of the decay energy
constexpr amrex::ParticleReal gamma_Bestar = (1._prt + 0.5_prt*E_decay/(m_alpha*c_sq));
constexpr amrex::ParticleReal gamma_Bestar_sq_minus_one = gamma_Bestar*gamma_Bestar - 1._prt;
const amrex::ParticleReal p_Bestar = m_alpha*PhysConst::c*std::sqrt(gamma_Bestar_sq_minus_one);
// Compute momentum of second alpha in Beryllium rest frame, assuming isotropic distribution
amrex::ParticleReal px_Bestar, py_Bestar, pz_Bestar;
ParticleUtils::RandomizeVelocity(px_Bestar, py_Bestar, pz_Bestar, p_Bestar, engine);
// Next step is to convert momentum of second alpha to lab frame
amrex::ParticleReal px_alpha2, py_alpha2, pz_alpha2;
// Preliminary calculation: compute Beryllium velocity v_Be
const amrex::ParticleReal p_Be_sq = px_Be*px_Be + py_Be*py_Be + pz_Be*pz_Be;
const amrex::ParticleReal g_Be = std::sqrt(1._prt + p_Be_sq / (mBe_sq*c_sq));
const amrex::ParticleReal mg_Be = m_beryllium*g_Be;
const amrex::ParticleReal v_Bex = px_Be / mg_Be;
const amrex::ParticleReal v_Bey = py_Be / mg_Be;
const amrex::ParticleReal v_Bez = pz_Be / mg_Be;
const amrex::ParticleReal v_Be_sq = v_Bex*v_Bex + v_Bey*v_Bey + v_Bez*v_Bez;
// Convert momentum of second alpha to lab frame, using equation (13) of F. Perez et al.,
// Phys.Plasmas.19.083104 (2012)
if ( v_Be_sq > std::numeric_limits<amrex::ParticleReal>::min() )
{
const amrex::ParticleReal vcDps = v_Bex*px_Bestar + v_Bey*py_Bestar + v_Bez*pz_Bestar;
const amrex::ParticleReal factor0 = (g_Be-1._prt)/v_Be_sq;
const amrex::ParticleReal factor = factor0*vcDps + m_alpha*gamma_Bestar*g_Be;
px_alpha2 = px_Bestar + v_Bex * factor;
py_alpha2 = py_Bestar + v_Bey * factor;
pz_alpha2 = pz_Bestar + v_Bez * factor;
}
else // If Beryllium velocity is zero, we are already in the lab frame
{
px_alpha2 = px_Bestar;
py_alpha2 = py_Bestar;
pz_alpha2 = pz_Bestar;
}
// Compute momentum of third alpha in lab frame, using total momentum conservation
const amrex::ParticleReal px_alpha3 = px_Be - px_alpha2;
const amrex::ParticleReal py_alpha3 = py_Be - py_alpha2;
const amrex::ParticleReal pz_alpha3 = pz_Be - pz_alpha2;
// Fill alpha species momentum data with the computed momentum (note that we actually
// create 6 alphas, 3 at the position of the proton and 3 at the position of the boron, so
// each computed momentum is used twice)
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start] = px_alpha1/m_alpha;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start] = py_alpha1/m_alpha;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start] = pz_alpha1/m_alpha;
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 1] = px_alpha1/m_alpha;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 1] = py_alpha1/m_alpha;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 1] = pz_alpha1/m_alpha;
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 2] = px_alpha2/m_alpha;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 2] = py_alpha2/m_alpha;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 2] = pz_alpha2/m_alpha;
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 3] = px_alpha2/m_alpha;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 3] = py_alpha2/m_alpha;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 3] = pz_alpha2/m_alpha;
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 4] = px_alpha3/m_alpha;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 4] = py_alpha3/m_alpha;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 4] = pz_alpha3/m_alpha;
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 5] = px_alpha3/m_alpha;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 5] = py_alpha3/m_alpha;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 5] = pz_alpha3/m_alpha;
}
}
#endif // PROTON_BORON_FUSION_INITIALIZE_MOMENTUM_H
|
