mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Extended plots to elements and nuclides
This commit is contained in:
Adam Nelson
parent
8de82dcbd3
commit
76991a223f
5 changed files with 597 additions and 209 deletions
|
|
@ -24,6 +24,7 @@ from openmc.statepoint import *
|
|||
from openmc.summary import *
|
||||
from openmc.particle_restart import *
|
||||
from openmc.mixin import *
|
||||
from openmc.plot_data import *
|
||||
|
||||
try:
|
||||
from openmc.opencg_compatible import *
|
||||
|
|
|
|||
|
|
@ -5,10 +5,12 @@ import os
|
|||
|
||||
from six import string_types
|
||||
from xml.etree import ElementTree as ET
|
||||
import numpy as np
|
||||
|
||||
import openmc
|
||||
from openmc.checkvalue import check_type, check_length
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.data import NATURAL_ABUNDANCE, atomic_mass
|
||||
from openmc.plot_data import *
|
||||
|
||||
|
||||
class Element(object):
|
||||
|
|
@ -80,8 +82,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
|
||||
|
|
@ -212,7 +214,7 @@ class Element(object):
|
|||
# its natural nuclides
|
||||
else:
|
||||
for nuclide in natural_nuclides:
|
||||
abundances[nuclide] = NATURAL_ABUNDNACE[nuclide]
|
||||
abundances[nuclide] = NATURAL_ABUNDANCE[nuclide]
|
||||
|
||||
# Modify mole fractions if enrichment provided
|
||||
if enrichment is not None:
|
||||
|
|
@ -257,3 +259,202 @@ class Element(object):
|
|||
isotopes.append((nuc, percent*abundance, percent_type))
|
||||
|
||||
return isotopes
|
||||
|
||||
def plot_xs(self, types, divisor_types=None, temperature=294.,
|
||||
Erange=(1.E-5, 20.E6), enrichment=None, cross_sections=None):
|
||||
"""Creates a figure of continuous-energy microscopic cross sections
|
||||
for this element
|
||||
|
||||
Parameters
|
||||
----------
|
||||
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.
|
||||
Erange : tuple of floats
|
||||
Energy range (in eV) to plot the cross section within
|
||||
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).
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
fig : matplotlib.figure.Figure
|
||||
Matplotlib Figure of the generated macroscopic cross section
|
||||
|
||||
"""
|
||||
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
E, data = self.calculate_xs(types, temperature, cross_sections)
|
||||
|
||||
if divisor_types:
|
||||
cv.check_length('divisor types', divisor_types, len(types),
|
||||
len(types))
|
||||
Ediv, data_div = self.calculate_xs(divisor_types, temperature,
|
||||
cross_sections)
|
||||
|
||||
# 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 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]
|
||||
data = data_new
|
||||
|
||||
# Generate the plot
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(111)
|
||||
iE_max = np.searchsorted(E, Erange[1])
|
||||
min_data = np.finfo(np.float64).max
|
||||
max_data = np.finfo(np.float64).min
|
||||
for i in range(len(data)):
|
||||
if np.sum(data[i, :]) > 0.:
|
||||
ax.loglog(E, data[i, :], label=types[i])
|
||||
min_data = min(min_data, np.min(data[i, :iE_max]))
|
||||
max_data = max(max_data, np.max(data[i, :iE_max]))
|
||||
|
||||
ax.set_xlabel('Energy [eV]')
|
||||
ax.set_ylabel('Elemental Cross Section [1/cm]')
|
||||
ax.legend(loc='best')
|
||||
ax.set_xlim(Erange)
|
||||
ax.set_ylim(min_data, max_data)
|
||||
if self.name is not None:
|
||||
title = 'Cross Section for ' + self.name
|
||||
ax.set_title(title)
|
||||
|
||||
return fig
|
||||
|
||||
def calculate_xs(self, types, temperature=294., enrichment=None,
|
||||
cross_sections=None):
|
||||
"""Calculates continuous-energy macroscopic cross sections of a
|
||||
requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
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.
|
||||
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).
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
unionE : 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
|
||||
|
||||
"""
|
||||
|
||||
# Check types
|
||||
if cross_sections is not None:
|
||||
cv.check_type('cross_sections', cross_sections, str)
|
||||
cv.check_iterable_type('types', types, str)
|
||||
|
||||
# Parse the types
|
||||
mts = []
|
||||
ops = []
|
||||
yields = []
|
||||
for line in types:
|
||||
if line in PLOT_TYPES:
|
||||
mts.append(PLOT_TYPES_MT[line])
|
||||
yields.append(PLOT_TYPES_YIELD[line])
|
||||
ops.append(PLOT_TYPES_OP[line])
|
||||
else:
|
||||
# Not a built-in type, we have to parse it ourselves
|
||||
raise NotImplementedError()
|
||||
|
||||
cv.check_type('temperature', temperature, Real)
|
||||
|
||||
# Expand elements in to nuclides with atomic densities
|
||||
nuclides = self.expand(100., 'ao', enrichment=enrichment,
|
||||
cross_sections=cross_sections)
|
||||
|
||||
# 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")
|
||||
|
||||
# For ease of processing split out nuc and nuc_density
|
||||
nucs = []
|
||||
nuc_fractions = []
|
||||
for nuclide in nuclides.items():
|
||||
nuc_name, nuc_data = nuclide
|
||||
nuc, nuc_fraction, nuc_fraction_type = nuc_data
|
||||
nucs.append(nuc)
|
||||
nuc_fractions.append(nuc_fraction)
|
||||
|
||||
# Identify the nuclides which need S(a,b) data
|
||||
sabs = {}
|
||||
for nuc in nucs:
|
||||
sabs[nuc.name] = None
|
||||
|
||||
if sab_name:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(sab_name)
|
||||
for nuc in sab.nuclides:
|
||||
sabs[nuc] = library.get_by_materials(sab_name)['path']
|
||||
|
||||
# Now we can create the data sets to be plotted
|
||||
xs = []
|
||||
E = []
|
||||
n = -1
|
||||
for nuclide in nuclides.items():
|
||||
n += 1
|
||||
# import pdb; pdb.set_trace()
|
||||
nuc_obj = openmc.Nuclide(nuclide[0].name)
|
||||
sab_tab = sabs[nucs[n].name]
|
||||
temp_E, temp_xs = nuc_obj.calculate_xs(types, temperature, sab_tab,
|
||||
cross_sections)
|
||||
E.append(temp_E)
|
||||
xs.append(temp_xs)
|
||||
|
||||
# Condense the data for every nuclide
|
||||
# First create a union energy grid
|
||||
unionE = E[0]
|
||||
for n in range(1, len(E)):
|
||||
unionE = np.union1d(unionE, E[n])
|
||||
|
||||
# Now we can combine all the nuclidic data
|
||||
data = np.zeros((len(mts), len(unionE)))
|
||||
for l in range(len(mts)):
|
||||
if types[l] == 'unity':
|
||||
data[l, :] = 1.
|
||||
else:
|
||||
for n in range(len(nuclides)):
|
||||
data[l, :] += nuc_fractions[n] * xs[n][l](unionE)
|
||||
|
||||
return unionE, data
|
||||
|
|
|
|||
|
|
@ -3,6 +3,7 @@ from copy import deepcopy
|
|||
from numbers import Real, Integral
|
||||
import warnings
|
||||
from xml.etree import ElementTree as ET
|
||||
import os
|
||||
|
||||
from six import string_types
|
||||
import numpy as np
|
||||
|
|
@ -11,6 +12,7 @@ import openmc
|
|||
import openmc.data
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.clean_xml import sort_xml_elements, clean_xml_indentation
|
||||
from openmc.plot_data import *
|
||||
|
||||
|
||||
# A static variable for auto-generated Material IDs
|
||||
|
|
@ -27,62 +29,6 @@ def reset_auto_material_id():
|
|||
DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum',
|
||||
'macro']
|
||||
|
||||
# Supported keywords for material xs plotting
|
||||
_PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
|
||||
'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity']
|
||||
|
||||
# MTs to combine to generate associated plot_types
|
||||
_PLOT_TYPES_MT = {'total': (1,), 'scatter': (1, 27), 'elastic': (2,),
|
||||
'inelastic': (1, 27, 2), 'fission': (18,),
|
||||
'absorption': (27,), 'capture': (101,),
|
||||
'nu-fission': (18,),
|
||||
'nu-scatter': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, 29, 30,
|
||||
32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152,
|
||||
153, 154, 156, 157, 158, 159, 160, 161, 162,
|
||||
163, 164, 165, 166, 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),
|
||||
'unity': (0,)}
|
||||
# Operations to use when combining MTs the first np.add is used in reference
|
||||
# to zero
|
||||
_PLOT_TYPES_OP = {'total': (), 'scatter': (np.subtract,), 'elastic': (),
|
||||
'inelastic': (np.subtract, np.subtract),
|
||||
'fission': (), 'absorption': (),
|
||||
'capture': (), 'nu-fission': (),
|
||||
'nu-scatter': (np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add),
|
||||
'unity': ()}
|
||||
# Whether or not to multiply the reaction by the yield as well
|
||||
_PLOT_TYPES_YIELD = {'total': (False,), 'scatter': (False, False),
|
||||
'elastic': (False,), 'inelastic': (False, False, False),
|
||||
'fission': (False,), 'absorption': (False,),
|
||||
'capture': (False,), 'nu-fission': (True,),
|
||||
'nu-scatter': (True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True),
|
||||
'unity': (False,)}
|
||||
|
||||
|
||||
class Material(object):
|
||||
"""A material composed of a collection of nuclides/elements that can be
|
||||
|
|
@ -637,14 +583,14 @@ class Material(object):
|
|||
|
||||
return nuclides
|
||||
|
||||
def get_nuclide_atom_densities(self, library):
|
||||
def get_nuclide_atom_densities(self, cross_sections=None):
|
||||
"""Returns all nuclides in the material and their atomic densities in
|
||||
units of atom/b-cm
|
||||
|
||||
Parameters
|
||||
----------
|
||||
library : openmc.data.DataLibrary
|
||||
Library of data to use for plotting.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -656,7 +602,20 @@ class Material(object):
|
|||
|
||||
import scipy.constants as sc
|
||||
|
||||
cv.check_type('library', library, openmc.data.DataLibrary)
|
||||
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")
|
||||
|
||||
# Expand elements in to nuclides
|
||||
nuclides = self.get_nuclide_densities()
|
||||
|
|
@ -728,28 +687,28 @@ class Material(object):
|
|||
|
||||
return nuclides
|
||||
|
||||
def plot_xs(self, library, types, divisor_types=None, temperature=294.,
|
||||
Erange=(1.E-5, 20.E6)):
|
||||
def plot_xs(self, types, divisor_types=None, temperature=294.,
|
||||
Erange=(1.E-5, 20.E6), cross_sections=None):
|
||||
"""Creates a figure of continuous-energy macroscopic cross sections
|
||||
for this material
|
||||
|
||||
Parameters
|
||||
----------
|
||||
library : openmc.data.DataLibrary
|
||||
Library of data to use for plotting.
|
||||
types : Iterable of values of _PLOT_TYPES
|
||||
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.
|
||||
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.
|
||||
Erange : tuple of floats
|
||||
Erange : 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.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -760,20 +719,20 @@ class Material(object):
|
|||
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
E, data = self.calculate_xs(library, types, temperature)
|
||||
E, data = self.calculate_xs(types, temperature, cross_sections)
|
||||
|
||||
if divisor_types:
|
||||
cv.check_length('divisor types', divisor_types, len(types),
|
||||
len(types))
|
||||
Ediv, data_div = self.calculate_xs(library, divisor_types,
|
||||
temperature)
|
||||
Ediv, data_div = self.calculate_xs(divisor_types, temperature,
|
||||
cross_sections)
|
||||
|
||||
# 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)))
|
||||
# import pdb; pdb.set_trace()
|
||||
|
||||
for l in range(len(types)):
|
||||
data_new[l, :] = \
|
||||
np.divide(np.interp(E, Enum, data[l, :]),
|
||||
|
|
@ -805,21 +764,21 @@ class Material(object):
|
|||
|
||||
return fig
|
||||
|
||||
def calculate_xs(self, library, types, temperature=294.):
|
||||
def calculate_xs(self, types, temperature=294., cross_sections=None):
|
||||
"""Calculates continuous-energy macroscopic cross sections of a
|
||||
requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
library : openmc.data.DataLibrary
|
||||
Library of data to use for plotting.
|
||||
types : Iterable of values of _PLOT_TYPES
|
||||
The type of cross sections to include in the plot.
|
||||
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.
|
||||
|
||||
Returns
|
||||
-------
|
||||
|
|
@ -831,10 +790,9 @@ class Material(object):
|
|||
|
||||
"""
|
||||
|
||||
import scipy.constants as sc
|
||||
|
||||
# Check types
|
||||
cv.check_type('library', library, openmc.data.DataLibrary)
|
||||
if cross_sections is not None:
|
||||
cv.check_type('cross_sections', cross_sections, str)
|
||||
cv.check_iterable_type('types', types, str)
|
||||
|
||||
# Parse the types
|
||||
|
|
@ -842,25 +800,35 @@ class Material(object):
|
|||
ops = []
|
||||
yields = []
|
||||
for line in types:
|
||||
if line in _PLOT_TYPES:
|
||||
mts.append(_PLOT_TYPES_MT[line])
|
||||
yields.append(_PLOT_TYPES_YIELD[line])
|
||||
ops.append(_PLOT_TYPES_OP[line])
|
||||
if line in PLOT_TYPES:
|
||||
mts.append(PLOT_TYPES_MT[line])
|
||||
yields.append(PLOT_TYPES_YIELD[line])
|
||||
ops.append(PLOT_TYPES_OP[line])
|
||||
else:
|
||||
# Not a built-in type, we have to parse it ourselves
|
||||
raise NotImplementedError()
|
||||
|
||||
# Convert temperature to format needed for access in the library
|
||||
if self.temperature is not None:
|
||||
strT = "{}K".format(int(round(self.temperature)))
|
||||
T = self.temperature
|
||||
else:
|
||||
cv.check_type('temperature', temperature, Real)
|
||||
strT = "{}K".format(int(round(temperature)))
|
||||
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")
|
||||
|
||||
# Expand elements in to nuclides with atomic densities
|
||||
nuclides = self.get_nuclide_atom_densities(library)
|
||||
nuclides = self.get_nuclide_atom_densities(cross_sections)
|
||||
|
||||
# For ease of processing split out nuc and nuc_density
|
||||
nucs = []
|
||||
|
|
@ -888,118 +856,13 @@ class Material(object):
|
|||
n = -1
|
||||
for nuclide in nuclides.items():
|
||||
n += 1
|
||||
lib = library.get_by_materials(nuclide[0])
|
||||
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 = np.finfo(np.float64).max
|
||||
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])))
|
||||
|
||||
# Create an energy grid composed of either the S(a,b) and
|
||||
# nuclide's grid, or just the nuclide's grid, depending on if
|
||||
# the S(a,b) data is available for this nuclide
|
||||
sab_tab = sabs[nucs[n].name]
|
||||
if sab_tab:
|
||||
sab = openmc.data.ThermalScattering.from_hdf5(sab_tab)
|
||||
# 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])))
|
||||
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)
|
||||
E.append(grid)
|
||||
else:
|
||||
E.append(nuc.energy[nucT])
|
||||
xs.append([])
|
||||
for i, mt_set in enumerate(mts):
|
||||
# Get the reaction xs data from the nuclide
|
||||
funcs = []
|
||||
op = ops[i]
|
||||
for mt, yield_check in zip(mt_set, yields[i]):
|
||||
if mt == 2:
|
||||
if sab_tab:
|
||||
# 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.Piecewise(
|
||||
[sab_sum, nuc[mt].xs[nucT]],
|
||||
[sab_Emax])
|
||||
funcs.append(pw_funcs)
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt in nuc:
|
||||
if yield_check:
|
||||
found_it = False
|
||||
for prod in nuc[mt].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)
|
||||
found_it = True
|
||||
break
|
||||
if not found_it:
|
||||
for prod in nuc[mt].products:
|
||||
if prod.particle == 'neutron' and \
|
||||
prod.emission_mode == 'prompt':
|
||||
func = openmc.data.Combination(
|
||||
[nuc[mt].xs[nucT],
|
||||
prod.yield_], [np.multiply])
|
||||
funcs.append(func)
|
||||
found_it = True
|
||||
break
|
||||
if not found_it:
|
||||
# Assume the yield is 1
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt == 0:
|
||||
funcs.append(lambda x: 1.)
|
||||
else:
|
||||
funcs.append(lambda x: 0.)
|
||||
xs[-1].append(openmc.data.Combination(funcs, op))
|
||||
else:
|
||||
raise ValueError(nuclide[0] + " not in library")
|
||||
# import pdb; pdb.set_trace()
|
||||
nuc_obj = openmc.Nuclide(nuclide[0].name)
|
||||
sab_tab = sabs[nucs[n].name]
|
||||
temp_E, temp_xs = nuc_obj.calculate_xs(types, T, sab_tab,
|
||||
cross_sections)
|
||||
E.append(temp_E)
|
||||
xs.append(temp_xs)
|
||||
|
||||
# Condense the data for every nuclide
|
||||
# First create a union energy grid
|
||||
|
|
|
|||
|
|
@ -1,10 +1,13 @@
|
|||
from numbers import Integral
|
||||
from numbers import Integral, Real
|
||||
import sys
|
||||
import warnings
|
||||
import os
|
||||
|
||||
from six import string_types
|
||||
|
||||
from openmc.checkvalue import check_type
|
||||
import openmc.checkvalue as cv
|
||||
from openmc.plot_data import *
|
||||
import openmc.data
|
||||
|
||||
|
||||
class Nuclide(object):
|
||||
|
|
@ -72,7 +75,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:
|
||||
|
|
@ -93,3 +96,266 @@ class Nuclide(object):
|
|||
raise ValueError(msg)
|
||||
|
||||
self._scattering = scattering
|
||||
|
||||
def plot_xs(self, types, divisor_types=None, temperature=294.,
|
||||
Erange=(1.E-5, 20.E6), sab_name=None, cross_sections=None):
|
||||
"""Creates a figure of continuous-energy cross sections for this
|
||||
nuclide
|
||||
|
||||
Parameters
|
||||
----------
|
||||
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.
|
||||
Erange : 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.
|
||||
cross_sections : str, optional
|
||||
Location of cross_sections.xml file. Default is None.
|
||||
|
||||
Returns
|
||||
-------
|
||||
fig : matplotlib.figure.Figure
|
||||
Matplotlib Figure of the generated macroscopic cross section
|
||||
|
||||
"""
|
||||
|
||||
from matplotlib import pyplot as plt
|
||||
|
||||
E, data = self.calculate_xs(types, temperature, sab_name,
|
||||
cross_sections)
|
||||
|
||||
if divisor_types:
|
||||
cv.check_length('divisor types', divisor_types, len(types),
|
||||
len(types))
|
||||
Ediv, data_div = self.calculate_xs(divisor_types, temperature,
|
||||
sab_name, cross_sections)
|
||||
|
||||
# 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 = []
|
||||
|
||||
for l in range(len(types)):
|
||||
data_new.append(openmc.data.Combination([data[l], data_div[l]],
|
||||
[np.divide]))
|
||||
if divisor_types[l] != 'unity':
|
||||
types[l] = types[l] + ' / ' + divisor_types[l]
|
||||
data = data_new
|
||||
|
||||
# Generate the plot
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(111)
|
||||
iE_max = np.searchsorted(E, Erange[1])
|
||||
min_data = np.finfo(np.float64).max
|
||||
max_data = np.finfo(np.float64).min
|
||||
for i in range(len(data)):
|
||||
to_plot = data[i](E)
|
||||
if np.sum(to_plot) > 0.:
|
||||
ax.loglog(E, to_plot, label=types[i])
|
||||
min_data = min(min_data, np.min(to_plot[:iE_max]))
|
||||
max_data = max(max_data, np.max(to_plot[:iE_max]))
|
||||
|
||||
ax.set_xlabel('Energy [eV]')
|
||||
ax.set_ylabel('Microscopic Cross Section [b]')
|
||||
ax.legend(loc='best')
|
||||
ax.set_xlim(Erange)
|
||||
ax.set_ylim(min_data, max_data)
|
||||
if self.name is not None:
|
||||
title = 'Microscopic Cross Section for ' + self.name
|
||||
ax.set_title(title)
|
||||
|
||||
return fig
|
||||
|
||||
def calculate_xs(self, types, temperature=294., sab_name=None,
|
||||
cross_sections=None):
|
||||
"""Calculates continuous-energy cross sections of a requested type
|
||||
|
||||
Parameters
|
||||
----------
|
||||
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.
|
||||
|
||||
Returns
|
||||
-------
|
||||
E : 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
|
||||
|
||||
"""
|
||||
|
||||
# Check types
|
||||
if cross_sections is not None:
|
||||
cv.check_type('cross_sections', cross_sections, str)
|
||||
cv.check_iterable_type('types', types, str)
|
||||
if sab_name:
|
||||
cv.check_type('sab_name', sab_name, str)
|
||||
|
||||
# Parse the types
|
||||
mts = []
|
||||
ops = []
|
||||
yields = []
|
||||
for line in types:
|
||||
if line in openmc.plot_data.PLOT_TYPES:
|
||||
mts.append(PLOT_TYPES_MT[line])
|
||||
yields.append(PLOT_TYPES_YIELD[line])
|
||||
ops.append(PLOT_TYPES_OP[line])
|
||||
else:
|
||||
# Not a built-in type, we have to parse it ourselves
|
||||
raise NotImplementedError()
|
||||
|
||||
# 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")
|
||||
|
||||
# Convert temperature to format needed for access in the library
|
||||
cv.check_type('temperature', temperature, Real)
|
||||
strT = "{}K".format(int(round(temperature)))
|
||||
T = temperature
|
||||
|
||||
# Now we can create the data sets to be plotted
|
||||
E = []
|
||||
xs = []
|
||||
lib = library.get_by_materials(self.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 = np.finfo(np.float64).max
|
||||
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)
|
||||
E = grid
|
||||
else:
|
||||
E = nuc.energy[nucT]
|
||||
|
||||
for i, mt_set in enumerate(mts):
|
||||
# Get the reaction xs data from the nuclide
|
||||
funcs = []
|
||||
op = ops[i]
|
||||
for mt, yield_check in zip(mt_set, yields[i]):
|
||||
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 yield_check:
|
||||
found_it = False
|
||||
for prod in nuc[mt].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)
|
||||
found_it = True
|
||||
break
|
||||
if not found_it:
|
||||
for prod in nuc[mt].products:
|
||||
if prod.particle == 'neutron' and \
|
||||
prod.emission_mode == 'prompt':
|
||||
func = openmc.data.Combination(
|
||||
[nuc[mt].xs[nucT],
|
||||
prod.yield_], [np.multiply])
|
||||
funcs.append(func)
|
||||
found_it = True
|
||||
break
|
||||
if not found_it:
|
||||
# Assume the yield is 1
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
else:
|
||||
funcs.append(nuc[mt].xs[nucT])
|
||||
elif mt == 0:
|
||||
funcs.append(lambda x: 1.)
|
||||
else:
|
||||
funcs.append(lambda x: 0.)
|
||||
xs.append(openmc.data.Combination(funcs, op))
|
||||
else:
|
||||
raise ValueError(nuclide[0] + " not in library")
|
||||
|
||||
return E, xs
|
||||
|
|
|
|||
57
openmc/plot_data.py
Normal file
57
openmc/plot_data.py
Normal file
|
|
@ -0,0 +1,57 @@
|
|||
import numpy as np
|
||||
|
||||
# Supported keywords for material xs plotting
|
||||
PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission',
|
||||
'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity']
|
||||
|
||||
# MTs to combine to generate associated plot_types
|
||||
PLOT_TYPES_MT = {'total': (2, 3,), 'scatter': (1, 27), 'elastic': (2,),
|
||||
'inelastic': (1, 27, 2), 'fission': (18,),
|
||||
'absorption': (27,), 'capture': (101,),
|
||||
'nu-fission': (18,),
|
||||
'nu-scatter': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, 29, 30,
|
||||
32, 33, 34, 35, 36, 37, 41, 42, 44, 45, 152,
|
||||
153, 154, 156, 157, 158, 159, 160, 161, 162,
|
||||
163, 164, 165, 166, 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),
|
||||
'unity': (0,)}
|
||||
# Operations to use when combining MTs the first np.add is used in reference
|
||||
# to zero
|
||||
PLOT_TYPES_OP = {'total': (np.add,), 'scatter': (np.subtract,), 'elastic': (),
|
||||
'inelastic': (np.subtract, np.subtract),
|
||||
'fission': (), 'absorption': (),
|
||||
'capture': (), 'nu-fission': (),
|
||||
'nu-scatter': (np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add,
|
||||
np.add, np.add, np.add, np.add, np.add),
|
||||
'unity': ()}
|
||||
# Whether or not to multiply the reaction by the yield as well
|
||||
PLOT_TYPES_YIELD = {'total': (False, False), 'scatter': (False, False),
|
||||
'elastic': (False,), 'inelastic': (False, False, False),
|
||||
'fission': (False,), 'absorption': (False,),
|
||||
'capture': (False,), 'nu-fission': (True,),
|
||||
'nu-scatter': (True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True, True, True, True, True,
|
||||
True),
|
||||
'unity': (False,)}
|
||||
Loading…
Add table
Add a link
Reference in a new issue