crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-26 05:05:30 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

Move count_target (distribcell) to C++

Browse source
This commit is contained in:
Sterling Harper 2018-05-11 22:43:27 -04:00
parent 5419403e10
commit df6e2a2cef
3 changed files with 79 additions and 177 deletions
Download patch file Download diff file Expand all files Collapse all files

172
src/geometry.F90
View file

@ -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

61
src/geometry_aux.cpp
View file

@ -1,3 +1,6 @@
//! \file geometry_aux.cpp
//! Auxilary functions for geometry initialization and general data handling.
#include <sstream>
#include <unordered_set>
@ -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

23
src/input_xml.F90
View file

@ -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