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

Substantially changed logic for fission bank algorithm to improve performance and handle multiple messages.

Browse source
This commit is contained in:
Paul Romano 2012-02-07 14:37:16 -05:00
parent 97b838fab4
commit c3b91a3cbf
8 changed files with 248 additions and 160 deletions
Download patch file Download diff file Expand all files Collapse all files

2
src/DEPENDENCIES
View file

@ -118,6 +118,8 @@ intercycle.o: global.o
intercycle.o: output.o
intercycle.o: particle_header.o
intercycle.o: random_lcg.o
intercycle.o: search.o
intercycle.o: string.o
intercycle.o: tally_header.o
intercycle.o: timing.o

2
src/global.F90
View file

@ -117,6 +117,7 @@ module global
integer(8) :: bank_first ! index of first particle in bank
integer(8) :: bank_last ! index of last particle in bank
integer(8) :: work ! number of particles per processor
integer(8) :: maxwork ! maximum number of particles per processor
! cycle keff
real(8) :: keff = ONE
@ -150,7 +151,6 @@ module global
type(Timer) :: time_ic_tallies ! timer for intercycle accumulate tallies
type(Timer) :: time_ic_sample ! timer for intercycle sampling
type(Timer) :: time_ic_sendrecv ! timer for intercycle SEND/RECV
type(Timer) :: time_ic_rebuild ! timer for intercycle source bank rebuild
type(Timer) :: time_inactive ! timer for inactive cycles
type(Timer) :: time_active ! timer for active cycles
type(Timer) :: time_compute ! timer for computation

3
src/hdf5_interface.F90
View file

