From 4e50b6e8c7330fbea96438cf35e807b33f70b705 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 5 Nov 2016 16:55:34 -0400 Subject: [PATCH 01/31] Ability to read in a cross_sections.xml from python and initial capability for plotting macroscopic xs --- openmc/data/library.py | 36 ++++++++++++++ openmc/material.py | 109 +++++++++++++++++++++++++++++++++++++++++ 2 files changed, 145 insertions(+) diff --git a/openmc/data/library.py b/openmc/data/library.py index 0485af064..3e01c2464 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -22,6 +22,12 @@ class DataLibrary(EqualityMixin): def __init__(self): self.libraries = [] + def __getitem__(self, material): + for library in self.libraries: + if material in library['materials']: + return library + return None + def register_file(self, filename): """Register a file with the data library. @@ -78,3 +84,33 @@ class DataLibrary(EqualityMixin): tree = ET.ElementTree(root) tree.write(path, xml_declaration=True, encoding='utf-8', method='xml') + + @classmethod + def from_xml(cls, path): + """Read cross section data library from an XML file. + + Parameters + ---------- + path : str + Path to XML file to read. + + """ + + data = cls() + + tree = ET.parse(path) + root = tree.getroot() + if root.find('directory') is not None: + directory = root.find('directory').text + else: + directory = os.path.dirname(path) + + for lib_element in root.findall('library'): + filename = os.path.join(directory, lib_element.attrib['path']) + filetype = lib_element.attrib['type'] + materials = lib_element.attrib['materials'].split() + library = {'path': filename, 'type': filetype, + 'materials': materials} + data.libraries.append(library) + + return data diff --git a/openmc/material.py b/openmc/material.py index c265c9069..271acfd8f 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -6,6 +6,9 @@ from xml.etree import ElementTree as ET import sys from six import string_types +from scipy.interpolate import interp1d +import numpy as np +import h5py import openmc import openmc.data @@ -27,6 +30,13 @@ 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'] +# MTs to sum to generate associated plot_types +_PLOT_TYPES_MT = {'total': (1,), 'scatter': (2, 4), 'elastic': (2,), + 'inelastic': (1,), 'fission': (18,), 'absorption': (102,)} + class Material(object): """A material composed of a collection of nuclides/elements that can be @@ -581,6 +591,105 @@ class Material(object): return nuclides + def plot_xs(self, library, types, labels=None): + """Creates a figure of macroscopic cross sections for this material + + Parameters + ---------- + library : openmc.data.DataLibrary + Library of data to use for plotting. + types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption'} or list thereof + The type of cross sections to include in the plot. This can either + be an MT number, a tuple of MT numbers (indicating they are to be + summed before plotting) or a string describing a common type. + These values can be either a single value, or an iterable + in the case of multiple sets per plot. + labels : str or list of str, optional + Labels to use in the plot legend for each type requested. + Single string (if there is only one type) or a list of strings + with a length matching the length of types. + + Returns + ------- + fig : matplotlib.figure.Figure + Matplotlib Figure of the generated macroscopic cross section + + """ + + # Check types + # ## Need to add error messages in some ELSEs + # ## Need to add in label handling + if isinstance(types, Integral): + cv.check_greater_than('types', types, 0) + elif isinstance(types, list): + for line in types: + if isinstance(line, Integral): + cv.check_greater_than('line in types', line, 0) + elif isinstance(line, tuple): + for entry in line: + if isinstance(entry, Integral): + cv.check_greater_than('entry in line in types', + entry, 0) + elif line in _PLOT_TYPES: + # Replace with the MTs to sum to simplify downstream code + line = _PLOT_TYPES_MT[line] + + # Convert temperature to format for accessing in the library + if self.temperature is not None: + T = "{}K".format(int(round(self.temperature))) + else: + T = None + + if isinstance(types, (Integral, str)): + types_ = [(types,)] + elif isinstance(types, tuple): + types_ = types + + # Now we can create the data sets to be plotted + xs = [] + E = [] + percent = [] + + for n, nuclide in enumerate(self.nuclides): + nuc, pct, units = nuclide + # import pdb; pdb.set_trace() + percent.append(pct) + lib = library[nuc] + if lib is not None: + nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) + # ## Need to check temperatures + nucT = T + # ## + E.append(nuc.energy[nucT]) + xs.append([]) + for l, line in enumerate(types_): + # Get the reaction from the nuclide + for mt in line: + funcs = [] + if mt in nuc: + funcs.append(nuc[mt].xs[nucT]) + # ## What if funcs is empty? + xs[-1].append(openmc.data.Sum(funcs)) + else: + raise ValueError(nuclide + " not in library") + + # 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(types_), len(self.nuclides))) + for l in range(len(types_)): + for n in range(len(self.nuclides)): + # ## + density = percent[n] + # ## + data[l, n] += density * xs[n][l](unionE) + + return data + def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") xml_element.set("name", nuclide[0].name) From ddf642e1956e3524084cedf12cfa014bdaf5c69b Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 6 Nov 2016 08:26:16 -0500 Subject: [PATCH 02/31] Incorporated ability to include densities and to plot s(a,b) data when in the material --- openmc/data/function.py | 57 ++++++++ openmc/material.py | 317 ++++++++++++++++++++++++++++++++++------ 2 files changed, 327 insertions(+), 47 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index e9ecf82c5..634bf5664 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -432,6 +432,63 @@ class Sum(EqualityMixin): self._functions = functions +class Piecewise(EqualityMixin): + """Piecewise composition of multiple functions. + + This class allows you to create a callable object which is composed + of other callable objects each over a set domain. + + Parameters + ---------- + functions : Iterable of Callable + Functions which are to be combined in a piecewise fashion + breakpoints : Iterable of float + The breakpoints between each function. The functions + in the functions attribute must be provided in order of + increasing domains. + + Attributes + ---------- + functions : Iterable of Callable + Functions which are to be combined in a piecewise fashion + breakpoints : Iterable of float + The breakpoints between each function + + """ + + def __init__(self, functions, breakpoints): + self.functions = functions + self.breakpoints = breakpoints + + def __call__(self, x): + i = np.searchsorted(self.breakpoints, x) + if isinstance(x, Iterable): + ans = np.empty_like(x) + for j in range(len(i)): + ans[j] = self.functions[i[j]](x[j]) + return ans + else: + return self.functions[i](x) + + @property + def functions(self): + return self._functions + + @property + def breakpoints(self): + return self._breakpoints + + @functions.setter + def functions(self, functions): + cv.check_type('functions', functions, Iterable, Callable) + self._functions = functions + + @breakpoints.setter + def breakpoints(self, breakpoints): + cv.check_iterable_type('breakpoints', breakpoints, Real) + self._breakpoints = breakpoints + + class ResonancesWithBackground(EqualityMixin): """Cross section in resolved resonance region. diff --git a/openmc/material.py b/openmc/material.py index 271acfd8f..46734a42f 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -6,9 +6,9 @@ from xml.etree import ElementTree as ET import sys from six import string_types -from scipy.interpolate import interp1d import numpy as np -import h5py +import scipy.constants as sc +from matplotlib import pyplot as plt import openmc import openmc.data @@ -32,10 +32,27 @@ DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum', # Supported keywords for material xs plotting _PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', - 'absorption'] + 'absorption', 'non-fission capture', 'n-alpha'] # MTs to sum to generate associated plot_types -_PLOT_TYPES_MT = {'total': (1,), 'scatter': (2, 4), 'elastic': (2,), - 'inelastic': (1,), 'fission': (18,), 'absorption': (102,)} +_PLOT_TYPES_MT = {'total': (2, 3,), + '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, 185, 186, 187, 188, 189, + 190, 194, 195, 196, 198, 199, 200, 875, 891), + 'elastic': (2,), + 'inelastic': (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, 185, 186, 187, 188, 189, + 190, 194, 195, 196, 198, 199, 200, 875, 891), + 'fission': (18,), 'absorption': (27,), + 'non-fission capture': (101,), + 'n-alpha': (107,)} class Material(object): @@ -580,7 +597,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: @@ -591,23 +608,26 @@ class Material(object): return nuclides - def plot_xs(self, library, types, labels=None): + def plot_xs(self, library, types, temperature=294., Erange=(1.E-5, 20.E6)): """Creates a figure of macroscopic cross sections for this material Parameters ---------- library : openmc.data.DataLibrary Library of data to use for plotting. - types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption'} or list thereof + types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'non-fission capture', 'n-alpha'} or list thereof The type of cross sections to include in the plot. This can either be an MT number, a tuple of MT numbers (indicating they are to be summed before plotting) or a string describing a common type. These values can be either a single value, or an iterable in the case of multiple sets per plot. - labels : str or list of str, optional - Labels to use in the plot legend for each type requested. - Single string (if there is only one type) or a list of strings - with a length matching the length of 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 Returns ------- @@ -616,63 +636,269 @@ class Material(object): """ + E, data, labels = self.calculate_xs(library, types, temperature) + + # 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=labels[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('Macroscopic 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 = 'Macroscopic Cross Section for ' + self.name + ax.set_title(title) + + return fig + + def calculate_xs(self, library, types, temperature=294.): + """Calculates macroscopic cross sections of a requested type + + Parameters + ---------- + library : openmc.data.DataLibrary + Library of data to use for plotting. + types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'n-alpha'} or list thereof + The type of cross sections to include in the plot. This can either + be an MT number, a tuple of MT numbers (indicating they are to be + summed before plotting) or a string describing a common type. + These values can be either a single value, or an iterable + in the case of multiple sets per plot. + 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. + + 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 + labels: Iterable of string-type + Name of cross section type for every type requested + + """ + # Check types - # ## Need to add error messages in some ELSEs - # ## Need to add in label handling if isinstance(types, Integral): cv.check_greater_than('types', types, 0) + labels = [openmc.data.REACTION_NAME[types]] + elif types in _PLOT_TYPES: + labels = [types] + # Replace with the MTs to sum to simplify downstream code + types = _PLOT_TYPES_MT[types] elif isinstance(types, list): - for line in types: - if isinstance(line, Integral): - cv.check_greater_than('line in types', line, 0) - elif isinstance(line, tuple): - for entry in line: + labels = [] + for t in range(len(types)): + if isinstance(types[t], Integral): + cv.check_greater_than('type in types', types[t], 0) + labels.append(openmc.data.REACTION_NAME[types[t]]) + elif isinstance(types[t], tuple): + labels.append('') + for e, entry in enumerate(types[t]): if isinstance(entry, Integral): - cv.check_greater_than('entry in line in types', + cv.check_greater_than('entry in type in types', entry, 0) - elif line in _PLOT_TYPES: + if e == len(types[t]) - 1: + labels[-1] += openmc.data.REACTION_NAME + else: + labels[-1] += openmc.data.REACTION_NAME + ' + ' + else: + raise ValueError("Invalid entry, " + "{}, in types".format(str(entry))) + elif types[t] in _PLOT_TYPES: + labels.append(types[t]) # Replace with the MTs to sum to simplify downstream code - line = _PLOT_TYPES_MT[line] + types[t] = _PLOT_TYPES_MT[types[t]] + else: + raise ValueError("Invalid type, " + "{}, in types".format(str(types[t]))) - # Convert temperature to format for accessing in the library + # Convert temperature to format needed for access in the library if self.temperature is not None: - T = "{}K".format(int(round(self.temperature))) + strT = "{}K".format(int(round(self.temperature))) + T = self.temperature else: - T = None + # ## What about default temperature? + cv.check_type('temperature', temperature, Real) + strT = "{}K".format(int(round(temperature))) + T = temperature if isinstance(types, (Integral, str)): types_ = [(types,)] elif isinstance(types, tuple): + types_ = [types] + else: types_ = types + # 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_data = nuclide + nuc, nuc_density, nuc_density_type = nuc_data + 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. + + # Pre-determine the nuclides which need s(a,b) data + sabs = {} + for nuc in nucs: + sabs[nuc.name] = None + + for sab_name in self._sab: + sab = openmc.data.ThermalScattering.from_hdf5( + library[sab_name]['path']) + for nuc in sab.nuclides: + sabs[nuc] = library[sab_name]['path'] + # Now we can create the data sets to be plotted xs = [] E = [] - percent = [] - - for n, nuclide in enumerate(self.nuclides): - nuc, pct, units = nuclide - # import pdb; pdb.set_trace() - percent.append(pct) - lib = library[nuc] + awrs = [] + n = -1 + for nuclide in nuclides.items(): + n += 1 + lib = library[nuclide[0]] + # nuc, nuc_data = nuclide + # lib = library[nuc] if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) - # ## Need to check temperatures - nucT = T - # ## - E.append(nuc.energy[nucT]) + awrs.append(nuc.atomic_weight_ratio) + if not percent_in_atom: + nuc_densities[n] = -nuc_densities[n] / awrs[-1] + # 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 l, line in enumerate(types_): - # Get the reaction from the nuclide + # Get the reaction xs data from the nuclide + funcs = [] for mt in line: - funcs = [] - if mt in nuc: + 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: funcs.append(nuc[mt].xs[nucT]) - # ## What if funcs is empty? xs[-1].append(openmc.data.Sum(funcs)) else: raise ValueError(nuclide + " not in library") + # 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 * \ + sc.Avogadro / sc.value('neutron mass in u') * 1.E-24 + nuc_densities = density * nuc_densities + # Condense the data for every nuclide # First create a union energy grid unionE = E[0] @@ -680,15 +906,12 @@ class Material(object): unionE = np.union1d(unionE, E[n]) # Now we can combine all the nuclidic data - data = np.zeros((len(types_), len(self.nuclides))) + data = np.zeros((len(types_), len(unionE))) for l in range(len(types_)): - for n in range(len(self.nuclides)): - # ## - density = percent[n] - # ## - data[l, n] += density * xs[n][l](unionE) + for n in range(len(nuclides)): + data[l, :] += nuc_densities[n] * xs[n][l](unionE) - return data + return unionE, data, labels def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") From 32f70cd21ac9a00337879dedb9260dc6def7e8a9 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 6 Nov 2016 08:47:12 -0500 Subject: [PATCH 03/31] Separated out routine to calculate the atom densities from the xs calculation routine --- openmc/material.py | 135 ++++++++++++++++++++++++++++++--------------- 1 file changed, 91 insertions(+), 44 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 46734a42f..39549dbd1 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -608,6 +608,94 @@ class Material(object): return nuclides + def get_nuclide_atom_densities(self, library): + """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. + + Returns + ------- + nuclides : dict + Dictionary whose keys are nuclide names and values are 2-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_data = nuclide + nuc, nuc_density, nuc_density_type = nuc_data + 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 = [] + n = -1 + for nuclide in nuclides.items(): + n += 1 + lib = library[nuclide[0]] + if lib is not None: + nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) + awrs.append(nuc.atomic_weight_ratio) + if not percent_in_atom: + nuc_densities[n] = -nuc_densities[n] / awrs[-1] + else: + raise ValueError(nuclide + " not in library") + + # 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 * \ + sc.Avogadro / sc.value('neutron mass in u') * 1.E-24 + nuc_densities = density * nuc_densities + + result = OrderedDict() + n = -1 + for nuc in nucs: + n += 1 + result[nuc] = (nuc, nuc_densities[n]) + + nuclides = result + return nuclides + def plot_xs(self, library, types, temperature=294., Erange=(1.E-5, 20.E6)): """Creates a figure of macroscopic cross sections for this material @@ -732,7 +820,6 @@ class Material(object): strT = "{}K".format(int(round(self.temperature))) T = self.temperature else: - # ## What about default temperature? cv.check_type('temperature', temperature, Real) strT = "{}K".format(int(round(temperature))) T = temperature @@ -744,41 +831,18 @@ class Material(object): else: types_ = types - # 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 + # Expand elements in to nuclides with atomic densities + nuclides = self.get_nuclide_atom_densities(library) # 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_data = nuclide - nuc, nuc_density, nuc_density_type = nuc_data + nuc, nuc_density = nuc_data 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. # Pre-determine the nuclides which need s(a,b) data sabs = {} @@ -794,7 +858,6 @@ class Material(object): # Now we can create the data sets to be plotted xs = [] E = [] - awrs = [] n = -1 for nuclide in nuclides.items(): n += 1 @@ -803,9 +866,6 @@ class Material(object): # lib = library[nuc] if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) - awrs.append(nuc.atomic_weight_ratio) - if not percent_in_atom: - nuc_densities[n] = -nuc_densities[n] / awrs[-1] # Obtain the nearest temperature if strT in nuc.temperatures: nucT = strT @@ -886,19 +946,6 @@ class Material(object): else: raise ValueError(nuclide + " not in library") - # 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 * \ - sc.Avogadro / sc.value('neutron mass in u') * 1.E-24 - nuc_densities = density * nuc_densities - # Condense the data for every nuclide # First create a union energy grid unionE = E[0] From 273e0a32cef1f3e828aba4b3e5caf1c297bf7a71 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 6 Nov 2016 12:42:31 -0500 Subject: [PATCH 04/31] Minor cleanup --- openmc/material.py | 30 +++++++++++++++--------------- 1 file changed, 15 insertions(+), 15 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 39549dbd1..587abf3e1 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -32,7 +32,8 @@ DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum', # Supported keywords for material xs plotting _PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', - 'absorption', 'non-fission capture', 'n-alpha'] + 'absorption', 'capture', 'n-alpha'] + # MTs to sum to generate associated plot_types _PLOT_TYPES_MT = {'total': (2, 3,), 'scatter': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, 29, @@ -51,7 +52,7 @@ _PLOT_TYPES_MT = {'total': (2, 3,), 180, 181, 183, 184, 185, 186, 187, 188, 189, 190, 194, 195, 196, 198, 199, 200, 875, 891), 'fission': (18,), 'absorption': (27,), - 'non-fission capture': (101,), + 'capture': (101,), 'n-alpha': (107,)} @@ -620,11 +621,13 @@ class Material(object): Returns ------- nuclides : dict - Dictionary whose keys are nuclide names and values are 2-tuples of + Dictionary whose keys are nuclide names and values are tuples of (nuclide, density in atom/b-cm) """ + cv.check_type('library', library, openmc.data.DataLibrary) + # Expand elements in to nuclides nuclides = self.get_nuclide_densities() @@ -672,7 +675,7 @@ class Material(object): if not percent_in_atom: nuc_densities[n] = -nuc_densities[n] / awrs[-1] else: - raise ValueError(nuclide + " not in library") + raise ValueError(nuclide[0] + " not in library") # Now that we have the awr, lets finish calculating densities sum_percent = np.sum(nuc_densities) @@ -687,13 +690,12 @@ class Material(object): sc.Avogadro / sc.value('neutron mass in u') * 1.E-24 nuc_densities = density * nuc_densities - result = OrderedDict() + nuclides = OrderedDict() n = -1 for nuc in nucs: n += 1 - result[nuc] = (nuc, nuc_densities[n]) + nuclides[nuc] = (nuc, nuc_densities[n]) - nuclides = result return nuclides def plot_xs(self, library, types, temperature=294., Erange=(1.E-5, 20.E6)): @@ -703,7 +705,7 @@ class Material(object): ---------- library : openmc.data.DataLibrary Library of data to use for plotting. - types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'non-fission capture', 'n-alpha'} or list thereof + types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'capture', 'n-alpha'} or list thereof The type of cross sections to include in the plot. This can either be an MT number, a tuple of MT numbers (indicating they are to be summed before plotting) or a string describing a common type. @@ -781,6 +783,7 @@ class Material(object): """ # Check types + cv.check_type('library', library, openmc.data.DataLibrary) if isinstance(types, Integral): cv.check_greater_than('types', types, 0) labels = [openmc.data.REACTION_NAME[types]] @@ -834,17 +837,16 @@ class Material(object): # Expand elements in to nuclides with atomic densities nuclides = self.get_nuclide_atom_densities(library) - # For ease of processing split out nuc, nuc_density, - # and nuc_density_type in to separate arrays + # For ease of processing split out nuc and nuc_density nucs = [] nuc_densities = [] for nuclide in nuclides.items(): - nuc, nuc_data = nuclide + nuc_name, nuc_data = nuclide nuc, nuc_density = nuc_data nucs.append(nuc) nuc_densities.append(nuc_density) - # Pre-determine the nuclides which need s(a,b) data + # Identify the nuclides which need S(a,b) data sabs = {} for nuc in nucs: sabs[nuc.name] = None @@ -862,8 +864,6 @@ class Material(object): for nuclide in nuclides.items(): n += 1 lib = library[nuclide[0]] - # nuc, nuc_data = nuclide - # lib = library[nuc] if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) # Obtain the nearest temperature @@ -944,7 +944,7 @@ class Material(object): funcs.append(nuc[mt].xs[nucT]) xs[-1].append(openmc.data.Sum(funcs)) else: - raise ValueError(nuclide + " not in library") + raise ValueError(nuclide[0] + " not in library") # Condense the data for every nuclide # First create a union energy grid From 739dd53c9bb14ea96a93eb4651a400eedf294d75 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 6 Nov 2016 15:48:36 -0500 Subject: [PATCH 05/31] Yay! Now can handle data operations aside from summing, allowing for cheaper calculations of multi-MT datasets (i.e., its easier to find scatter when its just total-absorption), and I put in the ability to get nu-fission and nu-scatter by way of the product yield. Whew. --- openmc/data/function.py | 55 +++++++++++ openmc/material.py | 206 ++++++++++++++++++++++------------------ 2 files changed, 171 insertions(+), 90 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 634bf5664..489910c15 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -397,6 +397,61 @@ class Polynomial(np.polynomial.Polynomial, Function1D): return cls(dataset.value) +class Combination(EqualityMixin): + """Combination of multiple functions with a user-defined operator + + This class allows you to create a callable object which represents the + combination of other callable objects by way of a user-defined operator. + + Parameters + ---------- + functions : Iterable of Callable + Functions which are to be added together + operations : Iterable of numpy.ufunc + Operations to perform between functions; note that the standard order + of operations (i.e., PEMDAS) will not be followed + + + Attributes + ---------- + functions : Iterable of Callable + Functions which are to be added together + operations : Iterable of numpy.ufunc + Operations to perform between functions + + """ + + def __init__(self, functions, operations): + self.functions = functions + self.operations = operations + + def __call__(self, x): + ans = np.zeros_like(x) + for i, operation in enumerate(self.operations): + ans = operation(ans, self.functions[i](x)) + return ans + + @property + def functions(self): + return self._functions + + @functions.setter + def functions(self, functions): + cv.check_type('functions', functions, Iterable, Callable) + self._functions = functions + + @property + def operations(self): + return self._operations + + @operations.setter + def operations(self, operations): + cv.check_type('operations', operations, Iterable, np.ufunc) + length = len(self.functions) + cv.check_length('operations', operations, length, length_max=length) + self._operations = operations + + class Sum(EqualityMixin): """Sum of multiple functions. diff --git a/openmc/material.py b/openmc/material.py index 587abf3e1..0a6402b78 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -3,7 +3,6 @@ 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 @@ -30,30 +29,60 @@ 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 +# Supported keywords for material xs plotting _PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', - 'absorption', 'capture', 'n-alpha'] + 'absorption', 'capture', 'nu-fission', 'nu-scatter'] -# MTs to sum to generate associated plot_types -_PLOT_TYPES_MT = {'total': (2, 3,), - '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, 185, 186, 187, 188, 189, - 190, 194, 195, 196, 198, 199, 200, 875, 891), - 'elastic': (2,), - 'inelastic': (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, 185, 186, 187, 188, 189, - 190, 194, 195, 196, 198, 199, 200, 875, 891), - 'fission': (18,), 'absorption': (27,), - 'capture': (101,), - 'n-alpha': (107,)} +# 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)} +# 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, np.subtract), + 'elastic': (np.add,), + 'inelastic': (np.add, np.subtract, np.subtract), + 'fission': (np.add,), 'absorption': (np.add,), + 'capture': (np.add,), 'nu-fission': (np.add,), + '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, + np.add)} +# 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)} class Material(object): @@ -705,18 +734,16 @@ class Material(object): ---------- library : openmc.data.DataLibrary Library of data to use for plotting. - types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'capture', 'n-alpha'} or list thereof - The type of cross sections to include in the plot. This can either - be an MT number, a tuple of MT numbers (indicating they are to be - summed before plotting) or a string describing a common type. - These values can be either a single value, or an iterable - in the case of multiple sets per plot. + types : Iterable of {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption'} + The type of cross sections to include in the plot. In addition to + the above, a custom string can be used which includes MT numbers + and operations to perform such as '2 + 3' 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 Energy range (in eV) to plot the cross section within Returns @@ -726,7 +753,7 @@ class Material(object): """ - E, data, labels = self.calculate_xs(library, types, temperature) + E, data = self.calculate_xs(library, types, temperature) # Generate the plot fig = plt.figure() @@ -736,7 +763,7 @@ class Material(object): 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=labels[i]) + 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])) @@ -758,13 +785,11 @@ class Material(object): ---------- library : openmc.data.DataLibrary Library of data to use for plotting. - types : int, tuples of int, {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'n-alpha'} or list thereof - The type of cross sections to include in the plot. This can either - be an MT number, a tuple of MT numbers (indicating they are to be - summed before plotting) or a string describing a common type. - These values can be either a single value, or an iterable - in the case of multiple sets per plot. - temperature: float, optional + types : Iterable of {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption'} + The type of cross sections to include in the plot. In addition to + the above, a custom string can be used which includes MT numbers + and operations to perform such as '2 + 3' + 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 @@ -772,51 +797,30 @@ class Material(object): Returns ------- - unionE: numpy.array + unionE : numpy.array Energies at which cross sections are calculated, in units of eV - data: numpy.ndarray + data : numpy.ndarray Macroscopic cross sections calculated at the energy grid described by unionE - labels: Iterable of string-type - Name of cross section type for every type requested """ # Check types cv.check_type('library', library, openmc.data.DataLibrary) - if isinstance(types, Integral): - cv.check_greater_than('types', types, 0) - labels = [openmc.data.REACTION_NAME[types]] - elif types in _PLOT_TYPES: - labels = [types] - # Replace with the MTs to sum to simplify downstream code - types = _PLOT_TYPES_MT[types] - elif isinstance(types, list): - labels = [] - for t in range(len(types)): - if isinstance(types[t], Integral): - cv.check_greater_than('type in types', types[t], 0) - labels.append(openmc.data.REACTION_NAME[types[t]]) - elif isinstance(types[t], tuple): - labels.append('') - for e, entry in enumerate(types[t]): - if isinstance(entry, Integral): - cv.check_greater_than('entry in type in types', - entry, 0) - if e == len(types[t]) - 1: - labels[-1] += openmc.data.REACTION_NAME - else: - labels[-1] += openmc.data.REACTION_NAME + ' + ' - else: - raise ValueError("Invalid entry, " - "{}, in types".format(str(entry))) - elif types[t] in _PLOT_TYPES: - labels.append(types[t]) - # Replace with the MTs to sum to simplify downstream code - types[t] = _PLOT_TYPES_MT[types[t]] - else: - raise ValueError("Invalid type, " - "{}, in types".format(str(types[t]))) + 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() # Convert temperature to format needed for access in the library if self.temperature is not None: @@ -827,13 +831,6 @@ class Material(object): strT = "{}K".format(int(round(temperature))) T = temperature - if isinstance(types, (Integral, str)): - types_ = [(types,)] - elif isinstance(types, tuple): - types_ = [types] - else: - types_ = types - # Expand elements in to nuclides with atomic densities nuclides = self.get_nuclide_atom_densities(library) @@ -925,10 +922,11 @@ class Material(object): else: E.append(nuc.energy[nucT]) xs.append([]) - for l, line in enumerate(types_): + for i, mt_set in enumerate(mts): # Get the reaction xs data from the nuclide funcs = [] - for mt in line: + 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 @@ -941,8 +939,36 @@ class Material(object): else: funcs.append(nuc[mt].xs[nucT]) elif mt in nuc: - funcs.append(nuc[mt].xs[nucT]) - xs[-1].append(openmc.data.Sum(funcs)) + 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.add, 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.add, 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]) + else: + funcs.append(lambda x: 0.) + xs[-1].append(openmc.data.Combination(funcs, op)) else: raise ValueError(nuclide[0] + " not in library") @@ -953,12 +979,12 @@ class Material(object): unionE = np.union1d(unionE, E[n]) # Now we can combine all the nuclidic data - data = np.zeros((len(types_), len(unionE))) - for l in range(len(types_)): + data = np.zeros((len(mts), len(unionE))) + for l in range(len(mts)): for n in range(len(nuclides)): data[l, :] += nuc_densities[n] * xs[n][l](unionE) - return unionE, data, labels + return unionE, data def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") From 3f610d03ca66fe318560875aa816e7dccfed052f Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 6 Nov 2016 19:43:32 -0500 Subject: [PATCH 06/31] Added the ability to divide by a type, so now we can easily (and not with complicated routines within calculate_xs) get nu-scatter/scatter, nu-fission/fission and something like a moderating ratio (just needs xi). --- openmc/material.py | 59 ++++++++++++++++++++++++++++++++++------------ 1 file changed, 44 insertions(+), 15 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 0a6402b78..1e0ff5661 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -31,7 +31,7 @@ DENSITY_UNITS = ['g/cm3', 'g/cc', 'kg/cm3', 'atom/b-cm', 'atom/cm3', 'sum', # Supported keywords for material xs plotting _PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', - 'absorption', 'capture', 'nu-fission', 'nu-scatter'] + '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,), @@ -44,7 +44,8 @@ _PLOT_TYPES_MT = {'total': (1,), 'scatter': (1, 27), 'elastic': (2,), 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)} + 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.add, np.subtract), @@ -64,7 +65,8 @@ _PLOT_TYPES_OP = {'total': (np.add,), 'scatter': (np.add, np.subtract), 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': (np.add,)} # 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), @@ -82,7 +84,8 @@ _PLOT_TYPES_YIELD = {'total': (False,), 'scatter': (False, False), True, True, True, True, True, True, True, True, True, True, True, True, True, True, True, - True)} + True), + 'unity': (False,)} class Material(object): @@ -727,17 +730,20 @@ class Material(object): return nuclides - def plot_xs(self, library, types, temperature=294., Erange=(1.E-5, 20.E6)): + def plot_xs(self, library, types, divisor_types=None, temperature=294., + Erange=(1.E-5, 20.E6)): """Creates a figure of macroscopic cross sections for this material Parameters ---------- library : openmc.data.DataLibrary Library of data to use for plotting. - types : Iterable of {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption'} - The type of cross sections to include in the plot. In addition to - the above, a custom string can be used which includes MT numbers - and operations to perform such as '2 + 3' + 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 @@ -755,6 +761,26 @@ class Material(object): E, data = self.calculate_xs(library, types, temperature) + if divisor_types: + cv.check_length('divisor types', divisor_types, len(types), + len(types)) + Ediv, data_div = self.calculate_xs(library, divisor_types, + temperature) + + # 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, :]), + 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) @@ -785,10 +811,8 @@ class Material(object): ---------- library : openmc.data.DataLibrary Library of data to use for plotting. - types : Iterable of {'total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption'} - The type of cross sections to include in the plot. In addition to - the above, a custom string can be used which includes MT numbers - and operations to perform such as '2 + 3' + types : Iterable of values of _PLOT_TYPES + The type of cross sections to include in the plot. temperature : float, optional Temperature in Kelvin to plot. If not specified, a default temperature of 294K will be plotted. Note that the nearest @@ -966,6 +990,8 @@ class Material(object): 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)) @@ -981,8 +1007,11 @@ class Material(object): # Now we can combine all the nuclidic data data = np.zeros((len(mts), len(unionE))) for l in range(len(mts)): - for n in range(len(nuclides)): - data[l, :] += nuc_densities[n] * xs[n][l](unionE) + if types[l] == 'unity': + data[l, :] = 1. + else: + for n in range(len(nuclides)): + data[l, :] += nuc_densities[n] * xs[n][l](unionE) return unionE, data From 0f9c4455a5bf5a8af1037abcf4cae52ad47c1726 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 7 Nov 2016 18:58:26 -0500 Subject: [PATCH 07/31] Cleanup of function --- openmc/data/function.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/openmc/data/function.py b/openmc/data/function.py index 489910c15..e87f0a9c4 100644 --- a/openmc/data/function.py +++ b/openmc/data/function.py @@ -401,7 +401,8 @@ class Combination(EqualityMixin): """Combination of multiple functions with a user-defined operator This class allows you to create a callable object which represents the - combination of other callable objects by way of a user-defined operator. + combination of other callable objects by way of a series of user-defined + operators connecting each of the callable objects. Parameters ---------- @@ -409,7 +410,7 @@ class Combination(EqualityMixin): Functions which are to be added together operations : Iterable of numpy.ufunc Operations to perform between functions; note that the standard order - of operations (i.e., PEMDAS) will not be followed + of operations (i.e., PEMDAS) will not be followed. Attributes @@ -491,7 +492,7 @@ class Piecewise(EqualityMixin): """Piecewise composition of multiple functions. This class allows you to create a callable object which is composed - of other callable objects each over a set domain. + of multiple other callable objects, each applying to a specific interval Parameters ---------- @@ -500,7 +501,7 @@ class Piecewise(EqualityMixin): breakpoints : Iterable of float The breakpoints between each function. The functions in the functions attribute must be provided in order of - increasing domains. + increasing intervals. Attributes ---------- From 131fdce5b34836412f7570cf5ee3ec3fac903048 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 7 Nov 2016 19:47:29 -0500 Subject: [PATCH 08/31] fixing failing tests --- openmc/material.py | 12 ++++++++---- 1 file changed, 8 insertions(+), 4 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index 1e0ff5661..a1281fef1 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -6,8 +6,6 @@ from xml.etree import ElementTree as ET from six import string_types import numpy as np -import scipy.constants as sc -from matplotlib import pyplot as plt import openmc import openmc.data @@ -732,7 +730,8 @@ class Material(object): def plot_xs(self, library, types, divisor_types=None, temperature=294., Erange=(1.E-5, 20.E6)): - """Creates a figure of macroscopic cross sections for this material + """Creates a figure of continuous-energy macroscopic cross sections + for this material Parameters ---------- @@ -759,6 +758,8 @@ class Material(object): """ + from matplotlib import pyplot as plt + E, data = self.calculate_xs(library, types, temperature) if divisor_types: @@ -805,7 +806,8 @@ class Material(object): return fig def calculate_xs(self, library, types, temperature=294.): - """Calculates macroscopic cross sections of a requested type + """Calculates continuous-energy macroscopic cross sections of a + requested type Parameters ---------- @@ -829,6 +831,8 @@ class Material(object): """ + import scipy.constants as sc + # Check types cv.check_type('library', library, openmc.data.DataLibrary) cv.check_iterable_type('types', types, str) From 8de82dcbd32f26ae198f635247c5012623e161dc Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Tue, 8 Nov 2016 20:29:33 -0500 Subject: [PATCH 09/31] Modifications to support new function routine --- openmc/material.py | 29 ++++++++++++++--------------- 1 file changed, 14 insertions(+), 15 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index a1281fef1..37dbe0d2d 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -46,11 +46,10 @@ _PLOT_TYPES_MT = {'total': (1,), 'scatter': (1, 27), 'elastic': (2,), '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.add, np.subtract), - 'elastic': (np.add,), - 'inelastic': (np.add, np.subtract, np.subtract), - 'fission': (np.add,), 'absorption': (np.add,), - 'capture': (np.add,), 'nu-fission': (np.add,), +_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, @@ -62,9 +61,8 @@ _PLOT_TYPES_OP = {'total': (np.add,), 'scatter': (np.add, np.subtract), 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': (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), @@ -656,6 +654,8 @@ class Material(object): """ + import scipy.constants as sc + cv.check_type('library', library, openmc.data.DataLibrary) # Expand elements in to nuclides @@ -698,7 +698,7 @@ class Material(object): n = -1 for nuclide in nuclides.items(): n += 1 - lib = library[nuclide[0]] + lib = library.get_by_materials(nuclide[0]) if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) awrs.append(nuc.atomic_weight_ratio) @@ -878,9 +878,9 @@ class Material(object): for sab_name in self._sab: sab = openmc.data.ThermalScattering.from_hdf5( - library[sab_name]['path']) + library.get_by_materials(sab_name)['path']) for nuc in sab.nuclides: - sabs[nuc] = library[sab_name]['path'] + sabs[nuc] = library.get_by_materials(sab_name)['path'] # Now we can create the data sets to be plotted xs = [] @@ -888,7 +888,7 @@ class Material(object): n = -1 for nuclide in nuclides.items(): n += 1 - lib = library[nuclide[0]] + lib = library.get_by_materials(nuclide[0]) if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) # Obtain the nearest temperature @@ -974,7 +974,7 @@ class Material(object): prod.emission_mode == 'total': func = openmc.data.Combination( [nuc[mt].xs[nucT], prod.yield_], - [np.add, np.multiply]) + [np.multiply]) funcs.append(func) found_it = True break @@ -984,8 +984,7 @@ class Material(object): prod.emission_mode == 'prompt': func = openmc.data.Combination( [nuc[mt].xs[nucT], - prod.yield_], - [np.add, np.multiply]) + prod.yield_], [np.multiply]) funcs.append(func) found_it = True break From 76991a223f62c666b3f7fefcf1e9a94f61d74af3 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 9 Nov 2016 18:10:41 -0500 Subject: [PATCH 10/31] Extended plots to elements and nuclides --- openmc/__init__.py | 1 + openmc/element.py | 209 +++++++++++++++++++++++++++++++++- openmc/material.py | 267 +++++++++++-------------------------------- openmc/nuclide.py | 272 +++++++++++++++++++++++++++++++++++++++++++- openmc/plot_data.py | 57 ++++++++++ 5 files changed, 597 insertions(+), 209 deletions(-) create mode 100644 openmc/plot_data.py diff --git a/openmc/__init__.py b/openmc/__init__.py index 48f43be91..98de9b70d 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -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 * diff --git a/openmc/element.py b/openmc/element.py index d56c3b350..90bde474d 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -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 diff --git a/openmc/material.py b/openmc/material.py index 37dbe0d2d..c357da721 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -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 diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 55afd4e9c..31b658a1c 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -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 diff --git a/openmc/plot_data.py b/openmc/plot_data.py new file mode 100644 index 000000000..03e929b52 --- /dev/null +++ b/openmc/plot_data.py @@ -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,)} From 7524637e2351670b4d8d3ca56c291dfafd0e16fc Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 05:20:38 -0500 Subject: [PATCH 11/31] Added kwargs to plots, switched to taking a file string for the cross_sections.xml file instead of an initialized DataLibrary, and avoided precision issues when calculating the inelastic xs --- openmc/element.py | 16 +++++++--------- openmc/material.py | 19 ++++++++----------- openmc/nuclide.py | 16 +++++++--------- openmc/plot_data.py | 36 +++++++++++++++++++++++++++++++++--- 4 files changed, 55 insertions(+), 32 deletions(-) diff --git a/openmc/element.py b/openmc/element.py index 90bde474d..80777fcfe 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -261,7 +261,8 @@ class Element(object): return isotopes def plot_xs(self, types, divisor_types=None, temperature=294., - Erange=(1.E-5, 20.E6), enrichment=None, cross_sections=None): + Erange=(1.E-5, 20.E6), enrichment=None, cross_sections=None, + **kwargs): """Creates a figure of continuous-energy microscopic cross sections for this element @@ -286,6 +287,9 @@ class Element(object): (natural composition). cross_sections : str, optional Location of cross_sections.xml file. Default is None. + **kwargs + All keyword arguments are passed to + :func:`matplotlib.pyplot.figure`. Returns ------- @@ -319,22 +323,16 @@ class Element(object): data = data_new # Generate the plot - fig = plt.figure() + fig = plt.figure(**kwargs) 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) @@ -426,7 +424,7 @@ class Element(object): 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'] + sabs[nuc] = library.get_by_material(sab_name)['path'] # Now we can create the data sets to be plotted xs = [] diff --git a/openmc/material.py b/openmc/material.py index c357da721..36f728ad1 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -657,7 +657,7 @@ class Material(object): n = -1 for nuclide in nuclides.items(): n += 1 - lib = library.get_by_materials(nuclide[0]) + lib = library.get_by_material(nuclide[0]) if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) awrs.append(nuc.atomic_weight_ratio) @@ -688,7 +688,7 @@ class Material(object): return nuclides def plot_xs(self, types, divisor_types=None, temperature=294., - Erange=(1.E-5, 20.E6), cross_sections=None): + Erange=(1.E-5, 20.E6), cross_sections=None, **kwargs): """Creates a figure of continuous-energy macroscopic cross sections for this material @@ -709,6 +709,9 @@ class Material(object): Energy range (in eV) to plot the cross section within cross_sections : str, optional Location of cross_sections.xml file. Default is None. + **kwargs + All keyword arguments are passed to + :func:`matplotlib.pyplot.figure`. Returns ------- @@ -742,22 +745,16 @@ class Material(object): data = data_new # Generate the plot - fig = plt.figure() + fig = plt.figure(**kwargs) 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('Macroscopic 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 = 'Macroscopic Cross Section for ' + self.name ax.set_title(title) @@ -846,9 +843,9 @@ class Material(object): for sab_name in self._sab: sab = openmc.data.ThermalScattering.from_hdf5( - library.get_by_materials(sab_name)['path']) + library.get_by_material(sab_name)['path']) for nuc in sab.nuclides: - sabs[nuc] = library.get_by_materials(sab_name)['path'] + sabs[nuc] = library.get_by_material(sab_name)['path'] # Now we can create the data sets to be plotted xs = [] diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 31b658a1c..938b6a666 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -98,7 +98,8 @@ class Nuclide(object): 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): + Erange=(1.E-5, 20.E6), sab_name=None, cross_sections=None, + **kwargs): """Creates a figure of continuous-energy cross sections for this nuclide @@ -121,6 +122,9 @@ class Nuclide(object): 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. + **kwargs + All keyword arguments are passed to + :func:`matplotlib.pyplot.figure`. Returns ------- @@ -154,23 +158,17 @@ class Nuclide(object): data = data_new # Generate the plot - fig = plt.figure() + fig = plt.figure(**kwargs) 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) @@ -245,7 +243,7 @@ class Nuclide(object): # Now we can create the data sets to be plotted E = [] xs = [] - lib = library.get_by_materials(self.name) + lib = library.get_by_material(self.name) if lib is not None: nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) # Obtain the nearest temperature diff --git a/openmc/plot_data.py b/openmc/plot_data.py index 03e929b52..ff3d799a2 100644 --- a/openmc/plot_data.py +++ b/openmc/plot_data.py @@ -6,7 +6,14 @@ PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', # 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,), + 'inelastic': (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), + 'fission': (18,), 'absorption': (27,), 'capture': (101,), 'nu-fission': (18,), 'nu-scatter': (2, 4, 11, 16, 17, 22, 23, 24, 25, 28, 29, 30, @@ -20,7 +27,18 @@ PLOT_TYPES_MT = {'total': (2, 3,), 'scatter': (1, 27), 'elastic': (2,), # 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), + 'inelastic': (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), 'fission': (), 'absorption': (), 'capture': (), 'nu-fission': (), 'nu-scatter': (np.add, np.add, np.add, np.add, np.add, @@ -38,7 +56,19 @@ PLOT_TYPES_OP = {'total': (np.add,), 'scatter': (np.subtract,), 'elastic': (), '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), + 'elastic': (False,), + 'inelastic': (False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False, + False, False, False, False, False), 'fission': (False,), 'absorption': (False,), 'capture': (False,), 'nu-fission': (True,), 'nu-scatter': (True, True, True, True, True, From 4fe275b38163131947c75aab8ed4b63980a25183 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 06:22:38 -0500 Subject: [PATCH 12/31] Code and resultant plot cleanup; also ensured every plot type explicitly included MT=2 so that S(a,b) data could properly be spliced in --- openmc/element.py | 76 ++++++++++++++++++++++++--------------------- openmc/material.py | 38 +++++++++++------------ openmc/nuclide.py | 13 ++++++-- openmc/plot_data.py | 74 +++++++++++++------------------------------ 4 files changed, 91 insertions(+), 110 deletions(-) diff --git a/openmc/element.py b/openmc/element.py index 80777fcfe..b1cd19e9b 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -1,4 +1,5 @@ from collections import OrderedDict +from numbers import Real import re import sys import os @@ -256,13 +257,13 @@ 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 def plot_xs(self, types, divisor_types=None, temperature=294., - Erange=(1.E-5, 20.E6), enrichment=None, cross_sections=None, - **kwargs): + Erange=(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 @@ -281,12 +282,14 @@ class Element(object): 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. 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. **kwargs All keyword arguments are passed to :func:`matplotlib.pyplot.figure`. @@ -300,13 +303,14 @@ class Element(object): from matplotlib import pyplot as plt - E, data = self.calculate_xs(types, temperature, cross_sections) + 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, - cross_sections) + sab_name, cross_sections) # Create a new union grid, interpolate data and data_div on to that # grid, and then do the actual division @@ -326,11 +330,20 @@ class Element(object): fig = plt.figure(**kwargs) ax = fig.add_subplot(111) 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.: - ax.loglog(E, data[i, :], label=types[i]) + plot_func(E, data[i, :], label=types[i]) ax.set_xlabel('Energy [eV]') - ax.set_ylabel('Elemental Cross Section [1/cm]') + if divisor_types: + ax.set_ylabel('Elemental Cross Section') + else: + ax.set_ylabel('Elemental Cross Section [b]') ax.legend(loc='best') ax.set_xlim(Erange) if self.name is not None: @@ -339,8 +352,8 @@ class Element(object): return fig - def calculate_xs(self, types, temperature=294., enrichment=None, - cross_sections=None): + def calculate_xs(self, types, temperature=294., sab_name=None, + cross_sections=None, enrichment=None): """Calculates continuous-energy macroscopic cross sections of a requested type @@ -353,12 +366,14 @@ class Element(object): 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). - cross_sections : str, optional - Location of cross_sections.xml file. Default is None. Returns ------- @@ -390,10 +405,6 @@ class Element(object): 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: @@ -407,20 +418,17 @@ class Element(object): raise ValueError("cross_sections or OPENMC_CROSS_SECTIONS " "environmental variable must be set") + # Expand elements in to nuclides with atomic densities + nuclides = self.expand(100., 'ao', enrichment=enrichment, + cross_sections=cross_sections) + # 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) + nuc_fractions = [nuclide[1] for nuclide in nuclides] - # Identify the nuclides which need S(a,b) data + # Identify the nuclides which have S(a,b) data sabs = {} - for nuc in nucs: - sabs[nuc.name] = None - + for nuclide in nuclides: + sabs[nuclide[0].name] = None if sab_name: sab = openmc.data.ThermalScattering.from_hdf5(sab_name) for nuc in sab.nuclides: @@ -429,14 +437,10 @@ class Element(object): # 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) + for nuclide in nuclides: + sab_tab = sabs[nuclide[0].name] + temp_E, temp_xs = nuclide[0].calculate_xs(types, temperature, + sab_tab, cross_sections) E.append(temp_E) xs.append(temp_xs) diff --git a/openmc/material.py b/openmc/material.py index 3b6c784b3..ec68913ca 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -773,11 +773,20 @@ class Material(object): fig = plt.figure(**kwargs) ax = fig.add_subplot(111) 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.: - ax.loglog(E, data[i, :], label=types[i]) + plot_func(E, data[i, :], label=types[i]) ax.set_xlabel('Energy [eV]') - ax.set_ylabel('Macroscopic Cross Section [1/cm]') + if divisor_types: + ax.set_ylabel('Macroscopic Cross Section') + else: + ax.set_ylabel('Macroscopic Cross Section [1/cm]') ax.legend(loc='best') ax.set_xlim(Erange) if self.name is not None: @@ -853,19 +862,12 @@ class Material(object): nuclides = self.get_nuclide_atom_densities(cross_sections) # For ease of processing split out nuc and nuc_density - nucs = [] - nuc_densities = [] - for nuclide in nuclides.items(): - nuc_name, nuc_data = nuclide - nuc, nuc_density = nuc_data - nucs.append(nuc) - nuc_densities.append(nuc_density) + nuc_densities = [nuclide[1][1] for nuclide in nuclides.items()] - # Identify the nuclides which need S(a,b) data + # Identify the nuclides which have S(a,b) data sabs = {} - for nuc in nucs: - sabs[nuc.name] = None - + for nuclide in nuclides.items(): + sabs[nuclide[0].name] = None for sab_name in self._sab: sab = openmc.data.ThermalScattering.from_hdf5( library.get_by_material(sab_name)['path']) @@ -875,14 +877,10 @@ class Material(object): # 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, T, sab_tab, - cross_sections) + sab_tab = sabs[nuclide[0].name] + temp_E, temp_xs = nuclide[0].calculate_xs(types, T, sab_tab, + cross_sections) E.append(temp_E) xs.append(temp_xs) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 938b6a666..fda44c43a 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -162,11 +162,20 @@ class Nuclide(object): ax = fig.add_subplot(111) 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.: - ax.loglog(E, to_plot, label=types[i]) + plot_func(E, to_plot, label=types[i]) ax.set_xlabel('Energy [eV]') - ax.set_ylabel('Microscopic Cross Section [b]') + if divisor_types: + ax.set_ylabel('Microscopic Cross Section') + else: + ax.set_ylabel('Microscopic Cross Section [b]') ax.legend(loc='best') ax.set_xlim(Erange) if self.name is not None: diff --git a/openmc/plot_data.py b/openmc/plot_data.py index ff3d799a2..ac59c16ea 100644 --- a/openmc/plot_data.py +++ b/openmc/plot_data.py @@ -5,7 +5,15 @@ 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,), +PLOT_TYPES_MT = {'total': (2, 3,), + '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), + 'elastic': (2,), 'inelastic': (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, @@ -26,62 +34,24 @@ PLOT_TYPES_MT = {'total': (2, 3,), 'scatter': (1, 27), 'elastic': (2,), '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.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), +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, 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), + 'nu-scatter': (np.add,) * (len(PLOT_TYPES_MT['nu-scatter']) - 1), 'unity': ()} + # Whether or not to multiply the reaction by the yield as well -PLOT_TYPES_YIELD = {'total': (False, False), 'scatter': (False, False), +PLOT_TYPES_YIELD = {'total': (False, False), + 'scatter': (False,) * len(PLOT_TYPES_MT['scatter']), 'elastic': (False,), - 'inelastic': (False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False, - False, False, False, False, False), + 'inelastic': (False,) * len(PLOT_TYPES_MT['inelastic']), '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), + 'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']), 'unity': (False,)} + +# Types of plots to plot linearly in y +PLOT_TYPES_LINEAR = ['nu-fission / fission', 'nu-scatter / scatter'] From 240028d13dc620e4e52ed1e1f16c0cebb019158f Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 14:50:24 -0500 Subject: [PATCH 13/31] Added ability to provide an axis to the plot_xs routines, and added a type for slowing-down power so we can plot moderating ratios. Why not right? --- openmc/element.py | 21 ++++++++++++++++----- openmc/material.py | 26 ++++++++++++++++++-------- openmc/nuclide.py | 27 +++++++++++++++++++++------ openmc/plot_data.py | 25 +++++++++++++++++++++---- 4 files changed, 76 insertions(+), 23 deletions(-) diff --git a/openmc/element.py b/openmc/element.py index b1cd19e9b..d90cdf4fc 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -261,7 +261,7 @@ class Element(object): return isotopes - def plot_xs(self, types, divisor_types=None, temperature=294., + def plot_xs(self, types, divisor_types=None, temperature=294., axis=None, Erange=(1.E-5, 20.E6), sab_name=None, cross_sections=None, enrichment=None, **kwargs): """Creates a figure of continuous-energy microscopic cross sections @@ -280,6 +280,9 @@ class Element(object): 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. Erange : tuple of floats Energy range (in eV) to plot the cross section within sab_name : str, optional @@ -297,7 +300,10 @@ class Element(object): Returns ------- fig : matplotlib.figure.Figure - Matplotlib Figure of the generated macroscopic cross section + If axis is None, then a Matplotlib Figure of the generated + macroscopic cross section will be returned. Otherwise, a value of + None will be returned as the figure and axes have already been + generated. """ @@ -327,8 +333,12 @@ class Element(object): data = data_new # Generate the plot - fig = plt.figure(**kwargs) - ax = fig.add_subplot(111) + if axis is None: + fig = plt.figure(**kwargs) + ax = fig.add_subplot(111) + else: + fig = None + ax = axis 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 @@ -337,11 +347,12 @@ class Element(object): 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]) ax.set_xlabel('Energy [eV]') if divisor_types: - ax.set_ylabel('Elemental Cross Section') + ax.set_ylabel('Elemental Data') else: ax.set_ylabel('Elemental Cross Section [b]') ax.legend(loc='best') diff --git a/openmc/material.py b/openmc/material.py index ec68913ca..ed54bffc3 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -666,8 +666,7 @@ class Material(object): nuc_densities = [] nuc_density_types = [] for nuclide in nuclides.items(): - nuc, nuc_data = nuclide - nuc, nuc_density, nuc_density_type = nuc_data + nuc, nuc_density, nuc_density_type = nuclide[1] nucs.append(nuc) nuc_densities.append(nuc_density) nuc_density_types.append(nuc_density_type) @@ -713,7 +712,8 @@ class Material(object): return nuclides def plot_xs(self, types, divisor_types=None, temperature=294., - Erange=(1.E-5, 20.E6), cross_sections=None, **kwargs): + axis=None, Erange=(1.E-5, 20.E6), cross_sections=None, + **kwargs): """Creates a figure of continuous-energy macroscopic cross sections for this material @@ -730,6 +730,9 @@ class Material(object): 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. Erange : tuple of floats, optional Energy range (in eV) to plot the cross section within cross_sections : str, optional @@ -740,8 +743,11 @@ class Material(object): Returns ------- - fig : matplotlib.figure.Figure - Matplotlib Figure of the generated macroscopic cross section + fig : matplotlib.figure.Figure or None + If axis is None, then a Matplotlib Figure of the generated + macroscopic cross section will be returned. Otherwise, a value of + None will be returned as the figure and axes have already been + generated. """ @@ -770,8 +776,12 @@ class Material(object): data = data_new # Generate the plot - fig = plt.figure(**kwargs) - ax = fig.add_subplot(111) + if axis is None: + fig = plt.figure(**kwargs) + ax = fig.add_subplot(111) + else: + fig = None + ax = axis 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 @@ -784,7 +794,7 @@ class Material(object): ax.set_xlabel('Energy [eV]') if divisor_types: - ax.set_ylabel('Macroscopic Cross Section') + ax.set_ylabel('Macroscopic Data') else: ax.set_ylabel('Macroscopic Cross Section [1/cm]') ax.legend(loc='best') diff --git a/openmc/nuclide.py b/openmc/nuclide.py index fda44c43a..b549da896 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -97,7 +97,7 @@ class Nuclide(object): self._scattering = scattering - def plot_xs(self, types, divisor_types=None, temperature=294., + def plot_xs(self, types, divisor_types=None, temperature=294., axis=None, Erange=(1.E-5, 20.E6), sab_name=None, cross_sections=None, **kwargs): """Creates a figure of continuous-energy cross sections for this @@ -116,6 +116,9 @@ class Nuclide(object): 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. Erange : tuple of floats Energy range (in eV) to plot the cross section within sab_name : str, optional @@ -129,7 +132,10 @@ class Nuclide(object): Returns ------- fig : matplotlib.figure.Figure - Matplotlib Figure of the generated macroscopic cross section + If axis is None, then a Matplotlib Figure of the generated + macroscopic cross section will be returned. Otherwise, a value of + None will be returned as the figure and axes have already been + generated. """ @@ -158,8 +164,12 @@ class Nuclide(object): data = data_new # Generate the plot - fig = plt.figure(**kwargs) - ax = fig.add_subplot(111) + if axis is None: + fig = plt.figure(**kwargs) + ax = fig.add_subplot(111) + else: + fig = None + ax = axis for i in range(len(data)): to_plot = data[i](E) # Set to loglog or semilogx depending on if we are plotting a data @@ -173,7 +183,7 @@ class Nuclide(object): ax.set_xlabel('Energy [eV]') if divisor_types: - ax.set_ylabel('Microscopic Cross Section') + ax.set_ylabel('Microscopic Data') else: ax.set_ylabel('Microscopic Cross Section [b]') ax.legend(loc='best') @@ -357,8 +367,13 @@ class Nuclide(object): funcs.append(nuc[mt].xs[nucT]) else: funcs.append(nuc[mt].xs[nucT]) - elif mt == 0: + 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)) diff --git a/openmc/plot_data.py b/openmc/plot_data.py index ac59c16ea..ed87f2f00 100644 --- a/openmc/plot_data.py +++ b/openmc/plot_data.py @@ -2,7 +2,12 @@ import numpy as np # Supported keywords for material xs plotting PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', - 'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity'] + 'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity', + 'slowing-down power'] + +# Special MT values +UNITY_MT = -1 +XI_MT = -2 # MTs to combine to generate associated plot_types PLOT_TYPES_MT = {'total': (2, 3,), @@ -31,7 +36,15 @@ PLOT_TYPES_MT = {'total': (2, 3,), 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 183, 184, 190, 194, 196, 198, 199, 200, 875, 891), - 'unity': (0,)} + 'unity': (UNITY_MT,), + 'slowing-down power': (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, XI_MT)} # Operations to use when combining MTs the first np.add is used in reference # to zero PLOT_TYPES_OP = {'total': (np.add,), @@ -41,7 +54,9 @@ PLOT_TYPES_OP = {'total': (np.add,), 'fission': (), 'absorption': (), 'capture': (), 'nu-fission': (), 'nu-scatter': (np.add,) * (len(PLOT_TYPES_MT['nu-scatter']) - 1), - 'unity': ()} + 'unity': (), + 'slowing-down power': + (np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,)} # Whether or not to multiply the reaction by the yield as well PLOT_TYPES_YIELD = {'total': (False, False), @@ -51,7 +66,9 @@ PLOT_TYPES_YIELD = {'total': (False, False), 'fission': (False,), 'absorption': (False,), 'capture': (False,), 'nu-fission': (True,), 'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']), - 'unity': (False,)} + 'unity': (False,), + 'slowing-down power': + (True,) * len(PLOT_TYPES_MT['slowing-down power'])} # Types of plots to plot linearly in y PLOT_TYPES_LINEAR = ['nu-fission / fission', 'nu-scatter / scatter'] From 2a50aff588fed0224699456a9e319bfb298a4610 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 15:20:40 -0500 Subject: [PATCH 14/31] Pre-PR cleanup and addition of ability to read MT in types --- openmc/element.py | 18 ++---------------- openmc/material.py | 19 ++----------------- openmc/nuclide.py | 13 +++++++++---- 3 files changed, 13 insertions(+), 37 deletions(-) diff --git a/openmc/element.py b/openmc/element.py index d90cdf4fc..6fe988195 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -400,20 +400,6 @@ class Element(object): 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) # If cross_sections is None, get the cross sections from the @@ -462,8 +448,8 @@ class Element(object): 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)): + data = np.zeros((len(types), len(unionE))) + for l in range(len(types)): if types[l] == 'unity': data[l, :] = 1. else: diff --git a/openmc/material.py b/openmc/material.py index ed54bffc3..1889a4765 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -834,21 +834,6 @@ class Material(object): # 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() - if self.temperature is not None: T = self.temperature else: @@ -901,8 +886,8 @@ class Material(object): 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)): + data = np.zeros((len(types), len(unionE))) + for l in range(len(types)): if types[l] == 'unity': data[l, :] = 1. else: diff --git a/openmc/nuclide.py b/openmc/nuclide.py index b549da896..439aedf48 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -200,8 +200,10 @@ class Nuclide(object): Parameters ---------- - types : Iterable of values of PLOT_TYPES - The type of cross sections to calculate + 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 @@ -224,7 +226,6 @@ class Nuclide(object): # 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) @@ -239,7 +240,11 @@ class Nuclide(object): ops.append(PLOT_TYPES_OP[line]) else: # Not a built-in type, we have to parse it ourselves - raise NotImplementedError() + cv.check_type('MT in types', line, Integral) + cv.check_greater_than('MT in types', line, 0) + mts.append((line,)) + yields.append((False,)) + ops.append(()) # If cross_sections is None, get the cross sections from the # OPENMC_CROSS_SECTIONS environment variable From d4f4166ba7ec48be7328a565f4fc4bc467557331 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 20:37:03 -0500 Subject: [PATCH 15/31] resolved most of the comments from @samuelshaner --- openmc/data/library.py | 20 ++++++++-- openmc/element.py | 65 ++++++++++++++------------------ openmc/material.py | 85 +++++++++++++++--------------------------- openmc/nuclide.py | 58 +++++++++++++--------------- openmc/plot_data.py | 14 ++++--- 5 files changed, 110 insertions(+), 132 deletions(-) diff --git a/openmc/data/library.py b/openmc/data/library.py index 9f6d93146..3beb25543 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -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: diff --git a/openmc/element.py b/openmc/element.py index 6fe988195..905cdc6c9 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -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 diff --git a/openmc/material.py b/openmc/material.py index 1889a4765..cd80b2089 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -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") diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 439aedf48..c42d4f8d5 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -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 diff --git a/openmc/plot_data.py b/openmc/plot_data.py index ed87f2f00..1efeda00c 100644 --- a/openmc/plot_data.py +++ b/openmc/plot_data.py @@ -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'} From 995ba2669482d79941dee5018c8dd32934df30c0 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 21:58:31 -0500 Subject: [PATCH 16/31] Moved plotting and calculate_xs routines to plotter, condensed code as needed. --- openmc/__init__.py | 2 +- openmc/element.py | 192 -------------- openmc/material.py | 178 ------------- openmc/nuclide.py | 287 --------------------- openmc/plot_data.py | 78 ------ openmc/plotter.py | 604 ++++++++++++++++++++++++++++++++++++++++++++ 6 files changed, 605 insertions(+), 736 deletions(-) delete mode 100644 openmc/plot_data.py create mode 100644 openmc/plotter.py diff --git a/openmc/__init__.py b/openmc/__init__.py index 98de9b70d..0685a80b1 100644 --- a/openmc/__init__.py +++ b/openmc/__init__.py @@ -24,7 +24,7 @@ from openmc.statepoint import * from openmc.summary import * from openmc.particle_restart import * from openmc.mixin import * -from openmc.plot_data import * +from openmc.plotter import * try: from openmc.opencg_compatible import * diff --git a/openmc/element.py b/openmc/element.py index 905cdc6c9..edfbee136 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -1,17 +1,13 @@ from collections import OrderedDict -from numbers import Real import re -import sys import os from six import string_types from xml.etree import ElementTree as ET -import numpy as np import openmc import openmc.checkvalue as cv from openmc.data import NATURAL_ABUNDANCE, atomic_mass -from openmc.plot_data import * class Element(object): @@ -260,191 +256,3 @@ class Element(object): isotopes.append((nuc, percent * abundance, percent_type)) return isotopes - - def plot_xs(self, types, divisor_types=None, temperature=294., axis=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 - - 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. - 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 - 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). - **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 - macroscopic 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 - - 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 = 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)): - if np.sum(data[i, :]) > 0.: - print('max = {:.2E}'.format(np.max(data[i, :]))) - plot_func(E, data[i, :], label=types[i]) - - ax.set_xlabel('Energy [eV]') - if divisor_types: - ax.set_ylabel('Elemental Data') - else: - ax.set_ylabel('Elemental Cross Section [b]') - ax.legend(loc='best') - ax.set_xlim(energy_range) - 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., sab_name=None, - cross_sections=None, enrichment=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. - 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.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 energy_grid - - """ - - # Check types - if cross_sections is not None: - cv.check_type('cross_sections', cross_sections, str) - cv.check_iterable_type('types', types, str) - cv.check_type('temperature', temperature, Real) - - # 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, - cross_sections=cross_sections) - - # For ease of processing split out nuc and nuc_density - nuc_fractions = [nuclide[1] for nuclide in nuclides] - - # Identify the nuclides which have S(a,b) data - sabs = {} - for nuclide in nuclides: - sabs[nuclide[0].name] = None - 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: - sab_tab = sabs[nuclide[0].name] - temp_E, temp_xs = nuclide[0].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 - energy_grid = E[0] - for n in range(1, len(E)): - energy_grid = np.union1d(energy_grid, E[n]) - - # 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 n in range(len(nuclides)): - data[line, :] += nuc_fractions[n] * xs[n][line](energy_grid) - - return energy_grid, data diff --git a/openmc/material.py b/openmc/material.py index cd80b2089..16d1fd147 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -3,7 +3,6 @@ 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 @@ -12,7 +11,6 @@ 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 @@ -697,182 +695,6 @@ class Material(object): return nuclides - def plot_xs(self, types, divisor_types=None, temperature=294., - 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 - - 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. - 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, optional - Energy range (in eV) to plot the cross section within - cross_sections : str, optional - Location of cross_sections.xml file. Default is None. - **kwargs - All keyword arguments are passed to - :func:`matplotlib.pyplot.figure`. - - Returns - ------- - fig : matplotlib.figure.Figure or None - If axis is None, then a Matplotlib Figure of the generated - macroscopic 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 - - 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 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)): - if np.sum(data[i, :]) > 0.: - plot_func(E, data[i, :], label=types[i]) - - ax.set_xlabel('Energy [eV]') - if divisor_types: - ax.set_ylabel('Macroscopic Data') - else: - ax.set_ylabel('Macroscopic Cross Section [1/cm]') - ax.legend(loc='best') - ax.set_xlim(energy_range) - if self.name is not None: - title = 'Macroscopic Cross Section for ' + self.name - ax.set_title(title) - - return fig - - def calculate_xs(self, types, temperature=294., 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. - cross_sections : str, optional - Location of cross_sections.xml file. Default is None. - - Returns - ------- - 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 energy_grid - - """ - - # Check types - if cross_sections is not None: - cv.check_type('cross_sections', cross_sections, str) - if self.temperature is not None: - T = self.temperature - else: - cv.check_type('temperature', temperature, Real) - T = temperature - - # 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) - - # For ease of processing split out nuc and nuc_density - nuc_densities = [nuclide[1][1] for nuclide in nuclides.items()] - - # Identify the nuclides which have S(a,b) data - sabs = {} - for nuclide in nuclides.items(): - sabs[nuclide[0].name] = None - for sab_name in self._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'] - - # Now we can create the data sets to be plotted - xs = [] - E = [] - for nuclide in nuclides.items(): - sab_tab = sabs[nuclide[0].name] - temp_E, temp_xs = nuclide[0].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 - energy_grid = E[0] - for n in range(1, len(E)): - energy_grid = np.union1d(energy_grid, E[n]) - - # 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 n in range(len(nuclides)): - data[line, :] += nuc_densities[n] * xs[n][line](energy_grid) - - return energy_grid, data - def _get_nuclide_xml(self, nuclide, distrib=False): xml_element = ET.Element("nuclide") xml_element.set("name", nuclide[0].name) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index c42d4f8d5..03062ee0e 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -1,13 +1,8 @@ -from numbers import Integral, Real -import sys import warnings -import os from six import string_types import openmc.checkvalue as cv -from openmc.plot_data import * -import openmc.data class Nuclide(object): @@ -96,285 +91,3 @@ class Nuclide(object): raise ValueError(msg) self._scattering = scattering - - def plot_xs(self, types, divisor_types=None, temperature=294., axis=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 - - 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. - 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 - 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. - **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 - macroscopic 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 - - 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 line in range(len(types)): - data_new.append(openmc.data.Combination([data[line], - data_div[line]], - [np.divide])) - 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)): - to_plot = data[i](E) - if np.sum(to_plot) > 0.: - plot_func(E, to_plot, label=types[i]) - - ax.set_xlabel('Energy [eV]') - if divisor_types: - ax.set_ylabel('Microscopic Data') - else: - ax.set_ylabel('Microscopic Cross Section [b]') - ax.legend(loc='best') - ax.set_xlim(energy_range) - 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 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.array - 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 - if cross_sections is not None: - cv.check_type('cross_sections', cross_sections, 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 - cv.check_type('MT in types', line, Integral) - cv.check_greater_than('MT in types', line, 0) - mts.append((line,)) - yields.append((False,)) - ops.append(()) - - # 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) - 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(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) - 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, 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 == 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(nuclide[0] + " not in library") - - return energy_grid, xs diff --git a/openmc/plot_data.py b/openmc/plot_data.py deleted file mode 100644 index 1efeda00c..000000000 --- a/openmc/plot_data.py +++ /dev/null @@ -1,78 +0,0 @@ -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', 'damage'] - -# Special MT values -UNITY_MT = -1 -XI_MT = -2 - -# MTs to combine to generate associated plot_types -PLOT_TYPES_MT = {'total': (2, 3,), - '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), - 'elastic': (2,), - 'inelastic': (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), - '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': (UNITY_MT,), - 'slowing-down power': (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, 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': ()} - -# Whether or not to multiply the reaction by the yield as well -PLOT_TYPES_YIELD = {'total': (False, False), - 'scatter': (False,) * len(PLOT_TYPES_MT['scatter']), - 'elastic': (False,), - 'inelastic': (False,) * len(PLOT_TYPES_MT['inelastic']), - 'fission': (False,), 'absorption': (False,), - 'capture': (False,), 'nu-fission': (True,), - 'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']), - 'unity': (False,), - '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', - 'nu-fission / absorption', 'fission / absorption'} diff --git a/openmc/plotter.py b/openmc/plotter.py new file mode 100644 index 000000000..0e604c8c1 --- /dev/null +++ b/openmc/plotter.py @@ -0,0 +1,604 @@ +from numbers import Integral, Real + +import numpy as np +from matplotlib import pyplot as plt + +import openmc.checkvalue as cv +import openmc.data + +# Supported keywords for material 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 +PLOT_TYPES_MT = {'total': (2, 3,), + '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), + 'elastic': (2,), + 'inelastic': (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), + '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': (UNITY_MT,), + 'slowing-down power': (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, 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': ()} + +# Whether or not to multiply the reaction by the yield as well +PLOT_TYPES_YIELD = {'total': (False, False), + 'scatter': (False,) * len(PLOT_TYPES_MT['scatter']), + 'elastic': (False,), + 'inelastic': (False,) * len(PLOT_TYPES_MT['inelastic']), + 'fission': (False,), 'absorption': (False,), + 'capture': (False,), 'nu-fission': (True,), + 'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']), + 'unity': (False,), + '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', + 'nu-fission / absorption', 'fission / absorption'} + + +def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, + energy_range=(1.E-5, 20.E6), 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. + 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; 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. + + """ + + 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) + + 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) + + # 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) + if data_type == 'nuclide': + data_new = [] + else: + data_new = np.zeros((len(types), len(E))) + + for line in range(len(types)): + if data_type == 'nuclide': + data_new.append(openmc.data.Combination([data[line], + data_div[line]], + [np.divide])) + else: + 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)): + if data_type == 'nuclide': + to_plot = data[i](E) + else: + to_plot = data[i, :] + if np.sum(to_plot) > 0.: + plot_func(E, to_plot, label=types[i]) + + ax.set_xlabel('Energy [eV]') + if divisor_types: + ax.set_ylabel('Data') + else: + ax.set_ylabel('Cross Section [b]') + 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) + + 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.array + 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, str) + if enrichment: + cv.check_type('enrichment', enrichment, Real) + + if isinstance(this, openmc.Nuclide): + energy_grid, data = _calculate_xs_nuclide(this, types, temperature, + sab_name, cross_sections) + elif isinstance(this, openmc.Element): + energy_grid, data = _calculate_xs_element(this, types, temperature, + sab_name, cross_sections, + enrichment) + elif isinstance(this, openmc.Material): + energy_grid, data = _calculate_xs_material(this, types, temperature, + cross_sections) + else: + raise TypeError("Invalid type") + + return energy_grid, data + + +def _calculate_xs_element(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 + 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. + 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.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 energy_grid + + """ + + # Load the library + library = openmc.data.DataLibrary.from_xml(cross_sections) + + # Expand elements in to nuclides with atomic densities + nuclides = this.expand(100., 'ao', enrichment=enrichment, + cross_sections=cross_sections) + + # For ease of processing split out nuc and nuc_density + nuc_fractions = [nuclide[1] for nuclide in nuclides] + + # Identify the nuclides which have S(a,b) data + sabs = {} + for nuclide in nuclides: + sabs[nuclide[0].name] = None + 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: + sab_tab = sabs[nuclide[0].name] + temp_E, temp_xs = calculate_xs(nuclide[0], 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 + energy_grid = E[0] + for n in range(1, len(E)): + energy_grid = np.union1d(energy_grid, E[n]) + + # 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 n in range(len(nuclides)): + data[line, :] += nuc_fractions[n] * xs[n][line](energy_grid) + + 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.array + 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 + + """ + + # 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 + cv.check_type('MT in types', line, Integral) + cv.check_greater_than('MT in types', line, 0) + mts.append((line,)) + yields.append((False,)) + ops.append(()) + + # 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 = 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) + 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, 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 == 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_material(this, types, temperature=294., cross_sections=None): + """Calculates continuous-energy macroscopic cross sections of a + requested type + + Parameters + ---------- + this : openmc.Material + 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. + cross_sections : str, optional + Location of cross_sections.xml file. Default is None. + + Returns + ------- + 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 energy_grid + + """ + + if this.temperature is not None: + T = this.temperature + else: + T = temperature + + # Load the library + library = openmc.data.DataLibrary.from_xml(cross_sections) + + # Expand elements in to nuclides with atomic densities + nuclides = this.get_nuclide_atom_densities(cross_sections) + + # For ease of processing split out nuc and nuc_density + nuc_densities = [nuclide[1][1] for nuclide in nuclides.items()] + + # Identify the nuclides which have S(a,b) data + sabs = {} + for nuclide in nuclides.items(): + sabs[nuclide[0].name] = None + 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'] + + # Now we can create the data sets to be plotted + xs = [] + E = [] + for nuclide in nuclides.items(): + sab_tab = sabs[nuclide[0].name] + temp_E, temp_xs = calculate_xs(nuclide[0], 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 + energy_grid = E[0] + for n in range(1, len(E)): + energy_grid = np.union1d(energy_grid, E[n]) + + # 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 n in range(len(nuclides)): + data[line, :] += nuc_densities[n] * xs[n][line](energy_grid) + + return energy_grid, data From ba5a6714a23e9fd6f1444a9b60e50f4c2537461a Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 22:13:14 -0500 Subject: [PATCH 17/31] Forgot enrichment --- openmc/plotter.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 0e604c8c1..4201abf9d 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -141,13 +141,14 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, else: raise TypeError("Invalid type for plotting") - E, data = calculate_xs(this, types, temperature, sab_name, cross_sections) + 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) + 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 From 00268a3c55b3070d55d6ec4ce253b12326cde647 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 11 Nov 2016 22:18:12 -0500 Subject: [PATCH 18/31] Fixed Travis failure: moved matplotlib import to plot_xs instead of at the module level --- openmc/plotter.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 4201abf9d..09b7b1bf8 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -1,7 +1,6 @@ from numbers import Integral, Real import numpy as np -from matplotlib import pyplot as plt import openmc.checkvalue as cv import openmc.data @@ -132,6 +131,8 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, """ + from matplotlib import pyplot as plt + if isinstance(this, openmc.Nuclide): data_type = 'nuclide' elif isinstance(this, openmc.Element): From 51babd69782c7282e92420dd5279a99a49f55b92 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 13 Nov 2016 21:54:12 -0500 Subject: [PATCH 19/31] Minor ylabel changes --- openmc/plotter.py | 17 ++++++++++++++--- 1 file changed, 14 insertions(+), 3 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 09b7b1bf8..31fb57f5e 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -61,7 +61,7 @@ PLOT_TYPES_OP = {'total': (np.add,), 'capture': (), 'nu-fission': (), 'nu-scatter': (np.add,) * (len(PLOT_TYPES_MT['nu-scatter']) - 1), 'unity': (), - 'slowing-down power': + 'slowing-down power': (np.add,) * (len(PLOT_TYPES_MT['slowing-down power']) - 2) + (np.multiply,), 'damage': ()} @@ -197,9 +197,20 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, ax.set_xlabel('Energy [eV]') if divisor_types: - ax.set_ylabel('Data') + if data_type == 'nuclide': + ylabel = 'Nuclidic Microscopic Data' + elif data_type == 'element': + ylabel = 'Elemental Microscopic Data' + elif data_type == 'material': + ylabel = 'Macroscopic Data' else: - ax.set_ylabel('Cross Section [b]') + 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') ax.set_xlim(energy_range) if this.name is not None: From a428ccb2315ba6b362475890ef8f838f5d0077f9 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 13 Nov 2016 21:54:36 -0500 Subject: [PATCH 20/31] And one minor editorial :-( --- openmc/plotter.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 31fb57f5e..b8eda28e6 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -5,7 +5,7 @@ import numpy as np import openmc.checkvalue as cv import openmc.data -# Supported keywords for material xs plotting +# Supported keywords for xs plotting PLOT_TYPES = ['total', 'scatter', 'elastic', 'inelastic', 'fission', 'absorption', 'capture', 'nu-fission', 'nu-scatter', 'unity', 'slowing-down power', 'damage'] From 625903e968b40ffe84445962bc1fe54cfe0a2c35 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 14 Nov 2016 20:48:13 -0500 Subject: [PATCH 21/31] One more minor typo --- openmc/data/library.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/data/library.py b/openmc/data/library.py index 3beb25543..34cd380a5 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -114,7 +114,7 @@ class DataLibrary(EqualityMixin): data = cls() - # If cross_sections is None, get the cross sections from the + # 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') From b8d9fe461f22e59437af10023a20eaf8cd92bbe4 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 16 Nov 2016 21:47:07 -0500 Subject: [PATCH 22/31] Resolving @paulromano comments --- openmc/data/data.py | 6 ++++ openmc/material.py | 17 ++++----- openmc/plotter.py | 85 +++++++++++++++------------------------------ 3 files changed, 40 insertions(+), 68 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index a3e715812..92be7cade 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -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 diff --git a/openmc/material.py b/openmc/material.py index 16d1fd147..ece75be19 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -664,17 +664,12 @@ class Material(object): sum_percent = 0. awrs = [] - n = -1 - for nuclide in nuclides.items(): - n += 1 - lib = library.get_by_material(nuclide[0]) - if lib is not None: - nuc = openmc.data.IncidentNeutron.from_hdf5(lib['path']) - awrs.append(nuc.atomic_weight_ratio) - if not percent_in_atom: - nuc_densities[n] = -nuc_densities[n] / awrs[-1] + 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] + " not in library") + raise ValueError(nuclide[0] + " is invalid") # Now that we have the awr, lets finish calculating densities sum_percent = np.sum(nuc_densities) @@ -686,7 +681,7 @@ class Material(object): sum_percent += x * awrs[n] sum_percent = 1. / sum_percent density = -density * sum_percent * \ - sc.Avogadro / sc.value('neutron mass in u') * 1.E-24 + openmc.data.AVOGADRO / openmc.data.NEUTRON_MASS * 1.E-24 nuc_densities = density * nuc_densities nuclides = OrderedDict() diff --git a/openmc/plotter.py b/openmc/plotter.py index b8eda28e6..f0275dda6 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -15,42 +15,18 @@ UNITY_MT = -1 XI_MT = -2 # MTs to combine to generate associated plot_types -PLOT_TYPES_MT = {'total': (2, 3,), - '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), - 'elastic': (2,), - 'inelastic': (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), - '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': (UNITY_MT,), - 'slowing-down power': (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, XI_MT), - 'damage': (444,)} +_INELASTIC = [mt for mt in openmc.data.SUM_RULES[3] if mt not in [5, 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,), @@ -84,8 +60,7 @@ PLOT_TYPES_LINEAR = {'nu-fission / fission', 'nu-scatter / scatter', def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, - energy_range=(1.E-5, 20.E6), sab_name=None, cross_sections=None, - enrichment=None, **kwargs): + sab_name=None, cross_sections=None, enrichment=None, **kwargs): """Creates a figure of continuous-energy cross sections for this item Parameters @@ -106,8 +81,6 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, 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 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 @@ -212,10 +185,10 @@ def plot_xs(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) + # Set to the most likely expected range + ax.set_xlim((1.E-5, 20.E6)) 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 @@ -246,7 +219,7 @@ def calculate_xs(this, types, temperature=294., sab_name=None, Returns ------- - energy_grid : numpy.array + 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 @@ -302,7 +275,7 @@ def _calculate_xs_element(this, types, temperature=294., sab_name=None, Returns ------- - energy_grid : numpy.array + energy_grid : numpy.ndarray Energies at which cross sections are calculated, in units of eV data : numpy.ndarray Macroscopic cross sections calculated at the energy grid described @@ -381,7 +354,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, Returns ------- - energy_grid : numpy.array + 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 @@ -405,6 +378,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, mts.append((line,)) yields.append((False,)) ops.append(()) + print(mts) # Load the library library = openmc.data.DataLibrary.from_xml(cross_sections) @@ -427,7 +401,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, 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 + min_delta = float('inf') closest_t = -1 for t in data_Ts: if abs(data_Ts[t] - T) < min_delta: @@ -496,7 +470,6 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, 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': @@ -504,9 +477,8 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, [nuc[mt].xs[nucT], prod.yield_], [np.multiply]) funcs.append(func) - found_it = True break - if not found_it: + else: for prod in nuc[mt].products: if prod.particle == 'neutron' and \ prod.emission_mode == 'prompt': @@ -514,11 +486,10 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, [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: + # Assume the yield is 1 + funcs.append(nuc[mt].xs[nucT]) else: funcs.append(nuc[mt].xs[nucT]) elif mt == UNITY_MT: @@ -557,7 +528,7 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): Returns ------- - energy_grid : numpy.array + energy_grid : numpy.ndarray Energies at which cross sections are calculated, in units of eV data : numpy.ndarray Macroscopic cross sections calculated at the energy grid described @@ -602,8 +573,8 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): # Condense the data for every nuclide # First create a union energy grid energy_grid = E[0] - for n in range(1, len(E)): - energy_grid = np.union1d(energy_grid, E[n]) + 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))) From 1a62b59dc3f6ed2c4a3f746553d4a2eb103ecb8b Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Thu, 17 Nov 2016 17:17:38 -0500 Subject: [PATCH 23/31] Resolved @samuelshaner comments --- openmc/material.py | 5 ----- openmc/plotter.py | 4 ++-- 2 files changed, 2 insertions(+), 7 deletions(-) diff --git a/openmc/material.py b/openmc/material.py index ece75be19..84d790feb 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -623,11 +623,6 @@ class Material(object): """ - import scipy.constants as sc - - # Load the library - library = openmc.data.DataLibrary.from_xml(cross_sections) - # Expand elements in to nuclides nuclides = self.get_nuclide_densities() diff --git a/openmc/plotter.py b/openmc/plotter.py index f0275dda6..626450456 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -165,6 +165,7 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, to_plot = data[i](E) else: to_plot = data[i, :] + to_plot = np.nan_to_num(to_plot) if np.sum(to_plot) > 0.: plot_func(E, to_plot, label=types[i]) @@ -287,7 +288,7 @@ def _calculate_xs_element(this, types, temperature=294., sab_name=None, library = openmc.data.DataLibrary.from_xml(cross_sections) # Expand elements in to nuclides with atomic densities - nuclides = this.expand(100., 'ao', enrichment=enrichment, + nuclides = this.expand(1., 'ao', enrichment=enrichment, cross_sections=cross_sections) # For ease of processing split out nuc and nuc_density @@ -378,7 +379,6 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, mts.append((line,)) yields.append((False,)) ops.append(()) - print(mts) # Load the library library = openmc.data.DataLibrary.from_xml(cross_sections) From b4c556db828f23e52f8ff1d4a3c188dfc2e8b58d Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Fri, 18 Nov 2016 19:55:24 -0500 Subject: [PATCH 24/31] incorporated MT=5, simplified yield usage in the plotter --- openmc/data/library.py | 3 +- openmc/material.py | 7 +---- openmc/plotter.py | 65 +++++++++++++++++++++--------------------- 3 files changed, 36 insertions(+), 39 deletions(-) diff --git a/openmc/data/library.py b/openmc/data/library.py index 34cd380a5..c179f78f8 100644 --- a/openmc/data/library.py +++ b/openmc/data/library.py @@ -1,5 +1,6 @@ import os import xml.etree.ElementTree as ET +from six import string_types import h5py @@ -124,7 +125,7 @@ class DataLibrary(EqualityMixin): raise ValueError("Either path or OPENMC_CROSS_SECTIONS " "environmental variable must be set") - check_type('path', path, str) + check_type('path', path, string_types) tree = ET.parse(path) root = tree.getroot() diff --git a/openmc/material.py b/openmc/material.py index 84d790feb..08b001e47 100644 --- a/openmc/material.py +++ b/openmc/material.py @@ -606,15 +606,10 @@ class Material(object): return nuclides - def get_nuclide_atom_densities(self, cross_sections=None): + def get_nuclide_atom_densities(self): """Returns all nuclides in the material and their atomic densities in units of atom/b-cm - Parameters - ---------- - cross_sections : str, optional - Location of cross_sections.xml file. Default is None. - Returns ------- nuclides : dict diff --git a/openmc/plotter.py b/openmc/plotter.py index 626450456..e20363ca4 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -1,4 +1,5 @@ from numbers import Integral, Real +from six import string_types import numpy as np @@ -15,7 +16,7 @@ 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 not in [5, 27]] +_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], @@ -41,19 +42,6 @@ PLOT_TYPES_OP = {'total': (np.add,), (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), - 'scatter': (False,) * len(PLOT_TYPES_MT['scatter']), - 'elastic': (False,), - 'inelastic': (False,) * len(PLOT_TYPES_MT['inelastic']), - 'fission': (False,), 'absorption': (False,), - 'capture': (False,), 'nu-fission': (True,), - 'nu-scatter': (True,) * len(PLOT_TYPES_MT['nu-scatter']), - 'unity': (False,), - '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', 'nu-fission / absorption', 'fission / absorption'} @@ -230,7 +218,7 @@ def calculate_xs(this, types, temperature=294., sab_name=None, # Check types cv.check_type('temperature', temperature, Real) if sab_name: - cv.check_type('sab_name', sab_name, str) + cv.check_type('sab_name', sab_name, string_types) if enrichment: cv.check_type('enrichment', enrichment, Real) @@ -370,15 +358,18 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, for line in types: if line in PLOT_TYPES: mts.append(PLOT_TYPES_MT[line]) - yields.append(PLOT_TYPES_YIELD[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,)) - yields.append((False,)) ops.append(()) + yields.append(False) # Load the library library = openmc.data.DataLibrary.from_xml(cross_sections) @@ -456,7 +447,7 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, # Get the reaction xs data from the nuclide funcs = [] op = ops[i] - for mt, yield_check in zip(mt_set, yields[i]): + for mt in mt_set: if mt == 2: if sab_name: # Then we need to do a piece-wise function of @@ -468,28 +459,38 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, funcs.append(pw_funcs) else: funcs.append(nuc[mt].xs[nucT]) + elif mt == 5 and mt in nuc: + # Only consider the (n,misc) products if neutrons are + # included in the outgoing channel since (n,misc) is only + # explicitly needed for scatter cross sections. + for prod in nuc[mt].products: + if prod.particle == 'neutron' and \ + prod.emission_mode in ('total', 'prompt'): + if yields[i]: + func = openmc.data.Combination( + [nuc[mt].xs[nucT], prod.yield_], + [np.multiply]) + else: + func = nuc[mt].xs[nucT] + + funcs.append(func) + break + else: + funcs.append(lambda x: 0.) + elif mt in nuc: - if yield_check: + if yields[i]: for prod in nuc[mt].products: if prod.particle == 'neutron' and \ - prod.emission_mode == 'total': + prod.emission_mode in ('total', 'prompt'): func = openmc.data.Combination( [nuc[mt].xs[nucT], prod.yield_], [np.multiply]) funcs.append(func) break else: - 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) - break - else: - # Assume the yield is 1 - funcs.append(nuc[mt].xs[nucT]) + # Assume the yield is 1 + funcs.append(nuc[mt].xs[nucT]) else: funcs.append(nuc[mt].xs[nucT]) elif mt == UNITY_MT: @@ -545,7 +546,7 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): library = openmc.data.DataLibrary.from_xml(cross_sections) # Expand elements in to nuclides with atomic densities - nuclides = this.get_nuclide_atom_densities(cross_sections) + nuclides = this.get_nuclide_atom_densities() # For ease of processing split out nuc and nuc_density nuc_densities = [nuclide[1][1] for nuclide in nuclides.items()] From e764a0852fcc809de030eb84f82fcdc7ce10dbd7 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 19 Nov 2016 21:02:50 -0500 Subject: [PATCH 25/31] Fixed the stupid for n in range... --- openmc/plotter.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index e20363ca4..72d4a43fa 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -304,8 +304,8 @@ def _calculate_xs_element(this, types, temperature=294., sab_name=None, # Condense the data for every nuclide # First create a union energy grid energy_grid = E[0] - for n in range(1, len(E)): - energy_grid = np.union1d(energy_grid, E[n]) + 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))) From 9afb0c8d68c01b513fbe5c1d9240e290f67104fe Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 19 Nov 2016 22:52:30 -0500 Subject: [PATCH 26/31] Made the _calculate_xs routines consistent in what they return to give us a little simplification --- openmc/plotter.py | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 72d4a43fa..97add26ff 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -149,13 +149,9 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, plot_func = ax.loglog # Plot the data for i in range(len(data)): - if data_type == 'nuclide': - to_plot = data[i](E) - else: - to_plot = data[i, :] - to_plot = np.nan_to_num(to_plot) - if np.sum(to_plot) > 0.: - plot_func(E, to_plot, 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: @@ -223,11 +219,17 @@ def calculate_xs(this, types, temperature=294., sab_name=None, cv.check_type('enrichment', enrichment, Real) if isinstance(this, openmc.Nuclide): - energy_grid, data = _calculate_xs_nuclide(this, types, temperature, - sab_name, cross_sections) + 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_element(this, types, temperature, - sab_name, cross_sections, + cross_sections, sab_name, enrichment) elif isinstance(this, openmc.Material): energy_grid, data = _calculate_xs_material(this, types, temperature, @@ -299,7 +301,10 @@ def _calculate_xs_element(this, types, temperature=294., sab_name=None, temp_E, temp_xs = calculate_xs(nuclide[0], types, temperature, sab_tab, cross_sections) E.append(temp_E) - xs.append(temp_xs) + # Since the energy grids are different, store the cross sections as + # a tabulated function so they can be calculated on any grid needed. + xs.append([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 @@ -345,9 +350,8 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, ------- 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 + data : Iterable of Callable + Requested cross section functions """ @@ -569,7 +573,10 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): temp_E, temp_xs = calculate_xs(nuclide[0], types, T, sab_tab, cross_sections) E.append(temp_E) - xs.append(temp_xs) + # Since the energy grids are different, store the cross sections as + # a tabulated function so they can be calculated on any grid needed. + xs.append([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 From 24798a0ecd6e51ad2e11f8cd5839ef2847282e3f Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 19 Nov 2016 23:29:59 -0500 Subject: [PATCH 27/31] Combined the elemental and material calculate_xs routines --- openmc/plotter.py | 178 ++++++++++++++++------------------------------ 1 file changed, 63 insertions(+), 115 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 97add26ff..cf2c97b92 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -228,11 +228,11 @@ 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_element(this, types, temperature, - cross_sections, sab_name, - enrichment) + 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_material(this, types, temperature, + energy_grid, data = _calculate_xs_elem_mat(this, types, temperature, cross_sections) else: raise TypeError("Invalid type") @@ -240,90 +240,6 @@ def calculate_xs(this, types, temperature=294., sab_name=None, return energy_grid, data -def _calculate_xs_element(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 - 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. - 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 - Macroscopic cross sections calculated at the energy grid described - by energy_grid - - """ - - # Load the library - library = openmc.data.DataLibrary.from_xml(cross_sections) - - # 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 nuc and nuc_density - nuc_fractions = [nuclide[1] for nuclide in nuclides] - - # Identify the nuclides which have S(a,b) data - sabs = {} - for nuclide in nuclides: - sabs[nuclide[0].name] = None - 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: - sab_tab = sabs[nuclide[0].name] - temp_E, temp_xs = calculate_xs(nuclide[0], types, temperature, 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.append([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 n in range(len(nuclides)): - data[line, :] += nuc_fractions[n] * xs[n][line](energy_grid) - - 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 @@ -513,14 +429,14 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, return energy_grid, xs -def _calculate_xs_material(this, types, temperature=294., cross_sections=None): - """Calculates continuous-energy macroscopic cross sections of a - requested type +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 - Material object to source data from + 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 @@ -530,53 +446,83 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): 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 - Macroscopic cross sections calculated at the energy grid described - by energy_grid + Cross sections calculated at the energy grid described by energy_grid """ - if this.temperature is not None: - T = this.temperature + 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) - # Expand elements in to nuclides with atomic densities - nuclides = this.get_nuclide_atom_densities() + 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 nuc and nuc_density - nuc_densities = [nuclide[1][1] for nuclide in nuclides.items()] + # 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 + 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 + 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].name] = None - 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'] + 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 = [] + xs = {} E = [] + # for nuclide in nuclides: for nuclide in nuclides.items(): - sab_tab = sabs[nuclide[0].name] - temp_E, temp_xs = calculate_xs(nuclide[0], types, T, sab_tab, - cross_sections) + 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.append([openmc.data.Tabulated1D(temp_E, temp_xs[line]) - for line in range(len(types))]) + 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 @@ -590,7 +536,9 @@ def _calculate_xs_material(this, types, temperature=294., cross_sections=None): if types[line] == 'unity': data[line, :] = 1. else: - for n in range(len(nuclides)): - data[line, :] += nuc_densities[n] * xs[n][line](energy_grid) + for nuclide in nuclides.items(): + name = nuclide[0] + data[line, :] += (nuc_fractions[name] * + xs[name][line](energy_grid)) return energy_grid, data From 27c07b02d63b1f5b6bd0b7d4093af6c50055db6e Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sat, 19 Nov 2016 23:32:49 -0500 Subject: [PATCH 28/31] cleaned up some comments --- openmc/plotter.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index cf2c97b92..45783d6cf 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -476,21 +476,23 @@ def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None, 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 @@ -512,7 +514,6 @@ def _calculate_xs_elem_mat(this, types, temperature=294., cross_sections=None, # Now we can create the data sets to be plotted xs = {} E = [] - # for nuclide in nuclides: for nuclide in nuclides.items(): name = nuclide[0] nuc = nuclide[1] From 14e3684a6a4013cb27900d474a60a18424bf0ed8 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 20 Nov 2016 06:05:36 -0500 Subject: [PATCH 29/31] Fixed issue with plotting nuclides with divisor types --- openmc/plotter.py | 16 ++++------------ 1 file changed, 4 insertions(+), 12 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index 45783d6cf..d6c9527bd 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -116,20 +116,12 @@ def plot_xs(this, types, divisor_types=None, temperature=294., axis=None, # grid, and then do the actual division Enum = E[:] E = np.union1d(Enum, Ediv) - if data_type == 'nuclide': - data_new = [] - else: - data_new = np.zeros((len(types), len(E))) + data_new = np.zeros((len(types), len(E))) for line in range(len(types)): - if data_type == 'nuclide': - data_new.append(openmc.data.Combination([data[line], - data_div[line]], - [np.divide])) - else: - data_new[line, :] = \ - np.divide(np.interp(E, Enum, data[line, :]), - np.interp(E, Ediv, data_div[line, :])) + 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 From 20ce145e366a50240fbec7e552f6fb60ed46c3d7 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 21 Nov 2016 20:30:54 -0500 Subject: [PATCH 30/31] Read in total_nu as a derived product for MT=18 if it only exists as summed, and more correctly treated yields --- openmc/data/neutron.py | 3 +++ openmc/plotter.py | 50 +++++++++++++++++++++++------------------- 2 files changed, 30 insertions(+), 23 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index f478aafc7..581adc417 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -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]) diff --git a/openmc/plotter.py b/openmc/plotter.py index d6c9527bd..f576dd8a6 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -1,5 +1,6 @@ from numbers import Integral, Real from six import string_types +from itertools import chain import numpy as np @@ -371,38 +372,41 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, funcs.append(pw_funcs) else: funcs.append(nuc[mt].xs[nucT]) - elif mt == 5 and mt in nuc: - # Only consider the (n,misc) products if neutrons are - # included in the outgoing channel since (n,misc) is only - # explicitly needed for scatter cross sections. - for prod in nuc[mt].products: - if prod.particle == 'neutron' and \ - prod.emission_mode in ('total', 'prompt'): - if yields[i]: - func = openmc.data.Combination( - [nuc[mt].xs[nucT], prod.yield_], - [np.multiply]) - else: - func = nuc[mt].xs[nucT] - - funcs.append(func) - break - else: - funcs.append(lambda x: 0.) - elif mt in nuc: if yields[i]: - for prod in nuc[mt].products: + for prod in chain(nuc[mt].products, + nuc[mt].derived_products): if prod.particle == 'neutron' and \ - prod.emission_mode in ('total', 'prompt'): + prod.emission_mode == 'total': func = openmc.data.Combination( [nuc[mt].xs[nucT], prod.yield_], [np.multiply]) funcs.append(func) break else: - # Assume the yield is 1 - funcs.append(nuc[mt].xs[nucT]) + # Total doesn't exist so we have to create from + # prompt and delayed + 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( + [nuc[mt].xs[nucT], prod.yield_, + func], [np.multiply, np.add]) + else: + func = openmc.data.Combination( + [nuc[mt].xs[nucT], prod.yield_], + [np.multiply]) + if func: + funcs.append(func) + 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: From 42baabd155152d6d343b6252baf86ff9aa054784 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Mon, 21 Nov 2016 20:35:06 -0500 Subject: [PATCH 31/31] Whoops, improperly weighted the data --- openmc/plotter.py | 15 ++++++++------- 1 file changed, 8 insertions(+), 7 deletions(-) diff --git a/openmc/plotter.py b/openmc/plotter.py index f576dd8a6..aa4317d30 100644 --- a/openmc/plotter.py +++ b/openmc/plotter.py @@ -374,6 +374,8 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, 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 \ @@ -385,7 +387,8 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, break else: # Total doesn't exist so we have to create from - # prompt and delayed + # prompt and delayed. This is used for scatter + # multiplication. func = None for prod in chain(nuc[mt].products, nuc[mt].derived_products): @@ -393,14 +396,12 @@ def _calculate_xs_nuclide(this, types, temperature=294., sab_name=None, prod.emission_mode != 'total': if func: func = openmc.data.Combination( - [nuc[mt].xs[nucT], prod.yield_, - func], [np.multiply, np.add]) + [prod.yield_, func], [np.add]) else: - func = openmc.data.Combination( - [nuc[mt].xs[nucT], prod.yield_], - [np.multiply]) + func = prod.yield_ if func: - funcs.append(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