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