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