Some major simplifications: combined plot_xs and plot_mgxs so now the same interface works for both, and added support in a very generic way for all the different types of data one would expect to plot for MGXS plots (i.e., all the types in the library itself including choices of certain orders/delayed groups)

This commit is contained in:
Adam Nelson 2016-11-20 16:11:04 -05:00
parent 39dcd86e54
commit f76ee0867c

View file

@ -13,11 +13,14 @@ PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
# Supported keywoards for multi-group cross section plotting
PLOT_TYPES_MGXS = ['total', 'absorption', 'scatter', 'fission',
'kappa-fission', 'chi', 'chi-prompt', 'nu-fission',
'prompt-nu-fission', 'inverse-velocity', 'unity']
# Add on values for scattering moments
PLOT_TYPES_MGXS += ['scatter-' + str(i)
for i in range(0, openmc.mgxs.MAX_LEGENDRE + 1)]
'kappa-fission', 'nu-fission', 'prompt-nu-fission',
'deleyed-nu-fission', 'chi', 'chi-prompt', 'chi-delayed',
'inverse-velocity', 'beta', 'decay rate', 'unity']
# Create a dictionary which can be used to convert PLOT_TYPES_MGXS to the
# openmc.XSdata attribute name needed to access the data
_PLOT_ATTR = {line: line.replace(' ', '_').replace('-', '_')
for line in PLOT_TYPES_MGXS}
_PLOT_ATTR['scatter'] = 'scatter_matrix'
# Special MT values
UNITY_MT = -1
@ -56,7 +59,9 @@ PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter',
def plot_xs(this, types, divisor_types=None, temperature=294., axis=None,
sab_name=None, cross_sections=None, enrichment=None, **kwargs):
sab_name=None, ce_cross_sections=None, mg_cross_sections=None,
enrichment=None, plot_CE=True, orders=None, divisor_orders=None,
**kwargs):
"""Creates a figure of continuous-energy cross sections for this item
Parameters
@ -80,138 +85,24 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None,
sab_name : str, optional
Name of S(a,b) library to apply to MT=2 data when applicable; only used
for items which are instances of openmc.Element or openmc.Nuclide
cross_sections : str, optional
ce_cross_sections : str, optional
Location of cross_sections.xml file. Default is None.
enrichment : float, optional
Enrichment for U235 in weight percent. For example, input 4.95 for
4.95 weight percent enriched U. Default is None. This is only used for
items which are instances of openmc.Element
**kwargs
All keyword arguments are passed to
:func:`matplotlib.pyplot.figure`.
Returns
-------
fig : matplotlib.figure.Figure
If axis is None, then a Matplotlib Figure of the generated
cross section will be returned. Otherwise, a value of
None will be returned as the figure and axes have already been
generated.
"""
from matplotlib import pyplot as plt
if isinstance(this, openmc.Nuclide):
data_type = 'nuclide'
elif isinstance(this, openmc.Element):
data_type = 'element'
elif isinstance(this, openmc.Material):
data_type = 'material'
else:
raise TypeError("Invalid type for plotting")
E, data = calculate_xs(this, types, temperature, sab_name, cross_sections,
enrichment)
if divisor_types:
cv.check_length('divisor types', divisor_types, len(types),
len(types))
Ediv, data_div = calculate_xs(this, divisor_types, temperature,
sab_name, cross_sections, enrichment)
# Create a new union grid, interpolate data and data_div on to that
# grid, and then do the actual division
Enum = E[:]
E = np.union1d(Enum, Ediv)
data_new = np.zeros((len(types), len(E)))
for line in range(len(types)):
data_new[line, :] = \
np.divide(np.interp(E, Enum, data[line, :]),
np.interp(E, Ediv, data_div[line, :]))
if divisor_types[line] != 'unity':
types[line] = types[line] + ' / ' + divisor_types[line]
data = data_new
# Generate the plot
if axis is None:
fig = plt.figure(**kwargs)
ax = fig.add_subplot(111)
else:
fig = None
ax = axis
# Set to loglog or semilogx depending on if we are plotting a data
# type which we expect to vary linearly
if set(types).issubset(PLOT_TYPES_LINEAR):
plot_func = ax.semilogx
else:
plot_func = ax.loglog
# Plot the data
for i in range(len(data)):
data[i, :] = np.nan_to_num(data[i, :])
if np.sum(data[i, :]) > 0.:
plot_func(E, data[i, :], label=types[i])
ax.set_xlabel('Energy [eV]')
if divisor_types:
if data_type == 'nuclide':
ylabel = 'Nuclidic Microscopic Data'
elif data_type == 'element':
ylabel = 'Elemental Microscopic Data'
elif data_type == 'material':
ylabel = 'Macroscopic Data'
else:
if data_type == 'nuclide':
ylabel = 'Microscopic Cross Section [b]'
elif data_type == 'element':
ylabel = 'Elemental Cross Section [b]'
elif data_type == 'material':
ylabel = 'Macroscopic Cross Section [1/cm]'
ax.set_ylabel(ylabel)
ax.legend(loc='best')
# Set to the most likely expected range
ax.set_xlim((1.E-5, 20.E6))
if this.name is not None:
ax.set_title('Cross Section for ' + this.name)
return fig
def plot_mgxs(this, types, divisor_types=None, temperature=294., axis=None,
energy_range=(1.E-5, 20.E6), cross_sections=None,
ce_cross_sections=None, enrichment=None, **kwargs):
"""Creates a figure of multi-group cross sections for this item
Parameters
----------
this : {openmc.Element, openmc.Nuclide, openmc.Material, openmc.Macroscopic}
Object to source data from
types : Iterable of values of PLOT_TYPES
The type of cross sections to include in the plot.
divisor_types : Iterable of values of PLOT_TYPES, optional
Cross section types which will divide those produced by types
before plotting. A type of 'unity' can be used to effectively not
divide some types.
temperature : float, optional
Temperature in Kelvin to plot. If not specified, a default
temperature of 294K will be plotted. Note that the nearest
temperature in the library for each nuclide will be used as opposed
to using any interpolation.
axis : matplotlib.axes, optional
A previously generated axis to use for plotting. If not specified,
a new axis and figure will be generated.
energy_range : tuple of floats
Energy range (in eV) to plot the cross section within
cross_sections : str, optional
mg_cross_sections : str, optional
Location of MGXS HDF5 Library file. Default is None.
cross_sections : str, optional
Location of continuous-energy cross_sections.xml file. Default is None.
This is used only for expanding an openmc.Element object passed as this
enrichment : float, optional
Enrichment for U235 in weight percent. For example, input 4.95 for
4.95 weight percent enriched U. Default is None. This is only used for
items which are instances of openmc.Element
plot_CE : bool
Denotes whether or not continuous-energy will be plotted. Defaults to
plotting the continuous-energy data.
orders : Iterable of Integral
The scattering order or delayed group index to use for the
corresponding entry in types. Defaults to the 0th order for scattering
and the total delayed neutron data. This only applies to plots of
multi-group data.
divisor_orders : Iterable of Integral
Same as orders, but for divisor_types
**kwargs
All keyword arguments are passed to
:func:`matplotlib.pyplot.figure`.
@ -228,6 +119,8 @@ def plot_mgxs(this, types, divisor_types=None, temperature=294., axis=None,
from matplotlib import pyplot as plt
cv.check_type("plot_CE", plot_CE, bool)
if isinstance(this, openmc.Nuclide):
data_type = 'nuclide'
elif isinstance(this, openmc.Element):
@ -239,21 +132,49 @@ def plot_mgxs(this, types, divisor_types=None, temperature=294., axis=None,
else:
raise TypeError("Invalid type for plotting")
E, data = calculate_mgxs(this, types, temperature, cross_sections,
ce_cross_sections, enrichment)
if plot_CE:
# Calculate for the CE cross sections
E, data = calculate_cexs(this, types, temperature, sab_name,
ce_cross_sections, enrichment)
if divisor_types:
cv.check_length('divisor types', divisor_types, len(types),
len(types))
Ediv, data_div = calculate_cexs(this, divisor_types, temperature,
sab_name, ce_cross_sections,
enrichment)
if divisor_types:
cv.check_length('divisor types', divisor_types, len(types),
len(types))
Ediv, data_div = calculate_mgxs(this, divisor_types, temperature,
cross_sections, ce_cross_sections,
enrichment)
# Create a new union grid, interpolate data and data_div on to that
# grid, and then do the actual division
Enum = E[:]
E = np.union1d(Enum, Ediv)
data_new = np.zeros((len(types), len(E)))
# Perform the division
for line in range(len(types)):
data[line, :] = np.divide(data[line, :], data_div[line, :])
if divisor_types[line] != 'unity':
types[line] = types[line] + ' / ' + divisor_types[line]
for line in range(len(types)):
data_new[line, :] = \
np.divide(np.interp(E, Enum, data[line, :]),
np.interp(E, Ediv, data_div[line, :]))
if divisor_types[line] != 'unity':
types[line] = types[line] + ' / ' + divisor_types[line]
data = data_new
else:
# Calculate for MG cross sections
E, data = calculate_mgxs(this, types, orders, temperature,
mg_cross_sections, ce_cross_sections,
enrichment)
if divisor_types:
cv.check_length('divisor types', divisor_types, len(types),
len(types))
Ediv, data_div = calculate_mgxs(this, divisor_types,
divisor_orders, temperature,
mg_cross_sections,
ce_cross_sections, enrichment)
# Perform the division
for line in range(len(types)):
data[line, :] = np.divide(data[line, :], data_div[line, :])
if divisor_types[line] != 'unity':
types[line] = types[line] + ' / ' + divisor_types[line]
# Generate the plot
if axis is None:
@ -268,9 +189,12 @@ def plot_mgxs(this, types, divisor_types=None, temperature=294., axis=None,
plot_func = ax.semilogx
else:
plot_func = ax.loglog
# Plot the data
for i in range(len(data)):
plot_func(E, data[i], label=types[i])
data[i, :] = np.nan_to_num(data[i, :])
if np.sum(data[i, :]) > 0.:
plot_func(E, data[i, :], label=types[i])
ax.set_xlabel('Energy [eV]')
if divisor_types:
@ -289,16 +213,14 @@ def plot_mgxs(this, types, divisor_types=None, temperature=294., axis=None,
ylabel = 'Macroscopic Cross Section [1/cm]'
ax.set_ylabel(ylabel)
ax.legend(loc='best')
ax.set_xlim(energy_range)
if this.name is not None:
title = 'Cross Section for ' + this.name
ax.set_title(title)
ax.set_title('Cross Section for ' + this.name)
return fig
def calculate_xs(this, types, temperature=294., sab_name=None,
cross_sections=None, enrichment=None):
def calculate_cexs(this, types, temperature=294., sab_name=None,
cross_sections=None, enrichment=None):
"""Calculates continuous-energy cross sections of a requested type
Parameters
@ -338,8 +260,8 @@ def calculate_xs(this, types, temperature=294., sab_name=None,
cv.check_type('enrichment', enrichment, Real)
if isinstance(this, openmc.Nuclide):
energy_grid, xs = _calculate_xs_nuclide(this, types, temperature,
sab_name, cross_sections)
energy_grid, xs = _calculate_cexs_nuclide(this, types, temperature,
sab_name, cross_sections)
# Convert xs (Iterable of Callable) to a grid of cross section values
# calculated on @ the points in energy_grid for consistency with the
# element and material functions.
@ -347,20 +269,20 @@ def calculate_xs(this, types, temperature=294., sab_name=None,
for line in range(len(types)):
data[line, :] = xs[line](energy_grid)
elif isinstance(this, openmc.Element):
energy_grid, data = _calculate_xs_elem_mat(this, types, temperature,
cross_sections, sab_name,
enrichment)
energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature,
cross_sections, sab_name,
enrichment)
elif isinstance(this, openmc.Material):
energy_grid, data = _calculate_xs_elem_mat(this, types, temperature,
cross_sections)
energy_grid, data = _calculate_cexs_elem_mat(this, types, temperature,
cross_sections)
else:
raise TypeError("Invalid type")
return energy_grid, data
def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None,
cross_sections=None):
def _calculate_cexs_nuclide(this, types, temperature=294., sab_name=None,
cross_sections=None):
"""Calculates continuous-energy cross sections of a requested type
Parameters
@ -548,8 +470,9 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None,
return energy_grid, xs
def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None,
sab_name=None, enrichment=None):
def _calculate_cexs_elem_mat(this, types, temperature=294.,
cross_sections=None, sab_name=None,
enrichment=None):
"""Calculates continuous-energy cross sections of a requested type
Parameters
@ -637,7 +560,8 @@ def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None,
name = nuclide[0]
nuc = nuclide[1]
sab_tab = sabs[name]
temp_E, temp_xs = calculate_xs(nuc, types, T, sab_tab, cross_sections)
temp_E, temp_xs = calculate_cexs(nuc, types, T, sab_tab,
cross_sections)
E.append(temp_E)
# Since the energy grids are different, store the cross sections as
# a tabulated function so they can be calculated on any grid needed.
@ -664,8 +588,9 @@ def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None,
return energy_grid, data
def calculate_mgxs(this, types, temperature=294., cross_sections=None,
ce_cross_sections=None, enrichment=None):
def calculate_mgxs(this, types, orders=None, temperature=294.,
cross_sections=None, ce_cross_sections=None,
enrichment=None):
"""Calculates continuous-energy cross sections of a requested type
Parameters
@ -674,6 +599,10 @@ def calculate_mgxs(this, types, temperature=294., cross_sections=None,
Object to source data from
types : Iterable of values of PLOT_TYPES
The type of cross sections to calculate
orders : Iterable of Integral
The scattering order or delayed group index to use for the
corresponding entry in types. Defaults to the 0th order for scattering
and the total delayed neutron data.
temperature : float, optional
Temperature in Kelvin to plot. If not specified, a default
temperature of 294K will be plotted. Note that the nearest
@ -702,15 +631,18 @@ def calculate_mgxs(this, types, temperature=294., cross_sections=None,
cv.check_type('temperature', temperature, Real)
if enrichment:
cv.check_type('enrichment', enrichment, Real)
cv.check_iterable_type('types', types, string_types)
cv.check_type("cross_sections", cross_sections, str)
library = openmc.MGXSLibrary.from_hdf5(cross_sections)
if isinstance(this, (openmc.Nuclide, openmc.Macroscopic)):
mgxs = _calculate_mgxs_nuc_macro(this, types, library, temperature)
mgxs = _calculate_mgxs_nuc_macro(this, types, library, orders,
temperature)
elif isinstance(this, (openmc.Element, openmc.Material)):
mgxs = _calculate_mgxs_elem_mat(this, types, library, temperature,
ce_cross_sections, enrichment)
mgxs = _calculate_mgxs_elem_mat(this, types, library, orders,
temperature, ce_cross_sections,
enrichment)
else:
raise TypeError("Invalid type")
@ -727,7 +659,8 @@ def calculate_mgxs(this, types, temperature=294., cross_sections=None,
return np.flipud(energy_grid), data
def _calculate_mgxs_nuc_macro(this, types, library, temperature=294.):
def _calculate_mgxs_nuc_macro(this, types, library, orders=None,
temperature=294.):
"""Determines the multi-group cross sections of a nuclide or macroscopic
object
@ -740,6 +673,10 @@ def _calculate_mgxs_nuc_macro(this, types, library, temperature=294.):
in openmc.PLOT_TYPES_MGXS
library : openmc.MGXSLibrary
MGXS Library containing the data of interest
orders : Iterable of Integral
The scattering order or delayed group index to use for the
corresponding entry in types. Defaults to the 0th order for scattering
and the total delayed neutron data.
temperature : float, optional
Temperature in Kelvin to plot. If not specified, a default
temperature of 294K will be plotted. Note that the nearest
@ -753,10 +690,17 @@ def _calculate_mgxs_nuc_macro(this, types, library, temperature=294.):
"""
# Check the parameters
for line in types:
# Check the parameters and grab order/delayed groups
if orders:
cv.check_iterable_type('orders', orders, Integral,
min_depth=len(types), max_depth=len(types))
else:
orders = [None] * len(types)
for i, line in enumerate(types):
cv.check_type("line", line, str)
cv.check_value("line", line, PLOT_TYPES_MGXS)
if orders[i]:
cv.check_greater_than("order value", orders[i], 0, equality=True)
xsdata = library[this.name]
@ -767,26 +711,46 @@ def _calculate_mgxs_nuc_macro(this, types, library, temperature=294.):
# Get the data
data = np.zeros((len(types), library.energy_groups.num_groups))
for i, line in enumerate(types):
if line == 'unity':
if 'fission' in line and not xsdata.fissionable:
continue
elif line == 'unity':
data[i, :] = 1.
elif line.startswith('scatter'):
# We have to remove the outgoing dependence
attr = line.replace(' ', '_').replace('-', '_')
matrix = xsdata.scatter_matrix[t]
# Sum over outgoing groups
vector = np.sum(matrix, axis=1)
# Now get the actual order of interest
if line == 'scatter':
order = 0
else:
order = int(line.split('-')[1])
if order < xsdata.xs_shapes["[G][G'][Order]"][-1]:
data[i, :] = vector[:, order]
else:
data[i, :] = 0.
else:
attr = line.replace(' ', '_').replace('-', '_')
data[i, :] = getattr(xsdata, attr)[t]
attr = _PLOT_ATTR[line]
temp_data = getattr(xsdata, attr)[t]
if temp_data.shape in (xsdata.xs_shapes["[G']"],
xsdata.xs_shapes["[G]"]):
data[i, :] = temp_data
elif temp_data.shape == xsdata.xs_shapes["[G][G']"]:
data[i, :] = np.sum(temp_data, axis=1)
elif temp_data.shape == xsdata.xs_shapes["[DG]"]:
if orders[i]:
if orders[i] < len(temp_data.shape[0]):
data[i, :] = temp_data[orders[i]]
else:
data[i, :] = np.sum(temp_data[:])
elif temp_data.shape in (xsdata.xs_shapes["[G'][DG]"],
xsdata.xs_shapes["[G][DG]"]):
if orders[i]:
if orders[i] < len(temp_data.shape[1]):
data[i, :] = temp_data[:, orders[i]]
else:
data[i, :] = np.sum(temp_data[:, :], axis=1)
elif temp_data.shape == xsdata.xs_shapes["[G][G'][DG]"]:
temp_data = np.sum(temp_data, axis=1)
if orders[i]:
if orders[i] < len(temp_data.shape[1]):
data[i, :] = temp_data[:, orders[i]]
else:
data[i, :] = np.sum(temp_data[:, :], axis=1)
elif temp_data.shape == xsdata.xs_shapes["[G][G'][Order]"]:
temp_data = np.sum(temp_data, axis=1)
if orders[i]:
order = orders[i]
else:
order = 0
if order < temp_data.shape[1]:
data[i, :] = temp_data[:, order]
else:
raise ValueError("{} not present in provided MGXS "
"library".format(this.name))
@ -794,8 +758,9 @@ def _calculate_mgxs_nuc_macro(this, types, library, temperature=294.):
return data
def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
ce_cross_sections=None, enrichment=None):
def _calculate_mgxs_elem_mat(this, types, library, orders=None,
temperature=294., ce_cross_sections=None,
enrichment=None):
"""Determines the multi-group cross sections of an element or material
object
@ -808,6 +773,10 @@ def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
in openmc.PLOT_TYPES_MGXS
library : openmc.MGXSLibrary
MGXS Library containing the data of interest
orders : Iterable of Integral
The scattering order or delayed group index to use for the
corresponding entry in types. Defaults to the 0th order for scattering
and the total delayed neutron data.
temperature : float, optional
Temperature in Kelvin to plot. If not specified, a default
temperature of 294K will be plotted. Note that the nearest
@ -828,11 +797,6 @@ def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
"""
# Check the parameters
for line in types:
cv.check_type("line", line, str)
cv.check_value("line", line, PLOT_TYPES_MGXS)
if isinstance(this, openmc.Material):
if this.temperature is not None:
T = this.temperature
@ -861,7 +825,7 @@ def _calculate_mgxs_elem_mat(this, types, library, temperature=294.,
nuc_data = []
for nuclide in nuclides.items():
nuc_data.append(_calculate_mgxs_nuc_macro(nuclide[0], types, library,
T))
orders, T))
# Combine across the nuclides
data = np.zeros((len(types), library.energy_groups.num_groups))