@ -620,7 +620,6 @@ contains
call hdf5_make_double(timing_group, "time_tallies", time_ic_tallies % elapsed)
call hdf5_make_double(timing_group, "time_sample", time_ic_sample % elapsed)
call hdf5_make_double(timing_group, "time_sendrecv", time_ic_sendrecv % elapsed)
call hdf5_make_double(timing_group, "time_rebuild", time_ic_rebuild % elapsed)
call hdf5_make_double(timing_group, "time_inactive", time_inactive % elapsed)
call hdf5_make_double(timing_group, "time_active", time_active % elapsed)
call hdf5_make_double(timing_group, "time_total", time_total % elapsed)
@ -642,8 +641,6 @@ contains
"description", "Time between cycles sampling source sites (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_sendrecv", &
"description", "Time between cycles SEND/RECVing source sites (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_rebuild", &
"description", "Time between cycles reconstructing source bank (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_inactive", &
"description", "Total time in inactive cycles (s)", hdf5_err)
call h5ltset_attribute_string_f(timing_group, "time_active", &

277
src/intercycle.F90
View file

@ -7,6 +7,8 @@ module intercycle
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use string, only: to_str
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
@ -24,30 +26,27 @@ contains
subroutine synchronize_bank()
integer :: i, j, k ! loop indices
integer(8) :: start ! starting index in local fission bank
integer(8) :: finish ! ending index in local fission bank
integer(8) :: total ! total sites in global fission bank
integer(8) :: index_temp ! index in temporary source bank
integer :: send_to_left ! total # of bank sites to send/recv to or from left
integer :: send_to_right ! total # of bank sites to send/recv to or from right
integer(8) :: sites_needed ! # of sites to be sampled
real(8) :: p_sample ! probability of sampling a site
type(Bank), allocatable :: &
temp_sites(:), & ! local array of extra sites on each node
left_bank(:), & ! bank sites to send/recv to or from left node
right_bank(:) ! bank sites to send/recv to or fram right node
integer :: i ! loop indices
integer :: j ! loop indices
integer(8) :: start ! starting index in global bank
integer(8) :: finish ! ending index in global bank
integer(8) :: total ! total sites in global fission bank
integer(8) :: index_temp ! index in temporary source bank
integer(8) :: sites_needed ! # of sites to be sampled
real(8) :: p_sample ! probability of sampling a site
type(Bank), save, allocatable :: &
& temp_sites(:) ! local array of extra sites on each node
#ifdef MPI
integer :: status(MPI_STATUS_SIZE) ! message status
integer :: request ! communication request for sending sites
integer :: request_left ! communication request for recv sites from left
integer :: request_right ! communication request for recv sites from right
integer :: n ! number of sites to send/recv
integer :: neighbor ! processor to send/recv data from
integer :: request(20) ! communication request for send/recving sites
integer :: n_request ! number of communication requests
integer(8) :: index_local ! index in local source bank
integer(8), save, allocatable :: &
& bank_position(:) ! starting positions in global source bank
#endif
message = "Collecting number of fission sites..."
call write_message(8)
! In order to properly understand the fission bank algorithm, you need to
! think of the fission and source bank as being one global array divided
! over multiple processors. At the start, each processor has a random amount
@ -73,11 +72,12 @@ contains
#endif
! If there are not that many particles per cycle, it's possible that no
! fission sites were created at all. There is no way to proceed in that
! circumstance so issue a fatal error and stop the simulation.
! fission sites were created at all on a single processor. Rather than add
! extra logic to treat this circumstance, we really want to ensure the user
! runs enough particles to avoid this in the first place.
if (total == 0) then
message = "No fission sites banked!"
if (n_bank == 0) then
message = "No fission sites banked on procesor " // to_str(rank)
call fatal_error()
end if
@ -99,18 +99,14 @@ contains
end if
p_sample = real(sites_needed,8)/real(total,8)
message = "Sampling fission sites..."
call write_message(8)
call timer_start(time_ic_sample)
! ==========================================================================
! SAMPLE N_PARTICLES FROM FISSION BANK AND PLACE IN TEMP_SITES
! Allocate temporary source bank
! TODO: consider only allocating this once?
index_temp = 0_8
allocate(temp_sites(2*work))
if (.not. allocated(temp_sites)) allocate(temp_sites(3*work))
do i = 1, int(n_bank,4)
@ -132,47 +128,42 @@ contains
end if
end do
! At this point, the sampling of source sites is done and now we need to figure out where to send source sites
! Now that we've sampled sites, check where the boundaries of data are for
! the source bank
! At this point, the sampling of source sites is done and now we need to
! figure out where to send source sites. Since it is possible that one
! processor's share of the source bank spans more than just the immediate
! neighboring processors, we have to perform an ALLGATHER to determine the
! indices for all processors
#ifdef MPI
! First do an exclusive scan to get the starting indices for
start = 0_8
call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, &
MPI_COMM_WORLD, mpi_err)
finish = start + index_temp
total = finish
call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, &
MPI_COMM_WORLD, mpi_err)
! Allocate space for bank_position if this hasn't been done yet
if (.not. allocated(bank_position)) allocate(bank_position(n_procs))
call MPI_ALLGATHER(start, 1, MPI_INTEGER8, bank_position, 1, &
MPI_INTEGER8, MPI_COMM_WORLD, mpi_err)
#else
start = 0_8
finish = index_temp
total = index_temp
#endif
! Determine how many sites to send to adjacent nodes
send_to_left = int(bank_first - 1_8 - start, 4)
send_to_right = int(finish - bank_last, 4)
! Check to make sure number of sites is not more than size of bank
if (abs(send_to_left) > work .or. abs(send_to_right) > work) then
message = "Tried sending sites to neighboring process greater than " &
// "the size of the source bank."
call fatal_error()
end if
! Now that the sampling is complete, we need to ensure that we have exactly
! n_particles source sites. The way this is done in a reproducible manner is
! to adjust only the source sites on the last processor.
if (rank == n_procs - 1) then
if (total > n_particles) then
if (finish > n_particles) then
! If we have extra sites sampled, we will simply discard the extra
! ones on the last processor
if (rank == n_procs - 1) then
index_temp = index_temp - send_to_right
end if
index_temp = n_particles - start
elseif (total < n_particles) then
! If we have too few sites, grab sites from the very end of the
elseif (finish < n_particles) then
! If we have too few sites, repeat sites from the very end of the
! fission bank
sites_needed = n_particles - total
sites_needed = n_particles - finish
do i = 1, int(sites_needed,4)
index_temp = index_temp + 1
temp_sites(index_temp) = fission_bank(n_bank - sites_needed + i)
@ -180,117 +171,107 @@ contains
end if
! the last processor should not be sending sites to right
send_to_right = 0
finish = bank_last
end if
call timer_stop(time_ic_sample)
call timer_start(time_ic_sendrecv)
#ifdef MPI
message = "Sending fission sites..."
call write_message(8)
! ==========================================================================
! SEND BANK SITES TO NEIGHBORS
allocate(left_bank(abs(send_to_left)))
allocate(right_bank(abs(send_to_right)))
if (send_to_right > 0) then
i = index_temp - send_to_right + 1
call MPI_ISEND(temp_sites(i), send_to_right, MPI_BANK, rank+1, 0, &
MPI_COMM_WORLD, request, mpi_err)
else if (send_to_right < 0) then
call MPI_IRECV(right_bank, -send_to_right, MPI_BANK, rank+1, 1, &
MPI_COMM_WORLD, request_right, mpi_err)
index_local = 1
n_request = 0
SEND_SITES: do while (start < finish)
! Determine the index of the processor which has the first part of the
! source_bank for the local processor
neighbor = start / maxwork
! Determine the number of sites to send
n = min((neighbor + 1)*maxwork, finish) - start
! Initiate an asynchronous send of source sites to the neighboring
! process
if (neighbor /= rank) then
n_request = n_request + 1
call MPI_ISEND(temp_sites(index_local), n, MPI_BANK, neighbor, &
rank, MPI_COMM_WORLD, request(n_request), mpi_err)
end if
! Increment all indices
start = start + n
index_local = index_local + n
neighbor = neighbor + 1
end do SEND_SITES
! ==========================================================================
! RECEIVE BANK SITES FROM NEIGHBORS OR TEMPORARY BANK
start = bank_first - 1
index_local = 1
! Determine what process has the source sites that will need to be stored at
! the beginning of this processor's source bank.
if (start >= bank_position(n_procs)) then
neighbor = n_procs - 1
else
neighbor = binary_search(bank_position, n_procs, start) - 1
end if
if (send_to_left < 0) then
call MPI_IRECV(left_bank, -send_to_left, MPI_BANK, rank-1, 0, &
MPI_COMM_WORLD, request_left, mpi_err)
else if (send_to_left > 0) then
call MPI_ISEND(temp_sites(1), send_to_left, MPI_BANK, rank-1, 1, &
MPI_COMM_WORLD, request, mpi_err)
end if
RECV_SITES: do while (start < bank_last)
! Determine how many sites need to be received
if (neighbor == n_procs - 1) then
n = min(n_particles, (rank+1)*maxwork) - start
else
n = min(bank_position(neighbor+2), min(n_particles, &
(rank+1)*maxwork)) - start
end if
if (neighbor /= rank) then
! If the source sites are not on this processor, initiate an
! asynchronous receive for the source sites
n_request = n_request + 1
call MPI_IRECV(source_bank(index_local), n, MPI_BANK, &
neighbor, neighbor, MPI_COMM_WORLD, request(n_request), mpi_err)
else
! If the source sites are on this procesor, we can simply copy them
! from the temp_sites bank
index_temp = start - bank_position(rank+1) + 1
source_bank(index_local:index_local+n-1) = &
temp_sites(index_temp:index_temp+n-1)
end if
! Increment all indices
start = start + n
index_local = index_local + n
neighbor = neighbor + 1
end do RECV_SITES
! Send we initiated a series of asynchronous ISENDs and IRECVs, now we have
! to ensure that the data has actually been communicated before moving on to
! the next cycle
call MPI_WAITALL(n_request, request, MPI_STATUSES_IGNORE, mpi_err)
! Deallocate space for bank_position on the last cycle
if (current_cycle == n_cycles) deallocate(bank_position)
#else
source_bank = temp_sites(1:n_particles)
#endif
call timer_stop(time_ic_sendrecv)
call timer_start(time_ic_rebuild)
message = "Constructing source bank..."
call write_message(8)
! ==========================================================================
! RECONSTRUCT SOURCE BANK
if (send_to_left < 0 .and. send_to_right >= 0) then
i = -send_to_left ! size of first block
j = int(index_temp,4) - send_to_right ! size of second block
source_bank(i+1:i+j) = temp_sites(1:j)
#ifdef MPI
call MPI_WAIT(request_left, status, mpi_err)
#endif
call copy_from_bank(left_bank, 1, i)
else if (send_to_left >= 0 .and. send_to_right < 0) then
i = int(index_temp,4) - send_to_left ! size of first block
j = -send_to_right ! size of second block
source_bank(1:i) = temp_sites(send_to_left+1:send_to_left+i)
#ifdef MPI
call MPI_WAIT(request_right, status, mpi_err)
#endif
call copy_from_bank(right_bank, i+1, j)
else if (send_to_left >= 0 .and. send_to_right >= 0) then
i = int(index_temp,4) - send_to_left - send_to_right
source_bank = temp_sites(send_to_left+1:send_to_left+i)
else if (send_to_left < 0 .and. send_to_right < 0) then
i = -send_to_left
j = int(index_temp,4)
k = -send_to_right
source_bank(i+1:i+j) = temp_sites(1:j)
#ifdef MPI
call MPI_WAIT(request_left, status, mpi_err)
#endif
call copy_from_bank(left_bank, 1, i)
#ifdef MPI
call MPI_WAIT(request_right, status, mpi_err)
#endif
call copy_from_bank(right_bank, i+j+1, k)
end if
call timer_stop(time_ic_rebuild)
#ifdef MPI
deallocate(left_bank)
deallocate(right_bank)
#endif
deallocate(temp_sites)
! Deallocate space for the temporary source bank on the last cycle
if (current_cycle == n_cycles) deallocate(temp_sites)
end subroutine synchronize_bank
!===============================================================================
! COPY_FROM_BANK
!===============================================================================
subroutine copy_from_bank(temp_bank, i_start, n_sites)
integer, intent(in) :: n_sites ! # of bank sites to copy
type(Bank), intent(in) :: temp_bank(n_sites)
integer, intent(in) :: i_start ! starting index in source_bank
integer :: i ! index in temp_bank
integer :: i_source ! index in source_bank
type(Bank), pointer :: src => null()
do i = 1, n_sites
i_source = i_start + i - 1
src => source_bank(i_source)
! Copy attributes
src % xyz = temp_bank(i) % xyz
src % uvw = temp_bank(i) % uvw
src % E = temp_bank(i) % E
end do
end subroutine copy_from_bank
!===============================================================================
! REDUCE_TALLIES collects all the results from tallies onto one processor
!===============================================================================

