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