Chombo + EB  3.2
BlockWriteI.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 #ifndef _BLOCKWRITEI_H_
12 #define _BLOCKWRITEI_H_
13 
14 #include "LevelData.H"
15 #include "HDF5Portable.H"
16 #include "CH_HDF5.H"
17 #include <string>
18 #include <map>
19 #include "RealVect.H"
20 #include "CH_Timer.H"
21 #include "LoadBalance.H"
22 #include "LayoutIterator.H"
23 #include "Vector.H"
24 #include "memtrack.H"
25 #include "FluxBox.H"
26 #include "SPMD.H"
27 
28 //==================================================================//
29 template <class T> void
30 blockLocalOffsets(Vector<long long>& a_localSizes, //=numpts*ncomp
31  //for each box
32  long long & a_localTotalSize, //=numpts*ncomp
33  Vector<Box> & a_localBoxes,
34  const BoxLayoutData<T>& a_data,
35  const Interval& a_comps,
36  const IntVect& a_outputGhost)
37 {
38  if (!(T::preAllocatable()==0))
39  {
40  MayDay::Error("non static preallocatable data not covered yet.");
41  }
42 
43  Vector<int> thisSize(1);
44  const BoxLayout& layout = a_data.boxLayout();
45  T dummy;
46  unsigned int curIndex = 0;
47  DataIterator it(layout.dataIterator());
48 
49  a_localTotalSize = 0;
50  a_localSizes.resize(it.size());
51  a_localBoxes.resize(it.size());
52  for (it.reset(); it.ok(); ++it)
53  {
54  long long curSize = 0;
55  const Box& curBox = layout[it()];
56  Box region = curBox;
57  region.grow(a_outputGhost);
58  dataSize(dummy, thisSize, region, a_comps);
59 
60  curSize = thisSize[0];
61  a_localSizes[curIndex] = curSize;
62  a_localBoxes[curIndex] = curBox;
63 
64  a_localTotalSize += thisSize[0];
65  curIndex += 1;
66  }
67 }
68 //==================================================================//
69 template <class T> void
70 blockWriteToBuffer(void* a_buffer,
71  const Vector<long long>& a_sizes, //ncomp*npts
72  const BoxLayoutData<T>& a_data,
73  const Interval& a_comps,
74  const IntVect& a_outputGhost)
75 {
76  int* curPtr = (int*)(a_buffer);
77  int curIndex = 0;
78  for (DataIterator it = a_data.dataIterator(); it.ok(); ++it)
79  {
80  Vector<void*> buffers(1, (void*) curPtr);
81  const T& data = a_data[it()];
82  //unsigned int index = a_data.boxLayout().index(it());
83  Box box = a_data.box(it());
84  box.grow(a_outputGhost);
85  write(data, buffers, box, a_comps); //write T to buffer
86 
87  curPtr += a_sizes[curIndex];
88  curIndex += 1;
89  }
90 
91 }
92 //==================================================================//
93 int
94 gatherBoxesAndOffsets(long long& a_offsetThisProc,
95  long long& a_allProcSize,
96  Vector<long long>& a_globalOffsets,
97  Vector<Box>& a_globalBoxes,
98  const Vector<long long>& a_localBoxSizes,
99  const Vector<Box>& a_localBoxes,
100  const long long& a_localAllBoxSize) //= numpts * ncomp
101 
102 {
103  Vector<Vector<Box > > allBoxes;
104  Vector<Vector<long long> > allBoxSizes;
105  Vector<long long> allSizes;
106  int iprocdest = 0;
107  gather(allBoxes, a_localBoxes, iprocdest);
108  gather(allBoxSizes, a_localBoxSizes, iprocdest);
109  gather(allSizes, a_localAllBoxSize, iprocdest);
110  Vector<long long> offsetsProcs(numProc());
111  if (procID() == iprocdest)
112  {
113  a_allProcSize = 0;
114  long long numBoxes = 0;
115  for (int iproc = 0; iproc < numProc(); iproc++)
116  {
117  a_allProcSize += allSizes[iproc];
118  numBoxes += allBoxes[iproc].size();
119  }
120  a_globalBoxes. resize(numBoxes);
121  a_globalOffsets.resize(numBoxes+1); //need the last one to hold
122  //the size
123  long long curOffset = 0;
124  int curIndex = 0;
125  for (int iproc = 0; iproc < numProc(); iproc++)
126  {
127  for (int ivec = 0; ivec < allBoxSizes[iproc].size(); ivec++)
128  {
129  if (ivec == 0)
130  {
131  offsetsProcs[iproc] = curOffset;
132  }
133  long long boxSize = allBoxSizes[iproc][ivec];
134  Box curBox = allBoxes[iproc][ivec];
135 
136  a_globalBoxes [curIndex] = curBox;
137  a_globalOffsets[curIndex] = curOffset;
138 
139  curIndex += 1;
140  curOffset += boxSize;
141  }
142  }
143  if (curOffset != a_allProcSize)
144  {
145  return -43; //offsets should add up to size
146  }
147  a_globalOffsets[curIndex] = curOffset; //make last one the total size
148 
149  }
150  a_offsetThisProc = offsetsProcs[procID()];
151  broadcast(a_allProcSize, iprocdest);
152  broadcast(a_globalOffsets, iprocdest);
153  broadcast(a_globalBoxes, iprocdest);
154  return 0;
155 }
156 //==================================================================//
157 int
159  void* a_buffer,
160  const std::string& a_name,
161  Vector<Box>& a_boxes,
162  Vector<long long>& a_sizes, //ncomp*npts
163  const Vector<hid_t>& a_types,
164  const BoxLayout& a_layout,
165  const long long& a_thisprocsize)
166 {
167 
168  herr_t err;
169  //need the offset into the entire data set for my proc
170  long long offsetthisproc, sumprocsize;
171  Vector<Box> globalBoxes;
172  Vector<long long> globalOffsets;
173  //write boxes and offsets
174  //get all procs boxes
175  int ret = gatherBoxesAndOffsets(offsetthisproc, sumprocsize, globalOffsets, globalBoxes, a_sizes, a_boxes, a_thisprocsize);
176  if (ret!=0)
177  {
178  return ret;
179  }
180 
181  if (globalBoxes.size() != a_layout.size())
182  {
183  return -12;
184  }
185 
186  {//offset write
187  char offsetname[1024];
188  sprintf(offsetname, "%s:offsets=0",a_name.c_str());
189 
190  hsize_t flatdims[1];
191  flatdims[0] = globalOffsets.size();
192  hid_t offsetspace = H5Screate_simple(1, flatdims, NULL);
193  hid_t offsetdata = H5Dcreate(a_handle.groupID(), offsetname,
194  H5T_NATIVE_LLONG, offsetspace, H5P_DEFAULT);
195 
196  if (procID() == 0)
197  {
198  hid_t memdataspace = H5Screate_simple(1, flatdims, NULL);
199  err = H5Dwrite(offsetdata, H5T_NATIVE_LLONG, memdataspace, offsetspace,
200  H5P_DEFAULT, &(globalOffsets[0]));
201  if (err < 0) return err;
202  H5Sclose(memdataspace);
203  }
204 
205  H5Sclose(offsetspace);
206  H5Dclose(offsetdata);
207  }
208  { //data write
209  hsize_t hs_procsize[1];
210  hsize_t hs_allprocsize[1];
211  ch_offset_t ch_offset[1];
212 
213  ch_offset[0] = offsetthisproc;
214  hs_procsize[0] = a_thisprocsize; //size of buffer on this proc
215  hs_allprocsize[0] = sumprocsize; //size of buffer on all procs
216 
217  char dataname[1024];
218  sprintf(dataname, "%s:datatype=0",a_name.c_str());
219 
220  hid_t dataspace = H5Screate_simple(1, hs_allprocsize, NULL);
221  hid_t memdataspace = H5Screate_simple(1, hs_procsize, NULL);
222 
223  hid_t dataset = H5Dcreate(a_handle.groupID(), dataname,
224  a_types[0], dataspace, H5P_DEFAULT);
225 
226  //select where in the file it will be written
227  err = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET,
228  ch_offset, NULL,
229  hs_procsize, NULL);
230  if (err < 0) return err;
231 
232  err = H5Dwrite(dataset, a_types[0], memdataspace, dataspace,
233  H5P_DEFAULT, a_buffer);
234  if (err < 0) return err;
235 
236  H5Sclose(memdataspace);
237  H5Sclose(dataspace);
238  H5Dclose(dataset);
239  }
240  return 0;
241 }
242 
243 template <class T>
244 int blockWrite(HDF5Handle& a_handle,
245  const LevelData<T>& a_data,
246  const std::string& a_name,
247  const IntVect& a_outputGhost,
248  const Interval& a_in_comps)
249 {
250  HDF5HeaderData info;
251  info.m_intvect["ghost"] = a_data.ghostVect();
252  IntVect og(a_outputGhost);
253  og.min(a_data.ghostVect());
254  info.m_intvect["outputGhost"] = og;
255  std::string group = a_handle.getGroup();
256  a_handle.setGroup(group+"/"+a_name+"_attributes");
257  info.writeToFile(a_handle);
258  a_handle.setGroup(group);
259  return blockWrite(a_handle, (const BoxLayoutData<T>&)a_data, a_name, og, a_in_comps);
260 }
261 
262 //==================================================================//
263 template <class T> int
265  const BoxLayoutData<T>& a_data,
266  const std::string& a_name,
267  const IntVect& a_outputGhost,
268  const Interval& a_comps)
269 {
270  CH_TIME("block write Level");
271  int ret = 0;
272 
273  const BoxLayout& layout = a_data.boxLayout();
274  T dummy; // used for preallocatable methods for dumb compilers.
275 
276  Vector<hid_t> types;
277  dataTypes(types, dummy);
278  if (types.size() != 1)
279  {
280  MayDay::Error("premature generality, I only deal with types of size 1");
281  }
282 
283  //all this stuff is local to this proc
284  Vector<long long> sizes; //=numpts*ncomp for each box
285  Vector<Box> boxes;
286  long long totalsize; //=numpts*ncomp
287 
288  //calculate offset stuff on this proc and allocate and output to
289  //local buffer
290  blockLocalOffsets(sizes, totalsize, boxes, a_data, a_comps, a_outputGhost);
291 
292  //allocate buffer --- this totalsize is the local size
293  size_t type_size = H5Tget_size(types[0]);
294 
295  //the argument is the number of chars in the buffer length
296  void* buffer = mallocMT(totalsize*type_size);
297 
298  //write to the buffer
299  blockWriteToBuffer(buffer, sizes, a_data, a_comps, a_outputGhost);
300 
301  //write attributes
302  HDF5HeaderData info;
303  info.m_int["comps"] = a_comps.size();
304  info.m_string["objectType"] = name(dummy); //this is a function call
305  std::string group = a_handle.getGroup();
306 
307  a_handle.setGroup(group+"/"+a_name+"_attributes");
308  info.writeToFile(a_handle);
309  a_handle.setGroup(group);
310 
311  //exchange offset information
312  //output stuff to file
313  ret = blockWriteBufferToFile(a_handle,
314  buffer,
315  a_name,
316  boxes,
317  sizes,
318  types,
319  layout,
320  totalsize);
321  if (ret != 0)
322  {
323  return ret;
324  }
325 // ret = blockWriteBufferToFile(a_handle, buffer, a_name, boxes, offsets, layout, totalsize);
326 
327  //free the buffer memory
328  freeMT(buffer);
329  buffer = NULL;
330 
331  return ret;
332 }
333 
334 template <class T> int
336  const int& a_level,
337  const LevelData<T>& a_data,
338  const Real& a_dx,
339  const Real& a_dt,
340  const Real& a_time,
341  const Box& a_domain,
342  const int& a_refRatio,
343  const IntVect& a_outputGhost,
344  const Interval& a_comps)
345 {
346  int error;
347  char levelName[10];
348  std::string currentGroup = a_handle.getGroup();
349  sprintf(levelName, "/level_%i",a_level);
350  error = a_handle.setGroup(currentGroup + levelName);
351  if (error != 0) return 1;
352 
353  HDF5HeaderData meta;
354  meta.m_real["dx"] = a_dx;
355  meta.m_real["dt"] = a_dt;
356  meta.m_real["time"] = a_time;
357  meta.m_box["prob_domain"] = a_domain;
358  meta.m_int["ref_ratio"] = a_refRatio;
359 
360  error = meta.writeToFile(a_handle);
361  if (error != 0) return 2;
362 
363  error = write(a_handle, a_data.boxLayout());
364  if (error != 0) return 3;
365 
366  error = blockWrite(a_handle, a_data, std::string("data"), a_outputGhost, a_comps);
367  if (error != 0) return 4;
368 
369  a_handle.setGroup(currentGroup);
370 
371  return 0;
372 }
373 //==================================================================//
374 template <class T> int
376  const int& a_level,
377  LevelData<T>& a_data,
378  Real& a_dx,
379  Real& a_dt,
380  Real& a_time,
381  Box& a_domain,
382  int& a_refRatio,
383  const Interval& a_comps,
384  bool a_setGhost)
385 {
386  HDF5HeaderData header;
387  header.readFromFile(a_handle);
388  //unused
389  // int nComp = header.m_int["num_components"];
390 
391  int error;
392  char levelName[10];
393  std::string currentGroup = a_handle.getGroup();
394  sprintf(levelName, "/level_%i",a_level);
395  error = a_handle.setGroup(currentGroup + levelName);
396  if (error != 0) return 1;
397 
398  HDF5HeaderData meta;
399  error = meta.readFromFile(a_handle);
400  if (error != 0) return 2;
401  a_dx = meta.m_real["dx"];
402  a_dt = meta.m_real["dt"];
403  a_time = meta.m_real["time"];
404  a_domain = meta.m_box["prob_domain"];
405  a_refRatio = meta.m_int["ref_ratio"];
406 
407  //read in the boxes and load balance
408  Vector<Box> boxes;
409  error = read(a_handle, boxes);
410  Vector<int> procIDs;
411  LoadBalance(procIDs, boxes);
412 
413  //make the layout
414  DisjointBoxLayout layout(boxes, procIDs);
415 
416  layout.close();
417  if (error != 0) return 3;
418  //do the very complicated read
419  error = blockRead(a_handle, a_data, "data", layout, a_comps);
420 
421  if (error != 0) return 4;
422 
423  //return the handle to its input state
424  a_handle.setGroup(currentGroup);
425 
426  return 0;
427 }
428 
429 //////////
430 template <class T> int
431 blockRead(HDF5Handle& a_handle, LevelData<T>& a_data, const std::string& a_name,
432  const DisjointBoxLayout& a_layout, const Interval& a_comps, bool a_redefineData)
433 {
434  if (a_redefineData)
435  {
436  HDF5HeaderData info;
437  std::string group = a_handle.getGroup();
438  if (a_handle.setGroup(group+"/"+a_name+"_attributes"))
439  {
440  std::string message = "error opening "
441  +a_handle.getGroup()+"/"+a_name+"_attributes" ;
442  MayDay::Warning(message.c_str());
443  return 1;
444  }
445  info.readFromFile(a_handle);
446  a_handle.setGroup(group);
447  int ncomp = info.m_int["comps"];
448  IntVect ghost = info.m_intvect["ghost"];
449  if (a_comps.end() > 0 && ncomp < a_comps.end())
450  {
451  MayDay::Error("attempt to read component interval that is not available");
452  }
453  if (a_comps.size() == 0)
454  a_data.define(a_layout, ncomp, ghost);
455  else
456  a_data.define(a_layout, a_comps.size(), ghost);
457  }
458  return blockRead(a_handle, (BoxLayoutData<T>&)a_data, a_name, a_layout, a_comps, false);
459 
460 }
461 
462 template <class T>
463 int
464 blockRead(HDF5Handle& a_handle, BoxLayoutData<T>& a_data, const std::string& a_name,
465  const BoxLayout& a_layout, const Interval& a_comps, bool a_redefineData)
466 {
467 
468  int ret = read(a_handle, a_data, a_name, a_layout, a_comps, a_redefineData);
469  return ret;
470 }
471 
472 #endif
int blockWriteBufferToFile(HDF5Handle &a_handle, void *a_buffer, const std::string &a_name, Vector< Box > &a_boxes, Vector< long long > &a_sizes, const Vector< hid_t > &a_types, const BoxLayout &a_layout, const long long &a_thisprocsize)
Definition: BlockWriteI.H:158
map< std::string, int > m_int
Definition: CH_HDF5.H:548
IntVect & min(const IntVect &p)
Definition: IntVect.H:1115
#define freeMT(a_a)
Definition: memtrack.H:160
data to be added to HDF5 files.
Definition: CH_HDF5.H:517
A not-necessarily-disjoint collective of boxes.
Definition: BoxLayout.H:145
map< std::string, Real > m_real
Definition: CH_HDF5.H:545
#define mallocMT(a_a)
Definition: memtrack.H:159
void blockLocalOffsets(Vector< long long > &a_localSizes, long long &a_localTotalSize, Vector< Box > &a_localBoxes, const BoxLayoutData< T > &a_data, const Interval &a_comps, const IntVect &a_outputGhost)
Definition: BlockWriteI.H:30
int gatherBoxesAndOffsets(long long &a_offsetThisProc, long long &a_allProcSize, Vector< long long > &a_globalOffsets, Vector< Box > &a_globalBoxes, const Vector< long long > &a_localBoxSizes, const Vector< Box > &a_localBoxes, const long long &a_localAllBoxSize)
Definition: BlockWriteI.H:94
Definition: DataIterator.H:190
unsigned int size() const
Returns the total number of boxes in the BoxLayout.
int LoadBalance(Vector< Vector< int > > &a_procAssignments, Real &a_effRatio, const Vector< Vector< Box > > &a_Grids, const Vector< Vector< long > > &a_ComputeLoads, const Vector< int > &a_RefRatios, int a_nProc=numProc())
unsigned int numProc()
number of parallel processes
int size() const
Definition: Interval.H:75
void dataTypes(Vector< hid_t > &a_types, const BaseFab< int > &dummy)
Definition: CH_HDF5.H:799
void dataSize(const BaseFab< int > &item, Vector< int > &a_sizes, const Box &box, const Interval &comps)
Definition: CH_HDF5.H:830
int blockWriteLevel(HDF5Handle &a_handle, const int &a_level, const LevelData< T > &a_data, const Real &a_dx, const Real &a_dt, const Real &a_time, const Box &a_domain, const int &a_refRatio, const IntVect &a_outputGhost, const Interval &a_comps)
Definition: BlockWriteI.H:335
Definition: EBInterface.H:45
void resize(unsigned int isize)
Definition: Vector.H:346
virtual void close()
map< std::string, IntVect > m_intvect
Definition: CH_HDF5.H:554
void gather(Vector< T > &a_outVec, const T &a_input, int a_dest)
Definition: SPMDI.H:197
#define CH_TIME(name)
Definition: CH_Timer.H:82
Structure for passing component ranges in code.
Definition: Interval.H:23
const char * name(const FArrayBox &a_dummySpecializationArg)
Definition: CH_HDF5.H:864
new code
Definition: BoxLayoutData.H:170
Data on a BoxLayout.
Definition: BoxLayoutData.H:97
double Real
Definition: REAL.H:33
virtual void define(const DisjointBoxLayout &dp, int comps, const IntVect &ghost=IntVect::Zero, const DataFactory< T > &a_factory=DefaultDataFactory< T >())
Definition: LevelDataI.H:83
A BoxLayout that has a concept of disjointedness.
Definition: DisjointBoxLayout.H:30
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.
int write(HDF5Handle &a_handle, const BoxLayout &a_layout, const std::string &name="boxes")
writes BoxLayout to HDF5 file.
const BoxLayout & boxLayout() const
Definition: LayoutData.H:107
DataIterator dataIterator() const
Parallel iterator.
const hid_t & groupID() const
int blockWrite(HDF5Handle &a_handle, const LevelData< T > &a_data, const std::string &a_name, const IntVect &a_outputGhost, const Interval &a_in_comps)
Definition: BlockWriteI.H:244
hssize_t ch_offset_t
Definition: CH_HDF5.H:693
static void Warning(const char *const a_msg=m_nullString)
Print out message to cerr and continue.
A Rectangular Domain on an Integer Lattice.
Definition: Box.H:465
void blockWriteToBuffer(void *a_buffer, const Vector< long long > &a_sizes, const BoxLayoutData< T > &a_data, const Interval &a_comps, const IntVect &a_outputGhost)
Definition: BlockWriteI.H:70
int writeToFile(HDF5Handle &file) const
map< std::string, Box > m_box
Definition: CH_HDF5.H:557
Handle to a particular group in an HDF file.
Definition: CH_HDF5.H:292
An integer Vector in SpaceDim-dimensional space.
Definition: CHArray.H:42
const std::string & getGroup() const
DataIterator dataIterator() const
Definition: LayoutDataI.H:78
size_t size() const
Definition: Vector.H:192
Box & grow(int i)
grow functions
Definition: Box.H:2247
virtual bool ok() const
return true if this iterator is still in its Layout
Definition: LayoutIterator.H:117
int setGroup(const std::string &groupAbsPath)
void broadcast(T &a_inAndOut, int a_src)
broadcast to every process
Definition: SPMDI.H:207
int end() const
Definition: Interval.H:102
int blockReadLevel(HDF5Handle &a_handle, const int &a_level, LevelData< T > &a_data, Real &a_dx, Real &a_dt, Real &a_time, Box &a_domain, int &a_refRatio, const Interval &a_comps, bool a_setGhost)
Definition: BlockWriteI.H:375
map< std::string, std::string > m_string
Definition: CH_HDF5.H:551
int procID()
local process ID
Box box(const DataIndex &a_index) const
Definition: LayoutDataI.H:66
int readFromFile(HDF5Handle &file)
void read(T &item, Vector< Vector< char > > &a_allocatedBuffers, const Box &box, const Interval &comps)
Definition: CH_HDF5.H:49
int blockRead(HDF5Handle &a_handle, LevelData< T > &a_data, const std::string &a_name, const DisjointBoxLayout &a_layout, const Interval &a_comps, bool a_redefineData)
read LevelData named a_name from location specified by a_handle.
Definition: BlockWriteI.H:431
const IntVect & ghostVect() const
Definition: LevelData.H:170