diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 5f34d4700c..e334f2bf11 100644 --- a/docs/source/io_formats/statepoint.rst +++ b/docs/source/io_formats/statepoint.rst @@ -201,15 +201,28 @@ if run_mode == 'k-eigenvalue': **/tallies/tally /filter /type** (*char[]*) Type of the j-th filter. Can be 'universe', 'material', 'cell', 'cellborn', - 'surface', 'mesh', 'energy', 'energyout', or 'distribcell'. + 'surface', 'mesh', 'energy', 'energyout', 'distribcell', 'mu', 'polar', + 'azimuthal', 'delayedgroup', or 'energyfunction'. **/tallies/tally /filter /n_bins** (*int*) - Number of bins for the j-th filter. + Number of bins for the j-th filter. Not present for 'energyfunction' + filters. **/tallies/tally /filter /bins** (*int[]* or *double[]*) - Value for each filter bin of this type. + Value for each filter bin of this type. Not present for 'energyfunction' + filters. + +**/tallies/tally /filter /energy** (*double[]*) + + Energy grid points for energyfunction interpolation. Only used for + 'energyfunction' filters. + +**/tallies/tally /filter /y** (*double[]*) + + Interpolant values for energyfunction interpolation. Only used for + 'energyfunction' filters. **/tallies/tally /nuclides** (*char[][]*) diff --git a/docs/source/pythonapi/index.rst b/docs/source/pythonapi/index.rst index 857826e097..1125a4b32d 100644 --- a/docs/source/pythonapi/index.rst +++ b/docs/source/pythonapi/index.rst @@ -146,6 +146,7 @@ Constructing Tallies openmc.AzimuthalFilter openmc.DistribcellFilter openmc.DelayedGroupFilter + openmc.EnergyFunctionFilter openmc.Mesh openmc.Trigger openmc.Tally diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index 4a81615a10..09eab32621 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1481,9 +1481,12 @@ The ```` element accepts the following sub-elements: *Default*: "" :filter: - Specify a filter that restricts contributions to the tally to particles - within certain regions of phase space. This element and its - attributes/sub-elements are described below. + Specify a filter that modifies tally behavior. Most tallies (e.g. ``cell``, + ``energy``, and ``material``) restrict the tally so that only particles + within certain regions of phase space contribute to the tally. Others + (e.g. ``delayedgroup`` and ``energyfunction``) can apply some other function + to the scored values. This element and its attributes/sub-elements are + described below. .. note:: You may specify zero, one, or multiple filters to apply to the tally. To @@ -1494,7 +1497,7 @@ The ```` element accepts the following sub-elements: :type: The type of the filter. Accepted options are "cell", "cellborn", "material", "universe", "energy", "energyout", "mesh", "distribcell", - and "delayedgroup". + "delayedgroup", and "energyfunction". :bins: For each filter type, the corresponding ``bins`` entry is given as @@ -1629,6 +1632,23 @@ The ```` element accepts the following sub-elements: + :energyfunction: + ``energyfunction`` filters do not use the ``bins`` entry. Instead + they use ``energy`` and ``y``. + + :energy: + ``energyfunction`` filters multiply tally scores by an arbitrary + function. The function is described by a piecewise linear-linear set of + (energy, y) values. This entry specifies the energy values. The function + will be evaluated as zero outside of the bounds of this energy grid. + (Only used for ``energyfunction`` filters) + + :y: + ``energyfunction`` filters multiply tally scores by an arbitrary + function. The function is described by a piecewise linear-linear set of + (energy, y) values. This entry specifies the y values. (Only used + for ``energyfunction`` filters) + :nuclides: If specified, the scores listed will be for particular nuclides, not the summation of reactions from all nuclides. The format for nuclides should be diff --git a/openmc/filter.py b/openmc/filter.py index 77e409c404..dd3e21e5c8 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1,6 +1,7 @@ from abc import ABCMeta from collections import Iterable, OrderedDict import copy +import hashlib from numbers import Real, Integral from xml.etree import ElementTree as ET @@ -13,7 +14,7 @@ import openmc.checkvalue as cv _FILTER_TYPES = ['universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', 'mu', 'polar', 'azimuthal', - 'distribcell', 'delayedgroup'] + 'distribcell', 'delayedgroup', 'energyfunction'] _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in', 3: 'x-max out', 4: 'x-max in', @@ -25,9 +26,38 @@ _CURRENT_NAMES = {1: 'x-min out', 2: 'x-min in', class FilterMeta(ABCMeta): def __new__(cls, name, bases, namespace, **kwargs): + # Check the class name. if not name.endswith('Filter'): raise ValueError("All filter class names must end with 'Filter'") + + # Create a 'short_name' attribute that removes the 'Filter' suffix. namespace['short_name'] = name[:-6] + + # Subclass methods can sort of inherit the docstring of parent class + # methods. If a function is defined without a docstring, most (all?) + # Python interpreters will search through the parent classes to see if + # there is a docstring for a function with the same name, and they will + # use that docstring. However, Sphinx does not have that functionality. + # This chunk of code handles this docstring inheritance manually so that + # the autodocumentation will pick it up. + if name != 'Filter': + # Look for newly-defined functions that were also in Filter. + for func_name in namespace: + if func_name in Filter.__dict__: + # Inherit the docstring from Filter if not defined. + if isinstance(namespace[func_name], + (classmethod, staticmethod)): + new_doc = namespace[func_name].__func__.__doc__ + old_doc = Filter.__dict__[func_name].__func__.__doc__ + if new_doc is None and old_doc is not None: + namespace[func_name].__func__.__doc__ = old_doc + else: + new_doc = namespace[func_name].__doc__ + old_doc = Filter.__dict__[func_name].__doc__ + if new_doc is None and old_doc is not None: + namespace[func_name].__doc__ = old_doc + + # Make the class. return super(FilterMeta, cls).__new__(cls, name, bases, namespace, **kwargs) @@ -93,7 +123,7 @@ class Filter(object): def __repr__(self): string = type(self).__name__ + '\n' - string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins) + string += '{: <16}=\t{}\n'.format('\tBins', self.bins) return string @classmethod @@ -154,12 +184,8 @@ class Filter(object): @bins.setter def bins(self, bins): - # Make sure the bins are a numpy array. - bins = np.array(bins) - - # If the bin is 0D numpy array, promote to 1D. - if bins.shape == (): - bins.shape = (1,) + # Format the bins as a 1D numpy array. + bins = np.atleast_1d(bins) # Check the bin values. self.check_bins(bins) @@ -194,8 +220,15 @@ class Filter(object): pass - def to_xml(self): - """Return XML Element representing the Filter.""" + def to_xml_element(self): + """Return XML Element representing the Filter. + + Returns + ------- + ElementTree.Element + + """ + element = ET.Element('filter') element.set('type', self.short_name.lower()) element.set('bins', ' '.join(str(b) for b in self.bins)) @@ -1610,3 +1643,238 @@ class DelayedGroupFilter(IntegralFilter): filter's bins. """ + + +class EnergyFunctionFilter(Filter): + """Multiplies tally scores by an arbitrary function of incident energy. + + The arbitrary function is described by a piecewise linear-linear + interpolation of energy and y values. Values outside of the given energy + range will be evaluated as zero. + + Parameters + ---------- + energy : Iterable of Real + A grid of energy values in eV. + y : iterable of Real + A grid of interpolant values in eV. + + Attributes + ---------- + energy : Iterable of Real + A grid of energy values in eV. + y : iterable of Real + A grid of interpolant values in eV. + num_bins : Integral + The number of filter bins (always 1 for this filter) + stride : Integral + The number of filter, nuclide and score bins within each of this + filter's bins. + + """ + + def __init__(self, energy, y): + self.energy = energy + self.y = y + self._stride = None + + def __eq__(self, other): + if type(self) is not type(other): + return False + elif not all(self.energy == other.energy): + return False + elif not all(self.y == other.y): + return False + else: + return True + + def __gt__(self, other): + if type(self) is not type(other): + if self.short_name in _FILTER_TYPES and \ + other.short_name in _FILTER_TYPES: + delta = _FILTER_TYPES.index(self.short_name) - \ + _FILTER_TYPES.index(other.short_name) + return delta > 0 + else: + return False + else: + return False + + def __lt__(self, other): + if type(self) is not type(other): + if self.short_name in _FILTER_TYPES and \ + other.short_name in _FILTER_TYPES: + delta = _FILTER_TYPES.index(self.short_name) - \ + _FILTER_TYPES.index(other.short_name) + return delta < 0 + else: + return False + else: + return False + + def __hash__(self): + # For some reason, it seems the __hash__ method is not inherited when we + # overwrite __repr__. + return hash(repr(self)) + + def __repr__(self): + string = type(self).__name__ + '\n' + string += '{: <16}=\t{}\n'.format('\tEnergy', self.energy) + string += '{: <16}=\t{}\n'.format('\tInterpolant', self.y) + return string + + @classmethod + def from_hdf5(cls, group, **kwargs): + if group['type'].value.decode() != cls.short_name.lower(): + raise ValueError("Expected HDF5 data for filter type '" + + cls.short_name.lower() + "' but got '" + + group['type'].value.decode() + " instead") + + energy = group['energy'].value + y = group['y'].value + + return cls(energy, y) + + @classmethod + def from_tabulated1d(cls, tab1d): + """Construct a filter from a Tabulated1D object. + + Parameters + ---------- + tab1d : openmc.data.Tabulated1D + A linear-linear Tabulated1D object with only a single interpolation + region. + + Returns + ------- + EnergyFunctionFilter + + """ + cv.check_type('EnergyFunctionFilter tab1d', tab1d, + openmc.data.Tabulated1D) + if tab1d.n_regions > 1: + raise ValueError('Only Tabulated1Ds with a single interpolation ' + 'region are supported') + if tab1d.interpolation[0] != 2: + raise ValueError('Only linear-linar Tabulated1Ds are supported') + + return cls(tab1d.x, tab1d.y) + + @property + def energy(self): + return self._energy + + @property + def y(self): + return self._y + + @property + def bins(self): + raise RuntimeError('EnergyFunctionFilters have no bins.') + + @property + def num_bins(self): + return 1 + + @energy.setter + def energy(self, energy): + # Format the bins as a 1D numpy array. + energy = np.atleast_1d(energy) + + # Make sure the values are Real and positive. + cv.check_type('filter energy grid', energy, Iterable, Real) + for E in energy: + cv.check_greater_than('filter energy grid', E, 0, equality=True) + + self._energy = energy + + @y.setter + def y(self, y): + # Format the bins as a 1D numpy array. + y = np.atleast_1d(y) + + # Make sure the values are Real. + cv.check_type('filter interpolant values', y, Iterable, Real) + + self._y = y + + @bins.setter + def bins(self, bins): + raise RuntimeError('EnergyFunctionFilters have no bins.') + + def to_xml_element(self): + element = ET.Element('filter') + element.set('type', self.short_name.lower()) + element.set('energy', ' '.join(str(e) for e in self.energy)) + element.set('y', ' '.join(str(y) for y in self.y)) + return element + + def can_merge(self, other): + return False + + def is_subset(self, other): + return self == other + + def get_bin_index(self, filter_bin): + # This filter only has one bin. Always return 0. + return 0 + + def get_bin(self, bin_index): + """This function is invalid for EnergyFunctionFilters.""" + raise RuntimeError('EnergyFunctionFilters have no get_bin() method') + + def get_pandas_dataframe(self, data_size, **kwargs): + """Builds a Pandas DataFrame for the Filter's bins. + + This method constructs a Pandas DataFrame object for the filter with + columns annotated by filter bin information. This is a helper method for + :meth:`Tally.get_pandas_dataframe`. + + Parameters + ---------- + data_size : Integral + The total number of bins in the tally corresponding to this filter + + Returns + ------- + pandas.DataFrame + A Pandas DataFrame with a column that is filled with a hash of this + filter. EnergyFunctionFilters have only 1 bin so the purpose of this + DataFrame column is to differentiate the filter from other + EnergyFunctionFilters. The number of rows in the DataFrame is the + same as the total number of bins in the corresponding tally. + + Raises + ------ + ImportError + When Pandas is not installed + + See also + -------- + Tally.get_pandas_dataframe(), CrossFilter.get_pandas_dataframe() + + """ + + import pandas as pd + df = pd.DataFrame() + + # There is no clean way of sticking all the energy, y data into a + # DataFrame so instead we'll just make a column with the filter name + # and fill it with a hash of the __repr__. We want a hash that is + # reproducible after restarting the interpreter so we'll use hashlib.md5 + # rather than the intrinsic hash(). + hash_fun = hashlib.md5() + hash_fun.update(repr(self).encode('utf-8')) + out = hash_fun.hexdigest() + + # The full 16 bytes make for a really wide column. Just 7 bytes (14 + # hex characters) of the digest are probably sufficient. + out = out[:14] + + filter_bins = np.repeat(out, self.stride) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + df = pd.concat([df, pd.DataFrame( + {self.short_name.lower(): filter_bins})]) + + return df diff --git a/openmc/tallies.py b/openmc/tallies.py index 9584387de0..a9a1390dfa 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -186,11 +186,8 @@ class Tally(object): string += '{: <16}=\t{}\n'.format('\tDerivative ID', str(self.derivative.id)) - string += '{: <16}=\n'.format('\tFilters') - - for self_filter in self.filters: - string += '{: <16}\t\t{}\t{}\n'.format('', - type(self_filter).__name__, self_filter.bins) + filters = ', '.join(type(f).__name__ for f in self.filters) + string += '{: <16}=\t{}\n'.format('\tFilters', filters) string += '{: <16}=\t'.format('\tNuclides') @@ -1039,7 +1036,7 @@ class Tally(object): # Optional Tally filters for self_filter in self.filters: - element.append(self_filter.to_xml()) + element.append(self_filter.to_xml_element()) # Optional Nuclides if len(self.nuclides) > 0: @@ -1330,6 +1327,10 @@ class Tally(object): elif isinstance(self_filter, openmc.DistribcellFilter): bins = [b for b in range(self_filter.num_bins)] + # EnergyFunctionFilters don't have bins so just add a None + elif isinstance(self_filter, openmc.EnergyFunctionFilter): + bins = [None] + # Create list of IDs for bins for all other filter types else: bins = self_filter.bins @@ -1703,146 +1704,6 @@ class Tally(object): data = np.reshape(data, new_shape) return data - def export_results(self, filename='tally-results', directory='.', - format='hdf5', append=True): - """Exports tallly results to an HDF5 or Python pickle binary file. - - Parameters - ---------- - filename : str - The name of the file for the results (default is 'tally-results') - directory : str - The name of the directory for the results (default is '.') - format : str - The format for the exported file - HDF5 ('hdf5', default) and - Python pickle ('pkl') files are supported - append : bool - Whether or not to append the results to the file (default is True) - - Raises - ------ - KeyError - When this method is called before the Tally is populated with data. - - """ - - # Ensure that the tally has data - if self._sum is None or self._sum_sq is None and not self.derived: - msg = 'The Tally ID="{0}" has no data to export'.format(self.id) - raise KeyError(msg) - - if not isinstance(filename, string_types): - msg = 'Unable to export the results for Tally ID="{0}" to ' \ - 'filename="{1}" since it is not a ' \ - 'string'.format(self.id, filename) - raise ValueError(msg) - - elif not isinstance(directory, string_types): - msg = 'Unable to export the results for Tally ID="{0}" to ' \ - 'directory="{1}" since it is not a ' \ - 'string'.format(self.id, directory) - raise ValueError(msg) - - elif format not in ['hdf5', 'pkl', 'csv']: - msg = 'Unable to export the results for Tally ID="{0}" to format ' \ - '"{1}" since it is not supported'.format(self.id, format) - raise ValueError(msg) - - elif not isinstance(append, bool): - msg = 'Unable to export the results for Tally ID="{0}" since the ' \ - 'append parameter is not True/False'.format(self.id) - raise ValueError(msg) - - # Make directory if it does not exist - if not os.path.exists(directory): - os.makedirs(directory) - - # HDF5 binary file - if format == 'hdf5': - filename = directory + '/' + filename + '.h5' - - if append: - tally_results = h5py.File(filename, 'a') - else: - tally_results = h5py.File(filename, 'w') - - # Create an HDF5 group within the file for this particular Tally - tally_group = tally_results.create_group('Tally-{0}'.format(self.id)) - - # Add basic Tally data to the HDF5 group - tally_group.create_dataset('id', data=self.id) - tally_group.create_dataset('name', data=self.name) - tally_group.create_dataset('estimator', data=self.estimator) - tally_group.create_dataset('scores', data=np.array(self.scores)) - - # Add a string array of the nuclides to the HDF5 group - nuclides = [] - - for nuclide in self.nuclides: - nuclides.append(nuclide.name) - - tally_group.create_dataset('nuclides', data=np.array(nuclides)) - - # Create an HDF5 sub-group for the Filters - filter_group = tally_group.create_group('filters') - - for self_filter in self.filters: - filter_group.create_dataset(self_filter.type, - filter=self_filter.bins) - - # Add all results to the main HDF5 group for the Tally - tally_group.create_dataset('sum', data=self.sum) - tally_group.create_dataset('sum_sq', data=self.sum_sq) - tally_group.create_dataset('mean', data=self.mean) - tally_group.create_dataset('std_dev', data=self.std_dev) - - # Close the Tally results HDF5 file - tally_results.close() - - # Python pickle binary file - elif format == 'pkl': - # Load the dictionary from the Pickle file - filename = directory + '/' + filename + '.pkl' - - if os.path.exists(filename) and append: - tally_results = pickle.load(open(filename, 'rb')) - else: - tally_results = {} - - # Create a nested dictionary within the file for this particular Tally - tally_results['Tally-{0}'.format(self.id)] = {} - tally_group = tally_results['Tally-{0}'.format(self.id)] - - # Add basic Tally data to the nested dictionary - tally_group['id'] = self.id - tally_group['name'] = self.name - tally_group['estimator'] = self.estimator - tally_group['scores'] = np.array(self.scores) - - # Add a string array of the nuclides to the HDF5 group - nuclides = [] - - for nuclide in self.nuclides: - nuclides.append(nuclide.name) - - tally_group['nuclides'] = np.array(nuclides) - - # Create a nested dictionary for the Filters - tally_group['filters'] = {} - filter_group = tally_group['filters'] - - for self_filter in self.filters: - filter_group[self_filter.type] = self_filter.bins - - # Add all results to the main sub-dictionary for the Tally - tally_group['sum'] = self.sum - tally_group['sum_sq'] = self.sum_sq - tally_group['mean'] = self.mean - tally_group['std_dev'] = self.std_dev - - # Pickle the Tally results to a file - pickle.dump(tally_results, open(filename, 'wb')) - def hybrid_product(self, other, binary_op, filter_product=None, nuclide_product=None, score_product=None): """Combines filters, scores and nuclides with another tally. @@ -2260,11 +2121,15 @@ class Tally(object): filters = [type(filter1), type(filter2)] if isinstance(filter1, openmc.DistribcellFilter): filter1_bins = [b for b in range(filter1.num_bins)] + elif isinstance(filter1, openmc.EnergyFunctionFilter): + filter1_bins = [None] else: filter1_bins = [filter1.get_bin(i) for i in range(filter1.num_bins)] if isinstance(filter2, openmc.DistribcellFilter): filter2_bins = [b for b in range(filter2.num_bins)] + elif isinstance(filter2, openmc.EnergyFunctionFilter): + filter2_bins = [None] else: filter2_bins = [filter2.get_bin(i) for i in range(filter2.num_bins)] @@ -3096,6 +2961,8 @@ class Tally(object): if isinstance(find_filter, openmc.DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) + elif isinstance(find_filter, openmc.EnergyFunctionFilter): + filter_bins = [None] else: num_bins = find_filter.num_bins filter_bins = \ @@ -3251,6 +3118,8 @@ class Tally(object): if isinstance(find_filter, openmc.DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) + elif isinstance(find_filter, openmc.EnergyFunctionFilter): + filter_bins = [None] else: num_bins = find_filter.num_bins filter_bins = \ diff --git a/src/constants.F90 b/src/constants.F90 index 044ea1c6bc..2c27ce3ee1 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -342,21 +342,22 @@ module constants integer, parameter :: NO_BIN_FOUND = -1 ! Tally filter and map types - integer, parameter :: N_FILTER_TYPES = 13 + integer, parameter :: N_FILTER_TYPES = 14 integer, parameter :: & - FILTER_UNIVERSE = 1, & - FILTER_MATERIAL = 2, & - FILTER_CELL = 3, & - FILTER_CELLBORN = 4, & - FILTER_SURFACE = 5, & - FILTER_MESH = 6, & - FILTER_ENERGYIN = 7, & - FILTER_ENERGYOUT = 8, & - FILTER_DISTRIBCELL = 9, & - FILTER_MU = 10, & - FILTER_POLAR = 11, & - FILTER_AZIMUTHAL = 12, & - FILTER_DELAYEDGROUP = 13 + FILTER_UNIVERSE = 1, & + FILTER_MATERIAL = 2, & + FILTER_CELL = 3, & + FILTER_CELLBORN = 4, & + FILTER_SURFACE = 5, & + FILTER_MESH = 6, & + FILTER_ENERGYIN = 7, & + FILTER_ENERGYOUT = 8, & + FILTER_DISTRIBCELL = 9, & + FILTER_MU = 10, & + FILTER_POLAR = 11, & + FILTER_AZIMUTHAL = 12, & + FILTER_DELAYEDGROUP = 13, & + FILTER_ENERGYFUNCTION = 14 ! Mesh types integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1bbf68f48b..dd8be29ea9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -2973,18 +2973,21 @@ contains temp_str = to_lower(temp_str) ! Determine number of bins - if (check_for_node(node_filt, "bins")) then - if (temp_str == 'energy' .or. temp_str == 'energyout' .or. & - temp_str == 'mu' .or. temp_str == 'polar' .or. & - temp_str == 'azimuthal') then - n_words = get_arraysize_double(node_filt, "bins") - else - n_words = get_arraysize_integer(node_filt, "bins") + select case(temp_str) + case ("energy", "energyout", "mu", "polar", "azimuthal") + if (.not. check_for_node(node_filt, "bins")) then + call fatal_error("Bins not set in filter on tally " & + // trim(to_str(t % id))) end if - else - call fatal_error("Bins not set in filter on tally " & - // trim(to_str(t % id))) - end if + n_words = get_arraysize_double(node_filt, "bins") + case ("mesh", "universe", "material", "cell", "distribcell", & + "cellborn", "surface", "delayedgroup") + if (.not. check_for_node(node_filt, "bins")) then + call fatal_error("Bins not set in filter on tally " & + // trim(to_str(t % id))) + end if + n_words = get_arraysize_integer(node_filt, "bins") + end select ! Determine type of filter select case (temp_str) @@ -3282,6 +3285,39 @@ contains ! Set the filter index in the tally find_filter array t % find_filter(FILTER_AZIMUTHAL) = j + case ('energyfunction') + ! Allocate and declare the filter type. + allocate(EnergyFunctionFilter::t % filters(j) % obj) + select type (filt => t % filters(j) % obj) + type is (EnergyFunctionFilter) + filt % n_bins = 1 + ! Make sure this is continuous-energy mode. + if (.not. run_CE) then + call fatal_error("EnergyFunction filters are only supported for & + &continuous-energy transport calculations") + end if + + ! Allocate and store energy grid. + if (.not. check_for_node(node_filt, "energy")) then + call fatal_error("Energy grid not specified for EnergyFunction & + &filter on tally " // trim(to_str(t % id))) + end if + n_words = get_arraysize_double(node_filt, "energy") + allocate(filt % energy(n_words)) + call get_node_array(node_filt, "energy", filt % energy) + + ! Allocate and store interpolant values. + if (.not. check_for_node(node_filt, "y")) then + call fatal_error("y values not specified for EnergyFunction & + &filter on tally " // trim(to_str(t % id))) + end if + n_words = get_arraysize_double(node_filt, "y") + allocate(filt % y(n_words)) + call get_node_array(node_filt, "y", filt % y) + end select + ! Set the filter index in the tally find_filter array + t % find_filter(FILTER_ENERGYFUNCTION) = j + case default ! Specified tally filter is invalid, raise error call fatal_error("Unknown filter type '" & diff --git a/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 2143c8e472..57f5f37478 100644 --- a/src/relaxng/tallies.rnc +++ b/src/relaxng/tallies.rnc @@ -25,8 +25,8 @@ element tallies { | attribute variable { ( "nuclide_density" ) } ) & (element nuclide { xsd:string { maxLength = "12" } } - | attribute nuclide { xsd:string { maxLength = "12" } } ) | - ) + | attribute nuclide { xsd:string { maxLength = "12" } } ) + ) | (element variable { ( "temperature") } | attribute variable { ( "temperature" ) } ) ) @@ -39,14 +39,24 @@ element tallies { (element estimator { ( "analog" | "tracklength" | "collision" ) } | attribute estimator { ( "analog" | "tracklength" | "collision" ) })? & element filter { - (element type { ( "cell" | "cellborn" | "material" | "universe" | - "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | - "polar" | "azimuthal" | "delayedgroup") } | - attribute type { ( "cell" | "cellborn" | "material" | "universe" | - "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | - "polar" | "azimuthal" | "delayedgroup") }) & - (element bins { list { xsd:double+ } } | - attribute bins { list { xsd:double+ } }) + ( + (element type { ( "cell" | "cellborn" | "material" | "universe" | + "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | + "polar" | "azimuthal" | "delayedgroup" | "energyfunction") } | + attribute type { ( "cell" | "cellborn" | "material" | "universe" | + "surface" | "distribcell" | "mesh" | "energy" | "energyout" | "mu" | + "polar" | "azimuthal" | "delayedgroup" | "energyfunction") }) & + (element bins { list { xsd:double+ } } | + attribute bins { list { xsd:double+ } }) + ) | + ( + (element type { ("energyfunction") } | + attribute type { ("energyfunction") }) & + (element energy { list { xsd:double+ } } | + attribute energy { list { xsd:double+ } }) & + (element y { list { xsd:double+ } } | + attribute y { list { xsd:double+ } }) + ) }* & element nuclides { list { xsd:string { maxLength = "12" }+ } diff --git a/src/relaxng/tallies.rng b/src/relaxng/tallies.rng index fc60b8373a..fde06bac6c 100644 --- a/src/relaxng/tallies.rng +++ b/src/relaxng/tallies.rng @@ -196,60 +196,106 @@ - - - - - cell - cellborn - material - universe - surface - distribcell - mesh - energy - energyout - mu - polar - azimuthal - delayedgroup - - - - - cell - cellborn - material - universe - surface - distribcell - mesh - energy - energyout - mu - polar - azimuthal - delayedgroup - - - - - - - - - - - - - - - - - - - - + + + + + + cell + cellborn + material + universe + surface + distribcell + mesh + energy + energyout + mu + polar + azimuthal + delayedgroup + energyfunction + + + + + cell + cellborn + material + universe + surface + distribcell + mesh + energy + energyout + mu + polar + azimuthal + delayedgroup + energyfunction + + + + + + + + + + + + + + + + + + + + + + + + energyfunction + + + energyfunction + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 37265cc5ff..4b372defe0 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -1,7 +1,7 @@ module tally_filter use algorithm, only: binary_search - use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION + use constants, only: ONE, NO_BIN_FOUND, FP_PRECISION, ERROR_REAL use dict_header, only: DictIntInt use geometry_header, only: BASE_UNIVERSE, RectLattice, HexLattice use global @@ -191,6 +191,20 @@ module tally_filter procedure :: text_label => text_label_azimuthal end type AzimuthalFilter +!=============================================================================== +! EnergyFunctionFilter multiplies tally scores by an arbitrary function of +! incident energy described by a piecewise linear-linear interpolation. +!=============================================================================== + type, extends(TallyFilter) :: EnergyFunctionFilter + real(8), allocatable :: energy(:) + real(8), allocatable :: y(:) + + contains + procedure :: get_next_bin => get_next_bin_energyfunction + procedure :: to_statepoint => to_statepoint_energyfunction + procedure :: text_label => text_label_energyfunction + end type EnergyFunctionFilter + contains !=============================================================================== @@ -229,6 +243,8 @@ contains logical :: end_in_mesh ! ending coordinates inside mesh? type(RegularMesh), pointer :: m + weight = ERROR_REAL + ! Get a pointer to the mesh. m => meshes(this % mesh) @@ -237,10 +253,10 @@ contains ! valid mesh bin. if (current_bin == NO_BIN_FOUND) then call get_mesh_bin(m, p % coord(1) % xyz, next_bin) + weight = ONE else next_bin = NO_BIN_FOUND end if - weight = ONE else ! A track can span multiple mesh bins so we need to handle a lot of @@ -471,13 +487,14 @@ contains ! Starting one coordinate level deeper, find the next bin. next_bin = NO_BIN_FOUND + weight = ERROR_REAL do i = start, p % n_coord if (this % map % has_key(p % coord(i) % universe)) then next_bin = this % map % get_key(p % coord(i) % universe) + weight = ONE exit end if end do - weight = ONE end subroutine get_next_bin_universe subroutine to_statepoint_universe(this, filter_group) @@ -532,12 +549,13 @@ contains real(8), intent(out) :: weight next_bin = NO_BIN_FOUND + weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then if (this % map % has_key(p % material)) then next_bin = this % map % get_key(p % material) + weight = ONE end if end if - weight = ONE end subroutine get_next_bin_material subroutine to_statepoint_material(this, filter_group) @@ -607,13 +625,14 @@ contains ! Starting one coordinate level deeper, find the next bin. next_bin = NO_BIN_FOUND + weight = ERROR_REAL do i = start, p % n_coord if (this % map % has_key(p % coord(i) % cell)) then next_bin = this % map % get_key(p % coord(i) % cell) + weight = ONE exit end if end do - weight = ONE end subroutine get_next_bin_cell subroutine to_statepoint_cell(this, filter_group) @@ -667,10 +686,8 @@ contains integer, intent(out) :: next_bin real(8), intent(out) :: weight - logical :: cell_found integer :: distribcell_index, offset, i - cell_found = .false. if (current_bin == NO_BIN_FOUND) then distribcell_index = cells(this % cell) % distribcell_index offset = 0 @@ -693,15 +710,13 @@ contains end if if (this % cell == p % coord(i) % cell) then next_bin = offset + 1 - cell_found = .true. - exit + weight = ONE + return end if end do - else - next_bin = NO_BIN_FOUND end if - if (.not. cell_found) next_bin = NO_BIN_FOUND - weight = ONE + next_bin = NO_BIN_FOUND + weight = ERROR_REAL end subroutine get_next_bin_distribcell subroutine to_statepoint_distribcell(this, filter_group) @@ -756,12 +771,13 @@ contains real(8), intent(out) :: weight next_bin = NO_BIN_FOUND + weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then if (this % map % has_key(p % cell_born)) then next_bin = this % map % get_key(p % cell_born) + weight = ONE end if end if - weight = ONE end subroutine get_next_bin_cellborn subroutine to_statepoint_cellborn(this, filter_group) @@ -818,15 +834,16 @@ contains integer :: i next_bin = NO_BIN_FOUND + weight = ERROR_REAL if (current_bin == NO_BIN_FOUND) then do i = 1, this % n_bins if (p % surface == this % surfaces(i)) then next_bin = i + weight = ONE exit end if end do end if - weight = ONE end subroutine get_next_bin_surface subroutine to_statepoint_surface(this, filter_group) @@ -891,6 +908,7 @@ contains ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch next_bin = num_energy_groups - next_bin + 1 + weight = ONE else ! Make sure the correct energy is used. @@ -903,16 +921,18 @@ contains ! Check if energy of the particle is within energy bins. if (E < this % bins(1) .or. E > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find incoming energy bin. next_bin = binary_search(this % bins, n + 1, E) + weight = ONE end if end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_energy subroutine to_statepoint_energy(this, filter_group) @@ -960,21 +980,24 @@ contains ! Tallies are ordered in increasing groups, group indices ! however are the opposite, so switch next_bin = num_energy_groups - next_bin + 1 + weight = ONE else ! Check if energy of the particle is within energy bins. if (p % E < this % bins(1) .or. p % E > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find incoming energy bin. next_bin = binary_search(this % bins, n + 1, p % E) + weight = ONE end if end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_energyout subroutine to_statepoint_energyout(this, filter_group) @@ -1012,10 +1035,11 @@ contains if (current_bin == NO_BIN_FOUND) then next_bin = 1 + weight = ONE else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_dg subroutine to_statepoint_dg(this, filter_group) @@ -1054,15 +1078,17 @@ contains ! Check if energy of the particle is within energy bins. if (p % mu < this % bins(1) .or. p % mu > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find incoming energy bin. next_bin = binary_search(this % bins, n + 1, p % mu) + weight = ONE end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_mu subroutine to_statepoint_mu(this, filter_group) @@ -1115,15 +1141,17 @@ contains ! Check if particle is within polar angle bins. if (theta < this % bins(1) .or. theta > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find polar angle bin. next_bin = binary_search(this % bins, n + 1, theta) + weight = ONE end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_polar subroutine to_statepoint_polar(this, filter_group) @@ -1176,15 +1204,17 @@ contains ! Check if particle is within azimuthal angle bins. if (phi < this % bins(1) .or. phi > this % bins(n + 1)) then next_bin = NO_BIN_FOUND + weight = ERROR_REAL else ! Search to find azimuthal angle bin. next_bin = binary_search(this % bins, n + 1, phi) + weight = ONE end if else next_bin = NO_BIN_FOUND + weight = ERROR_REAL end if - weight = ONE end subroutine get_next_bin_azimuthal subroutine to_statepoint_azimuthal(this, filter_group) @@ -1209,6 +1239,84 @@ contains // trim(to_str(E1)) // ")" end function text_label_azimuthal +!=============================================================================== +! EnergyFunctionFilter methods +!=============================================================================== + subroutine get_next_bin_energyfunction(this, p, estimator, current_bin, & + next_bin, weight) + class(EnergyFunctionFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + integer, value, intent(in) :: current_bin + integer, intent(out) :: next_bin + real(8), intent(out) :: weight + + integer :: n, indx + real(8) :: E, f + + select type(this) + type is (EnergyFunctionFilter) + if (current_bin == NO_BIN_FOUND) then + n = size(this % energy) + + ! Make sure the correct energy is used. + if (estimator == ESTIMATOR_TRACKLENGTH) then + E = p % E + else + E = p % last_E + end if + + ! Check if energy of the particle is within energy bins. + if (E < this % energy(1) .or. E > this % energy(n)) then + next_bin = NO_BIN_FOUND + weight = ERROR_REAL + + else + ! Search to find incoming energy bin. + indx = binary_search(this % energy, n, E) + + ! Compute an interpolation factor between nearest bins. + f = (E - this % energy(indx)) & + / (this % energy(indx+1) - this % energy(indx)) + + ! Interpolate on the lin-lin grid. + next_bin = 1 + weight = (ONE - f) * this % y(indx) + f * this % y(indx+1) + end if + + else + next_bin = NO_BIN_FOUND + weight = ERROR_REAL + end if + end select + end subroutine get_next_bin_energyfunction + + subroutine to_statepoint_energyfunction(this, filter_group) + class(EnergyFunctionFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + select type(this) + type is (EnergyFunctionFilter) + call write_dataset(filter_group, "type", "energyfunction") + call write_dataset(filter_group, "energy", this % energy) + call write_dataset(filter_group, "y", this % y) + end select + end subroutine to_statepoint_energyfunction + + function text_label_energyfunction(this, bin) result(label) + class(EnergyFunctionFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + select type(this) + type is (EnergyFunctionFilter) + write(label, FMT="(A, ES8.1, A, ES8.1, A, ES8.1, A, ES8.1, A)") & + "Energy Function f([", this % energy(1), ", ..., ", & + this % energy(size(this % energy)), "]) = [", this % y(1), & + ", ..., ", this % y(size(this % y)), "]" + end select + end function text_label_energyfunction + !=============================================================================== ! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and ! a target offset to build a string which is the path from the base universe to diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat new file mode 100644 index 0000000000..c9961071ae --- /dev/null +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -0,0 +1 @@ +2a08e03d9e0f40486fa26ba00800d264f43906a38c851da48f0fd76bacc0b862e8f45f16c970b7c62db03d96e3acb1b58e74179f554a851ccaaf4492e97fc163 \ No newline at end of file diff --git a/tests/test_filter_energyfun/results_true.dat b/tests/test_filter_energyfun/results_true.dat new file mode 100644 index 0000000000..d3dc3ad6a5 --- /dev/null +++ b/tests/test_filter_energyfun/results_true.dat @@ -0,0 +1,2 @@ + energyfunction nuclide score mean std. dev. +0 f4ae05cee02fae Am241 (n,gamma) 1.00e-01 6.22e-03 diff --git a/tests/test_filter_energyfun/test_filter_energyfun.py b/tests/test_filter_energyfun/test_filter_energyfun.py new file mode 100644 index 0000000000..b7fa68fa4d --- /dev/null +++ b/tests/test_filter_energyfun/test_filter_energyfun.py @@ -0,0 +1,64 @@ +#!/usr/bin/env python + +import os +import sys +import glob +import hashlib +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + + +class FilterEnergyFunHarness(PyAPITestHarness): + def _build_inputs(self): + # Build the default material, geometry, settings inputs. + self._input_set.build_default_materials_and_geometry() + self._input_set.build_default_settings() + + # Add Am241 to the fuel. + self._input_set.materials[1].add_nuclide('Am241', 1e-7) + + # Define Am242m / Am242 branching ratio from ENDF/B-VII.1 data. + x = [1e-5, 3.69e-1, 1e3, 1e5, 6e5, 1e6, 2e6, 4e6, 3e7] + y = [0.1, 0.1, 0.1333, 0.158, 0.18467, 0.25618, 0.4297, 0.48, 0.48] + + # Make an EnergyFunctionFilter directly from the x and y lists. + filt1 = openmc.EnergyFunctionFilter(x, y) + + # Also make a filter with the .from_tabulated1d constructor. Make sure + # the filters are identical. + tab1d = openmc.data.Tabulated1D(x, y) + filt2 = openmc.EnergyFunctionFilter.from_tabulated1d(tab1d) + assert filt1 == filt2, 'Error with the .from_tabulated1d constructor' + + # Make tallies. + tallies = [openmc.Tally(), openmc.Tally()] + for t in tallies: + t.scores = ['(n,gamma)'] + t.nuclides = ['Am241'] + tallies[1].filters = [filt1] + self._input_set.tallies = openmc.Tallies(tallies) + + # Export inputs to xml. + self._input_set.export() + + def _get_results(self): + # Read the statepoint file. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + sp = openmc.StatePoint(statepoint) + + # Use tally arithmetic to compute the branching ratio. + br_tally = sp.tallies[10001] / sp.tallies[10000] + + # Output the tally in a Pandas DataFrame. + return br_tally.get_pandas_dataframe().to_string() + '\n' + + def _cleanup(self): + super(FilterEnergyFunHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + + +if __name__ == '__main__': + harness = FilterEnergyFunHarness('statepoint.10.h5', True) + harness.main()