From 429667e8b8a3e623360f30aa2f52121a0da45e83 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 25 Nov 2016 02:34:49 -0500 Subject: [PATCH 01/10] Implement LinLinEnergyFilter in Fortran --- openmc/filter.py | 76 +++++++++++++++++++++++++++++++++++++ src/constants.F90 | 5 ++- src/input_xml.F90 | 58 ++++++++++++++++++++++------ src/tally_filter.F90 | 90 ++++++++++++++++++++++++++++++++++++++++++++ 4 files changed, 216 insertions(+), 13 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 88a88f99d..b9e4c4286 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1436,3 +1436,79 @@ class DelayedGroupFilter(IntegralFilter): filter's bins. """ + + +class LinLinEnergyFilter(Filter): + """ + + 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 + + @property + def energy(self): + return self._energy + + @property + def y(self): + return self._y + + @energy.setter + def energy(self, energy): + # Make sure the energy grid is a numpy array. + energy = np.array(energy) + + # If the grid is 0D numpy array, promote to 1D. + if energy.shape == (): + energy.shape = (1,) + + # 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): + # Make sure the values are in a numpy array. + y = np.array(y) + + # If the array is 0D, promote to 1D. + if y.shape == (): + y.shape = (1,) + + # Make sure the values are Real. + cv.check_type('filter interpolant values', y, Iterable, Real) + + self._y = y + + def to_xml(self): + """Return XML Element representing the Filter.""" + 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 diff --git a/src/constants.F90 b/src/constants.F90 index 044ea1c6b..bd9a0508b 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -342,7 +342,7 @@ 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, & @@ -356,7 +356,8 @@ module constants FILTER_MU = 10, & FILTER_POLAR = 11, & FILTER_AZIMUTHAL = 12, & - FILTER_DELAYEDGROUP = 13 + FILTER_DELAYEDGROUP = 13, & + FILTER_LINLINENERGY = 14 ! Mesh types integer, parameter :: & diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1bbf68f48..9279a48b8 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 ('linlinenergy') + ! Allocate and declare the filter type. + allocate(LinLinEnergyFilter::t % filters(j) % obj) + select type (filt => t % filters(j) % obj) + type is (LinLinEnergyFilter) + filt % n_bins = 1 + ! Make sure this is continuous-energy mode. + if (.not. run_CE) then + call fatal_error("LinLinEnergy 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 LinLinEnergy & + &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 LinLinEnergy & + &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_LINLINENERGY) = j + case default ! Specified tally filter is invalid, raise error call fatal_error("Unknown filter type '" & diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 37265cc5f..064e9d646 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -191,6 +191,19 @@ module tally_filter procedure :: text_label => text_label_azimuthal end type AzimuthalFilter +!=============================================================================== +! LinLinEnergyFilter +!=============================================================================== + type, extends(TallyFilter) :: LinLinEnergyFilter + real(8), allocatable :: energy(:) + real(8), allocatable :: y(:) + + contains + procedure :: get_next_bin => get_next_bin_linlinenergy + procedure :: to_statepoint => to_statepoint_linlinenergy + procedure :: text_label => text_label_linlinenergy + end type LinLinEnergyFilter + contains !=============================================================================== @@ -1209,6 +1222,83 @@ contains // trim(to_str(E1)) // ")" end function text_label_azimuthal +!=============================================================================== +! LinLinEnergyFilter methods +!=============================================================================== + subroutine get_next_bin_linlinenergy(this, p, estimator, current_bin, & + next_bin, weight) + class(LinLinEnergyFilter), 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 (LinLinEnergyFilter) + 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 = ONE + + 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 = ONE + end if + end select + end subroutine get_next_bin_linlinenergy + + subroutine to_statepoint_linlinenergy(this, filter_group) + class(LinLinEnergyFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + select type(this) + type is (LinLinEnergyFilter) + call write_dataset(filter_group, "type", "linlinenergy") + call write_dataset(filter_group, "energy", this % energy) + call write_dataset(filter_group, "y", this % y) + end select + end subroutine to_statepoint_linlinenergy + + function text_label_linlinenergy(this, bin) result(label) + class(LinLinEnergyFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + select type(this) + type is (LinLinEnergyFilter) + label = "Lin-Lin Energy Interpolation [" & + // trim(to_str(this % energy(1))) // ", ..., " & + // trim(to_str(this % energy(size(this % energy)))) // "]" + end select + end function text_label_linlinenergy + !=============================================================================== ! 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 From ce0a28de136f6293d584220c9742e10f8fc65cf8 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 4 Dec 2016 16:42:11 -0500 Subject: [PATCH 02/10] Implement LinLinEnergyFilter in PyAPI --- openmc/filter.py | 129 +++++++++++++++++++++++++++++++++++++++++++++- openmc/tallies.py | 19 +++++-- 2 files changed, 142 insertions(+), 6 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index b9e4c4286..318a56fac 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -94,7 +94,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 @@ -1467,6 +1467,63 @@ class LinLinEnergyFilter(Filter): 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) + @property def energy(self): return self._energy @@ -1475,6 +1532,14 @@ class LinLinEnergyFilter(Filter): def y(self): return self._y + @property + def bins(self): + raise RuntimeError('LinLinEnergyFilters have no bins.') + + @property + def num_bins(self): + return 1 + @energy.setter def energy(self, energy): # Make sure the energy grid is a numpy array. @@ -1505,6 +1570,10 @@ class LinLinEnergyFilter(Filter): self._y = y + @bins.setter + def bins(self, bins): + raise RuntimeError('LinLinEnergyFilters have no bins.') + def to_xml(self): """Return XML Element representing the Filter.""" element = ET.Element('filter') @@ -1512,3 +1581,61 @@ class LinLinEnergyFilter(Filter): 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 LinLinEnergyFilters.""" + raise RuntimeError('LinLinEnergyFilters 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. LinLinEnergyFilters have only 1 bin so the purpose of this + DataFrame column is to differentiate the filter from other + LinLinEnergyFilters. 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() + + """ + + # Initialize Pandas DataFrame + import pandas as pd + df = pd.DataFrame() + + filter_bins = np.repeat(hash(self), 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 d5d6df03f..b7d33b13f 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') @@ -1329,6 +1326,10 @@ class Tally(object): elif isinstance(self_filter, openmc.DistribcellFilter): bins = np.arange(self_filter.num_bins) + # LinLinEnergyFilters don't have bins so just add a None + elif isinstance(self_filter, openmc.LinLinEnergyFilter): + bins = [None] + # Create list of IDs for bins for all other filter types else: bins = self_filter.bins @@ -2259,11 +2260,15 @@ class Tally(object): filters = [type(filter1), type(filter2)] if isinstance(filter1, openmc.DistribcellFilter): filter1_bins = np.arange(filter1.num_bins) + elif isinstance(filter1, openmc.LinLinEnergyFilter): + filter1_bins = [None] else: filter1_bins = [filter1.get_bin(i) for i in range(filter1.num_bins)] if isinstance(filter2, openmc.DistribcellFilter): filter2_bins = np.arange(filter2.num_bins) + elif isinstance(filter2, openmc.LinLinEnergyFilter): + filter2_bins = [None] else: filter2_bins = [filter2.get_bin(i) for i in range(filter2.num_bins)] @@ -3095,6 +3100,8 @@ class Tally(object): if isinstance(find_filter, openmc.DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) + elif isinstance(find_filter, openmc.LinLinEnergyFilter): + filter_bins = [None] else: num_bins = find_filter.num_bins filter_bins = \ @@ -3242,6 +3249,8 @@ class Tally(object): if isinstance(find_filter, openmc.DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) + elif isinstance(find_filter, openmc.LinLinEnergyFilter): + filter_bins = [None] else: num_bins = find_filter.num_bins filter_bins = \ From 09a35ea1cdce15e0a09ed7cf5f6e94c0da95f2b5 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 4 Dec 2016 16:44:20 -0500 Subject: [PATCH 03/10] Remove Tally.export_results() method --- openmc/tallies.py | 140 ---------------------------------------------- 1 file changed, 140 deletions(-) diff --git a/openmc/tallies.py b/openmc/tallies.py index b7d33b13f..8eb654b4a 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1703,146 +1703,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. From 09689c4194173cf80d5782d71b0f96688860f4c3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 4 Dec 2016 17:37:18 -0500 Subject: [PATCH 04/10] Change LinLinEnergy to more generic EnergyFunction --- openmc/filter.py | 23 +++++++++------- openmc/tallies.py | 12 ++++----- src/constants.F90 | 28 +++++++++---------- src/input_xml.F90 | 14 +++++----- src/tally_filter.F90 | 64 +++++++++++++++++++++++--------------------- 5 files changed, 73 insertions(+), 68 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 318a56fac..ca761db27 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -14,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', @@ -1438,8 +1438,11 @@ class DelayedGroupFilter(IntegralFilter): """ -class LinLinEnergyFilter(Filter): - """ +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. Parameters ---------- @@ -1534,7 +1537,7 @@ class LinLinEnergyFilter(Filter): @property def bins(self): - raise RuntimeError('LinLinEnergyFilters have no bins.') + raise RuntimeError('EnergyFunctionFilters have no bins.') @property def num_bins(self): @@ -1572,7 +1575,7 @@ class LinLinEnergyFilter(Filter): @bins.setter def bins(self, bins): - raise RuntimeError('LinLinEnergyFilters have no bins.') + raise RuntimeError('EnergyFunctionFilters have no bins.') def to_xml(self): """Return XML Element representing the Filter.""" @@ -1593,8 +1596,8 @@ class LinLinEnergyFilter(Filter): return 0 def get_bin(self, bin_index): - """This function is invalid for LinLinEnergyFilters.""" - raise RuntimeError('LinLinEnergyFilters have no get_bin() method') + """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. @@ -1612,10 +1615,10 @@ class LinLinEnergyFilter(Filter): ------- pandas.DataFrame A Pandas DataFrame with a column that is filled with a hash of this - filter. LinLinEnergyFilters have only 1 bin so the purpose of this + filter. EnergyFunctionFilters have only 1 bin so the purpose of this DataFrame column is to differentiate the filter from other - LinLinEnergyFilters. The number of rows in the DataFrame is the same - as the total number of bins in the corresponding tally. + EnergyFunctionFilters. The number of rows in the DataFrame is the + same as the total number of bins in the corresponding tally. Raises ------ diff --git a/openmc/tallies.py b/openmc/tallies.py index 8eb654b4a..fe690e2d6 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1326,8 +1326,8 @@ class Tally(object): elif isinstance(self_filter, openmc.DistribcellFilter): bins = np.arange(self_filter.num_bins) - # LinLinEnergyFilters don't have bins so just add a None - elif isinstance(self_filter, openmc.LinLinEnergyFilter): + # 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 @@ -2120,14 +2120,14 @@ class Tally(object): filters = [type(filter1), type(filter2)] if isinstance(filter1, openmc.DistribcellFilter): filter1_bins = np.arange(filter1.num_bins) - elif isinstance(filter1, openmc.LinLinEnergyFilter): + 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 = np.arange(filter2.num_bins) - elif isinstance(filter2, openmc.LinLinEnergyFilter): + elif isinstance(filter2, openmc.EnergyFunctionFilter): filter2_bins = [None] else: filter2_bins = [filter2.get_bin(i) for i in range(filter2.num_bins)] @@ -2960,7 +2960,7 @@ class Tally(object): if isinstance(find_filter, openmc.DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) - elif isinstance(find_filter, openmc.LinLinEnergyFilter): + elif isinstance(find_filter, openmc.EnergyFunctionFilter): filter_bins = [None] else: num_bins = find_filter.num_bins @@ -3109,7 +3109,7 @@ class Tally(object): if isinstance(find_filter, openmc.DistribcellFilter): filter_bins = np.arange(find_filter.num_bins) - elif isinstance(find_filter, openmc.LinLinEnergyFilter): + elif isinstance(find_filter, openmc.EnergyFunctionFilter): filter_bins = [None] else: num_bins = find_filter.num_bins diff --git a/src/constants.F90 b/src/constants.F90 index bd9a0508b..2c27ce3ee 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -344,20 +344,20 @@ module constants ! Tally filter and map types 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_LINLINENERGY = 14 + 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 9279a48b8..dd8be29ea 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -3285,21 +3285,21 @@ contains ! Set the filter index in the tally find_filter array t % find_filter(FILTER_AZIMUTHAL) = j - case ('linlinenergy') + case ('energyfunction') ! Allocate and declare the filter type. - allocate(LinLinEnergyFilter::t % filters(j) % obj) + allocate(EnergyFunctionFilter::t % filters(j) % obj) select type (filt => t % filters(j) % obj) - type is (LinLinEnergyFilter) + type is (EnergyFunctionFilter) filt % n_bins = 1 ! Make sure this is continuous-energy mode. if (.not. run_CE) then - call fatal_error("LinLinEnergy filters are only supported for & + 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 LinLinEnergy & + 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") @@ -3308,7 +3308,7 @@ contains ! Allocate and store interpolant values. if (.not. check_for_node(node_filt, "y")) then - call fatal_error("y values not specified for LinLinEnergy & + 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") @@ -3316,7 +3316,7 @@ contains 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_LINLINENERGY) = j + t % find_filter(FILTER_ENERGYFUNCTION) = j case default ! Specified tally filter is invalid, raise error diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 064e9d646..f2680f4a9 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -192,17 +192,18 @@ module tally_filter end type AzimuthalFilter !=============================================================================== -! LinLinEnergyFilter +! EnergyFunctionFilter multiplies tally scores by an arbitrary function of +! incident energy described by a piecewise linear-linear interpolation. !=============================================================================== - type, extends(TallyFilter) :: LinLinEnergyFilter + type, extends(TallyFilter) :: EnergyFunctionFilter real(8), allocatable :: energy(:) real(8), allocatable :: y(:) contains - procedure :: get_next_bin => get_next_bin_linlinenergy - procedure :: to_statepoint => to_statepoint_linlinenergy - procedure :: text_label => text_label_linlinenergy - end type LinLinEnergyFilter + procedure :: get_next_bin => get_next_bin_energyfunction + procedure :: to_statepoint => to_statepoint_energyfunction + procedure :: text_label => text_label_energyfunction + end type EnergyFunctionFilter contains @@ -1223,22 +1224,22 @@ contains end function text_label_azimuthal !=============================================================================== -! LinLinEnergyFilter methods +! EnergyFunctionFilter methods !=============================================================================== - subroutine get_next_bin_linlinenergy(this, p, estimator, current_bin, & + subroutine get_next_bin_energyfunction(this, p, estimator, current_bin, & next_bin, weight) - class(LinLinEnergyFilter), 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 + 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 (LinLinEnergyFilter) + type is (EnergyFunctionFilter) if (current_bin == NO_BIN_FOUND) then n = size(this % energy) @@ -1272,32 +1273,33 @@ contains weight = ONE end if end select - end subroutine get_next_bin_linlinenergy + end subroutine get_next_bin_energyfunction - subroutine to_statepoint_linlinenergy(this, filter_group) - class(LinLinEnergyFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group + subroutine to_statepoint_energyfunction(this, filter_group) + class(EnergyFunctionFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group select type(this) - type is (LinLinEnergyFilter) - call write_dataset(filter_group, "type", "linlinenergy") + 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_linlinenergy + end subroutine to_statepoint_energyfunction - function text_label_linlinenergy(this, bin) result(label) - class(LinLinEnergyFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label + 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 (LinLinEnergyFilter) - label = "Lin-Lin Energy Interpolation [" & - // trim(to_str(this % energy(1))) // ", ..., " & - // trim(to_str(this % energy(size(this % energy)))) // "]" + 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_linlinenergy + end function text_label_energyfunction !=============================================================================== ! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and From 9c7d9432a2471b0fb493f0a9afb9012790c1a651 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 4 Dec 2016 20:05:18 -0500 Subject: [PATCH 05/10] Add EnergyFunctionFilter test --- openmc/filter.py | 17 +++++- tests/test_filter_energyfun/inputs_true.dat | 1 + tests/test_filter_energyfun/results_true.dat | 2 + .../test_filter_energyfun.py | 56 +++++++++++++++++++ 4 files changed, 74 insertions(+), 2 deletions(-) create mode 100644 tests/test_filter_energyfun/inputs_true.dat create mode 100644 tests/test_filter_energyfun/results_true.dat create mode 100644 tests/test_filter_energyfun/test_filter_energyfun.py diff --git a/openmc/filter.py b/openmc/filter.py index ca761db27..e16499660 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1631,11 +1631,24 @@ class EnergyFunctionFilter(Filter): """ - # Initialize Pandas DataFrame import pandas as pd df = pd.DataFrame() - filter_bins = np.repeat(hash(self), self.stride) + # 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(). + from hashlib import md5 + hash_fun = 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 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( diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat new file mode 100644 index 000000000..26f723392 --- /dev/null +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -0,0 +1 @@ +4dac56412e341cb64fa50d67e4dcf389b3b2cafb9bd8401eadc49654e050bad7988ca8239d5d9ee3e2953e1d36a2ca91a56557c7a34ccc72d520b41f462d9f1b \ 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 000000000..d3dc3ad6a --- /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 000000000..34141b038 --- /dev/null +++ b/tests/test_filter_energyfun/test_filter_energyfun.py @@ -0,0 +1,56 @@ +#!/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) + + # Make an EnergyFunctionFilter for the Am242m / Am242 branching ratio. + 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] + filt = openmc.EnergyFunctionFilter(x, y) + + # Make tallies. + tallies = [openmc.Tally() for i in range(2)] + for t in tallies: + t.scores = ['102'] + t.nuclides = ['Am241'] + tallies[1].filters = [filt] + 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.*', True) + harness.main() From f1eeb091a7584aa05377d3d7b353b86dfc5a1e7b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 4 Dec 2016 21:11:15 -0500 Subject: [PATCH 06/10] Add EnergyFunctionFilters to documentation --- docs/source/io_formats/statepoint.rst | 19 +++- docs/source/pythonapi/index.rst | 1 + docs/source/usersguide/input.rst | 27 ++++- src/relaxng/tallies.rnc | 30 +++-- src/relaxng/tallies.rng | 154 +++++++++++++++++--------- 5 files changed, 160 insertions(+), 71 deletions(-) diff --git a/docs/source/io_formats/statepoint.rst b/docs/source/io_formats/statepoint.rst index 5f34d4700..e334f2bf1 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 e364452b1..d66b36aa1 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 4a81615a1..ea407e725 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,22 @@ 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. (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/src/relaxng/tallies.rnc b/src/relaxng/tallies.rnc index 2143c8e47..57f5f3747 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 fc60b8373..fde06bac6 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 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 310fa134be0881f33aa6c91aa14584049dd62591 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 5 Dec 2016 16:51:55 -0500 Subject: [PATCH 07/10] Address #768 comments --- docs/source/usersguide/input.rst | 5 +- openmc/filter.py | 24 ++++--- openmc/tallies.py | 2 +- src/tally_filter.F90 | 62 ++++++++++++------- tests/test_filter_energyfun/inputs_true.dat | 2 +- .../test_filter_energyfun.py | 9 +-- 6 files changed, 60 insertions(+), 44 deletions(-) diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index ea407e725..09eab3262 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -1639,8 +1639,9 @@ The ```` element accepts the following sub-elements: :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. (Only used - for ``energyfunction`` filters) + (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 diff --git a/openmc/filter.py b/openmc/filter.py index ac6e961cf..6bd7e2623 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 @@ -158,8 +159,7 @@ class Filter(object): bins = np.array(bins) # If the bin is 0D numpy array, promote to 1D. - if bins.shape == (): - bins.shape = (1,) + bins = np.atleast_1d(bins) # Check the bin values. self.check_bins(bins) @@ -194,7 +194,7 @@ class Filter(object): pass - def to_xml(self): + def to_xml_element(self): """Return XML Element representing the Filter.""" element = ET.Element('filter') element.set('type', self.short_name.lower()) @@ -1616,7 +1616,8 @@ 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. + interpolation of energy and y values. Values outside of the given energy + range will be evaluated as zero. Parameters ---------- @@ -1723,8 +1724,7 @@ class EnergyFunctionFilter(Filter): energy = np.array(energy) # If the grid is 0D numpy array, promote to 1D. - if energy.shape == (): - energy.shape = (1,) + energy = np.atleast_1d(energy) # Make sure the values are Real and positive. cv.check_type('filter energy grid', energy, Iterable, Real) @@ -1739,8 +1739,7 @@ class EnergyFunctionFilter(Filter): y = np.array(y) # If the array is 0D, promote to 1D. - if y.shape == (): - y.shape = (1,) + y = np.atleast_1d(y) # Make sure the values are Real. cv.check_type('filter interpolant values', y, Iterable, Real) @@ -1751,7 +1750,7 @@ class EnergyFunctionFilter(Filter): def bins(self, bins): raise RuntimeError('EnergyFunctionFilters have no bins.') - def to_xml(self): + def to_xml_element(self): """Return XML Element representing the Filter.""" element = ET.Element('filter') element.set('type', self.short_name.lower()) @@ -1813,13 +1812,12 @@ class EnergyFunctionFilter(Filter): # 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(). - from hashlib import md5 - hash_fun = md5() + 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 of the - # digest are probably sufficient. + # 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) diff --git a/openmc/tallies.py b/openmc/tallies.py index 35fb32175..a9a1390df 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1036,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: diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index f2680f4a9..4b372defe 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 @@ -243,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) @@ -251,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 @@ -485,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) @@ -546,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) @@ -621,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) @@ -681,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 @@ -707,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) @@ -770,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) @@ -832,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) @@ -905,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. @@ -917,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) @@ -974,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) @@ -1026,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) @@ -1068,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) @@ -1129,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) @@ -1190,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) @@ -1253,7 +1269,7 @@ contains ! 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 = ONE + weight = ERROR_REAL else ! Search to find incoming energy bin. @@ -1270,7 +1286,7 @@ contains else next_bin = NO_BIN_FOUND - weight = ONE + weight = ERROR_REAL end if end select end subroutine get_next_bin_energyfunction diff --git a/tests/test_filter_energyfun/inputs_true.dat b/tests/test_filter_energyfun/inputs_true.dat index 26f723392..c9961071a 100644 --- a/tests/test_filter_energyfun/inputs_true.dat +++ b/tests/test_filter_energyfun/inputs_true.dat @@ -1 +1 @@ -4dac56412e341cb64fa50d67e4dcf389b3b2cafb9bd8401eadc49654e050bad7988ca8239d5d9ee3e2953e1d36a2ca91a56557c7a34ccc72d520b41f462d9f1b \ No newline at end of file +2a08e03d9e0f40486fa26ba00800d264f43906a38c851da48f0fd76bacc0b862e8f45f16c970b7c62db03d96e3acb1b58e74179f554a851ccaaf4492e97fc163 \ No newline at end of file diff --git a/tests/test_filter_energyfun/test_filter_energyfun.py b/tests/test_filter_energyfun/test_filter_energyfun.py index 34141b038..36dffb871 100644 --- a/tests/test_filter_energyfun/test_filter_energyfun.py +++ b/tests/test_filter_energyfun/test_filter_energyfun.py @@ -18,15 +18,16 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Add Am241 to the fuel. self._input_set.materials[1].add_nuclide('Am241', 1e-7) - # Make an EnergyFunctionFilter for the Am242m / Am242 branching ratio. + # Make an EnergyFunctionFilter for the 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] filt = openmc.EnergyFunctionFilter(x, y) # Make tallies. - tallies = [openmc.Tally() for i in range(2)] + tallies = [openmc.Tally(), openmc.Tally()] for t in tallies: - t.scores = ['102'] + t.scores = ['(n,gamma)'] t.nuclides = ['Am241'] tallies[1].filters = [filt] self._input_set.tallies = openmc.Tallies(tallies) @@ -52,5 +53,5 @@ class FilterEnergyFunHarness(PyAPITestHarness): if __name__ == '__main__': - harness = FilterEnergyFunHarness('statepoint.10.*', True) + harness = FilterEnergyFunHarness('statepoint.10.h5', True) harness.main() From 4dee578b6cd7f74b6e3dba20080362fd6e7cfe6b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 5 Dec 2016 17:56:32 -0500 Subject: [PATCH 08/10] Manually inherit Filter subclass method docstrings --- openmc/filter.py | 28 ++++++++++++++++++++++++++++ 1 file changed, 28 insertions(+) diff --git a/openmc/filter.py b/openmc/filter.py index 6bd7e2623..6eba9ec0f 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -26,9 +26,37 @@ _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): + 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) From 7b7840a5fcd3990323b86897c176e03d0b3b29a4 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 7 Dec 2016 15:38:18 -0500 Subject: [PATCH 09/10] EnergyFunctionFilter.from_tabulated1d constructor --- openmc/filter.py | 25 +++++++++++++++++++ .../test_filter_energyfun.py | 15 ++++++++--- 2 files changed, 36 insertions(+), 4 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 6eba9ec0f..100b0f8c9 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -1730,6 +1730,31 @@ class EnergyFunctionFilter(Filter): 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 diff --git a/tests/test_filter_energyfun/test_filter_energyfun.py b/tests/test_filter_energyfun/test_filter_energyfun.py index 36dffb871..b7fa68fa4 100644 --- a/tests/test_filter_energyfun/test_filter_energyfun.py +++ b/tests/test_filter_energyfun/test_filter_energyfun.py @@ -18,18 +18,25 @@ class FilterEnergyFunHarness(PyAPITestHarness): # Add Am241 to the fuel. self._input_set.materials[1].add_nuclide('Am241', 1e-7) - # Make an EnergyFunctionFilter for the Am242m / Am242 branching ratio - # (from ENDF/B-VII.1 data). + # 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] - filt = openmc.EnergyFunctionFilter(x, y) + + # 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 = [filt] + tallies[1].filters = [filt1] self._input_set.tallies = openmc.Tallies(tallies) # Export inputs to xml. From 2bf43e166a8f0f67303ce803eef961b36914bb2a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 9 Dec 2016 14:03:11 -0500 Subject: [PATCH 10/10] Address #768 comments --- openmc/filter.py | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 100b0f8c9..dd3e21e5c 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -45,7 +45,8 @@ class FilterMeta(ABCMeta): 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): + 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: @@ -183,10 +184,7 @@ 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. + # Format the bins as a 1D numpy array. bins = np.atleast_1d(bins) # Check the bin values. @@ -223,7 +221,14 @@ class Filter(object): pass def to_xml_element(self): - """Return XML Element representing the Filter.""" + """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)) @@ -1773,10 +1778,7 @@ class EnergyFunctionFilter(Filter): @energy.setter def energy(self, energy): - # Make sure the energy grid is a numpy array. - energy = np.array(energy) - - # If the grid is 0D numpy array, promote to 1D. + # Format the bins as a 1D numpy array. energy = np.atleast_1d(energy) # Make sure the values are Real and positive. @@ -1788,10 +1790,7 @@ class EnergyFunctionFilter(Filter): @y.setter def y(self, y): - # Make sure the values are in a numpy array. - y = np.array(y) - - # If the array is 0D, promote to 1D. + # Format the bins as a 1D numpy array. y = np.atleast_1d(y) # Make sure the values are Real. @@ -1804,7 +1803,6 @@ class EnergyFunctionFilter(Filter): raise RuntimeError('EnergyFunctionFilters have no bins.') def to_xml_element(self): - """Return XML Element representing the Filter.""" element = ET.Element('filter') element.set('type', self.short_name.lower()) element.set('energy', ' '.join(str(e) for e in self.energy))