WarpX
Loading...
Searching...
No Matches
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
10#include "Particles/Algorithms/KineticEnergy.H"
11#include "Particles/Collision/BinaryCollision/Bremsstrahlung/BremsstrahlungFunc.H"
12#include "Particles/Collision/BinaryCollision/Coulomb/PairWiseCoulombCollisionFunc.H"
13#include "Particles/Collision/BinaryCollision/DSMC/DSMCFunc.H"
14#include "Particles/Collision/BinaryCollision/NuclearFusion/NuclearFusionFunc.H"
15#include "Particles/Collision/BinaryCollision/LinearBreitWheeler/LinearBreitWheelerCollisionFunc.H"
16#include "Particles/Collision/BinaryCollision/LinearCompton/LinearComptonCollisionFunc.H"
17#include "Particles/Collision/BinaryCollision/ParticleCreationFunc.H"
18#include "Particles/Collision/BinaryCollision/ShuffleFisherYates.H"
19#include "Particles/Collision/CollisionBase.H"
20#include "Particles/ParticleCreation/SmartCopy.H"
21#include "Particles/ParticleCreation/SmartUtils.H"
22#include "Particles/Pusher/GetAndSetPosition.H"
23#include "Particles/MultiParticleContainer.H"
24#include "Particles/WarpXParticleContainer.H"
25#include "Utils/ParticleUtils.H"
26#include "Utils/TextMsg.H"
27#include "Utils/WarpXAlgorithmSelection.H"
28#include "WarpX.H"
29
30#include "Particles/MultiParticleContainer_fwd.H"
31#include "Particles/WarpXParticleContainer_fwd.H"
32
33#include <AMReX.H>
34#include <AMReX_Algorithm.H>
35#include <AMReX_BLassert.H>
36#include <AMReX_Config.H>
37#include <AMReX_DenseBins.H>
38#include <AMReX_Extension.H>
39#include <AMReX_Geometry.H>
40#include <AMReX_GpuAtomic.H>
41#include <AMReX_GpuContainers.H>
42#include <AMReX_GpuControl.H>
43#include <AMReX_GpuDevice.H>
44#include <AMReX_GpuLaunch.H>
45#include <AMReX_GpuQualifiers.H>
46#include <AMReX_LayoutData.H>
47#include <AMReX_MFIter.H>
48#include <AMReX_PODVector.H>
49#include <AMReX_ParmParse.H>
50#include <AMReX_Particles.H>
51#include <AMReX_ParticleTile.H>
52#include <AMReX_Random.H>
53#include <AMReX_REAL.H>
54#include <AMReX_Scan.H>
55#include <AMReX_Utility.H>
56#include <AMReX_Vector.H>
57
58#include <AMReX_BaseFwd.H>
59
60#include <cmath>
61#include <string>
62
72template <typename CollisionFunctor,
73 typename CopyTransformFunctor = NoParticleCreationFunc>
74class BinaryCollision final
75 : public CollisionBase
76{
77 // Define shortcuts for frequently-used type names
78 using ParticleType = WarpXParticleContainer::ParticleType;
79 using ParticleTileType = WarpXParticleContainer::ParticleTileType;
80 using ParticleTileDataType = ParticleTileType::ParticleTileDataType;
81 using ParticleBins = amrex::DenseBins<ParticleTileDataType>;
82 using index_type = ParticleBins::index_type;
83
84public:
92 BinaryCollision (const std::string& collision_name, MultiParticleContainer const * const mypc)
93 : CollisionBase(collision_name)
94 {
95 using namespace amrex::literals;
96 if(m_species_names.size() != 2) {
97 WARPX_ABORT_WITH_MESSAGE("Binary collision " + collision_name + " must have exactly two species.");
98 }
99
100 const CollisionType collision_type = BinaryCollisionUtils::get_collision_type(collision_name, mypc);
101
102 m_isSameSpecies = (m_species_names[0] == m_species_names[1]);
103
104 m_binary_collision_functor = CollisionFunctor(collision_name, mypc, m_isSameSpecies);
105
106 m_use_global_debye_length = m_binary_collision_functor.use_global_debye_length();
107
108 const amrex::ParmParse pp_collision_name(collision_name);
109 pp_collision_name.queryarr("product_species", m_product_species);
110
111
112 if (collision_type == CollisionType::PairwiseCoulomb) {
113 // Input parameter set for all pairwise Coulomb collisions
114 const amrex::ParmParse pp_collisions("collisions");
115 pp_collisions.query("correct_energy_momentum", m_correct_energy_momentum);
116 pp_collisions.query("energy_correction_sort_by_weight", m_energy_correction_sort_by_weight);
117 pp_collisions.query("np_warning_threshold", m_np_warning_threshold);
118 pp_collisions.query("beta_weight_exponent", m_beta_weight_exponent);
119 pp_collisions.query("energy_fraction", m_energy_fraction);
120
121 // Input parameter set for this collision
122 pp_collision_name.query("correct_energy_momentum", m_correct_energy_momentum);
123 pp_collision_name.query("energy_correction_sort_by_weight", m_energy_correction_sort_by_weight);
124 pp_collision_name.query("np_warning_threshold", m_np_warning_threshold);
125 pp_collision_name.query("beta_weight_exponent", m_beta_weight_exponent);
126 pp_collision_name.query("energy_fraction", m_energy_fraction);
127 }
128
129 // If DSMC the colliding species are also product species.
130 // Therefore, we insert the colliding species at the beginning of `m_product_species`.
131 if (collision_type == CollisionType::DSMC) {
132 // If the scattering process is ionization ensure that the
133 // explicitly specified "target" species, i.e., the species that
134 // undergoes ionization, is second in the species list for this
135 // collision set. The reason for this is that during the collision
136 // operation, an outgoing particle of the first species type will
137 // be created.
138 std::string target_species;
139 pp_collision_name.query("ionization_target_species", target_species);
140 if (!target_species.empty()) {
141 if (m_species_names[0] == target_species) {
142 std::swap(m_species_names[0], m_species_names[1]);
143 } else if (m_species_names[1] != target_species) {
144 WARPX_ABORT_WITH_MESSAGE("DSMC: Ionization target species, " + target_species + " must be one of the colliding species.");
145 }
146 }
147
148 m_product_species.insert( m_product_species.begin(), m_species_names.begin(), m_species_names.end() );
149 // Note that, if a species is both a colliding species and a product species
150 // (e.g., e- in e- + H -> 2 e- + H+) then it will be listed twice in `m_product_species`.
151 // (e.g., m_product_species = [e-, H, e-, H+])
152 // The code for ionization in `SplitAndScatterFunc` handles this case correctly.
153 }
154 m_have_product_species = !m_product_species.empty();
155
156 if ((std::is_same_v<CopyTransformFunctor, NoParticleCreationFunc>) && (m_have_product_species)) {
157 WARPX_ABORT_WITH_MESSAGE( "Binary collision " + collision_name +
158 " does not produce species. Thus, `product_species` should not be specified in the input script." );
159 }
160 m_copy_transform_functor = CopyTransformFunctor(collision_name, mypc);
161 }
162
163 ~BinaryCollision () override = default;
164
165 BinaryCollision ( BinaryCollision const &) = default;
166 BinaryCollision& operator= ( BinaryCollision const & ) = default;
167
168 BinaryCollision ( BinaryCollision&& ) = delete;
169 BinaryCollision& operator= ( BinaryCollision&& ) = delete;
170
178 void doCollisions (amrex::Real cur_time, amrex::Real dt, MultiParticleContainer* mypc) override
179 {
180 amrex::ignore_unused(cur_time);
181
182 auto& species1 = mypc->GetParticleContainerFromName(m_species_names[0]);
183 auto& species2 = mypc->GetParticleContainerFromName(m_species_names[1]);
184
185 // In case of particle creation, create the necessary vectors
186 const int n_product_species = m_product_species.size();
187 amrex::Vector<WarpXParticleContainer*> product_species_vector;
188 amrex::Vector<SmartCopyFactory> copy_factory_species1;
189 amrex::Vector<SmartCopyFactory> copy_factory_species2;
190 amrex::Vector<SmartCopy> copy_species1;
191 amrex::Vector<SmartCopy> copy_species2;
192 for (int i = 0; i < n_product_species; i++)
193 {
194 auto& product = mypc->GetParticleContainerFromName(m_product_species[i]);
195 product.defineAllParticleTiles();
196 product_species_vector.push_back(&product);
197 // Although the copy factories are not explicitly reused past this point, we need to
198 // store them in vectors so that the data that they own, which is used by the smart
199 // copy functors, does not go out of scope at the end of this for loop.
200 copy_factory_species1.push_back(SmartCopyFactory(species1, product));
201 copy_factory_species2.push_back(SmartCopyFactory(species2, product));
202 copy_species1.push_back(copy_factory_species1[i].getSmartCopy());
203 copy_species2.push_back(copy_factory_species2[i].getSmartCopy());
204 }
205#ifdef AMREX_USE_GPU
206 amrex::Gpu::DeviceVector<SmartCopy> device_copy_species1(n_product_species);
207 amrex::Gpu::DeviceVector<SmartCopy> device_copy_species2(n_product_species);
208 amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, copy_species1.begin(),
209 copy_species1.end(), device_copy_species1.begin());
210 amrex::Gpu::copyAsync(amrex::Gpu::hostToDevice, copy_species2.begin(),
211 copy_species2.end(), device_copy_species2.begin());
212 amrex::Gpu::streamSynchronize();
213 auto *copy_species1_data = device_copy_species1.data();
214 auto *copy_species2_data = device_copy_species2.data();
215#else
216 auto *copy_species1_data = copy_species1.data();
217 auto *copy_species2_data = copy_species2.data();
218#endif
219 if (m_have_product_species){
220 species1.defineAllParticleTiles();
221 if (!m_isSameSpecies) { species2.defineAllParticleTiles(); }
222 }
223
224 // Enable tiling
225 amrex::MFItInfo info;
226 if (amrex::Gpu::notInLaunchRegion()) { info.EnableTiling(species1.tile_size); }
227
228 // Loop over refinement levels
229 for (int lev = 0; lev <= species1.finestLevel(); ++lev){
230
231 amrex::LayoutData<amrex::Real>* cost = WarpX::getCosts(lev);
232
233 // Loop over all grids/tiles at this level
234#ifdef AMREX_USE_OMP
235 info.SetDynamic(true);
236#pragma omp parallel if (amrex::Gpu::notInLaunchRegion())
237#endif
238 for (amrex::MFIter mfi = species1.MakeMFIter(lev, info); mfi.isValid(); ++mfi){
239 if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
240 {
241 amrex::Gpu::synchronize();
242 }
243 auto wt = static_cast<amrex::Real>(amrex::second());
244
245 doCollisionsWithinTile( dt, lev, mfi, species1, species2, product_species_vector,
246 copy_species1_data, copy_species2_data);
247
248 if (cost && WarpX::load_balance_costs_update_algo == LoadBalanceCostsUpdateAlgo::Timers)
249 {
250 amrex::Gpu::synchronize();
251 wt = static_cast<amrex::Real>(amrex::second()) - wt;
252 amrex::HostDevice::Atomic::Add( &(*cost)[mfi.index()], wt);
253 }
254 }
255
256 if (m_have_product_species) {
257 // The fact that there are product species indicates that particles of
258 // the colliding species (`species1` and `species2`) may be removed
259 // (i.e., marked as invalid) in the process of creating new product particles.
260 species1.deleteInvalidParticles();
261 if (!m_isSameSpecies) { species2.deleteInvalidParticles(); }
262 }
263 }
264 }
265
278 void doCollisionsWithinTile (
279 amrex::Real dt, int const lev, amrex::MFIter const& mfi,
280 WarpXParticleContainer& species_1,
281 WarpXParticleContainer& species_2,
282 amrex::Vector<WarpXParticleContainer*> product_species_vector,
283 SmartCopy* copy_species1, SmartCopy* copy_species2)
284 {
285 using namespace ParticleUtils;
286 using namespace amrex::literals;
287
288 ABLASTR_PROFILE("BinaryCollision::doCollisionsWithinTile");
289
290 const auto& binary_collision_functor = m_binary_collision_functor.executor();
291 const bool have_product_species = m_have_product_species;
292
293 // Store product species data in vectors
294 const int n_product_species = m_product_species.size();
295 amrex::Vector<ParticleTileType*> tile_products;
296 amrex::Vector<GetParticlePosition<PIdx>> get_position_products;
297 amrex::Vector<index_type> products_np;
298 amrex::Vector<amrex::ParticleReal> products_mass;
299 constexpr int getpos_offset = 0;
300 for (int i = 0; i < n_product_species; i++)
301 {
302 ParticleTileType& ptile_product = product_species_vector[i]->ParticlesAt(lev, mfi);
303 tile_products.push_back(&ptile_product);
304 get_position_products.push_back(GetParticlePosition<PIdx>(ptile_product,
305 getpos_offset));
306 products_np.push_back(ptile_product.numParticles());
307 products_mass.push_back(product_species_vector[i]->getMass());
308 }
309 auto *tile_products_data = tile_products.data();
310
311 bool const use_global_debye_length = m_use_global_debye_length;
312 amrex::Real * global_debye_length_data = nullptr;
313 if (use_global_debye_length) {
314 WarpX & warpx = WarpX::GetInstance();
315 amrex::MultiFab & global_debye_length = *warpx.m_fields.get(warpx::fields::FieldType::global_debye_length, lev);
316 amrex::FArrayBox & global_debye_length_fab = global_debye_length[mfi];
317 global_debye_length_data = global_debye_length_fab.dataPtr();
318 }
319
320 amrex::Geometry const& geom_lev = WarpX::GetInstance().Geom(lev);
321 // dV is level-specific: cell volume at this refinement level (smaller on fine levels).
322 amrex::ParticleReal const dV = AMREX_D_TERM(geom_lev.CellSize(0), *geom_lev.CellSize(1), *geom_lev.CellSize(2));
323#if defined(WARPX_DIM_RZ) || defined(WARPX_DIM_RCYLINDER) || defined(WARPX_DIM_RSPHERE)
324 amrex::Box const& cbx = mfi.tilebox(amrex::IntVect::TheZeroVector()); //Cell-centered box
325#if defined(WARPX_DIM_RZ)
326 auto const lo = lbound(cbx);
327 auto const hi = ubound(cbx);
328 int const nr = hi.x - lo.x + 1;
329#endif
330 amrex::XDim3 const xyzmin = WarpX::LowerCorner(cbx, lev, 0._rt);
331 amrex::Real const rmin = xyzmin.x;
332 auto const dr = geom_lev.CellSize(0);
333#endif
334
335 auto volume_factor = [=] AMREX_GPU_DEVICE(int i_cell) noexcept {
336#if defined(WARPX_DIM_RZ)
337 // Return the radial factor for the volume element, dV
338 // DenseBins uses x-fastest ordering: i_cell = iz * nr + ir
339 int const ri = i_cell % nr;
340 // rr is radius at the cell center
341 amrex::ParticleReal const rr = rmin + (ri + 0.5_prt)*dr;
342 return 2.0_prt*static_cast<amrex::ParticleReal>(MathConst::pi)*rr;
343#elif defined(WARPX_DIM_RCYLINDER)
344 int const ri = i_cell;
345 // rr is radius at the cell center
346 amrex::ParticleReal const rr = rmin + (ri + 0.5_prt)*dr;
347 return 2.0_prt*static_cast<amrex::ParticleReal>(MathConst::pi)*rr;
348#elif defined(WARPX_DIM_RSPHERE)
349 // This needs double checking
350 int const ri = i_cell;
351 // rr is radius at the cell center
352 amrex::ParticleReal const rr = rmin + (ri + 0.5_prt)*dr;
353 return 4.0_prt*static_cast<amrex::ParticleReal>(MathConst::pi)*rr*rr;
354#else
355 // No factor is needed for Cartesian
356 amrex::ignore_unused(i_cell);
357 return 1.0_prt;
358#endif
359 };
360
361 if ( m_isSameSpecies ) // species_1 == species_2
362 {
363 // Extract particles in the tile that `mfi` points to
364 ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
365
366 // Find the particles that are in each cell of this tile
367 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findParticlesInEachCell", prof_findParticlesInEachCell);
368 ParticleBins bins_1 = findParticlesInEachCell( geom_lev, mfi, ptile_1 );
369 ABLASTR_PROFILE_VAR_STOP(prof_findParticlesInEachCell);
370
371 // Loop over cells, and collide the particles in each cell
372
373 // Extract low-level data
374 auto const n_cells = static_cast<int>(bins_1.numBins());
375 // - Species 1
376 auto np1 = ptile_1.numParticles();
377 const auto soa_1 = ptile_1.getParticleTileData();
378 index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr();
379 index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr();
380 index_type const* AMREX_RESTRICT bins_1_ptr = bins_1.binsPtr();
381 const amrex::ParticleReal q1 = species_1.getCharge();
382 const amrex::ParticleReal m1 = species_1.getMass();
383 auto get_position_1 = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
384
385 /*
386 The following calculations are only required when creating product particles
387 */
388 const int n_cells_products = have_product_species ? n_cells: 0;
389 amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
390 index_type* AMREX_RESTRICT p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
391
392 // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
393 // For a single species, the number of pair in a cell is half the number of particles
394 // in that cell, rounded up to the next higher integer.
395 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::computeNumberOfPairs", prof_computeNumberOfPairs);
396 amrex::ParallelFor( n_cells_products,
397 [=] AMREX_GPU_DEVICE (int i_cell) noexcept
398 {
399 const auto n_part_in_cell = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
400 // Particular case: if there's only 1 particle in a cell, then there's no pair
401 p_n_pairs_in_each_cell[i_cell] = (n_part_in_cell == 1)? 0: (n_part_in_cell+1)/2;
402 }
403 );
404
405 // Start indices of the pairs in a cell. Will be used for particle creation.
406 amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
407 const index_type n_total_pairs = (n_cells_products == 0) ? 0:
408 amrex::Scan::ExclusiveSum(n_cells_products,
409 p_n_pairs_in_each_cell, pair_offsets.data());
410 index_type* AMREX_RESTRICT p_pair_offsets = pair_offsets.dataPtr();
411
412 amrex::Gpu::DeviceVector<index_type> n_ind_pairs_in_each_cell(n_cells+1);
413 index_type* AMREX_RESTRICT p_n_ind_pairs_in_each_cell = n_ind_pairs_in_each_cell.dataPtr();
414
415 amrex::ParallelFor( n_cells+1,
416 [=] AMREX_GPU_DEVICE (int i_cell) noexcept
417 {
418 const auto n_part_in_cell = (i_cell < n_cells)? cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell]: 0;
419 // number of independent collisions in each cell
420 p_n_ind_pairs_in_each_cell[i_cell] = n_part_in_cell/2;
421 }
422 );
423
424 // start indices of independent collisions.
425 amrex::Gpu::DeviceVector<index_type> coll_offsets(n_cells+1);
426 // number of total independent collision pairs
427 const auto n_independent_pairs = (int) amrex::Scan::ExclusiveSum(n_cells+1,
428 p_n_ind_pairs_in_each_cell, coll_offsets.data(), amrex::Scan::RetSum{true});
429 index_type* AMREX_RESTRICT p_coll_offsets = coll_offsets.dataPtr();
430 ABLASTR_PROFILE_VAR_STOP(prof_computeNumberOfPairs);
431
432 // mask: for DSMC, the index+1 of the selected scattering process (0 = no collision);
433 // for other collision types, 1 if particle creation occurs, 0 otherwise
434 amrex::Gpu::DeviceVector<index_type> mask(n_total_pairs, index_type(0));
435 index_type* AMREX_RESTRICT p_mask = mask.dataPtr();
436 // Will be filled with the index of the first particle of a given pair
437 amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
438 index_type* AMREX_RESTRICT p_pair_indices_1 = pair_indices_1.dataPtr();
439 // Will be filled with the index of the second particle of a given pair
440 amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
441 index_type* AMREX_RESTRICT p_pair_indices_2 = pair_indices_2.dataPtr();
442 // How much weight should be given to the produced particles (and removed from the
443 // reacting particles)
444 amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
445 amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight =
446 pair_reaction_weight.dataPtr();
447 // Extra data needed when products are created
448 int const n_product_data = (binary_collision_functor.m_need_product_data ? n_total_pairs : 0);
449 amrex::Gpu::DeviceVector<amrex::ParticleReal> product_data(n_product_data);
450 amrex::ParticleReal* AMREX_RESTRICT p_product_data = product_data.dataPtr();
451 /*
452 End of calculations only required when creating product particles
453 */
454
455 amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
456 amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
457 amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
458 amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
459
460 // create vectors to store density and temperature on cell level
461 amrex::Gpu::DeviceVector<amrex::ParticleReal> n1_vec;
462 amrex::Gpu::DeviceVector<amrex::ParticleReal> T1_vec;
463 amrex::Gpu::DeviceVector<amrex::ParticleReal> vx1_vec;
464 amrex::Gpu::DeviceVector<amrex::ParticleReal> vy1_vec;
465 amrex::Gpu::DeviceVector<amrex::ParticleReal> vz1_vec;
466 amrex::Gpu::DeviceVector<amrex::ParticleReal> vs1_vec;
467
468 if (binary_collision_functor.m_computeSpeciesDensities) {
469 n1_vec.resize(n_cells, 0.0_prt);
470 }
471 if (binary_collision_functor.m_computeSpeciesTemperatures) {
472 T1_vec.resize(n_cells, 0.0_prt);
473 vx1_vec.resize(n_cells, 0.0_prt);
474 vy1_vec.resize(n_cells, 0.0_prt);
475 vz1_vec.resize(n_cells, 0.0_prt);
476 vs1_vec.resize(n_cells, 0.0_prt);
477 }
478 amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr();
479 amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr();
480 amrex::ParticleReal* AMREX_RESTRICT vx1_in_each_cell = vx1_vec.dataPtr();
481 amrex::ParticleReal* AMREX_RESTRICT vy1_in_each_cell = vy1_vec.dataPtr();
482 amrex::ParticleReal* AMREX_RESTRICT vz1_in_each_cell = vz1_vec.dataPtr();
483 amrex::ParticleReal* AMREX_RESTRICT vs1_in_each_cell = vs1_vec.dataPtr();
484
485 // Create vectors to store energy and momentum in each cell.
486 // This is used to correct the energy and momentum in the
487 // cell after the collisions.
488 amrex::Gpu::DeviceVector<amrex::ParticleReal> ww_weighted_sum_vec;
489 amrex::Gpu::DeviceVector<amrex::ParticleReal> KE_vec;
490 amrex::Gpu::DeviceVector<amrex::ParticleReal> px_vec;
491 amrex::Gpu::DeviceVector<amrex::ParticleReal> py_vec;
492 amrex::Gpu::DeviceVector<amrex::ParticleReal> pz_vec;
493 if (m_correct_energy_momentum) {
494 ww_weighted_sum_vec.resize(n_cells, 0.0_prt);
495 KE_vec.resize(n_cells, 0.0_prt);
496 px_vec.resize(n_cells, 0.0_prt);
497 py_vec.resize(n_cells, 0.0_prt);
498 pz_vec.resize(n_cells, 0.0_prt);
499 }
500 amrex::ParticleReal* AMREX_RESTRICT ww_weighted_sum_in_each_cell = ww_weighted_sum_vec.dataPtr();
501 amrex::ParticleReal* AMREX_RESTRICT KE_in_each_cell = KE_vec.dataPtr();
502 amrex::ParticleReal* AMREX_RESTRICT px_in_each_cell = px_vec.dataPtr();
503 amrex::ParticleReal* AMREX_RESTRICT py_in_each_cell = py_vec.dataPtr();
504 amrex::ParticleReal* AMREX_RESTRICT pz_in_each_cell = pz_vec.dataPtr();
505
506 // Save the momentum before the collision since sometimes the correction fails
507 // and the only solution is to restore the pre-collision values.
508 // The implicit solver already has ux_n etc so use those if available,
509 // otherwise create temporaries.
510 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1x_before;
511 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1y_before;
512 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1z_before;
513 amrex::ParticleReal* AMREX_RESTRICT u1x_before_ptr = nullptr;
514 amrex::ParticleReal* AMREX_RESTRICT u1y_before_ptr = nullptr;
515 amrex::ParticleReal* AMREX_RESTRICT u1z_before_ptr = nullptr;
516 if (m_correct_energy_momentum) {
517 // Check if ux_n was added.
518 std::vector<std::string> const & real_names1 = species_1.GetRealSoANames();
519 auto const pos1 = std::find(real_names1.begin(), real_names1.end(), "ux_n");
520 if (pos1 != real_names1.end()) {
521 int const n_builtin_real = PIdx::nattribs;
522 auto const u1x_ni = static_cast<int>(std::distance(real_names1.begin(), pos1));
523 int const u1x_runtime_ni = u1x_ni - n_builtin_real;
524 u1x_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni];
525 u1y_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 1];
526 u1z_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 2];
527 } else {
528 u1x_before.resize(np1);
529 u1y_before.resize(np1);
530 u1z_before.resize(np1);
531 u1x_before_ptr = u1x_before.dataPtr();
532 u1y_before_ptr = u1y_before.dataPtr();
533 u1z_before_ptr = u1z_before.dataPtr();
534 }
535
536 }
537
538 amrex::ParticleReal const beta_weight_exponent = m_beta_weight_exponent;
539
540 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures", prof_findDensityTemperatures);
541
542 if (m_correct_energy_momentum ||
543 binary_collision_functor.m_computeSpeciesDensities ||
544 binary_collision_functor.m_computeSpeciesTemperatures) {
545
546 bool const correct_energy_momentum = m_correct_energy_momentum;
547
548 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::atomics", prof_findDensityTemperatures_atomics);
549 // Loop over particles and compute quantities needed for energy conservation and the collsion parameters
550 amrex::ParallelFor( np1,
551 [=] AMREX_GPU_DEVICE (int ip) noexcept
552 {
553 if (correct_energy_momentum) {
554 u1x_before_ptr[ip] = u1x[ip];
555 u1y_before_ptr[ip] = u1y[ip];
556 u1z_before_ptr[ip] = u1z[ip];
557
558 const int i_cell = bins_1_ptr[ip];
559
560 // Compute the local energy and momentum.
561 amrex::ParticleReal const w1_scaled = std::pow(w1[ip], beta_weight_exponent);
562 amrex::Gpu::Atomic::AddNoRet(&ww_weighted_sum_in_each_cell[i_cell], w1_scaled*m1);
563 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], w1[ip]*u1x[ip]*m1);
564 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], w1[ip]*u1y[ip]*m1);
565 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], w1[ip]*u1z[ip]*m1);
566 const amrex::ParticleReal KE = Algorithms::KineticEnergy(u1x[ip], u1y[ip], u1z[ip], m1);
567 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], w1[ip]*KE);
568 }
569
570 // compute local density [1/m^3]
571 if (binary_collision_functor.m_computeSpeciesDensities) {
572 amrex::Gpu::Atomic::AddNoRet(&n1_in_each_cell[bins_1_ptr[ip]],
573 w1[ip]/(dV*volume_factor(bins_1_ptr[ip])));
574 }
575
576 // compute local temperature [Joules]
577 if (binary_collision_functor.m_computeSpeciesTemperatures) {
578 const amrex::ParticleReal us = u1x[ip]*u1x[ip] + u1y[ip]*u1y[ip] + u1z[ip]*u1z[ip];
579 auto constexpr inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
580 const amrex::ParticleReal gm = std::sqrt( 1.0_prt + us*inv_c2);
581 amrex::Gpu::Atomic::AddNoRet(&T1_in_each_cell[bins_1_ptr[ip]], w1[ip]);
582 amrex::Gpu::Atomic::AddNoRet(&vx1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1x[ip]/gm);
583 amrex::Gpu::Atomic::AddNoRet(&vy1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1y[ip]/gm);
584 amrex::Gpu::Atomic::AddNoRet(&vz1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1z[ip]/gm);
585 amrex::Gpu::Atomic::AddNoRet(&vs1_in_each_cell[bins_1_ptr[ip]], w1[ip]*us/gm/gm);
586 }
587 }
588 );
589 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_atomics);
590
591 }
592
593 // Finish temperature calculation
594 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::finishTemperature", prof_findDensityTemperatures_finish);
595 amrex::ParallelFor( n_cells, [=] AMREX_GPU_DEVICE (int i_cell) noexcept
596 {
597 // The particles from species1 that are in the cell `i_cell` are
598 // given by the `indices_1[cell_start_1:cell_stop_1]`
599 index_type const cell_start_1 = cell_offsets_1[i_cell];
600 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
601
602 // Do not need if there is only one particle in the cell
603 if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
604
605 // finish temperature calculation if needed
606 if (binary_collision_functor.m_computeSpeciesTemperatures) {
607 const amrex::ParticleReal invsum = 1._prt/T1_in_each_cell[i_cell];
608 auto vx1 = vx1_in_each_cell[i_cell] * invsum;
609 auto vy1 = vy1_in_each_cell[i_cell] * invsum;
610 auto vz1 = vz1_in_each_cell[i_cell] * invsum;
611 auto vs1 = vs1_in_each_cell[i_cell] * invsum;
612
613 T1_in_each_cell[i_cell] = m1/(3._prt)*(vs1 -(vx1*vx1+vy1*vy1+vz1*vz1));
614 }
615 }
616 );
617 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_finish);
618
619 // shuffle
620 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::shuffle", prof_findDensityTemperatures_shuffle);
621 amrex::ParallelForRNG( n_cells,
622 [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
623 {
624 // The particles from species1 that are in the cell `i_cell` are
625 // given by the `indices_1[cell_start_1:cell_stop_1]`
626 index_type const cell_start_1 = cell_offsets_1[i_cell];
627 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
628
629 // Do not shuffle if there is only one particle in the cell
630 if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
631
632 ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
633 }
634 );
635 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_shuffle);
636 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures);
637
638 // Loop over independent particle pairs
639 // To speed up binary collisions on GPU, we try to expose as much parallelism
640 // as possible (while avoiding race conditions): Instead of looping with one GPU
641 // thread per cell, we loop with one GPU thread per "independent pairs" (i.e. pairs
642 // that do not touch the same macroparticles, so that there is no race condition),
643 // where the number of independent pairs is determined by the lower number of
644 // macroparticles of either species, within each cell.
645 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::LoopOverCollisions", prof_loopOverCollisions);
646 amrex::ParallelForRNG( n_independent_pairs,
647 [=] AMREX_GPU_DEVICE (int i_coll, amrex::RandomEngine const& engine) noexcept
648 {
649 // to avoid type mismatch errors
650 auto ui_coll = (index_type)i_coll;
651
652 // Use a bisection algorithm to find the index of the cell in which this pair is located
653 const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
654
655 // The particles from species1 that are in the cell `i_cell` are
656 // given by the `indices_1[cell_start_1:cell_stop_1]`
657 index_type const cell_start_1 = cell_offsets_1[i_cell];
658 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
659 index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2;
660
661 // collision number of the cell
662 const index_type coll_idx = ui_coll - p_coll_offsets[i_cell];
663
664 // Same but for the pairs
665 index_type const cell_start_pair = have_product_species?
666 p_pair_offsets[i_cell] : 0;
667
668 // Get the local density and temperature for this cell
669 amrex::ParticleReal n1 = 0.0;
670 amrex::ParticleReal T1 = 0.0;
671 if (binary_collision_functor.m_computeSpeciesDensities) {
672 n1 = n1_in_each_cell[i_cell];
673 }
674 if (binary_collision_functor.m_computeSpeciesTemperatures) {
675 T1 = T1_in_each_cell[i_cell];
676 }
677
678 amrex::Real global_lamdb = 0.;
679 if (use_global_debye_length) {
680 global_lamdb = global_debye_length_data[i_cell];
681 }
682
683 // Call the function in order to perform collisions
684 // If there are product species, mask, p_pair_indices_1/2, and
685 // p_pair_reaction_weight and p_product_data are filled here
686 binary_collision_functor(
687 cell_start_1, cell_half_1,
688 cell_half_1, cell_stop_1,
689 indices_1, indices_1,
690 soa_1, soa_1, get_position_1, get_position_1,
691 n1, n1, T1, T1, global_lamdb,
692 q1, q1, m1, m1, dt, dV*volume_factor(i_cell), coll_idx,
693 cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
694 p_pair_reaction_weight, p_product_data, engine);
695 }
696 );
697 ABLASTR_PROFILE_VAR_STOP(prof_loopOverCollisions);
698 // Create the new product particles and define their initial values
699 // num_added: how many particles of each product species have been created
700 //NOLINTNEXTLINE(readability-suspicious-call-argument)
701 const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
702 ptile_1, ptile_1,
703 product_species_vector,
704 tile_products_data,
705 m1, m1,
706 products_mass, p_mask, products_np,
707 copy_species1, copy_species2,
708 p_pair_indices_1, p_pair_indices_2,
709 p_pair_reaction_weight, p_product_data);
710
711 for (int i = 0; i < n_product_species; i++)
712 {
713 setNewParticleIDs(*(tile_products_data[i]), static_cast<int>(products_np[i]), num_added[i]);
714 }
715
716 if (m_correct_energy_momentum) {
717 // Loop over cells and calculate the change in energy and momentum.
718 amrex::ParticleReal const energy_fraction = m_energy_fraction;
719
720 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::correctEnergyMomentum", prof_correctEnergyMomentum);
721 amrex::ParallelFor( np1,
722 [=] AMREX_GPU_DEVICE (int i1) noexcept
723 {
724
725 const int i_cell = bins_1_ptr[i1];
726
727 // Subract the momentum from the initial values.
728 // Whatever is left over, will be distributed among the particles.
729 amrex::ParticleReal const px = w1[i1]*u1x[i1]*m1;
730 amrex::ParticleReal const py = w1[i1]*u1y[i1]*m1;
731 amrex::ParticleReal const pz = w1[i1]*u1z[i1]*m1;
732 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], -px);
733 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], -py);
734 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], -pz);
735 }
736 );
737
738 amrex::ParallelFor( np1,
739 [=] AMREX_GPU_DEVICE (int i1) noexcept
740 {
741
742 const int i_cell = bins_1_ptr[i1];
743
744 amrex::ParticleReal const ww_weighted_sum = ww_weighted_sum_in_each_cell[i_cell];
745 amrex::ParticleReal const w_factor = std::pow(w1[i1], beta_weight_exponent - 1._prt);
746 u1x[i1] += w_factor*px_in_each_cell[i_cell]/ww_weighted_sum;
747 u1y[i1] += w_factor*py_in_each_cell[i_cell]/ww_weighted_sum;
748 u1z[i1] += w_factor*pz_in_each_cell[i_cell]/ww_weighted_sum;
749
750 amrex::ParticleReal const KE_after = Algorithms::KineticEnergy(u1x[i1], u1y[i1], u1z[i1], m1);
751 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], -w1[i1]*KE_after);
752 }
753 );
754
755 amrex::Gpu::Buffer<amrex::Long> failed_corrections({0});
756 amrex::Long* failed_corrections_ptr = failed_corrections.data();
757 const int np_warning_threshold = m_np_warning_threshold;
758
759 // The particles are sorted by weight in decreasing order. This gives more of the correction
760 // to the heavier particles, having a smallar effect on their momenta, and typically
761 // requiring fewer particles to absorb or give off the extra energy.
762 const int energy_correction_sort_by_weight_flag = m_energy_correction_sort_by_weight ? sort : no_sort;
763 auto heapSortDecreasing = ParticleUtils::HeapSortDecreasing();
764
765 amrex::ParallelFor(amrex::TypeList<amrex::CompileTimeOptions<no_sort,sort>>{},
766 {energy_correction_sort_by_weight_flag},
767 n_cells,
768 [=] AMREX_GPU_DEVICE (int i_cell, auto energy_correction_sort_by_weight_control) noexcept
769 {
770
771 index_type const cell_start_1 = cell_offsets_1[i_cell];
772 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
773
774 // No correction needed if there is only one particle in the cell
775 if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
776
777 amrex::ParticleReal deltaEp1 = -KE_in_each_cell[i_cell];
778
779 if (deltaEp1 != 0.) {
780
781 // The ignore_unused is needed to fix the error "An extended __device__
782 // lambda cannot first-capture variable in constexpr-if context"
783 amrex::ignore_unused(heapSortDecreasing, indices_1, w1);
784 if constexpr (energy_correction_sort_by_weight_control == sort) {
785 int const numCell1 = (cell_stop_1 - cell_start_1);
786 heapSortDecreasing(indices_1, w1, cell_start_1, numCell1);
787 }
788
789 // Adjust species 1 particles to absorb deltaEp1.
790 const bool correction_failed =
791 ParticleUtils::ModifyEnergyPairwise(u1x, u1y, u1z, w1, indices_1,
792 cell_start_1, cell_stop_1, m1,
793 energy_fraction, deltaEp1);
794 if (correction_failed) {
795 // If the correction failed, give up on the collision and restore the
796 // momentum to what it was beforehand.
797 for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
798 u1x[ indices_1[i1] ] = u1x_before_ptr[ indices_1[i1] ];
799 u1y[ indices_1[i1] ] = u1y_before_ptr[ indices_1[i1] ];
800 u1z[ indices_1[i1] ] = u1z_before_ptr[ indices_1[i1] ];
801 }
802 int const numCell1 = (cell_stop_1 - cell_start_1);
803 if (numCell1 > np_warning_threshold) {
804 amrex::Gpu::Atomic::Add(failed_corrections_ptr, amrex::Long(1));
805 }
806 }
807 }
808
809 }
810 );
811
812 amrex::Long const num_failed_corrections = *(failed_corrections.copyToHost());
813 if (num_failed_corrections > 0) {
814 ablastr::warn_manager::WMRecordWarning("BinaryCollision::doCollisionsWithinTile",
815 "The energy correction failed for " + std::to_string(num_failed_corrections) + " cells " +
816 "for Coulomb collisions for species " + m_species_names[0] + ". " +
817 "The collisions in those cells was cancelled.");
818 }
819
820 ABLASTR_PROFILE_VAR_STOP(prof_correctEnergyMomentum);
821 }
822
823 }
824 else // species_1 != species_2
825 {
826 // Extract particles in the tile that `mfi` points to
827 ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
828 ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
829
830 // Find the particles that are in each cell of this tile
831 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findParticlesInEachCell", prof_findParticlesInEachCell);
832 ParticleBins bins_1 = findParticlesInEachCell( geom_lev, mfi, ptile_1 );
833 ParticleBins bins_2 = findParticlesInEachCell( geom_lev, mfi, ptile_2 );
834 ABLASTR_PROFILE_VAR_STOP(prof_findParticlesInEachCell);
835
836 // Loop over cells, and collide the particles in each cell
837
838 // Extract low-level data
839 auto const n_cells = static_cast<int>(bins_1.numBins());
840 // - Species 1
841 auto np1 = ptile_1.numParticles();
842 const auto soa_1 = ptile_1.getParticleTileData();
843 index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr();
844 index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr();
845 index_type const* AMREX_RESTRICT bins_1_ptr = bins_1.binsPtr();
846 const amrex::ParticleReal q1 = species_1.getCharge();
847 const amrex::ParticleReal m1 = species_1.getMass();
848 auto get_position_1 = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
849 // - Species 2
850 auto np2 = ptile_2.numParticles();
851 const auto soa_2 = ptile_2.getParticleTileData();
852 index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr();
853 index_type const* AMREX_RESTRICT cell_offsets_2 = bins_2.offsetsPtr();
854 index_type const* AMREX_RESTRICT bins_2_ptr = bins_2.binsPtr();
855 const amrex::ParticleReal q2 = species_2.getCharge();
856 const amrex::ParticleReal m2 = species_2.getMass();
857 auto get_position_2 = GetParticlePosition<PIdx>(ptile_2, getpos_offset);
858
859 /*
860 The following calculations are only required when creating product particles
861 */
862 const int n_cells_products = have_product_species ? n_cells: 0;
863 amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
864 index_type* AMREX_RESTRICT p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
865
866 // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
867 // For different species, the number of pairs in a cell is the number of particles of
868 // the species that has the most particles in that cell
869 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::computeNumberOfPairs", prof_computeNumberOfPairs);
870 amrex::ParallelFor( n_cells_products,
871 [=] AMREX_GPU_DEVICE (int i_cell) noexcept
872 {
873 const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
874 const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
875 // Particular case: no pair if a species has no particle in that cell
876 if (n_part_in_cell_1 == 0 || n_part_in_cell_2 == 0) {
877 p_n_pairs_in_each_cell[i_cell] = 0;
878 } else {
879 p_n_pairs_in_each_cell[i_cell] =
880 amrex::max(n_part_in_cell_1,n_part_in_cell_2);
881 }
882 }
883 );
884
885 // Start indices of the pairs in a cell. Will be used for particle creation
886 amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
887 const index_type n_total_pairs = (n_cells_products == 0) ? 0:
888 amrex::Scan::ExclusiveSum(n_cells_products,
889 p_n_pairs_in_each_cell, pair_offsets.data());
890 index_type* AMREX_RESTRICT p_pair_offsets = pair_offsets.dataPtr();
891
892 amrex::Gpu::DeviceVector<index_type> n_ind_pairs_in_each_cell(n_cells+1);
893 index_type* AMREX_RESTRICT p_n_ind_pairs_in_each_cell = n_ind_pairs_in_each_cell.dataPtr();
894
895 amrex::ParallelFor( n_cells+1,
896 [=] AMREX_GPU_DEVICE (int i_cell) noexcept
897 {
898 if (i_cell < n_cells)
899 {
900 const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
901 const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
902 p_n_ind_pairs_in_each_cell[i_cell] = amrex::min(n_part_in_cell_1, n_part_in_cell_2);
903 }
904 else
905 {
906 p_n_ind_pairs_in_each_cell[i_cell] = 0;
907 }
908 }
909 );
910
911 // start indices of independent collisions.
912 amrex::Gpu::DeviceVector<index_type> coll_offsets(n_cells+1);
913 // number of total independent collision pairs
914 const auto n_independent_pairs = (int) amrex::Scan::ExclusiveSum(n_cells+1,
915 p_n_ind_pairs_in_each_cell, coll_offsets.data(), amrex::Scan::RetSum{true});
916 index_type* AMREX_RESTRICT p_coll_offsets = coll_offsets.dataPtr();
917 ABLASTR_PROFILE_VAR_STOP(prof_computeNumberOfPairs);
918
919 // mask: for DSMC, the index+1 of the selected scattering process (0 = no collision);
920 // for other collision types, 1 if particle creation occurs, 0 otherwise
921 amrex::Gpu::DeviceVector<index_type> mask(n_total_pairs, index_type(0));
922 index_type* AMREX_RESTRICT p_mask = mask.dataPtr();
923 // Will be filled with the index of the first particle of a given pair
924 amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
925 index_type* AMREX_RESTRICT p_pair_indices_1 = pair_indices_1.dataPtr();
926 // Will be filled with the index of the second particle of a given pair
927 amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
928 index_type* AMREX_RESTRICT p_pair_indices_2 = pair_indices_2.dataPtr();
929 // How much weight should be given to the produced particles (and removed from the
930 // reacting particles)
931 amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
932 amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight =
933 pair_reaction_weight.dataPtr();
934 // Extra data needed when products are created
935 int const n_product_data = (binary_collision_functor.m_need_product_data ? n_total_pairs : 0);
936 amrex::Gpu::DeviceVector<amrex::ParticleReal> product_data(n_product_data);
937 amrex::ParticleReal* AMREX_RESTRICT p_product_data = product_data.dataPtr();
938 /*
939 End of calculations only required when creating product particles
940 */
941
942 amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
943 amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
944 amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
945 amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
946
947 amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
948 amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
949 amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
950 amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
951
952 // create vectors to store density and temperature on cell level
953 amrex::Gpu::DeviceVector<amrex::ParticleReal> n1_vec, n2_vec;
954 amrex::Gpu::DeviceVector<amrex::ParticleReal> T1_vec, T2_vec;
955 amrex::Gpu::DeviceVector<amrex::ParticleReal> vx1_vec, vx2_vec;
956 amrex::Gpu::DeviceVector<amrex::ParticleReal> vy1_vec, vy2_vec;
957 amrex::Gpu::DeviceVector<amrex::ParticleReal> vz1_vec, vz2_vec;
958 amrex::Gpu::DeviceVector<amrex::ParticleReal> vs1_vec, vs2_vec;
959
960 if (binary_collision_functor.m_computeSpeciesDensities) {
961 n1_vec.resize(n_cells, 0.0_prt);
962 n2_vec.resize(n_cells, 0.0_prt);
963 }
964 if (binary_collision_functor.m_computeSpeciesTemperatures) {
965 T1_vec.resize(n_cells, 0.0_prt);
966 T2_vec.resize(n_cells, 0.0_prt);
967 vx1_vec.resize(n_cells, 0.0_prt);
968 vx2_vec.resize(n_cells, 0.0_prt);
969 vy1_vec.resize(n_cells, 0.0_prt);
970 vy2_vec.resize(n_cells, 0.0_prt);
971 vz1_vec.resize(n_cells, 0.0_prt);
972 vz2_vec.resize(n_cells, 0.0_prt);
973 vs1_vec.resize(n_cells, 0.0_prt);
974 vs2_vec.resize(n_cells, 0.0_prt);
975 }
976 amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr();
977 amrex::ParticleReal* AMREX_RESTRICT n2_in_each_cell = n2_vec.dataPtr();
978 amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr();
979 amrex::ParticleReal* AMREX_RESTRICT T2_in_each_cell = T2_vec.dataPtr();
980 amrex::ParticleReal* AMREX_RESTRICT vx1_in_each_cell = vx1_vec.dataPtr();
981 amrex::ParticleReal* AMREX_RESTRICT vx2_in_each_cell = vx2_vec.dataPtr();
982 amrex::ParticleReal* AMREX_RESTRICT vy1_in_each_cell = vy1_vec.dataPtr();
983 amrex::ParticleReal* AMREX_RESTRICT vy2_in_each_cell = vy2_vec.dataPtr();
984 amrex::ParticleReal* AMREX_RESTRICT vz1_in_each_cell = vz1_vec.dataPtr();
985 amrex::ParticleReal* AMREX_RESTRICT vz2_in_each_cell = vz2_vec.dataPtr();
986 amrex::ParticleReal* AMREX_RESTRICT vs1_in_each_cell = vs1_vec.dataPtr();
987 amrex::ParticleReal* AMREX_RESTRICT vs2_in_each_cell = vs2_vec.dataPtr();
988
989 // Create vectors to store energy and momentum in each cell.
990 // This is used to correct the energy and momentum in the
991 // cell after the collisions.
992 amrex::Gpu::DeviceVector<amrex::ParticleReal> ww_weighted_sum_vec;
993 amrex::Gpu::DeviceVector<amrex::ParticleReal> KE_vec;
994 amrex::Gpu::DeviceVector<amrex::ParticleReal> px_vec;
995 amrex::Gpu::DeviceVector<amrex::ParticleReal> py_vec;
996 amrex::Gpu::DeviceVector<amrex::ParticleReal> pz_vec;
997 if (m_correct_energy_momentum) {
998 ww_weighted_sum_vec.resize(n_cells, 0.0_prt);
999 KE_vec.resize(n_cells, 0.0_prt);
1000 px_vec.resize(n_cells, 0.0_prt);
1001 py_vec.resize(n_cells, 0.0_prt);
1002 pz_vec.resize(n_cells, 0.0_prt);
1003 }
1004 amrex::ParticleReal* AMREX_RESTRICT ww_weighted_sum_in_each_cell = ww_weighted_sum_vec.dataPtr();
1005 amrex::ParticleReal* AMREX_RESTRICT KE_in_each_cell = KE_vec.dataPtr();
1006 amrex::ParticleReal* AMREX_RESTRICT px_in_each_cell = px_vec.dataPtr();
1007 amrex::ParticleReal* AMREX_RESTRICT py_in_each_cell = py_vec.dataPtr();
1008 amrex::ParticleReal* AMREX_RESTRICT pz_in_each_cell = pz_vec.dataPtr();
1009
1010 // Save the momentum before the collision since sometimes the correction fails
1011 // and the only solution is to restore the pre-collision values.
1012 // The implicit solver already has ux_n etc so use those if available,
1013 // otherwise create temporaries.
1014 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1x_before;
1015 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1y_before;
1016 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1z_before;
1017 amrex::ParticleReal* AMREX_RESTRICT u1x_before_ptr = nullptr;
1018 amrex::ParticleReal* AMREX_RESTRICT u1y_before_ptr = nullptr;
1019 amrex::ParticleReal* AMREX_RESTRICT u1z_before_ptr = nullptr;
1020 if (m_correct_energy_momentum) {
1021 // Check if ux_n was added.
1022 std::vector<std::string> const & real_names1 = species_1.GetRealSoANames();
1023 auto const pos1 = std::find(real_names1.begin(), real_names1.end(), "ux_n");
1024 if (pos1 != real_names1.end()) {
1025 int const n_builtin_real = PIdx::nattribs;
1026 auto const u1x_ni = static_cast<int>(std::distance(real_names1.begin(), pos1));
1027 int const u1x_runtime_ni = u1x_ni - n_builtin_real;
1028 u1x_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni];
1029 u1y_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 1];
1030 u1z_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 2];
1031 } else {
1032 u1x_before.resize(np1);
1033 u1y_before.resize(np1);
1034 u1z_before.resize(np1);
1035 u1x_before_ptr = u1x_before.dataPtr();
1036 u1y_before_ptr = u1y_before.dataPtr();
1037 u1z_before_ptr = u1z_before.dataPtr();
1038 }
1039 }
1040
1041 amrex::Gpu::DeviceVector<amrex::ParticleReal> u2x_before;
1042 amrex::Gpu::DeviceVector<amrex::ParticleReal> u2y_before;
1043 amrex::Gpu::DeviceVector<amrex::ParticleReal> u2z_before;
1044 amrex::ParticleReal* AMREX_RESTRICT u2x_before_ptr = nullptr;
1045 amrex::ParticleReal* AMREX_RESTRICT u2y_before_ptr = nullptr;
1046 amrex::ParticleReal* AMREX_RESTRICT u2z_before_ptr = nullptr;
1047 if (m_correct_energy_momentum) {
1048 // Check if ux_n was added.
1049 std::vector<std::string> const & real_names2 = species_2.GetRealSoANames();
1050 auto const pos2 = std::find(real_names2.begin(), real_names2.end(), "ux_n");
1051 if (pos2 != real_names2.end()) {
1052 int const n_builtin_real = PIdx::nattribs;
1053 auto const u2x_ni = static_cast<int>(std::distance(real_names2.begin(), pos2));
1054 int const u2x_runtime_ni = u2x_ni - n_builtin_real;
1055 u2x_before_ptr = soa_2.m_runtime_rdata[u2x_runtime_ni];
1056 u2y_before_ptr = soa_2.m_runtime_rdata[u2x_runtime_ni + 1];
1057 u2z_before_ptr = soa_2.m_runtime_rdata[u2x_runtime_ni + 2];
1058 } else {
1059 u2x_before.resize(np2);
1060 u2y_before.resize(np2);
1061 u2z_before.resize(np2);
1062 u2x_before_ptr = u2x_before.dataPtr();
1063 u2y_before_ptr = u2y_before.dataPtr();
1064 u2z_before_ptr = u2z_before.dataPtr();
1065 }
1066 }
1067
1068 amrex::ParticleReal const beta_weight_exponent = m_beta_weight_exponent;
1069
1070 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures", prof_findDensityTemperatures);
1071
1072 if (m_correct_energy_momentum ||
1073 binary_collision_functor.m_computeSpeciesDensities ||
1074 binary_collision_functor.m_computeSpeciesTemperatures) {
1075
1076 bool const correct_energy_momentum = m_correct_energy_momentum;
1077
1078 // Loop over particles and compute quantities needed for energy conservation and the collsion parameters
1079 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::atomics", prof_findDensityTemperatures_atomics);
1080 amrex::ParallelFor( std::max(np1, np2),
1081 [=] AMREX_GPU_DEVICE (int ip) noexcept
1082 {
1083 if (correct_energy_momentum) {
1084 if (ip < np1) {
1085 const int i1 = ip;
1086 const int i_cell = bins_1_ptr[i1];
1087
1088 // Save the momenta before the collision
1089 u1x_before_ptr[i1] = u1x[i1];
1090 u1y_before_ptr[i1] = u1y[i1];
1091 u1z_before_ptr[i1] = u1z[i1];
1092
1093 // Compute the local energy and momentum.
1094 amrex::ParticleReal const w1_scaled = std::pow(w1[i1], beta_weight_exponent);
1095 amrex::Gpu::Atomic::AddNoRet(&ww_weighted_sum_in_each_cell[i_cell], w1_scaled*m1);
1096 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], w1[i1]*u1x[i1]*m1);
1097 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], w1[i1]*u1y[i1]*m1);
1098 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], w1[i1]*u1z[i1]*m1);
1099 const amrex::ParticleReal KE1 = Algorithms::KineticEnergy(u1x[i1], u1y[i1], u1z[i1], m1);
1100 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], w1[i1]*KE1);
1101 }
1102
1103 if (ip < np2) {
1104 const int i2 = ip;
1105 const int i_cell = bins_2_ptr[i2];
1106
1107 // Save the momenta before the collision
1108 u2x_before_ptr[i2] = u2x[i2];
1109 u2y_before_ptr[i2] = u2y[i2];
1110 u2z_before_ptr[i2] = u2z[i2];
1111
1112 // Compute the local energy and momentum.
1113 amrex::ParticleReal const w2_scaled = std::pow(w2[i2], beta_weight_exponent);
1114 amrex::Gpu::Atomic::AddNoRet(&ww_weighted_sum_in_each_cell[i_cell], w2_scaled*m2);
1115 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], w2[i2]*u2x[i2]*m2);
1116 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], w2[i2]*u2y[i2]*m2);
1117 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], w2[i2]*u2z[i2]*m2);
1118 const amrex::ParticleReal KE2 = Algorithms::KineticEnergy(u2x[i2], u2y[i2], u2z[i2], m2);
1119 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], w2[i2]*KE2);
1120 }
1121 }
1122
1123 // compute local densities [1/m^3]
1124 if (binary_collision_functor.m_computeSpeciesDensities) {
1125 if (ip < np1) {
1126 amrex::Gpu::Atomic::AddNoRet(&n1_in_each_cell[bins_1_ptr[ip]],
1127 w1[ip]/(dV*volume_factor(bins_1_ptr[ip])));
1128 }
1129 if (ip < np2) {
1130 amrex::Gpu::Atomic::AddNoRet(&n2_in_each_cell[bins_2_ptr[ip]],
1131 w2[ip]/(dV*volume_factor(bins_2_ptr[ip])));
1132 }
1133 }
1134
1135 // compute local temperatures [Joules]
1136 if (binary_collision_functor.m_computeSpeciesTemperatures) {
1137 if (ip < np1) {
1138 const amrex::ParticleReal us = u1x[ip]*u1x[ip] + u1y[ip]*u1y[ip] + u1z[ip]*u1z[ip];
1139 auto constexpr inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
1140 const amrex::ParticleReal gm = std::sqrt( 1.0_prt + us*inv_c2);
1141 amrex::Gpu::Atomic::AddNoRet(&T1_in_each_cell[bins_1_ptr[ip]], w1[ip]);
1142 amrex::Gpu::Atomic::AddNoRet(&vx1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1x[ip]/gm);
1143 amrex::Gpu::Atomic::AddNoRet(&vy1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1y[ip]/gm);
1144 amrex::Gpu::Atomic::AddNoRet(&vz1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1z[ip]/gm);
1145 amrex::Gpu::Atomic::AddNoRet(&vs1_in_each_cell[bins_1_ptr[ip]], w1[ip]*us/gm/gm);
1146 }
1147 if (ip < np2) {
1148 const amrex::ParticleReal us = u2x[ip]*u2x[ip] + u2y[ip]*u2y[ip] + u2z[ip]*u2z[ip];
1149 auto constexpr inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
1150 const amrex::ParticleReal gm = std::sqrt( 1.0_prt + us*inv_c2);
1151 amrex::Gpu::Atomic::AddNoRet(&T2_in_each_cell [bins_2_ptr[ip]], w2[ip]);
1152 amrex::Gpu::Atomic::AddNoRet(&vx2_in_each_cell[bins_2_ptr[ip]], w2[ip]*u2x[ip]/gm);
1153 amrex::Gpu::Atomic::AddNoRet(&vy2_in_each_cell[bins_2_ptr[ip]], w2[ip]*u2y[ip]/gm);
1154 amrex::Gpu::Atomic::AddNoRet(&vz2_in_each_cell[bins_2_ptr[ip]], w2[ip]*u2z[ip]/gm);
1155 amrex::Gpu::Atomic::AddNoRet(&vs2_in_each_cell[bins_2_ptr[ip]], w2[ip]*us/gm/gm);
1156 }
1157 }
1158 }
1159 );
1160 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_atomics);
1161
1162 }
1163
1164 // Finish temperature calculation - we launch 2*n_cells threads to compute both species simultaneously
1165 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::finishTemperature", prof_findDensityTemperatures_finish);
1166 amrex::ParallelFor( 2*n_cells, [=] AMREX_GPU_DEVICE (int i) noexcept
1167 {
1168 const int i_cell = i < n_cells ? i : i - n_cells;
1169
1170 // The particles from species1 that are in the cell `i_cell` are
1171 // given by the `indices_1[cell_start_1:cell_stop_1]`
1172 index_type const cell_start_1 = cell_offsets_1[i_cell];
1173 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1174
1175 // Same for species 2
1176 index_type const cell_start_2 = cell_offsets_2[i_cell];
1177 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1178
1179 // Do not need if one species is missing in the cell
1180 if ( cell_stop_1 - cell_start_1 < 1 ||
1181 cell_stop_2 - cell_start_2 < 1 ) { return; }
1182
1183 // finish temperature calculation if needed
1184 if (binary_collision_functor.m_computeSpeciesTemperatures) {
1185 if (i < n_cells) {
1186 const amrex::ParticleReal invsum1 = 1._prt/T1_in_each_cell[i_cell];
1187 auto vx1 = vx1_in_each_cell[i_cell] * invsum1;
1188 auto vy1 = vy1_in_each_cell[i_cell] * invsum1;
1189 auto vz1 = vz1_in_each_cell[i_cell] * invsum1;
1190 auto vs1 = vs1_in_each_cell[i_cell] * invsum1;
1191
1192 T1_in_each_cell[i_cell] = m1/(3._prt)*(vs1 -(vx1*vx1+vy1*vy1+vz1*vz1));
1193 } else {
1194 const amrex::ParticleReal invsum2 = 1._prt/T2_in_each_cell[i_cell];
1195 auto vx2 = vx2_in_each_cell[i_cell] * invsum2;
1196 auto vy2 = vy2_in_each_cell[i_cell] * invsum2;
1197 auto vz2 = vz2_in_each_cell[i_cell] * invsum2;
1198 auto vs2 = vs2_in_each_cell[i_cell] * invsum2;
1199
1200 T2_in_each_cell[i_cell] = m2/(3._prt)*(vs2 -(vx2*vx2+vy2*vy2+vz2*vz2));
1201 }
1202 }
1203 }
1204 );
1205 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_finish);
1206
1207 // shuffle - we launch 2*n_cells threads to compute both species simultaneously
1208 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::shuffle", prof_findDensityTemperatures_shuffle);
1209 amrex::ParallelForRNG( 2*n_cells,
1210 [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept
1211 {
1212 const int i_cell = i < n_cells ? i : i - n_cells;
1213
1214 // The particles from species1 that are in the cell `i_cell` are
1215 // given by the `indices_1[cell_start_1:cell_stop_1]`
1216 index_type const cell_start_1 = cell_offsets_1[i_cell];
1217 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1218
1219 // Same for species 2
1220 index_type const cell_start_2 = cell_offsets_2[i_cell];
1221 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1222
1223 // Do not collide if one species is missing in the cell
1224 if ( cell_stop_1 - cell_start_1 < 1 ||
1225 cell_stop_2 - cell_start_2 < 1 ) { return; }
1226
1227 if (i < n_cells) {
1228 ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
1229 } else {
1230 ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine);
1231 }
1232 }
1233 );
1234 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_shuffle);
1235
1236 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures);
1237 // Loop over independent particle pairs
1238 // To speed up binary collisions on GPU, we try to expose as much parallelism
1239 // as possible (while avoiding race conditions): Instead of looping with one GPU
1240 // thread per cell, we loop with one GPU thread per "independent pairs" (i.e. pairs
1241 // that do not touch the same macroparticles, so that there is no race condition),
1242 // where the number of independent pairs is determined by the lower number of
1243 // macroparticles of either species, within each cell.
1244 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::LoopOverCollisions", prof_loopOverCollisions);
1245 amrex::ParallelForRNG( n_independent_pairs,
1246 [=] AMREX_GPU_DEVICE (int i_coll, amrex::RandomEngine const& engine) noexcept
1247 {
1248 // to avoid type mismatch errors
1249 auto ui_coll = (index_type)i_coll;
1250
1251 // Use a bisection algorithm to find the index of the cell in which this pair is located
1252 const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
1253
1254 // The particles from species1 that are in the cell `i_cell` are
1255 // given by the `indices_1[cell_start_1:cell_stop_1]`
1256 index_type const cell_start_1 = cell_offsets_1[i_cell];
1257 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1258 // Same for species 2
1259 index_type const cell_start_2 = cell_offsets_2[i_cell];
1260 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1261
1262 // collision number of the cell
1263 const index_type coll_idx = ui_coll - p_coll_offsets[i_cell];
1264
1265 // Same but for the pairs
1266 index_type const cell_start_pair = have_product_species?
1267 p_pair_offsets[i_cell]: 0;
1268
1269 // ux from species1 can be accessed like this:
1270 // ux_1[ indices_1[i] ], where i is between
1271 // cell_start_1 (inclusive) and cell_start_2 (exclusive)
1272
1273 // Get the local densities and temperatures for this cell
1274 amrex::ParticleReal n1 = 0.0, n2 = 0.0;
1275 amrex::ParticleReal T1 = 0.0, T2 = 0.0;
1276 if (binary_collision_functor.m_computeSpeciesDensities) {
1277 n1 = n1_in_each_cell[i_cell];
1278 n2 = n2_in_each_cell[i_cell];
1279 }
1280 if (binary_collision_functor.m_computeSpeciesTemperatures) {
1281 T1 = T1_in_each_cell[i_cell];
1282 T2 = T2_in_each_cell[i_cell];
1283 }
1284
1285 amrex::Real global_lamdb = 0.;
1286 if (use_global_debye_length) {
1287 global_lamdb = global_debye_length_data[i_cell];
1288 }
1289
1290 // Call the function in order to perform collisions
1291 // If there are product species, p_mask, p_pair_indices_1/2, and
1292 // p_pair_reaction_weight and p_product_data are filled here
1293 binary_collision_functor(
1294 cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
1295 indices_1, indices_2,
1296 soa_1, soa_2, get_position_1, get_position_2,
1297 n1, n2, T1, T2, global_lamdb,
1298 q1, q2, m1, m2, dt, dV*volume_factor(i_cell), coll_idx,
1299 cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
1300 p_pair_reaction_weight, p_product_data, engine);
1301 }
1302 );
1303 ABLASTR_PROFILE_VAR_STOP(prof_loopOverCollisions);
1304
1305 // Create the new product particles and define their initial values
1306 // num_added: how many particles of each product species have been created
1307 //NOLINTNEXTLINE(readability-suspicious-call-argument)
1308 const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
1309 ptile_1, ptile_2,
1310 product_species_vector,
1311 tile_products_data,
1312 m1, m2,
1313 products_mass, p_mask, products_np,
1314 copy_species1, copy_species2,
1315 p_pair_indices_1, p_pair_indices_2,
1316 p_pair_reaction_weight, p_product_data);
1317
1318 for (int i = 0; i < n_product_species; i++)
1319 {
1320 setNewParticleIDs(*(tile_products_data[i]), static_cast<int>(products_np[i]), num_added[i]);
1321 }
1322
1323 if (m_correct_energy_momentum) {
1324 // Loop over cells and calculate the change in energy and momentum.
1325 amrex::ParticleReal const energy_fraction = m_energy_fraction;
1326
1327 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::correctEnergyMomentum", prof_correctEnergyMomentum);
1328 amrex::ParallelFor( std::max(np1, np2),
1329 [=] AMREX_GPU_DEVICE (int ip) noexcept
1330 {
1331
1332 if (ip < np1) {
1333 const int i1 = ip;
1334 const int i_cell = bins_1_ptr[i1];
1335
1336 // Subract the momentum from the initial values.
1337 // Whatever is left over, will be distributed among the particles.
1338 amrex::ParticleReal const px = w1[i1]*u1x[i1]*m1;
1339 amrex::ParticleReal const py = w1[i1]*u1y[i1]*m1;
1340 amrex::ParticleReal const pz = w1[i1]*u1z[i1]*m1;
1341 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], -px);
1342 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], -py);
1343 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], -pz);
1344 }
1345
1346 if (ip < np2) {
1347 const int i2 = ip;
1348 const int i_cell = bins_2_ptr[i2];
1349
1350 // Subract the momentum from the initial values.
1351 // Whatever is left over, will be distributed among the particles.
1352 amrex::ParticleReal const px = w2[i2]*u2x[i2]*m2;
1353 amrex::ParticleReal const py = w2[i2]*u2y[i2]*m2;
1354 amrex::ParticleReal const pz = w2[i2]*u2z[i2]*m2;
1355 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], -px);
1356 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], -py);
1357 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], -pz);
1358 }
1359 }
1360 );
1361
1362 amrex::ParallelFor( std::max(np1, np2),
1363 [=] AMREX_GPU_DEVICE (int ip) noexcept
1364 {
1365 if (ip < np1) {
1366 const int i1 = ip;
1367 const int i_cell = bins_1_ptr[i1];
1368
1369 amrex::ParticleReal const ww_weighted_sum = ww_weighted_sum_in_each_cell[i_cell];
1370 amrex::ParticleReal const w_factor = std::pow(w1[i1], beta_weight_exponent - 1._prt);
1371 u1x[i1] += w_factor*px_in_each_cell[i_cell]/ww_weighted_sum;
1372 u1y[i1] += w_factor*py_in_each_cell[i_cell]/ww_weighted_sum;
1373 u1z[i1] += w_factor*pz_in_each_cell[i_cell]/ww_weighted_sum;
1374 }
1375
1376 if (ip < np2) {
1377 const int i2 = ip;
1378 const int i_cell = bins_2_ptr[i2];
1379
1380 amrex::ParticleReal const ww_weighted_sum = ww_weighted_sum_in_each_cell[i_cell];
1381 amrex::ParticleReal const w_factor = std::pow(w2[i2], beta_weight_exponent - 1._prt);
1382 u2x[i2] += w_factor*px_in_each_cell[i_cell]/ww_weighted_sum;
1383 u2y[i2] += w_factor*py_in_each_cell[i_cell]/ww_weighted_sum;
1384 u2z[i2] += w_factor*pz_in_each_cell[i_cell]/ww_weighted_sum;
1385 }
1386 }
1387 );
1388
1389 amrex::Gpu::Buffer<amrex::Long> failed_corrections({0});
1390 amrex::Long* failed_corrections_ptr = failed_corrections.data();
1391 const int np_warning_threshold = m_np_warning_threshold;
1392
1393 // The particles are sorted by weight in decreasing order. This gives more of the correction
1394 // to the heavier particles, having a smallar effect on their momenta, and typically
1395 // requiring fewer particles to absorb or give off the extra energy.
1396 const int energy_correction_sort_by_weight_flag = m_energy_correction_sort_by_weight ? sort : no_sort;
1397 auto heapSortDecreasing = ParticleUtils::HeapSortDecreasing();
1398
1399 amrex::ParallelFor(amrex::TypeList<amrex::CompileTimeOptions<no_sort,sort>>{},
1400 {energy_correction_sort_by_weight_flag},
1401 n_cells,
1402 [=] AMREX_GPU_DEVICE (int i_cell, auto energy_correction_sort_by_weight_control) noexcept
1403 {
1404 // The particles from species1 that are in the cell `i_cell` are
1405 // given by the `indices_1[cell_start_1:cell_stop_1]`.
1406 index_type const cell_start_1 = cell_offsets_1[i_cell];
1407 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1408 // Same for species 2
1409 index_type const cell_start_2 = cell_offsets_2[i_cell];
1410 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1411
1412 // Do not do the rebalance if one species is missing in the cell.
1413 if ( cell_stop_1 - cell_start_1 < 1 ||
1414 cell_stop_2 - cell_start_2 < 1 ) { return; }
1415
1416 amrex::ParticleReal const psq = px_in_each_cell[i_cell]*px_in_each_cell[i_cell] +
1417 py_in_each_cell[i_cell]*py_in_each_cell[i_cell] +
1418 pz_in_each_cell[i_cell]*pz_in_each_cell[i_cell];
1419 if (psq > 0.) {
1420 amrex::ParticleReal w1_sum = 0.0_prt;
1421 amrex::ParticleReal w2_sum = 0.0_prt;
1422 amrex::ParticleReal KE1_after = 0._prt;
1423 amrex::ParticleReal KE2_after = 0._prt;
1424 for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
1425 w1_sum += w1[ indices_1[i1] ];
1426 KE1_after += w1[ indices_1[i1] ]*Algorithms::KineticEnergy(u1x[ indices_1[i1] ],
1427 u1y[ indices_1[i1] ],
1428 u1z[ indices_1[i1] ], m1);
1429 }
1430 for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) {
1431 w2_sum += w2[ indices_2[i2] ];
1432 KE2_after += w2[ indices_2[i2] ]*Algorithms::KineticEnergy(u2x[ indices_2[i2] ],
1433 u2y[ indices_2[i2] ],
1434 u2z[ indices_2[i2] ], m2);
1435 }
1436
1437 const amrex::ParticleReal KE_after = KE1_after + KE2_after;
1438 const amrex::ParticleReal deltaE = KE_after - KE_in_each_cell[i_cell];
1439 amrex::ParticleReal deltaEp1, deltaEp2;
1440 int const numCell1 = (cell_stop_1 - cell_start_1);
1441 int const numCell2 = (cell_stop_2 - cell_start_2);
1442 if (numCell1 == 1) {
1443 deltaEp1 = 0.0;
1444 deltaEp2 = deltaE;
1445 } else if (numCell2 == 1) {
1446 deltaEp1 = deltaE;
1447 deltaEp2 = 0.0;
1448 } else {
1449 const amrex::ParticleReal Etotdenom = w1_sum*KE1_after/numCell1 + w2_sum*KE2_after/numCell2;
1450 deltaEp1 = w1_sum*KE1_after/Etotdenom*deltaE/numCell1;
1451 deltaEp2 = w2_sum*KE2_after/Etotdenom*deltaE/numCell2;
1452 }
1453
1454 // The ignore_unused is needed to fix the error "An extended __device__
1455 // lambda cannot first-capture variable in constexpr-if context"
1456 amrex::ignore_unused(heapSortDecreasing);
1457
1458 bool correction1_failed = false;
1459 bool correction2_failed = false;
1460 if (deltaEp1 != 0. && numCell1 > 1) {
1461 if constexpr (energy_correction_sort_by_weight_control == sort) {
1462 heapSortDecreasing(indices_1, w1, cell_start_1, numCell1);
1463 }
1464 // Adjust species 1 particles to absorb deltaEp1.
1465 correction1_failed =
1466 ParticleUtils::ModifyEnergyPairwise(u1x, u1y, u1z, w1, indices_1,
1467 cell_start_1, cell_stop_1, m1,
1468 energy_fraction, deltaEp1);
1469 }
1470
1471 if (deltaEp2 != 0. && numCell2 > 1) {
1472 if constexpr (energy_correction_sort_by_weight_control == sort) {
1473 heapSortDecreasing(indices_2, w2, cell_start_2, numCell2);
1474 }
1475 // Adjust species 2 particles to absorb deltaEp2.
1476 correction2_failed =
1477 ParticleUtils::ModifyEnergyPairwise(u2x, u2y, u2z, w2, indices_2,
1478 cell_start_2, cell_stop_2, m2,
1479 energy_fraction, deltaEp2);
1480 }
1481
1482 if (correction1_failed || correction2_failed ||
1483 (numCell1 == 1 && numCell2 == 1)) {
1484 // If the correction failed, or if each species has only one particle
1485 // and the correction can not be applied, give up on the collision and
1486 // restore the momentum to what it was beforehand.
1487 for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
1488 u1x[ indices_1[i1] ] = u1x_before_ptr[ indices_1[i1] ];
1489 u1y[ indices_1[i1] ] = u1y_before_ptr[ indices_1[i1] ];
1490 u1z[ indices_1[i1] ] = u1z_before_ptr[ indices_1[i1] ];
1491 }
1492 for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) {
1493 u2x[ indices_2[i2] ] = u2x_before_ptr[ indices_2[i2] ];
1494 u2y[ indices_2[i2] ] = u2y_before_ptr[ indices_2[i2] ];
1495 u2z[ indices_2[i2] ] = u2z_before_ptr[ indices_2[i2] ];
1496 }
1497 if ((correction1_failed && numCell1 > np_warning_threshold) ||
1498 (correction2_failed && numCell2 > np_warning_threshold)) {
1499 amrex::Gpu::Atomic::Add(failed_corrections_ptr, amrex::Long(1));
1500 }
1501 }
1502
1503 }
1504 }
1505 );
1506
1507 amrex::Long const num_failed_corrections = *(failed_corrections.copyToHost());
1508 if (num_failed_corrections > 0) {
1509 ablastr::warn_manager::WMRecordWarning("BinaryCollision::doCollisionsWithinTile",
1510 "The energy correction failed for " + std::to_string(num_failed_corrections) + " cells " +
1511 "for Coulomb collisions between species " + m_species_names[0] + " and " + m_species_names[1] + ". " +
1512 "The collisions in those cells was cancelled.");
1513 }
1514
1515 ABLASTR_PROFILE_VAR_STOP(prof_correctEnergyMomentum);
1516 }
1517
1518 } // end if ( m_isSameSpecies)
1519
1520 }
1521
1522private:
1523
1524 bool m_isSameSpecies;
1525 bool m_have_product_species;
1526
1527 bool m_correct_energy_momentum = false;
1528 bool m_energy_correction_sort_by_weight = false;
1529 int m_np_warning_threshold = 20;
1530
1531 amrex::ParticleReal m_beta_weight_exponent = 1.;
1532 amrex::ParticleReal m_energy_fraction = 0.05;
1533
1534 enum energy_correction_sort_by_weight_flags : int { no_sort, sort };
1535
1536 amrex::Vector<std::string> m_product_species;
1537 // functor that performs collisions within a cell
1538 CollisionFunctor m_binary_collision_functor;
1539 // functor that creates new particles and initializes their parameters
1540 CopyTransformFunctor m_copy_transform_functor;
1541
1542};
1543
1544#endif // WARPX_PARTICLES_COLLISION_BINARYCOLLISION_H_
AMReX_Algorithm.H
AMReX_BLassert.H
AMReX_BaseFwd.H
AMReX_DenseBins.H
AMReX_Extension.H
AMREX_RESTRICT
#define AMREX_RESTRICT
AMReX_Geometry.H
AMReX_GpuAtomic.H
AMReX_GpuContainers.H
AMReX_GpuControl.H
AMReX_GpuDevice.H
AMReX_GpuLaunch.H
AMReX_GpuQualifiers.H
AMREX_GPU_DEVICE
#define AMREX_GPU_DEVICE
cbx
Box cbx
mask
Array4< int const > mask
AMReX_LayoutData.H
AMReX_MFIter.H
AMReX_PODVector.H
AMReX_ParmParse.H
AMReX_ParticleTile.H
AMReX_Particles.H
AMReX_REAL.H
AMReX_Random.H
AMREX_D_TERM
#define AMREX_D_TERM(a, b, c)
AMReX_Scan.H
AMReX_Utility.H
AMReX_Vector.H
CollisionType
CollisionType
Definition BinaryCollisionUtils.H:19
CollisionType::PairwiseCoulomb
@ PairwiseCoulomb
Definition BinaryCollisionUtils.H:26
CollisionType::DSMC
@ DSMC
Definition BinaryCollisionUtils.H:25
MultiParticleContainer_fwd.H
ABLASTR_PROFILE_VAR
#define ABLASTR_PROFILE_VAR(fname, vname)
Definition ProfilerWrapper.H:14
ABLASTR_PROFILE
#define ABLASTR_PROFILE(fname)
Definition ProfilerWrapper.H:13
ABLASTR_PROFILE_VAR_STOP
#define ABLASTR_PROFILE_VAR_STOP(vname)
Definition ProfilerWrapper.H:17
ShuffleFisherYates
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
setNewParticleIDs
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
TextMsg.H
WARPX_ABORT_WITH_MESSAGE
#define WARPX_ABORT_WITH_MESSAGE(MSG)
Definition TextMsg.H:15
LoadBalanceCostsUpdateAlgo::Timers
@ Timers
load balance according to in-code timer-based weights (i.e., with costs)
Definition WarpXAlgorithmSelection.H:118
WarpXParticleContainer_fwd.H
BinaryCollision::m_have_product_species
bool m_have_product_species
Definition BinaryCollision.H:1525
BinaryCollision::index_type
ParticleBins::index_type index_type
Definition BinaryCollision.H:82
BinaryCollision::ParticleBins
amrex::DenseBins< ParticleTileDataType > ParticleBins
Definition BinaryCollision.H:81
BinaryCollision::m_copy_transform_functor
CopyTransformFunctor m_copy_transform_functor
Definition BinaryCollision.H:1540
BinaryCollision::m_energy_correction_sort_by_weight
bool m_energy_correction_sort_by_weight
Definition BinaryCollision.H:1528
BinaryCollision::m_isSameSpecies
bool m_isSameSpecies
Definition BinaryCollision.H:1524
BinaryCollision::ParticleTileType
WarpXParticleContainer::ParticleTileType ParticleTileType
Definition BinaryCollision.H:79
BinaryCollision::ParticleType
WarpXParticleContainer::ParticleType ParticleType
Definition BinaryCollision.H:78
BinaryCollision::~BinaryCollision
~BinaryCollision() override=default
BinaryCollision::doCollisions
void doCollisions(amrex::Real cur_time, amrex::Real dt, MultiParticleContainer *mypc) override
Definition BinaryCollision.H:178
BinaryCollision::ParticleTileDataType
ParticleTileType::ParticleTileDataType ParticleTileDataType
Definition BinaryCollision.H:80
BinaryCollision::doCollisionsWithinTile
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:278
BinaryCollision::m_np_warning_threshold
int m_np_warning_threshold
Definition BinaryCollision.H:1529
BinaryCollision::m_beta_weight_exponent
amrex::ParticleReal m_beta_weight_exponent
Definition BinaryCollision.H:1531
BinaryCollision::BinaryCollision
BinaryCollision(const std::string &collision_name, MultiParticleContainer const *const mypc)
Constructor of the BinaryCollision class.
Definition BinaryCollision.H:92
BinaryCollision::m_energy_fraction
amrex::ParticleReal m_energy_fraction
Definition BinaryCollision.H:1532
BinaryCollision::m_correct_energy_momentum
bool m_correct_energy_momentum
Definition BinaryCollision.H:1527
BinaryCollision::m_product_species
amrex::Vector< std::string > m_product_species
Definition BinaryCollision.H:1536
BinaryCollision::energy_correction_sort_by_weight_flags
energy_correction_sort_by_weight_flags
Definition BinaryCollision.H:1534
BinaryCollision::no_sort
@ no_sort
Definition BinaryCollision.H:1534
BinaryCollision::sort
@ sort
Definition BinaryCollision.H:1534
BinaryCollision::m_binary_collision_functor
CollisionFunctor m_binary_collision_functor
Definition BinaryCollision.H:1538
BinaryCollision::BinaryCollision
BinaryCollision(BinaryCollision const &)=default
BinaryCollision::operator=
BinaryCollision & operator=(BinaryCollision const &)=default
BinaryCollision::BinaryCollision
BinaryCollision(BinaryCollision &&)=delete
CollisionBase::m_species_names
amrex::Vector< std::string > m_species_names
Definition CollisionBase.H:46
CollisionBase::use_global_debye_length
bool use_global_debye_length() const
Definition CollisionBase.H:41
CollisionBase::CollisionBase
CollisionBase(const std::string &collision_name)
Definition CollisionBase.cpp:14
CollisionBase::m_use_global_debye_length
bool m_use_global_debye_length
Definition CollisionBase.H:50
MultiParticleContainer
Definition MultiParticleContainer.H:69
MultiParticleContainer::GetParticleContainerFromName
WarpXParticleContainer & GetParticleContainerFromName(const std::string &name) const
Definition MultiParticleContainer.cpp:406
NoParticleCreationFunc
This class does nothing and is used as second template parameter for binary collisions that do not cr...
Definition ParticleCreationFunc.H:303
SmartCopyFactory
A factory for creating SmartCopy functors.
Definition SmartCopy.H:133
WarpX
Definition WarpX.H:87
WarpX::load_balance_costs_update_algo
static auto load_balance_costs_update_algo
Definition WarpX.H:198
WarpX::GetInstance
static WarpX & GetInstance()
Definition WarpX.cpp:299
WarpX::getCosts
static amrex::LayoutData< amrex::Real > * getCosts(int lev)
Definition WarpX.cpp:3390
WarpX::LowerCorner
static amrex::XDim3 LowerCorner(const amrex::Box &bx, int lev, amrex::Real time_shift_delta)
Return the lower corner of the box in real units.
Definition WarpX.cpp:3223
WarpXParticleContainer
Definition WarpXParticleContainer.H:195
WarpXParticleContainer::defineAllParticleTiles
void defineAllParticleTiles() noexcept
Definition WarpXParticleContainer.cpp:2802
WarpXParticleContainer::getCharge
amrex::ParticleReal getCharge() const
Definition WarpXParticleContainer.H:581
WarpXParticleContainer::getMass
amrex::ParticleReal getMass() const
Definition WarpXParticleContainer.H:583
amrex::AmrMesh::Geom
const Vector< Geometry > & Geom() const noexcept
amrex::BaseFab::dataPtr
T * dataPtr(int n=0) noexcept
amrex::CoordSys::CellSize
const Real * CellSize() const noexcept
amrex::DenseBins::offsetsPtr
index_type * offsetsPtr() noexcept
amrex::DenseBins::permutationPtr
index_type * permutationPtr() noexcept
amrex::DenseBins::numBins
Long numBins() const noexcept
amrex::DenseBins::binsPtr
index_type * binsPtr() noexcept
amrex::IntVectND< 3 >::TheZeroVector
__host__ static __device__ constexpr IntVectND< dim > TheZeroVector() noexcept
amrex::MFIter::tilebox
Box tilebox() const noexcept
amrex::MFIter::isValid
bool isValid() const noexcept
amrex::PODVector::begin
iterator begin() noexcept
amrex::PODVector::resize
void resize(size_type a_new_size, GrowthStrategy strategy=GrowthStrategy::Poisson)
amrex::PODVector::dataPtr
T * dataPtr() noexcept
amrex::PODVector::data
T * data() noexcept
amrex::ParmParse::queryarr
int queryarr(std::string_view name, std::vector< int > &ref, int start_ix=FIRST, int num_val=ALL) const
amrex::ParmParse::query
int query(std::string_view name, bool &ref, int ival=FIRST) const
amrex::ParticleContainer_impl::ParticlesAt
const ParticleTileType & ParticlesAt(int lev, int grid, int tile) const
amrex::ParticleContainer_impl::GetRealSoANames
std::vector< std::string > GetRealSoANames() const
amrex::ParticleContainer_impl< SoAParticle< T_NArrayReal, T_NArrayInt >, T_NArrayReal, T_NArrayInt, Allocator, CellAssignor >::ParticleTileType
std::conditional_t< is_rtsoa_pc, ParticleTileRT< typename ParticleType::RealType, typename ParticleType::IntType >, ParticleTile< ParticleType, NArrayReal, NArrayInt, Allocator > > ParticleTileType
amrex::ParticleContainer_impl< SoAParticle< T_NArrayReal, T_NArrayInt >, T_NArrayReal, T_NArrayInt, Allocator, CellAssignor >::ParticleType
T_ParticleType ParticleType
amrex::Real
amrex_real Real
amrex::ParticleReal
amrex_particle_real ParticleReal
amrex::Long
amrex_long Long
amrex::Scan::ExclusiveSum
T ExclusiveSum(N n, T const *in, T *out, RetSum a_ret_sum=retSum)
amrex::ubound
__host__ __device__ Dim3 ubound(Array4< T > const &a) noexcept
amrex::lbound
__host__ __device__ Dim3 lbound(Array4< T > const &a) noexcept
amrex::Gpu::DeviceVector
PODVector< T, ArenaAllocator< T > > DeviceVector
amrex::bisect
__host__ __device__ T bisect(T lo, T hi, F f, T tol=1e-12, int max_iter=100)
amrex::min
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
amrex::max
__host__ __device__ constexpr const T & max(const T &a, const T &b) noexcept
Algorithms::KineticEnergy
AMREX_GPU_HOST_DEVICE AMREX_INLINE T KineticEnergy(const amrex::ParticleReal ux, const amrex::ParticleReal uy, const amrex::ParticleReal uz, const amrex::ParticleReal mass)
Computes the kinetic energy of a particle. This method should not be used with photons.
Definition KineticEnergy.H:36
BinaryCollisionUtils::get_collision_type
CollisionType get_collision_type(const std::string &collision_name, MultiParticleContainer const *const mypc)
Definition BinaryCollisionUtils.cpp:27
ParticleUtils
Definition ParticleUtils.cpp:24
ParticleUtils::ModifyEnergyPairwise
AMREX_GPU_HOST_DEVICE AMREX_INLINE bool ModifyEnergyPairwise(amrex::ParticleReal *const AMREX_RESTRICT uxp, amrex::ParticleReal *const AMREX_RESTRICT uyp, amrex::ParticleReal *const AMREX_RESTRICT uzp, const amrex::ParticleReal *const AMREX_RESTRICT w, T_index const *const AMREX_RESTRICT indices, T_index const cell_start, T_index const cell_stop, amrex::ParticleReal const mass, amrex::ParticleReal const energy_fraction, amrex::ParticleReal &deltaE)
Definition ParticleUtils.H:385
ParticleUtils::findParticlesInEachCell
amrex::DenseBins< ParticleTileDataType > findParticlesInEachCell(amrex::Geometry const &geom_lev, amrex::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:35
ablastr::constant::SI::inv_c2_v
constexpr auto inv_c2_v
inverse of the square of the vacuum speed of light [s^2/m^2] (variable template)
Definition constant.H:153
ablastr::constant::math::pi
constexpr auto pi
ratio of a circle's circumference to its diameter
Definition constant.H:29
ablastr::warn_manager::WMRecordWarning
void WMRecordWarning(const std::string &topic, const std::string &text, const WarnPriority &priority=WarnPriority::medium)
Helper function to abbreviate the call to WarnManager::GetInstance().RecordWarning (recording a warni...
Definition WarnManager.cpp:320
amrex::Gpu::Atomic::AddNoRet
__host__ __device__ AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
amrex::Gpu::Atomic::Add
__host__ __device__ AMREX_FORCE_INLINE T Add(T *sum, T value) noexcept
amrex::Gpu::synchronize
void synchronize() noexcept
amrex::Gpu::copyAsync
void copyAsync(HostToDevice, InIter begin, InIter end, OutIter result) noexcept
amrex::Gpu::hostToDevice
static constexpr HostToDevice hostToDevice
amrex::Gpu::streamSynchronize
void streamSynchronize() noexcept
amrex::Gpu::notInLaunchRegion
bool notInLaunchRegion() noexcept
amrex::HostDevice::Atomic::Add
__host__ __device__ AMREX_FORCE_INLINE void Add(T *const sum, T const value) noexcept
amrex::ignore_unused
__host__ __device__ void ignore_unused(const Ts &...)
amrex::ParallelFor
void ParallelFor(TypeList< CTOs... > ctos, std::array< int, sizeof...(CTOs)> const &runtime_options, T N, F &&f)
amrex::Box
BoxND< 3 > Box
amrex::second
double second() noexcept
amrex::ParallelForRNG
AMREX_ATTRIBUTE_FLATTEN_FOR void ParallelForRNG(T n, L const &f) noexcept
amrex::int
const int[]
warpx::fields::FieldType::global_debye_length
@ global_debye_length
Definition Fields.H:97
warpx
Definition FieldBoundaries.cpp:18
GetParticlePosition
Functor that can be used to extract the positions of the macroparticles inside a ParallelFor kernel.
Definition GetAndSetPosition.H:75
PIdx::nattribs
@ nattribs
Definition WarpXParticleContainer.H:70
PIdx::uz
@ uz
Definition WarpXParticleContainer.H:70
PIdx::w
@ w
Definition WarpXParticleContainer.H:70
PIdx::uy
@ uy
Definition WarpXParticleContainer.H:70
PIdx::ux
@ ux
Definition WarpXParticleContainer.H:70
ParticleUtils::HeapSortDecreasing
A GPU compartible sort of a index vector vector based on another GPU vector's values,...
Definition ParticleUtils.H:542
SmartCopy
This is a functor for performing a "smart copy" that works in both host and device code.
Definition SmartCopy.H:34
amrex::MFItInfo::SetDynamic
MFItInfo & SetDynamic(bool f) noexcept
amrex::MFItInfo::EnableTiling
MFItInfo & EnableTiling(const IntVect &ts=FabArrayBase::mfiter_tile_size) noexcept
Generated by doxygen 1.17.0