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
|
/* Copyright 2021 Neil Zaim
*
* This file is part of WarpX.
*
* License: BSD-3-Clause-LBNL
*/
#ifndef SINGLE_NUCLEAR_FUSION_EVENT_H_
#define SINGLE_NUCLEAR_FUSION_EVENT_H_
#include "BoschHaleFusionCrossSection.H"
#include "ProtonBoronFusionCrossSection.H"
#include "Particles/Collision/BinaryCollision/BinaryCollisionUtils.H"
#include "Utils/WarpXConst.H"
#include <AMReX_Algorithm.H>
#include <AMReX_Math.H>
#include <AMReX_Random.H>
#include <AMReX_REAL.H>
#include <cmath>
/**
* \brief This function computes whether the collision between two particles result in a
* nuclear fusion event, using the algorithm described in Higginson et al., Journal of
* Computational Physics 388, 439-453 (2019). If nuclear fusion occurs, the mask is set to true
* for that given pair of particles and the weight of the produced particles is stored in
* p_pair_reaction_weight.
*
* @tparam index_type type of the index argument
* @param[in] u1x,u1y,u1z momenta of the first colliding particle
* @param[in] u2x,u2y,u2z momenta of the second colliding particle
* @param[in] m1,m2 masses
* @param[in] w1,w2 effective weight of the colliding particles
* @param[in] dt is the time step length between two collision calls.
* @param[in] dV is the volume of the corresponding cell.
* @param[in] pair_index is the index of the colliding pair
* @param[out] p_mask is a mask that will be set to true if fusion occurs for that pair
* @param[out] p_pair_reaction_weight stores the weight of the product particles
* @param[in] fusion_multiplier factor used to increase the number of fusion events by
* decreasing the weight of the produced particles
* @param[in] multiplier_ratio factor used to take into account unsampled pairs (i.e. the fact
* that a particle only collides with one or few particles of the other species)
* @param[in] probability_threshold probability threshold above which we decrease the fusion
* multiplier
* @param[in] probability_target_value if the probability threshold is exceeded, this is used
* to determine by how much the fusion multiplier is reduced
* @param[in] fusion_type the physical fusion process to model
* @param[in] engine the random engine.
*/
template <typename index_type>
AMREX_GPU_HOST_DEVICE AMREX_INLINE
void SingleNuclearFusionEvent (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 w1, amrex::ParticleReal w2,
const amrex::Real& dt, const amrex::ParticleReal& dV, const int& pair_index,
index_type* AMREX_RESTRICT p_mask,
amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
const amrex::ParticleReal& fusion_multiplier,
const int& multiplier_ratio,
const amrex::ParticleReal& probability_threshold,
const amrex::ParticleReal& probability_target_value,
const NuclearFusionType& fusion_type,
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 center of mass frame
using namespace amrex::literals;
using namespace amrex::Math;
const amrex::ParticleReal w_min = amrex::min(w1, w2);
const amrex::ParticleReal w_max = amrex::max(w1, w2);
constexpr auto one_pr = amrex::ParticleReal(1.);
constexpr auto inv_four_pr = amrex::ParticleReal(1./4.);
constexpr amrex::ParticleReal c_sq = PhysConst::c * PhysConst::c;
constexpr amrex::ParticleReal inv_csq = one_pr / ( c_sq );
const amrex::ParticleReal m1_sq = m1*m1;
const amrex::ParticleReal m2_sq = m2*m2;
// Compute Lorentz factor gamma in the lab frame
const amrex::ParticleReal g1 = std::sqrt( one_pr + (u1x*u1x+u1y*u1y+u1z*u1z)*inv_csq );
const amrex::ParticleReal g2 = std::sqrt( one_pr + (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 = powi<2>(p1x + p2x) +
powi<2>(p1y+p2y) +
powi<2>(p1z+p2z);
// Total energy in the lab frame
const amrex::ParticleReal E_lab = (m1 * g1 + m2 * g2) * c_sq;
// Total energy squared 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;
// Kinetic energy in the center of mass frame
const amrex::ParticleReal E_star = std::sqrt(E_star_sq);
const amrex::ParticleReal E_kin_star = E_star - (m1 + m2)*c_sq;
// Compute fusion cross section as a function of kinetic energy in the center of mass frame
auto fusion_cross_section = amrex::ParticleReal(0.);
if (fusion_type == NuclearFusionType::ProtonBoronToAlphas)
{
fusion_cross_section = ProtonBoronFusionCrossSection(E_kin_star);
}
else if ((fusion_type == NuclearFusionType::DeuteriumTritiumToNeutronHelium)
|| (fusion_type == NuclearFusionType::DeuteriumDeuteriumToProtonTritium)
|| (fusion_type == NuclearFusionType::DeuteriumDeuteriumToNeutronHelium))
{
fusion_cross_section = BoschHaleFusionCrossSection(E_kin_star, fusion_type, m1, m2);
}
// Square of the norm of the momentum of one of the particles 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
// The expression below is specifically written in a form that avoids returning
// small negative numbers due to machine precision errors, for low-energy particles
const amrex::ParticleReal E_ratio = E_star/((m1 + m2)*c_sq);
const amrex::ParticleReal p_star_sq = m1*m2*c_sq * ( powi<2>(E_ratio) - one_pr )
+ powi<2>(m1 - m2)*c_sq*inv_four_pr * powi<2>( E_ratio - 1._prt/E_ratio );
// Lorentz factors in the center of mass frame
const amrex::ParticleReal g1_star = std::sqrt(one_pr + p_star_sq / (m1_sq*c_sq));
const amrex::ParticleReal g2_star = std::sqrt(one_pr + p_star_sq / (m2_sq*c_sq));
// relative velocity in the center of mass frame
const amrex::ParticleReal v_rel = std::sqrt(p_star_sq) * (one_pr/(m1*g1_star) +
one_pr/(m2*g2_star));
// Fusion cross section and relative velocity are computed in the center of mass frame.
// On the other hand, the particle densities (weight over volume) in the lab frame are used. To
// take into account this discrepancy, we need to multiply the fusion probability by the ratio
// between the Lorentz factors in the COM frame and the Lorentz factors in the lab frame
// (see Perez et al., Phys.Plasmas.19.083104 (2012))
const amrex::ParticleReal lab_to_COM_factor = g1_star*g2_star/(g1*g2);
// First estimate of probability to have fusion reaction
amrex::ParticleReal probability_estimate = multiplier_ratio * fusion_multiplier *
lab_to_COM_factor * w_max * fusion_cross_section * v_rel * dt / dV;
// Effective fusion multiplier
amrex::ParticleReal fusion_multiplier_eff = fusion_multiplier;
// If the fusion probability is too high and the fusion multiplier greater than one, we risk to
// systematically underestimate the fusion yield. In this case, we reduce the fusion multiplier
// to reduce the fusion probability
if (probability_estimate > probability_threshold)
{
// We aim for a fusion probability of probability_target_value but take into account
// the constraint that the fusion_multiplier cannot be smaller than one
fusion_multiplier_eff = amrex::max(fusion_multiplier *
probability_target_value / probability_estimate , one_pr);
probability_estimate *= fusion_multiplier_eff/fusion_multiplier;
}
// Compute actual fusion probability that is always between zero and one
// In principle this is obtained by computing 1 - exp(-probability_estimate)
// However, the computation of this quantity can fail numerically when probability_estimate is
// too small (e.g. exp(-probability_estimate) returns 1 and the computation returns 0).
// In this case, we simply use "probability_estimate" instead of 1 - exp(-probability_estimate)
// The threshold exp_threshold at which we switch between the two formulas is determined by the
// fact that computing the exponential is only useful if it can resolve the x^2/2 term of its
// Taylor expansion, i.e. the square of probability_estimate should be greater than the
// machine epsilon.
#ifdef AMREX_SINGLE_PRECISION_PARTICLES
constexpr auto exp_threshold = amrex::ParticleReal(1.e-3);
#else
constexpr auto exp_threshold = amrex::ParticleReal(5.e-8);
#endif
const amrex::ParticleReal probability = (probability_estimate < exp_threshold) ?
probability_estimate: one_pr - std::exp(-probability_estimate);
// Get a random number
const amrex::ParticleReal random_number = amrex::Random(engine);
// If we have a fusion event, set the mask the true and fill the product weight array
if (random_number < probability)
{
p_mask[pair_index] = true;
p_pair_reaction_weight[pair_index] = w_min/fusion_multiplier_eff;
}
else
{
p_mask[pair_index] = false;
}
}
#endif // SINGLE_NUCLEAR_FUSION_EVENT_H_
|
