Added capability to due batching and multiple generations per batch. Not tested yet.

This commit is contained in:
Paul Romano 2012-03-15 17:50:01 -04:00
parent 2662a1b43b
commit 1e70be6f9b
10 changed files with 143 additions and 116 deletions

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

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

View file

@ -4,9 +4,10 @@
<options rootname="settings" />
<typedef name="criticality_xml">
<component name="cycles" type="integer" />
<component name="batches" type="integer" default="0" />
<component name="inactive" type="integer" />
<component name="particles" type="word" length="25" default="''" />
<component name="generations_per_batch" type="integer" default="1" />
</typedef>
<typedef name="source_xml">