Mgxs may now be used in external linked programs

This commit is contained in:
Gavin Ridley 2019-11-07 13:27:16 -05:00
parent 9c06a128d9
commit a2bcb07e03
8 changed files with 163 additions and 89 deletions

View file

@ -345,7 +345,7 @@ read_dataset(hid_t obj_id, const char* name, Position& r, bool indep=false)
}
template <typename T, std::size_t N>
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<T, N>& 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 <typename T, std::size_t N>
void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
inline void read_nd_vector(hid_t obj_id, const char* name, xt::xtensor<T, N>& result,
bool must_have=false)
{
if (object_exists(obj_id, name)) {

View file

@ -20,6 +20,14 @@ struct MgxsInterface
int num_energy_groups;
int num_delayed_groups;
// List of available names in the HDF5 file
std::vector<std::string> xs_names;
std::vector<std::string> xs_to_read;
std::vector<std::vector<double>> xs_temps_to_read;
// Name of the HDF5 file which contains mgxs
std::string cross_sections_path;
std::vector<Mgxs> nuclides_MG;
std::vector<Mgxs> macro_xs;
@ -27,12 +35,20 @@ struct MgxsInterface
std::vector<double> energy_bin_avg;
std::vector<double> rev_energy_bins;
// temperatues of each available nuclide
std::vector<std::vector<double>> 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<std::string> xs_to_read,
const std::vector<std::vector<double>> xs_temps);
void setNuclidesToRead(std::vector<std::string> arg_xs_to_read);
void setNuclideTemperaturesToRead(std::vector<std::vector<double>> 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<double>& temperature);
@ -41,13 +57,27 @@ struct MgxsInterface
std::vector<std::vector<double>> 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
//==============================================================================

View file

@ -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
//!

View file

@ -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

View file

@ -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();
}

View file

@ -288,7 +288,8 @@ Mgxs::Mgxs(hid_t xs_id, const std::vector<double>& 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<double>& 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];

View file

@ -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<std::string> xs_to_read,
const std::vector<std::vector<double>> 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<std::string> arg_xs_to_read)
{
// Check to remove all duplicates
xs_to_read = arg_xs_to_read;
}
void MgxsInterface::setNuclideTemperaturesToRead(std::vector<std::vector<double>> 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<std::vector<double>> nuc_temps(data::nuclide_map.size());
std::vector<std::vector<double>> 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<std::string> already_read;
// Build vector of nuclide names
std::vector<std::string> 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<xs_to_read.size(); ++i_nuc)
add_mgxs(file_id, xs_to_read[i_nuc], xs_temps_to_read[i_nuc]);
file_close(file_id);
create_macro_xs();
}
//==============================================================================
@ -113,7 +112,7 @@ MgxsInterface::add_mgxs(hid_t file_id, const std::string& name,
+ "provided MGXS Library");
}
nuclides_MG.emplace_back(xs_grp, temperature, num_energy_groups,
nuclides_MG.emplace_back(xsgavin.keith.ridley@gmail.com_grp, temperature, num_energy_groups,
num_delayed_groups);
close_group(xs_grp);
}
@ -183,18 +182,21 @@ std::vector<std::vector<double>> 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<int>(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<std::vector<double>> these_nuc_temps(data::nuclide_map.size());
std::vector<std::vector<double>> dummy;
get_temperatures(these_nuc_temps, dummy);
// Build vector of nuclide names which are to be read
std::vector<std::string> nuclide_names(data::nuclide_map.size());
for (const auto& kv : data::nuclide_map) {
nuclide_names[kv.second] = kv.first;
}
std::unordered_set<std::string> 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;
}
}
}
}
//==============================================================================

View file

@ -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<double, 2> 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<double, 2> 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<double, 2> 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<double, 3> 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<double, 3> 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<double, 3> 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<int, 2> gmin({n_ang, n_g}, 0.);
read_nd_vector(scatt_grp, "g_min", gmin, true);
xt::xtensor<int, 2> gmax({n_ang, n_g}, 0.);