DOLFINx 0.10.0.0
DOLFINx C++ interface
Loading...
Searching...
No Matches
HDF5Interface.h
1// Copyright (C) 2012 Chris N. Richardson and Garth N. Wells
2//
3// This file is part of DOLFINx (https://www.fenicsproject.org)
4//
5// SPDX-License-Identifier: LGPL-3.0-or-later
6
7#pragma once
8
9#include <array>
10#include <cassert>
11#include <chrono>
12#include <cstdint>
13#include <dolfinx/common/log.h>
14#include <filesystem>
15#include <hdf5.h>
16#include <mpi.h>
17#include <numeric>
18#include <string>
19#include <vector>
20
21namespace dolfinx::io::hdf5
22{
24template <typename T>
25hid_t hdf5_type()
26{
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>)
42 {
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.");
46 }
47 else
48 {
49 throw std::runtime_error("Cannot get HDF5 primitive data type. No "
50 "specialised function for this data type.");
51 }
52}
53
59hid_t open_file(MPI_Comm comm, const std::filesystem::path& filename,
60 const std::string& mode, bool use_mpi_io);
61
64void close_file(hid_t handle);
65
68void flush_file(hid_t handle);
69
73std::filesystem::path get_filename(hid_t handle);
74
79bool has_dataset(hid_t handle, const std::string& dataset_path);
80
85hid_t open_dataset(hid_t handle, const std::string& path);
86
91std::vector<std::int64_t> get_dataset_shape(hid_t handle,
92 const std::string& dataset_path);
93
103void set_mpi_atomicity(hid_t handle, bool atomic);
104
109bool get_mpi_atomicity(hid_t handle);
110
114void add_group(hid_t handle, const std::string& dataset_path);
115
126template <typename T>
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,
130 bool use_chunking)
131{
132 // Data rank
133 const int rank = global_size.size();
134 assert(rank != 0);
135 if (rank > 2)
136 {
137 throw std::runtime_error("Cannot write dataset to HDF5 file"
138 "Only rank 1 and rank 2 dataset are supported");
139 }
140
141 // Get HDF5 data type
142 const hid_t h5type = hdf5::hdf5_type<T>();
143
144 // Hyperslab selection parameters
145 std::vector<hsize_t> count(global_size.begin(), global_size.end());
146 count[0] = range[1] - range[0];
147
148 // Data offsets
149 std::vector<hsize_t> offset(rank, 0);
150 offset[0] = range[0];
151
152 // Dataset dimensions
153 const std::vector<hsize_t> dimsf(global_size.begin(), global_size.end());
154
155 // Create a global data space
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");
159
160 // Set chunking parameters
161 hid_t chunking_properties;
162 if (use_chunking)
163 {
164 // Set chunk size and limit to 1kB min/1MB max
165 hsize_t chunk_size = dimsf[0] / 2;
166 if (chunk_size > 1048576)
167 chunk_size = 1048576;
168 if (chunk_size < 1024)
169 chunk_size = 1024;
170
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);
174 }
175 else
176 chunking_properties = H5P_DEFAULT;
177
178 // Check that group exists and recursively create if required
179 const std::string group_name(dataset_path, 0, dataset_path.rfind('/'));
180 add_group(file_handle, group_name);
181
182 // Create global dataset (using dataset_path)
183 const hid_t dset_id
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.");
188
189 // Close global data space
190 if (herr_t status = H5Sclose(filespace0); status < 0)
191 throw std::runtime_error("Failed to close HDF5 global data space.");
192
193 // Create a local 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.");
197
198 // Create a file dataspace within the global space - a hyperslab
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);
202 if (status < 0)
203 throw std::runtime_error("Failed to create HDF5 dataspace.");
204
205 // Set parallel access
206 const hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
207 if (use_mpi_io)
208 {
209 if (herr_t status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
210 status < 0)
211 {
212 throw std::runtime_error(
213 "Failed to set HDF5 data transfer property list.");
214 }
215 }
216
217 // Write local dataset into selected hyperslab
218 if (H5Dwrite(dset_id, h5type, memspace, filespace1, plist_id, data) < 0)
219 {
220 throw std::runtime_error(
221 "Failed to write HDF5 local dataset into hyperslab.");
222 }
223
224 if (use_chunking)
225 {
226 // Close chunking properties
227 if (H5Pclose(chunking_properties) < 0)
228 throw std::runtime_error("Failed to close HDF5 chunking properties.");
229 }
230
231 // Close dataset collectively
232 if (H5Dclose(dset_id) < 0)
233 throw std::runtime_error("Failed to close HDF5 dataset.");
234
235 // Close hyperslab
236 if (H5Sclose(filespace1) < 0)
237 throw std::runtime_error("Failed to close HDF5 hyperslab.");
238
239 // Close local dataset
240 if (H5Sclose(memspace) < 0)
241 throw std::runtime_error("Failed to close local HDF5 dataset.");
242
243 // Release file-access template
244 if (H5Pclose(plist_id) < 0)
245 throw std::runtime_error("Failed to release HDF5 file-access template.");
246}
247
257template <typename T>
258std::vector<T> read_dataset(hid_t dset_id, std::array<std::int64_t, 2> range,
259 bool allow_cast)
260{
261 auto timer_start = std::chrono::system_clock::now();
262
263 if (!allow_cast)
264 {
265 // Check that HDF5 dataset type and the type T are the same
266
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.");
272 else if (!eq)
273 {
274 H5Tclose(dtype);
275 throw std::runtime_error("Wrong type for reading from HDF5. Use \"h5ls "
276 "-v\" to inspect the types in the HDF5 file.");
277 }
278 }
279
280 // Open dataspace
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.");
284
285 // Get rank of data set
286 int rank = H5Sget_simple_extent_ndims(dataspace);
287 if (rank < 1)
288 throw std::runtime_error("Failed to get rank of data space.");
289 else if (rank > 2)
290 spdlog::warn("io::hdf5::read_dataset untested for rank > 2.");
291
292 // Allocate data for shape
293 std::vector<hsize_t> shape(rank);
294
295 // Get size in each dimension
296 if (int ndims = H5Sget_simple_extent_dims(dataspace, shape.data(), nullptr);
297 ndims != rank)
298 {
299 throw std::runtime_error("Failed to get dimensionality of dataspace.");
300 }
301
302 // Hyperslab selection
303 std::vector<hsize_t> offset(rank, 0);
304 std::vector<hsize_t> count = shape;
305 if (range[0] != -1 and range[1] != -1)
306 {
307 offset[0] = range[0];
308 count[0] = range[1] - range[0];
309 }
310 else
311 offset[0] = 0;
312
313 // Select a block in the dataset beginning at offset[], with
314 // size=count[]
315 if (herr_t status
316 = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset.data(), nullptr,
317 count.data(), nullptr);
318 status < 0)
319 {
320 throw std::runtime_error("Failed to select HDF5 hyperslab.");
321 }
322
323 // Create a memory dataspace
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.");
327
328 // Create local data to read into
329 std::vector<T> data(
330 std::reduce(count.begin(), count.end(), 1, std::multiplies{}));
331
332 // Read data on each process
333 hid_t h5type = hdf5::hdf5_type<T>();
334 if (herr_t status
335 = H5Dread(dset_id, h5type, memspace, dataspace, H5P_DEFAULT, data.data());
336 status < 0)
337 {
338 throw std::runtime_error("Failed to read HDF5 data.");
339 }
340
341 // Close dataspace
342 if (herr_t status = H5Sclose(dataspace); status < 0)
343 throw std::runtime_error("Failed to close HDF5 dataspace.");
344
345 // Close memspace
346 if (herr_t status = H5Sclose(memspace); status < 0)
347 throw std::runtime_error("Failed to close HDF5 memory space.");
348
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);
353
354 return data;
355}
356} // namespace dolfinx::io::hdf5
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64