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

Move read_ce_cross_sections to C++

Browse source
This commit is contained in:
Paul Romano 2019-01-23 10:35:28 -06:00
parent fe9f774d31
commit 9a2298f66d
27 changed files with 454 additions and 497 deletions
Download patch file Download diff file Expand all files Collapse all files

1
CMakeLists.txt
View file

@ -335,7 +335,6 @@ add_library(libopenmc SHARED
src/random_lcg.F90
src/reaction_header.F90
src/relaxng
src/sab_header.F90
src/set_header.F90
src/settings.F90
src/simulation_header.F90

2
include/openmc/cross_sections.h
View file

@ -57,7 +57,7 @@ extern std::map<LibraryKey, std::size_t> library_map;
//! Read cross sections file (either XML or multigroup H5) and populate data
//! libraries
extern "C" void read_cross_sections_xml();
void read_cross_sections_xml();
//! Read cross_sections.xml and populate data libraries
void read_ce_cross_sections_xml();

33
include/openmc/geometry_aux.h
View file

@ -6,7 +6,7 @@
#include <cstdint>
#include <string>
#include <vector>
namespace openmc {
@ -14,13 +14,38 @@ namespace openmc {
//! Replace Universe, Lattice, and Material IDs with indices.
//==============================================================================
extern "C" void adjust_indices();
void adjust_indices();
//==============================================================================
//! Assign defaults to cells with undefined temperatures.
//==============================================================================
extern "C" void assign_temperatures();
void assign_temperatures();
//==============================================================================
//! \brief Obtain a list of temperatures that each nuclide/thermal scattering
//! 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)
//!
//! \param[out] nuc_temps Vector of temperatures for each nuclide
//! \param[out] thermal_temps Vector of tempratures for each thermal scattering
//! table
//==============================================================================
void get_temperatures(std::vector<std::vector<double>>& nuc_temps,
std::vector<std::vector<double>>& thermal_temps);
//==============================================================================
//! \brief Perform final setup for geometry
//!
//! \param[out] nuc_temps Vector of temperatures for each nuclide
//! \param[out] thermal_temps Vector of tempratures for each thermal scattering
//! table
//==============================================================================
void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
std::vector<std::vector<double>>& thermal_temps);
//==============================================================================
//! Figure out which Universe is the root universe.
@ -36,7 +61,7 @@ extern "C" int32_t find_root_universe();
//! Populate all data structures needed for distribcells.
//==============================================================================
extern "C" void prepare_distribcell();
void prepare_distribcell();
//==============================================================================
//! Recursively search through the geometry and count cell instances.

1
include/openmc/initialize.h
View file

@ -11,6 +11,7 @@ int parse_command_line(int argc, char* argv[]);
#ifdef OPENMC_MPI
void initialize_mpi(MPI_Comm intracomm);
#endif
void read_input_xml();
}

3
include/openmc/material.h
View file

@ -93,6 +93,9 @@ 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
//==============================================================================

6
include/openmc/nuclide.h
View file

@ -169,6 +169,12 @@ private:
static int XS_PHOTON_PROD;
};
//==============================================================================
// Non-member functions
//==============================================================================
void check_data_version(hid_t file_id);
//==============================================================================
// Global variables
//==============================================================================

2
include/openmc/output.h
View file

@ -37,7 +37,7 @@ extern "C" void print_particle(Particle* p);
//! Display plot information.
//==============================================================================
extern "C" void print_plot();
void print_plot();
//==============================================================================
//! Display information regarding cell overlap checking.

2
include/openmc/settings.h
View file

@ -97,7 +97,7 @@ extern double weight_survive; //!< Survival weight after Russian roulette
//! Read settings from XML file
//! \param[in] root XML node for <settings>
extern "C" void read_settings_xml();
void read_settings_xml();
extern "C" void read_settings_xml_f(pugi::xml_node_struct* root_ptr);

1
include/openmc/timer.h
View file

@ -20,6 +20,7 @@ extern Timer time_bank_sendrecv;
extern Timer time_finalize;
extern Timer time_inactive;
extern Timer time_initialize;
extern Timer time_read_xs;
extern Timer time_tallies;
extern Timer time_total;
extern Timer time_transport;

5
src/api.F90
View file

@ -87,6 +87,9 @@ contains
subroutine free_memory_cmfd() bind(C)
end subroutine free_memory_cmfd
subroutine sab_clear() bind(C)
end subroutine
end interface
call free_memory_geometry()
@ -97,7 +100,7 @@ contains
call free_memory_nuclide()
call free_memory_photon()
call free_memory_settings()
call free_memory_sab()
call sab_clear()
call free_memory_source()
call free_memory_mesh()
call free_memory_tally()

3
src/finalize.cpp
View file

