From 76991a223f62c666b3f7fefcf1e9a94f61d74af3 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Wed, 9 Nov 2016 18:10:41 -0500 Subject: [PATCH] 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 48f43be91d..98de9b70d4 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 d56c3b3507..90bde474d2 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 37dbe0d2da..c357da7218 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 55afd4e9c3..31b658a1cf 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 0000000000..03e929b52a --- /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,)}