From 68625ca314a722bf02e994b9aba5f4d74c9a412e Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 25 Aug 2016 11:17:37 -0500 Subject: [PATCH] Improve storage of temperature-dependent cross sections. - Got rid of XS class - Renamed Reaction.T_data to Reaction.xs - Got rid of Reaction.add_temperature() - Simplified IncidentNeutron.add_temperature_from_ace() --- openmc/data/neutron.py | 118 +++++++++------------------ openmc/data/reaction.py | 172 +++++++++++++++------------------------- 2 files changed, 102 insertions(+), 188 deletions(-) diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index dbf2da23c..6020f0962 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -1,6 +1,7 @@ from __future__ import division, unicode_literals import sys from collections import OrderedDict, Iterable, Mapping +from itertools import chain from numbers import Integral, Real from warnings import warn @@ -160,8 +161,6 @@ class IncidentNeutron(EqualityMixin): self.metastable = metastable self.atomic_weight_ratio = atomic_weight_ratio self.kTs = kTs - self.temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" - for kT in kTs] self.energy = {} self._fission_energy = None self.reactions = OrderedDict() @@ -221,6 +220,10 @@ class IncidentNeutron(EqualityMixin): def urr(self): return self._urr + @property + def temperatures(self): + return ["{}K".format(int(round(kT / K_BOLTZMANN))) for kT in self.kTs] + @name.setter def name(self, name): cv.check_type('name', name, basestring) @@ -296,66 +299,31 @@ class IncidentNeutron(EqualityMixin): """ - if isinstance(ace_or_filename, Table): - ace = ace_or_filename - else: - ace = get_table(ace_or_filename) + data = IncidentNeutron.from_ace(ace_or_filename, metastable_scheme) - # Obtain the information needed to check if this ACE file is for the - # same nuclide. - zaid, xs = ace.name.split('.') - name, element, Z, mass_number, metastable = \ - _get_metadata(int(zaid), metastable_scheme) + # Check if temprature already exists + strT = data.temperatures[0] + if strT in self.temperatures: + warn('Cross sections at T={} already exist.'.format(strT)) + return - # If this ACE data matches the data within self then get the data - if ace.temperature not in self.kTs: - if name == self.name: - # Add temperature and kTs - strT = str(int(round(ace.temperature / K_BOLTZMANN))) + "K" - self.temperatures.append(strT) - self.kTs.append(ace.temperature) - # Read energy grid - n_energy = ace.nxs[3] - energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy] - self.energy[strT] = energy - total_xs = \ - Tabulated1D(energy, ace.xss[ace.jxs[1] + - n_energy:ace.jxs[1] + - 2 * n_energy]) - abs_xs = Tabulated1D(energy, ace.xss[ace.jxs[1] + 2 * - n_energy:ace.jxs[1] + - 3 * n_energy]) + # Check that name matches + if data.name != self.name: + raise ValueError('Data provided for an incorrect nuclide.') - self.summed_reactions[1].add_temperature(strT, 0, total_xs) - if 27 in self.summed_reactions: - self.summed_reactions[27].add_temperature(strT, 0, abs_xs) + # Add temperature + self.kTs += data.kTs - # Read each reaction and get the xs data out of it - n_reaction = ace.nxs[4] + 1 - for i in range(n_reaction): - rx = Reaction.from_ace(ace, i) + # Add energy grid + self.energy[strT] = data.energy[strT] - xsdata = list(rx.T_data.values())[0] - self.reactions[rx.mt].add_temperature(strT, - xsdata.threshold_idx, - xsdata.xs) - - # Obtain data for the summed photon reactions - for mt in self.summed_reactions: - if mt not in [1, 27]: - # Create summed appropriate cross section - mts = self.get_reaction_components(mt) - xsdata = Sum([self.reactions[mt_i].T_data[strT].xs - for mt_i in mts]) - - self.summed_reactions[mt].add_temperature(strT, 0, - xsdata) - else: - raise ValueError('Data provided for an incorrect nuclide') - - else: - raise Warning('Temperature data set already within ' - 'IncidentNeutron object') + # Add normal and summed reactions + for mt in chain(data.reactions, data.summed_reactions): + if mt not in self: + raise ValueError("Tried to add cross sections for MT={} at T={}" + " but this reaction doesn't exist.".format( + mt, strT)) + self[mt].xs[strT] = data[mt].xs[strT] def get_reaction_components(self, mt): """Determine what reactions make up summed reaction. @@ -474,7 +442,6 @@ class IncidentNeutron(EqualityMixin): kTs = [] for temp in kTg: kTs.append(kTg[temp].value) - temperatures = [str(int(round(kT / K_BOLTZMANN))) + "K" for kT in kTs] data = cls(name, atomic_number, mass_number, metastable, atomic_weight_ratio, kTs) @@ -500,15 +467,12 @@ class IncidentNeutron(EqualityMixin): # MTs never depend on lower MTs. for mt_sum in sorted(SUM_RULES, reverse=True): if mt_sum not in data: - for it, T in enumerate(data.temperatures): - xs_components = \ - [data[mt].T_data[T].xs for mt in SUM_RULES[mt_sum] - if mt in data] + data.summed_reactions[mt_sum] = rx = Reaction(mt_sum) + for T in data.temperatures: + xs_components = [data[mt].xs[T] for mt in SUM_RULES[mt_sum] + if mt in data] if len(xs_components) > 0: - if it == 0: - data.summed_reactions[mt_sum] = Reaction(mt_sum) - data.summed_reactions[mt_sum].add_temperature( - T, 0, Sum(xs_components)) + rx.xs[T] = Sum(xs_components) # Read unresolved resonance probability tables if 'urr' in group: @@ -561,8 +525,6 @@ class IncidentNeutron(EqualityMixin): # Assign temperature to the running list kTs = [ace.temperature] - temperatures = [str(int(round(ace.temperature / K_BOLTZMANN))) + "K"] - # If mass number hasn't been specified, make an educated guess zaid, xs = ace.name.split('.') name, element, Z, mass_number, metastable = \ @@ -571,22 +533,24 @@ class IncidentNeutron(EqualityMixin): data = cls(name, Z, mass_number, metastable, ace.atomic_weight_ratio, kTs) + # Get string of temperature to use as a dictionary key + strT = data.temperatures[0] + # Read energy grid n_energy = ace.nxs[3] energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy] - data.energy[temperatures[0]] = energy + data.energy[strT] = energy total_xs = ace.xss[ace.jxs[1] + n_energy:ace.jxs[1] + 2 * n_energy] absorption_xs = ace.xss[ace.jxs[1] + 2*n_energy:ace.jxs[1] + 3 * n_energy] # Create summed reactions (total and absorption) total = Reaction(1) - total.add_temperature(temperatures[-1], 0, - Tabulated1D(energy, total_xs)) + total.xs[strT] = Tabulated1D(energy, total_xs) data.summed_reactions[1] = total + absorption = Reaction(27) - absorption.add_temperature(temperatures[-1], 0, - Tabulated1D(energy, absorption_xs)) + absorption.xs[strT] = Tabulated1D(energy, absorption_xs) data.summed_reactions[27] = absorption # Read each reaction @@ -616,12 +580,8 @@ class IncidentNeutron(EqualityMixin): warn('Photon production is present for MT={} but no ' 'reaction components exist.'.format(mt)) continue - threshold_idx = \ - np.amin([data.reactions[mt_i].T_data[temperatures[-1]].xs.x[0] - for mt_i in mts]) - xsvals = Sum([data.reactions[mt_i].T_data[temperatures[-1]].xs - for mt_i in mts]) - rx.add_temperature(temperatures[-1], threshold_idx, xsvals) + rx.xs[strT] = Sum([data.reactions[mt_i].xs[strT] + for mt_i in mts]) # Determine summed cross section rx.products += _get_photon_products(ace, rx) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index 9dc0d26a4..45d668b92 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1,5 +1,5 @@ from __future__ import division, unicode_literals -from collections import Iterable, Callable +from collections import Iterable, Callable, MutableMapping from copy import deepcopy from numbers import Real, Integral from warnings import warn @@ -11,7 +11,7 @@ from openmc.mixin import EqualityMixin from openmc.stats import Uniform from .angle_distribution import AngleDistribution from .angle_energy import AngleEnergy -from .function import Tabulated1D, Polynomial +from .function import Tabulated1D, Polynomial, Function1D from .data import REACTION_NAME, K_BOLTZMANN from .product import Product from .uncorrelated import UncorrelatedAngleEnergy @@ -211,7 +211,7 @@ def _get_photon_products(ace, rx): # Get photon production cross section photon_prod_xs = ace.xss[idx + 2:idx + 2 + n_energy] - neutron_xs = list(rx.T_data.values())[0].xs(energy) + neutron_xs = list(rx.xs.values())[0](energy) idx = np.where(neutron_xs > 0.) # Calculate photon yield @@ -251,78 +251,6 @@ def _get_photon_products(ace, rx): return photons -class XS(EqualityMixin): - """An energy-dependent cross-section for a reaction channel - - Parameters - ---------- - threshold_idx : int - The index on the energy grid corresponding to the threshold of this - reaction. - xs : dict of callable - Microscopic cross section for this reaction as a function of incident - energy; these cross sections are provided in a dictionary where the key - is the temperature of the cross section set. - - Attributes - ---------- - threshold : float - Threshold of the reaction in MeV - threshold_idx : int - The index on the energy grid corresponding to the threshold of this - reaction. - xs : dict of callable - Microscopic cross section for this reaction as a function of incident - energy; these cross sections are provided in a dictionary where the key - is the temperature of the cross section set. - - """ - - def __init__(self, threshold_idx, xs): - self._threshold_idx = threshold_idx - self._xs = xs - - @property - def threshold_idx(self): - return self._threshold_idx - - @property - def threshold(self): - return self.xs.x[0] - - @property - def xs(self): - return self._xs - - @threshold_idx.setter - def threshold_idx(self, threshold_idx): - cv.check_type('threshold_idx', threshold_idx, Integral) - cv.check_greater_than('threshold_idx', threshold_idx, 0, equality=True) - self._threshold_idx = threshold_idx - - @xs.setter - def xs(self, xs): - cv.check_type('reaction cross section', xs, Callable) - if isinstance(xs, Tabulated1D): - for y in xs.y: - cv.check_greater_than('reaction cross section', y, 0.0, True) - self._xs = xs - - def to_hdf5(self, group): - """Write XS to an HDF5 group - - Parameters - ---------- - group : h5py.Group - HDF5 group to write to - - """ - - group.attrs['threshold_idx'] = self.threshold_idx + 1 - if self.xs is not None: - group.create_dataset('xs', data=self.xs.y) - - class Reaction(EqualityMixin): """A nuclear reaction @@ -333,8 +261,7 @@ class Reaction(EqualityMixin): Parameters ---------- mt : int - The ENDF MT number for this reaction. On occasion, MCNP uses MT numbers - that don't correspond exactly to the ENDF specification. + The ENDF MT number for this reaction. Attributes ---------- @@ -346,14 +273,9 @@ class Reaction(EqualityMixin): The ENDF MT number for this reaction. q_value : float The Q-value of this reaction in MeV. - table : openmc.data.ace.Table - The ACE table which contains this reaction. threshold : float Threshold of the reaction in MeV - threshold_idx : int - The index on the energy grid corresponding to the threshold of this - reaction. - T_data : dict of openmc.data.XS + xs : dict of str to openmc.data.Function1D Microscopic cross section for this reaction as a function of incident energy; these cross sections are provided in a dictionary where the key is the temperature of the cross section set. @@ -366,12 +288,13 @@ class Reaction(EqualityMixin): """ def __init__(self, mt): - self.center_of_mass = True + self._center_of_mass = True + self._q_value = 0. + self._xs = {} + self._products = [] + self._derived_products = [] + self.mt = mt - self.q_value = 0. - self.T_data = {} - self.products = [] - self.derived_products = [] def __repr__(self): if self.mt in REACTION_NAME: @@ -391,6 +314,14 @@ class Reaction(EqualityMixin): def products(self): return self._products + @property + def derived_products(self): + return self._derived_products + + @property + def xs(self): + return self._xs + @center_of_mass.setter def center_of_mass(self, center_of_mass): cv.check_type('center of mass', center_of_mass, (bool, np.bool_)) @@ -406,9 +337,19 @@ class Reaction(EqualityMixin): cv.check_type('reaction products', products, Iterable, Product) self._products = products - def add_temperature(self, temperature, threshold_idx, xs): - cv.check_type('temperature', temperature, str) - self.T_data[temperature] = XS(threshold_idx, xs) + @derived_products.setter + def derived_products(self, derived_products): + cv.check_type('reaction derived products', derived_products, + Iterable, Product) + self._derived_products = derived_products + + @xs.setter + def xs(self, xs): + cv.check_type('reaction cross section dictionary', xs, MutableMapping) + for key, value in xs.items(): + cv.check_type('reaction cross section temperature', key, basestring) + cv.check_type('reaction cross section', value, Function1D) + self._xs = xs def to_hdf5(self, group): """Write reaction to an HDF5 group @@ -427,27 +368,30 @@ class Reaction(EqualityMixin): group.attrs['label'] = np.string_(self.mt) group.attrs['Q_value'] = self.q_value group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 - for T in self.T_data: + for T in self.xs: Tgroup = group.create_group(T) - self.T_data[T].to_hdf5(Tgroup) + if self.xs[T] is not None: + dset = Tgroup.create_dataset('xs', data=self.xs[T].y) + if hasattr(self.xs[T], '_threshold_idx'): + threshold_idx = self.xs[T]._threshold_idx + 1 + else: + threshold_idx = 1 + dset.attrs['threshold_idx'] = threshold_idx for i, p in enumerate(self.products): pgroup = group.create_group('product_{}'.format(i)) p.to_hdf5(pgroup) @classmethod - def from_hdf5(cls, group, energy, temperatures): + def from_hdf5(cls, group, energy): """Generate reaction from an HDF5 group Parameters ---------- group : h5py.Group HDF5 group to write to - energy : Iterable of float - Array of energies at which cross sections are tabulated at - temperatures : Iterable of str - List of string representations the temperatures of the target - nuclide in the data set. The temperatures are strings of the - temperature, rounded to the nearest integer; e.g., '294K' + energy : dict + Dictionary whose keys are temperatures (e.g., '300K') and values are + arrays of energies at which cross sections are tabulated at. Returns ------- @@ -461,13 +405,21 @@ class Reaction(EqualityMixin): rx.q_value = group.attrs['Q_value'] rx.center_of_mass = bool(group.attrs['center_of_mass']) - # Read cross section data - for T in temperatures: - Tgroup = group[T] - if 'xs' in Tgroup: - threshold_idx = Tgroup.attrs['threshold_idx'] - 1 - xs = Tabulated1D(energy[T][threshold_idx:], Tgroup['xs'].value) - rx.T_data[T] = XS(threshold_idx, xs) + # Read cross section at each temperature + for T, Tgroup in group.items(): + if T.endswith('K'): + if 'xs' in Tgroup: + # Make sure temperature has associated energy grid + if T not in energy: + raise ValueError( + 'Could not create reaction cross section for MT={} ' + 'at T={} because no corresponding energy grid ' + 'exists.'.format(mt, T)) + xs = Tgroup['xs'].value + threshold_idx = Tgroup['xs'].attrs['threshold_idx'] - 1 + tabulated_xs = Tabulated1D(energy[T][threshold_idx:], xs) + tabulated_xs._threshold_idx = threshold_idx + rx.xs[T] = tabulated_xs # Determine number of products n_product = 0 @@ -522,7 +474,8 @@ class Reaction(EqualityMixin): xs[xs < 0.0] = 0.0 tabulated_xs = Tabulated1D(energy, xs) - rx.T_data[strT] = XS(threshold_idx, tabulated_xs) + tabulated_xs._threshold_idx = threshold_idx + rx.xs[strT] = tabulated_xs # ================================================================== # YIELD AND ANGLE-ENERGY DISTRIBUTION @@ -582,7 +535,8 @@ class Reaction(EqualityMixin): elastic_xs[elastic_xs < 0.0] = 0.0 tabulated_xs = Tabulated1D(grid, elastic_xs) - rx.T_data[strT] = XS(0, tabulated_xs) + tabulated_xs._threshold_idx = 0 + rx.xs[strT] = tabulated_xs # No energy distribution for elastic scattering neutron = Product('neutron')