1
src/output.F90
View file

@ -937,7 +937,6 @@ contains
write(ou,100) " Accumulating tallies", time_ic_tallies % elapsed
write(ou,100) " Sampling source sites", time_ic_sample % elapsed
write(ou,100) " SEND/RECV source sites", time_ic_sendrecv % elapsed
write(ou,100) " Reconstruct source bank", time_ic_rebuild % elapsed
write(ou,100) "Total time in inactive cycles", time_inactive % elapsed
write(ou,100) "Total time in active cycles", time_active % elapsed
write(ou,100) "Total time elapsed", time_total % elapsed

8
src/physics.F90
View file

@ -1331,7 +1331,7 @@ contains
! determine index on incoming energy grid and interpolation factor
lc = 2 + 2*NR
i = binary_search(edist % data(lc+1), NE, E_in)
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
(edist%data(lc+i+1) - edist%data(lc+i))
@ -1407,7 +1407,7 @@ contains
i = NE - 1
r = ONE
else
i = binary_search(edist % data(lc+1), NE, E_in)
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
(edist%data(lc+i+1) - edist%data(lc+i))
end if
@ -1633,7 +1633,7 @@ contains
i = NE - 1
r = ONE
else
i = binary_search(edist % data(lc+1), NE, E_in)
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
(edist%data(lc+i+1) - edist%data(lc+i))
end if
@ -1776,7 +1776,7 @@ contains
i = NE - 1
r = ONE
else
i = binary_search(edist % data(lc+1), NE, E_in)
i = binary_search(edist % data(lc+1:lc+NE), NE, E_in)
r = (E_in - edist%data(lc+i)) / &
(edist%data(lc+i+1) - edist%data(lc+i))
end if

