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

Merge branch 'angles' into develop

Browse source
This commit is contained in:
Paul Romano 2016-11-28 20:41:22 -06:00
commit 879d41d701
16 changed files with 815 additions and 451 deletions
Download patch file Download diff file Expand all files Collapse all files

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

File diff suppressed because one or more lines are too long

10
examples/python/pincell_multigroup/build-xml.py
View file

@ -1,3 +1,5 @@
import numpy as np
import openmc
import openmc.mgxs
@ -26,7 +28,7 @@ uo2_xsdata.set_total(
0.5644058])
uo2_xsdata.set_absorption([8.0248E-03, 3.7174E-03, 2.6769E-02, 9.6236E-02,
3.0020E-02, 1.1126E-01, 2.8278E-01])
uo2_xsdata.set_scatter_matrix(
scatter_matrix = np.array(
[[[0.1275370, 0.0423780, 0.0000094, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.3244560, 0.0016314, 0.0000000, 0.0000000, 0.0000000, 0.0000000],
[0.0000000, 0.0000000, 0.4509400, 0.0026792, 0.0000000, 0.0000000, 0.0000000],
@ -34,6 +36,8 @@ uo2_xsdata.set_scatter_matrix(
[0.0000000, 0.0000000, 0.0000000, 0.0001253, 0.2714010, 0.0102550, 0.0000000],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0012968, 0.2658020, 0.0168090],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0085458, 0.2730800]]])
scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
uo2_xsdata.set_scatter_matrix(scatter_matrix)
uo2_xsdata.set_fission([7.21206E-03, 8.19301E-04, 6.45320E-03,
1.85648E-02, 1.78084E-02, 8.30348E-02,
2.16004E-01])
@ -50,7 +54,7 @@ h2o_xsdata.set_total([0.15920605, 0.412969593, 0.59030986, 0.58435,
h2o_xsdata.set_absorption([6.0105E-04, 1.5793E-05, 3.3716E-04,
1.9406E-03, 5.7416E-03, 1.5001E-02,
3.7239E-02])
h2o_xsdata.set_scatter_matrix(
scatter_matrix = np.array(
[[[0.0444777, 0.1134000, 0.0007235, 0.0000037, 0.0000001, 0.0000000, 0.0000000],
[0.0000000, 0.2823340, 0.1299400, 0.0006234, 0.0000480, 0.0000074, 0.0000010],
[0.0000000, 0.0000000, 0.3452560, 0.2245700, 0.0169990, 0.0026443, 0.0005034],
@ -58,6 +62,8 @@ h2o_xsdata.set_scatter_matrix(
[0.0000000, 0.0000000, 0.0000000, 0.0000714, 0.1391380, 0.5118200, 0.0612290],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0022157, 0.6999130, 0.5373200],
[0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.0000000, 0.1324400, 2.4807000]]])
scatter_matrix = np.rollaxis(scatter_matrix, 0, 3)
h2o_xsdata.set_scatter_matrix(scatter_matrix)
mg_cross_sections_file = openmc.MGXSLibrary(groups)
mg_cross_sections_file.add_xsdatas([uo2_xsdata, h2o_xsdata])

BIN
examples/xml/pincell_multigroup/mgxs.h5
View file

Binary file not shown.

325
openmc/filter.py
View file

@ -99,9 +99,13 @@ class Filter(object):
@classmethod
def _recursive_subclasses(cls):
"""Return all subclasses and their subclasses, etc."""
subs = cls.__subclasses__()
subsubs = [grand for s in subs for grand in s.__subclasses__()]
return subs + subsubs
all_subclasses = []
for subclass in cls.__subclasses__():
all_subclasses.append(subclass)
all_subclasses.extend(subclass._recursive_subclasses())
return all_subclasses
@classmethod
def from_hdf5(cls, group, **kwargs):
@ -773,7 +777,114 @@ class MeshFilter(Filter):
return df
class EnergyFilter(Filter):
class RealFilter(Filter):
"""Tally modifier that describes phase-space and other characteristics
Parameters
----------
bins : Iterable of Real
A grid of bin values.
Attributes
----------
bins : Iterable of Real
A grid of bin values.
num_bins : Integral
The number of filter bins
stride : Integral
The number of filter, nuclide and score bins within each of this
filter's bins.
"""
def __gt__(self, other):
if type(self) is type(other):
# Compare largest/smallest bin edges in filters
# This logic is used when merging tallies with real filters
return self.bins[0] >= other.bins[-1]
else:
return super(RealFilter, self).__gt__(other)
@property
def num_bins(self):
return len(self.bins) - 1
@num_bins.setter
def num_bins(self, num_bins):
cv.check_type('filter num_bins', num_bins, Integral)
cv.check_greater_than('filter num_bins', num_bins, 0, equality=True)
self._num_bins = num_bins
def can_merge(self, other):
if type(self) is not type(other):
return False
if self.bins[0] == other.bins[-1]:
# This low edge coincides with other's high edge
return True
elif self.bins[-1] == other.bins[0]:
# This high edge coincides with other's low edge
return True
else:
return False
def merge(self, other):
if not self.can_merge(other):
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(type(self), type(other))
raise ValueError(msg)
# Merge unique filter bins
merged_bins = np.concatenate((self.bins, other.bins))
merged_bins = np.unique(merged_bins)
# Create a new filter with these bins
return type(self)(sorted(merged_bins))
def is_subset(self, other):
"""Determine if another filter is a subset of this filter.
If all of the bins in the other filter are included as bins in this
filter, then it is a subset of this filter.
Parameters
----------
other : openmc.Filter
The filter to query as a subset of this filter
Returns
-------
bool
Whether or not the other filter is a subset of this filter
"""
if type(self) is not type(other):
return False
elif len(self.bins) != len(other.bins):
return False
else:
return np.allclose(self.bins, other.bins)
def get_bin_index(self, filter_bin):
i = np.where(self.bins == filter_bin[1])[0]
if len(i) == 0:
msg = 'Unable to get the bin index for Filter since "{0}" ' \
'is not one of the bins'.format(filter_bin)
raise ValueError(msg)
else:
return i[0] - 1
def get_bin(self, bin_index):
cv.check_type('bin_index', bin_index, Integral)
cv.check_greater_than('bin_index', bin_index, 0, equality=True)
cv.check_less_than('bin_index', bin_index, self.num_bins)
# Construct 2-tuple of lower, upper bins for real-valued filters
return (self.bins[bin_index], self.bins[bin_index + 1])
class EnergyFilter(RealFilter):
"""Bins tally events based on incident particle energy.
Parameters
@ -793,23 +904,16 @@ class EnergyFilter(Filter):
"""
def __gt__(self, other):
if type(self) is type(other):
# Compare largest/smallest energy bin edges in energy filters
# This logic is used when merging tallies with energy filters
return self.bins[0] >= other.bins[-1]
def get_bin_index(self, filter_bin):
# Use lower energy bound to find index for RealFilters
deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
min_delta = np.min(deltas)
if min_delta < 1E-3:
return deltas.argmin() - 1
else:
return super(EnergyFilter, self).__gt__(other)
@property
def num_bins(self):
return len(self.bins) - 1
@num_bins.setter
def num_bins(self, num_bins):
cv.check_type('filter num_bins', num_bins, Integral)
cv.check_greater_than('filter num_bins', num_bins, 0, equality=True)
self._num_bins = num_bins
msg = 'Unable to get the bin index for Filter since "{0}" ' \
'is not one of the bins'.format(filter_bin)
raise ValueError(msg)
def check_bins(self, bins):
for edge in bins:
@ -831,59 +935,6 @@ class EnergyFilter(Filter):
'increasing'.format(bins, type(self))
raise ValueError(msg)
def can_merge(self, other):
if type(self) is not type(other):
return False
if self.bins[0] == other.bins[-1]:
# This low energy edge coincides with other's high energy edge
return True
elif self.bins[-1] == other.bins[0]:
# This high energy edge coincides with other's low energy edge
return True
else:
return False
def merge(self, other):
if not self.can_merge(other):
msg = 'Unable to merge "{0}" with "{1}" ' \
'filters'.format(type(self), other.type)
raise ValueError(msg)
# Merge unique filter bins
merged_bins = np.concatenate((self.bins, other.bins))
merged_bins = np.unique(merged_bins)
# Create a new filter with these bins
return type(self)(sorted(merged_bins))
def is_subset(self, other):
if type(self) is not type(other):
return False
elif len(self.bins) != len(other.bins):
return False
else:
return np.allclose(self.bins, other.bins)
def get_bin_index(self, filter_bin):
# Use lower energy bound to find index for energy Filters
deltas = np.abs(self.bins - filter_bin[1]) / filter_bin[1]
min_delta = np.min(deltas)
if min_delta < 1E-3:
return deltas.argmin() - 1
else:
msg = 'Unable to get the bin index for Filter since "{0}" ' \
'is not one of the bins'.format(filter_bin)
raise ValueError(msg)
def get_bin(self, bin_index):
cv.check_type('bin_index', bin_index, Integral)
cv.check_greater_than('bin_index', bin_index, 0, equality=True)
cv.check_less_than('bin_index', bin_index, self.num_bins)
# Construct 2-tuple of lower, upper energies for energy(out) filters
return (self.bins[bin_index], self.bins[bin_index+1])
def get_pandas_dataframe(self, data_size, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
@ -1227,7 +1278,7 @@ class DistribcellFilter(Filter):
return df
class MuFilter(Filter):
class MuFilter(RealFilter):
"""Bins tally events based on particle scattering angle.
Parameters
@ -1255,8 +1306,82 @@ class MuFilter(Filter):
"""
def check_bins(self, bins):
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, type(self))
raise ValueError(msg)
elif edge < -1.:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is less than -1'.format(edge, type(self))
raise ValueError(msg)
elif edge > 1.:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is greater than 1'.format(edge, type(self))
raise ValueError(msg)
class PolarFilter(Filter):
# Check that bin edges are monotonically increasing
for index in range(1, len(bins)):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
This method constructs a Pandas DataFrame object for the filter with
columns annotated by filter bin information. This is a helper method
for :meth:`Tally.get_pandas_dataframe`.
Parameters
----------
data_size : Integral
The total number of bins in the tally corresponding to this filter
Returns
-------
pandas.DataFrame
A Pandas DataFrame with one column of the lower energy bound and one
column of upper energy bound for each filter bin. The number of
rows in the DataFrame is the same as the total number of bins in the
corresponding tally, with the filter bin appropriately tiled to map
to the corresponding tally bins.
Raises
------
ImportError
When Pandas is not installed
See also
--------
Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe()
"""
# Initialize Pandas DataFrame
import pandas as pd
df = pd.DataFrame()
# Extract the lower and upper energy bounds, then repeat and tile
# them as necessary to account for other filters.
lo_bins = np.repeat(self.bins[:-1], self.stride)
hi_bins = np.repeat(self.bins[1:], self.stride)
tile_factor = data_size / len(lo_bins)
lo_bins = np.tile(lo_bins, tile_factor)
hi_bins = np.tile(hi_bins, tile_factor)
# Add the new energy columns to the DataFrame.
df.loc[:, self.short_name.lower() + ' low'] = lo_bins
df.loc[:, self.short_name.lower() + ' high'] = hi_bins
return df
class PolarFilter(RealFilter):
"""Bins tally events based on the incident particle's direction.
Parameters
@ -1284,6 +1409,30 @@ class PolarFilter(Filter):
"""
def check_bins(self, bins):
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, type(self))
raise ValueError(msg)
elif edge < 0.:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is less than 0'.format(edge, type(self))
raise ValueError(msg)
elif edge > np.pi:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is greater than pi'.format(edge, type(self))
raise ValueError(msg)
# Check that bin edges are monotonically increasing
for index in range(1, len(bins)):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, **kwargs):
"""Builds a Pandas DataFrame for the Filter's bins.
@ -1335,7 +1484,7 @@ class PolarFilter(Filter):
return df
class AzimuthalFilter(Filter):
class AzimuthalFilter(RealFilter):
"""Bins tally events based on the incident particle's direction.
Parameters
@ -1363,6 +1512,30 @@ class AzimuthalFilter(Filter):
"""
def check_bins(self, bins):
for edge in bins:
if not isinstance(edge, Real):
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is a non-integer or floating point ' \
'value'.format(edge, type(self))
raise ValueError(msg)
elif edge < -np.pi:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is less than -pi'.format(edge, type(self))
raise ValueError(msg)
elif edge > np.pi:
msg = 'Unable to add bin edge "{0}" to a "{1}" ' \
'since it is greater than pi'.format(edge, type(self))
raise ValueError(msg)
# Check that bin edges are monotonically increasing
for index in range(1, len(bins)):
if bins[index] < bins[index-1]:
msg = 'Unable to add bin edges "{0}" to a "{1}" Filter ' \
'since they are not monotonically ' \
'increasing'.format(bins, type(self))
raise ValueError(msg)
def get_pandas_dataframe(self, data_size, distribcell_paths=True):
"""Builds a Pandas DataFrame for the Filter's bins.

138
openmc/mgxs/library.py
View file

@ -59,15 +59,23 @@ class Library(object):
The spatial domain(s) for which MGXS in the Library are computed
correction : {'P0', None}
Apply the P0 correction to scattering matrices if set to 'P0'
scatter_format : {'legendre', 'histogram'}
Representation of the angular scattering distribution (default is
'legendre')
legendre_order : int
The highest legendre moment in the scattering matrices (default is 0)
The highest Legendre moment in the scattering matrix; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0)
histogram_bins : int
The number of equally-spaced bins for the histogram representation of
the angular scattering distribution; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16)
energy_groups : openmc.mgxs.EnergyGroups
Energy group structure for energy condensation
num_delayed_groups : int
Number of delayed groups
estimator : str or None
The tally estimator used to compute multi-group cross sections. If None,
the default for each MGXS type is used.
The tally estimator used to compute multi-group cross sections.
If None, the default for each MGXS type is used.
tally_trigger : openmc.Trigger
An (optional) tally precision trigger given to each tally used to
compute the cross section
@ -101,7 +109,9 @@ class Library(object):
self._energy_groups = None
self._num_delayed_groups = 0
self._correction = 'P0'
self._scatter_format = 'legendre'
self._legendre_order = 0
self._histogram_bins = 16
self._tally_trigger = None
self._all_mgxs = OrderedDict()
self._sp_filename = None
@ -130,7 +140,9 @@ class Library(object):
clone._domain_type = self.domain_type
clone._domains = copy.deepcopy(self.domains)
clone._correction = self.correction
clone._scatter_format = self.scatter_format
clone._legendre_order = self.legendre_order
clone._histogram_bins = self.histogram_bins
clone._energy_groups = copy.deepcopy(self.energy_groups, memo)
clone._num_delayed_groups = self.num_delayed_groups
clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo)
@ -209,10 +221,18 @@ class Library(object):
def correction(self):
return self._correction
@property
def scatter_format(self):
return self._scatter_format
@property
def legendre_order(self):
return self._legendre_order
@property
def histogram_bins(self):
return self._histogram_bins
@property
def tally_trigger(self):
return self._tally_trigger
@ -267,7 +287,7 @@ class Library(object):
def by_nuclide(self, by_nuclide):
cv.check_type('by_nuclide', by_nuclide, bool)
if by_nuclide == True and self.domain_type == 'mesh':
if by_nuclide and self.domain_type == 'mesh':
raise ValueError('Unable to create MGXS library by nuclide with '
'mesh domain')
@ -277,7 +297,7 @@ class Library(object):
def domain_type(self, domain_type):
cv.check_value('domain type', domain_type, openmc.mgxs.DOMAIN_TYPES)
if self.by_nuclide == True and domain_type == 'mesh':
if self.by_nuclide and domain_type == 'mesh':
raise ValueError('Unable to create MGXS library by nuclide with '
'mesh domain')
@ -337,26 +357,72 @@ class Library(object):
def correction(self, correction):
cv.check_value('correction', correction, ('P0', None))
if correction == 'P0' and self.legendre_order > 0:
warn('The P0 correction will be ignored since the scattering '
'order "{}" is greater than zero'.format(self.legendre_order))
if self.scatter_format == 'legendre':
if correction == 'P0' and self.legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warn(msg)
elif self.scatter_format == 'histogram':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self._correction = correction
@scatter_format.setter
def scatter_format(self, scatter_format):
cv.check_value('scatter_format', scatter_format,
openmc.mgxs.MU_TREATMENTS)
if scatter_format == 'histogram' and self.correction == 'P0':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self.correction = None
self._scatter_format = scatter_format
@legendre_order.setter
def legendre_order(self, legendre_order):
cv.check_type('legendre_order', legendre_order, Integral)
cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
cv.check_greater_than('legendre_order', legendre_order, 0,
equality=True)
cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the scattering ' \
'order {} is greater than zero'.format(self.legendre_order)
warn(msg, RuntimeWarning)
self.correction = None
if self.scatter_format == 'legendre':
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warn(msg, RuntimeWarning)
self.correction = None
elif self.scatter_format == 'histogram':
msg = 'The legendre order will be ignored since the ' \
'scatter format is set to histogram'
warn(msg)
self._legendre_order = legendre_order
@histogram_bins.setter
def histogram_bins(self, histogram_bins):
cv.check_type('histogram_bins', histogram_bins, Integral)
cv.check_greater_than('histogram_bins', histogram_bins, 0)
if self.scatter_format == 'legendre':
msg = 'The histogram bins will be ignored since the ' \
'scatter format is set to legendre'
warn(msg)
elif self.scatter_format == 'histogram':
if self.correction == 'P0':
msg = 'The P0 correction will be ignored since ' \
'a histogram representation of the scattering '\
'kernel is requested'
warn(msg, RuntimeWarning)
self.correction = None
self._histogram_bins = histogram_bins
@tally_trigger.setter
def tally_trigger(self, tally_trigger):
cv.check_type('tally trigger', tally_trigger, openmc.Trigger)
@ -420,7 +486,7 @@ class Library(object):
mgxs.delayed_groups = None
else:
delayed_groups \
= list(range(1,self.num_delayed_groups+1))
= list(range(1, self.num_delayed_groups + 1))
mgxs.delayed_groups = delayed_groups
# If a tally trigger was specified, add it to the MGXS
@ -430,7 +496,9 @@ class Library(object):
# Specify whether to use a transport ('P0') correction
if isinstance(mgxs, openmc.mgxs.ScatterMatrixXS):
mgxs.correction = self.correction
mgxs.scatter_format = self.scatter_format
mgxs.legendre_order = self.legendre_order
mgxs.histogram_bins = self.histogram_bins
self.all_mgxs[domain.id][mgxs_type] = mgxs
@ -809,13 +877,13 @@ class Library(object):
return pickle.load(open(full_filename, 'rb'))
def get_xsdata(self, domain, xsdata_name, nuclide='total', xs_type='macro',
order=None, subdomain=None):
subdomain=None):
"""Generates an openmc.XSdata object describing a multi-group cross section
dataset for writing to an openmc.MGXSLibrary object.
Note that this method does not build an XSdata
object with nested temperature tables. The temperature of each
XSdata object will be left at the default value of 300K.
XSdata object will be left at the default value of 294K.
Parameters
----------
@ -830,10 +898,6 @@ class Library(object):
Provide the macro or micro cross section in units of cm^-1 or
barns. Defaults to 'macro'. If the Library object is not tallied by
nuclide this will be set to 'macro' regardless.
order : int
Scattering order for this data entry. Default is None,
which will set the XSdata object to use the order of the
Library.
subdomain : iterable of int
This parameter is not used unless using a mesh domain. In that
case, the subdomain is an [i,j,k] index (1-based indexing) of the
@ -863,10 +927,6 @@ class Library(object):
cv.check_type('xsdata_name', xsdata_name, string_types)
cv.check_type('nuclide', nuclide, string_types)
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
cv.check_type('order', order, (type(None), Integral))
if order is not None:
cv.check_greater_than('order', order, 0, equality=True)
cv.check_less_than('order', order, 10, equality=True)
if subdomain is not None:
cv.check_iterable_type('subdomain', subdomain, Integral,
max_depth=3)
@ -888,16 +948,7 @@ class Library(object):
xsdata = openmc.XSdata(name, self.energy_groups)
xsdata.num_delayed_groups = self.num_delayed_groups
if order is None:
# Set the order to the Library's order (the defualt behavior)
xsdata.order = self.legendre_order
else:
# Set the order of the xsdata object to the minimum of
# the provided order or the Library's order.
xsdata.order = min(order, self.legendre_order)
# Right now only 'legendre' data and isotropic weighting is supported
self.scatter_format = 'legendre'
# Right now only isotropic weighting is supported
self.representation = 'isotropic'
if nuclide != 'total':
@ -1037,10 +1088,17 @@ class Library(object):
# accounted for approximately by using an adjusted
# absorption cross section.
if 'total' in self.mgxs_types or 'transport' in self.mgxs_types:
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] \
= np.subtract(xsdata._total[i], np.sum(
xsdata._scatter_matrix[i][0, :, :], axis=1))
if xsdata.scatter_format == 'legendre':
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] = \
np.subtract(xsdata._total[i], np.sum(
xsdata._scatter_matrix[i][0, :, :], axis=1))
elif xsdata.scatter_format == 'histogram':
for i in range(len(xsdata.temperatures)):
xsdata._absorption[i] = \
np.subtract(xsdata._total[i], np.sum(np.sum(
xsdata._scatter_matrix[i][:, :, :], axis=0),
axis=1))
return xsdata
@ -1321,7 +1379,7 @@ class Library(object):
'scattering matrix is not provided.')
# Total or transport can be present, but if using
# self.correction=="P0", then we should use transport.
if (((self.correction is "P0") and
if (((self.correction == "P0") and
('nu-transport' not in self.mgxs_types))):
error_flag = True
warn('A "nu-transport" MGXS type is required since a "P0" '

257
openmc/mgxs/mgxs.py
View file

@ -59,6 +59,9 @@ _DOMAINS = (openmc.Cell,
openmc.Material,
openmc.Mesh)
# Supported ScatterMatrixXS and NuScatterMatrixXS angular distribution types
MU_TREATMENTS = ('legendre', 'histogram')
@add_metaclass(ABCMeta)
class MGXS(object):
@ -1526,9 +1529,19 @@ class MGXS(object):
else:
df = df.drop('score', axis=1)
# Determine if change-in-angle bins are included in the MGXS to
# properly tile the group boundaries
if 'mu low' in df:
# Find the length of the mu filters indirectly from the number
# of times the mu bins repeats.
num_mu = int(df.shape[0] /
df[df['mu low'] == df['mu low'][0]].shape[0])
else:
num_mu = 1
# Override energy groups bounds with indices
all_groups = np.arange(self.num_groups, 0, -1, dtype=np.int)
all_groups = np.repeat(all_groups, len(query_nuclides))
all_groups = np.repeat(all_groups, len(query_nuclides) * num_mu)
if 'energy low [eV]' in df and 'energyout low [eV]' in df:
df.rename(columns={'energy low [eV]': 'group in'},
inplace=True)
@ -3207,8 +3220,8 @@ class NuScatterXS(MGXS):
class ScatterMatrixXS(MatrixMGXS):
r"""A scattering matrix multi-group cross section for one or more Legendre
moments.
r"""A scattering matrix multi-group cross section with the cosine of the
change-in-angle represented as one or more Legendre moments or a histogram.
This class can be used for both OpenMC input generation and tally data
post-processing to compute spatially-homogenized and energy-integrated
@ -3226,7 +3239,7 @@ class ScatterMatrixXS(MatrixMGXS):
For a spatial domain :math:`V`, incoming energy group
:math:`[E_{g'},E_{g'-1}]`, and outgoing energy group :math:`[E_g,E_{g-1}]`,
the scattering moments are calculated as:
the Legendre scattering moments are calculated as:
.. math::
@ -3266,9 +3279,18 @@ class ScatterMatrixXS(MatrixMGXS):
Attributes
----------
correction : 'P0' or None
Apply the P0 correction to scattering matrices if set to 'P0'
Apply the P0 correction to scattering matrices if set to 'P0'; this is
used only if :attr:`ScatterMatrixXS.scatter_format` is 'legendre'
scatter_format : {'legendre', or 'histogram'}
Representation of the angular scattering distribution (default is
'legendre')
legendre_order : int
The highest Legendre moment in the scattering matrix (default is 0)
The highest Legendre moment in the scattering matrix; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'legendre'. (default is 0)
histogram_bins : int
The number of equally-spaced bins for the histogram representation of
the angular scattering distribution; this is used if
:attr:`ScatterMatrixXS.scatter_format` is 'histogram'. (default is 16)
name : str, optional
Name of the multi-group cross section
rxn_type : str
@ -3335,7 +3357,9 @@ class ScatterMatrixXS(MatrixMGXS):
groups, by_nuclide, name)
self._rxn_type = 'scatter'
self._correction = 'P0'
self._scatter_format = 'legendre'
self._legendre_order = 0
self._histogram_bins = 16
self._hdf5_key = 'scatter matrix'
self._estimator = 'analog'
self._valid_estimators = ['analog']
@ -3343,26 +3367,39 @@ class ScatterMatrixXS(MatrixMGXS):
def __deepcopy__(self, memo):
clone = super(ScatterMatrixXS, self).__deepcopy__(memo)
clone._correction = self.correction
clone._scatter_format = self.scatter_format
clone._legendre_order = self.legendre_order
clone._histogram_bins = self.histogram_bins
return clone
@property
def correction(self):
return self._correction
@property
def scatter_format(self):
return self._scatter_format
@property
def legendre_order(self):
return self._legendre_order
@property
def histogram_bins(self):
return self._histogram_bins
@property
def scores(self):
scores = ['flux']
if self.correction == 'P0' and self.legendre_order == 0:
scores += ['{}-0'.format(self.rxn_type),
'{}-1'.format(self.rxn_type)]
else:
scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)]
if self.scatter_format == 'legendre':
if self.correction == 'P0' and self.legendre_order == 0:
scores += ['{}-0'.format(self.rxn_type),
'{}-1'.format(self.rxn_type)]
else:
scores += ['{}-P{}'.format(self.rxn_type, self.legendre_order)]
elif self.scatter_format == 'histogram':
scores += [self.rxn_type]
return scores
@ -3372,10 +3409,15 @@ class ScatterMatrixXS(MatrixMGXS):
energy = openmc.EnergyFilter(group_edges)
energyout = openmc.EnergyoutFilter(group_edges)
if self.correction == 'P0' and self.legendre_order == 0:
filters = [[energy], [energy, energyout], [energyout]]
else:
filters = [[energy], [energy, energyout]]
if self.scatter_format == 'legendre':
if self.correction == 'P0' and self.legendre_order == 0:
filters = [[energy], [energy, energyout], [energyout]]
else:
filters = [[energy], [energy, energyout]]
elif self.scatter_format == 'histogram':
bins = np.linspace(-1., 1., num=self.histogram_bins + 1,
endpoint=True)
filters = [[energy], [energy, energyout, openmc.MuFilter(bins)]]
return filters
@ -3383,20 +3425,24 @@ class ScatterMatrixXS(MatrixMGXS):
def rxn_rate_tally(self):
if self._rxn_rate_tally is None:
if self.scatter_format == 'legendre':
# If using P0 correction subtract scatter-1 from the diagonal
if self.correction == 'P0' and self.legendre_order == 0:
scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
energy_filter = copy.deepcopy(energy_filter)
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
self._rxn_rate_tally = scatter_p0 - scatter_p1
# If using P0 correction subtract scatter-1 from the diagonal
if self.correction == 'P0' and self.legendre_order == 0:
scatter_p0 = self.tallies['{}-0'.format(self.rxn_type)]
scatter_p1 = self.tallies['{}-1'.format(self.rxn_type)]
energy_filter = scatter_p0.find_filter(openmc.EnergyFilter)
energy_filter = copy.deepcopy(energy_filter)
scatter_p1 = scatter_p1.diagonalize_filter(energy_filter)
self._rxn_rate_tally = scatter_p0 - scatter_p1
# Extract scattering moment reaction rate Tally
else:
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
self._rxn_rate_tally = self.tallies[tally_key]
# Extract scattering moment reaction rate Tally
else:
tally_key = '{}-P{}'.format(self.rxn_type,
self.legendre_order)
self._rxn_rate_tally = self.tallies[tally_key]
elif self.scatter_format == 'histogram':
# Extract scattering rate distribution tally
self._rxn_rate_tally = self.tallies[self.rxn_type]
self._rxn_rate_tally.sparse = self.sparse
@ -3406,27 +3452,52 @@ class ScatterMatrixXS(MatrixMGXS):
def correction(self, correction):
cv.check_value('correction', correction, ('P0', None))
if correction == 'P0' and self.legendre_order > 0:
msg = 'The P0 correction will be ignored since the scattering ' \
'order {} is greater than zero'.format(self.legendre_order)
if self.scatter_format == 'legendre':
if correction == 'P0' and self.legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warnings.warn(msg)
elif self.scatter_format == 'histogram':
msg = 'The P0 correction will be ignored since the ' \
'scatter format is set to histogram'
warnings.warn(msg)
self._correction = correction
@scatter_format.setter
def scatter_format(self, scatter_format):
cv.check_value('scatter_format', scatter_format, MU_TREATMENTS)
self._scatter_format = scatter_format
@legendre_order.setter
def legendre_order(self, legendre_order):
cv.check_type('legendre_order', legendre_order, Integral)
cv.check_greater_than('legendre_order', legendre_order, 0, equality=True)
cv.check_greater_than('legendre_order', legendre_order, 0,
equality=True)
cv.check_less_than('legendre_order', legendre_order, 10, equality=True)
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the scattering ' \
'order {} is greater than zero'.format(self.legendre_order)
warnings.warn(msg, RuntimeWarning)
self.correction = None
if self.scatter_format == 'legendre':
if self.correction == 'P0' and legendre_order > 0:
msg = 'The P0 correction will be ignored since the ' \
'scattering order {} is greater than '\
'zero'.format(self.legendre_order)
warnings.warn(msg, RuntimeWarning)
self.correction = None
elif self.scatter_format == 'histogram':
msg = 'The legendre order will be ignored since the ' \
'scatter format is set to histogram'
warnings.warn(msg)
self._legendre_order = legendre_order
@histogram_bins.setter
def histogram_bins(self, histogram_bins):
cv.check_type('histogram_bins', histogram_bins, Integral)
cv.check_greater_than('histogram_bins', histogram_bins, 0)
self._histogram_bins = histogram_bins
def load_from_statepoint(self, statepoint):
"""Extracts tallies in an OpenMC StatePoint with the data needed to
compute multi-group cross sections.
@ -3456,12 +3527,16 @@ class ScatterMatrixXS(MatrixMGXS):
self._rxn_rate_tally = None
self._loaded_sp = False
# Expand scores to match the format in the statepoint
# e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
if self.correction != 'P0' or self.legendre_order != 0:
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
self.tallies[tally_key].scores = \
[self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)]
if self.scatter_format == 'legendre':
# Expand scores to match the format in the statepoint
# e.g., "scatter-P2" -> "scatter-0", "scatter-1", "scatter-2"
if self.correction != 'P0' or self.legendre_order != 0:
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
self.tallies[tally_key].scores = \
[self.rxn_type + '-{}'.format(i)
for i in range(self.legendre_order + 1)]
elif self.scatter_format == 'histogram':
self.tallies[self.rxn_type].scores = [self.rxn_type]
super(ScatterMatrixXS, self).load_from_statepoint(statepoint)
@ -3508,7 +3583,7 @@ class ScatterMatrixXS(MatrixMGXS):
slice_xs._xs_tally = None
# Slice the Legendre order if needed
if legendre_order != 'same':
if legendre_order != 'same' and self.scatter_format == 'legendre':
cv.check_type('legendre_order', legendre_order, Integral)
cv.check_less_than('legendre_order', legendre_order,
self.legendre_order, equality=True)
@ -3517,7 +3592,8 @@ class ScatterMatrixXS(MatrixMGXS):
# Slice the scattering tally
tally_key = '{}-P{}'.format(self.rxn_type, self.legendre_order)
expand_scores = \
[self.rxn_type + '-{}'.format(i) for i in range(self.legendre_order+1)]
[self.rxn_type + '-{}'.format(i)
for i in range(self.legendre_order + 1)]
slice_xs.tallies[tally_key] = \
slice_xs.tallies[tally_key].get_slice(scores=expand_scores)
@ -3548,7 +3624,8 @@ class ScatterMatrixXS(MatrixMGXS):
This method constructs a 5D NumPy array for the requested
multi-group cross section data for one or more subdomains
(1st dimension), energy groups in (2nd dimension), energy groups out
(3rd dimension), nuclides (4th dimension), and moments (5th dimension).
(3rd dimension), nuclides (4th dimension), and moments/histograms
(5th dimension).
NOTE: The scattering moments are not multiplied by the :math:`(2l+1)/2`
prefactor in the expansion of the scattering source into Legendre
@ -3637,7 +3714,7 @@ class ScatterMatrixXS(MatrixMGXS):
filter_bins.append((self.energy_groups.get_group_bounds(group),))
# Construct CrossScore for requested scattering moment
if moment != 'all':
if moment != 'all' and self.scatter_format == 'legendre':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
@ -3687,9 +3764,19 @@ class ScatterMatrixXS(MatrixMGXS):
else:
num_out_groups = len(out_groups)
if self.scatter_format == 'histogram':
num_mu_bins = self.histogram_bins
else:
num_mu_bins = 1
# Reshape tally data array with separate axes for domain and energy
num_subdomains = int(xs.shape[0] / (num_in_groups * num_out_groups))
new_shape = (num_subdomains, num_in_groups, num_out_groups)
num_subdomains = int(xs.shape[0] /
(num_mu_bins * num_in_groups * num_out_groups))
if self.scatter_format == 'histogram':
new_shape = (num_subdomains, num_in_groups, num_out_groups,
num_mu_bins)
else:
new_shape = (num_subdomains, num_in_groups, num_out_groups)
new_shape += xs.shape[1:]
xs = np.reshape(xs, new_shape)
@ -3700,11 +3787,22 @@ class ScatterMatrixXS(MatrixMGXS):
# Reverse data if user requested increasing energy groups since
# tally data is stored in order of increasing energies
if order_groups == 'increasing':
xs = xs[:, ::-1, ::-1, :]
xs = xs[:, ::-1, ::-1, ...]
if squeeze:
xs = np.squeeze(xs)
xs = np.atleast_2d(xs)
# We want to squeeze out everything but the in_groups, out_groups,
# and, if needed, num_mu_bins dimension. These must not be squeezed
# so 1-group problems have the correct shape.
if self.scatter_format == 'histogram':
axes = (5, 4, 0)
else:
axes = (4, 3, 0)
# Squeeze will return a ValueError if the axis has a size greater
# than 1, so try each axis in axes one at a time, catching the
# ValueError as needed.
for axis in axes:
if xs.shape[axis] == 1:
xs = np.squeeze(xs, axis=axis)
return xs
@ -3755,28 +3853,41 @@ class ScatterMatrixXS(MatrixMGXS):
df = super(ScatterMatrixXS, self).get_pandas_dataframe(
groups, nuclides, xs_type, distribcell_paths)
# Add a moment column to dataframe
if self.legendre_order > 0:
# Insert a column corresponding to the Legendre moments
moments = ['P{}'.format(i) for i in range(self.legendre_order+1)]
moments = np.tile(moments, int(df.shape[0] / len(moments)))
df['moment'] = moments
if self.scatter_format == 'legendre':
# Add a moment column to dataframe
if self.legendre_order > 0:
# Insert a column corresponding to the Legendre moments
moments = ['P{}'.format(i)
for i in range(self.legendre_order + 1)]
moments = np.tile(moments, int(df.shape[0] / len(moments)))
df['moment'] = moments
# Place the moment column before the mean column
columns = df.columns.tolist()
mean_index = [i for i, s in enumerate(columns) if 'mean' in s][0]
if self.domain_type == 'mesh':
df = df[columns[:mean_index] + [('moment', '')] + columns[mean_index:-1]]
else:
df = df[columns[:mean_index] + ['moment'] + columns[mean_index:-1]]
# Place the moment column before the mean column
columns = df.columns.tolist()
mean_index \
= [i for i, s in enumerate(columns) if 'mean' in s][0]
if self.domain_type == 'mesh':
df = df[columns[:mean_index] + [('moment', '')] +
columns[mean_index:-1]]
else:
df = df[columns[:mean_index] + ['moment'] +
columns[mean_index:-1]]
# Select rows corresponding to requested scattering moment
if moment != 'all':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
df = df[df['moment'] == 'P{}'.format(moment)]
# Select rows corresponding to requested scattering moment
if moment != 'all':
cv.check_type('moment', moment, Integral)
cv.check_greater_than('moment', moment, 0, equality=True)
cv.check_less_than(
'moment', moment, self.legendre_order, equality=True)
df = df[df['moment'] == 'P{}'.format(moment)]
elif self.scatter_format == 'histogram':
# Replace the mu low and mu high columns with a single mu bin
del df['mu high']
df.rename(columns={'mu low': 'mu bins'}, inplace=True)
bins = [i + 1 for i in range(self.histogram_bins)]
bins = np.tile(bins, int(df.shape[0] / len(bins)))
df['mu bins'] = bins
return df
@ -3827,7 +3938,7 @@ class ScatterMatrixXS(MatrixMGXS):
cv.check_value('xs_type', xs_type, ['macro', 'micro'])
if self.correction != 'P0':
if self.correction != 'P0' and self.scatter_format == 'legendre':
rxn_type = '{0} (P{1})'.format(self.rxn_type, moment)
else:
rxn_type = self.rxn_type

208
openmc/mgxs_library.py
View file

@ -15,8 +15,8 @@ from openmc.checkvalue import check_type, check_value, check_greater_than, \
# Supported incoming particle MGXS angular treatment representations
_REPRESENTATIONS = ['isotropic', 'angle']
_SCATTER_TYPES = ['tabular', 'legendre', 'histogram']
_XS_SHAPES = ["[Order][G][G']", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
"[G'][DG]", "[G][G'][DG]"]
_XS_SHAPES = ["[G][G'][Order]", "[G]", "[G']", "[G][G']", "[DG]", "[G][DG]",
"[G'][DG]"]
class XSdata(object):
@ -134,7 +134,7 @@ class XSdata(object):
Note that some cross sections can be input in more than one shape so they
are listed multiple times:
[Order][G][G']: scatter_matrix
[G][G'][Order]: scatter_matrix
[G]: total, absorption, fission, kappa_fission, nu_fission,
prompt_nu_fission, delayed_nu_fission, inverse_velocity
@ -303,10 +303,9 @@ class XSdata(object):
self._xs_shapes["[G][G'][DG]"] = (self.energy_groups.num_groups,
self.energy_groups.num_groups,
self.num_delayed_groups)
self._xs_shapes["[Order][G][G']"] \
= (self.num_orders, self.energy_groups.num_groups,
self.energy_groups.num_groups)
self._xs_shapes["[G][G'][Order]"] \
= (self.energy_groups.num_groups,
self.energy_groups.num_groups, self.num_orders)
# If representation is by angle prepend num polar and num azim
if self.representation == 'angle':
@ -724,7 +723,7 @@ class XSdata(object):
"""
# Get the accepted shapes for this xs
shapes = [self.xs_shapes["[Order][G][G']"]]
shapes = [self.xs_shapes["[G][G'][Order]"]]
# Convert to a numpy array so we can easily get the shape for checking
scatter = np.asarray(scatter)
@ -1521,33 +1520,41 @@ class XSdata(object):
check_type('temperature', temperature, Real)
check_value('temperature', temperature, self.temperatures)
if self.scatter_format != 'legendre':
msg = 'Anisotropic scattering representations other than ' \
'Legendre expansions have not yet been implemented in ' \
'openmc.mgxs.'
raise ValueError(msg)
# Set the value of scatter_format based on the same value within
# scatter
self.scatter_format = scatter.scatter_format
# If the user has not defined XSdata.order, then we will set
# the order based on the data within scatter.
# Otherwise, we will check to see that XSdata.order to match
# Otherwise, we will check to see that XSdata.order matches
# the order of scatter
if self.order is None:
self.order = scatter.legendre_order
else:
check_value('legendre_order', scatter.legendre_order,
[self.order])
if self.scatter_format == 'legendre':
if self.order is None:
self.order = scatter.legendre_order
else:
check_value('legendre_order', scatter.legendre_order,
[self.order])
elif self.scatter_format == 'histogram':
if self.order is None:
self.order = scatter.histogram_bins
else:
check_value('histogram_bins', scatter.histogram_bins,
[self.order])
i = np.where(self.temperatures == temperature)[0][0]
if self.representation == 'isotropic':
# Get the scattering orders in the outermost dimension
self._scatter_matrix[i] = np.zeros((self.num_orders,
self.energy_groups.num_groups,
self.energy_groups.num_groups))
for moment in range(self.num_orders):
self._scatter_matrix[i][moment, :, :] = \
if self.scatter_format == 'legendre':
# Get the scattering orders in the outermost dimension
self._scatter_matrix[i] = \
np.zeros(self.xs_shapes["[G][G'][Order]"])
for moment in range(self.num_orders):
self._scatter_matrix[i][:, :, moment] = \
scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
moment=moment, subdomains=subdomain)
else:
self._scatter_matrix[i] = \
scatter.get_xs(nuclides=nuclide, xs_type=xs_type,
moment=moment, subdomains=subdomain)
subdomains=subdomain)
elif self.representation == 'angle':
msg = 'Angular-Dependent MGXS have not yet been implemented'
raise ValueError(msg)
@ -1625,6 +1632,10 @@ class XSdata(object):
scatt = scatter.get_xs(nuclides=nuclide,
xs_type=xs_type, moment=0,
subdomains=subdomain)
if scatter.scatter_format == 'histogram':
scatt = np.sum(scatt, axis=0)
if nuscatter.scatter_format == 'histogram':
nuscatt = np.sum(nuscatt, axis=0)
self._multiplicity_matrix[i] = np.divide(nuscatt, scatt)
elif self.representation == 'angle':
msg = 'Angular-Dependent MGXS have not yet been implemented'
@ -1656,7 +1667,7 @@ class XSdata(object):
if self.num_polar is not None:
grp.attrs['num_polar'] = self.num_polar
grp.attrs['scatter_shape'] = np.string_("[Order][G][G']")
grp.attrs['scatter_shape'] = np.string_("[G][G'][Order]")
if self.scatter_format is not None:
grp.attrs['scatter_format'] = np.string_(self.scatter_format)
if self.order is not None:
@ -1735,97 +1746,88 @@ class XSdata(object):
# Get the sparse scattering data to print to the library
G = self.energy_groups.num_groups
if self.representation == 'isotropic':
g_out_bounds = np.zeros((G, 2), dtype=np.int)
for g_in in range(G):
nz = np.nonzero(self._scatter_matrix[i][0, g_in, :])
g_out_bounds[g_in, 0] = nz[0][0]
g_out_bounds[g_in, 1] = nz[0][-1]
# Now create the flattened scatter matrix array
matrix = self._scatter_matrix[i]
flat_scatt = []
for g_in in range(G):
for g_out in range(g_out_bounds[g_in, 0],
g_out_bounds[g_in, 1] + 1):
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",
data=np.array(flat_scatt))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
matrix = self._multiplicity_matrix[i][:, :]
flat_mult = []
for g_in in range(G):
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",
data=np.array(flat_mult))
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :] += 1
scatt_grp.create_dataset("g_min", data=g_out_bounds[:, 0])
scatt_grp.create_dataset("g_max", data=g_out_bounds[:, 1])
Np = 1
Na = 1
elif self.representation == 'angle':
Np = self.num_polar
Na = self.num_azimuthal
g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
for p in range(Np):
for a in range(Na):
for g_in in range(G):
matrix = self._scatter_matrix[i][p, a, 0, g_in, :]
nz = np.nonzero(matrix)
g_out_bounds = np.zeros((Np, Na, G, 2), dtype=np.int)
for p in range(Np):
for a in range(Na):
for g_in in range(G):
if self.scatter_format == 'legendre':
if self.representation == 'isotropic':
matrix = \
self._scatter_matrix[i][g_in, :, 0]
elif self.representation == 'angle':
matrix = \
self._scatter_matrix[i][p, a, g_in, :, 0]
elif self.scatter_format == 'histogram':
if self.representation == 'isotropic':
matrix = \
np.sum(self._scatter_matrix[i][g_in, :, :],
axis=1)
elif self.representation == 'angle':
matrix = \
np.sum(self._scatter_matrix[i][p, a, g_in, :, :],
axis=1)
nz = np.nonzero(matrix)
# It is possible that there only zeros in matrix
# and therefore nz will be empty, in that case set
# g_out_bounds to 0s
if len(nz[0]) == 0:
g_out_bounds[p, a, g_in, :] = 0
else:
g_out_bounds[p, a, g_in, 0] = nz[0][0]
g_out_bounds[p, a, g_in, 1] = nz[0][-1]
# Now create the flattened scatter matrix array
flat_scatt = []
for p in range(Np):
for a in range(Na):
if self.representation == 'isotropic':
matrix = self._scatter_matrix[i][:, :, :]
elif self.representation == 'angle':
matrix = self._scatter_matrix[i][p, a, :, :, :]
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
for l in range(len(matrix[g_in, g_out, :])):
flat_scatt.append(matrix[g_in, g_out, l])
# And write it.
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:
# Now create the flattened scatter matrix array
flat_scatt = []
flat_mult = []
for p in range(Np):
for a in range(Na):
matrix = self._scatter_matrix[i][p, a, :, :, :]
if self.representation == 'isotropic':
matrix = self._multiplicity_matrix[i][:, :]
elif self.representation == 'angle':
matrix = self._multiplicity_matrix[i][p, a, :, :]
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
g_out_bounds[p, a, g_in, 1] + 1):
for l in range(len(matrix[:, g_in, g_out])):
flat_scatt.append(matrix[l, g_in, g_out])
flat_mult.append(matrix[g_in, g_out])
# And write it.
scatt_grp = xs_grp.create_group('scatter_data')
scatt_grp.create_dataset("scatter_matrix",
data=np.array(flat_scatt))
scatt_grp.create_dataset("multiplicity_matrix",
data=np.array(flat_mult))
# Repeat for multiplicity
if self._multiplicity_matrix[i] is not None:
# Now create the flattened scatter matrix array
flat_mult = []
for p in range(Np):
for a in range(Na):
matrix = self._multiplicity_matrix[i][p, a, :, :]
for g_in in range(G):
for g_out in range(g_out_bounds[p, a, g_in, 0],
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",
data=np.array(flat_mult))
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :, :, :] += 1
# And finally, adjust g_out_bounds for 1-based group counting
# and write it.
g_out_bounds[:, :, :, :] += 1
if self.representation == 'isotropic':
scatt_grp.create_dataset("g_min", data=g_out_bounds[0, 0, :, 0])
scatt_grp.create_dataset("g_max", data=g_out_bounds[0, 0, :, 1])
elif self.representation == 'angle':
scatt_grp.create_dataset("g_min", data=g_out_bounds[:, :, :, 0])
scatt_grp.create_dataset("g_max", data=g_out_bounds[:, :, :, 1])

