1 #ifndef OPENMC_HDF5_INTERFACE_H
2 #define OPENMC_HDF5_INTERFACE_H
3
4 #include <algorithm> // for min
5 #include <complex>
6 #include <cstddef>
7 #include <cstring> // for strlen
8 #include <string>
9 #include <sstream>
10 #include <type_traits>
11
12 #include "hdf5.h"
13 #include "hdf5_hl.h"
14 #include "xtensor/xadapt.hpp"
15 #include "xtensor/xarray.hpp"
16
17 #include "openmc/array.h"
18 #include "openmc/error.h"
19 #include "openmc/position.h"
20 #include "openmc/vector.h"
21
22 namespace openmc {
23
24 //==============================================================================
25 // Low-level internal functions
26 //==============================================================================
27
28 void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id,
29 void* buffer);
30
31 void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name,
32 hid_t mem_type_id, const void* buffer);
33
34 void read_dataset_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id,
35 hid_t mem_space_id, bool indep, void* buffer);
36
37 void write_dataset_lowlevel(hid_t group_id, int ndim, const hsize_t* dims,
38 const char* name, hid_t mem_type_id, hid_t mem_space_id, bool indep,
39 const void* buffer);
40
41 bool using_mpio_device(hid_t obj_id);
42
43 //==============================================================================
44 // Normal functions that are used to read/write files
45 //==============================================================================
46
47 hid_t create_group(hid_t parent_id, const std::string& name);
48
create_group(hid_t parent_id,const std::stringstream & name)49 inline hid_t create_group(hid_t parent_id, const std::stringstream& name)
50 {return create_group(parent_id, name.str());}
51
52
53 hid_t file_open(const std::string& filename, char mode, bool parallel=false);
54
55 void write_string(hid_t group_id, const char* name, const std::string& buffer,
56 bool indep);
57
58 vector<hsize_t> attribute_shape(hid_t obj_id, const char* name);
59 vector<std::string> dataset_names(hid_t group_id);
60 void ensure_exists(hid_t obj_id, const char* name, bool attribute=false);
61 vector<std::string> group_names(hid_t group_id);
62 vector<hsize_t> object_shape(hid_t obj_id);
63 std::string object_name(hid_t obj_id);
64
65 //==============================================================================
66 // Fortran compatibility functions
67 //==============================================================================
68
69 extern "C" {
70 bool attribute_exists(hid_t obj_id, const char* name);
71 size_t attribute_typesize(hid_t obj_id, const char* name);
72 hid_t create_group(hid_t parent_id, const char* name);
73 void close_dataset(hid_t dataset_id);
74 void close_group(hid_t group_id);
75 int dataset_ndims(hid_t dset);
76 size_t dataset_typesize(hid_t obj_id, const char* name);
77 hid_t file_open(const char* filename, char mode, bool parallel);
78 void file_close(hid_t file_id);
79 void get_name(hid_t obj_id, char* name);
80 int get_num_datasets(hid_t group_id);
81 int get_num_groups(hid_t group_id);
82 void get_datasets(hid_t group_id, char* name[]);
83 void get_groups(hid_t group_id, char* name[]);
84 void get_shape(hid_t obj_id, hsize_t* dims);
85 void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims);
86 bool object_exists(hid_t object_id, const char* name);
87 hid_t open_dataset(hid_t group_id, const char* name);
88 hid_t open_group(hid_t group_id, const char* name);
89 void read_attr_double(hid_t obj_id, const char* name, double* buffer);
90 void read_attr_int(hid_t obj_id, const char* name, int* buffer);
91 void read_attr_string(hid_t obj_id, const char* name, size_t slen,
92 char* buffer);
93 void read_complex(hid_t obj_id, const char* name,
94 std::complex<double>* buffer, bool indep);
95 void read_double(hid_t obj_id, const char* name, double* buffer, bool indep);
96 void read_int(hid_t obj_id, const char* name, int* buffer, bool indep);
97 void read_llong(hid_t obj_id, const char* name, long long* buffer,
98 bool indep);
99 void read_string(hid_t obj_id, const char* name, size_t slen, char* buffer,
100 bool indep);
101
102
103 void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
104 double* results);
105 void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims,
106 const char* name, const double* buffer);
107 void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims,
108 const char* name, const int* buffer);
109 void write_attr_string(hid_t obj_id, const char* name, const char* buffer);
110 void write_double(hid_t group_id, int ndim, const hsize_t* dims,
111 const char* name, const double* buffer, bool indep);
112 void write_int(hid_t group_id, int ndim, const hsize_t* dims,
113 const char* name, const int* buffer, bool indep);
114 void write_llong(hid_t group_id, int ndim, const hsize_t* dims,
115 const char* name, const long long* buffer, bool indep);
116 void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen,
117 const char* name, char const* buffer, bool indep);
118 void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score,
119 const double* results);
120 } // extern "C"
121
122 //==============================================================================
123 // Template struct used to map types to HDF5 datatype IDs, which are stored
124 // using the type hid_t. By having a single static data member, the template can
125 // be specialized for each type we know of. The specializations appear in the
126 // .cpp file since they are definitions.
127 //==============================================================================
128
129 template<typename T>
130 struct H5TypeMap { static const hid_t type_id; };
131
132 //==============================================================================
133 // Templates/overloads for read_attribute
134 //==============================================================================
135
136 // Scalar version
137 template<typename T>
read_attribute(hid_t obj_id,const char * name,T & buffer)138 void read_attribute(hid_t obj_id, const char* name, T& buffer)
139 {
140 read_attr(obj_id, name, H5TypeMap<T>::type_id, &buffer);
141 }
142
143 // array version
144 template<typename T, std::size_t N>
read_attribute(hid_t obj_id,const char * name,array<T,N> & buffer)145 inline void read_attribute(hid_t obj_id, const char* name, array<T, N>& buffer)
146 {
147 read_attr(obj_id, name, H5TypeMap<T>::type_id, buffer.data());
148 }
149
150 // vector version
151 template<typename T>
read_attribute(hid_t obj_id,const char * name,vector<T> & vec)152 void read_attribute(hid_t obj_id, const char* name, vector<T>& vec)
153 {
154 // Get shape of attribute array
155 auto shape = attribute_shape(obj_id, name);
156
157 // Allocate new array to read data into
158 std::size_t size = 1;
159 for (const auto x : shape)
160 size *= x;
161 vec.resize(size);
162
163 // Read data from attribute
164 read_attr(obj_id, name, H5TypeMap<T>::type_id, vec.data());
165 }
166
167 // Generic array version
168 template<typename T>
read_attribute(hid_t obj_id,const char * name,xt::xarray<T> & arr)169 void read_attribute(hid_t obj_id, const char* name, xt::xarray<T>& arr)
170 {
171 // Get shape of attribute array
172 auto shape = attribute_shape(obj_id, name);
173
174 // Allocate new array to read data into
175 std::size_t size = 1;
176 for (const auto x : shape)
177 size *= x;
178 vector<T> buffer(size);
179
180 // Read data from attribute
181 read_attr(obj_id, name, H5TypeMap<T>::type_id, buffer.data());
182
183 // Adapt array into xarray
184 arr = xt::adapt(buffer, shape);
185 }
186
187 // overload for std::string
188 inline void
read_attribute(hid_t obj_id,const char * name,std::string & str)189 read_attribute(hid_t obj_id, const char* name, std::string& str)
190 {
191 // Create buffer to read data into
192 auto n = attribute_typesize(obj_id, name);
193 char* buffer = new char[n];
194
195 // Read attribute and set string
196 read_attr_string(obj_id, name, n, buffer);
197 str = std::string{buffer, n};
198 delete[] buffer;
199 }
200
201 // overload for vector<std::string>
read_attribute(hid_t obj_id,const char * name,vector<std::string> & vec)202 inline void read_attribute(
203 hid_t obj_id, const char* name, vector<std::string>& vec)
204 {
205 auto dims = attribute_shape(obj_id, name);
206 auto m = dims[0];
207
208 // Allocate a C char array to get strings
209 auto n = attribute_typesize(obj_id, name);
210 char* buffer = new char[m*n];
211
212 // Read char data in attribute
213 read_attr_string(obj_id, name, n, buffer);
214
215 for (decltype(m) i = 0; i < m; ++i) {
216 // Determine proper length of string -- strlen doesn't work because
217 // buffer[i] might not have any null characters
218 std::size_t k = 0;
219 for (; k < n; ++k) if (buffer[i*n + k] == '\0') break;
220
221 // Create string based on (char*, size_t) constructor
222 vec.emplace_back(&buffer[i*n], k);
223 }
224 delete[] buffer;
225 }
226
227 //==============================================================================
228 // Templates/overloads for read_dataset and related methods
229 //==============================================================================
230
231 // Template for scalars. We need to be careful that the compiler does not use
232 // this version of read_dataset for vectors, arrays, or other non-scalar types.
233 // enable_if_t allows us to conditionally remove the function from overload
234 // resolution when the type T doesn't meet a certain criterion.
235 template<typename T> inline
236 std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
237 read_dataset(hid_t obj_id, const char* name, T& buffer, bool indep=false)
238 {
239 read_dataset_lowlevel(obj_id, name, H5TypeMap<T>::type_id, H5S_ALL, indep,
240 &buffer);
241 }
242
243 // overload for std::string
244 inline void
245 read_dataset(hid_t obj_id, const char* name, std::string& str, bool indep=false)
246 {
247 // Create buffer to read data into
248 auto n = dataset_typesize(obj_id, name);
249 char* buffer = new char[n];
250
251 // Read attribute and set string
252 read_string(obj_id, name, n, buffer, indep);
253 str = std::string{buffer, n};
254 }
255
256 // array version
257 template<typename T, std::size_t N>
258 inline void read_dataset(
259 hid_t dset, const char* name, array<T, N>& buffer, bool indep = false)
260 {
261 read_dataset_lowlevel(dset, name, H5TypeMap<T>::type_id, H5S_ALL, indep,
262 buffer.data());
263 }
264
265 // vector version
266 template<typename T>
267 void read_dataset(hid_t dset, vector<T>& vec, bool indep = false)
268 {
269 // Get shape of dataset
270 vector<hsize_t> shape = object_shape(dset);
271
272 // Resize vector to appropriate size
273 vec.resize(shape[0]);
274
275 // Read data into vector
276 read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
277 vec.data());
278 }
279
280 template<typename T>
281 void read_dataset(
282 hid_t obj_id, const char* name, vector<T>& vec, bool indep = false)
283 {
284 hid_t dset = open_dataset(obj_id, name);
285 read_dataset(dset, vec, indep);
286 close_dataset(dset);
287 }
288
289 template <typename T>
290 void read_dataset(hid_t dset, xt::xarray<T>& arr, bool indep=false)
291 {
292 // Get shape of dataset
293 vector<hsize_t> shape = object_shape(dset);
294
295 // Allocate space in the array to read data into
296 std::size_t size = 1;
297 for (const auto x : shape)
298 size *= x;
299 arr.resize(shape);
300
301 // Read data from attribute
302 read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
303 arr.data());
304 }
305
306 template<>
307 void read_dataset(hid_t dset, xt::xarray<std::complex<double>>& arr,
308 bool indep);
309
310 template <typename T>
311 void read_dataset(hid_t obj_id, const char* name, xt::xarray<T>& arr,
312 bool indep=false)
313 {
314 // Open dataset and read array
315 hid_t dset = open_dataset(obj_id, name);
316 read_dataset(dset, arr, indep);
317 close_dataset(dset);
318 }
319
320
321 template <typename T, std::size_t N>
322 void read_dataset(hid_t obj_id, const char* name, xt::xtensor<T, N>& arr,
323 bool indep=false)
324 {
325 // Open dataset and read array
326 hid_t dset = open_dataset(obj_id, name);
327
328 // Get shape of dataset
329 vector<hsize_t> hsize_t_shape = object_shape(dset);
330 close_dataset(dset);
331
332 // cast from hsize_t to size_t
333 vector<size_t> shape(hsize_t_shape.size());
334 for (int i = 0; i < shape.size(); i++) {
335 shape[i] = static_cast<size_t>(hsize_t_shape[i]);
336 }
337
338 // Allocate new xarray to read data into
339 xt::xarray<T> xarr(shape);
340
341 // Read data from the dataset
342 read_dataset(obj_id, name, xarr);
343
344 // Copy into xtensor
345 arr = xarr;
346 }
347
348 // overload for Position
349 inline void
350 read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false)
351 {
352 array<double, 3> x;
353 read_dataset(obj_id, name, x, indep);
354 r.x = x[0];
355 r.y = x[1];
356 r.z = x[2];
357 }
358
359 template <typename T, std::size_t N>
360 inline void read_dataset_as_shape(hid_t obj_id, const char* name,
361 xt::xtensor<T, N>& arr, bool indep=false)
362 {
363 hid_t dset = open_dataset(obj_id, name);
364
365 // Allocate new array to read data into
366 std::size_t size = 1;
367 for (const auto x : arr.shape())
368 size *= x;
369 vector<T> buffer(size);
370
371 // Read data from attribute
372 read_dataset_lowlevel(dset, nullptr, H5TypeMap<T>::type_id, H5S_ALL, indep,
373 buffer.data());
374
375 // Adapt into xarray
376 arr = xt::adapt(buffer, arr.shape());
377
378 close_dataset(dset);
379 }
380
381
382 template <typename T, std::size_t N>
383 inline void read_nd_vector(hid_t obj_id, const char* name,
384 xt::xtensor<T, N>& result, bool must_have=false)
385 {
386 if (object_exists(obj_id, name)) {
387 read_dataset_as_shape(obj_id, name, result, true);
388 } else if (must_have) {
389 fatal_error(std::string("Must provide " + std::string(name) + "!"));
390 }
391 }
392
393 //==============================================================================
394 // Templates/overloads for write_attribute
395 //==============================================================================
396
397 template<typename T> inline void
write_attribute(hid_t obj_id,const char * name,T buffer)398 write_attribute(hid_t obj_id, const char* name, T buffer)
399 {
400 write_attr(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id, &buffer);
401 }
402
403 inline void
write_attribute(hid_t obj_id,const char * name,const char * buffer)404 write_attribute(hid_t obj_id, const char* name, const char* buffer)
405 {
406 write_attr_string(obj_id, name, buffer);
407 }
408
409 inline void
write_attribute(hid_t obj_id,const char * name,const std::string & buffer)410 write_attribute(hid_t obj_id, const char* name, const std::string& buffer)
411 {
412 write_attr_string(obj_id, name, buffer.c_str());
413 }
414
415 template<typename T, std::size_t N>
write_attribute(hid_t obj_id,const char * name,const array<T,N> & buffer)416 inline void write_attribute(
417 hid_t obj_id, const char* name, const array<T, N>& buffer)
418 {
419 hsize_t dims[] {N};
420 write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
421 }
422
423 template<typename T>
write_attribute(hid_t obj_id,const char * name,const vector<T> & buffer)424 inline void write_attribute(
425 hid_t obj_id, const char* name, const vector<T>& buffer)
426 {
427 hsize_t dims[] {buffer.size()};
428 write_attr(obj_id, 1, dims, name, H5TypeMap<T>::type_id, buffer.data());
429 }
430
431 inline void
write_attribute(hid_t obj_id,const char * name,Position r)432 write_attribute(hid_t obj_id, const char* name, Position r)
433 {
434 array<double, 3> buffer {r.x, r.y, r.z};
435 write_attribute(obj_id, name, buffer);
436 }
437
438
439
440 //==============================================================================
441 // Templates/overloads for write_dataset
442 //==============================================================================
443
444 // Template for scalars (ensured by SFINAE)
445 template<typename T> inline
446 std::enable_if_t<std::is_scalar<std::decay_t<T>>::value>
write_dataset(hid_t obj_id,const char * name,T buffer)447 write_dataset(hid_t obj_id, const char* name, T buffer)
448 {
449 write_dataset_lowlevel(obj_id, 0, nullptr, name, H5TypeMap<T>::type_id,
450 H5S_ALL, false, &buffer);
451 }
452
453 inline void
write_dataset(hid_t obj_id,const char * name,const char * buffer)454 write_dataset(hid_t obj_id, const char* name, const char* buffer)
455 {
456 write_string(obj_id, name, buffer, false);
457 }
458
459 template<typename T, std::size_t N>
write_dataset(hid_t obj_id,const char * name,const array<T,N> & buffer)460 inline void write_dataset(
461 hid_t obj_id, const char* name, const array<T, N>& buffer)
462 {
463 hsize_t dims[] {N};
464 write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap<T>::type_id,
465 H5S_ALL, false, buffer.data());
466 }
467
write_dataset(hid_t obj_id,const char * name,const vector<std::string> & buffer)468 inline void write_dataset(
469 hid_t obj_id, const char* name, const vector<std::string>& buffer)
470 {
471 auto n {buffer.size()};
472 hsize_t dims[] {n};
473
474 // Determine length of longest string, including \0
475 size_t m = 1;
476 for (const auto& s : buffer) {
477 m = std::max(m, s.size() + 1);
478 }
479
480 // Copy data into contiguous buffer
481 char* temp = new char[n*m];
482 std::fill(temp, temp + n*m, '\0');
483 for (decltype(n) i = 0; i < n; ++i) {
484 std::copy(buffer[i].begin(), buffer[i].end(), temp + i*m);
485 }
486
487 // Write 2D data
488 write_string(obj_id, 1, dims, m, name, temp, false);
489
490 // Free temp array
491 delete[] temp;
492 }
493
494 template<typename T>
write_dataset(hid_t obj_id,const char * name,const vector<T> & buffer)495 inline void write_dataset(
496 hid_t obj_id, const char* name, const vector<T>& buffer)
497 {
498 hsize_t dims[] {buffer.size()};
499 write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap<T>::type_id,
500 H5S_ALL, false, buffer.data());
501 }
502
503 // Template for xarray, xtensor, etc.
504 template<typename D> inline void
write_dataset(hid_t obj_id,const char * name,const xt::xcontainer<D> & arr)505 write_dataset(hid_t obj_id, const char* name, const xt::xcontainer<D>& arr)
506 {
507 using T = typename D::value_type;
508 auto s = arr.shape();
509 vector<hsize_t> dims {s.cbegin(), s.cend()};
510 write_dataset_lowlevel(obj_id, dims.size(), dims.data(), name,
511 H5TypeMap<T>::type_id, H5S_ALL, false, arr.data());
512 }
513
514 inline void
write_dataset(hid_t obj_id,const char * name,Position r)515 write_dataset(hid_t obj_id, const char* name, Position r)
516 {
517 array<double, 3> buffer {r.x, r.y, r.z};
518 write_dataset(obj_id, name, buffer);
519 }
520
521 inline void
write_dataset(hid_t obj_id,const char * name,std::string buffer)522 write_dataset(hid_t obj_id, const char* name, std::string buffer)
523 {
524 write_string(obj_id, name, buffer.c_str(), false);
525 }
526
527 } // namespace openmc
528 #endif // OPENMC_HDF5_INTERFACE_H
529