sequenceio.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_HDF5_SEQUENCEIO_HH
#define DUNE_FUFEM_HDF5_SEQUENCEIO_HH
 
#include <array>
#include <functional>
#include <numeric>
#include <vector>
 
#if __has_include(<hdf5.h>)
 
#include <hdf5.h>
#include <hdf5_hl.h>
 
#include <dune/common/fvector.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/matrix.hh>
 
#include "file.hh"
#include "frombuffer.hh"
#include "tobuffer.hh"
#include "typetraits.hh"
 
namespace HDF5 {
template <int spatialDimensions, typename ctype = double, typename T = hsize_t>
class SequenceIO {
  int static const dimensions = 1 + spatialDimensions;
 
public:
  template <typename... Args>
  SequenceIO(Grouplike &file, std::string datasetname, Args... args)
      : file_(file), capacity_{{T(args)...}} {
    static_assert(sizeof...(args) == spatialDimensions,
                  "wrong number of arguments");
 
    if (file_.hasDataset(datasetname)) {
      dset_ = file_.openDataset(datasetname);
      checkDatatype();
      checkDimensions();
      checkChunkSize();
    } else {
      if (file_.isReadOnly())
        DUNE_THROW(Dune::Exception, "Dataset not found: " << datasetname);
 
      auto const initial_dims = maxExtent(0);
      auto const max_dims = maxExtent(H5S_UNLIMITED);
      hid_t file_space =
          H5Screate_simple(dimensions, initial_dims.data(), max_dims.data());
 
      hid_t plist = H5Pcreate(H5P_DATASET_CREATE);
      H5Pset_layout(plist, H5D_CHUNKED);
      auto const chunk_dims = maxExtent(1);
      H5Pset_chunk(plist, dimensions, chunk_dims.data());
 
      dset_ = file_.createDataset(datasetname, file_space, plist,
                                  TypeTraits<ctype>::getType());
 
      H5Pclose(plist);
      H5Sclose(file_space);
    }
  }
 
  ~SequenceIO() { H5Dclose(dset_); }
 
  void add(size_t timeStep, std::vector<ctype> const &buffer) {
    if (buffer.size() != entrySize())
      DUNE_THROW(Dune::Exception, "Buffer size incorrect");
 
    auto dims = getExtent();
    if (timeStep >= dims[0]) {
      dims[0] = timeStep + 1;
      H5Dset_extent(dset_, dims.data());
    }
    auto const start = minExtent(timeStep);
    H5DOwrite_chunk(dset_, H5P_DEFAULT, 0, start.data(),
                    entrySize() * sizeof(ctype), buffer.data());
  }
 
#if H5_VERSION_GE(1, 10, 0)
  void append(std::vector<ctype> const &buffer, size_t skip = 0) {
    if (buffer.size() != entrySize())
      DUNE_THROW(Dune::Exception, "Buffer size incorrect");
 
    H5DOappend(dset_, H5P_DEFAULT, skip, 1, TypeTraits<ctype>::getType(),
               buffer.data());
  }
#endif
 
  void read(size_t timeStep, std::vector<ctype> &buffer,
            std::array<T, spatialDimensions> &capacity) {
    buffer.resize(entrySize());
 
    hid_t file_space = H5Dget_space(dset_);
 
    std::array<T, dimensions> dims;
    H5Sget_simple_extent_dims(file_space, dims.data(), NULL);
    if (timeStep >= dims[0])
      DUNE_THROW(Dune::Exception, "Tried to read past final entry");
 
    auto const start = minExtent(timeStep);
    auto const count = maxExtent(1);
    H5Sselect_hyperslab(file_space, H5S_SELECT_SET, start.data(), NULL,
                        count.data(), NULL);
 
    hid_t mem_space = H5Screate_simple(dimensions, count.data(), NULL);
    H5Dread(dset_, TypeTraits<ctype>::getType(), mem_space, file_space,
            H5P_DEFAULT, buffer.data());
    H5Sclose(mem_space);
 
    H5Sclose(file_space);
    capacity = capacity_;
  }
 
private:
  std::array<T, dimensions> getExtent() const {
    hid_t file_space = H5Dget_space(dset_);
    std::array<T, dimensions> dims;
    H5Sget_simple_extent_dims(file_space, dims.data(), NULL);
    H5Sclose(file_space);
    return dims;
  }
 
