diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index ceb2ea472..e5d327bb8 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 925a2cc01..3c6fd8c1b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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 diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 4529e0780..87cb5b43e 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -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", & diff --git a/src/intercycle.F90 b/src/intercycle.F90 index b5df72bc0..e95cda141 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -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 !=============================================================================== diff --git a/src/output.F90 b/src/output.F90 index fb390f005..814c0f2e3 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/physics.F90 b/src/physics.F90 index 9e6bffef3..0830b17e9 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -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 diff --git a/src/search.F90 b/src/search.F90 index 46bbd89f0..3cdb51bfc 100644 --- a/src/search.F90 +++ b/src/search.F90 @@ -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 diff --git a/src/source.F90 b/src/source.F90 index 7fd2a6229..55809a1db 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -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