Block-Structured AMR Software Framework
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
AMReX_ParticleContainerI.H
Go to the documentation of this file.
2
3#include <algorithm>
4#include <iterator>
5#include <set>
6#include <stdexcept>
7#include <string>
8#include <type_traits>
9#include <vector>
10
11namespace amrex {
12
13template <typename ParticleType, int NArrayReal, int NArrayInt,
14 template<class> class Allocator, class CellAssignor>
15void
17{
18 num_real_comm_comps = 0;
19 int comm_comps_start = 0;
20 if constexpr (!ParticleType::is_soa_particle) {
21 comm_comps_start += AMREX_SPACEDIM + NStructReal;
22 }
23 for (int i = comm_comps_start; i < comm_comps_start + NumRealComps(); ++i) {
24 if (h_redistribute_real_comp[i]) {++num_real_comm_comps;}
25 }
26
27 num_int_comm_comps = 0;
28 comm_comps_start = 2 + NStructInt;
29 for (int i = comm_comps_start; i < comm_comps_start + NumIntComps(); ++i) {
30 if (h_redistribute_int_comp[i]) {++num_int_comm_comps;}
31 }
32
33 if constexpr (ParticleType::is_soa_particle) {
34 particle_size = sizeof(uint64_t); // idcpu
35 } else {
36 particle_size = sizeof(ParticleType);
37 }
38 superparticle_size = particle_size +
39 num_real_comm_comps*sizeof(ParticleReal) + num_int_comm_comps*sizeof(int);
40}
41
42template <typename ParticleType, int NArrayReal, int NArrayInt,
43 template<class> class Allocator, class CellAssignor>
44void
46{
47 levelDirectoriesCreated = false;
48 usePrePost = false;
49 doUnlink = true;
50
51 // by default communicate all components
52 if constexpr (ParticleType::is_soa_particle)
53 {
54 h_redistribute_real_comp.resize(NArrayReal, true);
55 } else {
56 h_redistribute_real_comp.resize(AMREX_SPACEDIM + NStructReal + NArrayReal, true);
57 }
58 h_redistribute_int_comp.resize(2 + NStructInt + NArrayInt, true);
59
60 SetParticleSize();
61
62 // add default names for SoA Real and Int compile-time arguments
63 m_soa_rdata_names.clear();
64 for (int i=0; i<NArrayReal; ++i)
65 {
66 m_soa_rdata_names.push_back(getDefaultCompNameReal<ParticleType>(i));
67 }
68 m_soa_idata_names.clear();
69 for (int i=0; i<NArrayInt; ++i)
70 {
71 m_soa_idata_names.push_back(getDefaultCompNameInt<ParticleType>(i));
72 }
73
74 static bool initialized = false;
75 if ( ! initialized)
76 {
77 static_assert(sizeof(ParticleType)%sizeof(RealType) == 0,
78 "sizeof ParticleType is not a multiple of sizeof RealType");
79
80 ParmParse pp("particles");
81 pp.queryAdd("do_tiling", do_tiling);
82 Vector<int> tilesize(AMREX_SPACEDIM);
83 if (pp.queryarr("tile_size", tilesize, 0, AMREX_SPACEDIM)) {
84 for (int i=0; i<AMREX_SPACEDIM; ++i) { tile_size[i] = tilesize[i]; }
85 }
86
87 static_assert(std::is_standard_layout<ParticleType>::value,
88 "Particle type must be standard layout");
89 // && std::is_trivial<ParticleType>::value,
90 // "Particle type must be standard layout and trivial.");
91
92 pp.query("use_prepost", usePrePost);
93 pp.query("do_unlink", doUnlink);
94 pp.queryAdd("do_mem_efficient_sort", memEfficientSort);
95 pp.queryAdd("use_comms_arena", use_comms_arena);
96
97 initialized = true;
98 }
99}
100
101template<
102 typename ParticleType,
103 int NArrayReal,
104 int NArrayInt,
105 template<class> class Allocator,
106 class CellAssignor
107>
108void
110 std::vector<std::string> const & rdata_name, std::vector<std::string> const & idata_name
111)
112{
113 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rdata_name.size() == NArrayReal, "rdata_name must be equal to NArrayReal");
114 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(idata_name.size() == NArrayInt, "idata_name must be equal to NArrayInt");
115
116 // ensure names for components are unique
117 std::set<std::string> const unique_r_names(rdata_name.begin(), rdata_name.end());
118 std::set<std::string> const unique_i_names(idata_name.begin(), idata_name.end());
119 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(rdata_name.size() == unique_r_names.size(), "SetSoACompileTimeNames: Provided names in rdata_name are not unique!");
120 AMREX_ALWAYS_ASSERT_WITH_MESSAGE(idata_name.size() == unique_i_names.size(), "SetSoACompileTimeNames: Provided names in idata_name are not unique!");
121
122 for (int i=0; i<NArrayReal; ++i)
123 {
124 m_soa_rdata_names.at(i) = rdata_name.at(i);
125 }
126 for (int i=0; i<NArrayInt; ++i)
127 {
128 m_soa_idata_names.at(i) = idata_name.at(i);
129 }
130}
131
132template<
133 typename ParticleType,
134 int NArrayReal,
135 int NArrayInt,
136 template<class> class Allocator,
137 class CellAssignor
138>
139bool
141{
142 return std::find(m_soa_rdata_names.begin(), m_soa_rdata_names.end(), name) != std::end(m_soa_rdata_names);
143}
144
145template <typename ParticleType, int NArrayReal, int NArrayInt,
146 template<class> class Allocator, class CellAssignor>
147bool
149{
150 return std::find(m_soa_idata_names.begin(), m_soa_idata_names.end(), name) != std::end(m_soa_idata_names);
151}
152
153template<
154 typename ParticleType,
155 int NArrayReal,
156 int NArrayInt,
157 template<class> class Allocator,
158 class CellAssignor
159>
160int
162{
163 const auto it = std::find(m_soa_rdata_names.begin(), m_soa_rdata_names.end(), name);
164
165 if (it == m_soa_rdata_names.end())
166 {
167 throw std::runtime_error("GetRealCompIndex: Component " + name + " does not exist!");
168 }
169 else
170 {
171 return std::distance(m_soa_rdata_names.begin(), it);
172 }
173}
174
175template<
176 typename ParticleType,
177 int NArrayReal,
178 int NArrayInt,
179 template<class> class Allocator,
180 class CellAssignor
181>
182int
184{
185 const auto it = std::find(m_soa_idata_names.begin(), m_soa_idata_names.end(), name);
186
187 if (it == m_soa_idata_names.end())
188 {
189 throw std::runtime_error("GetIntCompIndex: Component " + name + " does not exist!");
190 }
191 else
192 {
193 return std::distance(m_soa_idata_names.begin(), it);
194 }
195}
196
197template <typename ParticleType, int NArrayReal, int NArrayInt,
198 template<class> class Allocator, class CellAssignor>
199template <typename P, typename Assignor>
202{
203 const Geometry& geom = Geom(lev);
204 const auto& domain = geom.Domain();
205 const auto& plo = geom.ProbLoArray();
206 const auto& dxi = geom.InvCellSizeArray();
207
208 return Assignor{}(p, plo, dxi, domain);
209}
210
211template <typename ParticleType, int NArrayReal, int NArrayInt,
212 template<class> class Allocator, class CellAssignor>
213template <typename P>
214bool
216::Where (const P& p,
217 ParticleLocData& pld,
218 int lev_min,
219 int lev_max,
220 int nGrow,
221 int local_grid) const
222{
223
224 AMREX_ASSERT(m_gdb != nullptr);
225
226 if (lev_max == -1) {
227 lev_max = finestLevel();
228 }
229
230 AMREX_ASSERT(lev_max <= finestLevel());
231
232 AMREX_ASSERT(nGrow == 0 || (nGrow >= 0 && lev_min == lev_max));
233
234 std::vector< std::pair<int, Box> > isects;
235
236 for (int lev = lev_max; lev >= lev_min; lev--) {
237 const IntVect& iv = Index(p, lev);
238 if (lev == pld.m_lev) {
239 // The fact that we are here means this particle does not belong to any finer grids.
240 if (pld.m_grid >= 0) {
241 if (pld.m_grown_gridbox.contains(iv)) {
242 pld.m_cell = iv;
243 if (!pld.m_tilebox.contains(iv)) {
244 pld.m_tile = getTileIndex(iv, pld.m_gridbox, do_tiling, tile_size, pld.m_tilebox);
245 }
246 return true;
247 }
248 }
249 }
250
251 int grid;
252 const BoxArray& ba = ParticleBoxArray(lev);
254
255 if (local_grid < 0) {
256 bool findfirst = (nGrow == 0) ? true : false;
257 ba.intersections(Box(iv, iv), isects, findfirst, nGrow);
258 grid = isects.empty() ? -1 : isects[0].first;
259 if (nGrow > 0 && isects.size() > 1) {
260 for (auto & isect : isects) {
261 Box bx = ba.getCellCenteredBox(isect.first);
262 for (int dir = 0; dir < AMREX_SPACEDIM; ++dir) {
263 Box gbx = bx;
265 gr[dir] = nGrow;
266 gbx.grow(gr);
267 if (gbx.contains(iv)) {
268 grid = isect.first;
269 }
270 }
271 }
272 }
273 } else {
274 grid = (*redistribute_mask_ptr)[local_grid](iv, 0);
275 }
276
277 if (grid >= 0) {
278 const Box& bx = ba.getCellCenteredBox(grid);
279 pld.m_lev = lev;
280 pld.m_grid = grid;
281 pld.m_tile = getTileIndex(iv, bx, do_tiling, tile_size, pld.m_tilebox);
282 pld.m_cell = iv;
283 pld.m_gridbox = bx;
284 pld.m_grown_gridbox = amrex::grow(bx, nGrow);
285 return true;
286 }
287 }
288
289 return false;
290}
291
292template <typename ParticleType, int NArrayReal, int NArrayInt,
293 template<class> class Allocator, class CellAssignor>
294template <typename P>
295bool
298 ParticleLocData& pld,
299 int lev_min,
300 int lev_max,
301 int local_grid) const
302{
303
304 AMREX_ASSERT(m_gdb != nullptr);
305
306 if (!Geom(0).isAnyPeriodic()) { return false; }
307
308 if (lev_max == -1) {
309 lev_max = finestLevel();
310 }
311
312 AMREX_ASSERT(lev_max <= finestLevel());
313
314 // Create a copy "dummy" particle to check for periodic outs.
315 Particle<0, 0> p_prime;
316 AMREX_D_TERM(p_prime.pos(0) = p.pos(0);,
317 p_prime.pos(1) = p.pos(1);,
318 p_prime.pos(2) = p.pos(2));
319 if (PeriodicShift(p_prime)) {
320 std::vector< std::pair<int,Box> > isects;
321 for (int lev = lev_max; lev >= lev_min; lev--) {
322
323 int grid;
324 IntVect iv;
325 const BoxArray& ba = ParticleBoxArray(lev);
327
328 if (local_grid < 0) {
329 iv = Index<amrex::Particle<0, 0>, DefaultAssignor>(p_prime, lev);
330 ba.intersections(Box(iv, iv), isects, true, 0);
331 grid = isects.empty() ? -1 : isects[0].first;
332 } else {
333 iv = Index<amrex::Particle<0, 0>, DefaultAssignor>(p_prime, lev);
334 if (ba[local_grid].contains(iv))
335 {
336 grid = local_grid;
337 }
338 else
339 {
340 ba.intersections(Box(iv, iv), isects, true, 0);
341 grid = isects.empty() ? -1 : isects[0].first;
342 if(grid == -1)
343 {
344 grid = (*redistribute_mask_ptr)[local_grid](Index(p, lev), 0);
345 }
346 }
347 }
348
349 if (grid >= 0) {
350 AMREX_D_TERM(p.pos(0) = p_prime.pos(0);,
351 p.pos(1) = p_prime.pos(1);,
352 p.pos(2) = p_prime.pos(2););
353
354 const Box& bx = ba.getCellCenteredBox(grid);
355
356 pld.m_lev = lev;
357 pld.m_grid = grid;
358 pld.m_tile = getTileIndex(iv, bx, do_tiling, tile_size, pld.m_tilebox);
359 pld.m_cell = iv;
360 pld.m_gridbox = bx;
361 pld.m_grown_gridbox = bx;
362 return true;
363 }
364 }
365 }
366
367 return false;
369
370
371template <typename ParticleType, int NArrayReal, int NArrayInt,
372 template<class> class Allocator, class CellAssignor>
373template <typename P>
374bool
376::PeriodicShift (P& p) const
377{
378 const auto& geom = Geom(0);
379 const auto plo = geom.ProbLoArray();
380 const auto phi = geom.ProbHiArray();
381 const auto rlo = geom.ProbLoArrayInParticleReal();
382 const auto rhi = geom.ProbHiArrayInParticleReal();
383 const auto is_per = geom.isPeriodicArray();
384
385 return enforcePeriodic(p, plo, phi, rlo, rhi, is_per);
386}
387
388template <typename ParticleType, int NArrayReal, int NArrayInt,
389 template<class> class Allocator, class CellAssignor>
393 bool /*update*/,
394 bool verbose,
395 ParticleLocData pld) const
396{
397 AMREX_ASSERT(m_gdb != nullptr);
398
399 bool ok = Where(p, pld);
400
401 if (!ok && Geom(0).isAnyPeriodic())
402 {
403 // Attempt to shift the particle back into the domain if it
404 // crossed a periodic boundary.
405 PeriodicShift(p);
406 ok = Where(p, pld);
407 }
408
409 if (!ok) {
410 // invalidate the particle.
411 if (verbose) {
412 amrex::AllPrint()<< "Invalidating out-of-domain particle: " << p << '\n';
413 }
414
415 AMREX_ASSERT(p.id() > 0);
416
417 p.id() = -p.id();
418 }
419
420 return pld;
421}
422
423template <typename ParticleType, int NArrayReal, int NArrayInt,
424 template<class> class Allocator, class CellAssignor>
425void
431
432template <typename ParticleType, int NArrayReal, int NArrayInt,
433 template<class> class Allocator, class CellAssignor>
434void
436{
438 int nlevs = std::max(0, finestLevel()+1);
439 m_particles.resize(nlevs);
440}
441
442template <typename ParticleType, int NArrayReal, int NArrayInt,
443 template<class> class Allocator, class CellAssignor>
444template <typename P>
445void
447 int lev_min, int lev_max, int nGrow, int local_grid) const
448{
449 bool success;
450 if (Geom(0).outsideRoundoffDomain(AMREX_D_DECL(p.pos(0), p.pos(1), p.pos(2))))
451 {
452 // Note that EnforcePeriodicWhere may shift the particle if it is successful.
453 success = EnforcePeriodicWhere(p, pld, lev_min, lev_max, local_grid);
454 if (!success && lev_min == 0)
456 // The particle has left the domain; invalidate it.
457 p.id() = -p.id();
458 success = true;
459 }
460 }
461 else
462 {
463 success = Where(p, pld, lev_min, lev_max, 0, local_grid);
464 }
465
466 if (!success)
467 {
468 success = (nGrow > 0) && Where(p, pld, lev_min, lev_min, nGrow);
469 pld.m_grown_gridbox = pld.m_gridbox; // reset grown box for subsequent calls.
470 }
471
472 if (!success)
473 {
474 amrex::Abort("ParticleContainer::locateParticle(): invalid particle.");
475 }
476}
477
478template <typename ParticleType, int NArrayReal, int NArrayInt,
479 template<class> class Allocator, class CellAssignor>
480Long
482{
483 Long nparticles = 0;
484 for (int lev = 0; lev <= finestLevel(); lev++) {
485 nparticles += NumberOfParticlesAtLevel(lev,only_valid,true);
486 }
487 if (!only_local) {
489 }
490 return nparticles;
491}
492
493template <typename ParticleType, int NArrayReal, int NArrayInt,
494 template<class> class Allocator, class CellAssignor>
497{
498 AMREX_ASSERT(lev >= 0 && lev < int(m_particles.size()));
499
500 LayoutData<Long> np_per_grid_local(ParticleBoxArray(lev),
501 ParticleDistributionMap(lev));
502
503 for (ParConstIterType pti(*this, lev); pti.isValid(); ++pti)
504 {
505 int gid = pti.index();
506 if (only_valid)
507 {
508 const auto& ptile = ParticlesAt(lev, pti);
509 const int np = ptile.numParticles();
510 auto const ptd = ptile.getConstParticleTileData();
512 ReduceOps<ReduceOpSum> reduce_op;
513 ReduceData<int> reduce_data(reduce_op);
514 using ReduceTuple = typename decltype(reduce_data)::Type;
515
516 reduce_op.eval(np, reduce_data,
517 [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple
518 {
519 return (ptd.id(i) > 0) ? 1 : 0;
520 });
521
522 int np_valid = amrex::get<0>(reduce_data.value(reduce_op));
523 np_per_grid_local[gid] += np_valid;
524 } else
525 {
526 np_per_grid_local[gid] += pti.numParticles();
528 }
529
530 Vector<Long> nparticles(np_per_grid_local.size(), 0);
531 if (only_local)
533 for (ParConstIterType pti(*this, lev); pti.isValid(); ++pti)
534 {
535 nparticles[pti.index()] = np_per_grid_local[pti.index()];
536 }
537 }
538 else
539 {
540 ParallelDescriptor::GatherLayoutDataToVector(np_per_grid_local, nparticles,
542 ParallelDescriptor::Bcast(nparticles.data(), nparticles.size(),
544 }
545
546 return nparticles;
547}
548
549template <typename ParticleType, int NArrayReal, int NArrayInt,
550 template<class> class Allocator, class CellAssignor>
553 Long nparticles = 0;
555 if (level < 0 || level >= int(m_particles.size())) { return nparticles; }
557 if (only_valid) {
559 ReduceData<unsigned long long> reduce_data(reduce_op);
560 using ReduceTuple = typename decltype(reduce_data)::Type;
561
562 for (const auto& kv : GetParticles(level)) {
563 const auto& ptile = kv.second;
564 auto const ptd = ptile.getConstParticleTileData();
565
566 reduce_op.eval(ptile.numParticles(), reduce_data,
567 [=] AMREX_GPU_DEVICE (int i) -> ReduceTuple
568 {
569 return (ptd.id(i) > 0) ? 1 : 0;
570 });
571 }
573 nparticles = static_cast<Long>(amrex::get<0>(reduce_data.value(reduce_op)));
574 }
575 else {
576 for (const auto& kv : GetParticles(level)) {
577 const auto& ptile = kv.second;
578 nparticles += ptile.numParticles();
579 }
580 }
581
582 if (!only_local) {
584 }
585
586 return nparticles;
587}
588
589//
590// This includes both valid and invalid particles since invalid particles still take up space.
591//
593template <typename ParticleType, int NArrayReal, int NArrayInt,
594 template<class> class Allocator, class CellAssignor>
595std::array<Long, 3>
598{
599 Long cnt = 0;
600
601 for (unsigned lev = 0; lev < m_particles.size(); lev++) {
602 const auto& pmap = m_particles[lev];
603 for (const auto& kv : pmap) {
604 const auto& ptile = kv.second;
605 cnt += ptile.numParticles();
606 }
607 }
608
609 Long mn = cnt, mx = mn;
610
611 const int IOProc = ParallelContext::IOProcessorNumberSub();
612 const Long sz = sizeof(ParticleType)+NumRealComps()*sizeof(ParticleReal)+NumIntComps()*sizeof(int);
614#ifdef AMREX_LAZY
615 Lazy::QueueReduction( [=] () mutable {
616#endif
620
621 amrex::Print() << "ParticleContainer spread across MPI nodes - bytes (num particles): [Min: "
622 << mn*sz
623 << " (" << mn << ")"
624 << ", Max: "
625 << mx*sz
626 << " (" << mx << ")"
627 << ", Total: "
628 << cnt*sz
629 << " (" << cnt << ")]\n";
630#ifdef AMREX_LAZY
631 });
632#endif
633
634 return {mn*sz, mx*sz, cnt*sz};
635}
636
637template <typename ParticleType, int NArrayReal, int NArrayInt,
638 template<class> class Allocator, class CellAssignor>
639std::array<Long, 3>
640ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
641::PrintCapacity () const
642{
643 Long cnt = 0;
644
645 for (unsigned lev = 0; lev < m_particles.size(); lev++) {
646 const auto& pmap = m_particles[lev];
647 for (const auto& kv : pmap) {
648 const auto& ptile = kv.second;
649 cnt += ptile.capacity();
651 }
652
653 Long mn = cnt, mx = mn;
654
655 const int IOProc = ParallelContext::IOProcessorNumberSub();
657#ifdef AMREX_LAZY
658 Lazy::QueueReduction( [=] () mutable {
659#endif
663
664 amrex::Print() << "ParticleContainer spread across MPI nodes - bytes: [Min: "
665 << mn
666 << ", Max: "
667 << mx
668 << ", Total: "
669 << cnt
670 << "]\n";
671#ifdef AMREX_LAZY
672 });
673#endif
674
675 return {mn, mx, cnt};
676}
677
678template <typename ParticleType, int NArrayReal, int NArrayInt,
679 template<class> class Allocator, class CellAssignor>
680void
682{
683 for (unsigned lev = 0; lev < m_particles.size(); lev++) {
684 auto& pmap = m_particles[lev];
685 for (auto& kv : pmap) {
686 auto& ptile = kv.second;
687 ptile.shrink_to_fit();
688 }
689 }
690}
691
698template <typename ParticleType, int NArrayReal, int NArrayInt,
699 template<class> class Allocator, class CellAssignor>
700void
702{
703 BL_PROFILE("ParticleContainer::Increment");
704
705 AMREX_ASSERT(OK());
706 if (m_particles.empty()) { return; }
707 AMREX_ASSERT(lev >= 0 && lev < int(m_particles.size()));
708 AMREX_ASSERT(numParticlesOutOfRange(*this, 0) == 0);
709
710 const auto& geom = Geom(lev);
711 const auto plo = geom.ProbLoArray();
712 const auto dxi = geom.InvCellSizeArray();
713 const auto domain = geom.Domain();
714 amrex::ParticleToMesh(*this, mf, lev,
715 [=] AMREX_GPU_DEVICE (const typename ParticleTileType::ConstParticleTileDataType& ptd, int ip,
716 amrex::Array4<amrex::Real> const& count)
718 const auto p = make_particle<ConstParticleType>{}(ptd, ip);
719 CellAssignor assignor;
720 IntVect iv = assignor(p, plo, dxi, domain);
721 amrex::Gpu::Atomic::AddNoRet(&count(iv), 1.0_rt);
722 }, false);
723}
724
725template <typename ParticleType, int NArrayReal, int NArrayInt,
726 template<class> class Allocator, class CellAssignor>
727Long
729{
730 BL_PROFILE("ParticleContainer::IncrementWithTotal(lev)");
731 Increment(mf, lev);
732 return TotalNumberOfParticles(true, local);
733}
734
735template <typename ParticleType, int NArrayReal, int NArrayInt,
736 template<class> class Allocator, class CellAssignor>
737void
739{
740 BL_PROFILE("ParticleContainer::RemoveParticlesAtLevel()");
741 if (level >= int(this->m_particles.size())) { return; }
742
743 if (!this->m_particles[level].empty())
744 {
745 ParticleLevel().swap(this->m_particles[level]);
746 }
747}
748
749template <typename ParticleType, int NArrayReal, int NArrayInt,
750 template<class> class Allocator, class CellAssignor>
751void
753{
754 BL_PROFILE("ParticleContainer::RemoveParticlesNotAtFinestLevel()");
755 AMREX_ASSERT(this->finestLevel()+1 == int(this->m_particles.size()));
756
757 Long cnt = 0;
758
759 for (unsigned lev = 0; lev < m_particles.size() - 1; ++lev) {
760 auto& pmap = m_particles[lev];
761 if (!pmap.empty()) {
762 for (auto& kv : pmap) {
763 const auto& pbx = kv.second;
764 cnt += pbx.numParticles();
765 }
766 ParticleLevel().swap(pmap);
767 }
768 }
769
770 //
771 // Print how many particles removed on each processor if any were removed.
772 //
773 if (this->m_verbose > 1 && cnt > 0) {
774 amrex::AllPrint() << "Processor " << ParallelContext::MyProcSub() << " removed " << cnt
775 << " particles not in finest level\n";
776 }
777}
778
780{
781
785
787 const GpuArray<Real, AMREX_SPACEDIM> & dxi, const Box& domain)
788 : m_assign_buffer_grid(assign_buffer_grid), m_plo(plo), m_dxi(dxi), m_domain(domain)
789 {}
790
791 template <typename SrcData>
793 int operator() (const SrcData& src, int src_i) const noexcept
794 {
795 auto iv = getParticleCell(src.m_aos[src_i], m_plo, m_dxi, m_domain);
796 return (m_assign_buffer_grid(iv)!=-1);
797 }
798};
799
801{
802 template <typename DstData, typename SrcData>
804 void operator() (DstData& dst, const SrcData& src,
805 int src_i, int dst_i) const noexcept
806 {
807 copyParticle(dst, src, src_i, dst_i);
808
809 (dst.m_aos[dst_i]).id() = LongParticleIds::VirtualParticleID;
810 (dst.m_aos[dst_i]).cpu() = 0;
811 }
812};
813
814
815template <typename ParticleType, int NArrayReal, int NArrayInt,
816 template<class> class Allocator, class CellAssignor>
817void
818ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
819::CreateVirtualParticles (int level, AoS& virts) const
820{
821 ParticleTileType ptile;
822 CreateVirtualParticles(level, ptile);
823 ptile.GetArrayOfStructs().swap(virts);
824}
825
826template <typename ParticleType, int NArrayReal, int NArrayInt,
827 template<class> class Allocator, class CellAssignor>
828void
831{
832 BL_PROFILE("ParticleContainer::CreateVirtualParticles()");
833 AMREX_ASSERT(level > 0);
834 AMREX_ASSERT(virts.empty());
835
836 if (level >= static_cast<int>(m_particles.size())) {
837 return;
838 }
839
840 std::string aggregation_type = AggregationType();
841 int aggregation_buffer = AggregationBuffer();
842
843 if (aggregation_type == "None")
844 {
845 auto virts_offset = virts.numParticles();
846 for(ParConstIterType pti(*this, level); pti.isValid(); ++pti)
847 {
848 const auto& src_tile = ParticlesAt(level, pti);
849
850 auto np = src_tile.numParticles();
851 virts.resize(virts_offset+np);
852 transformParticles(virts, src_tile, 0, virts_offset, np, TransformerVirt());
853 virts_offset += np;
854 }
855 }
856 if (aggregation_type == "Cell")
857 {
858 //Components would be based on
859 int nComp = AMREX_SPACEDIM + NStructReal + NArrayReal;
860 // NArrayReal, NStructInt, NArrayInt behavior as before
861 int nGhost = 0;
862 MultiFab mf(ParticleBoxArray(level), ParticleDistributionMap(level), nComp, nGhost);
863
864 nComp = 1 + NStructInt + NArrayInt;
865 iMultiFab imf(ParticleBoxArray(level), ParticleDistributionMap(level), nComp, nGhost);
866
867 const auto& geom = Geom(level);
868 const auto plo = geom.ProbLoArray();
869 const auto dxi = geom.InvCellSizeArray();
870 const auto domain = geom.Domain();
871
872 BoxList bl_buffer;
873 bl_buffer.complementIn(Geom(level).Domain(), ParticleBoxArray(level));
874 BoxArray buffer(std::move(bl_buffer));
875 buffer.grow(aggregation_buffer);
876
878 locator.build(buffer, geom);
879 AssignGrid<DenseBinIteratorFactory<Box>> assign_buffer_grid = locator.getGridAssignor();
880
881 amrex::ParticleToMesh(*this, mf, level,
882 [=] AMREX_GPU_DEVICE (const ParticleType& p,
883 amrex::Array4<amrex::Real> const& partData,
886 {
887 auto iv = getParticleCell(p, plo_loc, dxi_loc, domain);
888 if(assign_buffer_grid(iv)==-1)
889 {
890 //Ordering may make this not unique
891 for (int i = 0; i < NArrayReal; ++i)
892 {
893 amrex::Gpu::Atomic::AddNoRet(&partData(iv,AMREX_SPACEDIM+NStructReal+i), partData(iv,AMREX_SPACEDIM)!=0.0 ? static_cast<Real>(0) : static_cast<Real>(p.rdata(NStructReal+i)));
894 }
895 //Add the rdata(0)-weighted sum of position
896 for (int i = 0; i < AMREX_SPACEDIM; ++i)
897 {
898 amrex::Gpu::Atomic::AddNoRet(&partData(iv,i), static_cast<Real>((p.rdata(0)*p.pos(i))));
899 }
900 //Add the rdata(0)-weighted sum of other rdata fields
901 for (int i = 1; i < NStructReal; ++i)
902 {
903 amrex::Gpu::Atomic::AddNoRet(&partData(iv,AMREX_SPACEDIM+i), static_cast<Real>((p.rdata(0)*p.rdata(i))));
904 }
905 //Add the rdata(0) sum
906 for (int i = 0; i < 1; ++i)
907 {
908 amrex::Gpu::Atomic::AddNoRet(&partData(iv,AMREX_SPACEDIM+i), static_cast<Real>(p.rdata(0)));
909 }
910 }
911
912 }); //skipping extra false argument, doing mf.setVal(0) at beginning
913
914 amrex::ParticleToMesh(*this, imf, level,
915 [=] AMREX_GPU_DEVICE (const ParticleType& p,
916 amrex::Array4<int> const& partData,
919 {
920
921 auto iv = getParticleCell(p, plo_loc, dxi_loc, domain);
922 if(assign_buffer_grid(iv)==-1)
923 {
924 //if this cell has no particle id info, do a straight copy to store idata
925 if(partData(iv,0)==0)
926 {
927 //Add 1 to indicate at least 1 particle at cell iv
928 amrex::Gpu::Atomic::AddNoRet(&partData(iv,0), 1);
929 for (int i = 0; i < NStructInt; ++i)
930 {
931 amrex::Gpu::Atomic::AddNoRet(&partData(iv,1+i), p.idata(i));
932 }
933 for (int i = 0; i < NArrayInt; ++i)
934 {
935 amrex::Gpu::Atomic::AddNoRet(&partData(iv,1+NStructInt+i), p.idata(NStructInt+i));
936 }
937 }
938 }
939 });
940
941 //There may be a better way to ensure virts is the right length
942 virts.resize(imf.sum(0));
943
944 int last_offset = 0;
945 for (MFIter mfi(mf); mfi.isValid(); ++mfi)
946 {
947 const auto bx = mfi.tilebox();
948 const auto partData = mf.array(mfi);
949 const auto imf_arr = imf.array(mfi);
950
951 Gpu::DeviceVector<int> offsets(bx.numPts());
952 auto *offsets_ptr = offsets.dataPtr();
953 int next_offset = Scan::ExclusiveSum((int) bx.numPts(),(imf_arr.ptr(bx.smallEnd(),0)),(offsets.dataPtr()),Scan::retSum);
954 auto dst = virts.getParticleTileData();
955 ParallelFor(bx, [=] AMREX_GPU_DEVICE(int i, int j, int k)
956 {
957 if(imf_arr(i,j,k,0)!=0)
958 {
960 p.cpu() = 0;
962
963 //Set rdata(0) first so we can normalize the weighted fields
964 p.rdata(0) = static_cast<ParticleReal>(partData(i,j,k,AMREX_SPACEDIM));
965 //Set pos with the normalized weighted field
966 for (int n = 0; n < AMREX_SPACEDIM; ++n)
967 {
968 p.pos(n) = static_cast<ParticleReal>(partData(i,j,k,n) / p.rdata(0));
969 }
970 //Set rdata(n>0) with the normalized weighted field for NStructReal
971 for (int n = 1; n < NStructReal; ++n)
972 {
973 p.rdata(n) = static_cast<ParticleReal>(partData(i,j,k,AMREX_SPACEDIM+n) / p.rdata(0));
974 }
975 //Set rdata(n>0) with the normalized weighted field for NArrayReal
976 for (int n = 0; n < NArrayReal; ++n)
977 {
978 p.rdata(NStructReal+n) = static_cast<ParticleReal>(partData(i,j,k,AMREX_SPACEDIM+NStructReal+n));
979 }
980 //Set idata with the "first" particles idata field for NStructInt
981 for (int n = 0; n < NStructInt; ++n)
982 {
983 p.idata(n) = imf_arr(i,j,k,1+n);
984 }
985 //Set idata with the "first" particles idata field for NArrayInt
986 for (int n = 0; n < NArrayInt; ++n)
987 {
988 p.idata(NStructInt+n) = imf_arr(i,j,k,1+NStructInt+n);
989 }
990 //Set the new particle in dst tile
991 dst.setSuperParticle(p, last_offset+offsets_ptr[((i-imf_arr.begin.x)+(j-imf_arr.begin.y)*imf_arr.jstride+(k-imf_arr.begin.z)*imf_arr.kstride)]);
992 }
993
994 });
995 last_offset+=next_offset;
997 }
998
999 // last_offset should equal virts.numParticles()
1000 auto virts_offset = last_offset;
1001 for(ParConstIterType pti(*this, level); pti.isValid(); ++pti)
1002 {
1003 const auto& src_tile = ParticlesAt(level, pti);
1004
1005 auto np = src_tile.numParticles();
1006 virts.resize(virts_offset+np);
1007 virts_offset += filterAndTransformParticles(virts, src_tile, FilterVirt(assign_buffer_grid,plo,dxi,domain), TransformerVirt(),0,virts_offset);
1009 }
1010 virts.resize(virts_offset);
1012 }
1013}
1014
1016{
1017
1020
1026 : m_lev_min(level), m_lev_max(level+1), m_nGrow(nGrow), m_gid(gid), m_assign_grid(assign_grid)
1027 {}
1028
1029 template <typename SrcData>
1031 int operator() (const SrcData& src, int src_i) const noexcept
1032 {
1033 const auto tup_min = (m_assign_grid)(src.m_aos[src_i], m_lev_min, m_lev_max, m_nGrow, DefaultAssignor{});
1034 const auto tup_max = (m_assign_grid)(src.m_aos[src_i], m_lev_max, m_lev_max, m_nGrow, DefaultAssignor{});
1035 const auto p_boxes = amrex::get<0>(tup_min);
1036 const auto p_boxes_max = amrex::get<0>(tup_max);
1037 const auto p_levs_max = amrex::get<1>(tup_max);
1038 return p_boxes_max >=0 && p_boxes == m_gid && p_levs_max == m_lev_max;
1039 }
1040};
1041
1043{
1044
1045 template <typename DstData, typename SrcData>
1047 void operator() (DstData& dst, const SrcData& src,
1048 int src_i, int dst_i) const noexcept
1049 {
1050 copyParticle(dst, src, src_i, dst_i);
1051
1052 (dst.m_aos[dst_i]).id() = LongParticleIds::GhostParticleID;
1053 (dst.m_aos[dst_i]).cpu() = 0;
1054 }
1055};
1056
1057template <typename ParticleType, int NArrayReal, int NArrayInt,
1058 template<class> class Allocator, class CellAssignor>
1059void
1060ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
1061::CreateGhostParticles (int level, int nGrow, AoS& ghosts) const
1062{
1063 ParticleTileType ptile;
1064 CreateGhostParticles(level, nGrow, ptile);
1065 ptile.GetArrayOfStructs().swap(ghosts);
1066}
1067
1068template <typename ParticleType, int NArrayReal, int NArrayInt,
1069 template<class> class Allocator, class CellAssignor>
1070void
1072::CreateGhostParticles (int level, int nGrow, ParticleTileType& ghosts) const
1073{
1074 BL_PROFILE("ParticleContainer::CreateGhostParticles()");
1075 AMREX_ASSERT(ghosts.empty());
1076 AMREX_ASSERT(level < finestLevel());
1077
1078 if (level >= static_cast<int>(m_particles.size())) {
1079 return;
1080 }
1081
1082 if (! m_particle_locator.isValid(GetParGDB())) {
1083 m_particle_locator.build(GetParGDB());
1084 }
1085
1086 m_particle_locator.setGeometry(GetParGDB());
1087 AmrAssignGrid<DenseBinIteratorFactory<Box>> assign_grid = m_particle_locator.getGridAssignor();
1088 auto ghost_offset = ghosts.numParticles();
1089 for(ParConstIterType pti(*this, level); pti.isValid(); ++pti)
1090 {
1091 const auto& src_tile = ParticlesAt(level, pti);
1092 int gid = pti.index();
1093
1094 auto np = src_tile.numParticles();
1095 ghosts.resize(ghost_offset+np);
1096 ghost_offset += filterAndTransformParticles(ghosts, src_tile, AssignGridFilter(assign_grid,gid,level,nGrow), TransformerGhost(),0,ghost_offset);
1097 }
1098 ghosts.resize(ghost_offset);
1099 Gpu::streamSynchronize();
1100}
1101
1102template <typename ParticleType, int NArrayReal, int NArrayInt,
1103 template<class> class Allocator, class CellAssignor>
1104void
1107{
1108 BL_PROFILE("ParticleContainer::clearParticles()");
1109
1110 for (int lev = 0; lev < static_cast<int>(m_particles.size()); ++lev)
1111 {
1112 for (auto& kv : m_particles[lev]) { kv.second.resize(0); }
1113 particle_detail::clearEmptyEntries(m_particles[lev]);
1114 }
1115}
1116
1117template <typename ParticleType, int NArrayReal, int NArrayInt,
1118 template<class> class Allocator, class CellAssignor>
1119template <class PCType, std::enable_if_t<IsParticleContainer<PCType>::value, int> foo>
1120void
1127
1128template <typename ParticleType, int NArrayReal, int NArrayInt,
1129 template<class> class Allocator, class CellAssignor>
1130template <class PCType, std::enable_if_t<IsParticleContainer<PCType>::value, int> foo>
1131void
1133addParticles (const PCType& other, bool local)
1134{
1136 addParticles(other, [] AMREX_GPU_HOST_DEVICE (const PData& /*data*/, int /*i*/) { return 1; }, local);
1137}
1138
1139template <typename ParticleType, int NArrayReal, int NArrayInt,
1140 template<class> class Allocator, class CellAssignor>
1141template <class F, class PCType,
1142 std::enable_if_t<IsParticleContainer<PCType>::value, int> foo,
1143 std::enable_if_t<! std::is_integral_v<F>, int> bar>
1144void
1146copyParticles (const PCType& other, F&& f, bool local)
1147{
1148 BL_PROFILE("ParticleContainer::copyParticles");
1149 clearParticles();
1150 addParticles(other, std::forward<F>(f), local);
1151}
1152
1153template <typename ParticleType, int NArrayReal, int NArrayInt,
1154 template<class> class Allocator, class CellAssignor>
1155template <class F, class PCType,
1156 std::enable_if_t<IsParticleContainer<PCType>::value, int> foo,
1157 std::enable_if_t<! std::is_integral_v<F>, int> bar>
1158void
1160addParticles (const PCType& other, F const& f, bool local)
1161{
1162 BL_PROFILE("ParticleContainer::addParticles");
1163
1164 for (int lev = 0; lev < other.numLevels(); ++lev)
1165 {
1166 const auto& plevel_other = other.GetParticles(lev);
1167 for(MFIter mfi = other.MakeMFIter(lev); mfi.isValid(); ++mfi)
1168 {
1169 auto index = std::make_pair(mfi.index(), mfi.LocalTileIndex());
1170 if(plevel_other.find(index) == plevel_other.end()) { continue; }
1171
1172 auto& ptile = DefineAndReturnParticleTile(lev, mfi.index(), mfi.LocalTileIndex());
1173 const auto& ptile_other = plevel_other.at(index);
1174 auto np = ptile_other.numParticles();
1175 if (np == 0) { continue; }
1176
1177 auto dst_index = ptile.numParticles();
1178 ptile.resize(dst_index + np);
1179
1180 auto count = amrex::filterParticles(ptile, ptile_other, f, 0, dst_index, np);
1181
1182 ptile.resize(dst_index + count);
1183 }
1184 }
1185
1186 if (! local) { Redistribute(); }
1187}
1188
1189//
1190// This redistributes valid particles and discards invalid ones.
1191//
1192template <typename ParticleType, int NArrayReal, int NArrayInt,
1193 template<class> class Allocator, class CellAssignor>
1194void
1196::Redistribute (int lev_min, int lev_max, int nGrow, int local, bool remove_negative)
1197{
1198 BL_PROFILE_SYNC_START_TIMED("SyncBeforeComms: Redist");
1199
1200#ifdef AMREX_USE_GPU
1201 if ( Gpu::inLaunchRegion() )
1202 {
1203 RedistributeGPU(lev_min, lev_max, nGrow, local, remove_negative);
1204 }
1205 else
1206 {
1207 RedistributeCPU(lev_min, lev_max, nGrow, local, remove_negative);
1208 }
1209#else
1210 RedistributeCPU(lev_min, lev_max, nGrow, local, remove_negative);
1211#endif
1212
1214}
1215
1216template <typename ParticleType, int NArrayReal, int NArrayInt,
1217 template<class> class Allocator, class CellAssignor>
1218template <class index_type>
1219void
1221::ReorderParticles (int lev, const MFIter& mfi, const index_type* permutations)
1222{
1223 auto& ptile = ParticlesAt(lev, mfi);
1224 const size_t np = ptile.numParticles();
1225 const size_t np_total = np + ptile.numNeighborParticles();
1226
1227 if (memEfficientSort) {
1228#if defined(AMREX_USE_CUDA) && defined(_WIN32)
1229 if (!ParticleType::is_soa_particle) {
1230#else
1231 if constexpr (!ParticleType::is_soa_particle) {
1232#endif
1233 static_assert(sizeof(ParticleType)%4 == 0 && sizeof(uint32_t) == 4);
1234 using tmp_t = std::conditional_t<sizeof(ParticleType)%8 == 0,
1235 uint64_t, uint32_t>;
1236 constexpr std::size_t nchunks = sizeof(ParticleType) / sizeof(tmp_t);
1238 auto* ptmp = tmp.data();
1239 auto* paos = (tmp_t*)(ptile.getParticleTileData().m_aos);
1240 for (std::size_t ichunk = 0; ichunk < nchunks; ++ichunk) {
1241 // Do not need to reorder neighbor particles
1243 {
1244 ptmp[i] = paos[permutations[i]*nchunks+ichunk];
1245 });
1247 {
1248 paos[i*nchunks+ichunk] = ptmp[i];
1249 });
1250 }
1251 Gpu::streamSynchronize();
1252 } else {
1253 typename SoA::IdCPU tmp_idcpu(np_total);
1254
1255 auto src = ptile.GetStructOfArrays().GetIdCPUData().data();
1256 uint64_t* dst = tmp_idcpu.data();
1257 AMREX_HOST_DEVICE_FOR_1D( np_total, i,
1258 {
1259 dst[i] = i < np ? src[permutations[i]] : src[i];
1260 });
1261
1262 Gpu::streamSynchronize();
1263
1264 ptile.GetStructOfArrays().GetIdCPUData().swap(tmp_idcpu);
1265 }
1266
1267 { // Create a scope for the temporary vector below
1268 RealVector tmp_real(np_total);
1269 for (int comp = 0; comp < NArrayReal + m_num_runtime_real; ++comp) {
1270 auto src = ptile.GetStructOfArrays().GetRealData(comp).data();
1271 ParticleReal* dst = tmp_real.data();
1272 AMREX_HOST_DEVICE_FOR_1D( np_total, i,
1273 {
1274 dst[i] = i < np ? src[permutations[i]] : src[i];
1275 });
1276
1277 Gpu::streamSynchronize();
1278
1279 ptile.GetStructOfArrays().GetRealData(comp).swap(tmp_real);
1280 }
1281 }
1282
1283 IntVector tmp_int(np_total);
1284 for (int comp = 0; comp < NArrayInt + m_num_runtime_int; ++comp) {
1285 auto src = ptile.GetStructOfArrays().GetIntData(comp).data();
1286 int* dst = tmp_int.data();
1287 AMREX_HOST_DEVICE_FOR_1D( np_total , i,
1288 {
1289 dst[i] = i < np ? src[permutations[i]] : src[i];
1290 });
1291
1292 Gpu::streamSynchronize();
1293
1294 ptile.GetStructOfArrays().GetIntData(comp).swap(tmp_int);
1295 }
1296 } else {
1297 ParticleTileType ptile_tmp;
1298 ptile_tmp.define(m_num_runtime_real, m_num_runtime_int, &m_soa_rdata_names, &m_soa_idata_names);
1299 ptile_tmp.resize(np_total);
1300 // copy re-ordered particles
1301 gatherParticles(ptile_tmp, ptile, np, permutations);
1302 // copy neighbor particles
1303 amrex::copyParticles(ptile_tmp, ptile, np, np, np_total-np);
1304 ptile.swap(ptile_tmp);
1305 }
1306}
1307
1308template <typename ParticleType, int NArrayReal, int NArrayInt,
1309 template<class> class Allocator, class CellAssignor>
1310void
1315
1316template <typename ParticleType, int NArrayReal, int NArrayInt,
1317 template<class> class Allocator, class CellAssignor>
1318void
1321{
1322 BL_PROFILE("ParticleContainer::SortParticlesByBin()");
1323
1324 if (bin_size == IntVect::TheZeroVector()) { return; }
1325
1326 for (int lev = 0; lev < numLevels(); ++lev)
1327 {
1328 const Geometry& geom = Geom(lev);
1329 const auto dxi = geom.InvCellSizeArray();
1330 const auto plo = geom.ProbLoArray();
1331 const auto domain = geom.Domain();
1332
1333 for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi)
1334 {
1335 auto& ptile = ParticlesAt(lev, mfi);
1336 const size_t np = ptile.numParticles();
1337
1338 const Box& box = mfi.validbox();
1339
1340 int ntiles = numTilesInBox(box, true, bin_size);
1341
1342 m_bins.build(np, ptile.getParticleTileData(), ntiles,
1343 GetParticleBin{plo, dxi, domain, bin_size, box});
1344 ReorderParticles(lev, mfi, m_bins.permutationPtr());
1346 }
1347}
1348
1349template <typename ParticleType, int NArrayReal, int NArrayInt,
1350 template<class> class Allocator, class CellAssignor>
1351void
1352ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
1353::SortParticlesForDeposition (IntVect idx_type)
1354{
1355 BL_PROFILE("ParticleContainer::SortParticlesForDeposition()");
1356
1357 for (int lev = 0; lev < numLevels(); ++lev)
1358 {
1359 const Geometry& geom = Geom(lev);
1360
1361 for(MFIter mfi = MakeMFIter(lev); mfi.isValid(); ++mfi)
1362 {
1363 const auto& ptile = ParticlesAt(lev, mfi);
1364 const size_t np = ptile.numParticles();
1365
1366 const Box& box = mfi.validbox();
1367
1368 using index_type = typename decltype(m_bins)::index_type;
1370 PermutationForDeposition<index_type>(perm, np, ptile, box, geom, idx_type);
1371 ReorderParticles(lev, mfi, perm.dataPtr());
1372 }
1373 }
1374}
1375
1376//
1377// The GPU implementation of Redistribute
1378//
1379template <typename ParticleType, int NArrayReal, int NArrayInt,
1380 template<class> class Allocator, class CellAssignor>
1381void
1383::RedistributeGPU (int lev_min, int lev_max, int nGrow, int local, bool remove_negative)
1384{
1385#ifdef AMREX_USE_GPU
1386
1387 if (local) { AMREX_ASSERT(numParticlesOutOfRange(*this, lev_min, lev_max, local) == 0); }
1388
1389 // sanity check
1390 AMREX_ALWAYS_ASSERT(do_tiling == false);
1391
1392 BL_PROFILE("ParticleContainer::RedistributeGPU()");
1393 BL_PROFILE_VAR_NS("Redistribute_partition", blp_partition);
1394
1395 int theEffectiveFinestLevel = m_gdb->finestLevel();
1396 while (!m_gdb->LevelDefined(theEffectiveFinestLevel)) { theEffectiveFinestLevel--; }
1397
1398 if (int(m_particles.size()) < theEffectiveFinestLevel+1) {
1399 if (Verbose()) {
1400 amrex::Print() << "ParticleContainer::Redistribute() resizing containers from "
1401 << m_particles.size() << " to "
1402 << theEffectiveFinestLevel + 1 << '\n';
1403 }
1404 m_particles.resize(theEffectiveFinestLevel+1);
1405 m_dummy_mf.resize(theEffectiveFinestLevel+1);
1406 }
1407
1408 for (int lev = 0; lev < theEffectiveFinestLevel+1; ++lev) { RedefineDummyMF(lev); }
1409
1410 int finest_lev_particles;
1411 if (lev_max == -1) {
1412 lev_max = theEffectiveFinestLevel;
1413 finest_lev_particles = m_particles.size() - 1;
1414 } else {
1415 finest_lev_particles = lev_max;
1416 }
1417 AMREX_ASSERT(lev_max <= finestLevel());
1418
1419 this->defineBufferMap();
1420
1421 if (! m_particle_locator.isValid(GetParGDB())) { m_particle_locator.build(GetParGDB()); }
1422 m_particle_locator.setGeometry(GetParGDB());
1423 auto assign_grid = m_particle_locator.getGridAssignor();
1424
1425 BL_PROFILE_VAR_START(blp_partition);
1427 int num_levels = finest_lev_particles + 1;
1428 op.setNumLevels(num_levels);
1429 Vector<std::map<int, int> > new_sizes(num_levels);
1430 const auto plo = Geom(0).ProbLoArray();
1431 const auto phi = Geom(0).ProbHiArray();
1432 const auto rlo = Geom(0).ProbLoArrayInParticleReal();
1433 const auto rhi = Geom(0).ProbHiArrayInParticleReal();
1434 const auto is_per = Geom(0).isPeriodicArray();
1435 for (int lev = lev_min; lev <= finest_lev_particles; ++lev)
1436 {
1437 auto& plev = m_particles[lev];
1438 for (auto& kv : plev)
1439 {
1440 int gid = kv.first.first;
1441 int tid = kv.first.second;
1442 auto index = std::make_pair(gid, tid);
1443
1444 auto& src_tile = plev[index];
1445 const size_t np = src_tile.numParticles();
1446
1447 int num_stay = partitionParticlesByDest(src_tile, assign_grid,
1448 std::forward<CellAssignor>(CellAssignor{}),
1449 BufferMap(),
1450 plo, phi, rlo, rhi, is_per, lev, gid, tid,
1451 lev_min, lev_max, nGrow, remove_negative);
1452
1453 int num_move = np - num_stay;
1454 new_sizes[lev][gid] = num_stay;
1455 op.resize(gid, lev, num_move);
1456
1457 auto p_boxes = op.m_boxes[lev][gid].dataPtr();
1458 auto p_levs = op.m_levels[lev][gid].dataPtr();
1459 auto p_src_indices = op.m_src_indices[lev][gid].dataPtr();
1460 auto p_periodic_shift = op.m_periodic_shift[lev][gid].dataPtr();
1461 auto ptd = src_tile.getParticleTileData();
1462
1463 AMREX_FOR_1D ( num_move, i,
1464 {
1465 const auto p = make_particle<ParticleType>{}(ptd,i + num_stay);
1466
1467 if (p.id() < 0)
1469 p_boxes[i] = -1;
1470 p_levs[i] = -1;
1471 }
1472 else
1473 {
1474 const auto tup = assign_grid(p, lev_min, lev_max, nGrow,
1475 std::forward<CellAssignor>(CellAssignor{}));
1476 p_boxes[i] = amrex::get<0>(tup);
1477 p_levs[i] = amrex::get<1>(tup);
1478 }
1479 p_periodic_shift[i] = IntVect(AMREX_D_DECL(0,0,0));
1480 p_src_indices[i] = i+num_stay;
1481 });
1482 }
1483 }
1484 BL_PROFILE_VAR_STOP(blp_partition);
1485
1486 ParticleCopyPlan plan;
1487
1488 plan.build(*this, op, h_redistribute_int_comp,
1489 h_redistribute_real_comp, local);
1490
1491 // by default, this uses The_Arena();
1494
1495 if (use_comms_arena) {
1496 snd_buffer.setArena(The_Comms_Arena());
1497 rcv_buffer.setArena(The_Comms_Arena());
1498 }
1499
1500 packBuffer(*this, op, plan, snd_buffer);
1501
1502 // clear particles from container
1503 for (int lev = lev_min; lev <= lev_max; ++lev)
1504 {
1505 auto& plev = m_particles[lev];
1506 for (auto& kv : plev)
1507 {
1508 int gid = kv.first.first;
1509 int tid = kv.first.second;
1510 auto index = std::make_pair(gid, tid);
1511 auto& tile = plev[index];
1512 tile.resize(new_sizes[lev][gid]);
1513 }
1514 }
1516 for (int lev = lev_min; lev <= lev_max; lev++)
1517 {
1518 particle_detail::clearEmptyEntries(m_particles[lev]);
1519 }
1520
1521 if (int(m_particles.size()) > theEffectiveFinestLevel+1) {
1522 if (m_verbose > 0) {
1523 amrex::Print() << "ParticleContainer::Redistribute() resizing m_particles from "
1524 << m_particles.size() << " to " << theEffectiveFinestLevel+1 << '\n';
1525 }
1526 AMREX_ASSERT(int(m_particles.size()) >= 2);
1527
1528 m_particles.resize(theEffectiveFinestLevel + 1);
1529 m_dummy_mf.resize(theEffectiveFinestLevel + 1);
1530 }
1531
1532 if (ParallelDescriptor::UseGpuAwareMpi())
1533 {
1534 plan.buildMPIFinish(BufferMap());
1535 communicateParticlesStart(*this, plan, snd_buffer, rcv_buffer);
1536 unpackBuffer(*this, plan, snd_buffer, RedistributeUnpackPolicy());
1538 unpackRemotes(*this, plan, rcv_buffer, RedistributeUnpackPolicy());
1539 }
1540 else
1541 {
1542 Gpu::Device::streamSynchronize();
1543 Gpu::PinnedVector<char> pinned_snd_buffer;
1544 Gpu::PinnedVector<char> pinned_rcv_buffer;
1545
1546 if (snd_buffer.arena()->isPinned()) {
1547 plan.buildMPIFinish(BufferMap());
1548 Gpu::Device::streamSynchronize();
1549 communicateParticlesStart(*this, plan, snd_buffer, pinned_rcv_buffer);
1550 } else {
1551 pinned_snd_buffer.resize(snd_buffer.size());
1552 Gpu::dtoh_memcpy_async(pinned_snd_buffer.dataPtr(), snd_buffer.dataPtr(), snd_buffer.size());
1553 plan.buildMPIFinish(BufferMap());
1554 Gpu::Device::streamSynchronize();
1555 communicateParticlesStart(*this, plan, pinned_snd_buffer, pinned_rcv_buffer);
1556 }
1557
1558 rcv_buffer.resize(pinned_rcv_buffer.size());
1559 unpackBuffer(*this, plan, snd_buffer, RedistributeUnpackPolicy());
1561 Gpu::htod_memcpy_async(rcv_buffer.dataPtr(), pinned_rcv_buffer.dataPtr(), pinned_rcv_buffer.size());
1562 unpackRemotes(*this, plan, rcv_buffer, RedistributeUnpackPolicy());
1563 }
1564
1565 Gpu::Device::streamSynchronize();
1566 AMREX_ASSERT(numParticlesOutOfRange(*this, lev_min, lev_max, nGrow) == 0);
1567#else
1568 amrex::ignore_unused(lev_min,lev_max,nGrow,local,remove_negative);
1569#endif
1570}
1571
1572//
1573// The CPU implementation of Redistribute
1574//
1575template <typename ParticleType, int NArrayReal, int NArrayInt,
1576 template<class> class Allocator, class CellAssignor>
1577void
1578ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
1579::RedistributeCPU (int lev_min, int lev_max, int nGrow, int local, bool remove_negative)
1580{
1581 BL_PROFILE("ParticleContainer::RedistributeCPU()");
1582
1583 const int MyProc = ParallelContext::MyProcSub();
1584 auto strttime = amrex::second();
1585
1586 if (local > 0) { BuildRedistributeMask(0, local); }
1587
1588 // On startup there are cases where Redistribute() could be called
1589 // with a given finestLevel() where that AmrLevel has yet to be defined.
1590 int theEffectiveFinestLevel = m_gdb->finestLevel();
1591
1592 while (!m_gdb->LevelDefined(theEffectiveFinestLevel)) {
1593 theEffectiveFinestLevel--;
1594 }
1595
1596 if (int(m_particles.size()) < theEffectiveFinestLevel+1) {
1597 if (Verbose()) {
1598 amrex::Print() << "ParticleContainer::Redistribute() resizing containers from "
1599 << m_particles.size() << " to "
1600 << theEffectiveFinestLevel + 1 << '\n';
1601 }
1602 m_particles.resize(theEffectiveFinestLevel+1);
1603 m_dummy_mf.resize(theEffectiveFinestLevel+1);
1604 }
1605
1606 // It is important to do this even if we don't have more levels because we may have changed the
1607 // grids at this level in a regrid.
1608 for (int lev = 0; lev < theEffectiveFinestLevel+1; ++lev) {
1609 RedefineDummyMF(lev);
1610 }
1611
1612 int finest_lev_particles;
1613 if (lev_max == -1) {
1614 lev_max = theEffectiveFinestLevel;
1615 finest_lev_particles = m_particles.size() - 1;
1616 } else {
1617 finest_lev_particles = lev_max;
1618 }
1619 AMREX_ASSERT(lev_max <= finestLevel());
1620
1621 // This will hold the valid particles that go to another process
1622 std::map<int, Vector<char> > not_ours;
1623
1624 int num_threads = OpenMP::get_max_threads();
1625
1626 // these are temporary buffers for each thread
1627 std::map<int, Vector<Vector<char> > > tmp_remote;
1630 tmp_local.resize(theEffectiveFinestLevel+1);
1631 soa_local.resize(theEffectiveFinestLevel+1);
1632
1633 // we resize these buffers outside the parallel region
1634 for (int lev = lev_min; lev <= lev_max; lev++) {
1635 for (MFIter mfi(*m_dummy_mf[lev], this->do_tiling ? this->tile_size : IntVect::TheZeroVector());
1636 mfi.isValid(); ++mfi) {
1637 auto index = std::make_pair(mfi.index(), mfi.LocalTileIndex());
1638 tmp_local[lev][index].resize(num_threads);
1639 soa_local[lev][index].resize(num_threads);
1640 for (int t = 0; t < num_threads; ++t) {
1641 soa_local[lev][index][t].define(m_num_runtime_real, m_num_runtime_int, &m_soa_rdata_names, &m_soa_idata_names);
1642 }
1643 }
1644 }
1645 if (local) {
1646 for (int i = 0; i < neighbor_procs.size(); ++i) {
1647 tmp_remote[neighbor_procs[i]].resize(num_threads);
1648 }
1649 } else {
1650 for (int i = 0; i < ParallelContext::NProcsSub(); ++i) {
1651 tmp_remote[i].resize(num_threads);
1652 }
1653 }
1654
1655 // first pass: for each tile in parallel, in each thread copies the particles that
1656 // need to be moved into it's own, temporary buffer.
1657 for (int lev = lev_min; lev <= finest_lev_particles; lev++) {
1658 auto& pmap = m_particles[lev];
1659
1660 Vector<std::pair<int, int> > grid_tile_ids;
1661 Vector<ParticleTileType*> ptile_ptrs;
1662 for (auto& kv : pmap)
1663 {
1664 grid_tile_ids.push_back(kv.first);
1665 ptile_ptrs.push_back(&(kv.second));
1666 }
1667
1668#ifdef AMREX_USE_OMP
1669#pragma omp parallel for
1670#endif
1671 for (int pmap_it = 0; pmap_it < static_cast<int>(ptile_ptrs.size()); ++pmap_it)
1672 {
1673 int thread_num = OpenMP::get_thread_num();
1674 int grid = grid_tile_ids[pmap_it].first;
1675 int tile = grid_tile_ids[pmap_it].second;
1676 auto& soa = ptile_ptrs[pmap_it]->GetStructOfArrays();
1677 auto& aos = ptile_ptrs[pmap_it]->GetArrayOfStructs();
1678
1679 // AMREX_ASSERT_WITH_MESSAGE((NumRealComps() == 0 && NumIntComps() == 0)
1680 // || aos.size() == soa.size(),
1681 // "The AoS and SoA data on this tile are different sizes - "
1682 // "perhaps particles have not been initialized correctly?");
1683 unsigned npart = ptile_ptrs[pmap_it]->numParticles();
1684 ParticleLocData pld;
1685
1686 if constexpr (!ParticleType::is_soa_particle){
1687
1688 if (npart != 0) {
1689 Long last = npart - 1;
1690 Long pindex = 0;
1691 while (pindex <= last) {
1692 ParticleType& p = aos[pindex];
1693
1694 if ((remove_negative == false) && (p.id() < 0)) {
1695 ++pindex;
1696 continue;
1697 }
1698
1699 if (p.id() < 0)
1700 {
1701 aos[pindex] = aos[last];
1702 for (int comp = 0; comp < NumRealComps(); comp++) {
1703 soa.GetRealData(comp)[pindex] = soa.GetRealData(comp)[last];
1704 }
1705 for (int comp = 0; comp < NumIntComps(); comp++) {
1706 soa.GetIntData(comp)[pindex] = soa.GetIntData(comp)[last];
1707 }
1708 correctCellVectors(last, pindex, grid, aos[pindex]);
1709 --last;
1710 continue;
1711 }
1712
1713 locateParticle(p, pld, lev_min, lev_max, nGrow, local ? grid : -1);
1714
1715 particlePostLocate(p, pld, lev);
1716
1717 if (p.id() < 0)
1718 {
1719 aos[pindex] = aos[last];
1720 for (int comp = 0; comp < NumRealComps(); comp++) {
1721 soa.GetRealData(comp)[pindex] = soa.GetRealData(comp)[last];
1722 }
1723 for (int comp = 0; comp < NumIntComps(); comp++) {
1724 soa.GetIntData(comp)[pindex] = soa.GetIntData(comp)[last];
1725 }
1726 correctCellVectors(last, pindex, grid, aos[pindex]);
1727 --last;
1728 continue;
1729 }
1730
1731 const int who = ParallelContext::global_to_local_rank(ParticleDistributionMap(pld.m_lev)[pld.m_grid]);
1732 if (who == MyProc) {
1733 if (pld.m_lev != lev || pld.m_grid != grid || pld.m_tile != tile) {
1734 // We own it but must shift it to another place.
1735 auto index = std::make_pair(pld.m_grid, pld.m_tile);
1736 AMREX_ASSERT(tmp_local[pld.m_lev][index].size() == num_threads);
1737 tmp_local[pld.m_lev][index][thread_num].push_back(p);
1738 for (int comp = 0; comp < NumRealComps(); ++comp) {
1739 RealVector& arr = soa_local[pld.m_lev][index][thread_num].GetRealData(comp);
1740 arr.push_back(soa.GetRealData(comp)[pindex]);
1741 }
1742 for (int comp = 0; comp < NumIntComps(); ++comp) {
1743 IntVector& arr = soa_local[pld.m_lev][index][thread_num].GetIntData(comp);
1744 arr.push_back(soa.GetIntData(comp)[pindex]);
1745 }
1746
1747 p.id() = -p.id(); // Invalidate the particle
1748 }
1749 }
1750 else {
1751 auto& particles_to_send = tmp_remote[who][thread_num];
1752 auto old_size = particles_to_send.size();
1753 auto new_size = old_size + superparticle_size;
1754 particles_to_send.resize(new_size);
1755 std::memcpy(&particles_to_send[old_size], &p, particle_size);
1756 char* dst = &particles_to_send[old_size] + particle_size;
1757 int array_comp_start = 0;
1758 if constexpr (!ParticleType::is_soa_particle) {
1759 array_comp_start = AMREX_SPACEDIM + NStructReal;
1760 }
1761 for (int comp = 0; comp < NumRealComps(); comp++) {
1762 if (h_redistribute_real_comp[array_comp_start + comp]) {
1763 std::memcpy(dst, &soa.GetRealData(comp)[pindex], sizeof(ParticleReal));
1764 dst += sizeof(ParticleReal);
1765 }
1766 }
1767 array_comp_start = 2 + NStructInt;
1768 for (int comp = 0; comp < NumIntComps(); comp++) {
1769 if (h_redistribute_int_comp[array_comp_start + comp]) {
1770 std::memcpy(dst, &soa.GetIntData(comp)[pindex], sizeof(int));
1771 dst += sizeof(int);
1772 }
1773 }
1774
1775 p.id() = -p.id(); // Invalidate the particle
1776 }
1777
1778 if (p.id() < 0)
1779 {
1780 aos[pindex] = aos[last];
1781 for (int comp = 0; comp < NumRealComps(); comp++) {
1782 soa.GetRealData(comp)[pindex] = soa.GetRealData(comp)[last];
1783 }
1784 for (int comp = 0; comp < NumIntComps(); comp++) {
1785 soa.GetIntData(comp)[pindex] = soa.GetIntData(comp)[last];
1786 }
1787 correctCellVectors(last, pindex, grid, aos[pindex]);
1788 --last;
1789 continue;
1790 }
1791
1792 ++pindex;
1793 }
1794
1795 aos().erase(aos().begin() + last + 1, aos().begin() + npart);
1796 for (int comp = 0; comp < NumRealComps(); comp++) {
1797 RealVector& rdata = soa.GetRealData(comp);
1798 rdata.erase(rdata.begin() + last + 1, rdata.begin() + npart);
1799 }
1800 for (int comp = 0; comp < NumIntComps(); comp++) {
1801 IntVector& idata = soa.GetIntData(comp);
1802 idata.erase(idata.begin() + last + 1, idata.begin() + npart);
1803 }
1804 }
1805
1806 } else { // soa particle
1807
1808 auto particle_tile = ptile_ptrs[pmap_it];
1809 if (npart != 0) {
1810 Long last = npart - 1;
1811 Long pindex = 0;
1812 auto ptd = particle_tile->getParticleTileData();
1813 while (pindex <= last) {
1814 ParticleType p(ptd,pindex);
1815
1816 if ((remove_negative == false) && (p.id() < 0)) {
1817 ++pindex;
1818 continue;
1819 }
1820
1821 if (p.id() < 0){
1822 soa.GetIdCPUData()[pindex] = soa.GetIdCPUData()[last];
1823 for (int comp = 0; comp < NumRealComps(); comp++) {
1824 soa.GetRealData(comp)[pindex] = soa.GetRealData(comp)[last];
1825 }
1826 for (int comp = 0; comp < NumIntComps(); comp++) {
1827 soa.GetIntData(comp)[pindex] = soa.GetIntData(comp)[last];
1828 }
1829 correctCellVectors(last, pindex, grid, ptd[pindex]);
1830 --last;
1831 continue;
1832 }
1833
1834 locateParticle(p, pld, lev_min, lev_max, nGrow, local ? grid : -1);
1835
1836 particlePostLocate(p, pld, lev);
1837
1838 if (p.id() < 0) {
1839 soa.GetIdCPUData()[pindex] = soa.GetIdCPUData()[last];
1840 for (int comp = 0; comp < NumRealComps(); comp++) {
1841 soa.GetRealData(comp)[pindex] = soa.GetRealData(comp)[last];
1842 }
1843 for (int comp = 0; comp < NumIntComps(); comp++) {
1844 soa.GetIntData(comp)[pindex] = soa.GetIntData(comp)[last];
1845 }
1846 correctCellVectors(last, pindex, grid, ptd[pindex]);
1847 --last;
1848 continue;
1849 }
1850
1851 const int who = ParallelContext::global_to_local_rank(ParticleDistributionMap(pld.m_lev)[pld.m_grid]);
1852 if (who == MyProc) {
1853 if (pld.m_lev != lev || pld.m_grid != grid || pld.m_tile != tile) {
1854 // We own it but must shift it to another place.
1855 auto index = std::make_pair(pld.m_grid, pld.m_tile);
1856 AMREX_ASSERT(soa_local[pld.m_lev][index].size() == num_threads);
1857 {
1858 auto& arr = soa_local[pld.m_lev][index][thread_num].GetIdCPUData();
1859 arr.push_back(soa.GetIdCPUData()[pindex]);
1860 }
1861 for (int comp = 0; comp < NumRealComps(); ++comp) {
1862 RealVector& arr = soa_local[pld.m_lev][index][thread_num].GetRealData(comp);
1863 arr.push_back(soa.GetRealData(comp)[pindex]);
1864 }
1865 for (int comp = 0; comp < NumIntComps(); ++comp) {
1866 IntVector& arr = soa_local[pld.m_lev][index][thread_num].GetIntData(comp);
1867 arr.push_back(soa.GetIntData(comp)[pindex]);
1868 }
1869
1870 p.id() = -p.id(); // Invalidate the particle
1871 }
1872 }
1873 else {
1874 auto& particles_to_send = tmp_remote[who][thread_num];
1875 auto old_size = particles_to_send.size();
1876 auto new_size = old_size + superparticle_size;
1877 particles_to_send.resize(new_size);
1878
1879 char* dst = &particles_to_send[old_size];
1880 {
1881 std::memcpy(dst, &soa.GetIdCPUData()[pindex], sizeof(uint64_t));
1882 dst += sizeof(uint64_t);
1883 }
1884 int array_comp_start = 0;
1885 if constexpr (!ParticleType::is_soa_particle) {
1886 array_comp_start = AMREX_SPACEDIM + NStructReal;
1887 }
1888 for (int comp = 0; comp < NumRealComps(); comp++) {
1889 if (h_redistribute_real_comp[array_comp_start + comp]) {
1890 std::memcpy(dst, &soa.GetRealData(comp)[pindex], sizeof(ParticleReal));
1891 dst += sizeof(ParticleReal);
1892 }
1893 }
1894 array_comp_start = 2 + NStructInt;
1895 for (int comp = 0; comp < NumIntComps(); comp++) {
1896 if (h_redistribute_int_comp[array_comp_start + comp]) {
1897 std::memcpy(dst, &soa.GetIntData(comp)[pindex], sizeof(int));
1898 dst += sizeof(int);
1899 }
1900 }
1901 p.id() = -p.id(); // Invalidate the particle
1902 }
1903
1904 if (p.id() < 0){
1905 soa.GetIdCPUData()[pindex] = soa.GetIdCPUData()[last];
1906 for (int comp = 0; comp < NumRealComps(); comp++) {
1907 soa.GetRealData(comp)[pindex] = soa.GetRealData(comp)[last];
1908 }
1909 for (int comp = 0; comp < NumIntComps(); comp++) {
1910 soa.GetIntData(comp)[pindex] = soa.GetIntData(comp)[last];
1911 }
1912 correctCellVectors(last, pindex, grid, ptd[pindex]);
1913 --last;
1914 continue;
1915 }
1916
1917 ++pindex;
1918 }
1919
1920 {
1921 auto& iddata = soa.GetIdCPUData();
1922 iddata.erase(iddata.begin() + last + 1, iddata.begin() + npart);
1923 }
1924 for (int comp = 0; comp < NumRealComps(); comp++) {
1925 RealVector& rdata = soa.GetRealData(comp);
1926 rdata.erase(rdata.begin() + last + 1, rdata.begin() + npart);
1927 }
1928 for (int comp = 0; comp < NumIntComps(); comp++) {
1929 IntVector& idata = soa.GetIntData(comp);
1930 idata.erase(idata.begin() + last + 1, idata.begin() + npart);
1931 }
1932 }
1933 }
1934 }
1935 }
1936
1937 for (int lev = lev_min; lev <= lev_max; lev++) {
1938 particle_detail::clearEmptyEntries(m_particles[lev]);
1939 }
1940
1941 // Second pass - for each tile in parallel, collect the particles we are owed from all thread's buffers.
1942 for (int lev = lev_min; lev <= lev_max; lev++) {
1943 typename std::map<std::pair<int, int>, Vector<ParticleVector > >::iterator pmap_it;
1944
1945 if constexpr(!ParticleType::is_soa_particle) {
1946 Vector<std::pair<int, int> > grid_tile_ids;
1947 Vector<Vector<ParticleVector>* > pvec_ptrs;
1948
1949 // we need to create any missing map entries in serial here
1950 for (pmap_it=tmp_local[lev].begin(); pmap_it != tmp_local[lev].end(); pmap_it++)
1951 {
1952 DefineAndReturnParticleTile(lev, pmap_it->first.first, pmap_it->first.second);
1953 grid_tile_ids.push_back(pmap_it->first);
1954 pvec_ptrs.push_back(&(pmap_it->second));
1955 }
1956
1957#ifdef AMREX_USE_OMP
1958#pragma omp parallel for
1959#endif
1960 for (int pit = 0; pit < static_cast<int>(pvec_ptrs.size()); ++pit)
1961 {
1962 auto index = grid_tile_ids[pit];
1963 auto& ptile = ParticlesAt(lev, index.first, index.second);
1964 auto& aos = ptile.GetArrayOfStructs();
1965 auto& soa = ptile.GetStructOfArrays();
1966 auto& aos_tmp = *(pvec_ptrs[pit]);
1967 auto& soa_tmp = soa_local[lev][index];
1968 for (int i = 0; i < num_threads; ++i) {
1969 aos.insert(aos.end(), aos_tmp[i].begin(), aos_tmp[i].end());
1970 aos_tmp[i].erase(aos_tmp[i].begin(), aos_tmp[i].end());
1971 for (int comp = 0; comp < NumRealComps(); ++comp) {
1972 RealVector& arr = soa.GetRealData(comp);
1973 RealVector& tmp = soa_tmp[i].GetRealData(comp);
1974 arr.insert(arr.end(), tmp.begin(), tmp.end());
1975 tmp.erase(tmp.begin(), tmp.end());
1976 }
1977 for (int comp = 0; comp < NumIntComps(); ++comp) {
1978 IntVector& arr = soa.GetIntData(comp);
1979 IntVector& tmp = soa_tmp[i].GetIntData(comp);
1980 arr.insert(arr.end(), tmp.begin(), tmp.end());
1981 tmp.erase(tmp.begin(), tmp.end());
1982 }
1983 }
1984 }
1985 } else { // soa particle
1986 Vector<std::pair<int, int> > grid_tile_ids;
1987
1988 // we need to create any missing map entries in serial here
1989 for (auto soa_map_it=soa_local[lev].begin(); soa_map_it != soa_local[lev].end(); soa_map_it++)
1990 {
1991 DefineAndReturnParticleTile(lev, soa_map_it->first.first, soa_map_it->first.second);
1992 grid_tile_ids.push_back(soa_map_it->first);
1993 }
1994
1995#ifdef AMREX_USE_OMP
1996#pragma omp parallel for
1997#endif
1998 for (int pit = 0; pit < static_cast<int>(grid_tile_ids.size()); ++pit) // NOLINT(modernize-loop-convert)
1999 {
2000 auto index = grid_tile_ids[pit];
2001 auto& ptile = ParticlesAt(lev, index.first, index.second);
2002 auto& soa = ptile.GetStructOfArrays();
2003 auto& soa_tmp = soa_local[lev][index];
2004 for (int i = 0; i < num_threads; ++i) {
2005 {
2006 auto& arr = soa.GetIdCPUData();
2007 auto& tmp = soa_tmp[i].GetIdCPUData();
2008 arr.insert(arr.end(), tmp.begin(), tmp.end());
2009 tmp.erase(tmp.begin(), tmp.end());
2010 }
2011 for (int comp = 0; comp < NumRealComps(); ++comp) {
2012 RealVector& arr = soa.GetRealData(comp);
2013 RealVector& tmp = soa_tmp[i].GetRealData(comp);
2014 arr.insert(arr.end(), tmp.begin(), tmp.end());
2015 tmp.erase(tmp.begin(), tmp.end());
2016 }
2017 for (int comp = 0; comp < NumIntComps(); ++comp) {
2018 IntVector& arr = soa.GetIntData(comp);
2019 IntVector& tmp = soa_tmp[i].GetIntData(comp);
2020 arr.insert(arr.end(), tmp.begin(), tmp.end());
2021 tmp.erase(tmp.begin(), tmp.end());
2022 }
2023 }
2024 }
2025 }
2026 }
2027
2028 for (auto& map_it : tmp_remote) {
2029 int who = map_it.first;
2030 not_ours[who];
2031 }
2032
2033 Vector<int> dest_proc_ids;
2034 Vector<Vector<Vector<char> >* > pbuff_ptrs;
2035 for (auto& kv : tmp_remote)
2036 {
2037 dest_proc_ids.push_back(kv.first);
2038 pbuff_ptrs.push_back(&(kv.second));
2039 }
2040
2041#ifdef AMREX_USE_OMP
2042#pragma omp parallel for
2043#endif
2044 for (int pmap_it = 0; pmap_it < static_cast<int>(pbuff_ptrs.size()); ++pmap_it)
2045 {
2046 int who = dest_proc_ids[pmap_it];
2047 Vector<Vector<char> >& tmp = *(pbuff_ptrs[pmap_it]);
2048 for (int i = 0; i < num_threads; ++i) {
2049 not_ours[who].insert(not_ours[who].end(), tmp[i].begin(), tmp[i].end());
2050 tmp[i].erase(tmp[i].begin(), tmp[i].end());
2051 }
2052 }
2053
2054 particle_detail::clearEmptyEntries(not_ours);
2055
2056 if (int(m_particles.size()) > theEffectiveFinestLevel+1) {
2057 // Looks like we lost an AmrLevel on a regrid.
2058 if (m_verbose > 0) {
2059 amrex::Print() << "ParticleContainer::Redistribute() resizing m_particles from "
2060 << m_particles.size() << " to " << theEffectiveFinestLevel+1 << '\n';
2061 }
2062 AMREX_ASSERT(int(m_particles.size()) >= 2);
2063
2064 m_particles.resize(theEffectiveFinestLevel + 1);
2065 m_dummy_mf.resize(theEffectiveFinestLevel + 1);
2066 }
2067
2068 if (ParallelContext::NProcsSub() == 1) {
2069 AMREX_ASSERT(not_ours.empty());
2070 }
2071 else {
2072 RedistributeMPI(not_ours, lev_min, lev_max, nGrow, local);
2073 }
2074
2075 AMREX_ASSERT(OK(lev_min, lev_max, nGrow));
2076
2077 if (m_verbose > 0) {
2078 auto stoptime = amrex::second() - strttime;
2079
2080 ByteSpread();
2081
2082#ifdef AMREX_LAZY
2083 Lazy::QueueReduction( [=] () mutable {
2084#endif
2085 ParallelReduce::Max(stoptime, ParallelContext::IOProcessorNumberSub(),
2086 ParallelContext::CommunicatorSub());
2087
2088 amrex::Print() << "ParticleContainer::Redistribute() time: " << stoptime << "\n\n";
2089#ifdef AMREX_LAZY
2090 });
2091#endif
2092 }
2093}
2094
2095template <typename ParticleType, int NArrayReal, int NArrayInt,
2096 template<class> class Allocator, class CellAssignor>
2097void
2099RedistributeMPI (std::map<int, Vector<char> >& not_ours,
2100 int lev_min, int lev_max, int nGrow, int local)
2101{
2102 BL_PROFILE("ParticleContainer::RedistributeMPI()");
2103 BL_PROFILE_VAR_NS("RedistributeMPI_locate", blp_locate);
2104 BL_PROFILE_VAR_NS("RedistributeMPI_copy", blp_copy);
2105
2106#ifdef AMREX_USE_MPI
2107
2108 using buffer_type = unsigned long long;
2109
2110 std::map<int, Vector<buffer_type> > mpi_snd_data;
2111 for (const auto& kv : not_ours)
2112 {
2113 auto nbt = (kv.second.size() + sizeof(buffer_type)-1)/sizeof(buffer_type);
2114 mpi_snd_data[kv.first].resize(nbt);
2115 std::memcpy((char*) mpi_snd_data[kv.first].data(), kv.second.data(), kv.second.size());
2116 }
2117
2118 const int NProcs = ParallelContext::NProcsSub();
2119 const int NNeighborProcs = neighbor_procs.size();
2120
2121 // We may now have particles that are rightfully owned by another CPU.
2122 Vector<Long> Snds(NProcs, 0), Rcvs(NProcs, 0); // bytes!
2123
2124 Long NumSnds = 0;
2125 if (local > 0)
2126 {
2127 AMREX_ALWAYS_ASSERT(lev_min == 0);
2128 AMREX_ALWAYS_ASSERT(lev_max == 0);
2129 BuildRedistributeMask(0, local);
2130 NumSnds = doHandShakeLocal(not_ours, neighbor_procs, Snds, Rcvs);
2131 }
2132 else
2133 {
2134 NumSnds = doHandShake(not_ours, Snds, Rcvs);
2135 }
2136
2137 const int SeqNum = ParallelDescriptor::SeqNum();
2138
2139 if ((! local) && NumSnds == 0) {
2140 return; // There's no parallel work to do.
2141 }
2142
2143 if (local)
2144 {
2145 Long tot_snds_this_proc = 0;
2146 Long tot_rcvs_this_proc = 0;
2147 for (int i = 0; i < NNeighborProcs; ++i) {
2148 tot_snds_this_proc += Snds[neighbor_procs[i]];
2149 tot_rcvs_this_proc += Rcvs[neighbor_procs[i]];
2150 }
2151 if ( (tot_snds_this_proc == 0) && (tot_rcvs_this_proc == 0) ) {
2152 return; // There's no parallel work to do.
2153 }
2154 }
2155
2156 Vector<int> RcvProc;
2157 Vector<std::size_t> rOffset; // Offset (in bytes) in the receive buffer
2158
2159 std::size_t TotRcvInts = 0;
2160 std::size_t TotRcvBytes = 0;
2161 for (int i = 0; i < NProcs; ++i) {
2162 if (Rcvs[i] > 0) {
2163 RcvProc.push_back(i);
2164 rOffset.push_back(TotRcvInts);
2165 TotRcvBytes += Rcvs[i];
2166 auto nbt = (Rcvs[i] + sizeof(buffer_type)-1)/sizeof(buffer_type);
2167 TotRcvInts += nbt;
2168 }
2169 }
2170
2171 const auto nrcvs = static_cast<int>(RcvProc.size());
2172 Vector<MPI_Status> stats(nrcvs);
2173 Vector<MPI_Request> rreqs(nrcvs);
2174
2175 // Allocate data for rcvs as one big chunk.
2176 Vector<unsigned long long> recvdata(TotRcvInts);
2177
2178 // Post receives.
2179 for (int i = 0; i < nrcvs; ++i) {
2180 const auto Who = RcvProc[i];
2181 const auto offset = rOffset[i];
2182 const auto Cnt = (Rcvs[Who] + sizeof(buffer_type)-1)/sizeof(buffer_type);
2183 AMREX_ASSERT(Cnt > 0);
2184 AMREX_ASSERT(Cnt < size_t(std::numeric_limits<int>::max()));
2185 AMREX_ASSERT(Who >= 0 && Who < NProcs);
2186
2187 rreqs[i] = ParallelDescriptor::Arecv(&recvdata[offset], Cnt, Who, SeqNum,
2188 ParallelContext::CommunicatorSub()).req();
2189 }
2190
2191 // Send.
2192 for (const auto& kv : mpi_snd_data) {
2193 const auto Who = kv.first;
2194 const auto Cnt = kv.second.size();
2195
2196 AMREX_ASSERT(Cnt > 0);
2197 AMREX_ASSERT(Who >= 0 && Who < NProcs);
2198 AMREX_ASSERT(Cnt < std::numeric_limits<int>::max());
2199
2200 ParallelDescriptor::Send(kv.second.data(), Cnt, Who, SeqNum,
2201 ParallelContext::CommunicatorSub());
2202 }
2203
2204 if (nrcvs > 0) {
2205 ParallelDescriptor::Waitall(rreqs, stats);
2206
2207 BL_PROFILE_VAR_START(blp_locate);
2208
2209 int npart = TotRcvBytes / superparticle_size;
2210
2211 Vector<int> rcv_levs(npart);
2212 Vector<int> rcv_grid(npart);
2213 Vector<int> rcv_tile(npart);
2214
2215 int ipart = 0;
2216 ParticleLocData pld;
2217 for (int j = 0; j < nrcvs; ++j)
2218 {
2219 const auto offset = rOffset[j];
2220 const auto Who = RcvProc[j];
2221 const auto Cnt = Rcvs[Who] / superparticle_size;
2222 for (int i = 0; i < int(Cnt); ++i)
2223 {
2224 char* pbuf = ((char*) &recvdata[offset]) + i*superparticle_size;
2225
2227
2228 if constexpr (ParticleType::is_soa_particle) {
2229 std::memcpy(&p.m_idcpu, pbuf, sizeof(uint64_t));
2230
2231 ParticleReal pos[AMREX_SPACEDIM];
2232 std::memcpy(&pos[0], pbuf + sizeof(uint64_t), AMREX_SPACEDIM*sizeof(ParticleReal));
2233 AMREX_D_TERM(p.pos(0) = pos[0];,
2234 p.pos(1) = pos[1];,
2235 p.pos(2) = pos[2]);
2236 } else {
2237 std::memcpy(&p, pbuf, sizeof(ParticleType));
2238 }
2239
2240 bool success = Where(p, pld, lev_min, lev_max, 0);
2241 if (!success)
2242 {
2243 success = (nGrow > 0) && Where(p, pld, lev_min, lev_min, nGrow);
2244 pld.m_grown_gridbox = pld.m_gridbox; // reset grown box for subsequent calls.
2245 }
2246 if (!success)
2247 {
2248 amrex::Abort("RedistributeMPI_locate:: invalid particle.");
2249 }
2250
2251 rcv_levs[ipart] = pld.m_lev;
2252 rcv_grid[ipart] = pld.m_grid;
2253 rcv_tile[ipart] = pld.m_tile;
2254
2255 ++ipart;
2256 }
2257 }
2258
2259 BL_PROFILE_VAR_STOP(blp_locate);
2260
2261 BL_PROFILE_VAR_START(blp_copy);
2262
2263#ifndef AMREX_USE_GPU
2264 ipart = 0;
2265 for (int i = 0; i < nrcvs; ++i)
2266 {
2267 const auto offset = rOffset[i];
2268 const auto Who = RcvProc[i];
2269 const auto Cnt = Rcvs[Who] / superparticle_size;
2270 for (int j = 0; j < int(Cnt); ++j)
2271 {
2272 auto& ptile = m_particles[rcv_levs[ipart]][std::make_pair(rcv_grid[ipart],
2273 rcv_tile[ipart])];
2274 char* pbuf = ((char*) &recvdata[offset]) + j*superparticle_size;
2275
2276 if constexpr (ParticleType::is_soa_particle) {
2277 uint64_t idcpudata;
2278 std::memcpy(&idcpudata, pbuf, sizeof(uint64_t));
2279 pbuf += sizeof(uint64_t);
2280 ptile.GetStructOfArrays().GetIdCPUData().push_back(idcpudata);
2281 } else {
2282 ParticleType p;
2283 std::memcpy(&p, pbuf, sizeof(ParticleType));
2284 pbuf += sizeof(ParticleType);
2285 ptile.push_back(p);
2286 }
2287
2288 int array_comp_start = 0;
2289 if constexpr (!ParticleType::is_soa_particle) {
2290 array_comp_start = AMREX_SPACEDIM + NStructReal;
2291 }
2292 for (int comp = 0; comp < NumRealComps(); ++comp) {
2293 if (h_redistribute_real_comp[array_comp_start + comp]) {
2294 ParticleReal rdata;
2295 std::memcpy(&rdata, pbuf, sizeof(ParticleReal));
2296 pbuf += sizeof(ParticleReal);
2297 ptile.push_back_real(comp, rdata);
2298 } else {
2299 ptile.push_back_real(comp, 0.0);
2300 }
2301 }
2302
2303 array_comp_start = 2 + NStructInt;
2304 for (int comp = 0; comp < NumIntComps(); ++comp) {
2305 if (h_redistribute_int_comp[array_comp_start + comp]) {
2306 int idata;
2307 std::memcpy(&idata, pbuf, sizeof(int));
2308 pbuf += sizeof(int);
2309 ptile.push_back_int(comp, idata);
2310 } else {
2311 ptile.push_back_int(comp, 0);
2312 }
2313 }
2314 ++ipart;
2315 }
2316 }
2317
2318#else
2320 host_particles.reserve(15);
2321 host_particles.resize(finestLevel()+1);
2322
2324 std::vector<Gpu::HostVector<ParticleReal> > > > host_real_attribs;
2325 host_real_attribs.reserve(15);
2326 host_real_attribs.resize(finestLevel()+1);
2327
2329 std::vector<Gpu::HostVector<int> > > > host_int_attribs;
2330 host_int_attribs.reserve(15);
2331 host_int_attribs.resize(finestLevel()+1);
2332
2334 host_idcpu.reserve(15);
2335 host_idcpu.resize(finestLevel()+1);
2336
2337 ipart = 0;
2338 for (int i = 0; i < nrcvs; ++i)
2339 {
2340 const auto offset = rOffset[i];
2341 const auto Who = RcvProc[i];
2342 const auto Cnt = Rcvs[Who] / superparticle_size;
2343 for (auto j = decltype(Cnt)(0); j < Cnt; ++j)
2344 {
2345 int lev = rcv_levs[ipart];
2346 std::pair<int, int> ind(std::make_pair(rcv_grid[ipart], rcv_tile[ipart]));
2347
2348 char* pbuf = ((char*) &recvdata[offset]) + j*superparticle_size;
2349
2350 host_real_attribs[lev][ind].resize(NumRealComps());
2351 host_int_attribs[lev][ind].resize(NumIntComps());
2352
2353 if constexpr (ParticleType::is_soa_particle) {
2354 uint64_t idcpudata;
2355 std::memcpy(&idcpudata, pbuf, sizeof(uint64_t));
2356 pbuf += sizeof(uint64_t);
2357 host_idcpu[lev][ind].push_back(idcpudata);
2358 } else {
2359 ParticleType p;
2360 std::memcpy(&p, pbuf, sizeof(ParticleType));
2361 pbuf += sizeof(ParticleType);
2362 host_particles[lev][ind].push_back(p);
2363 }
2364
2365 host_real_attribs[lev][ind].resize(NumRealComps());
2366 host_int_attribs[lev][ind].resize(NumIntComps());
2367
2368 // add the real...
2369 int array_comp_start = 0;
2370 if constexpr (!ParticleType::is_soa_particle) {
2371 array_comp_start = AMREX_SPACEDIM + NStructReal;
2372 }
2373 for (int comp = 0; comp < NumRealComps(); ++comp) {
2374 if (h_redistribute_real_comp[array_comp_start + comp]) {
2375 Real rdata;
2376 std::memcpy(&rdata, pbuf, sizeof(Real));
2377 pbuf += sizeof(Real);
2378 host_real_attribs[lev][ind][comp].push_back(rdata);
2379 } else {
2380 host_real_attribs[lev][ind][comp].push_back(0.0);
2381 }
2382 }
2383
2384 // ... and int array data
2385 array_comp_start = 2 + NStructInt;
2386 for (int comp = 0; comp < NumIntComps(); ++comp) {
2387 if (h_redistribute_int_comp[array_comp_start + comp]) {
2388 int idata;
2389 std::memcpy(&idata, pbuf, sizeof(int));
2390 pbuf += sizeof(int);
2391 host_int_attribs[lev][ind][comp].push_back(idata);
2392 } else {
2393 host_int_attribs[lev][ind][comp].push_back(0);
2394 }
2395 }
2396 ++ipart;
2397 }
2398 }
2399
2400 for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
2401 {
2402 for (auto& kv : host_particles[host_lev]) {
2403 auto grid = kv.first.first;
2404 auto tile = kv.first.second;
2405 const auto& src_tile = kv.second;
2406
2407 auto& dst_tile = GetParticles(host_lev)[std::make_pair(grid,tile)];
2408 auto old_size = dst_tile.size();
2409 auto new_size = old_size + src_tile.size();
2410 dst_tile.resize(new_size);
2411
2412 if constexpr (ParticleType::is_soa_particle) {
2413 Gpu::copyAsync(Gpu::hostToDevice,
2414 host_idcpu[host_lev][std::make_pair(grid,tile)].begin(),
2415 host_idcpu[host_lev][std::make_pair(grid,tile)].end(),
2416 dst_tile.GetStructOfArrays().GetIdCPUData().begin() + old_size);
2417 } else {
2418 Gpu::copyAsync(Gpu::hostToDevice,
2419 src_tile.begin(), src_tile.end(),
2420 dst_tile.GetArrayOfStructs().begin() + old_size);
2421 }
2422
2423 for (int i = 0; i < NumRealComps(); ++i) {
2424 Gpu::copyAsync(Gpu::hostToDevice,
2425 host_real_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
2426 host_real_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
2427 dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
2428 }
2429
2430 for (int i = 0; i < NumIntComps(); ++i) {
2431 Gpu::copyAsync(Gpu::hostToDevice,
2432 host_int_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
2433 host_int_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
2434 dst_tile.GetStructOfArrays().GetIntData(i).begin() + old_size);
2435 }
2436 }
2437 }
2438
2439 Gpu::Device::streamSynchronize();
2440#endif
2441
2442 BL_PROFILE_VAR_STOP(blp_copy);
2443 }
2444#else
2445 amrex::ignore_unused(not_ours,lev_min,lev_max,nGrow,local);
2446#endif
2447}
2448
2449template <typename ParticleType, int NArrayReal, int NArrayInt,
2450 template<class> class Allocator, class CellAssignor>
2451bool
2453{
2454 BL_PROFILE("ParticleContainer::OK()");
2455
2456 if (lev_max == -1) {
2457 lev_max = finestLevel();
2458}
2459
2460 return (numParticlesOutOfRange(*this, lev_min, lev_max, nGrow) == 0);
2461}
2462
2463template <typename ParticleType, int NArrayReal, int NArrayInt,
2464 template<class> class Allocator, class CellAssignor>
2465void
2467::AddParticlesAtLevel (AoS& particles, int level, int nGrow)
2468{
2469 ParticleTileType ptile;
2470 ptile.GetArrayOfStructs().swap(particles);
2471 AddParticlesAtLevel(ptile, level, nGrow);
2472}
2473
2474template <typename ParticleType, int NArrayReal, int NArrayInt,
2475 template<class> class Allocator, class CellAssignor>
2476void
2478::AddParticlesAtLevel (ParticleTileType& particles, int level, int nGrow)
2479{
2480 BL_PROFILE("ParticleContainer::AddParticlesAtLevel()");
2481
2482 if (int(m_particles.size()) < level+1)
2483 {
2484 if (Verbose())
2485 {
2486 amrex::Print() << "ParticleContainer::AddParticlesAtLevel resizing m_particles from "
2487 << m_particles.size()
2488 << " to "
2489 << level+1 << '\n';
2490 }
2491 m_particles.resize(level+1);
2492 m_dummy_mf.resize(level+1);
2493 for (int lev = 0; lev < level+1; ++lev) {
2494 RedefineDummyMF(lev);
2495 }
2496 }
2497
2498 auto& ptile = DefineAndReturnParticleTile(level, 0, 0);
2499 int old_np = ptile.size();
2500 int num_to_add = particles.size();
2501 int new_np = old_np + num_to_add;
2502 ptile.resize(new_np);
2503 amrex::copyParticles(ptile, particles, 0, old_np, num_to_add);
2504 Redistribute(level, level, nGrow);
2505 particles.resize(0);
2506}
2507
2508// This is the single-level version for cell-centered density
2509template <typename ParticleType, int NArrayReal, int NArrayInt,
2510 template<class> class Allocator, class CellAssignor>
2511void
2513AssignCellDensitySingleLevel (int rho_index,
2514 MultiFab& mf_to_be_filled,
2515 int lev,
2516 int ncomp,
2517 int particle_lvl_offset) const
2518{
2519 BL_PROFILE("ParticleContainer::AssignCellDensitySingleLevel()");
2520
2521 if (rho_index != 0) { amrex::Abort("AssignCellDensitySingleLevel only works if rho_index = 0"); }
2522
2523 MultiFab* mf_pointer;
2524
2525 if (OnSameGrids(lev, mf_to_be_filled)) {
2526 // If we are already working with the internal mf defined on the
2527 // particle_box_array, then we just work with this.
2528 mf_pointer = &mf_to_be_filled;
2529 }
2530 else {
2531 // If mf_to_be_filled is not defined on the particle_box_array, then we need
2532 // to make a temporary here and copy into mf_to_be_filled at the end.
2533 mf_pointer = new MultiFab(ParticleBoxArray(lev),
2534 ParticleDistributionMap(lev),
2535 ncomp, mf_to_be_filled.nGrow());
2536 }
2537
2538 // We must have ghost cells for each FAB so that a particle in one grid can spread
2539 // its effect to an adjacent grid by first putting the value into ghost cells of its
2540 // own grid. The mf->SumBoundary call then adds the value from one grid's ghost cell
2541 // to another grid's valid region.
2542 if (mf_pointer->nGrow() < 1) {
2543 amrex::Error("Must have at least one ghost cell when in AssignCellDensitySingleLevel");
2544 }
2545
2546 const auto strttime = amrex::second();
2547
2548 const auto dxi = Geom(lev).InvCellSizeArray();
2549 const auto plo = Geom(lev).ProbLoArray();
2550 const auto pdxi = Geom(lev + particle_lvl_offset).InvCellSizeArray();
2551
2552 if (Geom(lev).isAnyPeriodic() && ! Geom(lev).isAllPeriodic())
2553 {
2554 amrex::Error("AssignCellDensitySingleLevel: problem must be periodic in no or all directions");
2555 }
2556
2557 mf_pointer->setVal(0);
2558
2560#ifdef AMREX_USE_OMP
2561#pragma omp parallel if (Gpu::notInLaunchRegion())
2562#endif
2563 {
2564 FArrayBox local_rho;
2565 for (ParConstIter pti(*this, lev); pti.isValid(); ++pti) {
2566 const Long np = pti.numParticles();
2567 auto ptd = pti.GetParticleTile().getConstParticleTileData();
2568 FArrayBox& fab = (*mf_pointer)[pti];
2569 auto rhoarr = fab.array();
2570#ifdef AMREX_USE_OMP
2571 Box tile_box;
2572 if (Gpu::notInLaunchRegion())
2573 {
2574 tile_box = pti.tilebox();
2575 tile_box.grow(mf_pointer->nGrow());
2576 local_rho.resize(tile_box,ncomp);
2577 local_rho.setVal<RunOn::Host>(0.0);
2578 rhoarr = local_rho.array();
2579 }
2580#endif
2581
2582 if (particle_lvl_offset == 0)
2583 {
2585 {
2586 auto p = make_particle<ParticleType>{}(ptd,i);
2587 amrex_deposit_cic(p, ncomp, rhoarr, plo, dxi);
2588 });
2589 }
2590 else
2591 {
2593 {
2594 auto p = make_particle<ParticleType>{}(ptd,i);
2595 amrex_deposit_particle_dx_cic(p, ncomp, rhoarr, plo, dxi, pdxi);
2596 });
2597 }
2598
2599#ifdef AMREX_USE_OMP
2600 if (Gpu::notInLaunchRegion())
2601 {
2602 fab.atomicAdd<RunOn::Host>(local_rho, tile_box, tile_box, 0, 0, ncomp);
2603 }
2604#endif
2605 }
2606 }
2607
2608 mf_pointer->SumBoundary(Geom(lev).periodicity());
2609
2610 // If ncomp > 1, first divide the momenta (component n)
2611 // by the mass (component 0) in order to get velocities.
2612 // Be careful not to divide by zero.
2613 for (int n = 1; n < ncomp; n++)
2614 {
2615 for (MFIter mfi(*mf_pointer); mfi.isValid(); ++mfi)
2616 {
2617 (*mf_pointer)[mfi].protected_divide<RunOn::Device>((*mf_pointer)[mfi],0,n,1);
2618 }
2619 }
2620
2621 // Only multiply the first component by (1/vol) because this converts mass
2622 // to density. If there are additional components (like velocity), we don't
2623 // want to divide those by volume.
2624 const Real* dx = Geom(lev).CellSize();
2625 const Real vol = AMREX_D_TERM(dx[0], *dx[1], *dx[2]);
2626
2627 mf_pointer->mult(Real(1.0)/vol, 0, 1, mf_pointer->nGrow());
2628
2629 // If mf_to_be_filled is not defined on the particle_box_array, then we need
2630 // to copy here from mf_pointer into mf_to_be_filled.
2631 if (mf_pointer != &mf_to_be_filled)
2632 {
2633 mf_to_be_filled.ParallelCopy(*mf_pointer,0,0,ncomp,0,0);
2634 delete mf_pointer;
2635 }
2636
2637 if (m_verbose > 1)
2638 {
2639 auto stoptime = amrex::second() - strttime;
2640
2641 ParallelReduce::Max(stoptime, ParallelContext::IOProcessorNumberSub(),
2642 ParallelContext::CommunicatorSub());
2643
2644 amrex::Print() << "ParticleContainer::AssignCellDensitySingleLevel) time: "
2645 << stoptime << '\n';
2646 }
2647}
2648
2649template <typename ParticleType, int NArrayReal, int NArrayInt,
2650 template<class> class Allocator, class CellAssignor>
2651void
2653ResizeRuntimeRealComp (int new_size, bool communicate)
2654{
2655 int old_size = m_num_runtime_real;
2656
2657 m_runtime_comps_defined = (new_size > 0);
2658 m_num_runtime_real = new_size;
2659 int cur_size = h_redistribute_real_comp.size();
2660 h_redistribute_real_comp.resize(cur_size-old_size+new_size, communicate);
2661 SetParticleSize();
2662
2663 for (int lev = 0; lev < numLevels(); ++lev) {
2664 for (ParIterType pti(*this,lev); pti.isValid(); ++pti) {
2665 auto& tile = DefineAndReturnParticleTile(lev, pti);
2666 auto np = tile.numParticles();
2667 if (np > 0 && new_size > old_size) {
2668 auto& soa = tile.GetStructOfArrays();
2669 soa.resize(np);
2670 }
2671 }
2672 }
2673}
2674
2675template <typename ParticleType, int NArrayReal, int NArrayInt,
2676 template<class> class Allocator, class CellAssignor>
2677void
2679ResizeRuntimeIntComp (int new_size, bool communicate)
2680{
2681 int old_size = m_num_runtime_int;
2682
2683 m_runtime_comps_defined = (new_size > 0);
2684 m_num_runtime_int = new_size;
2685 int cur_size = h_redistribute_int_comp.size();
2686 h_redistribute_int_comp.resize(cur_size-old_size+new_size, communicate);
2687 SetParticleSize();
2688
2689 for (int lev = 0; lev < numLevels(); ++lev) {
2690 for (ParIterType pti(*this,lev); pti.isValid(); ++pti) {
2691 auto& tile = DefineAndReturnParticleTile(lev, pti);
2692 auto np = tile.numParticles();
2693 if (np > 0 && new_size > old_size) {
2694 auto& soa = tile.GetStructOfArrays();
2695 soa.resize(np);
2696 }
2697 }
2698 }
2699}
2700
2701}
#define BL_PROFILE_VAR_START(vname)
Definition AMReX_BLProfiler.H:562
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define BL_PROFILE_VAR_STOP(vname)
Definition AMReX_BLProfiler.H:563
#define BL_PROFILE_SYNC_STOP()
Definition AMReX_BLProfiler.H:645
#define BL_PROFILE_SYNC_START_TIMED(fname)
Definition AMReX_BLProfiler.H:644
#define BL_PROFILE_VAR_NS(fname, vname)
Definition AMReX_BLProfiler.H:561
#define AMREX_ALWAYS_ASSERT_WITH_MESSAGE(EX, MSG)
Definition AMReX_BLassert.H:49
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_ALWAYS_ASSERT(EX)
Definition AMReX_BLassert.H:50
#define AMREX_FORCE_INLINE
Definition AMReX_Extension.H:119
#define AMREX_FOR_1D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:97
#define AMREX_HOST_DEVICE_FOR_1D(...)
Definition AMReX_GpuLaunchMacrosC.nolint.H:105
#define AMREX_GPU_DEVICE
Definition AMReX_GpuQualifiers.H:18
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
amrex::ParmParse pp
Input file parser instance for the given namespace.
Definition AMReX_HypreIJIface.cpp:15
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:129
#define AMREX_D_DECL(a, b, c)
Definition AMReX_SPACE.H:104
Print on all processors of the default communicator.
Definition AMReX_Print.H:117
BaseFab< T > & atomicAdd(const BaseFab< T > &x) noexcept
Atomic FAB addition (a[i] <- a[i] + b[i]).
Definition AMReX_BaseFab.H:2954
AMREX_FORCE_INLINE Array4< T const > array() const noexcept
Definition AMReX_BaseFab.H:379
void setVal(T const &x, const Box &bx, int dcomp, int ncomp) noexcept
The setVal functions set sub-regions in the BaseFab to a constant value. This most general form speci...
Definition AMReX_BaseFab.H:1869
A collection of Boxes stored in an Array.
Definition AMReX_BoxArray.H:550
IndexType ixType() const noexcept
Return index type of this BoxArray.
Definition AMReX_BoxArray.H:837
BoxArray & grow(int n)
Grow each Box in the BoxArray by the specified amount.
std::vector< std::pair< int, Box > > intersections(const Box &bx) const
Return intersections of Box and BoxArray.
Box getCellCenteredBox(int index) const noexcept
Return cell-centered box at element index of this BoxArray.
Definition AMReX_BoxArray.H:730
A class for managing a List of Boxes that share a common IndexType. This class implements operations ...
Definition AMReX_BoxList.H:52
BoxList & complementIn(const Box &b, const BoxList &bl)
AMREX_GPU_HOST_DEVICE BoxND & grow(int i) noexcept
Definition AMReX_Box.H:627
AMREX_GPU_HOST_DEVICE bool contains(const IntVectND< dim > &p) const noexcept
Returns true if argument is contained within BoxND.
Definition AMReX_Box.H:204
GpuArray< Real, AMREX_SPACEDIM > InvCellSizeArray() const noexcept
Definition AMReX_CoordSys.H:87
A Fortran Array of REALs.
Definition AMReX_FArrayBox.H:229
void resize(const Box &b, int N=1, Arena *ar=nullptr)
For debugging purposes we hide BaseFab version and do some extra work.
Definition AMReX_FArrayBox.cpp:178
int size() const noexcept
Return the number of FABs in the FabArray.
Definition AMReX_FabArrayBase.H:109
int nGrow(int direction=0) const noexcept
Return the grow factor that defines the region of definition.
Definition AMReX_FabArrayBase.H:77
void ParallelCopy(const FabArray< FAB > &src, const Periodicity &period=Periodicity::NonPeriodic(), CpOp op=FabArrayBase::COPY)
Definition AMReX_FabArray.H:840
void SumBoundary(const Periodicity &period=Periodicity::NonPeriodic())
Sum values in overlapped cells. The destination is limited to valid cells.
Definition AMReX_FabArray.H:3174
Array4< typename FabArray< FAB >::value_type const > array(const MFIter &mfi) const noexcept
Definition AMReX_FabArray.H:560
void setVal(value_type val)
Set all components in the entire region of each FAB to val.
Definition AMReX_FabArray.H:2412
Rectangular problem domain geometry.
Definition AMReX_Geometry.H:73
const Box & Domain() const noexcept
Returns our rectangular domain.
Definition AMReX_Geometry.H:210
GpuArray< Real, AMREX_SPACEDIM > ProbLoArray() const noexcept
Definition AMReX_Geometry.H:186
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool cellCentered() const noexcept
True if the IndexTypeND is CELL based in all directions.
Definition AMReX_IndexType.H:101
AMREX_GPU_HOST_DEVICE static AMREX_FORCE_INLINE constexpr IntVectND< dim > TheZeroVector() noexcept
This static member function returns a reference to a constant IntVectND object, all of whose dim argu...
Definition AMReX_IntVect.H:670
a one-thingy-per-box distributed object
Definition AMReX_LayoutData.H:13
Definition AMReX_MFIter.H:57
int LocalTileIndex() const noexcept
The current local tile index in the current grid;.
Definition AMReX_MFIter.H:150
bool isValid() const noexcept
Is the iterator valid i.e. is it associated with a FAB?
Definition AMReX_MFIter.H:141
Box validbox() const noexcept
Return the valid Box in which the current tile resides.
Definition AMReX_MFIter.H:132
int index() const noexcept
The index into the underlying BoxArray of the current FAB.
Definition AMReX_MFIter.H:144
A collection (stored as an array) of FArrayBox objects.
Definition AMReX_MultiFab.H:38
void mult(Real val, int comp, int num_comp, int nghost=0)
Scales the value of each cell in the specified subregion of the MultiFab by the scalar val (a[i] <- a...
Definition AMReX_MultiFab.cpp:1460
Definition AMReX_PODVector.H:262
size_type size() const noexcept
Definition AMReX_PODVector.H:591
T * dataPtr() noexcept
Definition AMReX_PODVector.H:613
void resize(size_type a_new_size)
Definition AMReX_PODVector.H:641
T * data() noexcept
Definition AMReX_PODVector.H:609
Definition AMReX_ParIter.H:142
Definition AMReX_ParIter.H:113
Parse Parameters From Command Line and Input Files.
Definition AMReX_ParmParse.H:320
int queryarr(const char *name, std::vector< int > &ref, int start_ix=FIRST, int num_val=ALL) const
Same as queryktharr() but searches for last occurrence of name.
Definition AMReX_ParmParse.cpp:1378
int queryAdd(const char *name, T &ref)
If name is found, the value in the ParmParse database will be stored in the ref argument....
Definition AMReX_ParmParse.H:1014
int query(const char *name, bool &ref, int ival=FIRST) const
Same as querykth() but searches for the last occurrence of name.
Definition AMReX_ParmParse.cpp:1311
virtual void reserveData()
Definition AMReX_ParticleContainerBase.cpp:41
virtual void resizeData()
Definition AMReX_ParticleContainerBase.cpp:46
A distributed container for Particles sorted onto the levels, grids, and tiles of a block-structured ...
Definition AMReX_ParticleContainer.H:146
IntVect Index(const P &p, int lev) const
Definition AMReX_ParticleContainerI.H:201
std::map< std::pair< int, int >, ParticleTileType > ParticleLevel
Definition AMReX_ParticleContainer.H:187
void SetSoACompileTimeNames(std::vector< std::string > const &rdata_name, std::vector< std::string > const &idata_name)
Definition AMReX_ParticleContainerI.H:109
typename ParticleTileType::AoS AoS
Definition AMReX_ParticleContainer.H:188
typename SoA::RealVector RealVector
Definition AMReX_ParticleContainer.H:191
void locateParticle(P &p, ParticleLocData &pld, int lev_min, int lev_max, int nGrow, int local_grid=-1) const
Definition AMReX_ParticleContainerI.H:446
void SetParticleSize()
Definition AMReX_ParticleContainerI.H:16
void RemoveParticlesAtLevel(int level)
The Following methods are for managing Virtual and Ghost Particles.
Definition AMReX_ParticleContainerI.H:738
void Increment(MultiFab &mf, int level)
Definition AMReX_ParticleContainerI.H:701
bool HasIntComp(std::string const &name)
Definition AMReX_ParticleContainerI.H:148
void ShrinkToFit()
Definition AMReX_ParticleContainerI.H:681
bool HasRealComp(std::string const &name)
Definition AMReX_ParticleContainerI.H:140
void resizeData() override
This resizes the vector of dummy MultiFabs used by the ParticleContainer for the current number of le...
Definition AMReX_ParticleContainerI.H:435
Long IncrementWithTotal(MultiFab &mf, int level, bool local=false)
Definition AMReX_ParticleContainerI.H:728
Long NumberOfParticlesAtLevel(int level, bool only_valid=true, bool only_local=false) const
Returns # of particles at specified the level.
Definition AMReX_ParticleContainerI.H:551
typename Particle< NStructReal, NStructInt >::RealType RealType
The type of the Real data.
Definition AMReX_ParticleContainer.H:173
int GetRealCompIndex(std::string const &name)
Definition AMReX_ParticleContainerI.H:161
int GetIntCompIndex(std::string const &name)
Definition AMReX_ParticleContainerI.H:183
void RemoveParticlesNotAtFinestLevel()
Definition AMReX_ParticleContainerI.H:752
Vector< Long > NumberOfParticlesInGrid(int level, bool only_valid=true, bool only_local=false) const
Definition AMReX_ParticleContainerI.H:496
typename SoA::IntVector IntVector
Definition AMReX_ParticleContainer.H:192
ParticleLocData Reset(ParticleType &prt, bool update, bool verbose=true, ParticleLocData pld=ParticleLocData()) const
Updates a particle's location (Where), tries to periodic shift any particles that have left the domai...
Definition AMReX_ParticleContainerI.H:392
Long TotalNumberOfParticles(bool only_valid=true, bool only_local=false) const
Returns # of particles at all levels.
Definition AMReX_ParticleContainerI.H:481
void reserveData() override
This reserves data in the vector of dummy MultiFabs used by the ParticleContainer for the maximum num...
Definition AMReX_ParticleContainerI.H:426
T_ParticleType ParticleType
Definition AMReX_ParticleContainer.H:148
Definition AMReX_ParticleLocator.H:104
AssignGrid< BinIteratorFactory > getGridAssignor() const noexcept
Definition AMReX_ParticleLocator.H:183
void build(const BoxArray &ba, const Geometry &geom)
Definition AMReX_ParticleLocator.H:111
This class provides the user with a few print options.
Definition AMReX_Print.H:35
Definition AMReX_Reduce.H:249
Definition AMReX_Reduce.H:364
std::enable_if_t< IsFabArray< MF >::value > eval(MF const &mf, IntVect const &nghost, D &reduce_data, F &&f)
Definition AMReX_Reduce.H:441
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition AMReX_Vector.H:27
Long size() const noexcept
Definition AMReX_Vector.H:50
Definition AMReX_iMultiFab.H:32
Long sum(int comp, int nghost=0, bool local=false) const
Returns the sum in component comp.
Definition AMReX_iMultiFab.cpp:392
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void AddNoRet(T *sum, T value) noexcept
Definition AMReX_GpuAtomic.H:281
void streamSynchronize() noexcept
Definition AMReX_GpuDevice.H:237
void QueueReduction(Func f)
Definition AMReX_Lazy.cpp:7
constexpr Long VirtualParticleID
Definition AMReX_Particle.H:19
void Sum(T &v, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:204
MPI_Comm CommunicatorSub() noexcept
sub-communicator for current frame
Definition AMReX_ParallelContext.H:70
int MyProcSub() noexcept
my sub-rank in current frame
Definition AMReX_ParallelContext.H:76
int IOProcessorNumberSub() noexcept
IO sub-rank in current frame.
Definition AMReX_ParallelContext.H:78
void Bcast(void *, int, MPI_Datatype, int, MPI_Comm)
Definition AMReX_ParallelDescriptor.cpp:1285
void GatherLayoutDataToVector(const LayoutData< T > &sendbuf, Vector< T > &recvbuf, int root)
Gather LayoutData values to a vector on root.
Definition AMReX_ParallelDescriptor.H:1211
void Min(KeyValuePair< K, V > &vi, int root, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:263
void Max(KeyValuePair< K, V > &vi, int root, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:232
void Sum(T &v, int root, MPI_Comm comm)
Definition AMReX_ParallelReduce.H:320
T ExclusiveSum(N n, T const *in, T *out, RetSum a_ret_sum=retSum)
Definition AMReX_Scan.H:1229
static constexpr RetSum retSum
Definition AMReX_Scan.H:29
Definition AMReX_Amr.cpp:49
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void amrex_deposit_cic(P const &p, int nc, amrex::Array4< amrex::Real > const &rho, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi)
Definition AMReX_Particle_mod_K.H:13
int nComp(FabArrayBase const &fa)
void communicateParticlesStart(const PC &pc, ParticleCopyPlan &plan, const SndBuffer &snd_buffer, RcvBuffer &rcv_buffer)
Definition AMReX_ParticleCommunication.H:500
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)
Definition AMReX_CTOParallelForImpl.H:191
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE IntVect getParticleCell(P const &p, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi) noexcept
Returns the cell index for a given particle using the provided lower bounds and cell sizes.
Definition AMReX_ParticleUtil.H:374
int partitionParticlesByDest(PTile &ptile, const PLocator &ploc, CellAssignor const &assignor, const ParticleBufferMap &pmap, const GpuArray< Real, AMREX_SPACEDIM > &plo, const GpuArray< Real, AMREX_SPACEDIM > &phi, const GpuArray< ParticleReal, AMREX_SPACEDIM > &rlo, const GpuArray< ParticleReal, AMREX_SPACEDIM > &rhi, const GpuArray< int, AMREX_SPACEDIM > &is_per, int lev, int gid, int, int lev_min, int lev_max, int nGrow, bool remove_negative)
Definition AMReX_ParticleUtil.H:625
void unpackRemotes(PC &pc, const ParticleCopyPlan &plan, Buffer &rcv_buffer, UnpackPolicy const &policy)
Definition AMReX_ParticleCommunication.H:600
void copyParticles(DstTile &dst, const SrcTile &src) noexcept
Copy particles from src to dst. This version copies all the particles, writing them to the beginning ...
Definition AMReX_ParticleTransformation.H:158
Long doHandShake(const std::map< int, Vector< char > > &not_ours, Vector< Long > &Snds, Vector< Long > &Rcvs)
Definition AMReX_ParticleMPIUtil.cpp:25
Index filterAndTransformParticles(DstTile &dst, const SrcTile &src, Index *mask, F const &f, Index src_start, Index dst_start) noexcept
Conditionally copy particles from src to dst based on the value of mask. A transformation will also b...
Definition AMReX_ParticleTransformation.H:449
bool initialized
Definition AMReX_DistributionMapping.cpp:32
Long doHandShakeLocal(const std::map< int, Vector< char > > &not_ours, const Vector< int > &neighbor_procs, Vector< Long > &Snds, Vector< Long > &Rcvs)
Definition AMReX_ParticleMPIUtil.cpp:50
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 end(BoxND< dim > const &box) noexcept
Definition AMReX_Box.H:1890
double second() noexcept
Definition AMReX_Utility.cpp:922
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int numTilesInBox(const Box &box, const bool a_do_tiling, const IntVect &a_tile_size)
Definition AMReX_ParticleUtil.H:270
void communicateParticlesFinish(const ParticleCopyPlan &plan)
Definition AMReX_ParticleCommunication.cpp:384
Arena * The_Comms_Arena()
Definition AMReX_Arena.cpp:676
Index filterParticles(DstTile &dst, const SrcTile &src, const Index *mask) noexcept
Conditionally copy particles from src to dst based on the value of mask.
Definition AMReX_ParticleTransformation.H:325
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE BoxND< dim > grow(const BoxND< dim > &b, int i) noexcept
Grow BoxND in all directions by given amount.
Definition AMReX_Box.H:1211
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE bool enforcePeriodic(P &p, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &phi, amrex::GpuArray< amrex::ParticleReal, AMREX_SPACEDIM > const &rlo, amrex::GpuArray< amrex::ParticleReal, AMREX_SPACEDIM > const &rhi, amrex::GpuArray< int, AMREX_SPACEDIM > const &is_per) noexcept
Definition AMReX_ParticleUtil.H:459
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition AMReX.H:127
void ParticleToMesh(PC const &pc, const Vector< MultiFab * > &mf, int lev_min, int lev_max, F &&f, bool zero_out_input=true, bool vol_weight=true)
Definition AMReX_AmrParticles.H:156
void Error(const std::string &msg)
Print out message to cerr and exit via amrex::Abort().
Definition AMReX.cpp:224
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void amrex_deposit_particle_dx_cic(P const &p, int nc, amrex::Array4< amrex::Real > const &rho, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &plo, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &dxi, amrex::GpuArray< amrex::Real, AMREX_SPACEDIM > const &pdxi)
Definition AMReX_Particle_mod_K.H:118
void gatherParticles(PTile &dst, const PTile &src, N np, const Index *inds)
Gather particles copies particles into contiguous order from an arbitrary order. Specifically,...
Definition AMReX_ParticleTransformation.H:666
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE Dim3 begin(BoxND< dim > const &box) noexcept
Definition AMReX_Box.H:1881
int numParticlesOutOfRange(Iterator const &pti, int nGrow)
Returns the number of particles that are more than nGrow cells from the box correspond to the input i...
Definition AMReX_ParticleUtil.H:34
int Verbose() noexcept
Definition AMReX.cpp:169
void Abort(const std::string &msg)
Print out message to cerr and exit via abort().
Definition AMReX.cpp:230
const int[]
Definition AMReX_BLProfiler.cpp:1664
int verbose
Definition AMReX_DistributionMapping.cpp:36
void transformParticles(DstTile &dst, const SrcTile &src, F &&f) noexcept
Apply the function f to all the particles in src, writing the result to dst. This version does all th...
Definition AMReX_ParticleTransformation.H:210
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int getTileIndex(const IntVect &iv, const Box &box, const bool a_do_tiling, const IntVect &a_tile_size, Box &tbx)
Definition AMReX_ParticleUtil.H:222
void unpackBuffer(PC &pc, const ParticleCopyPlan &plan, const Buffer &snd_buffer, UnpackPolicy const &policy)
Definition AMReX_ParticleCommunication.H:435
void packBuffer(const PC &pc, const ParticleCopyOp &op, const ParticleCopyPlan &plan, Buffer &snd_buffer)
Definition AMReX_ParticleCommunication.H:336
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void copyParticle(const ParticleTileData< T_ParticleType, NAR, NAI > &dst, const ConstParticleTileData< T_ParticleType, NAR, NAI > &src, int src_i, int dst_i) noexcept
A general single particle copying routine that can run on the GPU.
Definition AMReX_ParticleTransformation.H:31
Definition AMReX_ParticleLocator.H:216
Definition AMReX_Array4.H:61
Definition AMReX_ParticleContainerI.H:1016
amrex::AmrAssignGrid< amrex::DenseBinIteratorFactory< amrex::Box > > m_assign_grid
Definition AMReX_ParticleContainerI.H:1019
int m_lev_max
Definition AMReX_ParticleContainerI.H:1018
AssignGridFilter(amrex::AmrAssignGrid< amrex::DenseBinIteratorFactory< amrex::Box > > assign_grid, int gid, int level, int nGrow)
This filters based on matching grids.
Definition AMReX_ParticleContainerI.H:1025
int m_nGrow
Definition AMReX_ParticleContainerI.H:1018
int m_lev_min
Definition AMReX_ParticleContainerI.H:1018
AMREX_GPU_HOST_DEVICE int operator()(const SrcData &src, int src_i) const noexcept
Definition AMReX_ParticleContainerI.H:1031
int m_gid
Definition AMReX_ParticleContainerI.H:1018
Definition AMReX_ParticleLocator.H:14
Definition AMReX_ParticleTile.H:501
Definition AMReX_ParticleUtil.H:432
Definition AMReX_DenseBins.H:32
Definition AMReX_ParticleContainerI.H:780
GpuArray< Real, AMREX_SPACEDIM > m_dxi
Definition AMReX_ParticleContainerI.H:783
Box m_domain
Definition AMReX_ParticleContainerI.H:784
FilterVirt(const amrex::AssignGrid< amrex::DenseBinIteratorFactory< amrex::Box > > &assign_buffer_grid, const GpuArray< Real, AMREX_SPACEDIM > &plo, const GpuArray< Real, AMREX_SPACEDIM > &dxi, const Box &domain)
Definition AMReX_ParticleContainerI.H:786
GpuArray< Real, AMREX_SPACEDIM > m_plo
Definition AMReX_ParticleContainerI.H:783
AMREX_GPU_HOST_DEVICE int operator()(const SrcData &src, int src_i) const noexcept
Definition AMReX_ParticleContainerI.H:793
amrex::AssignGrid< amrex::DenseBinIteratorFactory< amrex::Box > > m_assign_buffer_grid
Definition AMReX_ParticleContainerI.H:782
Definition AMReX_ParticleUtil.H:341
Definition AMReX_Array.H:34
uint64_t m_idcpu
Definition AMReX_Particle.H:252
Definition AMReX_ParticleCommunication.H:58
void setNumLevels(int num_levels)
Definition AMReX_ParticleCommunication.cpp:14
Vector< std::map< int, Gpu::DeviceVector< IntVect > > > m_periodic_shift
Definition AMReX_ParticleCommunication.H:62
Vector< std::map< int, Gpu::DeviceVector< int > > > m_boxes
Definition AMReX_ParticleCommunication.H:59
Vector< std::map< int, Gpu::DeviceVector< int > > > m_levels
Definition AMReX_ParticleCommunication.H:60
Vector< std::map< int, Gpu::DeviceVector< int > > > m_src_indices
Definition AMReX_ParticleCommunication.H:61
void resize(int gid, int lev, int size)
Definition AMReX_ParticleCommunication.cpp:22
Definition AMReX_ParticleCommunication.H:81
void buildMPIFinish(const ParticleBufferMap &map)
Definition AMReX_ParticleCommunication.cpp:213
void build(const PC &pc, const ParticleCopyOp &op, const Vector< int > &int_comp_mask, const Vector< int > &real_comp_mask, bool local)
Definition AMReX_ParticleCommunication.H:131
A struct used for storing a particle's position in the AMR hierarchy.
Definition AMReX_ParticleContainer.H:92
Box m_grown_gridbox
Definition AMReX_ParticleContainer.H:99
IntVect m_cell
Definition AMReX_ParticleContainer.H:96
int m_grid
Definition AMReX_ParticleContainer.H:94
int m_tile
Definition AMReX_ParticleContainer.H:95
int m_lev
Definition AMReX_ParticleContainer.H:93
Box m_tilebox
Definition AMReX_ParticleContainer.H:98
Box m_gridbox
Definition AMReX_ParticleContainer.H:97
Definition AMReX_ParticleTile.H:702
std::size_t size() const
Returns the total number of particles (real and neighbor)
Definition AMReX_ParticleTile.H:832
ParticleTileDataType getParticleTileData()
Definition AMReX_ParticleTile.H:1137
int numParticles() const
Returns the number of real particles (excluding neighbors)
Definition AMReX_ParticleTile.H:845
AoS & GetArrayOfStructs()
Definition AMReX_ParticleTile.H:820
void resize(std::size_t count)
Definition AMReX_ParticleTile.H:911
bool empty() const
Definition AMReX_ParticleTile.H:826
The struct used to store particles.
Definition AMReX_Particle.H:295
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE ParticleIDWrapper id() &
Definition AMReX_Particle.H:315
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RealVect pos() const &
Definition AMReX_Particle.H:338
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE RealType & rdata(int index) &
Definition AMReX_Particle.H:356
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE ParticleCPUWrapper cpu() &
Definition AMReX_Particle.H:312
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE int & idata(int index) &
Definition AMReX_Particle.H:427
Definition AMReX_ParticleContainerI.H:1043
Definition AMReX_ParticleContainerI.H:801
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE void operator()(DstData &dst, const SrcData &src, int src_i, int dst_i) const noexcept
Definition AMReX_ParticleContainerI.H:804
Definition AMReX_MakeParticle.H:18