Chaste  Release::2018.1
Hdf5DataWriter.cpp
1 /*
2 
3 Copyright (c) 2005-2018, University of Oxford.
4 All rights reserved.
5 
6 University of Oxford means the Chancellor, Masters and Scholars of the
7 University of Oxford, having an administrative office at Wellington
8 Square, Oxford OX1 2JD, UK.
9 
10 This file is part of Chaste.
11 
12 Redistribution and use in source and binary forms, with or without
13 modification, are permitted provided that the following conditions are met:
14  * Redistributions of source code must retain the above copyright notice,
15  this list of conditions and the following disclaimer.
16  * Redistributions in binary form must reproduce the above copyright notice,
17  this list of conditions and the following disclaimer in the documentation
18  and/or other materials provided with the distribution.
19  * Neither the name of the University of Oxford nor the names of its
20  contributors may be used to endorse or promote products derived from this
21  software without specific prior written permission.
22 
23 THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
24 AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
25 IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
26 ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
27 LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
28 CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE
29 GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
30 HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT
31 LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
32 OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
33 
34 */
35 
36 /*
37  * Implementation file for Hdf5DataWriter class.
38  *
39  */
40 #include <set>
41 #include <cstring> //For strcmp etc. Needed in gcc-4.4
42 #include <boost/scoped_array.hpp>
43 
44 #include "Hdf5DataWriter.hpp"
45 
46 #include "Exception.hpp"
47 #include "OutputFileHandler.hpp"
48 #include "PetscTools.hpp"
49 #include "Version.hpp"
50 #include "MathsCustomFunctions.hpp"
51 
53  const std::string& rDirectory,
54  const std::string& rBaseName,
55  bool cleanDirectory,
56  bool extendData,
57  std::string datasetName,
58  bool useCache)
59  : AbstractHdf5Access(rDirectory, rBaseName, datasetName),
60  mrVectorFactory(rVectorFactory),
61  mCleanDirectory(cleanDirectory),
62  mUseExistingFile(extendData),
63  mIsInDefineMode(true),
64  mIsFixedDimensionSet(false),
65  mEstimatedUnlimitedLength(1u),
66  mFileFixedDimensionSize(0u),
67  mDataFixedDimensionSize(0u),
68  mLo(mrVectorFactory.GetLow()),
69  mHi(mrVectorFactory.GetHigh()),
70  mNumberOwned(0u),
71  mOffset(0u),
72  mNeedExtend(false),
73  mUseMatrixForIncompleteData(false),
74  mCurrentTimeStep(0),
75  mSinglePermutation(nullptr),
76  mDoublePermutation(nullptr),
77  mSingleIncompleteOutputMatrix(nullptr),
78  mDoubleIncompleteOutputMatrix(nullptr),
79  mUseOptimalChunkSizeAlgorithm(true),
80  mNumberOfChunks(0),
81  mChunkTargetSize(0x20000), // 128 K
82  mAlignment(0), // No alignment
83  mUseCache(useCache),
84  mCacheFirstTimeStep(0u)
85 {
86  mChunkSize[0] = 0;
87  mChunkSize[1] = 0;
88  mChunkSize[2] = 0;
89  mFixedChunkSize[0] = 0;
90  mFixedChunkSize[1] = 0;
91  mFixedChunkSize[2] = 0;
92 
94  {
95  EXCEPTION("You are asking to delete a file and then extend it, change arguments to constructor.");
96  }
97 
98  if (!mUseExistingFile && mDatasetName != "Data")
99  {
100  //User is trying to add a new dataset, but they are not extending a existing file
101  EXCEPTION("Adding new data only makes sense when extending an existing file");
102  }
103 
104  if (mUseExistingFile)
105  {
106  // Variables should already be defined if we are extending.
107  mIsInDefineMode = false;
108 
109  // If the file exists already, open it - this call will check it exists.
110  OpenFile();
111 
112  // If the dataset we are interested in doesn't exist then close the file
113  // We will go on to define variables and open the file as usual (except for it pre-existing).
115  {
116  //std::cout << "Dataset: " << mDatasetName << " doesn't exist in the file.\n";
117  H5Fclose(mFileId);
118  mIsInDefineMode = true;
119  }
120  // If dataset does exist then leave file open and try to extend it.
121  else
122  {
123  // Where to find the file
124  assert(mCleanDirectory==false);
125 
126  mVariablesDatasetId = H5Dopen(mFileId, mDatasetName.c_str(), H5P_DEFAULT);
127  hid_t variables_dataspace = H5Dget_space(mVariablesDatasetId);
128  //unsigned variables_dataset_rank = H5Sget_simple_extent_ndims(variables_dataspace);
129  hsize_t dataset_max_sizes[DATASET_DIMS];
130  H5Sget_simple_extent_dims(variables_dataspace, mDatasetDims, dataset_max_sizes); // put dims into mDatasetDims
131  H5Sclose(variables_dataspace);
132 
133  // Check that an unlimited dimension is defined
134  if (dataset_max_sizes[0] != H5S_UNLIMITED)
135  {
136  H5Dclose(mVariablesDatasetId);
137  H5Fclose(mFileId);
138  EXCEPTION("Tried to open a datafile for extending which doesn't have an unlimited dimension.");
139  }
141 
142  // Sanity check other dimension sizes
143  for (unsigned i=1; i<DATASET_DIMS; i++) // Zero is excluded since it is unlimited
144  {
145  assert(mDatasetDims[i] == dataset_max_sizes[i]);
146  }
148  mIsFixedDimensionSet = true;
149  mVariables.reserve(mDatasetDims[2]);
150 
151  // Figure out what the variables are
152  hid_t attribute_id = H5Aopen_name(mVariablesDatasetId, "Variable Details");
153  hid_t attribute_type = H5Aget_type(attribute_id);
154  hid_t attribute_space = H5Aget_space(attribute_id);
155  hsize_t attr_dataspace_dim;
156  H5Sget_simple_extent_dims(attribute_space, &attr_dataspace_dim, nullptr);
157  unsigned num_columns = H5Sget_simple_extent_npoints(attribute_space);
158  assert(num_columns == mDatasetDims[2]); // I think...
159 
160  char* string_array = (char *)malloc(sizeof(char)*MAX_STRING_SIZE*(int)num_columns);
161  H5Aread(attribute_id, attribute_type, string_array);
162 
163  // Loop over columns/variables
164  for (unsigned index=0; index<num_columns; index++)
165  {
166  // Read the string from the raw vector
167  std::string column_name_unit(&string_array[MAX_STRING_SIZE*index]);
168 
169  // Find location of unit name
170  size_t name_length = column_name_unit.find('(');
171  size_t unit_length = column_name_unit.find(')') - name_length - 1;
172 
173  // Create the variable
174  DataWriterVariable var;
175  var.mVariableName = column_name_unit.substr(0, name_length);
176  var.mVariableUnits = column_name_unit.substr(name_length+1, unit_length);
177  mVariables.push_back(var);
178  }
179 
180  // Free memory, release ids
181  free(string_array);
182  H5Tclose(attribute_type);
183  H5Sclose(attribute_space);
184  H5Aclose(attribute_id);
185 
186  // Now deal with time
188 
189  // How many time steps have been written so far?
190  hid_t timestep_dataspace = H5Dget_space(mUnlimitedDatasetId);
191  hsize_t num_timesteps;
192  H5Sget_simple_extent_dims(timestep_dataspace, &num_timesteps, nullptr);
193  H5Sclose(timestep_dataspace);
194  mCurrentTimeStep = (long)num_timesteps - 1;
196 
197  // Incomplete data?
198  attribute_id = H5Aopen_name(mVariablesDatasetId, "IsDataComplete");
199  if (attribute_id < 0)
200  {
201  // LCOV_EXCL_START
202  // Old format, before we added the attribute.
203  EXCEPTION("Extending old-format files isn't supported.");
204  // LCOV_EXCL_STOP
205  }
206  else
207  {
208  attribute_type = H5Aget_type(attribute_id);
209  attribute_space = H5Aget_space(attribute_id);
210  unsigned is_data_complete;
211  H5Aread(attribute_id, H5T_NATIVE_UINT, &is_data_complete);
212  mIsDataComplete = (is_data_complete == 1) ? true : false;
213  H5Tclose(attribute_type);
214  H5Sclose(attribute_space);
215  H5Aclose(attribute_id);
216  }
217  if (mIsDataComplete)
218  {
220  mOffset = mLo;
222  }
223  else
224  {
225  // Read which nodes appear in the file (mIncompleteNodeIndices)
226  attribute_id = H5Aopen_name(mVariablesDatasetId, "NodeMap");
227  attribute_type = H5Aget_type(attribute_id);
228  attribute_space = H5Aget_space(attribute_id);
229 
230  // Get the dataset/dataspace dimensions
231  unsigned num_node_indices = H5Sget_simple_extent_npoints(attribute_space);
232 
233  // Read data from hyperslab in the file into the hyperslab in memory
234  mIncompleteNodeIndices.clear();
235  mIncompleteNodeIndices.resize(num_node_indices);
236  H5Aread(attribute_id, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);
237 
238  // Release ids
239  H5Tclose(attribute_type);
240  H5Sclose(attribute_space);
241  H5Aclose(attribute_id);
242 
243  // Set up what data we can
249  mDataFixedDimensionSize = UINT_MAX;
250  H5Dclose(mVariablesDatasetId);
251  H5Dclose(mUnlimitedDatasetId);
252  H5Fclose(mFileId);
253  EXCEPTION("Unable to extend an incomplete data file at present.");
254  }
255 
256  // Record chunk dimensions (useful for cached write mode)
257  hid_t dcpl = H5Dget_create_plist(mVariablesDatasetId); // get dataset creation property list
258  H5Pget_chunk(dcpl, DATASET_DIMS, mChunkSize );
259  if (mUseCache)
260  {
261  // Reserve space. Enough for one chunk in the time dimension.
263  }
264 
265  // Done
267  }
268  }
269 }
270 
272 {
273  Close();
274 
275  if (mSinglePermutation)
276  {
278  }
279  if (mDoublePermutation)
280  {
282  }
284  {
286  }
288  {
290  }
291 }
292 
294 {
295  OutputFileHandler output_file_handler(mDirectory, mCleanDirectory);
296  std::string file_name = mDirectory.GetAbsolutePath() + mBaseName + ".h5";
297 
298  if (mUseExistingFile)
299  {
300  FileFinder h5_file(file_name, RelativeTo::Absolute);
301  if (!h5_file.Exists())
302  {
303  EXCEPTION("Hdf5DataWriter could not open " + file_name + " , as it does not exist.");
304  }
305  }
306 
307  // Set up a property list saying how we'll open the file
308  hid_t fapl = H5Pcreate(H5P_FILE_ACCESS);
309  H5Pset_fapl_mpio(fapl, PETSC_COMM_WORLD, MPI_INFO_NULL);
310 
311  // Set size of each dimension in main dataset.
312  mDatasetDims[0] = mEstimatedUnlimitedLength; // While developing we got a non-documented "only the first dimension can be extendible" error.
313  mDatasetDims[1] = mFileFixedDimensionSize; // or should this be mDataFixedDimensionSize?
314  mDatasetDims[2] = mVariables.size();
315 
316  // Open the file and free the property list
317  std::string attempting_to;
318  if (mUseExistingFile)
319  {
320  mFileId = H5Fopen(file_name.c_str(), H5F_ACC_RDWR, fapl);
321  attempting_to = "open";
322  }
323  else
324  {
325  hid_t fcpl = H5Pcreate(H5P_FILE_CREATE);
326  /*
327  * Align objects to multiples of some number of bytes. Useful for
328  * striped filesystem e.g. Lustre. Ideally this should match the chunk
329  * size (see SetTargetChunkSize).
330  */
331  if (mAlignment != 0)
332  {
333  H5Pset_alignment(fapl, 0, mAlignment);
334  }
335 
336  /*
337  * The stripe size can be set on the directory just before creation of the H5
338  * file, and set back to normal just after, so that the H5 file inherits these
339  * settings, by uncommenting the lines below. Adjust the command for your
340  * specific filesystem!
341  *
342  * A simpler solution (where supported) might be to use an environment variable, e.g.
343  * export MPICH_MPIIO_HINTS="*.h5:striping_factor=48:striping_unit=1048576"
344  */
345  /*
346  std::string command;
347  // Turn on striping for directory, which the newly created HDF5 file will inherit.
348  if (PetscTools::AmMaster())
349  {
350  // Increase stripe count
351  command = "lfs setstripe --size 1M --count -1 "; // -1 means use all OSTs
352  command.append(mDirectory.GetAbsolutePath());
353  std::cout << command << std::endl;
354  system(command.c_str());
355  }
356  */
357 
358  // Create file
359  mFileId = H5Fcreate(file_name.c_str(), H5F_ACC_TRUNC, fcpl, fapl);
360 
361  /*
362  // Turn off striping so other output files stay unstriped.
363  if (PetscTools::AmMaster())
364  {
365  // Use one stripe
366  command = "lfs setstripe --size 1M --count 1 ";
367  command.append(mDirectory.GetAbsolutePath());
368  std::cout << command << std::endl;
369  system(command.c_str());
370  }
371  PetscTools::Barrier(); // Don't let other processes run away
372  */
373 
374  attempting_to = "create";
375  H5Pclose(fcpl);
376  }
377 
378  H5Pclose(fapl);
379 
380  if (mFileId < 0)
381  {
383  EXCEPTION("Hdf5DataWriter could not " << attempting_to << " " << file_name <<
384  " , H5F" << attempting_to << " error code = " << mFileId);
385  }
386 }
387 
388 void Hdf5DataWriter::DefineFixedDimension(long dimensionSize)
389 {
390  if (!mIsInDefineMode)
391  {
392  EXCEPTION("Cannot define variables when not in Define mode");
393  }
394  if (dimensionSize < 1)
395  {
396  EXCEPTION("Fixed dimension must be at least 1 long");
397  }
399  {
400  EXCEPTION("Fixed dimension already set");
401  }
402 
403  // Work out the ownership details
407  mOffset = mLo;
408  mFileFixedDimensionSize = dimensionSize;
409  mDataFixedDimensionSize = dimensionSize;
410  mIsFixedDimensionSet = true;
411 }
412 
413 void Hdf5DataWriter::DefineFixedDimension(const std::vector<unsigned>& rNodesToOuput, long vecSize)
414 {
415  unsigned vector_size = rNodesToOuput.size();
416 
417  for (unsigned index=0; index < vector_size-1; index++)
418  {
419  if (rNodesToOuput[index] >= rNodesToOuput[index+1])
420  {
421  EXCEPTION("Input should be monotonic increasing");
422  }
423  }
424 
425  if ((int) rNodesToOuput.back() >= vecSize)
426  {
427  EXCEPTION("Vector size doesn't match nodes to output");
428  }
429 
430  DefineFixedDimension(vecSize);
431 
432  mFileFixedDimensionSize = vector_size;
433  mIsDataComplete = false;
434  mIncompleteNodeIndices = rNodesToOuput;
436 }
437 
438 void Hdf5DataWriter::DefineFixedDimensionUsingMatrix(const std::vector<unsigned>& rNodesToOuput, long vecSize)
439 {
440  unsigned vector_size = rNodesToOuput.size();
441 
442  for (unsigned index=0; index < vector_size-1; index++)
443  {
444  if (rNodesToOuput[index] >= rNodesToOuput[index+1])
445  {
446  EXCEPTION("Input should be monotonic increasing");
447  }
448  }
449 
450  if ((int) rNodesToOuput.back() >= vecSize)
451  {
452  EXCEPTION("Vector size doesn't match nodes to output");
453  }
454 
455  DefineFixedDimension(vecSize);
456 
457  mFileFixedDimensionSize = vector_size;
458  mIsDataComplete = false;
459  mIncompleteNodeIndices = rNodesToOuput;
462 
463  // Make sure we've not done it already
464  assert(mSingleIncompleteOutputMatrix == nullptr);
465  assert(mDoubleIncompleteOutputMatrix == nullptr);
468 
469  // Only do local rows
470  for (unsigned row_index = mOffset; row_index < mOffset + mNumberOwned; row_index++)
471  {
472  // Put zero on the diagonal
473  MatSetValue(mSingleIncompleteOutputMatrix, row_index, row_index, 0.0, INSERT_VALUES);
474 
475  // Put one at (i,j)
476  MatSetValue(mSingleIncompleteOutputMatrix, row_index, rNodesToOuput[row_index], 1.0, INSERT_VALUES);
477 
478  unsigned bi_index = 2*row_index;
479  unsigned perm_index = 2*rNodesToOuput[row_index];
480 
481  // Put zeroes on the diagonal
482  MatSetValue(mDoubleIncompleteOutputMatrix, bi_index, bi_index, 0.0, INSERT_VALUES);
483  MatSetValue(mDoubleIncompleteOutputMatrix, bi_index+1, bi_index+1, 0.0, INSERT_VALUES);
484 
485  // Put ones at (i,j)
486  MatSetValue(mDoubleIncompleteOutputMatrix, bi_index, perm_index, 1.0, INSERT_VALUES);
487  MatSetValue(mDoubleIncompleteOutputMatrix, bi_index+1, perm_index+1, 1.0, INSERT_VALUES);
488  }
489  MatAssemblyBegin(mSingleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
490  MatAssemblyBegin(mDoubleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
491  MatAssemblyEnd(mSingleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
492  MatAssemblyEnd(mDoubleIncompleteOutputMatrix, MAT_FINAL_ASSEMBLY);
493 
494 // MatView(mSingleIncompleteOutputMatrix, PETSC_VIEWER_STDOUT_WORLD);
495 }
496 
498 {
499  mOffset = 0;
500  mNumberOwned = 0;
501  for (unsigned i=0; i<mIncompleteNodeIndices.size(); i++)
502  {
503  if (mIncompleteNodeIndices[i] < mLo)
504  {
505  mOffset++;
506  }
507  else if (mIncompleteNodeIndices[i] < mHi)
508  {
509  mNumberOwned++;
510  }
511  }
512 }
513 
514 int Hdf5DataWriter::DefineVariable(const std::string& rVariableName,
515  const std::string& rVariableUnits)
516 {
517  if (!mIsInDefineMode)
518  {
519  EXCEPTION("Cannot define variables when not in Define mode");
520  }
521 
522  CheckVariableName(rVariableName);
523  CheckUnitsName(rVariableUnits);
524 
525  // Check for the variable being already defined
526  for (unsigned index=0; index<mVariables.size(); index++)
527  {
528  if (mVariables[index].mVariableName == rVariableName)
529  {
530  EXCEPTION("Variable name already exists");
531  }
532  }
533 
534  DataWriterVariable new_variable;
535  new_variable.mVariableName = rVariableName;
536  new_variable.mVariableUnits = rVariableUnits;
537  int variable_id;
538 
539  // Add the variable to the variable vector
540  mVariables.push_back(new_variable);
541 
542  // Use the index of the variable vector as the variable ID.
543  // This is ok since there is no way to remove variables.
544  variable_id = mVariables.size() - 1;
545 
546  return variable_id;
547 }
548 
549 void Hdf5DataWriter::CheckVariableName(const std::string& rName)
550 {
551  if (rName.length() == 0)
552  {
553  EXCEPTION("Variable name not allowed: may not be blank.");
554  }
555  CheckUnitsName(rName);
556 }
557 
558 void Hdf5DataWriter::CheckUnitsName(const std::string& rName)
559 {
560  for (unsigned i=0; i<rName.length(); i++)
561  {
562  if (!isalnum(rName[i]) && !(rName[i]=='_'))
563  {
564  std::string error = "Variable name/units '" + rName + "' not allowed: may only contain alphanumeric characters or '_'.";
565  EXCEPTION(error);
566  }
567  }
568 }
569 
571 {
572  return mIsInDefineMode;
573 }
574 
576 {
577  // Check that at least one variable has been defined
578  if (mVariables.size() < 1)
579  {
580  EXCEPTION("Cannot end define mode. No variables have been defined.");
581  }
582 
583  // Check that a fixed dimension has been defined
584  if (mIsFixedDimensionSet == false)
585  {
586  EXCEPTION("Cannot end define mode. One fixed dimension should be defined.");
587  }
588 
589  OpenFile();
590 
591  mIsInDefineMode = false;
592 
593  /*
594  * Create "Data" dataset
595  */
596 
597  // Set max dims of dataset
598  hsize_t dataset_max_dims[DATASET_DIMS];
600  {
601  dataset_max_dims[0] = H5S_UNLIMITED;
602  }
603  else
604  {
605  dataset_max_dims[0] = 1;
606  }
607  dataset_max_dims[1] = mDatasetDims[1];
608  dataset_max_dims[2] = mDatasetDims[2];
609 
610  // Set chunk dimensions for the new dataset
611  SetChunkSize();
612 
613  if (mUseCache)
614  {
615  // Reserve space. Enough for one chunk in the time dimension.
617  }
618 
619  // Create chunked dataset and clean up
620  hid_t cparms = H5Pcreate (H5P_DATASET_CREATE);
621  H5Pset_chunk( cparms, DATASET_DIMS, mChunkSize);
622  hid_t filespace = H5Screate_simple(DATASET_DIMS, mDatasetDims, dataset_max_dims);
623  mVariablesDatasetId = H5Dcreate(mFileId, mDatasetName.c_str(), H5T_NATIVE_DOUBLE, filespace,
624  H5P_DEFAULT, cparms, H5P_DEFAULT);
625  SetMainDatasetRawChunkCache(); // Set large cache (even though parallel drivers don't currently use it!)
626  H5Sclose(filespace);
627  H5Pclose(cparms);
628 
629  // Create dataspace for the name, unit attribute
630  const unsigned MAX_STRING_SIZE = 100;
631  hsize_t columns[1] = {mVariables.size()};
632  hid_t colspace = H5Screate_simple(1, columns, nullptr);
633 
634  // Create attribute for variable names
635  char* col_data = (char*) malloc(mVariables.size() * sizeof(char) * MAX_STRING_SIZE);
636 
637  char* col_data_offset = col_data;
638  for (unsigned var=0; var<mVariables.size(); var++)
639  {
640  std::string full_name = mVariables[var].mVariableName + "(" + mVariables[var].mVariableUnits + ")";
641  strcpy (col_data_offset, full_name.c_str());
642  col_data_offset += sizeof(char) * MAX_STRING_SIZE;
643  }
644 
645  // Create the type 'string'
646  hid_t string_type = H5Tcopy(H5T_C_S1);
647  H5Tset_size(string_type, MAX_STRING_SIZE );
648  hid_t attr = H5Acreate(mVariablesDatasetId, "Variable Details", string_type, colspace,
649  H5P_DEFAULT, H5P_DEFAULT);
650 
651  // Write to the attribute
652  H5Awrite(attr, string_type, col_data);
653 
654  // Close dataspace & attribute
655  free(col_data);
656  H5Sclose(colspace);
657  H5Aclose(attr);
658 
659  // Create "boolean" attribute telling the data to be incomplete or not
660  columns[0] = 1;
661  colspace = H5Screate_simple(1, columns, nullptr);
662  attr = H5Acreate(mVariablesDatasetId, "IsDataComplete", H5T_NATIVE_UINT, colspace,
663  H5P_DEFAULT, H5P_DEFAULT);
664 
665  // Write to the attribute - note that the native boolean is not predictable
666  unsigned is_data_complete = mIsDataComplete ? 1 : 0;
667  H5Awrite(attr, H5T_NATIVE_UINT, &is_data_complete);
668 
669  H5Sclose(colspace);
670  H5Aclose(attr);
671 
672  if (!mIsDataComplete)
673  {
674  // We need to write a map
675  // Create "unsigned" attribute with the map
676  columns[0] = mFileFixedDimensionSize;
677  colspace = H5Screate_simple(1, columns, nullptr);
678  attr = H5Acreate(mVariablesDatasetId, "NodeMap", H5T_NATIVE_UINT, colspace,
679  H5P_DEFAULT, H5P_DEFAULT);
680 
681  // Write to the attribute
682  H5Awrite(attr, H5T_NATIVE_UINT, &mIncompleteNodeIndices[0]);
683 
684  H5Sclose(colspace);
685  H5Aclose(attr);
686  }
687 
688  /*
689  * Create "Time" dataset
690  */
692  {
693  hsize_t time_dataset_dims[1] = {mEstimatedUnlimitedLength};
694  hsize_t time_dataset_max_dims[1] = {H5S_UNLIMITED};
695 
696  /*
697  * Modify dataset creation properties to enable chunking. Set the chunk size in the "Time"
698  * dataset to 128 doubles, i.e. 1 KiB. See #2336.
699  */
700  hsize_t time_chunk_dims[1] = {128u};
701  hid_t time_cparms = H5Pcreate (H5P_DATASET_CREATE);
702  H5Pset_chunk( time_cparms, 1, time_chunk_dims);
703 
704  hid_t time_filespace = H5Screate_simple(1, time_dataset_dims, time_dataset_max_dims);
705 
706  // Create the unlimited dimension dataset
707 
708  // * Files post r18257 (inc. Release 3.2 onwards) use "<DatasetName>_Unlimited" for "<DatasetName>"'s
709  // unlimited variable,
710  // - a new attribute "Name" has been added to the Unlimited Dataset to allow it to assign
711  // any name to the unlimited variable. Which can then be easily read by Hdf5DataReader.
712 
713  mUnlimitedDatasetId = H5Dcreate(mFileId, (mDatasetName + "_Unlimited").c_str(), H5T_NATIVE_DOUBLE, time_filespace,
714  H5P_DEFAULT, time_cparms, H5P_DEFAULT);
715 
716  // Create the dataspace for the attribute
717  hsize_t one = 1;
718  hid_t one_column_space = H5Screate_simple(1, &one, nullptr);
719 
720  // Create an attribute for the time unit
721  hid_t unit_attr = H5Acreate(mUnlimitedDatasetId, "Unit", string_type, one_column_space,
722  H5P_DEFAULT, H5P_DEFAULT);
723 
724  // Create an attribute for the time name
725  hid_t name_attr = H5Acreate(mUnlimitedDatasetId, "Name", string_type, one_column_space,
726  H5P_DEFAULT, H5P_DEFAULT);
727 
728  // Copy the unit to a string MAX_STRING_SIZE long and write it
729  char unit_string[MAX_STRING_SIZE];
730  strcpy(unit_string, mUnlimitedDimensionUnit.c_str());
731  H5Awrite(unit_attr, string_type, unit_string);
732 
733  // Copy the unit to a string MAX_STRING_SIZE long and write it
734  char name_string[MAX_STRING_SIZE];
735  strcpy(name_string, mUnlimitedDimensionName.c_str());
736  H5Awrite(name_attr, string_type, name_string);
737 
738  // Close the filespace
739  H5Pclose(time_cparms);
740  H5Sclose(one_column_space);
741  H5Aclose(unit_attr);
742  H5Aclose(name_attr);
743  H5Sclose(time_filespace);
744  }
745 
746  /*
747  * Create the provenance attribute
748  */
749 
750  // Create a longer type for 'string'
751  const unsigned MAX_PROVENANCE_STRING_SIZE = 1023;
752  hid_t long_string_type = H5Tcopy(H5T_C_S1);
753  H5Tset_size(long_string_type, MAX_PROVENANCE_STRING_SIZE );
754  hsize_t prov_columns[1] = {1};
755  hid_t provenance_space = H5Screate_simple(1, prov_columns, nullptr);
756  char* provenance_data = (char*) malloc(sizeof(char) * MAX_PROVENANCE_STRING_SIZE);
757  assert( ChasteBuildInfo::GetProvenanceString().length() < MAX_PROVENANCE_STRING_SIZE);
758 
759  strcpy(provenance_data, ChasteBuildInfo::GetProvenanceString().c_str());
760  hid_t prov_attr = H5Acreate(mVariablesDatasetId, "Chaste Provenance", long_string_type, provenance_space,
761  H5P_DEFAULT, H5P_DEFAULT);
762 
763  // Write to the attribute
764  H5Awrite(prov_attr, long_string_type, provenance_data);
765 
766  // Close dataspace & attribute
767  free(provenance_data);
768  H5Sclose(provenance_space);
769  H5Aclose(prov_attr);
770 }
771 
772 void Hdf5DataWriter::PutVector(int variableID, Vec petscVector)
773 {
774  if (mIsInDefineMode)
775  {
776  EXCEPTION("Cannot write data while in define mode.");
777  }
778 
779  int vector_size;
780  VecGetSize(petscVector, &vector_size);
781 
782  if ((unsigned) vector_size != mDataFixedDimensionSize)
783  {
784  EXCEPTION("Vector size doesn't match fixed dimension");
785  }
786 
787  if (mUseCache)
788  {
789  if (mDatasetDims[2] != 1 )
790  {
791  //Covered by TestHdf5DataWriterMultipleColumnsCachedFails
792  EXCEPTION("Cached writes must write all variables at once.");
793  }
795  {
796  //Covered by TestHdf5DataWriterSingleColumnsCachedFails
797  EXCEPTION("Cached writes require an unlimited dimension.");
798  }
799  }
800 
801  // Make sure that everything is actually extended to the correct dimension.
802  PossiblyExtend();
803 
804  Vec output_petsc_vector;
805 
806  // Decide what to write
807  if (mSinglePermutation == nullptr)
808  {
809  // No permutation - just write
810  output_petsc_vector = petscVector;
811  }
812  else
813  {
814  assert(mIsDataComplete);
815  // Make a vector with the same pattern (doesn't copy the data)
816  VecDuplicate(petscVector, &output_petsc_vector);
817 
818  // Apply the permutation matrix
819  MatMult(mSinglePermutation, petscVector, output_petsc_vector);
820  }
821 
822  // Define memspace and hyperslab
823  hid_t memspace, hyperslab_space;
824  if (mNumberOwned != 0)
825  {
826  hsize_t v_size[1] = {mNumberOwned};
827  memspace = H5Screate_simple(1, v_size, nullptr);
828 
829  hsize_t count[DATASET_DIMS] = {1, mNumberOwned, 1};
830  hsize_t offset_dims[DATASET_DIMS] = {mCurrentTimeStep, mOffset, (unsigned)(variableID)};
831 
832  hyperslab_space = H5Dget_space(mVariablesDatasetId);
833  H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset_dims, nullptr, count, nullptr);
834  }
835  else
836  {
837  memspace = H5Screate(H5S_NULL);
838  hyperslab_space = H5Screate(H5S_NULL);
839  }
840 
841  // Create property list for collective dataset
842  hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
843  H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
844 
845  double* p_petsc_vector;
846  VecGetArray(output_petsc_vector, &p_petsc_vector);
847 
848  if (mIsDataComplete)
849  {
850  if (mUseCache)
851  {
852  //Covered by TestHdf5DataWriterSingleColumnCached
853  mDataCache.insert(mDataCache.end(), p_petsc_vector, p_petsc_vector+mNumberOwned);
854  }
855  else
856  {
857  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector);
858  }
859  }
860  else
861  {
863  {
864  // Make a vector of the required size
866 
867  // Fill the vector by multiplying complete data by incomplete output matrix
868  MatMult(mSingleIncompleteOutputMatrix, petscVector, output_petsc_vector);
869 
870  double* p_petsc_vector_incomplete;
871  VecGetArray(output_petsc_vector, &p_petsc_vector_incomplete);
872 
873  if (mUseCache)
874  {
875  //Covered by TestHdf5DataWriterSingleIncompleteUsingMatrixCached
876  mDataCache.insert(mDataCache.end(), p_petsc_vector_incomplete, p_petsc_vector_incomplete+mNumberOwned);
877  }
878  else
879  {
880  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector_incomplete);
881  }
882  }
883  else
884  {
885  // Make a local copy of the data you own
886  boost::scoped_array<double> local_data(new double[mNumberOwned]);
887  for (unsigned i=0; i<mNumberOwned; i++)
888  {
889  local_data[i] = p_petsc_vector[ mIncompleteNodeIndices[mOffset+i]-mLo ];
890 
891  }
892  if (mUseCache)
893  {
894  //Covered by TestHdf5DataWriterFullFormatIncompleteCached
895  mDataCache.insert(mDataCache.end(), local_data.get(), local_data.get()+mNumberOwned);
896  }
897  else
898  {
899  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, local_data.get());
900  }
901  }
902  }
903 
904  VecRestoreArray(output_petsc_vector, &p_petsc_vector);
905 
906  H5Sclose(memspace);
907  H5Sclose(hyperslab_space);
908  H5Pclose(property_list_id);
909 
910  if (petscVector != output_petsc_vector)
911  {
912  // Free local vector
913  PetscTools::Destroy(output_petsc_vector);
914  }
915 }
916 
917 void Hdf5DataWriter::PutStripedVector(std::vector<int> variableIDs, Vec petscVector)
918 {
919  if (mIsInDefineMode)
920  {
921  EXCEPTION("Cannot write data while in define mode.");
922  }
923 
924  if (variableIDs.size() <= 1)
925  {
926  EXCEPTION("The PutStripedVector method requires at least two variables ID. If only one is needed, use PutVector method instead");
927  }
928 
929  const unsigned NUM_STRIPES=variableIDs.size();
930 
931  if (mUseCache)
932  {
933  if (NUM_STRIPES != mDatasetDims[2])
934  {
935  //Covered by TestHdf5DataWriterFullFormatStripedIncompleteCached
936  EXCEPTION("Cached writes must write all variables at once.");
937  }
939  {
940  //Covered by TestHdf5DataWriterStripedNoTimeCachedFails
941  EXCEPTION("Cached writes require an unlimited dimension.");
942  }
943  }
944 
945  int firstVariableID=variableIDs[0];
946 
947  // Currently the method only works with consecutive columns, can be extended if needed
948  for (unsigned i = 1; i < variableIDs.size(); i++)
949  {
950  if (variableIDs[i]-variableIDs[i-1] != 1)
951  {
952  EXCEPTION("Columns should be consecutive. Try reordering them.");
953  }
954  }
955 
956  int vector_size;
957  VecGetSize(petscVector, &vector_size);
958 
959  if ((unsigned) vector_size != NUM_STRIPES*mDataFixedDimensionSize)
960  {
961  EXCEPTION("Vector size doesn't match fixed dimension");
962  }
963 
964  // Make sure that everything is actually extended to the correct dimension
965  PossiblyExtend();
966 
967  Vec output_petsc_vector;
968 
969  // Decide what to write
970  if (mDoublePermutation == nullptr)
971  {
972  // No permutation - just write
973  output_petsc_vector = petscVector;
974  }
975  else
976  {
977  assert(mIsDataComplete);
978  // Make a vector with the same pattern (doesn't copy the data)
979  VecDuplicate(petscVector, &output_petsc_vector);
980 
981  // Apply the permutation matrix
982  MatMult(mDoublePermutation, petscVector, output_petsc_vector);
983  }
984  // Define memspace and hyperslab
985  hid_t memspace, hyperslab_space;
986  if (mNumberOwned != 0)
987  {
988  hsize_t v_size[1] = {mNumberOwned*NUM_STRIPES};
989  memspace = H5Screate_simple(1, v_size, nullptr);
990 
991  hsize_t start[DATASET_DIMS] = {mCurrentTimeStep, mOffset, (unsigned)(firstVariableID)};
992  hsize_t stride[DATASET_DIMS] = {1, 1, 1};//we are imposing contiguous variables, hence the stride is 1 (3rd component)
993  hsize_t block_size[DATASET_DIMS] = {1, mNumberOwned, 1};
994  hsize_t number_blocks[DATASET_DIMS] = {1, 1, NUM_STRIPES};
995 
996  hyperslab_space = H5Dget_space(mVariablesDatasetId);
997  H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, start, stride, number_blocks, block_size);
998  }
999  else
1000  {
1001  memspace = H5Screate(H5S_NULL);
1002  hyperslab_space = H5Screate(H5S_NULL);
1003  }
1004 
1005  // Create property list for collective dataset write, and write! Finally.
1006  hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
1007  H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
1008 
1009  double* p_petsc_vector;
1010  VecGetArray(output_petsc_vector, &p_petsc_vector);
1011 
1012  if (mIsDataComplete)
1013  {
1014  if (mUseCache)
1015  {
1016  // Covered by TestHdf5DataWriterStripedCached
1017  mDataCache.insert(mDataCache.end(), p_petsc_vector, p_petsc_vector+mNumberOwned*NUM_STRIPES);
1018  }
1019  else
1020  {
1021  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector);
1022  }
1023  }
1024  else
1025  {
1026  if (variableIDs.size() < 3) // incomplete data and striped vector is supported only for NUM_STRIPES=2...for the moment
1027  {
1029  {
1030  // Make a vector of the required size
1031  output_petsc_vector = PetscTools::CreateVec(2*mFileFixedDimensionSize, 2*mNumberOwned);
1032 
1033  // Fill the vector by multiplying complete data by incomplete output matrix
1034  MatMult(mDoubleIncompleteOutputMatrix, petscVector, output_petsc_vector);
1035 
1036  double* p_petsc_vector_incomplete;
1037  VecGetArray(output_petsc_vector, &p_petsc_vector_incomplete);
1038 
1039  if (mUseCache)
1040  {
1041  //Covered by TestHdf5DataWriterFullFormatStripedIncompleteUsingMatrixCached
1042  mDataCache.insert(mDataCache.end(), p_petsc_vector_incomplete, p_petsc_vector_incomplete+2*mNumberOwned);
1043  }
1044  else
1045  {
1046  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, p_petsc_vector_incomplete);
1047  }
1048  }
1049  else
1050  {
1051  // Make a local copy of the data you own
1052  boost::scoped_array<double> local_data(new double[mNumberOwned*NUM_STRIPES]);
1053  for (unsigned i=0; i<mNumberOwned; i++)
1054  {
1055  unsigned local_node_number = mIncompleteNodeIndices[mOffset+i] - mLo;
1056  local_data[NUM_STRIPES*i] = p_petsc_vector[ local_node_number*NUM_STRIPES ];
1057  local_data[NUM_STRIPES*i+1] = p_petsc_vector[ local_node_number*NUM_STRIPES + 1];
1058  }
1059 
1060  if (mUseCache)
1061  {
1062  //Covered by TestHdf5DataWriterFullFormatStripedIncompleteCached
1063  mDataCache.insert(mDataCache.end(), local_data.get(), local_data.get()+2*mNumberOwned);
1064  }
1065  else
1066  {
1067  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, local_data.get());
1068  }
1069  }
1070  }
1071  else
1072  {
1073  EXCEPTION("The PutStripedVector functionality for incomplete data is supported for only 2 stripes");
1074  }
1075  }
1076 
1077  VecRestoreArray(output_petsc_vector, &p_petsc_vector);
1078 
1079  H5Sclose(memspace);
1080  H5Sclose(hyperslab_space);
1081  H5Pclose(property_list_id);
1082 
1083  if (petscVector != output_petsc_vector)
1084  {
1085  // Free local vector
1086  PetscTools::Destroy(output_petsc_vector);
1087  }
1088 }
1089 
1091 {
1092  return mUseCache;
1093 }
1094 
1096 {
1097  // The HDF5 writes are collective which means that if a process has nothing to write from
1098  // its cache then it must still proceed in step with the other processes.
1099  bool any_nonempty_caches = PetscTools::ReplicateBool( !mDataCache.empty() );
1100  if (!any_nonempty_caches)
1101  {
1102  // Nothing to do
1103  return;
1104  }
1105 
1106 // PRINT_3_VARIABLES(mCacheFirstTimeStep, mOffset, 0)
1107 // PRINT_3_VARIABLES(mCurrentTimeStep-mCacheFirstTimeStep, mNumberOwned, mDatasetDims[2])
1108 // PRINT_VARIABLE(mDataCache.size())
1109 
1110  // Define memspace and hyperslab
1111  hid_t memspace, hyperslab_space;
1112  if (mNumberOwned != 0)
1113  {
1114  hsize_t v_size[1] = {mDataCache.size()};
1115  memspace = H5Screate_simple(1, v_size, nullptr);
1116 
1119  assert((mCurrentTimeStep-mCacheFirstTimeStep)*mNumberOwned*mDatasetDims[2] == mDataCache.size()); // Got size right?
1120 
1121  hyperslab_space = H5Dget_space(mVariablesDatasetId);
1122  H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, start, nullptr, count, nullptr);
1123  }
1124  else
1125  {
1126  memspace = H5Screate(H5S_NULL);
1127  hyperslab_space = H5Screate(H5S_NULL);
1128  }
1129 
1130  // Create property list for collective dataset write
1131  hid_t property_list_id = H5Pcreate(H5P_DATASET_XFER);
1132  H5Pset_dxpl_mpio(property_list_id, H5FD_MPIO_COLLECTIVE);
1133 
1134  // Write!
1135  H5Dwrite(mVariablesDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, property_list_id, &mDataCache[0]);
1136 
1137  // Tidy up
1138  H5Sclose(memspace);
1139  H5Sclose(hyperslab_space);
1140  H5Pclose(property_list_id);
1141 
1142  mCacheFirstTimeStep = mCurrentTimeStep; // Update where we got to
1143  mDataCache.clear(); // Clear out cache
1144 }
1145 
1147 {
1148  if (mIsInDefineMode)
1149  {
1150  EXCEPTION("Cannot write data while in define mode.");
1151  }
1152 
1154  {
1155  EXCEPTION("PutUnlimitedVariable() called but no unlimited dimension has been set");
1156  }
1157 
1158  // Make sure that everything is actually extended to the correct dimension.
1159  PossiblyExtend();
1160 
1161  // This datum is only written by the master
1162  if (!PetscTools::AmMaster())
1163  {
1164  return;
1165  }
1166 
1167  hsize_t size[1] = {1};
1168  hid_t memspace = H5Screate_simple(1, size, nullptr);
1169 
1170  // Select hyperslab in the file.
1171  hsize_t count[1] = {1};
1172  hsize_t offset[1] = {mCurrentTimeStep};
1173  hid_t hyperslab_space = H5Dget_space(mUnlimitedDatasetId);
1174  H5Sselect_hyperslab(hyperslab_space, H5S_SELECT_SET, offset, nullptr, count, nullptr);
1175 
1176  H5Dwrite(mUnlimitedDatasetId, H5T_NATIVE_DOUBLE, memspace, hyperslab_space, H5P_DEFAULT, &value);
1177 
1178  H5Sclose(hyperslab_space);
1179  H5Sclose(memspace);
1180 }
1181 
1183 {
1184  if (mIsInDefineMode)
1185  {
1186  return; // Nothing to do...
1187  }
1188 
1189  if (mUseCache)
1190  {
1191  WriteCache();
1192  }
1193 
1194  H5Dclose(mVariablesDatasetId);
1196  {
1197  H5Dclose(mUnlimitedDatasetId);
1198  }
1199  H5Fclose(mFileId);
1200 
1201  // Cope with being called twice (e.g. if a user calls Close then the destructor)
1202  mIsInDefineMode = true;
1203 }
1204 
1205 void Hdf5DataWriter::DefineUnlimitedDimension(const std::string& rVariableName,
1206  const std::string& rVariableUnits,
1207  unsigned estimatedLength)
1208 {
1210  {
1211  EXCEPTION("Unlimited dimension already set. Cannot be defined twice");
1212  }
1213 
1214  if (!mIsInDefineMode)
1215  {
1216  EXCEPTION("Cannot define variables when not in Define mode");
1217  }
1218 
1219  this->mIsUnlimitedDimensionSet = true;
1220  this->mUnlimitedDimensionName = rVariableName;
1221  this->mUnlimitedDimensionUnit = rVariableUnits;
1222  mEstimatedUnlimitedLength = estimatedLength;
1223 }
1224 
1226 {
1228  {
1229  EXCEPTION("Trying to advance along an unlimited dimension without having defined any");
1230  }
1231 
1232  mCurrentTimeStep++;
1233 
1234  /*
1235  * Write when stepping over a chunk boundary. Note: NOT the same as write
1236  * out when the chunk size == the cache size, because we might have started
1237  * part-way through a chunk.
1238  */
1239  if (mUseCache && (mCurrentTimeStep % mChunkSize[0] == 0))
1240  {
1241  WriteCache();
1242  }
1243 
1244  /*
1245  * Extend the dataset (only reached when adding to an existing dataset,
1246  * or if mEstimatedUnlimitedLength hasn't been set and has defaulted to 1).
1247  */
1248  if (mCurrentTimeStep >= (long unsigned) mEstimatedUnlimitedLength)
1249  {
1250  mDatasetDims[0]++;
1251  mNeedExtend = true;
1252  }
1253 }
1254 
1256 {
1257  if (mNeedExtend)
1258  {
1259  H5Dset_extent( mVariablesDatasetId, mDatasetDims );
1260  H5Dset_extent( mUnlimitedDatasetId, mDatasetDims );
1261  }
1262  mNeedExtend = false;
1263 }
1264 
1266 {
1267  // Set internal counter to 0
1268  mCurrentTimeStep = 0;
1269  // Set dataset to 1 x nodes x vars
1270  mDatasetDims[0] = 1;
1271  mNeedExtend = 1;
1272  PossiblyExtend(); // Abusing the notation here, this is probably a contraction.
1273 }
1274 
1275 int Hdf5DataWriter::GetVariableByName(const std::string& rVariableName)
1276 {
1277  int id = -1;
1278 
1279  // Check for the variable name in the existing variables
1280  for (unsigned index=0; index<mVariables.size(); index++)
1281  {
1282  if (mVariables[index].mVariableName == rVariableName)
1283  {
1284  id = index;
1285  break;
1286  }
1287  }
1288  if (id == -1)
1289  {
1290  EXCEPTION("Variable does not exist in hdf5 definitions.");
1291  }
1292  return id;
1293 }
1294 
1295 bool Hdf5DataWriter::ApplyPermutation(const std::vector<unsigned>& rPermutation, bool unsafeExtendingMode)
1296 {
1297  if (unsafeExtendingMode == false && !mIsInDefineMode)
1298  {
1299  EXCEPTION("Cannot define permutation when not in Define mode");
1300  }
1301 
1302  if (rPermutation.empty())
1303  {
1304  return false;
1305  }
1306 
1307  if (rPermutation.size() != mFileFixedDimensionSize ||
1308  rPermutation.size() != mDataFixedDimensionSize)
1309  {
1310  EXCEPTION("Permutation doesn't match the expected problem size");
1311  }
1312 
1313  // Permutation checker
1314  std::set<unsigned> permutation_pigeon_hole;
1315 
1316  // Fill up the pigeon holes
1317  bool identity_map = true;
1318  for (unsigned i=0; i<mDataFixedDimensionSize; i++)
1319  {
1320  permutation_pigeon_hole.insert(rPermutation[i]);
1321  if (identity_map && i != rPermutation[i])
1322  {
1323  identity_map = false;
1324  }
1325  }
1326  if (identity_map)
1327  {
1328  // Do nothing
1329  return false;
1330  }
1331 
1332  /*
1333  * The pigeon-hole principle states that each index appears exactly once
1334  * so if any don't appear then either one appears twice or something out of
1335  * scope has appeared.
1336  */
1337  for (unsigned i=0; i<mDataFixedDimensionSize; i++)
1338  {
1339  if (permutation_pigeon_hole.count(i) != 1u)
1340  {
1341  EXCEPTION("Permutation vector doesn't contain a valid permutation");
1342  }
1343  }
1344  // Make sure we've not done it already
1345  assert(mSinglePermutation == nullptr);
1346  assert(mDoublePermutation == nullptr);
1347  PetscTools::SetupMat(mSinglePermutation, mDataFixedDimensionSize, mDataFixedDimensionSize, 2, mHi - mLo, mHi - mLo);
1348  PetscTools::SetupMat(mDoublePermutation, 2*mDataFixedDimensionSize, 2*mDataFixedDimensionSize, 4, 2*(mHi - mLo), 2*(mHi - mLo));
1349 #if (PETSC_VERSION_MAJOR == 3) //PETSc 3.x.x
1350  MatSetOption(mSinglePermutation, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
1351  MatSetOption(mDoublePermutation, MAT_IGNORE_OFF_PROC_ENTRIES, PETSC_TRUE);
1352 #else
1353  MatSetOption(mSinglePermutation, MAT_IGNORE_OFF_PROC_ENTRIES);
1354  MatSetOption(mDoublePermutation, MAT_IGNORE_OFF_PROC_ENTRIES);
1355 #endif
1356  // Only do local rows
1357  for (unsigned row_index=mLo; row_index<mHi; row_index++)
1358  {
1359  // Put zero on the diagonal
1360  MatSetValue(mSinglePermutation, row_index, row_index, 0.0, INSERT_VALUES);
1361 
1362  // Put one at (i,j)
1363  MatSetValue(mSinglePermutation, row_index, rPermutation[row_index], 1.0, INSERT_VALUES);
1364 
1365  unsigned bi_index = 2*row_index;
1366  unsigned perm_index = 2*rPermutation[row_index];
1367 
1368  // Put zeroes on the diagonal
1369  MatSetValue(mDoublePermutation, bi_index, bi_index, 0.0, INSERT_VALUES);
1370  MatSetValue(mDoublePermutation, bi_index+1, bi_index+1, 0.0, INSERT_VALUES);
1371 
1372  // Put ones at (i,j)
1373  MatSetValue(mDoublePermutation, bi_index, perm_index, 1.0, INSERT_VALUES);
1374  MatSetValue(mDoublePermutation, bi_index+1, perm_index+1, 1.0, INSERT_VALUES);
1375  }
1376  MatAssemblyBegin(mSinglePermutation, MAT_FINAL_ASSEMBLY);
1377  MatAssemblyBegin(mDoublePermutation, MAT_FINAL_ASSEMBLY);
1378  MatAssemblyEnd(mSinglePermutation, MAT_FINAL_ASSEMBLY);
1379  MatAssemblyEnd(mDoublePermutation, MAT_FINAL_ASSEMBLY);
1380  return true;
1381 }
1382 
1383 void Hdf5DataWriter::SetFixedChunkSize(const unsigned& rTimestepsPerChunk,
1384  const unsigned& rNodesPerChunk,
1385  const unsigned& rVariablesPerChunk)
1386 {
1387  assert(mIsInDefineMode);
1388 
1390  mFixedChunkSize[0] = rTimestepsPerChunk;
1391  mFixedChunkSize[1] = rNodesPerChunk;
1392  mFixedChunkSize[2] = rVariablesPerChunk;
1393 }
1394 
1396 {
1397  // Number of chunks for istore_k optimisation
1398  hsize_t num_chunks = 1;
1399  for (unsigned i=0; i<DATASET_DIMS; ++i)
1400  {
1401  num_chunks *= CeilDivide(mDatasetDims[i], mChunkSize[i]);
1402  }
1403  return num_chunks;
1404 }
1405 
1406 void Hdf5DataWriter::CalculateChunkDims( unsigned targetSize, unsigned* pChunkSizeInBytes, bool* pAllOneChunk )
1407 {
1408  bool all_one_chunk = true;
1409  unsigned chunk_size_in_bytes = 8u; // 8 bytes/double
1410  unsigned divisors[DATASET_DIMS];
1411  // Loop over dataset dimensions, dividing each dimension into the integer number of chunks that results
1412  // in the number of entries closest to the targetSize. This means the chunks will span the dataset with
1413  // as little waste as possible.
1414  for (unsigned i=0; i<DATASET_DIMS; ++i)
1415  {
1416  // What do I divide the dataset by to get targetSize entries per chunk?
1417  divisors[i] = CeilDivide(mDatasetDims[i], targetSize);
1418  // If I divide my dataset into divisors pieces, how big is each chunk?
1419  mChunkSize[i] = CeilDivide(mDatasetDims[i], divisors[i]);
1420  // If my chunks are this big, how many bytes is that?
1421  chunk_size_in_bytes *= mChunkSize[i];
1422  // Check if all divisors==1, which means we have one big chunk.
1423  all_one_chunk = all_one_chunk && divisors[i]==1u;
1424  }
1425  // Update pointers
1426  *pChunkSizeInBytes = chunk_size_in_bytes;
1427  *pAllOneChunk = all_one_chunk;
1428 }
1429 
1431 {
1433  {
1434  const unsigned target_size_in_bytes = mChunkTargetSize;
1435 
1436  unsigned target_size = 0;
1437  unsigned chunk_size_in_bytes;
1438  bool all_one_chunk;
1439 
1440  // While the chunks are too small, make mChunkSize[i]s larger, unless
1441  // we end up with a chunk that spans the dataset.
1442  do
1443  {
1444  target_size++;
1445  CalculateChunkDims(target_size, &chunk_size_in_bytes, &all_one_chunk);
1446  }
1447  while ( chunk_size_in_bytes < target_size_in_bytes && !all_one_chunk);
1448 
1449  // Go one step back if the target size has been exceeded
1450  if (chunk_size_in_bytes > target_size_in_bytes && !all_one_chunk)
1451  {
1452  target_size--;
1453  CalculateChunkDims(target_size, &chunk_size_in_bytes, &all_one_chunk);
1454  }
1455  }
1456  else
1457  {
1458  /*
1459  * User-provided chunk dims.
1460  */
1461  for (unsigned i=0; i<DATASET_DIMS; ++i)
1462  {
1463  mChunkSize[i] = mFixedChunkSize[i];
1464  }
1465  }
1466 
1468 
1469  /*
1470  if (PetscTools::AmMaster())
1471  {
1472  std::cout << "Hdf5DataWriter dataset contains " << mNumberOfChunks << " chunks of " << chunk_size_in_bytes << " B." << std::endl;
1473  }
1474  */
1475 }
1476 
1478 {
1479  if (!mIsInDefineMode)
1480  {
1481  /* Must be in define mode, i.e. creating a new dataset (and possibly a
1482  * new HDF5 file) to set the dataset chunk dims. */
1483  EXCEPTION("Cannot set chunk target size when not in define mode.");
1484  }
1485  mChunkTargetSize = targetSize;
1486 }
1487 
1489 {
1490  /* Note: calling this method after OpenFile() is pointless as that's where
1491  * H5Pset_alignment happens, so the obvious way to protect this method is
1492  * to check mFileId to assert H5Fopen has not been called yet.
1493  * Unfortunately it's uninitialised so is not always safe to check!
1494  *
1495  * Fortunately OpenFile() is only called in one of two ways:
1496  * 1. In the constructor in combination with mUseExistingFile. Subsequently
1497  * calling this method will end up throwing in the first block below
1498  * (since mUseExistingFile is const).
1499  * 2. By EndDefineMode(). Subsequently calling this method will end up in
1500  * the second block below.
1501  */
1502  if (mUseExistingFile)
1503  {
1504  // Must be called on new HDF5 files.
1505  EXCEPTION("Alignment parameter can only be set for new HDF5 files.");
1506  }
1507  if (!mIsInDefineMode)
1508  {
1509  // Creating a new file but EndDefineMode() called already.
1510  EXCEPTION("Cannot set alignment parameter when not in define mode.");
1511  }
1512 
1513  mAlignment = alignment;
1514 }
void ComputeIncompleteOffset()
hsize_t CalculateNumberOfChunks()
bool ApplyPermutation(const std::vector< unsigned > &rPermutation, bool unsafeExtendingMode=false)
bool mUseOptimalChunkSizeAlgorithm
std::string mUnlimitedDimensionName
static bool ReplicateBool(bool flag)
Definition: PetscTools.cpp:186
void CheckVariableName(const std::string &rName)
const bool mUseExistingFile
int GetVariableByName(const std::string &rVariableName)
unsigned CeilDivide(unsigned numerator, unsigned denominator)
hsize_t mChunkTargetSize
std::string GetAbsolutePath() const
Definition: FileFinder.cpp:224
std::string mUnlimitedDimensionUnit
#define EXCEPTION(message)
Definition: Exception.hpp:143
void CalculateChunkDims(unsigned targetSize, unsigned *pChunkSizeInBytes, bool *pAllOneChunk)
static bool AmMaster()
Definition: PetscTools.cpp:120
static Vec CreateVec(int size, int localSize=PETSC_DECIDE, bool ignoreOffProcEntries=true)
Definition: PetscTools.cpp:214
void DefineFixedDimensionUsingMatrix(const std::vector< unsigned > &rNodesToOuput, long vecSize)
void SetTargetChunkSize(hsize_t targetSize)
bool mUseMatrixForIncompleteData
bool DoesDatasetExist(const std::string &rDatasetName)
hsize_t mDatasetDims[DATASET_DIMS]
const bool mCleanDirectory
void SetFixedChunkSize(const unsigned &rTimestepsPerChunk, const unsigned &rNodesPerChunk, const unsigned &rVariablesPerChunk)
unsigned mEstimatedUnlimitedLength
long unsigned mCacheFirstTimeStep
void AdvanceAlongUnlimitedDimension()
std::vector< unsigned > mIncompleteNodeIndices
DistributedVectorFactory & mrVectorFactory
hsize_t mFixedChunkSize[DATASET_DIMS]
int DefineVariable(const std::string &rVariableName, const std::string &rVariableUnits)
static void SetupMat(Mat &rMat, int numRows, int numColumns, unsigned rowPreallocation, int numLocalRows=PETSC_DECIDE, int numLocalColumns=PETSC_DECIDE, bool ignoreOffProcEntries=true, bool newAllocationError=true)
Definition: PetscTools.cpp:268
Mat mDoubleIncompleteOutputMatrix
void PutVector(int variableID, Vec petscVector)
static void Destroy(Vec &rVec)
Definition: PetscTools.hpp:352
unsigned mNumberOwned
bool Exists() const
Definition: FileFinder.cpp:183
std::vector< double > mDataCache
void DefineFixedDimension(long dimensionSize)
Hdf5DataWriter(DistributedVectorFactory &rVectorFactory, const std::string &rDirectory, const std::string &rBaseName, bool cleanDirectory=true, bool extendData=false, std::string datasetName="Data", bool useCache=false)
void SetAlignment(hsize_t alignment)
static std::string GetProvenanceString()
virtual ~Hdf5DataWriter()
Mat mSingleIncompleteOutputMatrix
void CheckUnitsName(const std::string &rName)
void PutStripedVector(std::vector< int > variableIDs, Vec petscVector)
std::vector< DataWriterVariable > mVariables
void PutUnlimitedVariable(double value)
static const unsigned DATASET_DIMS
hsize_t mNumberOfChunks
unsigned mDataFixedDimensionSize
long unsigned mCurrentTimeStep
hsize_t mChunkSize[DATASET_DIMS]
unsigned mFileFixedDimensionSize
virtual void EndDefineMode()
void DefineUnlimitedDimension(const std::string &rVariableName, const std::string &rVariableUnits, unsigned estimatedLength=1)