Make sure distribcell paths match those in summary file

This commit is contained in:
Paul Romano 2017-03-06 08:54:35 -06:00
parent d22f10d60e
commit 9b7c2fd413
4 changed files with 189 additions and 177 deletions

View file

@ -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)

View file

@ -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:

View file

@ -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

View file

@ -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