The TS-GitLab will have to shut down towards the end of the year — please think about migrating your projects to GitLab.com or GitHub.
(This is still a very early message, meant to keep you informed. There will be more urgent ones in the future.)

util_h5file.hh 6.77 KB
Newer Older
Dion Haefner's avatar
Dion Haefner committed
1 2 3
#ifndef DUNE_DORIE_UTIL_H5FILE_HH
#define DUNE_DORIE_UTIL_H5FILE_HH

4 5 6 7
#include <string>
#include <cassert>
#include <iostream>
#include <vector>
8
#include <exception>
9

Santiago Ospina's avatar
Santiago Ospina committed
10 11
#include <hdf5.h>

Dion Haefner's avatar
Dion Haefner committed
12
namespace Dune {
Lukas Riedel's avatar
Lukas Riedel committed
13 14 15 16 17 18 19
namespace Dorie {

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

private:
20
    //! Verbosity of this class
Lukas Riedel's avatar
Lukas Riedel committed
21
    const int _verbose;
22
    //! ID of the H5 file
Lukas Riedel's avatar
Lukas Riedel committed
23
    hid_t _file_id;
24
    //! ID of the currently opened group
Lukas Riedel's avatar
Lukas Riedel committed
25 26 27 28 29 30 31 32
    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):
Lukas Riedel's avatar
Lukas Riedel committed
33 34 35
        _verbose(verbose),
        _file_id(-1),
        _group_id(-1)
Dion Haefner's avatar
Dion Haefner committed
36
    {
37 38 39 40 41
        // 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);

Dion Haefner's avatar
Dion Haefner committed
42
        // open the file for reading
Lukas Riedel's avatar
Lukas Riedel committed
43 44
        _file_id = H5Fopen(file_path.c_str(), H5F_ACC_RDONLY, h5_plist_id);
        assert(_file_id > -1);
Dion Haefner's avatar
Dion Haefner committed
45

46
        // release properties
Lukas Riedel's avatar
Lukas Riedel committed
47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64
        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);
65

Lukas Riedel's avatar
Lukas Riedel committed
66 67 68 69 70 71 72 73 74 75 76
        // 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="./")
    {
Lukas Riedel's avatar
Lukas Riedel committed
77 78 79 80 81 82 83
        // 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
Lukas Riedel's avatar
Lukas Riedel committed
84 85 86 87 88 89
        _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);
        }
    }

90 91 92
    /// Read data from a multi-dimensional dataset.
    /** Data is read collectively from all parallel processes.
     *  The output data arrays are resized automatically.
Lukas Riedel's avatar
Lukas Riedel committed
93
     *
94 95 96 97 98 99 100 101 102 103 104 105
     *  \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].
Lukas Riedel's avatar
Lukas Riedel committed
106
     */
107 108 109 110 111 112
    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
                     )
Lukas Riedel's avatar
Lukas Riedel committed
113 114 115 116 117 118 119 120
    {
        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);
Dion Haefner's avatar
Dion Haefner committed
121 122

        // open the dataset
Lukas Riedel's avatar
Lukas Riedel committed
123 124 125
        hid_t dataset_id = H5Dopen(_group_id,
                                   dataset_name.c_str(),
                                   H5P_DEFAULT);
126 127 128
        if (dataset_id < 0) {
            throw std::runtime_error("Cannot open dataset " + dataset_path);
        }
Dion Haefner's avatar
Dion Haefner committed
129 130

        // get the dataspace
Lukas Riedel's avatar
Lukas Riedel committed
131 132
        hid_t dataspace_id = H5Dget_space(dataset_id);
        assert(dataspace_id > -1);
Dion Haefner's avatar
Dion Haefner committed
133 134

        // get the dimension (2-d or 3-d)
Lukas Riedel's avatar
Lukas Riedel committed
135
        const hsize_t arr_dim = H5Sget_simple_extent_ndims(dataspace_id);
Dion Haefner's avatar
Dion Haefner committed
136 137

        // get the size of the problem
138
        std::vector<hsize_t> dims(arr_dim, 0);
Lukas Riedel's avatar
Lukas Riedel committed
139 140 141
        herr_t status = H5Sget_simple_extent_dims(dataspace_id,
                                                  dims.data(),
                                                  0);
Dion Haefner's avatar
Dion Haefner committed
142 143
        assert(status > -1);

144 145 146
        const auto local_size = std::accumulate(dims.begin(),
                                                dims.end(),
                                                1,
147
                                                std::multiplies<hsize_t>());
148 149 150
        const std::vector<hsize_t> offset(arr_dim, 0);

        // reverse order of dimensions!
151 152
        shape.resize(arr_dim);
        std::copy(dims.rbegin(), dims.rend(), shape.begin());
Dion Haefner's avatar
Dion Haefner committed
153 154

        // create the memory space
155
        hid_t memspace_id = H5Screate_simple(arr_dim, dims.data(), NULL);
Dion Haefner's avatar
Dion Haefner committed
156 157

        //select the hyperslab
158
        status = H5Sselect_hyperslab(memspace_id, H5S_SELECT_SET, offset.data(), NULL, dims.data(), NULL);
159
        assert(status > -1);
Dion Haefner's avatar
Dion Haefner committed
160 161

        //resize the return data
162
        data.resize(local_size);
Dion Haefner's avatar
Dion Haefner committed
163 164

        // set up the collective transfer properties list
165 166 167
        hid_t h5_plist_id = H5Pcreate(H5P_DATASET_XFER);
        H5Pset_dxpl_mpio(h5_plist_id, H5FD_MPIO_COLLECTIVE);
        assert(h5_plist_id > -1);
Dion Haefner's avatar
Dion Haefner committed
168

169 170 171 172 173 174 175
        // read the actual data
        status = H5Dread(dataset_id,
                         data_type,
                         memspace_id,
                         dataspace_id,
                         h5_plist_id,
                         data.data()
Dion Haefner's avatar
Dion Haefner committed
176
        );
177 178 179 180 181

        if (status < 0) {
            throw std::runtime_error("Error reading data from dataset " 
                                     + dataset_path);
        }
Dion Haefner's avatar
Dion Haefner committed
182 183

        // close the identifiers
184 185 186 187 188
        status = H5Dclose(dataset_id);
        status = H5Sclose(dataspace_id);
        status = H5Sclose(memspace_id);
        status = H5Pclose(h5_plist_id);
        assert(status > -1);
Lukas Riedel's avatar
Lukas Riedel committed
189
    }
Dion Haefner's avatar
Dion Haefner committed
190

Lukas Riedel's avatar
Lukas Riedel committed
191
private:
Dion Haefner's avatar
Dion Haefner committed
192

Lukas Riedel's avatar
Lukas Riedel committed
193 194 195 196 197 198 199 200 201 202 203 204 205
    /// 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);
        }
Dion Haefner's avatar
Dion Haefner committed
206

Lukas Riedel's avatar
Lukas Riedel committed
207 208 209
        return dataset_path;
    }
};
Dion Haefner's avatar
Dion Haefner committed
210

Lukas Riedel's avatar
Lukas Riedel committed
211 212
} // namespace Dorie
} // namespace Dune
Dion Haefner's avatar
Dion Haefner committed
213 214

#endif