crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-29 06:35:48 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Resolution of @paulromano comments

Browse source
This commit is contained in:
Adam Nelson 2016-10-10 15:19:47 -04:00
parent a52921e40a
commit 935d72af0d
14 changed files with 254 additions and 528 deletions
Download patch file Download diff file Expand all files Collapse all files

6
.gitignore vendored
View file

@ -50,8 +50,6 @@ tests/ctestscript.run
# HDF5 files
*.h5
!tests/1d_mgxs.h5
!examples/xml/pincell_multigroup/mgxs.h5
# Build files
src/CMakeCache.txt
@ -91,7 +89,3 @@ docs/source/pythonapi/examples/mgxs
docs/source/pythonapi/examples/tracks
docs/source/pythonapi/examples/fission-rates
docs/source/pythonapi/examples/plots
# VSCode workspace directory
.vscode
.vscode/*

52
docs/source/io_formats/mgxs_library.rst
View file

@ -38,8 +38,8 @@ MGXS Library Specification
The data within <library name> contains the temperature-dependent multi-group
data for the nuclide or material that it represents.
:Attributes: - **awr** (*double*) -- The atomic weight ratio (optional, i.e. it
is not meaningful for material-wise data)
:Attributes: - **atomic_weight_ratio** (*double*) -- The atomic weight ratio (optional,
i.e. it is not meaningful for material-wise data)
- **fissionable** (*int*) -- Whether the dataset is fissionable
(1) or not (0).
- **representation** (*char[]*) -- The method used to generate and
@ -49,26 +49,27 @@ data for the nuclide or material that it represents.
data was generated with angular dependent fluxes and thus the
data is angle-dependent. Valid values are either "isotropic" or
"angle".
- **num-azimuthal** (*int*) -- Number of equal width angular bins
- **num_azimuthal** (*int*) -- Number of equal width angular bins
that the azimuthal angular domain is subdivided if the
`representation` attribute is "angle". This parameter is
ignored otherwise.
- **num-polar** (*int*) -- Number of equal width angular bins
- **num_polar** (*int*) -- Number of equal width angular bins
that the polar angular domain is subdivided if the
`representation` attribute is "angle". This parameter is
ignored otherwise.
- **scatter-format** (*char[]*) -- The representation of the
- **scatter_format** (*char[]*) -- The representation of the
scattering angular distribution. The options are either
"legendre", "histogram", or "tabular". If not provided, the
default of "legendre" will be assumed.
- **order** (*int*) -- Either the Legendre order, number of bins,
or number of points (depending on the value of `scatter-format`)
or number of points (depending on the value of `scatter_format`)
used to describe the angular distribution associated with each
group-to-group transfer probability.
**/<library name>/kTs/**
:Datasets: - **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
:Datasets:
- **<TTT>K** (*double*) -- kT values (in MeV) for each Temperature
TTT (in Kelvin), rounded to the nearest integer
**/<library name>/<TTT>K/**
@ -109,44 +110,44 @@ Temperature-dependent data, provided for temperature <TTT>K.
the `nu-fission` data must represent the fission neutron energy
spectra as well and thus will have one additional dimension
for the outgoing energy group. In this case, `nu-fission` has the
same dimensionality as `multiplicity matrix`.
- **inverse velocities** (*double[]*) -- Average inverse velocity
same dimensionality as `multiplicity_matrix`.
- **inverse_velocities** (*double[]*) -- Average inverse velocity
for each of the groups in the library. This dataset is optional.
**/<library name>/<TTT>K/scatter data/**
**/<library name>/<TTT>K/scatter_data/**
Data specific to neutron scattering for the temperature <TTT>K
:Datasets: - **g_min** (*int[]* or *int[][][]) --
:Datasets: - **g_min** (*int[]* or *int[][][]*) --
Minimum (most energetic) outgoing groups with non-zero values of
the scattering matrix. These group numbers use the standard
ordering where the fastest neutron energy group is group 1 while
the slowest neutron energy group is group G.
The dimensionality of `g_out bounds` is:
`g_min[g_in]`, or `g_min[num-polar][num-azimuthal][g_in]`.
The dimensionality of `g_min` is:
`g_min[g_in]`, or `g_min[num_polar][num_azimuthal][g_in]`.
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle".
- **g_max** (*int[]* or *int[][][]) --
- **g_max** (*int[]* or *int[][][]*) --
Maximum (least energetic) outgoing groups with non-zero values of
the scattering matrix. These group numbers use the standard
ordering where the fastest neutron energy group is group 1 while
the slowest neutron energy group is group G.
The dimensionality of `g_out bounds` is:
`g_max[g_in]`, or `g_max[num-polar][num-azimuthal][g_in]`.
The dimensionality of `g_max` is:
`g_max[g_in]`, or `g_max[num_polar][num_azimuthal][g_in]`.
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle".
- **scatter matrix** (*double[]*) -- Flattened representation of the
- **scatter_matrix** (*double[]*) -- Flattened representation of the
scattering moment matrices. The pre-flattened array is shaped as
follows (in row-major format):
`scatter matrix[order(+1)][g_in][g_out]`, or
`scatter matrix[num-polar][num-azimuthal][order(+1)][g_in][g_out]`
`scatter_matrix[order(+1)][g_in][g_out]`, or
`scatter_matrix[num_polar][num_azimuthal][order(+1)][g_in][g_out]`
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle". Note that if the value of
`scatter-format` is "legendre", the order dimension will be one
`scatter_format` is "legendre", the order dimension will be one
larger than the value of `order`, otherwise it will match `order`.
Finally, the g_out dimension has a dimensionality of
`g_out bounds`[0] to `g_out bounds`[1].
- **multiplicity matrix** (*double[]*) -- Flattened representation of
`g_min` to `g_max`.
- **multiplicity_matrix** (*double[]*) -- Flattened representation of
the scattering moment matrices. This dataset provides the code with
a scaling factor to account for neutrons being produced in (n,xn)
reactions. This is assumed isotropic and therefore is not repeated
@ -154,9 +155,8 @@ Data specific to neutron scattering for the temperature <TTT>K
optional, if it is not provided no multiplication (i.e., values of
1.0) will be assumed.
The pre-flattened array is shaped as follows (in row-major format):
`multiplicity matrix[g_in][g_out]`, or
`multiplicity matrix[num-polar][num-azimuthal][g_in][g_out]`
`multiplicity_matrix[g_in][g_out]`, or
`multiplicity_matrix[num_polar][num_azimuthal][g_in][g_out]`
The former is used when `representation` is "isotropic", and the
latter when `representation` is "angle". Finally, the g_out
dimension has a dimensionality of `g_out bounds`[0] to
`g_out bounds`[1].
dimension has a dimensionality of `g_min` to `g_max`.

6
openmc/mesh.py
View file

@ -298,12 +298,12 @@ 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
# Values of +/-1e10 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., boundary_type='reflective')]
zplanes = [openmc.ZPlane(z0=-1e10., boundary_type='reflective'),
openmc.ZPlane(z0=1e10., boundary_type='reflective')]
else:
zplanes = [openmc.ZPlane(z0=self.lower_left[2],
boundary_type=bc[4]),

2
openmc/mgxs/library.py
View file

@ -912,7 +912,7 @@ class Library(object):
self.representation = 'isotropic'
if nuclide is not 'total':
xsdata.awr = self._nuclides[nuclide][1]
xsdata.atomic_weight_ratio = self._nuclides[nuclide][1]
if subdomain is None:
subdomain = 'all'

95
openmc/mgxs_library.py
View file

@ -31,7 +31,7 @@ class XSdata(object):
representation : {'isotropic', 'angle'}, optional
Method used in generating the MGXS (isotropic or angle-dependent flux
weighting). Defaults to 'isotropic'
temperatures : numpy.ndarray
temperatures : Iterable of float
Temperatures (in units of Kelvin) of the provided datasets. Defaults
to a single temperature at 294K.
@ -39,7 +39,7 @@ class XSdata(object):
----------
name : str
Unique identifier for the xsdata object
awr : float
aromic_weight_ratio : float
Atomic weight ratio of an isotope. That is, the ratio of the mass
of the isotope to the mass of a single neutron.
temperatures : numpy.ndarray
@ -60,24 +60,26 @@ class XSdata(object):
weighting).
num_azimuthal : int
Number of equal width angular bins that the azimuthal angular domain is
subdivided into. This only applies when ``representation`` is "angle".
subdivided into. This only applies when :attr:`XSdata.representation`
is "angle".
num_polar : int
Number of equal width angular bins that the polar angular domain is
subdivided into. This only applies when ``representation`` is "angle".
subdivided into. This only applies when :attr:`XSdata.representation`
is "angle".
use_chi : bool
Whether or not a chi vector or nu-fission matrix was used.
vector_shape : iterable of int
Dimensionality of vector multi-group cross sections (e.g., the total
cross section). The return result depends on the value of
``representation``.
:attr:`XSdata.representation`.
matrix_shape : iterable of int
Dimensionality of matrix multi-group cross sections (e.g., the
fission matrix cross section). The return result depends on the
value of ``representation``.
value of :attr:`XSdata.representation`.
pn_matrix_shape : iterable of int
Dimensionality of scattering matrix data (e.g., the
scattering matrix cross section). The return result depends on the
value of ``representation``.
value of :attr:`XSdata.representation`.
total : dict of numpy.ndarray
Group-wise total cross section.
absorption : dict of numpy.ndarray
@ -101,20 +103,22 @@ class XSdata(object):
approximation that the fission spectra does not depend on incoming
energy. If the user does not wish to make this approximation, then
this should not be provided and this information included in the
``nu_fission`` attribute instead.
:attr:`XSdata.nu_fission` attribute instead.
nu_fission : dict of numpy.ndarray
Group-wise fission production cross section vector (i.e., if ``chi`` is
provided), or is the group-wise fission production matrix.
Group-wise fission production cross section vector (i.e., if
:attr:`XSdata.chi` is provided), or is the group-wise fission production
matrix.
inverse_velocities : dict of numpy.ndarray
Inverse of velocities, in units of sec/cm.
Notes
-----
The parameters containing cross section data have dimensionalities which
depend upon the value of ``representation`` as well as the number of
Legendre or other angular dimensions as described by ``order``. The
``vector_shape``, ``matrix_shape``, and ``pn_matrix_shape`` properties are
provided to obtain the dimensionality of the data for each temperature.
depend upon the value of :attr:`XSdata.representation` as well as the
number of Legendre or other angular dimensions as described by
:attr:`XSdata.order`. The :attr:`XSdata.vector_shape`,
:attr:`XSdata.matrix_shape`, and :attr:`XSdata.pn_matrix_shape` properties
are provided to obtain the dimensionality of the data for each temperature.
The following are cross sections which should use each of these properties:
@ -135,7 +139,7 @@ class XSdata(object):
self.energy_groups = energy_groups
self.temperatures = temperatures
self.representation = representation
self._awr = None
self._atomic_weight_ratio = None
self._fissionable = False
self._scatter_format = 'legendre'
self._order = None
@ -165,8 +169,8 @@ class XSdata(object):
return self._representation
@property
def awr(self):
return self._awr
def atomic_weight_ratio(self):
return self._atomic_weight_ratio
@property
def fissionable(self):
@ -286,12 +290,12 @@ class XSdata(object):
check_value('representation', representation, _REPRESENTATIONS)
self._representation = representation
@awr.setter
def awr(self, awr):
# Check validity of type and that the awr value is > 0
check_type('awr', awr, Real)
check_greater_than('awr', awr, 0.0)
self._awr = awr
@atomic_weight_ratio.setter
def atomic_weight_ratio(self, atomic_weight_ratio):
# Check validity of type and that the atomic_weight_ratio value is > 0
check_type('atomic_weight_ratio', atomic_weight_ratio, Real)
check_greater_than('atomic_weight_ratio', atomic_weight_ratio, 0.0)
self._atomic_weight_ratio = atomic_weight_ratio
@fissionable.setter
def fissionable(self, fissionable):
@ -300,14 +304,9 @@ class XSdata(object):
@temperatures.setter
def temperatures(self, temperatures):
check_type('temperatures', temperatures, Iterable,
expected_iter_type=Real)
# Convert to a numpy array so we can easily get the shape for checking
nptemperatures = np.asarray(temperatures)
check_iterable_type('temperatures', temperatures, Real)
check_value('temperatures dimensionality', nptemperatures.ndim, [1])
self._temperatures = nptemperatures
self._temperatures = temperatures
@scatter_format.setter
def scatter_format(self, scatter_format):
@ -1113,10 +1112,18 @@ class XSdata(object):
self._multiplicity_matrix[i] = \
np.nan_to_num(self._multiplicity_matrix[i])
def _get_xsdata_group(self, file):
def to_hdf5(self, file):
"""Write XSdata to an HDF5 file
Parameters
----------
file : h5py.File
HDF5 File (a root Group) to write to
"""
grp = file.create_group(self.name)
if self.awr is not None:
grp.attrs['awr'] = self.awr
if self.atomic_weight_ratio is not None:
grp.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio
if self.fissionable is not None:
grp.attrs['fissionable'] = self.fissionable
if self.representation is not None:
@ -1124,11 +1131,11 @@ class XSdata(object):
dtype='S')
if self.representation == 'angle':
if self.num_azimuthal is not None:
grp.attrs['num-azimuthal'] = self.num_azimuthal
grp.attrs['num_azimuthal'] = self.num_azimuthal
if self.num_polar is not None:
grp.attrs['num-polar'] = self.num_polar
grp.attrs['num_polar'] = self.num_polar
if self.scatter_format is not None:
grp.attrs['scatter-format'] = np.array(self.scatter_format,
grp.attrs['scatter_format'] = np.array(self.scatter_format,
dtype='S')
if self.order is not None:
grp.attrs['order'] = self.order
@ -1184,8 +1191,8 @@ class XSdata(object):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
# And write it.
scatt_grp = xs_grp.create_group('scatter data')
scatt_grp.create_dataset("scatter matrix",
scatt_grp = xs_grp.create_group('scatter_data')
scatt_grp.create_dataset("scatter_matrix",
data=np.array(flat_scatt))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
@ -1196,7 +1203,7 @@ class XSdata(object):
for g_out in range(g_out_bounds[g_in, 0],
g_out_bounds[g_in, 1] + 1):
flat_mult.append(matrix[g_in, g_out])
scatt_grp.create_dataset("multiplicity matrix",
scatt_grp.create_dataset("multiplicity_matrix",
data=np.array(flat_mult))
# And finally, adjust g_out_bounds for 1-based group counting
@ -1227,8 +1234,8 @@ class XSdata(object):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
# And write it.
scatt_grp = xs_grp.create_group('scatter data')
scatt_grp.create_dataset("scatter matrix",
scatt_grp = xs_grp.create_group('scatter_data')
scatt_grp.create_dataset("scatter_matrix",
data=np.array(flat_scatt))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
@ -1242,7 +1249,7 @@ class XSdata(object):
g_out_bounds[p, a, g_in, 1] + 1):
flat_mult.append(matrix[g_in, g_out])
# And write it.
scatt_grp.create_dataset("multiplicity matrix",
scatt_grp.create_dataset("multiplicity_matrix",
data=np.array(flat_mult))
# And finally, adjust g_out_bounds for 1-based group counting
@ -1253,7 +1260,7 @@ class XSdata(object):
# Add the kinetics data
if self._inverse_velocities[i] is not None:
xs_grp.create_dataset("inverse-velocities",
xs_grp.create_dataset("inverse_velocities",
data=self._inverse_velocities[i])
@ -1364,6 +1371,6 @@ class MGXSLibrary(object):
file.attrs['group structure'] = self.energy_groups.group_edges
for xsdata in self._xsdatas:
xsdata._get_xsdata_group(file)
xsdata.to_hdf5(file)
file.close()

232
scripts/openmc-update-mgxs.py
View file

@ -1,232 +0,0 @@
#!/usr/bin/env python
"""Update OpenMC's deprecated multi-group cross section XML files to the latest
HDF5-based format.
Usage information can be obtained by running 'openmc-update-mgxs --help':
usage: openmc-update-mgxs [-h] in out
Update mgxs.xml files to the latest format. This will remove 'outside'
attributes/elements from lattices and replace them with 'outer' attributes. For
'cell' elements, any 'surfaces' attributes/elements will be renamed
'region'. Note that this script will not delete the given files; it will append
'.original' to the given files and write new ones.
positional arguments:
in Input mgxs xml file
out Output mgxs hdf5 file
optional arguments:
-h, --help show this help message and exit
"""
from __future__ import print_function
from shutil import move
import warnings
import xml.etree.ElementTree as ET
import argparse
import h5py
import numpy as np
import openmc.mgxs_library
description = """\
Update OpenMC's deprecated multi-group cross section XML files to the latest
HDF5-based format."""
def parse_args():
"""Read the input files from the commandline."""
# Create argument parser
parser = argparse.ArgumentParser(description=description,
formatter_class=argparse.RawTextHelpFormatter)
parser.add_argument('-i', '--input', type=argparse.FileType('r'),
help='input XML file')
parser.add_argument('-o', '--output', nargs='?', default='',
help='output file, in HDF5 format')
parser.add_argument('-c', '--compression', type=int,
help='HDF5 Compression Level')
args = vars(parser.parse_args())
if args['output'] == '':
filename = args['input'].name
extension = filename[filename.rfind('.'):]
if extension == '.xml':
filename = filename[:filename.rfind('.')] + '.h5'
args['output'] = filename
# Parse and return commandline arguments.
return args
def get_data(element, entry):
try:
value = element.find(entry).text
except:
if entry in element.attrib:
value = element.attrib[entry]
else:
value = None
if value is not None:
value = value.strip()
return value
if __name__ == '__main__':
args = parse_args()
# Parse the XML data.
tree = ET.parse(args['input'])
root = tree.getroot()
# Get old metadata
temp = tree.find('group_structure').text.strip()
temp = np.array(temp.split())
group_structure = temp.astype(np.float)
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
temp = tree.find('inverse_velocities')
if temp is not None:
temp = temp.text.strip()
temp = np.array(temp.split())
inverse_velocities = temp.astype(np.float)
else:
inverse_velocities = None
xsd = []
names = []
# Now move on to the cross section data itself
for xsdata_elem in root.iter('xsdata'):
name = get_data(xsdata_elem, 'name')
temperature = get_data(xsdata_elem, 'kT')
if temperature is not None:
temperature = \
float(temperature) / openmc.data.K_BOLTZMANN
else:
temperature = 294.
temperatures = [temperature]
awr = get_data(xsdata_elem, 'awr')
if awr is not None:
awr = float(awr)
representation = get_data(xsdata_elem, 'representation')
if representation is None:
representation = 'isotropic'
if representation == 'angle':
n_azi = int(get_data(xsdata_elem, 'num_azimuthal'))
n_pol = int(get_data(xsdata_elem, 'num_polar'))
scatter_format = get_data(xsdata_elem, 'scatt_type')
if scatter_format is None:
scatter_format = 'legendre'
order = int(get_data(xsdata_elem, 'order'))
tab_leg = get_data(xsdata_elem, 'tabular_legendre')
if tab_leg is not None:
warnings.Warning('The tabular_legendre option has moved to the '
'settings.xml file and must be added manually')
# Either add the data to a previously existing xsdata (if it is
# for the same 'name' but a different temperature), or create a
# new one.
try:
# It is in our list, so store that entry
i = names.index(name)
except:
# It is not in our list, so add it
i = -1
xsd.append(openmc.XSdata(name, energy_groups,
temperatures=temperatures,
representation=representation))
if awr is not None:
xsd[-1].awr = awr
if representation == 'angle':
xsd[-1].num_azimuthal = n_azi
xsd[-1].num_polar = n_pol
xsd[-1].scatter_format = scatter_format
xsd[-1].order = order
names.append(name)
if scatter_format == 'legendre':
order_dim = order + 1
else:
order_dim = order
if i != -1:
xsd[i].add_temperature(temperature)
temp = get_data(xsdata_elem, 'total')
if temp is not None:
temp = np.array(temp.split())
total = temp.astype(np.float)
total = np.reshape(total, xsd[i].vector_shape)
xsd[i].set_total(total, temperature)
if inverse_velocities is not None:
xsd[i].set_inverse_velocities(inverse_velocities, temperature)
temp = get_data(xsdata_elem, 'absorption')
temp = np.array(temp.split())
absorption = temp.astype(np.float)
absorption = np.reshape(absorption, xsd[i].vector_shape)
xsd[i].set_absorption(absorption, temperature)
temp = get_data(xsdata_elem, 'scatter')
temp = np.array(temp.split())
temp = temp.astype(np.float)
scatter = np.reshape(temp, xsd[i].pn_matrix_shape)
xsd[i].set_scatter_matrix(scatter, temperature)
temp = get_data(xsdata_elem, 'multiplicity')
if temp is not None:
temp = np.array(temp.split())
temp = temp.astype(np.float)
multiplicity = np.reshape(temp, xsd[i].matrix_shape)
xsd[i].set_multiplicity_matrix(multiplicity, temperature)
temp = get_data(xsdata_elem, 'fission')
if temp is not None:
temp = np.array(temp.split())
fission = temp.astype(np.float)
fission = np.reshape(fission, xsd[i].vector_shape)
xsd[i].set_fission(fission, temperature)
temp = get_data(xsdata_elem, 'kappa_fission')
if temp is not None:
temp = np.array(temp.split())
kappa_fission = temp.astype(np.float)
kappa_fission = np.reshape(kappa_fission, xsd[i].vector_shape)
xsd[i].set_kappa_fission(kappa_fission, temperature)
temp = get_data(xsdata_elem, 'chi')
if temp is not None:
temp = np.array(temp.split())
chi = temp.astype(np.float)
chi = np.reshape(chi, xsd[i].vector_shape)
xsd[i].set_chi(chi, temperature)
else:
chi = None
temp = get_data(xsdata_elem, 'nu_fission')
if temp is not None:
temp = np.array(temp.split())
temp = temp.astype(np.float)
if chi is not None:
nu_fission = np.reshape(temp, xsd[i].vector_shape)
else:
nu_fission = np.reshape(temp, xsd[i].matrix_shape)
xsd[i].set_nu_fission(nu_fission, temperature)
# Build library as we go, but first we have enough to initialize it
lib = openmc.MGXSLibrary(energy_groups)
lib.add_xsdatas(xsd)
lib.export_to_hdf5(args['output'])

78
src/geometry_header.F90
View file

@ -1,6 +1,12 @@
module geometry_header
use constants, only: HALF, TWO, THREE, INFINITY
use algorithm, only: find
use constants, only: HALF, TWO, THREE, INFINITY, K_BOLTZMANN, &
MATERIAL_VOID, NONE
use dict_header, only: DictCharInt, DictIntInt
use material_header, only: Material
use stl_vector, only: VectorReal
use string, only: to_lower
implicit none
@ -319,4 +325,74 @@ contains
end if
end function get_local_hex
!===============================================================================
! 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(cells, materials, material_dict, nuclide_dict, &
n_nucs, nuc_temps, sab_dict, n_sabs, sab_temps)
type(Cell), allocatable, intent(in) :: cells(:)
type(Material), allocatable, intent(in) :: materials(:)
type(DictIntInt), intent(in) :: material_dict
type(DictCharInt), intent(in) :: nuclide_dict
integer, intent(in) :: n_nucs
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(DictCharInt), optional, intent(in) :: sab_dict
integer, optional, intent(in) :: n_sabs
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_nucs))
if (present(n_sabs) .and. present(sab_temps)) allocate(sab_temps(n_sabs))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
else
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(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_key(to_lower(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
if (present(sab_temps) .and. present(sab_dict) .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_key(to_lower(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
end subroutine get_temperatures
end module geometry_header

33
src/hdf5_interface.F90
View file

@ -85,17 +85,16 @@ module hdf5_interface
public :: write_dataset
public :: read_dataset
public :: check_attribute
public :: attribute_exists
public :: write_attribute
public :: read_attribute
public :: file_create
public :: file_open
public :: file_close
public :: create_group
public :: check_group
public :: object_exists
public :: open_group
public :: close_group
public :: check_dataset
public :: open_dataset
public :: close_dataset
public :: get_shape
@ -253,7 +252,7 @@ contains
! CHECK_ATTRIBUTE Checks to see if an attribute exists in the object
!===============================================================================
function check_attribute(object_id, name) result(exists)
function attribute_exists(object_id, name) result(exists)
integer(HID_T), intent(in) :: object_id
character(*), intent(in) :: name ! name of group
logical :: exists
@ -263,13 +262,13 @@ contains
! Check if attribute exists
call h5aexists_by_name_f(object_id, '.', trim(name), exists, hdf5_err)
end function check_attribute
end function attribute_exists
!===============================================================================
! CHECK_GROUP Checks to see if a group exists in the object
!===============================================================================
function check_group(object_id, name) result(exists)
function object_exists(object_id, name) result(exists)
integer(HID_T), intent(in) :: object_id
character(*), intent(in) :: name ! name of group
logical :: exists
@ -279,21 +278,7 @@ contains
! Check if group exists
call h5ltpath_valid_f(object_id, trim(name), .true., exists, hdf5_err)
end function check_group
!===============================================================================
! CHECK_DATASET Checks to see if a dataset exists in the object
!===============================================================================
function check_dataset(object_id, name) result(exists)
integer(HID_T), intent(in) :: object_id
character(*), intent(in) :: name ! name of group
logical :: exists
! Wrap check_group since the method used there will work here too
exists = check_group(object_id, name)
end function check_dataset
end function object_exists
!===============================================================================
! GET_DATASETS Gets a list of all the datasets in a given location.
@ -367,7 +352,7 @@ contains
integer :: hdf5_err ! HDF5 error code
! Check if group exists
exists = check_group(group_id, name)
exists = object_exists(group_id, name)
! open group if it exists
if (exists) then
@ -390,7 +375,7 @@ contains
logical :: exists ! does the group exist
! Check if group exists
exists = check_group(group_id, name)
exists = object_exists(group_id, name)
! create group
if (exists) then
@ -428,7 +413,7 @@ contains
integer :: hdf5_err ! HDF5 error code
! Check if group exists
exists = check_group(group_id, name)
exists = object_exists(group_id, name)
! open group if it exists
if (exists) then

87
src/input_xml.F90
View file

@ -1,17 +1,15 @@
module input_xml
use hdf5
use algorithm, only: find
use cmfd_input, only: configure_cmfd
use constants
use dict_header, only: DictIntInt, ElemKeyValueCI
use dict_header, only: DictIntInt, DictCharInt, ElemKeyValueCI
use distribution_multivariate
use distribution_univariate
use endf, only: reaction_name
use energy_grid, only: grid_method, n_log_bins
use error, only: fatal_error, warning
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice
use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, &
get_temperatures
use global
use hdf5_interface
use list_header, only: ListChar, ListInt, ListReal
@ -61,6 +59,9 @@ contains
call time_read_xs % stop()
end if
! Normalize atom/weight percents
if (run_mode /= MODE_PLOTTING) call normalize_ao()
end subroutine read_input_xml
!===============================================================================
@ -2097,7 +2098,9 @@ contains
call assign_temperatures(material_temps)
! Determine desired temperatures for each nuclide and S(a,b) table
call get_temperatures(nuc_temps, sab_temps)
call get_temperatures(cells, materials, material_dict, nuclide_dict, &
n_nuclides_total, nuc_temps, sab_dict, &
n_sab_tables, sab_temps)
! Read continuous-energy cross sections
if (run_CE .and. run_mode /= MODE_PLOTTING) then
@ -2106,9 +2109,6 @@ contains
call time_read_xs % stop()
end if
! Normalize atom/weight percents
if (run_mode /= MODE_PLOTTING) call normalize_ao()
! Clear dictionary
call library_dict % clear()
end subroutine read_materials
@ -4631,7 +4631,7 @@ contains
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
if (check_attribute(file_id, "groups")) then
if (attribute_exists(file_id, "groups")) then
! Get neutron group count
call read_attribute(energy_groups, file_id, "groups")
else
@ -4640,7 +4640,7 @@ contains
allocate(rev_energy_bins(energy_groups + 1))
allocate(energy_bins(energy_groups + 1))
if (check_attribute(file_id, "group structure")) then
if (attribute_exists(file_id, "group structure")) then
! Get neutron group structure
call read_attribute(energy_bins, file_id, "group structure")
else
@ -5621,7 +5621,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = ONE
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
end if
! if given weight percent, convert all values so that they are divided
@ -5646,7 +5646,7 @@ contains
if (run_CE) then
awr = nuclides(mat % nuclide(j)) % awr
else
awr = ONE
awr = nuclides_MG(mat % nuclide(j)) % obj % awr
end if
x = mat % atom_density(j)
sum_percent = sum_percent + x*awr
@ -5910,67 +5910,6 @@ contains
end do
end subroutine assign_temperatures
!===============================================================================
! 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, sab_temps)
type(VectorReal), allocatable, intent(out) :: nuc_temps(:)
type(VectorReal), 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_total))
allocate(sab_temps(n_sab_tables))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
else
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(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_key(to_lower(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
if (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_key(to_lower(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
end subroutine get_temperatures
!===============================================================================
! READ_0K_ELASTIC_SCATTERING
!===============================================================================

56
src/mgxs_data.F90
View file

@ -3,8 +3,9 @@ module mgxs_data
use constants
use algorithm, only: find
use error, only: fatal_error
use geometry_header, only: get_temperatures
use global
use hdf5
use hdf5_interface
use material_header, only: Material
use mgxs_header
use output, only: write_message
@ -49,7 +50,8 @@ contains
call write_message("Loading Cross Section Data...", 5)
! Get temperatures
call get_temperatures(temps)
call get_temperatures(cells, materials, material_dict, nuclide_dict, &
n_nuclides_total, temps)
! Open file for reading
file_id = file_open(path_cross_sections, 'r', parallel=.true.)
@ -94,7 +96,7 @@ contains
call write_message("Loading " // trim(name) // " Data...", 5)
! Check to make sure cross section set exists in the library
if (check_group(file_id, trim(name))) then
if (object_exists(file_id, trim(name))) then
xsdata_group = open_group(file_id, trim(name))
else
call fatal_error("Data for '" // trim(name) // "' does not exist in "&
@ -102,7 +104,7 @@ contains
end if
! First find out the data representation
if (check_attribute(xsdata_group, "representation")) then
if (attribute_exists(xsdata_group, "representation")) then
call read_attribute(temp_str, xsdata_group, "representation")
if (trim(temp_str) == 'isotropic') then
representation = MGXS_ISOTROPIC
@ -236,51 +238,5 @@ contains
end subroutine get_mat_kTs
!===============================================================================
! GET_TEMPERATURES returns a list of temperatures that each MGXS 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(temps)
type(VectorReal), allocatable, intent(out) :: temps(:)
integer :: i, j, k
integer :: i_nuclide ! index in nuclides array
integer :: i_material
real(8) :: temperature ! temperature in Kelvin
allocate(temps(n_nuclides_total))
do i = 1, size(cells)
do j = 1, size(cells(i) % material)
! Skip any non-material cells and void materials
if (cells(i) % material(j) == NONE .or. &
cells(i) % material(j) == MATERIAL_VOID) cycle
! Get temperature of cell (rounding to nearest integer)
if (size(cells(i) % sqrtkT) > 1) then
temperature = cells(i) % sqrtkT(j)**2 / K_BOLTZMANN
else
temperature = cells(i) % sqrtkT(1)**2 / K_BOLTZMANN
end if
i_material = material_dict % get_key(cells(i) % material(j))
associate (mat => materials(i_material))
NUC_NAMES_LOOP: do k = 1, size(mat % names)
! Get index in temps array
i_nuclide = nuclide_dict % get_key(to_lower(mat % names(k)))
! Add temperature if it hasn't already been added
if (find(temps(i_nuclide), temperature) == -1) then
call temps(i_nuclide) % push_back(temperature)
end if
end do NUC_NAMES_LOOP
end associate
end do
end do
end subroutine get_temperatures
end module mgxs_data

114
src/mgxs_header.F90
View file

@ -229,8 +229,8 @@ module mgxs_header
! Get rid of leading '/'
this % name = trim(this % name(2:))
if (check_attribute(xs_id, "awr")) then
call read_attribute(this % awr, xs_id, "awr")
if (attribute_exists(xs_id, "atomic_weight_ratio")) then
call read_attribute(this % awr, xs_id, "atomic_weight_ratio")
else
this % awr = -ONE
end if
@ -323,8 +323,8 @@ module mgxs_header
end select
! Load the remaining metadata
if (check_attribute(xs_id, "scatter-format")) then
call read_attribute(temp_str, xs_id, "scatter-format")
if (attribute_exists(xs_id, "scatter_format")) then
call read_attribute(temp_str, xs_id, "scatter_format")
temp_str = trim(temp_str)
if (to_lower(temp_str) == 'legendre') then
this % scatter_format = ANGLE_LEGENDRE
@ -333,19 +333,19 @@ module mgxs_header
else if (to_lower(temp_str) == 'tabular') then
this % scatter_format = ANGLE_TABULAR
else
call fatal_error("Invalid scatter-format option!")
call fatal_error("Invalid scatter_format option!")
end if
else
this % scatter_format = ANGLE_LEGENDRE
end if
if (check_attribute(xs_id, "fissionable")) then
if (attribute_exists(xs_id, "fissionable")) then
call read_attribute(this % fissionable, xs_id, "fissionable")
else
call fatal_error("Fissionable element must be set!")
end if
! Get the library's value for the order
if (check_attribute(xs_id, "order")) then
if (attribute_exists(xs_id, "order")) then
call read_attribute(order_dim, xs_id, "order")
else
call fatal_error("Order must be provided!")
@ -366,16 +366,16 @@ module mgxs_header
! information therein
select type(this)
type is (MgxsAngle)
if (check_attribute(xs_id, "num-polar")) then
call read_attribute(this % n_pol, xs_id, "num-polar")
if (attribute_exists(xs_id, "num_polar")) then
call read_attribute(this % n_pol, xs_id, "num_polar")
else
call fatal_error("num-polar must be provided!")
call fatal_error("num_polar must be provided!")
end if
if (check_attribute(xs_id, "num-azimuthal")) then
call read_attribute(this % n_azi, xs_id, "num-azimuthal")
if (attribute_exists(xs_id, "num_azimuthal")) then
call read_attribute(this % n_azi, xs_id, "num_azimuthal")
else
call fatal_error("num-azimuthal must be provided!")
call fatal_error("num_azimuthal must be provided!")
end if
! Set angle data to use equally-spaced bins
@ -433,7 +433,7 @@ module mgxs_header
if (this % fissionable) then
allocate(xs % nu_fission(groups))
allocate(xs % chi(groups, groups))
if (check_dataset(xsdata_grp, "chi")) then
if (object_exists(xsdata_grp, "chi")) then
! Chi was provided, that means we need chi and nu-fission vectors
! Get chi
allocate(temp_arr(groups))
@ -448,7 +448,7 @@ module mgxs_header
deallocate(temp_arr)
! Get nu_fission (as a vector)
if (check_dataset(xsdata_grp, "nu-fission")) then
if (object_exists(xsdata_grp, "nu-fission")) then
call read_dataset(xs % nu_fission, xsdata_grp, "nu-fission")
else
call fatal_error("If fissionable, must provide nu-fission!")
@ -457,7 +457,7 @@ module mgxs_header
else
! chi isnt provided but is within nu_fission, existing as a matrix
! So, get nu_fission (as a matrix)
if (check_dataset(xsdata_grp, "nu-fission")) then
if (object_exists(xsdata_grp, "nu-fission")) then
allocate(temp_2d(groups, groups))
call read_dataset(temp_2d, xsdata_grp, "nu-fission")
else
@ -481,7 +481,7 @@ 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
if (check_dataset(xsdata_grp, "fission")) then
if (object_exists(xsdata_grp, "fission")) then
allocate(xs % fission(groups))
call read_dataset(xs % fission, xsdata_grp, "fission")
else
@ -490,7 +490,7 @@ module mgxs_header
end if
end if
if (get_kfiss) then
if (check_dataset(xsdata_grp, "kappa-fission")) then
if (object_exists(xsdata_grp, "kappa-fission")) then
allocate(xs % k_fission(groups))
call read_dataset(xs % k_fission, xsdata_grp, "kappa-fission")
else
@ -500,7 +500,7 @@ module mgxs_header
end if
end if
if (check_dataset(xsdata_grp, "absorption")) then
if (object_exists(xsdata_grp, "absorption")) then
allocate(xs % absorption(groups))
call read_dataset(xs % absorption, xsdata_grp, "absorption")
else
@ -508,21 +508,21 @@ module mgxs_header
end if
! Get scattering data
if (.not. check_group(xsdata_grp, "scatter data")) &
call fatal_error("Must provide 'scatter data'")
scatt_grp = open_group(xsdata_grp, 'scatter data')
if (.not. object_exists(xsdata_grp, "scatter_data")) &
call fatal_error("Must provide 'scatter_data'")
scatt_grp = open_group(xsdata_grp, 'scatter_data')
! First get the outgoing group boundary indices
if (check_dataset(scatt_grp, "g_min")) then
if (object_exists(scatt_grp, "g_min")) then
allocate(gmin(groups))
call read_dataset(gmin, scatt_grp, "g_min")
else
call fatal_error("'g_min' for the scatter matrix must be provided")
call fatal_error("'g_min' for the scatter_data must be provided")
end if
if (check_dataset(scatt_grp, "g_max")) then
if (object_exists(scatt_grp, "g_max")) then
allocate(gmax(groups))
call read_dataset(gmax, scatt_grp, "g_max")
else
call fatal_error("'g_max' for the scatter matrix must be provided")
call fatal_error("'g_max' for the scatter_data must be provided")
end if
! Now use this information to find the length of a container array
@ -533,9 +533,9 @@ module mgxs_header
end do
! Allocate flattened array
allocate(temp_arr(length))
if (.not. check_dataset(scatt_grp, 'scatter matrix')) &
call fatal_error("'scatter matrix' must be provided")
call read_dataset(temp_arr, scatt_grp, "scatter matrix")
if (.not. object_exists(scatt_grp, 'scatter_matrix')) &
call fatal_error("'scatter_matrix' must be provided")
call read_dataset(temp_arr, scatt_grp, "scatter_matrix")
! Compare the number of orders given with the maximum order of the
! problem. Strip off the supefluous orders if needed.
@ -611,7 +611,7 @@ module mgxs_header
deallocate(input_scatt)
! Now get the multiplication matrix
if (check_dataset(scatt_grp, 'multiplicity matrix')) then
if (object_exists(scatt_grp, 'multiplicity_matrix')) then
! Now use this information to find the length of a container array
! to hold the flattened data
length = 0
@ -620,7 +620,7 @@ module mgxs_header
end do
! Allocate flattened array
allocate(temp_arr(length))
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
call read_dataset(temp_arr, scatt_grp, "multiplicity_matrix")
! Convert temp_arr to a jagged array ((gin) % data(gout)) for passing
! to ScattData
@ -664,7 +664,7 @@ module mgxs_header
! Get, or infer, total xs data.
allocate(xs % total(groups))
if (check_dataset(xsdata_grp, "total")) then
if (object_exists(xsdata_grp, "total")) then
call read_dataset(xs % total, xsdata_grp, "total")
else
xs % total(:) = xs % absorption(:) + xs % scatter % scattxs(:)
@ -677,9 +677,9 @@ module mgxs_header
end do
! Get kinetics data
if (check_dataset(xsdata_grp, "inverse-velocities")) then
if (object_exists(xsdata_grp, "inverse_velocities")) then
allocate(xs % inv_vel(groups))
call read_dataset(xs % inv_vel, xsdata_grp, "inverse-velocities")
call read_dataset(xs % inv_vel, xsdata_grp, "inverse_velocities")
end if
! Close the groups we have opened and deallocate
@ -731,7 +731,7 @@ module mgxs_header
if (this % fissionable) then
allocate(xs % nu_fission(groups, this % n_azi, this % n_pol))
allocate(xs % chi(groups, groups, this % n_azi, this % n_pol))
if (check_dataset(xsdata_grp, "chi")) then
if (object_exists(xsdata_grp, "chi")) then
! Chi was provided, that means we need chi and nu-fission vectors
! Get chi
allocate(temp_arr(groups * this % n_azi * this % n_pol))
@ -764,7 +764,7 @@ module mgxs_header
deallocate(temp_arr)
! Get nu_fission (as a vector)
if (check_dataset(xsdata_grp, "nu-fission")) then
if (object_exists(xsdata_grp, "nu-fission")) then
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "nu-fission")
xs % nu_fission = reshape(temp_arr, (/groups, this % n_azi, &
@ -777,7 +777,7 @@ module mgxs_header
else
! chi isnt provided but is within nu_fission, existing as a matrix
! So, get nu_fission (as a matrix)
if (check_dataset(xsdata_grp, "nu-fission")) then
if (object_exists(xsdata_grp, "nu-fission")) then
allocate(temp_arr(groups * groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "nu-fission")
allocate(temp_4d(groups, groups, this % n_azi, this % n_pol))
@ -816,7 +816,7 @@ 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
if (check_dataset(xsdata_grp, "fission")) then
if (object_exists(xsdata_grp, "fission")) then
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "fission")
allocate(xs % fission(groups, this % n_azi, this % n_pol))
@ -829,7 +829,7 @@ module mgxs_header
end if
end if
if (get_kfiss) then
if (check_dataset(xsdata_grp, "kappa-fission")) then
if (object_exists(xsdata_grp, "kappa-fission")) then
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "kappa-fission")
allocate(xs % k_fission(groups, this % n_azi, this % n_pol))
@ -843,7 +843,7 @@ module mgxs_header
end if
end if
if (check_dataset(xsdata_grp, "absorption")) then
if (object_exists(xsdata_grp, "absorption")) then
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "absorption")
allocate(xs % absorption(groups, this % n_azi, this % n_pol))
@ -855,27 +855,27 @@ module mgxs_header
end if
! Get scattering data
if (.not. check_group(xsdata_grp, "scatter data")) &
call fatal_error("Must provide 'scatter data'")
scatt_grp = open_group(xsdata_grp, 'scatter data')
if (.not. object_exists(xsdata_grp, "scatter_data")) &
call fatal_error("Must provide 'scatter_data'")
scatt_grp = open_group(xsdata_grp, 'scatter_data')
! First get the outgoing group boundary indices
if (check_dataset(scatt_grp, "g_min")) then
if (object_exists(scatt_grp, "g_min")) then
allocate(int_arr(groups * this % n_azi * this % n_pol))
call read_dataset(int_arr, scatt_grp, "g_min")
allocate(gmin(groups, this % n_azi, this % n_pol))
gmin = reshape(int_arr, (/groups, this % n_azi, this % n_pol/))
deallocate(int_arr)
else
call fatal_error("'g_min' for the scatter matrix must be provided")
call fatal_error("'g_min' for the scatter_data must be provided")
end if
if (check_dataset(scatt_grp, "g_max")) then
if (object_exists(scatt_grp, "g_max")) then
allocate(int_arr(groups * this % n_azi * this % n_pol))
call read_dataset(int_arr, scatt_grp, "g_max")
allocate(gmax(groups, this % n_azi, this % n_pol))
gmax = reshape(int_arr, (/groups, this % n_azi, this % n_pol/))
deallocate(int_arr)
else
call fatal_error("'g_max' for the scatter matrix must be provided")
call fatal_error("'g_max' for the scatter_data must be provided")
end if
! Now use this information to find the length of a container array
@ -891,9 +891,9 @@ module mgxs_header
end do
! Allocate flattened array
allocate(temp_arr(length))
if (.not. check_dataset(scatt_grp, 'scatter matrix')) &
call fatal_error("'scatter matrix' must be provided")
call read_dataset(temp_arr, scatt_grp, "scatter matrix")
if (.not. object_exists(scatt_grp, 'scatter_matrix')) &
call fatal_error("'scatter_matrix' must be provided")
call read_dataset(temp_arr, scatt_grp, "scatter_matrix")
! Compare the number of orders given with the maximum order of the
! problem. Strip off the superfluous orders if needed.
@ -988,7 +988,7 @@ module mgxs_header
deallocate(input_scatt)
! Now get the multiplication matrix
if (check_dataset(scatt_grp, 'multiplicity matrix')) then
if (object_exists(scatt_grp, 'multiplicity_matrix')) then
! Now use this information to find the length of a container array
! to hold the flattened data
length = 0
@ -1001,7 +1001,7 @@ module mgxs_header
end do
! Allocate flattened array
allocate(temp_arr(length))
call read_dataset(temp_arr, scatt_grp, "multiplicity matrix")
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))
@ -1067,7 +1067,7 @@ module mgxs_header
end do
allocate(xs % total(groups, this % n_azi, this % n_pol))
if (check_dataset(xsdata_grp, "total")) then
if (object_exists(xsdata_grp, "total")) then
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "total")
xs % total = reshape(temp_arr, (/groups, this % n_azi, &
@ -1095,10 +1095,10 @@ module mgxs_header
end do
! Get kinetics data
if (check_dataset(xsdata_grp, "inverse-velocities")) then
if (object_exists(xsdata_grp, "inverse_velocities")) then
allocate(xs % inv_vel(groups, this % n_azi, this % n_pol))
allocate(temp_arr(groups * this % n_azi * this % n_pol))
call read_dataset(temp_arr, xsdata_grp, "inverse-velocities")
call read_dataset(temp_arr, xsdata_grp, "inverse_velocities")
xs % inv_vel = reshape(temp_arr, (/groups, this % n_azi, &
this % n_pol/))
deallocate(temp_arr)
@ -1372,7 +1372,7 @@ module mgxs_header
! dense matrix for this storage, with a reduction to the sparse
! format at the end.
! Get the multiplicity matrix
! Get the multiplicity_matrix
! To combine from nuclidic data we need to use the final relationship
! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))
@ -1639,7 +1639,7 @@ module mgxs_header
! dense matrix for this storage, with a reduction to the sparse
! format at the end.
! Get the multiplicity matrix
! Get the multiplicity_matrix
! To combine from nuclidic data we need to use the final relationship
! mult_{gg'} = sum_i(N_i*nuscatt_{i,g,g'}) /
! sum_i(N_i*(nuscatt_{i,g,g'} / mult_{i,g,g'}))

8
src/nuclide_header.F90
View file

@ -13,7 +13,7 @@ module nuclide_header
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, &
check_group, get_name, get_groups
object_exists, get_name, get_groups
use list_header, only: ListInt
use math, only: evaluate_legendre
use multipole_header, only: MultipoleArray
@ -350,7 +350,7 @@ module nuclide_header
call close_group(rxs_group)
! Read unresolved resonance probability tables if present
if (check_group(group_id, 'urr')) then
if (object_exists(group_id, 'urr')) then
this % urr_present = .true.
allocate(this % urr_data(n_temperature))
@ -391,7 +391,7 @@ module nuclide_header
end if
! Check for nu-total
if (check_group(group_id, 'total_nu')) then
if (object_exists(group_id, 'total_nu')) then
nu_group = open_group(group_id, 'total_nu')
! Read total nu data
@ -410,7 +410,7 @@ module nuclide_header
end if
! Read fission energy release data if present
if (check_group(group_id, 'fission_energy_release')) then
if (object_exists(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

6
src/sab_header.F90
View file

@ -9,7 +9,7 @@ module sab_header
use error, only: warning, fatal_error
use hdf5, only: HID_T, HSIZE_T, SIZE_T
use hdf5_interface, only: read_attribute, get_shape, open_group, close_group, &
open_dataset, read_dataset, close_dataset, get_datasets, check_group, &
open_dataset, read_dataset, close_dataset, get_datasets, object_exists, &
get_name
use secondary_correlated, only: CorrelatedAngleEnergy
use stl_vector, only: VectorInt, VectorReal
@ -209,7 +209,7 @@ contains
T_group = open_group(group_id, temp_str)
! Coherent elastic data
if (check_group(T_group, 'elastic')) then
if (object_exists(T_group, 'elastic')) then
! Read cross section data
elastic_group = open_group(T_group, 'elastic')
dset_id = open_dataset(elastic_group, 'xs')
@ -251,7 +251,7 @@ contains
end if
! Inelastic data
if (check_group(T_group, 'inelastic')) then
if (object_exists(T_group, 'inelastic')) then
! Read type of inelastic data
inelastic_group = open_group(T_group, 'inelastic')

7
src/secondary_uncorrelated.F90
View file

@ -8,7 +8,8 @@ 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, check_group
use hdf5_interface, only: read_attribute, open_group, close_group, &
object_exists
use random_lcg, only: prn
!===============================================================================
@ -60,14 +61,14 @@ contains
character(MAX_WORD_LEN) :: type
! Check if angle group is present & read
if (check_group(group_id, 'angle')) then
if (object_exists(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 & read
if (check_group(group_id, 'energy')) then
if (object_exists(group_id, 'energy')) then
energy_group = open_group(group_id, 'energy')
call read_attribute(type, energy_group, 'type')
select case (type)