Chombo + EB + MF  3.2
ParticleIOI.H
Go to the documentation of this file.
1 #ifdef CH_LANG_CC
2 /*
3  * _______ __
4  * / ___/ / ___ __ _ / / ___
5  * / /__/ _ \/ _ \/ V \/ _ \/ _ \
6  * \___/_//_/\___/_/_/_/_.__/\___/
7  * Please refer to Copyright.txt, in Chombo's root directory.
8  */
9 #endif
10 
11 // functions for I/O of particle data
12 
13 #ifndef _PARTICLEIOI_H_
14 #define _PARTICLEIOI_H_
15 
16 #include "ParticleIO.H"
17 
18 #include "NamespaceHeader.H"
19 
20 // save data in small chunks
21 #define _CHUNK (1024*1024)
22 
23 #ifdef CH_USE_HDF5
24 
25 template <class T>
27  const vector<T>& a_vect,
28  const hid_t& H5T_type,
29  const std::string& a_dataname)
30 {
31  // store vect info
32  hsize_t length = a_vect.size();
33  hid_t dataspace= H5Screate_simple(1, &length, NULL);
34 
35 #ifdef H516
36  hid_t dataset = H5Dcreate(a_handle.groupID(), a_dataname.c_str(),
37  H5T_type, dataspace,
38  H5P_DEFAULT);
39 #else
40  hid_t dataset = H5Dcreate2(a_handle.groupID(), a_dataname.c_str(),
41  H5T_type, dataspace,
42  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
43 #endif
44 
45  CH_assert(dataspace >= 0);
46  CH_assert(dataset >= 0);
47  if (procID() == 0)
48  {
49  hid_t memdataspace = H5Screate_simple(1, &length, NULL);
50  CH_assert(memdataspace >= 0);
51  int err = H5Dwrite(dataset, H5T_type, memdataspace,
52  dataspace, H5P_DEFAULT, &a_vect[0]);
53  CH_assert(err >= 0);
54  H5Sclose(memdataspace);
55  }
56  H5Sclose(dataspace);
57  H5Dclose(dataset);
58 }
59 
60 // read vect to header
61 template <class T>
63  vector<T>& a_vect,
64  const hid_t& H5T_type,
65  const std::string& a_dataname)
66 {
67  // read in offset info
68  // hsize_t length = sizeof(T) * a_vect.size();
69  hsize_t length = a_vect.size();
70 
71  hid_t dataspace= H5Screate_simple(1, &length, NULL);
72 #ifdef H516
73  hid_t dataset = H5Dopen(a_handle.groupID(), a_dataname.c_str());
74 #else
75  hid_t dataset = H5Dopen2(a_handle.groupID(), a_dataname.c_str(), H5P_DEFAULT);
76 #endif
77  CH_assert(dataspace >= 0);
78  CH_assert(dataset >= 0);
79 
80  hid_t memdataspace = H5Screate_simple(1, &length, NULL);
81  CH_assert(memdataspace >= 0);
82  int err = H5Dread(dataset, H5T_type, memdataspace,
83  dataspace, H5P_DEFAULT, &a_vect[0]);
84  CH_assert(err >= 0);
85 
86  H5Sclose(memdataspace);
87  H5Sclose(dataspace);
88  H5Dclose(dataset);
89 }
90 
91 // write whole hierarchy of particle data
92 template<class P> void
94  const ParticleData<P>& a_particles,
95  const std::string& a_dataType)
96 {
97  // timer
98  CH_TIMERS("writeParticlesToHDF");
99 
100  // grids
101  const BoxLayout& grids = a_particles.getBoxes();
102 
103  // loc part per box
104  vector<unsigned long long> locParticlesPerBox(grids.size(),0);
105 
106  // loc and tot number of particles
107  unsigned long long numLocalParticles = 0;
108 
109  // count number of particles per box
110  for (DataIterator di(grids); di.ok(); ++di)
111  {
112  const size_t numItems=a_particles[di].numItems();
113  numLocalParticles += (unsigned long long)numItems;
114  locParticlesPerBox[grids.index(di())] = (unsigned long long)numItems;
115  }
116 
117  vector<unsigned long long> particlesPerBox(grids.size());
118 
119 #ifdef CH_MPI
120  int result = MPI_Allreduce(&locParticlesPerBox[0], &particlesPerBox[0], locParticlesPerBox.size(),
121  MPI_UNSIGNED_LONG_LONG, MPI_SUM, Chombo_MPI::comm);
122  if (result != MPI_SUCCESS)
123  {
124  MayDay::Error("MPI communcation error in ParticleIO");
125  }
126 
127 #else
128 
129  for (int i=0; i<grids.size(); i++)
130  {
131  particlesPerBox[i]=(unsigned long long)locParticlesPerBox[i];
132  }
133 
134 #endif
135 
136  write_hdf_part_header(a_handle, grids, particlesPerBox, a_dataType);
137 
138  // size of individual objects
139  size_t objSize = P().size();
140 
141  // the number of components per particle
142  size_t numComps = objSize / sizeof(Real);
143 
144  // compute offsets and tot number of particles; order matters
145  vector<unsigned long long> offsets;
146  unsigned long long totNumParticles = 0;
147  for (int i=0; i<grids.size(); i++)
148  {
149  offsets.push_back(numComps * totNumParticles);
150  totNumParticles += particlesPerBox[i];
151  }
152 
153  // now store particle info
154  if (totNumParticles==0)
155  {
156  return;
157  }
158 
159  CH_TIMER("writeParticlesToHDF::checkpoint",t_cp);
160  CH_START(t_cp);
161 
162  // data size
163  hsize_t dataSize = totNumParticles * numComps;
164 
165  // a single dataspace
166  hid_t dataspace = H5Screate_simple(1, &dataSize, NULL);
167 
168  // create data type
169  hid_t H5T_type = H5T_NATIVE_DOUBLE;
170 
171  //
172  std::string dataname = a_dataType+":data";
173 
174  // open dataset
175 #ifdef H516
176  hid_t dataset = H5Dcreate(a_handle.groupID(), dataname.c_str(),
177  H5T_type, dataspace,
178  H5P_DEFAULT);
179 #else
180  hid_t dataset = H5Dcreate2(a_handle.groupID(), dataname.c_str(),
181  H5T_type, dataspace,
182  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
183 #endif
184 
185  if (numLocalParticles > 0)
186  {
187  // allocate data buffer where to linearize the particle data
188  const size_t chunkSize = objSize*_CHUNK;
189  char* chunk = new char[chunkSize];
190 
191  if (chunk == NULL) {
192  MayDay::Error("WritePart::Error: new returned NULL pointer ");
193  }
194 
195  for (DataIterator di(grids); di.ok(); ++di)
196  {
197  size_t offset = offsets[ grids.index(di()) ];
198 
199  // particle counter
200  int ip = 0;
201 
202  // linearize particle data
203  const List<P>& pList = a_particles[di].listItems();
204 
205  for (ListIterator<P> li(pList); li.ok(); ++li)
206  {
207  pList[li].linearOut((void*)chunk);
208  chunk += objSize;
209 
210  if (++ip % _CHUNK == 0)
211  {
212  // rewind pointer position, write buffered data
213  chunk -= chunkSize;
214  writeDataChunk(offset,dataspace,dataset,H5T_type,
215  chunkSize / sizeof(Real),chunk);
216  }
217  }
218 
219  // write our residual particles
220  const size_t nResidual = pList.length()%_CHUNK;
221  if (nResidual>0)
222  {
223  // rewind pointer position and write buffered data
224  chunk -= (nResidual*objSize);
225  writeDataChunk(offset,dataspace,dataset,H5T_type,
226  nResidual*objSize / sizeof(Real),chunk);
227  }
228  }
229 
230  // free buffer
231  delete[] chunk; chunk=NULL;
232  }
233 
234  // done
235  H5Sclose(dataspace);
236  H5Dclose(dataset);
237 
238  // stop timer
239  CH_STOP(t_cp);
240 }
241 
242 // read whole set of particle data
243 template<class P>
245  ParticleData<P>& a_particles,
246  const std::string& a_dataType)
247 {
248  // timer
249  CH_TIMERS("readParticlesFromHDF");
250 
251  // number of particles on each box
252  vector<unsigned long long> particlesPerBox;
253 
254  // grids
255  Vector<Box> grids;
256  read_hdf_part_header(a_handle,
257  grids,
258  particlesPerBox,
259  a_dataType,
260  a_handle.getGroup());
261 
262  // load balancing should be handled outside
263  const BoxLayout& bl = a_particles.getBoxes();
264 
265  // size of individual objects
266  size_t objSize = P().size();
267 
268  // the number of components per particle
269  size_t numComps = objSize / sizeof(Real);
270 
271  // compute offsets and tot number of particles; order matters
272  vector<unsigned long long> offsets;
273  unsigned long long totNumParticles = 0;
274  for (int i=0; i<grids.size(); i++)
275  {
276  offsets.push_back(numComps * totNumParticles);
277  totNumParticles += particlesPerBox[i];
278  }
279 
280  // now read in particle data
281  if (totNumParticles==0)
282  {
283  return;
284  }
285 
286  CH_TIMER("readParticlesFromHDF::checkpoint",t_cp);
287  CH_START(t_cp);
288 
289  // P object
290  P p;
291 
292  // data type
293  hid_t H5T_type = H5T_NATIVE_DOUBLE;
294 
295  //
296  std::string dataname = a_dataType+":data";
297 
298  // one single dataset
299 #ifdef H516
300  hid_t dataset = H5Dopen(a_handle.groupID(), dataname.c_str());
301 #else
302  hid_t dataset = H5Dopen2(a_handle.groupID(), dataname.c_str(), H5P_DEFAULT);
303 #endif
304 
305  hid_t dataspace= H5Dget_space(dataset);
306 
307  // read in data relative to this proc in small chuncks
308  const size_t chunkSize = objSize*_CHUNK;
309 
310  // allocate data buffer to read in linearized particle data chunks
311  char* chunk = new char[chunkSize];
312 
313  for (DataIterator di(bl); di.ok(); ++di)
314  {
315  size_t boxData= numComps * particlesPerBox[ bl.index(di()) ];
316  size_t offset = offsets[ bl.index(di()) ];
317 
318  List<P>& items= a_particles[di].listItems();
319 
320  size_t dataIn = 0;
321  while (dataIn < boxData)
322  {
323  // read data chunk of the right size
324  int size = ((boxData-dataIn) >= (chunkSize / sizeof(Real))) ? chunkSize / sizeof(Real) : boxData - dataIn;
325  readDataChunk(offset, dataspace, dataset, H5T_type, size, chunk);
326 
327  // list to store the particle data from the HDF5 file
328  for (int ip=0; ip<size/numComps; ip++)
329  {
330  p.linearIn(chunk);
331  chunk += objSize;
332 
333  // assign to data holder
334  items.add(p);
335  }
336  dataIn += size;
337  chunk -= size*sizeof(Real);
338  }
339  }
340 
341  // free alllocated memory
342  delete[] chunk; chunk=NULL;
343 
344  // done with this component
345  H5Sclose(dataspace);
346  H5Dclose(dataset);
347 
348  // stop timer
349  CH_STOP(t_cp);
350 }
351 
352 // write whole hierarchy of particle data
353 template<class P> void
355  const Vector<P>& a_particles,
356  const Box& a_domain,
357  const std::string& a_dataType)
358 {
359  // timer
360  CH_TIMERS("writeParticlesToHDF");
361 
362  // We make a fake BoxLayout with only one Box.
363  Vector<Box> boxes;
364  boxes.push_back(a_domain);
365  Vector<int> procIDs;
366  procIDs.push_back(0);
367  BoxLayout grids(boxes, procIDs);
368 
369  unsigned long long totNumParticles = a_particles.size();
370 
371  // In this case, there is only one box
372  vector<unsigned long long> particlesPerBox;
373  particlesPerBox.push_back(totNumParticles);
374 
375  write_hdf_part_header(a_handle, grids, particlesPerBox, a_dataType);
376 
377  // size of individual objects
378  size_t objSize = P().size();
379 
380  // the number of components per particle
381  size_t numComps = objSize / sizeof(Real);
382 
383  // don't do anything if there are no particles
384  if (totNumParticles==0)
385  {
386  return;
387  }
388 
389  // data size
390  hsize_t dataSize = totNumParticles * numComps;
391 
392  // a single dataspace
393  hid_t dataspace = H5Screate_simple(1, &dataSize, NULL);
394 
395  // create data type
396  hid_t H5T_type = H5T_NATIVE_DOUBLE;
397 
398  //
399  std::string dataname = a_dataType+":data";
400 
401  // open dataset
402 
403 #ifdef H516
404  hid_t dataset = H5Dcreate(a_handle.groupID(), dataname.c_str(),
405  H5T_type, dataspace,
406  H5P_DEFAULT);
407 #else
408  hid_t dataset = H5Dcreate2(a_handle.groupID(), dataname.c_str(),
409  H5T_type, dataspace,
410  H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
411 #endif
412 
413  // allocate data buffer in which to linearize the particle data
414  const size_t bufferSize = objSize*totNumParticles;
415  char* buffer = new char[bufferSize];
416 
417  // are we good to go?
418  if (buffer == NULL) {
419  MayDay::Error("WritePart::Error: new returned NULL pointer ");
420  }
421 
422  // remember, only one box
423  size_t offset = 0;
424 
425  // write out the particles to the buffer
426  for (int ip = 0; ip < totNumParticles; ip++)
427  {
428  a_particles[ip].linearOut((void*)buffer);
429  buffer += objSize;
430  }
431 
432  // rewind the buffer pointer
433  buffer -= bufferSize;
434 
435  // write the buffer to the HDF5
436  writeDataChunk(offset, dataspace, dataset, H5T_type,
437  bufferSize / sizeof(Real), buffer);
438 
439  // free buffer
440  delete[] buffer;
441  buffer=NULL;
442 
443  // done
444  H5Sclose(dataspace);
445  H5Dclose(dataset);
446 }
447 
448 #endif // HDF5
449 
450 #include "NamespaceFooter.H"
451 
452 #endif
#define CH_TIMERS(name)
Definition: CH_Timer.H:133
void writeDataChunk(size_t &offset, const hid_t &dataspace, const hid_t &dataset, const hid_t &H5T_type, const unsigned long dataLength, const void *const data)
Write chunk of data and upgrade offset. Not meant to be a user-facing function.
void readParticlesFromHDF(HDF5Handle &a_handle, ParticleData< P > &a_particles, const std::string &a_dataType)
Definition: ParticleIOI.H:244
void readDataChunk(size_t &offset, const hid_t &dataspace, const hid_t &dataset, const hid_t &H5T_type, const unsigned long dataLength, void *const data, const size_t stride=1, const size_t block=1)
Read chunk of data and upgrade offset. Not meant to be a user-facing function.
#define CH_assert(cond)
Definition: CHArray.H:37
A not-necessarily-disjoint collective of boxes.
Definition: BoxLayout.H:145
#define CH_START(tpointer)
Definition: CH_Timer.H:145
const BoxLayout & getBoxes() const
Get the BoxLayout on which this ParticleData.
Definition: ParticleDataI.H:193
void read_vect_from_header(HDF5Handle &a_handle, vector< T > &a_vect, const hid_t &H5T_type, const std::string &a_dataname)
A helper function that reads a vector of 'T's to the HDF5 file under a_dataname.
Definition: ParticleIOI.H:62
Definition: DataIterator.H:190
unsigned int size() const
Returns the total number of boxes in the BoxLayout.
void dataSize(const BaseFab< int > &item, Vector< int > &a_sizes, const Box &box, const Interval &comps)
Definition: CH_HDF5.H:830
#define CH_TIMER(name, tpointer)
Definition: CH_Timer.H:63
void push_back(const T &in)
Definition: Vector.H:295
void writeParticlesToHDF(HDF5Handle &a_handle, const ParticleData< P > &a_particles, const std::string &a_dataType)
Definition: ParticleIOI.H:93
Iterator over a List.
Definition: List.H:20
void add(const T &value)
Adds a copy of the value to the end of the List<T>.
Definition: ListImplem.H:49
double Real
Definition: REAL.H:33
unsigned int index(const LayoutIndex &index) const
Definition: BoxLayout.H:724
int length() const
Returns the number of objects in the List<T>.
Definition: ListImplem.H:56
static void Error(const char *const a_msg=m_nullString, int m_exitCode=CH_DEFAULT_ERROR_CODE)
Print out message to cerr and exit with the specified exit code.
void read_hdf_part_header(HDF5Handle &a_handle, Vector< Box > &a_grids, vector< unsigned long long > &a_particlesPerBox, const std::string &a_dataType, const std::string &a_path)
Read the particle header from the HDF5 file. Not meant to be a user-facing function.
const hid_t & groupID() const
void write_vect_to_header(HDF5Handle &a_handle, const vector< T > &a_vect, const hid_t &H5T_type, const std::string &a_dataname)
A helper function that writes a vector of 'T's to the HDF5 file under a_dataname. ...
Definition: ParticleIOI.H:26
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:469
void write_hdf_part_header(HDF5Handle &a_handle, const BoxLayout &a_grids, const vector< unsigned long long > &a_partPerBox, const std::string &a_dataType)
Write the particle header to the HDF5 file. Not meant to be a user-facing function.
#define CH_STOP(tpointer)
Definition: CH_Timer.H:150
Handle to a particular group in an HDF file.
Definition: CH_HDF5.H:292
bool ok() const
Return true if the iterator is not past the end of the list.
Definition: List.H:465
const std::string & getGroup() const
size_t size() const
Definition: Vector.H:192
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:117
#define _CHUNK
Definition: ParticleIOI.H:21
int procID()
local process ID
Definition: ParticleData.H:67