From 06751d53734647b3371afeb3352aafb387675489 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 24 Jan 2019 07:26:44 -0600 Subject: [PATCH] Move openmc_load_nuclide, openmc_material_add_nuclide, and read_mgxs to C++ --- include/openmc/capi.h | 2 +- include/openmc/cross_sections.h | 8 + include/openmc/hdf5_interface.h | 7 + include/openmc/material.h | 3 - include/openmc/mgxs.h | 35 +---- include/openmc/mgxs_interface.h | 9 +- include/openmc/nuclide.h | 3 +- include/openmc/scattdata.h | 6 +- include/openmc/xsdata.h | 38 ++--- src/cross_sections.cpp | 220 +++++++++++++++++++++++++++ src/geometry_header.F90 | 49 ------ src/initialize.cpp | 3 +- src/input_xml.F90 | 6 - src/material.cpp | 261 ++++++-------------------------- src/material_header.F90 | 67 -------- src/mgxs.cpp | 54 +++---- src/mgxs_data.F90 | 85 +---------- src/mgxs_interface.F90 | 30 +--- src/mgxs_interface.cpp | 97 +++++++++--- src/nuclide.cpp | 56 ++++++- src/nuclide_header.F90 | 97 ++---------- src/scattdata.cpp | 5 +- src/xsdata.cpp | 141 ++++++++--------- 23 files changed, 552 insertions(+), 730 deletions(-) diff --git a/include/openmc/capi.h b/include/openmc/capi.h index aad72a1386..84e5f5318f 100644 --- a/include/openmc/capi.h +++ b/include/openmc/capi.h @@ -59,7 +59,7 @@ extern "C" { int openmc_init(int argc, char* argv[], const void* intracomm); int openmc_legendre_filter_get_order(int32_t index, int* order); int openmc_legendre_filter_set_order(int32_t index, int order); - int openmc_load_nuclide(const char name[]); + int openmc_load_nuclide(const char* name); int openmc_material_add_nuclide(int32_t index, const char name[], double density); int openmc_material_get_densities(int32_t index, int** nuclides, double** densities, int* n); int openmc_material_get_id(int32_t index, int32_t* id); diff --git a/include/openmc/cross_sections.h b/include/openmc/cross_sections.h index 11af53ea64..d496f56d39 100644 --- a/include/openmc/cross_sections.h +++ b/include/openmc/cross_sections.h @@ -59,6 +59,14 @@ extern std::map library_map; //! libraries void read_cross_sections_xml(); + +//! Load nuclide and thermal scattering data from HDF5 files +//! +//! \param[in] nuc_temps Temperatures for each nuclide in [K] +//! \param[in] thermal_temps Temperatures for each thermal scattering table in [K] +void read_ce_cross_sections(const std::vector>& nuc_temps, + const std::vector>& thermal_temps); + //! Read cross_sections.xml and populate data libraries void read_ce_cross_sections_xml(); diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index b0c3acd4d4..d4b6f48e5f 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -137,6 +137,13 @@ 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, std::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, std::vector& vec) diff --git a/include/openmc/material.h b/include/openmc/material.h index 0aaf5db120..eda102d22a 100644 --- a/include/openmc/material.h +++ b/include/openmc/material.h @@ -93,9 +93,6 @@ private: void calculate_photon_xs(const Particle& p) const; }; -void read_ce_cross_sections(const std::vector>& nuc_temps, - const std::vector>& thermal_temps); - //============================================================================== // Fortran compatibility //============================================================================== diff --git a/include/openmc/mgxs.h b/include/openmc/mgxs.h index e80355bee7..9d680434ce 100644 --- a/include/openmc/mgxs.h +++ b/include/openmc/mgxs.h @@ -57,34 +57,25 @@ class Mgxs { //! @param in_fissionable Is this item fissionable or not. //! @param in_scatter_format Denotes whether Legendre, Tabular, or //! Histogram scattering is used. - //! @param in_num_groups Number of energy groups. - //! @param in_num_delayed_groups Number of delayed groups. //! @param in_is_isotropic Is this an isotropic or angular with respect to //! the incoming particle. //! @param in_polar Polar angle grid. //! @param in_azimuthal Azimuthal angle grid. void init(const std::string& in_name, double in_awr, const std::vector& in_kTs, - bool in_fissionable, int in_scatter_format, int in_num_groups, - int in_num_delayed_groups, bool in_is_isotropic, + bool in_fissionable, int in_scatter_format, bool in_is_isotropic, const std::vector& in_polar, const std::vector& in_azimuthal); //! \brief Initializes the Mgxs object metadata from the HDF5 file //! //! @param xs_id HDF5 group id for the cross section data. - //! @param in_num_groups Number of energy groups. - //! @param in_num_delayed_groups Number of delayed groups. //! @param temperature Temperatures to read. - //! @param tolerance Tolerance of temperature selection method. //! @param temps_to_read Resultant list of temperatures in the library //! to read which correspond to the requested temperatures. //! @param order_dim Resultant dimensionality of the scattering order. - //! @param method Method of choosing nearest temperatures. void - metadata_from_hdf5(hid_t xs_id, int in_num_groups, - int in_num_delayed_groups, const std::vector& temperature, - double tolerance, std::vector& temps_to_read, int& order_dim, - int& method); + metadata_from_hdf5(hid_t xs_id, const std::vector& temperature, + std::vector& temps_to_read, int& order_dim); //! \brief Performs the actual act of combining the microscopic data for a //! single temperature. @@ -118,21 +109,8 @@ class Mgxs { //! \brief Constructor that loads the Mgxs object from the HDF5 file //! //! @param xs_id HDF5 group id for the cross section data. - //! @param energy_groups Number of energy groups. - //! @param delayed_groups Number of delayed groups. //! @param temperature Temperatures to read. - //! @param tolerance Tolerance of temperature selection method. - //! @param max_order Maximum order requested by the user; - //! this is only used for Legendre scattering. - //! @param legendre_to_tabular Flag to denote if any Legendre provided - //! should be converted to a Tabular representation. - //! @param legendre_to_tabular_points If a conversion is requested, this - //! provides the number of points to use in the tabular representation. - //! @param method Method of choosing nearest temperatures. - Mgxs(hid_t xs_id, int energy_groups, - int delayed_groups, const std::vector& temperature, double tolerance, - int max_order, bool legendre_to_tabular, - int legendre_to_tabular_points, int& method); + Mgxs(hid_t xs_id, const std::vector& temperature); //! \brief Constructor that initializes and populates all data to build a //! macroscopic cross section from microscopic cross section. @@ -141,11 +119,8 @@ class Mgxs { //! @param mat_kTs temperatures (in units of eV) that data is needed. //! @param micros Microscopic objects to combine. //! @param atom_densities Atom densities of those microscopic quantities. - //! @param tolerance Tolerance of temperature selection method. - //! @param method Method of choosing nearest temperatures. Mgxs(const std::string& in_name, const std::vector& mat_kTs, - const std::vector& micros, const std::vector& atom_densities, - double tolerance, int& method); + const std::vector& micros, const std::vector& atom_densities); //! \brief Provides a cross section value given certain parameters //! diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index 17be1a18e9..41cc63c1ec 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -20,6 +20,7 @@ namespace data { extern std::vector nuclides_MG; extern std::vector macro_xs; extern "C" int num_energy_groups; +extern "C" int num_delayed_groups; extern std::vector energy_bins; extern std::vector energy_bin_avg; extern std::vector rev_energy_bins; @@ -30,11 +31,11 @@ extern std::vector rev_energy_bins; // Mgxs data loading interface methods //============================================================================== +void read_mgxs(); + extern "C" void -add_mgxs_c(hid_t file_id, const char* name, int energy_groups, - int delayed_groups, int n_temps, const double temps[], double tolerance, - int max_order, bool legendre_to_tabular, int legendre_to_tabular_points, - int& method); +add_mgxs_c(hid_t file_id, const std::string& name, + const std::vector& temperature); extern "C" bool query_fissionable_c(int n_nuclides, const int i_nuclides[]); diff --git a/include/openmc/nuclide.h b/include/openmc/nuclide.h index 1e621e5401..c34dbba18d 100644 --- a/include/openmc/nuclide.h +++ b/include/openmc/nuclide.h @@ -97,7 +97,7 @@ public: }; // Constructors - Nuclide(hid_t group, const double* temperature, int n, int i_nuclide); + Nuclide(hid_t group, const std::vector& temperature, int i_nuclide); //! Initialize logarithmic grid for energy searches void init_grid(); @@ -173,6 +173,7 @@ private: // Non-member functions //============================================================================== +//! Checks for the right version of nuclear data within HDF5 files void check_data_version(hid_t file_id); //============================================================================== diff --git a/include/openmc/scattdata.h b/include/openmc/scattdata.h index 8ea81f252b..1a3a252a51 100644 --- a/include/openmc/scattdata.h +++ b/include/openmc/scattdata.h @@ -140,8 +140,7 @@ class ScattDataLegendre: public ScattData { // Friend convert_legendre_to_tabular so it has access to protected // parameters friend void - convert_legendre_to_tabular(ScattDataLegendre& leg, - ScattDataTabular& tab, int n_mu); + convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab); public: @@ -223,8 +222,7 @@ class ScattDataTabular: public ScattData { // Friend convert_legendre_to_tabular so it has access to protected // parameters friend void - convert_legendre_to_tabular(ScattDataLegendre& leg, - ScattDataTabular& tab, int n_mu); + convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab); public: diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h index fa75b10aa9..a629d786af 100644 --- a/include/openmc/xsdata.h +++ b/include/openmc/xsdata.h @@ -24,50 +24,42 @@ class XsData { private: //! \brief Reads scattering data from the HDF5 file void - scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, - int scatter_format, int final_scatter_format, int order_data, - int max_order, int legendre_to_tabular_points); + scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, + int scatter_format, int final_scatter_format, int order_data); //! \brief Reads fission data from the HDF5 file void - fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, - size_t delayed_groups, bool is_isotropic); + fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic); //! \brief Reads fission data formatted as chi and nu-fission vectors from // the HDF5 file when beta is provided. void - fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups, bool is_isotropic); + fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic); //! \brief Reads fission data formatted as chi and nu-fission vectors from // the HDF5 file when beta is not provided. void - fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups); + fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang); //! \brief Reads fission data formatted as chi and nu-fission vectors from // the HDF5 file when no delayed data is provided. void - fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups); + fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang); //! \brief Reads fission data formatted as a nu-fission matrix from // the HDF5 file when beta is provided. void - fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups, bool is_isotropic); + fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic); //! \brief Reads fission data formatted as a nu-fission matrix from // the HDF5 file when beta is not provided. void - fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups); + fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang); //! \brief Reads fission data formatted as a nu-fission matrix from // the HDF5 file when no delayed data is provided. void - fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups); + fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang); public: @@ -106,8 +98,7 @@ class XsData { //! @param scatter_format The scattering representation of the file. //! @param n_pol Number of polar angles. //! @param n_azi Number of azimuthal angles. - XsData(size_t num_groups, size_t num_delayed_groups, bool fissionable, - int scatter_format, int n_pol, int n_azi); + XsData(bool fissionable, int scatter_format, int n_pol, int n_azi); //! \brief Loads the XsData object from the HDF5 file //! @@ -118,20 +109,13 @@ class XsData { //! this is different from scatter_format if converting a Legendre to //! a tabular representation. //! @param order_data The dimensionality of the scattering data in the file. - //! @param max_order Maximum order requested by the user; - //! this is only used for Legendre scattering. - //! @param legendre_to_tabular Flag to denote if any Legendre provided - //! should be converted to a Tabular representation. - //! @param legendre_to_tabular_points If a conversion is requested, this - //! provides the number of points to use in the tabular representation. //! @param is_isotropic Is this an isotropic or angular with respect to //! the incoming particle. //! @param n_pol Number of polar angles. //! @param n_azi Number of azimuthal angles. void from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, - int final_scatter_format, int order_data, int max_order, - int legendre_to_tabular_points, bool is_isotropic, int n_pol, + int final_scatter_format, int order_data, bool is_isotropic, int n_pol, int n_azi); //! \brief Combines the microscopic data to a macroscopic object. diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index 08b4cbad2e..5a33277a62 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -4,13 +4,21 @@ #include "openmc/container_util.h" #include "openmc/error.h" #include "openmc/file_utils.h" +#include "openmc/hdf5_interface.h" +#include "openmc/material.h" +#include "openmc/message_passing.h" +#include "openmc/nuclide.h" +#include "openmc/photon.h" #include "openmc/settings.h" +#include "openmc/simulation.h" #include "openmc/string_utils.h" +#include "openmc/thermal.h" #include "openmc/xml_interface.h" #include "pugixml.hpp" #include // for getenv +#include namespace openmc { @@ -152,6 +160,218 @@ void read_cross_sections_xml() } } +extern "C" void read_multipole_data(int i_nuclide); +extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr, + const double* temps, int n, int n_nuclide); + +void +read_ce_cross_sections(const std::vector>& nuc_temps, + const std::vector>& thermal_temps) +{ + std::unordered_set already_read; + + // Construct a vector of nuclide names because we haven't loaded nuclide data + // yet, but we need to know the name of the i-th nuclide + std::vector nuclide_names(data::nuclide_map.size()); + std::vector thermal_names(data::thermal_scatt_map.size()); + for (const auto& kv : data::nuclide_map) { + nuclide_names[kv.second] = kv.first; + } + for (const auto& kv : data::thermal_scatt_map) { + thermal_names[kv.second] = kv.first; + } + + // Read cross sections + for (const auto& mat : model::materials) { + for (int i_nuc : mat->nuclide_) { + // Find name of corresponding nuclide. Because we haven't actually loaded + // data, we don't have the name available, so instead we search through + // all key/value pairs in nuclide_map + std::string& name = nuclide_names[i_nuc]; + + if (already_read.find(name) == already_read.end()) { + LibraryKey key {Library::Type::neutron, name}; + int idx = data::library_map[key]; + std::string& filename = data::libraries[idx].path_; + + write_message("Reading " + name + " from " + filename, 6); + + // Open file and make sure version is sufficient + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Read nuclide data from HDF5 + hid_t group = open_group(file_id, name.c_str()); + int i_nuclide = data::nuclides.size(); + data::nuclides.push_back(std::make_unique( + group, nuc_temps[i_nuc], i_nuclide)); + + // Read from Fortran too + nuclide_from_hdf5(group, data::nuclides.back().get(), + &nuc_temps[i_nuc].front(), nuc_temps[i_nuc].size(), i_nuclide + 1); + + close_group(group); + file_close(file_id); + + // Determine if minimum/maximum energy for this nuclide is greater/less + // than the previous + if (data::nuclides[i_nuclide]->grid_.size() >= 1) { + // TODO: off-by-one + int neutron = static_cast(ParticleType::neutron) - 1; + data::energy_min[neutron] = std::max(data::energy_min[neutron], + data::nuclides[i_nuclide]->grid_[0].energy.front()); + data::energy_max[neutron] = std::min(data::energy_max[neutron], + data::nuclides[i_nuclide]->grid_[0].energy.back()); + } + + // Add name and alias to dictionary + already_read.insert(name); + + // Check if elemental data has been read, if needed + int pos = name.find_first_of("0123456789"); + std::string element = name.substr(pos); + if (settings::photon_transport) { + if (already_read.find(element) == already_read.end()) { + // Read photon interaction data from HDF5 photon library + LibraryKey key {Library::Type::photon, element}; + int idx = data::library_map[key]; + std::string& filename = data::libraries[idx].path_; + int i_element = data::element_map[element]; + write_message("Reading " + element + " from " + filename, 6); + + // Open file and make sure version is sufficient + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Read element data from HDF5 + hid_t group = open_group(file_id, element.c_str()); + data::elements.emplace_back(group, data::elements.size()); + + // Determine if minimum/maximum energy for this element is greater/less than + // the previous + const auto& elem {data::elements.back()}; + if (elem.energy_.size() >= 1) { + // TODO: off-by-one + int photon = static_cast(ParticleType::photon) - 1; + int n = elem.energy_.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], + std::exp(elem.energy_(1))); + data::energy_max[photon] = std::min(data::energy_max[photon], + std::exp(elem.energy_(n - 1))); + } + + close_group(group); + file_close(file_id); + + // Add element to set + already_read.insert(element); + } + } + + // Read multipole file into the appropriate entry on the nuclides array + if (settings::temperature_multipole) read_multipole_data(i_nuclide); + } + } + } + + for (auto& mat : model::materials) { + // Skip materials with no S(a,b) tables + if (mat->thermal_tables_.empty()) continue; + + for (const auto& table : mat->thermal_tables_) { + // Get name of S(a,b) table + int i_table = table.index_table; + std::string& name = thermal_names[i_table]; + + if (already_read.find(name) == already_read.end()) { + LibraryKey key {Library::Type::thermal, name}; + int idx = data::library_map[key]; + std::string& filename = data::libraries[idx].path_; + + write_message("Reading " + name + " from " + filename, 6); + + // Open file and make sure version matches + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Read thermal scattering data from HDF5 + hid_t group = open_group(file_id, name.c_str()); + data::thermal_scatt.push_back(std::make_unique( + group, thermal_temps[i_table])); + close_group(group); + file_close(file_id); + + // Add name to dictionary + already_read.insert(name); + } + } // thermal_tables_ + } // materials + + // Finish setting up materials (normalizing densities, etc.) + for (auto& mat : model::materials) { + mat->finalize(); + } + + // Set up logarithmic grid for nuclides + for (auto& nuc : data::nuclides) { + nuc->init_grid(); + } + int neutron = static_cast(ParticleType::neutron) - 1; + simulation::log_spacing = std::log(data::energy_max[neutron] / + data::energy_min[neutron]) / settings::n_log_bins; + + if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { + // Determine if minimum/maximum energy for bremsstrahlung is greater/less + // than the current minimum/maximum + if (data::ttb_e_grid.size() >= 1) { + // TODO: off-by-one + int photon = static_cast(ParticleType::photon) - 1; + int n_e = data::ttb_e_grid.size(); + data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1)); + data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1)); + } + + // Take logarithm of energies since they are log-log interpolated + data::ttb_e_grid = xt::log(data::ttb_e_grid); + } + + // Show which nuclide results in lowest energy for neutron transport + for (const auto& nuc : data::nuclides) { + // If a nuclide is present in a material that's not used in the model, its + // grid has not been allocated + if (nuc->grid_.size() > 0) { + double max_E = nuc->grid_[0].energy.back(); + int neutron = static_cast(ParticleType::neutron) - 1; + if (max_E == data::energy_max[neutron]) { + write_message("Maximum neutron transport energy: " + + std::to_string(data::energy_max[neutron]) + " eV for " + + nuc->name_, 7); + if (mpi::master && data::energy_max[neutron] < 20.0e6) { + warning("Maximum neutron energy is below 20 MeV. This may bias " + " the results."); + } + break; + } + } + } + + // If the user wants multipole, make sure we found a multipole library. + if (settings::temperature_multipole) { + bool mp_found = false; + for (const auto& nuc : data::nuclides) { + if (nuc->multipole_) { + mp_found = true; + break; + } + } + if (mpi::master && !mp_found) { + warning("Windowed multipole functionality is turned on, but no multipole " + "libraries were found. Make sure that windowed multipole data is " + "present in your cross_sections.xml file."); + } + } +} + void read_ce_cross_sections_xml() { // Check if cross_sections.xml exists diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 8dbef4bffb..132ec93abe 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -176,55 +176,6 @@ contains sqrtkT = cell_sqrtkT_c(this % ptr, i) end function cell_sqrtkT -!=============================================================================== -! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table -! appears at in the model. Later, this list is used to determine the actual -! temperatures to read (which may be different if interpolation is used) -!=============================================================================== - - subroutine get_temperatures(nuc_temps) - type(VectorReal), allocatable, intent(out) :: nuc_temps(:) - - integer :: i, j, k - integer :: i_nuclide ! index in nuclides array - integer :: i_material - real(8) :: temperature ! temperature in Kelvin - - allocate(nuc_temps(n_nuclides)) - - do i = 1, size(cells) - ! Skip non-material cells. - if (cells(i) % fill() /= C_NONE) cycle - - do j = 1, cells(i) % material_size() - ! Skip void materials - if (cells(i) % material(j) == MATERIAL_VOID) cycle - - ! Get temperature of cell (rounding to nearest integer) - if (cells(i) % sqrtkT_size() > 1) then - temperature = cells(i) % sqrtkT(j-1)**2 / K_BOLTZMANN - else - temperature = cells(i) % sqrtkT(0)**2 / K_BOLTZMANN - end if - - i_material = cells(i) % material(j) - - associate (mat => materials(i_material)) - NUC_NAMES_LOOP: do k = 1, size(mat % names) - ! Get index in nuc_temps array - i_nuclide = nuclide_dict % get(mat % names(k)) - - ! Add temperature if it hasn't already been added - if (find(nuc_temps(i_nuclide), temperature) == -1) then - call nuc_temps(i_nuclide) % push_back(temperature) - end if - end do NUC_NAMES_LOOP - end associate - end do - end do - - end subroutine get_temperatures - !=============================================================================== ! FREE_MEMORY_GEOMETRY deallocates global arrays defined in this module !=============================================================================== diff --git a/src/initialize.cpp b/src/initialize.cpp index 8b4a3bbe96..9f879509a2 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -17,8 +17,8 @@ #include "openmc/error.h" #include "openmc/geometry_aux.h" #include "openmc/hdf5_interface.h" -#include "openmc/material.h" // TODO: remove this #include "openmc/message_passing.h" +#include "openmc/mgxs_interface.h" #include "openmc/nuclide.h" #include "openmc/output.h" #include "openmc/random_lcg.h" @@ -36,7 +36,6 @@ extern "C" void print_version(); extern "C" void read_command_line(); extern "C" void read_geometry_xml(); extern "C" void read_materials_xml(); -extern "C" void read_mgxs(); extern "C" void read_plots_xml(); extern "C" void read_tallies_xml(); extern "C" void write_summary(); diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 83e8e4ff3d..fd21afdbb3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -15,7 +15,6 @@ module input_xml use material_header use mesh_header use message_passing - use mgxs_data, only: create_macro_xs, read_mgxs use mgxs_interface use nuclide_header use multipole_header @@ -97,11 +96,6 @@ module input_xml contains -!=============================================================================== -! READ_INPUT_XML calls each of the separate subroutines for reading settings, -! geometry, materials, and tallies. -!=============================================================================== - !=============================================================================== ! READ_SETTINGS_XML reads data from a settings.xml file and parses it, checking ! for errors and placing properly-formatted data in the right data structures diff --git a/src/material.cpp b/src/material.cpp index 9e58523184..3b3570b2f1 100644 --- a/src/material.cpp +++ b/src/material.cpp @@ -10,6 +10,7 @@ #include "xtensor/xoperation.hpp" #include "xtensor/xview.hpp" +#include "openmc/capi.h" #include "openmc/cross_sections.h" #include "openmc/container_util.h" #include "openmc/error.h" @@ -830,223 +831,57 @@ read_materials(pugi::xml_node* node) } } -extern "C" void read_multipole_data(int i_nuclide); -extern "C" void nuclide_from_hdf5(hid_t group, const Nuclide* ptr, - const double* temps, int n, int n_nuclide); - -void -read_ce_cross_sections(const std::vector>& nuc_temps, - const std::vector>& thermal_temps) -{ - std::unordered_set already_read; - - // Construct a vector of nuclide names because we haven't loaded nuclide data - // yet, but we need to know the name of the i-th nuclide - std::vector nuclide_names(data::nuclide_map.size()); - std::vector thermal_names(data::thermal_scatt_map.size()); - for (const auto& kv : data::nuclide_map) { - nuclide_names[kv.second] = kv.first; - } - for (const auto& kv : data::thermal_scatt_map) { - thermal_names[kv.second] = kv.first; - } - - // Read cross sections - for (const auto& mat : model::materials) { - for (int i_nuc : mat->nuclide_) { - // Find name of corresponding nuclide. Because we haven't actually loaded - // data, we don't have the name available, so instead we search through - // all key/value pairs in nuclide_map - std::string& name = nuclide_names[i_nuc]; - - if (already_read.find(name) == already_read.end()) { - LibraryKey key {Library::Type::neutron, name}; - int idx = data::library_map[key]; - std::string& filename = data::libraries[idx].path_; - - write_message("Reading " + name + " from " + filename, 6); - - // Open file and make sure version is sufficient - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Read nuclide data from HDF5 - hid_t group = open_group(file_id, name.c_str()); - int i_nuclide = data::nuclides.size(); - data::nuclides.push_back(std::make_unique( - group, &nuc_temps[i_nuc].front(), nuc_temps[i_nuc].size(), i_nuclide - )); - - // Read from Fortran too - nuclide_from_hdf5(group, data::nuclides.back().get(), - &nuc_temps[i_nuc].front(), nuc_temps[i_nuc].size(), i_nuclide + 1); - - close_group(group); - file_close(file_id); - - // Determine if minimum/maximum energy for this nuclide is greater/less - // than the previous - if (data::nuclides[i_nuclide]->grid_.size() >= 1) { - // TODO: off-by-one - int neutron = static_cast(ParticleType::neutron) - 1; - data::energy_min[neutron] = std::max(data::energy_min[neutron], - data::nuclides[i_nuclide]->grid_[0].energy.front()); - data::energy_max[neutron] = std::min(data::energy_max[neutron], - data::nuclides[i_nuclide]->grid_[0].energy.back()); - } - - // Add name and alias to dictionary - already_read.insert(name); - - // Check if elemental data has been read, if needed - int pos = name.find_first_of("0123456789"); - std::string element = name.substr(pos); - if (settings::photon_transport) { - if (already_read.find(element) == already_read.end()) { - // Read photon interaction data from HDF5 photon library - LibraryKey key {Library::Type::photon, element}; - int idx = data::library_map[key]; - std::string& filename = data::libraries[idx].path_; - int i_element = data::element_map[element]; - write_message("Reading " + element + " from " + filename, 6); - - // Open file and make sure version is sufficient - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Read element data from HDF5 - hid_t group = open_group(file_id, element.c_str()); - data::elements.emplace_back(group, data::elements.size()); - - // Determine if minimum/maximum energy for this element is greater/less than - // the previous - const auto& elem {data::elements.back()}; - if (elem.energy_.size() >= 1) { - // TODO: off-by-one - int photon = static_cast(ParticleType::photon) - 1; - int n = elem.energy_.size(); - data::energy_min[photon] = std::max(data::energy_min[photon], - std::exp(elem.energy_(1))); - data::energy_max[photon] = std::min(data::energy_max[photon], - std::exp(elem.energy_(n - 1))); - } - - close_group(group); - file_close(file_id); - - // Add element to set - already_read.insert(element); - } - } - - // Read multipole file into the appropriate entry on the nuclides array - if (settings::temperature_multipole) read_multipole_data(i_nuclide); - } - } - } - - for (auto& mat : model::materials) { - // Skip materials with no S(a,b) tables - if (mat->thermal_tables_.empty()) continue; - - for (const auto& table : mat->thermal_tables_) { - // Get name of S(a,b) table - int i_table = table.index_table; - std::string& name = thermal_names[i_table]; - - if (already_read.find(name) == already_read.end()) { - LibraryKey key {Library::Type::thermal, name}; - int idx = data::library_map[key]; - std::string& filename = data::libraries[idx].path_; - - write_message("Reading " + name + " from " + filename, 6); - - // Open file and make sure version matches - hid_t file_id = file_open(filename, 'r'); - check_data_version(file_id); - - // Read thermal scattering data from HDF5 - hid_t group = open_group(file_id, name.c_str()); - data::thermal_scatt.push_back(std::make_unique( - group, thermal_temps[i_table])); - close_group(group); - file_close(file_id); - - // Add name to dictionary - already_read.insert(name); - } - } // thermal_tables_ - } // materials - - // Finish setting up materials (normalizing densities, etc.) - for (auto& mat : model::materials) { - mat->finalize(); - } - - // Set up logarithmic grid for nuclides - for (auto& nuc : data::nuclides) { - nuc->init_grid(); - } - int neutron = static_cast(ParticleType::neutron) - 1; - simulation::log_spacing = std::log(data::energy_max[neutron] / - data::energy_min[neutron]) / settings::n_log_bins; - - if (settings::photon_transport && settings::electron_treatment == ELECTRON_TTB) { - // Determine if minimum/maximum energy for bremsstrahlung is greater/less - // than the current minimum/maximum - if (data::ttb_e_grid.size() >= 1) { - // TODO: off-by-one - int photon = static_cast(ParticleType::photon) - 1; - int n_e = data::ttb_e_grid.size(); - data::energy_min[photon] = std::max(data::energy_min[photon], data::ttb_e_grid(1)); - data::energy_max[photon] = std::min(data::energy_max[photon], data::ttb_e_grid(n_e - 1)); - } - - // Take logarithm of energies since they are log-log interpolated - data::ttb_e_grid = xt::log(data::ttb_e_grid); - } - - // Show which nuclide results in lowest energy for neutron transport - for (const auto& nuc : data::nuclides) { - // If a nuclide is present in a material that's not used in the model, its - // grid has not been allocated - if (nuc->grid_.size() > 0) { - double max_E = nuc->grid_[0].energy.back(); - int neutron = static_cast(ParticleType::neutron) - 1; - if (max_E == data::energy_max[neutron]) { - write_message("Maximum neutron transport energy: " + - std::to_string(data::energy_max[neutron]) + " eV for " + - nuc->name_, 7); - if (mpi::master && data::energy_max[neutron] < 20.0e6) { - warning("Maximum neutron energy is below 20 MeV. This may bias " - " the results."); - } - break; - } - } - } - - // If the user wants multipole, make sure we found a multipole library. - if (settings::temperature_multipole) { - bool mp_found = false; - for (const auto& nuc : data::nuclides) { - if (nuc->multipole_) { - mp_found = true; - break; - } - } - if (mpi::master && !mp_found) { - warning("Windowed multipole functionality is turned on, but no multipole " - "libraries were found. Make sure that windowed multipole data is " - "present in your cross_sections.xml file."); - } - } -} - //============================================================================== // C API //============================================================================== +extern "C" int +openmc_material_add_nuclide(int32_t index, const char* name, double density) +{ + int err = 0; + if (index >= 1 && index <= model::materials.size()) { + Material* m = model::materials[index - 1]; + + // Check if nuclide is already in material + for (int i = 0; i < m->nuclide_.size(); ++i) { + int i_nuc = m->nuclide_[i]; + if (data::nuclides[i_nuc]->name_ == name) { + double awr = data::nuclides[i_nuc]->awr_; + m->density_ += density - m->atom_density_(i); + m->density_gpcc_ += (density - m->atom_density_(i)) + * awr * MASS_NEUTRON / N_AVOGADRO; + m->atom_density_(i) = density; + return 0; + } + } + + // If nuclide wasn't found, extend nuclide/density arrays + err = openmc_load_nuclide(name); + + if (err == 0) { + // Append new nuclide/density + int i_nuc = data::nuclide_map[name]; + m->nuclide_.push_back(i_nuc); + + auto n = m->nuclide_.size(); + + // Create copy of atom_density_ array with one extra entry + xt::xtensor atom_density = xt::zeros({n}); + xt::view(atom_density, xt::range(0, n-1)) = m->atom_density_; + atom_density(n) = density; + m->atom_density_ = atom_density; + + m->density_ += density; + m->density_gpcc_ += density * data::nuclides[i_nuc]->awr_ + * MASS_NEUTRON / N_AVOGADRO; + } + } else { + set_errmsg("Index in materials array is out of bounds."); + return OPENMC_E_OUT_OF_BOUNDS; + } + return err; +} + extern "C" int openmc_material_get_volume(int32_t index, double* volume) { diff --git a/src/material_header.F90 b/src/material_header.F90 index 962b58b669..42ef35c63d 100644 --- a/src/material_header.F90 +++ b/src/material_header.F90 @@ -18,7 +18,6 @@ module material_header public :: free_memory_material public :: openmc_extend_materials public :: openmc_get_material_index - public :: openmc_material_add_nuclide public :: openmc_material_get_id public :: openmc_material_get_densities public :: openmc_material_get_volume @@ -273,72 +272,6 @@ contains end function openmc_get_material_index - function openmc_material_add_nuclide(index, name, density) result(err) bind(C) - ! Add a nuclide at a specified density in atom/b-cm to a material - integer(C_INT32_T), value, intent(in) :: index - character(kind=C_CHAR) :: name(*) - real(C_DOUBLE), value, intent(in) :: density - integer(C_INT) :: err - - integer :: j, k, n - real(8) :: awr - integer, allocatable :: new_nuclide(:) - real(8), allocatable :: new_density(:) - character(:), allocatable :: name_ - - name_ = to_f_string(name) - - err = E_UNASSIGNED - if (index >= 1 .and. index <= size(materials)) then - associate (m => materials(index)) - ! Check if nuclide is already in material - do j = 1, m % n_nuclides - k = m % nuclide(j) - if (nuclides(k) % name == name_) then - awr = nuclide_awr(k) - m % density = m % density + density - m % atom_density(j) - m % density_gpcc = m % density_gpcc + (density - & - m % atom_density(j)) * awr * MASS_NEUTRON / N_AVOGADRO - m % atom_density(j) = density - err = 0 - end if - end do - - ! If nuclide wasn't found, extend nuclide/density arrays - if (err /= 0) then - ! If nuclide hasn't been loaded, load it now - err = openmc_load_nuclide(name) - - if (err == 0) then - ! Extend arrays - n = m % n_nuclides - allocate(new_nuclide(n + 1)) - if (n > 0) new_nuclide(1:n) = m % nuclide - call move_alloc(FROM=new_nuclide, TO=m % nuclide) - - allocate(new_density(n + 1)) - if (n > 0) new_density(1:n) = m % atom_density - call move_alloc(FROM=new_density, TO=m % atom_density) - - ! Append new nuclide/density - k = nuclide_dict % get(name_) - m % nuclide(n + 1) = k - m % atom_density(n + 1) = density - m % density = m % density + density - m % density_gpcc = m % density_gpcc + & - density * nuclide_awr(k) * MASS_NEUTRON / N_AVOGADRO - m % n_nuclides = n + 1 - end if - end if - end associate - else - err = E_OUT_OF_BOUNDS - call set_errmsg("Index in materials array is out of bounds.") - end if - - end function openmc_material_add_nuclide - - function openmc_material_get_densities(index, nuclides, densities, n) & result(err) bind(C) ! returns an array of nuclide densities in a material diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 85107cbf92..43554f5735 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -16,7 +16,9 @@ #include "openmc/error.h" #include "openmc/math_functions.h" +#include "openmc/mgxs_interface.h" #include "openmc/random_lcg.h" +#include "openmc/settings.h" #include "openmc/string_utils.h" @@ -41,7 +43,7 @@ std::vector macro_xs; void Mgxs::init(const std::string& in_name, double in_awr, const std::vector& in_kTs, bool in_fissionable, int in_scatter_format, - int in_num_groups, int in_num_delayed_groups, bool in_is_isotropic, + bool in_is_isotropic, const std::vector& in_polar, const std::vector& in_azimuthal) { // Set the metadata @@ -51,8 +53,8 @@ Mgxs::init(const std::string& in_name, double in_awr, kTs = xt::adapt(in_kTs); fissionable = in_fissionable; scatter_format = in_scatter_format; - num_groups = in_num_groups; - num_delayed_groups = in_num_delayed_groups; + num_groups = data::num_energy_groups; + num_delayed_groups = data::num_delayed_groups; xs.resize(in_kTs.size()); is_isotropic = in_is_isotropic; n_pol = in_polar.size(); @@ -73,9 +75,8 @@ Mgxs::init(const std::string& in_name, double in_awr, //============================================================================== void -Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, - int in_num_delayed_groups, const std::vector& temperature, - double tolerance, std::vector& temps_to_read, int& order_dim, int& method) +Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector& temperature, + std::vector& temps_to_read, int& order_dim) { // get name char char_name[MAX_WORD_LEN]; @@ -114,20 +115,20 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, // If only one temperature is available, lets just use nearest temperature // interpolation - if ((num_temps == 1) && (method == TEMPERATURE_INTERPOLATION)) { + if ((num_temps == 1) && (settings::temperature_method == TEMPERATURE_INTERPOLATION)) { warning("Cross sections for " + strtrim(name) + " are only available " + "at one temperature. Reverting to the nearest temperature " + "method."); - method = TEMPERATURE_NEAREST; + settings::temperature_method = TEMPERATURE_NEAREST; } - switch(method) { + switch(settings::temperature_method) { case TEMPERATURE_NEAREST: // Determine actual temperatures to read for (const auto& T : temperature) { auto i_closest = xt::argmin(xt::abs(available_temps - T))[0]; double temp_actual = available_temps[i_closest]; - if (std::fabs(temp_actual - T) < tolerance) { + if (std::fabs(temp_actual - T) < settings::temperature_tolerance) { if (std::find(temps_to_read.begin(), temps_to_read.end(), std::round(temp_actual)) == temps_to_read.end()) { temps_to_read.push_back(std::round(temp_actual)); @@ -276,39 +277,34 @@ Mgxs::metadata_from_hdf5(hid_t xs_id, int in_num_groups, // Finally use this data to initialize the MGXS Object init(in_name, in_awr, in_kTs, in_fissionable, in_scatter_format, - in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar, + in_is_isotropic, in_polar, in_azimuthal); } //============================================================================== -Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups, - const std::vector& temperature, double tolerance, int max_order, - bool legendre_to_tabular, int legendre_to_tabular_points, int& method) +Mgxs::Mgxs(hid_t xs_id, const std::vector& temperature) { // Call generic data gathering routine (will populate the metadata) int order_data; std::vector temps_to_read; - metadata_from_hdf5(xs_id, energy_groups, delayed_groups, temperature, - tolerance, temps_to_read, order_data, method); + metadata_from_hdf5(xs_id, temperature, temps_to_read, order_data); // Set number of energy and delayed groups int final_scatter_format = scatter_format; - if (legendre_to_tabular) { + if (settings::legendre_to_tabular) { if (scatter_format == ANGLE_LEGENDRE) final_scatter_format = ANGLE_TABULAR; } // Load the more specific XsData information for (int t = 0; t < temps_to_read.size(); t++) { - xs[t] = XsData(energy_groups, delayed_groups, fissionable, - final_scatter_format, n_pol, n_azi); + xs[t] = XsData(fissionable, final_scatter_format, n_pol, n_azi); // Get the temperature as a string and then open the HDF5 group std::string temp_str = std::to_string(temps_to_read[t]) + "K"; hid_t xsdata_grp = open_group(xs_id, temp_str.c_str()); xs[t].from_hdf5(xsdata_grp, fissionable, scatter_format, - final_scatter_format, order_data, max_order, - legendre_to_tabular_points, is_isotropic, n_pol, n_azi); + final_scatter_format, order_data, is_isotropic, n_pol, n_azi); close_group(xsdata_grp); } // end temperature loop @@ -320,8 +316,7 @@ Mgxs::Mgxs(hid_t xs_id, int energy_groups, int delayed_groups, //============================================================================== Mgxs::Mgxs(const std::string& in_name, const std::vector& mat_kTs, - const std::vector& micros, const std::vector& atom_densities, - double tolerance, int& method) + const std::vector& micros, const std::vector& atom_densities) { // Get the minimum data needed to initialize: // Dont need awr, but lets just initialize it anyways @@ -341,13 +336,12 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector& mat_kTs, std::vector in_azimuthal = micros[0]->azimuthal; init(in_name, in_awr, mat_kTs, in_fissionable, in_scatter_format, - in_num_groups, in_num_delayed_groups, in_is_isotropic, in_polar, - in_azimuthal); + in_is_isotropic, in_polar, in_azimuthal); // Create the xs data for each temperature for (int t = 0; t < mat_kTs.size(); t++) { - xs[t] = XsData(in_num_groups, in_num_delayed_groups, in_fissionable, - in_scatter_format, in_polar.size(), in_azimuthal.size()); + xs[t] = XsData(in_fissionable, in_scatter_format, in_polar.size(), + in_azimuthal.size()); // Find the right temperature index to use double temp_desired = mat_kTs[t]; @@ -357,13 +351,13 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector& mat_kTs, std::vector micro_t(micros.size(), 0); std::vector micro_t_interp(micros.size(), 0.); for (int m = 0; m < micros.size(); m++) { - switch(method) { + switch(settings::temperature_method) { case TEMPERATURE_NEAREST: { micro_t[m] = xt::argmin(xt::abs(micros[m]->kTs - temp_desired))[0]; auto temp_actual = micros[m]->kTs[micro_t[m]]; - if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * tolerance) { + if (std::abs(temp_actual - temp_desired) >= K_BOLTZMANN * settings::temperature_tolerance) { std::stringstream msg; msg << "MGXS Library does not contain cross section for " << name << " at or near " << std::round(temp_desired / K_BOLTZMANN) << "K."; @@ -395,7 +389,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector& mat_kTs, // If we are doing nearest temperature interpolation, then we don't need // to do the 2nd temperature int num_interp_points = 2; - if (method == TEMPERATURE_NEAREST) num_interp_points = 1; + if (settings::temperature_method == TEMPERATURE_NEAREST) num_interp_points = 1; for (int interp_point = 0; interp_point < num_interp_points; interp_point++) { std::vector interp(micros.size()); std::vector temp_indices(micros.size()); diff --git a/src/mgxs_data.F90 b/src/mgxs_data.F90 index 80c8cd7a08..05e5571fad 100644 --- a/src/mgxs_data.F90 +++ b/src/mgxs_data.F90 @@ -6,7 +6,7 @@ module mgxs_data use algorithm, only: find use dict_header, only: DictCharInt use error, only: fatal_error, write_message - use geometry_header, only: get_temperatures, cells + use geometry_header, only: cells use hdf5_interface use material_header, only: Material, materials, n_materials use mgxs_interface @@ -19,89 +19,6 @@ module mgxs_data contains -!=============================================================================== -! READ_XS reads all the cross sections for the problem and stores them in -! nuclides and sab_tables arrays -!=============================================================================== - - subroutine read_mgxs() bind(C) - integer :: i ! index in materials array - integer :: j ! index over nuclides in material - integer :: i_nuclide ! index in nuclides array - character(20) :: name ! name of library to load - type(Material), pointer :: mat - type(SetChar) :: already_read - integer(HID_T) :: file_id - logical :: file_exists - type(VectorReal), allocatable, target :: temps(:) - character(MAX_WORD_LEN) :: word - integer, allocatable :: array(:) - - ! Check if MGXS Library exists - inquire(FILE=path_cross_sections, EXIST=file_exists) - if (.not. file_exists) then - - ! Could not find MGXS Library file - call fatal_error("Cross sections HDF5 file '" & - // trim(path_cross_sections) // "' does not exist!") - end if - - call write_message("Loading cross section data...", 5) - - ! Get temperatures - call get_temperatures(temps) - - ! Open file for reading - file_id = file_open(path_cross_sections, 'r', parallel=.true.) - - ! Read filetype - call read_attribute(word, file_id, "filetype") - if (word /= 'mgxs') then - call fatal_error("Provided MGXS Library is not a MGXS Library file.") - end if - - ! Read revision number for the MGXS Library file and make sure it matches - ! with the current version - call read_attribute(array, file_id, "version") - if (any(array /= VERSION_MGXS_LIBRARY)) then - call fatal_error("MGXS Library file version does not match current & - &version supported by OpenMC.") - end if - - ! ========================================================================== - ! READ ALL MGXS CROSS SECTION TABLES - - ! Loop over all files - MATERIAL_LOOP: do i = 1, n_materials - mat => materials(i) - - NUCLIDE_LOOP: do j = 1, mat % n_nuclides - name = trim(mat % names(j)) // C_NULL_CHAR - i_nuclide = mat % nuclide(j) - - if (.not. already_read % contains(name)) then - call add_mgxs_c(file_id, name, num_energy_groups, num_delayed_groups, & - temps(i_nuclide) % size(), temps(i_nuclide) % data, & - temperature_tolerance, max_order, & - logical(legendre_to_tabular, C_BOOL), & - legendre_to_tabular_points, temperature_method) - - call already_read % add(name) - end if - end do NUCLIDE_LOOP - - call mat % set_fissionable( & - logical(query_fissionable_c(mat % n_nuclides, mat % nuclide))) - - end do MATERIAL_LOOP - - call file_close(file_id) - - ! Avoid some valgrind leak errors - call already_read % clear() - - end subroutine read_mgxs - !=============================================================================== ! CREATE_MACRO_XS generates the macroscopic xs from the microscopic input data !=============================================================================== diff --git a/src/mgxs_interface.F90 b/src/mgxs_interface.F90 index 77f629b80e..b09c6a55ca 100644 --- a/src/mgxs_interface.F90 +++ b/src/mgxs_interface.F90 @@ -7,34 +7,6 @@ module mgxs_interface implicit none interface - - subroutine add_mgxs_c(file_id, name, energy_groups, delayed_groups, & - n_temps, temps, tolerance, max_order, legendre_to_tabular, & - legendre_to_tabular_points, method) bind(C) - use ISO_C_BINDING - import HID_T - implicit none - integer(HID_T), value, intent(in) :: file_id - character(kind=C_CHAR),intent(in) :: name(*) - integer(C_INT), value, intent(in) :: energy_groups - integer(C_INT), value, intent(in) :: delayed_groups - integer(C_INT), value, intent(in) :: n_temps - real(C_DOUBLE), intent(in) :: temps(1:n_temps) - real(C_DOUBLE), value, intent(in) :: tolerance - integer(C_INT), value, intent(in) :: max_order - logical(C_BOOL),value, intent(in) :: legendre_to_tabular - integer(C_INT), value, intent(in) :: legendre_to_tabular_points - integer(C_INT), intent(inout) :: method - end subroutine add_mgxs_c - - function query_fissionable_c(n_nuclides, i_nuclides) result(result) bind(C) - use ISO_C_BINDING - implicit none - integer(C_INT), value, intent(in) :: n_nuclides - integer(C_INT), intent(in) :: i_nuclides(1:n_nuclides) - logical(C_BOOL) :: result - end function query_fissionable_c - subroutine create_macro_xs_c(name, n_nuclides, i_nuclides, n_temps, temps, & atom_densities, tolerance, method) bind(C) use ISO_C_BINDING @@ -130,7 +102,7 @@ module mgxs_interface integer(C_INT), bind(C) :: num_energy_groups ! Number of delayed groups - integer(C_INT) :: num_delayed_groups + integer(C_INT), bind(C) :: num_delayed_groups ! Energy group structure with decreasing energy real(8), allocatable :: energy_bins(:) diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index 710d76a526..88e62ed860 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -1,10 +1,17 @@ #include "openmc/mgxs_interface.h" #include +#include #include "openmc/cross_sections.h" #include "openmc/error.h" +#include "openmc/file_utils.h" +#include "openmc/geometry_aux.h" +#include "openmc/hdf5_interface.h" +#include "openmc/material.h" #include "openmc/math_functions.h" +#include "openmc/nuclide.h" +#include "openmc/settings.h" namespace openmc { @@ -25,30 +32,87 @@ std::vector rev_energy_bins; // Mgxs data loading interface methods //============================================================================== -void -add_mgxs_c(hid_t file_id, const char* name, int energy_groups, - int delayed_groups, int n_temps, const double temps[], double tolerance, - int max_order, bool legendre_to_tabular, int legendre_to_tabular_points, - int& method) +void read_mgxs() { - // Convert temps to a vector for the from_hdf5 function - std::vector temperature(temps, temps + n_temps); + // Check if MGXS Library exists + if (!file_exists(settings::path_cross_sections)) { + // Could not find MGXS Library file + fatal_error("Cross sections HDF5 file '" + settings::path_cross_sections + + "' does not exist."); + } + write_message("Loading cross section data...", 5); + + // Get temperatures + std::vector> nuc_temps(data::nuclide_map.size()); + std::vector> dummy; + get_temperatures(nuc_temps, dummy); + + // Open file for reading + hid_t file_id = file_open(settings::path_cross_sections, 'r'); + + // Read filetype + std::string type; + read_attribute(file_id, "filetype", type); + if (type != "mgxs") { + fatal_error("Provided MGXS Library is not a MGXS Library file."); + } + + // Read revision number for the MGXS Library file and make sure it matches + // with the current version + std::array array; + read_attribute(file_id, "version", array); + if (array != VERSION_MGXS_LIBRARY) { + fatal_error("MGXS Library file version does not match current version " + "supported by OpenMC."); + } + + // ========================================================================== + // READ ALL MGXS CROSS SECTION TABLES + + std::unordered_set already_read; + + // Build vector of nuclide names + std::vector nuclide_names(data::nuclide_map.size()); + for (const auto& kv : data::nuclide_map) { + nuclide_names[kv.second] = kv.first; + } + + // Loop over all files + for (const auto& mat : model::materials) { + for (int i_nuc : mat->nuclide_) { + std::string& name = nuclide_names[i_nuc]; + + if (already_read.find(name) == already_read.end()) { + add_mgxs_c(file_id, name, nuc_temps[i_nuc]); + already_read.insert(name); + } + + if (data::nuclides_MG[i_nuc].fissionable) { + mat->fissionable_ = true; + } + } + } + + file_close(file_id); +} + +void +add_mgxs_c(hid_t file_id, const std::string& name, + const std::vector& temperature) +{ write_message("Loading " + std::string(name) + " data...", 6); // Check to make sure cross section set exists in the library hid_t xs_grp; - if (object_exists(file_id, name)) { - xs_grp = open_group(file_id, name); + if (object_exists(file_id, name.c_str())) { + xs_grp = open_group(file_id, name.c_str()); } else { fatal_error("Data for " + std::string(name) + " does not exist in " + "provided MGXS Library"); } - Mgxs mg(xs_grp, energy_groups, delayed_groups, temperature, tolerance, - max_order, legendre_to_tabular, legendre_to_tabular_points, method); - - data::nuclides_MG.push_back(mg); + data::nuclides_MG.emplace_back(xs_grp, temperature); close_group(xs_grp); } @@ -86,13 +150,10 @@ create_macro_xs_c(const char* mat_name, int n_nuclides, const int i_nuclides[], mgxs_ptr[n] = &data::nuclides_MG[i_nuclides[n] - 1]; } - Mgxs macro(mat_name, temperature, mgxs_ptr, atom_densities_vec, - tolerance, method); - data::macro_xs.emplace_back(macro); + data::macro_xs.emplace_back(mat_name, temperature, mgxs_ptr, atom_densities_vec); } else { // Preserve the ordering of materials by including a blank entry - Mgxs macro; - data::macro_xs.emplace_back(macro); + data::macro_xs.emplace_back(); } } diff --git a/src/nuclide.cpp b/src/nuclide.cpp index 7421708029..730284dc07 100644 --- a/src/nuclide.cpp +++ b/src/nuclide.cpp @@ -1,6 +1,8 @@ #include "openmc/nuclide.h" +#include "openmc/capi.h" #include "openmc/container_util.h" +#include "openmc/cross_sections.h" #include "openmc/endf.h" #include "openmc/error.h" #include "openmc/hdf5_interface.h" @@ -47,7 +49,7 @@ int Nuclide::XS_FISSION {2}; int Nuclide::XS_NU_FISSION {3}; int Nuclide::XS_PHOTON_PROD {4}; -Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) +Nuclide::Nuclide(hid_t group, const std::vector& temperature, int i_nuclide) : i_nuclide_{i_nuclide} { // Get name of nuclide from group, removing leading '/' @@ -82,6 +84,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) // temperature range was given, in which case all temperatures in the range // are loaded irrespective of what temperatures actually appear in the model std::vector temps_to_read; + int n = temperature.size(); double T_min = n > 0 ? settings::temperature_range[0] : 0.0; double T_max = n > 0 ? settings::temperature_range[1] : INFTY; if (T_max > 0.0) { @@ -95,8 +98,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) switch (settings::temperature_method) { case TEMPERATURE_NEAREST: // Find nearest temperatures - for (int i = 0; i < n; ++i) { - double T_desired = temperature[i]; + for (double T_desired : temperature) { // Determine closest temperature double min_delta_T = INFTY; @@ -130,9 +132,7 @@ Nuclide::Nuclide(hid_t group, const double* temperature, int n, int i_nuclide) case TEMPERATURE_INTERPOLATION: // If temperature interpolation or multipole is selected, get a list of // bounding temperatures for each actual temperature present in the model - for (int i = 0; i < n; ++i) { - double T_desired = temperature[i]; - + for (double T_desired : temperature) { bool found_pair = false; for (int j = 0; j < temps_available.size() - 1; ++j) { if (temps_available[j] <= T_desired && T_desired < temps_available[j + 1]) { @@ -888,6 +888,50 @@ void check_data_version(hid_t file_id) } } +//============================================================================== +// C API +//============================================================================== + +extern "C" void extend_nuclides(); + +extern "C" int openmc_load_nuclide(const char* name) +{ + if (data::nuclide_map.find(name) == data::nuclide_map.end()) { + const auto& it = data::library_map.find({Library::Type::neutron, name}); + if (it != data::library_map.end()) { + // Extend nuclides array on Fortran side + extend_nuclides(); + + // Get filename for library containing nuclide + int idx = it->second; + std::string& filename = data::libraries[idx].path_; + + // Open file and make sure version is sufficient + hid_t file_id = file_open(filename, 'r'); + check_data_version(file_id); + + // Read nuclide data from HDF5 + hid_t group = open_group(file_id, name); + std::vector temperature; + int i_nuclide = data::nuclides.size(); + data::nuclides.push_back(std::make_unique( + group, temperature, i_nuclide)); + close_group(group); + file_close(file_id); + + // Add entry to nuclide dictionary + data::nuclide_map[name] = i_nuclide; + + // Initialize nuclide grid + data::nuclides.back()->init_grid(); + } else { + set_errmsg("Nuclide '" + std::string{name} + "' is not present in library."); + return OPENMC_E_DATA; + } + } + return 0; +} + //============================================================================== // Fortran compatibility functions //============================================================================== diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 2ff858994e..8af1dbdc97 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -255,13 +255,10 @@ contains integer :: i integer :: i_closest integer :: n_temperature - integer(HID_T) :: nu_group integer(HID_T) :: energy_group, energy_dset integer(HID_T) :: kT_group integer(HID_T) :: rxs_group integer(HID_T) :: rx_group - integer(HID_T) :: fer_group ! fission_energy_release group - integer(HID_T) :: fer_dset integer(HSIZE_T) :: j integer(HSIZE_T) :: dims(1) character(MAX_WORD_LEN) :: temp_str @@ -632,31 +629,6 @@ contains call nuclide_calculate_elastic_xs_c(this % ptr) end subroutine nuclide_calculate_elastic_xs -!=============================================================================== -! CHECK_DATA_VERSION checks for the right version of nuclear data within HDF5 -! files -!=============================================================================== - - subroutine check_data_version(file_id) - integer(HID_T), intent(in) :: file_id - - integer, allocatable :: version(:) - - if (attribute_exists(file_id, 'version')) then - call read_attribute(version, file_id, 'version') - if (version(1) /= HDF5_VERSION(1)) then - call fatal_error("HDF5 data format uses version " // trim(to_str(& - version(1))) // "." // trim(to_str(version(2))) // " whereas & - &your installation of OpenMC expects version " // trim(to_str(& - HDF5_VERSION(1))) // ".x data.") - end if - else - call fatal_error("HDF5 data does not indicate a version. Your & - &installation of OpenMC expects version " // trim(to_str(& - HDF5_VERSION(1))) // ".x data.") - end if - end subroutine check_data_version - !=============================================================================== ! FREE_MEMORY_NUCLIDE deallocates global arrays defined in this module !=============================================================================== @@ -714,63 +686,6 @@ contains end function openmc_get_nuclide_index - function openmc_load_nuclide(name) result(err) bind(C) - ! Load a nuclide from the cross section library - character(kind=C_CHAR), intent(in) :: name(*) - integer(C_INT) :: err - - integer :: n - integer(HID_T) :: file_id - integer(HID_T) :: group_id - character(:), allocatable :: name_ - character(MAX_FILE_LEN) :: filename - real(8) :: minmax(2) = [ZERO, INFINITY] - type(VectorReal) :: temperature - type(Nuclide), allocatable :: new_nuclides(:) - - ! Copy array of C_CHARs to normal Fortran string - name_ = to_f_string(name) - - err = 0 - if (.not. nuclide_dict % has(name_)) then - if (library_present(LIBRARY_NEUTRON, name_)) then - ! allocate extra space in nuclides array - n = n_nuclides - allocate(new_nuclides(n + 1)) - new_nuclides(1:n) = nuclides(:) - call move_alloc(FROM=new_nuclides, TO=nuclides) - n = n + 1 - - filename = library_path(LIBRARY_NEUTRON, name_) - - ! Open file and make sure version is sufficient - file_id = file_open(filename, 'r') - call check_data_version(file_id) - - ! Read nuclide data from HDF5 - group_id = open_group(file_id, name_) - ! call nuclides(n) % from_hdf5(group_id, temperature, & - ! temperature_method, temperature_tolerance, minmax, & - ! master, n) - call close_group(group_id) - call file_close(file_id) - - ! Add entry to nuclide dictionary - call nuclide_dict % set(name_, n) - n_nuclides = n - - ! Initialize nuclide grid - call nuclides(n) % init_grid() - else - err = E_DATA - call set_errmsg("Nuclide '" // trim(name_) // "' is not present & - &in library.") - end if - end if - - end function openmc_load_nuclide - - function openmc_nuclide_name(index, name) result(err) bind(C) ! Return the name of a nuclide with a given index integer(C_INT), value, intent(in) :: index @@ -795,4 +710,16 @@ contains end if end function openmc_nuclide_name + subroutine extend_nuclides() bind(C) + integer :: n + type(Nuclide), allocatable :: new_nuclides(:) + + ! allocate extra space in nuclides array + n = n_nuclides + allocate(new_nuclides(n + 1)) + new_nuclides(1:n) = nuclides(:) + call move_alloc(FROM=new_nuclides, TO=nuclides) + n = n + 1 + end subroutine + end module nuclide_header diff --git a/src/scattdata.cpp b/src/scattdata.cpp index 21e7f9df7a..a8f63d0f69 100644 --- a/src/scattdata.cpp +++ b/src/scattdata.cpp @@ -10,6 +10,7 @@ #include "openmc/error.h" #include "openmc/math_functions.h" #include "openmc/random_lcg.h" +#include "openmc/settings.h" namespace openmc { @@ -850,10 +851,10 @@ ScattDataTabular::combine(const std::vector& those_scatts, //============================================================================== void -convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab, - int n_mu) +convert_legendre_to_tabular(ScattDataLegendre& leg, ScattDataTabular& tab) { // See if the user wants us to figure out how many points to use + int n_mu = settings::legendre_to_tabular_points; if (n_mu == C_NONE) { // then we will use 2 pts if its P0, or the default if a higher order // TODO use an error minimization algorithm that also picks n_mu diff --git a/src/xsdata.cpp b/src/xsdata.cpp index 8abc1edfd2..347793e5b7 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -13,7 +13,9 @@ #include "openmc/constants.h" #include "openmc/error.h" #include "openmc/math_functions.h" +#include "openmc/mgxs_interface.h" #include "openmc/random_lcg.h" +#include "openmc/settings.h" namespace openmc { @@ -22,10 +24,11 @@ namespace openmc { // XsData class methods //============================================================================== -XsData::XsData(size_t energy_groups, size_t num_delayed_groups, bool fissionable, - int scatter_format, int n_pol, int n_azi) +XsData::XsData(bool fissionable, int scatter_format, int n_pol, int n_azi) { size_t n_ang = n_pol * n_azi; + size_t n_dg = data::num_delayed_groups; + size_t n_g = data::num_energy_groups; // check to make sure scatter format is OK before we allocate if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR && @@ -33,7 +36,7 @@ XsData::XsData(size_t energy_groups, size_t num_delayed_groups, bool fissionable fatal_error("Invalid scatter_format!"); } // allocate all [temperature][angle][in group] quantities - std::vector shape = {n_ang, energy_groups}; + std::vector shape {n_ang, n_g}; total = xt::zeros(shape); absorption = xt::zeros(shape); inverse_velocity = xt::zeros(shape); @@ -45,20 +48,20 @@ XsData::XsData(size_t energy_groups, size_t num_delayed_groups, bool fissionable } // allocate decay_rate; [temperature][angle][delayed group] - shape[1] = num_delayed_groups; + shape[1] = n_dg; decay_rate = xt::zeros(shape); if (fissionable) { - shape = {n_ang, num_delayed_groups, energy_groups}; + shape = {n_ang, n_dg, n_g}; // allocate delayed_nu_fission; [temperature][angle][delay group][in group] delayed_nu_fission = xt::zeros(shape); // chi_prompt; [temperature][angle][in group][out group] - shape = {n_ang, energy_groups, energy_groups}; + shape = {n_ang, n_g, n_g}; chi_prompt = xt::zeros(shape); // chi_delayed; [temperature][angle][delay group][in group][out group] - shape = {n_ang, num_delayed_groups, energy_groups, energy_groups}; + shape = {n_ang, n_dg, n_g, n_g}; chi_delayed = xt::zeros(shape); } @@ -78,18 +81,15 @@ XsData::XsData(size_t energy_groups, size_t num_delayed_groups, bool fissionable void XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, - int final_scatter_format, int order_data, int max_order, - int legendre_to_tabular_points, bool is_isotropic, int n_pol, int n_azi) + int final_scatter_format, int order_data, bool is_isotropic, int n_pol, int n_azi) { // Reconstruct the dimension information so it doesn't need to be passed size_t n_ang = n_pol * n_azi; size_t energy_groups = total.shape()[1]; - size_t delayed_groups = decay_rate.shape()[1]; // Set the fissionable-specific data if (fissionable) { - fission_from_hdf5(xsdata_grp, n_ang, energy_groups, delayed_groups, - is_isotropic); + fission_from_hdf5(xsdata_grp, n_ang, is_isotropic); } // Get the non-fission-specific data read_nd_vector(xsdata_grp, "decay_rate", decay_rate); @@ -97,8 +97,8 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, read_nd_vector(xsdata_grp, "inverse-velocity", inverse_velocity); // Get scattering data - scatter_from_hdf5(xsdata_grp, n_ang, energy_groups, scatter_format, - final_scatter_format, order_data, max_order, legendre_to_tabular_points); + scatter_from_hdf5(xsdata_grp, n_ang, scatter_format, + final_scatter_format, order_data); // Check absorption to ensure it is not 0 since it is often the // denominator in tally methods @@ -123,12 +123,15 @@ XsData::from_hdf5(hid_t xsdata_grp, bool fissionable, int scatter_format, void XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups, bool is_isotropic) + bool is_isotropic) { // Data is provided as nu-fission and chi with a beta for delayed info + size_t n_g = data::num_energy_groups; + size_t n_dg = data::num_delayed_groups; + // Get chi - xt::xtensor temp_chi({n_ang, energy_groups}, 0.); + xt::xtensor temp_chi({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "chi", temp_chi, true); // Normalize chi by summing over the outgoing groups for each incoming angle @@ -141,7 +144,7 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, xt::all()); // Get nu-fission - xt::xtensor temp_nufiss({n_ang, energy_groups}, 0.); + xt::xtensor temp_nufiss({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "nu-fission", temp_nufiss, true); // Get beta (strategy will depend upon the number of dimensions in beta) @@ -151,7 +154,7 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, int ndim_target = 1; if (!is_isotropic) ndim_target += 2; if (beta_ndims == ndim_target) { - xt::xtensor temp_beta({n_ang, delayed_groups}, 0.); + xt::xtensor temp_beta({n_ang, n_dg}, 0.); read_nd_vector(xsdata_grp, "beta", temp_beta, true); // Set prompt_nu_fission = (1. - beta_total)*nu_fission @@ -162,8 +165,7 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, xt::view(temp_beta, xt::all(), xt::all(), xt::newaxis()) * xt::view(temp_nufiss, xt::all(), xt::newaxis(), xt::all()); } else if (beta_ndims == ndim_target + 1) { - xt::xtensor temp_beta({n_ang, delayed_groups, energy_groups}, - 0.); + xt::xtensor temp_beta({n_ang, n_dg, n_g}, 0.); read_nd_vector(xsdata_grp, "beta", temp_beta, true); // Set prompt_nu_fission = (1. - beta_total)*nu_fission @@ -176,20 +178,22 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, } void -XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups) +XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang) { // Data is provided separately as prompt + delayed nu-fission and chi + size_t n_g = data::num_energy_groups; + size_t n_dg = data::num_delayed_groups; + // Get chi-prompt - xt::xtensor temp_chi_p({n_ang, energy_groups}, 0.); + xt::xtensor temp_chi_p({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true); // Normalize chi by summing over the outgoing groups for each incoming angle temp_chi_p /= xt::view(xt::sum(temp_chi_p, {1}), xt::all(), xt::newaxis()); // Get chi-delayed - xt::xtensor temp_chi_d({n_ang, delayed_groups, energy_groups}, 0.); + xt::xtensor temp_chi_d({n_ang, n_dg, n_g}, 0.); read_nd_vector(xsdata_grp, "chi-delayed", temp_chi_d, true); // Normalize chi by summing over the outgoing groups for each incoming angle @@ -209,14 +213,15 @@ XsData::fission_vector_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, } void -XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups) +XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang) { // No beta is provided and there is no prompt/delay distinction. // Therefore, the code only considers the data as prompt. + size_t n_g = data::num_energy_groups; + // Get chi - xt::xtensor temp_chi({n_ang, energy_groups}, 0.); + xt::xtensor temp_chi({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "chi", temp_chi, true); // Normalize chi by summing over the outgoing groups for each incoming angle @@ -232,13 +237,15 @@ XsData::fission_vector_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, //============================================================================== void -XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups, bool is_isotropic) +XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic) { // Data is provided as nu-fission and chi with a beta for delayed info + size_t n_g = data::num_energy_groups; + size_t n_dg = data::num_delayed_groups; + // Get nu-fission matrix - xt::xtensor temp_matrix({n_ang, energy_groups, energy_groups}, 0.); + xt::xtensor temp_matrix({n_ang, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); // Get beta (strategy will depend upon the number of dimensions in beta) @@ -248,7 +255,7 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, int ndim_target = 1; if (!is_isotropic) ndim_target += 2; if (beta_ndims == ndim_target) { - xt::xtensor temp_beta({n_ang, delayed_groups}, 0.); + xt::xtensor temp_beta({n_ang, n_dg}, 0.); read_nd_vector(xsdata_grp, "beta", temp_beta, true); xt::xtensor temp_beta_sum({n_ang}, 0.); @@ -274,10 +281,10 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, xt::view(temp_matrix, xt::all(), xt::newaxis(), xt::all(), xt::all()); } else if (beta_ndims == ndim_target + 1) { - xt::xtensor temp_beta({n_ang, delayed_groups, energy_groups}, 0.); + xt::xtensor temp_beta({n_ang, n_dg, n_g}, 0.); read_nd_vector(xsdata_grp, "beta", temp_beta, true); - xt::xtensor temp_beta_sum({n_ang, energy_groups}, 0.); + xt::xtensor temp_beta_sum({n_ang, n_g}, 0.); temp_beta_sum = xt::sum(temp_beta, {1}); // prompt_nu_fission is the sum of this matrix over outgoing groups and @@ -308,13 +315,15 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, } void -XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups, size_t delayed_groups) +XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang) { // Data is provided separately as prompt + delayed nu-fission and chi + size_t n_g = data::num_energy_groups; + size_t n_dg = data::num_delayed_groups; + // Get the prompt nu-fission matrix - xt::xtensor temp_matrix_p({n_ang, energy_groups, energy_groups}, 0.); + xt::xtensor temp_matrix_p({n_ang, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true); // prompt_nu_fission is the sum over outgoing groups @@ -326,8 +335,7 @@ XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, xt::view(prompt_nu_fission, xt::all(), xt::all(), xt::newaxis()); // Get the delayed nu-fission matrix - xt::xtensor temp_matrix_d({n_ang, delayed_groups, energy_groups, - energy_groups}, 0.); + xt::xtensor temp_matrix_d({n_ang, n_dg, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "delayed-nu-fission", temp_matrix_d, true); // delayed_nu_fission is the sum over outgoing groups @@ -340,14 +348,15 @@ XsData::fission_matrix_no_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, } void -XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, - size_t energy_groups) +XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang) { // No beta is provided and there is no prompt/delay distinction. // Therefore, the code only considers the data as prompt. + size_t n_g = data::num_energy_groups; + // Get nu-fission matrix - xt::xtensor temp_matrix({n_ang, energy_groups, energy_groups}, 0.); + xt::xtensor temp_matrix({n_ang, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); // prompt_nu_fission is the sum over outgoing groups @@ -362,8 +371,7 @@ XsData::fission_matrix_no_delayed_from_hdf5(hid_t xsdata_grp, size_t n_ang, //============================================================================== void -XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, - size_t delayed_groups, bool is_isotropic) +XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic) { // Get the fission and kappa_fission data xs; these are optional read_nd_vector(xsdata_grp, "fission", fission); @@ -373,33 +381,29 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, // as a nu-fission matrix or a set of chi and nu-fission vectors if (object_exists(xsdata_grp, "chi") || object_exists(xsdata_grp, "chi-prompt")) { - if (delayed_groups == 0) { - fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang, energy_groups); + if (data::num_delayed_groups == 0) { + fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang); } else { if (object_exists(xsdata_grp, "beta")) { - fission_vector_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, - delayed_groups, is_isotropic); + fission_vector_beta_from_hdf5(xsdata_grp, n_ang, is_isotropic); } else { - fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, - delayed_groups); + fission_vector_no_beta_from_hdf5(xsdata_grp, n_ang); } } } else { - if (delayed_groups == 0) { - fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang, energy_groups); + if (data::num_delayed_groups == 0) { + fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang); } else { if (object_exists(xsdata_grp, "beta")) { - fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, - delayed_groups, is_isotropic); + fission_matrix_beta_from_hdf5(xsdata_grp, n_ang, is_isotropic); } else { - fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang, energy_groups, - delayed_groups); + fission_matrix_no_beta_from_hdf5(xsdata_grp, n_ang); } } } // Combine prompt_nu_fission and delayed_nu_fission into nu_fission - if (delayed_groups == 0) { + if (data::num_delayed_groups == 0) { nu_fission = prompt_nu_fission; } else { nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1}); @@ -409,9 +413,8 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, //============================================================================== void -XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, - int scatter_format, int final_scatter_format, int order_data, - int max_order, int legendre_to_tabular_points) +XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, + int scatter_format, int final_scatter_format, int order_data) { if (!object_exists(xsdata_grp, "scatter_data")) { fatal_error("Must provide scatter_data group!"); @@ -419,9 +422,10 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, hid_t scatt_grp = open_group(xsdata_grp, "scatter_data"); // Get the outgoing group boundary indices - xt::xtensor gmin({n_ang, energy_groups}, 0.); + size_t n_g = data::num_energy_groups; + xt::xtensor gmin({n_ang, n_g}, 0.); read_nd_vector(scatt_grp, "g_min", gmin, true); - xt::xtensor gmax({n_ang, energy_groups}, 0.); + xt::xtensor gmax({n_ang, n_g}, 0.); read_nd_vector(scatt_grp, "g_max", gmax, true); // Make gmin and gmax start from 0 vice 1 as they do in the library @@ -432,7 +436,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, // data. size_t length = order_data * xt::sum(gmax - gmin + 1)(); - double_4dvec input_scatt(n_ang, double_3dvec(energy_groups)); + double_4dvec input_scatt(n_ang, double_3dvec(n_g)); xt::xtensor temp_arr({length}, 0.); read_nd_vector(scatt_grp, "scatter_matrix", temp_arr, true); @@ -440,7 +444,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, // strip off the superfluous orders if needed int order_dim; if (scatter_format == ANGLE_LEGENDRE) { - order_dim = std::min(order_data - 1, max_order) + 1; + order_dim = std::min(order_data - 1, settings::max_order) + 1; } else { order_dim = order_data; } @@ -449,7 +453,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, // scatt data size_t temp_idx = 0; for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gin = 0; gin < n_g; gin++) { input_scatt[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); for (size_t i_gout = 0; i_gout < input_scatt[a][gin].size(); i_gout++) { input_scatt[a][gin][i_gout].resize(order_dim); @@ -463,7 +467,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, } // Get multiplication matrix - double_3dvec temp_mult(n_ang, double_2dvec(energy_groups)); + double_3dvec temp_mult(n_ang, double_2dvec(n_g)); if (object_exists(scatt_grp, "multiplicity_matrix")) { temp_arr.resize({length / order_data}); read_nd_vector(scatt_grp, "multiplicity_matrix", temp_arr); @@ -471,7 +475,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, // convert the flat temp_arr to a jagged array for passing to scatt data size_t temp_idx = 0; for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gin = 0; gin < n_g; gin++) { temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { temp_mult[a][gin][i_gout] = temp_arr[temp_idx++]; @@ -481,7 +485,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, } else { // Use a default: multiplicities are 1.0. for (size_t a = 0; a < n_ang; a++) { - for (size_t gin = 0; gin < energy_groups; gin++) { + for (size_t gin = 0; gin < n_g; gin++) { temp_mult[a][gin].resize(gmax(a, gin) - gmin(a, gin) + 1); for (size_t i_gout = 0; i_gout < temp_mult[a][gin].size(); i_gout++) { temp_mult[a][gin][i_gout] = 1.; @@ -504,8 +508,7 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, size_t energy_groups, // Now create a tabular version of legendre_scatt convert_legendre_to_tabular(legendre_scatt, - *static_cast(scatter[a].get()), - legendre_to_tabular_points); + *static_cast(scatter[a].get())); scatter_format = final_scatter_format; }