114
src/search.F90
View file

@ -6,6 +6,10 @@ module search
integer, parameter :: MAX_ITERATION = 64
interface binary_search
module procedure binary_search_real, binary_search_int4, binary_search_int8
end interface binary_search
contains
!===============================================================================
@ -13,7 +17,7 @@ contains
! value lies in the array. This is used extensively for energy grid searching
!===============================================================================
function binary_search(array, n, val) result(array_index)
function binary_search_real(array, n, val) result(array_index)
integer, intent(in) :: n
real(8), intent(in) :: array(n)
@ -64,6 +68,112 @@ contains
array_index = L
end function binary_search
end function binary_search_real
function binary_search_int4(array, n, val) result(array_index)
integer, intent(in) :: n
integer, intent(in) :: array(n)
integer, intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
real(8) :: testval
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
message = "Value outside of array during binary search"
call fatal_error()
end if
n_iteration = 0
do while (R - L > 1)
! Check boundaries
if (val > array(L) .and. val < array(L+1)) then
array_index = L
return
elseif (val > array(R-1) .and. val < array(R)) then
array_index = R - 1
return
end if
! Find values at midpoint
array_index = L + (R - L)/2
testval = array(array_index)
if (val >= testval) then
L = array_index
elseif (val < testval) then
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
message = "Reached maximum number of iterations on binary search."
call fatal_error()
end if
end do
array_index = L
end function binary_search_int4
function binary_search_int8(array, n, val) result(array_index)
integer, intent(in) :: n
integer(8), intent(in) :: array(n)
integer(8), intent(in) :: val
integer :: array_index
integer :: L
integer :: R
integer :: n_iteration
real(8) :: testval
L = 1
R = n
if (val < array(L) .or. val > array(R)) then
message = "Value outside of array during binary search"
call fatal_error()
end if
n_iteration = 0
do while (R - L > 1)
! Check boundaries
if (val > array(L) .and. val < array(L+1)) then
array_index = L
return
elseif (val > array(R-1) .and. val < array(R)) then
array_index = R - 1
return
end if
! Find values at midpoint
array_index = L + (R - L)/2
testval = array(array_index)
if (val >= testval) then
L = array_index
elseif (val < testval) then
R = array_index
end if
! check for large number of iterations
n_iteration = n_iteration + 1
if (n_iteration == MAX_ITERATION) then
message = "Reached maximum number of iterations on binary search."
call fatal_error()
end if
end do
array_index = L
end function binary_search_int8
end module search

1
src/source.F90
View file

@ -24,7 +24,6 @@ contains
integer :: i ! loop index over processors
integer(8) :: j ! loop index over bank sites
integer :: k ! dummy loop index
integer(8) :: maxwork ! maxinum # of particles per processor
integer(8) :: bytes ! size of fission/source bank
integer(8) :: id ! particle id
integer :: alloc_err ! allocation error code