diff --git a/CMakeLists.txt b/CMakeLists.txt index e29b5429d3..73c2ebcf8d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -447,9 +447,9 @@ add_library(libopenmc SHARED src/tallies/tally.cpp src/timer.cpp src/thermal.cpp + src/wmp.cpp src/xml_interface.cpp - src/xsdata.cpp - src/tallies/tally.cpp) + src/xsdata.cpp) set_target_properties(libopenmc PROPERTIES OUTPUT_NAME openmc diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index 5b5d2a9991..b0c3acd4d4 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -253,15 +253,18 @@ void read_dataset(hid_t dset, xt::xarray& arr, bool indep=false) std::size_t size = 1; for (const auto x : shape) size *= x; - std::vector buffer(size); + T* buffer = new T[size]; // Read data from attribute - read_dataset(dset, nullptr, H5TypeMap::type_id, buffer.data(), indep); + read_dataset(dset, nullptr, H5TypeMap::type_id, buffer, indep); // Adapt into xarray - arr = xt::adapt(buffer, shape); + arr = xt::adapt(buffer, size, xt::acquire_ownership(), shape); } +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) { diff --git a/include/openmc/wmp.h b/include/openmc/wmp.h new file mode 100644 index 0000000000..c8058d3694 --- /dev/null +++ b/include/openmc/wmp.h @@ -0,0 +1,35 @@ +#ifndef OPENMC_WMP_H +#define OPENMC_WMP_H + +#include "hdf5.h" +#include "xtensor/xtensor.hpp" + +#include + +namespace openmc { + +class WindowedMultipole { +public: + // Types, aliases + using cdouble = std::complex; + + // Constructors, destructors + WindowedMultipole(hid_t group); + + // Data members + std::string name_; //!< Name of nuclide + bool fissionable_; //!< Is the nuclide fissionable? + xt::xtensor data_; //!< Poles and residues + double sqrt_awr_; //!< Square root of atomic weight ratio + double E_min_; //!< Minimum energy in [eV] + double E_max_; //!< Maximum energy in [eV] + double spacing_; //!< Spacing in sqrt(E) space + int fit_order_; //!< Order of the fit + xt::xtensor windows_; //!< Indices of pole at start/end of window + xt::xtensor curvefit_; //!< Fitting function (reaction, coeff index, window index) + xt::xtensor broaden_poly_; //!< Whether to broaden curvefit +}; + +} // namespace openmc + +#endif // OPENMC_WMP_H diff --git a/src/hdf5_interface.cpp b/src/hdf5_interface.cpp index 865b6eb202..980a8688f2 100644 --- a/src/hdf5_interface.cpp +++ b/src/hdf5_interface.cpp @@ -487,6 +487,25 @@ read_dataset(hid_t obj_id, const char* name, hid_t mem_type_id, if (name) H5Dclose(dset); } +template<> +void read_dataset(hid_t dset, xt::xarray>& arr, bool indep) +{ + // Get shape of dataset + std::vector shape = object_shape(dset); + + // Allocate new array to read data into + std::size_t size = 1; + for (const auto x : shape) + size *= x; + std::vector> buffer(size); + + // Read data from attribute + read_complex(dset, nullptr, buffer.data(), indep); + + // Adapt into xarray + arr = xt::adapt(buffer, shape); +} + void read_double(hid_t obj_id, const char* name, double* buffer, bool indep) @@ -757,6 +776,8 @@ using_mpio_device(hid_t obj_id) // Specializations of the H5TypeMap template struct template<> +const hid_t H5TypeMap::type_id = H5T_NATIVE_INT8; +template<> const hid_t H5TypeMap::type_id = H5T_NATIVE_INT; template<> const hid_t H5TypeMap::type_id = H5T_NATIVE_ULONG; diff --git a/src/wmp.cpp b/src/wmp.cpp new file mode 100644 index 0000000000..5ad4218c80 --- /dev/null +++ b/src/wmp.cpp @@ -0,0 +1,47 @@ +#include "openmc/wmp.h" + +#include "openmc/hdf5_interface.h" + +namespace openmc { + +WindowedMultipole::WindowedMultipole(hid_t group) +{ + // Get name of nuclide from group, removing leading '/' + name_ = object_name(group).substr(1); + + // Read scalar values. + read_dataset(group, "spacing", spacing_); + read_dataset(group, "sqrtAWR", sqrt_awr_); + read_dataset(group, "E_min", E_min_); + read_dataset(group, "E_max", E_max_); + + // Read the "data" array. Use its shape to figure out the number of poles + // and residue types in this data. + read_dataset(group, "data", data_); + int n_residues = data_.shape()[1] - 1; + + // Check to see if this data includes fission residues. + fissionable_ = (n_residues == 3); + + // Read the "windows" array and use its shape to figure out the number of + // windows. + read_dataset(group, "windows", windows_); + int n_windows = windows_.shape()[0]; + + // Read the "broaden_poly" arrays. + read_dataset(group, "broaden_poly", broaden_poly_); + if (n_windows != broaden_poly_.shape()[0]) { + fatal_error("broaden_poly array shape is not consistent with the windows " + "array shape in WMP library for " + name_ + "."); + } + + // Read the "curvefit" array. + read_dataset(group, "curvefit", curvefit_); + if (n_windows != broaden_poly_.shape()[0]) { + fatal_error("curvefit array shape is not consistent with the windows " + "array shape in WMP library for " + name_ + "."); + } + fit_order_ = curvefit_.shape()[1] - 1; +} + +} // namespace openmc