Simplify distribcell find_offset logic

This commit is contained in:
Sterling Harper 2018-05-24 12:53:09 -04:00
parent 07ded4e7be
commit 9bf2eaa42e

View file

@ -139,17 +139,11 @@ contains
character(*), intent(inout) :: path ! Path to offset
integer :: map ! Index in maps vector
integer :: i, j ! Index over cells
integer :: i ! Index over cells
integer :: k, l, m ! Indices in lattice
integer :: old_k, old_l, old_m ! Previous indices in lattice
integer :: n_x, n_y, n_z ! Lattice cell array dimensions
integer :: n ! Number of cells to search
integer :: cell_index ! Index in cells array
integer :: lat_offset ! Offset from lattice
integer :: temp_offset ! Looped sum of offsets
integer :: i_univ ! index in universes array
logical :: this_cell = .false. ! Advance in this cell?
logical :: later_cell = .false. ! Fill cells after this one?
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
@ -157,8 +151,6 @@ contains
! Get the distribcell index for this cell
map = cells(i_cell) % distribcell_index
n = size(univ % cells)
! Write to the geometry stack
i_univ = universe_dict % get(univ % id)
if (i_univ == root_universe) then
@ -168,7 +160,7 @@ contains
end if
! Look through all cells in this universe
do i = 1, n
do i = 1, size(univ % cells)
! If the cell matches the goal and the offset matches final, write to the
! geometry stack
if (univ % cells(i) == i_cell .and. offset == target_offset) then
@ -178,234 +170,110 @@ contains
end if
end do
! Find the fill cell or lattice cell that we need to enter
do i = 1, n
later_cell = .false.
! Find the fill cell or lattice cell that contains the target offset.
do i = size(univ % cells), 2, -1
c => cells(univ % cells(i))
this_cell = .false.
! Material cells don't contain other cells so ignore them.
if (c % type() /= FILL_MATERIAL) then
! If we got here, we still think the target is in this universe
! or further down, but it's not this exact cell.
! Compare offset to next cell to see if we should enter this cell
if (i /= n) then
do j = i+1, n
c => cells(univ % cells(j))
! Skip normal cells which do not have offsets
if (c % type() == FILL_MATERIAL) cycle
! Break loop once we've found the next cell with an offset
exit
end do
! Ensure we didn't just end the loop by iteration
if (c % type() /= FILL_MATERIAL) then
! There are more cells in this universe that it could be in
later_cell = .true.
! Two cases, lattice or fill cell
if (c % type() == FILL_UNIVERSE) then
temp_offset = c % offset(map-1)
! Get the offset of the first lattice location
else
lat => lattices(c % fill() + 1) % obj
temp_offset = lat % offset(map-1, [0, 0, 0])
end if
! If the final offset is in the range of offset - temp_offset+offset
! then the goal is in this cell
if (target_offset < temp_offset + offset) then
this_cell = .true.
end if
end if
end if
if (n == 1 .and. c % type() /= FILL_MATERIAL) then
this_cell = .true.
end if
if (.not. later_cell) then
this_cell = .true.
end if
! Get pointer to THIS cell because target must be in this cell
if (this_cell) then
cell_index = univ % cells(i)
c => cells(cell_index)
path = trim(path) // "->c" // to_str(c%id())
! ====================================================================
! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL
if (c % type() == FILL_UNIVERSE) then
! Enter this cell to update the current offset
offset = c % offset(map-1) + offset
next_univ => universes(c % fill() + 1)
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
! ====================================================================
! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL
elseif (c % type() == FILL_LATTICE) then
! Set current lattice
temp_offset = offset + c % offset(map-1)
else
! Get the offset of the first lattice location
lat => lattices(c % fill() + 1) % obj
select type (lat)
! ==================================================================
! RECTANGULAR LATTICES
type is (RectLattice)
! Write to the geometry stack
path = trim(path) // "->l" // to_str(lat%id())
n_x = lat % n_cells(1)
n_y = lat % n_cells(2)
n_z = lat % n_cells(3)
old_m = 1
old_l = 1
old_k = 1
! Loop over lattice coordinates
do m = 1, n_z
do l = 1, n_y
do k = 1, n_x
if (target_offset >= lat % offset(map-1, [k-1, l-1, m-1]) &
+ 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-1, [k-1, l-1, m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
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
old_l = l
old_k = k
cycle
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map-1, [old_k-1, old_l-1, old_m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1)
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
end do
end do
end do
! ==================================================================
! HEXAGONAL LATTICES
type is (HexLattice)
! Write to the geometry stack
path = trim(path) // "->l" // to_str(lat%id())
n_z = lat % n_axial
n_y = 2 * lat % n_rings - 1
n_x = 2 * lat % n_rings - 1
old_m = 1
old_l = 1
old_k = 1
! Loop over lattice coordinates
do m = 1, n_z
do l = 1, n_y
do k = 1, n_x
! This array position is never used
if (k + l < lat % n_rings + 1) then
cycle
! This array position is never used
else if (k + l > 3*lat % n_rings - 1) then
cycle
end if
if (target_offset >= lat % offset(map-1, [k-1, l-1, m-1]) &
+ 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-1, [k-1, l-1, m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
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
old_l = l
old_k = k
cycle
end if
else
! Target is at this lattice position
lat_offset = lat % offset(map-1, [old_k-1, old_l-1, old_m-1])
offset = offset + lat_offset
next_univ => universes(lat % get([old_k-1, old_l-1, old_m-1])+1)
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
end do
end do
end do
end select
temp_offset = offset + lat % offset(map-1, [0, 0, 0])
end if
! The desired cell is the first cell that gives an offset smaller or
! equal to the target offset.
if (temp_offset <= target_offset) exit
end if
end do
! Add the cell to the path string.
c => cells(univ % cells(i))
path = trim(path) // "->c" // to_str(c%id())
if (c % type() == FILL_UNIVERSE) then
! Recurse into the fill cell.
offset = c % offset(map-1) + offset
next_univ => universes(c % fill() + 1)
call find_offset(i_cell, next_univ, target_offset, offset, path)
elseif (c % type() == FILL_LATTICE) then
! Recurse into the lattice cell.
lat => lattices(c % fill() + 1) % obj
select type (lat)
type is (RectLattice)
path = trim(path) // "->l" // to_str(lat%id())
n_x = lat % n_cells(1)
n_y = lat % n_cells(2)
n_z = lat % n_cells(3)
do m = n_z, 1, -1
do l = n_y, 1, -1
do k = n_x, 1, -1
temp_offset = offset + lat % offset(map-1, [k-1, l-1, m-1])
if (temp_offset <= target_offset) then
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
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
offset = temp_offset
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
end if
end do
end do
end do
type is (HexLattice)
path = trim(path) // "->l" // to_str(lat%id())
n_z = lat % n_axial
n_y = 2 * lat % n_rings - 1
n_x = 2 * lat % n_rings - 1
do m = n_z, 1, -1
do l = n_y, 1, -1
do k = n_x, 1, -1
if (k + l < lat % n_rings + 1) then
cycle
else if (k + l > 3*lat % n_rings - 1) then
cycle
end if
temp_offset = offset + lat % offset(map-1, [k-1, l-1, m-1])
if (temp_offset <= target_offset) then
next_univ => universes(lat % get([k-1, l-1, m-1])+1)
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
offset = temp_offset
call find_offset(i_cell, next_univ, target_offset, offset, path)
return
end if
end do
end do
end do
end select
end if
end subroutine find_offset
end module tally_filter_distribcell