resolved most of the comments from @samuelshaner

This commit is contained in:
Adam Nelson 2016-11-11 20:37:03 -05:00
parent 2a50aff588
commit d4f4166ba7
5 changed files with 110 additions and 132 deletions

View file

@ -5,6 +5,7 @@ import h5py
from openmc.mixin import EqualityMixin
from openmc.clean_xml import clean_xml_indentation
from openmc.checkvalue import check_type
class DataLibrary(EqualityMixin):
@ -95,13 +96,14 @@ class DataLibrary(EqualityMixin):
method='xml')
@classmethod
def from_xml(cls, path):
def from_xml(cls, path=None):
"""Read cross section data library from an XML file.
Parameters
----------
path : str
Path to XML file to read.
path : str, optional
Path to XML file to read. If not provided, the
`OPENMC_CROSS_SECTIONS` environment variable will be used.
Returns
-------
@ -112,6 +114,18 @@ class DataLibrary(EqualityMixin):
data = cls()
# If cross_sections is None, get the cross sections from the
# OPENMC_CROSS_SECTIONS environment variable
if path is None:
path = os.environ.get('OPENMC_CROSS_SECTIONS')
# Check to make sure there was an environmental variable.
if path is None:
raise ValueError("Either path or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
check_type('path', path, str)
tree = ET.parse(path)
root = tree.getroot()
if root.find('directory') is not None:

View file

@ -262,7 +262,7 @@ class Element(object):
return isotopes
def plot_xs(self, types, divisor_types=None, temperature=294., axis=None,
Erange=(1.E-5, 20.E6), sab_name=None, cross_sections=None,
energy_range=(1.E-5, 20.E6), sab_name=None, cross_sections=None,
enrichment=None, **kwargs):
"""Creates a figure of continuous-energy microscopic cross sections
for this element
@ -283,7 +283,7 @@ class Element(object):
axis : matplotlib.axes, optional
A previously generated axis to use for plotting. If not specified,
a new axis and figure will be generated.
Erange : tuple of floats
energy_range : tuple of floats
Energy range (in eV) to plot the cross section within
sab_name : str, optional
Name of S(a,b) library to apply to MT=2 data when applicable.
@ -324,12 +324,12 @@ class Element(object):
E = np.union1d(Enum, Ediv)
data_new = np.zeros((len(types), len(E)))
for l in range(len(types)):
data_new[l, :] = \
np.divide(np.interp(E, Enum, data[l, :]),
np.interp(E, Ediv, data_div[l, :]))
if divisor_types[l] != 'unity':
types[l] = types[l] + ' / ' + divisor_types[l]
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
@ -339,13 +339,14 @@ class Element(object):
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)):
# Set to loglog or semilogx depending on if we are plotting a data
# type which we expect to vary linearly
if types[i] in PLOT_TYPES_LINEAR:
plot_func = ax.semilogx
else:
plot_func = ax.loglog
if np.sum(data[i, :]) > 0.:
print('max = {:.2E}'.format(np.max(data[i, :])))
plot_func(E, data[i, :], label=types[i])
@ -356,7 +357,7 @@ class Element(object):
else:
ax.set_ylabel('Elemental Cross Section [b]')
ax.legend(loc='best')
ax.set_xlim(Erange)
ax.set_xlim(energy_range)
if self.name is not None:
title = 'Cross Section for ' + self.name
ax.set_title(title)
@ -388,11 +389,11 @@ class Element(object):
Returns
-------
unionE : numpy.array
energy_grid : numpy.array
Energies at which cross sections are calculated, in units of eV
data : numpy.ndarray
Macroscopic cross sections calculated at the energy grid described
by unionE
by energy_grid
"""
@ -402,18 +403,8 @@ class Element(object):
cv.check_iterable_type('types', types, str)
cv.check_type('temperature', temperature, Real)
# If cross_sections is None, get the cross sections from the
# OPENMC_CROSS_SECTIONS environment variable
if cross_sections is None:
cross_sections = os.environ.get('OPENMC_CROSS_SECTIONS')
# If a cross_sections library is present, check natural nuclides
# against the nuclides in the library
if cross_sections is not None:
library = openmc.data.DataLibrary.from_xml(cross_sections)
else:
raise ValueError("cross_sections or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
# Load the library
library = openmc.data.DataLibrary.from_xml(cross_sections)
# Expand elements in to nuclides with atomic densities
nuclides = self.expand(100., 'ao', enrichment=enrichment,
@ -443,17 +434,17 @@ class Element(object):
# Condense the data for every nuclide
# First create a union energy grid
unionE = E[0]
energy_grid = E[0]
for n in range(1, len(E)):
unionE = np.union1d(unionE, E[n])
energy_grid = np.union1d(energy_grid, E[n])
# Now we can combine all the nuclidic data
data = np.zeros((len(types), len(unionE)))
for l in range(len(types)):
if types[l] == 'unity':
data[l, :] = 1.
data = np.zeros((len(types), len(energy_grid)))
for line in range(len(types)):
if types[line] == 'unity':
data[line, :] = 1.
else:
for n in range(len(nuclides)):
data[l, :] += nuc_fractions[n] * xs[n][l](unionE)
data[line, :] += nuc_fractions[n] * xs[n][line](energy_grid)
return unionE, data
return energy_grid, data

View file

@ -627,20 +627,8 @@ class Material(object):
import scipy.constants as sc
cv.check_type('cross_sections', cross_sections, str)
# If cross_sections is None, get the cross sections from the
# OPENMC_CROSS_SECTIONS environment variable
if cross_sections is None:
cross_sections = os.environ.get('OPENMC_CROSS_SECTIONS')
# If a cross_sections library is present, check natural nuclides
# against the nuclides in the library
if cross_sections is not None:
library = openmc.data.DataLibrary.from_xml(cross_sections)
else:
raise ValueError("cross_sections or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
# Load the library
library = openmc.data.DataLibrary.from_xml(cross_sections)
# Expand elements in to nuclides
nuclides = self.get_nuclide_densities()
@ -704,15 +692,13 @@ class Material(object):
nuc_densities = density * nuc_densities
nuclides = OrderedDict()
n = -1
for nuc in nucs:
n += 1
for n, nuc in enumerate(nucs):
nuclides[nuc] = (nuc, nuc_densities[n])
return nuclides
def plot_xs(self, types, divisor_types=None, temperature=294.,
axis=None, Erange=(1.E-5, 20.E6), cross_sections=None,
axis=None, energy_range=(1.E-5, 20.E6), cross_sections=None,
**kwargs):
"""Creates a figure of continuous-energy macroscopic cross sections
for this material
@ -733,7 +719,7 @@ class Material(object):
axis : matplotlib.axes, optional
A previously generated axis to use for plotting. If not specified,
a new axis and figure will be generated.
Erange : tuple of floats, optional
energy_range : tuple of floats, optional
Energy range (in eV) to plot the cross section within
cross_sections : str, optional
Location of cross_sections.xml file. Default is None.
@ -767,12 +753,12 @@ class Material(object):
E = np.union1d(Enum, Ediv)
data_new = np.zeros((len(types), len(E)))
for l in range(len(types)):
data_new[l, :] = \
np.divide(np.interp(E, Enum, data[l, :]),
np.interp(E, Ediv, data_div[l, :]))
if divisor_types[l] != 'unity':
types[l] = types[l] + ' / ' + divisor_types[l]
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
@ -782,13 +768,14 @@ class Material(object):
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)):
# Set to loglog or semilogx depending on if we are plotting a data
# type which we expect to vary linearly
if types[i] in PLOT_TYPES_LINEAR:
plot_func = ax.semilogx
else:
plot_func = ax.loglog
if np.sum(data[i, :]) > 0.:
plot_func(E, data[i, :], label=types[i])
@ -798,7 +785,7 @@ class Material(object):
else:
ax.set_ylabel('Macroscopic Cross Section [1/cm]')
ax.legend(loc='best')
ax.set_xlim(Erange)
ax.set_xlim(energy_range)
if self.name is not None:
title = 'Macroscopic Cross Section for ' + self.name
ax.set_title(title)
@ -823,11 +810,11 @@ class Material(object):
Returns
-------
unionE : numpy.array
energy_grid : numpy.array
Energies at which cross sections are calculated, in units of eV
data : numpy.ndarray
Macroscopic cross sections calculated at the energy grid described
by unionE
by energy_grid
"""
@ -840,18 +827,8 @@ class Material(object):
cv.check_type('temperature', temperature, Real)
T = temperature
# If cross_sections is None, get the cross sections from the
# OPENMC_CROSS_SECTIONS environment variable
if cross_sections is None:
cross_sections = os.environ.get('OPENMC_CROSS_SECTIONS')
# If a cross_sections library is present, check natural nuclides
# against the nuclides in the library
if cross_sections is not None:
library = openmc.data.DataLibrary.from_xml(cross_sections)
else:
raise ValueError("cross_sections or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
# Load the library
library = openmc.data.DataLibrary.from_xml(cross_sections)
# Expand elements in to nuclides with atomic densities
nuclides = self.get_nuclide_atom_densities(cross_sections)
@ -881,20 +858,20 @@ class Material(object):
# Condense the data for every nuclide
# First create a union energy grid
unionE = E[0]
energy_grid = E[0]
for n in range(1, len(E)):
unionE = np.union1d(unionE, E[n])
energy_grid = np.union1d(energy_grid, E[n])
# Now we can combine all the nuclidic data
data = np.zeros((len(types), len(unionE)))
for l in range(len(types)):
if types[l] == 'unity':
data[l, :] = 1.
data = np.zeros((len(types), len(energy_grid)))
for line in range(len(types)):
if types[line] == 'unity':
data[line, :] = 1.
else:
for n in range(len(nuclides)):
data[l, :] += nuc_densities[n] * xs[n][l](unionE)
data[line, :] += nuc_densities[n] * xs[n][line](energy_grid)
return unionE, data
return energy_grid, data
def _get_nuclide_xml(self, nuclide, distrib=False):
xml_element = ET.Element("nuclide")

View file

@ -98,10 +98,9 @@ class Nuclide(object):
self._scattering = scattering
def plot_xs(self, types, divisor_types=None, temperature=294., axis=None,
Erange=(1.E-5, 20.E6), sab_name=None, cross_sections=None,
energy_range=(1.E-5, 20.E6), sab_name=None, cross_sections=None,
**kwargs):
"""Creates a figure of continuous-energy cross sections for this
nuclide
"""Creates a figure of continuous-energy cross sections for this nuclide
Parameters
----------
@ -119,7 +118,7 @@ class Nuclide(object):
axis : matplotlib.axes, optional
A previously generated axis to use for plotting. If not specified,
a new axis and figure will be generated.
Erange : tuple of floats
energy_range : tuple of floats
Energy range (in eV) to plot the cross section within
sab_name : str, optional
Name of S(a,b) library to apply to MT=2 data when applicable.
@ -156,11 +155,12 @@ class Nuclide(object):
E = np.union1d(Enum, Ediv)
data_new = []
for l in range(len(types)):
data_new.append(openmc.data.Combination([data[l], data_div[l]],
for line in range(len(types)):
data_new.append(openmc.data.Combination([data[line],
data_div[line]],
[np.divide]))
if divisor_types[l] != 'unity':
types[l] = types[l] + ' / ' + divisor_types[l]
if divisor_types[line] != 'unity':
types[line] = types[line] + ' / ' + divisor_types[line]
data = data_new
# Generate the plot
@ -170,14 +170,15 @@ class Nuclide(object):
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)):
to_plot = data[i](E)
# Set to loglog or semilogx depending on if we are plotting a data
# type which we expect to vary linearly
if types[i] in PLOT_TYPES_LINEAR:
plot_func = ax.semilogx
else:
plot_func = ax.loglog
if np.sum(to_plot) > 0.:
plot_func(E, to_plot, label=types[i])
@ -187,7 +188,7 @@ class Nuclide(object):
else:
ax.set_ylabel('Microscopic Cross Section [b]')
ax.legend(loc='best')
ax.set_xlim(Erange)
ax.set_xlim(energy_range)
if self.name is not None:
title = 'Microscopic Cross Section for ' + self.name
ax.set_title(title)
@ -216,10 +217,11 @@ class Nuclide(object):
Returns
-------
E : numpy.array
energy_grid : numpy.array
Energies at which cross sections are calculated, in units of eV
data : numpy.ndarray
Cross sections calculated at the energy grid described by unionE
Cross sections calculated at the energy grid described by
energy_grid
"""
@ -246,18 +248,8 @@ class Nuclide(object):
yields.append((False,))
ops.append(())
# If cross_sections is None, get the cross sections from the
# OPENMC_CROSS_SECTIONS environment variable
if cross_sections is None:
cross_sections = os.environ.get('OPENMC_CROSS_SECTIONS')
# If a cross_sections library is present, check natural nuclides
# against the nuclides in the library
if cross_sections is not None:
library = openmc.data.DataLibrary.from_xml(cross_sections)
else:
raise ValueError("cross_sections or OPENMC_CROSS_SECTIONS "
"environmental variable must be set")
# Load the library
library = openmc.data.DataLibrary.from_xml(cross_sections)
# Convert temperature to format needed for access in the library
cv.check_type('temperature', temperature, Real)
@ -265,7 +257,7 @@ class Nuclide(object):
T = temperature
# Now we can create the data sets to be plotted
E = []
energy_grid = []
xs = []
lib = library.get_by_material(self.name)
if lib is not None:
@ -325,9 +317,9 @@ class Nuclide(object):
if inelastic.x[-1] > sab_Emax:
sab_Emax = inelastic.x[-1]
sab_funcs.append(inelastic)
E = grid
energy_grid = grid
else:
E = nuc.energy[nucT]
energy_grid = nuc.energy[nucT]
for i, mt_set in enumerate(mts):
# Get the reaction xs data from the nuclide
@ -385,4 +377,4 @@ class Nuclide(object):
else:
raise ValueError(nuclide[0] + " not in library")
return E, xs
return energy_grid, xs

