From 9bf2eaa42e559e57cbb82e02384920fec3134d67 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Thu, 24 May 2018 12:53:09 -0400 Subject: [PATCH] Simplify distribcell find_offset logic --- src/tallies/tally_filter_distribcell.F90 | 328 +++++++---------------- 1 file changed, 98 insertions(+), 230 deletions(-) diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index fa0a6b296..e121cae5e 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -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