DOLFINx 0.9.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{
23
25template <typename T>
26hid_t hdf5_type()
27{
28 if constexpr (std::is_same_v<T, float>)
29 return H5T_NATIVE_FLOAT;
30 else if constexpr (std::is_same_v<T, double>)
31 return H5T_NATIVE_DOUBLE;
32 else if constexpr (std::is_same_v<T, std::int32_t>)
33 return H5T_NATIVE_INT32;
34 else if constexpr (std::is_same_v<T, std::uint32_t>)
35 return H5T_NATIVE_UINT32;
36 else if constexpr (std::is_same_v<T, std::int64_t>)
37 return H5T_NATIVE_INT64;
38 else if constexpr (std::is_same_v<T, std::uint64_t>)
39 return H5T_NATIVE_UINT64;
40 else if constexpr (std::is_same_v<T, std::size_t>)
41 {
42 throw std::runtime_error(
43 "Cannot determine size of std::size_t. std::size_t is not the same "
44 "size as long or int.");
45 }
46 else
47 {
48 throw std::runtime_error("Cannot get HDF5 primitive data type. No "
49 "specialised function for this data type.");
50 }
51}
52
58hid_t open_file(MPI_Comm comm, const std::filesystem::path& filename,
59 const std::string& mode, bool use_mpi_io);
60
63void close_file(hid_t handle);
64
67void flush_file(hid_t handle);
68
72std::filesystem::path get_filename(hid_t handle);
73
78bool has_dataset(hid_t handle, const std::string& dataset_path);
79
84hid_t open_dataset(hid_t handle, const std::string& path);
85
90std::vector<std::int64_t> get_dataset_shape(hid_t handle,
91 const std::string& dataset_path);
92
99void set_mpi_atomicity(hid_t handle, bool atomic);
100
105bool get_mpi_atomicity(hid_t handle);
106
110void add_group(hid_t handle, const std::string& dataset_path);
111
122template <typename T>
123void write_dataset(hid_t file_handle, const std::string& dataset_path,
124 const T* data, std::array<std::int64_t, 2> range,
125 const std::vector<int64_t>& global_size, bool use_mpi_io,
126 bool use_chunking)
127{
128 // Data rank
129 const int rank = global_size.size();
130 assert(rank != 0);
131 if (rank > 2)
132 {
133 throw std::runtime_error("Cannot write dataset to HDF5 file"
134 "Only rank 1 and rank 2 dataset are supported");
135 }
136
137 // Get HDF5 data type
138 const hid_t h5type = hdf5::hdf5_type<T>();
139
140 // Hyperslab selection parameters
141 std::vector<hsize_t> count(global_size.begin(), global_size.end());
142 count[0] = range[1] - range[0];
143
144 // Data offsets
145 std::vector<hsize_t> offset(rank, 0);
146 offset[0] = range[0];
147
148 // Dataset dimensions
149 const std::vector<hsize_t> dimsf(global_size.begin(), global_size.end());
150
151 // Create a global data space
152 const hid_t filespace0 = H5Screate_simple(rank, dimsf.data(), nullptr);
153 if (filespace0 == H5I_INVALID_HID)
154 throw std::runtime_error("Failed to create HDF5 data space");
155
156 // Set chunking parameters
157 hid_t chunking_properties;
158 if (use_chunking)
159 {
160 // Set chunk size and limit to 1kB min/1MB max
161 hsize_t chunk_size = dimsf[0] / 2;
162 if (chunk_size > 1048576)
163 chunk_size = 1048576;
164 if (chunk_size < 1024)
165 chunk_size = 1024;
166
167 hsize_t chunk_dims[2] = {chunk_size, dimsf[1]};
168 chunking_properties = H5Pcreate(H5P_DATASET_CREATE);
169 H5Pset_chunk(chunking_properties, rank, chunk_dims);
170 }
171 else
172 chunking_properties = H5P_DEFAULT;
173
174 // Check that group exists and recursively create if required
175 const std::string group_name(dataset_path, 0, dataset_path.rfind('/'));
176 add_group(file_handle, group_name);
177
178 // Create global dataset (using dataset_path)
179 const hid_t dset_id
180 = H5Dcreate2(file_handle, dataset_path.c_str(), h5type, filespace0,
181 H5P_DEFAULT, chunking_properties, H5P_DEFAULT);
182 if (dset_id == H5I_INVALID_HID)
183 throw std::runtime_error("Failed to create HDF5 global dataset.");
184
185 // Close global data space
186 if (herr_t status = H5Sclose(filespace0); status < 0)
187 throw std::runtime_error("Failed to close HDF5 global data space.");
188
189 // Create a local data space
190 const hid_t memspace = H5Screate_simple(rank, count.data(), nullptr);
191 if (memspace == H5I_INVALID_HID)
192 throw std::runtime_error("Failed to create HDF5 local data space.");
193
194 // Create a file dataspace within the global space - a hyperslab
195 const hid_t filespace1 = H5Dget_space(dset_id);
196 herr_t status = H5Sselect_hyperslab(filespace1, H5S_SELECT_SET, offset.data(),
197 nullptr, count.data(), nullptr);
198 if (status < 0)
199 throw std::runtime_error("Failed to create HDF5 dataspace.");
200
201 // Set parallel access
202 const hid_t plist_id = H5Pcreate(H5P_DATASET_XFER);
203 if (use_mpi_io)
204 {
205 if (herr_t status = H5Pset_dxpl_mpio(plist_id, H5FD_MPIO_COLLECTIVE);
206 status < 0)
207 {
208 throw std::runtime_error(
209 "Failed to set HDF5 data transfer property list.");
210 }
211 }
212
213 // Write local dataset into selected hyperslab
214 if (H5Dwrite(dset_id, h5type, memspace, filespace1, plist_id, data) < 0)
215 {
216 throw std::runtime_error(
217 "Failed to write HDF5 local dataset into hyperslab.");
218 }
219
220 if (use_chunking)
221 {
222 // Close chunking properties
223 if (H5Pclose(chunking_properties) < 0)
224 throw std::runtime_error("Failed to close HDF5 chunking properties.");
225 }
226
227 // Close dataset collectively
228 if (H5Dclose(dset_id) < 0)
229 throw std::runtime_error("Failed to close HDF5 dataset.");
230
231 // Close hyperslab
232 if (H5Sclose(filespace1) < 0)
233 throw std::runtime_error("Failed to close HDF5 hyperslab.");
234
235 // Close local dataset
236 if (H5Sclose(memspace) < 0)
237 throw std::runtime_error("Failed to close local HDF5 dataset.");
238
239 // Release file-access template
240 if (H5Pclose(plist_id) < 0)
241 throw std::runtime_error("Failed to release HDF5 file-access template.");
242}
243
253template <typename T>
254std::vector<T> read_dataset(hid_t dset_id, std::array<std::int64_t, 2> range,
255 bool allow_cast)
256{
257 auto timer_start = std::chrono::system_clock::now();
258
259 if (!allow_cast)
260 {
261 // Check that HDF5 dataset type and the type T are the same
262
263 hid_t dtype = H5Dget_type(dset_id);
264 if (dtype == H5I_INVALID_HID)
265 throw std::runtime_error("Failed to get HDF5 data type.");
266 if (htri_t eq = H5Tequal(dtype, hdf5::hdf5_type<T>()); eq < 0)
267 throw std::runtime_error("HDF5 datatype equality test failed.");
268 else if (!eq)
269 {
270 H5Tclose(dtype);
271 throw std::runtime_error("Wrong type for reading from HDF5. Use \"h5ls "
272 "-v\" to inspect the types in the HDF5 file.");
273 }
274 }
275
276 // Open dataspace
277 hid_t dataspace = H5Dget_space(dset_id);
278 if (dataspace == H5I_INVALID_HID)
279 throw std::runtime_error("Failed to open HDF5 data space.");
280
281 // Get rank of data set
282 int rank = H5Sget_simple_extent_ndims(dataspace);
283 if (rank < 1)
284 throw std::runtime_error("Failed to get rank of data space.");
285 else if (rank > 2)
286 spdlog::warn("io::hdf5::read_dataset untested for rank > 2.");
287
288 // Allocate data for shape
289 std::vector<hsize_t> shape(rank);
290
291 // Get size in each dimension
292 if (int ndims = H5Sget_simple_extent_dims(dataspace, shape.data(), nullptr);
293 ndims != rank)
294 {
295 throw std::runtime_error("Failed to get dimensionality of dataspace.");
296 }
297
298 // Hyperslab selection
299 std::vector<hsize_t> offset(rank, 0);
300 std::vector<hsize_t> count = shape;
301 if (range[0] != -1 and range[1] != -1)
302 {
303 offset[0] = range[0];
304 count[0] = range[1] - range[0];
305 }
306 else
307 offset[0] = 0;
308
309 // Select a block in the dataset beginning at offset[], with
310 // size=count[]
311 if (herr_t status
312 = H5Sselect_hyperslab(dataspace, H5S_SELECT_SET, offset.data(), nullptr,
313 count.data(), nullptr);
314 status < 0)
315 {
316 throw std::runtime_error("Failed to select HDF5 hyperslab.");
317 }
318
319 // Create a memory dataspace
320 hid_t memspace = H5Screate_simple(rank, count.data(), nullptr);
321 if (memspace == H5I_INVALID_HID)
322 throw std::runtime_error("Failed to create HDF5 dataspace.");
323
324 // Create local data to read into
325 std::vector<T> data(
326 std::reduce(count.begin(), count.end(), 1, std::multiplies{}));
327
328 // Read data on each process
329 hid_t h5type = hdf5::hdf5_type<T>();
330 if (herr_t status
331 = H5Dread(dset_id, h5type, memspace, dataspace, H5P_DEFAULT, data.data());
332 status < 0)
333 {
334 throw std::runtime_error("Failed to read HDF5 data.");
335 }
336
337 // Close dataspace
338 if (herr_t status = H5Sclose(dataspace); status < 0)
339 throw std::runtime_error("Failed to close HDF5 dataspace.");
340
341 // Close memspace
342 if (herr_t status = H5Sclose(memspace); status < 0)
343 throw std::runtime_error("Failed to close HDF5 memory space.");
344
345 auto timer_end = std::chrono::system_clock::now();
346 std::chrono::duration<double> dt = (timer_end - timer_start);
347 double data_rate = data.size() * sizeof(T) / (1e6 * dt.count());
348 spdlog::info("HDF5 Read data rate: {} MB/s", data_rate);
349
350 return data;
351}
352} // namespace dolfinx::io::hdf5
int rank(MPI_Comm comm)
Return process rank for the communicator.
Definition MPI.cpp:64