diff --git a/src/error.F90 b/src/error.F90 index 4f4c53d78a..59cc6d80e0 100644 --- a/src/error.F90 +++ b/src/error.F90 @@ -76,9 +76,10 @@ contains write(ERROR_UNIT,*) ! Write information on current cycle and particle - if (current_cycle > 0) then - write(ERROR_UNIT,'(1X,A,I11) ') 'Cycle: ', current_cycle - write(ERROR_UNIT,'(1X,A,I11)') 'Particle: ', p % id + if (current_batch > 0) then + write(ERROR_UNIT,'(1X,A,I12) ') 'Batch: ', current_batch + write(ERROR_UNIT,'(1X,A,I12) ') 'Generation:', current_gen + write(ERROR_UNIT,'(1X,A,I12)') 'Particle: ', p % id write(ERROR_UNIT,*) end if diff --git a/src/global.F90 b/src/global.F90 index bd4b7471c4..d7edf18215 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -125,11 +125,13 @@ module global ! ============================================================================ ! CRITICALITY SIMULATION VARIABLES - integer(8) :: n_particles = 10000 ! # of particles per cycle - integer :: n_cycles = 500 ! # of cycles - integer :: n_inactive = 50 ! # of inactive cycles - integer :: n_active ! # of active cycles - integer :: current_cycle = 0 ! current cycle + integer(8) :: n_particles ! # of particles per generation + integer :: n_batches ! # of batches + integer :: n_inactive ! # of inactive batches + integer :: n_active ! # of active batches + integer :: gen_per_batch = 1 ! # of generations per batch + integer :: current_batch = 0 ! current batch + integer :: current_gen = 0 ! current generation ! External source type(ExtSource), target :: external_source @@ -143,7 +145,7 @@ module global integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor - ! cycle keff + ! single-genreation keff real(8) :: keff = ONE real(8) :: keff_std @@ -217,7 +219,8 @@ module global ! Trace for single particle logical :: trace - integer :: trace_cycle + integer :: trace_batch + integer :: trace_gen integer(8) :: trace_particle contains diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index fafe30109c..cfc7383291 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -77,19 +77,22 @@ contains call hdf5_make_double(hdf5_output_file, "n_particles", real(n_particles,8)) ! Use H5LT interface to write n_cycles, n_inactive, and n_active - call hdf5_make_integer(hdf5_output_file, "n_cycles", n_cycles) + call hdf5_make_integer(hdf5_output_file, "n_batches", n_batches) call hdf5_make_integer(hdf5_output_file, "n_inactive", n_inactive) call hdf5_make_integer(hdf5_output_file, "n_active", n_active) + call hdf5_make_integer(hdf5_output_file, "gen_per_batch", gen_per_batch) ! Add description of each variable call h5ltset_attribute_string_f(hdf5_output_file, "n_particles", & "description", "Number of particles per cycle", hdf5_err) - call h5ltset_attribute_string_f(hdf5_output_file, "n_cycles", & - "description", "Total number of cycles", hdf5_err) + call h5ltset_attribute_string_f(hdf5_output_file, "n_batches", & + "description", "Total number of batches", hdf5_err) call h5ltset_attribute_string_f(hdf5_output_file, "n_inactive", & "description", "Number of inactive cycles", hdf5_err) call h5ltset_attribute_string_f(hdf5_output_file, "n_active", & "description", "Number of active cycles", hdf5_err) + call h5ltset_attribute_string_f(hdf5_output_file, "gen_per_batch", & + "description", "Number of generations per batch", hdf5_err) end if call hdf5_write_geometry() @@ -788,7 +791,7 @@ contains call h5ltset_attribute_string_f(timing_group, "time_transport", & "description", "Time in transport only (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_intercycle", & - "description", "Total time between cycles (s)", hdf5_err) + "description", "Total time between generations (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_tallies", & "description", "Time between cycles accumulating tallies (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_sample", & @@ -796,16 +799,16 @@ contains 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_inactive", & - "description", "Total time in inactive cycles (s)", hdf5_err) + "description", "Total time in inactive batches (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_active", & - "description", "Total time in active cycles (s)", hdf5_err) + "description", "Total time in active batches (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_finalize", & "description", "Total time for finalization (s)", hdf5_err) call h5ltset_attribute_string_f(timing_group, "time_total", & "description", "Total time elapsed (s)", hdf5_err) ! Write calculation rate - total_particles = n_particles * n_cycles + total_particles = n_particles * n_batches * gen_per_batch speed = real(total_particles) / (time_inactive % elapsed + & time_active % elapsed) call hdf5_make_double(timing_group, "neutrons_per_second", speed) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3f74d514a6..29c72309a9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -91,7 +91,7 @@ contains end if ! Criticality information - if (criticality % cycles > 0) then + if (criticality % batches > 0) then problem_type = PROB_CRITICALITY ! Check number of particles @@ -102,9 +102,13 @@ contains n_particles = str_to_int(criticality % particles) ! Copy cycle information - n_cycles = criticality % cycles - n_inactive = criticality % inactive - n_active = n_cycles - n_inactive + n_batches = criticality % batches + n_inactive = criticality % inactive + n_active = n_batches - n_inactive + gen_per_batch = criticality % generations_per_batch + else + message = "Need to specify number of batches with tag." + call fatal_error() end if ! Verbosity @@ -151,8 +155,9 @@ contains ! Particle trace if (associated(trace_)) then - trace_cycle = trace_(1) - trace_particle = trace_(2) + trace_batch = trace_(1) + trace_gen = trace_(2) + trace_particle = trace_(3) end if ! Shannon Entropy mesh diff --git a/src/intercycle.F90 b/src/intercycle.F90 index 5ecd79fe84..6d3ae131e8 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -10,7 +10,7 @@ module intercycle use random_lcg, only: prn, set_particle_seed, prn_skip use search, only: binary_search use string, only: to_str - use tally, only: accumulate_cycle_estimate + use tally, only: accumulate_batch_estimate use tally_header, only: TallyObject use timing, only: timer_start, timer_stop @@ -22,8 +22,8 @@ contains !=============================================================================== ! SYNCHRONIZE_BANK samples source sites from the fission sites that were -! accumulated during the cycle. This routine is what allows this Monte Carlo to -! scale to large numbers of processors where other codes cannot. +! accumulated during the generation. This routine is what allows this Monte +! Carlo to scale to large numbers of processors where other codes cannot. !=============================================================================== subroutine synchronize_bank() @@ -73,7 +73,7 @@ contains total = n_bank #endif - ! If there are not that many particles per cycle, it's possible that no + ! If there are not that many particles per generation, it's possible that no ! 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. @@ -83,12 +83,12 @@ contains call fatal_error() end if - ! 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 + ! Make sure all processors start at the same point for random sampling. Then + ! skip ahead in the sequence using the starting index in the 'global' + ! fission bank for each processor. - call set_particle_seed(int(current_cycle,8)) + call set_particle_seed(int((current_batch - 1)*gen_per_batch + & + current_gen,8)) call prn_skip(start) ! Determine how many fission sites we need to sample from the source bank @@ -257,20 +257,22 @@ contains ! 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 + ! the next generation 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) + ! Deallocate space for bank_position on the very last generation + if (current_batch == n_batches .and. current_gen == gen_per_batch) & + deallocate(bank_position) #else source_bank = temp_sites(1:n_particles) #endif call timer_stop(time_ic_sendrecv) - ! Deallocate space for the temporary source bank on the last cycle - if (current_cycle == n_cycles) deallocate(temp_sites) + ! Deallocate space for the temporary source bank on the last generation + if (current_batch == n_batches .and. current_gen == gen_per_batch) & + deallocate(temp_sites) end subroutine synchronize_bank @@ -379,27 +381,22 @@ contains end subroutine shannon_entropy !=============================================================================== -! CALCULATE_KEFF calculates the single cycle estimate of keff as well as the -! mean and standard deviation of the mean for active cycles and displays them +! CALCULATE_KEFF calculates the single batch estimate of keff as well as the +! mean and standard deviation of the mean for active batches and displays them !=============================================================================== subroutine calculate_keff() - integer :: n ! active cycle number - real(8) :: k_cycle ! single cycle estimate of keff + integer :: n ! active batch number + real(8) :: k_batch ! single batch estimate of keff #ifdef MPI real(8) :: global_temp(N_GLOBAL_TALLIES) #endif - message = "Calculate cycle keff..." + message = "Calculate batch keff..." call write_message(8) - ! Since the creation of bank sites was originally weighted by the last - ! cycle keff, we need to multiply by that keff to get the current cycle's - ! value #ifdef MPI - global_tallies(K_ANALOG) % value = n_bank * keff - ! Copy global tallies into array to be reduced global_temp = global_tallies(:) % value @@ -416,49 +413,52 @@ contains ! Reset value on other processors global_tallies(:) % value = ZERO end if -#else - global_tallies(K_ANALOG) % value = n_bank * keff #endif ! Collect statistics and print output if (master) then - k_cycle = global_tallies(K_ANALOG) % value/n_particles + k_batch = global_tallies(K_ANALOG) % value/(n_particles*gen_per_batch) - if (current_cycle > n_inactive) then - ! Active cycle number - n = current_cycle - n_inactive + if (tallies_on) then + ! Active batch number + n = current_batch - n_inactive - ! Accumulate single cycle realizations of k - call accumulate_cycle_estimate(global_tallies) + ! Accumulate single batch realizations of k + call accumulate_batch_estimate(global_tallies) ! Determine mean and standard deviation of mean keff = global_tallies(K_ANALOG) % sum/n keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - keff*keff)/n) - ! Display output for this cycle - if (current_cycle > n_inactive + 1) then + ! Display output for this batch + if (current_batch > n_inactive + 1) then if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=103) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, & entropy, keff, keff_std else - write(UNIT=OUTPUT_UNIT, FMT=101) current_cycle, k_cycle, & + write(UNIT=OUTPUT_UNIT, FMT=101) current_batch, k_batch, & keff, keff_std end if else if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy + write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else - write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle + write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if end if else - ! Display output for inactive cycle + ! Display output for inactive batch if (entropy_on) then - write(UNIT=OUTPUT_UNIT, FMT=102) current_cycle, k_cycle, entropy + write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else - write(UNIT=OUTPUT_UNIT, FMT=100) current_cycle, k_cycle + write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if - keff = k_cycle + + ! Set keff + keff = k_batch + + ! Reset tally values + global_tallies(:) % value = ZERO end if end if diff --git a/src/main.F90 b/src/main.F90 index af9824a6a0..31092f79f6 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -11,7 +11,7 @@ program main use random_lcg, only: set_particle_seed use source, only: get_source_particle use string, only: to_str - use tally, only: synchronize_tallies + use tally, only: synchronize_tallies, add_to_score use timing, only: timer_start, timer_stop #ifdef MPI @@ -54,51 +54,68 @@ contains ! Display column titles if (entropy_on) then - message = " Cycle k(cycle) Entropy Average k" + message = " Batch k(batch) Entropy Average k" call write_message(1) message = " ===== ======== ======= ===================" call write_message(1) else - message = " Cycle k(cycle) Average k" + message = " Batch k(batch) Average k" call write_message(1) message = " ===== ======== ===================" call write_message(1) end if ! ========================================================================== - ! LOOP OVER CYCLES - CYCLE_LOOP: do current_cycle = 1, n_cycles + ! LOOP OVER BATCHES + BATCH_LOOP: do current_batch = 1, n_batches - message = "Simulating cycle " // trim(to_str(current_cycle)) // "..." + message = "Simulating batch " // trim(to_str(current_batch)) // "..." call write_message(8) - ! Set all tallies to zero - n_bank = 0 - ! ======================================================================= - ! LOOP OVER HISTORIES + ! LOOP OVER GENERATIONS + GENERATION_LOOP: do current_gen = 1, gen_per_batch - ! Start timer for transport - call timer_start(time_transport) + ! Set all tallies to zero + n_bank = 0 - HISTORY_LOOP: do i = 1, work + ! ==================================================================== + ! LOOP OVER HISTORIES - ! grab source particle from bank - call get_source_particle(i) + ! Start timer for transport + call timer_start(time_transport) - ! transport particle - call transport() + HISTORY_LOOP: do i = 1, work - end do HISTORY_LOOP + ! grab source particle from bank + call get_source_particle(i) - ! Accumulate time for transport - call timer_stop(time_transport) + ! transport particle + call transport() - ! ======================================================================= - ! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC + end do HISTORY_LOOP - ! Start timer for inter-cycle synchronization - call timer_start(time_intercycle) + ! Accumulate time for transport + call timer_stop(time_transport) + + ! ==================================================================== + ! WRAP UP FISSION BANK AND COMPUTE TALLIES, KEFF, ETC + + ! Start timer for inter-cycle synchronization + call timer_start(time_intercycle) + + ! Distribute fission bank across processors evenly + call synchronize_bank() + + ! Add to analog estimate of keff -- since the creation of bank sites + ! was originally weighted by the last cycle keff, we need to multiply + ! by that keff to get the current cycle's value + call add_to_score(global_tallies(K_ANALOG), n_bank * keff) + + ! Stop timer for inter-cycle synchronization + call timer_stop(time_intercycle) + + end do GENERATION_LOOP ! Collect tallies if (tallies_on) then @@ -110,25 +127,17 @@ contains ! Calculate shannon entropy if (entropy_on) call shannon_entropy() - ! Distribute fission bank across processors evenly - call synchronize_bank() - ! Collect results and statistics call calculate_keff() - ! print cycle information - ! Turn tallies on once inactive cycles are complete - if (current_cycle == n_inactive) then + if (current_batch == n_inactive) then tallies_on = .true. call timer_stop(time_inactive) call timer_start(time_active) end if - ! Stop timer for inter-cycle synchronization - call timer_stop(time_intercycle) - - end do CYCLE_LOOP + end do BATCH_LOOP call timer_stop(time_active) diff --git a/src/output.F90 b/src/output.F90 index 1a9141c0e0..920d9643f7 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -830,8 +830,9 @@ contains call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) if (problem_type == PROB_CRITICALITY) then write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' - write(UNIT_SUMMARY,101) 'Number of Cycles:', n_cycles - write(UNIT_SUMMARY,101) 'Number of Inactive Cycles:', n_inactive + write(UNIT_SUMMARY,101) 'Number of Batches:', n_batches + write(UNIT_SUMMARY,101) 'Number of Inactive Batches:', n_inactive + write(UNIT_SUMMARY,101) 'Generations per Batch:', gen_per_batch elseif (problem_type == PROB_SOURCE) then write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' end if @@ -950,9 +951,9 @@ contains write(ou,100) "Total time in simulation", time_inactive % elapsed + & time_active % elapsed write(ou,100) " Time in transport only", time_transport % elapsed - write(ou,100) " Time in inactive cycles", time_inactive % elapsed - write(ou,100) " Time in active cycles", time_active % elapsed - write(ou,100) " Time between cycles", time_intercycle % elapsed + write(ou,100) " Time in inactive batches", time_inactive % elapsed + write(ou,100) " Time in active batches", time_active % elapsed + write(ou,100) " Time between generations", time_intercycle % elapsed 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 @@ -960,7 +961,7 @@ contains write(ou,100) "Total time elapsed", time_total % elapsed ! display calculate rate and final keff - total_particles = n_particles * n_cycles + total_particles = n_particles * n_batches * gen_per_batch speed = real(total_particles) / (time_inactive % elapsed + & time_active % elapsed) string = to_str(speed) diff --git a/src/source.F90 b/src/source.F90 index fd9cbf1c68..943f18c52c 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -145,12 +145,13 @@ contains p % id = bank_first + index_source - 1 ! set random number seed - particle_seed = (current_cycle - 1)*n_particles + p % id + particle_seed = ((current_batch - 1)*gen_per_batch + & + current_gen - 1)*n_particles + p % id call set_particle_seed(particle_seed) ! set particle trace trace = .false. - if (current_cycle == trace_cycle .and. & + if (current_batch == trace_batch .and. current_gen == trace_gen .and. & p % id == trace_particle) trace = .true. end subroutine get_source_particle diff --git a/src/tally.F90 b/src/tally.F90 index af6d93e738..555e873ad2 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1048,28 +1048,31 @@ contains end subroutine add_to_score !=============================================================================== -! ACCUMULATE_CYCLE_ESTIMATE +! ACCUMULATE_BATCH_ESTIMATE !=============================================================================== - elemental subroutine accumulate_cycle_estimate(score) + elemental subroutine accumulate_batch_estimate(score) type(TallyScore), intent(inout) :: score + real(8) :: val + ! Add the sum and square of the sum of contributions from each cycle ! within a cycle to the variables sum and sum_sq. This will later allow us ! to calculate a variance on the tallies - score % sum = score % sum + score % value/n_particles - score % sum_sq = score % sum_sq + (score % value/n_particles)**2 + val = score % value/(n_particles*gen_per_batch) + score % sum = score % sum + val + score % sum_sq = score % sum_sq + val*val - ! Reset the single cycle estimate + ! Reset the single batch estimate score % value = ZERO - end subroutine accumulate_cycle_estimate + end subroutine accumulate_batch_estimate !=============================================================================== ! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history -! within the cycle to a new random variable +! within the batch to a new random variable !=============================================================================== subroutine synchronize_tallies() @@ -1086,7 +1089,7 @@ contains t => tallies(i) ! Loop over all filter and scoring bins - call accumulate_cycle_estimate(t % scores) + call accumulate_batch_estimate(t % scores) end do end subroutine synchronize_tallies diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 23247fe991..69388420ca 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -4,9 +4,10 @@ - + +