View file

@ -3,7 +3,7 @@ import numpy as np
# Supported keywords for material xs plotting
PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity',
'slowing-down power']
'slowing-down power', 'damage']
# Special MT values
UNITY_MT = -1
@ -44,7 +44,8 @@ PLOT_TYPES_MT = {'total': (2, 3,),
167, 168, 169, 170, 171, 172, 173, 174,
175, 176, 177, 178, 179, 180, 181, 183,
184, 190, 194, 196, 198, 199, 200, 875,
891, XI_MT)}
891, XI_MT),
'damage': (444,)}
# Operations to use when combining MTs the first np.add is used in reference
# to zero
PLOT_TYPES_OP = {'total': (np.add,),
@ -56,7 +57,8 @@ PLOT_TYPES_OP = {'total': (np.add,),
'nu-scatter': (np.add,) * (len(PLOT_TYPES_MT['nu-scatter']) - 1),
'unity': (),
'slowing-down power':
(np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,)}
(np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,),
'damage': ()}
# Whether or not to multiply the reaction by the yield as well
PLOT_TYPES_YIELD = {'total': (False, False),
@ -68,7 +70,9 @@ PLOT_TYPES_YIELD = {'total': (False, False),
'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']),
'unity': (False,),
'slowing-down power':
(True,) * len(PLOT_TYPES_MT['slowing-down power'])}
(True,) * len(PLOT_TYPES_MT['slowing-down power']),
'damage': (False,)}
# Types of plots to plot linearly in y
PLOT_TYPES_LINEAR = ['nu-fission / fission', 'nu-scatter / scatter']
PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter',
'nu-fission / absorption', 'fission / absorption'}