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
|
/* 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 "TwoProductFusionUtil.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_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;
// The constexprs ma_sq and mBe_sq underflow in single precision because we use SI units,
// which can cause compilation to fail or generate a warning, so we're explicitly setting
// them as double. Note that nuclear fusion module does not currently work with single
// precision anyways.
constexpr double m_alpha = PhysConst::m_u * 4.002602_prt;
constexpr double m_beryllium = PhysConst::m_p * 7.94748_prt;
constexpr double mBe_sq = m_beryllium*m_beryllium;
constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c;
// Compute the resulting momenta of the alpha and beryllium particles
// produced in the fusion reaction
amrex::ParticleReal ux_alpha1 = 0.0, uy_alpha1 = 0.0, uz_alpha1 = 0.0;
amrex::ParticleReal ux_Be = 0.0, uy_Be = 0.0, uz_Be = 0.0;
TwoProductFusionComputeProductMomenta (
soa_1.m_rdata[PIdx::ux][idx_1],
soa_1.m_rdata[PIdx::uy][idx_1],
soa_1.m_rdata[PIdx::uz][idx_1], m1,
soa_2.m_rdata[PIdx::ux][idx_2],
soa_2.m_rdata[PIdx::uy][idx_2],
soa_2.m_rdata[PIdx::uz][idx_2], m2,
ux_alpha1, uy_alpha1, uz_alpha1, m_alpha,
ux_Be, uy_Be, uz_Be, m_beryllium,
E_fusion, engine);
// Compute momentum of beryllium in lab frame
const amrex::ParticleReal px_Be = m_beryllium * ux_Be;
const amrex::ParticleReal py_Be = m_beryllium * uy_Be;
const amrex::ParticleReal pz_Be = m_beryllium * uz_Be;
// 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] = ux_alpha1;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start] = uy_alpha1;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start] = uz_alpha1;
soa_alpha.m_rdata[PIdx::ux][idx_alpha_start + 1] = ux_alpha1;
soa_alpha.m_rdata[PIdx::uy][idx_alpha_start + 1] = uy_alpha1;
soa_alpha.m_rdata[PIdx::uz][idx_alpha_start + 1] = uz_alpha1;
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
|
