Block-Structured AMR Software Framework
AMReX_ParticleInit.H
Go to the documentation of this file.
1 #ifndef AMREX_PARTICLEINIT_H
2 #define AMREX_PARTICLEINIT_H
3 #include <AMReX_Config.H>
4 
5 /*
6  \brief Initialize particles from an Ascii file in the following format:
7 
8  8
9  0.1 0.1 16.2 1000.0 4.0 1.0 6.0
10  8.1 0.1 16.2 1000.0 -5 0.0 -7.0
11  16.1 0.1 16.2 1000.0 6.0 -8.0 2.0
12  24.1 0.1 16.2 1000.0 9.0 4.0 8.0
13  0.1 8.1 16.2 1000.0 -8.0 -3.0 -10.0
14  8.1 8.1 16.2 1000.0 2.0 1.0 0.0
15  16.1 8.1 16.2 1000.0 0.0 2.0 3.0
16  24.1 8.1 16.2 1000.0 -9.0 7.0 5.0
17 
18  The first line is the number of particles. The remaining lines give the particle
19  data to read, one particle on each line. The first AMREX_SPACEDIM components are
20  positions. The next extradata components are the additional real data to read in.
21 
22  Integer data is not currently supported by this function.
23 
24  Parameters:
25 
26  file: the name of the Ascii file with the particles
27  extradata: the number of real components to read in, beyond the AMREX_SPACEDIM
28  positions
29  Nrep: pointer to IntVect that lets you replicate the incoming particles
30  across the domain so that you only need to specify a sub-volume of
31  them. By default particles are not replicated.
32  */
33 template <typename ParticleType, int NArrayReal, int NArrayInt,
34  template<class> class Allocator, class CellAssignor>
35 void
36 ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
37 ::InitFromAsciiFile (const std::string& file, int extradata, const IntVect* Nrep)
38 {
39  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitFromAsciiFile()");
40  AMREX_ASSERT(!file.empty());
41  AMREX_ASSERT(extradata <= NStructReal + NumRealComps());
42 
43  const int MyProc = ParallelDescriptor::MyProc();
44  const int NProcs = ParallelDescriptor::NProcs();
45  const auto strttime = amrex::second();
46 
47  int NReaders = MaxReaders(); // num ranks to read with
48  int NRedist = 1; // number of times to redistribute while reading
49 
50  // try to use some sensible heuristics
51  if (NProcs <= 1024)
52  {
53  if (NReaders > 1)
54  {
55  NRedist = 2;
56  }
57  }
58  else if (NProcs <= 4096)
59  {
60  NReaders = std::max(NReaders,128);
61  NRedist = 4;
62  }
63  else if (NProcs <= 8192)
64  {
65  NReaders = std::max(NReaders,384);
66  NRedist = 32;
67  }
68  else if (NProcs <= 16384)
69  {
70  NReaders = std::max(NReaders,512);
71  NRedist = 48;
72  }
73 
74  resizeData();
75 
76  IntVect lNrep(AMREX_D_DECL(1,1,1));
77 
78  if (Nrep != nullptr) {
79  lNrep = *Nrep;
80  }
81 
82  Long how_many = 0;
83  Long how_many_read = 0;
84 
86  Vector<Gpu::HostVector<Real> > nreals(0);
87  if (extradata > NStructReal) { nreals.resize(extradata - NStructReal); }
88 
89  if (MyProc < NReaders)
90  {
91  VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
92 
93  std::ifstream ifs;
94 
95  ifs.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
96 
97  ifs.open(file.c_str(), std::ios::in);
98 
99  if (!ifs.good())
100  {
101  amrex::FileOpenFailed(file);
102  }
103 
104  int cnt = 0;
105 
106  ifs >> cnt >> std::ws;
107 
108  ParticleLocData pld;
109  ParticleType p, p_rep;
110  Vector<Real> r;
111 
112  if (extradata > NStructReal) { r.resize(extradata - NStructReal); }
113 
114  const int Chunk = cnt / NReaders;
115 
116  for (int i = 0; i < MyProc*Chunk; i++)
117  {
118  ifs.ignore(std::numeric_limits<std::streamsize>::max(), '\n');
119  ifs >> std::ws; // Eat newline.
120  }
121 
122  if (!ifs.good())
123  {
124  std::string msg("ParticleContainer::InitFromAsciiFile(");
125  msg += file;
126  msg += ") failed @ 1";
127  amrex::Error(msg.c_str());
128  }
129 
130  int MyCnt = Chunk;
131 
132  if (MyProc == (NReaders - 1))
133  {
134  MyCnt += cnt % NReaders;
135  }
136 
137  const Geometry& geom = Geom(0);
138 
139  for (int i = 0; i < MyCnt; i++)
140  {
141  AMREX_D_TERM(ifs >> p.pos(0);,
142  ifs >> p.pos(1);,
143  ifs >> p.pos(2););
144 
145  for (int n = 0; n < extradata; n++)
146  {
147  if (n < NStructReal)
148  {
149  ifs >> p.rdata(n);
150  }
151  else
152  {
153  ifs >> r[n - NStructReal];
154  }
155  }
156 
157  if (!ifs.good())
158  {
159  std::string msg("ParticleContainer::InitFromAsciiFile(");
160  msg += file; msg += ") failed @ 2";
161  amrex::Error(msg.c_str());
162  }
163 
164  if (!Where(p, pld))
165  {
166  PeriodicShift(p);
167 
168  if (!Where(p, pld))
169  {
170  if (m_verbose) {
171  amrex::AllPrint() << "BAD PARTICLE ID WOULD BE " << ParticleType::NextID() << '\n'
172  << "BAD PARTICLE POS "
173  << AMREX_D_TERM( p.pos(0),
174  << p.pos(1),
175  << p.pos(2))
176  << "\n";
177  }
178  amrex::Abort("ParticleContainer::InitFromAsciiFile(): invalid particle");
179  }
180  }
181 
182  // set these rather than reading them in
183  p.id() = ParticleType::NextID();
184  p.cpu() = MyProc;
185 
186  nparticles.push_back(p);
187  if(nreals.size() > extradata - NStructReal) {
188  for (int n = NStructReal; n < extradata; n++)
189  {
190  nreals[n-NStructReal].push_back(r[n-NStructReal]);
191  }
192  }
193 
194  how_many++;
195  how_many_read++;
196 
197  const Real DomSize[AMREX_SPACEDIM] =
198  { AMREX_D_DECL((geom.ProbHi(0)-geom.ProbLo(0))/lNrep[0],
199  (geom.ProbHi(1)-geom.ProbLo(1))/lNrep[1],
200  (geom.ProbHi(2)-geom.ProbLo(2))/lNrep[2]) };
201  int rep[AMREX_SPACEDIM];
202 
203 #if AMREX_SPACEDIM > 2
204  for (rep[2] = 1; rep[2] <= lNrep[2]; rep[2]++)
205  {
206 #endif
207 #if AMREX_SPACEDIM > 1
208  for (rep[1] = 1; rep[1] <= lNrep[1]; rep[1]++)
209  {
210 #endif
211  for (rep[0] = 1; rep[0] <= lNrep[0]; rep[0]++)
212  {
213  if (!(AMREX_D_TERM( (rep[0] == 1), && (rep[1] == 1), && (rep[2] == 1) ) ) )
214  {
215  // Shift the position.
216  for (int d=0; d<AMREX_SPACEDIM; ++d)
217  {
218  p_rep.pos(d) = static_cast<ParticleReal>(p.pos(d) + Real(rep[d]-1)*DomSize[d]);
219  }
220 
221  // Copy the extra data
222  for (int n = 0; n < extradata; n++)
223  {
224  if (n < NStructReal)
225  {
226  p_rep.rdata(n) = p.rdata(n);
227  }
228  }
229 
230  if (!Where(p_rep, pld))
231  {
232  PeriodicShift(p_rep);
233  if (!Where(p_rep, pld))
234  {
235  if (m_verbose) {
236  amrex::AllPrint() << "BAD REPLICATED PARTICLE ID WOULD BE " << ParticleType::NextID() << "\n";
237  }
238  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitFromAsciiFile(): invalid replicated particle");
239  }
240  }
241 
242  p_rep.id() = ParticleType::NextID();
243  p_rep.cpu() = MyProc;
244 
245  nparticles.push_back(p_rep);
246 
247  if (nreals.size() > extradata - NStructReal) {
248  for (int n = NStructReal; n < extradata; n++)
249  {
250  nreals[n-NStructReal].push_back(r[n-NStructReal]);
251  }
252  }
253 
254  how_many++;
255  }
256  }
257 #if AMREX_SPACEDIM > 1
258  }
259 #endif
260 #if AMREX_SPACEDIM > 2
261  }
262 #endif
263 
264  }
265 
266  }
267 
268  // Now Redistribute() each chunk separately to minimize memory bloat.
269  int NRedist_chunk = NReaders / NRedist;
270 
271  ParticleLocData pld;
272 
273  for (int nr = 0; nr < NRedist; nr++)
274  {
276  host_particles.reserve(15);
277  host_particles.resize(finestLevel()+1);
278 
279  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<ParticleReal>, NArrayReal > > > host_real_attribs;
280  host_real_attribs.reserve(15);
281  host_real_attribs.resize(finestLevel()+1);
282 
283  if (m_verbose > 0) {
284  amrex::Print() << "Redistributing from processor "
285  << nr*NRedist_chunk << " to "
286  << (nr+1)*NRedist_chunk-1 << '\n';
287  }
288 
289  for (int which = nr*NRedist_chunk; which < (nr+1)*NRedist_chunk; which++)
290  {
291  if (which == MyProc)
292  {
293  while (!nparticles.empty())
294  {
295  ParticleType& p = nparticles.back();
296  Where(p, pld);
297 
298  host_particles[pld.m_lev][std::make_pair(pld.m_grid, pld.m_tile)].push_back(p);
299 
300  if (nreals.size() > extradata - NStructReal && NArrayReal > 0)
301  {
302  for (int n = NStructReal; n < extradata; n++)
303  {
304  Real rdata = nreals[n-NStructReal].back();
305  host_real_attribs[pld.m_lev][std::make_pair(pld.m_grid, pld.m_tile)][n-NStructReal].push_back(rdata);
306  }
307  }
308 
309  nparticles.pop_back();
310 
311  if (nreals.size() > extradata - NStructReal)
312  {
313  for (int n = NStructReal; n < extradata; n++)
314  {
315  nreals[n-NStructReal].pop_back();
316  }
317  }
318  }
319  }
320  }
321 
322  for (int lev = 0; lev < static_cast<int>(host_particles.size()); ++lev)
323  {
324  for (auto& kv : host_particles[lev])
325  {
326 
327  auto grid = kv.first.first;
328  auto tile = kv.first.second;
329  const auto& src_tile = kv.second;
330 
331  auto& dst_tile = GetParticles(lev)[std::make_pair(grid,tile)];
332  auto old_size = dst_tile.GetArrayOfStructs().size();
333  auto new_size = old_size + src_tile.size();
334  dst_tile.resize(new_size);
335 
336  Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
337  dst_tile.GetArrayOfStructs().begin() + old_size);
338 
339  if((host_real_attribs[lev][std::make_pair(grid, tile)]).size() > (long unsigned int) NArrayReal) {
340  for (int i = 0; i < NArrayReal; ++i) {
342  host_real_attribs[lev][std::make_pair(grid,tile)][i].begin(),
343  host_real_attribs[lev][std::make_pair(grid,tile)][i].end(),
344  dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
345  }
346  }
347  }
348  }
350 
351  Redistribute();
352  }
353 
354  if (m_verbose > 0)
355  {
356  if (NRedist*NRedist_chunk < NReaders) {
357  amrex::Print() << "Redistributing from processor "
358  << NRedist*NRedist_chunk << " to "
359  << NReaders << '\n';
360  }
361  }
362  for (int which = NRedist*NRedist_chunk; which < NReaders; which++)
363  {
365  host_particles.reserve(15);
366  host_particles.resize(finestLevel()+1);
367 
368  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<ParticleReal>, NArrayReal > > > host_real_attribs;
369  host_real_attribs.reserve(15);
370  host_real_attribs.resize(finestLevel()+1);
371 
372  if (which == MyProc)
373  {
374  while (!nparticles.empty())
375  {
376  ParticleType& p = nparticles.back();
377  Where(p, pld);
378 
379  host_particles[pld.m_lev][std::make_pair(pld.m_grid, pld.m_tile)].push_back(p);
380  if((host_real_attribs[pld.m_lev][std::make_pair(pld.m_grid, pld.m_tile)]).size() > (long unsigned int) (extradata - NStructReal)) {
381  for (int n = NStructReal; n < extradata; n++)
382  {
383  Real rdata = nreals[n-NStructReal].back();
384  host_real_attribs[pld.m_lev][std::make_pair(pld.m_grid, pld.m_tile)][n-NStructReal].push_back(rdata);
385  }
386  }
387 
388  nparticles.pop_back();
389 
390  if (nreals.size() > extradata - NStructReal) {
391  for (int n = NStructReal; n < extradata; n++)
392  {
393  nreals[n-NStructReal].pop_back();
394  }
395  }
396  }
397  }
398 
399  for (int lev = 0; lev < static_cast<int>(host_particles.size()); ++lev)
400  {
401  for (auto& kv : host_particles[lev])
402  {
403  auto grid = kv.first.first;
404  auto tile = kv.first.second;
405  const auto& src_tile = kv.second;
406 
407  auto& dst_tile = GetParticles(lev)[std::make_pair(grid,tile)];
408  auto old_size = dst_tile.GetArrayOfStructs().size();
409  auto new_size = old_size + src_tile.size();
410  dst_tile.resize(new_size);
411 
412  Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
413  dst_tile.GetArrayOfStructs().begin() + old_size);
414 
415  for (int i = 0; i < NArrayReal; ++i) {
417  host_real_attribs[lev][std::make_pair(grid,tile)][i].begin(),
418  host_real_attribs[lev][std::make_pair(grid,tile)][i].end(),
419  dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
420  }
421  }
422  }
424 
425  Redistribute();
426 
427  }
428 
429  if (m_verbose > 0)
430  {
431  const int IOProcNumber = ParallelDescriptor::IOProcessorNumber();
432 
433  Long num_particles = how_many;
434 
435  ParallelDescriptor::ReduceLongSum(num_particles, IOProcNumber);
436 
437  if (AMREX_D_TERM(lNrep[0] == 1, && lNrep[1] == 1, && lNrep[2] == 1))
438  {
439  amrex::Print() << "Total number of particles: " << num_particles << '\n';
440  }
441  else
442  {
443  Long num_particles_read = how_many_read;
444 
445  ParallelDescriptor::ReduceLongSum(num_particles_read, IOProcNumber);
446 
447  amrex::Print() << "Replication the domain with vector "
448  << AMREX_D_TERM(lNrep[0] << " ", << lNrep[1] << " ", << lNrep[2]) << "\n"
449  << "Total number of particles read in : " << num_particles_read << '\n'
450  << "Total number of particles after replication: " << num_particles << '\n';
451  }
452  }
453 
454  AMREX_ASSERT(OK());
455 
456  if (m_verbose > 1)
457  {
458  ByteSpread();
459 
460  auto runtime = amrex::second() - strttime;
461 
463 
464  amrex::Print() << "InitFromAsciiFile() time: " << runtime << '\n';
465  }
466 }
467 
468 //
469 // The format of a binary particle init file:
470 //
471 // NP -- The number of particles in the file. A "Long".
472 // DM -- Our dimension. Either 1, 2, or 3. A "int".
473 // NX -- The amount of "extra" data. A "int".
474 // NP*(DM+NX) native floating-point numbers. A "float" or "double".
475 //
476 // Note that there is nothing separating all these values.
477 // They're packed into the binary file like sardines.
478 //
479 template <typename ParticleType, int NArrayReal, int NArrayInt,
480  template<class> class Allocator, class CellAssignor>
481 void
483 InitFromBinaryFile (const std::string& file,
484  int extradata)
485 {
486  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitFromBinaryFile()");
487  AMREX_ASSERT(!file.empty());
488  AMREX_ASSERT(extradata <= NStructReal);
489 
490  const int MyProc = ParallelDescriptor::MyProc();
491  const int NProcs = ParallelDescriptor::NProcs();
492  const int IOProc = ParallelDescriptor::IOProcessorNumber();
493  const auto strttime = amrex::second();
494  //
495  // The number of MPI processes that read from the file.
496  // We limit this to a rather small number since there's a limit
497  // to the number of independent I/O channels on most filesystems.
498  //
499  const int NReaders = MaxReaders();
500 
501  AMREX_ASSERT(NReaders <= NProcs);
502  //
503  // How many particles each NReaders reads before redistributing.
504  //
505  const Long NPartPerRedist = MaxParticlesPerRead();
506 
507  if (m_verbose > 0)
508  {
509  amrex::Print() << "Reading with " << NReaders << " readers\n"
510  << "Redistributing after every " << NPartPerRedist << " particles for each reader\n";
511  }
512  //
513  // tmp_particles should mirror how m_particles is built.
514  // At the end of this routine it'll be the new m_particles.
515  //
516  Vector<ParticleLevel> tmp_particles;
517 
518  tmp_particles.reserve(15); // So we don't ever have to do any copying on a resize.
519  tmp_particles.resize(finestLevel()+1);
520 
521  resizeData();
522 
523  //
524  // All the processors need to participate in Redistribute() though.
525  //
526  int NX = 0;
527  Long NP = 0;
528 
529  VisMF::IO_Buffer io_buffer(VisMF::IO_Buffer_Size);
530 
531  std::ifstream ifs;
532 
533  //
534  // The "set" of MPI processor numbers of the readers. We want
535  // to be able to easily check whether or not we're a reader.
536  //
537  std::set<int> readers;
538  //
539  // The same set but in ascending order.
540  //
541  Vector<int> rprocs(NReaders);
542 
543  if (NReaders == NProcs)
544  {
545  //
546  // Just set'm.
547  //
548  for (int i = 0; i < NProcs; i++) {
549  rprocs[i] = i;
550  }
551  }
552  else
553  {
554  //
555  // The I/O Proc builds a set of NReader integers in the range: [0,NProcs-1].
556  //
557  // It then broadcast'm to all MPI procs.
558  //
559  // We want these to be as evenly distributed over the full set of
560  // [0,NProcs-1] MPI processors as possible, so that when reading we
561  // minimize the number of readers per Node, and hence can use more
562  // of the available Node memory for reading.
563  //
565  {
566  do
567  {
568  auto n = int(amrex::Random() * Real(NProcs-1));
569 
570  AMREX_ASSERT(n >= 0);
571  AMREX_ASSERT(n < NProcs);
572 
573  readers.insert(n);
574  }
575  while (static_cast<int>(readers.size()) < NReaders);
576 
577  AMREX_ASSERT(static_cast<Long>(readers.size()) == rprocs.size());
578 
579  int i = 0;
580 
581  for (auto it = readers.cbegin(), End = readers.cend();
582  it != End;
583  ++it, ++i)
584  {
585  rprocs[i] = *it;
586  }
587  }
588 
589  ParallelDescriptor::Bcast(rprocs.dataPtr(), rprocs.size(), IOProc);
590  }
591 
592  if (readers.empty())
593  {
594  //
595  // Set readers for non I/O procs.
596  //
597  readers.insert(rprocs.begin(), rprocs.end());
598 
599  AMREX_ASSERT(static_cast<Long>(readers.size()) == rprocs.size());
600 
601  AMREX_ASSERT(rprocs.size() == NReaders);
602  }
603 
604  int RealSizeInFile = 0;
605 
606  if (readers.find(MyProc) != readers.end())
607  {
608  int DM = 0;
609 
610  ifs.rdbuf()->pubsetbuf(io_buffer.dataPtr(), io_buffer.size());
611 
612  ifs.open(file.c_str(), std::ios::in|std::ios::binary);
613 
614  if (!ifs.good()) {
615  amrex::FileOpenFailed(file);
616  }
617 
618  ifs.read((char*)&NP, sizeof(NP));
619  ifs.read((char*)&DM, sizeof(DM));
620  ifs.read((char*)&NX, sizeof(NX));
621  //
622  // NP MUST be positive!
623  //
624  if (NP <= 0) {
625  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitFromBinaryFile(): NP <= 0");
626  }
627  //
628  // DM must equal AMREX_SPACEDIM.
629  //
630  if (DM != AMREX_SPACEDIM) {
631  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitFromBinaryFile(): DM != AMREX_SPACEDIM");
632  }
633  //
634  // NX MUST be in [0,N].
635  //
636  if (NX < 0 || NX > NStructReal) {
637  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitFromBinaryFile(): NX < 0 || NX > N");
638  }
639  //
640  // Can't ask for more data than exists in the file!
641  //
642  if (extradata > NX) {
643  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitFromBinaryFile(): extradata > NX");
644  }
645  //
646  // Figure out whether we're dealing with floats or doubles.
647  //
648  // First get our current position.
649  //
650  const std::streamoff CURPOS = ifs.tellg();
651  //
652  // Seek to end of file.
653  //
654  ifs.seekg(0,std::ios::end);
655  //
656  // ENDPOS - CURPOS should bracket the particle data.
657  //
658  const std::streamoff ENDPOS = ifs.tellg();
659 
660  RealSizeInFile = int((ENDPOS - CURPOS) / (NP*(DM+NX)));
661 
662  AMREX_ASSERT(RealSizeInFile == sizeof(float) || RealSizeInFile == sizeof(double));
663  //
664  // Now set stream back to earlier position.
665  //
666  ifs.seekg(CURPOS, std::ios::beg);
667  //
668  // Skip to our place in the file.
669  //
670  int id = 0;
671  for ( ; id < NReaders; id++) {
672  if (rprocs[id] == MyProc) {
673  break;
674  }
675  }
676 
677  AMREX_ASSERT(id >= 0 && id < NReaders);
678 
679  const std::streamoff NSKIP = id * (NP/NReaders) * (DM+NX) * RealSizeInFile;
680 
681  if (NSKIP > 0)
682  {
683  ifs.seekg(NSKIP, std::ios::cur);
684  }
685 
686  if (!ifs.good())
687  {
688  std::string msg("ParticleContainer::InitFromBinaryFile(");
689  msg += file;
690  msg += ") failed @ 1";
691  amrex::Error(msg.c_str());
692  }
693  }
694  //
695  // Everyone needs to know NP -- the number of particles in the file.
696  //
698  //
699  // How many particles each reader gets to read.
700  //
701  Long MyCnt = NP / NReaders;
702 
703  if (MyProc == rprocs[0]) {
704  //
705  // Give any remainder to the first reader.
706  //
707  MyCnt += NP % NReaders;
708  }
709 
710  Long how_many_redists = NP / (NPartPerRedist*NReaders), how_many_read = 0;
711 
712  if (NP % (NPartPerRedist*NReaders)) { how_many_redists++; }
713 
714  Vector<float> fxtra, fignore;
715  Vector<double> dxtra, dignore;
716 
717  if (extradata > 0)
718  {
719  fxtra.resize(extradata);
720  dxtra.resize(extradata);
721  }
722 
723  if ((NX-extradata) > 0)
724  {
725  fignore.resize(NX-extradata);
726  dignore.resize(NX-extradata);
727  }
728 
729  ParticleLocData pld;
730 
731  for (int j = 0; j < how_many_redists; j++)
732  {
733 
735  host_particles.reserve(15);
736  host_particles.resize(finestLevel()+1);
737 
738  if (readers.find(MyProc) != readers.end())
739  {
740  ParticleType p;
741 
742  AMREX_ASSERT(MyCnt > how_many_read);
743 
744  const Long NRead = std::min((MyCnt-how_many_read), NPartPerRedist);
745 
746  for (Long i = 0; i < NRead; i++)
747  {
748  //
749  // We don't read in idata.id or idata.cpu. We'll set those later
750  // in a manner to guarantee the global uniqueness of the pair.
751  //
752  if (RealSizeInFile == sizeof(float))
753  {
754  float fpos[AMREX_SPACEDIM];
755 
756  ifs.read((char*)&fpos[0], AMREX_SPACEDIM*sizeof(float));
757 
758  AMREX_D_TERM(p.pos(0) = fpos[0];,
759  p.pos(1) = fpos[1];,
760  p.pos(2) = fpos[2];);
761 
762  }
763  else if (RealSizeInFile == sizeof(double))
764  {
765  double dpos[AMREX_SPACEDIM];
766 
767  ifs.read((char*)&dpos[0], AMREX_SPACEDIM*sizeof(double));
768 
769  AMREX_D_TERM(p.pos(0) = static_cast<ParticleReal>(dpos[0]);,
770  p.pos(1) = static_cast<ParticleReal>(dpos[1]);,
771  p.pos(2) = static_cast<ParticleReal>(dpos[2]););
772  }
773 
774  //
775  // Read in any "extradata".
776  //
777  if (extradata > 0)
778  {
779  if (RealSizeInFile == sizeof(float))
780  {
781  ifs.read((char*)fxtra.data(), std::streamsize(extradata*sizeof(float)));
782 
783  for (int ii = 0; ii < extradata; ii++) {
784  p.rdata(ii) = static_cast<ParticleReal>(fxtra[ii]);
785  }
786  }
787  else if (RealSizeInFile == sizeof(double))
788  {
789  ifs.read((char*)dxtra.data(), std::streamsize(extradata*sizeof(double)));
790 
791  for (int ii = 0; ii < extradata; ii++) {
792  p.rdata(ii) = static_cast<ParticleReal>(dxtra[ii]);
793  }
794  }
795  }
796  //
797  // Read any remaining data for this particle.
798  //
799  if ((NX-extradata) > 0)
800  {
801  if (RealSizeInFile == sizeof(float))
802  {
803  ifs.read((char*)fignore.data(), std::streamsize((NX-extradata)*sizeof(float)));
804  }
805  else if (RealSizeInFile == sizeof(double))
806  {
807  ifs.read((char*)dignore.data(), std::streamsize((NX-extradata)*sizeof(double)));
808  }
809  }
810 
811  if (!ifs.good())
812  {
813  std::string msg("ParticleContainer::InitFromBinaryFile(");
814  msg += file;
815  msg += ") failed @ 2";
816  amrex::Error(msg.c_str());
817  }
818 
819  if (!Where(p, pld))
820  {
821  PeriodicShift(p);
822 
823  if (!Where(p, pld))
824  {
825  if (m_verbose) {
826  amrex::AllPrint() << "BAD PARTICLE ID WOULD BE " << ParticleType::NextID() << '\n'
827  << "BAD PARTICLE POS "
828  << AMREX_D_TERM( p.pos(0),
829  << p.pos(1),
830  << p.pos(2))
831  << "\n";
832  }
833  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitFromBinaryFile(): invalid particle");
834  }
835  }
836 
837  p.id() = ParticleType::NextID();
838  p.cpu() = MyProc;
839 
840  host_particles[pld.m_lev][std::make_pair(pld.m_grid, pld.m_tile)].push_back(p);
841  }
842 
843  how_many_read += NRead;
844  }
845 
846  for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
847  {
848  for (auto& kv : host_particles[host_lev]) {
849  auto grid = kv.first.first;
850  auto tile = kv.first.second;
851  const auto& src_tile = kv.second;
852 
853  auto& dst_tile = GetParticles(host_lev)[std::make_pair(grid,tile)];
854  auto old_size = dst_tile.GetArrayOfStructs().size();
855  auto new_size = old_size + src_tile.size();
856  dst_tile.resize(new_size);
857 
858  Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
859  dst_tile.GetArrayOfStructs().begin() + old_size);
860  }
861  }
863 
864  Redistribute();
865 
866  //
867  // Move particles in m_particles into tmp_particles so that
868  // we don't keep trying to redistribute particles that have
869  // already been redistributed correctly.
870  //
871  for (int lev = 0; lev < m_particles.size(); lev++)
872  {
873  auto& pmap = m_particles[lev];
874  auto& tmp_pmap = tmp_particles[lev];
875 
876  for (auto& kv : pmap) {
877  auto& aos = kv.second.GetArrayOfStructs()();
878  auto& tmp_aos = tmp_pmap[kv.first].GetArrayOfStructs()();
879 
880  tmp_aos.insert(tmp_aos.end(), aos.begin(), aos.end());
881  ParticleVector().swap(aos);
882  }
883 
884  ParticleLevel().swap(pmap);
885  }
886  }
887  //
888  // Make tmp_particles the new m_particles.
889  //
890  tmp_particles.swap(m_particles);
891  //
892  // Add up all the particles read in to get the total number of particles.
893  //
894  if (m_verbose > 0)
895  {
896  Long num_particles_read = how_many_read;
897 
899 
900  amrex::Print() << "\nTotal number of particles: " << num_particles_read << '\n';
901  }
902 
903  AMREX_ASSERT(OK());
904 
905  if (m_verbose > 1)
906  {
907  ByteSpread();
908 
909  auto runtime = amrex::second() - strttime;
910 
912 
913  amrex::Print() << "InitFromBinaryFile() time: " << runtime << '\n';
914  }
915 
917 }
918 
919 //
920 // This function expects to read a file containing the pathnames of
921 // binary particles files needing to be read in for input. It expects
922 // one file name per line.
923 //
924 
925 template <typename ParticleType, int NArrayReal, int NArrayInt,
926  template<class> class Allocator, class CellAssignor>
927 void
929 InitFromBinaryMetaFile (const std::string& metafile,
930  int extradata)
931 {
932  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitFromBinaryMetaFile()");
933  const auto strttime = amrex::second();
934 
935  std::ifstream ifs(metafile.c_str(), std::ios::in);
936 
937  std::string file;
938 
939  for (;;)
940  {
941  std::getline(ifs,file);
942 
943  if (!ifs.good()) { break; }
944 
945  if (m_verbose > 1) {
946  amrex::Print() << "InitFromBinaryMetaFile: processing file: " << file << '\n';
947  }
948 
949  InitFromBinaryFile(file, extradata);
950  }
951 
952  if (m_verbose > 1)
953  {
954  ByteSpread();
955 
956  auto runtime = amrex::second() - strttime;
957 
959 
960  amrex::Print() << "InitFromBinaryMetaFile() time: " << runtime << '\n';
961  }
962 }
963 
964 template <typename ParticleType, int NArrayReal, int NArrayInt,
965  template<class> class Allocator, class CellAssignor>
966 void
968 InitRandom (Long icount,
969  ULong iseed,
970  const ParticleInitData& pdata,
971  bool serialize,
972  RealBox containing_bx)
973 {
974  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitRandom()");
975  AMREX_ASSERT(iseed > 0);
976  AMREX_ASSERT(icount > 0);
977 
978  AMREX_ASSERT(m_gdb != nullptr);
979 
980  const int MyProc = ParallelDescriptor::MyProc();
981  const int NProcs = ParallelDescriptor::NProcs();
982  const int IOProc = ParallelDescriptor::IOProcessorNumber();
983  const auto strttime = amrex::second();
984  const Geometry& geom = Geom(0);
985 
986  Real r, x, len[AMREX_SPACEDIM] = { AMREX_D_DECL(geom.ProbLength(0),
987  geom.ProbLength(1),
988  geom.ProbLength(2)) };
989 
990  // We will enforce that the particles are within the containing_bx.
991  // If containing_bx is not passed in, it defaults to the full domain.
992  if (!containing_bx.ok()) { containing_bx = geom.ProbDomain(); }
993 
994  // containing_bx is assumed to lie within the domain.
995  if (!geom.ProbDomain().contains(containing_bx))
996  {
997  containing_bx.setLo(geom.ProbLo());
998  containing_bx.setHi(geom.ProbHi());
999  }
1000 
1001  const Real* xlo = containing_bx.lo();
1002  const Real* xhi = containing_bx.hi();
1003 
1004  amrex::InitRandom(iseed+MyProc);
1005 
1006  if (serialize)
1007  {
1008  if(icount*AMREX_SPACEDIM >= std::numeric_limits<int>::max())
1009  {
1010  amrex::Abort(
1011  "InitRandom has serialize=true, but this would cause too much "
1012  "particle data to be sent from IOProc. Set serialize=false, "
1013  "or use fewer than " +
1014  std::to_string(
1015  amrex::Math::ceil(
1016  amrex::Real(std::numeric_limits<int>::max()) /
1017  amrex::Real(AMREX_SPACEDIM)
1018  )
1019  ) + " particles."
1020  );
1021  }
1022  //
1023  // We'll let IOProc generate the particles so we get the same
1024  // positions no matter how many CPUs we have. This is here
1025  // mainly for debugging purposes. It's not really useful for
1026  // very large numbers of particles.
1027  //
1028  //
1029  Vector<typename ParticleType::RealType> pos(icount*AMREX_SPACEDIM);
1030 
1032  {
1033  for (Long j = 0; j < icount; j++)
1034  {
1035  for (int i = 0; i < AMREX_SPACEDIM; i++)
1036  {
1037  do
1038  {
1039  r = amrex::Random();
1040  x = geom.ProbLo(i) + (r * len[i]);
1041  }
1042  while (static_cast<ParticleReal>(x) < static_cast<ParticleReal>(xlo[i]) || static_cast<ParticleReal>(x) >= static_cast<ParticleReal>(xhi[i]));
1043 
1044  pos[j*AMREX_SPACEDIM + i] = static_cast<ParticleReal>(x);
1045  }
1046  }
1047  }
1048 
1049  ParallelDescriptor::Bcast(pos.dataPtr(), icount*AMREX_SPACEDIM, IOProc);
1050 
1051  ParticleLocData pld;
1052 
1054  host_particles.reserve(15);
1055  host_particles.resize(finestLevel()+1);
1056 
1057  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<ParticleReal>, NArrayReal > > > host_real_attribs;
1058  host_real_attribs.reserve(15);
1059  host_real_attribs.resize(finestLevel()+1);
1060 
1061  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<int>, NArrayInt > > > host_int_attribs;
1062  host_int_attribs.reserve(15);
1063  host_int_attribs.resize(finestLevel()+1);
1064 
1066  host_idcpu.reserve(15);
1067  host_idcpu.resize(finestLevel()+1);
1068 
1069  for (Long j = 0; j < icount; j++)
1070  {
1071  Particle<0, 0> ptest;
1072 
1073  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1074  ptest.pos(i) = pos[j*AMREX_SPACEDIM + i];
1075  }
1076 
1077  if (!Where(ptest, pld)) {
1078  amrex::Abort("ParticleContainer::InitRandom(): invalid particle");
1079  }
1080 
1081  AMREX_ASSERT(pld.m_lev >= 0 && pld.m_lev <= finestLevel());
1082  std::pair<int, int> ind(pld.m_grid, pld.m_tile);
1083 
1084  const int who = ParticleDistributionMap(pld.m_lev)[pld.m_grid];
1085 
1086  if (who == MyProc) {
1087 
1088  if constexpr(!ParticleType::is_soa_particle)
1089  {
1090  ParticleType p;
1091  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1092  p.pos(i) = pos[j*AMREX_SPACEDIM + i];
1093  }
1094 
1095  // We own it. Add it at the appropriate level.
1096  p.id() = ParticleType::NextID();
1097  p.cpu() = MyProc;
1098 
1099  for (int i = 0; i < NStructInt; i++) {
1100  p.idata(i) = pdata.int_struct_data[i];
1101  }
1102 
1103  for (int i = 0; i < NStructReal; i++) {
1104  p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
1105  }
1106 
1107  // add the struct
1108  host_particles[pld.m_lev][ind].push_back(p);
1109 
1110  // add the real...
1111  for (int i = 0; i < NArrayReal; i++) {
1112  host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
1113  }
1114 
1115  // ... and int array data
1116  for (int i = 0; i < NArrayInt; i++) {
1117  host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
1118  }
1119  } else {
1120  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1121  host_real_attribs[pld.m_lev][ind][i].push_back(pos[j*AMREX_SPACEDIM+i]);
1122  }
1123 
1124  host_idcpu[pld.m_lev][ind].push_back(0);
1125  ParticleIDWrapper(host_idcpu[pld.m_lev][ind].back()) = ParticleType::NextID();
1126  ParticleCPUWrapper(host_idcpu[pld.m_lev][ind].back()) = ParallelDescriptor::MyProc();
1127 
1128  host_particles[pld.m_lev][ind];
1129 
1130  // add the real...
1131  for (int i = AMREX_SPACEDIM; i < NArrayReal; i++) {
1132  host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
1133  }
1134 
1135  // ... and int array data
1136  for (int i = 2; i < NArrayInt; i++) {
1137  host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
1138  }
1139  }
1140  }
1141  }
1142 
1143  for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
1144  {
1145  for (auto& kv : host_particles[host_lev]) {
1146  auto grid = kv.first.first;
1147  auto tile = kv.first.second;
1148  const auto& src_tile = kv.second;
1149 
1150  auto& dst_tile = GetParticles(host_lev)[std::make_pair(grid,tile)];
1151  auto old_size = dst_tile.size();
1152  auto new_size = old_size;
1153  if constexpr(!ParticleType::is_soa_particle)
1154  {
1155  new_size += src_tile.size();
1156  } else {
1157  new_size += host_real_attribs[host_lev][std::make_pair(grid,tile)][0].size();
1158  }
1159  dst_tile.resize(new_size);
1160 
1161  if constexpr(!ParticleType::is_soa_particle)
1162  {
1163  Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
1164  dst_tile.GetArrayOfStructs().begin() + old_size);
1165  } else {
1167  host_idcpu[host_lev][std::make_pair(grid,tile)].begin(),
1168  host_idcpu[host_lev][std::make_pair(grid,tile)].end(),
1169  dst_tile.GetStructOfArrays().GetIdCPUData().begin() + old_size);
1170  }
1171 
1172  for (int i = 0; i < NArrayReal; ++i) { // NOLINT(readability-misleading-indentation)
1174  host_real_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
1175  host_real_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
1176  dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
1177  }
1178 
1179  for (int i = 0; i < NArrayInt; ++i) {
1181  host_int_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
1182  host_int_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
1183  dst_tile.GetStructOfArrays().GetIntData(i).begin() + old_size);
1184  }
1185  }
1186  }
1188 
1189  AMREX_ASSERT(OK());
1190  }
1191  else {
1192  // We'll generate the particles in parallel.
1193  // Each CPU will key off the given seed to get independent streams of random numbers.
1194  Long M = icount / NProcs;
1195  // Processor 0 will get the slop.
1196  if (MyProc == 0) {
1197  M += (icount % NProcs);
1198  }
1199 
1200  ParticleLocData pld;
1201 
1203  host_particles.reserve(15);
1204  host_particles.resize(finestLevel()+1);
1205 
1206  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<ParticleReal>, NArrayReal > > > host_real_attribs;
1207  host_real_attribs.reserve(15);
1208  host_real_attribs.resize(finestLevel()+1);
1209 
1210  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<int>, NArrayInt > > > host_int_attribs;
1211  host_int_attribs.reserve(15);
1212  host_int_attribs.resize(finestLevel()+1);
1213 
1215  host_idcpu.reserve(15);
1216  host_idcpu.resize(finestLevel()+1);
1217 
1218  for (Long icnt = 0; icnt < M; icnt++) {
1219  Particle<0, 0> ptest;
1220  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1221  do {
1222  r = amrex::Random();
1223  x = geom.ProbLo(i) + (r * len[i]);
1224  }
1225  while (static_cast<ParticleReal>(x) < static_cast<ParticleReal>(xlo[i]) || static_cast<ParticleReal>(x) >= static_cast<ParticleReal>(xhi[i]));
1226 
1227  ptest.pos(i) = static_cast<ParticleReal>(x);
1228 
1229  AMREX_ASSERT(ptest.pos(i) < geom.ProbHi(i));
1230  }
1231 
1232  ptest.id() = ParticleType::NextID();
1233  ptest.cpu() = ParallelDescriptor::MyProc();
1234 
1235  // locate the particle
1236  if (!Where(ptest, pld))
1237  {
1238  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitRandom(): invalid particle");
1239  }
1240  AMREX_ASSERT(pld.m_lev >= 0 && pld.m_lev <= finestLevel());
1241  std::pair<int, int> ind(pld.m_grid, pld.m_tile);
1242 
1243  if constexpr(!ParticleType::is_soa_particle)
1244  {
1245  ParticleType p;
1246  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1247  p.pos(i) = ptest.pos(i);;
1248  }
1249 
1250  p.id() = ptest.id();
1251  p.cpu() = ptest.cpu();
1252 
1253  for (int i = 0; i < NStructReal; i++) {
1254  p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
1255  }
1256 
1257  for (int i = 0; i < NStructInt; i++) {
1258  p.idata(i) = pdata.int_struct_data[i];
1259  }
1260 
1261  // add the struct
1262  host_particles[pld.m_lev][ind].push_back(p);
1263 
1264  // add the real...
1265  for (int i = 0; i < NArrayReal; i++) {
1266  host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
1267  }
1268 
1269  // ... and int array data
1270  for (int i = 0; i < NArrayInt; i++) {
1271  host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
1272  }
1273  } else {
1274  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1275  host_real_attribs[pld.m_lev][ind][i].push_back(ptest.pos(i));
1276  }
1277 
1278  host_idcpu[pld.m_lev][ind].push_back(0);
1279  ParticleIDWrapper(host_idcpu[pld.m_lev][ind].back()) = ParticleType::NextID();
1280  ParticleCPUWrapper(host_idcpu[pld.m_lev][ind].back()) = ParallelDescriptor::MyProc();
1281 
1282  host_particles[pld.m_lev][ind];
1283 
1284  // add the real...
1285  for (int i = AMREX_SPACEDIM; i < NArrayReal; i++) {
1286  host_real_attribs[pld.m_lev][ind][i].push_back(static_cast<ParticleReal>(pdata.real_array_data[i]));
1287  }
1288 
1289  // ... and int array data
1290  for (int i = 2; i < NArrayInt; i++) {
1291  host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
1292  }
1293  }
1294  }
1295 
1296  for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
1297  {
1298  for (auto& kv : host_particles[host_lev]) {
1299  auto grid = kv.first.first;
1300  auto tile = kv.first.second;
1301  const auto& src_tile = kv.second;
1302 
1303  auto& dst_tile = GetParticles(host_lev)[std::make_pair(grid,tile)];
1304  auto old_size = dst_tile.size();
1305  auto new_size = old_size;
1306  if constexpr(!ParticleType::is_soa_particle)
1307  {
1308  new_size += src_tile.size();
1309  } else {
1310  new_size += host_real_attribs[host_lev][std::make_pair(grid,tile)][0].size();
1311  }
1312  dst_tile.resize(new_size);
1313 
1314  if constexpr(!ParticleType::is_soa_particle)
1315  {
1316  Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
1317  dst_tile.GetArrayOfStructs().begin() + old_size);
1318  } else {
1320  host_idcpu[host_lev][std::make_pair(grid,tile)].begin(),
1321  host_idcpu[host_lev][std::make_pair(grid,tile)].end(),
1322  dst_tile.GetStructOfArrays().GetIdCPUData().begin() + old_size);
1323  }
1324 
1325  for (int i = 0; i < NArrayReal; ++i) { // NOLINT(readability-misleading-indentation)
1327  host_real_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
1328  host_real_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
1329  dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
1330  }
1331 
1332  for (int i = 0; i < NArrayInt; ++i) {
1334  host_int_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
1335  host_int_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
1336  dst_tile.GetStructOfArrays().GetIntData(i).begin() + old_size);
1337  }
1338  }
1339  }
1341 
1342  // Let Redistribute() sort out where the particles belong.
1343  Redistribute();
1344  }
1345 
1346  if (m_verbose > 1)
1347  {
1348  auto stoptime = amrex::second() - strttime;
1349 
1350  ParallelDescriptor::ReduceRealMax(stoptime,IOProc);
1351 
1352  amrex::Print() << "ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitRandom() time: " << stoptime << '\n';
1353  }
1354 
1356 }
1357 
1358 template <typename ParticleType, int NArrayReal, int NArrayInt,
1359  template<class> class Allocator, class CellAssignor>
1360 void
1362 ::InitRandomPerBox (Long icount_per_box,
1363  ULong iseed,
1364  const ParticleInitData& pdata)
1365 {
1366  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitRandomPerBox()");
1367  AMREX_ASSERT(iseed > 0);
1368  AMREX_ASSERT(icount_per_box > 0);
1369 
1370  AMREX_ASSERT(m_gdb != nullptr);
1371 
1372  const int IOProc = ParallelDescriptor::IOProcessorNumber();
1373  const auto strttime = amrex::second();
1374  const Geometry& geom = Geom(0);
1375 
1376  ParticleLocData pld;
1377  ParticleType p;
1378 
1379  // This assumes level 0 since geom = m_gdb->Geom(0)
1380  const Real* dx = geom.CellSize();
1381 
1382  // We use exactly the same seed for every grid
1383  std::mt19937 mt(iseed);
1384  std::uniform_real_distribution<double> dist(0.0, 1.0);
1385 
1386  m_particles.resize(m_gdb->finestLevel()+1);
1387 
1388  for (int lev = 0; lev < m_particles.size(); lev++)
1389  {
1390  AMREX_ASSERT(m_particles[lev].empty());
1391  }
1392 
1393  // We'll generate the particles in parallel -- but no tiling here.
1394  for (MFIter mfi(*m_dummy_mf[0], false); mfi.isValid(); ++mfi)
1395  {
1396  Box grid = m_gdb->ParticleBoxArray(0)[mfi.index()];
1397  RealBox grid_box = RealBox(grid,dx,geom.ProbLo());
1398 
1399  for (Long icnt = 0; icnt < icount_per_box; icnt++) {
1400  for (Long jcnt = 0; jcnt < icount_per_box; jcnt++) {
1401  for (Long kcnt = 0; kcnt < icount_per_box; kcnt++)
1402  {
1403  AMREX_D_TERM(
1404  p.pos(0) = static_cast<ParticleReal>(grid_box.lo(0) + (dist(mt) + double(icnt)) / double(icount_per_box) * grid_box.length(0));,
1405  p.pos(1) = static_cast<ParticleReal>(grid_box.lo(1) + (dist(mt) + double(jcnt)) / double(icount_per_box) * grid_box.length(1));,
1406  p.pos(2) = static_cast<ParticleReal>(grid_box.lo(2) + (dist(mt) + double(kcnt)) / double(icount_per_box) * grid_box.length(2));
1407  );
1408 
1409  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1410  AMREX_ASSERT(p.pos(i) < grid_box.hi(i));
1411  }
1412 
1413  // the real struct data
1414  for (int i = 0; i < NStructReal; i++) {
1415  p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
1416  }
1417 
1418  // the int struct data
1419  p.id() = ParticleType::NextID();
1420  p.cpu() = ParallelDescriptor::MyProc();
1421 
1422  for (int i = 0; i < NStructInt; i++) {
1423  p.idata(i) = pdata.int_struct_data[i];
1424  }
1425 
1426  // locate the particle
1427  if (!Where(p, pld)) {
1428  amrex::Abort("ParticleContainer::InitRandomPerBox(): invalid particle");
1429  }
1430  AMREX_ASSERT(pld.m_lev >= 0 && pld.m_lev <= finestLevel());
1431  std::pair<int, int> ind(pld.m_grid, pld.m_tile);
1432 
1433  // add the struct
1434  m_particles[pld.m_lev][ind].push_back(p);
1435 
1436  // add the real...
1437  for (int i = 0; i < NArrayReal; i++) {
1438  m_particles[pld.m_lev][ind].push_back_real(i, static_cast<ParticleReal>(pdata.real_array_data[i]));
1439  }
1440 
1441  // ... and int array data
1442  for (int i = 0; i < NArrayInt; i++) {
1443  m_particles[pld.m_lev][ind].push_back_int(i, pdata.int_array_data[i]);
1444  }
1445 
1446  } } }
1447  }
1448 
1449  if (m_verbose > 1)
1450  {
1451  auto stoptime = amrex::second() - strttime;
1452 
1453  ParallelDescriptor::ReduceRealMax(stoptime,IOProc);
1454 
1455  amrex::Print() << "ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitRandomPerBox() time: " << stoptime << '\n';
1456  }
1457 }
1458 
1459 template <typename ParticleType, int NArrayReal, int NArrayInt,
1460  template<class> class Allocator, class CellAssignor>
1461 void
1463 InitOnePerCell (Real x_off, Real y_off, Real z_off, const ParticleInitData& pdata)
1464 {
1465  amrex::ignore_unused(y_off,z_off);
1466 
1467  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitOnePerCell()");
1468 
1469  AMREX_ASSERT(m_gdb != nullptr);
1470 
1471  // Note that x_off, y_off and z_off are the offsets from the lower left corner
1472  // in each cell as measured in units of dx, so they must be in [0,1]
1473  AMREX_ASSERT(x_off >= 0. && y_off >= 0. && z_off >= 0.);
1474  AMREX_ASSERT(x_off <= 1. && y_off <= 1. && z_off <= 1.);
1475 
1476  const int IOProc = ParallelDescriptor::IOProcessorNumber();
1477  const auto strttime = amrex::second();
1478  const Geometry& geom = Geom(0); // generates particles on level 0;
1479 
1480  const Real* dx = geom.CellSize();
1481 
1482  ParticleType p;
1483 
1484  // We'll generate the particles in parallel -- but no tiling of the grid here.
1485  for (MFIter mfi(*m_dummy_mf[0], false); mfi.isValid(); ++mfi) {
1486  Box grid = ParticleBoxArray(0)[mfi.index()];
1487  auto ind = std::make_pair(mfi.index(), mfi.LocalTileIndex());
1488  RealBox grid_box (grid,dx,geom.ProbLo());
1490  for (IntVect beg = grid.smallEnd(), end=grid.bigEnd(), cell = grid.smallEnd(); cell <= end; grid.next(cell))
1491  {
1492  // the real struct data
1493  AMREX_D_TERM(p.pos(0) = static_cast<ParticleReal>(grid_box.lo(0) + (x_off + cell[0]-beg[0])*dx[0]);,
1494  p.pos(1) = static_cast<ParticleReal>(grid_box.lo(1) + (y_off + cell[1]-beg[1])*dx[1]);,
1495  p.pos(2) = static_cast<ParticleReal>(grid_box.lo(2) + (z_off + cell[2]-beg[2])*dx[2]););
1496 
1497  for (int d = 0; d < AMREX_SPACEDIM; ++d) {
1498  AMREX_ASSERT(p.pos(d) < grid_box.hi(d));
1499  }
1500 
1501  for (int i = 0; i < NStructReal; i++) {
1502  p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
1503  }
1504 
1505  // the int struct data
1506  p.id() = ParticleType::NextID();
1507  p.cpu() = ParallelDescriptor::MyProc();
1508 
1509  for (int i = 0; i < NStructInt; i++) {
1510  p.idata(i) = pdata.int_struct_data[i];
1511  }
1512 
1513  // add the struct
1514  ptile_tmp.push_back(p);
1515 
1516  // add the real...
1517  for (int i = 0; i < NArrayReal; i++) {
1518  ptile_tmp.push_back_real(i, static_cast<ParticleReal>(pdata.real_array_data[i]));
1519  }
1520 
1521  // ... and int array data
1522  for (int i = 0; i < NArrayInt; i++) {
1523  ptile_tmp.push_back_int(i, pdata.int_array_data[i]);
1524  }
1525  }
1526 
1527  m_particles[0][ind].resize(ptile_tmp.numParticles());
1528  amrex::copyParticles(m_particles[0][ind], ptile_tmp);
1530  }
1531 
1532  Redistribute();
1533 
1534  if (m_verbose > 1) {
1535  auto stoptime = amrex::second() - strttime;
1536 
1537  ParallelDescriptor::ReduceRealMax(stoptime,IOProc);
1538 
1539  amrex::Print() << "ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitOnePerCell() time: " << stoptime << '\n';
1540  }
1541 }
1542 
1543 template <typename ParticleType, int NArrayReal, int NArrayInt,
1544  template<class> class Allocator, class CellAssignor>
1545 void
1547 InitNRandomPerCell (int n_per_cell, const ParticleInitData& pdata)
1548 {
1549  BL_PROFILE("ParticleContainer<NSR, NSI, NAR, NAI>::InitNRandomPerCell()");
1550 
1551  AMREX_ASSERT(m_gdb != nullptr);
1552 
1553  const int IOProc = ParallelDescriptor::IOProcessorNumber();
1554  const auto strttime = amrex::second();
1555  const Geometry& geom = Geom(0);
1556 
1557  // This assumes level 0 since geom = Geom(0)
1558  const Real* dx = geom.CellSize();
1559 
1560  resizeData();
1561 
1562  for (int lev = 0; lev < m_particles.size(); lev++) {
1563  AMREX_ASSERT(m_particles[lev].empty());
1564  }
1565 
1566  ParticleLocData pld;
1567  ParticleType p;
1568  Real r;
1569 
1570  // We'll generate the particles in parallel -- but no tiling here.
1571  for (MFIter mfi(*m_dummy_mf[0], false); mfi.isValid(); ++mfi)
1572  {
1573  Box grid = ParticleBoxArray(0)[mfi.index()];
1574  RealBox grid_box (grid,dx,geom.ProbLo());
1575 
1577  host_particles.reserve(15);
1578  host_particles.resize(finestLevel()+1);
1579 
1580  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<ParticleReal>, NArrayReal > > > host_real_attribs;
1581  host_real_attribs.reserve(15);
1582  host_real_attribs.resize(finestLevel()+1);
1583 
1584  Vector<std::map<std::pair<int, int>, std::array<Gpu::HostVector<int>, NArrayInt > > > host_int_attribs;
1585  host_int_attribs.reserve(15);
1586  host_int_attribs.resize(finestLevel()+1);
1587 
1588  for (IntVect beg = grid.smallEnd(), end=grid.bigEnd(),
1589  cell = grid.smallEnd(); cell <= end; grid.next(cell)) {
1590 
1591  for (int n = 0; n < n_per_cell; n++)
1592  {
1593  // the real struct data
1594  for (int i = 0; i < AMREX_SPACEDIM; i++) {
1595  constexpr int max_iter = 10;
1596  int iter = 0;
1597  while (iter < max_iter) {
1598  r = amrex::Random();
1599  p.pos(i) = static_cast<ParticleReal>(grid_box.lo(i) + (r + Real(cell[i]-beg[i]))*dx[i]);
1600  if (p.pos(i) < grid_box.hi(i)) { break; }
1601  iter++;
1602  }
1603  AMREX_ASSERT(p.pos(i) < grid_box.hi(i));
1604  }
1605 
1606  for (int i = 0; i < NStructReal; i++) {
1607  p.rdata(i) = static_cast<ParticleReal>(pdata.real_struct_data[i]);
1608  }
1609 
1610  // the int struct data
1611  p.id() = ParticleType::NextID();
1612  p.cpu() = ParallelDescriptor::MyProc();
1613 
1614  for (int i = 0; i < NStructInt; i++) {
1615  p.idata(i) = pdata.int_struct_data[i];
1616  }
1617 
1618  // locate the particle
1619  if (!Where(p, pld)) {
1620  amrex::Abort("ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitNRandomPerCell(): invalid particle");
1621  }
1622  AMREX_ASSERT(pld.m_lev >= 0 && pld.m_lev <= finestLevel());
1623  std::pair<int, int> ind(pld.m_grid, pld.m_tile);
1624 
1625  // add the struct
1626  host_particles[pld.m_lev][ind].push_back(p);
1627 
1628  // add the real...
1629  for (int i = 0; i < NArrayReal; i++) {
1630  host_real_attribs[pld.m_lev][ind][i].push_back(pdata.real_array_data[i]);
1631  }
1632 
1633  // ... and int array data
1634  for (int i = 0; i < NArrayInt; i++) {
1635  host_int_attribs[pld.m_lev][ind][i].push_back(pdata.int_array_data[i]);
1636  }
1637  }
1638  }
1639 
1640  for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
1641  {
1642  for (auto& kv : host_particles[host_lev]) {
1643  auto gid = kv.first.first;
1644  auto tid = kv.first.second;
1645  const auto& src_tid = kv.second;
1646 
1647  auto& dst_tile = GetParticles(host_lev)[std::make_pair(gid,tid)];
1648  auto old_size = dst_tile.GetArrayOfStructs().size();
1649  auto new_size = old_size + src_tid.size();
1650  dst_tile.resize(new_size);
1651 
1652  Gpu::copyAsync(Gpu::hostToDevice, src_tid.begin(), src_tid.end(),
1653  dst_tile.GetArrayOfStructs().begin() + old_size);
1654 
1655  for (int i = 0; i < NArrayReal; ++i)
1656  {
1658  host_real_attribs[host_lev][std::make_pair(gid,tid)][i].begin(),
1659  host_real_attribs[host_lev][std::make_pair(gid,tid)][i].end(),
1660  dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
1661  }
1662 
1663  for (int i = 0; i < NArrayInt; ++i)
1664  {
1666  host_int_attribs[host_lev][std::make_pair(gid,tid)][i].begin(),
1667  host_int_attribs[host_lev][std::make_pair(gid,tid)][i].end(),
1668  dst_tile.GetStructOfArrays().GetIntData(i).begin() + old_size);
1669  }
1670  }
1671  }
1673  }
1674 
1675  if (m_verbose > 1)
1676  {
1677  auto stoptime = amrex::second() - strttime;
1678 
1679  ParallelDescriptor::ReduceRealMax(stoptime,IOProc);
1680 
1681  amrex::Print() << "ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt>::InitNRandomPerCell() time: " << stoptime << '\n';
1682  }
1683 
1685 }
1686 #endif /*AMREX_PARTICLEINIT_H*/
#define BL_PROFILE(a)
Definition: AMReX_BLProfiler.H:551
#define AMREX_ASSERT(EX)
Definition: AMReX_BLassert.H:38
#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
AMREX_GPU_HOST_DEVICE const IntVectND< dim > & smallEnd() const &noexcept
Get the smallend of the BoxND.
Definition: AMReX_Box.H:105
AMREX_GPU_HOST_DEVICE void next(IntVectND< dim > &) const noexcept
Step through the rectangle. It is a runtime error to give a point not inside rectangle....
Definition: AMReX_Box.H:1059
AMREX_GPU_HOST_DEVICE const IntVectND< dim > & bigEnd() const &noexcept
Get the bigend.
Definition: AMReX_Box.H:116
const Real * CellSize() const noexcept
Returns the cellsize for each coordinate direction.
Definition: AMReX_CoordSys.H:71
Rectangular problem domain geometry.
Definition: AMReX_Geometry.H:73
const Real * ProbLo() const noexcept
Returns the lo end of the problem domain in each dimension.
Definition: AMReX_Geometry.H:178
const RealBox & ProbDomain() const noexcept
Returns the problem domain.
Definition: AMReX_Geometry.H:170
const Real * ProbHi() const noexcept
Returns the hi end of the problem domain in each dimension.
Definition: AMReX_Geometry.H:180
Real ProbLength(int dir) const noexcept
Returns length of problem domain in specified dimension.
Definition: AMReX_Geometry.H:208
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
Definition: AMReX_PODVector.H:246
T & back() noexcept
Definition: AMReX_PODVector.H:589
void pop_back() noexcept
Definition: AMReX_PODVector.H:571
bool empty() const noexcept
Definition: AMReX_PODVector.H:579
void push_back(const T &a_value)
Definition: AMReX_PODVector.H:556
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 AoS::ParticleVector ParticleVector
Definition: AMReX_ParticleContainer.H:192
T_ParticleType ParticleType
Definition: AMReX_ParticleContainer.H:147
This class provides the user with a few print options.
Definition: AMReX_Print.H:35
A Box with real dimensions. A RealBox is OK iff volume >= 0.
Definition: AMReX_RealBox.H:21
AMREX_GPU_HOST_DEVICE bool ok() const noexcept
Is the RealBox OK; i.e. does it have non-negative volume?
Definition: AMReX_RealBox.H:81
void setLo(const Real *a_lo) noexcept
Sets lo side.
Definition: AMReX_RealBox.H:64
AMREX_GPU_HOST_DEVICE bool contains(const Real *point, Real eps=0.0) const noexcept
Is the specified point contained in the RealBox?
Definition: AMReX_RealBox.H:102
AMREX_GPU_HOST_DEVICE const Real * hi() const &noexcept
Returns hide side.
Definition: AMReX_RealBox.H:51
AMREX_GPU_HOST_DEVICE const Real * lo() const &noexcept
Returns lo side.
Definition: AMReX_RealBox.H:46
void setHi(const Real *a_hi) noexcept
Sets hi side.
Definition: AMReX_RealBox.H:72
AMREX_GPU_HOST_DEVICE Real length(int dir) const noexcept
Returns length in specified direction.
Definition: AMReX_RealBox.H:62
This class is a thin wrapper around std::vector. Unlike vector, Vector::operator[] provides bound che...
Definition: AMReX_Vector.H:27
Long size() const noexcept
Definition: AMReX_Vector.H:50
T * dataPtr() noexcept
get access to the underlying data pointer
Definition: AMReX_Vector.H:46
AMREX_GPU_HOST_DEVICE Long size(T const &b) noexcept
integer version
Definition: AMReX_GpuRange.H:26
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
int NProcs()
Process ID in MPI_COMM_WORLD.
Definition: AMReX_MPMD.cpp:122
int MyProc()
Definition: AMReX_MPMD.cpp:117
void ReduceLongMax(Long &)
Long max reduction.
Definition: AMReX_ParallelDescriptor.cpp:1224
void ReduceLongSum(Long &)
Long sum reduction.
Definition: AMReX_ParallelDescriptor.cpp:1223
void Bcast(void *, int, MPI_Datatype, int, MPI_Comm)
Definition: AMReX_ParallelDescriptor.cpp:1282
int IOProcessorNumber() noexcept
Definition: AMReX_ParallelDescriptor.H:266
bool IOProcessor() noexcept
Is this CPU the I/O Processor? To get the rank number, call IOProcessorNumber()
Definition: AMReX_ParallelDescriptor.H:275
void ReduceRealMax(Vector< std::reference_wrapper< Real > > const &)
Definition: AMReX_ParallelDescriptor.cpp:1215
@ min
Definition: AMReX_ParallelReduce.H:18
@ max
Definition: AMReX_ParallelReduce.H:17
static constexpr int M
Definition: AMReX_OpenBC.H:13
void FileOpenFailed(const std::string &file)
Output a message and abort when couldn't open the file.
Definition: AMReX_Utility.cpp:131
void InitRandom(ULong cpu_seed, int nprocs, ULong gpu_seed)
Set the seed of the random number generator.
Definition: AMReX_Random.cpp:89
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
Real Random()
Generate a psuedo-random double from uniform distribution.
Definition: AMReX_Random.cpp:123
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 void ignore_unused(const Ts &...)
This shuts up the compiler about unused variables.
Definition: AMReX.H:111
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 Dim3 begin(BoxND< dim > const &box) noexcept
Definition: AMReX_Box.H:1881
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
Definition: AMReX_Particle.H:140
Definition: AMReX_Particle.H:37
A struct used to pass initial data into the various Init methods. This struct is used to pass initial...
Definition: AMReX_ParticleContainer.H:116
std::array< int, NStructInt > int_struct_data
Definition: AMReX_ParticleContainer.H:118
std::array< int, NArrayInt > int_array_data
Definition: AMReX_ParticleContainer.H:120
std::array< double, NArrayReal > real_array_data
Definition: AMReX_ParticleContainer.H:119
std::array< double, NStructReal > real_struct_data
Definition: AMReX_ParticleContainer.H:117
A struct used for storing a particle's position in the AMR hierarchy.
Definition: AMReX_ParticleContainer.H:91
int m_grid
Definition: AMReX_ParticleContainer.H:93
int m_tile
Definition: AMReX_ParticleContainer.H:94
int m_lev
Definition: AMReX_ParticleContainer.H:92
Definition: AMReX_ParticleTile.H:693
void push_back(const ParticleType &p)
Definition: AMReX_ParticleTile.H:914
int numParticles() const
Returns the number of real particles (excluding neighbors)
Definition: AMReX_ParticleTile.H:836
void push_back_int(int comp, int v)
Definition: AMReX_ParticleTile.H:1010
void push_back_real(int comp, ParticleReal v)
Definition: AMReX_ParticleTile.H:958
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 ParticleCPUWrapper cpu() &
Definition: AMReX_Particle.H:312