From 9b7c2fd4137b08ee837d54b323229eb3e22e1e44 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Mar 2017 08:54:35 -0600 Subject: [PATCH] Make sure distribcell paths match those in summary file --- openmc/lattice.py | 127 +++++++++++++++++++++++------------------- openmc/universe.py | 10 ++-- src/geometry.F90 | 129 +++++++++++++++++++------------------------ src/tally_filter.F90 | 100 ++++++++++++++++++--------------- 4 files changed, 189 insertions(+), 177 deletions(-) diff --git a/openmc/lattice.py b/openmc/lattice.py index 298b5842a3..7994048b17 100644 --- a/openmc/lattice.py +++ b/openmc/lattice.py @@ -378,6 +378,51 @@ class Lattice(object): return all_universes + def get_universe(self, idx): + """Return universe corresponding to a lattice element index + + Parameters + ---------- + idx : Iterable of int + Lattice element indices in the :math:`(x,y,z)` coordinate system + + Returns + ------- + openmc.Universe + Universe with given indices + + """ + idx_u = self.get_universe_index(idx) + if self.ndim == 2: + return self.universes[idx_u[0]][idx_u[1]] + else: + return self.universes[idx_u[0]][idx_u[1]][idx_u[2]] + + def find(self, point): + """Find cells/universes/lattices which contain a given point + + Parameters + ---------- + point : 3-tuple of float + Cartesian coordinatesof the point + + Returns + ------- + list + Sequence of universes, cells, and lattices which are traversed to + find the given point + + """ + idx, p = self.find_element(point) + if self.is_valid_index(idx): + u = self.get_universe(idx) + else: + if self.outer is not None: + u = self.outer + else: + return [] + return [(self, idx)] + u.find(p) + class RectLattice(Lattice): """A lattice consisting of rectangular prisms. @@ -516,6 +561,15 @@ class RectLattice(Lattice): return list(np.broadcast(*np.ogrid[ :self.shape[2], :self.shape[1], :self.shape[0]])) + @property + def _natural_indices(self): + if self.ndim == 2: + nx, ny = self.shape + return np.broadcast(*np.ogrid[:nx, :ny]) + else: + nx, ny, nz = self.shape + return np.broadcast(*np.ogrid[:nx, :ny, :nz]) + @property def lower_left(self): return self._lower_left @@ -681,32 +735,6 @@ class RectLattice(Lattice): 0 <= idx[1] < self.shape[1] and 0 <= idx[2] < self.shape[2]) - def find(self, point): - """Find cells/universes/lattices which contain a given point - - Parameters - ---------- - point : 3-tuple of float - Cartesian coordinatesof the point - - Returns - ------- - list - Sequence of universes, cells, and lattices which are traversed to - find the given point - - """ - idx, p = self.find_element(point) - if self.is_valid_index(idx): - idx_u = self.get_universe_index(idx) - u = self.universes[idx_u] - else: - if self.outer is not None: - u = self.outer - else: - return [] - return [(self, idx)] + u.find(p) - def create_xml_subelement(self, xml_element): # Determine if XML element already contains subelement for this Lattice @@ -922,6 +950,23 @@ class HexLattice(Lattice): for r in range(self.num_rings) for i in range(max(6*(self.num_rings - 1 - r), 1))] + @property + def _natural_indices(self): + r = self.num_rings + if self.num_axial is None: + for a in range(-r + 1, r): + for x in range(-r + 1, r): + idx = (x, a) + if self.is_valid_index(idx): + yield idx + else: + for z in range(self.num_axial): + for a in range(-r + 1, r): + for x in range(-r + 1, r): + idx = (x, a, z) + if self.is_valid_index(idx): + yield idx + @property def ndim(self): return 2 if isinstance(self.universes[0][0], openmc.Universe) else 3 @@ -1144,36 +1189,6 @@ class HexLattice(Lattice): else: return g < self.num_rings and 0 <= idx[2] < self.num_axial - def find(self, point): - """Find cells/universes/lattices which contain a given point - - Parameters - ---------- - point : 3-tuple of float - Cartesian coordinatesof the point - - Returns - ------- - list - Sequence of universes, cells, and lattices which are traversed to - find the given point - - """ - idx, p = self.find_element(point) - if self.is_valid_index(idx): - idx_u = self.get_universe_index(idx) - if self.num_axial is None: - u = self.universes[idx_u[0]][idx_u[1]] - else: - u = self.universes[idx_u[0]][idx_u[1]][idx_u[2]] - else: - if self.outer is not None: - u = self.outer - else: - return [] - - return [(self, idx)] + u.find(p) - def create_xml_subelement(self, xml_element): # Determine if XML element already contains subelement for this Lattice path = './hex_lattice[@id=\'{0}\']'.format(self._id) diff --git a/openmc/universe.py b/openmc/universe.py index fd941a224e..3e9963e492 100644 --- a/openmc/universe.py +++ b/openmc/universe.py @@ -538,12 +538,10 @@ class Universe(object): latt = cell.fill # Count instances in each universe in the lattice - for index in latt.indices: - latt_path = '{}->l{}{}->'.format(cell_path, latt.id, index) - if latt.ndim == 3: - univ = latt.universes[index[0]][index[1]][index[2]] - else: - univ = latt.universes[index[0]][index[1]] + for index in latt._natural_indices: + latt_path = '{}->l{}({})->'.format( + cell_path, latt.id, ",".join(str(x) for x in index)) + univ = latt.get_universe(index) univ._determine_paths(latt_path) else: diff --git a/src/geometry.F90 b/src/geometry.F90 index df963bc2ff..7365f5af2f 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -1087,9 +1087,9 @@ contains ! routine is called once upon initialization. !=============================================================================== - subroutine calc_offsets(goal, map, univ, counts, found) + subroutine calc_offsets(univ_id, map, univ, counts, found) - integer, intent(in) :: goal ! target universe ID + integer, intent(in) :: univ_id ! target universe ID integer, intent(in) :: map ! map index in vector of maps type(Universe), intent(in) :: univ ! universe searching in integer, intent(inout) :: counts(:,:) ! target count @@ -1126,7 +1126,7 @@ contains ! Count contents of this cell next_univ => universes(c % fill) - offset = offset + count_target(next_univ, counts, found, goal, map) + offset = offset + count_target(next_univ, counts, found, univ_id, map) ! Move into the next universe next_univ => universes(c % fill) @@ -1150,7 +1150,7 @@ contains lat % offset(map, j, k, m) = offset next_univ => universes(lat % universes(j, k, m)) offset = offset + & - count_target(next_univ, counts, found, goal, map) + count_target(next_univ, counts, found, univ_id, map) end do end do end do @@ -1171,7 +1171,7 @@ contains lat % offset(map, j, k, m) = offset next_univ => universes(lat % universes(j, k, m)) offset = offset + & - count_target(next_univ, counts, found, goal, map) + count_target(next_univ, counts, found, univ_id, map) end if end do end do @@ -1188,12 +1188,12 @@ contains ! universe ID beginning with the universe given. !=============================================================================== - recursive function count_target(univ, counts, found, goal, map) result(count) + recursive function count_target(univ, counts, found, univ_id, map) result(count) type(Universe), intent(in) :: univ ! universe to search through integer, intent(inout) :: counts(:,:) ! target count logical, intent(inout) :: found(:,:) ! target found - integer, intent(in) :: goal ! target universe ID + integer, intent(in) :: univ_id ! target universe ID integer, intent(in) :: map ! current map integer :: i ! index over cells @@ -1213,7 +1213,7 @@ contains ! If this is the target, it can't contain itself. ! Count = 1, then quit - if (univ % id == goal) then + if (univ % id == univ_id) then count = 1 counts(universe_dict % get_key(univ % id), map) = 1 found(universe_dict % get_key(univ % id), map) = .true. @@ -1241,12 +1241,12 @@ contains next_univ => universes(c % fill) ! Found target - stop since target cannot contain itself - if (next_univ % id == goal) then + if (next_univ % id == univ_id) then count = count + 1 return end if - count = count + count_target(next_univ, counts, found, goal, map) + count = count + count_target(next_univ, counts, found, univ_id, map) c => cells(cell_index) ! ==================================================================== @@ -1267,13 +1267,13 @@ contains next_univ => universes(lat % universes(j, k, m)) ! Found target - stop since target cannot contain itself - if (next_univ % id == goal) then + if (next_univ % id == univ_id) then count = count + 1 cycle end if count = count + & - count_target(next_univ, counts, found, goal, map) + count_target(next_univ, counts, found, univ_id, map) end do end do @@ -1295,13 +1295,13 @@ contains next_univ => universes(lat % universes(j, k, m)) ! Found target - stop since target cannot contain itself - if (next_univ % id == goal) then + if (next_univ % id == univ_id) then count = count + 1 cycle end if count = count + & - count_target(next_univ, counts, found, goal, map) + count_target(next_univ, counts, found, univ_id, map) end if end do end do @@ -1326,81 +1326,66 @@ contains type(Universe), intent(in) :: univ ! universe to search through - integer :: i ! index over cells - integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search - integer :: cell_index ! index in cells array - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next universe to loop through - class(Lattice), pointer :: lat ! pointer to current lattice + integer :: i ! index over cells + integer :: j, k, m ! indices in lattice + integer :: n ! number of cells to search n = size(univ % cells) do i = 1, n - cell_index = univ % cells(i) - ! get pointer to cell - c => cells(cell_index) - c % instances = c % instances + 1 + associate (c => cells(univ % cells(i))) + c % instances = c % instances + 1 - ! ==================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type == FILL_MATERIAL) then + ! ==================================================================== + ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL + if (c % type == FILL_UNIVERSE) then - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type == FILL_UNIVERSE) then + call count_instance(universes(c % fill)) - next_univ => universes(c % fill) + ! ==================================================================== + ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL + elseif (c % type == FILL_LATTICE) then - call count_instance(next_univ) - c => cells(cell_index) + ! Set current lattice + associate (lat => lattices(c % fill) % obj) - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type == FILL_LATTICE) then + select type (lat) + type is (RectLattice) - ! Set current lattice - lat => lattices(c % fill) % obj - - select type (lat) - - type is (RectLattice) - - ! Loop over lattice coordinates - do j = 1, lat % n_cells(1) - do k = 1, lat % n_cells(2) - do m = 1, lat % n_cells(3) - next_univ => universes(lat % universes(j, k, m)) - call count_instance(next_univ) + ! Loop over lattice coordinates + do j = 1, lat % n_cells(1) + do k = 1, lat % n_cells(2) + do m = 1, lat % n_cells(3) + call count_instance(universes(lat % universes(j, k, m))) + end do + end do end do - end do - end do - type is (HexLattice) + type is (HexLattice) - ! Loop over lattice coordinates - do m = 1, lat % n_axial - do k = 1, 2*lat % n_rings - 1 - do j = 1, 2*lat % n_rings - 1 - ! This array location is never used - if (j + k < lat % n_rings + 1) then - cycle - ! This array location is never used - else if (j + k > 3*lat % n_rings - 1) then - cycle - else - next_univ => universes(lat % universes(j, k, m)) - call count_instance(next_univ) - end if + ! Loop over lattice coordinates + do m = 1, lat % n_axial + do k = 1, 2*lat % n_rings - 1 + do j = 1, 2*lat % n_rings - 1 + ! This array location is never used + if (j + k < lat % n_rings + 1) then + cycle + ! This array location is never used + else if (j + k > 3*lat % n_rings - 1) then + cycle + else + call count_instance(universes(lat % universes(j, k, m))) + end if + end do + end do end do - end do - end do - end select - - end if + end select + end associate + end if + end associate end do end subroutine count_instance diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index f88362be33..35340fa5c2 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -1370,11 +1370,11 @@ contains ! the target cell with the given offset !=============================================================================== - recursive subroutine find_offset(goal, univ, final, offset, path) + recursive subroutine find_offset(i_cell, univ, target_offset, offset, path) - integer, intent(in) :: goal ! The target cell index + integer, intent(in) :: i_cell ! The target cell index type(Universe), intent(in) :: univ ! Universe to begin search - integer, intent(in) :: final ! Target offset + integer, intent(in) :: target_offset ! Target offset integer, intent(inout) :: offset ! Current offset character(*), intent(inout) :: path ! Path to offset @@ -1395,27 +1395,23 @@ contains class(Lattice), pointer :: lat ! Pointer to current lattice ! Get the distribcell index for this cell - map = cells(goal) % distribcell_index + map = cells(i_cell) % distribcell_index n = size(univ % cells) ! Write to the geometry stack i_univ = universe_dict % get_key(univ % id) if (i_univ == root_universe) then - path = trim(path) // to_str(univ%id) + path = trim(path) // "u" // to_str(univ%id) else - path = trim(path) // "->" // to_str(univ%id) + path = trim(path) // "->u" // to_str(univ%id) end if ! Look through all cells in this universe do i = 1, n ! 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 + if (univ % cells(i) == i_cell .and. offset == target_offset) return end do ! Find the fill cell or lattice cell that we need to enter @@ -1423,8 +1419,7 @@ contains later_cell = .false. - cell_index = univ % cells(i) - c => cells(cell_index) + c => cells(univ % cells(i)) this_cell = .false. @@ -1435,13 +1430,10 @@ contains do j = i+1, n - cell_index = univ % cells(j) - c => cells(cell_index) + c => cells(univ % cells(j)) ! Skip normal cells which do not have offsets - if (c % type == FILL_MATERIAL) then - cycle - end if + if (c % type == FILL_MATERIAL) cycle ! Break loop once we've found the next cell with an offset exit @@ -1465,7 +1457,7 @@ contains ! If the final offset is in the range of offset - temp_offset+offset ! then the goal is in this cell - if (final < temp_offset + offset) then + if (target_offset < temp_offset + offset) then this_cell = .true. end if end if @@ -1485,7 +1477,7 @@ contains cell_index = univ % cells(i) c => cells(cell_index) - path = trim(path) // "->" // to_str(c%id) + path = trim(path) // "->c" // to_str(c%id) ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL @@ -1495,7 +1487,7 @@ contains offset = c % offset(map) + offset next_univ => universes(c % fill) - call find_offset(goal, next_univ, final, offset, path) + call find_offset(i_cell, next_univ, target_offset, offset, path) return ! ==================================================================== @@ -1512,7 +1504,7 @@ contains type is (RectLattice) ! Write to the geometry stack - path = trim(path) // "->" // to_str(lat%id) + path = trim(path) // "->l" // to_str(lat%id) n_x = lat % n_cells(1) n_y = lat % n_cells(2) @@ -1526,16 +1518,21 @@ contains do l = 1, n_y do m = 1, n_z - if (final >= lat % offset(map, k, l, m) + offset) then + if (target_offset >= lat % offset(map, k, l, m) + offset) then if (k == n_x .and. l == n_y .and. m == n_z) then ! This is last lattice cell, so target must be here lat_offset = lat % offset(map, k, l, m) offset = offset + lat_offset next_univ => universes(lat % universes(k, l, m)) - path = trim(path) // "(" // trim(to_str(k)) // & - "," // trim(to_str(l)) // "," // & - trim(to_str(m)) // ")" - call find_offset(goal, next_univ, final, offset, path) + if (lat % is_3d) then + path = trim(path) // "(" // trim(to_str(k-1)) // & + "," // trim(to_str(l-1)) // "," // & + trim(to_str(m-1)) // ")" + else + path = trim(path) // "(" // trim(to_str(k-1)) // & + "," // trim(to_str(l-1)) // ")" + end if + call find_offset(i_cell, next_univ, target_offset, offset, path) return else old_m = m @@ -1548,10 +1545,15 @@ contains lat_offset = lat % offset(map, old_k, old_l, old_m) offset = offset + lat_offset next_univ => universes(lat % universes(old_k, old_l, old_m)) - path = trim(path) // "(" // trim(to_str(old_k)) // & - "," // trim(to_str(old_l)) // "," // & - trim(to_str(old_m)) // ")" - call find_offset(goal, next_univ, final, offset, path) + if (lat % is_3d) then + path = trim(path) // "(" // trim(to_str(old_k-1)) // & + "," // trim(to_str(old_l-1)) // "," // & + trim(to_str(old_m-1)) // ")" + else + path = trim(path) // "(" // trim(to_str(old_k-1)) // & + "," // trim(to_str(old_l-1)) // ")" + end if + call find_offset(i_cell, next_univ, target_offset, offset, path) return end if @@ -1564,7 +1566,7 @@ contains type is (HexLattice) ! Write to the geometry stack - path = trim(path) // "->" // to_str(lat%id) + path = trim(path) // "->l" // to_str(lat%id) n_z = lat % n_axial n_y = 2 * lat % n_rings - 1 @@ -1586,17 +1588,23 @@ contains cycle end if - if (final >= lat % offset(map, k, l, m) + offset) then + if (target_offset >= lat % offset(map, k, l, m) + offset) then if (k == lat % n_rings .and. l == n_y .and. m == n_z) then ! This is last lattice cell, so target must be here lat_offset = lat % offset(map, k, l, m) offset = offset + lat_offset next_univ => universes(lat % universes(k, l, m)) - path = trim(path) // "(" // & - trim(to_str(k - lat % n_rings)) // "," // & - trim(to_str(l - lat % n_rings)) // "," // & - trim(to_str(m)) // ")" - call find_offset(goal, next_univ, final, offset, path) + if (lat % is_3d) then + path = trim(path) // "(" // & + trim(to_str(k - lat % n_rings)) // "," // & + trim(to_str(l - lat % n_rings)) // "," // & + trim(to_str(m - 1)) // ")" + else + path = trim(path) // "(" // & + trim(to_str(k - lat % n_rings)) // "," // & + trim(to_str(l - lat % n_rings)) // ")" + end if + call find_offset(i_cell, next_univ, target_offset, offset, path) return else old_m = m @@ -1609,11 +1617,17 @@ contains lat_offset = lat % offset(map, old_k, old_l, old_m) offset = offset + lat_offset next_univ => universes(lat % universes(old_k, old_l, old_m)) - path = trim(path) // "(" // & - trim(to_str(old_k - lat % n_rings)) // "," // & - trim(to_str(old_l - lat % n_rings)) // "," // & - trim(to_str(old_m)) // ")" - call find_offset(goal, next_univ, final, offset, path) + if (lat % is_3d) then + path = trim(path) // "(" // & + trim(to_str(old_k - lat % n_rings)) // "," // & + trim(to_str(old_l - lat % n_rings)) // "," // & + trim(to_str(old_m - 1)) // ")" + else + path = trim(path) // "(" // & + trim(to_str(old_k - lat % n_rings)) // "," // & + trim(to_str(old_l - lat % n_rings)) // ")" + end if + call find_offset(i_cell, next_univ, target_offset, offset, path) return end if