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