WarpX
ParticleCreationFunc.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_PARTICLE_CREATION_FUNC_H_
9 #define WARPX_PARTICLE_CREATION_FUNC_H_
10 
11 #include "BinaryCollisionUtils.H"
12 
18 
19 #include <AMReX_DenseBins.H>
20 #include <AMReX_GpuAtomic.H>
21 #include <AMReX_GpuDevice.H>
22 #include <AMReX_GpuContainers.H>
23 #include <AMReX_INT.H>
24 #include <AMReX_Random.H>
25 #include <AMReX_REAL.H>
26 #include <AMReX_Vector.H>
27 
33 {
34  // Define shortcuts for frequently-used type names
41 
42 public:
46  ParticleCreationFunc () = default;
47 
54  ParticleCreationFunc (const std::string& collision_name, MultiParticleContainer const * mypc);
55 
96  const index_type& n_total_pairs,
97  const SoaData_type& soa_1, const SoaData_type& soa_2,
98  const amrex::Vector<WarpXParticleContainer*>& pc_products,
99  ParticleTileType** AMREX_RESTRICT tile_products,
100  const amrex::ParticleReal& m1, const amrex::ParticleReal& m2,
101  const amrex::Vector<amrex::ParticleReal>& products_mass,
102  const index_type* AMREX_RESTRICT p_mask,
103  const amrex::Vector<index_type>& products_np,
104  const SmartCopy* AMREX_RESTRICT copy_species1,
105  const SmartCopy* AMREX_RESTRICT copy_species2,
106  const index_type* AMREX_RESTRICT p_pair_indices_1,
107  const index_type* AMREX_RESTRICT p_pair_indices_2,
108  const amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight
109  ) const
110  {
111  using namespace amrex::literals;
112 
113  if (n_total_pairs == 0) { return amrex::Vector<int>(m_num_product_species, 0); }
114 
115  // Compute offset array and allocate memory for the produced species
116  amrex::Gpu::DeviceVector<index_type> offsets(n_total_pairs);
117  const auto total = amrex::Scan::ExclusiveSum(n_total_pairs, p_mask, offsets.data());
118  const index_type* AMREX_RESTRICT p_offsets = offsets.dataPtr();
120  for (int i = 0; i < m_num_product_species; i++)
121  {
122  // How many particles of product species i are created.
123  // Factor 2 is here because we currently create one product species at the position of
124  // each source particle of the binary collision. E.g., if a binary collision produces
125  // one electron, we create two electrons, one at the position of each particle that
126  // collided. This allows for exact charge conservation.
127  const index_type num_added = total * m_num_products_host[i] * 2;
128  num_added_vec[i] = static_cast<int>(num_added);
129  tile_products[i]->resize(products_np[i] + num_added);
130  }
131 
132  amrex::ParticleReal* AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
133  amrex::ParticleReal* AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
134  uint64_t* AMREX_RESTRICT idcpu1 = soa_1.m_idcpu;
135  uint64_t* AMREX_RESTRICT idcpu2 = soa_2.m_idcpu;
136 
137  // Create necessary GPU vectors, that will be used in the kernel below
138  amrex::Vector<SoaData_type> soa_products;
139  for (int i = 0; i < m_num_product_species; i++)
140  {
141  soa_products.push_back(tile_products[i]->getParticleTileData());
142  }
143 #ifdef AMREX_USE_GPU
147  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, soa_products.begin(),
148  soa_products.end(),
149  device_soa_products.begin());
150  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, products_np.begin(),
151  products_np.end(),
152  device_products_np.begin());
153  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, products_mass.begin(),
154  products_mass.end(),
155  device_products_mass.begin());
157  SoaData_type* AMREX_RESTRICT soa_products_data = device_soa_products.data();
158  const index_type* AMREX_RESTRICT products_np_data = device_products_np.data();
159  const amrex::ParticleReal* AMREX_RESTRICT products_mass_data = device_products_mass.data();
160 #else
161  SoaData_type* AMREX_RESTRICT soa_products_data = soa_products.data();
162  const index_type* AMREX_RESTRICT products_np_data = products_np.data();
163  const amrex::ParticleReal* AMREX_RESTRICT products_mass_data = products_mass.data();
164 #endif
165 
166  const int t_num_product_species = m_num_product_species;
167  const int* AMREX_RESTRICT p_num_products_device = m_num_products_device.data();
168  const CollisionType t_collision_type = m_collision_type;
169 
170  amrex::ParallelForRNG(n_total_pairs,
171  [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept
172  {
173  if (p_mask[i])
174  {
175  for (int j = 0; j < t_num_product_species; j++)
176  {
177  for (int k = 0; k < p_num_products_device[j]; k++)
178  {
179  // Factor 2 is here because we create one product species at the position
180  // of each source particle
181  const auto product_index = products_np_data[j] +
182  2*(p_offsets[i]*p_num_products_device[j] + k);
183  // Create product particle at position of particle 1
184  copy_species1[j](soa_products_data[j], soa_1, static_cast<int>(p_pair_indices_1[i]),
185  static_cast<int>(product_index), engine);
186  // Create another product particle at position of particle 2
187  copy_species2[j](soa_products_data[j], soa_2, static_cast<int>(p_pair_indices_2[i]),
188  static_cast<int>(product_index + 1), engine);
189 
190  // Set the weight of the new particles to p_pair_reaction_weight[i]/2
191  soa_products_data[j].m_rdata[PIdx::w][product_index] =
192  p_pair_reaction_weight[i]/amrex::ParticleReal(2.);
193  soa_products_data[j].m_rdata[PIdx::w][product_index + 1] =
194  p_pair_reaction_weight[i]/amrex::ParticleReal(2.);
195  }
196  }
197 
198  // Remove p_pair_reaction_weight[i] from the colliding particles' weights
199  amrex::Gpu::Atomic::AddNoRet(&w1[p_pair_indices_1[i]],
200  -p_pair_reaction_weight[i]);
201  amrex::Gpu::Atomic::AddNoRet(&w2[p_pair_indices_2[i]],
202  -p_pair_reaction_weight[i]);
203 
204  // Note: Particle::atomicSetID should also be provided as a standalone helper function in AMReX
205  // to replace the following lambda.
206  auto const atomicSetIdInvalid = [] AMREX_GPU_DEVICE (uint64_t & idcpu)
207  {
208 #if defined(AMREX_USE_OMP)
209 #pragma omp atomic write
211 #else
213  (unsigned long long *)&idcpu,
214  (unsigned long long)amrex::ParticleIdCpus::Invalid
215  );
216 #endif
217  };
218 
219  // If the colliding particle weight decreases to zero, remove particle by
220  // setting its id to invalid
221  if (w1[p_pair_indices_1[i]] <= std::numeric_limits<amrex::ParticleReal>::min())
222  {
223  atomicSetIdInvalid(idcpu1[p_pair_indices_1[i]]);
224 
225  }
226  if (w2[p_pair_indices_2[i]] <= std::numeric_limits<amrex::ParticleReal>::min())
227  {
228  atomicSetIdInvalid(idcpu2[p_pair_indices_2[i]]);
229  }
230 
231  // Initialize the product particles' momentum, using a function depending on the
232  // specific collision type
233  if (t_collision_type == CollisionType::ProtonBoronToAlphasFusion)
234  {
235  const index_type product_start_index = products_np_data[0] + 2*p_offsets[i]*
236  p_num_products_device[0];
237  ProtonBoronFusionInitializeMomentum(soa_1, soa_2, soa_products_data[0],
238  p_pair_indices_1[i], p_pair_indices_2[i],
239  product_start_index, m1, m2, engine);
240  }
241  else if ((t_collision_type == CollisionType::DeuteriumTritiumToNeutronHeliumFusion)
244  {
245  amrex::ParticleReal fusion_energy = 0.0_prt;
247  fusion_energy = 17.5893e6_prt * PhysConst::q_e;
248  }
249  else if (t_collision_type == CollisionType::DeuteriumDeuteriumToProtonTritiumFusion) {
250  fusion_energy = 4.032667e6_prt * PhysConst::q_e;
251  }
252  else if (t_collision_type == CollisionType::DeuteriumDeuteriumToNeutronHeliumFusion) {
253  fusion_energy = 3.268911e6_prt * PhysConst::q_e;
254  }
255  TwoProductFusionInitializeMomentum(soa_1, soa_2,
256  soa_products_data[0], soa_products_data[1],
257  p_pair_indices_1[i], p_pair_indices_2[i],
258  products_np_data[0] + 2*p_offsets[i]*p_num_products_device[0],
259  products_np_data[1] + 2*p_offsets[i]*p_num_products_device[1],
260  m1, m2, products_mass_data[0], products_mass_data[1], fusion_energy, engine);
261  }
262 
263  }
264  });
265 
266  // Initialize the user runtime components
267  for (int i = 0; i < m_num_product_species; i++)
268  {
269  const auto start_index = int(products_np[i]);
270  const auto stop_index = int(products_np[i] + num_added_vec[i]);
272  0, 0,
273  pc_products[i]->getUserRealAttribs(), pc_products[i]->getUserIntAttribs(),
274  pc_products[i]->getParticleComps(), pc_products[i]->getParticleiComps(),
275  pc_products[i]->getUserRealAttribParser(),
276  pc_products[i]->getUserIntAttribParser(),
277 #ifdef WARPX_QED
278  false, // do not initialize QED quantities, since they were initialized
279  // when calling the SmartCopy functors
280  pc_products[i]->get_breit_wheeler_engine_ptr(),
281  pc_products[i]->get_quantum_sync_engine_ptr(),
282 #endif
283  pc_products[i]->getIonizationInitialLevel(),
284  start_index, stop_index);
285  }
286 
288 
289  return num_added_vec;
290  }
291 
292 private:
293  // How many different type of species the collision produces
295  // Vectors of size m_num_product_species storing how many particles of a given species are
296  // produced by a collision event. These vectors are duplicated (one version for host and one
297  // for device) which is necessary with GPUs but redundant on CPU.
301 };
302 
303 
309 {
316 
317 public:
319 
320  NoParticleCreationFunc (const std::string& /*collision_name*/,
321  MultiParticleContainer const * const /*mypc*/) {}
322 
325  const index_type& /*n_total_pairs*/,
326  const SoaData_type& /*soa_1*/, const SoaData_type& /*soa_2*/,
328  ParticleTileType** /*tile_products*/,
329  const amrex::ParticleReal& /*m1*/, const amrex::ParticleReal& /*m2*/,
330  const amrex::Vector<amrex::ParticleReal>& /*products_mass*/,
331  const index_type* /*p_mask*/, const amrex::Vector<index_type>& /*products_np*/,
332  const SmartCopy* /*copy_species1*/, const SmartCopy* /*copy_species2*/,
333  const index_type* /*p_pair_indices_1*/, const index_type* /*p_pair_indices_2*/,
334  const amrex::ParticleReal* /*p_pair_reaction_weight*/
335  ) const
336  {
337  return {};
338  }
339 };
340 
341 #endif // WARPX_PARTICLE_CREATION_FUNC_H_
#define AMREX_RESTRICT
#define AMREX_INLINE
#define AMREX_GPU_DEVICE
CollisionType
Definition: BinaryCollisionUtils.H:17
@ ProtonBoronToAlphasFusion
@ DeuteriumDeuteriumToProtonTritiumFusion
@ DeuteriumDeuteriumToNeutronHeliumFusion
@ DeuteriumTritiumToNeutronHeliumFusion
Definition: MultiParticleContainer.H:66
This class does nothing and is used as second template parameter for binary collisions that do not cr...
Definition: ParticleCreationFunc.H:309
typename WarpXParticleContainer::ParticleType ParticleType
Definition: ParticleCreationFunc.H:310
NoParticleCreationFunc(const std::string &, MultiParticleContainer const *const)
Definition: ParticleCreationFunc.H:320
typename ParticleBins::index_type index_type
Definition: ParticleCreationFunc.H:314
NoParticleCreationFunc()=default
typename WarpXParticleContainer::ParticleTileType ParticleTileType
Definition: ParticleCreationFunc.H:311
typename ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition: ParticleCreationFunc.H:312
typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType SoaData_type
Definition: ParticleCreationFunc.H:315
This functor creates particles produced from a binary collision and sets their initial properties (po...
Definition: ParticleCreationFunc.H:33
CollisionType m_collision_type
Definition: ParticleCreationFunc.H:300
AMREX_INLINE amrex::Vector< int > operator()(const index_type &n_total_pairs, const SoaData_type &soa_1, const SoaData_type &soa_2, const amrex::Vector< WarpXParticleContainer * > &pc_products, ParticleTileType **AMREX_RESTRICT tile_products, const amrex::ParticleReal &m1, const amrex::ParticleReal &m2, const amrex::Vector< amrex::ParticleReal > &products_mass, const index_type *AMREX_RESTRICT p_mask, const amrex::Vector< index_type > &products_np, const SmartCopy *AMREX_RESTRICT copy_species1, const SmartCopy *AMREX_RESTRICT copy_species2, const index_type *AMREX_RESTRICT p_pair_indices_1, const index_type *AMREX_RESTRICT p_pair_indices_2, const amrex::ParticleReal *AMREX_RESTRICT p_pair_reaction_weight) const
operator() of the ParticleCreationFunc class. It creates new particles from binary collisions....
Definition: ParticleCreationFunc.H:95
int m_num_product_species
Definition: ParticleCreationFunc.H:294
amrex::Gpu::DeviceVector< int > m_num_products_device
Definition: ParticleCreationFunc.H:298
typename ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition: ParticleCreationFunc.H:37
typename WarpXParticleContainer::ParticleType ParticleType
Definition: ParticleCreationFunc.H:35
typename ParticleBins::index_type index_type
Definition: ParticleCreationFunc.H:39
typename WarpXParticleContainer::ParticleTileType ParticleTileType
Definition: ParticleCreationFunc.H:36
ParticleCreationFunc()=default
Default constructor of the ParticleCreationFunc class.
typename WarpXParticleContainer::ParticleTileType::ParticleTileDataType SoaData_type
Definition: ParticleCreationFunc.H:40
amrex::Gpu::HostVector< int > m_num_products_host
Definition: ParticleCreationFunc.H:299
unsigned int index_type
T * data() noexcept
iterator begin() noexcept
T * dataPtr() noexcept
ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > ParticleTileType
void DefaultInitializeRuntimeAttributes(PTile &ptile, const int n_external_attr_real, const int n_external_attr_int, const std::vector< std::string > &user_real_attribs, const std::vector< std::string > &user_int_attribs, const std::map< std::string, int > &particle_comps, const std::map< std::string, int > &particle_icomps, const std::vector< amrex::Parser * > &user_real_attrib_parser, const std::vector< amrex::Parser * > &user_int_attrib_parser, const bool do_qed_comps, BreitWheelerEngine *p_bw_engine, QuantumSynchrotronEngine *p_qs_engine, const int ionization_initial_level, int start, int stop)
Default initialize runtime attributes in a tile. This routine does not initialize the first n_externa...
Definition: DefaultInitialization.H:118
typename ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition: ParticleUtils.cpp:33
static constexpr auto q_e
elementary charge [C]
Definition: constant.H:50
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T Exch(T *address, T val) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
void synchronize() noexcept
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
static constexpr HostToDevice hostToDevice
void streamSynchronize() noexcept
constexpr std::uint64_t Invalid
T ExclusiveSum(N n, T const *in, T *out, RetSum a_ret_sum=retSum)
const int[]
AMREX_ATTRIBUTE_FLATTEN_FOR void ParallelForRNG(T n, L const &f) noexcept
i
Definition: check_interp_points_and_weights.py:174
@ w
weight
Definition: NamedComponentParticleContainer.H:33
This is a functor for performing a "smart copy" that works in both host and device code.
Definition: SmartCopy.H:34
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType