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: equal to 1 if particle creation occurs for a given pair, 0 otherwise
433 amrex::Gpu::DeviceVector<index_type> mask(n_total_pairs);
434 index_type* AMREX_RESTRICT p_mask = mask.dataPtr();
435 // Will be filled with the index of the first particle of a given pair
436 amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
437 index_type* AMREX_RESTRICT p_pair_indices_1 = pair_indices_1.dataPtr();
438 // Will be filled with the index of the second particle of a given pair
439 amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
440 index_type* AMREX_RESTRICT p_pair_indices_2 = pair_indices_2.dataPtr();
441 // How much weight should be given to the produced particles (and removed from the
442 // reacting particles)
443 amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
444 amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight =
445 pair_reaction_weight.dataPtr();
446 // Extra data needed when products are created
447 int const n_product_data = (binary_collision_functor.m_need_product_data ? n_total_pairs : 0);
448 amrex::Gpu::DeviceVector<amrex::ParticleReal> product_data(n_product_data);
449 amrex::ParticleReal* AMREX_RESTRICT p_product_data = product_data.dataPtr();
450 /*
451 End of calculations only required when creating product particles
452 */
453
454 amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
455 amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
456 amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
457 amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
458
459 // create vectors to store density and temperature on cell level
460 amrex::Gpu::DeviceVector<amrex::ParticleReal> n1_vec;
461 amrex::Gpu::DeviceVector<amrex::ParticleReal> T1_vec;
462 amrex::Gpu::DeviceVector<amrex::ParticleReal> vx1_vec;
463 amrex::Gpu::DeviceVector<amrex::ParticleReal> vy1_vec;
464 amrex::Gpu::DeviceVector<amrex::ParticleReal> vz1_vec;
465 amrex::Gpu::DeviceVector<amrex::ParticleReal> vs1_vec;
466
467 if (binary_collision_functor.m_computeSpeciesDensities) {
468 n1_vec.resize(n_cells, 0.0_prt);
469 }
470 if (binary_collision_functor.m_computeSpeciesTemperatures) {
471 T1_vec.resize(n_cells, 0.0_prt);
472 vx1_vec.resize(n_cells, 0.0_prt);
473 vy1_vec.resize(n_cells, 0.0_prt);
474 vz1_vec.resize(n_cells, 0.0_prt);
475 vs1_vec.resize(n_cells, 0.0_prt);
476 }
477 amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr();
478 amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr();
479 amrex::ParticleReal* AMREX_RESTRICT vx1_in_each_cell = vx1_vec.dataPtr();
480 amrex::ParticleReal* AMREX_RESTRICT vy1_in_each_cell = vy1_vec.dataPtr();
481 amrex::ParticleReal* AMREX_RESTRICT vz1_in_each_cell = vz1_vec.dataPtr();
482 amrex::ParticleReal* AMREX_RESTRICT vs1_in_each_cell = vs1_vec.dataPtr();
483
484 // Create vectors to store energy and momentum in each cell.
485 // This is used to correct the energy and momentum in the
486 // cell after the collisions.
487 amrex::Gpu::DeviceVector<amrex::ParticleReal> ww_weighted_sum_vec;
488 amrex::Gpu::DeviceVector<amrex::ParticleReal> KE_vec;
489 amrex::Gpu::DeviceVector<amrex::ParticleReal> px_vec;
490 amrex::Gpu::DeviceVector<amrex::ParticleReal> py_vec;
491 amrex::Gpu::DeviceVector<amrex::ParticleReal> pz_vec;
492 if (m_correct_energy_momentum) {
493 ww_weighted_sum_vec.resize(n_cells, 0.0_prt);
494 KE_vec.resize(n_cells, 0.0_prt);
495 px_vec.resize(n_cells, 0.0_prt);
496 py_vec.resize(n_cells, 0.0_prt);
497 pz_vec.resize(n_cells, 0.0_prt);
498 }
499 amrex::ParticleReal* AMREX_RESTRICT ww_weighted_sum_in_each_cell = ww_weighted_sum_vec.dataPtr();
500 amrex::ParticleReal* AMREX_RESTRICT KE_in_each_cell = KE_vec.dataPtr();
501 amrex::ParticleReal* AMREX_RESTRICT px_in_each_cell = px_vec.dataPtr();
502 amrex::ParticleReal* AMREX_RESTRICT py_in_each_cell = py_vec.dataPtr();
503 amrex::ParticleReal* AMREX_RESTRICT pz_in_each_cell = pz_vec.dataPtr();
504
505 // Save the momentum before the collision since sometimes the correction fails
506 // and the only solution is to restore the pre-collision values.
507 // The implicit solver already has ux_n etc so use those if available,
508 // otherwise create temporaries.
509 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1x_before;
510 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1y_before;
511 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1z_before;
512 amrex::ParticleReal* AMREX_RESTRICT u1x_before_ptr = nullptr;
513 amrex::ParticleReal* AMREX_RESTRICT u1y_before_ptr = nullptr;
514 amrex::ParticleReal* AMREX_RESTRICT u1z_before_ptr = nullptr;
515 if (m_correct_energy_momentum) {
516 // Check if ux_n was added.
517 std::vector<std::string> const & real_names1 = species_1.GetRealSoANames();
518 auto const pos1 = std::find(real_names1.begin(), real_names1.end(), "ux_n");
519 if (pos1 != real_names1.end()) {
520 int const n_builtin_real = PIdx::nattribs;
521 auto const u1x_ni = static_cast<int>(std::distance(real_names1.begin(), pos1));
522 int const u1x_runtime_ni = u1x_ni - n_builtin_real;
523 u1x_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni];
524 u1y_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 1];
525 u1z_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 2];
526 } else {
527 u1x_before.resize(np1);
528 u1y_before.resize(np1);
529 u1z_before.resize(np1);
530 u1x_before_ptr = u1x_before.dataPtr();
531 u1y_before_ptr = u1y_before.dataPtr();
532 u1z_before_ptr = u1z_before.dataPtr();
533 }
534
535 }
536
537 amrex::ParticleReal const beta_weight_exponent = m_beta_weight_exponent;
538
539 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures", prof_findDensityTemperatures);
540
541 if (m_correct_energy_momentum ||
542 binary_collision_functor.m_computeSpeciesDensities ||
543 binary_collision_functor.m_computeSpeciesTemperatures) {
544
545 bool const correct_energy_momentum = m_correct_energy_momentum;
546
547 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::atomics", prof_findDensityTemperatures_atomics);
548 // Loop over particles and compute quantities needed for energy conservation and the collsion parameters
549 amrex::ParallelFor( np1,
550 [=] AMREX_GPU_DEVICE (int ip) noexcept
551 {
552 if (correct_energy_momentum) {
553 u1x_before_ptr[ip] = u1x[ip];
554 u1y_before_ptr[ip] = u1y[ip];
555 u1z_before_ptr[ip] = u1z[ip];
556
557 const int i_cell = bins_1_ptr[ip];
558
559 // Compute the local energy and momentum.
560 amrex::ParticleReal const w1_scaled = std::pow(w1[ip], beta_weight_exponent);
561 amrex::Gpu::Atomic::AddNoRet(&ww_weighted_sum_in_each_cell[i_cell], w1_scaled*m1);
562 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], w1[ip]*u1x[ip]*m1);
563 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], w1[ip]*u1y[ip]*m1);
564 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], w1[ip]*u1z[ip]*m1);
565 const amrex::ParticleReal KE = Algorithms::KineticEnergy(u1x[ip], u1y[ip], u1z[ip], m1);
566 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], w1[ip]*KE);
567 }
568
569 // compute local density [1/m^3]
570 if (binary_collision_functor.m_computeSpeciesDensities) {
571 amrex::Gpu::Atomic::AddNoRet(&n1_in_each_cell[bins_1_ptr[ip]],
572 w1[ip]/(dV*volume_factor(bins_1_ptr[ip])));
573 }
574
575 // compute local temperature [Joules]
576 if (binary_collision_functor.m_computeSpeciesTemperatures) {
577 const amrex::ParticleReal us = u1x[ip]*u1x[ip] + u1y[ip]*u1y[ip] + u1z[ip]*u1z[ip];
578 auto constexpr inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
579 const amrex::ParticleReal gm = std::sqrt( 1.0_prt + us*inv_c2);
580 amrex::Gpu::Atomic::AddNoRet(&T1_in_each_cell[bins_1_ptr[ip]], w1[ip]);
581 amrex::Gpu::Atomic::AddNoRet(&vx1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1x[ip]/gm);
582 amrex::Gpu::Atomic::AddNoRet(&vy1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1y[ip]/gm);
583 amrex::Gpu::Atomic::AddNoRet(&vz1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1z[ip]/gm);
584 amrex::Gpu::Atomic::AddNoRet(&vs1_in_each_cell[bins_1_ptr[ip]], w1[ip]*us/gm/gm);
585 }
586 }
587 );
588 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_atomics);
589
590 }
591
592 // Finish temperature calculation
593 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::finishTemperature", prof_findDensityTemperatures_finish);
594 amrex::ParallelFor( n_cells, [=] AMREX_GPU_DEVICE (int i_cell) noexcept
595 {
596 // The particles from species1 that are in the cell `i_cell` are
597 // given by the `indices_1[cell_start_1:cell_stop_1]`
598 index_type const cell_start_1 = cell_offsets_1[i_cell];
599 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
600
601 // Do not need if there is only one particle in the cell
602 if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
603
604 // finish temperature calculation if needed
605 if (binary_collision_functor.m_computeSpeciesTemperatures) {
606 const amrex::ParticleReal invsum = 1._prt/T1_in_each_cell[i_cell];
607 auto vx1 = vx1_in_each_cell[i_cell] * invsum;
608 auto vy1 = vy1_in_each_cell[i_cell] * invsum;
609 auto vz1 = vz1_in_each_cell[i_cell] * invsum;
610 auto vs1 = vs1_in_each_cell[i_cell] * invsum;
611
612 T1_in_each_cell[i_cell] = m1/(3._prt)*(vs1 -(vx1*vx1+vy1*vy1+vz1*vz1));
613 }
614 }
615 );
616 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_finish);
617
618 // shuffle
619 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::shuffle", prof_findDensityTemperatures_shuffle);
620 amrex::ParallelForRNG( n_cells,
621 [=] AMREX_GPU_DEVICE (int i_cell, amrex::RandomEngine const& engine) noexcept
622 {
623 // The particles from species1 that are in the cell `i_cell` are
624 // given by the `indices_1[cell_start_1:cell_stop_1]`
625 index_type const cell_start_1 = cell_offsets_1[i_cell];
626 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
627
628 // Do not shuffle if there is only one particle in the cell
629 if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
630
631 ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
632 }
633 );
634 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_shuffle);
635 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures);
636
637 // Loop over independent particle pairs
638 // To speed up binary collisions on GPU, we try to expose as much parallelism
639 // as possible (while avoiding race conditions): Instead of looping with one GPU
640 // thread per cell, we loop with one GPU thread per "independent pairs" (i.e. pairs
641 // that do not touch the same macroparticles, so that there is no race condition),
642 // where the number of independent pairs is determined by the lower number of
643 // macroparticles of either species, within each cell.
644 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::LoopOverCollisions", prof_loopOverCollisions);
645 amrex::ParallelForRNG( n_independent_pairs,
646 [=] AMREX_GPU_DEVICE (int i_coll, amrex::RandomEngine const& engine) noexcept
647 {
648 // to avoid type mismatch errors
649 auto ui_coll = (index_type)i_coll;
650
651 // Use a bisection algorithm to find the index of the cell in which this pair is located
652 const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
653
654 // The particles from species1 that are in the cell `i_cell` are
655 // given by the `indices_1[cell_start_1:cell_stop_1]`
656 index_type const cell_start_1 = cell_offsets_1[i_cell];
657 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
658 index_type const cell_half_1 = (cell_start_1+cell_stop_1)/2;
659
660 // collision number of the cell
661 const index_type coll_idx = ui_coll - p_coll_offsets[i_cell];
662
663 // Same but for the pairs
664 index_type const cell_start_pair = have_product_species?
665 p_pair_offsets[i_cell] : 0;
666
667 // Get the local density and temperature for this cell
668 amrex::ParticleReal n1 = 0.0;
669 amrex::ParticleReal T1 = 0.0;
670 if (binary_collision_functor.m_computeSpeciesDensities) {
671 n1 = n1_in_each_cell[i_cell];
672 }
673 if (binary_collision_functor.m_computeSpeciesTemperatures) {
674 T1 = T1_in_each_cell[i_cell];
675 }
676
677 amrex::Real global_lamdb = 0.;
678 if (use_global_debye_length) {
679 global_lamdb = global_debye_length_data[i_cell];
680 }
681
682 // Call the function in order to perform collisions
683 // If there are product species, mask, p_pair_indices_1/2, and
684 // p_pair_reaction_weight and p_product_data are filled here
685 binary_collision_functor(
686 cell_start_1, cell_half_1,
687 cell_half_1, cell_stop_1,
688 indices_1, indices_1,
689 soa_1, soa_1, get_position_1, get_position_1,
690 n1, n1, T1, T1, global_lamdb,
691 q1, q1, m1, m1, dt, dV*volume_factor(i_cell), coll_idx,
692 cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
693 p_pair_reaction_weight, p_product_data, engine);
694 }
695 );
696 ABLASTR_PROFILE_VAR_STOP(prof_loopOverCollisions);
697 // Create the new product particles and define their initial values
698 // num_added: how many particles of each product species have been created
699 //NOLINTNEXTLINE(readability-suspicious-call-argument)
700 const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
701 ptile_1, ptile_1,
702 product_species_vector,
703 tile_products_data,
704 m1, m1,
705 products_mass, p_mask, products_np,
706 copy_species1, copy_species2,
707 p_pair_indices_1, p_pair_indices_2,
708 p_pair_reaction_weight, p_product_data);
709
710 for (int i = 0; i < n_product_species; i++)
711 {
712 setNewParticleIDs(*(tile_products_data[i]), static_cast<int>(products_np[i]), num_added[i]);
713 }
714
715 if (m_correct_energy_momentum) {
716 // Loop over cells and calculate the change in energy and momentum.
717 amrex::ParticleReal const energy_fraction = m_energy_fraction;
718
719 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::correctEnergyMomentum", prof_correctEnergyMomentum);
720 amrex::ParallelFor( np1,
721 [=] AMREX_GPU_DEVICE (int i1) noexcept
722 {
723
724 const int i_cell = bins_1_ptr[i1];
725
726 // Subract the momentum from the initial values.
727 // Whatever is left over, will be distributed among the particles.
728 amrex::ParticleReal const px = w1[i1]*u1x[i1]*m1;
729 amrex::ParticleReal const py = w1[i1]*u1y[i1]*m1;
730 amrex::ParticleReal const pz = w1[i1]*u1z[i1]*m1;
731 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], -px);
732 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], -py);
733 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], -pz);
734 }
735 );
736
737 amrex::ParallelFor( np1,
738 [=] AMREX_GPU_DEVICE (int i1) noexcept
739 {
740
741 const int i_cell = bins_1_ptr[i1];
742
743 amrex::ParticleReal const ww_weighted_sum = ww_weighted_sum_in_each_cell[i_cell];
744 amrex::ParticleReal const w_factor = std::pow(w1[i1], beta_weight_exponent - 1._prt);
745 u1x[i1] += w_factor*px_in_each_cell[i_cell]/ww_weighted_sum;
746 u1y[i1] += w_factor*py_in_each_cell[i_cell]/ww_weighted_sum;
747 u1z[i1] += w_factor*pz_in_each_cell[i_cell]/ww_weighted_sum;
748
749 amrex::ParticleReal const KE_after = Algorithms::KineticEnergy(u1x[i1], u1y[i1], u1z[i1], m1);
750 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], -w1[i1]*KE_after);
751 }
752 );
753
754 amrex::Gpu::Buffer<amrex::Long> failed_corrections({0});
755 amrex::Long* failed_corrections_ptr = failed_corrections.data();
756 const int np_warning_threshold = m_np_warning_threshold;
757
758 // The particles are sorted by weight in decreasing order. This gives more of the correction
759 // to the heavier particles, having a smallar effect on their momenta, and typically
760 // requiring fewer particles to absorb or give off the extra energy.
761 const int energy_correction_sort_by_weight_flag = m_energy_correction_sort_by_weight ? sort : no_sort;
762 auto heapSortDecreasing = ParticleUtils::HeapSortDecreasing();
763
764 amrex::ParallelFor(amrex::TypeList<amrex::CompileTimeOptions<no_sort,sort>>{},
765 {energy_correction_sort_by_weight_flag},
766 n_cells,
767 [=] AMREX_GPU_DEVICE (int i_cell, auto energy_correction_sort_by_weight_control) noexcept
768 {
769
770 index_type const cell_start_1 = cell_offsets_1[i_cell];
771 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
772
773 // No correction needed if there is only one particle in the cell
774 if ( cell_stop_1 - cell_start_1 <= 1 ) { return; }
775
776 amrex::ParticleReal deltaEp1 = -KE_in_each_cell[i_cell];
777
778 if (deltaEp1 != 0.) {
779
780 // The ignore_unused is needed to fix the error "An extended __device__
781 // lambda cannot first-capture variable in constexpr-if context"
782 amrex::ignore_unused(heapSortDecreasing, indices_1, w1);
783 if constexpr (energy_correction_sort_by_weight_control == sort) {
784 int const numCell1 = (cell_stop_1 - cell_start_1);
785 heapSortDecreasing(indices_1, w1, cell_start_1, numCell1);
786 }
787
788 // Adjust species 1 particles to absorb deltaEp1.
789 const bool correction_failed =
790 ParticleUtils::ModifyEnergyPairwise(u1x, u1y, u1z, w1, indices_1,
791 cell_start_1, cell_stop_1, m1,
792 energy_fraction, deltaEp1);
793 if (correction_failed) {
794 // If the correction failed, give up on the collision and restore the
795 // momentum to what it was beforehand.
796 for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
797 u1x[ indices_1[i1] ] = u1x_before_ptr[ indices_1[i1] ];
798 u1y[ indices_1[i1] ] = u1y_before_ptr[ indices_1[i1] ];
799 u1z[ indices_1[i1] ] = u1z_before_ptr[ indices_1[i1] ];
800 }
801 int const numCell1 = (cell_stop_1 - cell_start_1);
802 if (numCell1 > np_warning_threshold) {
803 amrex::Gpu::Atomic::Add(failed_corrections_ptr, amrex::Long(1));
804 }
805 }
806 }
807
808 }
809 );
810
811 amrex::Long const num_failed_corrections = *(failed_corrections.copyToHost());
812 if (num_failed_corrections > 0) {
813 ablastr::warn_manager::WMRecordWarning("BinaryCollision::doCollisionsWithinTile",
814 "The energy correction failed for " + std::to_string(num_failed_corrections) + " cells " +
815 "for Coulomb collisions for species " + m_species_names[0] + ". " +
816 "The collisions in those cells was cancelled.");
817 }
818
819 ABLASTR_PROFILE_VAR_STOP(prof_correctEnergyMomentum);
820 }
821
822 }
823 else // species_1 != species_2
824 {
825 // Extract particles in the tile that `mfi` points to
826 ParticleTileType& ptile_1 = species_1.ParticlesAt(lev, mfi);
827 ParticleTileType& ptile_2 = species_2.ParticlesAt(lev, mfi);
828
829 // Find the particles that are in each cell of this tile
830 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findParticlesInEachCell", prof_findParticlesInEachCell);
831 ParticleBins bins_1 = findParticlesInEachCell( geom_lev, mfi, ptile_1 );
832 ParticleBins bins_2 = findParticlesInEachCell( geom_lev, mfi, ptile_2 );
833 ABLASTR_PROFILE_VAR_STOP(prof_findParticlesInEachCell);
834
835 // Loop over cells, and collide the particles in each cell
836
837 // Extract low-level data
838 auto const n_cells = static_cast<int>(bins_1.numBins());
839 // - Species 1
840 auto np1 = ptile_1.numParticles();
841 const auto soa_1 = ptile_1.getParticleTileData();
842 index_type* AMREX_RESTRICT indices_1 = bins_1.permutationPtr();
843 index_type const* AMREX_RESTRICT cell_offsets_1 = bins_1.offsetsPtr();
844 index_type const* AMREX_RESTRICT bins_1_ptr = bins_1.binsPtr();
845 const amrex::ParticleReal q1 = species_1.getCharge();
846 const amrex::ParticleReal m1 = species_1.getMass();
847 auto get_position_1 = GetParticlePosition<PIdx>(ptile_1, getpos_offset);
848 // - Species 2
849 auto np2 = ptile_2.numParticles();
850 const auto soa_2 = ptile_2.getParticleTileData();
851 index_type* AMREX_RESTRICT indices_2 = bins_2.permutationPtr();
852 index_type const* AMREX_RESTRICT cell_offsets_2 = bins_2.offsetsPtr();
853 index_type const* AMREX_RESTRICT bins_2_ptr = bins_2.binsPtr();
854 const amrex::ParticleReal q2 = species_2.getCharge();
855 const amrex::ParticleReal m2 = species_2.getMass();
856 auto get_position_2 = GetParticlePosition<PIdx>(ptile_2, getpos_offset);
857
858 /*
859 The following calculations are only required when creating product particles
860 */
861 const int n_cells_products = have_product_species ? n_cells: 0;
862 amrex::Gpu::DeviceVector<index_type> n_pairs_in_each_cell(n_cells_products);
863 index_type* AMREX_RESTRICT p_n_pairs_in_each_cell = n_pairs_in_each_cell.dataPtr();
864
865 // Compute how many pairs in each cell and store in n_pairs_in_each_cell array
866 // For different species, the number of pairs in a cell is the number of particles of
867 // the species that has the most particles in that cell
868 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::computeNumberOfPairs", prof_computeNumberOfPairs);
869 amrex::ParallelFor( n_cells_products,
870 [=] AMREX_GPU_DEVICE (int i_cell) noexcept
871 {
872 const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
873 const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
874 // Particular case: no pair if a species has no particle in that cell
875 if (n_part_in_cell_1 == 0 || n_part_in_cell_2 == 0) {
876 p_n_pairs_in_each_cell[i_cell] = 0;
877 } else {
878 p_n_pairs_in_each_cell[i_cell] =
879 amrex::max(n_part_in_cell_1,n_part_in_cell_2);
880 }
881 }
882 );
883
884 // Start indices of the pairs in a cell. Will be used for particle creation
885 amrex::Gpu::DeviceVector<index_type> pair_offsets(n_cells_products);
886 const index_type n_total_pairs = (n_cells_products == 0) ? 0:
887 amrex::Scan::ExclusiveSum(n_cells_products,
888 p_n_pairs_in_each_cell, pair_offsets.data());
889 index_type* AMREX_RESTRICT p_pair_offsets = pair_offsets.dataPtr();
890
891 amrex::Gpu::DeviceVector<index_type> n_ind_pairs_in_each_cell(n_cells+1);
892 index_type* AMREX_RESTRICT p_n_ind_pairs_in_each_cell = n_ind_pairs_in_each_cell.dataPtr();
893
894 amrex::ParallelFor( n_cells+1,
895 [=] AMREX_GPU_DEVICE (int i_cell) noexcept
896 {
897 if (i_cell < n_cells)
898 {
899 const auto n_part_in_cell_1 = cell_offsets_1[i_cell+1] - cell_offsets_1[i_cell];
900 const auto n_part_in_cell_2 = cell_offsets_2[i_cell+1] - cell_offsets_2[i_cell];
901 p_n_ind_pairs_in_each_cell[i_cell] = amrex::min(n_part_in_cell_1, n_part_in_cell_2);
902 }
903 else
904 {
905 p_n_ind_pairs_in_each_cell[i_cell] = 0;
906 }
907 }
908 );
909
910 // start indices of independent collisions.
911 amrex::Gpu::DeviceVector<index_type> coll_offsets(n_cells+1);
912 // number of total independent collision pairs
913 const auto n_independent_pairs = (int) amrex::Scan::ExclusiveSum(n_cells+1,
914 p_n_ind_pairs_in_each_cell, coll_offsets.data(), amrex::Scan::RetSum{true});
915 index_type* AMREX_RESTRICT p_coll_offsets = coll_offsets.dataPtr();
916 ABLASTR_PROFILE_VAR_STOP(prof_computeNumberOfPairs);
917
918 // mask: equal to 1 if particle creation occurs for a given pair, 0 otherwise
919 amrex::Gpu::DeviceVector<index_type> mask(n_total_pairs, index_type(0));
920 index_type* AMREX_RESTRICT p_mask = mask.dataPtr();
921 // Will be filled with the index of the first particle of a given pair
922 amrex::Gpu::DeviceVector<index_type> pair_indices_1(n_total_pairs);
923 index_type* AMREX_RESTRICT p_pair_indices_1 = pair_indices_1.dataPtr();
924 // Will be filled with the index of the second particle of a given pair
925 amrex::Gpu::DeviceVector<index_type> pair_indices_2(n_total_pairs);
926 index_type* AMREX_RESTRICT p_pair_indices_2 = pair_indices_2.dataPtr();
927 // How much weight should be given to the produced particles (and removed from the
928 // reacting particles)
929 amrex::Gpu::DeviceVector<amrex::ParticleReal> pair_reaction_weight(n_total_pairs);
930 amrex::ParticleReal* AMREX_RESTRICT p_pair_reaction_weight =
931 pair_reaction_weight.dataPtr();
932 // Extra data needed when products are created
933 int const n_product_data = (binary_collision_functor.m_need_product_data ? n_total_pairs : 0);
934 amrex::Gpu::DeviceVector<amrex::ParticleReal> product_data(n_product_data);
935 amrex::ParticleReal* AMREX_RESTRICT p_product_data = product_data.dataPtr();
936 /*
937 End of calculations only required when creating product particles
938 */
939
940 amrex::ParticleReal * const AMREX_RESTRICT w1 = soa_1.m_rdata[PIdx::w];
941 amrex::ParticleReal * const AMREX_RESTRICT u1x = soa_1.m_rdata[PIdx::ux];
942 amrex::ParticleReal * const AMREX_RESTRICT u1y = soa_1.m_rdata[PIdx::uy];
943 amrex::ParticleReal * const AMREX_RESTRICT u1z = soa_1.m_rdata[PIdx::uz];
944
945 amrex::ParticleReal * const AMREX_RESTRICT w2 = soa_2.m_rdata[PIdx::w];
946 amrex::ParticleReal * const AMREX_RESTRICT u2x = soa_2.m_rdata[PIdx::ux];
947 amrex::ParticleReal * const AMREX_RESTRICT u2y = soa_2.m_rdata[PIdx::uy];
948 amrex::ParticleReal * const AMREX_RESTRICT u2z = soa_2.m_rdata[PIdx::uz];
949
950 // create vectors to store density and temperature on cell level
951 amrex::Gpu::DeviceVector<amrex::ParticleReal> n1_vec, n2_vec;
952 amrex::Gpu::DeviceVector<amrex::ParticleReal> T1_vec, T2_vec;
953 amrex::Gpu::DeviceVector<amrex::ParticleReal> vx1_vec, vx2_vec;
954 amrex::Gpu::DeviceVector<amrex::ParticleReal> vy1_vec, vy2_vec;
955 amrex::Gpu::DeviceVector<amrex::ParticleReal> vz1_vec, vz2_vec;
956 amrex::Gpu::DeviceVector<amrex::ParticleReal> vs1_vec, vs2_vec;
957
958 if (binary_collision_functor.m_computeSpeciesDensities) {
959 n1_vec.resize(n_cells, 0.0_prt);
960 n2_vec.resize(n_cells, 0.0_prt);
961 }
962 if (binary_collision_functor.m_computeSpeciesTemperatures) {
963 T1_vec.resize(n_cells, 0.0_prt);
964 T2_vec.resize(n_cells, 0.0_prt);
965 vx1_vec.resize(n_cells, 0.0_prt);
966 vx2_vec.resize(n_cells, 0.0_prt);
967 vy1_vec.resize(n_cells, 0.0_prt);
968 vy2_vec.resize(n_cells, 0.0_prt);
969 vz1_vec.resize(n_cells, 0.0_prt);
970 vz2_vec.resize(n_cells, 0.0_prt);
971 vs1_vec.resize(n_cells, 0.0_prt);
972 vs2_vec.resize(n_cells, 0.0_prt);
973 }
974 amrex::ParticleReal* AMREX_RESTRICT n1_in_each_cell = n1_vec.dataPtr();
975 amrex::ParticleReal* AMREX_RESTRICT n2_in_each_cell = n2_vec.dataPtr();
976 amrex::ParticleReal* AMREX_RESTRICT T1_in_each_cell = T1_vec.dataPtr();
977 amrex::ParticleReal* AMREX_RESTRICT T2_in_each_cell = T2_vec.dataPtr();
978 amrex::ParticleReal* AMREX_RESTRICT vx1_in_each_cell = vx1_vec.dataPtr();
979 amrex::ParticleReal* AMREX_RESTRICT vx2_in_each_cell = vx2_vec.dataPtr();
980 amrex::ParticleReal* AMREX_RESTRICT vy1_in_each_cell = vy1_vec.dataPtr();
981 amrex::ParticleReal* AMREX_RESTRICT vy2_in_each_cell = vy2_vec.dataPtr();
982 amrex::ParticleReal* AMREX_RESTRICT vz1_in_each_cell = vz1_vec.dataPtr();
983 amrex::ParticleReal* AMREX_RESTRICT vz2_in_each_cell = vz2_vec.dataPtr();
984 amrex::ParticleReal* AMREX_RESTRICT vs1_in_each_cell = vs1_vec.dataPtr();
985 amrex::ParticleReal* AMREX_RESTRICT vs2_in_each_cell = vs2_vec.dataPtr();
986
987 // Create vectors to store energy and momentum in each cell.
988 // This is used to correct the energy and momentum in the
989 // cell after the collisions.
990 amrex::Gpu::DeviceVector<amrex::ParticleReal> ww_weighted_sum_vec;
991 amrex::Gpu::DeviceVector<amrex::ParticleReal> KE_vec;
992 amrex::Gpu::DeviceVector<amrex::ParticleReal> px_vec;
993 amrex::Gpu::DeviceVector<amrex::ParticleReal> py_vec;
994 amrex::Gpu::DeviceVector<amrex::ParticleReal> pz_vec;
995 if (m_correct_energy_momentum) {
996 ww_weighted_sum_vec.resize(n_cells, 0.0_prt);
997 KE_vec.resize(n_cells, 0.0_prt);
998 px_vec.resize(n_cells, 0.0_prt);
999 py_vec.resize(n_cells, 0.0_prt);
1000 pz_vec.resize(n_cells, 0.0_prt);
1001 }
1002 amrex::ParticleReal* AMREX_RESTRICT ww_weighted_sum_in_each_cell = ww_weighted_sum_vec.dataPtr();
1003 amrex::ParticleReal* AMREX_RESTRICT KE_in_each_cell = KE_vec.dataPtr();
1004 amrex::ParticleReal* AMREX_RESTRICT px_in_each_cell = px_vec.dataPtr();
1005 amrex::ParticleReal* AMREX_RESTRICT py_in_each_cell = py_vec.dataPtr();
1006 amrex::ParticleReal* AMREX_RESTRICT pz_in_each_cell = pz_vec.dataPtr();
1007
1008 // Save the momentum before the collision since sometimes the correction fails
1009 // and the only solution is to restore the pre-collision values.
1010 // The implicit solver already has ux_n etc so use those if available,
1011 // otherwise create temporaries.
1012 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1x_before;
1013 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1y_before;
1014 amrex::Gpu::DeviceVector<amrex::ParticleReal> u1z_before;
1015 amrex::ParticleReal* AMREX_RESTRICT u1x_before_ptr = nullptr;
1016 amrex::ParticleReal* AMREX_RESTRICT u1y_before_ptr = nullptr;
1017 amrex::ParticleReal* AMREX_RESTRICT u1z_before_ptr = nullptr;
1018 if (m_correct_energy_momentum) {
1019 // Check if ux_n was added.
1020 std::vector<std::string> const & real_names1 = species_1.GetRealSoANames();
1021 auto const pos1 = std::find(real_names1.begin(), real_names1.end(), "ux_n");
1022 if (pos1 != real_names1.end()) {
1023 int const n_builtin_real = PIdx::nattribs;
1024 auto const u1x_ni = static_cast<int>(std::distance(real_names1.begin(), pos1));
1025 int const u1x_runtime_ni = u1x_ni - n_builtin_real;
1026 u1x_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni];
1027 u1y_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 1];
1028 u1z_before_ptr = soa_1.m_runtime_rdata[u1x_runtime_ni + 2];
1029 } else {
1030 u1x_before.resize(np1);
1031 u1y_before.resize(np1);
1032 u1z_before.resize(np1);
1033 u1x_before_ptr = u1x_before.dataPtr();
1034 u1y_before_ptr = u1y_before.dataPtr();
1035 u1z_before_ptr = u1z_before.dataPtr();
1036 }
1037 }
1038
1039 amrex::Gpu::DeviceVector<amrex::ParticleReal> u2x_before;
1040 amrex::Gpu::DeviceVector<amrex::ParticleReal> u2y_before;
1041 amrex::Gpu::DeviceVector<amrex::ParticleReal> u2z_before;
1042 amrex::ParticleReal* AMREX_RESTRICT u2x_before_ptr = nullptr;
1043 amrex::ParticleReal* AMREX_RESTRICT u2y_before_ptr = nullptr;
1044 amrex::ParticleReal* AMREX_RESTRICT u2z_before_ptr = nullptr;
1045 if (m_correct_energy_momentum) {
1046 // Check if ux_n was added.
1047 std::vector<std::string> const & real_names2 = species_2.GetRealSoANames();
1048 auto const pos2 = std::find(real_names2.begin(), real_names2.end(), "ux_n");
1049 if (pos2 != real_names2.end()) {
1050 int const n_builtin_real = PIdx::nattribs;
1051 auto const u2x_ni = static_cast<int>(std::distance(real_names2.begin(), pos2));
1052 int const u2x_runtime_ni = u2x_ni - n_builtin_real;
1053 u2x_before_ptr = soa_2.m_runtime_rdata[u2x_runtime_ni];
1054 u2y_before_ptr = soa_2.m_runtime_rdata[u2x_runtime_ni + 1];
1055 u2z_before_ptr = soa_2.m_runtime_rdata[u2x_runtime_ni + 2];
1056 } else {
1057 u2x_before.resize(np2);
1058 u2y_before.resize(np2);
1059 u2z_before.resize(np2);
1060 u2x_before_ptr = u2x_before.dataPtr();
1061 u2y_before_ptr = u2y_before.dataPtr();
1062 u2z_before_ptr = u2z_before.dataPtr();
1063 }
1064 }
1065
1066 amrex::ParticleReal const beta_weight_exponent = m_beta_weight_exponent;
1067
1068 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures", prof_findDensityTemperatures);
1069
1070 if (m_correct_energy_momentum ||
1071 binary_collision_functor.m_computeSpeciesDensities ||
1072 binary_collision_functor.m_computeSpeciesTemperatures) {
1073
1074 bool const correct_energy_momentum = m_correct_energy_momentum;
1075
1076 // Loop over particles and compute quantities needed for energy conservation and the collsion parameters
1077 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::atomics", prof_findDensityTemperatures_atomics);
1078 amrex::ParallelFor( std::max(np1, np2),
1079 [=] AMREX_GPU_DEVICE (int ip) noexcept
1080 {
1081 if (correct_energy_momentum) {
1082 if (ip < np1) {
1083 const int i1 = ip;
1084 const int i_cell = bins_1_ptr[i1];
1085
1086 // Save the momenta before the collision
1087 u1x_before_ptr[i1] = u1x[i1];
1088 u1y_before_ptr[i1] = u1y[i1];
1089 u1z_before_ptr[i1] = u1z[i1];
1090
1091 // Compute the local energy and momentum.
1092 amrex::ParticleReal const w1_scaled = std::pow(w1[i1], beta_weight_exponent);
1093 amrex::Gpu::Atomic::AddNoRet(&ww_weighted_sum_in_each_cell[i_cell], w1_scaled*m1);
1094 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], w1[i1]*u1x[i1]*m1);
1095 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], w1[i1]*u1y[i1]*m1);
1096 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], w1[i1]*u1z[i1]*m1);
1097 const amrex::ParticleReal KE1 = Algorithms::KineticEnergy(u1x[i1], u1y[i1], u1z[i1], m1);
1098 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], w1[i1]*KE1);
1099 }
1100
1101 if (ip < np2) {
1102 const int i2 = ip;
1103 const int i_cell = bins_2_ptr[i2];
1104
1105 // Save the momenta before the collision
1106 u2x_before_ptr[i2] = u2x[i2];
1107 u2y_before_ptr[i2] = u2y[i2];
1108 u2z_before_ptr[i2] = u2z[i2];
1109
1110 // Compute the local energy and momentum.
1111 amrex::ParticleReal const w2_scaled = std::pow(w2[i2], beta_weight_exponent);
1112 amrex::Gpu::Atomic::AddNoRet(&ww_weighted_sum_in_each_cell[i_cell], w2_scaled*m2);
1113 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], w2[i2]*u2x[i2]*m2);
1114 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], w2[i2]*u2y[i2]*m2);
1115 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], w2[i2]*u2z[i2]*m2);
1116 const amrex::ParticleReal KE2 = Algorithms::KineticEnergy(u2x[i2], u2y[i2], u2z[i2], m2);
1117 amrex::Gpu::Atomic::AddNoRet(&KE_in_each_cell[i_cell], w2[i2]*KE2);
1118 }
1119 }
1120
1121 // compute local densities [1/m^3]
1122 if (binary_collision_functor.m_computeSpeciesDensities) {
1123 if (ip < np1) {
1124 amrex::Gpu::Atomic::AddNoRet(&n1_in_each_cell[bins_1_ptr[ip]],
1125 w1[ip]/(dV*volume_factor(bins_1_ptr[ip])));
1126 }
1127 if (ip < np2) {
1128 amrex::Gpu::Atomic::AddNoRet(&n2_in_each_cell[bins_2_ptr[ip]],
1129 w2[ip]/(dV*volume_factor(bins_2_ptr[ip])));
1130 }
1131 }
1132
1133 // compute local temperatures [Joules]
1134 if (binary_collision_functor.m_computeSpeciesTemperatures) {
1135 if (ip < np1) {
1136 const amrex::ParticleReal us = u1x[ip]*u1x[ip] + u1y[ip]*u1y[ip] + u1z[ip]*u1z[ip];
1137 auto constexpr inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
1138 const amrex::ParticleReal gm = std::sqrt( 1.0_prt + us*inv_c2);
1139 amrex::Gpu::Atomic::AddNoRet(&T1_in_each_cell[bins_1_ptr[ip]], w1[ip]);
1140 amrex::Gpu::Atomic::AddNoRet(&vx1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1x[ip]/gm);
1141 amrex::Gpu::Atomic::AddNoRet(&vy1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1y[ip]/gm);
1142 amrex::Gpu::Atomic::AddNoRet(&vz1_in_each_cell[bins_1_ptr[ip]], w1[ip]*u1z[ip]/gm);
1143 amrex::Gpu::Atomic::AddNoRet(&vs1_in_each_cell[bins_1_ptr[ip]], w1[ip]*us/gm/gm);
1144 }
1145 if (ip < np2) {
1146 const amrex::ParticleReal us = u2x[ip]*u2x[ip] + u2y[ip]*u2y[ip] + u2z[ip]*u2z[ip];
1147 auto constexpr inv_c2 = PhysConst::inv_c2_v<amrex::ParticleReal>;
1148 const amrex::ParticleReal gm = std::sqrt( 1.0_prt + us*inv_c2);
1149 amrex::Gpu::Atomic::AddNoRet(&T2_in_each_cell [bins_2_ptr[ip]], w2[ip]);
1150 amrex::Gpu::Atomic::AddNoRet(&vx2_in_each_cell[bins_2_ptr[ip]], w2[ip]*u2x[ip]/gm);
1151 amrex::Gpu::Atomic::AddNoRet(&vy2_in_each_cell[bins_2_ptr[ip]], w2[ip]*u2y[ip]/gm);
1152 amrex::Gpu::Atomic::AddNoRet(&vz2_in_each_cell[bins_2_ptr[ip]], w2[ip]*u2z[ip]/gm);
1153 amrex::Gpu::Atomic::AddNoRet(&vs2_in_each_cell[bins_2_ptr[ip]], w2[ip]*us/gm/gm);
1154 }
1155 }
1156 }
1157 );
1158 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_atomics);
1159
1160 }
1161
1162 // Finish temperature calculation - we launch 2*n_cells threads to compute both species simultaneously
1163 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::finishTemperature", prof_findDensityTemperatures_finish);
1164 amrex::ParallelFor( 2*n_cells, [=] AMREX_GPU_DEVICE (int i) noexcept
1165 {
1166 const int i_cell = i < n_cells ? i : i - n_cells;
1167
1168 // The particles from species1 that are in the cell `i_cell` are
1169 // given by the `indices_1[cell_start_1:cell_stop_1]`
1170 index_type const cell_start_1 = cell_offsets_1[i_cell];
1171 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1172
1173 // Same for species 2
1174 index_type const cell_start_2 = cell_offsets_2[i_cell];
1175 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1176
1177 // Do not need if one species is missing in the cell
1178 if ( cell_stop_1 - cell_start_1 < 1 ||
1179 cell_stop_2 - cell_start_2 < 1 ) { return; }
1180
1181 // finish temperature calculation if needed
1182 if (binary_collision_functor.m_computeSpeciesTemperatures) {
1183 if (i < n_cells) {
1184 const amrex::ParticleReal invsum1 = 1._prt/T1_in_each_cell[i_cell];
1185 auto vx1 = vx1_in_each_cell[i_cell] * invsum1;
1186 auto vy1 = vy1_in_each_cell[i_cell] * invsum1;
1187 auto vz1 = vz1_in_each_cell[i_cell] * invsum1;
1188 auto vs1 = vs1_in_each_cell[i_cell] * invsum1;
1189
1190 T1_in_each_cell[i_cell] = m1/(3._prt)*(vs1 -(vx1*vx1+vy1*vy1+vz1*vz1));
1191 } else {
1192 const amrex::ParticleReal invsum2 = 1._prt/T2_in_each_cell[i_cell];
1193 auto vx2 = vx2_in_each_cell[i_cell] * invsum2;
1194 auto vy2 = vy2_in_each_cell[i_cell] * invsum2;
1195 auto vz2 = vz2_in_each_cell[i_cell] * invsum2;
1196 auto vs2 = vs2_in_each_cell[i_cell] * invsum2;
1197
1198 T2_in_each_cell[i_cell] = m2/(3._prt)*(vs2 -(vx2*vx2+vy2*vy2+vz2*vz2));
1199 }
1200 }
1201 }
1202 );
1203 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_finish);
1204
1205 // shuffle - we launch 2*n_cells threads to compute both species simultaneously
1206 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::findDensityTemperatures::shuffle", prof_findDensityTemperatures_shuffle);
1207 amrex::ParallelForRNG( 2*n_cells,
1208 [=] AMREX_GPU_DEVICE (int i, amrex::RandomEngine const& engine) noexcept
1209 {
1210 const int i_cell = i < n_cells ? i : i - n_cells;
1211
1212 // The particles from species1 that are in the cell `i_cell` are
1213 // given by the `indices_1[cell_start_1:cell_stop_1]`
1214 index_type const cell_start_1 = cell_offsets_1[i_cell];
1215 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1216
1217 // Same for species 2
1218 index_type const cell_start_2 = cell_offsets_2[i_cell];
1219 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1220
1221 // Do not collide if one species is missing in the cell
1222 if ( cell_stop_1 - cell_start_1 < 1 ||
1223 cell_stop_2 - cell_start_2 < 1 ) { return; }
1224
1225 if (i < n_cells) {
1226 ShuffleFisherYates(indices_1, cell_start_1, cell_stop_1, engine);
1227 } else {
1228 ShuffleFisherYates(indices_2, cell_start_2, cell_stop_2, engine);
1229 }
1230 }
1231 );
1232 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures_shuffle);
1233
1234 ABLASTR_PROFILE_VAR_STOP(prof_findDensityTemperatures);
1235 // Loop over independent particle pairs
1236 // To speed up binary collisions on GPU, we try to expose as much parallelism
1237 // as possible (while avoiding race conditions): Instead of looping with one GPU
1238 // thread per cell, we loop with one GPU thread per "independent pairs" (i.e. pairs
1239 // that do not touch the same macroparticles, so that there is no race condition),
1240 // where the number of independent pairs is determined by the lower number of
1241 // macroparticles of either species, within each cell.
1242 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::LoopOverCollisions", prof_loopOverCollisions);
1243 amrex::ParallelForRNG( n_independent_pairs,
1244 [=] AMREX_GPU_DEVICE (int i_coll, amrex::RandomEngine const& engine) noexcept
1245 {
1246 // to avoid type mismatch errors
1247 auto ui_coll = (index_type)i_coll;
1248
1249 // Use a bisection algorithm to find the index of the cell in which this pair is located
1250 const int i_cell = amrex::bisect( p_coll_offsets, 0, n_cells, ui_coll );
1251
1252 // The particles from species1 that are in the cell `i_cell` are
1253 // given by the `indices_1[cell_start_1:cell_stop_1]`
1254 index_type const cell_start_1 = cell_offsets_1[i_cell];
1255 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1256 // Same for species 2
1257 index_type const cell_start_2 = cell_offsets_2[i_cell];
1258 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1259
1260 // collision number of the cell
1261 const index_type coll_idx = ui_coll - p_coll_offsets[i_cell];
1262
1263 // Same but for the pairs
1264 index_type const cell_start_pair = have_product_species?
1265 p_pair_offsets[i_cell]: 0;
1266
1267 // ux from species1 can be accessed like this:
1268 // ux_1[ indices_1[i] ], where i is between
1269 // cell_start_1 (inclusive) and cell_start_2 (exclusive)
1270
1271 // Get the local densities and temperatures for this cell
1272 amrex::ParticleReal n1 = 0.0, n2 = 0.0;
1273 amrex::ParticleReal T1 = 0.0, T2 = 0.0;
1274 if (binary_collision_functor.m_computeSpeciesDensities) {
1275 n1 = n1_in_each_cell[i_cell];
1276 n2 = n2_in_each_cell[i_cell];
1277 }
1278 if (binary_collision_functor.m_computeSpeciesTemperatures) {
1279 T1 = T1_in_each_cell[i_cell];
1280 T2 = T2_in_each_cell[i_cell];
1281 }
1282
1283 amrex::Real global_lamdb = 0.;
1284 if (use_global_debye_length) {
1285 global_lamdb = global_debye_length_data[i_cell];
1286 }
1287
1288 // Call the function in order to perform collisions
1289 // If there are product species, p_mask, p_pair_indices_1/2, and
1290 // p_pair_reaction_weight and p_product_data are filled here
1291 binary_collision_functor(
1292 cell_start_1, cell_stop_1, cell_start_2, cell_stop_2,
1293 indices_1, indices_2,
1294 soa_1, soa_2, get_position_1, get_position_2,
1295 n1, n2, T1, T2, global_lamdb,
1296 q1, q2, m1, m2, dt, dV*volume_factor(i_cell), coll_idx,
1297 cell_start_pair, p_mask, p_pair_indices_1, p_pair_indices_2,
1298 p_pair_reaction_weight, p_product_data, engine);
1299 }
1300 );
1301 ABLASTR_PROFILE_VAR_STOP(prof_loopOverCollisions);
1302
1303 // Create the new product particles and define their initial values
1304 // num_added: how many particles of each product species have been created
1305 //NOLINTNEXTLINE(readability-suspicious-call-argument)
1306 const amrex::Vector<int> num_added = m_copy_transform_functor(n_total_pairs,
1307 ptile_1, ptile_2,
1308 product_species_vector,
1309 tile_products_data,
1310 m1, m2,
1311 products_mass, p_mask, products_np,
1312 copy_species1, copy_species2,
1313 p_pair_indices_1, p_pair_indices_2,
1314 p_pair_reaction_weight, p_product_data);
1315
1316 for (int i = 0; i < n_product_species; i++)
1317 {
1318 setNewParticleIDs(*(tile_products_data[i]), static_cast<int>(products_np[i]), num_added[i]);
1319 }
1320
1321 if (m_correct_energy_momentum) {
1322 // Loop over cells and calculate the change in energy and momentum.
1323 amrex::ParticleReal const energy_fraction = m_energy_fraction;
1324
1325 ABLASTR_PROFILE_VAR("BinaryCollision::doCollisionsWithinTile::correctEnergyMomentum", prof_correctEnergyMomentum);
1326 amrex::ParallelFor( std::max(np1, np2),
1327 [=] AMREX_GPU_DEVICE (int ip) noexcept
1328 {
1329
1330 if (ip < np1) {
1331 const int i1 = ip;
1332 const int i_cell = bins_1_ptr[i1];
1333
1334 // Subract the momentum from the initial values.
1335 // Whatever is left over, will be distributed among the particles.
1336 amrex::ParticleReal const px = w1[i1]*u1x[i1]*m1;
1337 amrex::ParticleReal const py = w1[i1]*u1y[i1]*m1;
1338 amrex::ParticleReal const pz = w1[i1]*u1z[i1]*m1;
1339 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], -px);
1340 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], -py);
1341 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], -pz);
1342 }
1343
1344 if (ip < np2) {
1345 const int i2 = ip;
1346 const int i_cell = bins_2_ptr[i2];
1347
1348 // Subract the momentum from the initial values.
1349 // Whatever is left over, will be distributed among the particles.
1350 amrex::ParticleReal const px = w2[i2]*u2x[i2]*m2;
1351 amrex::ParticleReal const py = w2[i2]*u2y[i2]*m2;
1352 amrex::ParticleReal const pz = w2[i2]*u2z[i2]*m2;
1353 amrex::Gpu::Atomic::AddNoRet(&px_in_each_cell[i_cell], -px);
1354 amrex::Gpu::Atomic::AddNoRet(&py_in_each_cell[i_cell], -py);
1355 amrex::Gpu::Atomic::AddNoRet(&pz_in_each_cell[i_cell], -pz);
1356 }
1357 }
1358 );
1359
1360 amrex::ParallelFor( std::max(np1, np2),
1361 [=] AMREX_GPU_DEVICE (int ip) noexcept
1362 {
1363 if (ip < np1) {
1364 const int i1 = ip;
1365 const int i_cell = bins_1_ptr[i1];
1366
1367 amrex::ParticleReal const ww_weighted_sum = ww_weighted_sum_in_each_cell[i_cell];
1368 amrex::ParticleReal const w_factor = std::pow(w1[i1], beta_weight_exponent - 1._prt);
1369 u1x[i1] += w_factor*px_in_each_cell[i_cell]/ww_weighted_sum;
1370 u1y[i1] += w_factor*py_in_each_cell[i_cell]/ww_weighted_sum;
1371 u1z[i1] += w_factor*pz_in_each_cell[i_cell]/ww_weighted_sum;
1372 }
1373
1374 if (ip < np2) {
1375 const int i2 = ip;
1376 const int i_cell = bins_2_ptr[i2];
1377
1378 amrex::ParticleReal const ww_weighted_sum = ww_weighted_sum_in_each_cell[i_cell];
1379 amrex::ParticleReal const w_factor = std::pow(w2[i2], beta_weight_exponent - 1._prt);
1380 u2x[i2] += w_factor*px_in_each_cell[i_cell]/ww_weighted_sum;
1381 u2y[i2] += w_factor*py_in_each_cell[i_cell]/ww_weighted_sum;
1382 u2z[i2] += w_factor*pz_in_each_cell[i_cell]/ww_weighted_sum;
1383 }
1384 }
1385 );
1386
1387 amrex::Gpu::Buffer<amrex::Long> failed_corrections({0});
1388 amrex::Long* failed_corrections_ptr = failed_corrections.data();
1389 const int np_warning_threshold = m_np_warning_threshold;
1390
1391 // The particles are sorted by weight in decreasing order. This gives more of the correction
1392 // to the heavier particles, having a smallar effect on their momenta, and typically
1393 // requiring fewer particles to absorb or give off the extra energy.
1394 const int energy_correction_sort_by_weight_flag = m_energy_correction_sort_by_weight ? sort : no_sort;
1395 auto heapSortDecreasing = ParticleUtils::HeapSortDecreasing();
1396
1397 amrex::ParallelFor(amrex::TypeList<amrex::CompileTimeOptions<no_sort,sort>>{},
1398 {energy_correction_sort_by_weight_flag},
1399 n_cells,
1400 [=] AMREX_GPU_DEVICE (int i_cell, auto energy_correction_sort_by_weight_control) noexcept
1401 {
1402 // The particles from species1 that are in the cell `i_cell` are
1403 // given by the `indices_1[cell_start_1:cell_stop_1]`.
1404 index_type const cell_start_1 = cell_offsets_1[i_cell];
1405 index_type const cell_stop_1 = cell_offsets_1[i_cell+1];
1406 // Same for species 2
1407 index_type const cell_start_2 = cell_offsets_2[i_cell];
1408 index_type const cell_stop_2 = cell_offsets_2[i_cell+1];
1409
1410 // Do not do the rebalance if one species is missing in the cell.
1411 if ( cell_stop_1 - cell_start_1 < 1 ||
1412 cell_stop_2 - cell_start_2 < 1 ) { return; }
1413
1414 amrex::ParticleReal const psq = px_in_each_cell[i_cell]*px_in_each_cell[i_cell] +
1415 py_in_each_cell[i_cell]*py_in_each_cell[i_cell] +
1416 pz_in_each_cell[i_cell]*pz_in_each_cell[i_cell];
1417 if (psq > 0.) {
1418 amrex::ParticleReal w1_sum = 0.0_prt;
1419 amrex::ParticleReal w2_sum = 0.0_prt;
1420 amrex::ParticleReal KE1_after = 0._prt;
1421 amrex::ParticleReal KE2_after = 0._prt;
1422 for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
1423 w1_sum += w1[ indices_1[i1] ];
1424 KE1_after += w1[ indices_1[i1] ]*Algorithms::KineticEnergy(u1x[ indices_1[i1] ],
1425 u1y[ indices_1[i1] ],
1426 u1z[ indices_1[i1] ], m1);
1427 }
1428 for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) {
1429 w2_sum += w2[ indices_2[i2] ];
1430 KE2_after += w2[ indices_2[i2] ]*Algorithms::KineticEnergy(u2x[ indices_2[i2] ],
1431 u2y[ indices_2[i2] ],
1432 u2z[ indices_2[i2] ], m2);
1433 }
1434
1435 const amrex::ParticleReal KE_after = KE1_after + KE2_after;
1436 const amrex::ParticleReal deltaE = KE_after - KE_in_each_cell[i_cell];
1437 amrex::ParticleReal deltaEp1, deltaEp2;
1438 int const numCell1 = (cell_stop_1 - cell_start_1);
1439 int const numCell2 = (cell_stop_2 - cell_start_2);
1440 if (numCell1 == 1) {
1441 deltaEp1 = 0.0;
1442 deltaEp2 = deltaE;
1443 } else if (numCell2 == 1) {
1444 deltaEp1 = deltaE;
1445 deltaEp2 = 0.0;
1446 } else {
1447 const amrex::ParticleReal Etotdenom = w1_sum*KE1_after/numCell1 + w2_sum*KE2_after/numCell2;
1448 deltaEp1 = w1_sum*KE1_after/Etotdenom*deltaE/numCell1;
1449 deltaEp2 = w2_sum*KE2_after/Etotdenom*deltaE/numCell2;
1450 }
1451
1452 // The ignore_unused is needed to fix the error "An extended __device__
1453 // lambda cannot first-capture variable in constexpr-if context"
1454 amrex::ignore_unused(heapSortDecreasing);
1455
1456 bool correction1_failed = false;
1457 bool correction2_failed = false;
1458 if (deltaEp1 != 0. && numCell1 > 1) {
1459 if constexpr (energy_correction_sort_by_weight_control == sort) {
1460 heapSortDecreasing(indices_1, w1, cell_start_1, numCell1);
1461 }
1462 // Adjust species 1 particles to absorb deltaEp1.
1463 correction1_failed =
1464 ParticleUtils::ModifyEnergyPairwise(u1x, u1y, u1z, w1, indices_1,
1465 cell_start_1, cell_stop_1, m1,
1466 energy_fraction, deltaEp1);
1467 }
1468
1469 if (deltaEp2 != 0. && numCell2 > 1) {
1470 if constexpr (energy_correction_sort_by_weight_control == sort) {
1471 heapSortDecreasing(indices_2, w2, cell_start_2, numCell2);
1472 }
1473 // Adjust species 2 particles to absorb deltaEp2.
1474 correction2_failed =
1475 ParticleUtils::ModifyEnergyPairwise(u2x, u2y, u2z, w2, indices_2,
1476 cell_start_2, cell_stop_2, m2,
1477 energy_fraction, deltaEp2);
1478 }
1479
1480 if (correction1_failed || correction2_failed ||
1481 (numCell1 == 1 && numCell2 == 1)) {
1482 // If the correction failed, or if each species has only one particle
1483 // and the correction can not be applied, give up on the collision and
1484 // restore the momentum to what it was beforehand.
1485 for (index_type i1=cell_start_1; i1<cell_stop_1; ++i1) {
1486 u1x[ indices_1[i1] ] = u1x_before_ptr[ indices_1[i1] ];
1487 u1y[ indices_1[i1] ] = u1y_before_ptr[ indices_1[i1] ];
1488 u1z[ indices_1[i1] ] = u1z_before_ptr[ indices_1[i1] ];
1489 }
1490 for (index_type i2=cell_start_2; i2<cell_stop_2; ++i2) {
1491 u2x[ indices_2[i2] ] = u2x_before_ptr[ indices_2[i2] ];
1492 u2y[ indices_2[i2] ] = u2y_before_ptr[ indices_2[i2] ];
1493 u2z[ indices_2[i2] ] = u2z_before_ptr[ indices_2[i2] ];
1494 }
1495 if ((correction1_failed && numCell1 > np_warning_threshold) ||
1496 (correction2_failed && numCell2 > np_warning_threshold)) {
1497 amrex::Gpu::Atomic::Add(failed_corrections_ptr, amrex::Long(1));
1498 }
1499 }
1500
1501 }
1502 }
1503 );
1504
1505 amrex::Long const num_failed_corrections = *(failed_corrections.copyToHost());
1506 if (num_failed_corrections > 0) {
1507 ablastr::warn_manager::WMRecordWarning("BinaryCollision::doCollisionsWithinTile",
1508 "The energy correction failed for " + std::to_string(num_failed_corrections) + " cells " +
1509 "for Coulomb collisions between species " + m_species_names[0] + " and " + m_species_names[1] + ". " +
1510 "The collisions in those cells was cancelled.");
1511 }
1512
1513 ABLASTR_PROFILE_VAR_STOP(prof_correctEnergyMomentum);
1514 }
1515
1516 } // end if ( m_isSameSpecies)
1517
1518 }
1519
1520private:
1521
1522 bool m_isSameSpecies;
1523 bool m_have_product_species;
1524
1525 bool m_correct_energy_momentum = false;
1526 bool m_energy_correction_sort_by_weight = false;
1527 int m_np_warning_threshold = 20;
1528
1529 amrex::ParticleReal m_beta_weight_exponent = 1.;
1530 amrex::ParticleReal m_energy_fraction = 0.05;
1531
1532 enum energy_correction_sort_by_weight_flags : int { no_sort, sort };
1533
1534 amrex::Vector<std::string> m_product_species;
1535 // functor that performs collisions within a cell
1536 CollisionFunctor m_binary_collision_functor;
1537 // functor that creates new particles and initializes their parameters
1538 CopyTransformFunctor m_copy_transform_functor;
1539
1540};
1541
1542#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:17
CollisionType::PairwiseCoulomb
@ PairwiseCoulomb
Definition BinaryCollisionUtils.H:24
CollisionType::DSMC
@ DSMC
Definition BinaryCollisionUtils.H:23
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
Definition WarpXAlgorithmSelection.H:119
WarpXParticleContainer_fwd.H
BinaryCollision::m_have_product_species
bool m_have_product_species
Definition BinaryCollision.H:1523
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:1538
BinaryCollision::m_energy_correction_sort_by_weight
bool m_energy_correction_sort_by_weight
Definition BinaryCollision.H:1526
BinaryCollision::m_isSameSpecies
bool m_isSameSpecies
Definition BinaryCollision.H:1522
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:1527
BinaryCollision::m_beta_weight_exponent
amrex::ParticleReal m_beta_weight_exponent
Definition BinaryCollision.H:1529
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:1530
BinaryCollision::m_correct_energy_momentum
bool m_correct_energy_momentum
Definition BinaryCollision.H:1525
BinaryCollision::m_product_species
amrex::Vector< std::string > m_product_species
Definition BinaryCollision.H:1534
BinaryCollision::energy_correction_sort_by_weight_flags
energy_correction_sort_by_weight_flags
Definition BinaryCollision.H:1532
BinaryCollision::no_sort
@ no_sort
Definition BinaryCollision.H:1532
BinaryCollision::sort
@ sort
Definition BinaryCollision.H:1532
BinaryCollision::m_binary_collision_functor
CollisionFunctor m_binary_collision_functor
Definition BinaryCollision.H:1536
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:404
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:3387
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:3220
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::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
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:25
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::bisect
__host__ __device__ T bisect(T lo, T hi, F f, T tol=1e-12, int max_iter=100)
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::min
__host__ __device__ constexpr const T & min(const T &a, const T &b) noexcept
amrex::second
double second() noexcept
amrex::max
__host__ __device__ constexpr const T & max(const T &a, const T &b) 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.13.2