3
openmc/tallies.py
View file

@ -1320,7 +1320,8 @@ class Tally(object):
# Create list of 2-tuples for energy boundary bins
elif isinstance(self_filter, (openmc.EnergyFilter,
openmc.EnergyoutFilter)):
openmc.EnergyoutFilter, openmc.MuFilter,
openmc.PolarFilter, openmc.AzimuthalFilter)):
bins = []
for k in range(self_filter.num_bins):
bins.append((self_filter.bins[k], self_filter.bins[k+1]))

40
scripts/openmc-update-mgxs
View file

@ -23,12 +23,10 @@ optional arguments:
from __future__ import print_function
import os
from shutil import move
import warnings
import xml.etree.ElementTree as ET
import argparse
import h5py
import numpy as np
import openmc.mgxs_library
@ -51,7 +49,7 @@ def parse_args():
if args['output'] == '':
filename = args['input'].name
extension = filenameos.path.splitext()
extension = os.path.splitext(filename)
if extension == '.xml':
filename = filename[:filename.rfind('.')] + '.h5'
args['output'] = filename
@ -84,6 +82,8 @@ if __name__ == '__main__':
temp = tree.find('group_structure').text.strip()
temp = np.array(temp.split())
group_structure = temp.astype(np.float)
# Convert from MeV to eV
group_structure *= 1.e6
energy_groups = openmc.mgxs.EnergyGroups(group_structure)
temp = tree.find('inverse-velocity')
if temp is not None:
@ -103,7 +103,7 @@ if __name__ == '__main__':
temperature = get_data(xsdata_elem, 'kT')
if temperature is not None:
temperature = \
float(temperature) / openmc.data.K_BOLTZMANN
float(temperature) / openmc.data.K_BOLTZMANN * 1.E6
else:
temperature = 294.
temperatures = [temperature]
@ -163,7 +163,7 @@ if __name__ == '__main__':
if temp is not None:
temp = np.array(temp.split(), dtype=float)
total = temp.astype(np.float)
total.shape = xsd[i].vector_shape
total.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_total(total, temperature)
if inverse_velocity is not None:
@ -172,41 +172,55 @@ if __name__ == '__main__':
temp = get_data(xsdata_elem, 'absorption')
temp = np.array(temp.split())
absorption = temp.astype(np.float)
absorption.shape = xsd[i].vector_shape
absorption.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_absorption(absorption, temperature)
temp = get_data(xsdata_elem, 'scatter')
temp = np.array(temp.split())
scatter = temp.astype(np.float)
scatter.shape = xsd[i].pn_matrix_shape
# This is now a flattened-array of something that started with a
# shape of [Order][G][G']; we need to unflatten and then switch the
# ordering
in_shape = (order_dim, energy_groups.num_groups,
energy_groups.num_groups)
if representation == 'angle':
in_shape = (n_pol, n_azi) + in_shape
scatter.shape = in_shape
scatter = np.swapaxes(scatter, 2, 3)
scatter = np.swapaxes(scatter, 3, 4)
else:
scatter.shape = in_shape
scatter = np.swapaxes(scatter, 0, 1)
scatter = np.swapaxes(scatter, 1, 2)
xsd[i].set_scatter_matrix(scatter, temperature)
temp = get_data(xsdata_elem, 'multiplicity')
if temp is not None:
temp = np.array(temp.split())
multiplicity = temp.astype(np.float)
multiplicity.shape = xsd[i].matrix_shape
multiplicity.shape = xsd[i].xs_shapes["[G][G']"]
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.shape = xsd[i].vector_shape
fission.shape = xsd[i].xs_shapes['[G]']
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.shape = xsd[i].vector_shape
kappa_fission.shape = xsd[i].xs_shapes['[G]']
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.shape = xsd[i].vector_shape
chi.shape = xsd[i].xs_shapes['[G]']
xsd[i].set_chi(chi, temperature)
else:
chi = None
@ -216,9 +230,9 @@ if __name__ == '__main__':
temp = np.array(temp.split())
nu_fission = temp.astype(np.float)
if chi is not None:
nu_fission.shape = xsd[i].vector_shape
nu_fission.shape = xsd[i].xs_shapes['[G]']
else:
nu_fission.shape = xsd[i].matrix_shape
nu_fission.shape = xsd[i].xs_shapes["[G][G']"]
xsd[i].set_nu_fission(nu_fission, temperature)
# Build library as we go, but first we have enough to initialize it

