Block-Structured AMR Software Framework
 
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Modules Pages
Loading...
Searching...
No Matches
AMReX_ParticleHDF5.H
Go to the documentation of this file.
1#ifndef AMREX_PARTICLEHDF5_H
2#define AMREX_PARTICLEHDF5_H
3#include <AMReX.H>
4#include <AMReX_Config.H>
5
7
8#ifdef AMREX_USE_HDF5
9#include "hdf5.h"
10
11#ifdef AMREX_USE_HDF5_ZFP
12#include "H5Zzfp_lib.h"
13#include "H5Zzfp_props.h"
14#endif
15
16#ifdef AMREX_USE_HDF5_SZ
17#include "H5Z_SZ.h"
18#endif
19
20namespace amrex {
21
22template <typename ParticleType, int NArrayReal, int NArrayInt,
23 template<class> class Allocator, class CellAssignor>
24void
25ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
26::CheckpointHDF5 (const std::string& dir,
27 const std::string& name, bool /*is_checkpoint*/,
28 const Vector<std::string>& real_comp_names,
29 const Vector<std::string>& int_comp_names,
30 const std::string& compression) const
31{
32 Vector<int> write_real_comp;
33 Vector<std::string> tmp_real_comp_names;
34 for (int i = 0; i < NStructReal + NumRealComps(); ++i )
35 {
36 write_real_comp.push_back(1);
37 if (real_comp_names.size() == 0)
38 {
39 std::stringstream ss;
40 ss << "real_comp" << i;
41 tmp_real_comp_names.push_back(ss.str());
42 }
43 else
44 {
45 tmp_real_comp_names.push_back(real_comp_names[i]);
46 }
47 }
48
49 Vector<int> write_int_comp;
50 Vector<std::string> tmp_int_comp_names;
51 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
52 {
53 write_int_comp.push_back(1);
54 if (int_comp_names.size() == 0)
55 {
56 std::stringstream ss;
57 ss << "int_comp" << i;
58 tmp_int_comp_names.push_back(ss.str());
59 }
60 else
61 {
62 tmp_int_comp_names.push_back(int_comp_names[i]);
63 }
64 }
65
66 WriteHDF5ParticleData(dir, name, write_real_comp, write_int_comp,
67 tmp_real_comp_names, tmp_int_comp_names,
68 compression,
69 [=] AMREX_GPU_HOST_DEVICE (const SuperParticleType& p) -> int
70 {
71 return p.id() > 0;
72 }, true);
73}
74
75template <typename ParticleType, int NArrayReal, int NArrayInt,
76 template<class> class Allocator, class CellAssignor>
77void
78ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
79::CheckpointHDF5 (const std::string& dir, const std::string& name,
80 const std::string& compression) const
81{
82 Vector<int> write_real_comp;
83 Vector<std::string> real_comp_names;
84 for (int i = 0; i < NStructReal + NumRealComps(); ++i )
85 {
86 write_real_comp.push_back(1);
87 std::stringstream ss;
88 ss << "real_comp" << i;
89 real_comp_names.push_back(ss.str());
90 }
91
92 Vector<int> write_int_comp;
93 Vector<std::string> int_comp_names;
94 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
95 {
96 write_int_comp.push_back(1);
97 std::stringstream ss;
98 ss << "int_comp" << i;
99 int_comp_names.push_back(ss.str());
100 }
101
102 WriteHDF5ParticleData(dir, name, write_real_comp, write_int_comp,
103 real_comp_names, int_comp_names,
104 compression,
105 [=] AMREX_GPU_HOST_DEVICE (const SuperParticleType& p) -> int
106 {
107 return p.id() > 0;
108 });
109}
110
111template <typename ParticleType, int NArrayReal, int NArrayInt,
112 template<class> class Allocator, class CellAssignor>
113void
114ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
115::WritePlotFileHDF5 (const std::string& dir, const std::string& name,
116 const std::string& compression) const
117{
118 Vector<int> write_real_comp;
119 Vector<std::string> real_comp_names;
120 for (int i = 0; i < NStructReal + NumRealComps(); ++i )
121 {
122 write_real_comp.push_back(1);
123 std::stringstream ss;
124 ss << "real_comp" << i;
125 real_comp_names.push_back(ss.str());
126 }
127
128 Vector<int> write_int_comp;
129 Vector<std::string> int_comp_names;
130 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
131 {
132 write_int_comp.push_back(1);
133 std::stringstream ss;
134 ss << "int_comp" << i;
135 int_comp_names.push_back(ss.str());
136 }
137
138 WriteHDF5ParticleData(dir, name, write_real_comp, write_int_comp,
139 real_comp_names, int_comp_names,
140 compression,
141 [=] AMREX_GPU_HOST_DEVICE (const SuperParticleType& p)
142 {
143 return p.id() > 0;
144 });
145}
146
147template <typename ParticleType, int NArrayReal, int NArrayInt,
148 template<class> class Allocator, class CellAssignor>
149void
150ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
151::WritePlotFileHDF5 (const std::string& dir, const std::string& name,
152 const Vector<std::string>& real_comp_names,
153 const Vector<std::string>& int_comp_names,
154 const std::string& compression) const
155{
156 AMREX_ASSERT(real_comp_names.size() == NStructReal + NumRealComps());
157 AMREX_ASSERT( int_comp_names.size() == NStructInt + NumIntComps() );
158
159 Vector<int> write_real_comp;
160 for (int i = 0; i < NStructReal + NumRealComps(); ++i) { write_real_comp.push_back(1); }
161
162 Vector<int> write_int_comp;
163 for (int i = 0; i < NStructInt + NumIntComps(); ++i) { write_int_comp.push_back(1); }
164
165 WriteHDF5ParticleData(dir, name,
166 write_real_comp, write_int_comp,
167 real_comp_names, int_comp_names,
168 compression,
169 [=] AMREX_GPU_HOST_DEVICE (const SuperParticleType& p)
170 {
171 return p.id() > 0;
172 });
173}
174
175template <typename ParticleType, int NArrayReal, int NArrayInt,
176 template<class> class Allocator, class CellAssignor>
177void
178ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
179::WritePlotFileHDF5 (const std::string& dir, const std::string& name,
180 const Vector<std::string>& real_comp_names,
181 const std::string& compression) const
182{
183 AMREX_ASSERT(real_comp_names.size() == NStructReal + NumRealComps());
184
185 Vector<int> write_real_comp;
186 for (int i = 0; i < NStructReal + NumRealComps(); ++i) write_real_comp.push_back(1);
187
188 Vector<int> write_int_comp;
189 for (int i = 0; i < NStructInt + NumIntComps(); ++i) write_int_comp.push_back(1);
190
191 Vector<std::string> int_comp_names;
192 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
193 {
194 std::stringstream ss;
195 ss << "int_comp" << i;
196 int_comp_names.push_back(ss.str());
197 }
198
199 WriteHDF5ParticleData(dir, name,
200 write_real_comp, write_int_comp,
201 real_comp_names, int_comp_names,
202 compression,
203 [=] AMREX_GPU_HOST_DEVICE (const ParticleType& p)
204 {
205 return p.id() > 0;
206 });
207}
208
209template <typename ParticleType, int NArrayReal, int NArrayInt,
210 template<class> class Allocator, class CellAssignor>
211void
212ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
213::WritePlotFileHDF5 (const std::string& dir,
214 const std::string& name,
215 const Vector<int>& write_real_comp,
216 const Vector<int>& write_int_comp,
217 const std::string& compression) const
218{
219 AMREX_ASSERT(write_real_comp.size() == NStructReal + NumRealComps());
220 AMREX_ASSERT(write_int_comp.size() == NStructInt + NArrayInt );
221
222 Vector<std::string> real_comp_names;
223 for (int i = 0; i < NStructReal + NumRealComps(); ++i )
224 {
225 std::stringstream ss;
226 ss << "real_comp" << i;
227 real_comp_names.push_back(ss.str());
228 }
229
230 Vector<std::string> int_comp_names;
231 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
232 {
233 std::stringstream ss;
234 ss << "int_comp" << i;
235 int_comp_names.push_back(ss.str());
236 }
237
238 WriteHDF5ParticleData(dir, name, write_real_comp, write_int_comp,
239 real_comp_names, int_comp_names,
240 compression,
241 [=] AMREX_GPU_HOST_DEVICE (const SuperParticleType& p)
242 {
243 return p.id() > 0;
244 });
245}
246
247template <typename ParticleType, int NArrayReal, int NArrayInt,
248 template<class> class Allocator, class CellAssignor>
249void
250ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>::
251WritePlotFileHDF5 (const std::string& dir, const std::string& name,
252 const Vector<int>& write_real_comp,
253 const Vector<int>& write_int_comp,
254 const Vector<std::string>& real_comp_names,
255 const Vector<std::string>& int_comp_names,
256 const std::string& compression) const
257{
258 BL_PROFILE("ParticleContainer::WritePlotFile()");
259
260 WriteHDF5ParticleData(dir, name,
261 write_real_comp, write_int_comp,
262 real_comp_names, int_comp_names,
263 compression,
264 [=] AMREX_GPU_HOST_DEVICE (const SuperParticleType& p)
265 {
266 return p.id() > 0;
267 });
268}
269
270template <typename ParticleType, int NArrayReal, int NArrayInt,
271 template<class> class Allocator, class CellAssignor>
272template <class F, std::enable_if_t<!std::is_same_v<F, Vector<std::string>>>*>
273void
274ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
275::WritePlotFileHDF5 (const std::string& dir, const std::string& name,
276 const std::string& compression, F&& f) const
277{
278 Vector<int> write_real_comp;
279 Vector<std::string> real_comp_names;
280 for (int i = 0; i < NStructReal + NumRealComps(); ++i )
281 {
282 write_real_comp.push_back(1);
283 std::stringstream ss;
284 ss << "real_comp" << i;
285 real_comp_names.push_back(ss.str());
286 }
287
288 Vector<int> write_int_comp;
289 Vector<std::string> int_comp_names;
290 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
291 {
292 write_int_comp.push_back(1);
293 std::stringstream ss;
294 ss << "int_comp" << i;
295 int_comp_names.push_back(ss.str());
296 }
297
298 WriteHDF5ParticleData(dir, name, write_real_comp, write_int_comp,
299 real_comp_names, int_comp_names, compression,
300 std::forward<F>(f));
301}
302
303template <typename ParticleType, int NArrayReal, int NArrayInt,
304 template<class> class Allocator, class CellAssignor>
305template <class F>
306void
307ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
308::WritePlotFileHDF5 (const std::string& dir, const std::string& name,
309 const Vector<std::string>& real_comp_names,
310 const Vector<std::string>& int_comp_names,
311 const std::string& compression, F&& f) const
312{
313 AMREX_ASSERT(real_comp_names.size() == NStructReal + NumRealComps());
314 AMREX_ASSERT( int_comp_names.size() == NStructInt + NArrayInt );
315
316 Vector<int> write_real_comp;
317 for (int i = 0; i < NStructReal + NumRealComps(); ++i) { write_real_comp.push_back(1); }
318
319 Vector<int> write_int_comp;
320 for (int i = 0; i < NStructInt + NumIntComps(); ++i) { write_int_comp.push_back(1); }
321
322 WriteHDF5ParticleData(dir, name,
323 write_real_comp, write_int_comp,
324 real_comp_names, int_comp_names,
325 compression, std::forward<F>(f));
326}
327
328template <typename ParticleType, int NArrayReal, int NArrayInt,
329 template<class> class Allocator, class CellAssignor>
330template <class F, std::enable_if_t<!std::is_same_v<F, Vector<std::string>>>*>
331void
332ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
333::WritePlotFileHDF5 (const std::string& dir, const std::string& name,
334 const Vector<std::string>& real_comp_names,
335 const std::string& compression, F&& f) const
336{
337 AMREX_ASSERT(real_comp_names.size() == NStructReal + NumRealComps());
338
339 Vector<int> write_real_comp;
340 for (int i = 0; i < NStructReal + NumRealComps(); ++i) { write_real_comp.push_back(1); }
341
342 Vector<int> write_int_comp;
343 for (int i = 0; i < NStructInt + NumIntComps(); ++i) { write_int_comp.push_back(1); }
344
345 Vector<std::string> int_comp_names;
346 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
347 {
348 std::stringstream ss;
349 ss << "int_comp" << i;
350 int_comp_names.push_back(ss.str());
351 }
352
353 WriteHDF5ParticleData(dir, name,
354 write_real_comp, write_int_comp,
355 real_comp_names, int_comp_names,
356 compression, std::forward<F>(f));
357}
358
359template <typename ParticleType, int NArrayReal, int NArrayInt,
360 template<class> class Allocator, class CellAssignor>
361template <class F>
362void
363ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
364::WritePlotFileHDF5 (const std::string& dir,
365 const std::string& name,
366 const Vector<int>& write_real_comp,
367 const Vector<int>& write_int_comp,
368 const std::string& compression, F&& f) const
369{
370 AMREX_ASSERT(write_real_comp.size() == NStructReal + NumRealComps());
371 AMREX_ASSERT(write_int_comp.size() == NStructInt + NumIntComps() );
372
373 Vector<std::string> real_comp_names;
374 for (int i = 0; i < NStructReal + NumRealComps(); ++i )
375 {
376 std::stringstream ss;
377 ss << "real_comp" << i;
378 real_comp_names.push_back(ss.str());
379 }
380
381 Vector<std::string> int_comp_names;
382 for (int i = 0; i < NStructInt + NumIntComps(); ++i )
383 {
384 std::stringstream ss;
385 ss << "int_comp" << i;
386 int_comp_names.push_back(ss.str());
387 }
388
389 WriteHDF5ParticleData(dir, name, write_real_comp, write_int_comp,
390 real_comp_names, int_comp_names,
391 compression, std::forward<F>(f));
392}
393
394template <typename ParticleType, int NArrayReal, int NArrayInt,
395 template<class> class Allocator, class CellAssignor>
396template <class F>
397void
398ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>::
399WritePlotFileHDF5 (const std::string& dir, const std::string& name,
400 const Vector<int>& write_real_comp,
401 const Vector<int>& write_int_comp,
402 const Vector<std::string>& real_comp_names,
403 const Vector<std::string>& int_comp_names,
404 const std::string& compression, F&& f) const
405{
406 BL_PROFILE("ParticleContainer::WritePlotFile()");
407
408 WriteHDF5ParticleData(dir, name,
409 write_real_comp, write_int_comp,
410 real_comp_names, int_comp_names,
411 compression, std::forward<F>(f));
412}
413
414template <typename ParticleType, int NArrayReal, int NArrayInt,
415 template<class> class Allocator, class CellAssignor>
416template <class F>
417void
418ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
419::WriteHDF5ParticleData (const std::string& dir, const std::string& name,
420 const Vector<int>& write_real_comp,
421 const Vector<int>& write_int_comp,
422 const Vector<std::string>& real_comp_names,
423 const Vector<std::string>& int_comp_names,
424 const std::string& compression,
425 F&& f, bool is_checkpoint) const
426{
427 /* HDF5 async is implemented in WriteHDF5ParticleDataSync, enabled by compile flag */
428 /* if (AsyncOut::UseAsyncOut()) { */
429 /* WriteHDF5ParticleDataAsync(*this, dir, name, */
430 /* write_real_comp, write_int_comp, */
431 /* real_comp_names, int_comp_names); */
432 /* } else */
433 /* { */
434 WriteHDF5ParticleDataSync(*this, dir, name,
435 write_real_comp, write_int_comp,
436 real_comp_names, int_comp_names,
437 compression,
438 std::forward<F>(f), is_checkpoint);
439 /* } */
440}
441
442template <typename ParticleType, int NArrayReal, int NArrayInt,
443 template<class> class Allocator, class CellAssignor>
444void
445ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
446::CheckpointPreHDF5 ()
447{
448 if( ! usePrePost) {
449 return;
450 }
451
452 BL_PROFILE("ParticleContainer::CheckpointPre()");
453
454 const int IOProcNumber = ParallelDescriptor::IOProcessorNumber();
455 Long nparticles = 0;
456 Long maxnextid = ParticleType::NextID();
457
458 for (int lev = 0; lev < m_particles.size(); lev++) {
459 const auto& pmap = m_particles[lev];
460 for (const auto& kv : pmap) {
461 const auto& aos = kv.second.GetArrayOfStructs();
462 for (int k = 0; k < aos.numParticles(); ++k) {
463 const ParticleType& p = aos[k];
464 if (p.id() > 0) {
465 //
466 // Only count (and checkpoint) valid particles.
467 //
468 nparticles++;
469 }
470 }
471 }
472 }
473 ParallelDescriptor::ReduceLongSum(nparticles, IOProcNumber);
474
475 ParticleType::NextID(maxnextid);
476 ParallelDescriptor::ReduceLongMax(maxnextid, IOProcNumber);
477
478 nparticlesPrePost = nparticles;
479 maxnextidPrePost = maxnextid;
480
481 nParticlesAtLevelPrePost.clear();
482 nParticlesAtLevelPrePost.resize(finestLevel() + 1, 0);
483 for(int lev(0); lev <= finestLevel(); ++lev) {
484 nParticlesAtLevelPrePost[lev] = NumberOfParticlesAtLevel(lev);
485 }
486
487 whichPrePost.clear();
488 whichPrePost.resize(finestLevel() + 1);
489 countPrePost.clear();
490 countPrePost.resize(finestLevel() + 1);
491 wherePrePost.clear();
492 wherePrePost.resize(finestLevel() + 1);
493
494 filePrefixPrePost.clear();
495 filePrefixPrePost.resize(finestLevel() + 1);
496}
497
498
499template <typename ParticleType, int NArrayReal, int NArrayInt,
500 template<class> class Allocator, class CellAssignor>
501void
502ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
503::CheckpointPostHDF5 ()
504{
505 if( ! usePrePost) {
506 return;
507 }
508
509 BL_PROFILE("ParticleContainer::CheckpointPostHDF5()");
510
511 const int IOProcNumber = ParallelDescriptor::IOProcessorNumber();
512 std::ofstream HdrFile;
513 HdrFile.open(HdrFileNamePrePost.c_str(), std::ios::out | std::ios::app);
514
515 for(int lev(0); lev <= finestLevel(); ++lev) {
516 ParallelDescriptor::ReduceIntSum (whichPrePost[lev].dataPtr(), whichPrePost[lev].size(), IOProcNumber);
517 ParallelDescriptor::ReduceIntSum (countPrePost[lev].dataPtr(), countPrePost[lev].size(), IOProcNumber);
518 ParallelDescriptor::ReduceLongSum(wherePrePost[lev].dataPtr(), wherePrePost[lev].size(), IOProcNumber);
519
520
522 for(int j(0); j < whichPrePost[lev].size(); ++j) {
523 HdrFile << whichPrePost[lev][j] << ' ' << countPrePost[lev][j] << ' ' << wherePrePost[lev][j] << '\n';
524 }
525
526 const bool gotsome = (nParticlesAtLevelPrePost[lev] > 0);
527 if(gotsome && doUnlink) {
528// BL_PROFILE_VAR("PC<NNNN>::Checkpoint:unlink", unlink_post);
529 // Unlink any zero-length data files.
530 Vector<Long> cnt(nOutFilesPrePost,0);
531
532 for(int i(0), N = countPrePost[lev].size(); i < N; ++i) {
533 cnt[whichPrePost[lev][i]] += countPrePost[lev][i];
534 }
535
536 for(int i(0), N = cnt.size(); i < N; ++i) {
537 if(cnt[i] == 0) {
538 std::string FullFileName = NFilesIter::FileName(i, filePrefixPrePost[lev]);
539 FileSystem::Remove(FullFileName);
540 }
541 }
542 }
543 }
544 }
545
547 HdrFile.flush();
548 HdrFile.close();
549 if( ! HdrFile.good()) {
550 amrex::Abort("ParticleContainer::CheckpointPostHDF5(): problem writing HdrFile");
551 }
552 }
553}
554
555template <typename ParticleType, int NArrayReal, int NArrayInt,
556 template<class> class Allocator, class CellAssignor>
557void
558ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
559::WritePlotFilePreHDF5 ()
560{
562}
563
564
565template <typename ParticleType, int NArrayReal, int NArrayInt,
566 template<class> class Allocator, class CellAssignor>
567void
568ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
569::WritePlotFilePostHDF5 ()
570{
572}
573
574
575template <typename ParticleType, int NArrayReal, int NArrayInt,
576 template<class> class Allocator, class CellAssignor>
577void
578ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
579::WriteParticlesHDF5 (int lev, hid_t grp,
580 Vector<int>& which, Vector<int>& count, Vector<Long>& where,
581 const Vector<int>& write_real_comp,
582 const Vector<int>& write_int_comp,
583 const std::string& compression,
584 const Vector<std::map<std::pair<int, int>, IntVector>>& particle_io_flags,
585 bool is_checkpoint) const
586{
587 BL_PROFILE("ParticleContainer::WriteParticlesHDF5()");
588
589 // For a each grid, the tiles it contains
590 std::map<int, Vector<int> > tile_map;
591
592 int ret;
593 hid_t dxpl_col, dxpl_ind, dcpl_int, dcpl_real;
594 dxpl_col = H5Pcreate(H5P_DATASET_XFER);
595 dxpl_ind = H5Pcreate(H5P_DATASET_XFER);
596#ifdef AMREX_USE_MPI
597 H5Pset_dxpl_mpio(dxpl_col, H5FD_MPIO_COLLECTIVE);
598#endif
599 dcpl_int = H5Pcreate(H5P_DATASET_CREATE);
600 dcpl_real = H5Pcreate(H5P_DATASET_CREATE);
601
602 std::string mode_env, value_env;
603 double comp_value = -1.0;
604 std::string::size_type pos = compression.find('@');
605 if (pos != std::string::npos) {
606 mode_env = compression.substr(0, pos);
607 value_env = compression.substr(pos+1);
608 if (!value_env.empty()) {
609 comp_value = atof(value_env.c_str());
610 }
611 }
612
613 H5Pset_alloc_time(dcpl_int, H5D_ALLOC_TIME_LATE);
614 H5Pset_alloc_time(dcpl_real, H5D_ALLOC_TIME_LATE);
615
616 if (!mode_env.empty() && mode_env != "None") {
617 const char *chunk_env = NULL;
618 hsize_t chunk_dim = 1024;
619 chunk_env = getenv("HDF5_CHUNK_SIZE");
620 if (chunk_env != NULL) {
621 chunk_dim = atoi(chunk_env);
622 }
623
624 H5Pset_chunk(dcpl_int, 1, &chunk_dim);
625 H5Pset_chunk(dcpl_real, 1, &chunk_dim);
626
627#ifdef AMREX_USE_HDF5_ZFP
628 pos = mode_env.find("ZFP");
629 if (pos != std::string::npos) {
630 ret = H5Z_zfp_initialize();
631 if (ret < 0) { amrex::Abort("ZFP initialize failed!"); }
632 }
633#endif
634
635 if (mode_env == "ZLIB") {
636 H5Pset_shuffle(dcpl_int);
637 H5Pset_shuffle(dcpl_real);
638 H5Pset_deflate(dcpl_int, (int)comp_value);
639 H5Pset_deflate(dcpl_real, (int)comp_value);
640 }
641#ifdef AMREX_USE_HDF5_SZ
642 else if (mode_env == "SZ") {
643 ret = H5Z_SZ_Init((char*)value_env.c_str());
644 if (ret < 0) {
645 std::cout << "SZ config file:" << value_env.c_str() << std::endl;
646 amrex::Abort("SZ initialize failed, check SZ config file!");
647 }
648 }
649#endif
650#ifdef AMREX_USE_HDF5_ZFP
651 else if (mode_env == "ZFP_RATE") {
652 H5Pset_zfp_rate(dcpl_int, comp_value);
653 H5Pset_zfp_rate(dcpl_real, comp_value);
654 }
655 else if (mode_env == "ZFP_PRECISION") {
656 H5Pset_zfp_precision(dcpl_int, (unsigned int)comp_value);
657 H5Pset_zfp_precision(dcpl_real, (unsigned int)comp_value);
658 }
659 else if (mode_env == "ZFP_ACCURACY") {
660 H5Pset_zfp_accuracy(dcpl_int, comp_value);
661 H5Pset_zfp_accuracy(dcpl_real, comp_value);
662 }
663 else if (mode_env == "ZFP_REVERSIBLE") {
664 H5Pset_zfp_reversible(dcpl_int);
665 H5Pset_zfp_reversible(dcpl_real);
666 }
667#endif
668 if (ParallelDescriptor::MyProc() == 0) {
669 std::cout << "\nHDF5 particle checkpoint using " << mode_env << ", " <<
670 value_env << ", " << chunk_dim << std::endl;
671 }
672 }
673
674 for (const auto& kv : m_particles[lev])
675 {
676 const int grid = kv.first.first;
677 const int tile = kv.first.second;
678 tile_map[grid].push_back(tile);
679 const auto& pflags = particle_io_flags[lev].at(kv.first);
680
681 // Only write out valid particles.
682 count[grid] += particle_detail::countFlags(pflags);
683 }
684
685 MFInfo info;
686 info.SetAlloc(false);
687 MultiFab state(ParticleBoxArray(lev), ParticleDistributionMap(lev), 1,0,info);
688
689 int my_mfi_cnt = 0;
690 ULong my_mfi_int_total_size = 0, my_mfi_real_total_size = 0, int_size, real_size;
691 Vector<int> all_mfi_cnt(ParallelDescriptor::NProcs());
692 Vector<int> my_mfi_real_size;
693 Vector<int> my_mfi_int_size;
694 Vector<int> my_nparticles;
695 Vector<ULong> all_mfi_real_total_size(ParallelDescriptor::NProcs());
696 Vector<ULong> all_mfi_int_total_size(ParallelDescriptor::NProcs());
697 hid_t offset_id, offset_space, real_mem_space, real_dset_space, real_dset_id;
698 hid_t int_mem_space, int_dset_id, int_dset_space;
699 hsize_t total_mfi = 0, total_real_size = 0, total_int_size = 0, real_file_offset = 0, int_file_offset = 0;
700 hsize_t my_int_offset, my_int_count, my_real_offset, my_real_count;
701
702 // Count total number of components written
703 int real_comp_count = AMREX_SPACEDIM; // position values
704
705 for (int i = 0; i < NStructReal + NumRealComps(); ++i ) {
706 if (write_real_comp[i]) { ++real_comp_count; }
707 }
708
709 int int_comp_count = 2; // cpu and id values
710
711 for (int i = 0; i < NStructInt + NumIntComps(); ++i ) {
712 if (write_int_comp[i]) { ++int_comp_count; }
713 }
714
715 // Get the size for each mf so we know the amount of data from each rank
716 for (MFIter mfi(state); mfi.isValid(); ++mfi) {
717 const int grid = mfi.index();
718 if (count[grid] == 0)
719 continue;
720
721 int_size = count[grid] * int_comp_count;
722 my_mfi_int_size.push_back(int_size);
723 my_nparticles.push_back(count[grid]);
724 my_mfi_int_total_size += int_size;
725
726
727 real_size = count[grid] * real_comp_count;
728 my_mfi_real_size.push_back(real_size);
729 my_mfi_real_total_size += real_size;
730 my_mfi_cnt++;
731 }
732
733 #ifdef AMREX_USE_MPI
734 // Collect the number of mf and total size of mf from each rank
735 MPI_Allgather(&my_mfi_cnt, 1, ParallelDescriptor::Mpi_typemap<int>::type(), &(all_mfi_cnt[0]), 1,
737
738 for (int i = 0; i < ParallelDescriptor::NProcs(); i++)
739 total_mfi += all_mfi_cnt[i];
740
741 // Create the int data
742 MPI_Allgather(&my_mfi_int_total_size, 1, ParallelDescriptor::Mpi_typemap<ULong>::type(),
744 #else
745 all_mfi_cnt[0] = my_mfi_cnt;
746 all_mfi_int_total_size[0] = my_mfi_int_total_size;
747 #endif
748
749 int_file_offset = 0;
750 for (int i = 0; i < ParallelDescriptor::MyProc(); i++)
751 int_file_offset += all_mfi_int_total_size[i];
752 my_int_offset = int_file_offset;
753 my_int_count = 0;
754
755 for (int i = 0; i < ParallelDescriptor::NProcs(); i++)
756 total_int_size += all_mfi_int_total_size[i];
757
758 // SZ int compression seems to have issues at the moment
759/* #ifdef AMREX_USE_HDF5_SZ */
760/* if (mode_env == "SZ") { */
761/* size_t cd_nelmts; */
762/* unsigned int* cd_values = NULL; */
763/* unsigned filter_config; */
764/* SZ_metaDataToCdArray(&cd_nelmts, &cd_values, SZ_INT32, 0, 0, 0, 0, total_int_size); */
765/* H5Pset_filter(dcpl_int, H5Z_FILTER_SZ, H5Z_FLAG_MANDATORY, cd_nelmts, cd_values); */
766/* } */
767/* #endif */
768
769 hsize_t chunk_size;
770 if (H5Pget_layout(dcpl_int) == H5D_CHUNKED) {
771 if (H5Pget_chunk(dcpl_int, 1, &chunk_size) > -1) {
772 if (chunk_size > total_int_size) {
773 H5Pset_chunk(dcpl_int, 1, &total_int_size);
774 }
775 }
776 }
777
778 int_dset_space = H5Screate_simple(1, &total_int_size, NULL);
779#ifdef AMREX_USE_HDF5_ASYNC
780 int_dset_id = H5Dcreate_async(grp, "data:datatype=0", H5T_NATIVE_INT, int_dset_space, H5P_DEFAULT, dcpl_int, H5P_DEFAULT, es_par_g);
781#else
782 int_dset_id = H5Dcreate(grp, "data:datatype=0", H5T_NATIVE_INT, int_dset_space, H5P_DEFAULT, dcpl_int, H5P_DEFAULT);
783#endif
784
785 H5Sclose(int_dset_space);
786
787 // Create the real data
788 #ifdef AMREX_USE_MPI
789 MPI_Allgather(&my_mfi_real_total_size, 1, ParallelDescriptor::Mpi_typemap<ULong>::type(),
791 #else
792 all_mfi_real_total_size[0] = my_mfi_real_total_size;
793 #endif
794
795 for (int i = 0; i < ParallelDescriptor::NProcs(); i++)
796 total_real_size += all_mfi_real_total_size[i];
797
798#ifdef AMREX_USE_HDF5_SZ
799 if (mode_env == "SZ") {
800 size_t cd_nelmts;
801 unsigned int* cd_values = NULL;
802 unsigned filter_config;
803 if (sizeof(typename ParticleType::RealType) == 4) {
804 SZ_metaDataToCdArray(&cd_nelmts, &cd_values, SZ_FLOAT, 0, 0, 0, 0, total_real_size);
805 H5Pset_filter(dcpl_real, H5Z_FILTER_SZ, H5Z_FLAG_MANDATORY, cd_nelmts, cd_values);
806 }
807 else {
808 SZ_metaDataToCdArray(&cd_nelmts, &cd_values, SZ_DOUBLE, 0, 0, 0, 0, total_real_size);
809 H5Pset_filter(dcpl_real, H5Z_FILTER_SZ, H5Z_FLAG_MANDATORY, cd_nelmts, cd_values);
810 }
811 }
812#endif
813
814 if (H5Pget_layout(dcpl_real) == H5D_CHUNKED) {
815 if (H5Pget_chunk(dcpl_real, 1, &chunk_size) > -1) {
816 if (chunk_size > total_real_size) {
817 H5Pset_chunk(dcpl_real, 1, &total_real_size);
818 }
819 }
820 }
821
822 real_dset_space = H5Screate_simple(1, &total_real_size, NULL);
823 if (sizeof(typename ParticleType::RealType) == 4) {
824#ifdef AMREX_USE_HDF5_ASYNC
825 real_dset_id = H5Dcreate_async(grp, "data:datatype=1", H5T_NATIVE_FLOAT, real_dset_space,
826 H5P_DEFAULT, dcpl_real, H5P_DEFAULT, es_par_g);
827#else
828 real_dset_id = H5Dcreate(grp, "data:datatype=1", H5T_NATIVE_FLOAT, real_dset_space,
829 H5P_DEFAULT, dcpl_real, H5P_DEFAULT);
830#endif
831 }
832 else {
833#ifdef AMREX_USE_HDF5_ASYNC
834 real_dset_id = H5Dcreate_async(grp, "data:datatype=1", H5T_NATIVE_DOUBLE, real_dset_space,
835 H5P_DEFAULT, dcpl_real, H5P_DEFAULT, es_par_g);
836#else
837 real_dset_id = H5Dcreate(grp, "data:datatype=1", H5T_NATIVE_DOUBLE, real_dset_space,
838 H5P_DEFAULT, dcpl_real, H5P_DEFAULT);
839#endif
840 }
841 H5Sclose(real_dset_space);
842
843 real_file_offset = 0;
844 for (int i = 0; i < ParallelDescriptor::MyProc(); i++) {
845 real_file_offset += all_mfi_real_total_size[i];
846 }
847 my_real_offset = real_file_offset;
848 my_real_count = 0;
849
850 int max_mfi_count = 0, write_count = 0;
851 for (int i = 0; i < ParallelDescriptor::NProcs(); i++) {
852 if (max_mfi_count < all_mfi_cnt[i]) {
853 max_mfi_count = all_mfi_cnt[i];
854 }
855 }
856
857
858 for (MFIter mfi(state); mfi.isValid(); ++mfi)
859 {
860 const int grid = mfi.index();
861
862 if (count[grid] == 0) { continue; }
863
864 Vector<int> istuff;
865 Vector<ParticleReal> rstuff;
866 particle_detail::packIOData(istuff, rstuff, *this, lev, grid,
867 write_real_comp, write_int_comp,
868 particle_io_flags, tile_map[grid], count[grid],
869 is_checkpoint);
870
871 my_int_count = istuff.size();
872 int_mem_space = H5Screate_simple(1, &my_int_count, NULL);
873 /* std::cout << "Rank " << ParallelDescriptor::MyProc() << ": my_int_offset = " << */
874 /* my_int_offset << ", my_int_count = " << my_int_count << ", total_int_size = " << total_int_size << std::endl; */
875 int_dset_space = H5Screate_simple(1, &total_int_size, NULL);
876 H5Sselect_hyperslab (int_dset_space, H5S_SELECT_SET, &my_int_offset, NULL, &my_int_count, NULL);
877
878#ifdef AMREX_USE_HDF5_ASYNC
879 ret = H5Dwrite_async(int_dset_id, H5T_NATIVE_INT, int_mem_space, int_dset_space, dxpl_col, istuff.dataPtr(), es_par_g);
880#else
881 ret = H5Dwrite(int_dset_id, H5T_NATIVE_INT, int_mem_space, int_dset_space, dxpl_col, istuff.dataPtr());
882#endif
883 if (ret < 0) { amrex::Abort("H5Dwrite int_dset failed!"); }
884
885 H5Sclose(int_dset_space);
886 H5Sclose(int_mem_space);
887
888 my_int_offset += my_int_count;
889
890 my_real_count = rstuff.size();
891 real_mem_space = H5Screate_simple(1, &my_real_count, NULL);
892 /* std::cout << "Rank " << ParallelDescriptor::MyProc() << ": my_real_offset = " << */
893 /* my_real_offset << ", my_real_count = " << my_real_count << ", total_real_size = " << total_real_size << '\n'; */
894 real_dset_space = H5Screate_simple(1, &total_real_size, NULL);
895 H5Sselect_hyperslab (real_dset_space, H5S_SELECT_SET, &my_real_offset, NULL, &my_real_count, NULL);
896 if (sizeof(typename ParticleType::RealType) == 4) {
897#ifdef AMREX_USE_HDF5_ASYNC
898 ret = H5Dwrite_async(real_dset_id, H5T_NATIVE_FLOAT, real_mem_space, real_dset_space, dxpl_col, rstuff.dataPtr(), es_par_g);
899#else
900 ret = H5Dwrite(real_dset_id, H5T_NATIVE_FLOAT, real_mem_space, real_dset_space, dxpl_col, rstuff.dataPtr());
901#endif
902 } else {
903#ifdef AMREX_USE_HDF5_ASYNC
904 ret = H5Dwrite_async(real_dset_id, H5T_NATIVE_DOUBLE, real_mem_space, real_dset_space, dxpl_col, rstuff.dataPtr(), es_par_g);
905#else
906 ret = H5Dwrite(real_dset_id, H5T_NATIVE_DOUBLE, real_mem_space, real_dset_space, dxpl_col, rstuff.dataPtr());
907#endif
908 }
909
910 if (ret < 0) { amrex::Abort("H5Dwrite real_dset failed!"); }
911
912 H5Sclose(real_mem_space);
913 H5Sclose(real_dset_space);
914
915 my_real_offset += my_real_count;
916 write_count++;
917
918 } // end for (mfi)
919
920 // Dummy writes so that every rank participates to every possible H5Dwrite (collective)
921 while (write_count < max_mfi_count) {
922 int_dset_space = H5Screate_simple(1, &total_int_size, NULL);
923 real_dset_space = H5Screate_simple(1, &total_real_size, NULL);
924 H5Sselect_none(int_dset_space);
925 H5Sselect_none(real_dset_space);
926
927#ifdef AMREX_USE_HDF5_ASYNC
928 H5Dwrite_async(int_dset_id, H5T_NATIVE_INT, int_dset_space, int_dset_space, dxpl_col, NULL, es_par_g);
929 if (sizeof(typename ParticleType::RealType) == 4) {
930 H5Dwrite_async(real_dset_id, H5T_NATIVE_FLOAT, real_dset_space, real_dset_space, dxpl_col, NULL, es_par_g);
931 } else {
932 H5Dwrite_async(real_dset_id, H5T_NATIVE_DOUBLE, real_dset_space, real_dset_space, dxpl_col, NULL, es_par_g);
933 }
934#else
935 H5Dwrite(int_dset_id, H5T_NATIVE_INT, int_dset_space, int_dset_space, dxpl_col, NULL);
936 if (sizeof(typename ParticleType::RealType) == 4) {
937 H5Dwrite(real_dset_id, H5T_NATIVE_FLOAT, real_dset_space, real_dset_space, dxpl_col, NULL);
938 } else {
939 H5Dwrite(real_dset_id, H5T_NATIVE_DOUBLE, real_dset_space, real_dset_space, dxpl_col, NULL);
940 }
941#endif
942
943 H5Sclose(int_dset_space);
944 H5Sclose(real_dset_space);
945
946 write_count++;
947 }
948
949#ifdef AMREX_USE_HDF5_ASYNC
950 H5Dclose_async(real_dset_id, es_par_g);
951 H5Dclose_async(int_dset_id, es_par_g);
952#else
953 H5Dclose(real_dset_id);
954 H5Dclose(int_dset_id);
955#endif
956
957 // Create and write the size dataset
958 offset_space = H5Screate_simple(1, &total_mfi, NULL);
959#ifdef AMREX_USE_HDF5_ASYNC
960 offset_id = H5Dcreate_async(grp, "nparticles_grid", H5T_NATIVE_INT, offset_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT, es_par_g);
961#else
962 offset_id = H5Dcreate(grp, "nparticles_grid", H5T_NATIVE_INT, offset_space, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
963#endif
964
965 my_int_offset = 0;
966 for (int i = 0; i < ParallelDescriptor::MyProc(); i++) {
967 my_int_offset += all_mfi_cnt[i];
968 }
969 my_int_count = my_mfi_cnt;
970 int_mem_space = H5Screate_simple(1, &my_int_count, NULL);
971 /* std::cout << "Rank " << ParallelDescriptor::MyProc() << ": my_int_offset = " << */
972 /* my_int_offset << ", my_int_count = " << my_int_count << ", total_mfi = " << total_mfi << '\n'; */
973 H5Sselect_hyperslab (offset_space, H5S_SELECT_SET, &my_int_offset, NULL, &my_int_count, NULL);
974
975#ifdef AMREX_USE_HDF5_ASYNC
976 ret = H5Dwrite_async(offset_id, H5T_NATIVE_INT, int_mem_space, offset_space, dxpl_col, &(my_nparticles[0]), es_par_g);
977#else
978 ret = H5Dwrite(offset_id, H5T_NATIVE_INT, int_mem_space, offset_space, dxpl_col, &(my_nparticles[0]));
979#endif
980 if (ret < 0) { amrex::Abort("H5Dwrite offset failed!"); }
981
982 H5Pclose(dcpl_int);
983 H5Pclose(dcpl_real);
984 H5Pclose(dxpl_col);
985 H5Pclose(dxpl_ind);
986
987 H5Sclose(int_mem_space);
988 H5Sclose(offset_space);
989
990#ifdef AMREX_USE_HDF5_ASYNC
991 H5Dclose_async(offset_id, es_par_g);
992#else
993 H5Dclose(offset_id);
994#endif
995
996 /* std::cout << "Rank " << ParallelDescriptor::MyProc() << ": done WriteParticlesHDF5" << std::endl; */
997 return;
998} // End WriteParticlesHDF5
999
1000template <typename ParticleType, int NArrayReal, int NArrayInt,
1001 template<class> class Allocator, class CellAssignor>
1002void
1003ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
1004::RestartHDF5 (const std::string& dir, const std::string& file, bool /*is_checkpoint*/)
1005{
1006 RestartHDF5(dir, file);
1007}
1008
1009template <typename ParticleType, int NArrayReal, int NArrayInt,
1010 template<class> class Allocator, class CellAssignor>
1011void
1012ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
1013::RestartHDF5 (const std::string& dir, const std::string& file)
1014{
1015 BL_PROFILE("ParticleContainer::RestartHDF5()");
1016 AMREX_ASSERT(!dir.empty());
1017 AMREX_ASSERT(!file.empty());
1018
1019 const auto strttime = amrex::second();
1020
1021 std::string fullname = dir;
1022 if (!fullname.empty() && fullname[fullname.size()-1] != '/') {
1023 fullname += '/';
1024 }
1025
1026 fullname += file;
1027 fullname += ".h5";
1028
1029 hid_t fid, dset, grp, fapl, attr, atype, dspace, int_dset, real_dset;
1030 int ret;
1031
1032 fapl = H5Pcreate (H5P_FILE_ACCESS);
1033#ifdef AMREX_USE_MPI
1035#else
1036 SetHDF5fapl(fapl);
1037#endif
1038
1039 fid = H5Fopen(fullname.c_str(), H5F_ACC_RDONLY, fapl);
1040 if (fid < 0) {
1041 std::string msg("ParticleContainer::RestartHDF5(): unable to open file: ");
1042 msg += fullname;
1043 amrex::Abort(msg.c_str());
1044 }
1045
1046 attr = H5Aopen(fid, "version_name", H5P_DEFAULT);
1047 if (attr < 0) {
1048 std::string msg("ParticleContainer::RestartHDF5(): unable to open version attribute ");
1049 amrex::Abort(msg.c_str());
1050 }
1051
1052 atype = H5Aget_type(attr);
1053 if (atype < 0) {
1054 std::string msg("ParticleContainer::RestartHDF5(): unable to get type of attribute ");
1055 amrex::Abort(msg.c_str());
1056 }
1057
1058 size_t attr_len = H5Tget_size(atype);
1059 if (attr_len == 0) {
1060 std::string msg("ParticleContainer::RestartHDF5(): unable to get size of attribute ");
1061 amrex::Abort(msg.c_str());
1062 }
1063
1064 std::string version;
1065 version.resize(attr_len+1);
1066
1067 ret = H5Aread(attr, atype, &version[0]);
1068 H5Tclose(atype);
1069
1070 H5Aclose(attr);
1071
1072 AMREX_ASSERT(!version.empty());
1073
1074 // What do our version strings mean?
1075 // "Version_One_Dot_Zero" -- hard-wired to write out in double precision.
1076 // "Version_One_Dot_One" -- can write out either as either single or double precision.
1077 // Appended to the latter version string are either "_single" or "_double" to
1078 // indicate how the particles were written.
1079 // "Version_Two_Dot_Zero" -- this is the AMReX particle file format
1080 // "Version_Two_Dot_One" -- expanded particle ids to allow for 2**39-1 per proc
1081 std::string how;
1082 bool convert_ids = false;
1083 if (version.find("Version_Two_Dot_One") != std::string::npos) {
1084 convert_ids = true;
1085 }
1086 if (version.find("_single") != std::string::npos) {
1087 how = "single";
1088 }
1089 else if (version.find("_double") != std::string::npos) {
1090 how = "double";
1091 }
1092 else {
1093 std::string msg("ParticleContainer::Restart(): bad version string: ");
1094 msg += version;
1095 amrex::Error(version.c_str());
1096 }
1097
1098 std::string gname = "Chombo_global";
1099 std::string aname = "SpaceDim";
1100 int dm;
1101 grp = H5Gopen(fid, gname.c_str(), H5P_DEFAULT);
1102 if (grp < 0) {
1103 std::string msg("ParticleContainer::RestartHDF5(): unable to open group ");
1104 msg += gname;
1105 amrex::Abort(msg.c_str());
1106 }
1107 ret = ReadHDF5Attr(grp, aname.c_str(), &dm, H5T_NATIVE_INT);
1108 if (ret < 0) {
1109 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1110 msg += aname;
1111 amrex::Abort(msg.c_str());
1112 }
1113 H5Gclose(grp);
1114
1115 aname = "num_component_real";
1116 int nr;
1117 ret = ReadHDF5Attr(fid, aname.c_str(), &nr, H5T_NATIVE_INT);
1118 if (ret < 0) {
1119 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1120 msg += aname;
1121 amrex::Abort(msg.c_str());
1122 }
1123 if (nr != NStructReal + NumRealComps())
1124 amrex::Abort("ParticleContainer::Restart(): nr != NStructReal + NumRealComps()");
1125
1126 aname = "num_component_int";
1127 int ni;
1128 ret = ReadHDF5Attr(fid, aname.c_str(), &ni, H5T_NATIVE_INT);
1129 if (ret < 0) {
1130 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1131 msg += aname;
1132 amrex::Abort(msg.c_str());
1133 }
1134 if (ni != NStructInt + NumIntComps()) {
1135 amrex::Abort("ParticleContainer::Restart(): ni != NStructInt");
1136 }
1137
1138 aname = "nparticles";
1139 Long nparticles;
1140 ret = ReadHDF5Attr(fid, aname.c_str(), &nparticles, H5T_NATIVE_LLONG);
1141 if (ret < 0) {
1142 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1143 msg += aname;
1144 amrex::Abort(msg.c_str());
1145 }
1146 AMREX_ASSERT(nparticles >= 0);
1147
1148 aname = "maxnextid";
1149 Long maxnextid;
1150 ret = ReadHDF5Attr(fid, aname.c_str(), &maxnextid, H5T_NATIVE_LLONG);
1151 if (ret < 0) {
1152 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1153 msg += aname;
1154 amrex::Abort(msg.c_str());
1155 }
1156 AMREX_ASSERT(maxnextid > 0);
1157 ParticleType::NextID(maxnextid);
1158
1159 aname = "finest_level";
1160 int finest_level_in_file;
1161 ret = ReadHDF5Attr(fid, aname.c_str(), &finest_level_in_file, H5T_NATIVE_INT);
1162 if (ret < 0) {
1163 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1164 msg += aname;
1165 amrex::Abort(msg.c_str());
1166 }
1167 AMREX_ASSERT(finest_level_in_file >= 0);
1168
1169 hid_t comp_dtype = H5Tcreate (H5T_COMPOUND, 2 * AMREX_SPACEDIM * sizeof(int));
1170 if (1 == AMREX_SPACEDIM) {
1171 H5Tinsert (comp_dtype, "lo_i", 0 * sizeof(int), H5T_NATIVE_INT);
1172 H5Tinsert (comp_dtype, "hi_i", 1 * sizeof(int), H5T_NATIVE_INT);
1173 }
1174 else if (2 == AMREX_SPACEDIM) {
1175 H5Tinsert (comp_dtype, "lo_i", 0 * sizeof(int), H5T_NATIVE_INT);
1176 H5Tinsert (comp_dtype, "lo_j", 1 * sizeof(int), H5T_NATIVE_INT);
1177 H5Tinsert (comp_dtype, "hi_i", 2 * sizeof(int), H5T_NATIVE_INT);
1178 H5Tinsert (comp_dtype, "hi_j", 3 * sizeof(int), H5T_NATIVE_INT);
1179 }
1180 else if (3 == AMREX_SPACEDIM) {
1181 H5Tinsert (comp_dtype, "lo_i", 0 * sizeof(int), H5T_NATIVE_INT);
1182 H5Tinsert (comp_dtype, "lo_j", 1 * sizeof(int), H5T_NATIVE_INT);
1183 H5Tinsert (comp_dtype, "lo_k", 2 * sizeof(int), H5T_NATIVE_INT);
1184 H5Tinsert (comp_dtype, "hi_i", 3 * sizeof(int), H5T_NATIVE_INT);
1185 H5Tinsert (comp_dtype, "hi_j", 4 * sizeof(int), H5T_NATIVE_INT);
1186 H5Tinsert (comp_dtype, "hi_k", 5 * sizeof(int), H5T_NATIVE_INT);
1187 }
1188
1189 // Determine whether this is a dual-grid restart or not.
1190 Vector<BoxArray> particle_box_arrays(finest_level_in_file + 1);
1191 bool dual_grid = false;
1192
1193 for (int lev = 0; lev <= finest_level_in_file; lev++)
1194 {
1195 if (lev > finestLevel())
1196 {
1197 dual_grid = true;
1198 break;
1199 }
1200
1201 gname = "level_" + std::to_string(lev);
1202 grp = H5Gopen(fid, gname.c_str(), H5P_DEFAULT);
1203 if (grp < 0) {
1204 std::string msg("ParticleContainer::RestartHDF5(): unable to open group ");
1205 msg += gname;
1206 amrex::Abort(msg.c_str());
1207 }
1208
1209 dset = H5Dopen(grp, "boxes", H5P_DEFAULT);
1210 dspace = H5Dget_space(dset);
1211 hsize_t ngrid;
1212 H5Sget_simple_extent_dims(dspace, &ngrid, NULL);
1213
1214 Vector<int> boxes(ngrid*AMREX_SPACEDIM*2);
1215 ret = H5Dread(dset, comp_dtype, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(boxes[0]));
1216 if (ret < 0) {
1217 std::string msg("ParticleContainer::RestartHDF5(): unable to read nparticles_grid dataset");
1218 amrex::Abort(msg.c_str());
1219 }
1220 H5Sclose(dspace);
1221 H5Dclose(dset);
1222 H5Gclose(grp);
1223
1224 /* particle_box_arrays[lev].readFrom(phdr_file); */
1225 particle_box_arrays[lev] = ParticleBoxArray(lev);
1226 for (int i = 0; i < (int)ngrid; i++) {
1227 Box tmp (IntVect{AMREX_D_DECL(boxes[i*AMREX_SPACEDIM], boxes[i*AMREX_SPACEDIM+1], boxes[i*AMREX_SPACEDIM+2])},
1228 IntVect{AMREX_D_DECL(boxes[i*AMREX_SPACEDIM*2], boxes[i*AMREX_SPACEDIM*2+1], boxes[i*AMREX_SPACEDIM*2+2])});
1229 particle_box_arrays[lev][i] = tmp;
1230 }
1231
1232 if (! particle_box_arrays[lev].CellEqual(ParticleBoxArray(lev)))
1233 dual_grid = true;
1234 }
1235
1236 if (dual_grid) {
1237 for (int lev = 0; lev <= finestLevel(); lev++) {
1238 SetParticleBoxArray(lev, particle_box_arrays[lev]);
1239 DistributionMapping pdm(particle_box_arrays[lev]);
1240 SetParticleDistributionMap(lev, pdm);
1241 }
1242 }
1243
1244 Vector<int> ngrids(finest_level_in_file+1);
1245 for (int lev = 0; lev <= finest_level_in_file; lev++) {
1246 gname = "level_" + std::to_string(lev);
1247 grp = H5Gopen(fid, gname.c_str(), H5P_DEFAULT);
1248 if (grp < 0) {
1249 std::string msg("ParticleContainer::RestartHDF5(): unable to open group ");
1250 msg += gname;
1251 amrex::Abort(msg.c_str());
1252 }
1253
1254 aname = "ngrids";
1255 ret = ReadHDF5Attr(grp, aname.c_str(), &ngrids[lev], H5T_NATIVE_INT);
1256 if (ret < 0) {
1257 std::string msg("ParticleContainer::RestartHDF5(): unable to read attribute ");
1258 msg += aname;
1259 amrex::Abort(msg.c_str());
1260 }
1261
1262 AMREX_ASSERT(ngrids[lev] > 0);
1263 if (lev <= finestLevel()) {
1264 AMREX_ASSERT(ngrids[lev] == int(ParticleBoxArray(lev).size()));
1265 }
1266
1267 H5Gclose(grp);
1268 }
1269
1270 resizeData();
1271
1272 if (finest_level_in_file > finestLevel()) {
1273 m_particles.resize(finest_level_in_file+1);
1274 }
1275
1276 for (int lev = 0; lev <= finest_level_in_file; lev++) {
1277
1278 gname = "level_" + std::to_string(lev);
1279 grp = H5Gopen(fid, gname.c_str(), H5P_DEFAULT);
1280 if (grp < 0) {
1281 std::string msg("ParticleContainer::RestartHDF5(): unable to open group ");
1282 msg += gname;
1283 amrex::Abort(msg.c_str());
1284 }
1285
1286 dset = H5Dopen(grp, "nparticles_grid", H5P_DEFAULT);
1287 Vector<int> count(ngrids[lev]);
1288 ret = H5Dread(dset, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(count[0]));
1289 if (ret < 0) {
1290 std::string msg("ParticleContainer::RestartHDF5(): unable to read nparticles_grid dataset");
1291 amrex::Abort(msg.c_str());
1292 }
1293 H5Dclose(dset);
1294
1295 Vector<hsize_t> offset(ngrids[lev]);
1296 offset[0] = 0;
1297 for (int i = 1; i < ngrids[lev]; i++) {
1298 offset[i] = offset[i-1] + count[i-1];
1299 }
1300
1301 int_dset = H5Dopen(grp, "data:datatype=0", H5P_DEFAULT);
1302 if (int_dset < 0) {
1303 std::string msg("ParticleContainer::RestartHDF5(): unable to open int dataset");
1304 amrex::Abort(msg.c_str());
1305 }
1306 real_dset = H5Dopen(grp, "data:datatype=1", H5P_DEFAULT);
1307 if (real_dset < 0) {
1308 std::string msg("ParticleContainer::RestartHDF5(): unable to open int dataset");
1309 amrex::Abort(msg.c_str());
1310 }
1311
1312 Vector<int> grids_to_read;
1313 if (lev <= finestLevel()) {
1314 for (MFIter mfi(*m_dummy_mf[lev]); mfi.isValid(); ++mfi) {
1315 grids_to_read.push_back(mfi.index());
1316 }
1317 } else {
1318
1319 // we lost a level on restart. we still need to read in particles
1320 // on finer levels, and put them in the right place via Redistribute()
1321
1322 const int rank = ParallelDescriptor::MyProc();
1323 const int NReaders = MaxReaders();
1324 if (rank >= NReaders) { return; }
1325
1326 const int Navg = ngrids[lev] / NReaders;
1327 const int Nleft = ngrids[lev] - Navg * NReaders;
1328
1329 int lo, hi;
1330 if (rank < Nleft) {
1331 lo = rank*(Navg + 1);
1332 hi = lo + Navg + 1;
1333 }
1334 else {
1335 lo = rank * Navg + Nleft;
1336 hi = lo + Navg;
1337 }
1338
1339 for (int i = lo; i < hi; ++i) {
1340 grids_to_read.push_back(i);
1341 }
1342 }
1343
1344 for(int igrid = 0; igrid < static_cast<int>(grids_to_read.size()); ++igrid) {
1345 const int grid = grids_to_read[igrid];
1346
1347 if (how == "single") {
1348 ReadParticlesHDF5<float>(offset[grid], count[grid], grid, lev, int_dset, real_dset, finest_level_in_file, convert_ids);
1349 }
1350 else if (how == "double") {
1351 ReadParticlesHDF5<double>(offset[grid], count[grid], grid, lev, int_dset, real_dset, finest_level_in_file, convert_ids);
1352 }
1353 else {
1354 std::string msg("ParticleContainer::Restart(): bad parameter: ");
1355 msg += how;
1356 amrex::Error(msg.c_str());
1357 }
1358 }
1359
1360 H5Dclose(int_dset);
1361 H5Dclose(real_dset);
1362 H5Gclose(grp);
1363 } // end for level
1364
1365 H5Fclose(fid);
1366
1367 Redistribute();
1368
1369 AMREX_ASSERT(OK());
1370
1371 if (m_verbose > 1) {
1372 auto stoptime = amrex::second() - strttime;
1374 amrex::Print() << "ParticleContainer::Restart() time: " << stoptime << '\n';
1375 }
1376}
1377
1378// Read a batch of particles from the checkpoint file
1379template <typename ParticleType, int NArrayReal, int NArrayInt,
1380 template<class> class Allocator, class CellAssignor>
1381template <class RTYPE>
1382void
1383ParticleContainer_impl<ParticleType, NArrayReal, NArrayInt, Allocator, CellAssignor>
1384::ReadParticlesHDF5 (hsize_t offset, hsize_t cnt, int grd, int lev,
1385 hid_t int_dset, hid_t real_dset, int finest_level_in_file,
1386 bool convert_ids)
1387{
1388 BL_PROFILE("ParticleContainer::ReadParticlesHDF5()");
1389 AMREX_ASSERT(cnt > 0);
1390 AMREX_ASSERT(lev < int(m_particles.size()));
1391
1392 hid_t int_dspace, int_fspace, real_dspace, real_fspace;
1393
1394 // First read in the integer data in binary. We do not store
1395 // the m_lev and m_grid data on disk. We can easily recreate
1396 // that given the structure of the checkpoint file.
1397 const int iChunkSize = 2 + NStructInt + NumIntComps();
1398 Vector<int> istuff(cnt*iChunkSize);
1399 int_fspace = H5Dget_space(int_dset);
1400 hsize_t int_cnt = cnt*iChunkSize;
1401 hsize_t int_offset = offset*iChunkSize;
1402 int_dspace = H5Screate_simple(1, &int_cnt, NULL);
1403 H5Sselect_hyperslab (int_fspace, H5S_SELECT_SET, &int_offset, NULL, &int_cnt, NULL);
1404 H5Dread(int_dset, H5T_NATIVE_INT, int_dspace, int_fspace, H5P_DEFAULT, istuff.dataPtr());
1405
1406 H5Sclose(int_fspace);
1407 H5Sclose(int_dspace);
1408
1409 // Then the real data in binary.
1410 const int rChunkSize = AMREX_SPACEDIM + NStructReal + NumRealComps();
1411 Vector<RTYPE> rstuff(cnt*rChunkSize);
1412 real_fspace = H5Dget_space(real_dset);
1413 hsize_t real_cnt = cnt*rChunkSize;
1414 hsize_t real_offset = offset*rChunkSize;
1415 real_dspace = H5Screate_simple(1, &real_cnt, NULL);
1416 H5Sselect_hyperslab (real_fspace, H5S_SELECT_SET, &real_offset, NULL, &real_cnt, NULL);
1417 if (sizeof(RTYPE) == 4) {
1418 H5Dread(real_dset, H5T_NATIVE_FLOAT, real_dspace, real_fspace, H5P_DEFAULT, rstuff.dataPtr());
1419 } else {
1420 H5Dread(real_dset, H5T_NATIVE_DOUBLE, real_dspace, real_fspace, H5P_DEFAULT, rstuff.dataPtr());
1421 }
1422
1423 H5Sclose(real_fspace);
1424 H5Sclose(real_dspace);
1425
1426 // Now reassemble the particles.
1427 int* iptr = istuff.dataPtr();
1428 RTYPE* rptr = rstuff.dataPtr();
1429
1430 ParticleType p;
1431 ParticleLocData pld;
1432
1433 Vector<std::map<std::pair<int, int>, Gpu::HostVector<ParticleType> > > host_particles;
1434 host_particles.reserve(15);
1435 host_particles.resize(finest_level_in_file+1);
1436
1437 Vector<std::map<std::pair<int, int>,
1438 std::vector<Gpu::HostVector<ParticleReal> > > > host_real_attribs;
1439 host_real_attribs.reserve(15);
1440 host_real_attribs.resize(finest_level_in_file+1);
1441
1442 Vector<std::map<std::pair<int, int>,
1443 std::vector<Gpu::HostVector<int> > > > host_int_attribs;
1444 host_int_attribs.reserve(15);
1445 host_int_attribs.resize(finest_level_in_file+1);
1446
1447 for (hsize_t i = 0; i < cnt; i++) {
1448 if (convert_ids) {
1449 std::int32_t xi, yi;
1450 std::uint32_t xu, yu;
1451 xi = iptr[0];
1452 yi = iptr[1];
1453 std::memcpy(&xu, &xi, sizeof(xi));
1454 std::memcpy(&yu, &yi, sizeof(yi));
1455 p.m_idcpu = ((std::uint64_t)xu) << 32 | yu;
1456 } else {
1457 p.id() = iptr[0];
1458 p.cpu() = iptr[1];
1459 }
1460 iptr += 2;
1461
1462 for (int j = 0; j < NStructInt; j++)
1463 {
1464 p.idata(j) = *iptr;
1465 ++iptr;
1466 }
1467
1468 AMREX_ASSERT(p.id() > 0);
1469
1470 AMREX_D_TERM(p.pos(0) = ParticleReal(rptr[0]);,
1471 p.pos(1) = ParticleReal(rptr[1]);,
1472 p.pos(2) = ParticleReal(rptr[2]););
1473
1474 rptr += AMREX_SPACEDIM;
1475
1476 for (int j = 0; j < NStructReal; j++)
1477 {
1478 p.rdata(j) = ParticleReal(*rptr);
1479 ++rptr;
1480 }
1481
1482 locateParticle(p, pld, 0, finestLevel(), 0);
1483
1484 std::pair<int, int> ind(grd, pld.m_tile);
1485
1486 host_real_attribs[lev][ind].resize(NumRealComps());
1487 host_int_attribs[lev][ind].resize(NumIntComps());
1488
1489 // add the struct
1490 host_particles[lev][ind].push_back(p);
1491
1492 // add the real...
1493 for (int icomp = 0; icomp < NumRealComps(); icomp++) {
1494 host_real_attribs[lev][ind][icomp].push_back(ParticleReal(*rptr));
1495 ++rptr;
1496 }
1497
1498 // ... and int array data
1499 for (int icomp = 0; icomp < NumIntComps(); icomp++) {
1500 host_int_attribs[lev][ind][icomp].push_back(*iptr);
1501 ++iptr;
1502 }
1503 }
1504
1505 for (int host_lev = 0; host_lev < static_cast<int>(host_particles.size()); ++host_lev)
1506 {
1507 for (auto& kv : host_particles[host_lev]) {
1508 auto grid = kv.first.first;
1509 auto tile = kv.first.second;
1510 const auto& src_tile = kv.second;
1511
1512 auto& dst_tile = DefineAndReturnParticleTile(host_lev, grid, tile);
1513 auto old_size = dst_tile.GetArrayOfStructs().size();
1514 auto new_size = old_size + src_tile.size();
1515 dst_tile.resize(new_size);
1516
1517 Gpu::copyAsync(Gpu::hostToDevice, src_tile.begin(), src_tile.end(),
1518 dst_tile.GetArrayOfStructs().begin() + old_size);
1519
1520 for (int i = 0; i < NumRealComps(); ++i) {
1522 host_real_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
1523 host_real_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
1524 dst_tile.GetStructOfArrays().GetRealData(i).begin() + old_size);
1525 }
1526
1527 for (int i = 0; i < NumIntComps(); ++i) {
1529 host_int_attribs[host_lev][std::make_pair(grid,tile)][i].begin(),
1530 host_int_attribs[host_lev][std::make_pair(grid,tile)][i].end(),
1531 dst_tile.GetStructOfArrays().GetIntData(i).begin() + old_size);
1532 }
1533 }
1534 }
1535
1537}
1538
1539}
1540
1541#endif /*AMREX_USE_HDF5*/
1542#endif /*AMREX_PARTICLEHDF5_H*/
#define BL_PROFILE(a)
Definition AMReX_BLProfiler.H:551
#define AMREX_ASSERT(EX)
Definition AMReX_BLassert.H:38
#define AMREX_GPU_HOST_DEVICE
Definition AMReX_GpuQualifiers.H:20
Array4< int const > offset
Definition AMReX_HypreMLABecLap.cpp:1089
void CheckpointPreHDF5()
void RestartHDF5(const std::string &dir, const std::string &file)
Restart from checkpoint.
void WriteHDF5ParticleData(const std::string &dir, const std::string &name, const Vector< int > &write_real_comp, const Vector< int > &write_int_comp, const Vector< std::string > &real_comp_names, const Vector< std::string > &int_comp_names, const std::string &compression, F &&f, bool is_checkpoint=false) const
Writes particle data to disk in the HDF5 format.
void CheckpointPostHDF5()
#define AMREX_D_TERM(a, b, c)
Definition AMReX_SPACE.H:129
const std::string & FileName() const
Definition AMReX_NFiles.H:160
This class provides the user with a few print options.
Definition AMReX_Print.H:35
bool Remove(std::string const &filename)
Definition AMReX_FileSystem.cpp:190
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
MPI_Comm Communicator() noexcept
Definition AMReX_ParallelDescriptor.H:210
void ReduceIntSum(int &)
Integer sum reduction.
Definition AMReX_ParallelDescriptor.cpp:1255
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
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
std::enable_if_t< RunOnGpu< typename PC::template AllocatorType< int > >::value, void > packIOData(Vector< int > &idata, Vector< ParticleReal > &rdata, const PC &pc, int lev, int grid, const Vector< int > &write_real_comp, const Vector< int > &write_int_comp, const Vector< std::map< std::pair< int, int >, typename PC::IntVector > > &particle_io_flags, const Vector< int > &tiles, int np, bool is_checkpoint)
Definition AMReX_WriteBinaryParticleData.H:173
std::enable_if_t< RunOnGpu< typename Container< int, Allocator >::allocator_type >::value, amrex::Long > countFlags(const Vector< std::map< std::pair< int, int >, Container< int, Allocator > > > &particle_io_flags, const PC &pc)
Definition AMReX_WriteBinaryParticleData.H:78
Definition AMReX_Amr.cpp:49
BoxND< AMREX_SPACEDIM > Box
Definition AMReX_BaseFwd.H:27
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
IntVectND< AMREX_SPACEDIM > IntVect
Definition AMReX_BaseFwd.H:30
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
static void SetHDF5fapl(hid_t fapl, MPI_Comm comm)
Definition AMReX_PlotFileUtilHDF5.cpp:111
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
void WriteHDF5ParticleDataSync(PC const &pc, const std::string &dir, const std::string &name, const Vector< int > &write_real_comp, const Vector< int > &write_int_comp, const Vector< std::string > &real_comp_names, const Vector< std::string > &int_comp_names, const std::string &compression, F &&f, bool is_checkpoint)
Definition AMReX_WriteBinaryParticleDataHDF5.H:117
static int ReadHDF5Attr(hid_t loc, const char *name, void *data, hid_t dtype)
Definition AMReX_WriteBinaryParticleDataHDF5.H:61