mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-27 13:45:36 -04:00
Merge remote-tracking branch 'upstream/develop' into angles
This commit is contained in:
Adam Nelson
commit
906270b4fc
8 changed files with 652 additions and 14 deletions
|
|
@ -25,6 +25,7 @@ from openmc.statepoint import *
|
|||
from openmc.summary import *
|
||||
from openmc.particle_restart import *
|
||||
from openmc.mixin import *
|
||||
from openmc.plotter import *
|
||||
|
||||
try:
|
||||
from openmc.opencg_compatible import *
|
||||
|
|
|
|||
|
|
@ -186,3 +186,9 @@ K_BOLTZMANN = 8.6173324e-5
|
|||
|
||||
# Used for converting units in ACE data
|
||||
EV_PER_MEV = 1.0e6
|
||||
|
||||
# Avogadro's constant from CODATA 2010
|
||||
AVOGADRO = 6.02214129E23
|
||||
|
||||
# Neutron mass from CODATA 2010 in units of amu
|
||||
NEUTRON_MASS = 1.008664916
|
||||
|
|
|
|||
|
|
@ -1,10 +1,12 @@
|
|||
import os
|
||||
import xml.etree.ElementTree as ET
|
||||
from six import string_types
|
||||
|
||||
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 +97,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 +115,18 @@ class DataLibrary(EqualityMixin):
|
|||
|
||||
data = cls()
|
||||
|
||||
# If path 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, string_types)
|
||||
|
||||
tree = ET.parse(path)
|
||||
root = tree.getroot()
|
||||
if root.find('directory') is not None:
|
||||
|
|
|
|||
|
|
@ -511,6 +511,9 @@ class IncidentNeutron(EqualityMixin):
|
|||
rxs = [data[mt] for mt in SUM_RULES[mt_sum] if mt in data]
|
||||
if len(rxs) > 0:
|
||||
data.summed_reactions[mt_sum] = rx = Reaction(mt_sum)
|
||||
if rx.mt == 18 and 'total_nu' in group:
|
||||
tgroup = group['total_nu']
|
||||
rx.derived_products.append(Product.from_hdf5(tgroup))
|
||||
for T in data.temperatures:
|
||||
rx.xs[T] = Sum([rx_i.xs[T] for rx_i in rxs])
|
||||
|
||||
|
|
|
|||
|
|
@ -1,13 +1,12 @@
|
|||
from collections import OrderedDict
|
||||
import re
|
||||
import sys
|
||||
import os
|
||||
|
||||
from six import string_types
|
||||
from xml.etree import ElementTree as ET
|
||||
|
||||
import openmc
|
||||
from openmc.checkvalue import check_type, check_length
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data import NATURAL_ABUNDANCE, atomic_mass
|
||||
|
||||
|
||||
|
|
@ -80,8 +79,8 @@ class Element(object):
|
|||
|
||||
@name.setter
|
||||
def name(self, name):
|
||||
check_type('element name', name, string_types)
|
||||
check_length('element name', name, 1, 2)
|
||||
cv.check_type('element name', name, string_types)
|
||||
cv.check_length('element name', name, 1, 2)
|
||||
self._name = name
|
||||
|
||||
@scattering.setter
|
||||
|
|
@ -254,6 +253,6 @@ class Element(object):
|
|||
for nuclide, abundance in abundances.items():
|
||||
nuc = openmc.Nuclide(nuclide)
|
||||
nuc.scattering = self.scattering
|
||||
isotopes.append((nuc, percent*abundance, percent_type))
|
||||
isotopes.append((nuc, percent * abundance, percent_type))
|
||||
|
||||
return isotopes
|
||||
|
|
|
|||
|
|
@ -3,9 +3,9 @@ from copy import deepcopy
|
|||
from numbers import Real, Integral
|
||||
import warnings
|
||||
from xml.etree import ElementTree as ET
|
||||
import sys
|
||||
|
||||
from six import string_types
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
import openmc.data
|
||||
|
|
@ -595,7 +595,7 @@ class Material(object):
|
|||
nuclides = OrderedDict()
|
||||
|
||||
for nuclide, density, density_type in self._nuclides:
|
||||
nuclides[nuclide.name] = (nuclide, density)
|
||||
nuclides[nuclide.name] = (nuclide, density, density_type)
|
||||
|
||||
for ele, ele_pct, ele_pct_type, enr in self._elements:
|
||||
|
||||
|
|
@ -606,6 +606,80 @@ class Material(object):
|
|||
|
||||
return nuclides
|
||||
|
||||
def get_nuclide_atom_densities(self):
|
||||
"""Returns all nuclides in the material and their atomic densities in
|
||||
units of atom/b-cm
|
||||
|
||||
Returns
|
||||
-------
|
||||
nuclides : dict
|
||||
Dictionary whose keys are nuclide names and values are tuples of
|
||||
(nuclide, density in atom/b-cm)
|
||||
|
||||
"""
|
||||
|
||||
# Expand elements in to nuclides
|
||||
nuclides = self.get_nuclide_densities()
|
||||
|
||||
sum_density = False
|
||||
if self.density_units == 'sum':
|
||||
sum_density = True
|
||||
density = 0.
|
||||
elif self.density_units == 'macro':
|
||||
density = self.density
|
||||
elif self.density_units == 'g/cc' or self.density_units == 'g/cm3':
|
||||
density = -self.density
|
||||
elif self.density_units == 'kg/m3':
|
||||
density = -0.001 * self.density
|
||||
elif self.density_units == 'atom/b-cm':
|
||||
density = self.density
|
||||
elif self.density_units == 'atom/cm3' or self.density_units == 'atom/cc':
|
||||
density = 1.E-24 * self.density
|
||||
|
||||
# For ease of processing split out nuc, nuc_density,
|
||||
# and nuc_density_type in to separate arrays
|
||||
nucs = []
|
||||
nuc_densities = []
|
||||
nuc_density_types = []
|
||||
for nuclide in nuclides.items():
|
||||
nuc, nuc_density, nuc_density_type = nuclide[1]
|
||||
nucs.append(nuc)
|
||||
nuc_densities.append(nuc_density)
|
||||
nuc_density_types.append(nuc_density_type)
|
||||
|
||||
if sum_density:
|
||||
density = np.sum(nuc_densities)
|
||||
percent_in_atom = np.all(nuc_density_types == 'ao')
|
||||
density_in_atom = density > 0.
|
||||
sum_percent = 0.
|
||||
|
||||
awrs = []
|
||||
for n, nuclide in enumerate(nuclides.items()):
|
||||
awr = openmc.data.atomic_mass(nuclide[0])
|
||||
if awr is not None:
|
||||
awrs.append(awr / openmc.data.NEUTRON_MASS)
|
||||
else:
|
||||
raise ValueError(nuclide[0] + " is invalid")
|
||||
|
||||
# Now that we have the awr, lets finish calculating densities
|
||||
sum_percent = np.sum(nuc_densities)
|
||||
nuc_densities = nuc_densities / sum_percent
|
||||
if not density_in_atom:
|
||||
sum_percent = 0.
|
||||
for n, nuc in enumerate(nucs):
|
||||
x = nuc_densities[n]
|
||||
sum_percent += x * awrs[n]
|
||||
sum_percent = 1. / sum_percent
|
||||
density = -density * sum_percent * \
|
||||
openmc.data.AVOGADRO / openmc.data.NEUTRON_MASS * 1.E-24
|
||||
nuc_densities = density * nuc_densities
|
||||
|
||||
nuclides = OrderedDict()
|
||||
for n, nuc in enumerate(nucs):
|
||||
nuclides[nuc] = (nuc, nuc_densities[n])
|
||||
|
||||
return nuclides
|
||||
|
||||
def _get_nuclide_xml(self, nuclide, distrib=False):
|
||||
xml_element = ET.Element("nuclide")
|
||||
xml_element.set("name", nuclide[0].name)
|
||||
|
|
|
|||
|
|
@ -1,10 +1,8 @@
|
|||
from numbers import Integral
|
||||
import sys
|
||||
import warnings
|
||||
|
||||
from six import string_types
|
||||
|
||||
from openmc.checkvalue import check_type
|
||||
import openmc.checkvalue as cv
|
||||
|
||||
|
||||
class Nuclide(object):
|
||||
|
|
@ -72,7 +70,7 @@ class Nuclide(object):
|
|||
|
||||
@name.setter
|
||||
def name(self, name):
|
||||
check_type('name', name, string_types)
|
||||
cv.check_type('name', name, string_types)
|
||||
self._name = name
|
||||
|
||||
if '-' in name:
|
||||
|
|
|
|||
542
openmc/plotter.py
Normal file
542
openmc/plotter.py
Normal file
|
|
@ -0,0 +1,542 @@
|
|||
from numbers import Integral, Real
|
||||
from six import string_types
|
||||
from itertools import chain
|
||||
|
||||
import numpy as np
|
||||
|
||||
import openmc.checkvalue as cv
|
||||
import openmc.data
|
||||
|
||||
# Supported keywords for xs plotting
|
||||
PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
|
||||
'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity',
|
||||
'slowing-down power', 'damage']
|
||||
|
||||
# Special MT values
|
||||
UNITY_MT = -1
|
||||
XI_MT = -2
|
||||
|
||||
# MTs to combine to generate associated plot_types
|
||||
_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt != 27]
|
||||
PLOT_TYPES_MT = {'total': openmc.data.SUM_RULES[1],
|
||||
'scatter': [2] + _INELASTIC,
|
||||
'elastic': [2],
|
||||
'inelastic': _INELASTIC,
|
||||
'fission': [18],
|
||||
'absorption': [27], 'capture': [101],
|
||||
'nu-fission': [18],
|
||||
'nu-scatter': [2] + _INELASTIC,
|
||||
'unity': [UNITY_MT],
|
||||
'slowing-down power': [2] + _INELASTIC + [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,),
|
||||
'scatter': (np.add,) * (len(PLOT_TYPES_MT['scatter']) - 1),
|
||||
'elastic': (),
|
||||
'inelastic': (np.add,) * (len(PLOT_TYPES_MT['inelastic']) - 1),
|
||||
'fission': (), 'absorption': (),
|
||||
'capture': (), 'nu-fission': (),
|
||||
'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,),
|
||||
'damage': ()}
|
||||
|
||||
# Types of plots to plot linearly in y
|
||||
PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter',
|
||||
'nu-fission / absorption', 'fission / absorption'}
|
||||
|
||||
|
||||
def plot_xs(this, types, divisor_types=None, temperature=294., axis=None,
|
||||
sab_name=None, cross_sections=None, enrichment=None, **kwargs):
|
||||
"""Creates a figure of continuous-energy cross sections for this item
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : openmc.Element, openmc.Nuclide, or openmc.Material
|
||||
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.
|
||||
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
|
||||
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 calculate_xs(this, types, temperature=294., sab_name=None,
|
||||
cross_sections=None, enrichment=None):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : openmc.Element, openmc.Nuclide, or openmc.Material
|
||||
Object to source data from
|
||||
types : Iterable of values of PLOT_TYPES
|
||||
The type of cross sections to calculate
|
||||
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.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable.
|
||||
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
|
||||
(natural composition).
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy_grid : numpy.ndarray
|
||||
Energies at which cross sections are calculated, in units of eV
|
||||
data : numpy.ndarray
|
||||
Cross sections calculated at the energy grid described by energy_grid
|
||||
|
||||
"""
|
||||
|
||||
# Check types
|
||||
cv.check_type('temperature', temperature, Real)
|
||||
if sab_name:
|
||||
cv.check_type('sab_name', sab_name, string_types)
|
||||
if enrichment:
|
||||
cv.check_type('enrichment', enrichment, Real)
|
||||
|
||||
if isinstance(this, openmc.Nuclide):
|
||||
energy_grid, xs = _calculate_xs_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.
|
||||
data = np.zeros((len(types), len(energy_grid)))
|
||||
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)
|
||||
elif isinstance(this, openmc.Material):
|
||||
energy_grid, data = _calculate_xs_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):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : openmc.Nuclide
|
||||
Nuclide object to source data from
|
||||
types : Iterable of str or Integral
|
||||
The type of cross sections to calculate; values can either be those
|
||||
in openmc.PLOT_TYPES or integers which correspond to reaction
|
||||
channel (MT) numbers.
|
||||
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.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy_grid : numpy.ndarray
|
||||
Energies at which cross sections are calculated, in units of eV
|
||||
data : Iterable of Callable
|
||||
Requested cross section functions
|
||||
|
||||
"""
|
||||
|
||||
# Parse the types
|
||||
mts = []
|
||||
ops = []
|
||||
yields = []
|
||||
for line in types:
|
||||
if line in PLOT_TYPES:
|
||||
mts.append(PLOT_TYPES_MT[line])
|
||||
if line.startswith('nu'):
|
||||
yields.append(True)
|
||||
else:
|
||||
yields.append(False)
|
||||
ops.append(PLOT_TYPES_OP[line])
|
||||
else:
|
||||
# Not a built-in type, we have to parse it ourselves
|
||||
cv.check_type('MT in types', line, Integral)
|
||||
cv.check_greater_than('MT in types', line, 0)
|
||||
mts.append((line,))
|
||||
ops.append(())
|
||||
yields.append(False)
|
||||
|
||||
# Load the library
|
||||
library = openmc.data.DataLibrary.from_xml(cross_sections)
|
||||
|
||||
# Convert temperature to format needed for access in the library
|
||||
strT = "{}K".format(int(round(temperature)))
|
||||
T = temperature
|
||||
|
||||
# Now we can create the data sets to be plotted
|
||||
energy_grid = []
|
||||
xs = []
|
||||
lib = library.get_by_material(this.name)
|
||||
if lib is not None:
|
||||
nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path'])
|
||||
# Obtain the nearest temperature
|
||||
if strT in nuc.temperatures:
|
||||
nucT = strT
|
||||
else:
|
||||
data_Ts = nuc.temperatures
|
||||
for t in range(len(data_Ts)):
|
||||
# Take off the "K" and convert to a float
|
||||
data_Ts[t] = float(data_Ts[t][:-1])
|
||||
min_delta = float('inf')
|
||||
closest_t = -1
|
||||
for t in data_Ts:
|
||||
if abs(data_Ts[t] - T) < min_delta:
|
||||
closest_t = t
|
||||
nucT = "{}K".format(int(round(data_Ts[closest_t])))
|
||||
|
||||
# Prep S(a,b) data if needed
|
||||
if sab_name:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
|
||||
# Obtain the nearest temperature
|
||||
if strT in sab.temperatures:
|
||||
sabT = strT
|
||||
else:
|
||||
data_Ts = sab.temperatures
|
||||
for t in range(len(data_Ts)):
|
||||
# Take off the "K" and convert to a float
|
||||
data_Ts[t] = float(data_Ts[t][:-1])
|
||||
min_delta = np.finfo(np.float64).max
|
||||
closest_t = -1
|
||||
for t in data_Ts:
|
||||
if abs(data_Ts[t] - T) < min_delta:
|
||||
closest_t = t
|
||||
sabT = "{}K".format(int(round(data_Ts[closest_t])))
|
||||
|
||||
# Create an energy grid composed the S(a,b) and
|
||||
# the nuclide's grid
|
||||
grid = nuc.energy[nucT]
|
||||
sab_Emax = 0.
|
||||
sab_funcs = []
|
||||
if sab.elastic_xs:
|
||||
elastic = sab.elastic_xs[sabT]
|
||||
if isinstance(elastic, openmc.data.CoherentElastic):
|
||||
grid = np.union1d(grid, elastic.bragg_edges)
|
||||
if elastic.bragg_edges[-1] > sab_Emax:
|
||||
sab_Emax = elastic.bragg_edges[-1]
|
||||
elif isinstance(elastic, openmc.data.Tabulated1D):
|
||||
grid = np.union1d(grid, elastic.x)
|
||||
if elastic.x[-1] > sab_Emax:
|
||||
sab_Emax = elastic.x[-1]
|
||||
sab_funcs.append(elastic)
|
||||
if sab.inelastic_xs:
|
||||
inelastic = sab.inelastic_xs[sabT]
|
||||
grid = np.union1d(grid, inelastic.x)
|
||||
if inelastic.x[-1] > sab_Emax:
|
||||
sab_Emax = inelastic.x[-1]
|
||||
sab_funcs.append(inelastic)
|
||||
energy_grid = grid
|
||||
else:
|
||||
energy_grid = nuc.energy[nucT]
|
||||
|
||||
for i, mt_set in enumerate(mts):
|
||||
# Get the reaction xs data from the nuclide
|
||||
funcs = []
|
||||
op = ops[i]
|
||||
for mt in mt_set:
|
||||
if mt == 2:
|
||||
if sab_name:
|
||||
# Then we need to do a piece-wise function of
|
||||
# The S(a,b) and non-thermal data
|
||||
sab_sum = openmc.data.Sum(sab_funcs)
|
||||
pw_funcs = openmc.data.Regions1D(
|
||||
[sab_sum, nuc[mt].xs[nucT]],
|
||||
[sab_Emax])
|
||||
funcs.append(pw_funcs)
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt in nuc:
|
||||
if yields[i]:
|
||||
# Get the total yield first if available. This will be
|
||||
# used primarily for fission.
|
||||
for prod in chain(nuc[mt].products,
|
||||
nuc[mt].derived_products):
|
||||
if prod.particle == 'neutron' and \
|
||||
prod.emission_mode == 'total':
|
||||
func = openmc.data.Combination(
|
||||
[nuc[mt].xs[nucT], prod.yield_],
|
||||
[np.multiply])
|
||||
funcs.append(func)
|
||||
break
|
||||
else:
|
||||
# Total doesn't exist so we have to create from
|
||||
# prompt and delayed. This is used for scatter
|
||||
# multiplication.
|
||||
func = None
|
||||
for prod in chain(nuc[mt].products,
|
||||
nuc[mt].derived_products):
|
||||
if prod.particle == 'neutron' and \
|
||||
prod.emission_mode != 'total':
|
||||
if func:
|
||||
func = openmc.data.Combination(
|
||||
[prod.yield_, func], [np.add])
|
||||
else:
|
||||
func = prod.yield_
|
||||
if func:
|
||||
funcs.append(openmc.data.Combination(
|
||||
[func, nuc[mt].xs[nucT]], [np.multiply]))
|
||||
else:
|
||||
# If func is still None, then there were no
|
||||
# products. In that case, assume the yield is
|
||||
# one as its not provided for some summed
|
||||
# reactions like MT=4
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt == UNITY_MT:
|
||||
funcs.append(lambda x: 1.)
|
||||
elif mt == XI_MT:
|
||||
awr = nuc.atomic_weight_ratio
|
||||
alpha = ((awr - 1.) / (awr + 1.))**2
|
||||
xi = 1. + alpha * np.log(alpha) / (1. - alpha)
|
||||
funcs.append(lambda x: xi)
|
||||
else:
|
||||
funcs.append(lambda x: 0.)
|
||||
xs.append(openmc.data.Combination(funcs, op))
|
||||
else:
|
||||
raise ValueError(this.name + " not in library")
|
||||
|
||||
return energy_grid, xs
|
||||
|
||||
|
||||
def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None,
|
||||
sab_name=None, enrichment=None):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
this : {openmc.Material, openmc.Element}
|
||||
Object to source data from
|
||||
types : Iterable of values of PLOT_TYPES
|
||||
The type of cross sections to calculate
|
||||
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.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
sab_name : str, optional
|
||||
Name of S(a,b) library to apply to MT=2 data when applicable.
|
||||
enrichment : float, optional
|
||||
Enrichment for U235 in weight percent. For example, input 4.95 for
|
||||
4.95 weight percent enriched U. Default is None
|
||||
(natural composition).
|
||||
|
||||
Returns
|
||||
-------
|
||||
energy_grid : numpy.ndarray
|
||||
Energies at which cross sections are calculated, in units of eV
|
||||
data : numpy.ndarray
|
||||
Cross sections calculated at the energy grid described by energy_grid
|
||||
|
||||
"""
|
||||
|
||||
if isinstance(this, openmc.Material):
|
||||
if this.temperature is not None:
|
||||
T = this.temperature
|
||||
else:
|
||||
T = temperature
|
||||
else:
|
||||
T = temperature
|
||||
|
||||
# Load the library
|
||||
library = openmc.data.DataLibrary.from_xml(cross_sections)
|
||||
|
||||
if isinstance(this, openmc.Material):
|
||||
# Expand elements in to nuclides with atomic densities
|
||||
nuclides = this.get_nuclide_atom_densities()
|
||||
# For ease of processing split out the nuclide and its fraction
|
||||
nuc_fractions = {nuclide[1][0].name: nuclide[1][1]
|
||||
for nuclide in nuclides.items()}
|
||||
# Create a dict of [nuclide name] = nuclide object to carry forward
|
||||
# with a common nuclides format between openmc.Material and
|
||||
# openmc.Element objects
|
||||
nuclides = {nuclide[1][0].name: nuclide[1][0]
|
||||
for nuclide in nuclides.items()}
|
||||
else:
|
||||
# Expand elements in to nuclides with atomic densities
|
||||
nuclides = this.expand(1., 'ao', enrichment=enrichment,
|
||||
cross_sections=cross_sections)
|
||||
# For ease of processing split out the nuclide and its fraction
|
||||
nuc_fractions = {nuclide[0].name: nuclide[1] for nuclide in nuclides}
|
||||
# Create a dict of [nuclide name] = nuclide object to carry forward
|
||||
# with a common nuclides format between openmc.Material and
|
||||
# openmc.Element objects
|
||||
nuclides = {nuclide[0].name: nuclide[0] for nuclide in nuclides}
|
||||
|
||||
# Identify the nuclides which have S(a,b) data
|
||||
sabs = {}
|
||||
for nuclide in nuclides.items():
|
||||
sabs[nuclide[0]] = None
|
||||
if isinstance(this, openmc.Material):
|
||||
for sab_name in this._sab:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(
|
||||
library.get_by_material(sab_name)['path'])
|
||||
for nuc in sab.nuclides:
|
||||
sabs[nuc] = library.get_by_material(sab_name)['path']
|
||||
else:
|
||||
if sab_name:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
|
||||
for nuc in sab.nuclides:
|
||||
sabs[nuc] = library.get_by_material(sab_name)['path']
|
||||
|
||||
# Now we can create the data sets to be plotted
|
||||
xs = {}
|
||||
E = []
|
||||
for nuclide in nuclides.items():
|
||||
name = nuclide[0]
|
||||
nuc = nuclide[1]
|
||||
sab_tab = sabs[name]
|
||||
temp_E, temp_xs = calculate_xs(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.
|
||||
xs[name] = [openmc.data.Tabulated1D(temp_E, temp_xs[line])
|
||||
for line in range(len(types))]
|
||||
|
||||
# Condense the data for every nuclide
|
||||
# First create a union energy grid
|
||||
energy_grid = E[0]
|
||||
for grid in E[1:]:
|
||||
energy_grid = np.union1d(energy_grid, grid)
|
||||
|
||||
# Now we can combine all the nuclidic data
|
||||
data = np.zeros((len(types), len(energy_grid)))
|
||||
for line in range(len(types)):
|
||||
if types[line] == 'unity':
|
||||
data[line, :] = 1.
|
||||
else:
|
||||
for nuclide in nuclides.items():
|
||||
name = nuclide[0]
|
||||
data[line, :] += (nuc_fractions[name] *
|
||||
xs[name][line](energy_grid))
|
||||
|
||||
return energy_grid, data
|
||||
Loading…
Add table
Add a link
Reference in a new issue