diff --git a/openmc/arithmetic.py b/openmc/arithmetic.py index 8574c4873e..5869cad802 100644 --- a/openmc/arithmetic.py +++ b/openmc/arithmetic.py @@ -1,5 +1,7 @@ import sys +import copy from numbers import Integral +from collections import Iterable import numpy as np @@ -14,7 +16,7 @@ if sys.version_info[0] >= 3: _TALLY_ARITHMETIC_OPS = ['+', '-', '*', '/', '^'] # Acceptable tally aggregation operations -_TALLY_AGGREGATE_OPS = ['sum', 'mean'] +_TALLY_AGGREGATE_OPS = ['sum', 'avg'] class CrossScore(object): @@ -186,6 +188,23 @@ class CrossNuclide(object): return existing def __repr__(self): + return self.name + + + @property + def left_nuclide(self): + return self._left_nuclide + + @property + def right_nuclide(self): + return self._right_nuclide + + @property + def binary_op(self): + return self._binary_op + + @property + def name(self): string = '' @@ -207,18 +226,6 @@ class CrossNuclide(object): return string - @property - def left_nuclide(self): - return self._left_nuclide - - @property - def right_nuclide(self): - return self._right_nuclide - - @property - def binary_op(self): - return self._binary_op - @left_nuclide.setter def left_nuclide(self, left_nuclide): cv.check_type('left_nuclide', left_nuclide, @@ -430,7 +437,7 @@ class CrossFilter(object): filter_index = left_index * self.right_filter.num_bins + right_index return filter_index - def get_pandas_dataframe(self, datasize, summary=None): + def get_pandas_dataframe(self, data_size, summary=None): """Builds a Pandas DataFrame for the CrossFilter's bins. This method constructs a Pandas DataFrame object for the CrossFilter @@ -445,7 +452,7 @@ class CrossFilter(object): Parameters ---------- - datasize : Integral + data_size : Integral The total number of bins in the tally corresponding to this filter summary : None or Summary An optional Summary object to be used to construct columns for @@ -472,19 +479,18 @@ class CrossFilter(object): # If left and right filters are identical, do not combine bins if self.left_filter == self.right_filter: - df = self.left_filter.get_pandas_dataframe(datasize, summary) + df = self.left_filter.get_pandas_dataframe(data_size, summary) # If left and right filters are different, combine their bins else: - left_df = self.left_filter.get_pandas_dataframe(datasize, summary) - right_df = self.right_filter.get_pandas_dataframe(datasize, summary) + left_df = self.left_filter.get_pandas_dataframe(data_size, summary) + right_df = self.right_filter.get_pandas_dataframe(data_size, summary) left_df = left_df.astype(str) right_df = right_df.astype(str) df = '(' + left_df + ' ' + self.binary_op + ' ' + right_df + ')' return df - class AggregateScore(object): """A special-purpose tally score used to encapsulate an aggregate of a subset or all of tally's scores for tally aggregation. @@ -494,7 +500,7 @@ class AggregateScore(object): scores : Iterable of str or CrossScore The scores included in the aggregation aggregate_op : str - The tally aggregation operator (e.g., 'sum', 'mean', etc.) used + The tally aggregation operator (e.g., 'sum', 'avg', etc.) used to aggregate across a tally's scores with this AggregateScore Attributes @@ -502,7 +508,7 @@ class AggregateScore(object): scores : Iterable of str or CrossScore The scores included in the aggregation aggregate_op : str - The tally aggregation operator (e.g., 'sum', 'mean', etc.) used + The tally aggregation operator (e.g., 'sum', 'avg', etc.) used to aggregate across a tally's scores with this AggregateScore """ @@ -556,10 +562,16 @@ class AggregateScore(object): def aggregate_op(self): return self._aggregate_op + @property + def name(self): + + # Append each score in the aggregate to the string + string = '(' + ', '.join(self.scores) + ')' + return string + @scores.setter def scores(self, scores): - cv.check_iterable_type('scores', scores, - (basestring, CrossScore, AggregateScore)) + cv.check_iterable_type('scores', scores, basestring) self._scores = scores @aggregate_op.setter @@ -578,7 +590,7 @@ class AggregateNuclide(object): nuclides : Iterable of str or Nuclide or CrossNuclide The nuclides included in the aggregation aggregate_op : str - The tally aggregation operator (e.g., 'sum', 'mean', etc.) used + The tally aggregation operator (e.g., 'sum', 'avg', etc.) used to aggregate across a tally's nuclides with this AggregateNuclide Attributes @@ -586,7 +598,7 @@ class AggregateNuclide(object): nuclides : Iterable of str or Nuclide or CrossNuclide The nuclides included in the aggregation aggregate_op : str - The tally aggregation operator (e.g., 'sum', 'mean', etc.) used + The tally aggregation operator (e.g., 'sum', 'avg', etc.) used to aggregate across a tally's nuclides with this AggregateNuclide """ @@ -644,10 +656,19 @@ class AggregateNuclide(object): def aggregate_op(self): return self._aggregate_op + @property + def name(self): + + # Append each nuclide in the aggregate to the string + names = [nuclide.name if isinstance(nuclide, Nuclide) else str(nuclide) + for nuclide in self.nuclides] + string = '(' + ', '.join(map(str, names)) + ')' + return string + @nuclides.setter def nuclides(self, nuclides): cv.check_iterable_type('nuclides', nuclides, - (basestring, Nuclide, CrossNuclide, AggregateNuclide)) + (basestring, Nuclide, CrossNuclide)) self._nuclides = nuclides @aggregate_op.setter @@ -668,7 +689,7 @@ class AggregateFilter(object): bins : Iterable of tuple The filter bins included in the aggregation aggregate_op : str - The tally aggregation operator (e.g., 'sum', 'mean', etc.) used + The tally aggregation operator (e.g., 'sum', 'avg', etc.) used to aggregate across a tally filter's bins with this AggregateFilter Attributes @@ -678,7 +699,7 @@ class AggregateFilter(object): aggregate_filter : filter The filter included in the aggregation aggregate_op : str - The tally aggregation operator (e.g., 'sum', 'mean', etc.) used + The tally aggregation operator (e.g., 'sum', 'avg', etc.) used to aggregate across a tally filter's bins with this AggregateFilter bins : Iterable of tuple The filter bins included in the aggregation @@ -715,6 +736,21 @@ class AggregateFilter(object): def __ne__(self, other): return not self == other + def __gt__(self, other): + if self.type != other.type: + if self.aggregate_filter.type in _FILTER_TYPES and \ + other.aggregate_filter.type in _FILTER_TYPES: + delta = _FILTER_TYPES.index(self.aggregate_filter.type) - \ + _FILTER_TYPES.index(other.aggregate_filter.type) + return delta > 0 + else: + return False + else: + return False + + def __lt__(self, other): + return not self > other + def __repr__(self): string = 'AggregateFilter\n' string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type) @@ -759,7 +795,7 @@ class AggregateFilter(object): @property def num_bins(self): - return 1 if self.aggregate_filter else 0 + return len(self.bins) if self.aggregate_filter else 0 @property def stride(self): @@ -776,14 +812,13 @@ class AggregateFilter(object): @aggregate_filter.setter def aggregate_filter(self, aggregate_filter): - cv.check_type('aggregate_filter', aggregate_filter, - (Filter, CrossFilter, AggregateFilter)) + cv.check_type('aggregate_filter', aggregate_filter, (Filter, CrossFilter)) self._aggregate_filter = aggregate_filter @bins.setter def bins(self, bins): - cv.check_iterable_type('bins', bins, (Integral, tuple)) - self._bins = bins + cv.check_iterable_type('bins', bins, Iterable) + self._bins = list(map(tuple, bins)) @aggregate_op.setter def aggregate_op(self, aggregate_op): @@ -823,15 +858,14 @@ class AggregateFilter(object): """ - if filter_bin not in self.bins and \ - filter_bin != self._aggregate_filter.bins: + if filter_bin not in self.bins: msg = 'Unable to get the bin index for AggregateFilter since ' \ '"{0}" is not one of the bins'.format(filter_bin) raise ValueError(msg) else: - return 0 + return self.bins.index(filter_bin) - def get_pandas_dataframe(self, datasize, summary=None): + def get_pandas_dataframe(self, data_size, summary=None): """Builds a Pandas DataFrame for the AggregateFilter's bins. This method constructs a Pandas DataFrame object for the AggregateFilter @@ -840,7 +874,7 @@ class AggregateFilter(object): Parameters ---------- - datasize : Integral + data_size : Integral The total number of bins in the tally corresponding to this filter summary : None or Summary An optional Summary object to be used to construct columns for @@ -868,14 +902,80 @@ class AggregateFilter(object): import pandas as pd - # Construct a sring representing the filter aggregation - aggregate_bin = '{0}('.format(self.aggregate_op) - aggregate_bin += ', '.join(map(str, self.bins)) + ')' + # Create NumPy array of the bin tuples for repeating / tiling + filter_bins = np.empty(self.num_bins, dtype=tuple) + for i, bin in enumerate(self.bins): + filter_bins[i] = bin - # Construct NumPy array of bin repeated for each element in dataframe - aggregate_bin_array = np.array([aggregate_bin]) - aggregate_bin_array = np.repeat(aggregate_bin_array, datasize) + # Repeat and tile bins as needed for DataFrame + filter_bins = np.repeat(filter_bins, self.stride) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) - # Construct Pandas DataFrame for the AggregateFilter - df = pd.DataFrame({self.type: aggregate_bin_array}) + # Create DataFrame with aggregated bins + df = pd.DataFrame({self.type: filter_bins}) return df + + def can_merge(self, other): + """Determine if AggregateFilter can be merged with another. + + Parameters + ---------- + other : AggregateFilter + Filter to compare with + + Returns + ------- + bool + Whether the filter can be merged + + """ + + if not isinstance(other, AggregateFilter): + return False + + # Filters must be of the same type + elif self.type != other.type: + return False + + # None of the bins in this filter should match in the other filter + for bin in self.bins: + if bin in other.bins: + return False + + # If all conditional checks passed then filters are mergeable + return True + + def merge(self, other): + """Merge this aggregatefilter with another. + + Parameters + ---------- + other : AggregateFilter + Filter to merge with + + Returns + ------- + merged_filter : AggregateFilter + Filter resulting from the merge + + """ + + if not self.can_merge(other): + msg = 'Unable to merge "{0}" with "{1}" ' \ + 'filters'.format(self.type, other.type) + raise ValueError(msg) + + # Create deep copy of filter to return as merged filter + merged_filter = copy.deepcopy(self) + + # Merge unique filter bins + merged_bins = self.bins + other.bins + + # Sort energy bin edges + if 'energy' in self.type: + merged_bins = sorted(merged_bins) + + # Assign merged bins to merged filter + merged_filter.bins = list(merged_bins) + return merged_filter diff --git a/openmc/element.py b/openmc/element.py index 9f04abfdab..dda110ea7a 100644 --- a/openmc/element.py +++ b/openmc/element.py @@ -57,6 +57,15 @@ class Element(object): def __ne__(self, other): return not self == other + def __gt__(self, other): + return repr(self) > repr(other) + + def __lt__(self, other): + return not self > other + + def __hash__(self): + return hash(repr(self)) + def __hash__(self): return hash(repr(self)) diff --git a/openmc/filter.py b/openmc/filter.py index 783f39043e..2ae8eeb626 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -77,6 +77,25 @@ class Filter(object): def __ne__(self, other): return not self == other + def __gt__(self, other): + if self.type != other.type: + if self.type in _FILTER_TYPES and other.type in _FILTER_TYPES: + delta = _FILTER_TYPES.index(self.type) - \ + _FILTER_TYPES.index(other.type) + return delta > 0 + else: + return False + else: + # Compare largest/smallest energy bin edges in energy filters + # This logic is used when merging tallies with energy filters + if 'energy' in self.type and 'energy' in other.type: + return self.bins[0] >= other.bins[-1] + else: + return max(self.bins) > max(other.bins) + + def __lt__(self, other): + return not self > other + def __hash__(self): return hash(repr(self)) @@ -246,20 +265,28 @@ class Filter(object): return False # Filters must be of the same type - elif self.type != other.type: + if self.type != other.type: return False # Distribcell filters cannot have more than one bin - elif self.type == 'distribcell': + if self.type == 'distribcell': return False # Mesh filters cannot have more than one bin elif self.type == 'mesh': return False - # Different energy bins are not mergeable + # Different energy bins structures must be mutually exclusive and + # share only one shared bin edge at the minimum or maximum energy elif 'energy' in self.type: - return False + # This low energy edge coincides with other's high energy edge + if self.bins[0] == other.bins[-1]: + return True + # This high energy edge coincides with other's low energy edge + elif self.bins[-1] == other.bins[0]: + return True + else: + return False else: return True @@ -288,9 +315,21 @@ class Filter(object): merged_filter = copy.deepcopy(self) # Merge unique filter bins - merged_bins = list(set(np.concatenate((self.bins, other.bins)))) - merged_filter.bins = merged_bins - merged_filter.num_bins = len(merged_bins) + merged_bins = np.concatenate((self.bins, other.bins)) + merged_bins = np.unique(merged_bins) + + # Sort energy bin edges + if 'energy' in self.type: + merged_bins = sorted(merged_bins) + + # Assign merged bins to merged filter + merged_filter.bins = list(merged_bins) + + # Count bins in the merged filter + if 'energy' in merged_filter.type: + merged_filter.num_bins = len(merged_bins) - 1 + else: + merged_filter.num_bins = len(merged_bins) return merged_filter @@ -521,14 +560,8 @@ class Filter(object): """ - # Attempt to import Pandas - try: - import pandas as pd - except ImportError: - msg = 'The Pandas Python package must be installed on your system' - raise ImportError(msg) - # Initialize Pandas DataFrame + import pandas as pd df = pd.DataFrame() # mesh filters @@ -707,7 +740,6 @@ class Filter(object): filter_bins = np.repeat(filter_bins, self.stride) tile_factor = data_size / len(filter_bins) filter_bins = np.tile(filter_bins, tile_factor) - filter_bins = filter_bins df = pd.DataFrame({self.type : filter_bins}) # If OpenCG level info DataFrame was created, concatenate diff --git a/openmc/geometry.py b/openmc/geometry.py index 9788671f5c..dac0bd90f1 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -65,8 +65,10 @@ class Geometry(object): # Find the distribcell index of the cell. cells = self.get_all_cells() - if path[-1] in cells: - distribcell_index = cells[path[-1]].distribcell_index + for cell in cells: + if cell.id == path[-1]: + distribcell_index = cell.distribcell_index + break else: raise RuntimeError('Could not find cell {} specified in a \ distribcell filter'.format(path[-1])) @@ -94,7 +96,16 @@ class Geometry(object): """ - return self._root_universe.get_all_cells() + all_cells = self._root_universe.get_all_cells() + cells = set() + + for cell in all_cells.values(): + if cell._type == 'normal': + cells.add(cell) + + cells = list(cells) + cells.sort(key=lambda x: x.id) + return cells def get_all_universes(self): """Return all universes defined @@ -106,7 +117,15 @@ class Geometry(object): """ - return self._root_universe.get_all_universes() + all_universes = self._root_universe.get_all_universes() + universes = set() + + for universe in all_universes.values(): + universes.add(universe) + + universes = list(universes) + universes.sort(key=lambda x: x.id) + return universes def get_all_nuclides(self): """Return all nuclides assigned to a material in the geometry @@ -150,10 +169,19 @@ class Geometry(object): return materials def get_all_material_cells(self): + """Return all cells filled by a material + + Returns + ------- + list of openmc.universe.Cell + Cells filled by Materials in the geometry + + """ + all_cells = self.get_all_cells() material_cells = set() - for cell_id, cell in all_cells.items(): + for cell in all_cells: if cell._type == 'normal': material_cells.add(cell) @@ -174,9 +202,9 @@ class Geometry(object): all_universes = self.get_all_universes() material_universes = set() - for universe_id, universe in all_universes.items(): - cells = universe._cells - for cell_id, cell in cells.items(): + for universe in all_universes: + cells = universe.cells + for cell in cells: if cell._type == 'normal': material_universes.add(universe) @@ -184,6 +212,227 @@ class Geometry(object): material_universes.sort(key=lambda x: x.id) return material_universes + def get_all_lattices(self): + """Return all lattices defined + + Returns + ------- + list of openmc.universe.Lattice + Lattices in the geometry + + """ + + cells = self.get_all_cells() + lattices = set() + + for cell in cells: + if isinstance(cell.fill, openmc.Lattice): + lattices.add(cell.fill) + + lattices = list(lattices) + lattices.sort(key=lambda x: x.id) + return lattices + + def get_materials_by_name(self, name, case_sensitive=False, matching=False): + """Return a list of materials with matching names. + + Parameters + ---------- + name : str + The name to match + case_sensitive : bool + Whether to distinguish upper and lower case letters in each + material's name (default is True) + matching : bool + Whether the names must match completely (default is True) + + Returns + ------- + list of openmc.material.Material + Materials matching the queried name + + """ + + if not case_sensitive: + name = name.lower() + + all_materials = self.get_all_materials() + materials = set() + + for material in all_materials: + material_name = material.name + if not case_sensitive: + material_name = material_name.lower() + + if material_name == name: + materials.add(material) + elif not matching and name in material_name: + materials.add(material) + + materials = list(materials) + materials.sort(key=lambda x: x.id) + return materials + + def get_cells_by_name(self, name, case_sensitive=False, matching=False): + """Return a list of cells with matching names. + + Parameters + ---------- + name : str + The name to search match + case_sensitive : bool + Whether to distinguish upper and lower case letters in each + cell's name (default is True) + matching : bool + Whether the names must match completely (default is True) + + Returns + ------- + list of openmc.universe.Cell + Cells matching the queried name + + """ + + if not case_sensitive: + name = name.lower() + + all_cells = self.get_all_cells() + cells = set() + + for cell in all_cells: + cell_name = cell.name + if not case_sensitive: + cell_name = cell_name.lower() + + if cell_name == name: + cells.add(cell) + elif not matching and name in cell_name: + cells.add(cell) + + cells = list(cells) + cells.sort(key=lambda x: x.id) + return cells + + def get_cells_by_fill_name(self, name, case_sensitive=False, matching=False): + """Return a list of cells with fills with matching names. + + Parameters + ---------- + name : str + The name to match + case_sensitive : bool + Whether to distinguish upper and lower case letters in each + cell's name (default is True) + matching : bool + Whether the names must match completely (default is True) + + Returns + ------- + list of openmc.universe.Cell + Cells with fills matching the queried name + + """ + + if not case_sensitive: + name = name.lower() + + all_cells = self.get_all_cells() + cells = set() + + for cell in all_cells: + cell_fill_name = cell.fill.name + if not case_sensitive: + cell_fill_name = cell_fill_name.lower() + + if cell_fill_name == name: + cells.add(cell) + elif not matching and name in cell_fill_name: + cells.add(cell) + + cells = list(cells) + cells.sort(key=lambda x: x.id) + return cells + + def get_universes_by_name(self, name, case_sensitive=False, matching=False): + """Return a list of universes with matching names. + + Parameters + ---------- + name : str + The name to match + case_sensitive : bool + Whether to distinguish upper and lower case letters in each + universe's name (default is True) + matching : bool + Whether the names must match completely (default is True) + + Returns + ------- + list of openmc.universe.Universe + Universes matching the queried name + + """ + + if not case_sensitive: + name = name.lower() + + all_universes = self.get_all_universes() + universes = set() + + for universe in all_universes: + universe_name = universe.name + if not case_sensitive: + universe_name = universe_name.lower() + + if universe_name == name: + universes.add(universe) + elif not matching and name in universe_name: + universes.add(universe) + + universes = list(universes) + universes.sort(key=lambda x: x.id) + return universes + + def get_lattices_by_name(self, name, case_sensitive=False, matching=False): + """Return a list of lattices with matching names. + + Parameters + ---------- + name : str + The name to match + case_sensitive : bool + Whether to distinguish upper and lower case letters in each + lattice's name (default is True) + matching : bool + Whether the names must match completely (default is True) + + Returns + ------- + list of openmc.universe.Lattice + Lattices matching the queried name + + """ + + if not case_sensitive: + name = name.lower() + + all_lattices = self.get_all_lattices() + lattices = set() + + for lattice in all_lattices: + lattice_name = lattice.name + if not case_sensitive: + lattice_name = lattice_name.lower() + + if lattice_name == name: + lattices.add(lattice) + elif not matching and name in lattice_name: + lattices.add(lattice) + + lattices = list(lattices) + lattices.sort(key=lambda x: x.id) + return lattices + class GeometryFile(object): """Geometry file used for an OpenMC simulation. Corresponds directly to the diff --git a/openmc/mgxs/groups.py b/openmc/mgxs/groups.py index 253df520ee..a1e03c3371 100644 --- a/openmc/mgxs/groups.py +++ b/openmc/mgxs/groups.py @@ -54,10 +54,12 @@ class EnergyGroups(object): def __eq__(self, other): if not isinstance(other, EnergyGroups): return False - elif (self.group_edges != other.group_edges).all(): + elif self.num_groups != other.num_groups: return False - else: + elif np.allclose(self.group_edges, other.group_edges): return True + else: + return False def __ne__(self, other): return not self == other @@ -236,3 +238,64 @@ class EnergyGroups(object): condensed_groups.group_edges = group_edges return condensed_groups + + def can_merge(self, other): + """Determine if energy groups can be merged with another. + + Parameters + ---------- + other : EnergyGroups + EnergyGroups to compare with + + Returns + ------- + bool + Whether the energy groups can be merged + + """ + + if not isinstance(other, EnergyGroups): + return False + + # If the energy group structures match then groups are mergeable + if self == other: + return True + + # This low energy edge coincides with other's high energy edge + if self.group_edges[0] == other.group_edges[-1]: + return True + # This high energy edge coincides with other's low energy edge + elif self.group_edges[-1] == other.group_edges[0]: + return True + else: + return False + + def merge(self, other): + """Merge this energy groups with another. + + Parameters + ---------- + other : EnergyGroups + EnergyGroups to merge with + + Returns + ------- + merged_groups : EnergyGroups + EnergyGroups resulting from the merge + + """ + + if not self.can_merge(other): + raise ValueError('Unable to merge energy groups') + + # Create deep copy to return as merged energy groups + merged_groups = copy.deepcopy(self) + + # Merge unique filter bins + merged_edges = np.concatenate((self.group_edges, other.group_edges)) + merged_edges = np.unique(merged_edges) + merged_edges = sorted(merged_edges) + + # Assign merged edges to merged groups + merged_groups.group_edges = list(merged_edges) + return merged_groups diff --git a/openmc/mgxs/library.py b/openmc/mgxs/library.py index e87b4bdaef..a38e42d245 100644 --- a/openmc/mgxs/library.py +++ b/openmc/mgxs/library.py @@ -116,11 +116,11 @@ class Library(object): clone._by_nuclide = self.by_nuclide clone._mgxs_types = self.mgxs_types clone._domain_type = self.domain_type - clone._domains = self.domains + clone._domains = copy.deepcopy(self.domains) clone._correction = self.correction clone._energy_groups = copy.deepcopy(self.energy_groups, memo) clone._tally_trigger = copy.deepcopy(self.tally_trigger, memo) - clone._all_mgxs = self.all_mgxs + clone._all_mgxs = copy.deepcopy(self.all_mgxs) clone._sp_filename = self._sp_filename clone._keff = self._keff clone._sparse = self.sparse diff --git a/openmc/mgxs/mgxs.py b/openmc/mgxs/mgxs.py index 09e0f95907..6be8255f4a 100644 --- a/openmc/mgxs/mgxs.py +++ b/openmc/mgxs/mgxs.py @@ -67,6 +67,10 @@ class MGXS(object): The energy group structure for energy condensation by_nuclide : bool If true, computes cross sections for each nuclide in domain + nuclides : Iterable of basestring + The user-specified nuclides to compute cross sections. If by_nuclide + is True but nuclides are not specified by the user, all nuclides in the + spatial domain will be used. name : str, optional Name of the multi-group cross section. Used as a label to identify tallies in OpenMC 'tallies.xml' file. @@ -111,6 +115,8 @@ class MGXS(object): sparse : bool Whether or not the MGXS' tallies use SciPy's LIL sparse matrix format for compressed data storage + derived : bool + Whether or not the MGXS is merged from one or more other MGXS """ @@ -123,6 +129,7 @@ class MGXS(object): self._name = '' self._rxn_type = None self._by_nuclide = None + self._nuclides = None self._domain = None self._domain_type = None self._energy_groups = None @@ -131,6 +138,7 @@ class MGXS(object): self._rxn_rate_tally = None self._xs_tally = None self._sparse = False + self._derived = False self.name = name self.by_nuclide = by_nuclide @@ -151,6 +159,7 @@ class MGXS(object): clone._name = self.name clone._rxn_type = self.rxn_type clone._by_nuclide = self.by_nuclide + clone._nuclides = copy.deepcopy(self._nuclides) clone._domain = self.domain clone._domain_type = self.domain_type clone._energy_groups = copy.deepcopy(self.energy_groups, memo) @@ -158,6 +167,7 @@ class MGXS(object): clone._rxn_rate_tally = copy.deepcopy(self._rxn_rate_tally, memo) clone._xs_tally = copy.deepcopy(self._xs_tally, memo) clone._sparse = self.sparse + clone._derived = self.derived clone._tallies = OrderedDict() for tally_type, tally in self.tallies.items(): @@ -232,8 +242,7 @@ class MGXS(object): @property def num_subdomains(self): - tally = list(self.tallies.values())[0] - domain_filter = tally.find_filter(self.domain_type) + domain_filter = self.xs_tally.find_filter(self.domain_type) return domain_filter.num_bins @property @@ -250,6 +259,10 @@ class MGXS(object): else: return 'sum' + @property + def derived(self): + return self._derived + @name.setter def name(self, name): cv.check_type('name', name, basestring) @@ -260,6 +273,11 @@ class MGXS(object): cv.check_type('by_nuclide', by_nuclide, bool) self._by_nuclide = by_nuclide + @nuclides.setter + def nuclides(self, nuclides): + cv.check_iterable_type('nuclides', nuclides, basestring) + self._nuclides = nuclides + @domain.setter def domain(self, domain): cv.check_type('domain', domain, tuple(_DOMAINS)) @@ -389,8 +407,14 @@ class MGXS(object): if self.domain is None: raise ValueError('Unable to get all nuclides without a domain') - nuclides = self.domain.get_all_nuclides() - return nuclides.keys() + # If the user defined nuclides, return them + if self._nuclides: + return self._nuclides + + # Otherwise, return all nuclides in the spatial domain + else: + nuclides = self.domain.get_all_nuclides() + return nuclides.keys() def get_nuclide_density(self, nuclide): """Get the atomic number density in units of atoms/b-cm for a nuclide @@ -553,7 +577,7 @@ class MGXS(object): # If computing xs for each nuclide, replace CrossNuclides with originals if self.by_nuclide: self.xs_tally._nuclides = [] - nuclides = self.domain.get_all_nuclides() + nuclides = self.get_all_nuclides() for nuclide in nuclides: self.xs_tally.add_nuclide(openmc.Nuclide(nuclide)) @@ -682,7 +706,7 @@ class MGXS(object): # Construct a collection of the domain filter bins if not isinstance(subdomains, basestring): - cv.check_iterable_type('subdomains', subdomains, Integral) + cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=2) for subdomain in subdomains: filters.append(self.domain_type) filter_bins.append((subdomain,)) @@ -855,19 +879,181 @@ class MGXS(object): # Clone this MGXS to initialize the subdomain-averaged version avg_xs = copy.deepcopy(self) - avg_xs._rxn_rate_tally = None - avg_xs._xs_tally = None - # Average each of the tallies across subdomains - for tally_type, tally in avg_xs.tallies.items(): - tally_avg = tally.summation(filter_type=self.domain_type, - filter_bins=subdomains) - avg_xs.tallies[tally_type] = tally_avg + if self.derived: + avg_xs._rxn_rate_tally = avg_xs.rxn_rate_tally.average( + filter_type=self.domain_type, filter_bins=subdomains) + else: + avg_xs._rxn_rate_tally = None + avg_xs._xs_tally = None - avg_xs._domain_type = 'sum({0})'.format(self.domain_type) + # Average each of the tallies across subdomains + for tally_type, tally in avg_xs.tallies.items(): + tally_avg = tally.average(filter_type=self.domain_type, + filter_bins=subdomains) + avg_xs.tallies[tally_type] = tally_avg + + avg_xs._domain_type = 'avg({0})'.format(self.domain_type) avg_xs.sparse = self.sparse return avg_xs + def get_slice(self, nuclides=[], groups=[]): + """Build a sliced MGXS for the specified nuclides and energy groups. + + This method constructs a new MGXS to encapsulate a subset of the data + represented by this MGXS. The subset of data to include in the tally + slice is determined by the nuclides and energy groups specified in + the input parameters. + + Parameters + ---------- + nuclides : list of str + A list of nuclide name strings + (e.g., ['U-235', 'U-238']; default is []) + groups : list of Integral + A list of energy group indices starting at 1 for the high energies + (e.g., [1, 2, 3]; default is []) + + Returns + ------- + MGXS + A new tally which encapsulates the subset of data requested for the + nuclide(s) and/or energy group(s) requested in the parameters. + + """ + + cv.check_iterable_type('nuclides', nuclides, basestring) + cv.check_iterable_type('energy_groups', groups, Integral) + + # Build lists of filters and filter bins to slice + if len(groups) == 0: + filters = [] + filter_bins = [] + else: + filter_bins = [] + for group in groups: + group_bounds = self.energy_groups.get_group_bounds(group) + filter_bins.append(group_bounds) + filter_bins = [tuple(filter_bins)] + filters = ['energy'] + + # Clone this MGXS to initialize the sliced version + slice_xs = copy.deepcopy(self) + slice_xs._rxn_rate_tally = None + slice_xs._xs_tally = None + + # Slice each of the tallies across nuclides and energy groups + for tally_type, tally in slice_xs.tallies.items(): + slice_nuclides = [nuc for nuc in nuclides if nuc in tally.nuclides] + if len(groups) != 0 and tally.contains_filter('energy'): + tally_slice = tally.get_slice(filters=filters, + filter_bins=filter_bins, nuclides=slice_nuclides) + else: + tally_slice = tally.get_slice(nuclides=slice_nuclides) + slice_xs.tallies[tally_type] = tally_slice + + # Assign sliced energy group structure to sliced MGXS + if groups: + new_group_edges = [] + for group in groups: + group_edges = self.energy_groups.get_group_bounds(group) + new_group_edges.extend(group_edges) + new_group_edges = np.unique(new_group_edges) + slice_xs.energy_groups.group_edges = sorted(new_group_edges) + + # Assign sliced nuclides to sliced MGXS + if nuclides: + slice_xs.nuclides = nuclides + + slice_xs.sparse = self.sparse + return slice_xs + + def can_merge(self, other): + """Determine if another MGXS can be merged with this one + + If results have been loaded from a statepoint, then MGXS are only + mergeable along one and only one of enegy groups or nuclides. + + Parameters + ---------- + other : MGXS + MGXS to check for merging + + """ + + if not isinstance(other, type(self)): + return False + + # Compare reaction type, energy groups, nuclides, domain type + if self.rxn_type != other.rxn_type: + return False + elif not self.energy_groups.can_merge(other.energy_groups): + return False + elif self.by_nuclide != other.by_nuclide: + return False + elif self.domain_type != other.domain_type: + return False + elif 'distribcell' not in self.domain_type and self.domain != other.domain: + return False + elif not self.xs_tally.can_merge(other.xs_tally): + return False + elif not self.rxn_rate_tally.can_merge(other.rxn_rate_tally): + return False + + # If all conditionals pass then MGXS are mergeable + return True + + def merge(self, other): + """Merge another MGXS with this one + + MGXS are only mergeable if their energy groups and nuclides are either + identical or mutually exclusive. If results have been loaded from a + statepoint, then MGXS are only mergeable along one and only one of + energy groups or nuclides. + + Parameters + ---------- + other : MGXS + MGXS to merge with this one + + Returns + ------- + merged_mgxs : MGXS + Merged MGXS + + """ + + if not self.can_merge(other): + raise ValueError('Unable to merge MGXS') + + # Create deep copy of tally to return as merged tally + merged_mgxs = copy.deepcopy(self) + merged_mgxs._derived = True + + # Merge energy groups + if self.energy_groups != other.energy_groups: + merged_groups = self.energy_groups.merge(other.energy_groups) + merged_mgxs.energy_groups = merged_groups + + # Merge nuclides + if self.nuclides != other.nuclides: + + # The nuclides must be mutually exclusive + for nuclide in self.nuclides: + if nuclide in other.nuclides: + msg = 'Unable to merge MGXS with shared nuclides' + raise ValueError(msg) + + # Concatenate lists of nuclides for the merged MGXS + merged_mgxs.nuclides = self.nuclides + other.nuclides + + # Null base tallies but merge reaction rate and cross section tallies + merged_mgxs._tallies = OrderedDict() + merged_mgxs._rxn_rate_tally = self.rxn_rate_tally.merge(other.rxn_rate_tally) + merged_mgxs._xs_tally = self.xs_tally.merge(other.xs_tally) + + return merged_mgxs + def print_xs(self, subdomains='all', nuclides='all', xs_type='macro'): """Print a string representation for the multi-group cross section. @@ -1022,6 +1208,9 @@ class MGXS(object): cv.check_iterable_type('subdomains', subdomains, Integral) elif self.domain_type == 'distribcell': subdomains = np.arange(self.num_subdomains, dtype=np.int) + elif self.domain_type == 'avg(distribcell)': + domain_filter = self.xs_tally.find_filter('avg(distribcell)') + subdomains = domain_filter.bins else: subdomains = [self.domain.id] @@ -1239,14 +1428,13 @@ class MGXS(object): df.rename(columns={'energy low [MeV]': 'group in'}, inplace=True) in_groups = np.tile(all_groups, self.num_subdomains) - in_groups = np.repeat(in_groups, self.num_groups) + in_groups = np.repeat(in_groups, df.shape[0] / in_groups.size) df['group in'] = in_groups del df['energy high [MeV]'] df.rename(columns={'energyout low [MeV]': 'group out'}, inplace=True) - out_groups = \ - np.tile(all_groups, self.num_subdomains * self.num_groups) + out_groups = np.tile(all_groups, df.shape[0] / all_groups.size) df['group out'] = out_groups del df['energyout high [MeV]'] columns = ['group in', 'group out'] @@ -1285,8 +1473,7 @@ class MGXS(object): # Sort the dataframe by domain type id (e.g., distribcell id) and # energy groups such that data is from fast to thermal - df.sort([self.domain_type] + columns, inplace=True) - + df.sort_values(by=[self.domain_type] + columns, inplace=True) return df @@ -1329,7 +1516,7 @@ class TotalXS(MGXS): @property def rxn_rate_tally(self): - if self._rxn_rate_tally is None: + if self._rxn_rate_tally is None : self._rxn_rate_tally = self.tallies['total'] self._rxn_rate_tally.sparse = self.sparse return self._rxn_rate_tally @@ -1735,6 +1922,58 @@ class ScatterMatrixXS(MGXS): cv.check_value('correction', correction, ('P0', None)) self._correction = correction + def get_slice(self, nuclides=[], in_groups=[], out_groups=[]): + """Build a sliced ScatterMatrix for the specified nuclides and + energy groups. + + This method constructs a new MGXS to encapsulate a subset of the data + represented by this MGXS. The subset of data to include in the tally + slice is determined by the nuclides and energy groups specified in + the input parameters. + + Parameters + ---------- + nuclides : list of str + A list of nuclide name strings + (e.g., ['U-235', 'U-238']; default is []) + in_groups : list of Integral + A list of incoming energy group indices starting at 1 for the high + energies (e.g., [1, 2, 3]; default is []) + out_groups : list of Integral + A list of outgoing energy group indices starting at 1 for the high + energies (e.g., [1, 2, 3]; default is []) + + Returns + ------- + MGXS + A new tally which encapsulates the subset of data requested for the + nuclide(s) and/or energy group(s) requested in the parameters. + + """ + + # Call super class method and null out derived tallies + slice_xs = super(ScatterMatrixXS, self).get_slice(nuclides, in_groups) + slice_xs._rxn_rate_tally = None + slice_xs._xs_tally = None + + # Slice outgoing energy groups if needed + if len(out_groups) != 0: + filter_bins = [] + for group in out_groups: + group_bounds = self.energy_groups.get_group_bounds(group) + filter_bins.append(group_bounds) + filter_bins = [tuple(filter_bins)] + + # Slice each of the tallies across energyout groups + for tally_type, tally in slice_xs.tallies.items(): + if tally.contains_filter('energyout'): + tally_slice = tally.get_slice(filters=['energyout'], + filter_bins=filter_bins) + slice_xs.tallies[tally_type] = tally_slice + + slice_xs.sparse = self.sparse + return slice_xs + def get_xs(self, in_groups='all', out_groups='all', subdomains='all', nuclides='all', xs_type='macro', order_groups='increasing', value='mean'): @@ -1789,7 +2028,7 @@ class ScatterMatrixXS(MGXS): # Construct a collection of the domain filter bins if not isinstance(subdomains, basestring): - cv.check_iterable_type('subdomains', subdomains, Integral) + cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=2) for subdomain in subdomains: filters.append(self.domain_type) filter_bins.append((subdomain,)) @@ -2078,6 +2317,113 @@ class Chi(MGXS): return self._xs_tally + def get_slice(self, nuclides=[], groups=[]): + """Build a sliced Chi for the specified nuclides and energy groups. + + This method constructs a new MGXS to encapsulate a subset of the data + represented by this MGXS. The subset of data to include in the tally + slice is determined by the nuclides and energy groups specified in + the input parameters. + + Parameters + ---------- + nuclides : list of str + A list of nuclide name strings + (e.g., ['U-235', 'U-238']; default is []) + groups : list of Integral + A list of energy group indices starting at 1 for the high energies + (e.g., [1, 2, 3]; default is []) + + Returns + ------- + MGXS + A new tally which encapsulates the subset of data requested for the + nuclide(s) and/or energy group(s) requested in the parameters. + + """ + + # Temporarily remove energy filter from nu-fission-in since its + # group structure will work in super MGXS.get_slice(...) method + nu_fission_in = self.tallies['nu-fission-in'] + energy_filter = nu_fission_in.find_filter('energy') + nu_fission_in.remove_filter(energy_filter) + + # Call super class method and null out derived tallies + slice_xs = super(Chi, self).get_slice(nuclides, groups) + slice_xs._rxn_rate_tally = None + slice_xs._xs_tally = None + + # Slice energy groups if needed + if len(groups) != 0: + filter_bins = [] + for group in groups: + group_bounds = self.energy_groups.get_group_bounds(group) + filter_bins.append(group_bounds) + filter_bins = [tuple(filter_bins)] + + # Slice nu-fission-out tally along energyout filter + nu_fission_out = slice_xs.tallies['nu-fission-out'] + tally_slice = nu_fission_out.get_slice(filters=['energyout'], + filter_bins=filter_bins) + slice_xs._tallies['nu-fission-out'] = tally_slice + + # Add energy filter back to nu-fission-in tallies + self.tallies['nu-fission-in'].add_filter(energy_filter) + slice_xs._tallies['nu-fission-in'].add_filter(energy_filter) + + slice_xs.sparse = self.sparse + return slice_xs + + def merge(self, other): + """Merge another Chi with this one + + If results have been loaded from a statepoint, then Chi are only + mergeable along one and only one of energy groups or nuclides. + + Parameters + ---------- + other : MGXS + MGXS to merge with this one + + Returns + ------- + merged_mgxs : MGXS + Merged MGXS + """ + + if not self.can_merge(other): + raise ValueError('Unable to merge Chi') + + # Create deep copy of tally to return as merged tally + merged_mgxs = copy.deepcopy(self) + merged_mgxs._derived = True + merged_mgxs._rxn_rate_tally = None + merged_mgxs._xs_tally = None + + # Merge energy groups + if self.energy_groups != other.energy_groups: + merged_groups = self.energy_groups.merge(other.energy_groups) + merged_mgxs.energy_groups = merged_groups + + # Merge nuclides + if self.nuclides != other.nuclides: + + # The nuclides must be mutually exclusive + for nuclide in self.nuclides: + if nuclide in other.nuclides: + msg = 'Unable to merge Chi with shared nuclides' + raise ValueError(msg) + + # Concatenate lists of nuclides for the merged MGXS + merged_mgxs.nuclides = self.nuclides + other.nuclides + + # Merge tallies + for tally_key in self.tallies: + merged_tally = self.tallies[tally_key].merge(other.tallies[tally_key]) + merged_mgxs.tallies[tally_key] = merged_tally + + return merged_mgxs + def get_xs(self, groups='all', subdomains='all', nuclides='all', xs_type='macro', order_groups='increasing', value='mean'): """Returns an array of the fission spectrum. @@ -2129,7 +2475,7 @@ class Chi(MGXS): # Construct a collection of the domain filter bins if not isinstance(subdomains, basestring): - cv.check_iterable_type('subdomains', subdomains, Integral) + cv.check_iterable_type('subdomains', subdomains, Integral, max_depth=2) for subdomain in subdomains: filters.append(self.domain_type) filter_bins.append((subdomain,)) diff --git a/openmc/nuclide.py b/openmc/nuclide.py index 01fb2aa459..8e97f1a1c6 100644 --- a/openmc/nuclide.py +++ b/openmc/nuclide.py @@ -60,6 +60,12 @@ class Nuclide(object): def __ne__(self, other): return not self == other + def __gt__(self, other): + return repr(self) > repr(other) + + def __lt__(self, other): + return not self > other + def __hash__(self): return hash(repr(self)) diff --git a/openmc/tallies.py b/openmc/tallies.py index d1666694fe..fb66eb77f8 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -302,7 +302,7 @@ class Tally(object): @property def sum(self): - if not self._sp_filename: + if not self._sp_filename or self.derived: return None if not self._results_read: @@ -674,72 +674,193 @@ class Tally(object): self._nuclides.remove(nuclide) - def can_merge(self, tally): - """Determine if another tally can be merged with this one + def _can_merge_filters(self, other): + """Determine if another tally's filters can be merged with this one's + + The types of filters between the two tallies must match identically. + The bins in all of the filters must match identically, or be mergeable + in only one filter. This is a helper method for the can_merge(...) + and merge(...) methods. Parameters ---------- - tally : Tally - Tally to check for merging + other : Tally + Tally to check for mergeable filters """ - if not isinstance(tally, Tally): + # Two tallys must have the same number of filters + if len(self.filters) != len(other.filters): return False - # Must have same estimator - if self.estimator != tally.estimator: - return False - - # Must have same nuclides - if len(self.nuclides) != len(tally.nuclides): - return False - - for nuclide in self.nuclides: - if nuclide not in tally.nuclides: - return False - - # Must have same or mergeable filters - if len(self.filters) != len(tally.filters): - return False - - # Check if only one tally contains a delayed group filter - tally1_dg = False - for filter1 in self.filters: - if filter1.type == 'delayedgroup': - tally1_dg = True - - tally2_dg = False - for filter2 in tally.filters: - if filter2.type == 'delayedgroup': - tally2_dg = True - # Return False if only one tally has a delayed group filter - if (tally1_dg or tally2_dg) and not (tally1_dg and tally2_dg): + tally1_dg = self.contains_filter('delayedgroup') + tally2_dg = other.contains_filter('delayedgroup') + if sum([tally1_dg, tally2_dg]) == 1: return False # Look to see if all filters are the same, or one or more can be merged for filter1 in self.filters: + merge_filters = False mergeable_filter = False - for filter2 in tally.filters: - if filter1 == filter2 or filter1.can_merge(filter2): + for filter2 in other.filters: + + # If filters match, they are mergeable + if filter1 == filter2: mergeable_filter = True break + # If filters are first mergeable filters encountered + elif filter1.can_merge(filter2) and not merge_filters: + merge_filters = True + mergeable_filter = True + break + + # If filters are the second mergeable filters encountered + elif filter1.can_merge(filter2) and merge_filters: + return False + # If no mergeable filter was found, the tallies are not mergeable if not mergeable_filter: return False - # Tallies are mergeable if all conditional checks passed + # Tally filters are mergeable if all conditional checks passed return True - def merge(self, tally): - """Merge another tally with this one + def _can_merge_nuclides(self, other): + """Determine if another tally's nuclides can be merged with this one's + + The nuclides between the two tallies must be mutually exclusive or + identically matching. This is a helper method for the can_merge(...) + and merge(...) methods. Parameters ---------- - tally : Tally + other : Tally + Tally to check for mergeable nuclides + + """ + + no_nuclides_match = True + all_nuclides_match = True + + # Search for each of this tally's nuclides in the other tally + for nuclide in self.nuclides: + if nuclide not in other.nuclides: + all_nuclides_match = False + else: + no_nuclides_match = False + + # Search for each of the other tally's nuclides in this tally + for nuclide in other.nuclides: + if nuclide not in self.nuclides: + all_nuclides_match = False + else: + no_nuclides_match = False + + # Either all nuclides should match, or none should + if no_nuclides_match or all_nuclides_match: + return True + else: + return False + + def _can_merge_scores(self, other): + """Determine if another tally's scores can be merged with this one's + + The scores between the two tallies must be mutually exclusive or + identically matching. This is a helper method for the can_merge(...) + and merge(...) methods. + + Parameters + ---------- + other : Tally + Tally to check for mergeable scores + + """ + + no_scores_match = True + all_scores_match = True + + # Search for each of this tally's scores in the other tally + for score in self.scores: + if score not in other.scores: + all_scores_match = False + else: + no_scores_match = False + + # Search for each of the other tally's scores in this tally + for score in other.scores: + if score not in self.scores: + all_scores_match = False + else: + no_scores_match = False + + # Nuclides cannot be specified on 'flux' scores + if 'flux' in self.scores or 'flux' in other.scores: + if self.nuclides != other.nuclides: + return False + + # Either all scores should match, or none should + if no_scores_match or all_scores_match: + return True + else: + return False + + def can_merge(self, other): + """Determine if another tally can be merged with this one + + If results have been loaded from a statepoint, then tallies are only + mergeable along one and only one of filter bins, nuclides or scores. + + Parameters + ---------- + other : Tally + Tally to check for merging + + """ + + if not isinstance(other, Tally): + return False + + # Must have same estimator + if self.estimator != other.estimator: + return False + + equal_filters = sorted(self.filters) == sorted(other.filters) + equal_nuclides = sorted(self.nuclides) == sorted(other.nuclides) + equal_scores = sorted(self.scores) == sorted(other.scores) + equality = [equal_filters, equal_nuclides, equal_scores] + + # If all filters, nuclides and scores match then tallies are mergeable + if equal_filters and equal_nuclides and equal_scores: + return True + + # Variables to indicate matching filter bins, nuclides and scores + merge_filters = self._can_merge_filters(other) + merge_nuclides = self._can_merge_nuclides(other) + merge_scores = self._can_merge_scores(other) + mergeability = [merge_filters, merge_nuclides, merge_scores] + + if not all(mergeability): + return False + + # If the tally results have been read from the statepoint, we can only + # at least two of filters, nuclides and scores must match + elif self._results_read and sum(equality) < 2: + return False + else: + return True + + def merge(self, other): + """Merge another tally with this one + + If results have been loaded from a statepoint, then tallies are only + mergeable along one and only one of filter bins, nuclides or scores. + + Parameters + ---------- + other : Tally Tally to merge with this one Returns @@ -749,9 +870,9 @@ class Tally(object): """ - if not self.can_merge(tally): - msg = 'Unable to merge tally ID="{0}" with ' + \ - '"{1}"'.format(tally.id, self.id) + if not self.can_merge(other): + msg = 'Unable to merge tally ID="{0}" with ' \ + '"{1}"'.format(other.id, self.id) raise ValueError(msg) # Create deep copy of tally to return as merged tally @@ -760,23 +881,127 @@ class Tally(object): # Differentiate Tally with a new auto-generated Tally ID merged_tally.id = None - # Merge filters - for i, filter1 in enumerate(merged_tally.filters): - for filter2 in tally.filters: - if filter1 != filter2 and filter1.can_merge(filter2): - merged_filter = filter1.merge(filter2) - merged_tally.filters[i] = merged_filter - break + # If the two tallies are equal, simply return copy + if self == other: + return merged_tally - # Add unique scores from second tally to merged tally - for score in tally.scores: - if score not in merged_tally.scores: - merged_tally.add_score(score) + # Create deep copy of other tally to use for array concatenation + other_copy = copy.deepcopy(other) - # Add triggers from second tally to merged tally - for trigger in tally.triggers: + # Identify if filters, nuclides and scores are mergeable and/or equal + merge_filters = self._can_merge_filters(other) + merge_nuclides = self._can_merge_nuclides(other) + merge_scores = self._can_merge_scores(other) + equal_filters = sorted(self.filters) == sorted(other.filters) + equal_nuclides = sorted(self.nuclides) == sorted(other.nuclides) + equal_scores = sorted(self.scores) == sorted(other.scores) + + # If two tallies can be merged along a filter's bins + if merge_filters and not equal_filters: + + # Search for mergeable filters + for i, filter1 in enumerate(self.filters): + for j, filter2 in enumerate(other.filters): + if filter1 != filter2 and filter1.can_merge(filter2): + other_copy._swap_filters(other_copy.filters[i], filter2) + merged_tally.filters[i] = filter1.merge(filter2) + join_right = filter1 < filter2 + merge_axis = i + break + + # If two tallies can be merged along nuclide bins + if merge_nuclides and not equal_nuclides: + merge_axis = self.num_filters + join_right = True + + # Add unique nuclides from other tally to merged tally + for nuclide in other.nuclides: + if nuclide not in merged_tally.nuclides: + merged_tally.add_nuclide(nuclide) + + # If two tallies can be merged along score bins + if merge_scores and not equal_scores: + merge_axis = self.num_filters + 1 + join_right = True + + # Add unique scores from other tally to merged tally + for score in other.scores: + if score not in merged_tally.scores: + merged_tally.add_score(score) + + # Add triggers from other tally to merged tally + for trigger in other.triggers: merged_tally.add_trigger(trigger) + # If results have not been read, then return tally for input generation + if self._results_read is None: + return merged_tally + # Otherwise, this is a derived tally which needs merged results arrays + else: + self._derived = True + + # Update filter strides in merged tally + merged_tally._update_filter_strides() + + # Concatenate sum arrays if present in both tallies + if self.sum is not None and other_copy.sum is not None: + self_sum = self.get_reshaped_data(value='sum') + other_sum = other_copy.get_reshaped_data(value='sum') + + if join_right: + merged_sum = \ + np.concatenate((self_sum, other_sum), axis=merge_axis) + else: + merged_sum = \ + np.concatenate((other_sum, self_sum), axis=merge_axis) + + merged_tally._sum = np.reshape(merged_sum, merged_tally.shape) + + # Concatenate sum_sq arrays if present in both tallies + if self.sum_sq is not None and other.sum_sq is not None: + self_sum_sq = self.get_reshaped_data(value='sum_sq') + other_sum_sq = other_copy.get_reshaped_data(value='sum_sq') + + if join_right: + merged_sum_sq = \ + np.concatenate((self_sum_sq, other_sum_sq), axis=merge_axis) + else: + merged_sum_sq = \ + np.concatenate((other_sum_sq, self_sum_sq), axis=merge_axis) + + merged_tally._sum_sq = np.reshape(merged_sum_sq, merged_tally.shape) + + # Concatenate mean arrays if present in both tallies + if self.mean is not None and other.mean is not None: + self_mean = self.get_reshaped_data(value='mean') + other_mean = other_copy.get_reshaped_data(value='mean') + + if join_right: + merged_mean = \ + np.concatenate((self_mean, other_mean), axis=merge_axis) + else: + merged_mean = \ + np.concatenate((other_mean, self_mean), axis=merge_axis) + + merged_tally._mean = np.reshape(merged_mean, merged_tally.shape) + + # Concatenate std. dev. arrays if present in both tallies + if self.std_dev is not None and other.std_dev is not None: + self_std_dev = self.get_reshaped_data(value='std_dev') + other_std_dev = other_copy.get_reshaped_data(value='std_dev') + + if join_right: + merged_std_dev = \ + np.concatenate((self_std_dev, other_std_dev), axis=merge_axis) + else: + merged_std_dev = \ + np.concatenate((other_std_dev, self_std_dev), axis=merge_axis) + + merged_tally._std_dev = np.reshape(merged_std_dev, merged_tally.shape) + + # Sparsify merged tally if both tallies are sparse + merged_tally.sparse = self.sparse and other.sparse + return merged_tally def get_tally_xml(self): @@ -847,6 +1072,32 @@ class Tally(object): return element + def contains_filter(self, filter_type): + """Looks for a filter in the tally that matches a specified type + + Parameters + ---------- + filter_type : str + Type of the filter, e.g. 'mesh' + + Returns + ------- + filter_found : bool + True if the tally contains a filter of the requested type; + otherwise false + + """ + + filter_found = False + + # Look through all of this Tally's Filters for the type requested + for test_filter in self.filters: + if test_filter.type == filter_type: + filter_found = True + break + + return filter_found + def find_filter(self, filter_type): """Return a filter in the tally that matches a specified type @@ -1302,14 +1553,8 @@ class Tally(object): 'Summary info'.format(self.id) raise KeyError(msg) - # Attempt to import Pandas - try: - import pandas as pd - except ImportError: - msg = 'The Pandas Python package must be installed on your system' - raise ImportError(msg) - # Initialize a pandas dataframe for the tally data + import pandas as pd df = pd.DataFrame() # Find the total length of the tally data array @@ -1326,23 +1571,36 @@ class Tally(object): # Include DataFrame column for nuclides if user requested it if nuclides: nuclides = [] + column_name = 'nuclide' for nuclide in self.nuclides: - # Write Nuclide name if Summary info was linked with StatePoint if isinstance(nuclide, Nuclide): nuclides.append(nuclide.name) + elif isinstance(nuclide, AggregateNuclide): + nuclides.append(nuclide.name) + column_name = '{0}(nuclide)'.format(nuclide.aggregate_op) else: nuclides.append(nuclide) # Tile the nuclide bins into a DataFrame column nuclides = np.repeat(nuclides, len(self.scores)) tile_factor = data_size / len(nuclides) - df['nuclide'] = np.tile(nuclides, int(tile_factor)) + df[column_name] = np.tile(nuclides, int(tile_factor)) # Include column for scores if user requested it if scores: + scores = [] + column_name = 'score' + + for score in self.scores: + if isinstance(score, (basestring, CrossScore)): + scores.append(score) + elif isinstance(score, AggregateScore): + scores.append(score.name) + column_name = '{0}(score)'.format(score.aggregate_op) + tile_factor = data_size / len(self.scores) - df['score'] = np.tile(self.scores, int(tile_factor)) + df[column_name] = np.tile(scores, int(tile_factor)) # Append columns with mean, std. dev. for each tally bin df['mean'] = self.mean.ravel() @@ -1951,19 +2209,12 @@ class Tally(object): """ - # Check that results have been read - if not self.derived and self.sum is None: - msg = 'Unable to use tally arithmetic with Tally ID="{0}" ' \ - 'since it does not contain any results.'.format(self.id) - raise ValueError(msg) - - cv.check_type('filter1', filter1, Filter) - cv.check_type('filter2', filter2, Filter) + cv.check_type('filter1', filter1, (Filter, CrossFilter, AggregateFilter)) + cv.check_type('filter2', filter2, (Filter, CrossFilter, AggregateFilter)) # Check that the filters exist in the tally and are not the same if filter1 == filter2: - msg = 'Unable to swap a filter with itself' - raise ValueError(msg) + return elif filter1 not in self.filters: msg = 'Unable to swap "{0}" filter1 in Tally ID="{1}" since it ' \ 'does not contain such a filter'.format(filter1.type, self.id) @@ -2684,7 +2935,13 @@ class Tally(object): 'since it does not contain any results.'.format(self.id) raise ValueError(msg) + # Create deep copy of tally to return as sliced tally new_tally = copy.deepcopy(self) + new_tally._derived = True + + # Differentiate Tally with a new auto-generated Tally ID + new_tally.id = None + new_tally.sparse = False if not self.derived and self.sum is not None: @@ -2769,7 +3026,7 @@ class Tally(object): def summation(self, scores=[], filter_type=None, filter_bins=[], nuclides=[], remove_filter=False): """Vectorized sum of tally data across scores, filter bins and/or - nuclides using tally addition. + nuclides using tally aggregation. This method constructs a new tally to encapsulate the sum of the data represented by the summation of the data in this tally. The tally data @@ -2811,7 +3068,7 @@ class Tally(object): tally_sum._derived = True tally_sum._estimator = self.estimator tally_sum._num_realizations = self.num_realizations - tally_sum.with_batch_statistics = self.with_batch_statistics + tally_sum._with_batch_statistics = self.with_batch_statistics tally_sum._with_summary = self.with_summary tally_sum._sp_filename = self._sp_filename tally_sum._results_read = self._results_read @@ -2852,7 +3109,7 @@ class Tally(object): # Add AggregateFilter to the tally sum if not remove_filter: filter_sum = \ - AggregateFilter(self_filter, filter_bins, 'sum') + AggregateFilter(self_filter, [tuple(filter_bins)], 'sum') tally_sum.add_filter(filter_sum) # Add a copy of each filter not summed across to the tally sum @@ -2914,6 +3171,157 @@ class Tally(object): tally_sum.sparse = self.sparse return tally_sum + def average(self, scores=[], filter_type=None, + filter_bins=[], nuclides=[], remove_filter=False): + """Vectorized average of tally data across scores, filter bins and/or + nuclides using tally aggregation. + + This method constructs a new tally to encapsulate the average of the + data represented by the average of the data in this tally. The tally + data average is determined by the scores, filter bins and nuclides + specified in the input parameters. + + Parameters + ---------- + scores : list of str + A list of one or more score strings to average across + (e.g., ['absorption', 'nu-fission']; default is []) + filter_type : str + A filter type string (e.g., 'cell', 'energy') corresponding to the + filter bins to average across + filter_bins : Iterable of Integral or tuple + A list of the filter bins corresponding to the filter_type parameter + Each bin in the list is the integer ID for 'material', 'surface', + 'cell', 'cellborn', and 'universe' Filters. Each bin is an integer + for the cell instance ID for 'distribcell' Filters. Each bin is a + 2-tuple of floats for 'energy' and 'energyout' filters corresponding + to the energy boundaries of the bin of interest. Each bin is an + (x,y,z) 3-tuple for 'mesh' filters corresponding to the mesh cell of + interest. + nuclides : list of str + A list of nuclide name strings to average across + (e.g., ['U-235', 'U-238']; default is []) + remove_filter : bool + If a filter is being averaged over, this bool indicates whether to + remove that filter in the returned tally. Default is False. + + Returns + ------- + Tally + A new tally which encapsulates the average of data requested. + """ + + # Create new derived Tally for average + tally_avg = Tally() + tally_avg._derived = True + tally_avg._estimator = self.estimator + tally_avg._num_realizations = self.num_realizations + tally_avg._with_batch_statistics = self.with_batch_statistics + tally_avg._with_summary = self.with_summary + tally_avg._sp_filename = self._sp_filename + tally_avg._results_read = self._results_read + + # Get tally data arrays reshaped with one dimension per filter + mean = self.get_reshaped_data(value='mean') + std_dev = self.get_reshaped_data(value='std_dev') + + # Average across any filter bins specified by the user + if filter_type in _FILTER_TYPES: + find_filter = self.find_filter(filter_type) + + # If user did not specify filter bins, average across all bins + if len(filter_bins) == 0: + bin_indices = np.arange(find_filter.num_bins) + + if filter_type == 'distribcell': + filter_bins = np.arange(find_filter.num_bins) + else: + num_bins = find_filter.num_bins + filter_bins = \ + [(find_filter.get_bin(i)) for i in range(num_bins)] + + # Only average across bins specified by the user + else: + bin_indices = \ + [find_filter.get_bin_index(bin) for bin in filter_bins] + + # Average across the bins in the user-specified filter + for i, self_filter in enumerate(self.filters): + if self_filter.type == filter_type: + mean = np.take(mean, indices=bin_indices, axis=i) + std_dev = np.take(std_dev, indices=bin_indices, axis=i) + mean = np.mean(mean, axis=i, keepdims=True) + std_dev = np.mean(std_dev**2, axis=i, keepdims=True) + std_dev /= len(bin_indices) + std_dev = np.sqrt(std_dev) + + # Add AggregateFilter to the tally avg + if not remove_filter: + filter_sum = \ + AggregateFilter(self_filter, [tuple(filter_bins)], 'avg') + tally_avg.add_filter(filter_sum) + + # Add a copy of each filter not averaged across to the tally avg + else: + tally_avg.add_filter(copy.deepcopy(self_filter)) + + # Add a copy of this tally's filters to the tally avg + else: + tally_avg._filters = copy.deepcopy(self.filters) + + # Sum across any nuclides specified by the user + if len(nuclides) != 0: + nuclide_bins = [self.get_nuclide_index(nuclide) for nuclide in nuclides] + axis_index = self.num_filters + mean = np.take(mean, indices=nuclide_bins, axis=axis_index) + std_dev = np.take(std_dev, indices=nuclide_bins, axis=axis_index) + mean = np.mean(mean, axis=axis_index, keepdims=True) + std_dev = np.mean(std_dev**2, axis=axis_index, keepdims=True) + std_dev /= len(nuclide_bins) + std_dev = np.sqrt(std_dev) + + # Add AggregateNuclide to the tally avg + nuclide_avg = AggregateNuclide(nuclides, 'avg') + tally_avg.add_nuclide(nuclide_avg) + + # Add a copy of this tally's nuclides to the tally avg + else: + tally_avg._nuclides = copy.deepcopy(self.nuclides) + + # Sum across any scores specified by the user + if len(scores) != 0: + score_bins = [self.get_score_index(score) for score in scores] + axis_index = self.num_filters + 1 + mean = np.take(mean, indices=score_bins, axis=axis_index) + std_dev = np.take(std_dev, indices=score_bins, axis=axis_index) + mean = np.sum(mean, axis=axis_index, keepdims=True) + std_dev = np.sum(std_dev**2, axis=axis_index, keepdims=True) + std_dev /= len(score_bins) + std_dev = np.sqrt(std_dev) + + # Add AggregateScore to the tally avg + score_sum = AggregateScore(scores, 'avg') + tally_avg.add_score(score_sum) + + # Add a copy of this tally's scores to the tally avg + else: + tally_avg._scores = copy.deepcopy(self.scores) + + # Update the tally avg's filter strides + tally_avg._update_filter_strides() + + # Reshape condensed data arrays with one dimension for all filters + mean = np.reshape(mean, tally_avg.shape) + std_dev = np.reshape(std_dev, tally_avg.shape) + + # Assign tally avg's data with the new arrays + tally_avg._mean = mean + tally_avg._std_dev = std_dev + + # If original tally was sparse, sparsify the tally average + tally_avg.sparse = self.sparse + return tally_avg + def diagonalize_filter(self, new_filter): """Diagonalize the tally data array along a new axis of filter bins. diff --git a/src/distribution_multivariate.F90 b/src/distribution_multivariate.F90 index 694fac3013..68b246df84 100644 --- a/src/distribution_multivariate.F90 +++ b/src/distribution_multivariate.F90 @@ -16,15 +16,15 @@ module distribution_multivariate type, abstract :: UnitSphereDistribution real(8) :: reference_uvw(3) contains - procedure(iSample), deferred :: sample + procedure(unitsphere_distribution_sample_), deferred :: sample end type UnitSphereDistribution abstract interface - function iSample(this) result(uvw) + function unitsphere_distribution_sample_(this) result(uvw) import UnitSphereDistribution class(UnitSphereDistribution), intent(in) :: this real(8) :: uvw(3) - end function iSample + end function unitsphere_distribution_sample_ end interface !=============================================================================== @@ -58,15 +58,15 @@ module distribution_multivariate type, abstract :: SpatialDistribution contains - procedure(iSampleSpatial), deferred :: sample + procedure(spatial_distribution_sample_), deferred :: sample end type SpatialDistribution abstract interface - function iSampleSpatial(this) result(xyz) + function spatial_distribution_sample_(this) result(xyz) import SpatialDistribution class(SpatialDistribution), intent(in) :: this real(8) :: xyz(3) - end function iSampleSpatial + end function spatial_distribution_sample_ end interface type, extends(SpatialDistribution) :: CartesianIndependent diff --git a/src/distribution_univariate.F90 b/src/distribution_univariate.F90 index f3e4fdee2e..6c053e5942 100644 --- a/src/distribution_univariate.F90 +++ b/src/distribution_univariate.F90 @@ -16,7 +16,7 @@ module distribution_univariate type, abstract :: Distribution contains - procedure(iSample), deferred :: sample + procedure(distribution_sample_), deferred :: sample end type Distribution type DistributionContainer @@ -24,11 +24,11 @@ module distribution_univariate end type DistributionContainer abstract interface - function iSample(this) result(x) + function distribution_sample_(this) result(x) import Distribution class(Distribution), intent(in) :: this real(8) :: x - end function iSample + end function distribution_sample_ end interface !=============================================================================== diff --git a/src/energy_distribution.F90 b/src/energy_distribution.F90 index 2cfc1b1841..8b2cc10c9c 100644 --- a/src/energy_distribution.F90 +++ b/src/energy_distribution.F90 @@ -16,16 +16,16 @@ module energy_distribution type, abstract :: EnergyDistribution contains - procedure(iSampleEnergy), deferred :: sample + procedure(energy_distribution_sample_), deferred :: sample end type EnergyDistribution abstract interface - function iSampleEnergy(this, E_in) result(E_out) + function energy_distribution_sample_(this, E_in) result(E_out) import EnergyDistribution class(EnergyDistribution), intent(in) :: this real(8), intent(in) :: E_in real(8) :: E_out - end function iSampleEnergy + end function energy_distribution_sample_ end interface type :: EnergyDistributionContainer diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 519b15bff4..1adda3ea34 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -31,12 +31,10 @@ module geometry_header integer :: outer ! universe to tile outside the lat logical :: is_3d ! Lattice has cells on z axis integer, allocatable :: offset(:,:,:,:) ! Distribcell offsets - - contains - - procedure(are_valid_indices_), deferred :: are_valid_indices - procedure(get_indices_), deferred :: get_indices - procedure(get_local_xyz_), deferred :: get_local_xyz + contains + procedure(lattice_are_valid_indices_), deferred :: are_valid_indices + procedure(lattice_get_indices_), deferred :: get_indices + procedure(lattice_get_local_xyz_), deferred :: get_local_xyz end type Lattice abstract interface @@ -45,33 +43,33 @@ module geometry_header ! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the ! bounds of the lattice. Returns false otherwise. - function are_valid_indices_(this, i_xyz) result(is_valid) + function lattice_are_valid_indices_(this, i_xyz) result(is_valid) import Lattice class(Lattice), intent(in) :: this integer, intent(in) :: i_xyz(3) logical :: is_valid - end function are_valid_indices_ + end function lattice_are_valid_indices_ !=============================================================================== ! GET_INDICES returns the indices in a lattice for the given global xyz. - function get_indices_(this, global_xyz) result(i_xyz) + function lattice_get_indices_(this, global_xyz) result(i_xyz) import Lattice class(Lattice), intent(in) :: this real(8), intent(in) :: global_xyz(3) integer :: i_xyz(3) - end function get_indices_ + end function lattice_get_indices_ !=============================================================================== ! GET_LOCAL_XYZ returns the translated local version of the given global xyz. - function get_local_xyz_(this, global_xyz, i_xyz) result(local_xyz) + function lattice_get_local_xyz_(this, global_xyz, i_xyz) result(local_xyz) import Lattice class(Lattice), intent(in) :: this real(8), intent(in) :: global_xyz(3) integer, intent(in) :: i_xyz(3) real(8) :: local_xyz(3) - end function get_local_xyz_ + end function lattice_get_local_xyz_ end interface !=============================================================================== diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9ef439bb04..03ef8dcbc4 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -131,15 +131,21 @@ contains if (run_CE) then call get_environment_variable("OPENMC_CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then - call fatal_error("No cross_sections.xml file was specified in & - &settings.xml or in the OPENMC_CROSS_SECTIONS environment & - &variable. OpenMC needs such a file to identify where to & - &find ACE cross section libraries. Please consult the user's & - &guide at http://mit-crpg.github.io/openmc for information on & - &how to set up ACE cross section libraries.") - else - path_cross_sections = trim(env_variable) + call get_environment_variable("CROSS_SECTIONS", env_variable) + if (len_trim(env_variable) == 0) then + call fatal_error("No cross_sections.xml file was specified in & + &settings.xml or in the OPENMC_CROSS_SECTIONS environment & + &variable. OpenMC needs such a file to identify where to & + &find ACE cross section libraries. Please consult the user's & + &guide at http://mit-crpg.github.io/openmc for information on & + &how to set up ACE cross section libraries.") + else + call warning("The CROSS_SECTIONS environment variable is & + &deprecated. Please update your environment to use & + &OPENMC_CROSS_SECTIONS instead.") + end if end if + path_cross_sections = trim(env_variable) else call get_environment_variable("OPENMC_MG_CROSS_SECTIONS", env_variable) if (len_trim(env_variable) == 0) then diff --git a/src/nuclide_header.F90 b/src/nuclide_header.F90 index 9bcce57a75..43fea77b65 100644 --- a/src/nuclide_header.F90 +++ b/src/nuclide_header.F90 @@ -33,17 +33,15 @@ module nuclide_header logical :: fissionable ! nuclide is fissionable? contains - procedure(print_nuclide_), deferred :: print ! Writes nuclide info + procedure(nuclide_print_), deferred :: print ! Writes nuclide info end type Nuclide abstract interface - - subroutine print_nuclide_(this, unit) + subroutine nuclide_print_(this, unit) import Nuclide class(Nuclide),intent(in) :: this integer, optional, intent(in) :: unit - end subroutine print_nuclide_ - + end subroutine nuclide_print_ end interface type, extends(Nuclide) :: NuclideCE diff --git a/src/particle_header.F90 b/src/particle_header.F90 index c1f02eca24..4ad4119b76 100644 --- a/src/particle_header.F90 +++ b/src/particle_header.F90 @@ -217,7 +217,7 @@ contains integer, intent(in) :: type logical, intent(in) :: run_CE - integer :: n + integer(8) :: n ! Check to make sure that the hard-limit on secondary particles is not ! exceeded. diff --git a/src/scattdata_header.F90 b/src/scattdata_header.F90 index b04115e570..f8fddbc6b3 100644 --- a/src/scattdata_header.F90 +++ b/src/scattdata_header.F90 @@ -20,22 +20,22 @@ module scattdata_header real(8), allocatable :: data(:,:,:) ! (Order/Nmu x Gout x Gin) contains - procedure(init_), deferred :: init ! Initializes ScattData - procedure(calc_f_), deferred :: calc_f ! Calculates f, given mu - procedure(sample_), deferred :: sample ! sample the scatter event + procedure(scattdata_init_), deferred :: init ! Initializes ScattData + procedure(scattdata_calc_f_), deferred :: calc_f ! Calculates f, given mu + procedure(scattdata_sample_), deferred :: sample ! sample the scatter event end type ScattData abstract interface - subroutine init_(this, order, energy, mult, coeffs) + subroutine scattdata_init_(this, order, energy, mult, coeffs) import ScattData class(ScattData), intent(inout) :: this ! Object to work on integer, intent(in) :: order ! Data Order real(8), intent(in) :: energy(:,:) ! Energy Transfer Matrix real(8), intent(in) :: mult(:,:) ! Scatter Prod'n Matrix real(8), intent(in) :: coeffs(:,:,:) ! Coefficients to use - end subroutine init_ + end subroutine scattdata_init_ - pure function calc_f_(this, gin, gout, mu) result(f) + pure function scattdata_calc_f_(this, gin, gout, mu) result(f) import ScattData class(ScattData), intent(in) :: this ! The ScattData to evaluate integer, intent(in) :: gin ! Incoming Energy Group @@ -43,16 +43,16 @@ module scattdata_header real(8), intent(in) :: mu ! Angle of interest real(8) :: f ! Return value of f(mu) - end function calc_f_ + end function scattdata_calc_f_ - subroutine sample_(this, gin, gout, mu, wgt) + subroutine scattdata_sample_(this, gin, gout, mu, wgt) import ScattData class(ScattData), intent(in) :: this ! Scattering Object to Use integer, intent(in) :: gin ! Incoming neutron group integer, intent(out) :: gout ! Sampled outgoin group real(8), intent(out) :: mu ! Sampled change in angle real(8), intent(inout) :: wgt ! Particle weight - end subroutine sample_ + end subroutine scattdata_sample_ end interface type, extends(ScattData) :: ScattDataLegendre @@ -486,4 +486,4 @@ contains end subroutine scattdatatabular_sample -end module scattdata_header \ No newline at end of file +end module scattdata_header diff --git a/src/secondary_header.F90 b/src/secondary_header.F90 index 7449a5793b..d9a18b6b15 100644 --- a/src/secondary_header.F90 +++ b/src/secondary_header.F90 @@ -14,17 +14,17 @@ module secondary_header type, abstract :: AngleEnergy contains - procedure(iSampleAngleEnergy), deferred :: sample + procedure(angleenergy_sample_), deferred :: sample end type AngleEnergy abstract interface - subroutine iSampleAngleEnergy(this, E_in, E_out, mu) + subroutine angleenergy_sample_(this, E_in, E_out, mu) import AngleEnergy class(AngleEnergy), intent(in) :: this real(8), intent(in) :: E_in real(8), intent(out) :: E_out real(8), intent(out) :: mu - end subroutine iSampleAngleEnergy + end subroutine angleenergy_sample_ end interface type :: AngleEnergyContainer @@ -54,6 +54,7 @@ contains real(8), intent(out) :: E_out ! sampled outgoing energy real(8), intent(out) :: mu ! sampled scattering cosine + integer :: i ! loop counter integer :: n ! number of angle-energy distributions real(8) :: prob ! cumulative probability real(8) :: c ! sampled cumulative probability diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 76562f4937..4686552176 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -19,34 +19,34 @@ module surface_header contains procedure :: sense procedure :: reflect - procedure(iEvaluate), deferred :: evaluate - procedure(iDistance), deferred :: distance - procedure(iNormal), deferred :: normal + procedure(surface_evaluate_), deferred :: evaluate + procedure(surface_distance_), deferred :: distance + procedure(surface_normal_), deferred :: normal end type Surface abstract interface - pure function iEvaluate(this, xyz) result(f) + pure function surface_evaluate_(this, xyz) result(f) import Surface class(Surface), intent(in) :: this real(8), intent(in) :: xyz(3) real(8) :: f - end function iEvaluate + end function surface_evaluate_ - pure function iDistance(this, xyz, uvw, coincident) result(d) + pure function surface_distance_(this, xyz, uvw, coincident) result(d) import Surface class(Surface), intent(in) :: this real(8), intent(in) :: xyz(3) real(8), intent(in) :: uvw(3) logical, intent(in) :: coincident real(8) :: d - end function iDistance + end function surface_distance_ - pure function iNormal(this, xyz) result(uvw) + pure function surface_normal_(this, xyz) result(uvw) import Surface class(Surface), intent(in) :: this real(8), intent(in) :: xyz(3) real(8) :: uvw(3) - end function iNormal + end function surface_normal_ end interface !=============================================================================== diff --git a/src/tally.F90 b/src/tally.F90 index f265bdb274..5a54d9846e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -28,12 +28,12 @@ module tally !$omp threadprivate(position) - procedure(score_general_intfc), pointer :: score_general => null() - procedure(get_scoring_bins_intfc), pointer :: get_scoring_bins => null() + procedure(score_general_), pointer :: score_general => null() + procedure(get_scoring_bins_), pointer :: get_scoring_bins => null() abstract interface - subroutine score_general_intfc(p, t, start_index, filter_index, i_nuclide, & - atom_density, flux) + subroutine score_general_(p, t, start_index, filter_index, i_nuclide, & + atom_density, flux) import Particle import TallyObject type(Particle), intent(in) :: p @@ -43,15 +43,14 @@ module tally integer, intent(in) :: filter_index ! for % results real(8), intent(in) :: flux ! flux estimate real(8), intent(in) :: atom_density ! atom/b-cm - end subroutine score_general_intfc + end subroutine score_general_ - subroutine get_scoring_bins_intfc(p, i_tally, found_bin) + subroutine get_scoring_bins_(p, i_tally, found_bin) import Particle type(Particle), intent(in) :: p integer, intent(in) :: i_tally logical, intent(out) :: found_bin - end subroutine get_scoring_bins_intfc - + end subroutine get_scoring_bins_ end interface contains diff --git a/tests/test_asymmetric_lattice/test_asymmetric_lattice.py b/tests/test_asymmetric_lattice/test_asymmetric_lattice.py index 5a1d47ef85..fdb21db33c 100644 --- a/tests/test_asymmetric_lattice/test_asymmetric_lattice.py +++ b/tests/test_asymmetric_lattice/test_asymmetric_lattice.py @@ -19,13 +19,10 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness): # Build full core geometry from underlying input set self._input_set.build_default_materials_and_geometry() - # Extract all universes from the full core geometry - geometry = self._input_set.geometry.geometry - all_univs = geometry.get_all_universes() - # Extract universes encapsulating fuel and water assemblies - water = all_univs[7] - fuel = all_univs[8] + geometry = self._input_set.geometry.geometry + water = geometry.get_universes_by_name('water assembly (hot)')[0] + fuel = geometry.get_universes_by_name('fuel assembly (hot)')[0] # Construct a 3x3 lattice of fuel assemblies core_lat = openmc.RectLattice(name='3x3 Core Lattice', lattice_id=202) @@ -102,9 +99,10 @@ class AsymmetricLatticeTestHarness(PyAPITestHarness): outstr += ', '.join(map(str, tally.std_dev.flatten())) + '\n' # Extract fuel assembly lattices from the summary - all_cells = su.openmc_geometry.get_all_cells() - fuel = all_cells[80].fill - core = all_cells[1].fill + core = su.get_cell_by_id(1) + fuel = su.get_cell_by_id(80) + fuel = fuel.fill + core = core.fill # Append a string of lattice distribcell offsets to the string outstr += ', '.join(map(str, fuel.offsets.flatten())) + '\n' diff --git a/tests/test_mgxs_library_distribcell/results_true.dat b/tests/test_mgxs_library_distribcell/results_true.dat index 4936da4cec..ad9b949de9 100644 --- a/tests/test_mgxs_library_distribcell/results_true.dat +++ b/tests/test_mgxs_library_distribcell/results_true.dat @@ -1,5 +1,5 @@ - sum(distribcell) group in nuclide mean std. dev. -0 sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ... 1 total 0.720213 1.424323 sum(distribcell) group in nuclide mean std. dev. -0 sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ... 1 total 0 0 sum(distribcell) group in group out nuclide mean std. dev. -0 sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ... 1 1 total 0.70466 1.403916 sum(distribcell) group out nuclide mean std. dev. -0 sum(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, ... 1 total 0 0 \ No newline at end of file + avg(distribcell) group in nuclide mean std. dev. +0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0.720213 1.424323 avg(distribcell) group in nuclide mean std. dev. +0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0 0 avg(distribcell) group in group out nuclide mean std. dev. +0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 1 total 0.70466 1.403916 avg(distribcell) group out nuclide mean std. dev. +0 (0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13,... 1 total 0 0 \ No newline at end of file diff --git a/tests/test_tally_slice_merge/inputs_true.dat b/tests/test_tally_slice_merge/inputs_true.dat new file mode 100644 index 0000000000..29f0f1d827 --- /dev/null +++ b/tests/test_tally_slice_merge/inputs_true.dat @@ -0,0 +1 @@ +8d1ab9e4add51b99045e990ac9c3dad9447e9720d811bc430d4bfdd7c2c035424bcb7750e4a4d0ec0460ea1ef4be46ac58372ed01d55f5d8cfeebbce75559066 \ No newline at end of file diff --git a/tests/test_tally_slice_merge/results_true.dat b/tests/test_tally_slice_merge/results_true.dat new file mode 100644 index 0000000000..f66f174277 --- /dev/null +++ b/tests/test_tally_slice_merge/results_true.dat @@ -0,0 +1,49 @@ + energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 21 U-235 fission 9.86e-02 9.19e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 21 U-235 nu-fission 2.40e-01 2.24e-02 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 21 U-238 fission 1.37e-07 1.28e-08 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 21 U-238 nu-fission 3.42e-07 3.20e-08 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 21 U-235 fission 2.79e-02 6.02e-04 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 21 U-235 nu-fission 6.82e-02 1.46e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 21 U-238 fission 1.66e-02 1.15e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 21 U-238 nu-fission 4.58e-02 3.34e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 27 U-235 fission 5.78e-02 4.82e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 27 U-235 nu-fission 1.41e-01 1.17e-02 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 27 U-238 fission 8.18e-08 7.06e-09 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 27 U-238 nu-fission 2.04e-07 1.76e-08 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 27 U-235 fission 1.76e-02 1.94e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 27 U-235 nu-fission 4.31e-02 4.74e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 27 U-238 fission 9.88e-03 1.93e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 6.25e-07 2.00e+01 27 U-238 nu-fission 2.71e-02 5.21e-03 energy low [MeV] energy high [MeV] cell nuclide score mean std. dev. +0 0.00e+00 6.25e-07 21 U-235 fission 9.86e-02 9.19e-03 +1 0.00e+00 6.25e-07 21 U-235 nu-fission 2.40e-01 2.24e-02 +2 0.00e+00 6.25e-07 21 U-238 fission 1.37e-07 1.28e-08 +3 0.00e+00 6.25e-07 21 U-238 nu-fission 3.42e-07 3.20e-08 +4 0.00e+00 6.25e-07 27 U-235 fission 5.78e-02 4.82e-03 +5 0.00e+00 6.25e-07 27 U-235 nu-fission 1.41e-01 1.17e-02 +6 0.00e+00 6.25e-07 27 U-238 fission 8.18e-08 7.06e-09 +7 0.00e+00 6.25e-07 27 U-238 nu-fission 2.04e-07 1.76e-08 +8 6.25e-07 2.00e+01 21 U-235 fission 2.79e-02 6.02e-04 +9 6.25e-07 2.00e+01 21 U-235 nu-fission 6.82e-02 1.46e-03 +10 6.25e-07 2.00e+01 21 U-238 fission 1.66e-02 1.15e-03 +11 6.25e-07 2.00e+01 21 U-238 nu-fission 4.58e-02 3.34e-03 +12 6.25e-07 2.00e+01 27 U-235 fission 1.76e-02 1.94e-03 +13 6.25e-07 2.00e+01 27 U-235 nu-fission 4.31e-02 4.74e-03 +14 6.25e-07 2.00e+01 27 U-238 fission 9.88e-03 1.93e-03 +15 6.25e-07 2.00e+01 27 U-238 nu-fission 2.71e-02 5.21e-03 sum(distribcell) energy low [MeV] energy high [MeV] nuclide score mean std. dev. +0 (0, 100, 2000, 30000) 0.00e+00 6.25e-07 U-235 fission 0.00e+00 0.00e+00 +1 (0, 100, 2000, 30000) 0.00e+00 6.25e-07 U-235 nu-fission 0.00e+00 0.00e+00 +2 (0, 100, 2000, 30000) 0.00e+00 6.25e-07 U-238 fission 0.00e+00 0.00e+00 +3 (0, 100, 2000, 30000) 0.00e+00 6.25e-07 U-238 nu-fission 0.00e+00 0.00e+00 +4 (0, 100, 2000, 30000) 6.25e-07 2.00e+01 U-235 fission 0.00e+00 0.00e+00 +5 (0, 100, 2000, 30000) 6.25e-07 2.00e+01 U-235 nu-fission 0.00e+00 0.00e+00 +6 (0, 100, 2000, 30000) 6.25e-07 2.00e+01 U-238 fission 0.00e+00 0.00e+00 +7 (0, 100, 2000, 30000) 6.25e-07 2.00e+01 U-238 nu-fission 0.00e+00 0.00e+00 +8 (500, 5000, 50000) 0.00e+00 6.25e-07 U-235 fission 0.00e+00 0.00e+00 +9 (500, 5000, 50000) 0.00e+00 6.25e-07 U-235 nu-fission 0.00e+00 0.00e+00 +10 (500, 5000, 50000) 0.00e+00 6.25e-07 U-238 fission 0.00e+00 0.00e+00 +11 (500, 5000, 50000) 0.00e+00 6.25e-07 U-238 nu-fission 0.00e+00 0.00e+00 +12 (500, 5000, 50000) 6.25e-07 2.00e+01 U-235 fission 0.00e+00 0.00e+00 +13 (500, 5000, 50000) 6.25e-07 2.00e+01 U-235 nu-fission 0.00e+00 0.00e+00 +14 (500, 5000, 50000) 6.25e-07 2.00e+01 U-238 fission 0.00e+00 0.00e+00 +15 (500, 5000, 50000) 6.25e-07 2.00e+01 U-238 nu-fission 0.00e+00 0.00e+00 \ No newline at end of file diff --git a/tests/test_tally_slice_merge/test_tally_slice_merge.py b/tests/test_tally_slice_merge/test_tally_slice_merge.py new file mode 100644 index 0000000000..79acf182d6 --- /dev/null +++ b/tests/test_tally_slice_merge/test_tally_slice_merge.py @@ -0,0 +1,165 @@ +#!/usr/bin/env python + +import os +import sys +import glob +import hashlib +import itertools +sys.path.insert(0, os.pardir) +from testing_harness import PyAPITestHarness +import openmc + + +class TallySliceMergeTestHarness(PyAPITestHarness): + def _build_inputs(self): + + # The summary.h5 file needs to be created to read in the tallies + self._input_set.settings.output = {'summary': True} + + # Initialize the tallies file + tallies_file = openmc.TalliesFile() + + # Define nuclides and scores to add to both tallies + self.nuclides = ['U-235', 'U-238'] + self.scores = ['fission', 'nu-fission'] + + # Define filters for energy and spatial domain + + low_energy = openmc.Filter(type='energy', bins=[0., 0.625e-6]) + high_energy = openmc.Filter(type='energy', bins=[0.625e-6, 20.]) + merged_energies = low_energy.merge(high_energy) + + cell_21 = openmc.Filter(type='cell', bins=[21]) + cell_27 = openmc.Filter(type='cell', bins=[27]) + distribcell_filter = openmc.Filter(type='distribcell', bins=[21]) + + self.cell_filters = [cell_21, cell_27] + self.energy_filters = [low_energy, high_energy] + + # Initialize cell tallies with filters, nuclides and scores + tallies = [] + for cell_filter in self.energy_filters: + for energy_filter in self.cell_filters: + for nuclide in self.nuclides: + for score in self.scores: + tally = openmc.Tally() + tally.estimator = 'tracklength' + tally.add_score(score) + tally.add_nuclide(nuclide) + tally.add_filter(cell_filter) + tally.add_filter(energy_filter) + tallies.append(tally) + + # Merge all cell tallies together + while len(tallies) != 1: + halfway = int(len(tallies) / 2) + zip_split = zip(tallies[:halfway], tallies[halfway:]) + tallies = list(map(lambda xy: xy[0].merge(xy[1]), zip_split)) + + # Specify a name for the tally + tallies[0].name = 'cell tally' + + # Initialize a distribcell tally + distribcell_tally = openmc.Tally(name='distribcell tally') + distribcell_tally.estimator = 'tracklength' + distribcell_tally.add_filter(distribcell_filter) + distribcell_tally.add_filter(merged_energies) + for score in self.scores: + distribcell_tally.add_score(score) + for nuclide in self.nuclides: + distribcell_tally.add_nuclide(nuclide) + + # Add tallies to a TalliesFile + tallies_file = openmc.TalliesFile() + tallies_file.add_tally(tallies[0]) + tallies_file.add_tally(distribcell_tally) + + # Export tallies to file + self._input_set.tallies = tallies_file + super(TallySliceMergeTestHarness, self)._build_inputs() + + def _get_results(self, hash_output=False): + """Digest info in the statepoint and return as a string.""" + + # Read the statepoint file. + statepoint = glob.glob(os.path.join(os.getcwd(), self._sp_name))[0] + sp = openmc.StatePoint(statepoint) + + # Read the summary file. + summary = glob.glob(os.path.join(os.getcwd(), 'summary.h5'))[0] + su = openmc.Summary(summary) + sp.link_with_summary(su) + + # Extract the cell tally + tallies = [sp.get_tally(name='cell tally')] + + # Slice the tallies by cell filter bins + cell_filter_prod = itertools.product(tallies, self.cell_filters) + tallies = map(lambda tf: tf[0].get_slice(filters=[tf[1].type], + filter_bins=[tf[1].get_bin(0)]), cell_filter_prod) + + # Slice the tallies by energy filter bins + energy_filter_prod = itertools.product(tallies, self.energy_filters) + tallies = map(lambda tf: tf[0].get_slice(filters=[tf[1].type], + filter_bins=[(tf[1].get_bin(0),)]), energy_filter_prod) + + # Slice the tallies by nuclide + nuclide_prod = itertools.product(tallies, self.nuclides) + tallies = map(lambda tn: tn[0].get_slice(nuclides=[tn[1]]), nuclide_prod) + + # Slice the tallies by score + score_prod = itertools.product(tallies, self.scores) + tallies = map(lambda ts: ts[0].get_slice(scores=[ts[1]]), score_prod) + tallies = list(tallies) + + # Initialize an output string + outstr = '' + + # Append sliced Tally Pandas DataFrames to output string + for tally in tallies: + df = tally.get_pandas_dataframe() + outstr += df.to_string() + + # Merge all tallies together + while len(tallies) != 1: + halfway = int(len(tallies) / 2) + zip_split = zip(tallies[:halfway], tallies[halfway:]) + tallies = list(map(lambda xy: xy[0].merge(xy[1]), zip_split)) + + # Append merged Tally Pandas DataFrame to output string + df = tallies[0].get_pandas_dataframe() + outstr += df.to_string() + + # Extract the distribcell tally + distribcell_tally = sp.get_tally(name='distribcell tally') + + # Sum up a few subdomains from the distribcell tally + sum1 = distribcell_tally.summation(filter_type='distribcell', + filter_bins=[0,100,2000,30000]) + # Sum up a few subdomains from the distribcell tally + sum2 = distribcell_tally.summation(filter_type='distribcell', + filter_bins=[500,5000,50000]) + + # Merge the distribcell tally slices + merge_tally = sum1.merge(sum2) + + # Append merged Tally Pandas DataFrame to output string + df = merge_tally.get_pandas_dataframe() + outstr += df.to_string() + + # Hash the results if necessary + if hash_output: + sha512 = hashlib.sha512() + sha512.update(outstr.encode('utf-8')) + outstr = sha512.hexdigest() + + return outstr + + def _cleanup(self): + super(TallySliceMergeTestHarness, self)._cleanup() + f = os.path.join(os.getcwd(), 'tallies.xml') + if os.path.exists(f): os.remove(f) + +if __name__ == '__main__': + harness = TallySliceMergeTestHarness('statepoint.10.h5', True) + harness.main()