diff --git a/src/geometry.F90 b/src/geometry.F90 index be2e2be84..ed5906032 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -17,14 +17,21 @@ module geometry interface function cell_contains_c(cell_ptr, xyz, uvw, on_surface) & bind(C, name="cell_contains") result(in_cell) - use ISO_C_BINDING - implicit none + import C_PTR, C_DOUBLE, C_INT32_T, C_BOOL type(C_PTR), intent(in), value :: cell_ptr real(C_DOUBLE), intent(in) :: xyz(3) real(C_DOUBLE), intent(in) :: uvw(3) integer(C_INT32_T), intent(in), value :: on_surface logical(C_BOOL) :: in_cell end function cell_contains_c + + function count_universe_instances(search_univ, target_univ_id) bind(C) & + result(count) + import C_INT32_T, C_INT + integer(C_INT32_T), intent(in), value :: search_univ + integer(C_INT32_T), intent(in), value :: target_univ_id + integer(C_INT) :: count + end function end interface contains @@ -532,27 +539,22 @@ contains ! routine is called once upon initialization. !=============================================================================== - subroutine calc_offsets(univ_id, map, univ, counts, found) + subroutine calc_offsets(univ_id, map, univ) 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 - logical, intent(inout) :: found(:,:) ! target found integer :: i ! index over cells integer :: j, k, m ! indices in lattice - integer :: n ! number of cells to search integer :: offset ! total offset for a given cell integer :: cell_index ! index in cells array type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next universe to cycle through class(Lattice), pointer :: lat ! pointer to current lattice - n = size(univ % cells) offset = 0 - do i = 1, n + do i = 1, size(univ % cells) cell_index = univ % cells(i) @@ -566,16 +568,8 @@ contains ! ==================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL elseif (c % type() == FILL_UNIVERSE) then - ! Set offset for the cell on this level c % offset(map) = offset - - ! Count contents of this cell - next_univ => universes(c % fill) - offset = offset + count_target(next_univ, counts, found, univ_id, map) - - ! Move into the next universe - next_univ => universes(c % fill) - c => cells(cell_index) + offset = offset + count_universe_instances(c % fill - 1, univ_id) ! ==================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL @@ -587,30 +581,24 @@ contains 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) lat % offset(map, j, k, m) = offset - next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) - offset = offset + & - count_target(next_univ, counts, found, univ_id, map) + offset = offset + count_universe_instances(& + lat % get([j-1, k-1, m-1]), univ_id) end do end do end do 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 if (lat % are_valid_indices([j, k, m])) then lat % offset(map, j, k, m) = offset - next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) - offset = offset + & - count_target(next_univ, counts, found, univ_id, map) + offset = offset + count_universe_instances(& + lat % get([j-1, k-1, m-1]), univ_id) end if end do end do @@ -622,134 +610,6 @@ contains end subroutine calc_offsets -!=============================================================================== -! COUNT_TARGET recursively totals the numbers of occurances of a given -! universe ID beginning with the universe given. -!=============================================================================== - - 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) :: univ_id ! target universe ID - integer, intent(in) :: map ! current map - - 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 - integer :: count ! number of times target located - type(Cell), pointer :: c ! pointer to current cell - type(Universe), pointer :: next_univ ! next univ to loop through - class(Lattice), pointer :: lat ! pointer to current lattice - - ! Don't research places already checked - if (found(universe_dict % get(univ % id), map)) then - count = counts(universe_dict % get(univ % id), map) - return - end if - - ! If this is the target, it can't contain itself. - ! Count = 1, then quit - if (univ % id == univ_id) then - count = 1 - counts(universe_dict % get(univ % id), map) = 1 - found(universe_dict % get(univ % id), map) = .true. - return - end if - - count = 0 - n = size(univ % cells) - - do i = 1, n - - cell_index = univ % cells(i) - - ! get pointer to cell - c => cells(cell_index) - - ! ==================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - if (c % type() == FILL_MATERIAL) then - - ! ==================================================================== - ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL - elseif (c % type() == FILL_UNIVERSE) then - - next_univ => universes(c % fill) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - return - end if - - count = count + count_target(next_univ, counts, found, univ_id, map) - c => cells(cell_index) - - ! ==================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - elseif (c % type() == FILL_LATTICE) then - - ! 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 % get([j-1, k-1, m-1]) + 1) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - cycle - end if - - count = count + & - count_target(next_univ, counts, found, univ_id, map) - - end do - end do - end do - - 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 - if (lat % are_valid_indices([j, k, m])) then - next_univ => universes(lat % get([j-1, k-1, m-1]) + 1) - - ! Found target - stop since target cannot contain itself - if (next_univ % id == univ_id) then - count = count + 1 - cycle - end if - - count = count + & - count_target(next_univ, counts, found, univ_id, map) - end if - end do - end do - end do - - end select - - end if - end do - - counts(universe_dict % get(univ % id), map) = count - found(universe_dict % get(univ % id), map) = .true. - - end function count_target - !=============================================================================== ! MAXIMUM_LEVELS determines the maximum number of nested coordinate levels in ! the geometry diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index ec4da6f12..906845cb9 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -1,3 +1,6 @@ +//! \file geometry_aux.cpp +//! Auxilary functions for geometry initialization and general data handling. + #include #include @@ -11,6 +14,10 @@ namespace openmc { +//============================================================================== +//! Replace Universe, Lattice, and Material IDs with indices. +//============================================================================== + extern "C" void adjust_indices_c() { @@ -58,6 +65,12 @@ adjust_indices_c() } } +//============================================================================== +//! Figure out which Universe is the root universe. +//! +//! This function looks for a universe that is not listed in a Cell::fill or in +//! a Lattice. +//! @return The index of the root universe. //============================================================================== extern "C" int32_t @@ -101,10 +114,17 @@ find_root_universe() return root_univ; } +//============================================================================== +//! Recursively search through the geometry and count cell instances. +//! +//! This function will update the Cell::n_instances value for each cell in the +//! geometry. +//! @param univ_indx The index of the universe to begin searching from (probably +//! the root universe). //============================================================================== extern "C" void -count_instances(int32_t univ_indx) +count_cell_instances(int32_t univ_indx) { for (int32_t cell_indx : universes_c[univ_indx]->cells) { Cell &c {*cells_c[cell_indx]}; @@ -112,16 +132,51 @@ count_instances(int32_t univ_indx) if (c.type == FILL_UNIVERSE) { // This cell contains another universe. Recurse into that universe. - count_instances(c.fill-1); // TODO: off-by-one + count_cell_instances(c.fill-1); // TODO: off-by-one } else if (c.type == FILL_LATTICE) { // This cell contains a lattice. Recurse into the lattice universes. Lattice &lat {*lattices_c[c.fill-1]}; // TODO: off-by-one for (auto it = lat.begin(); it != lat.end(); ++it) { - count_instances(*it); + count_cell_instances(*it); } } } } +//============================================================================== +//! Recursively search through universes and count the number of instances of +//! the target universe in the geometry tree. +//! @param search_univ The index of the universe to begin searching from. +//! @param target_univ_id The ID of the universe to be counted. +//============================================================================== + +extern "C" int +count_universe_instances(int32_t search_univ, int32_t target_univ_id) +{ + // If this is the target, it can't contain itself. + if (universes_c[search_univ]->id == target_univ_id) { + return 1; + } + + int count {0}; + for (int32_t cell_indx : universes_c[search_univ]->cells) { + Cell &c {*cells_c[cell_indx]}; + + if (c.type == FILL_UNIVERSE) { + int32_t next_univ = c.fill - 1; // TODO: off-by-one + count += count_universe_instances(next_univ, target_univ_id); + + } else if (c.type == FILL_LATTICE) { + Lattice &lat {*lattices_c[c.fill - 1]}; //TODO: off-by-one + for (auto it = lat.begin(); it != lat.end(); ++it) { + int32_t next_univ = *it; + count += count_universe_instances(next_univ, target_univ_id); + } + } + } + + return count; +} + } // namespace openmc diff --git a/src/input_xml.F90 b/src/input_xml.F90 index e43906c9a..44a161f47 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -50,7 +50,7 @@ module input_xml subroutine adjust_indices_c() bind(C) end subroutine adjust_indices_c - subroutine count_instances_c(univ_indx) bind(C, name='count_instances') + subroutine count_cell_instances(univ_indx) bind(C) import C_INT32_T integer(C_INT32_T), intent(in), value :: univ_indx end subroutine @@ -143,7 +143,7 @@ contains ! Perform some final operations to set up the geometry call adjust_indices() - call count_instances_c(root_universe-1) + call count_cell_instances(root_universe-1) ! After reading input and basic geometry setup is complete, build lists of ! neighboring cells for efficient tracking @@ -3829,8 +3829,6 @@ contains 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. @@ -3885,12 +3883,12 @@ contains end do ! Allocate offset maps at each level in the geometry - call allocate_offsets(univ_list, counts, found) + call allocate_offsets(univ_list) ! 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) + call calc_offsets(univ_list(i), i, universes(j)) end do end do @@ -3901,11 +3899,9 @@ contains ! memory for distribcell offset tables !=============================================================================== - recursive subroutine allocate_offsets(univ_list, counts, found) + recursive subroutine allocate_offsets(univ_list) 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 @@ -3943,15 +3939,6 @@ contains ! 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