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