crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 21:55:41 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Move openmc_load_nuclide, openmc_material_add_nuclide, and read_mgxs to C++

Browse source
This commit is contained in:
Paul Romano 2019-01-24 07:26:44 -06:00
parent 9a2298f66d
commit 06751d5373
23 changed files with 552 additions and 730 deletions
Download patch file Download diff file Expand all files Collapse all files

2
include/openmc/capi.h
View file

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

8
include/openmc/cross_sections.h
View file

@ -59,6 +59,14 @@ extern std::map<LibraryKey, std::size_t> 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<std::vector<double>>& nuc_temps,
const std::vector<std::vector<double>>& thermal_temps);
//! Read cross_sections.xml and populate data libraries
void read_ce_cross_sections_xml();

7
include/openmc/hdf5_interface.h
View file

@ -137,6 +137,13 @@ void read_attribute(hid_t obj_id, const char* name, T& buffer)
read_attr(obj_id, name, H5TypeMap<T>::type_id, &buffer);
}
// array version
template<typename T, std::size_t N> inline void
read_attribute(hid_t obj_id, const char* name, std::array<T, N>& buffer)
{
read_attr(obj_id, name, H5TypeMap<T>::type_id, buffer.data());
}
// vector version
template<typename T>
void read_attribute(hid_t obj_id, const char* name, std::vector<T>& vec)

3
include/openmc/material.h
View file

@ -93,9 +93,6 @@ private:
void calculate_photon_xs(const Particle& p) const;
};
void read_ce_cross_sections(const std::vector<std::vector<double>>& nuc_temps,
const std::vector<std::vector<double>>& thermal_temps);
//==============================================================================
// Fortran compatibility
//==============================================================================

35
include/openmc/mgxs.h
View file

@ -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<double>& 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<double>& in_polar, const std::vector<double>& 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<double>& temperature,
double tolerance, std::vector<int>& temps_to_read, int& order_dim,
int& method);
metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
std::vector<int>& 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<double>& 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<double>& 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<double>& mat_kTs,
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
double tolerance, int& method);
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities);
//! \brief Provides a cross section value given certain parameters
//!

9
include/openmc/mgxs_interface.h
View file

@ -20,6 +20,7 @@ namespace data {
extern std::vector<Mgxs> nuclides_MG;
extern std::vector<Mgxs> macro_xs;
extern "C" int num_energy_groups;
extern "C" int num_delayed_groups;
extern std::vector<double> energy_bins;
extern std::vector<double> energy_bin_avg;
extern std::vector<double> rev_energy_bins;
@ -30,11 +31,11 @@ extern std::vector<double> 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<double>& temperature);
extern "C" bool
query_fissionable_c(int n_nuclides, const int i_nuclides[]);

3
include/openmc/nuclide.h
View file

@ -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<double>& 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);
//==============================================================================

6
include/openmc/scattdata.h
View file

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

38
include/openmc/xsdata.h
View file

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

220
src/cross_sections.cpp
View file