6
src/mgxs_header.F90
View file

@ -355,7 +355,7 @@ module mgxs_header
if (attribute_exists(xs_id, "scatter_shape")) then
call read_attribute(temp_str, xs_id, "scatter_shape")
temp_str = trim(temp_str)
if (to_lower(temp_str) /= "[order][g][g']") then
if (to_lower(temp_str) /= "[g][g'][order]") then
call fatal_error("Invalid scatter_shape option!")
end if
end if
@ -1995,8 +1995,8 @@ module mgxs_header
order_dim = order + 1
end if
! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for passing
! to ScattData
! Convert temp_1d to a jagged array ((gin) % data(l, gout)) for
! passing to ScattData
allocate(input_scatt(energy_groups, this % n_azi, this % n_pol))
index = 1

4
src/scattdata_header.F90
View file

@ -274,7 +274,7 @@ contains
! Get this by summing the un-normalized P0 coefficient in matrix
! over all outgoing groups
do gin = 1, groups
this % scattxs(gin) = sum(matrix(gin) % data(1, :), dim=1)
this % scattxs(gin) = sum(matrix(gin) % data(:, :))
end do
allocate(energy(groups))
@ -282,7 +282,7 @@ contains
! while also normalizing matrix itself (making CDF of f(mu=1)=1)
do gin = 1, groups
allocate(energy(gin) % data(gmin(gin):gmax(gin)))
do gout = 1, groups
do gout = gmin(gin), gmax(gin)
norm = sum(matrix(gin) % data(:, gout))
energy(gin) % data(gout) = norm
if (norm /= ZERO) then

