12 #include <dolfinx/common/log.h>
35 static hid_t
open_file(MPI_Comm mpi_comm,
const std::string& filename,
36 const std::string& mode,
const bool use_mpi_io);
62 static void write_dataset(
const hid_t handle,
const std::string& dataset_path,
64 const std::array<std::int64_t, 2>& range,
65 const std::vector<std::int64_t>& global_size,
66 bool use_mpi_io,
bool use_chunking);
78 const std::string& dataset_path,
79 const std::array<std::int64_t, 2>& range);
85 static bool has_dataset(
const hid_t handle,
const std::string& dataset_path);
91 static std::vector<std::int64_t>
112 static void add_group(
const hid_t handle,
const std::string& dataset_path);
115 template <
typename T>
116 static hid_t hdf5_type()
118 throw std::runtime_error(
"Cannot get HDF5 primitive data type. "
119 "No specialised function for this data type");
125 inline hid_t HDF5Interface::hdf5_type<float>()
127 return H5T_NATIVE_FLOAT;
131 inline hid_t HDF5Interface::hdf5_type<double>()
133 return H5T_NATIVE_DOUBLE;
137 inline hid_t HDF5Interface::hdf5_type<int>()
139 return H5T_NATIVE_INT;
143 inline hid_t HDF5Interface::hdf5_type<std::int64_t>()
145 return H5T_NATIVE_INT64;
149 inline hid_t HDF5Interface::hdf5_type<std::size_t>()
151 if (
sizeof(std::size_t) ==
sizeof(
unsigned long))
152 return H5T_NATIVE_ULONG;
153 else if (
sizeof(std::size_t) ==
sizeof(
unsigned int))
154 return H5T_NATIVE_UINT;
156 throw std::runtime_error(
"Cannot determine size of std::size_t. "
157 "std::size_t is not the same size as long or int");
160 template <
typename T>
162 const hid_t file_handle,
const std::string& dataset_path,
const T* data,
163 const std::array<std::int64_t, 2>& range,
164 const std::vector<int64_t>& global_size,
bool use_mpi_io,
bool use_chunking)
167 const std::size_t rank = global_size.size();
172 throw std::runtime_error(
"Cannot write dataset to HDF5 file"
173 "Only rank 1 and rank 2 dataset are supported");
177 const hid_t h5type = hdf5_type<T>();
180 std::vector<hsize_t> count(global_size.begin(), global_size.end());
181 count[0] = range[1] - range[0];
184 std::vector<hsize_t> offset(rank, 0);
185 offset[0] = range[0];
188 const std::vector<hsize_t> dimsf(global_size.begin(), global_size.end());
194 const hid_t filespace0 = H5Screate_simple(rank, dimsf.data(),
nullptr);
195 assert(filespace0 != HDF5_FAIL);
198 hid_t chunking_properties;
202 hsize_t chunk_size = dimsf[0] / 2;
203 if (chunk_size > 1048576)
204 chunk_size = 1048576;
205 if (chunk_size < 1024)
208 hsize_t chunk_dims[2] = {chunk_size, dimsf[1]};
209 chunking_properties = H5Pcreate(H5P_DATASET_CREATE);
210 H5Pset_chunk(chunking_properties, rank, chunk_dims);
213 chunking_properties = H5P_DEFAULT;
216 const std::string group_name(dataset_path, 0, dataset_path.rfind(
'/'));
217 add_group(file_handle, group_name);
221 = H5Dcreate2(file_handle, dataset_path.c_str(), h5type, filespace0,
222 H5P_DEFAULT, chunking_properties, H5P_DEFAULT);
223 assert(dset_id != HDF5_FAIL);
226 status = H5Sclose(filespace0);
227 assert(status != HDF5_FAIL);
230 const hid_t memspace = H5Screate_simple(rank, count.data(),
nullptr);
231 assert(memspace != HDF5_FAIL);
234 const hid_t filespace1 = H5Dget_space(dset_id);
235 status = H5Sselect_hyperslab(filespace1, H5S_SELECT_SET, offset.data(),
236 nullptr, count.data(),
nullptr);
237 assert(status != HDF5_FAIL);
240 const hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
243 #ifdef H5_HAVE_PARALLEL
244 status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
245 assert(status != HDF5_FAIL);
247 throw std::runtime_error(
"HDF5 library has not been configured with MPI");
252 status = H5Dwrite(dset_id, h5type, memspace, filespace1, plist_id, data);
253 assert(status != HDF5_FAIL);
258 status = H5Pclose(chunking_properties);
259 assert(status != HDF5_FAIL);
263 status = H5Dclose(dset_id);
264 assert(status != HDF5_FAIL);
267 status = H5Sclose(filespace1);
268 assert(status != HDF5_FAIL);
271 status = H5Sclose(memspace);
272 assert(status != HDF5_FAIL);
275 status = H5Pclose(plist_id);
276 assert(status != HDF5_FAIL);
279 template <
typename T>
280 inline std::vector<T>
282 const std::string& dataset_path,
283 const std::array<std::int64_t, 2>& range)
287 = H5Dopen2(file_handle, dataset_path.c_str(), H5P_DEFAULT);
288 assert(dset_id != HDF5_FAIL);
291 const hid_t dataspace = H5Dget_space(dset_id);
292 assert(dataspace != HDF5_FAIL);
295 const int rank = H5Sget_simple_extent_ndims(dataspace);
299 LOG(WARNING) <<
"HDF5Interface::read_dataset untested for rank > 2.";
302 std::vector<hsize_t> shape(rank);
305 const int ndims = H5Sget_simple_extent_dims(dataspace, shape.data(),
nullptr);
306 assert(ndims == rank);
309 std::vector<hsize_t> offset(rank, 0);
310 std::vector<hsize_t> count = shape;
311 if (range[0] != -1 and range[1] != -1)
313 offset[0] = range[0];
314 count[0] = range[1] - range[0];
321 herr_t status = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset.data(),
322 nullptr, count.data(),
nullptr);
323 assert(status != HDF5_FAIL);
326 const hid_t memspace = H5Screate_simple(rank, count.data(),
nullptr);
327 assert(memspace != HDF5_FAIL);
330 std::size_t data_size = 1;
331 for (std::size_t i = 0; i < count.size(); ++i)
332 data_size *= count[i];
333 std::vector<T> data(data_size);
336 const hid_t h5type = hdf5_type<T>();
338 = H5Dread(dset_id, h5type, memspace, dataspace, H5P_DEFAULT, data.data());
339 assert(status != HDF5_FAIL);
342 status = H5Sclose(dataspace);
343 assert(status != HDF5_FAIL);
346 status = H5Sclose(memspace);
347 assert(status != HDF5_FAIL);
350 status = H5Dclose(dset_id);
351 assert(status != HDF5_FAIL);