WarpX
NuclearFusionFunc.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_NUCLEAR_FUSION_FUNC_H_
9 #define WARPX_NUCLEAR_FUSION_FUNC_H_
10 
12 
18 #include "Utils/TextMsg.H"
19 #include "WarpX.H"
20 
21 #include <AMReX_Algorithm.H>
22 #include <AMReX_DenseBins.H>
23 #include <AMReX_ParmParse.H>
24 #include <AMReX_Random.H>
25 #include <AMReX_REAL.H>
26 #include <AMReX_Vector.H>
27 
37 {
38  // Define shortcuts for frequently-used type names
45 
46 public:
50  NuclearFusionFunc () = default;
51 
59  NuclearFusionFunc (const std::string& collision_name, MultiParticleContainer const * const mypc,
60  const bool isSameSpecies):
61  m_fusion_multiplier{amrex::ParticleReal{1.0}}, // default fusion multiplier
62  m_probability_threshold{amrex::ParticleReal{0.02}}, // default fusion probability threshold
63  m_probability_target_value{amrex::ParticleReal{0.002}}, // default fusion probability target_value
65  m_isSameSpecies{isSameSpecies}
66  {
67 
68 #ifdef AMREX_SINGLE_PRECISION_PARTICLES
69  WARPX_ABORT_WITH_MESSAGE("Nuclear fusion module does not currently work with single precision");
70 #endif
71 
72  const amrex::ParmParse pp_collision_name(collision_name);
74  pp_collision_name, "fusion_multiplier", m_fusion_multiplier);
76  pp_collision_name, "fusion_probability_threshold", m_probability_threshold);
78  pp_collision_name, "fusion_probability_target_value",
80 
86  }
87 
88  struct Executor {
128  index_type const I1s, index_type const I1e,
129  index_type const I2s, index_type const I2e,
130  index_type const* AMREX_RESTRICT I1,
131  index_type const* AMREX_RESTRICT I2,
132  const SoaData_type& soa_1, const SoaData_type& soa_2,
133  GetParticlePosition<PIdx> /*get_position_1*/, GetParticlePosition<PIdx> /*get_position_2*/,
134  amrex::ParticleReal const /*q1*/, amrex::ParticleReal const /*q2*/,
135  amrex::ParticleReal const m1, amrex::ParticleReal const m2,
136  amrex::Real const dt, amrex::Real const dV, index_type coll_idx,
137  index_type const cell_start_pair, index_type* AMREX_RESTRICT p_mask,
138  index_type* AMREX_RESTRICT p_pair_indices_1, index_type* AMREX_RESTRICT p_pair_indices_2,
139  amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight,
140  amrex::RandomEngine const& engine) const
141  {
142  amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
143  amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
144  amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
145  amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
146 
147  amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
148  amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
149  amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
150  amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
151 
152  // Number of macroparticles of each species
153  const index_type NI1 = I1e - I1s;
154  const index_type NI2 = I2e - I2s;
155  const index_type max_N = amrex::max(NI1,NI2);
156  const index_type min_N = amrex::min(NI1,NI2);
157 
158  index_type pair_index = cell_start_pair + coll_idx;
159 
160  // Because the number of particles of each species is not always equal (NI1 != NI2
161  // in general), some macroparticles will be paired with multiple macroparticles of the
162  // other species and we need to decrease their weight accordingly.
163  // c1 corresponds to the minimum number of times a particle of species 1 will be paired
164  // with a particle of species 2. Same for c2.
165  // index_type(1): https://github.com/AMReX-Codes/amrex/pull/3684
166  const index_type c1 = amrex::max(NI2/NI1, index_type(1));
167  const index_type c2 = amrex::max(NI1/NI2, index_type(1));
168 
169  // multiplier ratio to take into account unsampled pairs
170  const auto multiplier_ratio = static_cast<int>(
171  m_isSameSpecies ? 2*max_N - 1 : max_N);
172 
173 #if (defined WARPX_DIM_RZ)
174  amrex::ParticleReal * const AMREX_RESTRICT theta1 = soa_1.m_rdata[PIdx::theta];
175  amrex::ParticleReal * const AMREX_RESTRICT theta2 = soa_2.m_rdata[PIdx::theta];
176 #endif
177  index_type i1 = I1s + coll_idx;
178  index_type i2 = I2s + coll_idx;
179  // we will start from collision number = coll_idx and then add
180  // stride (smaller set size) until we do all collisions (larger set size)
181  for (index_type k = coll_idx; k < max_N; k += min_N)
182  {
183  // c1k : how many times the current particle of species 1 is paired with a particle
184  // of species 2. Same for c2k.
185  const index_type c1k = (k%NI1 < max_N%NI1) ? c1 + 1: c1;
186  const index_type c2k = (k%NI2 < max_N%NI2) ? c2 + 1: c2;
187 
188 #if (defined WARPX_DIM_RZ)
189  /* In RZ geometry, macroparticles can collide with other macroparticles
190  * in the same *cylindrical* cell. For this reason, collisions between macroparticles
191  * are actually not local in space. In this case, the underlying assumption is that
192  * particles within the same cylindrical cell represent a cylindrically-symmetry
193  * momentum distribution function. Therefore, here, we temporarily rotate the
194  * momentum of one of the macroparticles in agreement with this cylindrical symmetry.
195  * (This is technically only valid if we use only the m=0 azimuthal mode in the simulation;
196  * there is a corresponding assert statement at initialization.) */
197  amrex::ParticleReal const theta = theta2[I2[i2]]-theta1[I1[i1]];
198  amrex::ParticleReal const u1xbuf = u1x[I1[i1]];
199  u1x[I1[i1]] = u1xbuf*std::cos(theta) - u1y[I1[i1]]*std::sin(theta);
200  u1y[I1[i1]] = u1xbuf*std::sin(theta) + u1y[I1[i1]]*std::cos(theta);
201 #endif
202 
204  u1x[ I1[i1] ], u1y[ I1[i1] ], u1z[ I1[i1] ],
205  u2x[ I2[i2] ], u2y[ I2[i2] ], u2z[ I2[i2] ],
206  m1, m2, w1[ I1[i1] ]/c1k, w2[ I2[i2] ]/c2k,
207  dt, dV, static_cast<int>(pair_index), p_mask, p_pair_reaction_weight,
208  m_fusion_multiplier, multiplier_ratio,
211  m_fusion_type, engine);
212 
213 #if (defined WARPX_DIM_RZ)
214  amrex::ParticleReal const u1xbuf_new = u1x[I1[i1]];
215  u1x[I1[i1]] = u1xbuf_new*std::cos(-theta) - u1y[I1[i1]]*std::sin(-theta);
216  u1y[I1[i1]] = u1xbuf_new*std::sin(-theta) + u1y[I1[i1]]*std::cos(-theta);
217 #endif
218 
219  p_pair_indices_1[pair_index] = I1[i1];
220  p_pair_indices_2[pair_index] = I2[i2];
221  if (max_N == NI1) {
222  i1 += min_N;
223  }
224  if (max_N == NI2) {
225  i2 += min_N;
226  }
227  pair_index += min_N;
228  }
229  }
230 
231  amrex::ParticleReal m_fusion_multiplier;
232  amrex::ParticleReal m_probability_threshold;
233  amrex::ParticleReal m_probability_target_value;
236  };
237 
238  [[nodiscard]] Executor const& executor () const { return m_exe; }
239 
240 private:
241  // Factor used to increase the number of fusion reaction by decreasing the weight of the
242  // produced particles
243  amrex::ParticleReal m_fusion_multiplier;
244  // If the fusion multiplier is too high and results in a fusion probability that approaches
245  // 1, there is a risk of underestimating the total fusion yield. In these cases, we reduce
246  // the fusion multiplier used in a given collision. m_probability_threshold is the fusion
247  // probability threshold above which we reduce the fusion multiplier.
248  // m_probability_target_value is the target probability used to determine by how much
249  // the fusion multiplier should be reduced.
250  amrex::ParticleReal m_probability_threshold;
251  amrex::ParticleReal m_probability_target_value;
254 
256 };
257 
258 #endif // WARPX_NUCLEAR_FUSION_FUNC_H_
#define AMREX_RESTRICT
#define AMREX_INLINE
#define AMREX_GPU_HOST_DEVICE
NuclearFusionType
Definition: BinaryCollisionUtils.H:26
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)
This function computes whether the collision between two particles result in a nuclear fusion event,...
Definition: SingleNuclearFusionEvent.H:54
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition: TextMsg.H:15
Definition: MultiParticleContainer.H:66
This functor does binary nuclear fusions on a single cell. Particles of the two reacting species are ...
Definition: NuclearFusionFunc.H:37
NuclearFusionFunc(const std::string &collision_name, MultiParticleContainer const *const mypc, const bool isSameSpecies)
Constructor of the NuclearFusionFunc class.
Definition: NuclearFusionFunc.H:59
NuclearFusionFunc()=default
Default constructor of the NuclearFusionFunc class.
Executor const & executor() const
Definition: NuclearFusionFunc.H:238
amrex::ParticleReal m_fusion_multiplier
Definition: NuclearFusionFunc.H:243
bool m_isSameSpecies
Definition: NuclearFusionFunc.H:253
WarpXParticleContainer::ParticleType ParticleType
Definition: NuclearFusionFunc.H:39
ParticleBins::index_type index_type
Definition: NuclearFusionFunc.H:43
amrex::ParticleReal m_probability_threshold
Definition: NuclearFusionFunc.H:250
Executor m_exe
Definition: NuclearFusionFunc.H:255
amrex::ParticleReal m_probability_target_value
Definition: NuclearFusionFunc.H:251
NuclearFusionType m_fusion_type
Definition: NuclearFusionFunc.H:252
unsigned int index_type
ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > ParticleTileType
Definition: BinaryCollisionUtils.cpp:18
NuclearFusionType get_nuclear_fusion_type(const std::string &collision_name, MultiParticleContainer const *const mypc)
Definition: BinaryCollisionUtils.cpp:40
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & max(const T &a, const T &b) noexcept
AMREX_GPU_HOST_DEVICE constexpr AMREX_FORCE_INLINE const T & min(const T &a, const T &b) noexcept
float dt
Definition: stencil.py:442
int queryWithParser(const amrex::ParmParse &a_pp, char const *const str, T &val)
Definition: ParserUtils.H:137
Definition: NuclearFusionFunc.H:88
NuclearFusionType m_fusion_type
Definition: NuclearFusionFunc.H:234
amrex::ParticleReal m_probability_target_value
Definition: NuclearFusionFunc.H:233
bool m_isSameSpecies
Definition: NuclearFusionFunc.H:235
amrex::ParticleReal m_probability_threshold
Definition: NuclearFusionFunc.H:232
AMREX_GPU_HOST_DEVICE AMREX_INLINE void operator()(index_type const I1s, index_type const I1e, index_type const I2s, index_type const I2e, index_type const *AMREX_RESTRICT I1, index_type const *AMREX_RESTRICT I2, const SoaData_type &soa_1, const SoaData_type &soa_2, GetParticlePosition< PIdx >, GetParticlePosition< PIdx >, amrex::ParticleReal const, amrex::ParticleReal const, amrex::ParticleReal const m1, amrex::ParticleReal const m2, amrex::Real const dt, amrex::Real const dV, index_type coll_idx, index_type const cell_start_pair, index_type *AMREX_RESTRICT p_mask, index_type *AMREX_RESTRICT p_pair_indices_1, index_type *AMREX_RESTRICT p_pair_indices_2, amrex::ParticleReal *AMREX_RESTRICT p_pair_reaction_weight, amrex::RandomEngine const &engine) const
Executor of the NuclearFusionFunc class. Performs nuclear fusions at the cell level using the algorit...
Definition: NuclearFusionFunc.H:127
amrex::ParticleReal m_fusion_multiplier
Definition: NuclearFusionFunc.H:231
@ theta
RZ needs all three position components.
Definition: NamedComponentParticleContainer.H:36
@ uz
Definition: NamedComponentParticleContainer.H:34
@ w
weight
Definition: NamedComponentParticleContainer.H:33
@ uy
Definition: NamedComponentParticleContainer.H:34
@ ux
Definition: NamedComponentParticleContainer.H:34
GpuArray< ParticleReal *, NAR > m_rdata
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType