From a2bcb07e0333073a16e2035b05e967ffea571d2e Mon Sep 17 00:00:00 2001 From: Gavin Ridley Date: Thu, 7 Nov 2019 13:27:16 -0500 Subject: [PATCH] Mgxs may now be used in external linked programs --- include/openmc/hdf5_interface.h | 4 +- include/openmc/mgxs_interface.h | 38 +++++++- include/openmc/xsdata.h | 8 +- src/cross_sections.cpp | 3 +- src/initialize.cpp | 5 +- src/mgxs.cpp | 5 +- src/mgxs_interface.cpp | 159 +++++++++++++++++++++----------- src/xsdata.cpp | 30 ++---- 8 files changed, 163 insertions(+), 89 deletions(-) diff --git a/include/openmc/hdf5_interface.h b/include/openmc/hdf5_interface.h index d0474c8344..ee3419fa22 100644 --- a/include/openmc/hdf5_interface.h +++ b/include/openmc/hdf5_interface.h @@ -345,7 +345,7 @@ read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false) } template -void read_dataset_as_shape(hid_t obj_id, const char* name, +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); @@ -367,7 +367,7 @@ void read_dataset_as_shape(hid_t obj_id, const char* name, template -void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor& result, +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)) { diff --git a/include/openmc/mgxs_interface.h b/include/openmc/mgxs_interface.h index c72587be2b..0c1091c4b3 100644 --- a/include/openmc/mgxs_interface.h +++ b/include/openmc/mgxs_interface.h @@ -20,6 +20,14 @@ struct MgxsInterface int num_energy_groups; int num_delayed_groups; + // List of available names in the HDF5 file + std::vector xs_names; + std::vector xs_to_read; + std::vector> xs_temps_to_read; + + // Name of the HDF5 file which contains mgxs + std::string cross_sections_path; + std::vector nuclides_MG; std::vector macro_xs; @@ -27,12 +35,20 @@ struct MgxsInterface std::vector energy_bin_avg; std::vector rev_energy_bins; + // temperatues of each available nuclide + std::vector> nuc_temps; + MgxsInterface() = default; - // Construct from path to cross sections file - MgxsInterface(const std::string& path_cross_sections); + // Construct from path to cross sections file, as well as a list + // of XS to read and the corresponding temperatures for each XS + MgxsInterface(const std::string& path_cross_sections, + const std::vector xs_to_read, + const std::vector> xs_temps); + void setNuclidesToRead(std::vector arg_xs_to_read); + void setNuclideTemperaturesToRead(std::vector> xs_temps); - void init(const std::string& path_cross_sections); + void init(); void add_mgxs(hid_t file_id, const std::string& name, const std::vector& temperature); @@ -41,13 +57,27 @@ struct MgxsInterface std::vector> get_mat_kTs(); - void read_mg_cross_sections_header(); + // Reads just the header of the cross sections file, to find + // min & max energies as well as the available XS + void readHeader(const std::string& path_cross_sections); }; namespace data { extern MgxsInterface mgInterface; } +// Puts available XS in MGXS file to globals so that when +// materials are read, the MGXS specified in a material can +// be ensured to be present in the available data. +void putMgxsHeaderDataToGlobals(); + +// Set which nuclides and temperatures are to be read on +// mgInterface through global data +void setMgInterfaceNuclidesAndTemps(); + +// After macro XS have been read, materials can be marked as fissionable +void markFissionableMgxsMaterials(); + //============================================================================== // Mgxs tracking/transport/tallying interface methods //============================================================================== diff --git a/include/openmc/xsdata.h b/include/openmc/xsdata.h index a8aed5d17c..f7ba802b87 100644 --- a/include/openmc/xsdata.h +++ b/include/openmc/xsdata.h @@ -22,6 +22,9 @@ namespace openmc { class XsData { private: + //! Number of energy and delayed neutron groups + size_t n_g, n_dg; + //! \brief Reads scattering data from the HDF5 file void scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, @@ -98,7 +101,10 @@ 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(bool fissionable, int scatter_format, int n_pol, int n_azi); + //! @param n_groups Number of energy groups. + //! @param n_d_groups Number of delayed neutron groups. + XsData(bool fissionable, int scatter_format, int n_pol, int n_azi, + size_t n_groups, size_t n_d_groups); //! \brief Loads the XsData object from the HDF5 file //! diff --git a/src/cross_sections.cpp b/src/cross_sections.cpp index e5315b9fd5..84596f6151 100644 --- a/src/cross_sections.cpp +++ b/src/cross_sections.cpp @@ -155,7 +155,8 @@ void read_cross_sections_xml() if (settings::run_CE) { read_ce_cross_sections_xml(); } else { - data::mgInterface.read_mg_cross_sections_header(); + data::mgInterface.readHeader(settings::path_cross_sections); + putMgxsHeaderDataToGlobals(); } // Establish mapping between (type, material) and index in libraries diff --git a/src/initialize.cpp b/src/initialize.cpp index 203b5fe912..2d40a43924 100644 --- a/src/initialize.cpp +++ b/src/initialize.cpp @@ -260,8 +260,9 @@ void read_input_xml() read_ce_cross_sections(nuc_temps, thermal_temps); } else { // Create material macroscopic data for MGXS - data::mgInterface.init(settings::path_cross_sections); - data::mgInterface.create_macro_xs(); + setMgInterfaceNuclidesAndTemps(); + data::mgInterface.init(); + markFissionableMgxsMaterials(); } simulation::time_read_xs.stop(); } diff --git a/src/mgxs.cpp b/src/mgxs.cpp index 83adfd0b40..84381fe77d 100644 --- a/src/mgxs.cpp +++ b/src/mgxs.cpp @@ -288,7 +288,8 @@ Mgxs::Mgxs(hid_t xs_id, const std::vector& temperature, // Load the more specific XsData information for (int t = 0; t < temps_to_read.size(); t++) { - xs[t] = XsData(fissionable, final_scatter_format, n_pol, n_azi); + xs[t] = XsData(fissionable, final_scatter_format, n_pol, n_azi, + num_groups, num_delayed_groups); // 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()); @@ -332,7 +333,7 @@ Mgxs::Mgxs(const std::string& in_name, const std::vector& mat_kTs, // Create the xs data for each temperature for (int t = 0; t < mat_kTs.size(); t++) { xs[t] = XsData(in_fissionable, in_scatter_format, in_polar.size(), - in_azimuthal.size()); + in_azimuthal.size(), num_groups, num_delayed_groups); // Find the right temperature index to use double temp_desired = mat_kTs[t]; diff --git a/src/mgxs_interface.cpp b/src/mgxs_interface.cpp index befa946828..861926c6b1 100644 --- a/src/mgxs_interface.cpp +++ b/src/mgxs_interface.cpp @@ -26,29 +26,48 @@ namespace data { MgxsInterface mgInterface; } -MgxsInterface::MgxsInterface(const std::string& path_cross_sections) +MgxsInterface::MgxsInterface(const std::string& path_cross_sections, + const std::vector xs_to_read, + const std::vector> xs_temps) { - init(path_cross_sections); + readHeader(path_cross_sections); + setNuclidesToRead(xs_to_read); + setNuclideTemperaturesToRead(xs_temps); + init(); } -void MgxsInterface::init(const std::string& path_cross_sections) +// Should these perhaps unnecessary setters be lumped into one? +void MgxsInterface::setNuclidesToRead(std::vector arg_xs_to_read) +{ + // Check to remove all duplicates + xs_to_read = arg_xs_to_read; +} +void MgxsInterface::setNuclideTemperaturesToRead(std::vector> xs_temps) { + xs_temps_to_read = xs_temps; + if (xs_to_read.size() != xs_temps.size()) + fatal_error("The list of macro XS temperatures to read does not " + "correspond in length to the number of XS names. "); +} + +void MgxsInterface::init() +{ + + // Check that at least some data was set to be read + if (xs_to_read.size() == 0) + warning("No MGXS nuclides were set to be read."); + // Check if MGXS Library exists - if (!file_exists(path_cross_sections)) { + if (!file_exists(cross_sections_path)) { // Could not find MGXS Library file - fatal_error("Cross sections HDF5 file '" + path_cross_sections + + fatal_error("Cross sections HDF5 file '" + cross_sections_path + "' 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(path_cross_sections, 'r'); + hid_t file_id = file_open(cross_sections_path, 'r'); // Read filetype std::string type; @@ -68,32 +87,12 @@ void MgxsInterface::init(const std::string& path_cross_sections) // ========================================================================== // 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(file_id, name, nuc_temps[i_nuc]); - already_read.insert(name); - } - - if (nuclides_MG[i_nuc].fissionable) { - mat->fissionable_ = true; - } - } - } + for (unsigned i_nuc=0; i_nuc> MgxsInterface::get_mat_kTs() //============================================================================== -void MgxsInterface::read_mg_cross_sections_header() +void MgxsInterface::readHeader(const std::string& path_cross_sections) { + // Save name of HDF5 file to be read to struct data + cross_sections_path = path_cross_sections; + // Check if MGXS Library exists - if (!file_exists(settings::path_cross_sections)) { + if (!file_exists(cross_sections_path)) { // Could not find MGXS Library file - fatal_error("Cross sections HDF5 file '" + settings::path_cross_sections + + fatal_error("Cross sections HDF5 file '" + cross_sections_path + "' does not exist."); } write_message("Reading cross sections HDF5 file...", 5); // Open file for reading - hid_t file_id = file_open(settings::path_cross_sections, 'r', true); + hid_t file_id = file_open(cross_sections_path, 'r', true); ensure_exists(file_id, "energy_groups", true); read_attribute(file_id, "energy_groups", num_energy_groups); @@ -210,35 +212,84 @@ void MgxsInterface::read_mg_cross_sections_header() // Reverse energy bins std::copy(rev_energy_bins.crbegin(), rev_energy_bins.crend(), - std::back_inserter(data::mgInterface.energy_bins)); + std::back_inserter(energy_bins)); // Create average energies - for (int i = 0; i < data::mgInterface.energy_bins.size() - 1; ++i) { - data::mgInterface.energy_bin_avg.push_back(0.5* - (data::mgInterface.energy_bins[i] + data::mgInterface.energy_bins[i+1])); + for (int i = 0; i < energy_bins.size() - 1; ++i) { + energy_bin_avg.push_back(0.5* + (energy_bins[i] + energy_bins[i+1])); } // Add entries into libraries for MG data - auto names = group_names(file_id); - if (names.empty()) { + xs_names = group_names(file_id); + if (xs_names.empty()) { fatal_error("At least one MGXS data set must be present in mgxs " "library file!"); } - for (auto& name : names) { - Library lib {}; - lib.type_ = Library::Type::neutron; - lib.materials_.push_back(name); - data::libraries.push_back(lib); - } + // Close MGXS HDF5 file + file_close(file_id); +} +void putMgxsHeaderDataToGlobals() +{ // Get the minimum and maximum energies int neutron = static_cast(Particle::Type::neutron); data::energy_min[neutron] = data::mgInterface.energy_bins.back(); data::energy_max[neutron] = data::mgInterface.energy_bins.front(); - // Close MGXS HDF5 file - file_close(file_id); + // Save available XS names to library list, so that when + // materials are read, the specified mgxs can be confirmed + // as present + for (auto& name : data::mgInterface.xs_names) { + Library lib {}; + lib.type_ = Library::Type::neutron; + lib.materials_.push_back(name); + data::libraries.push_back(lib); + } +} + +void setMgInterfaceNuclidesAndTemps() +{ + // Get temperatures from global data + std::vector> these_nuc_temps(data::nuclide_map.size()); + std::vector> dummy; + get_temperatures(these_nuc_temps, dummy); + + // Build vector of nuclide names which are to be read + std::vector nuclide_names(data::nuclide_map.size()); + for (const auto& kv : data::nuclide_map) { + nuclide_names[kv.second] = kv.first; + } + + std::unordered_set already_read; + + // 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()) { + data::mgInterface.xs_to_read.push_back(name); + data::mgInterface.xs_temps_to_read.push_back(these_nuc_temps[i_nuc]); + // DBG + std::cout << these_nuc_temps[i_nuc][0] << std::endl; + already_read.insert(name); + } + } + } +} + +void markFissionableMgxsMaterials() +{ + // Loop over all files + for (const auto& mat : model::materials) { + for (int i_nuc : mat->nuclide_) { + if (data::mgInterface.nuclides_MG[i_nuc].fissionable) { + mat->fissionable_ = true; + } + } + } } //============================================================================== diff --git a/src/xsdata.cpp b/src/xsdata.cpp index ab6fb52d08..0ee0bd0b02 100644 --- a/src/xsdata.cpp +++ b/src/xsdata.cpp @@ -24,11 +24,12 @@ namespace openmc { // XsData class methods //============================================================================== -XsData::XsData(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_groups, size_t n_d_groups) : + n_g(n_groups), + n_dg(n_d_groups) { size_t n_ang = n_pol * n_azi; - size_t n_dg = data::mgInterface.num_delayed_groups; - size_t n_g = data::mgInterface.num_energy_groups; // check to make sure scatter format is OK before we allocate if (scatter_format != ANGLE_HISTOGRAM && scatter_format != ANGLE_TABULAR && @@ -127,9 +128,6 @@ XsData::fission_vector_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, { // Data is provided as nu-fission and chi with a beta for delayed info - size_t n_g = data::mgInterface.num_energy_groups; - size_t n_dg = data::mgInterface.num_delayed_groups; - // Get chi xt::xtensor temp_chi({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "chi", temp_chi, true); @@ -182,9 +180,6 @@ 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::mgInterface.num_energy_groups; - size_t n_dg = data::mgInterface.num_delayed_groups; - // Get chi-prompt xt::xtensor temp_chi_p({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "chi-prompt", temp_chi_p, true); @@ -218,8 +213,6 @@ 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::mgInterface.num_energy_groups; - // Get chi xt::xtensor temp_chi({n_ang, n_g}, 0.); read_nd_vector(xsdata_grp, "chi", temp_chi, true); @@ -241,9 +234,6 @@ XsData::fission_matrix_beta_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_is { // Data is provided as nu-fission and chi with a beta for delayed info - size_t n_g = data::mgInterface.num_energy_groups; - size_t n_dg = data::mgInterface.num_delayed_groups; - // Get nu-fission matrix xt::xtensor temp_matrix({n_ang, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); @@ -319,9 +309,6 @@ 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::mgInterface.num_energy_groups; - size_t n_dg = data::mgInterface.num_delayed_groups; - // Get the prompt nu-fission matrix xt::xtensor temp_matrix_p({n_ang, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "prompt-nu-fission", temp_matrix_p, true); @@ -353,8 +340,6 @@ 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::mgInterface.num_energy_groups; - // Get nu-fission matrix xt::xtensor temp_matrix({n_ang, n_g, n_g}, 0.); read_nd_vector(xsdata_grp, "nu-fission", temp_matrix, true); @@ -381,7 +366,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic) // 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 (data::mgInterface.num_delayed_groups == 0) { + if (n_dg == 0) { fission_vector_no_delayed_from_hdf5(xsdata_grp, n_ang); } else { if (object_exists(xsdata_grp, "beta")) { @@ -391,7 +376,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic) } } } else { - if (data::mgInterface.num_delayed_groups == 0) { + if (n_dg == 0) { fission_matrix_no_delayed_from_hdf5(xsdata_grp, n_ang); } else { if (object_exists(xsdata_grp, "beta")) { @@ -403,7 +388,7 @@ XsData::fission_from_hdf5(hid_t xsdata_grp, size_t n_ang, bool is_isotropic) } // Combine prompt_nu_fission and delayed_nu_fission into nu_fission - if (data::mgInterface.num_delayed_groups == 0) { + if (n_dg == 0) { nu_fission = prompt_nu_fission; } else { nu_fission = prompt_nu_fission + xt::sum(delayed_nu_fission, {1}); @@ -422,7 +407,6 @@ XsData::scatter_from_hdf5(hid_t xsdata_grp, size_t n_ang, hid_t scatt_grp = open_group(xsdata_grp, "scatter_data"); // Get the outgoing group boundary indices - size_t n_g = data::mgInterface.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, n_g}, 0.);