BIN
tests/1d_mgxs.h5
View file

Binary file not shown.

2
tests/test_mg_tallies/results_true.dat
View file

@ -1 +1 @@
3d2ce1b8bdd558fe9f8560e8bb91455e1fedd80a47349344399e22f5b0fac1391f1721e7f12a42bb0c458fd6c87ebf9fa38fe9539fcf69292b21cebf8e5f8989
864328b2c4f3c4bfa9c80c756acbedac6fbdf3d502c03c2f2a3e7461b774729c47a55341c9515eff246e1bd445485f0d0a05a5712663ee499046afdbe0c77aef

2
tests/test_tallies/inputs_true.dat
View file

@ -1 +1 @@
b1f616bc0e342e3886b4b27eca22c4d1e55e445979254d5ba5f9ea01048f53d0e9a01c575745ae130eb34108080d2256ab9e4cb926ea0050891a44ae37ecdd88
ac118c5796593df0f53491920875bbe77921882253e7f24942e365e5812d726fb3d3bab86da02b3aad93f6f8f20a090125ffa51f9fdf9aef84009702f9cc363d

2
tests/test_tallies/results_true.dat
View file

@ -1 +1 @@
ceb1420d48e097185dd0798a495676c8e968c37c6e75fd6232137a4df8f6fae33d8232b95c68f6479815a014d20ccfc7611f96b22fcb7da159b42346e00cdf4e
a8172f7492fcc69d21f53c20523ffd0ad2882a753a3e1d263d135153cc6b3ee8d8e38f8823a4a2eb3288d730579892875d2d19f85f083dd659f91414988115ea

