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