@ -16,7 +16,6 @@ using namespace openmc;
// Functions defined in Fortran
extern "C" void free_memory();
extern "C" void reset_timers_f();
int openmc_finalize()
{
@ -25,7 +24,6 @@ int openmc_finalize()
// Reset timers
reset_timers();
reset_timers_f();
// Reset global variables
settings::assume_separate = false;
@ -127,7 +125,6 @@ int openmc_hard_reset()
// Reset all tallies and timers
openmc_reset();
reset_timers();
reset_timers_f();
// Reset total generations and keff guess
simulation::keff = 1.0;

69
src/geometry_aux.cpp
View file

@ -6,6 +6,7 @@
#include "openmc/cell.h"
#include "openmc/constants.h"
#include "openmc/container_util.h"
#include "openmc/error.h"
#include "openmc/geometry.h"
#include "openmc/lattice.h"
@ -112,6 +113,74 @@ assign_temperatures()
//==============================================================================
void
get_temperatures(std::vector<std::vector<double>>& nuc_temps,
std::vector<std::vector<double>>& thermal_temps)
{
for (const auto& cell : model::cells) {
// Skip non-material cells.
if (cell->fill_ != C_NONE) continue;
for (int j = 0; j < cell->material_.size(); ++j) {
// Skip void materials
int i_material = cell->material_[j];
if (i_material == MATERIAL_VOID) continue;
// Get temperature of cell (rounding to nearest integer)
double sqrtkT = cell->sqrtkT_.size() == 1 ?
cell->sqrtkT_[j] : cell->sqrtkT_[0];
double temperature = sqrtkT*sqrtkT / K_BOLTZMANN;
const auto& mat {model::materials[i_material]};
for (const auto& i_nuc : mat->nuclide_) {
// Add temperature if it hasn't already been added
if (!contains(nuc_temps[i_nuc], temperature)) {
nuc_temps[i_nuc].push_back(temperature);
}
}
if (mat->thermal_tables_.size() > 0) {
for (const auto& table : mat->thermal_tables_) {
// Get index in data::thermal_scatt array
int i_sab = table.index_table;
// Add temperature if it hasn't already been added
if (!contains(thermal_temps[i_sab], temperature)) {
thermal_temps[i_sab].push_back(temperature);
}
}
}
}
}
}
//==============================================================================
void finalize_geometry(std::vector<std::vector<double>>& nuc_temps,
std::vector<std::vector<double>>& thermal_temps)
{
// Perform some final operations to set up the geometry
adjust_indices();
count_cell_instances(model::root_universe);
// Assign temperatures to cells that don't have temperatures already assigned
assign_temperatures();
// Determine desired temperatures for each nuclide and S(a,b) table
get_temperatures(nuc_temps, thermal_temps);
// Check to make sure there are not too many nested coordinate levels in the
// geometry since the coordinate list is statically allocated for performance
// reasons
if (maximum_levels(model::root_universe) > MAX_COORD) {
fatal_error("Too many nested coordinate levels in the geometry. "
"Try increasing the maximum number of coordinate levels by "
"providing the CMake -Dmaxcoord= option.");
}
}
//==============================================================================
int32_t
find_root_universe()
{

17
src/geometry_header.F90
View file

@ -182,18 +182,15 @@ contains
! temperatures to read (which may be different if interpolation is used)
!===============================================================================
subroutine get_temperatures(nuc_temps, sab_temps)
subroutine get_temperatures(nuc_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), optional, allocatable, intent(out) :: sab_temps(:)
integer :: i, j, k
integer :: i_nuclide ! index in nuclides array
integer :: i_sab ! index in S(a,b) array
integer :: i_material
real(8) :: temperature ! temperature in Kelvin
allocate(nuc_temps(n_nuclides))
if (present(sab_temps)) allocate(sab_temps(n_sab_tables))
do i = 1, size(cells)
! Skip non-material cells.
@ -222,18 +219,6 @@ contains
call nuc_temps(i_nuclide) % push_back(temperature)
end if
end do NUC_NAMES_LOOP
if (present(sab_temps) .and. mat % n_sab > 0) then
SAB_NAMES_LOOP: do k = 1, size(mat % sab_names)
! Get index in nuc_temps array
i_sab = sab_dict % get(mat % sab_names(k))
! Add temperature if it hasn't already been added
if (find(sab_temps(i_sab), temperature) == -1) then
call sab_temps(i_sab) % push_back(temperature)
end if
end do SAB_NAMES_LOOP
end if
end associate
end do
end do

66
src/initialize.cpp
View file

@ -5,6 +5,7 @@
#include <cstring>
#include <sstream>
#include <string>
#include <vector>
#ifdef _OPENMP
#include <omp.h>
@ -12,20 +13,33 @@
#include "openmc/capi.h"
#include "openmc/constants.h"
#include "openmc/cross_sections.h"
#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/nuclide.h"
#include "openmc/output.h"
#include "openmc/random_lcg.h"
#include "openmc/settings.h"
#include "openmc/simulation.h"
#include "openmc/string_utils.h"
#include "openmc/thermal.h"
#include "openmc/timer.h"
// data/functions from Fortran side
extern "C" void create_macro_xs();
extern "C" void normalize_ao();
extern "C" void print_usage();
extern "C" void print_version();
extern "C" void read_command_line();
extern "C" void read_input_xml();
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();
// Paths to various files
extern "C" {
@ -248,4 +262,54 @@ parse_command_line(int argc, char* argv[])
return 0;
}
void read_input_xml()
{
read_settings_xml();
read_cross_sections_xml();
read_materials_xml();
read_geometry_xml();
// Convert user IDs -> indices, assign temperatures
double_2dvec nuc_temps(data::nuclide_map.size());
double_2dvec thermal_temps(data::thermal_scatt_map.size());
finalize_geometry(nuc_temps, thermal_temps);
if (settings::run_mode != RUN_MODE_PLOTTING) {
simulation::time_read_xs.start();
if (settings::run_CE) {
// Read continuous-energy cross sections
read_ce_cross_sections(nuc_temps, thermal_temps);
} else {
// Create material macroscopic data for MGXS
read_mgxs();
create_macro_xs();
}
simulation::time_read_xs.stop();
}
read_tallies_xml();
// Initialize distribcell_filters
prepare_distribcell();
if (settings::run_mode == RUN_MODE_PLOTTING) {
// Read plots.xml if it exists
read_plots_xml();
if (mpi::master && settings::verbosity >= 5) print_plot();
} else {
// Normalize atom/weight percents
normalize_ao();
// Write summary information
if (mpi::master && settings::output_summary) write_summary();
// Warn if overlap checking is on
if (mpi::master && settings::check_overlaps) {
warning("Cell overlap checking is ON.");
}
}
}
} // namespace openmc

352
src/input_xml.F90
View file

@ -35,7 +35,6 @@ module input_xml
use tally_derivative_header
use tally_filter_header
use tally_filter
use timer_header, only: time_read_xs
use trigger_header
use volume_header
use xml_interface
@ -44,20 +43,11 @@ module input_xml
save
interface
subroutine adjust_indices() bind(C)
end subroutine adjust_indices
subroutine assign_temperatures() bind(C)
end subroutine assign_temperatures
subroutine count_cell_instances(univ_indx) bind(C)
import C_INT32_T
integer(C_INT32_T), intent(in), value :: univ_indx
end subroutine count_cell_instances
subroutine prepare_distribcell() bind(C)
end subroutine prepare_distribcell
subroutine read_surfaces(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
@ -68,9 +58,6 @@ module input_xml
type(C_PTR) :: node_ptr
end subroutine read_cells
subroutine read_cross_sections_xml() bind(C)
end subroutine
subroutine read_lattices(node_ptr) bind(C)
import C_PTR
type(C_PTR) :: node_ptr
@ -100,9 +87,6 @@ module input_xml
type(C_PTR) :: node_ptr
end subroutine read_plots
subroutine print_plot() bind(C)
end subroutine print_plot
subroutine set_particle_energy_bounds(particle, E_min, E_max) bind(C)
import C_INT, C_DOUBLE
integer(C_INT), value :: particle
@ -118,81 +102,6 @@ contains
! geometry, materials, and tallies.
!===============================================================================
subroutine read_input_xml() bind(C)
type(VectorReal), allocatable :: nuc_temps(:) ! List of T to read for each nuclide
type(VectorReal), allocatable :: sab_temps(:) ! List of T to read for each S(a,b)
call read_settings_xml()
call read_cross_sections_xml()
call read_materials_xml()
call read_geometry_xml()
! Convert user IDs -> indices, assign temperatures
call finalize_geometry(nuc_temps, sab_temps)
if (run_mode /= MODE_PLOTTING) then
call time_read_xs % start()
if (run_CE) then
! Read continuous-energy cross sections
call read_ce_cross_sections(nuc_temps, sab_temps)
else
! Create material macroscopic data for MGXS
call read_mgxs()
call create_macro_xs()
end if
call time_read_xs % stop()
end if
call read_tallies_xml()
! Initialize distribcell_filters
call prepare_distribcell()
if (run_mode == MODE_PLOTTING) then
! Read plots.xml if it exists
call read_plots_xml()
if (master .and. verbosity >= 5) call print_plot()
else
! Normalize atom/weight percents
call normalize_ao()
! Write summary information
if (master .and. output_summary) call write_summary()
! Warn if overlap checking is on
if (master .and. check_overlaps) &
call warning("Cell overlap checking is ON.")
end if
end subroutine read_input_xml
subroutine finalize_geometry(nuc_temps, sab_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), optional, allocatable, intent(out) :: sab_temps(:)
! Perform some final operations to set up the geometry
call adjust_indices()
call count_cell_instances(root_universe)
! Assign temperatures to cells that don't have temperatures already assigned
call assign_temperatures()
! Determine desired temperatures for each nuclide and S(a,b) table
call get_temperatures(nuc_temps, sab_temps)
! Check to make sure there are not too many nested coordinate levels in the
! geometry since the coordinate list is statically allocated for performance
! reasons
if (maximum_levels(root_universe) > MAX_COORD) then
call fatal_error("Too many nested coordinate levels in the geometry. &
&Try increasing the maximum number of coordinate levels by &
&providing the CMake -Dmaxcoord= option.")
end if
end subroutine finalize_geometry
!===============================================================================
! 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
@ -281,7 +190,7 @@ contains
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_geometry_xml()
subroutine read_geometry_xml() bind(C)
integer :: i, n
integer :: univ_id
@ -483,12 +392,11 @@ contains
end do
end subroutine allocate_cells
subroutine read_materials_xml()
subroutine read_materials_xml() bind(C)
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: k ! loop index
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: index_nuclide ! index in nuclides
integer :: index_element ! index in elements
integer :: index_sab ! index in sab_tables
@ -510,12 +418,10 @@ contains
type(XMLNode) :: node_mat
type(XMLNode) :: node_dens
type(XMLNode) :: node_nuc
type(XMLNode) :: node_sab
type(XMLNode), allocatable :: node_mat_list(:)
type(XMLNode), allocatable :: node_nuc_list(:)
type(XMLNode), allocatable :: node_ele_list(:)
type(XMLNode), allocatable :: node_macro_list(:)
type(XMLNode), allocatable :: node_sab_list(:)
! Display output message
call write_message("Reading materials XML file...", 5)
@ -829,63 +735,6 @@ contains
call densities % clear()
call list_iso_lab % clear()
! =======================================================================
! READ AND PARSE <sab> TAG FOR S(a,b) DATA
if (run_CE) then
! Get pointer list to XML <sab>
call get_node_list(node_mat, "sab", node_sab_list)
n_sab = size(node_sab_list)
if (n_sab > 0) then
! Set number of S(a,b) tables
mat % n_sab = n_sab
! Allocate names and indices for nuclides and tables -- for now we
! allocate these as the number of S(a,b) tables listed. Since a single
! table might apply to multiple nuclides, they are resized later if a
! table is indeed applied to multiple nuclides.
allocate(mat % sab_names(n_sab))
allocate(mat % i_sab_tables(n_sab))
allocate(mat % sab_fracs(n_sab))
do j = 1, n_sab
! Get pointer to S(a,b) table
node_sab = node_sab_list(j)
! Determine name of S(a,b) table
if (.not. check_for_node(node_sab, "name")) then
call fatal_error("Need to specify <name> for S(a,b) table.")
end if
call get_node_value(node_sab, "name", name)
name = trim(name)
mat % sab_names(j) = name
! Read the fraction of nuclei affected by this S(a,b) table
if (check_for_node(node_sab, "fraction")) then
call get_node_value(node_sab, "fraction", mat % sab_fracs(j))
else
mat % sab_fracs(j) = ONE
end if
! Check that this nuclide is listed in the cross_sections.xml file
if (.not. library_present(LIBRARY_THERMAL, name)) then
call fatal_error("Could not find S(a,b) table " // trim(name) &
// " in cross_sections.xml file!")
end if
! If this S(a,b) table hasn't been encountered yet, we need to add its
! name and alias to the sab_dict
if (.not. sab_dict % has(name)) then
index_sab = index_sab + 1
mat % i_sab_tables(j) = index_sab
call sab_dict % set(name, index_sab)
else
mat % i_sab_tables(j) = sab_dict % get(name)
end if
end do
end if
end if
! Add material to dictionary
call material_dict % set(mat % id(), i)
end do
@ -893,7 +742,7 @@ contains
! Set total number of nuclides and S(a,b) tables
n_nuclides = index_nuclide
n_elements = index_element
n_sab_tables = index_sab
allocate(nuclides(n_nuclides))
! Close materials XML file
call doc % clear()
@ -905,7 +754,7 @@ contains
! for errors and placing properly-formatted data in the right data structures
!===============================================================================
subroutine read_tallies_xml()
subroutine read_tallies_xml() bind(C)
integer :: i ! loop over user-specified tallies
integer :: j ! loop over words
@ -1826,7 +1675,7 @@ contains
! READ_PLOTS_XML reads data from a plots.xml file
!===============================================================================
subroutine read_plots_xml()
subroutine read_plots_xml() bind(C)
logical :: file_exists ! does plots.xml file exist?
character(MAX_LINE_LEN) :: filename ! absolute path to plots.xml
@ -1941,7 +1790,7 @@ contains
! NORMALIZE_AO Normalize the nuclide atom percents
!===============================================================================
subroutine normalize_ao()
subroutine normalize_ao() bind(C)
integer :: i ! index in materials array
integer :: j ! index over nuclides in material
real(8) :: sum_percent ! summation
@ -1959,7 +1808,7 @@ contains
do j = 1, size(mat % nuclide)
! determine atomic weight ratio
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
awr = nuclide_awr(mat % nuclide(j))
else
awr = get_awr_c(mat % nuclide(j))
end if
@ -1984,7 +1833,7 @@ contains
sum_percent = ZERO
do j = 1, mat % n_nuclides
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
awr = nuclide_awr(mat % nuclide(j))
else
awr = get_awr_c(mat % nuclide(j))
end if
@ -2003,7 +1852,7 @@ contains
mat % density_gpcc = ZERO
do j = 1, mat % n_nuclides
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
awr = nuclide_awr(mat % nuclide(j))
else
awr = ONE
end if
@ -2015,192 +1864,13 @@ contains
end subroutine normalize_ao
subroutine read_ce_cross_sections(nuc_temps, sab_temps)
type(VectorReal), intent(in) :: nuc_temps(:)
type(VectorReal), intent(in) :: sab_temps(:)
integer :: i, j
integer :: i_nuclide
integer :: i_element
integer :: i_sab
integer(C_INT) :: n
integer(HID_T) :: file_id
integer(HID_T) :: group_id
real(C_DOUBLE) :: dummy
logical :: mp_found ! if windowed multipole libraries were found
character(MAX_WORD_LEN) :: name
character(MAX_FILE_LEN) :: filename
character(3) :: element
type(SetChar) :: already_read
type(SetChar) :: element_already_read
interface
subroutine photon_from_hdf5(group) bind(C)
import HID_T
integer(HID_T), value :: group
end subroutine
subroutine read_ce_cross_sections_c() bind(C)
end subroutine
end interface
allocate(nuclides(n_nuclides))
! Read cross sections
do i = 1, size(materials)
do j = 1, size(materials(i) % names)
name = materials(i) % names(j)
if (.not. already_read % contains(name)) then
filename = library_path(LIBRARY_NEUTRON, name)
i_nuclide = nuclide_dict % get(name)
call write_message('Reading ' // trim(name) // ' from ' // &
trim(filename), 6)
! 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(i_nuclide) % from_hdf5(group_id, nuc_temps(i_nuclide), &
temperature_method, temperature_tolerance, temperature_range, &
master, i_nuclide)
call close_group(group_id)
call file_close(file_id)
! Determine if minimum/maximum energy for this nuclide is greater/less
! than the previous
if (size(nuclides(i_nuclide) % grid) >= 1) then
energy_min(NEUTRON) = max(energy_min(NEUTRON), &
nuclides(i_nuclide) % grid(1) % energy(1))
energy_max(NEUTRON) = min(energy_max(NEUTRON), nuclides(i_nuclide) % &
grid(1) % energy(size(nuclides(i_nuclide) % grid(1) % energy)))
call set_particle_energy_bounds(NEUTRON, energy_min(NEUTRON), &
energy_max(NEUTRON))
end if
! Add name and alias to dictionary
call already_read % add(name)
! Check if elemental data has been read, if needed
element = name(1:scan(name, '0123456789') - 1)
if (photon_transport) then
if (.not. element_already_read % contains(element)) then
! Read photon interaction data from HDF5 photon library
filename = library_path(LIBRARY_PHOTON, element)
i_element = element_dict % get(element)
call write_message('Reading ' // trim(element) // ' from ' // &
trim(filename), 6)
! Open file and make sure version is sufficient
file_id = file_open(filename, 'r')
call check_data_version(file_id)
! Read element data from HDF5
group_id = open_group(file_id, element)
call photon_from_hdf5(group_id)
call close_group(group_id)
call file_close(file_id)
! Add element to set
call element_already_read % add(element)
end if
end if
! Read multipole file into the appropriate entry on the nuclides array
if (temperature_multipole) call read_multipole_data(i_nuclide)
end if
end do
end do
do i = 1, size(materials)
! Skip materials with no S(a,b) tables
if (.not. allocated(materials(i) % sab_names)) cycle
do j = 1, size(materials(i) % sab_names)
! Get name of S(a,b) table
name = materials(i) % sab_names(j)
if (.not. already_read % contains(name)) then
filename = library_path(LIBRARY_THERMAL, name)
i_sab = sab_dict % get(name)
call write_message('Reading ' // trim(name) // ' from ' // &
trim(filename), 6)
! Open file and make sure version matches
file_id = file_open(filename, 'r')
call check_data_version(file_id)
! Read S(a,b) data from HDF5
group_id = open_group(file_id, name)
n = sab_temps(i_sab) % size()
if (n > 0) then
call sab_from_hdf5(group_id, sab_temps(i_sab) % data(1), n)
else
! In this case, data(1) doesn't exist, so we just pass a dummy value
call sab_from_hdf5(group_id, dummy, n)
end if
call close_group(group_id)
call file_close(file_id)
! Add name to dictionary
call already_read % add(name)
end if
end do
end do
call read_ce_cross_sections_c()
! Show which nuclide results in lowest energy for neutron transport
do i = 1, size(nuclides)
! If a nuclide is present in a material that's not used in the model, its
! grid has not been allocated
if (size(nuclides(i) % grid) > 0) then
if (nuclides(i) % grid(1) % energy(size(nuclides(i) % grid(1) % energy)) &
== energy_max(NEUTRON)) then
call write_message("Maximum neutron transport energy: " // &
trim(to_str(energy_max(NEUTRON))) // " eV for " // &
trim(adjustl(nuclides(i) % name)), 7)
if (master .and. energy_max(NEUTRON) < 20.e6) call warning("Maximum &
&neutron energy is below 20 MeV. This may bias the results.")
exit
end if
end if
end do
! If the user wants multipole, make sure we found a multipole library.
if (temperature_multipole) then
mp_found = .false.
do i = 1, size(nuclides)
if (nuclides(i) % mp_present) then
mp_found = .true.
exit
end if
end do
if (master .and. .not. mp_found) call 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.")
end if
call already_read % clear()
call element_already_read % clear()
end subroutine read_ce_cross_sections
!===============================================================================
! READ_MULTIPOLE_DATA checks for the existence of a multipole library in the
! directory and loads it using multipole_read
!===============================================================================
subroutine read_multipole_data(i_table)
integer, intent(in) :: i_table ! index in nuclides/sab_tables
subroutine read_multipole_data(i_table) bind(C)
integer(C_INT), value :: i_table ! index in nuclides/sab_tables
logical :: file_exists ! Does multipole library exist?
character(7) :: readable ! Is multipole library readable?

186
src/material.cpp
View file

@ -14,6 +14,7 @@
#include "openmc/container_util.h"
#include "openmc/error.h"
#include "openmc/math_functions.h"
#include "openmc/message_passing.h"
#include "openmc/mgxs_interface.h"
#include "openmc/nuclide.h"
#include "openmc/photon.h"
@ -829,8 +830,155 @@ read_materials(pugi::xml_node* node)
}
}
extern "C" void read_ce_cross_sections_c()
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();
}
@ -857,6 +1005,42 @@ extern "C" void read_ce_cross_sections_c()
// 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.");
}
}
}
//==============================================================================