  size_t entrySize() {
    return std::accumulate(capacity_.begin(), capacity_.end(), 1,
                           std::multiplies<T>());
  }
 
  void checkChunkSize() {
    hid_t plist = H5Dget_create_plist(dset_);
    H5D_layout_t layout = H5Pget_layout(plist);
    if (layout != H5D_CHUNKED)
      DUNE_THROW(Dune::Exception, "Layout not chunked");
 
    std::array<T, dimensions> chunk_dims;
    int const rank = H5Pget_chunk(plist, dimensions, chunk_dims.data());
    if (rank != dimensions)
      DUNE_THROW(Dune::Exception, "Unexpected chunk rank");
    if (!std::equal(capacity_.begin(), capacity_.end(), chunk_dims.begin() + 1))
      DUNE_THROW(Dune::Exception, "Unexpected chunk dimensions");
    if (chunk_dims[0] != 1)
      DUNE_THROW(Dune::Exception, "Unexpected chunk length");
 
    H5Pclose(plist);
  }
 
  void checkDimensions() {
    hid_t file_space = H5Dget_space(dset_);
    std::array<T, dimensions> dims;
    int const rank = H5Sget_simple_extent_dims(file_space, dims.data(), NULL);
 
    if (rank != dimensions)
      DUNE_THROW(Dune::Exception, "Unexpected dataset rank");
    if (!std::equal(capacity_.begin(), capacity_.end(), dims.begin() + 1))
      DUNE_THROW(Dune::Exception, "Unexpected dataset dimensions");
 
    H5Sclose(file_space);
  }
 
  void checkDatatype() {
    if (!H5Tequal(H5Dget_type(dset_), TypeTraits<ctype>::getType()))
      DUNE_THROW(Dune::Exception, "Unexpected data type");
  }
 
  std::array<T, dimensions> minExtent(T timeSteps) {
    std::array<T, dimensions> ret;
    ret[0] = timeSteps;
    std::fill(ret.begin() + 1, ret.end(), 0);
    return ret;
  }
 
  std::array<T, dimensions> maxExtent(T timeSteps) {
    std::array<T, dimensions> ret;
    ret[0] = timeSteps;
    std::copy(capacity_.begin(), capacity_.end(), ret.begin() + 1);
    return ret;
  }
 
  Grouplike &file_;
  std::array<T, spatialDimensions> const capacity_;
 
  hid_t dset_;
};
 
template <int spatialDimensions, typename ctype, typename T, class Data>
void addEntry(SequenceIO<spatialDimensions, ctype, T> &writer, size_t timeStep,
              Data const &data) {
  std::vector<ctype> buffer;
  toBuffer(data, buffer);
  writer.add(timeStep, buffer);
}
 
#if H5_VERSION_GE(1, 10, 0)
template <int spatialDimensions, typename ctype, typename T, class Data>
void appendEntry(SequenceIO<spatialDimensions, ctype, T> &writer,
                 Data const &data) {
  std::vector<ctype> buffer;
  toBuffer(data, buffer);
  writer.append(buffer);
}
#endif
 
template <int spatialDimensions, typename ctype, typename T, class Data>
void readEntry(SequenceIO<spatialDimensions, ctype, T> &reader, size_t timeStep,
               Data &data) {
  std::vector<ctype> buffer;
  std::array<T, spatialDimensions> dimensions;
  reader.read(timeStep, buffer, dimensions);
  fromBuffer(buffer, dimensions, data);
}
}
 
#else
  #warning Including the hdf5/sequenceio.hh but hdf5.h is missing.
#endif // __has_include(<hdf5.h>)
 
#endif