diff --git a/CMakeLists.txt b/CMakeLists.txt index 669494e059..b0bcd1b1b3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -378,6 +378,7 @@ set(LIBOPENMC_FORTRAN_SRC src/tallies/tally.F90 src/tallies/tally_filter.F90 src/tallies/tally_filter_header.F90 + src/tallies/tally_filter_distribcell.F90 src/tallies/tally_filter_energy.F90 src/tallies/tally_filter_material.F90 src/tallies/tally_filter_mesh.F90 diff --git a/src/summary.F90 b/src/summary.F90 index e4263d8d0a..6ee8b2e781 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -16,7 +16,7 @@ module summary use surface_header use string, only: to_str use tally_header, only: TallyObject - use tally_filter, only: find_offset + use tally_filter_distribcell, only: find_offset implicit none private diff --git a/src/tallies/tally_filter.F90 b/src/tallies/tally_filter.F90 index d48939f8a6..dc3aa36b18 100644 --- a/src/tallies/tally_filter.F90 +++ b/src/tallies/tally_filter.F90 @@ -14,6 +14,7 @@ module tally_filter use tally_filter_header ! Inherit other filters + use tally_filter_distribcell use tally_filter_energy use tally_filter_material use tally_filter_mesh @@ -57,19 +58,6 @@ module tally_filter procedure :: text_label => text_label_cell_from end type CellFromFilter -!=============================================================================== -! DISTRIBCELLFILTER specifies which distributed geometric cells tally events -! reside in. -!=============================================================================== - type, extends(TallyFilter) :: DistribcellFilter - integer :: cell - contains - procedure :: get_all_bins => get_all_bins_distribcell - procedure :: to_statepoint => to_statepoint_distribcell - procedure :: text_label => text_label_distribcell - procedure :: initialize => initialize_distribcell - end type DistribcellFilter - !=============================================================================== ! CELLBORNFILTER specifies which cell the particle was born in. !=============================================================================== @@ -352,83 +340,6 @@ contains label = "Cell from " // to_str(cells(this % cells(bin)) % id) end function text_label_cell_from -!=============================================================================== -! DistribcellFilter methods -!=============================================================================== - subroutine get_all_bins_distribcell(this, p, estimator, match) - class(DistribcellFilter), intent(in) :: this - type(Particle), intent(in) :: p - integer, intent(in) :: estimator - type(TallyFilterMatch), intent(inout) :: match - - integer :: distribcell_index, offset, i - - distribcell_index = cells(this % cell) % distribcell_index - offset = 0 - do i = 1, p % n_coord - if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then - offset = offset + cells(p % coord(i) % cell) % & - offset(distribcell_index) - elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then - if (lattices(p % coord(i + 1) % lattice) % obj & - % are_valid_indices([& - p % coord(i + 1) % lattice_x, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z])) then - offset = offset + lattices(p % coord(i + 1) % lattice) % obj % & - offset(distribcell_index, & - p % coord(i + 1) % lattice_x, & - p % coord(i + 1) % lattice_y, & - p % coord(i + 1) % lattice_z) - end if - end if - if (this % cell == p % coord(i) % cell) then - call match % bins % push_back(offset + 1) - call match % weights % push_back(ONE) - return - end if - end do - end subroutine get_all_bins_distribcell - - subroutine to_statepoint_distribcell(this, filter_group) - class(DistribcellFilter), intent(in) :: this - integer(HID_T), intent(in) :: filter_group - - call write_dataset(filter_group, "type", "distribcell") - call write_dataset(filter_group, "n_bins", this % n_bins) - call write_dataset(filter_group, "bins", cells(this % cell) % id) - end subroutine to_statepoint_distribcell - - subroutine initialize_distribcell(this) - class(DistribcellFilter), intent(inout) :: this - - integer :: id - - ! Convert id to index. - id = this % cell - if (cell_dict % has_key(id)) then - this % cell = cell_dict % get_key(id) - else - call fatal_error("Could not find cell " // trim(to_str(id)) & - &// " specified on tally filter.") - end if - end subroutine initialize_distribcell - - function text_label_distribcell(this, bin) result(label) - class(DistribcellFilter), intent(in) :: this - integer, intent(in) :: bin - character(MAX_LINE_LEN) :: label - - integer :: offset - type(Universe), pointer :: univ - - univ => universes(root_universe) - offset = 0 - label = '' - call find_offset(this % cell, univ, bin-1, offset, label) - label = "Distributed Cell " // label - end function text_label_distribcell - !=============================================================================== ! CellbornFilter methods !=============================================================================== @@ -788,288 +699,6 @@ contains end select end function text_label_energyfunction -!=============================================================================== -! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and -! a target offset to build a string which is the path from the base universe to -! the target cell with the given offset -!=============================================================================== - - recursive subroutine find_offset(i_cell, univ, target_offset, offset, path) - - integer, intent(in) :: i_cell ! The target cell index - type(Universe), intent(in) :: univ ! Universe to begin search - integer, intent(in) :: target_offset ! 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 - 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 - - ! 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_key(univ % id) - if (i_univ == root_universe) then - path = trim(path) // "u" // to_str(univ%id) - else - 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) == i_cell .and. offset == target_offset) then - c => cells(univ % cells(i)) - path = trim(path) // "->c" // to_str(c % id) - return - end if - end do - - ! Find the fill cell or lattice cell that we need to enter - do i = 1, n - - later_cell = .false. - - c => cells(univ % cells(i)) - - this_cell = .false. - - ! 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) - - ! Get the offset of the first lattice location - else - lat => lattices(c % fill) % obj - temp_offset = lat % offset(map, 1, 1, 1) - 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) + offset - - next_univ => universes(c % fill) - 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 - lat => lattices(c % fill) % 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 k = 1, n_x - do l = 1, n_y - do m = 1, n_z - - 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)) - 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, old_k, old_l, old_m) - offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) - 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, 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)) - 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, old_k, old_l, old_m) - offset = offset + lat_offset - next_univ => universes(lat % universes(old_k, old_l, old_m)) - 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 - - end if - end if - end do - end subroutine find_offset - !=============================================================================== ! C API FUNCTIONS !=============================================================================== diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 new file mode 100644 index 0000000000..d0f8adf07b --- /dev/null +++ b/src/tallies/tally_filter_distribcell.F90 @@ -0,0 +1,391 @@ +module tally_filter_distribcell + + use, intrinsic :: ISO_C_BINDING + + use hdf5, only: HID_T + + use constants + use error + use geometry_header + use hdf5_interface + use particle_header, only: Particle + use string, only: to_str + use tally_filter_header + + implicit none + private + public :: find_offset + +!=============================================================================== +! DISTRIBCELLFILTER specifies which distributed geometric cells tally events +! reside in. +!=============================================================================== + + type, public, extends(TallyFilter) :: DistribcellFilter + integer :: cell + contains + procedure :: get_all_bins => get_all_bins_distribcell + procedure :: to_statepoint => to_statepoint_distribcell + procedure :: text_label => text_label_distribcell + procedure :: initialize => initialize_distribcell + end type DistribcellFilter + +contains + + subroutine get_all_bins_distribcell(this, p, estimator, match) + class(DistribcellFilter), intent(in) :: this + type(Particle), intent(in) :: p + integer, intent(in) :: estimator + type(TallyFilterMatch), intent(inout) :: match + + integer :: distribcell_index, offset, i + + distribcell_index = cells(this % cell) % distribcell_index + offset = 0 + do i = 1, p % n_coord + if (cells(p % coord(i) % cell) % type == FILL_UNIVERSE) then + offset = offset + cells(p % coord(i) % cell) % & + offset(distribcell_index) + elseif (cells(p % coord(i) % cell) % type == FILL_LATTICE) then + if (lattices(p % coord(i + 1) % lattice) % obj & + % are_valid_indices([& + p % coord(i + 1) % lattice_x, & + p % coord(i + 1) % lattice_y, & + p % coord(i + 1) % lattice_z])) then + offset = offset + lattices(p % coord(i + 1) % lattice) % obj % & + offset(distribcell_index, & + p % coord(i + 1) % lattice_x, & + p % coord(i + 1) % lattice_y, & + p % coord(i + 1) % lattice_z) + end if + end if + if (this % cell == p % coord(i) % cell) then + call match % bins % push_back(offset + 1) + call match % weights % push_back(ONE) + return + end if + end do + end subroutine get_all_bins_distribcell + + subroutine to_statepoint_distribcell(this, filter_group) + class(DistribcellFilter), intent(in) :: this + integer(HID_T), intent(in) :: filter_group + + call write_dataset(filter_group, "type", "distribcell") + call write_dataset(filter_group, "n_bins", this % n_bins) + call write_dataset(filter_group, "bins", cells(this % cell) % id) + end subroutine to_statepoint_distribcell + + subroutine initialize_distribcell(this) + class(DistribcellFilter), intent(inout) :: this + + integer :: id + + ! Convert id to index. + id = this % cell + if (cell_dict % has_key(id)) then + this % cell = cell_dict % get_key(id) + else + call fatal_error("Could not find cell " // trim(to_str(id)) & + &// " specified on tally filter.") + end if + end subroutine initialize_distribcell + + function text_label_distribcell(this, bin) result(label) + class(DistribcellFilter), intent(in) :: this + integer, intent(in) :: bin + character(MAX_LINE_LEN) :: label + + integer :: offset + type(Universe), pointer :: univ + + univ => universes(root_universe) + offset = 0 + label = '' + call find_offset(this % cell, univ, bin-1, offset, label) + label = "Distributed Cell " // label + end function text_label_distribcell + +!=============================================================================== +! FIND_OFFSET (for distribcell) uses a given map number, a target cell ID, and +! a target offset to build a string which is the path from the base universe to +! the target cell with the given offset +!=============================================================================== + + recursive subroutine find_offset(i_cell, univ, target_offset, offset, path) + + integer, intent(in) :: i_cell ! The target cell index + type(Universe), intent(in) :: univ ! Universe to begin search + integer, intent(in) :: target_offset ! 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 + 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 + + ! 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_key(univ % id) + if (i_univ == root_universe) then + path = trim(path) // "u" // to_str(univ%id) + else + 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) == i_cell .and. offset == target_offset) then + c => cells(univ % cells(i)) + path = trim(path) // "->c" // to_str(c % id) + return + end if + end do + + ! Find the fill cell or lattice cell that we need to enter + do i = 1, n + + later_cell = .false. + + c => cells(univ % cells(i)) + + this_cell = .false. + + ! 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) + + ! Get the offset of the first lattice location + else + lat => lattices(c % fill) % obj + temp_offset = lat % offset(map, 1, 1, 1) + 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) + offset + + next_univ => universes(c % fill) + 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 + lat => lattices(c % fill) % 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 k = 1, n_x + do l = 1, n_y + do m = 1, n_z + + 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)) + 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, old_k, old_l, old_m) + offset = offset + lat_offset + next_univ => universes(lat % universes(old_k, old_l, old_m)) + 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, 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)) + 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, old_k, old_l, old_m) + offset = offset + lat_offset + next_univ => universes(lat % universes(old_k, old_l, old_m)) + 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 + + end if + end if + end do + end subroutine find_offset + +end module tally_filter_distribcell