diff --git a/.gitignore b/.gitignore index f6ba21dd38..c92c103401 100644 --- a/.gitignore +++ b/.gitignore @@ -61,5 +61,4 @@ data/nndc # PyCharm project configuration files .idea -.idea/* - +.idea/* \ No newline at end of file diff --git a/src/geometry.F90 b/src/geometry.F90 index 2a67619ad1..0b31d0cf89 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -129,7 +129,6 @@ contains type(Particle), intent(inout) :: p logical, intent(inout) :: found integer, optional :: search_cells(:) - integer :: i ! index over cells integer :: i_xyz(3) ! indices in lattice integer :: n ! number of cells to search @@ -583,7 +582,6 @@ contains type(Particle), intent(inout) :: p integer, intent(in) :: lattice_translation(3) - integer :: i_xyz(3) ! indices in lattice logical :: found ! particle found in cell? class(Lattice), pointer :: lat diff --git a/src/utils/openmc/filter.py b/src/utils/openmc/filter.py index 56463c44db..089f9f18ac 100644 --- a/src/utils/openmc/filter.py +++ b/src/utils/openmc/filter.py @@ -104,7 +104,7 @@ class Filter(object): self._type = type elif not type in FILTER_TYPES.values(): - msg = 'Unable to set Filter type to {0} since it is not one ' \ + msg = 'Unable to set Filter type to "{0}" since it is not one ' \ 'of the supported types'.format(type) raise ValueError(msg) @@ -118,7 +118,7 @@ class Filter(object): self.num_bins = 0 elif self._type is None: - msg = 'Unable to set bins for Filter to {0} since ' \ + msg = 'Unable to set bins for Filter to "{0}" since ' \ 'the Filter type has not yet been set'.format(bins) raise ValueError(msg) @@ -136,12 +136,12 @@ class Filter(object): for edge in bins: if not is_integer(edge): - msg = 'Unable to add bin {0} to a {1} Filter since ' \ + msg = 'Unable to add bin "{0}" to a {1} Filter since ' \ 'it is a non-integer'.format(edge, self._type) raise ValueError(msg) elif edge < 0: - msg = 'Unable to add bin {0} to a {1} Filter since ' \ + msg = 'Unable to add bin "{0}" to a {1} Filter since ' \ 'it is a negative integer'.format(edge, self._type) raise ValueError(msg) @@ -151,22 +151,22 @@ class Filter(object): for edge in bins: if not is_integer(edge) and not is_float(edge): - msg = 'Unable to add bin edge {0} to {1} Filter since ' \ - 'it is a non-integer or floating point ' \ + msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \ + 'since it is a non-integer or floating point ' \ 'value'.format(edge, self._type) raise ValueError(msg) elif edge < 0.: - msg = 'Unable to add bin edge {0} to {1} Filter since it ' \ - 'is a negative value'.format(edge, self._type) + msg = 'Unable to add bin edge "{0}" to a {1} Filter ' \ + 'since it is a negative value'.format(edge, self._type) raise ValueError(msg) # Check that bin edges are monotonically increasing for index in range(len(bins)): if index > 0 and bins[index] < bins[index-1]: - msg = 'Unable to add bin edges {0} to {1} Filter since ' \ - 'they are not monotonically ' \ + msg = 'Unable to add bin edges "{0}" to a {1} Filter ' \ + 'since they are not monotonically ' \ 'increasing'.format(bins, self._type) raise ValueError(msg) @@ -175,17 +175,17 @@ class Filter(object): elif self._type == 'mesh': if not len(bins) == 1: - msg = 'Unable to add bins {0} to a mesh Filter since ' \ + msg = 'Unable to add bins "{0}" to a mesh Filter since ' \ 'only a single mesh can be used per tally'.format(bins) raise ValueError(msg) elif not is_integer(bins[0]): - msg = 'Unable to add bin {0} to mesh Filter since it ' \ + msg = 'Unable to add bin "{0}" to mesh Filter since it ' \ 'is a non-integer'.format(bins[0]) raise ValueError(msg) elif bins[0] < 0: - msg = 'Unable to add bin {0} to mesh Filter since it ' \ + msg = 'Unable to add bin "{0}" to mesh Filter since it ' \ 'is a negative integer'.format(bins[0]) raise ValueError(msg) @@ -198,7 +198,7 @@ class Filter(object): def num_bins(self, num_bins): if not is_integer(num_bins) or num_bins < 0: - msg = 'Unable to set the number of bins {0} for a {1} Filter ' \ + msg = 'Unable to set the number of bins "{0}" for a {1} Filter ' \ 'since it is not a positive ' \ 'integer'.format(num_bins, self._type) raise ValueError(msg) @@ -210,7 +210,7 @@ class Filter(object): def mesh(self, mesh): if not isinstance(mesh, Mesh): - msg = 'Unable to set Mesh to {0} for Filter since it is not a ' \ + msg = 'Unable to set Mesh to "{0}" for Filter since it is not a ' \ 'Mesh object'.format(mesh) raise ValueError(msg) @@ -223,7 +223,7 @@ class Filter(object): def offset(self, offset): if not is_integer(offset): - msg = 'Unable to set offset {0} for a {1} Filter since it is a ' \ + msg = 'Unable to set offset "{0}" for a {1} Filter since it is a ' \ 'non-integer value'.format(offset, self._type) raise ValueError(msg) @@ -234,12 +234,12 @@ class Filter(object): def stride(self, stride): if not is_integer(stride): - msg = 'Unable to set stride {0} for a {1} Filter since it is a ' \ + msg = 'Unable to set stride "{0}" for a {1} Filter since it is a ' \ 'non-integer value'.format(stride, self._type) raise ValueError(msg) if stride < 0: - msg = 'Unable to set stride {0} for a {1} Filter since it is a ' \ + msg = 'Unable to set stride "{0}" for a {1} Filter since it is a ' \ 'negative value'.format(stride, self._type) raise ValueError(msg) @@ -288,17 +288,67 @@ class Filter(object): return merged_filter - def get_bin_index(self, bin): + def get_bin_index(self, filter_bin): + """Returns the index in the Filter for some bin. + + Parameters + ---------- + filter_bin : int, tuple + The bin is the integer ID for 'material', 'surface', 'cell', + 'cellborn', and 'universe' Filters. The bin is an integer for + the cell instance ID for 'distribcell' Filters. The bin is + a 2-tuple of floats for 'energy' and 'energyout' filters + corresponding to the energy boundaries of the bin of interest. + The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding to + the mesh cell of interest. + + Returns + ------- + The index in the Tally data array for this filter bin. + """ + + + # FIXME: This does not work for distribcells!!! try: - index = self._bins.index(bin) + # Filter bins for a mesh are an (x,y,z) tuple + if self.type == 'mesh': + + # Convert (x,y,z) to a single bin -- this is similar to + # subroutine mesh_indices_to_bin in openmc/src/mesh.F90. + if (len(self.mesh.dimension) == 3): + nx, ny, nz = self.mesh.dimension + val = (filter_bin[0] - 1) * ny * nz + \ + (filter_bin[1] - 1) * nz + \ + (filter_bin[2] - 1) + else: + nx, ny = self.mesh.dimension + val = (filter_bin[0] - 1) * ny + \ + (filter_bin[1] - 1) + + filter_index = val + + # Use lower energy bound to find index for energy Filters + elif self.type in ['energy', 'energyout']: + val = self.bins.index(filter_bin[0]) + filter_index = val + + # Filter bins for distribcell are the "IDs" of each unique placement + # of the Cell in the Geometry (integers starting at 0) + elif self._type == 'distribcell': + filter_index = filter_bin + + # Use ID for all other Filters (e.g., material, cell, etc.) + else: + val = self.bins.index(filter_bin) + filter_index = val except ValueError: - msg = 'Unable to get the bin index for Filter since {0} ' \ - 'is not one of the bins'.format(bin) + msg = 'Unable to get the bin index for Filter since "{0}" ' \ + 'is not one of the bins'.format(filter_bin) raise ValueError(msg) - return index + return filter_index def __repr__(self): diff --git a/src/utils/openmc/statepoint.py b/src/utils/openmc/statepoint.py index 636fd9da2f..f5805f0e2c 100644 --- a/src/utils/openmc/statepoint.py +++ b/src/utils/openmc/statepoint.py @@ -332,9 +332,6 @@ class StatePoint(object): subbase = '{0}{1}/filter '.format(base, tally_key) - # Initialize the stride - stride = 1 - # Initialize all Filters for j in range(1, n_filters+1): @@ -359,7 +356,6 @@ class StatePoint(object): bins = self._get_double( n_bins+1, path='{0}{1}/bins'.format(subbase, j)) - # FIXME elif FILTER_TYPES[filter_type] in ['mesh', 'distribcell']: bins = self._get_int( path='{0}{1}/bins'.format(subbase, j))[0] @@ -371,7 +367,6 @@ class StatePoint(object): # Create Filter object filter = openmc.Filter(FILTER_TYPES[filter_type], bins) filter.offset = offset - filter.stride = stride filter.num_bins = n_bins if FILTER_TYPES[filter_type] == 'mesh': @@ -381,9 +376,6 @@ class StatePoint(object): # Add Filter to the Tally tally.add_filter(filter) - # Update the stride for the next Filter - stride *= n_bins - # Read Nuclide bins n_nuclides = self._get_int( path='{0}{1}/n_nuclides'.format(base, tally_key))[0] @@ -406,6 +398,14 @@ class StatePoint(object): n_user_scores = self._get_int( path='{0}{1}/n_user_score_bins'.format(base, tally_key))[0] + # Compute and set the filter strides + for i in range(n_filters): + filter = tally.filters[i] + filter.stride = n_score_bins * n_nuclides + + for j in range(i+1, n_filters): + filter.stride *= tally.filters[j].num_bins + # Read scattering moment order strings (e.g., P3, Y-1,2, etc.) moments = [] subbase = '{0}{1}/moments/'.format(base, tally_key) @@ -435,7 +435,7 @@ class StatePoint(object): tally.add_score(score) # Add Tally to the global dictionary of all Tallies - self._tallies[tally_key] = tally + self.tallies[tally_key] = tally def read_results(self): @@ -574,76 +574,108 @@ class StatePoint(object): # Calculate sample mean and standard deviation for user-defined Tallies - for tally_id, tally in self._tallies.items(): + for tally_id, tally in self.tallies.items(): tally.compute_std_dev(t_value) - def get_tally(self, score, filters, nuclides, - name='', estimator='tracklength'): + def get_tally(self, scores=[], filters=[], nuclides=[], + name=None, id=None, estimator=None): """Finds and returns a Tally object with certain properties. + This routine searches the list of Tallies and returns the first Tally + found it finds which satisfieds all of the input parameters. + NOTE: The input parameters do not need to match the complete Tally + specification and may only represent a subset of the Tallies properties. + Parameters ---------- - score : str - The score string + scores : list + A list of one or more score strings (default is []). filters : list - A list of Filter objects + A list of Filter objects (default is []). nuclides : list - A list of Nuclide objects + A list of Nuclide objects (default is []). name : str - The name specified for the Tally (default is '') + The name specified for the Tally (default is None). + + id : int + The id specified for the Tally (default is None). estimator: str - The type of estimator ('tracklength' (default) or 'analog') + The type of estimator ('tracklength', 'analog'; default is None). + + Returns + ------- + A Tally object. + + Raises + ------ + LookupError : An error when a Tally meeting all of the input + parameters cannot be found in the statepoint. """ - # Loop over the domain-to-tallies mapping to find the Tally tally = None # Iterate over all tallies to find the appropriate one - for tally_id, test_tally in self._tallies.items(): + for tally_id, test_tally in self.tallies.items(): - # Determine if the queried Tally name is the same as this Tally - if not name == test_tally._name: + # Determine if Tally has queried name + if name and name != test_tally.name: continue - # Determine if the queried Tally estimator is the same as this Tally - if not estimator == test_tally._estimator: + # Determine if Tally has queried id + if id and id != test_tally.id: continue - # Determine if the queried Tally scores are the same as this Tally - if not score in test_tally._scores: + # Determine if Tally has queried estimator + if estimator and not estimator == test_tally.estimator: continue - # Determine if queried Tally filters is same length as this Tally - if len(filters) != len(test_tally._filters): - continue + # Determine if Tally has the queried score(s) + if scores: + contains_scores = True - # Determine if the queried Tally filters are the same as this Tally - contains_filters = True + # Iterate over the scores requested by the user + for score in scores: + if not score in test_tally.scores: + contains_scores = False + break - # Iterate over the filters requested by the user - for filter in filters: - if not filter in test_tally._filters: - contains_filters = False - break + if not contains_scores: + continue - # Determine if the queried Nuclide is in this Tally - contains_nuclides = True + # Determine if Tally has the queried Filter(s) + if filters: + contains_filters = True - # Iterate over the Nuclides requested by the user - for nuclide in nuclides: - if not nuclide in test_tally._nuclides: - contains_nuclides = False - break + # Iterate over the Filters requested by the user + for filter in filters: + if not filter in test_tally.filters: + contains_filters = False + break - # If the Tally contained all Filters and Nuclides, return the Tally - if contains_filters and contains_nuclides: - tally = test_tally - break + if not contains_filters: + continue + + # Determine if Tally has the queried Nuclide(s) + if nuclides: + contains_nuclides = True + + # Iterate over the Nuclides requested by the user + for nuclide in nuclides: + if not nuclide in test_tally.nuclides: + contains_nuclides = False + break + + if not contains_nuclides: + continue + + # If the current Tally met user's request, break loop and return it + tally = test_tally + break # If we did not find the Tally, return an error message if tally is None: @@ -652,41 +684,37 @@ class StatePoint(object): return tally - def get_tally_id(self, score, filters, name='', estimator='tracklength'): - """Retrieve the Tally ID for a given list of filters and score(s). + def link_with_summary(self, summary): + """Links Tallies and Filters with Summary model information. + + This routine retrieves model information (materials, geometry) from a + Summary object populated with an HDF5 'summary.h5' file and inserts + it into the Tally objects. This can be helpful when viewing and + manipulating large scale Tally data. Parameters ---------- - score : str - The score string + summary : Summary + A Summary object. - filters : list - A list of Filter objects - - name : str - The name specified for the Tally (default is '') - - estimator: str - The type of estimator ('tracklength' (default) or 'analog') + Raises + ------ + ValueError : An error when the argument passed to the 'summary' + parameter is not an openmc.Summary object. """ - tally = self.get_tally(score, filters, name, estimator) - return tally._id - - - def link_with_summary(self, summary): - if not isinstance(summary, openmc.summary.Summary): msg = 'Unable to link statepoint with {0} which ' \ 'is not a Summary object'.format(summary) raise ValueError(msg) - for tally_id, tally in self._tallies.items(): + for tally_id, tally in self.tallies.items(): # Get the Tally name from the summary file tally.name = summary.tallies[tally_id].name + tally.with_summary = True - nuclide_zaids = copy.deepcopy(tally._nuclides) + nuclide_zaids = copy.deepcopy(tally.nuclides) for nuclide_zaid in nuclide_zaids: @@ -696,32 +724,32 @@ class StatePoint(object): else: tally.add_nuclide(summary.nuclides[nuclide_zaid]) - for filter in tally._filters: + for filter in tally.filters: - if filter._type == 'surface': + if filter.type == 'surface': surface_ids = [] - for bin in filter._bins: - surface_ids.append(summary.surfaces[bin]._id) + for bin in filter.bins: + surface_ids.append(summary.surfaces[bin].id) filter.bins = surface_ids - if filter._type in ['cell', 'distribcell']: + if filter.type in ['cell', 'distribcell']: distribcell_ids = [] - for bin in filter._bins: - distribcell_ids.append(summary.cells[bin]._id) + for bin in filter.bins: + distribcell_ids.append(summary.cells[bin].id) filter.bins = distribcell_ids - if filter._type == 'universe': + if filter.type == 'universe': universe_ids = [] - for bin in filter._bins: - universe_ids.append(summary.universes[bin]._id) + for bin in filter.bins: + universe_ids.append(summary.universes[bin].id) filter.bins = universe_ids - if filter._type == 'material': + if filter.type == 'material': material_ids = [] - for bin in filter._bins: - material_ids.append(summary.materials[bin]._id) + for bin in filter.bins: + material_ids.append(summary.materials[bin].id) filter.bins = material_ids - + self._with_summary = True diff --git a/src/utils/openmc/summary.py b/src/utils/openmc/summary.py index fec616c064..61fb89a1ed 100644 --- a/src/utils/openmc/summary.py +++ b/src/utils/openmc/summary.py @@ -491,6 +491,10 @@ class Summary(object): self.tallies = {} # Read the number of tallies + if not 'tallies' in self._f.keys(): + self.n_tallies = 0 + return + self.n_tallies = self._f['tallies/n_tallies'][0] # OpenMC Tally keys @@ -532,15 +536,16 @@ class Summary(object): subsubbase = '{0}/filter {1}'.format(subbase, j) - # Read filter type (e.g., "cell", "energy", etc.) integer code - filter_type = self._f['{0}/type'.format(subsubbase)][0] + # Read filter type (e.g., "cell", "energy", etc.) + filter_type_code = self._f['{0}/type'.format(subsubbase)][0] + filter_type = openmc.FILTER_TYPES[filter_type_code] # Read the filter bins num_bins = self._f['{0}/n_bins'.format(subsubbase)][0] bins = self._f['{0}/bins'.format(subsubbase)][...] # Create Filter object - filter = openmc.Filter(openmc.FILTER_TYPES[filter_type], bins) + filter = openmc.Filter(filter_type, bins) filter.num_bins = num_bins # Add Filter to the Tally diff --git a/src/utils/openmc/tallies.py b/src/utils/openmc/tallies.py index 2971a85373..59dd03dde4 100644 --- a/src/utils/openmc/tallies.py +++ b/src/utils/openmc/tallies.py @@ -1,8 +1,12 @@ import copy import os +import itertools from xml.etree import ElementTree as ET +import numpy as np + from openmc import Mesh, Filter, Trigger, Nuclide +from openmc.summary import Summary from openmc.clean_xml import * from openmc.checkvalue import * @@ -31,6 +35,7 @@ class Tally(object): self._num_score_bins = 0 self._num_realizations = 0 + self._with_summary = False self._sum = None self._sum_sq = None @@ -46,30 +51,30 @@ class Tally(object): if existing is None: clone = type(self).__new__(type(self)) - clone._id = self._id - clone._name = self._name - clone._estimator = self._estimator - clone._num_score_bins = self._num_score_bins - clone._num_realizations = self._num_realizations - clone._sum = copy.deepcopy(self._sum, memo) - clone._sum_sq = copy.deepcopy(self._sum_sq, memo) - clone._mean = copy.deepcopy(self._mean, memo) - clone._std_dev = copy.deepcopy(self._std_dev, memo) + clone.id = self.id + clone.name = self.name + clone.estimator = self.estimator + clone.num_score_bins = self.num_score_bins + clone.num_realizations = self.num_realizations + clone._sum = copy.deepcopy(self.sum, memo) + clone._sum_sq = copy.deepcopy(self.sum_sq, memo) + clone._mean = copy.deepcopy(self.mean, memo) + clone._std_dev = copy.deepcopy(self.std_dev, memo) - clone._filters = [] - for filter in self._filters: + clone.filters = [] + for filter in self.filters: clone.add_filter(copy.deepcopy(filter, memo)) - clone._nuclides = [] - for nuclide in self._nuclides: + clone.nuclides = [] + for nuclide in self.nuclides: clone.add_nuclide(copy.deepcopy(nuclide, memo)) - clone._scores = [] - for score in self._scores: + clone.scores = [] + for score in self.scores: clone.add_score(score) - clone._triggers = [] - for trigger in self._triggers: + clone.triggers = [] + for trigger in self.triggers: clone.add_trigger(trigger) memo[id(self)] = clone @@ -84,21 +89,30 @@ class Tally(object): def __eq__(self, tally2): # Check all filters - for filter in self._filters: - if not filter in tally2._filters: + if len(self.filters) != len(tally2.filters): + return False + + for filter in self.filters: + if not filter in tally2.filters: return False # Check all nuclides - for nuclide in self._nuclides: - if not nuclide in tally2._nuclides: + if len(self.nuclides) != len(tally2.nuclides): + return False + + for nuclide in self.nuclides: + if not nuclide in tally2.nuclides: return False # Check all scores - for score in self._scores: - if not score in tally2._scores: + if len(self.scores) != len(tally2.scores): + return False + + for score in self.scores: + if not score in tally2.scores: return False - if self._estimator != tally2._estimator: + if self.estimator != tally2.estimator: return False return True @@ -107,17 +121,17 @@ class Tally(object): def __hash__(self): hashable = [] - for filter in self._filters: - hashable.append((filter._type, tuple(filter._bins))) + for filter in self.filters: + hashable.append((filter.type, tuple(filter.bins))) - for nuclide in self._nuclides: - hashable.append(nuclide._name) + for nuclide in self.nuclides: + hashable.append(nuclide.name) - for score in self._scores: + for score in self.scores: hashable.append(score) - hashable.append(self._estimator) - hashable.append(self._name) + hashable.append(self.estimator) + hashable.append(self.name) return hash(tuple(hashable)) @@ -129,7 +143,7 @@ class Tally(object): new_tally = Tally() new_tally._mean = self._mean + other._mean - new_tally._std_dev = np.sqrt(self._std_dev**2 + other._std_dev**2) + new_tally._std_dev = np.sqrt(self.std_dev**2 + other.std_dev**2) @property @@ -211,6 +225,11 @@ class Tally(object): return self._num_realizations + @property + def with_summary(self): + return self._with_summary + + @property def sum(self): return self._sum @@ -236,7 +255,7 @@ class Tally(object): if not estimator in ['analog', 'tracklength']: msg = 'Unable to set the estimator for Tally ID={0} to {1} since ' \ - 'it is not a valid estimator type'.format(self._id, estimator) + 'it is not a valid estimator type'.format(self.id, estimator) raise ValueError(msg) self._estimator = estimator @@ -245,8 +264,8 @@ class Tally(object): def add_trigger(self, trigger): if not isinstance(trigger, Trigger): - msg = 'Unable to add a tally trigger for for Tally ID={0} to ' \ - '{1} since it is not a valid estimator type'.format(trigger) + msg = 'Unable to add a tally trigger for Tally ID={0} to ' \ + 'since "{1}" is not a Trigger'.format(self.id, trigger) raise ValueError(msg) self._triggers.append(trigger) @@ -279,7 +298,7 @@ class Tally(object): if not is_string(name): msg = 'Unable to set name for Tally ID={0} with a non-string ' \ - 'value {1}'.format(self._id, name) + 'value "{1}"'.format(self.id, name) raise ValueError(msg) else: @@ -291,8 +310,8 @@ class Tally(object): global filters if not isinstance(filter, Filter): - msg = 'Unable to add Filter {0} to Tally ID={1} since it is not ' \ - 'a Filter object'.format(filter, self._id) + msg = 'Unable to add Filter "{0}" to Tally ID={1} since it is ' \ + 'not a Filter object'.format(filter, self.id) raise ValueError(msg) self._filters.append(filter) @@ -305,12 +324,12 @@ class Tally(object): def add_score(self, score): if not is_string(score): - msg = 'Unable to add score {0} to Tally ID={1} since it is not a ' \ - 'string'.format(score, self._id) + msg = 'Unable to add score "{0}" to Tally ID={1} since it is ' \ + 'not a string'.format(score, self.id) raise ValueError(msg) # If the score is already in the Tally, don't add it again - if score in self._scores: + if score in self.scores: return else: self._scores.append(score) @@ -325,13 +344,13 @@ class Tally(object): def num_realizations(self, num_realizations): if not is_integer(num_realizations): - msg = 'Unable to set the number of realizations to {0} for ' \ + msg = 'Unable to set the number of realizations to "{0}" for ' \ 'Tally ID={1} since it is not an ' \ 'integer'.format(num_realizations) raise ValueError(msg) elif num_realizations < 0: - msg = 'Unable to set the number of realizations to {0} for ' \ + msg = 'Unable to set the number of realizations to "{0}" for ' \ 'Tally ID={1} since it is a negative ' \ 'value'.format(num_realizations) raise ValueError(msg) @@ -339,18 +358,29 @@ class Tally(object): self._num_realizations = num_realizations + @with_summary.setter + def with_summary(self, with_summary): + + if not isinstance(with_summary, bool): + msg = 'Unable to set with_summary to a non-boolean ' \ + 'value "{0}"'.format(with_summary) + raise ValueError(msg) + + self._with_summary = with_summary + + def set_results(self, sum, sum_sq): if not isinstance(sum, (tuple, list, np.ndarray)): - msg = 'Unable to set the sum to {0}for Tally ID={1} since ' \ + msg = 'Unable to set the sum to "{0}" for Tally ID={1} since ' \ 'it is not a Python tuple/list or NumPy ' \ - 'array'.format(sum, self._id) + 'array'.format(sum, self.id) raise ValueError(msg) if not isinstance(sum_sq, (tuple, list, np.ndarray)): - msg = 'Unable to set the sum to {0}for Tally ID={1} since ' \ + msg = 'Unable to set the sum to "{0}" for Tally ID={1} since ' \ 'it is not a Python tuple/list or NumPy ' \ - 'array'.format(sum_sq, self._id) + 'array'.format(sum_sq, self.id) raise ValueError(msg) self._sum = sum @@ -359,9 +389,9 @@ class Tally(object): def remove_score(self, score): - if not score in self._scores: - msg = 'Unable to remove score {0} from Tally ID={1} since the ' \ - 'Tally does not contain this score'.format(score, self._id) + if not score in self.scores: + msg = 'Unable to remove score "{0}" from Tally ID={1} since the ' \ + 'Tally does not contain this score'.format(score, self.id) ValueError(msg) self._scores.remove(score) @@ -369,9 +399,9 @@ class Tally(object): def remove_filter(self, filter): - if not filter in self._filters: - msg = 'Unable to remove filter {0} from Tally ID={1} since the ' \ - 'Tally does not contain this filter'.format(filter, self._id) + if not filter in self.filters: + msg = 'Unable to remove filter "{0}" from Tally ID={1} since the ' \ + 'Tally does not contain this filter'.format(filter, self.id) ValueError(msg) self._filters.remove(filter) @@ -379,9 +409,9 @@ class Tally(object): def remove_nuclide(self, nuclide): - if not nuclide in self._nuclides: - msg = 'Unable to remove nuclide {0} from Tally ID={1} since the ' \ - 'Tally does not contain this nuclide'.format(nuclide, self._id) + if not nuclide in self.nuclides: + msg = 'Unable to remove nuclide "{0}" from Tally ID={1} since the ' \ + 'Tally does not contain this nuclide'.format(nuclide, self.id) ValueError(msg) self._nuclides.remove(nuclide) @@ -390,36 +420,36 @@ class Tally(object): def compute_std_dev(self, t_value=1.0): # Calculate sample mean and standard deviation - self._mean = self._sum / self._num_realizations - self._std_dev = np.sqrt((self._sum_sq / self._num_realizations - \ - self._mean**2) / (self._num_realizations - 1)) + self._mean = self.sum / self.num_realizations + self._std_dev = np.sqrt((self.sum_sq / self.num_realizations - \ + self.mean**2) / (self.num_realizations - 1)) self._std_dev *= t_value def __repr__(self): string = 'Tally\n' - string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self._id) - string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self._name) + string += '{0: <16}{1}{2}\n'.format('\tID', '=\t', self.id) + string += '{0: <16}{1}{2}\n'.format('\tName', '=\t', self.name) string += '{0: <16}\n'.format('\tFilters') - for filter in self._filters: - string += '{0: <16}\t\t{1}\t{2}\n'.format('', filter._type, - filter._bins) + for filter in self.filters: + string += '{0: <16}\t\t{1}\t{2}\n'.format('', filter.type, + filter.bins) string += '{0: <16}{1}'.format('\tNuclides', '=\t') - for nuclide in self._nuclides: + for nuclide in self.nuclides: if isinstance(nuclide, Nuclide): - string += '{0} '.format(nuclide._name) + string += '{0} '.format(nuclide.name) else: string += '{0} '.format(nuclide) string += '\n' - string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self._scores) - string += '{0: <16}{1}{2}\n'.format('\tEstimator', '=\t', self._estimator) + string += '{0: <16}{1}{2}\n'.format('\tScores', '=\t', self.scores) + string += '{0: <16}{1}{2}\n'.format('\tEstimator', '=\t', self.estimator) return string @@ -430,26 +460,26 @@ class Tally(object): return False # Must have same estimator - if self._estimator != tally._estimator: + if self.estimator != tally.estimator: return False # Must have same nuclides - if len(self._nuclides) != len(tally._nuclides): + if len(self.nuclides) != len(tally.nuclides): return False - for nuclide in self._nuclides: - if not nuclide in tally._nuclides: + for nuclide in self.nuclides: + if not nuclide in tally.nuclides: return False # Must have same or mergeable filters - if len(self._filters) != len(tally._filters): + if len(self.filters) != len(tally.filters): return False # Look to see if all filters are the same, or one or more can be merged - for filter1 in self._filters: + for filter1 in self.filters: mergeable_filter = False - for filter2 in tally._filters: + for filter2 in tally.filters: if filter1 == filter2 or filter1.can_merge(filter2): mergeable_filter = True break @@ -475,19 +505,19 @@ class Tally(object): merged_tally.id = None # Merge filters - for i, filter1 in enumerate(merged_tally._filters): - for filter2 in tally._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 + merged_tally.filters[i] = merged_filter break # Add scores from second tally to merged tally - for score in tally._scores: + for score in tally.scores: merged_tally.add_score(score) # Add triggers from second tally to merged tally - for trigger in tally._triggers: + for trigger in tally.triggers: merged_tally.add_trigger(trigger) return merged_tally @@ -498,33 +528,33 @@ class Tally(object): element = ET.Element("tally") # Tally ID - element.set("id", str(self._id)) + element.set("id", str(self.id)) # Optional Tally name - if self._name != '': - element.set("name", self._name) + if self.name != '': + element.set("name", self.name) # Optional Tally filters - for filter in self._filters: + for filter in self.filters: subelement = ET.SubElement(element, "filter") - subelement.set("type", str(filter._type)) + subelement.set("type", str(filter.type)) - if not filter._bins is None: + if not filter.bins is None: bins = '' - for bin in filter._bins: + for bin in filter.bins: bins += '{0} '.format(bin) subelement.set("bins", bins.rstrip(' ')) # Optional Nuclides - if len(self._nuclides) > 0: + if len(self.nuclides) > 0: nuclides = '' - for nuclide in self._nuclides: + for nuclide in self.nuclides: if isinstance(nuclide, Nuclide): - nuclides += '{0} '.format(nuclide._name) + nuclides += '{0} '.format(nuclide.name) else: nuclides += '{0} '.format(nuclide) @@ -532,198 +562,677 @@ class Tally(object): subelement.text = nuclides.rstrip(' ') # Scores - if len(self._scores) == 0: + if len(self.scores) == 0: msg = 'Unable to get XML for Tally ID={0} since it does not ' \ - 'contain any scores'.format(self._id) + 'contain any scores'.format(self.id) raise ValueError(msg) else: scores = '' - for score in self._scores: + for score in self.scores: scores += '{0} '.format(score) subelement = ET.SubElement(element, "scores") subelement.text = scores.rstrip(' ') # Tally estimator type - if not self._estimator is None: + if not self.estimator is None: subelement = ET.SubElement(element, "estimator") - subelement.text = self._estimator + subelement.text = self.estimator # Optional Triggers - for trigger in self._triggers: + for trigger in self.triggers: trigger.get_trigger_xml(element) return element - def find_filter(self, filter_type, bins): + def find_filter(self, filter_type): filter = None - for test_filter in self._filters: - - # Determine if the Filter has the same type as the one requested - if test_filter._type != filter_type: - continue - - # Determine if the Filter has the same bin edges as the one requested - elif test_filter._bins != bins: - continue - - else: + # 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 = test_filter break - # If we found the Filter, return it - if not filter is None: - return filter - - # Otherwise, throw an Exception - else: - msg = 'Unable to find filter type {0} with bin edges {1} in ' \ - 'Tally ID={2}'.format(filter_type, bins, self._id) + # If we did not find the Filter, throw an Exception + if filter is None: + msg = 'Unable to find filter type "{0}" in ' \ + 'Tally ID={1}'.format(filter_type, self.id) raise ValueError(msg) + return filter + + + def get_filter_index(self, filter_type, filter_bin): + """Returns the index in the Tally's results array for a Filter bin + + Parameters + ---------- + filter_type : str + The type of Filter (e.g., 'cell', 'energy', etc.) + + filter_bin : int, list + The bin is an integer ID for 'material', 'surface', 'cell', + 'cellborn', and 'universe' Filters. The bin is an integer for + the cell instance ID for 'distribcell' Filters. The bin is + a 2-tuple of floats for 'energy' and 'energyout' filters + corresponding to the energy boundaries of the bin of interest. + The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding to + the mesh cell of interest. + + Returns + ------- + The index in the Tally data array for this filter bin. + """ + + # Find the equivalent Filter in this Tally's list of Filters + filter = self.find_filter(filter_type) + + # Get the index for the requested bin from the Filter and return it + filter_index = filter.get_bin_index(filter_bin) + return filter_index + + + def get_nuclide_index(self, nuclide): + """Returns the index in the Tally's results array for a Nuclide bin + + Parameters + ---------- + nuclide : str + The name of the Nuclide (e.g., 'H-1', 'U-238') + + Returns + ------- + The index in the Tally data array for this nuclide. + + Raises + ------ + KeyError : An error when the argument passed to the 'nuclide' + parameter cannot be found in the Tally. + """ + + nuclide_index = -1 + + # Look for the user-requested nuclide in all of the Tally's Nuclides + for i, test_nuclide in enumerate(self.nuclides): + + # If the Summary was linked, then values are Nuclide objects + if isinstance(test_nuclide, Nuclide): + if test_nuclide._name == nuclide: + nuclide_index = i + break + + # If the Summary has not been linked, then values are ZAIDs + else: + if test_nuclide == nuclide: + nuclide_index = i + break + + if nuclide_index == -1: + msg = 'Unable to get the nuclide index for Tally since "{0}" ' \ + 'is not one of the nuclides'.format(nuclide) + raise KeyError(msg) + else: + return nuclide_index + def get_score_index(self, score): - - try: - index = self._scores.index(score) - - except ValueError: - msg = 'Unable to get the score index for Tally since {0} ' \ - 'is not one of the bins'.format(bin) - raise ValueError(msg) - - return index - - - def get_value(self, score, filters, filter_bins, nuclide=None, value='mean'): - """Returns a tally score value given a list of filters to satisfy. + """Returns the index in the Tally's results array for a score bin Parameters ---------- score : str - The score string of interest + The score string (e.g., 'absorption', 'nu-fission') - filters : list - A list of the filters of interest + Returns + ------- + The index in the Tally data array for this score. - filter_bins : list - A list of the filter bins of interest. These are integers for - material, surface, cell, cellborn, distribcell, universe filters, - and floats for energy or energyout filters. The bins are tuples - of three integers (x,y,z) for mesh filters. The order of the bins - in the list is assumed to correspond to the order of the filters. - - nuclide : Nuclide - The Nuclide of interest - - value : str - A string for the type of value to return ('mean' (default), 'std_dev', - 'sum', or 'sum_sq' are accepted) + Raises + ------ + ValueError: An error when the argument passed to the 'score' + parameter cannot be found in the Tally. """ - # Determine the score index from the score string - score_index = self._scores.index(score) + try: + score_index = self.scores.index(score) - # Determine the nuclide index from the nuclide string/object - if not nuclide is None: - nuclide_index = self._nuclides.index(nuclide) + except ValueError: + msg = 'Unable to get the score index for Tally since "{0}" ' \ + 'is not one of the scores'.format(score) + raise ValueError(msg) + + return score_index + + + def get_values(self, scores=[], filters=[], filter_bins=[], + nuclides=[], value='mean'): + """Returns a tally score value given a list of filters to satisfy. + + This routine constructs a 3D NumPy array for the requested Tally data + indexed by filter bin, nuclide bin, and score index. The routine will + order the data in the array + + Parameters + ---------- + scores : list + A list of one or more score strings + (e.g., ['absorption', 'nu-fission']; default is []) + + filters : list + A list of filter type strings + (e.g., ['mesh', 'energy']; default is []) + + filter_bins : list + A list of the filter bins corresponding to the filter_types + parameter (e.g., [1, (0., 0.625e-6)]; default is []). 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. + The bin is a (x,y,z) 3-tuple for 'mesh' filters corresponding + to the mesh cell of interest. The order of the bins in the list + must correspond of the filter_types parameter. + + nuclides : list + A list of nuclide name strings + (e.g., ['U-235', 'U-238']; default is []) + + value : str + A string for the type of value to return - 'mean' (default), + 'std_dev', 'rel_err', 'sum', or 'sum_sq' are accepted + + Returns + ------- + A scalar or NumPy array of the Tally data indexed in the order + each filter, nuclide and score is listed in the parameters. + + Raises + ------ + ValueError : An error when this routine is called before the Tally + is populated with data by the StatePoint.read_results() routine. + """ + + # Ensure that StatePoint.read_results() was called first + if self._sum is None or self._sum_sq is None: + msg = 'The Tally ID={0} has no data to return. Call the ' \ + 'StatePoint.read_results() routine before using ' \ + 'Tally.get_values(...)'.format(self.id) + raise ValueError(msg) + + + # Compute batch statistics if not yet computed + self.compute_std_dev() + + ############################ FILTERS ######################### + # Determine the score indices from any of the requested scores + if filters: + + # Initialize empty list of indices for each bin in each Filter + filter_indices = [] + + # Loop over all of the Tally's Filters + for i, filter in enumerate(self.filters): + + # Initialize empty list of indices for this Filter's bins + filter_indices.append([]) + + user_filter = False + + # If a user-requested Filter, get the user-requested bins + for j, test_filter in enumerate(filters): + if filter.type == test_filter: + bins = filter_bins[j] + user_filter = True + break + + # If not a user-requested Filter, get all bins + if not user_filter: + + # Create list of 2- or 3-tuples tuples for mesh cell bins + if filter.type == 'mesh': + dimension = filter.mesh.dimension + xyz = map(lambda x: np.arange(1,x+1), dimension) + bins = list(itertools.product(*xyz)) + + # Create list of 2-tuples for energy boundary bins + elif filter.type in ['energy', 'energyout']: + bins = [] + for i in range(filter.num_bins): + bins.append((filter.bins[i], filter.bins[i+1])) + + # Create list of IDs for bins for all other Filter types + else: + bins = filter.bins + + # Add indices for each bin in this Filter to the list + for bin in bins: + filter_indices[i].append( + self.get_filter_index(filter.type, bin)) + + # Apply cross-product sum between all filter bin indices + filter_indices = map(sum, itertools.product(*filter_indices)) + + # If user did not specify any specific Filters, use them all else: - nuclide_index = 0 + filter_indices = np.arange(self.num_filter_bins) - # Initialize index for Filter in Tally.results[:,:,:] - filter_index = 0 + ############################ NUCLIDES ######################## + # Determine the score indices from any of the requested scores + if nuclides: + nuclide_indices = np.zeros(len(nuclides), dtype=np.int) + for i, nuclide in enumerate(nuclides): + nuclide_indices[i] = self.get_nuclide_index(nuclide) - # Iterate over specified Filters to compute filter index - for i, filter in enumerate(filters): + # If user did not specify any specific Nuclides, use them all + else: + nuclide_indices = np.arange(self.num_nuclides) - # Find the equivalent Filter in this Tally's list of Filters - test_filter = self.find_filter(filter._type, filter._bins) + ############################# SCORES ######################### + # Determine the score indices from any of the requested scores + if scores: + score_indices = np.zeros(len(scores), dtype=np.int) + for i, score in enumerate(scores): + score_indices[i] = self.get_score_index(score) - # Filter bins for a mesh are an (x,y,z) tuple - if filter._type == 'mesh': + # If user did not specify any specific scores, use them all + else: + score_indices = np.arange(self.num_scores) - # Get the dimensions of the corresponding mesh - nx, ny, nz = test_filter._mesh._dimension - - # Convert (x,y,z) to a single bin -- this is similar to - # subroutine mesh_indices_to_bin in openmc/src/mesh.F90. - val = ((filter_bins[i][0] - 1) * ny * nz + - (filter_bins[i][1] - 1) * nz + - (filter_bins[i][2] - 1)) - filter_index += val * test_filter._stride - - # Filter bins for distribcell are the "IDs" of each unique placement - # of the Cell in the Geometry (integers starting at 0) - elif filter._type == 'distribcell': - bin = filter_bins[i] - filter_index += bin * test_filter._stride - - else: - bin = filter_bins[i] - bin_index = test_filter.get_bin_index(bin) - filter_index += bin_index * test_filter._stride + # Construct cross-product of all three index types with each other + indices = np.ix_(filter_indices, nuclide_indices, score_indices) # Return the desired result from Tally if value == 'mean': - return self._mean[filter_index, nuclide_index, score_index] + data = self.mean[indices] elif value == 'std_dev': - return self._std_dev[filter_index, nuclide_index, score_index] + data = self.std_dev[indices] + elif value == 'rel_err': + data = self.std_dev[indices] / self.mean[indices] elif value == 'sum': - return self._sum[filter_index, nuclide_index, score_index] + data = self.sum[indices] elif value == 'sum_sq': - return self._sum_sq[filter_index, nuclide_index, score_index] + data = self.sum_sq[indices] else: - msg = 'Unable to return results from Tally ID={0} for score {1} ' \ - 'since the value {2} is not \'mean\', \'std_dev\', ' \ - '\'sum\', or \'sum_sq\''.format(self._id, score, value) + msg = 'Unable to return results from Tally ID={0} since the ' \ + 'the requested value "{1}" is not \'mean\', \'std_dev\', ' \ + '\rel_err\', \'sum\', or \'sum_sq\''.format(self.id, value) raise LookupError(msg) + return data.squeeze() + + + def get_pandas_dataframe(self, filters=True, nuclides=True, + scores=True, summary=None): + """Build a Pandas DataFrame for the Tally data. + + This routine constructs a Pandas DataFrame object for the Tally data + with columns annotated by filter, nuclide and score bin information. + This capability has been tested for Pandas >=v0.13.1. However, if p + possible, it is recommended to use the v0.16 or newer versions of + Pandas since this this routine uses the Multi-index Pandas feature. + + Parameters + ---------- + filters : bool + Include columns with filter bin information (default is True). + + nuclides : bool + Include columns with nuclide bin information (default is True). + + scores : bool + Include columns with score bin information (default is True). + + summary : None or Summary + An optional Summary object to be used to construct columns for + for distribcell tally filters (default is None). The geometric + information in the Summary object is embedded into a Multi-index + column with a geometric "path" to each distribcell intance. + NOTE: This option requires the OpenCG Python package. + + Returns + ------- + A Pandas DataFrame with each column annotated by filter, nuclide + and score bin information (if these parameters are True), and the + mean and standard deviation of the Tally's data. + + Raises + ------ + KeyError : An error when this routine is called before the Tally + is populated with data by the StatePoint.read_results() routine. + """ + + # Ensure that StatePoint.read_results() was called first + if self._sum is None or self._sum_sq is None: + msg = 'The Tally ID={0} has no data to return. Call the ' \ + 'StatePoint.read_results() routine before using ' \ + 'Tally.get_pandas_dataframe(...)'.format(self.id) + raise KeyError(msg) + + # If using Summary, ensure StatePoint.link_with_summary(...) was called + if summary and not self.with_summary: + msg = 'The Tally ID={0} has not been linked with the Summary. ' \ + 'Call the StatePoint.link_with_summary(...) routine ' \ + 'before using Tally.get_pandas_dataframe(...) with ' \ + 'Summary info'.format(self.id) + raise KeyError(msg) + + # Attempt to import the pandas package + try: + import pandas as pd + except ImportError: + msg = 'The pandas Python package must be installed on your system' + raise ImportError(msg) + + # Compute batch statistics if not yet computed + self.compute_std_dev() + + # Initialize a pandas dataframe for the tally data + df = pd.DataFrame() + + # Find the total length of the tally data array + data_size = self.sum.size + + # Build DataFrame columns for filters if user requested them + if filters: + + for filter in self.filters: + + # mesh filters + if filter.type == 'mesh': + + # Initialize dictionary to build Pandas Multi-index column + filter_dict = {} + + # Append Mesh ID as outermost index of mult-index + mesh_id = filter.mesh.id + mesh_key = 'mesh {0}'.format(mesh_id) + + # Find mesh dimensions - use 3D indices for simplicity + if (len(filter.mesh.dimension) == 3): + nx, ny, nz = filter.mesh.dimension + else: + nx, ny = filter.mesh.dimension + nz = 1 + + # Generate multi-index sub-column for x-axis + filter_bins = np.arange(1, nx+1) + repeat_factor = ny * nz * filter.stride + filter_bins = np.repeat(filter_bins, repeat_factor) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + filter_dict[(mesh_key, 'x')] = filter_bins + + # Generate multi-index sub-column for y-axis + filter_bins = np.arange(1, ny+1) + repeat_factor = nz * filter.stride + filter_bins = np.repeat(filter_bins, repeat_factor) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + filter_dict[(mesh_key, 'y')] = filter_bins + + # Generate multi-index sub-column for z-axis + filter_bins = np.arange(1, nz+1) + repeat_factor = filter.stride + filter_bins = np.repeat(filter_bins, repeat_factor) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + filter_dict[(mesh_key, 'z')] = filter_bins + + # Append the multi-index column to the DataFrame + df = pd.concat([df, pd.DataFrame(filter_dict)], axis=1) + + # distribcell filters + elif filter.type == 'distribcell': + + if isinstance(summary, Summary): + + # Attempt to import the OpenCG package + try: + import opencg + except ImportError: + msg = 'The OpenCG package must be installed ' \ + 'to use a Summary for distribcell dataframes' + raise ImportError(msg) + + # Create and extract the OpenCG geometry the Summary + summary.make_opencg_geometry() + opencg_geometry = summary.opencg_geometry + openmc_geometry = summary.openmc_geometry + + # Use OpenCG to compute the number of regions + opencg_geometry.initializeCellOffsets() + num_regions = opencg_geometry._num_regions + + # Initialize a dictionary mapping OpenMC distribcell + # offsets to OpenCG LocalCoords linked lists + offsets_to_coords = {} + + # Use OpenCG to compute LocalCoords linked list for + # each region and store in dictionary + for region in range(num_regions): + coords = opencg_geometry.findRegion(region) + path = opencg.get_path(coords) + cell_id = path[-1] + + # If this region is in Cell corresponding to the + # distribcell filter bin, store it in dictionary + if cell_id == filter.bins[0]: + offset = openmc_geometry.get_offset(path, + filter.offset) + offsets_to_coords[offset] = coords + + # Each distribcell offset is a DataFrame bin + # Unravel the paths into DataFrame columns + num_offsets = len(offsets_to_coords) + + # Initialize termination condition for while loop + levels_remain = True + counter = 0 + + # Iterate over each level in the CSG tree hierarchy + while levels_remain: + levels_remain = False + + # Initialize dictionary to build Pandas Multi-index + # column for this level in the CSG tree hierarchy + level_dict = {} + + # Initialize prefix Multi-index keys + counter += 1 + level_key = 'level {0}'.format(counter) + univ_key = (level_key, 'univ', 'id') + cell_key = (level_key, 'cell', 'id') + lat_id_key = (level_key, 'lat', 'id') + lat_x_key = (level_key, 'lat', 'x') + lat_y_key = (level_key, 'lat', 'y') + lat_z_key = (level_key, 'lat', 'z') + + # Allocate NumPy arrays for each CSG level and + # each Multi-index column in the DataFrame + level_dict[univ_key] = np.empty(num_offsets) + level_dict[cell_key] = np.empty(num_offsets) + level_dict[lat_id_key] = np.empty(num_offsets) + level_dict[lat_x_key] = np.empty(num_offsets) + level_dict[lat_y_key] = np.empty(num_offsets) + level_dict[lat_z_key] = np.empty(num_offsets) + + # Initialize Multi-index columns to NaN - this is + # necessary since some distribcell instances may + # have very different LocalCoords linked lists + level_dict[univ_key][:] = np.nan + level_dict[cell_key][:] = np.nan + level_dict[lat_id_key][:] = np.nan + level_dict[lat_x_key][:] = np.nan + level_dict[lat_y_key][:] = np.nan + level_dict[lat_z_key][:] = np.nan + + # Iterate over all regions (distribcell instances) + for offset in range(num_offsets): + coords = offsets_to_coords[offset] + + # If entire LocalCoords has been unraveled into + # Multi-index columns already, continue + if coords == None: + continue + + # Assign entry to Universe Multi-index column + if coords._type == 'universe': + univ_id = coords._universe._id + cell_id = coords._cell._id + level_dict[univ_key][offset] = univ_id + level_dict[cell_key][offset] = cell_id + + # Assign entry to Lattice Multi-index column + else: + lat_id = coords._lattice._id + lat_x = coords._lat_x + lat_y = coords._lat_y + lat_z = coords._lat_z + level_dict[lat_id_key][offset] = lat_id + level_dict[lat_x_key][offset] = lat_x + level_dict[lat_y_key][offset] = lat_y + level_dict[lat_z_key][offset] = lat_z + + # Move to next node in LocalCoords linked list + if coords._next == None: + offsets_to_coords[offset] = None + else: + offsets_to_coords[offset] = coords._next + levels_remain = True + + # Tile the Multi-index columns + for level_key, level_bins in level_dict.items(): + level_bins = \ + np.repeat(level_bins, filter.stride) + tile_factor = data_size / len(level_bins) + level_bins = np.tile(level_bins, tile_factor) + level_dict[level_key] = level_bins + + # Append the multi-index column to the DataFrame + df = pd.concat([df, pd.DataFrame(level_dict)], + axis=1) + + # Create DataFrame column for distribcell instances IDs + # NOTE: This is performed regardless of whether the user + # requests Summary geomeric information + filter_bins = np.arange(filter.num_bins) + filter_bins = np.repeat(filter_bins, filter.stride) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + df[filter.type] = filter_bins + + # energy, energyout filters + elif 'energy' in filter.type: + bins = filter.bins + num_bins = filter.num_bins + + # Create strings for + template = '{0:.1e} - {1:.1e}' + filter_bins = [] + for i in range(num_bins): + filter_bins.append(template.format(bins[i], bins[i+1])) + + # Tile the energy bins into a DataFrame column + filter_bins = np.repeat(filter_bins, filter.stride) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + df[filter.type + ' [MeV]'] = filter_bins + + # universe, material, surface, cell, and cellborn filters + else: + filter_bins = np.repeat(filter.bins, filter.stride) + tile_factor = data_size / len(filter_bins) + filter_bins = np.tile(filter_bins, tile_factor) + df[filter.type] = filter_bins + + # Include DataFrame column for nuclides if user requested it + if nuclides: + nuclides = [] + + for nuclide in self.nuclides: + # Write Nuclide name if Summary info was linked with StatePoint + if isinstance(nuclide, Nuclide): + nuclides.append(nuclide.name) + 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, tile_factor) + + # Include column for scores if user requested it + if scores: + tile_factor = data_size / len(self.scores) + df['score'] = np.tile(self.scores, tile_factor) + + # Append columns with mean, std. dev. for each tally bin + df['mean'] = self.mean.ravel() + df['std. dev.'] = self.std_dev.ravel() + + df.index.name = 'bin' + df = df.dropna(axis=1) + return df + def export_results(self, filename='tally-results', directory='.', format='hdf5', append=True): - """Returns a tally score value given a list of filters to satisfy. + """Exports tallly results to an HDF5 or Python pickle binary file. Parameters ---------- filename : str - The name of the file for the results (default is 'tally-results') + The name of the file for the results (default is 'tally-results') directory : str - The name of the directory for the results (default is '.') + The name of the directory for the results (default is '.') format : str - The format for the exported file - HDF5 ('hdf5', default), Python - pickle ('pkl'), comma-separated values ('csv') files are supported. + The format for the exported file - HDF5 ('hdf5', default) and + Python pickle ('pkl') files are supported append : bool - Whether or not to append the results to the file (default is True) + Whether or not to append the results to the file (default is True) + + Raises + ------ + KeyError : An error when this routine is called before the Tally + is populated with data by the StatePoint.read_results() routine. """ + # Ensure that StatePoint.read_results() was called first + if self._sum is None or self._sum_sq is None: + msg = 'The Tally ID={0} has no data to export. Call the ' \ + 'StatePoint.read_results() routine before using ' \ + 'Tally.export_results(...)'.format(self.id) + raise KeyError(msg) + if not is_string(filename): msg = 'Unable to export the results for Tally ID={0} to ' \ - 'filename={1} since it is not a ' \ - 'string'.format(self._id, filename) + 'filename="{1}" since it is not a ' \ + 'string'.format(self.id, filename) raise ValueError(msg) elif not is_string(directory): msg = 'Unable to export the results for Tally ID={0} to ' \ - 'directory={1} since it is not a ' \ - 'string'.format(self._id, directory) + 'directory="{1}" since it is not a ' \ + 'string'.format(self.id, directory) raise ValueError(msg) elif not format in ['hdf5', 'pkl', 'csv']: - msg = 'Unable to export the results for Tally ID={0} to ' \ - 'format {1} since it is not supported'.format(self._id, format) + msg = 'Unable to export the results for Tally ID={0} to format ' \ + '"{1}" since it is not supported'.format(self.id, format) raise ValueError(msg) elif not isinstance(append, (bool, np.bool)): msg = 'Unable to export the results for Tally ID={0} since the ' \ - 'append parameters is not True/False'.format(self._id, append) + 'append parameters is not True/False'.format(self.id, append) raise ValueError(msg) # Make directory if it does not exist @@ -744,19 +1253,19 @@ class Tally(object): tally_results = h5py.File(filename, 'w') # Create an HDF5 group within the file for this particular Tally - tally_group = tally_results.create_group('Tally-{0}'.format(self._id)) + tally_group = tally_results.create_group('Tally-{0}'.format(self.id)) # Add basic Tally data to the HDF5 group - tally_group.create_dataset('id', data=self._id) - tally_group.create_dataset('name', data=self._name) - tally_group.create_dataset('estimator', data=self._estimator) - tally_group.create_dataset('scores', data=np.array(self._scores)) + tally_group.create_dataset('id', data=self.id) + tally_group.create_dataset('name', data=self.name) + tally_group.create_dataset('estimator', data=self.estimator) + tally_group.create_dataset('scores', data=np.array(self.scores)) # Add a string array of the nuclides to the HDF5 group nuclides = [] - for nuclide in self._nuclides: - nuclides.append(nuclide._name) + for nuclide in self.nuclides: + nuclides.append(nuclide.name) tally_group.create_dataset('nuclides', data=np.array(nuclides)) @@ -764,14 +1273,14 @@ class Tally(object): # Create an HDF5 sub-group for the Filters filter_group = tally_group.create_group('filters') - for filter in self._filters: - filter_group.create_dataset(filter._type, data=filter._bins) + for filter in self.filters: + filter_group.create_dataset(filter.type, data=filter.bins) # Add all results to the main HDF5 group for the Tally - tally_group.create_dataset('sum', data=self._sum) - tally_group.create_dataset('sum_sq', data=self._sum_sq) - tally_group.create_dataset('mean', data=self._mean) - tally_group.create_dataset('std_dev', data=self._std_dev) + tally_group.create_dataset('sum', data=self.sum) + tally_group.create_dataset('sum_sq', data=self.sum_sq) + tally_group.create_dataset('mean', data=self.mean) + tally_group.create_dataset('std_dev', data=self.std_dev) # Close the Tally results HDF5 file tally_results.close() @@ -791,20 +1300,20 @@ class Tally(object): tally_results = {} # Create a nested dictionary within the file for this particular Tally - tally_results['Tally-{0}'.format(self._id)] = {} - tally_group = tally_results['Tally-{0}'.format(self._id)] + tally_results['Tally-{0}'.format(self.id)] = {} + tally_group = tally_results['Tally-{0}'.format(self.id)] # Add basic Tally data to the nested dictionary - tally_group['id'] = self._id - tally_group['name'] = self._name - tally_group['estimator'] = self._estimator - tally_group['scores'] = np.array(self._scores) + tally_group['id'] = self.id + tally_group['name'] = self.name + tally_group['estimator'] = self.estimator + tally_group['scores'] = np.array(self.scores) # Add a string array of the nuclides to the HDF5 group nuclides = [] - for nuclide in self._nuclides: - nuclides.append(nuclide._name) + for nuclide in self.nuclides: + nuclides.append(nuclide.name) tally_group['nuclides']= np.array(nuclides) @@ -812,14 +1321,14 @@ class Tally(object): tally_group['filters'] = {} filter_group = tally_group['filters'] - for filter in self._filters: - filter_group[filter._type] = filter._bins + for filter in self.filters: + filter_group[filter.type] = filter.bins # Add all results to the main sub-dictionary for the Tally - tally_group['sum'] = self._sum - tally_group['sum_sq'] = self._sum_sq - tally_group['mean'] = self._mean - tally_group['std_dev'] = self._std_dev + tally_group['sum'] = self.sum + tally_group['sum_sq'] = self.sum_sq + tally_group['mean'] = self.mean + tally_group['std_dev'] = self.std_dev # Pickle the Tally results to a file pickle.dump(tally_results, open(filename, 'wb')) diff --git a/src/utils/openmc/trigger.py b/src/utils/openmc/trigger.py index f2504cb92f..140fe93515 100644 --- a/src/utils/openmc/trigger.py +++ b/src/utils/openmc/trigger.py @@ -57,7 +57,7 @@ class Trigger(object): if not trigger_type in ['variance', 'std_dev', 'rel_err']: msg = 'Unable to create a tally trigger with ' \ - 'type {0}'.format(trigger_type) + 'type "{0}"'.format(trigger_type) raise ValueError(msg) self._trigger_type = trigger_type @@ -68,7 +68,7 @@ class Trigger(object): if not is_float(threshold): msg = 'Unable to set a tally trigger threshold with ' \ - 'threshold {0}'.format(threshold) + 'threshold "{0}"'.format(threshold) raise ValueError(msg) self._threshold = threshold @@ -77,7 +77,7 @@ class Trigger(object): def add_score(self, score): if not is_string(score): - msg = 'Unable to add score {0} to tally trigger since ' \ + msg = 'Unable to add score "{0}" to tally trigger since ' \ 'it is not a string'.format(score) raise ValueError(msg) diff --git a/tests/test_cmfd_feed/results.py b/tests/test_cmfd_feed/results.py index 125d13a8ea..869b785f36 100644 --- a/tests/test_cmfd_feed/results.py +++ b/tests/test_cmfd_feed/results.py @@ -18,15 +18,15 @@ else: sp.read_results() # extract tally results and convert to vector -tally1 = sp.get_tally('flux', [openmc.Filter('mesh', [1])], \ - [-1], estimator='tracklength') -tally2 = sp.get_tally('flux', [openmc.Filter('mesh', [2])], \ - [-1], estimator='analog') -tally3 = sp.get_tally('nu-fission', [openmc.Filter('mesh', [2])], \ - [-1], estimator='analog') -tally4 = sp.get_tally('current', [openmc.Filter('mesh', [2]), \ +tally1 = sp.get_tally(scores=['flux'], filters=[openmc.Filter('mesh', [1])], \ + estimator='tracklength') +tally2 = sp.get_tally(scores=['flux'], filters=[openmc.Filter('mesh', [2])], \ + estimator='analog') +tally3 = sp.get_tally(scores=['nu-fission'], \ + filters=[openmc.Filter('mesh', [2])], estimator='analog') +tally4 = sp.get_tally(scores=['current'], filters=[openmc.Filter('mesh', [2]), \ openmc.Filter('surface', [1,2,3,4,5,6])], \ - [-1], estimator='analog') + estimator='analog') results1 = np.zeros((tally1.sum.size + tally1.sum.size, )) results1[0::2] = tally1.sum.ravel() diff --git a/tests/test_cmfd_nofeed/results.py b/tests/test_cmfd_nofeed/results.py index e104fbd9ae..d22c2b2b42 100644 --- a/tests/test_cmfd_nofeed/results.py +++ b/tests/test_cmfd_nofeed/results.py @@ -18,15 +18,15 @@ else: sp.read_results() # extract tally results and convert to vector -tally1 = sp.get_tally('flux', [openmc.Filter('mesh', [1])], \ - [-1], estimator='tracklength') -tally2 = sp.get_tally('flux', [openmc.Filter('mesh', [2])], \ - [-1], estimator='analog') -tally3 = sp.get_tally('nu-fission', [openmc.Filter('mesh', [2])], \ - [-1], estimator='analog') -tally4 = sp.get_tally('current', [openmc.Filter('mesh', [2]), \ +tally1 = sp.get_tally(scores=['flux'], filters=[openmc.Filter('mesh', [1])], \ + estimator='tracklength') +tally2 = sp.get_tally(scores=['flux'], filters=[openmc.Filter('mesh', [2])], \ + estimator='analog') +tally3 = sp.get_tally(scores=['nu-fission'], \ + filters=[openmc.Filter('mesh', [2])], estimator='analog') +tally4 = sp.get_tally(scores=['current'], filters=[openmc.Filter('mesh', [2]), \ openmc.Filter('surface', [1,2,3,4,5,6])], \ - [-1], estimator='analog') + estimator='analog') results1 = np.zeros((tally1.sum.size + tally1.sum.size, )) results1[0::2] = tally1.sum.ravel() diff --git a/tests/test_filter_distribcell/results.py b/tests/test_filter_distribcell/results.py index 155e7ebf02..a1dc38ebd4 100644 --- a/tests/test_filter_distribcell/results.py +++ b/tests/test_filter_distribcell/results.py @@ -21,89 +21,74 @@ sp3.compute_stdev() # analyze sp1 # compare distrib sum to cell -sp1_t1 = sp1._tallies[1] -sp1_t2 = sp1._tallies[2] +sp1_t1 = sp1.tallies[1] +sp1_t2 = sp1.tallies[2] -distribcell_filter = sp1_t1.find_filter(filter_type='distribcell', bins=[2]) -sp1_t1_d1 = sp1_t1.get_value(score='total', filters=[distribcell_filter], - filter_bins=[0], value='mean') -sp1_t1_d2 = sp1_t1.get_value(score='total', filters=[distribcell_filter], - filter_bins=[1], value='mean') -sp1_t1_d3 = sp1_t1.get_value(score='total', filters=[distribcell_filter], - filter_bins=[2], value='mean') -sp1_t1_d4 = sp1_t1.get_value(score='total', filters=[distribcell_filter], - filter_bins=[3], value='mean') +sp1_t1_d1 = sp1_t1.get_values(scores=['total'], filters=['distribcell'], + filter_bins=[(0,)], value='mean') +sp1_t1_d2 = sp1_t1.get_values(scores=['total'], filters=['distribcell'], + filter_bins=[(1,)], value='mean') +sp1_t1_d3 = sp1_t1.get_values(scores=['total'], filters=['distribcell'], + filter_bins=[(2,)], value='mean') +sp1_t1_d4 = sp1_t1.get_values(scores=['total'], filters=['distribcell'], + filter_bins=[(3,)], value='mean') sp1_t1_sum = sp1_t1_d1 + sp1_t1_d2 + sp1_t1_d3 + sp1_t1_d4 -cell_filter = sp1_t2.find_filter(filter_type='cell', bins=[2]) -sp1_t2_c201 = sp1_t2.get_value(score='total', filters=[cell_filter], - filter_bins=[2], value='mean') +cell_filter = sp1_t2.find_filter(filter_type='cell') +sp1_t2_c201 = sp1_t2.get_values(scores=['total'], value='mean') # analyze sp2 -sp2_t1 = sp2._tallies[1] -cell_filter = sp2_t1.find_filter(filter_type='cell', bins=[2,4,6,8]) -sp2_t1_c201 = sp2_t1.get_value(score='total', filters=[cell_filter], - filter_bins=[2], value='mean') -sp2_t1_c203 = sp2_t1.get_value(score='total', filters=[cell_filter], - filter_bins=[4], value='mean') -sp2_t1_c205 = sp2_t1.get_value(score='total', filters=[cell_filter], - filter_bins=[6], value='mean') -sp2_t1_c207 = sp2_t1.get_value(score='total', filters=[cell_filter], - filter_bins=[8], value='mean') +sp2_t1 = sp2.tallies[1] +sp2_t1_c201 = sp2_t1.get_values(scores=['total'], filters=['cell'], + filter_bins=[(2,)], value='mean') +sp2_t1_c203 = sp2_t1.get_values(scores=['total'], filters=['cell'], + filter_bins=[(4,)], value='mean') +sp2_t1_c205 = sp2_t1.get_values(scores=['total'], filters=['cell'], + filter_bins=[(6,)], value='mean') +sp2_t1_c207 = sp2_t1.get_values(scores=['total'], filters=['cell'], + filter_bins=[(8,)], value='mean') # analyze sp3 -sp3_t1 = sp3._tallies[1] -sp3_t2 = sp3._tallies[2] -sp3_t3 = sp3._tallies[3] -sp3_t4 = sp3._tallies[4] -sp3_t5 = sp3._tallies[5] -sp3_t6 = sp3._tallies[6] +sp3_t1 = sp3.tallies[1] +sp3_t2 = sp3.tallies[2] +sp3_t3 = sp3.tallies[3] +sp3_t4 = sp3.tallies[4] +sp3_t5 = sp3.tallies[5] +sp3_t6 = sp3.tallies[6] -cell_filter = sp3_t1.find_filter(filter_type='cell', bins=[1]) -distribcell_filter = sp3_t2.find_filter(filter_type='distribcell', bins=[1]) -sp3_t1_c1 = sp3_t1.get_value(score='total', filters=[cell_filter], - filter_bins=[1], value='mean') -sp3_t2_d1 = sp3_t2.get_value(score='total', filters=[distribcell_filter], - filter_bins=[0], value='mean') +sp3_t1_c1 = sp3_t1.get_values(scores=['total'], filters=['cell'], + filter_bins=[(1,)], value='mean') +sp3_t2_d1 = sp3_t2.get_values(scores=['total'], filters=['distribcell'], + filter_bins=[(0,)], value='mean') -cell_filter = sp3_t3.find_filter(filter_type='cell', bins=[26]) -distribcell_filter = sp3_t4.find_filter(filter_type='distribcell', bins=[26]) -sp3_t3_c60 = sp3_t3.get_value(score='total', filters=[cell_filter], - filter_bins=[26], value='mean') -sp3_t4_d = 0 -for i in range(241): - sp3_t4_d += sp3_t4.get_value(score='total', filters=[distribcell_filter], - filter_bins=[i], value='mean') +sp3_t3_c60 = sp3_t3.get_values(scores=['total'], filters=['cell'], + filter_bins=[(26,)], value='mean') +sp3_t4_d = sum(sp3_t4.get_values(scores=['total'], value='mean')) -cell_filter = sp3_t5.find_filter(filter_type='cell', bins=[19]) -distribcell_filter = sp3_t6.find_filter(filter_type='distribcell', bins=[19]) -sp3_t5_c27 = sp3_t5.get_value(score='total', filters=[cell_filter], - filter_bins=[19], value='mean') -sp3_t6_d = 0 -for i in range(63624): - sp3_t6_d += sp3_t6.get_value(score='total', filters=[distribcell_filter], - filter_bins=[i], value='mean') +sp3_t5_c27 = sp3_t5.get_values(scores=['total'], filters=['cell'], + filter_bins=[(19,)], value='mean') +sp3_t6_d = sum(sp3_t6.get_values(scores=['total'], value='mean')) # set up output string outstr = '' outstr += "{0:12.6E}\n".format(sp1_t1_sum) -outstr += "{0:12.6E}\n".format(sp1_t2_c201) -outstr += "{0:12.6E}\n".format(sp2_t1_c201) -outstr += "{0:12.6E}\n".format(sp1_t1_d4) -outstr += "{0:12.6E}\n".format(sp2_t1_c203) -outstr += "{0:12.6E}\n".format(sp1_t1_d1) -outstr += "{0:12.6E}\n".format(sp2_t1_c205) -outstr += "{0:12.6E}\n".format(sp1_t1_d3) -outstr += "{0:12.6E}\n".format(sp2_t1_c207) -outstr += "{0:12.6E}\n".format(sp1_t1_d2) -outstr += "{0:12.6E}\n".format(sp3_t1_c1) -outstr += "{0:12.6E}\n".format(sp3_t2_d1) -outstr += "{0:12.6E}\n".format(sp3_t3_c60) +outstr += "{0:12.6E}\n".format(sp1_t2_c201[()]) +outstr += "{0:12.6E}\n".format(sp2_t1_c201[()]) +outstr += "{0:12.6E}\n".format(sp1_t1_d4[()]) +outstr += "{0:12.6E}\n".format(sp2_t1_c203[()]) +outstr += "{0:12.6E}\n".format(sp1_t1_d1[()]) +outstr += "{0:12.6E}\n".format(sp2_t1_c205[()]) +outstr += "{0:12.6E}\n".format(sp1_t1_d3[()]) +outstr += "{0:12.6E}\n".format(sp2_t1_c207[()]) +outstr += "{0:12.6E}\n".format(sp1_t1_d2[()]) +outstr += "{0:12.6E}\n".format(sp3_t1_c1[()]) +outstr += "{0:12.6E}\n".format(sp3_t2_d1[()]) +outstr += "{0:12.6E}\n".format(sp3_t3_c60[()]) outstr += "{0:12.6E}\n".format(sp3_t4_d) -outstr += "{0:12.6E}\n".format(sp3_t5_c27) +outstr += "{0:12.6E}\n".format(sp3_t5_c27[()]) outstr += "{0:12.6E}\n".format(sp3_t6_d) # write results to file