From 8762642dda47965b2ea29cdda5c5c6c585937d02 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Fri, 25 Dec 2015 12:23:53 -0800 Subject: [PATCH 01/18] changed mean and std_dev direct access to property getters in tallies.py --- openmc/summary.py | 1 - openmc/tallies.py | 36 ++++++++++++++++++------------------ 2 files changed, 18 insertions(+), 19 deletions(-) diff --git a/openmc/summary.py b/openmc/summary.py index 45742b0c8a..ff6b59d3c7 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -532,7 +532,6 @@ class Summary(object): # If this is a moment, use generic moment order pattern = r'-n$|-pn$|-yn$' score = re.sub(pattern, '-' + moments[j].decode(), score) - tally.add_score(score) # Read filter metadata diff --git a/openmc/tallies.py b/openmc/tallies.py index a8569b7272..cdf6c1f87e 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -1847,7 +1847,7 @@ class Tally(object): data = self.get_values( filters=filters, filter_bins=filter_bins, value='mean') indices = self.get_filter_indices(filters, filter_bins) - self._mean[indices, :, :] = data + self.mean[indices, :, :] = data # Adjust the std_dev data array to relect the new filter order if self.std_dev is not None: @@ -1856,7 +1856,7 @@ class Tally(object): data = self.get_values( filters=filters, filter_bins=filter_bins, value='std_dev') indices = self.get_filter_indices(filters, filter_bins) - self._std_dev[indices, :, :] = data + self.std_dev[indices, :, :] = data def _swap_nuclides(self, nuclide1, nuclide2): """Reverse the ordering of two nuclides in this tally @@ -1913,17 +1913,17 @@ class Tally(object): # Adjust the mean data array to relect the new nuclide order if self.mean is not None: - nuclide1_mean = self._mean[:, nuclide1_index, :].copy() - nuclide2_mean = self._mean[:, nuclide2_index, :].copy() - self._mean[:, nuclide2_index, :] = nuclide1_mean - self._mean[:, nuclide1_index, :] = nuclide2_mean + nuclide1_mean = self.mean[:, nuclide1_index, :].copy() + nuclide2_mean = self.mean[:, nuclide2_index, :].copy() + self.mean[:, nuclide2_index, :] = nuclide1_mean + self.mean[:, nuclide1_index, :] = nuclide2_mean # Adjust the std_dev data array to relect the new nuclide order if self.std_dev is not None: - nuclide1_std_dev = self._std_dev[:, nuclide1_index, :].copy() - nuclide2_std_dev = self._std_dev[:, nuclide2_index, :].copy() - self._std_dev[:, nuclide2_index, :] = nuclide1_std_dev - self._std_dev[:, nuclide1_index, :] = nuclide2_std_dev + nuclide1_std_dev = self.std_dev[:, nuclide1_index, :].copy() + nuclide2_std_dev = self.std_dev[:, nuclide2_index, :].copy() + self.std_dev[:, nuclide2_index, :] = nuclide1_std_dev + self.std_dev[:, nuclide1_index, :] = nuclide2_std_dev def _swap_scores(self, score1, score2): """Reverse the ordering of two scores in this tally @@ -1985,17 +1985,17 @@ class Tally(object): # Adjust the mean data array to relect the new nuclide order if self.mean is not None: - score1_mean = self._mean[:, :, score1_index].copy() - score2_mean = self._mean[:, :, score2_index].copy() - self._mean[:, :, score2_index] = score1_mean - self._mean[:, :, score1_index] = score2_mean + score1_mean = self.mean[:, :, score1_index].copy() + score2_mean = self.mean[:, :, score2_index].copy() + self.mean[:, :, score2_index] = score1_mean + self.mean[:, :, score1_index] = score2_mean # Adjust the std_dev data array to relect the new nuclide order if self.std_dev is not None: - score1_std_dev = self._std_dev[:, :, score1_index].copy() - score2_std_dev = self._std_dev[:, :, score2_index].copy() - self._std_dev[:, :, score2_index] = score1_std_dev - self._std_dev[:, :, score1_index] = score2_std_dev + score1_std_dev = self.std_dev[:, :, score1_index].copy() + score2_std_dev = self.std_dev[:, :, score2_index].copy() + self.std_dev[:, :, score2_index] = score1_std_dev + self.std_dev[:, :, score1_index] = score2_std_dev def __add__(self, other): """Adds this tally to another tally or scalar value. From 4e939d382024660c8213b8d3ae9a5576daaa5339 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Mon, 28 Dec 2015 10:59:27 -0800 Subject: [PATCH 02/18] changed RectLattice universes ordering from [x][y][z] to [z][y][x] --- openmc/universe.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmc/universe.py b/openmc/universe.py index 98367c381b..14a3007afb 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -1119,7 +1119,7 @@ class RectLattice(Lattice): for z in range(self._dimension[2]): for y in range(self._dimension[1]): for x in range(self._dimension[0]): - universe = self._universes[x][y][z] + universe = self._universes[z][y][x] # Append Universe ID to the Lattice XML subelement universe_ids += '{0} '.format(universe._id) @@ -1137,7 +1137,7 @@ class RectLattice(Lattice): else: for y in range(self._dimension[1]): for x in range(self._dimension[0]): - universe = self._universes[x][y] + universe = self._universes[y][x] # Append Universe ID to Lattice XML subelement universe_ids += '{0} '.format(universe._id) From 022d8d3a44c0151dfad75c2baf72aa6f199dd08f Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 28 Dec 2015 23:14:05 -0500 Subject: [PATCH 03/18] Replace filter % offset w/ cell % distribcell_ind --- openmc/filter.py | 14 ---------- openmc/statepoint.py | 4 --- src/geometry_header.F90 | 9 ++++--- src/initialize.F90 | 58 ++++++++++++++--------------------------- src/output.F90 | 23 ++++++++-------- src/state_point.F90 | 1 - src/summary.F90 | 3 ++- src/tally.F90 | 6 +++-- src/tally_header.F90 | 1 - 9 files changed, 44 insertions(+), 75 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index 04935b8edc..dce3d4924c 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -42,8 +42,6 @@ class Filter(object): The number of filter bins mesh : Mesh or None A Mesh object for 'mesh' type filters. - offset : Integral - A value used to index tally bins for 'distribcell' tallies. stride : Integral The number of filter, nuclide and score bins within each of this filter's bins. @@ -57,7 +55,6 @@ class Filter(object): self._num_bins = 0 self._bins = None self._mesh = None - self._offset = -1 self._stride = None if type is not None: @@ -93,7 +90,6 @@ class Filter(object): clone._bins = copy.deepcopy(self.bins, memo) clone._num_bins = self.num_bins clone._mesh = copy.deepcopy(self.mesh, memo) - clone._offset = self.offset clone._stride = self.stride memo[id(self)] = clone @@ -108,7 +104,6 @@ class Filter(object): string = 'Filter\n' string += '{0: <16}{1}{2}\n'.format('\tType', '=\t', self.type) string += '{0: <16}{1}{2}\n'.format('\tBins', '=\t', self.bins) - string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self.offset) return string @property @@ -134,10 +129,6 @@ class Filter(object): def mesh(self): return self._mesh - @property - def offset(self): - return self._offset - @property def stride(self): return self._stride @@ -226,11 +217,6 @@ class Filter(object): self.type = 'mesh' self.bins = self.mesh.id - @offset.setter - def offset(self, offset): - cv.check_type('filter offset', offset, Integral) - self._offset = offset - @stride.setter def stride(self, stride): cv.check_type('filter stride', stride, Integral) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 2edf9badd4..bc020cfbce 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -360,9 +360,6 @@ class StatePoint(object): # Read the Filter type filter_type = self._f['{0}{1}/type'.format(subbase, j)].value.decode() - # Read the Filter offset - offset = self._f['{0}{1}/offset'.format(subbase, j)].value - n_bins = self._f['{0}{1}/n_bins'.format(subbase, j)].value # Read the bin values @@ -370,7 +367,6 @@ class StatePoint(object): # Create Filter object filter = openmc.Filter(filter_type, bins) - filter.offset = offset filter.num_bins = n_bins if filter_type == 'mesh': diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 6ca13740ba..684de883c2 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -131,10 +131,13 @@ module geometry_header integer, allocatable :: offset (:) ! Distribcell offset for tally ! counter integer, allocatable :: region(:) ! Definition of spatial region as - ! Boolean expression of half-spaces + ! Boolean expression of half-spaces integer, allocatable :: rpn(:) ! Reverse Polish notation for region - ! expression - logical :: simple ! Is the region simple (intersections only) + ! expression + logical :: simple ! Is the region simple (intersections + ! only) + integer :: distribcell_ind ! Index corresponding to this cell in + ! distribcell arrays ! Rotation matrix and translation vector real(8), allocatable :: translation(:) diff --git a/src/initialize.F90 b/src/initialize.F90 index 9d63eb4e93..acfe92b921 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1050,12 +1050,12 @@ contains do i = 1, n_tallies t => tallies(i) - do j = 1, t%n_filters - filter => t%filters(j) + do j = 1, t % n_filters + filter => t % filters(j) - if (filter%type == FILTER_DISTRIBCELL) then - if (.not. cell_list%contains(filter%int_bins(1))) then - call cell_list%add(filter%int_bins(1)) + if (filter % type == FILTER_DISTRIBCELL) then + if (.not. cell_list % contains(filter % int_bins(1))) then + call cell_list % add(filter % int_bins(1)) end if end if @@ -1066,8 +1066,8 @@ contains ! to determine the number of offset tables to allocate do i = 1, n_universes univ => universes(i) - do j = 1, univ%n_cells - if (cell_list%contains(univ%cells(j))) then + do j = 1, univ % n_cells + if (cell_list % contains(univ % cells(j))) then n_maps = n_maps + 1 end if end do @@ -1089,56 +1089,38 @@ contains do i = 1, n_universes univ => universes(i) - do j = 1, univ%n_cells + do j = 1, univ % n_cells + if (.not. cell_list % contains(univ % cells(j))) cycle - if (cell_list%contains(univ%cells(j))) then + cells(univ % cells(j)) % distribcell_ind = k - ! Loop over all tallies - do l = 1, n_tallies - t => tallies(l) - - do m = 1, t%n_filters - filter => t%filters(m) - - ! Loop over only distribcell filters - ! If filter points to cell we just found, set offset index - if (filter%type == FILTER_DISTRIBCELL) then - if (filter%int_bins(1) == univ%cells(j)) then - filter%offset = k - end if - end if - - end do - end do - - univ_list(k) = univ%id - k = k + 1 - end if + univ_list(k) = univ % id + k = k + 1 end do end do ! Allocate the offset tables for lattices do i = 1, n_lattices - lat => lattices(i)%obj + lat => lattices(i) % obj select type(lat) type is (RectLattice) - allocate(lat%offset(n_maps, lat%n_cells(1), lat%n_cells(2), & - lat%n_cells(3))) + allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & + lat % n_cells(3))) type is (HexLattice) - allocate(lat%offset(n_maps, 2 * lat%n_rings - 1, & - 2 * lat%n_rings - 1, lat%n_axial)) + allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & + 2 * lat % n_rings - 1, lat % n_axial)) end select - lat%offset(:, :, :, :) = 0 + lat % offset(:, :, :, :) = 0 end do ! Allocate offset table for fill cells do i = 1, n_cells - if (cells(i)%material == NONE) then - allocate(cells(i)%offset(n_maps)) + if (cells(i) % material == NONE) then + allocate(cells(i) % offset(n_maps)) end if end do diff --git a/src/output.F90 b/src/output.F90 index c31b20b70b..8be6c09dab 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1378,8 +1378,7 @@ contains label = '' univ => universes(BASE_UNIVERSE) offset = 0 - call find_offset(t % filters(i_filter) % offset, & - t % filters(i_filter) % int_bins(1), & + call find_offset(t % filters(i_filter) % int_bins(1), & univ, bin-1, offset, label) case (FILTER_SURFACE) i = t % filters(i_filter) % int_bins(bin) @@ -1413,15 +1412,15 @@ contains ! with the given offset !=============================================================================== - recursive subroutine find_offset(map, goal, univ, final, offset, path) + recursive subroutine find_offset(goal, univ, final, offset, path) - integer, intent(in) :: map ! Index in maps vector - integer, intent(in) :: goal ! The target cell ID + integer, intent(in) :: goal ! The target cell index type(Universe), intent(in) :: univ ! Universe to begin search integer, intent(in) :: final ! Target offset integer, intent(inout) :: offset ! Current offset character(*), intent(inout) :: path ! Path to offset + integer :: map ! Index in maps vector integer :: i, j ! Index over cells integer :: k, l, m ! Indices in lattice integer :: old_k, old_l, old_m ! Previous indices in lattice @@ -1436,6 +1435,9 @@ contains type(Universe), pointer :: next_univ ! Next universe to loop through class(Lattice), pointer :: lat ! Pointer to current lattice + ! Get the distribcell index for this cell + map = cells(goal) % distribcell_ind + n = univ % n_cells ! Write to the geometry stack @@ -1537,7 +1539,7 @@ contains offset = c % offset(map) + offset next_univ => universes(c % fill) - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return ! ==================================================================== @@ -1577,7 +1579,7 @@ contains path = trim(path) // "(" // trim(to_str(k)) // & "," // trim(to_str(l)) // "," // & trim(to_str(m)) // ")" - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return else old_m = m @@ -1593,7 +1595,7 @@ contains path = trim(path) // "(" // trim(to_str(old_k)) // & "," // trim(to_str(old_l)) // "," // & trim(to_str(old_m)) // ")" - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return end if @@ -1638,8 +1640,7 @@ contains trim(to_str(k - lat % n_rings)) // "," // & trim(to_str(l - lat % n_rings)) // "," // & trim(to_str(m)) // ")" - call find_offset(map, goal, next_univ, final, offset, & - path) + call find_offset(goal, next_univ, final, offset, path) return else old_m = m @@ -1656,7 +1657,7 @@ contains trim(to_str(old_k - lat % n_rings)) // "," // & trim(to_str(old_l - lat % n_rings)) // "," // & trim(to_str(old_m)) // ")" - call find_offset(map, goal, next_univ, final, offset, path) + call find_offset(goal, next_univ, final, offset, path) return end if diff --git a/src/state_point.F90 b/src/state_point.F90 index f9f4b5b757..046c857fd2 100644 --- a/src/state_point.F90 +++ b/src/state_point.F90 @@ -268,7 +268,6 @@ contains call write_dataset(filter_group, "type", "delayedgroup") end select - call write_dataset(filter_group, "offset", tally%filters(j)%offset) call write_dataset(filter_group, "n_bins", tally%filters(j)%n_bins) if (tally % filters(j) % type == FILTER_ENERGYIN .or. & tally % filters(j) % type == FILTER_ENERGYOUT .or. & diff --git a/src/summary.F90 b/src/summary.F90 index cb005d3be5..11c4cd434b 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -196,6 +196,8 @@ contains end do call write_dataset(cell_group, "region", adjustl(region_spec)) + call write_dataset(cell_group, "distribcell_ind", c % distribcell_ind) + call close_group(cell_group) end do CELL_LOOP @@ -540,7 +542,6 @@ contains filter_group = create_group(tally_group, "filter " // trim(to_str(j))) ! Write number of bins for this filter - call write_dataset(filter_group, "offset", t%filters(j)%offset) call write_dataset(filter_group, "n_bins", t%filters(j)%n_bins) ! Write filter bins diff --git a/src/tally.F90 b/src/tally.F90 index 09e30a2426..2ccd59a58e 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1712,6 +1712,7 @@ contains integer :: j integer :: n ! number of bins for single filter integer :: offset ! offset for distribcell + integer :: distribcell_ind ! index in distribcell arrays real(8) :: E ! particle energy real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively type(TallyObject), pointer :: t @@ -1756,12 +1757,13 @@ contains case (FILTER_DISTRIBCELL) ! determine next distribcell bin + distribcell_ind = cells(t % filters(i) % int_bins(1)) % distribcell_ind matching_bins(i) = NO_BIN_FOUND offset = 0 do j = 1, p % n_coord if (cells(p % coord(j) % cell) % type == CELL_FILL) then offset = offset + cells(p % coord(j) % cell) % & - offset(t % filters(i) % offset) + offset(distribcell_ind) elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then if (lattices(p % coord(j + 1) % lattice) % obj & % are_valid_indices([& @@ -1769,7 +1771,7 @@ contains p % coord(j + 1) % lattice_y, & p % coord(j + 1) % lattice_z])) then offset = offset + lattices(p % coord(j + 1) % lattice) % obj % & - offset(t % filters(i) % offset, & + offset(distribcell_ind, & p % coord(j + 1) % lattice_x, & p % coord(j + 1) % lattice_y, & p % coord(j + 1) % lattice_z) diff --git a/src/tally_header.F90 b/src/tally_header.F90 index 18b9219522..dcbba3d89e 100644 --- a/src/tally_header.F90 +++ b/src/tally_header.F90 @@ -55,7 +55,6 @@ module tally_header type TallyFilter integer :: type = NONE integer :: n_bins = 0 - integer :: offset = 0 ! Only used for distribcell filters integer, allocatable :: int_bins(:) real(8), allocatable :: real_bins(:) ! Only used for energy filters end type TallyFilter From 56e8e7d6818937cdcb57fc208532025e52768a68 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Tue, 29 Dec 2015 16:25:29 -0500 Subject: [PATCH 04/18] Update PyAPI with distribcell changes --- openmc/filter.py | 2 +- openmc/geometry.py | 20 +++++++++++++------- openmc/summary.py | 5 +++++ openmc/universe.py | 37 ++++++++++++++++++++++++++----------- src/input_xml.F90 | 3 ++- src/output.F90 | 14 +++++--------- 6 files changed, 52 insertions(+), 29 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index dce3d4924c..ea92f5257c 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -609,7 +609,7 @@ class Filter(object): # If this region is in Cell corresponding to the # distribcell filter bin, store it in dictionary if cell_id == self.bins[0]: - offset = openmc_geometry.get_offset(path, self.offset) + offset = openmc_geometry.get_offset(path) offsets_to_coords[offset] = coords # Each distribcell offset is a DataFrame bin diff --git a/openmc/geometry.py b/openmc/geometry.py index e848e0cddf..16f6a45e5a 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -42,7 +42,7 @@ class Geometry(object): self._root_universe = root_universe - def get_offset(self, path, filter_offset): + def get_offset(self, path): """Returns the corresponding location in the results array for a given path and filter number. This is primarily intended to post-processing result when a distribcell filter is used. @@ -55,8 +55,6 @@ class Geometry(object): lattice passed through. For the case of the lattice, a tuple should be provided to indicate which coordinates in the lattice should be entered. This should be in the form: (lat_id, i_x, i_y, i_z) - filter_offset : int - An integer that specifies which offset map the filter is using Returns ------- @@ -65,14 +63,22 @@ class Geometry(object): """ + # Find the distribcell index of the cell. + cells = self.get_all_cells() + if path[-1] in cells: + distribcell_ind = cells[path[-1]].distribcell_ind + else: + raise RuntimeError('Could not find cell {} specified in a \ + distribcell filter'.format(path[-1])) + # Return memoize'd offset if possible - if (path, filter_offset) in self._offsets: - offset = self._offsets[(path, filter_offset)] + if (path, distribcell_ind) in self._offsets: + offset = self._offsets[(path, distribcell_ind)] # Begin recursive call to compute offset starting with the base Universe else: - offset = self._root_universe.get_offset(path, filter_offset) - self._offsets[(path, filter_offset)] = offset + offset = self._root_universe.get_offset(path, distribcell_ind) + self._offsets[(path, distribcell_ind)] = offset # Return the final offset return offset diff --git a/openmc/summary.py b/openmc/summary.py index bc6551e7cb..9e2145b343 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -271,6 +271,11 @@ class Summary(object): cell.region = Region.from_expression( region, {s.id: s for s in self.surfaces.values()}) + # Get the distribcell index + ind = self._f['geometry/cells'][key]['distribcell_ind'].value + if ind != 0: + cell.distribcell_ind = ind + # Add the Cell to the global dictionary of all Cells self.cells[index] = cell diff --git a/openmc/universe.py b/openmc/universe.py index 98367c381b..b346ffd00d 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -63,6 +63,8 @@ class Cell(object): that is used to translate (shift) the universe. offsets : ndarray Array of offsets used for distributed cell searches + distribcell_ind : int + Index of this cell in distribcell arrays """ @@ -76,6 +78,7 @@ class Cell(object): self._rotation = None self._translation = None self._offsets = None + self._distribcell_ind = None def __eq__(self, other): if not isinstance(other, Cell): @@ -122,6 +125,8 @@ class Cell(object): string += '{0: <16}{1}{2}\n'.format('\tTranslation', '=\t', self._translation) string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets) + string += '{0: <16}{1}{2}\n'.format('\tDistribcell index', '=\t', + self._distribcell_ind) return string @@ -164,6 +169,10 @@ class Cell(object): def offsets(self): return self._offsets + @property + def distribcell_ind(self): + return self._distribcell_ind + @id.setter def id(self, cell_id): if cell_id is None: @@ -231,6 +240,11 @@ class Cell(object): cv.check_type('cell region', region, Region) self._region = region + @distribcell_ind.setter + def distribcell_ind(self, ind): + cv.check_type('distribcell index', ind, Integral) + self._distribcell_ind = ind + def add_surface(self, surface, halfspace): """Add a half-space to the list of half-spaces whose intersection defines the cell. @@ -271,7 +285,7 @@ class Cell(object): else: self.region = Intersection(self.region, region) - def get_offset(self, path, filter_offset): + def get_offset(self, path, distribcell_ind): # Get the current element and remove it from the list cell_id = path[0] path = path[1:] @@ -282,12 +296,12 @@ class Cell(object): # If the Cell is filled by a Universe elif self._type == 'fill': - offset = self._offsets[filter_offset-1] - offset += self._fill.get_offset(path, filter_offset) + offset = self._offsets[distribcell_ind-1] + offset += self._fill.get_offset(path, distribcell_ind) # If the Cell is filled by a Lattice else: - offset = self._fill.get_offset(path, filter_offset) + offset = self._fill.get_offset(path, distribcell_ind) return offset @@ -591,7 +605,7 @@ class Universe(object): self._cells.clear() - def get_offset(self, path, filter_offset): + def get_offset(self, path, distribcell_ind): # Get the current element and remove it from the list path = path[1:] @@ -599,7 +613,7 @@ class Universe(object): cell_id = path[0] # Make a recursive call to the Cell within this Universe - offset = self._cells[cell_id].get_offset(path, filter_offset) + offset = self._cells[cell_id].get_offset(path, distribcell_ind) # Return the offset computed at all nested Universe levels return offset @@ -1059,21 +1073,22 @@ class RectLattice(Lattice): cv.check_greater_than('lattice pitch', dim, 0.0) self._pitch = pitch - def get_offset(self, path, filter_offset): + def get_offset(self, path, distribcell_ind): # Get the current element and remove it from the list i = path[0] path = path[1:] # For 2D Lattices if len(self._dimension) == 2: - offset = self._offsets[i[1]-1, i[2]-1, 0, filter_offset-1] - offset += self._universes[i[1]][i[2]].get_offset(path, filter_offset) + offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_ind-1] + offset += self._universes[i[1]][i[2]].get_offset(path, + distribcell_ind) # For 3D Lattices else: - offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, filter_offset-1] + offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_ind-1] offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_offset(path, - filter_offset) + distribcell_ind) return offset diff --git a/src/input_xml.F90 b/src/input_xml.F90 index abd6e29092..a4317ed8e3 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1049,8 +1049,9 @@ contains do i = 1, n_cells c => cells(i) - ! Initialize the number of cell instances - this is a base case for distribcells + ! Initialize distribcell instances and distribcell index c % instances = 0 + c % distribcell_ind = NONE ! Get pointer to i-th cell node call get_list_item(node_cell_list, i, node_cell) diff --git a/src/output.F90 b/src/output.F90 index 8be6c09dab..5b938e98d3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1449,17 +1449,13 @@ contains ! Look through all cells in this universe do i = 1, n - - cell_index = univ % cells(i) - c => cells(cell_index) - - ! If the cell ID matches the goal and the offset matches final, - ! write to the geometry stack - if (cell_dict % get_key(c % id) == goal .AND. offset == final) then - path = trim(path) // "->" // to_str(c%id) + ! If the cell matches the goal and the offset matches final, write to the + ! geometry stack + if (univ % cells(i) == goal .AND. offset == final) then + c => cells(univ % cells(i)) + path = trim(path) // "->" // to_str(c % id) return end if - end do ! Find the fill cell or lattice cell that we need to enter From c4b83ced54a66dd52d229507edf2214beddea1f3 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Wed, 30 Dec 2015 14:43:33 -0500 Subject: [PATCH 05/18] Fixed bugs in OpenCG compatiblity and Summary Python modules for new z-y-z Universe ordering --- openmc/opencg_compatible.py | 6 +++--- openmc/summary.py | 4 ++-- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index afa57c78d9..62192355de 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -905,7 +905,7 @@ def get_opencg_lattice(openmc_lattice): for z in range(dimension[2]): for y in range(dimension[1]): for x in range(dimension[0]): - universe_id = universes[x][dimension[1]-y-1][z].id + universe_id = universes[z][dimension[1]-y-1][x].id universe_array[z][y][x] = unique_universes[universe_id] opencg_lattice = opencg.Lattice(lattice_id, name) @@ -963,7 +963,7 @@ def get_openmc_lattice(opencg_lattice): outer = opencg_lattice.outside # Initialize an empty array for the OpenMC nested Universes in this Lattice - universe_array = np.ndarray(tuple(np.array(dimension)), + universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), dtype=openmc.Universe) # Create OpenMC Universes for each unique nested Universe in this Lattice @@ -977,7 +977,7 @@ def get_openmc_lattice(opencg_lattice): for y in range(dimension[1]): for x in range(dimension[0]): universe_id = universes[z][y][x].id - universe_array[x][y][z] = unique_universes[universe_id] + universe_array[z][y][x] = unique_universes[universe_id] # Reverse y-dimension in array to match ordering in OpenCG universe_array = universe_array[:, ::-1, :] diff --git a/openmc/summary.py b/openmc/summary.py index bc6551e7cb..3d0e37ddad 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -356,8 +356,8 @@ class Summary(object): self.get_universe_by_id(universe_ids[x, y, z]) # Transpose, reverse y-dimension for appropriate ordering - shape = universes.shape - universes = np.transpose(universes, (1, 0, 2)) + shape = universes.shape[::-1] + universes = np.transpose(universes, (2, 1, 0)) universes.shape = shape universes = universes[:, ::-1, :] lattice.universes = universes From a28c9281e796fbc719a62bc73fdaf8187c23d260 Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Wed, 30 Dec 2015 13:29:14 -0800 Subject: [PATCH 06/18] added documentation for updated summary file format and incremented the REVISION_SUMMARY constant --- docs/source/usersguide/output/summary.rst | 5 +++++ src/constants.F90 | 2 +- 2 files changed, 6 insertions(+), 1 deletion(-) diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index f87f60c4a0..ce1c1221a8 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -306,6 +306,11 @@ The current revision of the summary file format is 1. than the number of user-specified scores since each score might have multiple scoring bins, e.g., scatter-PN. +**/tallies/tally /moment_orders** (*char[][]*) + + Tallying moment orders for Legendre and spherical harmonic tally expansions + (*e.g.*, 'P2', 'Y1,2', etc.). + **/tallies/tally /score_bins** (*char[][]*) Scoring bins for the tally. diff --git a/src/constants.F90 b/src/constants.F90 index ba77f35aba..2db3ec3c6c 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -14,7 +14,7 @@ module constants integer, parameter :: REVISION_STATEPOINT = 14 integer, parameter :: REVISION_PARTICLE_RESTART = 1 integer, parameter :: REVISION_TRACK = 1 - integer, parameter :: REVISION_SUMMARY = 1 + integer, parameter :: REVISION_SUMMARY = 2 ! ============================================================================ ! ADJUSTABLE PARAMETERS From 33375dfa30f80de5fe69dcfb47425d21642e479a Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 30 Dec 2015 15:57:44 -0500 Subject: [PATCH 07/18] Small fixes for #538 --- openmc/filter.py | 2 +- openmc/geometry.py | 5 +++-- openmc/universe.py | 22 +++++++++++----------- src/initialize.F90 | 14 +++++++------- 4 files changed, 22 insertions(+), 21 deletions(-) diff --git a/openmc/filter.py b/openmc/filter.py index ea92f5257c..d61fb648f6 100644 --- a/openmc/filter.py +++ b/openmc/filter.py @@ -609,7 +609,7 @@ class Filter(object): # If this region is in Cell corresponding to the # distribcell filter bin, store it in dictionary if cell_id == self.bins[0]: - offset = openmc_geometry.get_offset(path) + offset = openmc_geometry.get_cell_instance(path) offsets_to_coords[offset] = coords # Each distribcell offset is a DataFrame bin diff --git a/openmc/geometry.py b/openmc/geometry.py index 16f6a45e5a..baadfce3c2 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -42,7 +42,7 @@ class Geometry(object): self._root_universe = root_universe - def get_offset(self, path): + def get_cell_instance(self, path): """Returns the corresponding location in the results array for a given path and filter number. This is primarily intended to post-processing result when a distribcell filter is used. @@ -77,7 +77,8 @@ class Geometry(object): # Begin recursive call to compute offset starting with the base Universe else: - offset = self._root_universe.get_offset(path, distribcell_ind) + offset = self._root_universe.get_cell_instance(path, + distribcell_ind) self._offsets[(path, distribcell_ind)] = offset # Return the final offset diff --git a/openmc/universe.py b/openmc/universe.py index b346ffd00d..e7b2151a32 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -285,7 +285,7 @@ class Cell(object): else: self.region = Intersection(self.region, region) - def get_offset(self, path, distribcell_ind): + def get_cell_instance(self, path, distribcell_ind): # Get the current element and remove it from the list cell_id = path[0] path = path[1:] @@ -296,12 +296,12 @@ class Cell(object): # If the Cell is filled by a Universe elif self._type == 'fill': - offset = self._offsets[distribcell_ind-1] - offset += self._fill.get_offset(path, distribcell_ind) + offset = self.offsets[distribcell_ind-1] + offset += self.fill.get_cell_instance(path, distribcell_ind) # If the Cell is filled by a Lattice else: - offset = self._fill.get_offset(path, distribcell_ind) + offset = self.fill.get_cell_instance(path, distribcell_ind) return offset @@ -605,7 +605,7 @@ class Universe(object): self._cells.clear() - def get_offset(self, path, distribcell_ind): + def get_cell_instance(self, path, distribcell_ind): # Get the current element and remove it from the list path = path[1:] @@ -613,7 +613,7 @@ class Universe(object): cell_id = path[0] # Make a recursive call to the Cell within this Universe - offset = self._cells[cell_id].get_offset(path, distribcell_ind) + offset = self.cells[cell_id].get_cell_instance(path, distribcell_ind) # Return the offset computed at all nested Universe levels return offset @@ -1073,7 +1073,7 @@ class RectLattice(Lattice): cv.check_greater_than('lattice pitch', dim, 0.0) self._pitch = pitch - def get_offset(self, path, distribcell_ind): + def get_cell_instance(self, path, distribcell_ind): # Get the current element and remove it from the list i = path[0] path = path[1:] @@ -1081,14 +1081,14 @@ class RectLattice(Lattice): # For 2D Lattices if len(self._dimension) == 2: offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_ind-1] - offset += self._universes[i[1]][i[2]].get_offset(path, - distribcell_ind) + offset += self._universes[i[1]][i[2]].get_cell_instance(path, + distribcell_ind) # For 3D Lattices else: offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_ind-1] - offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_offset(path, - distribcell_ind) + offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance( + path, distribcell_ind) return offset diff --git a/src/initialize.F90 b/src/initialize.F90 index acfe92b921..1134040e9a 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1086,16 +1086,16 @@ contains found(:,:) = .false. k = 1 + ! Search through universes for distributed cells and assign each one a + ! unique distribcell array index. do i = 1, n_universes univ => universes(i) - do j = 1, univ % n_cells - if (.not. cell_list % contains(univ % cells(j))) cycle - - cells(univ % cells(j)) % distribcell_ind = k - - univ_list(k) = univ % id - k = k + 1 + if (cell_list % contains(univ % cells(j))) then + cells(univ % cells(j)) % distribcell_ind = k + univ_list(k) = univ % id + k = k + 1 + end if end do end do From 562739c24554576997dad5d4172dd14764460213 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Wed, 30 Dec 2015 17:09:46 -0500 Subject: [PATCH 08/18] HDF5 summary file now report the lattice universes in original input ordering --- openmc/opencg_compatible.py | 7 +++++-- openmc/summary.py | 22 +++++++--------------- src/summary.F90 | 5 +++-- 3 files changed, 15 insertions(+), 19 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 62192355de..ad2cd06f83 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -725,11 +725,11 @@ def get_openmc_cell(opencg_cell): else: openmc_cell.fill = get_openmc_material(fill) - if opencg_cell.rotation: + if opencg_cell.rotation is not None: rotation = np.asarray(opencg_cell.rotation, dtype=np.float64) openmc_cell.rotation = rotation - if opencg_cell.translation: + if opencg_cell.translation is not None: translation = np.asarray(opencg_cell.translation, dtype=np.float64) openmc_cell.translation = translation @@ -908,6 +908,9 @@ def get_opencg_lattice(openmc_lattice): universe_id = universes[z][dimension[1]-y-1][x].id universe_array[z][y][x] = unique_universes[universe_id] + # Reverse y-dimension in array to match ordering in OpenCG + universe_array = universe_array[:, ::-1, :] + opencg_lattice = opencg.Lattice(lattice_id, name) opencg_lattice.dimension = dimension opencg_lattice.width = pitch diff --git a/openmc/summary.py b/openmc/summary.py index 3d0e37ddad..1c4d39b39b 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -329,11 +329,8 @@ class Summary(object): self._f['geometry/lattices'][key]['lower_left'][...] pitch = self._f['geometry/lattices'][key]['pitch'][...] outer = self._f['geometry/lattices'][key]['outer'].value - universe_ids = \ - self._f['geometry/lattices'][key]['universes'][...] - universe_ids = np.swapaxes(universe_ids, 0, 1) - universe_ids = np.swapaxes(universe_ids, 1, 2) + self._f['geometry/lattices'][key]['universes'][...] # Create the Lattice lattice = openmc.RectLattice(lattice_id=lattice_id, name=name) @@ -349,22 +346,17 @@ class Summary(object): universes = \ np.ndarray(tuple(universe_ids.shape), dtype=openmc.Universe) - for x in range(universe_ids.shape[0]): + for z in range(universe_ids.shape[0]): for y in range(universe_ids.shape[1]): - for z in range(universe_ids.shape[2]): - universes[x, y, z] = \ - self.get_universe_by_id(universe_ids[x, y, z]) + for x in range(universe_ids.shape[2]): + universes[z, y, x] = \ + self.get_universe_by_id(universe_ids[z, y, x]) - # Transpose, reverse y-dimension for appropriate ordering - shape = universes.shape[::-1] - universes = np.transpose(universes, (2, 1, 0)) - universes.shape = shape - universes = universes[:, ::-1, :] + # Set the universes for the lattice lattice.universes = universes if offsets is not None: - offsets = np.swapaxes(offsets, 0, 1) - offsets = np.swapaxes(offsets, 1, 2) + offsets = np.swapaxes(offsets, 0, 2) lattice.offsets = offsets # Add the Lattice to the global dictionary of all Lattices diff --git a/src/summary.F90 b/src/summary.F90 index cb005d3be5..f525797e24 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -362,9 +362,10 @@ contains allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), & &lat%n_cells(3))) do j = 1, lat%n_cells(1) - do k = 1, lat%n_cells(2) + do k = 0, lat%n_cells(2) - 1 do m = 1, lat%n_cells(3) - lattice_universes(j,k,m) = universes(lat%universes(j,k,m))%id + lattice_universes(j, k+1, m) = & + universes(lat%universes(j, lat%n_cells(2) - k, m))%id end do end do end do From 665f824d630a0160341b3fb79c8caf889ba723ab Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Fri, 1 Jan 2016 16:08:07 -0500 Subject: [PATCH 09/18] Removed unnecessary double flip of y lattice ordering in OpenCG compatibility module --- openmc/opencg_compatible.py | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index ad2cd06f83..9b86042ddc 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -905,12 +905,9 @@ def get_opencg_lattice(openmc_lattice): for z in range(dimension[2]): for y in range(dimension[1]): for x in range(dimension[0]): - universe_id = universes[z][dimension[1]-y-1][x].id + universe_id = universes[z][y][x].id universe_array[z][y][x] = unique_universes[universe_id] - # Reverse y-dimension in array to match ordering in OpenCG - universe_array = universe_array[:, ::-1, :] - opencg_lattice = opencg.Lattice(lattice_id, name) opencg_lattice.dimension = dimension opencg_lattice.width = pitch From 8daac97bc3d51d074b29b925b749390c2f1e497c Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Fri, 1 Jan 2016 16:43:12 -0500 Subject: [PATCH 10/18] OpenCG compatiblity module now properly deals with 3D lattices --- openmc/opencg_compatible.py | 17 ++++++++++++++--- openmc/summary.py | 7 +++++++ src/summary.F90 | 9 +++++++-- 3 files changed, 28 insertions(+), 5 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index 9b86042ddc..d778f44d85 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -881,16 +881,27 @@ def get_opencg_lattice(openmc_lattice): universes = openmc_lattice.universes outer = openmc_lattice.outer + # Convert 2D dimension to 3D for OpenCG + if len(dimension) == 2: + new_dimension = np.ones(3, dtype=np.int) + new_dimension[:2] = dimension + dimension = new_dimension + + # Convert 2D pitch to 3D for OpenCG if len(pitch) == 2: - new_pitch = np.ones(3, dtype=np.float64) * np.inf + new_pitch = np.ones(3, dtype=np.float64) * np.finfo(np.float64).max new_pitch[:2] = pitch pitch = new_pitch + # Convert 2D lower left to 3D for OpenCG if len(lower_left) == 2: - new_lower_left = np.ones(3, dtype=np.float64) + new_lower_left = np.ones(3, dtype=np.float64) * np.finfo(np.float64).min new_lower_left[:2] = lower_left lower_left = new_lower_left + # Convert 2D universes array to 3D for OpenCG + universes = np.atleast_3d(universes) + # Initialize an empty array for the OpenCG nested Universes in this Lattice universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), dtype=opencg.Universe) @@ -905,7 +916,7 @@ def get_opencg_lattice(openmc_lattice): for z in range(dimension[2]): for y in range(dimension[1]): for x in range(dimension[0]): - universe_id = universes[z][y][x].id + universe_id = universes[x][y][z].id universe_array[z][y][x] = unique_universes[universe_id] opencg_lattice = opencg.Lattice(lattice_id, name) diff --git a/openmc/summary.py b/openmc/summary.py index 1c4d39b39b..5405ac0e1e 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -352,6 +352,13 @@ class Summary(object): universes[z, y, x] = \ self.get_universe_by_id(universe_ids[z, y, x]) + # Use 2D NumPy array to store lattice universes for 2D lattices + if len(dimension) == 2: + print('squeezing!') + universes = np.squeeze(universes) + universes = np.atleast_2d(universes) + print(universes.shape) + # Set the universes for the lattice lattice.universes = universes diff --git a/src/summary.F90 b/src/summary.F90 index f525797e24..0f64f80ec2 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -355,8 +355,13 @@ contains call write_dataset(lattice_group, "type", "rectangular") ! Write lattice dimensions, lower left corner, and pitch - call write_dataset(lattice_group, "dimension", lat%n_cells) - call write_dataset(lattice_group, "lower_left", lat%lower_left) + if (lat % is_3d) then + call write_dataset(lattice_group, "dimension", lat % n_cells) + call write_dataset(lattice_group, "lower_left", lat % lower_left) + else + call write_dataset(lattice_group, "dimension", lat % n_cells(1:2)) + call write_dataset(lattice_group, "lower_left", lat % lower_left) + end if ! Write lattice universes. allocate(lattice_universes(lat%n_cells(1), lat%n_cells(2), & From 6efaf427619e5ecad51cff885aca50a45b0a3d8e Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sat, 2 Jan 2016 08:49:57 -0500 Subject: [PATCH 11/18] Fixed OpenCG compatibility from 2D to 3D lattice reshaping --- openmc/opencg_compatible.py | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/openmc/opencg_compatible.py b/openmc/opencg_compatible.py index d778f44d85..61485b5099 100644 --- a/openmc/opencg_compatible.py +++ b/openmc/opencg_compatible.py @@ -900,7 +900,8 @@ def get_opencg_lattice(openmc_lattice): lower_left = new_lower_left # Convert 2D universes array to 3D for OpenCG - universes = np.atleast_3d(universes) + if len(universes.shape) == 2: + universes.shape = (1,) + universes.shape # Initialize an empty array for the OpenCG nested Universes in this Lattice universe_array = np.ndarray(tuple(np.array(dimension)[::-1]), @@ -916,7 +917,7 @@ def get_opencg_lattice(openmc_lattice): for z in range(dimension[2]): for y in range(dimension[1]): for x in range(dimension[0]): - universe_id = universes[x][y][z].id + universe_id = universes[z][y][x].id universe_array[z][y][x] = unique_universes[universe_id] opencg_lattice = opencg.Lattice(lattice_id, name) From 0dddf0db56477346bde7210a2749084474c0812a Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sat, 2 Jan 2016 13:29:26 -0500 Subject: [PATCH 12/18] Now enforce lattice universes storage as a NumPy array --- openmc/universe.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmc/universe.py b/openmc/universe.py index 14a3007afb..60cbf4fa79 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -807,7 +807,7 @@ class Lattice(object): def universes(self, universes): cv.check_iterable_type('lattice universes', universes, Universe, min_depth=2, max_depth=3) - self._universes = universes + self._universes = np.asarray(universes) def get_unique_universes(self): """Determine all unique universes in the lattice From c48af6f71aad7829a426f8e06eb2757ad626ff88 Mon Sep 17 00:00:00 2001 From: "wbinventor@gmail.com" Date: Sat, 2 Jan 2016 17:36:07 -0500 Subject: [PATCH 13/18] Removed debug print statements from summary.py --- openmc/summary.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/openmc/summary.py b/openmc/summary.py index 5405ac0e1e..71074c8867 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -354,10 +354,8 @@ class Summary(object): # Use 2D NumPy array to store lattice universes for 2D lattices if len(dimension) == 2: - print('squeezing!') universes = np.squeeze(universes) universes = np.atleast_2d(universes) - print(universes.shape) # Set the universes for the lattice lattice.universes = universes From 7fea70650f576674f7663fbbbdee2f543956293b Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Sun, 3 Jan 2016 13:44:42 -0500 Subject: [PATCH 14/18] Improve geometry.get_cell_instance docstring --- openmc/geometry.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/openmc/geometry.py b/openmc/geometry.py index baadfce3c2..eb2b9db558 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -43,9 +43,9 @@ class Geometry(object): self._root_universe = root_universe def get_cell_instance(self, path): - """Returns the corresponding location in the results array for a given path and - filter number. This is primarily intended to post-processing result when - a distribcell filter is used. + """Return the instance number for the final cell in a geometry path. + + The instance is an index into tally distribcell filter arrays. Parameters ---------- @@ -58,8 +58,8 @@ class Geometry(object): Returns ------- - offset : int - Location in the results array for the path and filter + instance : int + Index in tally results array for distribcell filters """ From 493c6fd9ccb9d1e23516ba83b24249fb1d4ccba7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 4 Jan 2016 06:54:35 -0600 Subject: [PATCH 15/18] Fix type of return argument for h5pget_driver_f. --- src/hdf5_interface.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index fc7a462e6a..01e50d983a 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -1898,7 +1898,7 @@ contains logical :: mpio integer :: hdf5_err - integer :: driver + integer(HID_T) :: driver integer(HID_T) :: file_id integer(HID_T) :: fapl_id From 9f34183d49d8fd8a1e2f30caf49029f0ba115cba Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 4 Jan 2016 10:40:14 -0600 Subject: [PATCH 16/18] Revert to having Reaction%multiplicity_E as a pointer. Apparently having multiplicity_E as an allocatable scalar causes gfortran 5.2 to segfault when compiling ace.F90. Not sure why exactly... --- src/ace_header.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/ace_header.F90 b/src/ace_header.F90 index 985371ff18..d7c298979c 100644 --- a/src/ace_header.F90 +++ b/src/ace_header.F90 @@ -48,7 +48,7 @@ module ace_header integer :: MT ! ENDF MT value real(8) :: Q_value ! Reaction Q value integer :: multiplicity ! Number of secondary particles released - type(Tab1), allocatable :: multiplicity_E ! Energy-dependent neutron yield + type(Tab1), pointer :: multiplicity_E => null() ! Energy-dependent neutron yield integer :: threshold ! Energy grid index of threshold logical :: scatter_in_cm ! scattering system in center-of-mass? logical :: multiplicity_with_E = .false. ! Flag to indicate E-dependent multiplicity @@ -308,6 +308,8 @@ module ace_header class(Reaction), intent(inout) :: this ! The Reaction object to clear + if (associated(this % multiplicity_E)) deallocate(this % multiplicity_E) + if (associated(this % edist)) then call this % edist % clear() deallocate(this % edist) From 8730a243c3d1fb6102f335965b18179156f10d7c Mon Sep 17 00:00:00 2001 From: Sam Shaner Date: Mon, 4 Jan 2016 12:47:43 -0800 Subject: [PATCH 17/18] fixed typos in comments and descrptions in Python API --- openmc/statepoint.py | 2 +- openmc/summary.py | 2 +- openmc/tallies.py | 4 +++- 3 files changed, 5 insertions(+), 3 deletions(-) diff --git a/openmc/statepoint.py b/openmc/statepoint.py index 845937f594..b3cb4b6bdf 100644 --- a/openmc/statepoint.py +++ b/openmc/statepoint.py @@ -404,7 +404,7 @@ class StatePoint(object): 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.) + # Read scattering moment order strings (e.g., P3, Y1,2, etc.) moments = self._f['{0}{1}/moment_orders'.format( base, tally_key)].value diff --git a/openmc/summary.py b/openmc/summary.py index 00cf2d47ea..a347b5baa5 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -518,7 +518,7 @@ class Summary(object): # Create Tally object and assign basic properties tally = openmc.Tally(tally_id, tally_name) - # Read scattering moment order strings (e.g., P3, Y-1,2, etc.) + # Read scattering moment order strings (e.g., P3, Y1,2, etc.) moments = self._f['{0}/moment_orders'.format(subbase)].value # Read score metadata diff --git a/openmc/tallies.py b/openmc/tallies.py index cdf6c1f87e..190f147279 100644 --- a/openmc/tallies.py +++ b/openmc/tallies.py @@ -66,7 +66,9 @@ class Tally(object): triggers : list of openmc.trigger.Trigger List of tally triggers num_scores : Integral - Total number of user-specified scores + Total number of scores, accounting for the fact that a single + user-specified score, e.g. scatter-P3 or flux-Y2,2, might have multiple + bins num_filter_bins : Integral Total number of filter bins accounting for all filters num_bins : Integral From bae0ad713683bff2e6bc2bc234c0ec91f66709eb Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Mon, 4 Jan 2016 18:25:59 -0500 Subject: [PATCH 18/18] Small fixes for #538 --- docs/source/usersguide/output/statepoint.rst | 4 --- docs/source/usersguide/output/summary.rst | 4 +++ openmc/geometry.py | 10 +++--- openmc/summary.py | 4 +-- openmc/universe.py | 38 ++++++++++---------- src/geometry_header.F90 | 2 +- src/initialize.F90 | 2 +- src/input_xml.F90 | 2 +- src/output.F90 | 2 +- src/summary.F90 | 2 +- src/tally.F90 | 9 ++--- 11 files changed, 40 insertions(+), 39 deletions(-) diff --git a/docs/source/usersguide/output/statepoint.rst b/docs/source/usersguide/output/statepoint.rst index d3c1729af6..15bc79f739 100644 --- a/docs/source/usersguide/output/statepoint.rst +++ b/docs/source/usersguide/output/statepoint.rst @@ -185,10 +185,6 @@ if run_mode == 'k-eigenvalue': Type of the j-th filter. Can be 'universe', 'material', 'cell', 'cellborn', 'surface', 'mesh', 'energy', 'energyout', or 'distribcell'. -**/tallies/tally /filter /offset** (*int*) - - Filter offset (used for distribcell filter). - **/tallies/tally /filter /n_bins** (*int*) Number of bins for the j-th filter. diff --git a/docs/source/usersguide/output/summary.rst b/docs/source/usersguide/output/summary.rst index f87f60c4a0..302f340799 100644 --- a/docs/source/usersguide/output/summary.rst +++ b/docs/source/usersguide/output/summary.rst @@ -121,6 +121,10 @@ The current revision of the summary file format is 1. Region specification for the cell. +**/geometry/cells/cell /distribcell_index** (*int*) + + Index of this cell in distribcell filter arrays. + **/geometry/surfaces/surface /index** (*int*) Index in surfaces array used internally in OpenMC. diff --git a/openmc/geometry.py b/openmc/geometry.py index eb2b9db558..fc8e19f071 100644 --- a/openmc/geometry.py +++ b/openmc/geometry.py @@ -66,20 +66,20 @@ class Geometry(object): # Find the distribcell index of the cell. cells = self.get_all_cells() if path[-1] in cells: - distribcell_ind = cells[path[-1]].distribcell_ind + distribcell_index = cells[path[-1]].distribcell_index else: raise RuntimeError('Could not find cell {} specified in a \ distribcell filter'.format(path[-1])) # Return memoize'd offset if possible - if (path, distribcell_ind) in self._offsets: - offset = self._offsets[(path, distribcell_ind)] + if (path, distribcell_index) in self._offsets: + offset = self._offsets[(path, distribcell_index)] # Begin recursive call to compute offset starting with the base Universe else: offset = self._root_universe.get_cell_instance(path, - distribcell_ind) - self._offsets[(path, distribcell_ind)] = offset + distribcell_index) + self._offsets[(path, distribcell_index)] = offset # Return the final offset return offset diff --git a/openmc/summary.py b/openmc/summary.py index 9e2145b343..1c4c69aab7 100644 --- a/openmc/summary.py +++ b/openmc/summary.py @@ -272,9 +272,9 @@ class Summary(object): region, {s.id: s for s in self.surfaces.values()}) # Get the distribcell index - ind = self._f['geometry/cells'][key]['distribcell_ind'].value + ind = self._f['geometry/cells'][key]['distribcell_index'].value if ind != 0: - cell.distribcell_ind = ind + cell.distribcell_index = ind # Add the Cell to the global dictionary of all Cells self.cells[index] = cell diff --git a/openmc/universe.py b/openmc/universe.py index e7b2151a32..8c50ddbc75 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -63,7 +63,7 @@ class Cell(object): that is used to translate (shift) the universe. offsets : ndarray Array of offsets used for distributed cell searches - distribcell_ind : int + distribcell_index : int Index of this cell in distribcell arrays """ @@ -78,7 +78,7 @@ class Cell(object): self._rotation = None self._translation = None self._offsets = None - self._distribcell_ind = None + self._distribcell_index = None def __eq__(self, other): if not isinstance(other, Cell): @@ -126,7 +126,7 @@ class Cell(object): self._translation) string += '{0: <16}{1}{2}\n'.format('\tOffset', '=\t', self._offsets) string += '{0: <16}{1}{2}\n'.format('\tDistribcell index', '=\t', - self._distribcell_ind) + self._distribcell_index) return string @@ -170,8 +170,8 @@ class Cell(object): return self._offsets @property - def distribcell_ind(self): - return self._distribcell_ind + def distribcell_index(self): + return self._distribcell_index @id.setter def id(self, cell_id): @@ -240,10 +240,10 @@ class Cell(object): cv.check_type('cell region', region, Region) self._region = region - @distribcell_ind.setter - def distribcell_ind(self, ind): + @distribcell_index.setter + def distribcell_index(self, ind): cv.check_type('distribcell index', ind, Integral) - self._distribcell_ind = ind + self._distribcell_index = ind def add_surface(self, surface, halfspace): """Add a half-space to the list of half-spaces whose intersection defines the @@ -285,7 +285,7 @@ class Cell(object): else: self.region = Intersection(self.region, region) - def get_cell_instance(self, path, distribcell_ind): + def get_cell_instance(self, path, distribcell_index): # Get the current element and remove it from the list cell_id = path[0] path = path[1:] @@ -296,12 +296,12 @@ class Cell(object): # If the Cell is filled by a Universe elif self._type == 'fill': - offset = self.offsets[distribcell_ind-1] - offset += self.fill.get_cell_instance(path, distribcell_ind) + offset = self.offsets[distribcell_index-1] + offset += self.fill.get_cell_instance(path, distribcell_index) # If the Cell is filled by a Lattice else: - offset = self.fill.get_cell_instance(path, distribcell_ind) + offset = self.fill.get_cell_instance(path, distribcell_index) return offset @@ -605,7 +605,7 @@ class Universe(object): self._cells.clear() - def get_cell_instance(self, path, distribcell_ind): + def get_cell_instance(self, path, distribcell_index): # Get the current element and remove it from the list path = path[1:] @@ -613,7 +613,7 @@ class Universe(object): cell_id = path[0] # Make a recursive call to the Cell within this Universe - offset = self.cells[cell_id].get_cell_instance(path, distribcell_ind) + offset = self.cells[cell_id].get_cell_instance(path, distribcell_index) # Return the offset computed at all nested Universe levels return offset @@ -1073,22 +1073,22 @@ class RectLattice(Lattice): cv.check_greater_than('lattice pitch', dim, 0.0) self._pitch = pitch - def get_cell_instance(self, path, distribcell_ind): + def get_cell_instance(self, path, distribcell_index): # Get the current element and remove it from the list i = path[0] path = path[1:] # For 2D Lattices if len(self._dimension) == 2: - offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_ind-1] + offset = self._offsets[i[1]-1, i[2]-1, 0, distribcell_index-1] offset += self._universes[i[1]][i[2]].get_cell_instance(path, - distribcell_ind) + distribcell_index) # For 3D Lattices else: - offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_ind-1] + offset = self._offsets[i[1]-1, i[2]-1, i[3]-1, distribcell_index-1] offset += self._universes[i[1]-1][i[2]-1][i[3]-1].get_cell_instance( - path, distribcell_ind) + path, distribcell_index) return offset diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 684de883c2..3dad3ba395 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -136,7 +136,7 @@ module geometry_header ! expression logical :: simple ! Is the region simple (intersections ! only) - integer :: distribcell_ind ! Index corresponding to this cell in + integer :: distribcell_index ! Index corresponding to this cell in ! distribcell arrays ! Rotation matrix and translation vector diff --git a/src/initialize.F90 b/src/initialize.F90 index 1134040e9a..30654be9d2 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -1092,7 +1092,7 @@ contains univ => universes(i) do j = 1, univ % n_cells if (cell_list % contains(univ % cells(j))) then - cells(univ % cells(j)) % distribcell_ind = k + cells(univ % cells(j)) % distribcell_index = k univ_list(k) = univ % id k = k + 1 end if diff --git a/src/input_xml.F90 b/src/input_xml.F90 index a4317ed8e3..07d739b889 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1051,7 +1051,7 @@ contains ! Initialize distribcell instances and distribcell index c % instances = 0 - c % distribcell_ind = NONE + c % distribcell_index = NONE ! Get pointer to i-th cell node call get_list_item(node_cell_list, i, node_cell) diff --git a/src/output.F90 b/src/output.F90 index 5b938e98d3..66e6831371 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -1436,7 +1436,7 @@ contains class(Lattice), pointer :: lat ! Pointer to current lattice ! Get the distribcell index for this cell - map = cells(goal) % distribcell_ind + map = cells(goal) % distribcell_index n = univ % n_cells diff --git a/src/summary.F90 b/src/summary.F90 index 11c4cd434b..7a430cf741 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -196,7 +196,7 @@ contains end do call write_dataset(cell_group, "region", adjustl(region_spec)) - call write_dataset(cell_group, "distribcell_ind", c % distribcell_ind) + call write_dataset(cell_group, "distribcell_index", c % distribcell_index) call close_group(cell_group) end do CELL_LOOP diff --git a/src/tally.F90 b/src/tally.F90 index 2ccd59a58e..7210efc5c8 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1712,7 +1712,7 @@ contains integer :: j integer :: n ! number of bins for single filter integer :: offset ! offset for distribcell - integer :: distribcell_ind ! index in distribcell arrays + integer :: distribcell_index ! index in distribcell arrays real(8) :: E ! particle energy real(8) :: theta, phi ! Polar and Azimuthal Angles, respectively type(TallyObject), pointer :: t @@ -1757,13 +1757,14 @@ contains case (FILTER_DISTRIBCELL) ! determine next distribcell bin - distribcell_ind = cells(t % filters(i) % int_bins(1)) % distribcell_ind + distribcell_index = cells(t % filters(i) % int_bins(1)) & + % distribcell_index matching_bins(i) = NO_BIN_FOUND offset = 0 do j = 1, p % n_coord if (cells(p % coord(j) % cell) % type == CELL_FILL) then offset = offset + cells(p % coord(j) % cell) % & - offset(distribcell_ind) + offset(distribcell_index) elseif(cells(p % coord(j) % cell) % type == CELL_LATTICE) then if (lattices(p % coord(j + 1) % lattice) % obj & % are_valid_indices([& @@ -1771,7 +1772,7 @@ contains p % coord(j + 1) % lattice_y, & p % coord(j + 1) % lattice_z])) then offset = offset + lattices(p % coord(j + 1) % lattice) % obj % & - offset(distribcell_ind, & + offset(distribcell_index, & p % coord(j + 1) % lattice_x, & p % coord(j + 1) % lattice_y, & p % coord(j + 1) % lattice_z)