crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Updated ipython notebook, fixed minor corrections throughout, and added h5 wrapper functions in hdf5_interface so the various from_hdf5 routines dont have to call low-level functions directly

Browse source
This commit is contained in:
Adam Nelson 2016-09-10 11:28:08 -04:00
parent 61feb8cc6d
commit 255aba1fc0
10 changed files with 206 additions and 196 deletions
Download patch file Download diff file Expand all files Collapse all files

112
docs/source/pythonapi/examples/mgxs-part-iv.ipynb
View file

File diff suppressed because one or more lines are too long

2
docs/source/usersguide/input.rst
View file

@ -737,7 +737,7 @@ coefficients. This element has the following attributes/sub-elements:
*Default*: 33
.. note:: This element is not used in the multi-group :ref:`energy_mode`.
.. note:: This element is only used in the multi-group :ref:`energy_mode`.
.. _temperature_default:

4
openmc/mesh.py
View file

@ -291,6 +291,10 @@ class Mesh(object):
# Would prefer to have the z ranges be the max supported float, but
# these values are apparently different between python and Fortran.
# Choosing a safe and sane default.
# Values of +/-1000 are used here as there seems to be an
# inconsistency between what numpy uses as the max float and what
# Fortran expects for a real(8), so this avoids code complication
# and achieves the same goal.
zplanes = [openmc.ZPlane(z0=-1000.,
boundary_type='reflective'),
openmc.ZPlane(z0=1000.,

6
openmc/mgxs_library.py
View file

@ -14,14 +14,10 @@ if sys.version_info[0] >= 3:
basestring = str
# Supported incoming particle MGXS angular treatment representations
_REPRESENTATION_ISOTROPIC = 1
_REPRESENTATION_ANGLE = 2
_SCATTER_TYPE_TABULAR = 3
_SCATTER_TYPE_LEGENDRE = 4
_SCATTER_TYPE_HISTOGRAM = 5
_REPRESENTATIONS = ['isotropic', 'angle']
_SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
class XSdata(object):
"""A multi-group cross section data set providing all the
multi-group data necessary for a multi-group OpenMC calculation.

22
src/hdf5_interface.F90
View file

@ -102,6 +102,7 @@ module hdf5_interface
public :: write_attribute_string
public :: get_groups
public :: get_datasets
public :: get_name
contains
@ -333,6 +334,27 @@ contains
end subroutine get_datasets
!===============================================================================
! GET_NAME Obtains the name of the current group in group_id
!===============================================================================
function get_name(group_id, name_len_) result(name)
integer(HID_T), intent(in) :: group_id
integer(SIZE_T), optional, intent(in) :: name_len_
character(len=255) :: name ! name of group
integer(SIZE_T) :: name_len, name_file_len
integer :: hdf5_err ! HDF5 error code
if (present(name_len_)) then
name_len = name_len_
else
name_len = 255
end if
call h5iget_name_f(group_id, name, name_len, name_file_len, hdf5_err)
end function get_name
!===============================================================================
! OPEN_GROUP opens an existing HDF5 group
!===============================================================================

169
src/mgxs_header.F90
View file

@ -3,8 +3,7 @@ module mgxs_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5iget_name_f
use h5lt, only: h5ltpath_valid_f
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: find, sort
use constants, only: MAX_FILE_LEN, ZERO, ONE, TWO, PI
@ -212,8 +211,7 @@ module mgxs_header
type(VectorInt), intent(out) :: temps_to_read ! Temperatures to read
integer, intent(out) :: order_dim ! Scattering data order size
integer :: hdf5_err
integer(SIZE_T) :: name_len, name_file_len
integer(SIZE_T) :: name_len
integer(HID_T) :: kT_group
character(MAX_FILE_LEN), allocatable :: dset_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
@ -225,7 +223,7 @@ module mgxs_header
! Get name of dataset from group
name_len = len(this % name)
call h5iget_name_f(xs_id, this % name, name_len, name_file_len, hdf5_err)
this % name = get_name(xs_id, name_len)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -452,8 +450,8 @@ module mgxs_header
! If we have a need* for the fission and kappa-fission x/s, get them
! (*Need is defined as will be using it to tally)
if (get_fiss) then
allocate(xs % fission(groups))
if (check_dataset(xsdata_grp, "fission")) then
allocate(xs % fission(groups))
call read_dataset(xs % fission, xsdata_grp, "fission")
else
call fatal_error("Fission data missing, required due to fission&
@ -461,8 +459,8 @@ module mgxs_header
end if
end if
if (get_kfiss) then
allocate(xs % k_fission(groups))
if (check_dataset(xsdata_grp, "kappa-fission")) then
allocate(xs % k_fission(groups))
call read_dataset(xs % k_fission, xsdata_grp, "kappa-fission")
else
call fatal_error("kappa-fission data missing, required due to &
@ -471,8 +469,8 @@ module mgxs_header
end if
end if
allocate(xs % absorption(groups))
if (check_dataset(xsdata_grp, "absorption")) then
allocate(xs % absorption(groups))
call read_dataset(xs % absorption, xsdata_grp, "absorption")
else
call fatal_error("Must provide absorption!")
@ -581,30 +579,37 @@ module mgxs_header
deallocate(input_scatt)
! Now get the multiplication matrix
! Now use this information to find the length of a container array
! to hold the flattened data
length = 0
do gin = 1, groups
length = length + (gmax(gin) - gmin(gin) + 1)
end do
! Allocate flattened array
allocate(temp_arr(length))
if (.not. check_dataset(scatt_grp, 'multiplicity matrix')) &
call fatal_error("'multiplicity matrix' must be provided")
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing
! to ScattData
allocate(temp_mult(groups))
index = 1
do gin = 1, groups
allocate(temp_mult(gin) % data(gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
temp_mult(gin) % data(gout) = temp_arr(index)
index = index + 1
if (check_dataset(scatt_grp, 'multiplicity matrix')) then
! Now use this information to find the length of a container array
! to hold the flattened data
length = 0
do gin = 1, groups
length = length + (gmax(gin) - gmin(gin) + 1)
end do
end do
deallocate(temp_arr)
! Allocate flattened array
allocate(temp_arr(length))
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing
! to ScattData
allocate(temp_mult(groups))
index = 1
do gin = 1, groups
allocate(temp_mult(gin) % data(gmin(gin):gmax(gin)))
do gout = gmin(gin), gmax(gin)
temp_mult(gin) % data(gout) = temp_arr(index)
index = index + 1
end do
end do
deallocate(temp_arr)
else
! Default to multiplicities of 1.0
allocate(temp_mult(groups))
do gin = 1, groups
allocate(temp_mult(gin) % data(gmin(gin):gmax(gin)))
temp_mult(gin) % data = ONE
end do
end if
! Allocate and initialize our ScattData Object.
if (this % scatter_type == ANGLE_HISTOGRAM) then
@ -641,6 +646,7 @@ module mgxs_header
! Close the groups we have opened and deallocate
call close_group(xsdata_grp)
call close_group(scatt_grp)
deallocate(scatt_coeffs, temp_mult)
end associate ! xs
end do ! Temperatures
@ -689,7 +695,7 @@ module mgxs_header
if (check_dataset(xsdata_grp, "chi")) then
! Chi was provided, that means we need chi and nu-fission vectors
! Get chi
allocate(temp_arr(1 * groups * this % n_azi * this % n_pol))
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "chi")
! Initialize counter for temp_arr
l = 0
@ -920,39 +926,50 @@ module mgxs_header
deallocate(input_scatt)
! Now get the multiplication matrix
! Now use this information to find the length of a container array
! to hold the flattened data
length = 0
do ipol = 1, this % n_pol
do iazi = 1, this % n_azi
do gin = 1, groups
length = length + (gmax(gin, iazi, ipol) - gmin(gin, iazi, ipol) + 1)
end do
end do
end do
! Allocate flattened array
allocate(temp_arr(length))
if (.not. check_dataset(scatt_grp, 'multiplicity matrix')) &
call fatal_error("'multiplicity matrix' must be provided")
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing
! to ScattData
allocate(temp_mult(groups, this % n_azi, this % n_pol))
index = 1
do ipol = 1, this % n_pol
do iazi = 1, this % n_azi
do gin = 1, groups
allocate(temp_mult(gin, iazi, ipol) % data( &
gmin(gin, iazi, ipol):gmax(gin, iazi, ipol)))
do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol)
temp_mult(gin, iazi, ipol) % data(gout) = temp_arr(index)
index = index + 1
if (check_dataset(scatt_grp, 'multiplicity matrix')) then
! Now use this information to find the length of a container array
! to hold the flattened data
length = 0
do ipol = 1, this % n_pol
do iazi = 1, this % n_azi
do gin = 1, groups
length = length + (gmax(gin, iazi, ipol) - gmin(gin, iazi, ipol) + 1)
end do
end do
end do
end do
deallocate(temp_arr)
! Allocate flattened array
allocate(temp_arr(length))
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing
! to ScattData
allocate(temp_mult(groups, this % n_azi, this % n_pol))
index = 1
do ipol = 1, this % n_pol
do iazi = 1, this % n_azi
do gin = 1, groups
allocate(temp_mult(gin, iazi, ipol) % data( &
gmin(gin, iazi, ipol):gmax(gin, iazi, ipol)))
do gout = gmin(gin, iazi, ipol), gmax(gin, iazi, ipol)
temp_mult(gin, iazi, ipol) % data(gout) = temp_arr(index)
index = index + 1
end do
end do
end do
end do
deallocate(temp_arr)
else
allocate(temp_mult(groups, this % n_azi, this % n_pol))
! Default to multiplicities of 1.0
do ipol = 1, this % n_pol
do iazi = 1, this % n_azi
do gin = 1, groups
allocate(temp_mult(gin, iazi, ipol) % data( &
gmin(gin, iazi, ipol):gmax(gin, iazi, ipol)))
temp_mult(gin, iazi, ipol) % data = ONE
end do
end do
end do
end if
! Allocate and initialize our ScattData Object.
allocate(xs % scatter(this % n_azi, this % n_pol))
@ -1013,6 +1030,7 @@ module mgxs_header
! Close the groups we have opened and deallocate
call close_group(xsdata_grp)
call close_group(scatt_grp)
deallocate(scatt_coeffs, temp_mult)
end associate ! xs
end do ! Temperatures
@ -1220,22 +1238,22 @@ module mgxs_header
type is (MgxsIso)
! Perform our operations which depend upon the type
! Add contributions to total, absorption, and fission data (if necessary)
this % xs(t) % total(:) = this % xs(t) % total(:) + &
atom_density * nuc % xs(nuc_t) % total(:)
this % xs(t) % absorption(:) = this % xs(t) % absorption(:) + &
atom_density * nuc % xs(nuc_t) % absorption(:)
this % xs(t) % total = this % xs(t) % total + &
atom_density * nuc % xs(nuc_t) % total
this % xs(t) % absorption = this % xs(t) % absorption + &
atom_density * nuc % xs(nuc_t) % absorption
if (nuc % fissionable) then
this % xs(t) % chi(:, :) = this % xs(t) % chi(:, :) + &
atom_density * nuc % xs(nuc_t) % chi(:, :)
this % xs(t) % nu_fission(:) = this % xs(t) % nu_fission(:) + &
atom_density * nuc % xs(nuc_t) % nu_fission(:)
this % xs(t) % chi = this % xs(t) % chi + &
atom_density * nuc % xs(nuc_t) % chi
this % xs(t) % nu_fission = this % xs(t) % nu_fission + &
atom_density * nuc % xs(nuc_t) % nu_fission
if (allocated(nuc % xs(nuc_t) % fission)) then
this % xs(t) % fission(:) = this % xs(t) % fission(:) + &
atom_density * nuc % xs(nuc_t) % fission(:)
this % xs(t) % fission = this % xs(t) % fission + &
atom_density * nuc % xs(nuc_t) % fission
end if
if (allocated(nuc % xs(nuc_t) % k_fission)) then
this % xs(t) % k_fission(:) = this % xs(t) % k_fission(:) + &
atom_density * nuc % xs(nuc_t) % k_fission(:)
this % xs(t) % k_fission = this % xs(t) % k_fission + &
atom_density * nuc % xs(nuc_t) % k_fission
end if
end if
@ -1547,13 +1565,12 @@ module mgxs_header
end do
end do
! Create the Jagged arrays
! Now create our jagged data from the dense data
call jagged_from_dense_2D(scatt_coeffs(:, :, :, iazi, ipol), &
jagged_scatt)
call jagged_from_dense_1D(temp_mult(:, :, iazi, ipol), &
jagged_mult, gmin, gmax)
! Initialize the ScattData Object
! Initialize the ScattData Object
call this % xs(t) % scatter(iazi, ipol) % obj % init(gmin, &
gmax, jagged_mult, jagged_scatt)

37
src/nuclide_header.F90
View file

@ -3,9 +3,7 @@ module nuclide_header
use, intrinsic :: ISO_FORTRAN_ENV
use, intrinsic :: ISO_C_BINDING
use hdf5, only: HID_T, HSIZE_T, SIZE_T, h5iget_name_f, h5gget_info_f, &
h5lget_name_by_idx_f, H5_INDEX_NAME_F, H5_ITER_INC_F
use h5lt, only: h5ltpath_valid_f
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use algorithm, only: sort, find
use constants
@ -14,7 +12,8 @@ module nuclide_header
use endf_header, only: Function1D, Polynomial, Tabulated1D
use error, only: fatal_error, warning
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape, get_datasets
open_dataset, read_dataset, close_dataset, get_shape, get_datasets, &
check_group, get_name, get_groups
use list_header, only: ListInt
use math, only: evaluate_legendre
use multipole_header, only: MultipoleArray
@ -194,10 +193,6 @@ module nuclide_header
real(8), intent(in) :: tolerance
integer :: i
integer :: storage_type
integer :: max_corder
integer :: n_links
integer :: hdf5_err
integer :: i_closest
integer :: n_temperature
integer(HID_T) :: urr_group, nu_group
@ -208,21 +203,21 @@ module nuclide_header
integer(HID_T) :: total_nu
integer(HID_T) :: fer_group ! fission_energy_release group
integer(HID_T) :: fer_dset
integer(SIZE_T) :: name_len, name_file_len
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: j
integer(HSIZE_T) :: dims(1)
character(MAX_WORD_LEN) :: temp_str
character(MAX_FILE_LEN), allocatable :: dset_names(:)
character(MAX_FILE_LEN), allocatable :: grp_names(:)
real(8), allocatable :: temps_available(:) ! temperatures available
real(8) :: temp_desired
real(8) :: temp_actual
logical :: exists
type(VectorInt) :: MTs
type(VectorInt) :: temps_to_read
! Get name of nuclide from group
name_len = len(this % name)
call h5iget_name_f(group_id, this % name, name_len, name_file_len, hdf5_err)
this % name = get_name(group_id, name_len)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -305,12 +300,10 @@ module nuclide_header
! Get MT values based on group names
rxs_group = open_group(group_id, 'reactions')
call h5gget_info_f(rxs_group, storage_type, n_links, max_corder, hdf5_err)
do j = 0, n_links - 1
call h5lget_name_by_idx_f(rxs_group, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, &
j, temp_str, hdf5_err, name_len)
if (starts_with(temp_str, "reaction_")) then
call MTs % push_back(int(str_to_int(temp_str(10:12))))
call get_groups(rxs_group, grp_names)
do j = 1, size(grp_names)
if (starts_with(grp_names(j), "reaction_")) then
call MTs % push_back(int(str_to_int(grp_names(j)(10:12))))
end if
end do
@ -326,8 +319,7 @@ module nuclide_header
call close_group(rxs_group)
! Read unresolved resonance probability tables if present
call h5ltpath_valid_f(group_id, 'urr', .true., exists, hdf5_err)
if (exists) then
if (check_group(group_id, 'urr')) then
this % urr_present = .true.
allocate(this % urr_data(n_temperature))
@ -368,8 +360,7 @@ module nuclide_header
end if
! Check for nu-total
call h5ltpath_valid_f(group_id, 'total_nu', .true., exists, hdf5_err)
if (exists) then
if (check_group(group_id, 'total_nu')) then
nu_group = open_group(group_id, 'total_nu')
! Read total nu data
@ -388,9 +379,7 @@ module nuclide_header
end if
! Read fission energy release data if present
call h5ltpath_valid_f(group_id, 'fission_energy_release', .true., exists, &
hdf5_err)
if (exists) then
if (check_group(group_id, 'fission_energy_release')) then
fer_group = open_group(group_id, 'fission_energy_release')
! Check to see if this is polynomial or tabulated data

17
src/reaction_header.F90
View file

@ -5,7 +5,7 @@ module reaction_header
use constants, only: MAX_WORD_LEN
use hdf5_interface, only: read_attribute, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_shape
open_dataset, read_dataset, close_dataset, get_shape, get_groups
use product_header, only: ReactionProduct
use stl_vector, only: VectorInt
use string, only: to_str, starts_with
@ -42,17 +42,12 @@ contains
integer :: i
integer :: cm
integer :: n_product
integer :: storage_type
integer :: max_corder
integer :: n_links
integer :: hdf5_err
integer(HID_T) :: pgroup
integer(HID_T) :: xs, temp_group
integer(SIZE_T) :: name_len
integer(HSIZE_T) :: dims(1)
integer(HSIZE_T) :: j
character(MAX_WORD_LEN) :: name
character(MAX_WORD_LEN) :: temp_str ! temperature dataset name, e.g. '294K'
character(MAX_WORD_LEN), allocatable :: grp_names(:)
call read_attribute(this % Q_value, group_id, 'Q_value')
call read_attribute(this % MT, group_id, 'mt')
@ -74,12 +69,10 @@ contains
end do
! Determine number of products
call h5gget_info_f(group_id, storage_type, n_links, max_corder, hdf5_err)
n_product = 0
do j = 0, n_links - 1
call h5lget_name_by_idx_f(group_id, ".", H5_INDEX_NAME_F, H5_ITER_INC_F, &
j, name, hdf5_err, name_len)
if (starts_with(name, "product_")) n_product = n_product + 1
call get_groups(group_id, grp_names)
do j = 1, size(grp_names)
if (starts_with(grp_names(j), "product_")) n_product = n_product + 1
end do
! Read products

16
src/sab_header.F90
View file

@ -8,9 +8,9 @@ module sab_header
use distribution_univariate, only: Tabular
use error, only: warning, fatal_error
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use h5lt, only: h5ltpath_valid_f, h5iget_name_f
use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_datasets
open_dataset, read_dataset, close_dataset, get_datasets, check_group, &
get_name
use secondary_correlated, only: CorrelatedAngleEnergy
use stl_vector, only: VectorInt, VectorReal
use string, only: to_str, str_to_int
@ -91,8 +91,7 @@ contains
integer :: n_energy, n_energy_out, n_mu
integer :: i_closest
integer :: n_temperature
integer :: hdf5_err
integer(SIZE_T) :: name_len, name_file_len
integer(SIZE_T) :: name_len
integer(HID_T) :: T_group
integer(HID_T) :: elastic_group
integer(HID_T) :: inelastic_group
@ -102,7 +101,6 @@ contains
integer(HSIZE_T) :: dims3(3)
real(8), allocatable :: temp(:,:)
character(20) :: type
logical :: exists
type(CorrelatedAngleEnergy) :: correlated_dist
character(MAX_WORD_LEN) :: temp_str
@ -114,7 +112,7 @@ contains
! Get name of table from group
name_len = len(this % name)
call h5iget_name_f(group_id, this % name, name_len, name_file_len, hdf5_err)
this % name = get_name(group_id, name_len)
! Get rid of leading '/'
this % name = trim(this % name(2:))
@ -180,8 +178,7 @@ contains
T_group = open_group(group_id, temp_str)
! Coherent elastic data
call h5ltpath_valid_f(T_group, 'elastic', .true., exists, hdf5_err)
if (exists) then
if (check_group(T_group, 'elastic')) then
! Read cross section data
elastic_group = open_group(T_group, 'elastic')
dset_id = open_dataset(elastic_group, 'xs')
@ -223,8 +220,7 @@ contains
end if
! Inelastic data
call h5ltpath_valid_f(T_group, 'inelastic', .true., exists, hdf5_err)
if (exists) then
if (check_group(T_group, 'inelastic')) then
! Read type of inelastic data
inelastic_group = open_group(T_group, 'inelastic')

17
src/secondary_uncorrelated.F90
View file

@ -1,6 +1,5 @@
module secondary_uncorrelated
use h5lt, only: h5ltpath_valid_f
use hdf5, only: HID_T
use angle_distribution, only: AngleDistribution
@ -9,7 +8,7 @@ module secondary_uncorrelated
use energy_distribution, only: EnergyDistribution, LevelInelastic, &
ContinuousTabular, MaxwellEnergy, Evaporation, WattEnergy, DiscretePhoton
use error, only: warning
use hdf5_interface, only: read_attribute, open_group, close_group
use hdf5_interface, only: read_attribute, open_group, close_group, check_group
use random_lcg, only: prn
!===============================================================================
@ -56,25 +55,19 @@ contains
class(UncorrelatedAngleEnergy), intent(inout) :: this
integer(HID_T), intent(in) :: group_id
logical :: exists
integer :: hdf5_err
integer(HID_T) :: energy_group
integer(HID_T) :: angle_group
character(MAX_WORD_LEN) :: type
! Check if energy group is present
call h5ltpath_valid_f(group_id, 'angle', .true., exists, hdf5_err)
if (exists) then
! Check if angle group is present & read
if (check_group(group_id, 'angle')) then
angle_group = open_group(group_id, 'angle')
call this%angle%from_hdf5(angle_group)
call close_group(angle_group)
end if
! Check if energy group is present
call h5ltpath_valid_f(group_id, 'energy', .true., exists, hdf5_err)
if (exists) then
! Check if energy group is present & read
if (check_group(group_id, 'energy')) then
energy_group = open_group(group_id, 'energy')
call read_attribute(type, energy_group, 'type')
select case (type)