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