13#include <dolfinx/common/log.h>
21namespace dolfinx::io::hdf5
27 if constexpr (std::is_same_v<T, float>)
28 return H5T_NATIVE_FLOAT;
29 else if constexpr (std::is_same_v<T, double>)
30 return H5T_NATIVE_DOUBLE;
31 else if constexpr (std::is_same_v<T, std::int32_t>)
32 return H5T_NATIVE_INT32;
33 else if constexpr (std::is_same_v<T, std::uint32_t>)
34 return H5T_NATIVE_UINT32;
35 else if constexpr (std::is_same_v<T, std::int64_t>)
36 return H5T_NATIVE_INT64;
37 else if constexpr (std::is_same_v<T, std::uint64_t>)
38 return H5T_NATIVE_UINT64;
39 else if constexpr (std::is_same_v<T, std::uint8_t>)
40 return H5T_NATIVE_UINT8;
41 else if constexpr (std::is_same_v<T, std::size_t>)
43 throw std::runtime_error(
44 "Cannot determine size of std::size_t. std::size_t is not the same "
45 "size as long or int.");
49 throw std::runtime_error(
"Cannot get HDF5 primitive data type. No "
50 "specialised function for this data type.");
59hid_t open_file(MPI_Comm comm,
const std::filesystem::path& filename,
60 const std::string& mode,
bool use_mpi_io);
64void close_file(hid_t handle);
68void flush_file(hid_t handle);
73std::filesystem::path get_filename(hid_t handle);
79bool has_dataset(hid_t handle,
const std::string& dataset_path);
85hid_t open_dataset(hid_t handle,
const std::string& path);
91std::vector<std::int64_t> get_dataset_shape(hid_t handle,
92 const std::string& dataset_path);
103void set_mpi_atomicity(hid_t handle,
bool atomic);
109bool get_mpi_atomicity(hid_t handle);
114void add_group(hid_t handle,
const std::string& dataset_path);
127void write_dataset(hid_t file_handle,
const std::string& dataset_path,
128 const T* data, std::array<std::int64_t, 2> range,
129 const std::vector<int64_t>& global_size,
bool use_mpi_io,
133 const int rank = global_size.size();
137 throw std::runtime_error(
"Cannot write dataset to HDF5 file"
138 "Only rank 1 and rank 2 dataset are supported");
142 const hid_t h5type = hdf5::hdf5_type<T>();
145 std::vector<hsize_t> count(global_size.begin(), global_size.end());
146 count[0] = range[1] - range[0];
149 std::vector<hsize_t> offset(rank, 0);
150 offset[0] = range[0];
153 const std::vector<hsize_t> dimsf(global_size.begin(), global_size.end());
156 const hid_t filespace0 = H5Screate_simple(rank, dimsf.data(),
nullptr);
157 if (filespace0 == H5I_INVALID_HID)
158 throw std::runtime_error(
"Failed to create HDF5 data space");
161 hid_t chunking_properties;
165 hsize_t chunk_size = dimsf[0] / 2;
166 if (chunk_size > 1048576)
167 chunk_size = 1048576;
168 if (chunk_size < 1024)
171 hsize_t chunk_dims[2] = {chunk_size, dimsf[1]};
172 chunking_properties = H5Pcreate(H5P_DATASET_CREATE);
173 H5Pset_chunk(chunking_properties, rank, chunk_dims);
176 chunking_properties = H5P_DEFAULT;
179 const std::string group_name(dataset_path, 0, dataset_path.rfind(
'/'));
180 add_group(file_handle, group_name);
184 = H5Dcreate2(file_handle, dataset_path.c_str(), h5type, filespace0,
185 H5P_DEFAULT, chunking_properties, H5P_DEFAULT);
186 if (dset_id == H5I_INVALID_HID)
187 throw std::runtime_error(
"Failed to create HDF5 global dataset.");
190 if (herr_t status = H5Sclose(filespace0); status < 0)
191 throw std::runtime_error(
"Failed to close HDF5 global data space.");
194 const hid_t memspace = H5Screate_simple(rank, count.data(),
nullptr);
195 if (memspace == H5I_INVALID_HID)
196 throw std::runtime_error(
"Failed to create HDF5 local data space.");
199 const hid_t filespace1 = H5Dget_space(dset_id);
200 herr_t status = H5Sselect_hyperslab(filespace1, H5S_SELECT_SET, offset.data(),
201 nullptr, count.data(),
nullptr);
203 throw std::runtime_error(
"Failed to create HDF5 dataspace.");
206 const hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
209 if (herr_t status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
212 throw std::runtime_error(
213 "Failed to set HDF5 data transfer property list.");
218 if (H5Dwrite(dset_id, h5type, memspace, filespace1, plist_id, data) < 0)
220 throw std::runtime_error(
221 "Failed to write HDF5 local dataset into hyperslab.");
227 if (H5Pclose(chunking_properties) < 0)
228 throw std::runtime_error(
"Failed to close HDF5 chunking properties.");
232 if (H5Dclose(dset_id) < 0)
233 throw std::runtime_error(
"Failed to close HDF5 dataset.");
236 if (H5Sclose(filespace1) < 0)
237 throw std::runtime_error(
"Failed to close HDF5 hyperslab.");
240 if (H5Sclose(memspace) < 0)
241 throw std::runtime_error(
"Failed to close local HDF5 dataset.");
244 if (H5Pclose(plist_id) < 0)
245 throw std::runtime_error(
"Failed to release HDF5 file-access template.");
258std::vector<T> read_dataset(hid_t dset_id, std::array<std::int64_t, 2> range,
261 auto timer_start = std::chrono::system_clock::now();
267 hid_t dtype = H5Dget_type(dset_id);
268 if (dtype == H5I_INVALID_HID)
269 throw std::runtime_error(
"Failed to get HDF5 data type.");
270 if (htri_t eq = H5Tequal(dtype, hdf5::hdf5_type<T>()); eq < 0)
271 throw std::runtime_error(
"HDF5 datatype equality test failed.");
275 throw std::runtime_error(
"Wrong type for reading from HDF5. Use \"h5ls "
276 "-v\" to inspect the types in the HDF5 file.");
281 hid_t dataspace = H5Dget_space(dset_id);
282 if (dataspace == H5I_INVALID_HID)
283 throw std::runtime_error(
"Failed to open HDF5 data space.");
286 int rank = H5Sget_simple_extent_ndims(dataspace);
288 throw std::runtime_error(
"Failed to get rank of data space.");
290 spdlog::warn(
"io::hdf5::read_dataset untested for rank > 2.");
293 std::vector<hsize_t> shape(rank);
296 if (
int ndims = H5Sget_simple_extent_dims(dataspace, shape.data(),
nullptr);
299 throw std::runtime_error(
"Failed to get dimensionality of dataspace.");
303 std::vector<hsize_t> offset(rank, 0);
304 std::vector<hsize_t> count = shape;
305 if (range[0] != -1 and range[1] != -1)
307 offset[0] = range[0];
308 count[0] = range[1] - range[0];
316 = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset.data(),
nullptr,
317 count.data(),
nullptr);
320 throw std::runtime_error(
"Failed to select HDF5 hyperslab.");
324 hid_t memspace = H5Screate_simple(rank, count.data(),
nullptr);
325 if (memspace == H5I_INVALID_HID)
326 throw std::runtime_error(
"Failed to create HDF5 dataspace.");
330 std::reduce(count.begin(), count.end(), 1, std::multiplies{}));
333 hid_t h5type = hdf5::hdf5_type<T>();
335 = H5Dread(dset_id, h5type, memspace, dataspace, H5P_DEFAULT, data.data());
338 throw std::runtime_error(
"Failed to read HDF5 data.");
342 if (herr_t status = H5Sclose(dataspace); status < 0)
343 throw std::runtime_error(
"Failed to close HDF5 dataspace.");
346 if (herr_t status = H5Sclose(memspace); status < 0)
347 throw std::runtime_error(
"Failed to close HDF5 memory space.");
349 auto timer_end = std::chrono::system_clock::now();
350 std::chrono::duration<double> dt = (timer_end - timer_start);
351 double data_rate = data.size() *
sizeof(T) / (1e6 * dt.count());
352 spdlog::info(
"HDF5 Read data rate: {} MB/s", data_rate);
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64