From 97b838fab45ae4832f6b1899644b46ea071388bc Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 6 Feb 2012 14:40:34 -0500 Subject: [PATCH] Changed source bank to Bank type. Initializing of source particles handled cleaner. --- src/global.F90 | 5 +-- src/intercycle.F90 | 106 +++++++++++++++++++++++++-------------------- src/main.F90 | 22 +++------- src/source.F90 | 93 ++++++++++++++++++++------------------- 4 files changed, 115 insertions(+), 111 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index 93d2891879..925a2cc019 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -111,13 +111,12 @@ module global type(ExtSource), target :: external_source ! Source and fission bank - type(Particle), allocatable, target :: source_bank(:) - type(Bank), allocatable, target :: fission_bank(:) + type(Bank), allocatable, target :: source_bank(:) + type(Bank), allocatable, target :: fission_bank(:) integer(8) :: n_bank ! # of sites in fission bank 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) :: source_index ! index for source particles ! cycle keff real(8) :: keff = ONE diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 9a8d6f1217..b5df72bc01 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -28,9 +28,9 @@ contains 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_local ! index for source bank - integer :: send_to_left ! # of bank sites to send/recv to or from left - integer :: send_to_right ! # of bank sites to send/recv to or from right + 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 :: & @@ -48,8 +48,16 @@ contains 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 + ! of fission bank sites -- each processor needs to know the total number of + ! sites in order to figure out the probability for selecting + ! sites. Furthermore, each proc also needs to know where in the 'global' + ! fission bank its own sites starts in order to ensure reproducibility by + ! skipping ahead to the proper seed. + #ifdef MPI - ! Determine starting index for fission bank and total sites in fission bank start = 0_8 call MPI_EXSCAN(n_bank, start, 1, MPI_INTEGER8, MPI_SUM, & MPI_COMM_WORLD, mpi_err) @@ -64,20 +72,25 @@ contains total = n_bank #endif - ! Check if there are no fission sites + ! 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. + if (total == 0) then message = "No fission sites banked!" call fatal_error() end if - ! Make sure all processors start at the same point for random sampling - call set_particle_seed(int(current_cycle,8)) + ! Make sure all processors start at the same point for random sampling by + ! using the current_cycle as the starting seed. Then skip ahead in the + ! sequence using the starting index in the 'global' fission bank for each + ! processor Skip ahead however many random numbers are needed - ! Skip ahead however many random numbers are needed + call set_particle_seed(int(current_cycle,8)) call prn_skip(start) - allocate(temp_sites(2*work)) - index_local = 0_8 ! Index for local source_bank + ! Determine how many fission sites we need to sample from the source bank + ! and the probability for selecting a site. if (total < n_particles) then sites_needed = mod(n_particles,total) @@ -93,6 +106,12 @@ contains ! ========================================================================== ! 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)) + do i = 1, int(n_bank,4) ! If there are less than n_particles particles banked, automatically add @@ -100,34 +119,35 @@ contains ! 1000 and 300 were banked, this would add 3 source sites per banked site ! and the remaining 100 would be randomly sampled. if (total < n_particles) then - do j = 1,int(n_particles/total) - ! If index is within this node's range, add site to source - index_local = index_local + 1 - temp_sites(index_local) = fission_bank(i) + do j = 1, int(n_particles/total) + index_temp = index_temp + 1 + temp_sites(index_temp) = fission_bank(i) end do end if ! Randomly sample sites needed if (prn() < p_sample) then - index_local = index_local + 1 - temp_sites(index_local) = fission_bank(i) + index_temp = index_temp + 1 + temp_sites(index_temp) = fission_bank(i) 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 + #ifdef MPI start = 0_8 - call MPI_EXSCAN(index_local, start, 1, MPI_INTEGER8, MPI_SUM, & + call MPI_EXSCAN(index_temp, start, 1, MPI_INTEGER8, MPI_SUM, & MPI_COMM_WORLD, mpi_err) - finish = start + index_local + finish = start + index_temp total = finish call MPI_BCAST(total, 1, MPI_INTEGER8, n_procs - 1, & MPI_COMM_WORLD, mpi_err) #else start = 0_8 - finish = index_local - total = index_local + finish = index_temp + total = index_temp #endif ! Determine how many sites to send to adjacent nodes @@ -146,7 +166,7 @@ contains ! If we have extra sites sampled, we will simply discard the extra ! ones on the last processor if (rank == n_procs - 1) then - index_local = index_local - send_to_right + index_temp = index_temp - send_to_right end if elseif (total < n_particles) then @@ -154,8 +174,8 @@ contains ! fission bank sites_needed = n_particles - total do i = 1, int(sites_needed,4) - index_local = index_local + 1 - temp_sites(index_local) = fission_bank(n_bank - sites_needed + i) + index_temp = index_temp + 1 + temp_sites(index_temp) = fission_bank(n_bank - sites_needed + i) end do end if @@ -176,7 +196,7 @@ contains allocate(right_bank(abs(send_to_right))) if (send_to_right > 0) then - i = index_local - send_to_right + 1 + 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 @@ -202,29 +222,29 @@ contains ! ========================================================================== ! 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_local,4) - send_to_right ! size of second block - call copy_from_bank(temp_sites, i+1, j) + 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_local,4) - send_to_left ! size of first block + i = int(index_temp,4) - send_to_left ! size of first block j = -send_to_right ! size of second block - call copy_from_bank(temp_sites(1+send_to_left), 1, i) + 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_local,4) - send_to_left - send_to_right - call copy_from_bank(temp_sites(1+send_to_left), 1, i) + 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_local,4) + j = int(index_temp,4) k = -send_to_right - call copy_from_bank(temp_sites, i+1, j) + source_bank(i+1:i+j) = temp_sites(1:j) #ifdef MPI call MPI_WAIT(request_left, status, mpi_err) #endif @@ -235,9 +255,6 @@ contains call copy_from_bank(right_bank, i+j+1, k) end if - ! Reset source index - source_index = 0_8 - call timer_stop(time_ic_rebuild) #ifdef MPI @@ -260,21 +277,16 @@ contains integer :: i ! index in temp_bank integer :: i_source ! index in source_bank - type(Particle), pointer :: p + type(Bank), pointer :: src => null() do i = 1, n_sites i_source = i_start + i - 1 - p => source_bank(i_source) - - ! set defaults - call initialize_particle(p) - - p % coord % xyz = temp_bank(i) % xyz - p % coord % uvw = temp_bank(i) % uvw - p % last_xyz = temp_bank(i) % xyz - p % E = temp_bank(i) % E - p % last_E = temp_bank(i) % E + 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 diff --git a/src/main.F90 b/src/main.F90 index 5d61c0b3e7..db84cbc79d 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -43,7 +43,7 @@ contains subroutine run_problem() - integer(8) :: particle_seed ! unique index for particle + integer(8) :: i ! index over histories in single cycle type(Particle), pointer :: p => null() if (master) call header("BEGIN SIMULATION", level=1) @@ -51,6 +51,9 @@ contains tallies_on = .false. call timer_start(time_inactive) + ! Allocate particle + allocate(p) + ! Display column titles if (entropy_on) then message = " Cycle k(cycle) Entropy Average k" @@ -79,23 +82,10 @@ contains ! ======================================================================= ! LOOP OVER HISTORIES - HISTORY_LOOP: do + HISTORY_LOOP: do i = 1, work ! grab source particle from bank - p => get_source_particle() - if ( .not. associated(p) ) then - ! no particles left in source bank - exit HISTORY_LOOP - end if - - ! set random number seed - particle_seed = (current_cycle - 1)*n_particles + p % id - call set_particle_seed(particle_seed) - - ! set particle trace - trace = .false. - if (current_cycle == trace_cycle .and. & - p % id == trace_particle) trace = .true. + call get_source_particle(p, i) ! transport particle call transport(p) diff --git a/src/source.F90 b/src/source.F90 index 1ecea63eda..7fd2a6229e 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -21,20 +21,22 @@ contains subroutine initialize_source() - type(Particle), pointer :: p => null() 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 :: alloc_err ! allocation error code integer(8) :: bytes ! size of fission/source bank + integer(8) :: id ! particle id + integer :: alloc_err ! allocation error code real(8) :: r(3) ! sampled coordinates real(8) :: phi ! azimuthal angle real(8) :: mu ! cosine of polar angle real(8) :: E ! outgoing energy real(8) :: p_min(3) ! minimum coordinates of source real(8) :: p_max(3) ! maximum coordinates of source +#ifndef NO_F2008 type(Bank) :: bank_obj +#endif message = "Initializing source particles..." call write_message(6) @@ -68,47 +70,38 @@ contains call fatal_error() end if - ! Check external source type - if (external_source%type == SRC_BOX) then - p_min = external_source%values(1:3) - p_max = external_source%values(4:6) - else - message = "Unsupported external source type: " // & - to_str(external_source%type) - call fatal_error() - end if - ! Initialize first cycle source bank do i = 0, n_procs - 1 if (rank == i) then ! ID's of first and last source particles - bank_first = i*maxwork + 1 - bank_last = min((i+1)*maxwork, n_particles) + bank_first = rank*maxwork + 1 + bank_last = min((rank+1)*maxwork, n_particles) ! number of particles for this processor work = bank_last - bank_first + 1 - do j = bank_first, bank_last - p => source_bank(j - bank_first + 1) - - ! set defaults - call initialize_particle(p) + do j = 1, work + id = bank_first + j - 1 + source_bank(j) % id = id ! initialize random number seed - call set_particle_seed(int(j,8)) + call set_particle_seed(id) - ! sample position - r = (/ (prn(), k = 1,3) /) - p % id = j - p % coord0 % xyz = p_min + r*(p_max - p_min) - p % last_xyz = p % coord0 % xyz + ! sample position from external source + select case (external_source % type) + case (SRC_BOX) + p_min = external_source % values(1:3) + p_max = external_source % values(4:6) + r = (/ (prn(), k = 1,3) /) + source_bank(j) % xyz = p_min + r*(p_max - p_min) + end select ! sample angle phi = TWO*PI*prn() mu = TWO*prn() - ONE - p % coord0 % uvw(1) = mu - p % coord0 % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) - p % coord0 % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) + source_bank(j) % uvw(1) = mu + source_bank(j) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi) + source_bank(j) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi) ! sample energy from Watt fission energy spectrum for U-235 do @@ -118,40 +111,50 @@ contains end do ! set particle energy - p % E = E - p % last_E = E + source_bank(j) % E = E end do end if end do - ! Reset source index - source_index = 0_8 - end subroutine initialize_source !=============================================================================== ! GET_SOURCE_PARTICLE returns the next source particle !=============================================================================== - function get_source_particle() result(p) + subroutine get_source_particle(p, index_source) type(Particle), pointer :: p + integer(8), intent(in) :: index_source - ! increment index - source_index = source_index + 1 + integer(8) :: particle_seed ! unique index for particle + type(Bank), pointer :: src => null() - ! if at end of bank, return null pointer - if (source_index > work) then - p => null() - return - end if + ! set defaults + call initialize_particle(p) ! point to next source particle - p => source_bank(source_index) + src => source_bank(index_source) - ! set id - p % id = bank_first + source_index - 1 + ! copy attributes from source bank site + p % coord % xyz = src % xyz + p % coord % uvw = src % uvw + p % last_xyz = src % xyz + p % E = src % E + p % last_E = src % E - end function get_source_particle + ! set identifier for particle + p % id = bank_first + index_source - 1 + + ! set random number seed + particle_seed = (current_cycle - 1)*n_particles + p % id + call set_particle_seed(particle_seed) + + ! set particle trace + trace = .false. + if (current_cycle == trace_cycle .and. & + p % id == trace_particle) trace = .true. + + end subroutine get_source_particle end module source