#ifndef OPENMC_HDF5_INTERFACE_H #define OPENMC_HDF5_INTERFACE_H #include // for min #include #include #include // for strlen #include #include #include #include "hdf5.h" #include "hdf5_hl.h" #include "xtensor/xadapt.hpp" #include "xtensor/xarray.hpp" #include "openmc/array.h" #include "openmc/error.h" #include "openmc/position.h" #include "openmc/vector.h" namespace openmc { //============================================================================== // Low-level internal functions //============================================================================== void read_attr(hid_t obj_id, const char* name, hid_t mem_type_id, void* buffer); void write_attr(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, hid_t mem_type_id, const void* buffer); void read_dataset_lowlevel(hid_t obj_id, const char* name, hid_t mem_type_id, hid_t mem_space_id, bool indep, void* buffer); void write_dataset_lowlevel(hid_t group_id, int ndim, const hsize_t* dims, const char* name, hid_t mem_type_id, hid_t mem_space_id, bool indep, const void* buffer); bool using_mpio_device(hid_t obj_id); //============================================================================== // Normal functions that are used to read/write files //============================================================================== hid_t create_group(hid_t parent_id, const std::string& name); inline hid_t create_group(hid_t parent_id, const std::stringstream& name) { return create_group(parent_id, name.str()); } hid_t file_open(const std::string& filename, char mode, bool parallel = false); hid_t open_group(hid_t group_id, const std::string& name); void write_string( hid_t group_id, const char* name, const std::string& buffer, bool indep); vector attribute_shape(hid_t obj_id, const char* name); vector dataset_names(hid_t group_id); void ensure_exists(hid_t obj_id, const char* name, bool attribute = false); vector group_names(hid_t group_id); vector object_shape(hid_t obj_id); std::string object_name(hid_t obj_id); hid_t open_object(hid_t group_id, const std::string& name); void close_object(hid_t obj_id); //============================================================================== // Fortran compatibility functions //============================================================================== extern "C" { bool attribute_exists(hid_t obj_id, const char* name); size_t attribute_typesize(hid_t obj_id, const char* name); hid_t create_group(hid_t parent_id, const char* name); void close_dataset(hid_t dataset_id); void close_group(hid_t group_id); int dataset_ndims(hid_t dset); size_t dataset_typesize(hid_t obj_id, const char* name); hid_t file_open(const char* filename, char mode, bool parallel); void file_close(hid_t file_id); void get_name(hid_t obj_id, std::string& name); int get_num_datasets(hid_t group_id); int get_num_groups(hid_t group_id); void get_datasets(hid_t group_id, char* name[]); void get_groups(hid_t group_id, char* name[]); void get_shape(hid_t obj_id, hsize_t* dims); void get_shape_attr(hid_t obj_id, const char* name, hsize_t* dims); bool object_exists(hid_t object_id, const char* name); hid_t open_dataset(hid_t group_id, const char* name); hid_t open_group(hid_t group_id, const char* name); void read_attr_double(hid_t obj_id, const char* name, double* buffer); void read_attr_int(hid_t obj_id, const char* name, int* buffer); void read_attr_string( hid_t obj_id, const char* name, size_t slen, char* buffer); void read_complex( hid_t obj_id, const char* name, std::complex* buffer, bool indep); void read_double(hid_t obj_id, const char* name, double* buffer, bool indep); void read_int(hid_t obj_id, const char* name, int* buffer, bool indep); void read_llong(hid_t obj_id, const char* name, long long* buffer, bool indep); void read_string( hid_t obj_id, const char* name, size_t slen, char* buffer, bool indep); void read_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, hsize_t n_results, double* results); void write_attr_double(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, const double* buffer); void write_attr_int(hid_t obj_id, int ndim, const hsize_t* dims, const char* name, const int* buffer); void write_attr_string(hid_t obj_id, const char* name, const char* buffer); void write_double(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const double* buffer, bool indep); void write_int(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const int* buffer, bool indep); void write_llong(hid_t group_id, int ndim, const hsize_t* dims, const char* name, const long long* buffer, bool indep); void write_string(hid_t group_id, int ndim, const hsize_t* dims, size_t slen, const char* name, const char* buffer, bool indep); void write_tally_results(hid_t group_id, hsize_t n_filter, hsize_t n_score, hsize_t n_results, const double* results); } // extern "C" //============================================================================== // Template struct used to map types to HDF5 datatype IDs, which are stored // using the type hid_t. By having a single static data member, the template can // be specialized for each type we know of. The specializations appear in the // .cpp file since they are definitions. //============================================================================== template struct H5TypeMap { static const hid_t type_id; }; //============================================================================== // Templates/overloads for read_attribute //============================================================================== // Scalar version template void read_attribute(hid_t obj_id, const char* name, T& buffer) { read_attr(obj_id, name, H5TypeMap::type_id, &buffer); } // array version template inline void read_attribute(hid_t obj_id, const char* name, array& buffer) { read_attr(obj_id, name, H5TypeMap::type_id, buffer.data()); } // vector version template void read_attribute(hid_t obj_id, const char* name, vector& vec) { // Get shape of attribute array auto shape = attribute_shape(obj_id, name); // Allocate new array to read data into std::size_t size = 1; for (const auto x : shape) size *= x; vec.resize(size); // Read data from attribute read_attr(obj_id, name, H5TypeMap::type_id, vec.data()); } // Generic array version template void read_attribute(hid_t obj_id, const char* name, xt::xarray& arr) { // Get shape of attribute array auto shape = attribute_shape(obj_id, name); // Allocate new array to read data into std::size_t size = 1; for (const auto x : shape) size *= x; vector buffer(size); // Read data from attribute read_attr(obj_id, name, H5TypeMap::type_id, buffer.data()); // Adapt array into xarray arr = xt::adapt(buffer, shape); } // overload for std::string inline void read_attribute(hid_t obj_id, const char* name, std::string& str) { // Create buffer to read data into auto n = attribute_typesize(obj_id, name); char* buffer = new char[n]; // Read attribute and set string read_attr_string(obj_id, name, n, buffer); str = std::string {buffer, n}; delete[] buffer; } // overload for vector inline void read_attribute( hid_t obj_id, const char* name, vector& vec) { auto dims = attribute_shape(obj_id, name); auto m = dims[0]; // Allocate a C char array to get strings auto n = attribute_typesize(obj_id, name); char* buffer = new char[m * n]; // Read char data in attribute read_attr_string(obj_id, name, n, buffer); for (decltype(m) i = 0; i < m; ++i) { // Determine proper length of string -- strlen doesn't work because // buffer[i] might not have any null characters std::size_t k = 0; for (; k < n; ++k) if (buffer[i * n + k] == '\0') break; // Create string based on (char*, size_t) constructor vec.emplace_back(&buffer[i * n], k); } delete[] buffer; } //============================================================================== // Templates/overloads for read_dataset and related methods //============================================================================== // Template for scalars. We need to be careful that the compiler does not use // this version of read_dataset for vectors, arrays, or other non-scalar types. // enable_if_t allows us to conditionally remove the function from overload // resolution when the type T doesn't meet a certain criterion. template inline std::enable_if_t>::value> read_dataset( hid_t obj_id, const char* name, T& buffer, bool indep = false) { read_dataset_lowlevel( obj_id, name, H5TypeMap::type_id, H5S_ALL, indep, &buffer); } // overload for std::string inline void read_dataset( hid_t obj_id, const char* name, std::string& str, bool indep = false) { // Create buffer to read data into auto n = dataset_typesize(obj_id, name); std::vector buffer(n, '\0'); // Read attribute and set string read_string(obj_id, name, n, buffer.data(), indep); str = std::string {buffer.begin(), buffer.end()}; } // array version template inline void read_dataset( hid_t dset, const char* name, array& buffer, bool indep = false) { read_dataset_lowlevel( dset, name, H5TypeMap::type_id, H5S_ALL, indep, buffer.data()); } // vector version template void read_dataset(hid_t dset, vector& vec, bool indep = false) { // Get shape of dataset vector shape = object_shape(dset); // Resize vector to appropriate size vec.resize(shape[0]); // Read data into vector read_dataset_lowlevel( dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, vec.data()); } template void read_dataset( hid_t obj_id, const char* name, vector& vec, bool indep = false) { hid_t dset = open_dataset(obj_id, name); read_dataset(dset, vec, indep); close_dataset(dset); } template void read_dataset(hid_t dset, xt::xarray& arr, bool indep = false) { // Get shape of dataset vector shape = object_shape(dset); // Allocate space in the array to read data into std::size_t size = 1; for (const auto x : shape) size *= x; arr.resize(shape); // Read data from attribute read_dataset_lowlevel( dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, arr.data()); } template<> void read_dataset( hid_t dset, xt::xarray>& arr, bool indep); template void read_dataset( hid_t obj_id, const char* name, xt::xarray& arr, bool indep = false) { // Open dataset and read array hid_t dset = open_dataset(obj_id, name); read_dataset(dset, arr, indep); close_dataset(dset); } template void read_dataset( hid_t obj_id, const char* name, xt::xtensor& arr, bool indep = false) { // Open dataset and read array hid_t dset = open_dataset(obj_id, name); // Get shape of dataset vector hsize_t_shape = object_shape(dset); close_dataset(dset); // cast from hsize_t to size_t vector shape(hsize_t_shape.size()); for (int i = 0; i < shape.size(); i++) { shape[i] = static_cast(hsize_t_shape[i]); } // Allocate new xarray to read data into xt::xarray xarr(shape); // Read data from the dataset read_dataset(obj_id, name, xarr); // Copy into xtensor arr = xarr; } // overload for Position inline void read_dataset( hid_t obj_id, const char* name, Position& r, bool indep = false) { array x; read_dataset(obj_id, name, x, indep); r.x = x[0]; r.y = x[1]; r.z = x[2]; } template inline void read_dataset_as_shape( hid_t obj_id, const char* name, xt::xtensor& arr, bool indep = false) { hid_t dset = open_dataset(obj_id, name); // Allocate new array to read data into std::size_t size = 1; for (const auto x : arr.shape()) size *= x; vector buffer(size); // Read data from attribute read_dataset_lowlevel( dset, nullptr, H5TypeMap::type_id, H5S_ALL, indep, buffer.data()); // Adapt into xarray arr = xt::adapt(buffer, arr.shape()); close_dataset(dset); } template inline void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, bool must_have = false) { if (object_exists(obj_id, name)) { read_dataset_as_shape(obj_id, name, result, true); } else if (must_have) { fatal_error(std::string("Must provide " + std::string(name) + "!")); } } //============================================================================== // Templates/overloads for write_attribute //============================================================================== template inline void write_attribute(hid_t obj_id, const char* name, T buffer) { write_attr(obj_id, 0, nullptr, name, H5TypeMap::type_id, &buffer); } inline void write_attribute(hid_t obj_id, const char* name, const char* buffer) { write_attr_string(obj_id, name, buffer); } inline void write_attribute( hid_t obj_id, const char* name, const std::string& buffer) { write_attr_string(obj_id, name, buffer.c_str()); } template inline void write_attribute( hid_t obj_id, const char* name, const array& buffer) { hsize_t dims[] {N}; write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); } template inline void write_attribute( hid_t obj_id, const char* name, const vector& buffer) { hsize_t dims[] {buffer.size()}; write_attr(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data()); } inline void write_attribute(hid_t obj_id, const char* name, Position r) { array buffer {r.x, r.y, r.z}; write_attribute(obj_id, name, buffer); } //============================================================================== // Templates/overloads for write_dataset //============================================================================== // Template for scalars (ensured by SFINAE) template inline std::enable_if_t>::value> write_dataset( hid_t obj_id, const char* name, T buffer) { write_dataset_lowlevel( obj_id, 0, nullptr, name, H5TypeMap::type_id, H5S_ALL, false, &buffer); } inline void write_dataset(hid_t obj_id, const char* name, const char* buffer) { write_string(obj_id, name, buffer, false); } template inline void write_dataset( hid_t obj_id, const char* name, const array& buffer) { hsize_t dims[] {N}; write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap::type_id, H5S_ALL, false, buffer.data()); } inline void write_dataset( hid_t obj_id, const char* name, const vector& buffer) { auto n {buffer.size()}; hsize_t dims[] {n}; // Determine length of longest string, including \0 size_t m = 1; for (const auto& s : buffer) { m = std::max(m, s.size() + 1); } // Copy data into contiguous buffer char* temp = new char[n * m]; std::fill(temp, temp + n * m, '\0'); for (decltype(n) i = 0; i < n; ++i) { std::copy(buffer[i].begin(), buffer[i].end(), temp + i * m); } // Write 2D data write_string(obj_id, 1, dims, m, name, temp, false); // Free temp array delete[] temp; } template inline void write_dataset( hid_t obj_id, const char* name, const vector& buffer) { hsize_t dims[] {buffer.size()}; write_dataset_lowlevel(obj_id, 1, dims, name, H5TypeMap::type_id, H5S_ALL, false, buffer.data()); } // Template for xarray, xtensor, etc. template inline void write_dataset( hid_t obj_id, const char* name, const xt::xcontainer& arr) { using T = typename D::value_type; auto s = arr.shape(); vector dims {s.cbegin(), s.cend()}; write_dataset_lowlevel(obj_id, dims.size(), dims.data(), name, H5TypeMap::type_id, H5S_ALL, false, arr.data()); } inline void write_dataset(hid_t obj_id, const char* name, Position r) { array buffer {r.x, r.y, r.z}; write_dataset(obj_id, name, buffer); } inline void write_dataset(hid_t obj_id, const char* name, std::string buffer) { write_string(obj_id, name, buffer.c_str(), false); } } // namespace openmc #endif // OPENMC_HDF5_INTERFACE_H