From 3b65dde0fd3dace45bfa422d505d1503694d4e77 Mon Sep 17 00:00:00 2001 From: Adam Nelson Date: Sun, 21 Aug 2016 06:05:19 -0400 Subject: [PATCH] Folded multiple temperatures in to each h5 data file; now im getting 50% less storage required of endf70. yay! Next I think I can get rid of some s(a,b) information, will investigate that and then on to revising openMC so it can read and use this data --- openmc/data/data.py | 4 + openmc/data/neutron.py | 288 +++++++++++++++++++++--------- openmc/data/reaction.py | 141 +++++++++++---- openmc/data/thermal.py | 354 +++++++++++++++++++++++++++++-------- scripts/openmc-ace-to-hdf5 | 92 ++++++---- 5 files changed, 657 insertions(+), 222 deletions(-) diff --git a/openmc/data/data.py b/openmc/data/data.py index ebc2cba43..4eb7860c2 100644 --- a/openmc/data/data.py +++ b/openmc/data/data.py @@ -218,3 +218,7 @@ def atomic_mass(isotope): isotope = isotope[:isotope.find('_')] return _ATOMIC_MASS.get(isotope.lower()) + +def kT_to_K(kT): + K = kT / 8.6173324e-11 + return K \ No newline at end of file diff --git a/openmc/data/neutron.py b/openmc/data/neutron.py index c6621e1b8..8ea8e38e3 100644 --- a/openmc/data/neutron.py +++ b/openmc/data/neutron.py @@ -7,7 +7,7 @@ from warnings import warn import numpy as np import h5py -from .data import ATOMIC_SYMBOL, SUM_RULES +from .data import ATOMIC_SYMBOL, SUM_RULES, kT_to_K from .ace import Table, get_table from .fission_energy import FissionEnergyRelease from .function import Tabulated1D, Sum @@ -21,6 +21,64 @@ if sys.version_info[0] >= 3: basestring = str +def _get_metadata(zaid, metastable_scheme='nndc'): + """Method to obtain the complete element name, element, Z, mass_number, + and metastable state + + Parameters + ---------- + zaid : int + ZAID (1000*Z + A) obtained from a library + metastable_scheme : {'nndc', 'mcnp'} + Determine how ZAID identifiers are to be interpreted in the case of + a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not + encode metastable information, different conventions are used among + different libraries. In MCNP libraries, the convention is to add 400 + for a metastable nuclide except for Am242m, for which 95242 is + metastable and 95642 (or 1095242 in newer libraries) is the ground + state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m. + + Returns + ------- + + """ + + cv.check_type('zaid', zaid, int) + cv.check_value('metastable_scheme', metastable_scheme, ['nndc', 'mcnp']) + + Z = zaid // 1000 + mass_number = zaid % 1000 + + if metastable_scheme == 'mcnp': + if zaid > 1000000: + # New SZA format + Z = Z % 1000 + if zaid == 1095242: + metastable = 0 + else: + metastable = zaid // 1000000 + else: + if zaid == 95242: + metastable = 1 + elif zaid == 95642: + metastable = 0 + else: + metastable = 1 if mass_number > 300 else 0 + elif metastable_scheme == 'nndc': + metastable = 1 if mass_number > 300 else 0 + + while mass_number > 3 * Z: + mass_number -= 100 + + # Determine name + element = ATOMIC_SYMBOL[Z] + name = '{}{}'.format(element, mass_number) + if metastable > 0: + name += '_m{}'.format(metastable) + + return (name, element, Z, mass_number, metastable) + + class IncidentNeutron(EqualityMixin): """Continuous-energy neutron interaction data. @@ -40,8 +98,9 @@ class IncidentNeutron(EqualityMixin): Metastable state of the nucleus. A value of zero indicates ground state. atomic_weight_ratio : float Atomic mass ratio of the target nuclide. - temperature : float - Temperature of the target nuclide in MeV. + kTs : Iterable float + List of temperatures of the target nuclide in the data set. + The temperatures have units of MeV. Attributes ---------- @@ -51,8 +110,10 @@ class IncidentNeutron(EqualityMixin): Atomic symbol of the nuclide, e.g., 'Zr' atomic_weight_ratio : float Atomic weight ratio of the target nuclide. - energy : numpy.ndarray + energy : dict of numpy.ndarray The energy values (MeV) at which reaction cross-sections are tabulated. + They keys of the dict are the temperature string ('296.3K') for each + set of energies fission_energy : None or openmc.data.FissionEnergyRelease The energy released by fission, tabulated by component (e.g. prompt neutrons or beta particles) and dependent on incident neutron energy @@ -69,23 +130,28 @@ class IncidentNeutron(EqualityMixin): summed_reactions : collections.OrderedDict Contains summed cross sections, e.g., the total cross section. The keys are the MT values and the values are Reaction objects. - temperature : float - Temperature of the target nuclide in MeV. + temperatures : Iterable of str + List of string representations the temperatures of the target nuclide + in the data set. The temperatures are strings with 1 decimal place, + i.e., '293.6K' + kTs : Iterable of float + List of temperatures of the target nuclide in the data set. + The temperatures have units of MeV. urr : None or openmc.data.ProbabilityTables Unresolved resonance region probability tables """ def __init__(self, name, atomic_number, mass_number, metastable, - atomic_weight_ratio, temperature): + atomic_weight_ratio, kTs): self.name = name self.atomic_number = atomic_number self.mass_number = mass_number self.metastable = metastable self.atomic_weight_ratio = atomic_weight_ratio - self.temperature = temperature - - self._energy = None + self.kTs = kTs + self.temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs] + self.energy = {} self._fission_energy = None self.reactions = OrderedDict() self.summed_reactions = OrderedDict() @@ -128,18 +194,10 @@ class IncidentNeutron(EqualityMixin): def atomic_weight_ratio(self): return self._atomic_weight_ratio - @property - def energy(self): - return self._energy - @property def fission_energy(self): return self._fission_energy - @property - def temperature(self): - return self._temperature - @property def reactions(self): return self._reactions @@ -159,7 +217,7 @@ class IncidentNeutron(EqualityMixin): @property def atomic_symbol(self): - return atomic_symbol[self.atomic_number] + return ATOMIC_SYMBOL[self.atomic_number] @atomic_number.setter def atomic_number(self, atomic_number): @@ -185,17 +243,6 @@ class IncidentNeutron(EqualityMixin): cv.check_greater_than('atomic weight ratio', atomic_weight_ratio, 0.0) self._atomic_weight_ratio = atomic_weight_ratio - @temperature.setter - def temperature(self, temperature): - cv.check_type('temperature', temperature, Real) - cv.check_greater_than('temperature', temperature, 0.0, True) - self._temperature = temperature - - @energy.setter - def energy(self, energy): - cv.check_type('energy grid', energy, Iterable, Real) - self._energy = energy - @fission_energy.setter def fission_energy(self, fission_energy): cv.check_type('fission energy release', fission_energy, @@ -218,6 +265,87 @@ class IncidentNeutron(EqualityMixin): (ProbabilityTables, type(None))) self._urr = urr + def add_temperature_from_ace(self, ace_or_filename, metastable_scheme='nndc'): + """Add data to the IncidentNeutron object from an ACE file at a + different temperature. + + Parameters + ---------- + ace_or_filename : openmc.data.ace.Table or str + ACE table to read from. If given as a string, it is assumed to be + the filename for the ACE file. + metastable_scheme : {'nndc', 'mcnp'} + Determine how ZAID identifiers are to be interpreted in the case of + a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not + encode metastable information, different conventions are used among + different libraries. In MCNP libraries, the convention is to add 400 + for a metastable nuclide except for Am242m, for which 95242 is + metastable and 95642 (or 1095242 in newer libraries) is the ground + state. For NNDC libraries, ZAID is given as 1000*Z + A + 100*m. + + """ + + if isinstance(ace_or_filename, Table): + ace = ace_or_filename + else: + ace = get_table(ace_or_filename) + + # 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) + + # 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 = "{0:.1f}K".format(kT_to_K(ace.temperature)) + 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]) + + 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) + + # 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) + + 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') + def get_reaction_components(self, mt): """Determine what reactions make up summed reaction. @@ -271,10 +399,12 @@ class IncidentNeutron(EqualityMixin): g.attrs['A'] = self.mass_number g.attrs['metastable'] = self.metastable g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio - g.attrs['temperature'] = self.temperature + g.attrs['kTs'] = self.kTs # Write energy grid - g.create_dataset('energy', data=self.energy) + eg = g.create_group('energy') + for temperature in self.temperatures: + eg.create_dataset(temperature, data=self.energy[temperature]) # Write reaction data rxs_group = g.create_group('reactions') @@ -327,19 +457,22 @@ class IncidentNeutron(EqualityMixin): mass_number = group.attrs['A'] metastable = group.attrs['metastable'] atomic_weight_ratio = group.attrs['atomic_weight_ratio'] - temperature = group.attrs['temperature'] + kTs = group.attrs['kTs'].tolist() + temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs] data = cls(name, atomic_number, mass_number, metastable, - atomic_weight_ratio, temperature) + atomic_weight_ratio, kTs) # Read energy grid - data.energy = group['energy'].value + e_group = group['energy'] + for temperature in temperatures: + data.energy[temperature] = e_group[temperature].value # Read reaction data rxs_group = group['reactions'] for name, obj in sorted(rxs_group.items()): if name.startswith('reaction_'): - rx = Reaction.from_hdf5(obj, data.energy) + rx = Reaction.from_hdf5(obj, data.energy, temperatures) data.reactions[rx.mt] = rx # Read total nu data if available @@ -351,12 +484,15 @@ class IncidentNeutron(EqualityMixin): # MTs never depend on lower MTs. for mt_sum in sorted(SUM_RULES, reverse=True): if mt_sum not in data: - xs_components = [data[mt].xs for mt in SUM_RULES[mt_sum] - if mt in data] - if len(xs_components) > 0: - rxn = Reaction(mt_sum) - rxn.xs = Sum(xs_components) - data.summed_reactions[mt_sum] = rxn + 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] + 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)) # Read unresolved resonance probability tables if 'urr' in group: @@ -376,9 +512,9 @@ class IncidentNeutron(EqualityMixin): Parameters ---------- - ace : openmc.data.ace.Table or str - ACE table to read from. If given as a string, it is assumed to be - the filename for the ACE file. + ace_or_filename : openmc.data.ace.Table or str + ACE table to read from. If the value is a string, it is assumed to + be the filename for the ACE file. metastable_scheme : {'nndc', 'mcnp'} Determine how ZAID identifiers are to be interpreted in the case of a metastable nuclide. Because the normal ZAID (=1000*Z + A) does not @@ -394,6 +530,8 @@ class IncidentNeutron(EqualityMixin): Incident neutron continuous-energy data """ + + # First obtain the data for the first provided ACE table/file if isinstance(ace_or_filename, Table): ace = ace_or_filename else: @@ -401,54 +539,37 @@ class IncidentNeutron(EqualityMixin): # If mass number hasn't been specified, make an educated guess zaid, xs = ace.name.split('.') - zaid = int(zaid) - Z = zaid // 1000 - mass_number = zaid % 1000 + name, element, Z, mass_number, metastable = \ + _get_metadata(int(zaid), metastable_scheme) - if metastable_scheme == 'mcnp': - if zaid > 1000000: - # New SZA format - Z = Z % 1000 - if zaid == 1095242: - metastable = 0 - else: - metastable = zaid // 1000000 - else: - if zaid == 95242: - metastable = 1 - elif zaid == 95642: - metastable = 0 - else: - metastable = 1 if mass_number > 300 else 0 - elif metastable_scheme == 'nndc': - metastable = 1 if mass_number > 300 else 0 + # Assign temperature to the running list + kTs = [ace.temperature] + temperatures = ["{0:.1f}K".format(kT_to_K(ace.temperature))] - while mass_number > 3*Z: - mass_number -= 100 - - # Determine name for group - element = ATOMIC_SYMBOL[Z] - if metastable > 0: - name = '{}{}_m{}.{}'.format(element, mass_number, metastable, xs) - else: - name = '{}{}.{}'.format(element, mass_number, xs) + # If mass number hasn't been specified, make an educated guess + zaid, xs = ace.name.split('.') + name, element, Z, mass_number, metastable = \ + _get_metadata(int(zaid), metastable_scheme) data = cls(name, Z, mass_number, metastable, - ace.atomic_weight_ratio, ace.temperature) + ace.atomic_weight_ratio, kTs) # Read energy grid n_energy = ace.nxs[3] energy = ace.xss[ace.jxs[1]:ace.jxs[1] + n_energy] - data.energy = 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] + data.energy[temperatures[0]] = 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.xs = Tabulated1D(energy, total_xs) + total.add_temperature(temperatures[-1], 0, + Tabulated1D(energy, total_xs)) data.summed_reactions[1] = total absorption = Reaction(27) - absorption.xs = Tabulated1D(energy, absorption_xs) + absorption.add_temperature(temperatures[-1], 0, + Tabulated1D(energy, absorption_xs)) data.summed_reactions[27] = absorption # Read each reaction @@ -478,7 +599,12 @@ class IncidentNeutron(EqualityMixin): warn('Photon production is present for MT={} but no ' 'reaction components exist.'.format(mt)) continue - rx.xs = Sum([data.reactions[mt_i].xs for mt_i in mts]) + 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) # Determine summed cross section rx.products += _get_photon_products(ace, rx) diff --git a/openmc/data/reaction.py b/openmc/data/reaction.py index e7580ec47..741822bf0 100644 --- a/openmc/data/reaction.py +++ b/openmc/data/reaction.py @@ -1,7 +1,7 @@ from __future__ import division, unicode_literals from collections import Iterable, Callable from copy import deepcopy -from numbers import Real +from numbers import Real, Integral from warnings import warn import numpy as np @@ -12,7 +12,7 @@ from openmc.stats import Uniform from .angle_distribution import AngleDistribution from .angle_energy import AngleEnergy from .function import Tabulated1D, Polynomial -from .data import REACTION_NAME +from .data import REACTION_NAME, kT_to_K 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 = rx.xs(energy) + neutron_xs = list(rx.T_data.values())[0].xs(energy) idx = np.where(neutron_xs > 0.) # Calculate photon yield @@ -251,6 +251,78 @@ 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 @@ -281,9 +353,10 @@ class Reaction(EqualityMixin): threshold_idx : int The index on the energy grid corresponding to the threshold of this reaction. - xs : callable + T_data : dict of openmc.data.XS Microscopic cross section for this reaction as a function of incident - energy + energy; these cross sections are provided in a dictionary where the key + is the temperature of the cross section set. products : Iterable of openmc.data.Product Reaction products derived_products : Iterable of openmc.data.Product @@ -296,8 +369,7 @@ class Reaction(EqualityMixin): self.center_of_mass = True self.mt = mt self.q_value = 0. - self.threshold_idx = 0 - self._xs = None + self.T_data = {} self.products = [] self.derived_products = [] @@ -319,14 +391,6 @@ class Reaction(EqualityMixin): def products(self): return self._products - @property - def threshold(self): - return self.xs.x[0] - - @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_)) @@ -342,13 +406,9 @@ class Reaction(EqualityMixin): cv.check_type('reaction products', products, Iterable, Product) self._products = products - @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 add_temperature(self, temperature, threshold_idx, xs): + cv.check_type('temperature', temperature, str) + self.T_data[temperature] = XS(threshold_idx, xs) def to_hdf5(self, group): """Write reaction to an HDF5 group @@ -366,16 +426,16 @@ class Reaction(EqualityMixin): else: group.attrs['label'] = np.string_(self.mt) group.attrs['Q_value'] = self.q_value - group.attrs['threshold_idx'] = self.threshold_idx + 1 group.attrs['center_of_mass'] = 1 if self.center_of_mass else 0 - if self.xs is not None: - group.create_dataset('xs', data=self.xs.y) + for T in self.T_data: + Tgroup = group.create_group(T) + self.T_data[T].to_hdf5(Tgroup) 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): + def from_hdf5(cls, group, energy, temperatures): """Generate reaction from an HDF5 group Parameters @@ -384,6 +444,8 @@ class Reaction(EqualityMixin): HDF5 group to write to energy : Iterable of float Array of energies at which cross sections are tabulated at + temperatures : Iterable of float + Array of temperatures at which to obtain the cross sections Returns ------- @@ -391,16 +453,19 @@ class Reaction(EqualityMixin): Reaction data """ + mt = group.attrs['mt'] rx = cls(mt) rx.q_value = group.attrs['Q_value'] - rx.threshold_idx = group.attrs['threshold_idx'] - 1 rx.center_of_mass = bool(group.attrs['center_of_mass']) - # Read cross section - if 'xs' in group: - xs = group['xs'].value - rx.xs = Tabulated1D(energy[rx.threshold_idx:], xs) + # 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) # Determine number of products n_product = 0 @@ -421,6 +486,10 @@ class Reaction(EqualityMixin): n_grid = ace.nxs[3] grid = ace.xss[ace.jxs[1]:ace.jxs[1] + n_grid] + # Convert data temperature to a "300.0K" number for indexing + # temperature data + strT = "{0:.1f}K".format(kT_to_K(ace.temperature)) + if i_reaction > 0: mt = int(ace.xss[ace.jxs[3] + i_reaction - 1]) rx = cls(mt) @@ -435,11 +504,11 @@ class Reaction(EqualityMixin): loc = int(ace.xss[ace.jxs[6] + i_reaction - 1]) # Determine starting index on energy grid - rx.threshold_idx = int(ace.xss[ace.jxs[7] + loc - 1]) - 1 + threshold_idx = int(ace.xss[ace.jxs[7] + loc - 1]) - 1 # Determine number of energies in reaction n_energy = int(ace.xss[ace.jxs[7] + loc]) - energy = grid[rx.threshold_idx:rx.threshold_idx + n_energy] + energy = grid[threshold_idx:threshold_idx + n_energy] # Read reaction cross section xs = ace.xss[ace.jxs[7] + loc + 1:ace.jxs[7] + loc + 1 + n_energy] @@ -450,7 +519,8 @@ class Reaction(EqualityMixin): "to zero.".format(rx.mt, ace.name)) xs[xs < 0.0] = 0.0 - rx.xs = Tabulated1D(energy, xs) + tabulated_xs = Tabulated1D(energy, xs) + rx.T_data[strT] = XS(threshold_idx, tabulated_xs) # ================================================================== # YIELD AND ANGLE-ENERGY DISTRIBUTION @@ -509,7 +579,8 @@ class Reaction(EqualityMixin): "Setting to zero.".format(ace.name)) elastic_xs[elastic_xs < 0.0] = 0.0 - rx.xs = Tabulated1D(grid, elastic_xs) + tabulated_xs = Tabulated1D(grid, elastic_xs) + rx.T_data[strT] = XS(0, tabulated_xs) # No energy distribution for elastic scattering neutron = Product('neutron') diff --git a/openmc/data/thermal.py b/openmc/data/thermal.py index 1551e8ac5..27c61e924 100644 --- a/openmc/data/thermal.py +++ b/openmc/data/thermal.py @@ -8,6 +8,7 @@ import h5py import openmc.checkvalue as cv from openmc.mixin import EqualityMixin +from .data import kT_to_K from .ace import Table, get_table from .angle_energy import AngleEnergy from .function import Tabulated1D @@ -89,7 +90,7 @@ class CoherentElastic(EqualityMixin): if isinstance(E, Iterable): E = np.asarray(E) idx = np.searchsorted(self.bragg_edges, E) - return self.factors[idx]/E + return self.factors[idx] / E def __len__(self): @@ -159,8 +160,9 @@ class ThermalScattering(EqualityMixin): ZAID identifier of the table, e.g. lwtr.10t. atomic_weight_ratio : float Atomic mass ratio of the target nuclide. - temperature : float - Temperature of the target nuclide in eV. + kTs : Iterable of float + List of temperatures of the target nuclide in the data set. + The temperatures have units of MeV. Attributes ---------- @@ -174,22 +176,29 @@ class ThermalScattering(EqualityMixin): approximation name : str Name of the table, e.g. lwtr.20t. - temperature : float - Temperature of the target nuclide in eV. + temperatures : Iterable of str + List of string representations the temperatures of the target nuclide + in the data set. The temperatures are strings with 1 decimal place, + i.e., '293.6K' + kTs : Iterable of float + List of temperatures of the target nuclide in the data set. + The temperatures have units of MeV. zaids : Iterable of int ZAID identifiers that the thermal scattering data applies to """ - def __init__(self, name, atomic_weight_ratio, temperature): + def __init__(self, name, atomic_weight_ratio, kTs): self.name = name self.atomic_weight_ratio = atomic_weight_ratio - self.temperature = temperature - self.elastic_xs = None - self.elastic_mu_out = None - self.inelastic_xs = None - self.inelastic_e_out = None - self.inelastic_mu_out = None + self.kTs = kTs + self.temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs] + self.elastic_xs = {} + self.elastic_mu_out = {} + self.inelastic_xs = {} + self.inelastic_e_out = {} + self.inelastic_mu_out = {} + self.inelastic_dist = {} self.secondary_mode = None self.zaids = [] @@ -217,26 +226,191 @@ class ThermalScattering(EqualityMixin): # Write basic data g = f.create_group(self.name) g.attrs['atomic_weight_ratio'] = self.atomic_weight_ratio - g.attrs['temperature'] = self.temperature + g.attrs['kTs'] = self.kTs g.attrs['zaids'] = self.zaids + g.attrs['secondary_mode'] = np.string_(self.secondary_mode) - # Write thermal elastic scattering - if self.elastic_xs is not None: - elastic_group = g.create_group('elastic') - self.elastic_xs.to_hdf5(elastic_group, 'xs') - if self.elastic_mu_out is not None: - elastic_group.create_dataset('mu_out', data=self.elastic_mu_out) + for T in self.temperatures: + Tg = g.create_group(T) + # Write thermal elastic scattering + if self.elastic_xs: + elastic_group = Tg.create_group('elastic') - # Write thermal inelastic scattering - if self.inelastic_xs is not None: - inelastic_group = g.create_group('inelastic') - self.inelastic_xs.to_hdf5(inelastic_group, 'xs') - inelastic_group.attrs['secondary_mode'] = np.string_(self.secondary_mode) - if self.secondary_mode in ('equal', 'skewed'): - inelastic_group.create_dataset('energy_out', data=self.inelastic_e_out) - inelastic_group.create_dataset('mu_out', data=self.inelastic_mu_out) - elif self.secondary_mode == 'continuous': - self.inelastic_dist.to_hdf5(inelastic_group) + self.elastic_xs[T].to_hdf5(elastic_group, 'xs') + if self.elastic_mu_out: + elastic_group.create_dataset('mu_out', + data=self.elastic_mu_out[T]) + + # Write thermal inelastic scattering + if self.inelastic_xs: + inelastic_group = Tg.create_group('inelastic') + self.inelastic_xs[T].to_hdf5(inelastic_group, 'xs') + if self.secondary_mode in ('equal', 'skewed'): + inelastic_group.create_dataset('energy_out', + data=self.inelastic_e_out[T]) + inelastic_group.create_dataset('mu_out', + data=self.inelastic_mu_out[T]) + elif self.secondary_mode == 'continuous': + self.inelastic_dist[T].to_hdf5(inelastic_group) + + f.close() + + def add_temperature_from_ace(self, ace_or_filename, name=None): + """Add data to the ThermalScattering object from an ACE file at a + different temperature. + + Parameters + ---------- + ace_or_filename : openmc.data.ace.Table or str + ACE table to read from. If given as a string, it is assumed to be + the filename for the ACE file. + name : str + GND-conforming name of the material, e.g. c_H_in_H2O. If none is + passed, the appropriate name is guessed based on the name of the ACE + table. + + Returns + ------- + openmc.data.ThermalScattering + Thermal scattering data + + """ + if isinstance(ace_or_filename, Table): + ace = ace_or_filename + else: + ace = get_table(ace_or_filename) + + # Get new name that is GND-consistent + ace_name, xs = ace.name.split('.') + if name is None: + if ace_name.lower() in _THERMAL_NAMES: + name = _THERMAL_NAMES[ace_name.lower()] + else: + # Make an educated guess? This actually works well for JEFF-3.2 + # which stupidly uses names like lw00.32t, lw01.32t, etc. for + # different temperatures + matches = get_close_matches( + ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) + if len(matches) > 0: + name = _THERMAL_NAMES[matches[0]] + else: + # OK, we give up. Just use the ACE name. + name = 'c_' + ace.name + warn('Thermal scattering material "{}" is not recognized. ' + 'Assigning a name of {}.'.format(ace.name, name)) + + # 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 = "{0:.1f}K".format(kT_to_K(ace.temperature)) + self.temperatures.append(strT) + self.kTs.append(ace.temperature) + + # Incoherent inelastic scattering cross section + idx = ace.jxs[1] + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy] + xs = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + self.inelastic_xs[strT] = Tabulated1D(energy, xs) + + # Make sure secondary_mode is always equal. This should always + # be the case, but to reduce future debugging should something + # change, this will alert the developers to the issue. + if ace.nxs[7] == 0: + secondary_mode = 'equal' + elif ace.nxs[7] == 1: + secondary_mode = 'skewed' + elif ace.nxs[7] == 2: + secondary_mode = 'continuous' + + if secondary_mode != self.secondary_mode: + raise ValueError('Secondary Modes are inconsistent.') + + n_energy_out = ace.nxs[4] + if self.secondary_mode in ('equal', 'skewed'): + n_mu = ace.nxs[3] + idx = ace.jxs[3] + self.inelastic_e_out[strT] = \ + ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2): + n_mu + 2] + self.inelastic_e_out[strT].shape = \ + (n_energy, n_energy_out) + + self.inelastic_mu_out[strT] = \ + ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)] + self.inelastic_mu_out[strT].shape = \ + (n_energy, n_energy_out, n_mu + 2) + self.inelastic_mu_out[strT] = \ + self.inelastic_mu_out[strT][:, :, 1:] + else: + n_mu = ace.nxs[3] - 1 + idx = ace.jxs[3] + locc = ace.xss[idx:idx + n_energy].astype(int) + n_energy_out = \ + ace.xss[idx + n_energy:idx + 2 * n_energy].astype(int) + energy_out = [] + mu_out = [] + for i in range(n_energy): + idx = locc[i] + + # Outgoing energy distribution for incoming energy i + e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3): + n_mu + 3] + p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3): + n_mu + 3] + c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3): + n_mu + 3] + eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True) + eout_i.c = c + + # Outgoing angle distribution for each + # (incoming, outgoing) energy pair + mu_i = [] + for j in range(n_energy_out[i]): + mu = ace.xss[idx + 4:idx + 4 + n_mu] + p_mu = 1. / n_mu * np.ones(n_mu) + mu_ij = Discrete(mu, p_mu) + mu_ij.c = np.cumsum(p_mu) + mu_i.append(mu_ij) + idx += 3 + n_mu + + energy_out.append(eout_i) + mu_out.append(mu_i) + + # Create correlated angle-energy distribution + breakpoints = [n_energy] + interpolation = [2] + energy = self.inelastic_xs[strT].x + self.inelastic_dist[strT] = CorrelatedAngleEnergy( + breakpoints, interpolation, energy, energy_out, mu_out) + + # Incoherent/coherent elastic scattering cross section + idx = ace.jxs[4] + if idx != 0: + n_energy = int(ace.xss[idx]) + energy = ace.xss[idx + 1: idx + 1 + n_energy] + P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] + + if ace.nxs[5] == 4: + self.elastic_xs[strT] = CoherentElastic(energy, P) + else: + self.elastic_xs[strT] = Tabulated1D(energy, P) + + # Angular distribution + n_mu = ace.nxs[6] + if n_mu != -1: + idx = ace.jxs[6] + self.elastic_mu_out[strT] = \ + ace.xss[idx:idx + n_energy * n_mu] + self.elastic_mu_out[strT].shape = \ + (n_energy, n_mu) + + else: + raise ValueError('Data provided for an incorrect library') + else: + raise Warning('Temperature data set already within ' + 'IncidentNeutron object') @classmethod def from_hdf5(cls, group_or_filename): @@ -263,35 +437,46 @@ class ThermalScattering(EqualityMixin): name = group.name[1:] atomic_weight_ratio = group.attrs['atomic_weight_ratio'] - temperature = group.attrs['temperature'] - table = cls(name, atomic_weight_ratio, temperature) + kTs = group.attrs['kTs'].tolist() + temperatures = ["{0:.1f}K".format(kT_to_K(kT)) for kT in kTs] + + table = cls(name, atomic_weight_ratio, kTs) table.zaids = group.attrs['zaids'] + table.secondary_mode = group.attrs['secondary_mode'].decode() # Read thermal elastic scattering - if 'elastic' in group: - elastic_group = group['elastic'] + for T in temperatures: + Tgroup = group[T] + if 'elastic' in Tgroup: + elastic_group = Tgroup['elastic'] - # Cross section - elastic_xs_type = elastic_group['xs'].attrs['type'].decode() - if elastic_xs_type == 'tab1': - table.elastic_xs = Tabulated1D.from_hdf5(elastic_group['xs']) - elif elastic_xs_type == 'bragg': - table.elastic_xs = CoherentElastic.from_hdf5(elastic_group['xs']) + # Cross section + elastic_xs_type = elastic_group['xs'].attrs['type'].decode() + if elastic_xs_type == 'Tabulated1D': + table.elastic_xs[T] = \ + Tabulated1D.from_hdf5(elastic_group['xs']) + elif elastic_xs_type == 'bragg': + table.elastic_xs[T] = \ + CoherentElastic.from_hdf5(elastic_group['xs']) - # Angular distribution - if 'mu_out' in elastic_group: - table.elastic_mu_out = elastic_group['mu_out'].value + # Angular distribution + if 'mu_out' in elastic_group: + table.elastic_mu_out[T] = \ + elastic_group['mu_out'].value - # Read thermal inelastic scattering - if 'inelastic' in group: - inelastic_group = group['inelastic'] - table.secondary_mode = inelastic_group.attrs['secondary_mode'].decode() - table.inelastic_xs = Tabulated1D.from_hdf5(inelastic_group['xs']) - if table.secondary_mode in ('equal', 'skewed'): - table.inelastic_e_out = inelastic_group['energy_out'] - table.inelastic_mu_out = inelastic_group['mu_out'] - elif table.secondary_mode == 'continuous': - table.inelastic_dist = AngleEnergy.from_hdf5(inelastic_group) + # Read thermal inelastic scattering + if 'inelastic' in Tgroup: + inelastic_group = Tgroup['inelastic'] + table.inelastic_xs[T] = \ + Tabulated1D.from_hdf5(inelastic_group['xs']) + if table.secondary_mode in ('equal', 'skewed'): + table.inelastic_e_out[T] = \ + inelastic_group['energy_out'] + table.inelastic_mu_out[T] = \ + inelastic_group['mu_out'] + elif table.secondary_mode == 'continuous': + table.inelastic_dist[T] = \ + AngleEnergy.from_hdf5(inelastic_group) return table @@ -301,7 +486,7 @@ class ThermalScattering(EqualityMixin): Parameters ---------- - ace : openmc.data.ace.Table or str + ace_or_filename : openmc.data.ace.Table or str ACE table to read from. If given as a string, it is assumed to be the filename for the ACE file. name : str @@ -324,7 +509,7 @@ class ThermalScattering(EqualityMixin): ace_name, xs = ace.name.split('.') if name is None: if ace_name.lower() in _THERMAL_NAMES: - name = _THERMAL_NAMES[ace_name.lower()] + '.' + xs + name = _THERMAL_NAMES[ace_name.lower()] else: # Make an educated guess?? This actually works well for JEFF-3.2 # which stupidly uses names like lw00.32t, lw01.32t, etc. for @@ -332,21 +517,25 @@ class ThermalScattering(EqualityMixin): matches = get_close_matches( ace_name.lower(), _THERMAL_NAMES.keys(), cutoff=0.5) if len(matches) > 0: - name = _THERMAL_NAMES[matches[0]] + '.' + xs + name = _THERMAL_NAMES[matches[0]] else: # OK, we give up. Just use the ACE name. name = 'c_' + ace.name warn('Thermal scattering material "{}" is not recognized. ' 'Assigning a name of {}.'.format(ace.name, name)) - table = cls(name, ace.atomic_weight_ratio, ace.temperature) + # Assign temperature to the running list + kTs = [ace.temperature] + temperatures = ["{0:.1f}K".format(kT_to_K(ace.temperature))] + + table = cls(name, ace.atomic_weight_ratio, kTs) # Incoherent inelastic scattering cross section idx = ace.jxs[1] n_energy = int(ace.xss[idx]) energy = ace.xss[idx+1 : idx+1+n_energy] xs = ace.xss[idx+1+n_energy : idx+1+2*n_energy] - table.inelastic_xs = Tabulated1D(energy, xs) + table.inelastic_xs[temperatures[0]] = Tabulated1D(energy, xs) if ace.nxs[7] == 0: table.secondary_mode = 'equal' @@ -359,34 +548,45 @@ class ThermalScattering(EqualityMixin): if table.secondary_mode in ('equal', 'skewed'): n_mu = ace.nxs[3] idx = ace.jxs[3] - table.inelastic_e_out = ace.xss[idx:idx+n_energy*n_energy_out*(n_mu+2):n_mu+2] - table.inelastic_e_out.shape = (n_energy, n_energy_out) + table.inelastic_e_out[temperatures[0]] = \ + ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2): + n_mu + 2] + table.inelastic_e_out[temperatures[0]].shape = \ + (n_energy, n_energy_out) - table.inelastic_mu_out = ace.xss[idx:idx+n_energy*n_energy_out*(n_mu+2)] - table.inelastic_mu_out.shape = (n_energy, n_energy_out, n_mu+2) - table.inelastic_mu_out = table.inelastic_mu_out[:, :, 1:] + table.inelastic_mu_out[temperatures[0]] = \ + ace.xss[idx:idx + n_energy * n_energy_out * (n_mu + 2)] + table.inelastic_mu_out[temperatures[0]].shape = \ + (n_energy, n_energy_out, n_mu+2) + table.inelastic_mu_out[temperatures[0]] = \ + table.inelastic_mu_out[temperatures[0]][:, :, 1:] else: n_mu = ace.nxs[3] - 1 idx = ace.jxs[3] locc = ace.xss[idx:idx + n_energy].astype(int) - n_energy_out = ace.xss[idx + n_energy:idx + 2*n_energy].astype(int) + n_energy_out = \ + ace.xss[idx + n_energy:idx + 2 * n_energy].astype(int) energy_out = [] mu_out = [] for i in range(n_energy): idx = locc[i] # Outgoing energy distribution for incoming energy i - e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3):n_mu + 3] - p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3):n_mu + 3] - c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3):n_mu + 3] + e = ace.xss[idx + 1:idx + 1 + n_energy_out[i]*(n_mu + 3): + n_mu + 3] + p = ace.xss[idx + 2:idx + 2 + n_energy_out[i]*(n_mu + 3): + n_mu + 3] + c = ace.xss[idx + 3:idx + 3 + n_energy_out[i]*(n_mu + 3): + n_mu + 3] eout_i = Tabular(e, p, 'linear-linear', ignore_negative=True) eout_i.c = c - # Outgoing angle distribution for each (incoming, outgoing) energy pair + # Outgoing angle distribution for each + # (incoming, outgoing) energy pair mu_i = [] for j in range(n_energy_out[i]): mu = ace.xss[idx + 4:idx + 4 + n_mu] - p_mu = 1./n_mu*np.ones(n_mu) + p_mu = 1. / n_mu * np.ones(n_mu) mu_ij = Discrete(mu, p_mu) mu_ij.c = np.cumsum(p_mu) mu_i.append(mu_ij) @@ -398,28 +598,30 @@ class ThermalScattering(EqualityMixin): # Create correlated angle-energy distribution breakpoints = [n_energy] interpolation = [2] - energy = table.inelastic_xs.x - table.inelastic_dist = CorrelatedAngleEnergy( + energy = table.inelastic_xs[temperatures[0]].x + table.inelastic_dist[temperatures[0]] = CorrelatedAngleEnergy( breakpoints, interpolation, energy, energy_out, mu_out) # Incoherent/coherent elastic scattering cross section idx = ace.jxs[4] if idx != 0: n_energy = int(ace.xss[idx]) - energy = ace.xss[idx+1 : idx+1+n_energy] - P = ace.xss[idx+1+n_energy : idx+1+2*n_energy] + energy = ace.xss[idx + 1: idx + 1 + n_energy] + P = ace.xss[idx + 1 + n_energy: idx + 1 + 2 * n_energy] if ace.nxs[5] == 4: - table.elastic_xs = CoherentElastic(energy, P) + table.elastic_xs[temperatures[0]] = CoherentElastic(energy, P) else: - table.elastic_xs = Tabulated1D(energy, P) + table.elastic_xs[temperatures[0]] = Tabulated1D(energy, P) # Angular distribution n_mu = ace.nxs[6] if n_mu != -1: idx = ace.jxs[6] - table.elastic_mu_out = ace.xss[idx:idx + n_energy*n_mu] - table.elastic_mu_out.shape = (n_energy, n_mu) + table.elastic_mu_out[temperatures[0]] = \ + ace.xss[idx:idx + n_energy * n_mu] + table.elastic_mu_out[temperatures[0]].shape = \ + (n_energy, n_mu) # Get relevant ZAIDs pairs = np.fromiter(map(lambda p: p[0], ace.pairs), int) diff --git a/scripts/openmc-ace-to-hdf5 b/scripts/openmc-ace-to-hdf5 index 05e3aff32..a21fc9001 100755 --- a/scripts/openmc-ace-to-hdf5 +++ b/scripts/openmc-ace-to-hdf5 @@ -124,47 +124,79 @@ for filename in ace_libraries: continue lib = openmc.data.ace.Library(filename) + nuclides = {} for table in lib.tables: - if table.name.endswith('c'): + name, xs = table.name.split('.') + if xs.endswith('c'): # Continuous-energy neutron data - try: + if name not in nuclides: neutron = openmc.data.IncidentNeutron.from_ace( - table, args.metastable) - except Exception as e: - print('Failed to convert {}: {}'.format(table.name, e)) - continue + table, args.metastable) + # try: + # neutron = openmc.data.IncidentNeutron.from_ace( + # table, args.metastable) + # except Exception as e: + # print('Failed to convert {}: {}'.format(table.name, e)) + # continue - # Fission energy release data, if available - if args.fission_energy_release is not None: - fer = openmc.data.FissionEnergyRelease.from_compact_hdf5( - args.fission_energy_release, neutron) - if fer is not None: - neutron.fission_energy = fer + # Fission energy release data, if available + if args.fission_energy_release is not None: + fer = openmc.data.FissionEnergyRelease.from_compact_hdf5( + args.fission_energy_release, neutron) + if fer is not None: + neutron.fission_energy = fer - print('Converting {} (ACE) to {} (HDF5)'.format(table.name, - neutron.name)) + print('Converting {} (ACE) to {} (HDF5)'.format(table.name, + neutron.name)) - # Determine filename - outfile = os.path.join(args.destination, - neutron.name.replace('.', '_') + '.h5') - neutron.export_to_hdf5(outfile, 'w') + # Determine filename + outfile = os.path.join(args.destination, + neutron.name.replace('.', '_') + '.h5') + neutron.export_to_hdf5(outfile, 'w') - # Register with library - library.register_file(outfile) + # Register with library + library.register_file(outfile) - elif table.name.endswith('t'): + # Add nuclide to list + nuclides[name] = outfile + else: + # Then we only need to append the data + print('Converting {} (ACE) to {} (HDF5)'.format(table.name, + neutron.name)) + neutron = \ + openmc.data.IncidentNeutron.from_hdf5(nuclides[name]) + neutron.add_temperature_from_ace(table, args.metastable) + neutron.export_to_hdf5(outfile + '_1', 'w') + os.rename(outfile + '_1', outfile) + + elif xs.endswith('t'): # Thermal scattering data - thermal = openmc.data.ThermalScattering.from_ace(table) - print('Converting {} (ACE) to {} (HDF5)'.format(table.name, - thermal.name)) + if name not in nuclides: + thermal = openmc.data.ThermalScattering.from_ace(table) + print('Converting {} (ACE) to {} (HDF5)'.format(table.name, + thermal.name)) - # Determine filename - outfile = os.path.join(args.destination, - thermal.name.replace('.', '_') + '.h5') - thermal.export_to_hdf5(outfile, 'w') + # Determine filename + outfile = os.path.join(args.destination, + thermal.name.replace('.', '_') + '.h5') + thermal.export_to_hdf5(outfile, 'w') - # Register with library - library.register_file(outfile, 'thermal') + # Register with library + library.register_file(outfile, 'thermal') + + # Add data to list + nuclides[name] = outfile + + else: + # Then we only need to append the data + print('Converting {} (ACE) to {} (HDF5)'.format(table.name, + thermal.name)) + # if table.name == 'poly.11t': + # import pdb; pdb.set_trace() + thermal = openmc.data.ThermalScattering.from_hdf5(nuclides[name]) + thermal.add_temperature_from_ace(table) + thermal.export_to_hdf5(outfile + '_1', 'w') + os.rename(outfile + '_1', outfile) # Write cross_sections.xml libpath = os.path.join(args.destination, 'cross_sections.xml')