WarpX
BinaryCollision.H
Go to the documentation of this file.
1 /* Copyright 2020-2021 Yinjian Zhao, David Grote, Neil Zaim
2  *
3  * This file is part of WarpX.
4  *
5  * License: BSD-3-Clause-LBNL
6  */
7 #ifndef WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
8 #define WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
9 
21 #include "Utils/ParticleUtils.H"
22 #include "Utils/TextMsg.H"
24 #include "WarpX.H"
25 
28 
29 #include <AMReX.H>
30 #include <AMReX_Algorithm.H>
31 #include <AMReX_BLassert.H>
32 #include <AMReX_Config.H>
33 #include <AMReX_DenseBins.H>
34 #include <AMReX_Extension.H>
35 #include <AMReX_Geometry.H>
36 #include <AMReX_GpuAtomic.H>
37 #include <AMReX_GpuContainers.H>
38 #include <AMReX_GpuControl.H>
39 #include <AMReX_GpuDevice.H>
40 #include <AMReX_GpuLaunch.H>
41 #include <AMReX_GpuQualifiers.H>
42 #include <AMReX_LayoutData.H>
43 #include <AMReX_MFIter.H>
44 #include <AMReX_PODVector.H>
45 #include <AMReX_ParmParse.H>
46 #include <AMReX_Particles.H>
47 #include <AMReX_ParticleTile.H>
48 #include <AMReX_Random.H>
49 #include <AMReX_REAL.H>
50 #include <AMReX_Scan.H>
51 #include <AMReX_Utility.H>
52 #include <AMReX_Vector.H>
53 
54 #include <AMReX_BaseFwd.H>
55 
56 #include <cmath>
57 #include <string>
58 
68 template <typename CollisionFunctor,
69  typename CopyTransformFunctor = NoParticleCreationFunc>
70 class BinaryCollision final
71  : public CollisionBase
72 {
73  // Define shortcuts for frequently-used type names
80 
81 public:
89  BinaryCollision (std::string collision_name, MultiParticleContainer const * const mypc)
90  : CollisionBase(collision_name)
91  {
92  if(m_species_names.size() != 2) {
93  WARPX_ABORT_WITH_MESSAGE("Binary collision " + collision_name + " must have exactly two species.");
94  }
95 
96  const CollisionType collision_type = BinaryCollisionUtils::get_collision_type(collision_name, mypc);
97 
99 
100  m_binary_collision_functor = CollisionFunctor(collision_name, mypc, m_isSameSpecies);
101 
102  const amrex::ParmParse pp_collision_name(collision_name);
103  pp_collision_name.queryarr("product_species", m_product_species);
104 
105  // if DSMC the colliding species are also product species
106  // Therefore, we insert the colliding species at the beginning of `m_product_species`
107  if (collision_type == CollisionType::DSMC) {
108  m_product_species.insert( m_product_species.begin(), m_species_names.begin(), m_species_names.end() );
109  }
111 
113  WARPX_ABORT_WITH_MESSAGE( "Binary collision " + collision_name +
114  " does not produce species. Thus, `product_species` should not be specified in the input script." );
115  }
116  m_copy_transform_functor = CopyTransformFunctor(collision_name, mypc);
117  }
118 
119  ~BinaryCollision () override = default;
120 
121  BinaryCollision ( BinaryCollision const &) = default;
123 
126 
134  void doCollisions (amrex::Real cur_time, amrex::Real dt, MultiParticleContainer* mypc) override
135  {
136  amrex::ignore_unused(cur_time);
137 
139  auto& species2 = mypc->GetParticleContainerFromName(m_species_names[1]);
140 
141  // In case of particle creation, create the necessary vectors
142  const int n_product_species = m_product_species.size();
143  amrex::Vector<WarpXParticleContainer*> product_species_vector;
144  amrex::Vector<SmartCopyFactory> copy_factory_species1;
145  amrex::Vector<SmartCopyFactory> copy_factory_species2;
146  amrex::Vector<SmartCopy> copy_species1;
147  amrex::Vector<SmartCopy> copy_species2;
148  for (int i = 0; i < n_product_species; i++)
149  {
150  auto& product = mypc->GetParticleContainerFromName(m_product_species[i]);
151  product.defineAllParticleTiles();
152  product_species_vector.push_back(&product);
153  // Although the copy factories are not explicitly reused past this point, we need to
154  // store them in vectors so that the data that they own, which is used by the smart
155  // copy functors, does not go out of scope at the end of this for loop.
156  copy_factory_species1.push_back(SmartCopyFactory(species1, product));
157  copy_factory_species2.push_back(SmartCopyFactory(species2, product));
158  copy_species1.push_back(copy_factory_species1[i].getSmartCopy());
159  copy_species2.push_back(copy_factory_species2[i].getSmartCopy());
160  }
161 #ifdef AMREX_USE_GPU
162  amrex::Gpu::DeviceVector<SmartCopy> device_copy_species1(n_product_species);
163  amrex::Gpu::DeviceVector<SmartCopy> device_copy_species2(n_product_species);
164  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, copy_species1.begin(),
165  copy_species1.end(), device_copy_species1.begin());
166  amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, copy_species2.begin(),
167  copy_species2.end(), device_copy_species2.begin());
169  auto *copy_species1_data = device_copy_species1.data();
170  auto *copy_species2_data = device_copy_species2.data();
171 #else
172  auto *copy_species1_data = copy_species1.data();
173  auto *copy_species2_data = copy_species2.data();
174 #endif
176  species1.defineAllParticleTiles();
177  if (!m_isSameSpecies) { species2.defineAllParticleTiles(); }
178  }
179 
180  // Enable tiling
181  amrex::MFItInfo info;
182  if (amrex::Gpu::notInLaunchRegion()) { info.EnableTiling(species1.tile_size); }
183 
184  // Loop over refinement levels
185  for (int lev = 0; lev <= species1.finestLevel(); ++lev){
186 
188 
189  // Loop over all grids/tiles at this level
190 #ifdef AMREX_USE_OMP
191  info.SetDynamic(true);
192 #pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
193 #endif
194  for (amrex::MFIter mfi = species1.MakeMFIter(lev, info); mfi.isValid(); ++mfi){
196  {
198  }
199  auto wt = static_cast<amrex::Real>(amrex::second());
200 
201  doCollisionsWithinTile( dt, lev, mfi, species1, species2, product_species_vector,
202  copy_species1_data, copy_species2_data);
203 
205  {
207  wt = static_cast<amrex::Real>(amrex::second()) - wt;
208  amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
209  }
210  }
211 
213  // The fact that there are product species indicates that particles of
214  // the colliding species (`species1` and `species2`) may be removed
215  // (i.e., marked as invalid) in the process of creating new product particles.
216  species1.deleteInvalidParticles();
217  if (!m_isSameSpecies) { species2.deleteInvalidParticles(); }
218  }
219  }
220  }
221 
235  amrex::Real dt, int const lev, amrex::MFIter const& mfi,
236  WarpXParticleContainer& species_1,
237  WarpXParticleContainer& species_2,
238  amrex::Vector<WarpXParticleContainer*> product_species_vector,
239  SmartCopy* copy_species1, SmartCopy* copy_species2)
240  {
241  using namespace ParticleUtils;
242  using namespace amrex::literals;
243 
244  const auto& binary_collision_functor = m_binary_collision_functor.executor();
245  const bool have_product_species = m_have_product_species;
246 
247  // Store product species data in vectors
248  const int n_product_species = m_product_species.size();
249  amrex::Vector<ParticleTileType*> tile_products;
250  amrex::Vector<GetParticlePosition<PIdx>> get_position_products;
251  amrex::Vector<index_type> products_np;
253  constexpr int getpos_offset = 0;
254  for (int i = 0; i < n_product_species; i++)
255  {
256  ParticleTileType& ptile_product = product_species_vector[i]->ParticlesAt(lev, mfi);
257  tile_products.push_back(&ptile_product);
258  get_position_products.push_back(GetParticlePosition<PIdx>(ptile_product,
259  getpos_offset));
260  products_np.push_back(ptile_product.numParticles());
261  products_mass.push_back(product_species_vector[i]->getMass());
262  }
263  auto *tile_products_data = tile_products.data();
264 
265  if ( m_isSameSpecies ) // species_1 == species_2
266  {
267  // Extract particles in the tile that `mfi` points to
268  ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
269 
270  // Find the particles that are in each cell of this tile
271  ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
272 
273  // Loop over cells, and collide the particles in each cell
274 
275  // Extract low-level data
276  auto const n_cells = static_cast<int>(bins_1.numBins());
277  // - Species 1
278  const auto soa_1 = ptile_1.getParticleTileData();
279  index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr();
280  index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr();
281  const amrex::ParticleReal q1 = species_1.getCharge();
282  const amrex::ParticleReal m1 = species_1.getMass();
283  auto get_position_1 = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
284 
285  amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
286 #if defined WARPX_DIM_1D_Z
287  auto dV = geom.CellSize(0);
288 #elif defined WARPX_DIM_XZ
289  auto dV = geom.CellSize(0) * geom.CellSize(1);
290 #elif defined WARPX_DIM_RZ
291  amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
292  const auto lo = lbound(cbx);
293  const auto hi = ubound(cbx);
294  int const nz = hi.y-lo.y+1;
295  auto dr = geom.CellSize(0);
296  auto dz = geom.CellSize(1);
297 #elif defined(WARPX_DIM_3D)
298  auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
299 #endif
300 
301 
302  /*
303  The following calculations are only required when creating product particles
304  */
305  const int n_cells_products = have_product_species ? n_cells: 0;
306  amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
307  index_type* AMREX_RESTRICT p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
308 
309  // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
310  // For a single species, the number of pair in a cell is half the number of particles
311  // in that cell, rounded up to the next higher integer.
312  amrex::ParallelFor( n_cells_products,
313  [=] AMREX_GPU_DEVICE (int i_cell) noexcept
314  {
315  const auto n_part_in_cell = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
316  // Particular case: if there's only 1 particle in a cell, then there's no pair
317  p_n_pairs_in_each_cell[i_cell] = (n_part_in_cell == 1)? 0: (n_part_in_cell+1)/2;
318  }
319  );
320 
321  // Start indices of the pairs in a cell. Will be used for particle creation.
322  amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
323  const index_type n_total_pairs = (n_cells_products == 0) ? 0:
324  amrex::Scan::ExclusiveSum(n_cells_products,
325  p_n_pairs_in_each_cell, pair_offsets.data());
326  index_type* AMREX_RESTRICT p_pair_offsets = pair_offsets.dataPtr();
327 
328  amrex::Gpu::DeviceVector<index_type> n_ind_pairs_in_each_cell(n_cells+1);
329  index_type* AMREX_RESTRICT p_n_ind_pairs_in_each_cell = n_ind_pairs_in_each_cell.dataPtr();
330 
331  amrex::ParallelFor( n_cells+1,
332  [=] AMREX_GPU_DEVICE (int i_cell) noexcept
333  {
334  const auto n_part_in_cell = (i_cell < n_cells)? cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell]: 0;
335  // number of independent collisions in each cell
336  p_n_ind_pairs_in_each_cell[i_cell] = n_part_in_cell/2;
337  }
338  );
339 
340  // start indices of independent collisions.
341  amrex::Gpu::DeviceVector<index_type> coll_offsets(n_cells+1);
342  // number of total independent collision pairs
343  const auto n_independent_pairs = (int) amrex::Scan::ExclusiveSum(n_cells+1,
344  p_n_ind_pairs_in_each_cell, coll_offsets.data(), amrex::Scan::RetSum{true});
345  index_type* AMREX_RESTRICT p_coll_offsets = coll_offsets.dataPtr();
346 
347  // mask: equal to 1 if particle creation occurs for a given pair, 0 otherwise
350  // Will be filled with the index of the first particle of a given pair
351  amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
352  index_type* AMREX_RESTRICT p_pair_indices_1 = pair_indices_1.dataPtr();
353  // Will be filled with the index of the second particle of a given pair
354  amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
355  index_type* AMREX_RESTRICT p_pair_indices_2 = pair_indices_2.dataPtr();
356  // How much weight should be given to the produced particles (and removed from the
357  // reacting particles)
358  amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
359  amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight =
360  pair_reaction_weight.dataPtr();
361  /*
362  End of calculations only required when creating product particles
363  */
364 
365  // Loop over cells
366  amrex::ParallelForRNG( n_cells,
367  [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
368  {
369  // The particles from species1 that are in the cell `i_cell` are
370  // given by the `indices_1[cell_start_1:cell_stop_1]`
371  index_type const cell_start_1 = cell_offsets_1[i_cell];
372  index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
373  index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2;
374 
375  // Do not collide if there is only one particle in the cell
376  if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
377 
378  // shuffle
380  indices_1, cell_start_1, cell_half_1, engine );
381  }
382  );
383 
384  // Loop over independent particle pairs
385  // To speed up binary collisions on GPU, we try to expose as much parallelism
386  // as possible (while avoiding race conditions): Instead of looping with one GPU
387  // thread per cell, we loop with one GPU thread per "independent pairs" (i.e. pairs
388  // that do not touch the same macroparticles, so that there is no race condition),
389  // where the number of independent pairs is determined by the lower number of
390  // macroparticles of either species, within each cell.
391  amrex::ParallelForRNG( n_independent_pairs,
392  [=] AMREX_GPU_DEVICE (int i_coll, amrex::RandomEngine const& engine) noexcept
393  {
394  // to avoid type mismatch errors
395  auto ui_coll = (index_type)i_coll;
396 
397  // Use a bisection algorithm to find the index of the cell in which this pair is located
398  const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
399 
400  // The particles from species1 that are in the cell `i_cell` are
401  // given by the `indices_1[cell_start_1:cell_stop_1]`
402  index_type const cell_start_1 = cell_offsets_1[i_cell];
403  index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
404  index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2;
405 
406  // collision number of the cell
407  const index_type coll_idx = ui_coll - p_coll_offsets[i_cell];
408 
409  // Same but for the pairs
410  index_type const cell_start_pair = have_product_species?
411  p_pair_offsets[i_cell] : 0;
412 
413 #if defined WARPX_DIM_RZ
414  const int ri = (i_cell - i_cell%nz) / nz;
415  auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz;
416 #endif
417  // Call the function in order to perform collisions
418  // If there are product species, mask, p_pair_indices_1/2, and
419  // p_pair_reaction_weight are filled here
420  binary_collision_functor(
421  cell_start_1, cell_half_1,
422  cell_half_1, cell_stop_1,
423  indices_1, indices_1,
424  soa_1, soa_1, get_position_1, get_position_1,
425  q1, q1, m1, m1, dt, dV, coll_idx,
426  cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
427  p_pair_reaction_weight, engine);
428  }
429  );
430 
431  // Create the new product particles and define their initial values
432  // num_added: how many particles of each product species have been created
433  const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
434  soa_1, soa_1,
435  product_species_vector,
436  tile_products_data,
437  m1, m1,
438  products_mass, p_mask, products_np,
439  copy_species1, copy_species2,
440  p_pair_indices_1, p_pair_indices_2,
441  p_pair_reaction_weight);
442 
443  for (int i = 0; i < n_product_species; i++)
444  {
445  setNewParticleIDs(*(tile_products_data[i]), static_cast<int>(products_np[i]), num_added[i]);
446  }
447  }
448  else // species_1 != species_2
449  {
450  // Extract particles in the tile that `mfi` points to
451  ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
452  ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
453 
454  // Find the particles that are in each cell of this tile
455  ParticleBins bins_1 = findParticlesInEachCell( lev, mfi, ptile_1 );
456  ParticleBins bins_2 = findParticlesInEachCell( lev, mfi, ptile_2 );
457 
458  // Loop over cells, and collide the particles in each cell
459 
460  // Extract low-level data
461  auto const n_cells = static_cast<int>(bins_1.numBins());
462  // - Species 1
463  const auto soa_1 = ptile_1.getParticleTileData();
464  index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr();
465  index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr();
466  const amrex::ParticleReal q1 = species_1.getCharge();
467  const amrex::ParticleReal m1 = species_1.getMass();
468  auto get_position_1 = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
469  // - Species 2
470  const auto soa_2 = ptile_2.getParticleTileData();
471  index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr();
472  index_type const* AMREX_RESTRICT cell_offsets_2 = bins_2.offsetsPtr();
473  const amrex::ParticleReal q2 = species_2.getCharge();
474  const amrex::ParticleReal m2 = species_2.getMass();
475  auto get_position_2 = GetParticlePosition<PIdx>(ptile_2, getpos_offset);
476 
477  amrex::Geometry const& geom = WarpX::GetInstance().Geom(lev);
478 #if defined WARPX_DIM_1D_Z
479  auto dV = geom.CellSize(0);
480 #elif defined WARPX_DIM_XZ
481  auto dV = geom.CellSize(0) * geom.CellSize(1);
482 #elif defined WARPX_DIM_RZ
483  amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
484  const auto lo = lbound(cbx);
485  const auto hi = ubound(cbx);
486  const int nz = hi.y-lo.y+1;
487  auto dr = geom.CellSize(0);
488  auto dz = geom.CellSize(1);
489 #elif defined(WARPX_DIM_3D)
490  auto dV = geom.CellSize(0) * geom.CellSize(1) * geom.CellSize(2);
491 #endif
492 
493  /*
494  The following calculations are only required when creating product particles
495  */
496  const int n_cells_products = have_product_species ? n_cells: 0;
497  amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
498  index_type* AMREX_RESTRICT p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
499 
500  // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
501  // For different species, the number of pairs in a cell is the number of particles of
502  // the species that has the most particles in that cell
503  amrex::ParallelFor( n_cells_products,
504  [=] AMREX_GPU_DEVICE (int i_cell) noexcept
505  {
506  const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
507  const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
508  // Particular case: no pair if a species has no particle in that cell
509  if (n_part_in_cell_1 == 0 || n_part_in_cell_2 == 0) {
510  p_n_pairs_in_each_cell[i_cell] = 0;
511  } else {
512  p_n_pairs_in_each_cell[i_cell] =
513  amrex::max(n_part_in_cell_1,n_part_in_cell_2);
514  }
515  }
516  );
517 
518  // Start indices of the pairs in a cell. Will be used for particle creation
519  amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
520  const index_type n_total_pairs = (n_cells_products == 0) ? 0:
521  amrex::Scan::ExclusiveSum(n_cells_products,
522  p_n_pairs_in_each_cell, pair_offsets.data());
523  index_type* AMREX_RESTRICT p_pair_offsets = pair_offsets.dataPtr();
524 
525  amrex::Gpu::DeviceVector<index_type> n_ind_pairs_in_each_cell(n_cells+1);
526  index_type* AMREX_RESTRICT p_n_ind_pairs_in_each_cell = n_ind_pairs_in_each_cell.dataPtr();
527 
528  amrex::ParallelFor( n_cells+1,
529  [=] AMREX_GPU_DEVICE (int i_cell) noexcept
530  {
531  if (i_cell < n_cells)
532  {
533  const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
534  const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
535  p_n_ind_pairs_in_each_cell[i_cell] = amrex::min(n_part_in_cell_1, n_part_in_cell_2);
536  }
537  else
538  {
539  p_n_ind_pairs_in_each_cell[i_cell] = 0;
540  }
541  }
542  );
543 
544  // start indices of independent collisions.
545  amrex::Gpu::DeviceVector<index_type> coll_offsets(n_cells+1);
546  // number of total independent collision pairs
547  const auto n_independent_pairs = (int) amrex::Scan::ExclusiveSum(n_cells+1,
548  p_n_ind_pairs_in_each_cell, coll_offsets.data(), amrex::Scan::RetSum{true});
549  index_type* AMREX_RESTRICT p_coll_offsets = coll_offsets.dataPtr();
550 
551  // mask: equal to 1 if particle creation occurs for a given pair, 0 otherwise
554  // Will be filled with the index of the first particle of a given pair
555  amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
556  index_type* AMREX_RESTRICT p_pair_indices_1 = pair_indices_1.dataPtr();
557  // Will be filled with the index of the second particle of a given pair
558  amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
559  index_type* AMREX_RESTRICT p_pair_indices_2 = pair_indices_2.dataPtr();
560  // How much weight should be given to the produced particles (and removed from the
561  // reacting particles)
562  amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
563  amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight =
564  pair_reaction_weight.dataPtr();
565  /*
566  End of calculations only required when creating product particles
567  */
568 
569 
570  // Loop over cells
571  amrex::ParallelForRNG( n_cells,
572  [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
573  {
574  // The particles from species1 that are in the cell `i_cell` are
575  // given by the `indices_1[cell_start_1:cell_stop_1]`
576  index_type const cell_start_1 = cell_offsets_1[i_cell];
577  index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
578  // Same for species 2
579  index_type const cell_start_2 = cell_offsets_2[i_cell];
580  index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
581 
582  // ux from species1 can be accessed like this:
583  // ux_1[ indices_1[i] ], where i is between
584  // cell_start_1 (inclusive) and cell_start_2 (exclusive)
585 
586  // Do not collide if one species is missing in the cell
587  if ( cell_stop_1 - cell_start_1 < 1 ||
588  cell_stop_2 - cell_start_2 < 1 ) { return; }
589 
590  // shuffle
591  ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
592  ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine);
593  }
594  );
595 
596  // Loop over independent particle pairs
597  // To speed up binary collisions on GPU, we try to expose as much parallelism
598  // as possible (while avoiding race conditions): Instead of looping with one GPU
599  // thread per cell, we loop with one GPU thread per "independent pairs" (i.e. pairs
600  // that do not touch the same macroparticles, so that there is no race condition),
601  // where the number of independent pairs is determined by the lower number of
602  // macroparticles of either species, within each cell.
603  amrex::ParallelForRNG( n_independent_pairs,
604  [=] AMREX_GPU_DEVICE (int i_coll, amrex::RandomEngine const& engine) noexcept
605  {
606  // to avoid type mismatch errors
607  auto ui_coll = (index_type)i_coll;
608 
609  // Use a bisection algorithm to find the index of the cell in which this pair is located
610  const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
611 
612  // The particles from species1 that are in the cell `i_cell` are
613  // given by the `indices_1[cell_start_1:cell_stop_1]`
614  index_type const cell_start_1 = cell_offsets_1[i_cell];
615  index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
616  // Same for species 2
617  index_type const cell_start_2 = cell_offsets_2[i_cell];
618  index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
619 
620  // collision number of the cell
621  const index_type coll_idx = ui_coll - p_coll_offsets[i_cell];
622 
623  // Same but for the pairs
624  index_type const cell_start_pair = have_product_species?
625  p_pair_offsets[i_cell]: 0;
626 
627  // ux from species1 can be accessed like this:
628  // ux_1[ indices_1[i] ], where i is between
629  // cell_start_1 (inclusive) and cell_start_2 (exclusive)
630 
631 #if defined WARPX_DIM_RZ
632  const int ri = (i_cell - i_cell%nz) / nz;
633  auto dV = MathConst::pi*(2.0_prt*ri+1.0_prt)*dr*dr*dz;
634 #endif
635  // Call the function in order to perform collisions
636  // If there are product species, p_mask, p_pair_indices_1/2, and
637  // p_pair_reaction_weight are filled here
638  binary_collision_functor(
639  cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
640  indices_1, indices_2,
641  soa_1, soa_2, get_position_1, get_position_2,
642  q1, q2, m1, m2, dt, dV, coll_idx,
643  cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
644  p_pair_reaction_weight, engine);
645  }
646  );
647 
648  // Create the new product particles and define their initial values
649  // num_added: how many particles of each product species have been created
650  const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
651  soa_1, soa_2,
652  product_species_vector,
653  tile_products_data,
654  m1, m2,
655  products_mass, p_mask, products_np,
656  copy_species1, copy_species2,
657  p_pair_indices_1, p_pair_indices_2,
658  p_pair_reaction_weight);
659 
660  for (int i = 0; i < n_product_species; i++)
661  {
662  setNewParticleIDs(*(tile_products_data[i]), static_cast<int>(products_np[i]), num_added[i]);
663  }
664 
665  } // end if ( m_isSameSpecies)
666 
667  }
668 
669 private:
670 
674  // functor that performs collisions within a cell
675  CollisionFunctor m_binary_collision_functor;
676  // functor that creates new particles and initializes their parameters
677  CopyTransformFunctor m_copy_transform_functor;
678 
679 };
680 
681 #endif // WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
#define AMREX_RESTRICT
#define AMREX_GPU_DEVICE
Array4< int const > mask
CollisionType
Definition: BinaryCollisionUtils.H:17
AMREX_GPU_HOST_DEVICE AMREX_INLINE void ShuffleFisherYates(T_index *array, T_index const is, T_index const ie, amrex::RandomEngine const &engine)
Definition: ShuffleFisherYates.H:20
void setNewParticleIDs(PTile &ptile, amrex::Long old_size, amrex::Long num_added)
Sets the ids of newly created particles to the next values.
Definition: SmartUtils.H:52
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition: TextMsg.H:15
This class performs generic binary collisions.
Definition: BinaryCollision.H:72
bool m_have_product_species
Definition: BinaryCollision.H:672
ParticleBins::index_type index_type
Definition: BinaryCollision.H:79
CopyTransformFunctor m_copy_transform_functor
Definition: BinaryCollision.H:677
bool m_isSameSpecies
Definition: BinaryCollision.H:671
WarpXParticleContainer::ParticleType ParticleType
Definition: BinaryCollision.H:74
BinaryCollision(std::string collision_name, MultiParticleContainer const *const mypc)
Constructor of the BinaryCollision class.
Definition: BinaryCollision.H:89
~BinaryCollision() override=default
void doCollisions(amrex::Real cur_time, amrex::Real dt, MultiParticleContainer *mypc) override
Definition: BinaryCollision.H:134
void doCollisionsWithinTile(amrex::Real dt, int const lev, amrex::MFIter const &mfi, WarpXParticleContainer &species_1, WarpXParticleContainer &species_2, amrex::Vector< WarpXParticleContainer * > product_species_vector, SmartCopy *copy_species1, SmartCopy *copy_species2)
Definition: BinaryCollision.H:234
BinaryCollision & operator=(BinaryCollision const &)=default
amrex::Vector< std::string > m_product_species
Definition: BinaryCollision.H:673
CollisionFunctor m_binary_collision_functor
Definition: BinaryCollision.H:675
BinaryCollision(BinaryCollision const &)=default
BinaryCollision(BinaryCollision &&)=delete
Definition: CollisionBase.H:18
amrex::Vector< std::string > m_species_names
Definition: CollisionBase.H:36
Definition: MultiParticleContainer.H:66
WarpXParticleContainer & GetParticleContainerFromName(const std::string &name) const
Definition: MultiParticleContainer.cpp:392
This class does nothing and is used as second template parameter for binary collisions that do not cr...
Definition: ParticleCreationFunc.H:309
A factory for creating SmartCopy functors.
Definition: SmartCopy.H:133
static WarpX & GetInstance()
Definition: WarpX.cpp:239
static amrex::LayoutData< amrex::Real > * getCosts(int lev)
Definition: WarpX.cpp:3045
static short load_balance_costs_update_algo
Definition: WarpX.H:173
Definition: WarpXParticleContainer.H:111
void defineAllParticleTiles() noexcept
Definition: WarpXParticleContainer.cpp:1518
amrex::ParticleReal getCharge() const
Definition: WarpXParticleContainer.H:371
amrex::ParticleReal getMass() const
Definition: WarpXParticleContainer.H:373
const Vector< Geometry > & Geom() const noexcept
const Real * CellSize() const noexcept
index_type * offsetsPtr() noexcept
index_type * permutationPtr() noexcept
unsigned int index_type
Long numBins() const noexcept
AMREX_GPU_HOST_DEVICE static constexpr AMREX_FORCE_INLINE IntVect TheZeroVector() noexcept
Box tilebox() const noexcept
T * data() noexcept
iterator begin() noexcept
T * dataPtr() noexcept
int queryarr(const char *name, std::vector< int > &ref, int start_ix=FIRST, int num_val=ALL) const
ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > ParticleTileType
const ParticleTileType & ParticlesAt(int lev, int grid, int tile) const
CollisionType get_collision_type(const std::string &collision_name, MultiParticleContainer const *const mypc)
Definition: BinaryCollisionUtils.cpp:20
Definition: ParticleUtils.cpp:26
ParticleBins findParticlesInEachCell(int lev, MFIter const &mfi, ParticleTileType &ptile)
Find the particles and count the particles that are in each cell. More specifically this function ret...
Definition: ParticleUtils.cpp:40
void synchronize() noexcept
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
static constexpr HostToDevice hostToDevice
void streamSynchronize() noexcept
bool notInLaunchRegion() noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void Add(T *const sum, T const value) noexcept
T ExclusiveSum(N n, T const *in, T *out, RetSum a_ret_sum=retSum)
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
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 ubound(Array4< T > const &a) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE T bisect(T lo, T hi, F f, T tol=1e-12, int max_iter=100)
double second() noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 lbound(Array4< T > const &a) noexcept
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
const int[]
std::enable_if_t< std::is_integral_v< T > > ParallelFor(TypeList< CTOs... >, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
AMREX_ATTRIBUTE_FLATTEN_FOR void ParallelForRNG(T n, L const &f) noexcept
i
Definition: check_interp_points_and_weights.py:174
int dz
Definition: compute_domain.py:36
float dt
Definition: stencil.py:442
value
Definition: updateAMReX.py:141
string species1
Definition: video_yt.py:35
@ Timers
load balance according to in-code timer-based weights (i.e., with costs)
Definition: WarpXAlgorithmSelection.H:147
This is a functor for performing a "smart copy" that works in both host and device code.
Definition: SmartCopy.H:34
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int const * dataPtr() const noexcept
MFItInfo & EnableTiling(const IntVect &ts=FabArrayBase::mfiter_tile_size) noexcept
MFItInfo & SetDynamic(bool f) noexcept
ParticleTileData< StorageParticleType, NArrayReal, NArrayInt > ParticleTileDataType
ParticleTileDataType getParticleTileData()
int numParticles() const