Combined main batch/particle loop structure for eigenvalue and fixed source

simulations in one 'simulation' mode. The fixed_source module is now gone, and
the eigenvalue module consists specifically of eigenvalue-related subroutines.
This commit is contained in:
Paul Romano 2015-07-02 20:01:24 +07:00
parent 2c8c771998
commit 14c62e70b6
5 changed files with 331 additions and 433 deletions

View file

@ -4,275 +4,24 @@ module eigenvalue
use message_passing
#endif
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use error, only: fatal_error, warning
use global
use math, only: t_percentile
use mesh, only: count_bank_sites
use mesh_header, only: StructuredMesh
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use particle_header, only: Particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use source, only: get_source_particle, initialize_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
implicit none
private
public :: run_eigenvalue
real(8) :: keff_generation ! Single-generation k on each
! processor
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
real(8) :: keff_generation ! Single-generation k on each processor
real(8) :: k_sum(2) = ZERO ! Used to reduce sum and sum_sq
contains
!===============================================================================
! RUN_EIGENVALUE encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a k-eigenvalue calculation.
!===============================================================================
subroutine run_eigenvalue()
type(Particle) :: p
integer(8) :: i_work
if (.not. restart_run) call initialize_source()
if (master) call header("K EIGENVALUE SIMULATION", level=1)
! Display column titles
if(master) call print_columns()
! Turn on inactive timer
call time_inactive % start()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
call initialize_batch()
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
cycle BATCH_LOOP
end if
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
call initialize_generation()
! Start timer for transport
call time_transport % start()
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call get_source_particle(p, current_work)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_generation()
end do GENERATION_LOOP
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
! Clear particle
call p % clear()
end subroutine run_eigenvalue
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
call write_message("Simulating batch " // trim(to_str(current_batch)) &
&// "...", 8)
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
if (current_batch == n_inactive + 1) then
! Switch from inactive batch timer to active batch timer
call time_inactive % stop()
call time_active % start()
! Enable active batches (and tallies_on if it hasn't been enabled)
active_batches = .true.
tallies_on = .true.
! Add user tallies to active tallies list
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
end if
! check CMFD initialize batch
if (cmfd_run) call cmfd_init_batch()
end subroutine initialize_batch
!===============================================================================
! INITIALIZE_GENERATION
!===============================================================================
subroutine initialize_generation()
! set overall generation number
overall_gen = gen_per_batch*(current_batch - 1) + current_gen
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call count_source_for_ufs()
! Store current value of tracklength k
keff_generation = global_tallies(K_TRACKLENGTH) % value
end subroutine initialize_generation
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
subroutine finalize_generation()
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + global_tally_collision
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + global_tally_absorption
!$omp end critical
! reset private tallies
global_tally_tracklength = ZERO
global_tally_collision = ZERO
global_tally_leakage = ZERO
global_tally_absorption = ZERO
!$omp end parallel
#ifdef _OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()
#endif
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
call time_bank % stop()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Collect results and statistics
call calculate_generation_keff()
call calculate_average_keff()
! Write generation output
if (master .and. current_gen /= gen_per_batch) call print_generation()
end subroutine finalize_generation
!===============================================================================
! FINALIZE_BATCH handles synchronization and accumulation of tallies,
! calculation of Shannon entropy, getting single-batch estimate of keff, and
! turning on tallies when appropriate
!===============================================================================
subroutine finalize_batch()
! Collect tallies
call time_tallies % start()
call synchronize_tallies()
call time_tallies % stop()
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
n_realizations = 0
end if
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Display output
if (master) call print_batch_keff()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()
end if
! Write out source point if it's been specified for this batch
if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. &
source_write) then
call write_source_point()
end if
if (master .and. current_batch == n_max_batches) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
end if
end subroutine finalize_batch
!===============================================================================
! SYNCHRONIZE_BANK samples source sites from the fission sites that were
! accumulated during the generation. This routine is what allows this Monte
@ -832,37 +581,6 @@ contains
end subroutine count_source_for_ufs
!===============================================================================
! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a
! batch using data read from a state point file
!===============================================================================
subroutine replay_batch_history
! Write message at beginning
if (current_batch == 1) then
call write_message("Replaying history from state point...", 1)
end if
do current_gen = 1, gen_per_batch
overall_gen = overall_gen + 1
call calculate_average_keff()
! print out batch keff
if (current_gen < gen_per_batch) then
if (master) call print_generation()
else
if (master) call print_batch_keff()
end if
end do
! Write message at end
if (current_batch == restart_batch) then
call write_message("Resuming simulation...", 1)
end if
end subroutine replay_batch_history
#ifdef _OPENMP
!===============================================================================
! JOIN_BANK_FROM_THREADS

View file

@ -1,142 +0,0 @@
module fixed_source
#ifdef MPI
use message_passing
#endif
use constants, only: ZERO, MAX_LINE_LEN
use global
use output, only: write_message, header
use particle_header, only: Particle
use random_lcg, only: set_particle_seed
use source, only: initialize_source, get_source_particle
use state_point, only: write_state_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies
use trigger, only: check_triggers
use tracking, only: transport
implicit none
contains
subroutine run_fixedsource()
type(Particle) :: p
integer(8) :: i_work ! index over histories in single cycle
if (.not. restart_run) call initialize_source()
if (master) call header("FIXED SOURCE TRANSPORT SIMULATION", level=1)
! Turn timer and tallies on
tallies_on = .true.
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
call time_active % start()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
! In a restart run, skip any batches that have already been simulated
if (restart_run .and. current_batch <= restart_batch) then
if (current_batch > n_inactive) n_realizations = n_realizations + 1
cycle BATCH_LOOP
end if
call initialize_batch()
overall_gen = current_batch
! Start timer for transport
call time_transport % start()
! =======================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call get_source_particle(p, current_work)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
end subroutine run_fixedsource
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
call write_message("Simulating batch " // trim(to_str(current_batch)) &
&// "...", 1)
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
end subroutine initialize_batch
!===============================================================================
! FINALIZE_BATCH
!===============================================================================
subroutine finalize_batch()
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
!$omp end critical
! reset private tallies
global_tally_leakage = ZERO
!$omp end parallel
! Collect and accumulate tallies
call time_tallies % start()
call synchronize_tallies()
call time_tallies % stop()
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()
end if
end subroutine finalize_batch
end module fixed_source

View file

@ -1,13 +1,12 @@
program main
use constants
use eigenvalue, only: run_eigenvalue
use finalize, only: finalize_run
use fixed_source, only: run_fixedsource
use global
use initialize, only: initialize_run
use particle_restart, only: run_particle_restart
use plot, only: run_plot
use simulation, only: run_simulation
implicit none
@ -16,10 +15,8 @@ program main
! start problem based on mode
select case (run_mode)
case (MODE_FIXEDSOURCE)
call run_fixedsource()
case (MODE_EIGENVALUE)
call run_eigenvalue()
case (MODE_FIXEDSOURCE, MODE_EIGENVALUE)
call run_simulation()
case (MODE_PLOTTING)
call run_plot()
case (MODE_PARTICLE)

325
src/simulation.F90 Normal file
View file

@ -0,0 +1,325 @@
module simulation
#ifdef MPI
use mpi
#endif
use cmfd_execute, only: cmfd_init_batch, execute_cmfd
use constants, only: ZERO
use eigenvalue, only: count_source_for_ufs, calculate_average_keff, &
calculate_combined_keff, calculate_generation_keff, &
shannon_entropy, synchronize_bank, keff_generation
#ifdef _OPENMP
use eigenvalue, only: join_bank_from_threads
#endif
use global
use output, only: write_message, header, print_columns, &
print_batch_keff, print_generation
use particle_header, only: Particle
use source, only: get_source_particle, initialize_source
use state_point, only: write_state_point, write_source_point
use string, only: to_str
use tally, only: synchronize_tallies, setup_active_usertallies, &
reset_result
use trigger, only: check_triggers
use tracking, only: transport
implicit none
private
public :: run_simulation
contains
!===============================================================================
! RUN_EIGENVALUE encompasses all the main logic where iterations are performed
! over the batches, generations, and histories in a k-eigenvalue calculation.
!===============================================================================
subroutine run_simulation()
type(Particle) :: p
integer(8) :: i_work
if (.not. restart_run) call initialize_source()
! Display header
if (master) then
if (run_mode == MODE_FIXEDSOURCE) then
call header("FIXED SOURCE TRANSPORT SIMULATION", level=1)
elseif (run_mode == MODE_EIGENVALUE) then
call header("K EIGENVALUE SIMULATION", level=1)
call print_columns()
end if
end if
! Turn on inactive timer
call time_inactive % start()
! ==========================================================================
! LOOP OVER BATCHES
BATCH_LOOP: do current_batch = 1, n_max_batches
call initialize_batch()
! Handle restart runs
if (restart_run .and. current_batch <= restart_batch) then
call replay_batch_history()
cycle BATCH_LOOP
end if
! =======================================================================
! LOOP OVER GENERATIONS
GENERATION_LOOP: do current_gen = 1, gen_per_batch
call initialize_generation()
! Start timer for transport
call time_transport % start()
! ====================================================================
! LOOP OVER PARTICLES
!$omp parallel do schedule(static) firstprivate(p)
PARTICLE_LOOP: do i_work = 1, work
current_work = i_work
! grab source particle from bank
call get_source_particle(p, current_work)
! transport particle
call transport(p)
end do PARTICLE_LOOP
!$omp end parallel do
! Accumulate time for transport
call time_transport % stop()
call finalize_generation()
end do GENERATION_LOOP
call finalize_batch()
if (satisfy_triggers) exit BATCH_LOOP
end do BATCH_LOOP
call time_active % stop()
! ==========================================================================
! END OF RUN WRAPUP
if (master) call header("SIMULATION FINISHED", level=1)
! Clear particle
call p % clear()
end subroutine run_simulation
!===============================================================================
! INITIALIZE_BATCH
!===============================================================================
subroutine initialize_batch()
if (run_mode == MODE_FIXEDSOURCE) then
call write_message("Simulating batch " // trim(to_str(current_batch)) &
// "...", 1)
end if
! Reset total starting particle weight used for normalizing tallies
total_weight = ZERO
if (current_batch == n_inactive + 1) then
! Switch from inactive batch timer to active batch timer
call time_inactive % stop()
call time_active % start()
! Enable active batches (and tallies_on if it hasn't been enabled)
active_batches = .true.
tallies_on = .true.
! Add user tallies to active tallies list
!$omp parallel
call setup_active_usertallies()
!$omp end parallel
end if
! check CMFD initialize batch
if (run_mode == MODE_EIGENVALUE) then
if (cmfd_run) call cmfd_init_batch()
end if
end subroutine initialize_batch
!===============================================================================
! INITIALIZE_GENERATION
!===============================================================================
subroutine initialize_generation()
! set overall generation number
overall_gen = gen_per_batch*(current_batch - 1) + current_gen
if (run_mode == MODE_EIGENVALUE) then
! Reset number of fission bank sites
n_bank = 0
! Count source sites if using uniform fission source weighting
if (ufs) call count_source_for_ufs()
! Store current value of tracklength k
keff_generation = global_tallies(K_TRACKLENGTH) % value
end if
end subroutine initialize_generation
!===============================================================================
! FINALIZE_GENERATION
!===============================================================================
subroutine finalize_generation()
! Update global tallies with the omp private accumulation variables
!$omp parallel
!$omp critical
if (run_mode == MODE_EIGENVALUE) then
global_tallies(K_COLLISION) % value = &
global_tallies(K_COLLISION) % value + global_tally_collision
global_tallies(K_ABSORPTION) % value = &
global_tallies(K_ABSORPTION) % value + global_tally_absorption
global_tallies(K_TRACKLENGTH) % value = &
global_tallies(K_TRACKLENGTH) % value + global_tally_tracklength
end if
global_tallies(LEAKAGE) % value = &
global_tallies(LEAKAGE) % value + global_tally_leakage
!$omp end critical
! reset private tallies
if (run_mode == MODE_EIGENVALUE) then
global_tally_collision = 0
global_tally_absorption = 0
global_tally_tracklength = 0
end if
global_tally_leakage = 0
!$omp end parallel
if (run_mode == MODE_EIGENVALUE) then
#ifdef _OPENMP
! Join the fission bank from each thread into one global fission bank
call join_bank_from_threads()
#endif
! Distribute fission bank across processors evenly
call time_bank % start()
call synchronize_bank()
call time_bank % stop()
! Calculate shannon entropy
if (entropy_on) call shannon_entropy()
! Collect results and statistics
call calculate_generation_keff()
call calculate_average_keff()
! Write generation output
if (master .and. current_gen /= gen_per_batch) call print_generation()
end if
end subroutine finalize_generation
!===============================================================================
! FINALIZE_BATCH handles synchronization and accumulation of tallies,
! calculation of Shannon entropy, getting single-batch estimate of keff, and
! turning on tallies when appropriate
!===============================================================================
subroutine finalize_batch()
! Collect tallies
call time_tallies % start()
call synchronize_tallies()
call time_tallies % stop()
! Reset global tally results
if (.not. active_batches) then
call reset_result(global_tallies)
n_realizations = 0
end if
if (run_mode == MODE_EIGENVALUE) then
! Perform CMFD calculation if on
if (cmfd_on) call execute_cmfd()
! Display output
if (master) call print_batch_keff()
! Calculate combined estimate of k-effective
if (master) call calculate_combined_keff()
end if
! Check_triggers
if (master) call check_triggers()
#ifdef MPI
call MPI_BCAST(satisfy_triggers, 1, MPI_LOGICAL, 0, &
MPI_COMM_WORLD, mpi_err)
#endif
if (satisfy_triggers .or. &
(trigger_on .and. current_batch == n_max_batches)) then
call statepoint_batch % add(current_batch)
end if
! Write out state point if it's been specified for this batch
if (statepoint_batch % contains(current_batch)) then
call write_state_point()
end if
! Write out source point if it's been specified for this batch
if ((sourcepoint_batch % contains(current_batch) .or. source_latest) .and. &
source_write) then
call write_source_point()
end if
if (master .and. current_batch == n_max_batches .and. &
run_mode == MODE_EIGENVALUE) then
! Make sure combined estimate of k-effective is calculated at the last
! batch in case no state point is written
call calculate_combined_keff()
end if
end subroutine finalize_batch
!===============================================================================
! REPLAY_BATCH_HISTORY displays keff and entropy for each generation within a
! batch using data read from a state point file
!===============================================================================
subroutine replay_batch_history
! Write message at beginning
if (current_batch == 1) then
call write_message("Replaying history from state point...", 1)
end if
if (run_mode == MODE_EIGENVALUE) then
do current_gen = 1, gen_per_batch
overall_gen = overall_gen + 1
call calculate_average_keff()
! print out batch keff
if (current_gen < gen_per_batch) then
if (master) call print_generation()
else
if (master) call print_batch_keff()
end if
end do
end if
! Write message at end
if (current_batch == restart_batch) then
call write_message("Resuming simulation...", 1)
end if
end subroutine replay_batch_history
end module simulation

View file

@ -1,7 +1,7 @@
current batch:
7.000000E+00
current gen:
0.000000E+00
1.000000E+00
particle id:
9.280000E+02
run mode: