diff --git a/src/initialize.F90 b/src/initialize.F90 index c2d5e31860..7752b86437 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -12,8 +12,6 @@ module initialize use dict_header, only: DictIntInt, ElemKeyValueII use set_header, only: SetInt use error, only: fatal_error, warning - use geometry, only: calc_offsets, & - maximum_levels use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& root_universe use global @@ -89,18 +87,6 @@ contains ! Read XML input files call read_input_xml() - ! Initialize distribcell_filters - call prepare_distribcell() - - ! Check to make sure there are not too many nested coordinate levels in the - ! geometry since the coordinate list is statically allocated for performance - ! reasons - if (maximum_levels(universes(root_universe)) > MAX_COORD) then - call fatal_error("Too many nested coordinate levels in the geometry. & - &Try increasing the maximum number of coordinate levels by & - &providing the CMake -Dmaxcoord= option.") - end if - if (run_mode /= MODE_PLOTTING) then ! Set up tally procedure pointers call init_tally_routines() @@ -490,179 +476,4 @@ contains end subroutine allocate_banks -!=============================================================================== -! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the -! offsets for distribcells -!=============================================================================== - - subroutine prepare_distribcell() - - integer :: i, j ! Tally, filter loop counters - logical :: distribcell_active ! Does simulation use distribcell? - integer, allocatable :: univ_list(:) ! Target offsets - integer, allocatable :: counts(:,:) ! Target count - logical, allocatable :: found(:,:) ! Target found - - ! Assume distribcell is not needed until proven otherwise. - distribcell_active = .false. - - ! We need distribcell if any tallies have distribcell filters. - do i = 1, n_tallies - do j = 1, size(tallies(i) % obj % filter) - select type(filt => filters(tallies(i) % obj % filter(j)) % obj) - type is (DistribcellFilter) - distribcell_active = .true. - end select - end do - end do - - ! We also need distribcell if any distributed materials or distributed - ! temperatues are present. - if (.not. distribcell_active) then - do i = 1, n_cells - if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then - distribcell_active = .true. - exit - end if - end do - end if - - ! If distribcell isn't used in this simulation then no more work left to do. - if (.not. distribcell_active) return - - ! Make sure the number of materials and temperatures matches the number of - ! cell instances. - do i = 1, n_cells - associate (c => cells(i)) - if (size(c % material) > 1) then - if (size(c % material) /= c % instances) then - call fatal_error("Cell " // trim(to_str(c % id)) // " was & - &specified with " // trim(to_str(size(c % material))) & - // " materials but has " // trim(to_str(c % instances)) & - // " distributed instances. The number of materials must & - &equal one or the number of instances.") - end if - end if - if (size(c % sqrtkT) > 1) then - if (size(c % sqrtkT) /= c % instances) then - call fatal_error("Cell " // trim(to_str(c % id)) // " was & - &specified with " // trim(to_str(size(c % sqrtkT))) & - // " temperatures but has " // trim(to_str(c % instances)) & - // " distributed instances. The number of temperatures must & - &equal one or the number of instances.") - end if - end if - end associate - end do - - ! Allocate offset maps at each level in the geometry - call allocate_offsets(univ_list, counts, found) - - ! Calculate offsets for each target distribcell - do i = 1, n_maps - do j = 1, n_universes - call calc_offsets(univ_list(i), i, universes(j), counts, found) - end do - end do - - end subroutine prepare_distribcell - -!=============================================================================== -! ALLOCATE_OFFSETS determines the number of maps needed and allocates required -! memory for distribcell offset tables -!=============================================================================== - - recursive subroutine allocate_offsets(univ_list, counts, found) - - integer, intent(out), allocatable :: univ_list(:) ! Target offsets - integer, intent(out), allocatable :: counts(:,:) ! Target count - logical, intent(out), allocatable :: found(:,:) ! Target found - - integer :: i, j, k ! Loop counters - type(SetInt) :: cell_list ! distribells to track - - ! Begin gathering list of cells in distribcell tallies - n_maps = 0 - - ! List all cells referenced in distribcell filters. - do i = 1, n_tallies - do j = 1, size(tallies(i) % obj % filter) - select type(filt => filters(tallies(i) % obj % filter(j)) % obj) - type is (DistribcellFilter) - call cell_list % add(filt % cell) - end select - end do - end do - - ! List all cells with multiple (distributed) materials or temperatures. - do i = 1, n_cells - if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then - call cell_list % add(i) - end if - end do - - ! Compute the number of unique universes containing these distribcells - ! to determine the number of offset tables to allocate - do i = 1, n_universes - do j = 1, size(universes(i) % cells) - if (cell_list % contains(universes(i) % cells(j))) then - n_maps = n_maps + 1 - end if - end do - end do - - ! Allocate the list of offset tables for each unique universe - allocate(univ_list(n_maps)) - - ! Allocate list to accumulate target distribcell counts in each universe - allocate(counts(n_universes, n_maps)) - counts(:,:) = 0 - - ! Allocate list to track if target distribcells are found in each universe - allocate(found(n_universes, n_maps)) - found(:,:) = .false. - - - ! Search through universes for distributed cells and assign each one a - ! unique distribcell array index. - k = 1 - do i = 1, n_universes - do j = 1, size(universes(i) % cells) - if (cell_list % contains(universes(i) % cells(j))) then - cells(universes(i) % cells(j)) % distribcell_index = k - univ_list(k) = universes(i) % id - k = k + 1 - end if - end do - end do - - ! Allocate the offset tables for lattices - do i = 1, n_lattices - associate(lat => lattices(i) % obj) - select type(lat) - - type is (RectLattice) - allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & - lat % n_cells(3))) - type is (HexLattice) - allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & - 2 * lat % n_rings - 1, lat % n_axial)) - end select - - lat % offset(:, :, :, :) = 0 - end associate - end do - - ! Allocate offset table for fill cells - do i = 1, n_cells - if (cells(i) % type /= FILL_MATERIAL) then - allocate(cells(i) % offset(n_maps)) - end if - end do - - ! Free up memory - call cell_list % clear() - - end subroutine allocate_offsets - end module initialize diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e7ccc1f7cf..9f0b77b330 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -11,7 +11,8 @@ module input_xml use distribution_univariate use endf, only: reaction_name use error, only: fatal_error, warning - use geometry, only: count_instance, neighbor_lists + use geometry, only: calc_offsets, maximum_levels, count_instance, & + neighbor_lists use geometry_header, only: Cell, Lattice, RectLattice, HexLattice, & get_temperatures, root_universe use global @@ -68,6 +69,10 @@ contains end if call read_tallies_xml() + + ! Initialize distribcell_filters + call prepare_distribcell() + if (cmfd_run) call configure_cmfd() if (.not. run_CE) then @@ -105,6 +110,15 @@ contains ! Determine desired txemperatures for each nuclide and S(a,b) table call get_temperatures(nuc_temps, sab_temps) + ! Check to make sure there are not too many nested coordinate levels in the + ! geometry since the coordinate list is statically allocated for performance + ! reasons + if (maximum_levels(universes(root_universe)) > MAX_COORD) then + call fatal_error("Too many nested coordinate levels in the geometry. & + &Try increasing the maximum number of coordinate levels by & + &providing the CMake -Dmaxcoord= option.") + end if + end subroutine finalize_geometry !=============================================================================== @@ -5583,4 +5597,179 @@ contains end subroutine adjust_indices +!=============================================================================== +! PREPARE_DISTRIBCELL initializes any distribcell filters present and sets the +! offsets for distribcells +!=============================================================================== + + subroutine prepare_distribcell() + + integer :: i, j ! Tally, filter loop counters + logical :: distribcell_active ! Does simulation use distribcell? + integer, allocatable :: univ_list(:) ! Target offsets + integer, allocatable :: counts(:,:) ! Target count + logical, allocatable :: found(:,:) ! Target found + + ! Assume distribcell is not needed until proven otherwise. + distribcell_active = .false. + + ! We need distribcell if any tallies have distribcell filters. + do i = 1, n_tallies + do j = 1, size(tallies(i) % obj % filter) + select type(filt => filters(tallies(i) % obj % filter(j)) % obj) + type is (DistribcellFilter) + distribcell_active = .true. + end select + end do + end do + + ! We also need distribcell if any distributed materials or distributed + ! temperatues are present. + if (.not. distribcell_active) then + do i = 1, n_cells + if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then + distribcell_active = .true. + exit + end if + end do + end if + + ! If distribcell isn't used in this simulation then no more work left to do. + if (.not. distribcell_active) return + + ! Make sure the number of materials and temperatures matches the number of + ! cell instances. + do i = 1, n_cells + associate (c => cells(i)) + if (size(c % material) > 1) then + if (size(c % material) /= c % instances) then + call fatal_error("Cell " // trim(to_str(c % id)) // " was & + &specified with " // trim(to_str(size(c % material))) & + // " materials but has " // trim(to_str(c % instances)) & + // " distributed instances. The number of materials must & + &equal one or the number of instances.") + end if + end if + if (size(c % sqrtkT) > 1) then + if (size(c % sqrtkT) /= c % instances) then + call fatal_error("Cell " // trim(to_str(c % id)) // " was & + &specified with " // trim(to_str(size(c % sqrtkT))) & + // " temperatures but has " // trim(to_str(c % instances)) & + // " distributed instances. The number of temperatures must & + &equal one or the number of instances.") + end if + end if + end associate + end do + + ! Allocate offset maps at each level in the geometry + call allocate_offsets(univ_list, counts, found) + + ! Calculate offsets for each target distribcell + do i = 1, n_maps + do j = 1, n_universes + call calc_offsets(univ_list(i), i, universes(j), counts, found) + end do + end do + + end subroutine prepare_distribcell + +!=============================================================================== +! ALLOCATE_OFFSETS determines the number of maps needed and allocates required +! memory for distribcell offset tables +!=============================================================================== + + recursive subroutine allocate_offsets(univ_list, counts, found) + + integer, intent(out), allocatable :: univ_list(:) ! Target offsets + integer, intent(out), allocatable :: counts(:,:) ! Target count + logical, intent(out), allocatable :: found(:,:) ! Target found + + integer :: i, j, k ! Loop counters + type(SetInt) :: cell_list ! distribells to track + + ! Begin gathering list of cells in distribcell tallies + n_maps = 0 + + ! List all cells referenced in distribcell filters. + do i = 1, n_tallies + do j = 1, size(tallies(i) % obj % filter) + select type(filt => filters(tallies(i) % obj % filter(j)) % obj) + type is (DistribcellFilter) + call cell_list % add(filt % cell) + end select + end do + end do + + ! List all cells with multiple (distributed) materials or temperatures. + do i = 1, n_cells + if (size(cells(i) % material) > 1 .or. size(cells(i) % sqrtkT) > 1) then + call cell_list % add(i) + end if + end do + + ! Compute the number of unique universes containing these distribcells + ! to determine the number of offset tables to allocate + do i = 1, n_universes + do j = 1, size(universes(i) % cells) + if (cell_list % contains(universes(i) % cells(j))) then + n_maps = n_maps + 1 + end if + end do + end do + + ! Allocate the list of offset tables for each unique universe + allocate(univ_list(n_maps)) + + ! Allocate list to accumulate target distribcell counts in each universe + allocate(counts(n_universes, n_maps)) + counts(:,:) = 0 + + ! Allocate list to track if target distribcells are found in each universe + allocate(found(n_universes, n_maps)) + found(:,:) = .false. + + + ! Search through universes for distributed cells and assign each one a + ! unique distribcell array index. + k = 1 + do i = 1, n_universes + do j = 1, size(universes(i) % cells) + if (cell_list % contains(universes(i) % cells(j))) then + cells(universes(i) % cells(j)) % distribcell_index = k + univ_list(k) = universes(i) % id + k = k + 1 + end if + end do + end do + + ! Allocate the offset tables for lattices + do i = 1, n_lattices + associate(lat => lattices(i) % obj) + select type(lat) + + type is (RectLattice) + allocate(lat % offset(n_maps, lat % n_cells(1), lat % n_cells(2), & + lat % n_cells(3))) + type is (HexLattice) + allocate(lat % offset(n_maps, 2 * lat % n_rings - 1, & + 2 * lat % n_rings - 1, lat % n_axial)) + end select + + lat % offset(:, :, :, :) = 0 + end associate + end do + + ! Allocate offset table for fill cells + do i = 1, n_cells + if (cells(i) % type /= FILL_MATERIAL) then + allocate(cells(i) % offset(n_maps)) + end if + end do + + ! Free up memory + call cell_list % clear() + + end subroutine allocate_offsets + end module input_xml