4
src/material_header.F90
View file

@ -295,7 +295,7 @@ contains
do j = 1, m % n_nuclides
k = m % nuclide(j)
if (nuclides(k) % name == name_) then
awr = nuclides(k) % awr
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
@ -326,7 +326,7 @@ contains
m % atom_density(n + 1) = density
m % density = m % density + density
m % density_gpcc = m % density_gpcc + &
density * nuclides(k) % awr * MASS_NEUTRON / N_AVOGADRO
density * nuclide_awr(k) * MASS_NEUTRON / N_AVOGADRO
m % n_nuclides = n + 1
end if
end if

4
src/mgxs_data.F90
View file

@ -24,7 +24,7 @@ contains
! nuclides and sab_tables arrays
!===============================================================================
subroutine read_mgxs()
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
@ -106,7 +106,7 @@ contains
! CREATE_MACRO_XS generates the macroscopic xs from the microscopic input data
!===============================================================================
subroutine create_macro_xs()
subroutine create_macro_xs() bind(C)
integer :: i_mat ! index in materials array
type(Material), pointer :: mat ! current material
type(VectorReal), allocatable :: kTs(:)

31
src/nuclide.cpp
View file

@ -866,6 +866,28 @@ void Nuclide::calculate_urr_xs(int i_temp, double E) const
}
//==============================================================================
// Non-member functions
//==============================================================================
void check_data_version(hid_t file_id)
{
if (attribute_exists(file_id, "version")) {
std::vector<int> version;
read_attribute(file_id, "version", version);
if (version[0] != HDF5_VERSION[0]) {
fatal_error("HDF5 data format uses version " + std::to_string(version[0])
+ "." + std::to_string(version[1]) + " whereas your installation of "
"OpenMC expects version " + std::to_string(HDF5_VERSION[0])
+ ".x data.");
}
} else {
fatal_error("HDF5 data does not indicate a version. Your installation of "
"OpenMC expects version " + std::to_string(HDF5_VERSION[0]) +
".x data.");
}
}
//==============================================================================
// Fortran compatibility functions
//==============================================================================
@ -877,15 +899,10 @@ set_particle_energy_bounds(int particle, double E_min, double E_max)
data::energy_max[particle - 1] = E_max;
}
extern "C" Nuclide* nuclide_from_hdf5_c(hid_t group, const double* temperature, int n)
{
data::nuclides.push_back(std::make_unique<Nuclide>(group, temperature, n,
data::nuclides.size()));
return data::nuclides.back().get();
}
extern "C" void nuclide_init_grid_c(Nuclide* nuc) { nuc->init_grid(); }
extern "C" double nuclide_awr(int i_nuc) { return data::nuclides[i_nuc - 1]->awr_; }
extern "C" Reaction* nuclide_reaction(Nuclide* nuc, int i_rx)
{
return nuc->reactions_[i_rx-1].get();

72
src/nuclide_header.F90
View file

@ -79,7 +79,6 @@ module nuclide_header
type(C_PTR) :: ptr
contains
procedure :: from_hdf5 => nuclide_from_hdf5
procedure :: init_grid => nuclide_init_grid
procedure :: nu => nuclide_nu
procedure, private :: create_derived => nuclide_create_derived
@ -195,6 +194,12 @@ module nuclide_header
logical(C_BOOL) :: b
end function
function nuclide_awr(i_nuc) result(awr) bind(C)
import C_INT, C_DOUBLE
integer(C_INT), value :: i_nuc
real(C_DOUBLE) :: awr
end function
function nuclide_fission_q_prompt(ptr, E) result(q) bind(C)
import C_PTR, C_DOUBLE
type(C_PTR), value :: ptr
@ -238,16 +243,14 @@ contains
ptr = C_LOC(micro_xs(1))
end function
subroutine nuclide_from_hdf5(this, group_id, temperature, method, tolerance, &
minmax, master, i_nuclide)
class(Nuclide), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
type(VectorReal), intent(in), target :: temperature ! list of desired temperatures
integer, intent(inout) :: method
real(8), intent(in) :: tolerance
real(8), intent(in) :: minmax(2) ! range of temperatures
logical(C_BOOL), intent(in) :: master ! if this is the master proc
integer, intent(in) :: i_nuclide ! Nuclide index in nuclides
subroutine nuclide_from_hdf5(group_id, ptr, temps, n, i_nuclide) bind(C)
integer(HID_T), value :: group_id
type(C_PTR), value :: ptr
type(C_PTR), value :: temps
integer(C_INT), value :: n
integer(C_INT), value :: i_nuclide
real(C_DOUBLE), pointer :: temperature(:) ! list of desired temperatures
integer :: i
integer :: i_closest
@ -270,23 +273,11 @@ contains
type(VectorInt) :: MTs
type(VectorInt) :: temps_to_read
interface
function nuclide_from_hdf5_c(group, temperature, n) result(ptr) bind(C)
import HID_T, C_DOUBLE, C_INT, C_PTR
integer(HID_T), value :: group
type(C_PTR), value :: temperature
integer(C_INT), value :: n
type(C_PTR) :: ptr
end function
end interface
! Get array passed
call c_f_pointer(temps, temperature, [n])
! Read data on C++ side
if (temperature % size() > 0) then
this % ptr = nuclide_from_hdf5_c(group_id, C_LOC(temperature % data(1)), &
temperature % size())
else
this % ptr = nuclide_from_hdf5_c(group_id, C_NULL_PTR, 0)
end if
associate (this => nuclides(i_nuclide))
this % ptr = ptr
! Get name of nuclide from group
this % name = get_name(group_id)
@ -311,35 +302,35 @@ contains
call sort(temps_available)
! If only one temperature is available, revert to nearest temperature
if (size(temps_available) == 1 .and. method == TEMPERATURE_INTERPOLATION) then
if (size(temps_available) == 1 .and. temperature_method == TEMPERATURE_INTERPOLATION) then
if (master) then
call warning("Cross sections for " // trim(this % name) // " are only &
&available at one temperature. Reverting to nearest temperature &
&method.")
end if
method = TEMPERATURE_NEAREST
temperature_method = TEMPERATURE_NEAREST
end if
! Determine actual temperatures to read -- start by checking whether a
! temperature range was given, in which case all temperatures in the range
! are loaded irrespective of what temperatures actually appear in the model
if (minmax(2) > ZERO) then
if (temperature_range(2) > ZERO) then
do i = 1, size(temps_available)
temp_actual = temps_available(i)
if (minmax(1) <= temp_actual .and. temp_actual <= minmax(2)) then
if (temperature_range(1) <= temp_actual .and. temp_actual <= temperature_range(2)) then
call temps_to_read % push_back(nint(temp_actual))
end if
end do
end if
select case (method)
select case (temperature_method)
case (TEMPERATURE_NEAREST)
! Find nearest temperatures
do i = 1, temperature % size()
temp_desired = temperature % data(i)
do i = 1, n
temp_desired = temperature(i)
i_closest = minloc(abs(temps_available - temp_desired), dim=1)
temp_actual = temps_available(i_closest)
if (abs(temp_actual - temp_desired) < tolerance) then
if (abs(temp_actual - temp_desired) < temperature_tolerance) then
if (find(temps_to_read, nint(temp_actual)) == -1) then
call temps_to_read % push_back(nint(temp_actual))
@ -362,8 +353,8 @@ contains
case (TEMPERATURE_INTERPOLATION)
! If temperature interpolation or multipole is selected, get a list of
! bounding temperatures for each actual temperature present in the model
TEMP_LOOP: do i = 1, temperature % size()
temp_desired = temperature % data(i)
TEMP_LOOP: do i = 1, n
temp_desired = temperature(i)
do j = 1, size(temps_available) - 1
if (temps_available(j) <= temp_desired .and. &
@ -441,6 +432,7 @@ contains
! Finalize with the nuclide index
this % i_nuclide = i_nuclide
end associate
end subroutine nuclide_from_hdf5
@ -757,9 +749,9 @@ contains
! 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 nuclides(n) % from_hdf5(group_id, temperature, &
! temperature_method, temperature_tolerance, minmax, &
! master, n)
call close_group(group_id)
call file_close(file_id)

2
src/output.F90
View file

@ -322,7 +322,7 @@ contains
! display time elapsed for various sections
write(ou,100) "Total time for initialization", time_initialize_elapsed()
write(ou,100) " Reading cross sections", time_read_xs % elapsed
write(ou,100) " Reading cross sections", time_read_xs_elapsed()
write(ou,100) "Total time in simulation", time_inactive_elapsed() + &
time_active_elapsed()
write(ou,100) " Time in transport only", time_transport_elapsed()

14
src/physics.cpp
View file

@ -438,10 +438,9 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
}
// Get pointers to nuclide/density arrays
int* nuclides;
double* densities;
int n;
openmc_material_get_densities(p->material, &nuclides, &densities, &n);
// TODO: off-by-one
const auto& mat {model::materials[p->material - 1]};
int n = mat->nuclide_.size();
*i_nuc_mat = 0;
double prob = 0.0;
@ -452,10 +451,9 @@ void sample_nuclide(const Particle* p, int mt, int* i_nuclide, int* i_nuc_mat)
fatal_error("Did not sample any nuclide during collision.");
}
// Find atom density
// TODO: off-by-one
*i_nuclide = nuclides[*i_nuc_mat] - 1;
double atom_density = densities[*i_nuc_mat];
// Get atom density
*i_nuclide = mat->nuclide_[*i_nuc_mat];
double atom_density = mat->atom_density_[*i_nuc_mat];
// Determine microscopic cross section
double sigma;

54
src/sab_header.F90
View file

@ -1,54 +0,0 @@
module sab_header
use, intrinsic :: ISO_C_BINDING
use dict_header, only: DictCharInt
use hdf5_interface
implicit none
private
public :: free_memory_sab, sab_from_hdf5, sab_has_nuclide, sab_threshold
! S(a,b) tables
integer, public :: n_sab_tables
type(DictCharInt), public :: sab_dict
interface
subroutine sab_from_hdf5(group_id, temperature, n) bind(C)
import HID_T, C_DOUBLE, C_INT
integer(HID_T), value :: group_id
real(C_DOUBLE), intent(in) :: temperature
integer(C_INT), value :: n
end subroutine
subroutine sab_clear() bind(C)
end subroutine
function sab_has_nuclide(i_sab, name) result(val) bind(C)
import C_INT, C_CHAR, C_BOOL
integer(C_INT), value :: i_sab
character(kind=C_CHAR), intent(in) :: name(*)
logical(C_BOOL) :: val
end function
function sab_threshold(i_sab) result(threshold) bind(C)
import C_INT, C_DOUBLE
integer(C_INT), value :: i_sab
real(C_DOUBLE) :: threshold
end function
end interface
contains
!===============================================================================
! FREE_MEMORY_SAB deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_sab()
n_sab_tables = 0
call sab_clear()
call sab_dict % clear()
end subroutine free_memory_sab
end module sab_header

2
src/state_point.F90
View file

@ -372,7 +372,7 @@ contains
call write_dataset(runtime_group, "total initialization", &
time_initialize_elapsed())
call write_dataset(runtime_group, "reading cross sections", &
time_read_xs % get_value())
time_read_xs_elapsed())
call write_dataset(runtime_group, "simulation", &
time_inactive_elapsed() + time_active_elapsed())
call write_dataset(runtime_group, "transport", &

2
src/summary.F90
View file

@ -26,7 +26,7 @@ contains
! WRITE_SUMMARY
!===============================================================================
subroutine write_summary()
subroutine write_summary() bind(C)
interface
subroutine write_geometry(file_id) bind(C)

3
src/timer.cpp
View file

@ -15,6 +15,7 @@ Timer time_bank_sendrecv;
Timer time_finalize;
Timer time_inactive;
Timer time_initialize;
Timer time_read_xs;
Timer time_tallies;
Timer time_total;
Timer time_transport;
@ -64,6 +65,7 @@ extern "C" double time_bank_sendrecv_elapsed() { return simulation::time_bank_se
extern "C" double time_finalize_elapsed() { return simulation::time_finalize.elapsed(); }
extern "C" double time_inactive_elapsed() { return simulation::time_inactive.elapsed(); }
extern "C" double time_initialize_elapsed() { return simulation::time_initialize.elapsed(); }
extern "C" double time_read_xs_elapsed() { return simulation::time_read_xs.elapsed(); }
extern "C" double time_tallies_elapsed() { return simulation::time_tallies.elapsed(); }
extern "C" double time_total_elapsed() { return simulation::time_total.elapsed(); }
extern "C" double time_transport_elapsed() { return simulation::time_transport.elapsed(); }
@ -81,6 +83,7 @@ void reset_timers()
simulation::time_finalize.reset();
simulation::time_inactive.reset();
simulation::time_initialize.reset();
simulation::time_read_xs.reset();
simulation::time_tallies.reset();
simulation::time_total.reset();
simulation::time_transport.reset();

14
src/timer_header.F90
View file

@ -35,6 +35,10 @@ module timer_header
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_read_xs_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
end function
function time_tallies_elapsed() result(t) bind(C)
import C_DOUBLE
real(C_DOUBLE) :: t
@ -72,8 +76,6 @@ module timer_header
! ============================================================================
! TIMING VARIABLES
type(Timer) :: time_read_xs ! timer for reading cross sections
contains
!===============================================================================
@ -137,12 +139,4 @@ contains
self % elapsed = ZERO
end subroutine timer_reset
!===============================================================================
! RESET_TIMERS resets timers on the Fortran side
!===============================================================================
subroutine reset_timers_f() bind(C)
call time_read_xs % reset()
end subroutine
end module timer_header