util_h5file.hh

#ifndef DUNE_DORIE_UTIL_H5FILE_HH
#define DUNE_DORIE_UTIL_H5FILE_HH

#include <string>
#include <cassert>
#include <iostream>
#include <vector>
#include <exception>

#include <hdf5.h>

namespace Dune {
namespace Dorie {

/// A very simple HDF5 dataset reader. Represents an open H5 file.
class H5File
{

private:
    //! Verbosity of this class
    const int _verbose;
    //! ID of the H5 file
    hid_t _file_id;
    //! ID of the currently opened group
    hid_t _group_id;

public:
    /// Open the H5 file and its base group.
    /** \param file_path Path to the H5 file
     *  \param verbose Verbosity of the output
     */
    explicit H5File(const std::string& file_path, const int verbose):
        _verbose(verbose),
        _file_id(-1),
        _group_id(-1)
    {
        // set up property list for collective I/O
        hid_t h5_plist_id = H5Pcreate(H5P_FILE_ACCESS);
        H5Pset_fapl_mpio(h5_plist_id, MPI_COMM_WORLD, MPI_INFO_NULL);
        assert(h5_plist_id > -1);

        // open the file for reading
        _file_id = H5Fopen(file_path.c_str(), H5F_ACC_RDONLY, h5_plist_id);
        assert(_file_id > -1);

        // release properties
        herr_t status = H5Pclose(h5_plist_id);
        assert(status > -1);

        // check for errors regardless of build type
        if (_file_id < 0) {
            throw std::runtime_error("Cannot open H5 file " + file_path);
        }

        // open the base group
        open_group();
    }

    /// Close the H5 file when this object is destroyed
    ~H5File ()
    {
        // close the opened group
        herr_t status = H5Gclose(_group_id);
        assert(status > -1);

        // close the opened file
        status = H5Fclose(_file_id);
        assert(status > -1);
    }

    /// Open a group
    /** \param group_path Full internal path to the group.
     *      Defaults to the base group.
     */
    void open_group(const std::string& group_path="./")
    {
        // close the group if we already opened it
        if (_group_id >= 0) {
            herr_t status = H5Gclose(_group_id);
            assert(status > -1);
        }

        // open the new group
        _group_id = H5Gopen(_file_id, group_path.c_str(), H5P_DEFAULT);
        if (_group_id < 0) {
            throw std::runtime_error("Cannot open H5 group " + group_path);
        }
    }

    /// Read data from a multi-dimensional dataset.
    /** Data is read collectively from all parallel processes.
     *  The output data arrays are resized automatically.
     *
     *  \tparam data_t The data type to read the H5 array entries into.
     *      Must implicitly convert to the data H5 data type.
     *  \tparam ext_t The data type for array extensions.
     *      Must implicitly convert to hid_t.
     *  \param[in] dataset_path The full internal path of the dataset.
     *      May contain subgroups.
     *  \param[in] data_type The H5 internal data type of the array to read.
     *  \param[out] data The vector for reading the data into. Is resized
     *      automatically.
     *  \param[out] shape The vector containing the extensions of the dataset.
     *      This shape is the reverse information of the H5 internal shape.
     *      It indicated directions x, y[, z].
     */
    template<typename data_t, typename ext_t>
    void read_dataset(const std::string& dataset_path,
                      const hid_t data_type,
                      std::vector<data_t>& data,
                      std::vector<ext_t>& shape
                     )
    {
        if(_verbose > 1) {
            std::cout << "Reading H5 parameter field " << dataset_path
                      << std::endl;
        }

        // open a subgroup if the dataset path is nested
        const auto dataset_name = open_nested_dataset(dataset_path);

        // open the dataset
        hid_t dataset_id = H5Dopen(_group_id,
                                   dataset_name.c_str(),
                                   H5P_DEFAULT);
        if (dataset_id < 0) {
            throw std::runtime_error("Cannot open dataset " + dataset_path);
        }

        // get the dataspace
        hid_t dataspace_id = H5Dget_space(dataset_id);
        assert(dataspace_id > -1);

        // get the dimension (2-d or 3-d)
        const hsize_t arr_dim = H5Sget_simple_extent_ndims(dataspace_id);

        // get the size of the problem
        std::vector<hsize_t> dims(arr_dim, 0);
        herr_t status = H5Sget_simple_extent_dims(dataspace_id,
                                                  dims.data(),
                                                  0);
        assert(status > -1);

        const auto local_size = std::accumulate(dims.begin(),
                                                dims.end(),
                                                1,
                                                std::multiplies<hsize_t>());
        const std::vector<hsize_t> offset(arr_dim, 0);

        // reverse order of dimensions!
        shape.resize(arr_dim);
        std::copy(dims.rbegin(), dims.rend(), shape.begin());

        // create the memory space
        hid_t memspace_id = H5Screate_simple(arr_dim, dims.data(), NULL);

        //select the hyperslab
        status = H5Sselect_hyperslab(memspace_id, H5S_SELECT_SET, offset.data(), NULL, dims.data(), NULL);
        assert(status > -1);

        //resize the return data
        data.resize(local_size);

        // set up the collective transfer properties list
        hid_t h5_plist_id = H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(h5_plist_id, H5FD_MPIO_COLLECTIVE);
        assert(h5_plist_id > -1);

        // read the actual data
        status = H5Dread(dataset_id,
                         data_type,
                         memspace_id,
                         dataspace_id,
                         h5_plist_id,
                         data.data()
        );

        if (status < 0) {
            throw std::runtime_error("Error reading data from dataset " 
                                     + dataset_path);
        }

        // close the identifiers
        status = H5Dclose(dataset_id);
        status = H5Sclose(dataspace_id);
        status = H5Sclose(memspace_id);
        status = H5Pclose(h5_plist_id);
        assert(status > -1);
    }

private:

    /// Split the dataset path and open its group, if applicable
    /** \param dataset_path (Possibly nested) dataset path.
     *  \return The name of the actual dataset to open.
     */
    std::string open_nested_dataset(std::string dataset_path)
    {
        const auto pos = dataset_path.rfind('/');
        // split group from dataset name
        if (pos != std::string::npos) {
            const auto group = dataset_path.substr(0, pos);
            dataset_path = dataset_path.substr(pos+1);
            open_group(group);
        }

        return dataset_path;
    }
};

} // namespace Dorie
} // namespace Dune

#endif