Nektar++
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
FieldIOHdf5.cpp
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // File: FieldIOHdf5.cpp
4 //
5 // For more information, please see: http://www.nektar.info/
6 //
7 // The MIT License
8 //
9 // Copyright (c) 2006 Division of Applied Mathematics, Brown University (USA),
10 // Department of Aeronautics, Imperial College London (UK), and Scientific
11 // Computing and Imaging Institute, University of Utah (USA).
12 //
13 // License for the specific language governing rights and limitations under
14 // Permission is hereby granted, free of charge, to any person obtaining a
15 // copy of this software and associated documentation files (the "Software"),
16 // to deal in the Software without restriction, including without limitation
17 // the rights to use, copy, modify, merge, publish, distribute, sublicense,
18 // and/or sell copies of the Software, and to permit persons to whom the
19 // Software is furnished to do so, subject to the following conditions:
20 //
21 // The above copyright notice and this permission notice shall be included
22 // in all copies or substantial portions of the Software.
23 //
24 // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
25 // OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
26 // FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
27 // THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
28 // LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
29 // FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
30 // DEALINGS IN THE SOFTWARE.
31 //
32 // Description: I/O routines relating to Fields into HDF
33 //
34 ////////////////////////////////////////////////////////////////////////////////
35 
37 #include <boost/unordered_set.hpp>
38 
39 namespace berrc = boost::system::errc;
40 
41 namespace Nektar
42 {
43 namespace LibUtilities
44 {
45 namespace H5
46 {
47 
49 {
50  return PredefinedDataType::Native<int>();
51 }
52 
53 }
54 
55 std::string FieldIOHdf5::className =
57  "Hdf5", FieldIOHdf5::create, "HDF5-based output of field data.");
58 
59 /// Version of the Nektar++ HDF5 format, which is embedded into the main NEKTAR
60 /// group as an attribute.
61 const unsigned int FieldIOHdf5::FORMAT_VERSION = 1;
62 
63 // The following definitions allow us to consistently refer to indexes pulled
64 // out of the various datasets.
65 
66 /// A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the
67 /// position of the number of elements in decomposition (i.e. field definition).
68 const unsigned int FieldIOHdf5::ELEM_DCMP_IDX = 0;
69 /// A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the
70 /// position of the number of data points in decomposition (i.e. field
71 /// definition).
72 const unsigned int FieldIOHdf5::VAL_DCMP_IDX = 1;
73 /// A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the
74 /// position of the number of elements multiplied by the dimension of the
75 /// element, giving number of modes when variable polynomial order is defined.
76 const unsigned int FieldIOHdf5::ORDER_DCMP_IDX = 2;
77 /// A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the
78 /// position of the number of the number of y-planes for homogeneous
79 /// simulations.
80 const unsigned int FieldIOHdf5::HOMY_DCMP_IDX = 3;
81 /// A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the
82 /// position of the number of the number of z-planes for homogeneous
83 /// simulations.
84 const unsigned int FieldIOHdf5::HOMZ_DCMP_IDX = 4;
85 /// A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the
86 /// position of the number of the number of strips for homogeneous simulations.
87 const unsigned int FieldIOHdf5::HOMS_DCMP_IDX = 5;
88 /// The hash of the field definition information, which defines the name of the
89 /// attribute containing the field definition itself.
90 const unsigned int FieldIOHdf5::HASH_DCMP_IDX = 6;
91 /// A helper for FieldIOHdf5::v_Write. Describes the maximum number of items in
92 /// the decomposition per field definition.
93 const unsigned int FieldIOHdf5::MAX_DCMPS = FieldIOHdf5::HASH_DCMP_IDX + 1;
94 
95 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
96 /// elements in the cnt array.
97 const unsigned int FieldIOHdf5::ELEM_CNT_IDX = 0;
98 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
99 /// data points in the cnt array.
100 const unsigned int FieldIOHdf5::VAL_CNT_IDX = 1;
101 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
102 /// order points in the cnt array.
103 const unsigned int FieldIOHdf5::ORDER_CNT_IDX = 2;
104 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
105 /// homogeneous y-planes in the cnt array.
106 const unsigned int FieldIOHdf5::HOMY_CNT_IDX = 3;
107 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
108 /// homogeneous z-planes in the cnt array.
109 const unsigned int FieldIOHdf5::HOMZ_CNT_IDX = 4;
110 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
111 /// homogeneous strips in the cnt array.
112 const unsigned int FieldIOHdf5::HOMS_CNT_IDX = 5;
113 /// A helper for FieldIOHdf5::v_Write. Describes the maximum number of items in
114 /// the cnt array per field definition.
115 const unsigned int FieldIOHdf5::MAX_CNTS = FieldIOHdf5::HOMS_CNT_IDX + 1;
116 
117 /// A helper for FieldIOHdf5::v_Write. Describes the position of the element IDs
118 /// within the indexing set.
119 const unsigned int FieldIOHdf5::IDS_IDX_IDX = 0;
120 /// A helper for FieldIOHdf5::v_Write. Describes the position of the data size
121 /// within the indexing set.
122 const unsigned int FieldIOHdf5::DATA_IDX_IDX = 1;
123 /// A helper for FieldIOHdf5::v_Write. Describes the position of the element
124 /// order within the indexing set.
125 const unsigned int FieldIOHdf5::ORDER_IDX_IDX = 2;
126 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
127 /// y-planes within the indexing set.
128 const unsigned int FieldIOHdf5::HOMY_IDX_IDX = 3;
129 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
130 /// z-planes within the indexing set.
131 const unsigned int FieldIOHdf5::HOMZ_IDX_IDX = 4;
132 /// A helper for FieldIOHdf5::v_Write. Describes the position of the number of
133 /// homogeneous strips within the indexing set.
134 const unsigned int FieldIOHdf5::HOMS_IDX_IDX = 5;
135 /// A helper for FieldIOHdf5::v_Write. Describes the maximum number of items in
136 /// the indexing set.
137 const unsigned int FieldIOHdf5::MAX_IDXS = FieldIOHdf5::HOMS_IDX_IDX + 1;
138 
139 /**
140  * @brief Construct the FieldIO object for HDF5 output.
141  *
142  * @param pComm Communicator object.
143  * @param sharedFilesystem True if this system has a shared filesystem.
144  */
146  bool sharedFilesystem)
147  : FieldIO(pComm, sharedFilesystem)
148 {
149 }
150 
151 /**
152  * @brief Write a HDF5 file to @p outFile given the field definitions @p
153  * fielddefs, field data @p fielddata and metadata @p fieldmetadatamap.
154  *
155  * The writing strategy for HDF5 output is as follows:
156  *
157  * - Each rank determines the amount of data needed to be written into each
158  * dataset.
159  * - Each rank communicates its decomposition information to the root process.
160  * - The root processor initialises the output structure, writes the
161  * decomposition dataset and all the field definition information.
162  * - Other ranks may have field definitions that do not belong to the root
163  * process, in which case they open the file and append this (since
164  * attributes cannot be written in parallel).
165  * - Each of the other ranks writes their data contributions to the rest of
166  * the set.
167  *
168  * @param outFile Output filename.
169  * @param fielddefs Input field definitions.
170  * @param fielddata Input field data.
171  * @param fieldmetadatamap Field metadata.
172  */
173 void FieldIOHdf5::v_Write(const std::string &outFile,
174  std::vector<FieldDefinitionsSharedPtr> &fielddefs,
175  std::vector<std::vector<NekDouble> > &fielddata,
176  const FieldMetaDataMap &fieldmetadatamap,
177  const bool backup)
178 {
179  std::stringstream prfx;
180  prfx << m_comm->GetRank() << ": FieldIOHdf5::v_Write(): ";
181  double tm0 = 0.0, tm1 = 0.0;
182 
183  if (m_comm->GetRank() == 0)
184  {
185  tm0 = m_comm->Wtime();
186  }
187 
188  SetUpOutput(outFile, false, backup);
189 
190  // We make a number of assumptions in this code:
191  // 1. All element ids have the same type: unsigned int
192  // 2. All elements within a given field have the same number of values
193  // 3. All element values have the same type, NekDouble
194 
195  // Determine the root MPI process, i.e., the lowest ranked process handling
196  // nMaxFields fields, that will be responsible for writing our file.
197  ASSERTL1(fielddefs.size() == fielddata.size(),
198  prfx.str() + "fielddefs and fielddata have incompatible lengths.");
199 
200  size_t nFields = fielddefs.size();
201  size_t nMaxFields = nFields;
202  m_comm->AllReduce(nMaxFields, LibUtilities::ReduceMax);
203 
204  int root_rank = -1;
205  bool amRoot = false;
206  LibUtilities::CommSharedPtr max_fields_comm;
207 
208  if (m_comm->GetSize() > 1)
209  {
210  max_fields_comm = m_comm->CommCreateIf((nFields == nMaxFields) ? 1 : 0);
211  }
212  else
213  {
214  max_fields_comm = m_comm;
215  }
216 
217  if (max_fields_comm)
218  {
219  int rank = m_comm->GetRank();
220  root_rank = rank;
221  max_fields_comm->AllReduce(root_rank, LibUtilities::ReduceMin);
222  amRoot = (rank == root_rank);
223  if (!amRoot)
224  {
225  root_rank = -1;
226  }
227  }
228 
229  m_comm->AllReduce(root_rank, LibUtilities::ReduceMax);
230  ASSERTL1(root_rank >= 0 && root_rank < m_comm->GetSize(),
231  prfx.str() + "invalid root rank.");
232 
233  std::vector<uint64_t> decomps(nMaxFields * MAX_DCMPS, 0);
234  std::vector<uint64_t> all_hashes(nMaxFields * m_comm->GetSize(), 0);
235  std::vector<uint64_t> cnts(MAX_CNTS, 0);
236  std::vector<std::string> fieldNames(nFields);
237  std::vector<std::string> shapeStrings(nFields);
238  std::vector<std::vector<NekDouble> > homoLengths(nFields);
239  std::vector<std::vector<unsigned int> > homoSIDs(nFields),
240  homoYIDs(nFields), homoZIDs(nFields);
241  std::vector<std::vector<unsigned int> > numModesPerDirVar(nFields);
242  std::vector<std::string> numModesPerDirUni(nFields);
243 
244  int homDim = -1;
245  int varOrder = 0;
246 
247  for (int f = 0; f < nFields; ++f)
248  {
249  if (!fielddefs[f]->m_uniOrder)
250  {
251  varOrder = 1;
252  break;
253  }
254  }
255 
256  m_comm->AllReduce(varOrder, LibUtilities::ReduceMax);
257 
258  // Calculate the total number of elements handled by this MPI process and
259  // the total number of bytes required to store the elements. Base the name
260  // of each field on the hash of the field definition.
261  for (int f = 0; f < nFields; ++f)
262  {
263  ASSERTL1(fielddata[f].size() > 0,
264  prfx.str() +
265  "fielddata vector must contain at least one value.");
266  ASSERTL1(fielddata[f].size() ==
267  fielddefs[f]->m_fields.size() *
268  CheckFieldDefinition(fielddefs[f]),
269  prfx.str() + "fielddata vector has invalid size.");
270 
271  std::size_t nFieldElems = fielddefs[f]->m_elementIDs.size();
272  std::size_t nElemVals = fielddata[f].size();
273 
274  decomps[f * MAX_DCMPS + ELEM_DCMP_IDX] = nFieldElems;
275  decomps[f * MAX_DCMPS + VAL_DCMP_IDX] = nElemVals;
276 
277  cnts[ELEM_CNT_IDX] += nFieldElems;
278  cnts[VAL_CNT_IDX] += nElemVals;
279 
280  // Hash the field specification
281  std::stringstream hashStream;
282  std::size_t nSubFields = fielddefs[f]->m_fields.size();
283  for (int sf = 0; sf < nSubFields; ++sf)
284  {
285  hashStream << fielddefs[f]->m_fields[sf];
286  }
287 
288  nSubFields = fielddefs[f]->m_basis.size();
289  for (int sf = 0; sf < nSubFields; ++sf)
290  {
291  hashStream << fielddefs[f]->m_basis[sf];
292  }
293 
294  // Determine SHAPE attribute
295  std::stringstream shapeStringStream;
296  shapeStringStream << ShapeTypeMap[fielddefs[f]->m_shapeType];
297 
298  if (fielddefs[f]->m_numHomogeneousDir > 0)
299  {
300  if (homDim == -1)
301  {
302  homDim = fielddefs[f]->m_numHomogeneousDir;
303  }
304 
305  ASSERTL1(homDim == fielddefs[f]->m_numHomogeneousDir,
306  "HDF5 does not support variable homogeneous directions in "
307  "the same file.");
308 
309  shapeStringStream << "-HomogenousExp"
310  << fielddefs[f]->m_numHomogeneousDir << "D";
311  }
312 
313  if (fielddefs[f]->m_homoStrips)
314  {
315  shapeStringStream << "-Strips";
316  }
317 
318  shapeStrings[f] = shapeStringStream.str();
319  hashStream << shapeStringStream.str();
320 
321  // Determine HOMOGENEOUS attributes
322  if (fielddefs[f]->m_numHomogeneousDir)
323  {
324  nSubFields = fielddefs[f]->m_homogeneousLengths.size();
325  homoLengths[f].resize(nSubFields);
326  for (int sf = 0; sf < nSubFields; ++sf)
327  {
328  uint64_t len = fielddefs[f]->m_homogeneousLengths[sf];
329  hashStream << len;
330  homoLengths[f][sf] = len;
331  }
332 
333  nSubFields = fielddefs[f]->m_homogeneousYIDs.size();
334  if (nSubFields > 0)
335  {
336  homoYIDs[f].resize(nSubFields);
337  decomps[f * MAX_DCMPS + HOMY_DCMP_IDX] = nSubFields;
338  cnts[HOMY_CNT_IDX] += nSubFields;
339  for (int sf = 0; sf < nSubFields; ++sf)
340  {
341  homoYIDs[f][sf] = fielddefs[f]->m_homogeneousYIDs[sf];
342  }
343  }
344 
345  nSubFields = fielddefs[f]->m_homogeneousZIDs.size();
346  if (nSubFields > 0)
347  {
348  homoZIDs[f].resize(nSubFields);
349  decomps[f * MAX_DCMPS + HOMZ_DCMP_IDX] = nSubFields;
350  cnts[HOMZ_CNT_IDX] += nSubFields;
351  for (int sf = 0; sf < nSubFields; ++sf)
352  {
353  homoZIDs[f][sf] = fielddefs[f]->m_homogeneousZIDs[sf];
354  }
355  }
356 
357  nSubFields = fielddefs[f]->m_homogeneousSIDs.size();
358  if (nSubFields > 0)
359  {
360  homoSIDs[f].resize(nSubFields);
361  decomps[f * MAX_DCMPS + HOMS_DCMP_IDX] = nSubFields;
362  cnts[HOMS_CNT_IDX] += nSubFields;
363  for (int sf = 0; sf < nSubFields; ++sf)
364  {
365  homoSIDs[f][sf] = fielddefs[f]->m_homogeneousSIDs[sf];
366  }
367  }
368  }
369 
370  if (fielddefs[f]->m_uniOrder)
371  {
372  std::vector<unsigned int> elemModes(fielddefs[f]->m_basis.size());
373 
374  for (std::vector<int>::size_type i = 0;
375  i < fielddefs[f]->m_basis.size(); ++i)
376  {
377  elemModes[i] = fielddefs[f]->m_numModes[i];
378  }
379 
380  if (varOrder)
381  {
382  for (std::vector<int>::size_type i = 0; i < nFieldElems; ++i)
383  {
384  std::copy(elemModes.begin(), elemModes.end(),
385  std::back_inserter(numModesPerDirVar[f]));
386  }
387  decomps[f * MAX_DCMPS + ORDER_DCMP_IDX] =
388  nFieldElems * elemModes.size();
389  cnts[ORDER_CNT_IDX] += nFieldElems * elemModes.size();
390  }
391  else
392  {
393  stringstream numModesStringStream;
394  numModesStringStream << "UNIORDER:";
395  for (std::vector<int>::size_type i = 0;
396  i < elemModes.size(); i++)
397  {
398  if (i > 0)
399  {
400  numModesStringStream << ",";
401  }
402  numModesStringStream << elemModes[i];
403  }
404 
405  numModesPerDirUni[f] = numModesStringStream.str();
406  hashStream << numModesPerDirUni[f];
407  }
408  }
409  else
410  {
411  numModesPerDirVar[f] = fielddefs[f]->m_numModes;
412  decomps[f * MAX_DCMPS + ORDER_DCMP_IDX] =
413  fielddefs[f]->m_numModes.size();
414  cnts[ORDER_CNT_IDX] += fielddefs[f]->m_numModes.size();
415  }
416 
417  boost::hash<std::string> string_hasher;
418  std::stringstream fieldNameStream;
419  uint64_t fieldDefHash = string_hasher(hashStream.str());
420 
421  decomps[f * MAX_DCMPS + HASH_DCMP_IDX] = fieldDefHash;
422  all_hashes[m_comm->GetRank() * nMaxFields + f] = fieldDefHash;
423 
424  fieldNameStream << fieldDefHash;
425  fieldNames[f] = fieldNameStream.str();
426  }
427 
428  // Gather information from all MPI processes
429  std::vector<uint64_t> all_cnts = m_comm->Gather(root_rank, cnts);
430  std::vector<uint64_t> all_idxs(m_comm->GetSize() * MAX_IDXS, 0);
431  std::vector<uint64_t> all_decomps = m_comm->Gather(root_rank, decomps);
432  std::vector<uint64_t> all_dsetsize(MAX_CNTS, 0);
433 
434  // The root rank creates the file layout from scratch
435  if (amRoot)
436  {
437  H5::FileSharedPtr outfile = H5::File::Create(outFile, H5F_ACC_TRUNC);
438  ASSERTL1(outfile, prfx.str() + "cannot create HDF5 file.");
439  H5::GroupSharedPtr root = outfile->CreateGroup("NEKTAR");
440  ASSERTL1(root, prfx.str() + "cannot create root group.");
441  TagWriterSharedPtr info_writer(new H5TagWriter(root));
442  AddInfoTag(info_writer, fieldmetadatamap);
443 
444  // Record file format version as attribute in main group.
445  root->SetAttribute("FORMAT_VERSION", FORMAT_VERSION);
446 
447  // Calculate the indexes to be used by each MPI process when reading the
448  // IDS and DATA datasets
449  std::size_t nTotElems = 0, nTotVals = 0, nTotOrder = 0;
450  std::size_t nTotHomY = 0, nTotHomZ = 0, nTotHomS = 0;
451  int nRanks = m_comm->GetSize();
452  for (int r = 0; r < nRanks; ++r)
453  {
454  all_idxs[r * MAX_IDXS + IDS_IDX_IDX] = nTotElems;
455  all_idxs[r * MAX_IDXS + DATA_IDX_IDX] = nTotVals;
456  all_idxs[r * MAX_IDXS + ORDER_IDX_IDX] = nTotOrder;
457  all_idxs[r * MAX_IDXS + HOMY_IDX_IDX] = nTotHomY;
458  all_idxs[r * MAX_IDXS + HOMZ_IDX_IDX] = nTotHomZ;
459  all_idxs[r * MAX_IDXS + HOMS_IDX_IDX] = nTotHomS;
460 
461  nTotElems += all_cnts[r * MAX_CNTS + ELEM_CNT_IDX];
462  nTotVals += all_cnts[r * MAX_CNTS + VAL_CNT_IDX];
463  nTotOrder += all_cnts[r * MAX_CNTS + ORDER_CNT_IDX];
464  nTotHomY += all_cnts[r * MAX_CNTS + HOMY_CNT_IDX];
465  nTotHomZ += all_cnts[r * MAX_CNTS + HOMZ_CNT_IDX];
466  nTotHomS += all_cnts[r * MAX_CNTS + HOMS_CNT_IDX];
467  }
468 
469  all_dsetsize[ELEM_CNT_IDX ] = nTotElems;
470  all_dsetsize[VAL_CNT_IDX ] = nTotVals;
471  all_dsetsize[ORDER_CNT_IDX] = nTotOrder;
472  all_dsetsize[HOMY_CNT_IDX ] = nTotHomY;
473  all_dsetsize[HOMZ_CNT_IDX ] = nTotHomZ;
474  all_dsetsize[HOMS_CNT_IDX ] = nTotHomS;
475 
476  // Create DECOMPOSITION dataset: basic field info for each MPI process
477  H5::DataTypeSharedPtr decomps_type =
478  H5::DataType::OfObject(all_decomps[0]);
479  H5::DataSpaceSharedPtr decomps_space =
480  H5::DataSpace::OneD(all_decomps.size());
481  H5::DataSetSharedPtr decomps_dset =
482  root->CreateDataSet("DECOMPOSITION", decomps_type, decomps_space);
483  ASSERTL1(decomps_dset,
484  prfx.str() + "cannot create DECOMPOSITION dataset.");
485 
486  // Create IDS dataset: element ids
487  H5::DataTypeSharedPtr ids_type =
488  H5::DataType::OfObject(fielddefs[0]->m_elementIDs[0]);
489  H5::DataSpaceSharedPtr ids_space = H5::DataSpace::OneD(nTotElems);
490  H5::DataSetSharedPtr ids_dset =
491  root->CreateDataSet("ELEMENTIDS", ids_type, ids_space);
492  ASSERTL1(ids_dset, prfx.str() + "cannot create ELEMENTIDS dataset.");
493 
494  // Create DATA dataset: element data
495  H5::DataTypeSharedPtr data_type =
496  H5::DataType::OfObject(fielddata[0][0]);
497  H5::DataSpaceSharedPtr data_space = H5::DataSpace::OneD(nTotVals);
498  H5::DataSetSharedPtr data_dset =
499  root->CreateDataSet("DATA", data_type, data_space);
500  ASSERTL1(data_dset, prfx.str() + "cannot create DATA dataset.");
501 
502  // Create HOMOGENEOUSYIDS dataset: homogeneous y-plane IDs
503  if (nTotHomY > 0)
504  {
505  H5::DataTypeSharedPtr homy_type =
506  H5::DataType::OfObject(homoYIDs[0][0]);
507  H5::DataSpaceSharedPtr homy_space = H5::DataSpace::OneD(nTotHomY);
508  H5::DataSetSharedPtr homy_dset =
509  root->CreateDataSet("HOMOGENEOUSYIDS", homy_type, homy_space);
510  ASSERTL1(homy_dset,
511  prfx.str() + "cannot create HOMOGENEOUSYIDS dataset.");
512  }
513 
514  // Create HOMOGENEOUSYIDS dataset: homogeneous z-plane IDs
515  if (nTotHomZ > 0)
516  {
517  H5::DataTypeSharedPtr homz_type =
518  H5::DataType::OfObject(homoZIDs[0][0]);
519  H5::DataSpaceSharedPtr homz_space = H5::DataSpace::OneD(nTotHomZ);
520  H5::DataSetSharedPtr homz_dset =
521  root->CreateDataSet("HOMOGENEOUSZIDS", homz_type, homz_space);
522  ASSERTL1(data_dset,
523  prfx.str() + "cannot create HOMOGENEOUSZIDS dataset.");
524  }
525 
526  // Create HOMOGENEOUSSIDS dataset: homogeneous strip IDs
527  if (nTotHomS > 0)
528  {
529  H5::DataTypeSharedPtr homs_type =
530  H5::DataType::OfObject(homoSIDs[0][0]);
531  H5::DataSpaceSharedPtr homs_space = H5::DataSpace::OneD(nTotHomS);
532  H5::DataSetSharedPtr homs_dset =
533  root->CreateDataSet("HOMOGENEOUSSIDS", homs_type, homs_space);
534  ASSERTL1(homs_dset,
535  prfx.str() + "cannot create HOMOGENEOUSSIDS dataset.");
536  }
537 
538  // Create POLYORDERS dataset: elemental polynomial orders
539  if (varOrder)
540  {
541  H5::DataTypeSharedPtr order_type =
542  H5::DataType::OfObject(numModesPerDirVar[0][0]);
543  H5::DataSpaceSharedPtr order_space = H5::DataSpace::OneD(nTotOrder);
544  H5::DataSetSharedPtr order_dset =
545  root->CreateDataSet("POLYORDERS", order_type, order_space);
546  ASSERTL1(order_dset,
547  prfx.str() + "cannot create POLYORDERS dataset.");
548  }
549  }
550 
551  m_comm->Bcast(all_dsetsize, root_rank);
552 
553  // Datasets, root group and HDF5 file are all closed automatically since
554  // they are now out of scope. Now we need to determine which process will
555  // write the group representing the field description in the HDF5 file. This
556  // next block of code performs this by finding all unique hashes and then
557  // determining one process that will create (possibly more than one) group
558  // for that hash. An alternative would be to communicate the field
559  // information to the root processor, but this is a bit convoluted.
560 
561  // This set stores the unique hashes.
562  std::set<uint64_t> hashToProc;
563  // This map takes ranks to hashes this process will write.
564  std::map<int, std::vector<uint64_t> > writingProcs;
565 
566  // Gather all field hashes to every processor.
567  m_comm->AllReduce(all_hashes, LibUtilities::ReduceMax);
568 
569  for (int n = 0; n < m_comm->GetSize(); ++n)
570  {
571  for (int i = 0; i < nMaxFields; ++i)
572  {
573  uint64_t hash = all_hashes[n*nMaxFields + i];
574 
575  // Note hash can be zero if, on this process, nFields < nMaxFields.
576  if (hashToProc.find(hash) != hashToProc.end() || hash == 0)
577  {
578  continue;
579  }
580  hashToProc.insert(hash);
581  writingProcs[n].push_back(hash);
582  }
583  }
584 
585  // Having constructed the map, go ahead and write the attributes out.
586  map<int, std::vector<uint64_t> >::iterator sIt;
587  for (sIt = writingProcs.begin(); sIt != writingProcs.end(); sIt++)
588  {
589  int rank = sIt->first;
590 
591  // Write out this rank's groups.
592  if (m_comm->GetRank() == rank)
593  {
594  H5::PListSharedPtr serialProps = H5::PList::Default();
596 
597  // Reopen the file
598  H5::FileSharedPtr outfile =
599  H5::File::Open(outFile, H5F_ACC_RDWR, serialProps);
600  ASSERTL1(outfile, prfx.str() + "cannot open HDF5 file.");
601  H5::GroupSharedPtr root = outfile->OpenGroup("NEKTAR");
602  ASSERTL1(root, prfx.str() + "cannot open root group.");
603 
604  // Write a HDF5 group for each field
605  hashToProc.clear();
606  for (int i = 0; i < sIt->second.size(); ++i)
607  {
608  for (int f = 0; f < nFields; ++f)
609  {
610  if (sIt->second[i] !=
611  all_hashes[m_comm->GetRank() * nMaxFields + f] ||
612  hashToProc.find(sIt->second[i]) != hashToProc.end())
613  {
614  continue;
615  }
616 
617  hashToProc.insert(sIt->second[i]);
618 
619  // Just in case we've already written this
620  H5::GroupSharedPtr field_group =
621  root->CreateGroup(fieldNames[f]);
622  ASSERTL1(field_group,
623  prfx.str() + "cannot create field group.");
624  field_group->SetAttribute("FIELDS", fielddefs[f]->m_fields);
625  field_group->SetAttribute("BASIS", fielddefs[f]->m_basis);
626  field_group->SetAttribute("SHAPE", shapeStrings[f]);
627 
628  if (homoLengths[f].size() > 0)
629  {
630  field_group->SetAttribute("HOMOGENEOUSLENGTHS",
631  homoLengths[f]);
632  }
633 
634  // If the field has only uniform order, we write the order
635  // into the NUMMODESPERDIR attribute. Otherwise, we'll go
636  // ahead and assume everything is mixed and fix this in the
637  // read later if required.
638  if (!varOrder)
639  {
640  field_group->SetAttribute("NUMMODESPERDIR",
641  numModesPerDirUni[f]);
642  }
643  else
644  {
645  std::string numModesPerDir = "MIXORDER";
646  field_group->SetAttribute("NUMMODESPERDIR",
647  numModesPerDir);
648  }
649  }
650  }
651  }
652 
653  // We block to avoid more than one processor opening the file at a time.
654  m_comm->Block();
655  }
656 
657  // Write the DECOMPOSITION dataset
658  if (amRoot)
659  {
660  H5::PListSharedPtr serialProps = H5::PList::Default();
662 
663  // Reopen the file
664  H5::FileSharedPtr outfile =
665  H5::File::Open(outFile, H5F_ACC_RDWR, serialProps);
666  ASSERTL1(outfile, prfx.str() + "cannot open HDF5 file.");
667  H5::GroupSharedPtr root = outfile->OpenGroup("NEKTAR");
668  ASSERTL1(root, prfx.str() + "cannot open root group.");
669 
670  // Write the DECOMPOSITION dataset
671  H5::DataSetSharedPtr decomps_dset = root->OpenDataSet("DECOMPOSITION");
672  ASSERTL1(decomps_dset,
673  prfx.str() + "cannot open DECOMPOSITION dataset.");
674 
675  H5::DataSpaceSharedPtr decomps_fspace = decomps_dset->GetSpace();
676  ASSERTL1(decomps_fspace,
677  prfx.str() + "cannot open DECOMPOSITION filespace.");
678 
679  decomps_fspace->SelectRange(0, all_decomps.size());
680  decomps_dset->Write(all_decomps, decomps_fspace, writeSR);
681  }
682 
683  // Initialise the dataset indexes for all MPI processes
684  std::vector<uint64_t> idx = m_comm->Scatter(root_rank, all_idxs);
685  uint64_t ids_i = idx[IDS_IDX_IDX];
686  uint64_t data_i = idx[DATA_IDX_IDX];
687  uint64_t order_i = idx[ORDER_IDX_IDX];
688  uint64_t homy_i = idx[HOMY_IDX_IDX];
689  uint64_t homz_i = idx[HOMZ_IDX_IDX];
690  uint64_t homs_i = idx[HOMS_IDX_IDX];
691 
692  // Set properties for parallel file access (if we're in parallel)
693  H5::PListSharedPtr parallelProps = H5::PList::Default();
695  if (m_comm->GetSize() > 1)
696  {
697  // Use MPI/O to access the file
698  parallelProps = H5::PList::FileAccess();
699  parallelProps->SetMpio(m_comm);
700  // Use collective IO
701  writePL = H5::PList::DatasetXfer();
702  writePL->SetDxMpioCollective();
703  }
704 
705  // Reopen the file
706  H5::FileSharedPtr outfile =
707  H5::File::Open(outFile, H5F_ACC_RDWR, parallelProps);
708  ASSERTL1(outfile, prfx.str() + "cannot open HDF5 file.");
709  H5::GroupSharedPtr root = outfile->OpenGroup("NEKTAR");
710  ASSERTL1(root, prfx.str() + "cannot open root group.");
711 
712  m_comm->Block();
713 
714  // all HDF5 groups have now been created. Open the IDS dataset and
715  // associated data space
716  H5::DataSetSharedPtr ids_dset = root->OpenDataSet("ELEMENTIDS");
717  ASSERTL1(ids_dset, prfx.str() + "cannot open ELEMENTIDS dataset.");
718  H5::DataSpaceSharedPtr ids_fspace = ids_dset->GetSpace();
719  ASSERTL1(ids_fspace, prfx.str() + "cannot open ELEMENTIDS filespace.");
720 
721  // Open the DATA dataset and associated data space
722  H5::DataSetSharedPtr data_dset = root->OpenDataSet("DATA");
723  ASSERTL1(data_dset, prfx.str() + "cannot open DATA dataset.");
724  H5::DataSpaceSharedPtr data_fspace = data_dset->GetSpace();
725  ASSERTL1(data_fspace, prfx.str() + "cannot open DATA filespace.");
726 
727  // Open the optional datasets and data spaces.
728  H5::DataSetSharedPtr order_dset, homy_dset, homz_dset, homs_dset;
729  H5::DataSpaceSharedPtr order_fspace , homy_fspace, homz_fspace, homs_fspace;
730 
731  if (all_dsetsize[ORDER_CNT_IDX])
732  {
733  order_dset = root->OpenDataSet("POLYORDERS");
734  ASSERTL1(order_dset, prfx.str() + "cannot open POLYORDERS dataset.");
735  order_fspace = order_dset->GetSpace();
736  ASSERTL1(order_fspace, prfx.str() + "cannot open POLYORDERS filespace.");
737  }
738 
739  if (all_dsetsize[HOMY_CNT_IDX])
740  {
741  homy_dset = root->OpenDataSet("HOMOGENEOUSYIDS");
742  ASSERTL1(homy_dset, prfx.str() + "cannot open HOMOGENEOUSYIDS dataset.");
743  homy_fspace = homy_dset->GetSpace();
744  ASSERTL1(homy_fspace, prfx.str() + "cannot open HOMOGENEOUSYIDS filespace.");
745  }
746 
747  if (all_dsetsize[HOMZ_CNT_IDX])
748  {
749  homz_dset = root->OpenDataSet("HOMOGENEOUSZIDS");
750  ASSERTL1(homz_dset, prfx.str() + "cannot open HOMOGENEOUSZIDS dataset.");
751  homz_fspace = homz_dset->GetSpace();
752  ASSERTL1(homz_fspace, prfx.str() + "cannot open HOMOGENEOUSZIDS filespace.");
753  }
754 
755  if (all_dsetsize[HOMS_CNT_IDX])
756  {
757  homs_dset = root->OpenDataSet("HOMOGENEOUSSIDS");
758  ASSERTL1(homs_dset, prfx.str() + "cannot open HOMOGENEOUSSIDS dataset.");
759  homs_fspace = homs_dset->GetSpace();
760  ASSERTL1(homs_fspace, prfx.str() + "cannot open HOMOGENEOUSSIDS filespace.");
761  }
762 
763  // Write the data
764  for (int f = 0; f < nFields; ++f)
765  {
766  // write the element ids
767  std::size_t nFieldElems = fielddefs[f]->m_elementIDs.size();
768  ids_fspace->SelectRange(ids_i, nFieldElems);
769  ids_dset->Write(fielddefs[f]->m_elementIDs, ids_fspace, writePL);
770  ids_i += nFieldElems;
771 
772  // write the element values
773  std::size_t nFieldVals = fielddata[f].size();
774  data_fspace->SelectRange(data_i, nFieldVals);
775  data_dset->Write(fielddata[f], data_fspace, writePL);
776  data_i += nFieldVals;
777  }
778 
779  if (order_dset)
780  {
781  for (int f = 0; f < nFields; ++f)
782  {
783  std::size_t nOrders = numModesPerDirVar[f].size();
784  order_fspace->SelectRange(order_i, nOrders);
785  order_dset->Write(numModesPerDirVar[f], order_fspace, writePL);
786  order_i += nOrders;
787  }
788  }
789 
790  if (homy_dset)
791  {
792  for (int f = 0; f < nFields; ++f)
793  {
794  std::size_t nYIDs = homoYIDs[f].size();
795  homy_fspace->SelectRange(homy_i, nYIDs);
796  homy_dset->Write(homoYIDs[f], homy_fspace, writePL);
797  homy_i += nYIDs;
798  }
799  }
800 
801  if (homz_dset)
802  {
803  for (int f = 0; f < nFields; ++f)
804  {
805  std::size_t nZIDs = homoZIDs[f].size();
806  homz_fspace->SelectRange(homz_i, nZIDs);
807  homz_dset->Write(homoZIDs[f], homz_fspace, writePL);
808  homz_i += nZIDs;
809  }
810  }
811 
812  if (homs_dset)
813  {
814  for (int f = 0; f < nFields; ++f)
815  {
816  std::size_t nSIDs = homoSIDs[f].size();
817  homs_fspace->SelectRange(homs_i, nSIDs);
818  homs_dset->Write(homoSIDs[f], homs_fspace, writePL);
819  homs_i += nSIDs;
820  }
821  }
822 
823  for (int f = nFields; f < nMaxFields; ++f)
824  {
825  // this MPI process is handling fewer than nMaxFields fields
826  // so, since this is a collective operation
827  // just rewrite the element ids and values of the last field
828  ids_dset->Write(
829  fielddefs[nFields - 1]->m_elementIDs, ids_fspace, writePL);
830  data_dset->Write(fielddata[nFields - 1], data_fspace, writePL);
831  }
832 
833  m_comm->Block();
834 
835  // all data has been written
836  if (m_comm->GetRank() == 0)
837  {
838  tm1 = m_comm->Wtime();
839  cout << " (" << tm1 - tm0 << "s, HDF5)" << endl;
840  }
841 }
842 
843 /**
844  * @brief Import a HDF5 format file.
845  *
846  * @param finfilename Input filename
847  * @param fielddefs Field definitions of resulting field
848  * @param fielddata Field data of resulting field
849  * @param fieldinfomap Field metadata of resulting field
850  * @param ElementIDs If specified, contains the list of element IDs on
851  * this rank. The resulting field definitions will only
852  * contain data for the element IDs specified in this
853  * array.
854  */
855 void FieldIOHdf5::v_Import(const std::string &infilename,
856  std::vector<FieldDefinitionsSharedPtr> &fielddefs,
857  std::vector<std::vector<NekDouble> > &fielddata,
858  FieldMetaDataMap &fieldinfomap,
859  const Array<OneD, int> &ElementIDs)
860 {
861  std::stringstream prfx;
862  int nRanks = m_comm->GetSize();
863 
864  // Set properties for parallel file access (if we're in parallel)
865  H5::PListSharedPtr parallelProps = H5::PList::Default();
867 
868  if (nRanks > 1)
869  {
870  // Use MPI/O to access the file
871  parallelProps = H5::PList::FileAccess();
872  parallelProps->SetMpio(m_comm);
873  // Use collective IO
874  readPL = H5::PList::DatasetXfer();
875  readPL->SetDxMpioCollective();
876  }
877 
879  infilename, parallelProps);
880 
881  // Open the root group of the hdf5 file
883  boost::static_pointer_cast<H5DataSource>(dataSource);
884  ASSERTL1(h5, prfx.str() + "cannot open HDF5 file.");
885  H5::GroupSharedPtr root = h5->Get()->OpenGroup("NEKTAR");
886  ASSERTL1(root, prfx.str() + "cannot open root group.");
887 
888  // Check format version
889  unsigned int formatVersion;
890  H5::Group::AttrIterator attrIt = root->attr_begin();
891  H5::Group::AttrIterator attrEnd = root->attr_end();
892  for (; attrIt != attrEnd; ++attrIt)
893  {
894  if (*attrIt == "FORMAT_VERSION")
895  {
896  break;
897  }
898  }
899 
900  ASSERTL0(attrIt != attrEnd,
901  "Unable to determine Nektar++ HDF5 file version");
902  root->GetAttribute("FORMAT_VERSION", formatVersion);
903 
904  ASSERTL0(formatVersion <= FORMAT_VERSION,
905  "File format if " + infilename + " is higher than supported in "
906  "this version of Nektar++");
907 
908  // Open the datasets
909  H5::DataSetSharedPtr decomps_dset = root->OpenDataSet("DECOMPOSITION");
910  ASSERTL1(decomps_dset, prfx.str() + "cannot open DECOMPOSITION dataset.");
911 
912  H5::DataSetSharedPtr ids_dset = root->OpenDataSet("ELEMENTIDS");
913  ASSERTL1(ids_dset, prfx.str() + "cannot open ELEMENTIDS dataset.");
914 
915  H5::DataSetSharedPtr data_dset = root->OpenDataSet("DATA");
916  ASSERTL1(data_dset, prfx.str() + "cannot open DATA dataset.");
917 
918  // Open the dataset file spaces
919  H5::DataSpaceSharedPtr decomps_fspace = decomps_dset->GetSpace();
920  ASSERTL1(decomps_fspace,
921  prfx.str() + "cannot open DECOMPOSITION filespace.");
922 
923  H5::DataSpaceSharedPtr ids_fspace = ids_dset->GetSpace();
924  ASSERTL1(ids_fspace, prfx.str() + "cannot open ELEMENTIDS filespace.");
925 
926  H5::DataSpaceSharedPtr data_fspace = data_dset->GetSpace();
927  ASSERTL1(data_fspace, prfx.str() + "cannot open DATA filespace.");
928 
929  // Read entire IDS data set to get list of global IDs.
930  std::vector<uint64_t> ids;
931 
932  ids_dset->Read(ids, ids_fspace, readPL);
933 
934  boost::unordered_set<uint64_t> toread;
935  if (ElementIDs != NullInt1DArray)
936  {
937  for (uint64_t i = 0; i < ElementIDs.num_elements(); ++i)
938  {
939  toread.insert(ElementIDs[i]);
940  }
941  }
942 
943  std::vector<uint64_t> decomps;
944  decomps_dset->Read(decomps, decomps_fspace, readPL);
945 
946  size_t nDecomps = decomps.size() / MAX_DCMPS;
947  size_t cnt = 0, cnt2 = 0;
948 
949  // Mapping from each decomposition to offsets in the data array.
950  vector<OffsetHelper> decompsToOffsets (nDecomps);
951 
952  // Mapping from each group's hash to a vector of element IDs. Note this has
953  // to be unsigned int, since that's what we use in FieldDefinitions.
954  map<uint64_t, vector<unsigned int> > groupsToElmts;
955 
956  // Mapping from each group's hash to each of the decompositions.
957  map<uint64_t, set<uint64_t> > groupsToDecomps;
958 
959  // True if we are pulling element IDs from ElementIDs.
960  bool selective = toread.size() > 0;
961 
962  // Counters for data offsets
963  OffsetHelper running;
964 
965  for (size_t i = 0; i < nDecomps; ++i, cnt += MAX_DCMPS)
966  {
967  uint64_t nElmt = decomps[cnt + ELEM_DCMP_IDX];
968  uint64_t groupHash = decomps[cnt + HASH_DCMP_IDX];
969 
970  vector<uint64_t> tmp;
971 
972  if (selective)
973  {
974  for (size_t j = 0; j < nElmt; ++j)
975  {
976  uint64_t elmtId = ids[cnt2 + j];
977  if (toread.find(elmtId) != toread.end())
978  {
979  tmp.push_back(elmtId);
980  }
981  }
982  }
983  else
984  {
985  tmp.insert(
986  tmp.begin(), ids.begin() + cnt2, ids.begin() + cnt2 + nElmt);
987  }
988 
989  vector<unsigned int> tmp2(nElmt);
990  for (size_t j = 0; j < nElmt; ++j)
991  {
992  tmp2[j] = ids[cnt2+j];
993  }
994 
995  cnt2 += nElmt;
996 
997  if (tmp.size() > 0)
998  {
999  groupsToDecomps[groupHash].insert(i);
1000  }
1001 
1002  groupsToElmts[i] = tmp2;
1003  decompsToOffsets[i] = running;
1004 
1005  running.data += decomps[cnt + VAL_DCMP_IDX];
1006  running.order += decomps[cnt + ORDER_DCMP_IDX];
1007  running.homy += decomps[cnt + HOMY_DCMP_IDX];
1008  running.homz += decomps[cnt + HOMZ_DCMP_IDX];
1009  running.homs += decomps[cnt + HOMS_DCMP_IDX];
1010  }
1011 
1012  map<uint64_t, set<uint64_t> >::iterator gIt;
1013  for (gIt = groupsToDecomps.begin(); gIt != groupsToDecomps.end(); ++gIt)
1014  {
1015  // Select region from dataset for this decomposition.
1017  for (sIt = gIt->second.begin(); sIt != gIt->second.end(); ++sIt)
1018  {
1019  std::stringstream fieldNameStream;
1020  fieldNameStream << gIt->first;
1021 
1022  FieldDefinitionsSharedPtr fielddef =
1024  ImportFieldDef(readPL, root, decomps, *sIt, decompsToOffsets[*sIt],
1025  fieldNameStream.str(), fielddef);
1026 
1027  fielddef->m_elementIDs = groupsToElmts[*sIt];
1028  fielddefs.push_back(fielddef);
1029 
1030  if (fielddata != NullVectorNekDoubleVector)
1031  {
1032  std::vector<NekDouble> decompFieldData;
1034  readPL, data_dset, data_fspace,
1035  decompsToOffsets[*sIt].data, decomps, *sIt, fielddef,
1036  decompFieldData);
1037  fielddata.push_back(decompFieldData);
1038  }
1039  }
1040  }
1041 
1042  m_comm->Block();
1043 }
1044 
1045 /**
1046  * @brief Import field definitions from a HDF5 document.
1047  *
1048  * @param readPL Reading parameter list.
1049  * @param root Root group containing field definitions.
1050  * @param group Group name to process.
1051  * @param def On output contains field definitions.
1052  */
1054  H5::PListSharedPtr readPL,
1055  H5::GroupSharedPtr root,
1056  std::vector<uint64_t> &decomps,
1057  uint64_t decomp,
1058  OffsetHelper offset,
1059  std::string group,
1061 {
1062  std::stringstream prfx;
1063  prfx << m_comm->GetRank() << ": FieldIOHdf5::ImportFieldDefsHdf5(): ";
1064 
1065  H5::GroupSharedPtr field = root->OpenGroup(group);
1066  ASSERTL1(field, prfx.str() + "cannot open field group, " + group + '.');
1067 
1068  def->m_uniOrder = false;
1069 
1070  H5::Group::AttrIterator attrIt = field->attr_begin();
1071  H5::Group::AttrIterator attrEnd = field->attr_end();
1072  for (; attrIt != attrEnd; ++attrIt)
1073  {
1074  const std::string &attrName = *attrIt;
1075  if (attrName == "FIELDS")
1076  {
1077  field->GetAttribute(attrName, def->m_fields);
1078  }
1079  else if (attrName == "SHAPE")
1080  {
1081  std::string shapeString;
1082  field->GetAttribute(attrName, shapeString);
1083 
1084  // check to see if homogeneous expansion and if so
1085  // strip down the shapeString definition
1086  size_t loc;
1087  //---> this finds the first location of 'n'!
1088  if (shapeString.find("Strips") != string::npos)
1089  {
1090  def->m_homoStrips = true;
1091  }
1092 
1093  if ((loc = shapeString.find_first_of("-")) != string::npos)
1094  {
1095  if (shapeString.find("Exp1D") != string::npos)
1096  {
1097  def->m_numHomogeneousDir = 1;
1098  }
1099  else // HomogeneousExp1D
1100  {
1101  def->m_numHomogeneousDir = 2;
1102  }
1103 
1104  shapeString.erase(loc, shapeString.length());
1105  }
1106 
1107  // get the geometrical shape
1108  bool valid = false;
1109  for (unsigned int j = 0; j < SIZE_ShapeType; j++)
1110  {
1111  if (ShapeTypeMap[j] == shapeString)
1112  {
1113  def->m_shapeType = (ShapeType)j;
1114  valid = true;
1115  break;
1116  }
1117  }
1118 
1119  ASSERTL0(valid, prfx.str() + std::string(
1120  "unable to correctly parse the shape type: ")
1121  .append(shapeString).c_str());
1122  }
1123  else if (attrName == "BASIS")
1124  {
1125  field->GetAttribute(attrName, def->m_basis);
1126  // check the basis is in range
1127  std::vector<BasisType>::const_iterator bIt = def->m_basis.begin();
1128  std::vector<BasisType>::const_iterator bEnd = def->m_basis.end();
1129  for (; bIt != bEnd; ++bIt)
1130  {
1131  BasisType bt = *bIt;
1132  ASSERTL0(bt >= 0 && bt < SIZE_BasisType,
1133  prfx.str() +
1134  "unable to correctly parse the basis types.");
1135  }
1136  }
1137  else if (attrName == "HOMOGENEOUSLENGTHS")
1138  {
1139  field->GetAttribute(attrName, def->m_homogeneousLengths);
1140  }
1141  else if (attrName == "NUMMODESPERDIR")
1142  {
1143  std::string numModesPerDir;
1144  field->GetAttribute(attrName, numModesPerDir);
1145 
1146  if (strstr(numModesPerDir.c_str(), "UNIORDER:"))
1147  {
1148  def->m_uniOrder = true;
1149  bool valid = ParseUtils::GenerateOrderedVector(
1150  numModesPerDir.c_str() + 9, def->m_numModes);
1151  ASSERTL0(valid,
1152  prfx.str() +
1153  "unable to correctly parse the number of modes.");
1154  }
1155  }
1156  else if (attrName == "POINTSTYPE")
1157  {
1158  std::string pointsString;
1159  field->GetAttribute(attrName, pointsString);
1160  def->m_pointsDef = true;
1161 
1162  std::vector<std::string> pointsStrings;
1164  pointsString.c_str(), pointsStrings);
1165  ASSERTL0(valid,
1166  prfx.str() +
1167  "unable to correctly parse the points types.");
1168  for (std::vector<std::string>::size_type i = 0;
1169  i < pointsStrings.size();
1170  i++)
1171  {
1172  valid = false;
1173  for (unsigned int j = 0; j < SIZE_PointsType; j++)
1174  {
1175  if (kPointsTypeStr[j] == pointsStrings[i])
1176  {
1177  def->m_points.push_back((PointsType)j);
1178  valid = true;
1179  break;
1180  }
1181  }
1182 
1183  ASSERTL0(
1184  valid,
1185  prfx.str() +
1186  std::string(
1187  "unable to correctly parse the points type: ")
1188  .append(pointsStrings[i])
1189  .c_str());
1190  }
1191  }
1192  else if (attrName == "NUMPOINTSPERDIR")
1193  {
1194  std::string numPointsPerDir;
1195  field->GetAttribute(attrName, numPointsPerDir);
1196  def->m_numPointsDef = true;
1197 
1198  bool valid = ParseUtils::GenerateOrderedVector(
1199  numPointsPerDir.c_str(), def->m_numPoints);
1200  ASSERTL0(valid,
1201  prfx.str() +
1202  "unable to correctly parse the number of points.");
1203  }
1204  else
1205  {
1206  std::string errstr("unknown attribute: ");
1207  errstr += attrName;
1208  ASSERTL1(false, prfx.str() + errstr.c_str());
1209  }
1210  }
1211 
1212  if (def->m_numHomogeneousDir >= 1)
1213  {
1214  H5::DataSetSharedPtr homz_dset = root->OpenDataSet("HOMOGENEOUSZIDS");
1215  H5::DataSpaceSharedPtr homz_fspace = homz_dset->GetSpace();
1216  uint64_t dsize = decomps[decomp * MAX_DCMPS + HOMZ_DCMP_IDX];
1217  homz_fspace->SelectRange(offset.homz, dsize);
1218  homz_dset->Read(def->m_homogeneousZIDs, homz_fspace, readPL);
1219  }
1220 
1221  if (def->m_numHomogeneousDir >= 2)
1222  {
1223  H5::DataSetSharedPtr homy_dset = root->OpenDataSet("HOMOGENEOUSYIDS");
1224  H5::DataSpaceSharedPtr homy_fspace = homy_dset->GetSpace();
1225  uint64_t dsize = decomps[decomp * MAX_DCMPS + HOMY_DCMP_IDX];
1226  homy_fspace->SelectRange(offset.homy, dsize);
1227  homy_dset->Read(def->m_homogeneousYIDs, homy_fspace, readPL);
1228  }
1229 
1230  if (def->m_homoStrips)
1231  {
1232  H5::DataSetSharedPtr homs_dset = root->OpenDataSet("HOMOGENEOUSSIDS");
1233  H5::DataSpaceSharedPtr homs_fspace = homs_dset->GetSpace();
1234  uint64_t dsize = decomps[decomp * MAX_DCMPS + HOMS_DCMP_IDX];
1235  homs_fspace->SelectRange(offset.homs, dsize);
1236  homs_dset->Read(def->m_homogeneousSIDs, homs_fspace, readPL);
1237  }
1238 
1239  if (!def->m_uniOrder)
1240  {
1241  H5::DataSetSharedPtr order_dset = root->OpenDataSet("POLYORDERS");
1242  H5::DataSpaceSharedPtr order_fspace = order_dset->GetSpace();
1243  uint64_t dsize = decomps[decomp * MAX_DCMPS + ORDER_DCMP_IDX];
1244  order_fspace->SelectRange(offset.order, dsize);
1245  order_dset->Read(def->m_numModes, order_fspace, readPL);
1246  }
1247 }
1248 
1249 /**
1250  * @brief Import the field data from the HDF5 document.
1251  *
1252  * @param readPL Reading parameter list.
1253  * @param data_dset Pointer to the `DATA` dataset.
1254  * @param data_fspace Pointer to the `DATA` data space.
1255  * @param data_i Index in the `DATA` dataset to start reading from.
1256  * @param decomps Information from the `DECOMPOSITION` dataset.
1257  * @param decomp Index of the decomposition.
1258  * @param fielddef Field definitions for this file
1259  * @param fielddata On return contains resulting field data.
1260  */
1262  H5::PListSharedPtr readPL,
1263  H5::DataSetSharedPtr data_dset,
1264  H5::DataSpaceSharedPtr data_fspace,
1265  uint64_t data_i,
1266  std::vector<uint64_t> &decomps,
1267  uint64_t decomp,
1268  const FieldDefinitionsSharedPtr fielddef,
1269  std::vector<NekDouble> &fielddata)
1270 {
1271  std::stringstream prfx;
1272  prfx << m_comm->GetRank() << ": FieldIOHdf5::ImportFieldData(): ";
1273 
1274  uint64_t nElemVals = decomps[decomp * MAX_DCMPS + VAL_DCMP_IDX];
1275  uint64_t nFieldVals = nElemVals;
1276 
1277  data_fspace->SelectRange(data_i, nFieldVals);
1278  data_dset->Read(fielddata, data_fspace, readPL);
1279 
1280  int datasize = CheckFieldDefinition(fielddef);
1281  ASSERTL0(
1282  fielddata.size() == datasize * fielddef->m_fields.size(),
1283  prfx.str() +
1284  "input data is not the same length as header information.");
1285 }
1286 
1287 /**
1288  * @brief Import field metadata from @p filename and return the data source
1289  * which wraps @p filename.
1290  *
1291  * @param filename Input filename.
1292  * @param fieldmetadatamap Resulting field metadata from @p dataSource.
1293  */
1295  const std::string &filename,
1296  FieldMetaDataMap &fieldmetadatamap)
1297 {
1298  H5::PListSharedPtr parallelProps = H5::PList::Default();
1299  DataSourceSharedPtr ans = H5DataSource::create(filename, parallelProps);
1300  ImportHDF5FieldMetaData(ans, fieldmetadatamap);
1301  return ans;
1302 }
1303 
1304 /**
1305  * @brief Import field metadata from @p dataSource.
1306  *
1307  * @param dataSource Input datasource, which should be a H5DataSource.
1308  * @param fieldmetadatamap Resulting field metadata from @p dataSource.
1309  */
1311  FieldMetaDataMap &fieldmetadatamap)
1312 {
1313  H5DataSourceSharedPtr hdf =
1314  boost::static_pointer_cast<H5DataSource>(dataSource);
1315 
1316  H5::GroupSharedPtr master = hdf->Get()->OpenGroup("NEKTAR");
1317  // New metadata format only in HDF
1318  H5::GroupSharedPtr metadata = master->OpenGroup("Metadata");
1319 
1320  if (metadata)
1321  {
1322  H5::Group::AttrIterator param = metadata->attr_begin(),
1323  pEnd = metadata->attr_end();
1324  for (; param != pEnd; ++param)
1325  {
1326  std::string paramString = *param;
1327  if (paramString != "Provenance")
1328  {
1329  std::string paramBodyStr;
1330  metadata->GetAttribute(paramString, paramBodyStr);
1331  fieldmetadatamap[paramString] = paramBodyStr;
1332  }
1333  }
1334  }
1335 }
1336 
1337 }
1338 }
boost::shared_ptr< H5DataSource > H5DataSourceSharedPtr
Definition: FieldIOHdf5.h:89
#define ASSERTL0(condition, msg)
Definition: ErrorUtil.hpp:198
static bool GenerateOrderedStringVector(const char *const str, std::vector< std::string > &vec)
Definition: ParseUtils.hpp:143
boost::shared_ptr< Group > GroupSharedPtr
Definition: FieldIOHdf5.h:50
FieldIOHdf5(LibUtilities::CommSharedPtr pComm, bool sharedFilesystem)
Construct the FieldIO object for HDF5 output.
static bool GenerateOrderedVector(const char *const str, std::vector< unsigned int > &vec)
Definition: ParseUtils.hpp:97
static const unsigned int HOMY_CNT_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of homogeneous y-planes in th...
Definition: FieldIOHdf5.h:190
static boost::shared_ptr< DataType > AllocateSharedPtr()
Allocate a shared pointer from the memory pool.
static FileSharedPtr Create(const std::string &filename, unsigned mode, PListSharedPtr createPL=PList::Default(), PListSharedPtr accessPL=PList::Default())
Definition: H5.cpp:555
virtual void v_Write(const std::string &outFile, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata, const FieldMetaDataMap &fieldinfomap=NullFieldMetaDataMap, const bool backup=false)
Write a HDF5 file to outFile given the field definitions fielddefs, field data fielddata and metadata...
const std::string kPointsTypeStr[]
Definition: Foundations.hpp:69
static const unsigned int MAX_DCMPS
A helper for FieldIOHdf5::v_Write. Describes the maximum number of items in the decomposition per fie...
Definition: FieldIOHdf5.h:185
static const unsigned int HOMS_IDX_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of homogeneous strips within ...
Definition: FieldIOHdf5.h:200
boost::shared_ptr< FieldDefinitions > FieldDefinitionsSharedPtr
Definition: FieldIO.h:181
static DataTypeSharedPtr OfObject(const T &obj)
Definition: H5.h:378
static const unsigned int ORDER_IDX_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the element order within the indexing se...
Definition: FieldIOHdf5.h:197
boost::shared_ptr< DataSet > DataSetSharedPtr
Definition: H5.h:100
static DataTypeSharedPtr GetType()
Definition: H5.h:538
void ImportHDF5FieldMetaData(DataSourceSharedPtr dataSource, FieldMetaDataMap &fieldmetadatamap)
Import field metadata from dataSource.
boost::shared_ptr< File > FileSharedPtr
Definition: H5.h:98
boost::shared_ptr< DataType > DataTypeSharedPtr
Definition: H5.h:86
boost::shared_ptr< PList > PListSharedPtr
Definition: H5.h:102
FieldIOFactory & GetFieldIOFactory()
Returns the FieldIO factory.
Definition: FieldIO.cpp:74
std::map< std::string, std::string > FieldMetaDataMap
Definition: FieldIO.h:54
int CheckFieldDefinition(const FieldDefinitionsSharedPtr &fielddefs)
Check field definitions for correctness and return storage size.
Definition: FieldIO.cpp:545
const char *const ShapeTypeMap[]
Definition: ShapeType.hpp:66
boost::shared_ptr< DataSource > DataSourceSharedPtr
Definition: FieldIO.h:81
static const unsigned int HASH_DCMP_IDX
The hash of the field definition information, which defines the name of the attribute containing the ...
Definition: FieldIOHdf5.h:184
static DataSourceSharedPtr create(const std::string &fn, H5::PListSharedPtr parallelProps)
Static constructor for this data source.
Definition: FieldIOHdf5.h:79
boost::shared_ptr< Comm > CommSharedPtr
Pointer to a Communicator object.
Definition: Comm.h:55
void ImportFieldDef(H5::PListSharedPtr readPL, H5::GroupSharedPtr root, std::vector< uint64_t > &decomps, uint64_t decomp, OffsetHelper offset, std::string group, FieldDefinitionsSharedPtr def)
Import field definitions from a HDF5 document.
static const unsigned int FORMAT_VERSION
Version of the Nektar++ HDF5 format, which is embedded into the main NEKTAR group as an attribute...
Definition: FieldIOHdf5.h:176
boost::shared_ptr< DataSpace > DataSpaceSharedPtr
Definition: H5.h:88
static DataSpaceSharedPtr OneD(hsize_t size)
Definition: H5.cpp:403
static const unsigned int HOMY_IDX_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of y-planes within the indexi...
Definition: FieldIOHdf5.h:198
static const unsigned int MAX_IDXS
A helper for FieldIOHdf5::v_Write. Describes the maximum number of items in the indexing set...
Definition: FieldIOHdf5.h:201
static FileSharedPtr Open(const std::string &filename, unsigned mode, PListSharedPtr accessPL=PList::Default())
Definition: H5.cpp:563
static std::vector< std::vector< NekDouble > > NullVectorNekDoubleVector
static PListSharedPtr DatasetXfer()
Properties for raw data transfer.
Definition: H5.cpp:114
int GetSize(const ConstArray< OneD, NekDouble > &x)
Definition: StdExpUtil.cpp:111
static const unsigned int HOMZ_CNT_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of homogeneous z-planes in th...
Definition: FieldIOHdf5.h:191
static const unsigned int VAL_DCMP_IDX
A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the position of the number of ...
Definition: FieldIOHdf5.h:179
static const unsigned int IDS_IDX_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the element IDs within the indexing set...
Definition: FieldIOHdf5.h:195
std::string SetUpOutput(const std::string outname, bool perRank, bool backup=false)
Set up the filesystem ready for output.
Definition: FieldIO.cpp:400
LibUtilities::CommSharedPtr m_comm
Communicator to use when writing parallel format.
Definition: FieldIO.h:265
static std::string className
Name of class.
Definition: FieldIOHdf5.h:212
StandardMatrixTag boost::call_traits< LhsDataType >::const_reference rhs typedef NekMatrix< LhsDataType, StandardMatrixTag >::iterator iterator
virtual void v_Import(const std::string &infilename, std::vector< FieldDefinitionsSharedPtr > &fielddefs, std::vector< std::vector< NekDouble > > &fielddata=NullVectorNekDoubleVector, FieldMetaDataMap &fieldinfomap=NullFieldMetaDataMap, const Array< OneD, int > &ElementIDs=NullInt1DArray)
Import a HDF5 format file.
static const unsigned int MAX_CNTS
A helper for FieldIOHdf5::v_Write. Describes the maximum number of items in the cnt array per field d...
Definition: FieldIOHdf5.h:193
static const unsigned int ELEM_CNT_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of elements in the cnt array...
Definition: FieldIOHdf5.h:187
static const unsigned int ORDER_CNT_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of order points in the cnt ar...
Definition: FieldIOHdf5.h:189
static const unsigned int HOMY_DCMP_IDX
A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the position of the number of ...
Definition: FieldIOHdf5.h:181
static const unsigned int DATA_IDX_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the data size within the indexing set...
Definition: FieldIOHdf5.h:196
static const unsigned int HOMS_CNT_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of homogeneous strips in the ...
Definition: FieldIOHdf5.h:192
static const unsigned int HOMZ_DCMP_IDX
A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the position of the number of ...
Definition: FieldIOHdf5.h:182
static const unsigned int ELEM_DCMP_IDX
A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the position of the number of ...
Definition: FieldIOHdf5.h:178
static const unsigned int HOMZ_IDX_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of z-planes within the indexi...
Definition: FieldIOHdf5.h:199
void ImportFieldData(H5::PListSharedPtr readPL, H5::DataSetSharedPtr data_dset, H5::DataSpaceSharedPtr data_fspace, uint64_t data_i, std::vector< uint64_t > &decomps, uint64_t decomp, const FieldDefinitionsSharedPtr fielddef, std::vector< NekDouble > &fielddata)
Import the field data from the HDF5 document.
Class for operating on Nektar++ input/output files.
Definition: FieldIO.h:224
static PListSharedPtr FileAccess()
Properties for file access.
Definition: H5.cpp:96
H5::FileSharedPtr Get()
Get H5::FileSharedPtr reference to file.
Definition: FieldIOHdf5.h:67
boost::shared_ptr< TagWriter > TagWriterSharedPtr
Definition: FieldIO.h:71
static FieldIOSharedPtr create(LibUtilities::CommSharedPtr pComm, bool sharedFilesystem)
Creates an instance of this class.
Definition: FieldIOHdf5.h:204
static const unsigned int HOMS_DCMP_IDX
A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the position of the number of ...
Definition: FieldIOHdf5.h:183
static Array< OneD, int > NullInt1DArray
static const unsigned int ORDER_DCMP_IDX
A helper for FieldIOHdf5::v_Write and FieldIOHdf5::v_Import. Describes the position of the number of ...
Definition: FieldIOHdf5.h:180
#define ASSERTL1(condition, msg)
Assert Level 1 – Debugging which is used whether in FULLDEBUG or DEBUG compilation mode...
Definition: ErrorUtil.hpp:228
void AddInfoTag(TagWriterSharedPtr root, const FieldMetaDataMap &fieldmetadatamap)
Add provenance information to the field metadata map.
Definition: FieldIO.cpp:334
static const unsigned int VAL_CNT_IDX
A helper for FieldIOHdf5::v_Write. Describes the position of the number of data points in the cnt arr...
Definition: FieldIOHdf5.h:188
static PListSharedPtr Default()
Default options.
Definition: H5.cpp:78
virtual DataSourceSharedPtr v_ImportFieldMetaData(const std::string &filename, FieldMetaDataMap &fieldmetadatamap)
Import field metadata from filename and return the data source which wraps filename.
tKey RegisterCreatorFunction(tKey idKey, CreatorFunction classCreator, tDescription pDesc="")
Register a class with the factory.
Definition: NekFactory.hpp:215