@ -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 <cstdlib> // for getenv
#include <unordered_set>
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<std::vector<double>>& nuc_temps,
const std::vector<std::vector<double>>& thermal_temps)
{
std::unordered_set<std::string> 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<std::string> nuclide_names(data::nuclide_map.size());
std::vector<std::string> 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<Nuclide>(
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<int>(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<int>(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<ThermalScattering>(
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<int>(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<int>(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<int>(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

49
src/geometry_header.F90
View file

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

3
src/initialize.cpp
View file

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

6
src/input_xml.F90
View file

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

261
src/material.cpp
View file

@ -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<std::vector<double>>& nuc_temps,
const std::vector<std::vector<double>>& thermal_temps)
{
std::unordered_set<std::string> 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<std::string> nuclide_names(data::nuclide_map.size());
std::vector<std::string> 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<Nuclide>(
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<int>(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<int>(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<ThermalScattering>(
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<int>(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<int>(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<int>(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<double, 1> atom_density = xt::zeros<double>({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)
{

67
src/material_header.F90
View file

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

54
src/mgxs.cpp
View file

@ -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<Mgxs> macro_xs;
void
Mgxs::init(const std::string& in_name, double in_awr,
const std::vector<double>& 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<double>& in_polar, const std::vector<double>& 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<double>& temperature,
double tolerance, std::vector<int>& temps_to_read, int& order_dim, int& method)
Mgxs::metadata_from_hdf5(hid_t xs_id, const std::vector<double>& temperature,
std::vector<int>& 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<double>& 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<double>& temperature)
{
// Call generic data gathering routine (will populate the metadata)
int order_data;
std::vector<int> 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<double>& mat_kTs,
const std::vector<Mgxs*>& micros, const std::vector<double>& atom_densities,
double tolerance, int& method)
const std::vector<Mgxs*>& micros, const std::vector<double>& 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<double>& mat_kTs,
std::vector<double> 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<double>& mat_kTs,
std::vector<int> micro_t(micros.size(), 0);
std::vector<double> 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<double>& 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<double> interp(micros.size());
std::vector<double> temp_indices(micros.size());

85
src/mgxs_data.F90
View file

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

30
src/mgxs_interface.F90
View file

@ -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(:)

97
src/mgxs_interface.cpp
View file

@ -1,10 +1,17 @@
#include "openmc/mgxs_interface.h"
#include <string>
#include <unordered_set>
#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<double> 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<double> 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<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(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<int, 2> 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<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_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<double>& 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();
}
}

56
src/nuclide.cpp
View file

@ -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<double>& 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<int> 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<double> temperature;
int i_nuclide = data::nuclides.size();
data::nuclides.push_back(std::make_unique<Nuclide>(
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
//==============================================================================

97
src/nuclide_header.F90
View file

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

5
src/scattdata.cpp
View file

@ -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<ScattData*>& 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

141
src/xsdata.cpp
View file

@ -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<size_t> shape = {n_ang, energy_groups};
std::vector<size_t> shape {n_ang, n_g};
total = xt::zeros<double>(shape);
absorption = xt::zeros<double>(shape);
inverse_velocity = xt::zeros<double>(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<double>(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<double>(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<double>(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<double>(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<double, 2> temp_chi({n_ang, energy_groups}, 0.);
xt::xtensor<double, 2> 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<double, 2> temp_nufiss({n_ang, energy_groups}, 0.);
xt::xtensor<double, 2> 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<double, 2> temp_beta({n_ang, delayed_groups}, 0.);
xt::xtensor<double, 2> 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<double, 3> temp_beta({n_ang, delayed_groups, energy_groups},
0.);
xt::xtensor<double, 3> 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<double, 2> temp_chi_p({n_ang, energy_groups}, 0.);
xt::xtensor<double, 2> 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<double, 3> temp_chi_d({n_ang, delayed_groups, energy_groups}, 0.);
xt::xtensor<double, 3> 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<double, 2> temp_chi({n_ang, energy_groups}, 0.);
xt::xtensor<double, 2> 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<double, 3> temp_matrix({n_ang, energy_groups, energy_groups}, 0.);
xt::xtensor<double, 3> 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<double, 2> temp_beta({n_ang, delayed_groups}, 0.);
xt::xtensor<double, 2> temp_beta({n_ang, n_dg}, 0.);
read_nd_vector(xsdata_grp, "beta", temp_beta, true);
xt::xtensor<double, 1> 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<double, 3> temp_beta({n_ang, delayed_groups, energy_groups}, 0.);
xt::xtensor<double, 3> temp_beta({n_ang, n_dg, n_g}, 0.);
read_nd_vector(xsdata_grp, "beta", temp_beta, true);
xt::xtensor<double, 2> temp_beta_sum({n_ang, energy_groups}, 0.);
xt::xtensor<double, 2> 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<double, 3> temp_matrix_p({n_ang, energy_groups, energy_groups}, 0.);
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);
// 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<double, 4> temp_matrix_d({n_ang, delayed_groups, energy_groups,
energy_groups}, 0.);
xt::xtensor<double, 4> 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<double, 3> temp_matrix({n_ang, energy_groups, energy_groups}, 0.);
xt::xtensor<double, 3> 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<int, 2> gmin({n_ang, energy_groups}, 0.);
size_t n_g = data::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, energy_groups}, 0.);
xt::xtensor<int, 2> 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<double, 1> 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<ScattDataTabular*>(scatter[a].get()),
legendre_to_tabular_points);
*static_cast<ScattDataTabular*>(scatter[a].get()));
scatter_format = final_scatter_format;
}