52
tests/test_tallies/test_tallies.py
View file

@ -3,12 +3,14 @@
import os
import sys
sys.path.insert(0, os.pardir)
from testing_harness import PyAPITestHarness
from openmc.filter import *
from openmc import Mesh, Tally, Tallies
from openmc.source import Source
from openmc.stats import Box
class TalliesTestHarness(PyAPITestHarness):
def _build_inputs(self):
# Build default materials/geometry
@ -21,33 +23,27 @@ class TalliesTestHarness(PyAPITestHarness):
self._input_set.settings.source = Source(space=Box(
[-160, -160, -183], [160, 160, 183]))
azimuthal_bins = (-3.1416, -1.8850, -0.6283, 0.6283, 1.8850, 3.1416)
azimuthal_filter1 = AzimuthalFilter(azimuthal_bins)
azimuthal_bins = (-3.14159, -1.8850, -0.6283, 0.6283, 1.8850, 3.14159)
azimuthal_filter = AzimuthalFilter(azimuthal_bins)
azimuthal_tally1 = Tally()
azimuthal_tally1.filters = [azimuthal_filter1]
azimuthal_tally1.filters = [azimuthal_filter]
azimuthal_tally1.scores = ['flux']
azimuthal_tally1.estimator = 'tracklength'
azimuthal_tally2 = Tally()
azimuthal_tally2.filters = [azimuthal_filter1]
azimuthal_tally2.filters = [azimuthal_filter]
azimuthal_tally2.scores = ['flux']
azimuthal_tally2.estimator = 'analog'
azimuthal_filter2 = AzimuthalFilter(5)
azimuthal_tally3 = Tally()
azimuthal_tally3.filters = [azimuthal_filter2]
azimuthal_tally3.scores = ['flux']
azimuthal_tally3.estimator = 'tracklength'
mesh_2x2 = Mesh(mesh_id=1)
mesh_2x2.lower_left = [-182.07, -182.07]
mesh_2x2.upper_right = [182.07, 182.07]
mesh_2x2.dimension = [2, 2]
mesh_filter = MeshFilter(mesh_2x2)
azimuthal_tally4 = Tally()
azimuthal_tally4.filters = [azimuthal_filter2, mesh_filter]
azimuthal_tally4.scores = ['flux']
azimuthal_tally4.estimator = 'tracklength'
azimuthal_tally3 = Tally()
azimuthal_tally3.filters = [azimuthal_filter, mesh_filter]
azimuthal_tally3.scores = ['flux']
azimuthal_tally3.estimator = 'tracklength'
cellborn_tally = Tally()
cellborn_tally.filters = [CellbornFilter((10, 21, 22, 23))]
@ -76,20 +72,17 @@ class TalliesTestHarness(PyAPITestHarness):
material_tally.filters = [MaterialFilter((1, 2, 3, 4))]
material_tally.scores = ['total']
mu_bins = (-1.0, -0.5, 0.0, 0.5, 1.0)
mu_filter = MuFilter(mu_bins)
mu_tally1 = Tally()
mu_tally1.filters = [MuFilter((-1.0, -0.5, 0.0, 0.5, 1.0))]
mu_tally1.filters = [mu_filter]
mu_tally1.scores = ['scatter', 'nu-scatter']
mu_filter = MuFilter(5)
mu_tally2 = Tally()
mu_tally2.filters = [mu_filter]
mu_tally2.filters = [mu_filter, mesh_filter]
mu_tally2.scores = ['scatter', 'nu-scatter']
mu_tally3 = Tally()
mu_tally3.filters = [mu_filter, mesh_filter]
mu_tally3.scores = ['scatter', 'nu-scatter']
polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.1416)
polar_bins = (0.0, 0.6283, 1.2566, 1.8850, 2.5132, 3.14159)
polar_filter = PolarFilter(polar_bins)
polar_tally1 = Tally()
polar_tally1.filters = [polar_filter]
@ -101,17 +94,11 @@ class TalliesTestHarness(PyAPITestHarness):
polar_tally2.scores = ['flux']
polar_tally2.estimator = 'analog'
polar_filter2 = PolarFilter((5,))
polar_tally3 = Tally()
polar_tally3.filters = [polar_filter2]
polar_tally3.filters = [polar_filter, mesh_filter]
polar_tally3.scores = ['flux']
polar_tally3.estimator = 'tracklength'
polar_tally4 = Tally()
polar_tally4.filters = [polar_filter2, mesh_filter]
polar_tally4.scores = ['flux']
polar_tally4.estimator = 'tracklength'
universe_tally = Tally()
universe_tally.filters = [UniverseFilter((1, 2, 3, 4, 6, 8))]
universe_tally.scores = ['total']
@ -173,10 +160,9 @@ class TalliesTestHarness(PyAPITestHarness):
self._input_set.tallies = Tallies()
self._input_set.tallies += (
[azimuthal_tally1, azimuthal_tally2, azimuthal_tally3,
azimuthal_tally4, cellborn_tally, dg_tally, energy_tally,
energyout_tally, transfer_tally, material_tally, mu_tally1,
mu_tally2, mu_tally3, polar_tally1, polar_tally2, polar_tally3,
polar_tally4, universe_tally])
cellborn_tally, dg_tally, energy_tally, energyout_tally,
transfer_tally, material_tally, mu_tally1, mu_tally2,
polar_tally1, polar_tally2, polar_tally3, universe_tally])
self._input_set.tallies += score_tallies
self._input_set.tallies += flux_tallies
self._input_set.tallies += (scatter_tally1, scatter_tally2)