From 60ca9a92f7589b98890062c05eebcf26a4fd8df6 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 6 Apr 2012 12:44:12 -0400 Subject: [PATCH 1/8] Added no_reduce option and introduced n_realizations and total_weight global variables. No_reduce does not yet work properly. --- src/global.F90 | 7 ++++++ src/input_xml.F90 | 13 +++++++++++ src/main.F90 | 3 +++ src/source.F90 | 3 +++ src/tally.F90 | 29 +++++++++++++++++++----- src/xml-fortran/templates/settings_t.xml | 1 + 6 files changed, 50 insertions(+), 6 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index 5da8782b5..a420d5b77 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -117,6 +117,10 @@ module global integer :: n_tracklength_tallies = 0 ! # of track-length tallies integer :: n_current_tallies = 0 ! # of surface current tallies + ! Normalization for statistics + integer :: n_realizations ! # of independent realizations + real(8) :: total_weight ! total starting particle weight in realization + ! Flag for turning tallies on logical :: tallies_on @@ -174,6 +178,9 @@ module global integer :: mpi_err ! MPI error code integer :: MPI_BANK ! MPI datatype for fission bank + ! No reduction at end of batch + logical :: no_reduce = .false. + ! ============================================================================ ! TIMING VARIABLES diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1b7f423c4..91bc79efc 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -68,6 +68,7 @@ contains energy_grid_ = "union" seed_ = 0_8 write_source_ = "" + no_reduce_ = "" source_ % type = "" ! Parse settings.xml file @@ -293,6 +294,18 @@ contains ! Check if the user has specified to write binary source file if (trim(write_source_) == 'on') write_source = .true. + ! Check if the user has specified to not reduce tallies at the end of every + ! batch + if (trim(no_reduce_) == 'on') no_reduce = .true. + + ! Determine number of realizations + if (no_reduce) then + n_realizations = n_active * n_procs + else + n_realizations = n_active + end if + + end subroutine read_settings_xml !=============================================================================== diff --git a/src/main.F90 b/src/main.F90 index 9cd90d6d0..d5906aaf0 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -72,6 +72,9 @@ contains message = "Simulating batch " // trim(to_str(current_batch)) // "..." call write_message(8) + ! Reset total starting weight + total_weight = ZERO + ! ======================================================================= ! LOOP OVER GENERATIONS GENERATION_LOOP: do current_gen = 1, gen_per_batch diff --git a/src/source.F90 b/src/source.F90 index 0c8caeb49..fe3b80fc8 100644 --- a/src/source.F90 +++ b/src/source.F90 @@ -172,6 +172,9 @@ contains if (current_batch == trace_batch .and. current_gen == trace_gen .and. & p % id == trace_particle) trace = .true. + ! Add paricle's starting weight to count for normalizing tallies later + total_weight = total_weight + src % wgt + end subroutine get_source_particle !=============================================================================== diff --git a/src/tally.F90 b/src/tally.F90 index f6d18eaa1..a90289a0d 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1310,8 +1310,12 @@ contains type(TallyObject), pointer :: t => null() #ifdef MPI - call reduce_tallies() - if (.not. master) return + if (no_reduce) then + if (current_batch == n_batches) call reduce_tallies() + else + call reduce_tallies() + if (.not. master) return + end if #endif do i = 1, n_tallies @@ -1369,6 +1373,19 @@ contains end do + ! We also need to determine the total starting weight of particles from the + ! last realization + if (.not. no_reduce) then + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, total_weight, 1, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + else + call MPI_REDUCE(total_weight, total_weight, 1, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + end if + end if + + end subroutine reduce_tallies #endif @@ -1808,7 +1825,7 @@ contains ! within a cycle to the variables sum and sum_sq. This will later allow us ! to calculate a variance on the tallies - val = score % value/(n_particles*gen_per_batch) + val = score % value/total_weight score % sum = score % sum + val score % sum_sq = score % sum_sq + val*val @@ -1830,9 +1847,9 @@ contains ! have used Bessel's correction so that the estimator of the variance of the ! sample mean is unbiased. - score % sum = score % sum/n_active - score % sum_sq = sqrt((score % sum_sq/n_active - score % sum**2) / & - (n_active - 1)) + score % sum = score % sum/n_realizations + score % sum_sq = sqrt((score % sum_sq/n_realizations - score % sum * & + score % sum) / (n_realizations - 1)) end subroutine statistics_score diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml index 321565955..ba648329d 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -31,6 +31,7 @@ + From 855e9c3a2818eaff623e8144ee9654cee4c70bdf Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Apr 2012 10:52:46 -0400 Subject: [PATCH 2/8] Significantly changed synchronize_tallies and calculate_keff to get no_reduce option to work. --- src/global.F90 | 7 ++- src/intercycle.F90 | 113 +++++++++++++++++++++++++------------- src/tally.F90 | 133 ++++++++++++++++++++++++++++++++++++++------- 3 files changed, 192 insertions(+), 61 deletions(-) diff --git a/src/global.F90 b/src/global.F90 index a420d5b77..7853ec917 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -150,9 +150,10 @@ module global integer(8) :: work ! number of particles per processor integer(8) :: maxwork ! maximum number of particles per processor - ! single-genreation keff - real(8) :: keff = ONE - real(8) :: keff_std + ! Temporary k-effective values + real(8) :: k_batch ! single batch estimate of k + real(8) :: keff = ONE ! average k over active cycles + real(8) :: keff_std ! standard deviation of average k ! Shannon entropy logical :: entropy_on = .false. diff --git a/src/intercycle.F90 b/src/intercycle.F90 index c7482c2d0..8c9d55e05 100644 --- a/src/intercycle.F90 +++ b/src/intercycle.F90 @@ -355,54 +355,84 @@ contains subroutine calculate_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 + integer :: n ! active batch number + real(8) :: temp(2) ! used to reduce sum and sum_sq message = "Calculate batch keff..." call write_message(8) + ! ========================================================================= + ! SINGLE-BATCH ESTIMATE OF K-EFFECTIVE + + if (.not. tallies_on) k_batch = global_tallies(K_ANALOG) % value + #ifdef MPI - ! Copy global tallies into array to be reduced - global_temp = global_tallies(:) % value - + ! Reduce value of k_batch if running in parallel if (master) then - call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) - - ! Transfer values back to global_tallies on master - global_tallies(:) % value = global_temp + call MPI_REDUCE(MPI_IN_PLACE, k_batch, 1, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) else - call MPI_REDUCE(global_temp, global_temp, N_GLOBAL_TALLIES, & - MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) - - ! Reset value on other processors - global_tallies(:) % value = ZERO + call MPI_REDUCE(k_batch, k_batch, 1, MPI_REAL8, MPI_SUM, 0, & + MPI_COMM_WORLD, mpi_err) end if #endif - ! Collect statistics and print output - if (master) then - k_batch = global_tallies(K_ANALOG) % value/(n_particles*gen_per_batch) + ! Normalize single batch estimate of k + k_batch = k_batch/(n_particles * gen_per_batch) - if (tallies_on) then - ! Active batch number + if (tallies_on) then + ! ======================================================================= + ! ACTIVE BATCHES + + if (no_reduce) then + ! In this case, no reduce was ever done on global_tallies. Thus, we + ! need to reduce the values in sum and sum^2 to get the sample mean + ! and its standard deviation + + ! Define number of realizations + n = (current_batch - n_inactive) * n_procs + +#ifdef MPI + if (current_batch /= n_batches) then + call MPI_REDUCE(global_tallies(K_ANALOG) % sum, temp, 2, & + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + else + temp(1) = global_tallies(K_ANALOG) % sum + temp(2) = global_tallies(K_ANALOG) % sum_sq + end if +#else + temp(1) = global_tallies(K_ANALOG) % sum + temp(2) = global_tallies(K_ANALOG) % sum_sq +#endif + + ! Sample mean of k + keff = temp(1) / n + + if (n > 1) then + ! Standard deviation of the sample mean of k + keff_std = sqrt((temp(2)/n - keff*keff)/(n - 1)) + end if + else + ! In this case, global_tallies has already been reduced, so we don't + ! need to perform any more reductions and just take the values from + ! global_tallies directly + + ! Define number of realizations n = current_batch - n_inactive - ! Accumulate single batch realizations of k - call accumulate_score(global_tallies) + ! Sample mean of keff + keff = global_tallies(K_ANALOG) % sum / n - ! Determine sample mean of keff - keff = global_tallies(K_ANALOG) % sum/n - - ! Display output for this batch - if (current_batch > n_inactive + 1) then - ! Determine standard deviation of the sample mean of keff + if (n > 1) then + ! Standard deviation of the sample mean of k keff_std = sqrt((global_tallies(K_ANALOG) % sum_sq/n - & keff*keff)/(n - 1)) + end if + end if + ! Display output for this batch + if (master) then + if (current_batch > n_inactive + 1) then if (entropy_on) then write(UNIT=OUTPUT_UNIT, FMT=103) current_batch, k_batch, & entropy, keff, keff_std @@ -417,20 +447,25 @@ contains write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if end if - else - ! Display output for inactive batch + end if + else + ! ======================================================================= + ! INACTIVE BATCHES + + ! Display output for inactive batch + if (master) then if (entropy_on) then write(UNIT=OUTPUT_UNIT, FMT=102) current_batch, k_batch, entropy else write(UNIT=OUTPUT_UNIT, FMT=100) current_batch, k_batch end if - - ! Set keff - keff = k_batch - - ! Reset tally values - global_tallies(:) % value = ZERO end if + + ! Set keff + keff = k_batch + + ! Reset tally values + global_tallies(:) % value = ZERO end if #ifdef MPI diff --git a/src/tally.F90 b/src/tally.F90 index a90289a0d..dba6f741a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1307,38 +1307,46 @@ contains subroutine synchronize_tallies() integer :: i ! index in tallies array - type(TallyObject), pointer :: t => null() #ifdef MPI - if (no_reduce) then - if (current_batch == n_batches) call reduce_tallies() - else - call reduce_tallies() + if (.not. no_reduce) then + call reduce_tally_values() if (.not. master) return end if #endif + ! Accumulate scores for each tally do i = 1, n_tallies - t => tallies(i) - - ! Loop over all filter and scoring bins - call accumulate_score(t % scores) + call accumulate_score(tallies(i) % scores) end do + ! Before accumulating scores for global_tallies, we need to get the current + ! batch estimate of k_analog for displaying to output + k_batch = global_tallies(K_ANALOG) % value + + ! Accumulate scores for global tallies + call accumulate_score(global_tallies) + +#ifdef MPI + if (no_reduce .and. current_batch == n_batches) & + call reduce_tally_sums() +#endif + end subroutine synchronize_tallies !=============================================================================== -! REDUCE_TALLIES collects all the results from tallies onto one processor +! REDUCE_TALLY_VALUES collects all the results from tallies onto one processor !=============================================================================== #ifdef MPI - subroutine reduce_tallies() + subroutine reduce_tally_values() integer :: i ! loop index for tallies integer :: n ! number of filter bins integer :: m ! number of score bins integer :: n_bins ! total number of bins real(8), allocatable :: tally_temp(:,:) ! contiguous array of scores + real(8) :: global_temp(N_GLOBAL_TALLIES) type(TallyObject), pointer :: t => null() do i = 1, n_tallies @@ -1353,26 +1361,42 @@ contains tally_temp = t % scores(:,:) % value if (master) then - ! The MPI_IN_PLACE specifier allows the master to copy values into a - ! receive buffer without having a temporary variable - call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + ! The MPI_IN_PLACE specifier allows the master to copy values into + ! a receive buffer without having a temporary variable + call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, & + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Transfer values to value on master t % scores(:,:) % value = tally_temp else ! Receive buffer not significant at other processors - call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & - 0, MPI_COMM_WORLD, mpi_err) + call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, & + MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) ! Reset value on other processors t % scores(:,:) % value = 0 end if deallocate(tally_temp) - end do + ! Copy global tallies into array to be reduced + global_temp = global_tallies(:) % value + + if (master) then + call MPI_REDUCE(MPI_IN_PLACE, global_temp, N_GLOBAL_TALLIES, & + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + + ! Transfer values back to global_tallies on master + global_tallies(:) % value = global_temp + else + call MPI_REDUCE(global_temp, global_temp, N_GLOBAL_TALLIES, & + MPI_REAL8, MPI_SUM, 0, MPI_COMM_WORLD, mpi_err) + + ! Reset value on other processors + global_tallies(:) % value = ZERO + end if + ! We also need to determine the total starting weight of particles from the ! last realization if (.not. no_reduce) then @@ -1385,8 +1409,79 @@ contains end if end if + end subroutine reduce_tally_values +#endif - end subroutine reduce_tallies +!=============================================================================== +! REDUCE_TALLY_SUMS gathers data from the sum and sum_sq member variables of +! TallyScores rather than the data from values. This is used if no_reduce is +! turned on and scores are accumulated on every processor *before* being +! reduced. +!=============================================================================== + +#ifdef MPI + subroutine reduce_tally_sums() + + integer :: i ! loop index for tallies + integer :: n ! number of filter bins + integer :: m ! number of score bins + integer :: n_bins ! total number of bins + real(8), allocatable :: tally_temp(:,:,:) ! contiguous array of scores + real(8) :: global_temp(2,N_GLOBAL_TALLIES) + type(TallyObject), pointer :: t => null() + + do i = 1, n_tallies + t => tallies(i) + + n = t % n_total_bins + m = t % n_score_bins + n_bins = n*m*2 + + allocate(tally_temp(2,n,m)) + + tally_temp(1,:,:) = t % scores(:,:) % sum + tally_temp(2,:,:) = t % scores(:,:) % sum_sq + + if (master) then + ! The MPI_IN_PLACE specifier allows the master to copy values into a + ! receive buffer without having a temporary variable + call MPI_REDUCE(MPI_IN_PLACE, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + + ! Transfer values to value on master + t % scores(:,:) % sum = tally_temp(1,:,:) + t % scores(:,:) % sum_sq = tally_temp(2,:,:) + else + ! Receive buffer not significant at other processors + call MPI_REDUCE(tally_temp, tally_temp, n_bins, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + end if + + deallocate(tally_temp) + + end do + + n_bins = 2 * N_GLOBAL_TALLIES + + global_temp(1,:) = global_tallies(:) % sum + global_temp(2,:) = global_tallies(:) % sum_sq + + if (master) then + ! The MPI_IN_PLACE specifier allows the master to copy values into a + ! receive buffer without having a temporary variable + call MPI_REDUCE(MPI_IN_PLACE, global_temp, n_bins, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + + ! Transfer values to value on master + global_tallies(:) % sum = global_temp(1,:) + global_tallies(:) % sum_sq = global_temp(2,:) + else + ! Receive buffer not significant at other processors + call MPI_REDUCE(global_temp, global_temp, n_bins, MPI_REAL8, MPI_SUM, & + 0, MPI_COMM_WORLD, mpi_err) + end if + + end subroutine reduce_tally_sums #endif !=============================================================================== From f4f2a85b6f086c46648224637dc0a10014e77393 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Apr 2012 11:12:22 -0400 Subject: [PATCH 3/8] Added 'modes' to improve program logic. --- src/constants.F90 | 6 ++++-- src/finalize.F90 | 4 ++-- src/geometry.F90 | 6 +++--- src/global.F90 | 11 +++-------- src/hdf5_interface.F90 | 2 +- src/initialize.F90 | 8 ++++---- src/input_xml.F90 | 2 +- src/main.F90 | 19 ++++++++++--------- src/output.F90 | 7 ++++--- 9 files changed, 32 insertions(+), 33 deletions(-) diff --git a/src/constants.F90 b/src/constants.F90 index 90828d301..9bd5df3dc 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -325,9 +325,11 @@ module constants SRC_POINT = 2, & ! Source at a single point SRC_FILE = 3 ! Source from a file + ! Running modes integer, parameter :: & - PROB_SOURCE = 1, & ! External source problem - PROB_CRITICALITY = 2 ! Criticality problem + MODE_FIXEDSOURCE = 1, & ! Fixed source mode + MODE_CRITICALITY = 2, & ! Criticality mode + MODE_PLOTTING = 3 ! Plotting mode ! Unit numbers integer, parameter :: UNIT_SUMMARY = 11 ! unit # for writing summary file diff --git a/src/finalize.F90 b/src/finalize.F90 index e239c6db9..a29a46706 100644 --- a/src/finalize.F90 +++ b/src/finalize.F90 @@ -24,7 +24,7 @@ contains ! Start finalization timer call timer_start(time_finalize) - if (.not. plotting) then + if (run_mode /= MODE_PLOTTING) then ! Calculate statistics for tallies and write to tallies.out call tally_statistics() if (master) call write_tallies() @@ -36,7 +36,7 @@ contains ! stop timers and show timing statistics call timer_stop(time_finalize) call timer_stop(time_total) - if (master .and. (.not. plotting)) call print_runtime() + if (master .and. (run_mode /= MODE_PLOTTING)) call print_runtime() #ifdef HDF5 ! Write time statistics to HDF5 output diff --git a/src/geometry.F90 b/src/geometry.F90 index 0abda64c1..ee67ae772 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -258,7 +258,7 @@ contains call write_message() end if - if (surf % bc == BC_VACUUM .and. (.not. plotting)) then + if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then ! ======================================================================= ! PARTICLE LEAKS OUT OF PROBLEM @@ -289,7 +289,7 @@ contains end if return - elseif (surf % bc == BC_REFLECT .and. (.not. plotting)) then + elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then ! ======================================================================= ! PARTICLE REFLECTS FROM SURFACE @@ -447,7 +447,7 @@ contains call find_cell(found) ! Couldn't find next cell anywhere! - if ((.not. found) .and. (.not. plotting)) then + if ((.not. found) .and. (run_mode /= MODE_PLOTTING)) then message = "After particle " // trim(to_str(p % id)) // " crossed surface " & // trim(to_str(surfaces(abs(p%surface)) % id)) // " it could not be & &located in any cell and it did not leak." diff --git a/src/global.F90 b/src/global.F90 index 7853ec917..a4d26264b 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -205,11 +205,6 @@ module global real(8) :: weight_cutoff = 0.25 real(8) :: weight_survive = 1.0 - ! ============================================================================ - ! PLOTTING VARIABLES - - logical :: plotting = .false. - ! ============================================================================ ! HDF5 VARIABLES @@ -221,6 +216,9 @@ module global ! ============================================================================ ! MISCELLANEOUS VARIABLES + ! Mode to run in (fixed source, criticality, plotting, etc) + integer :: run_mode = MODE_CRITICALITY + character(MAX_FILE_LEN) :: path_input ! Path to input file character(MAX_FILE_LEN) :: path_cross_sections ! Path to cross_sections.xml @@ -230,9 +228,6 @@ module global ! Random number seed integer(8) :: seed = 1_8 - ! Problem type - integer :: problem_type = PROB_CRITICALITY - ! The verbosity controls how much information will be printed to the ! screen and in logs integer :: verbosity = 7 diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index cfc738329..af85517a8 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -71,7 +71,7 @@ contains subroutine hdf5_write_summary() ! Write criticality information - if (problem_type == PROB_CRITICALITY) then + if (run_mode == MODE_CRITICALITY) then ! Need to write integer(8)'s using double instead since there is no H5LT ! call for making a dataset of type long call hdf5_make_double(hdf5_output_file, "n_particles", real(n_particles,8)) diff --git a/src/initialize.F90 b/src/initialize.F90 index c018efa3f..22770b6a3 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -82,7 +82,7 @@ contains ! Read plots.xml if it exists -- this has to be done separate from the other ! XML files because we need the PRNG to be initialized first - if (plotting) call read_plots_xml() + if (run_mode == MODE_PLOTTING) call read_plots_xml() ! Set up universe structures call prepare_universes() @@ -94,7 +94,7 @@ contains ! neighboring cells for efficient tracking call neighbor_lists() - if (.not. plotting) then + if (run_mode /= MODE_PLOTTING) then ! Read cross section summary file to determine what files contain ! cross-sections call read_cross_sections_xml() @@ -124,7 +124,7 @@ contains ! stop timer for initialization if (master) then - if (plotting) then + if (run_mode == MODE_PLOTTING) then call print_geometry() call print_plot() else @@ -243,7 +243,7 @@ contains if (starts_with(argv(i), "-")) then select case (argv(i)) case ('-p', '-plot', '--plot') - plotting = .true. + run_mode = MODE_PLOTTING case ('-?', '-help', '--help') call print_usage() stop diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 91bc79efc..c2262f4ba 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -95,7 +95,7 @@ contains ! Criticality information if (criticality % batches > 0) then - problem_type = PROB_CRITICALITY + if (run_mode /= MODE_PLOTTING) run_mode = MODE_CRITICALITY ! Check number of particles if (len_trim(criticality % particles) == 0) then diff --git a/src/main.F90 b/src/main.F90 index d5906aaf0..dba439b0d 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -23,12 +23,13 @@ program main ! set up problem call initialize_run() - ! start problem - if (plotting) then + ! start problem based on mode + select case (run_mode) + case (MODE_CRITICALITY) + call run_criticality() + case (MODE_PLOTTING) call run_plot() - else - call run_problem() - end if + end select ! finalize run call finalize_run() @@ -36,11 +37,11 @@ program main contains !=============================================================================== -! RUN_PROBLEM encompasses all the main logic where iterations are performed over -! the cycles and histories. +! RUN_CRITICALITY encompasses all the main logic where iterations are performed +! over the cycles and histories. !=============================================================================== - subroutine run_problem() + subroutine run_criticality() integer(8) :: i ! index over histories in single cycle @@ -147,6 +148,6 @@ contains if (master) call header("SIMULATION FINISHED", level=1) - end subroutine run_problem + end subroutine run_criticality end program main diff --git a/src/output.F90 b/src/output.F90 index d486be1c7..d95a735db 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -873,14 +873,15 @@ contains ! Display problem summary call header("PROBLEM SUMMARY", unit=UNIT_SUMMARY) - if (problem_type == PROB_CRITICALITY) then + select case(run_mode) + case (MODE_CRITICALITY) write(UNIT_SUMMARY,100) 'Problem type:', 'Criticality' 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 + case (MODE_FIXEDSOURCE) write(UNIT_SUMMARY,100) 'Problem type:', 'External Source' - end if + end select write(UNIT_SUMMARY,101) 'Number of Particles:', n_particles ! Display geometry summary From 3ffbd80efd2738078737848825f03e07ec216a6a Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Apr 2012 11:27:05 -0400 Subject: [PATCH 4/8] Added criticality module. --- src/DEPENDENCIES | 19 +++--- src/OBJECTS | 1 + src/criticality.F90 | 130 +++++++++++++++++++++++++++++++++++++++++ src/main.F90 | 139 ++------------------------------------------ 4 files changed, 149 insertions(+), 140 deletions(-) create mode 100644 src/criticality.F90 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 7d70c783e..08f362792 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -22,6 +22,16 @@ cross_section.o: material_header.o cross_section.o: random_lcg.o cross_section.o: search.o +criticality.o: constants.o +criticality.o: global.o +criticality.o: intercycle.o +criticality.o: output.o +criticality.o: physics.o +criticality.o: source.o +criticality.o: string.o +criticality.o: tally.o +criticality.o: timing.o + datatypes.o: datatypes_header.o doppler.o: constants.o @@ -138,16 +148,11 @@ interpolation.o: search.o interpolation.o: string.o main.o: constants.o +main.o: criticality.o +main.o: finalize.o main.o: global.o main.o: initialize.o -main.o: intercycle.o -main.o: output.o -main.o: physics.o main.o: plot.o -main.o: source.o -main.o: string.o -main.o: tally.o -main.o: timing.o mesh.o: constants.o mesh.o: global.o diff --git a/src/OBJECTS b/src/OBJECTS index e78fcd0c6..69229e181 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -3,6 +3,7 @@ ace.o \ ace_header.o \ bank_header.o \ cross_section.o \ +criticality.o \ datatypes.o \ datatypes_header.o \ doppler.o \ diff --git a/src/criticality.F90 b/src/criticality.F90 new file mode 100644 index 000000000..92f3c95b0 --- /dev/null +++ b/src/criticality.F90 @@ -0,0 +1,130 @@ +module criticality + + use constants, only: ZERO + use global + use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank, & + count_source_for_ufs + use output, only: write_message, header + use physics, only: transport + use source, only: get_source_particle + use string, only: to_str + use tally, only: synchronize_tallies + use timing, only: timer_start, timer_stop + +contains + +!=============================================================================== +! RUN_CRITICALITY encompasses all the main logic where iterations are performed +! over the cycles and histories. +!=============================================================================== + + subroutine run_criticality() + + integer(8) :: i ! index over histories in single cycle + + if (master) call header("BEGIN SIMULATION", level=1) + + tallies_on = .false. + call timer_start(time_inactive) + + ! Allocate particle + allocate(p) + + ! Display column titles + if (entropy_on) then + message = " Batch k(batch) Entropy Average k" + call write_message(1) + message = " ===== ======== ======= ===================" + call write_message(1) + else + message = " Batch k(batch) Average k" + call write_message(1) + message = " ===== ======== ===================" + call write_message(1) + end if + + ! ========================================================================== + ! LOOP OVER BATCHES + BATCH_LOOP: do current_batch = 1, n_batches + + message = "Simulating batch " // trim(to_str(current_batch)) // "..." + call write_message(8) + + ! Reset total starting weight + total_weight = ZERO + + ! ======================================================================= + ! LOOP OVER GENERATIONS + GENERATION_LOOP: do current_gen = 1, gen_per_batch + + ! Set all tallies to zero + n_bank = 0 + + ! Count source sites if using uniform fission source weighting + if (ufs) call count_source_for_ufs() + + ! ==================================================================== + ! LOOP OVER HISTORIES + + ! Start timer for transport + call timer_start(time_transport) + + HISTORY_LOOP: do i = 1, work + + ! grab source particle from bank + call get_source_particle(i) + + ! transport particle + call transport() + + end do HISTORY_LOOP + + ! 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() + + ! Stop timer for inter-cycle synchronization + call timer_stop(time_intercycle) + + end do GENERATION_LOOP + + ! Collect tallies + if (tallies_on) then + call timer_start(time_ic_tallies) + call synchronize_tallies() + call timer_stop(time_ic_tallies) + end if + + ! Calculate shannon entropy + if (entropy_on) call shannon_entropy() + + ! Collect results and statistics + call calculate_keff() + + ! Turn tallies on once inactive cycles are complete + if (current_batch == n_inactive) then + tallies_on = .true. + call timer_stop(time_inactive) + call timer_start(time_active) + end if + + end do BATCH_LOOP + + call timer_stop(time_active) + + ! ========================================================================== + ! END OF RUN WRAPUP + + if (master) call header("SIMULATION FINISHED", level=1) + + end subroutine run_criticality + +end module criticality diff --git a/src/main.F90 b/src/main.F90 index dba439b0d..5de70d8f5 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1,22 +1,11 @@ program main - use constants - use global - use finalize, only: finalize_run - use initialize, only: initialize_run - use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank, & - count_source_for_ufs - use output, only: write_message, header - use plotter, only: run_plot - use physics, only: transport - use source, only: get_source_particle - use string, only: to_str - use tally, only: synchronize_tallies - use timing, only: timer_start, timer_stop - -#ifdef MPI - use mpi -#endif + use constants, only: MODE_CRITICALITY, MODE_PLOTTING + use criticality, only: run_criticality + use finalize, only: finalize_run + use global, only: run_mode + use initialize, only: initialize_run + use plotter, only: run_plot implicit none @@ -34,120 +23,4 @@ program main ! finalize run call finalize_run() -contains - -!=============================================================================== -! RUN_CRITICALITY encompasses all the main logic where iterations are performed -! over the cycles and histories. -!=============================================================================== - - subroutine run_criticality() - - integer(8) :: i ! index over histories in single cycle - - if (master) call header("BEGIN SIMULATION", level=1) - - tallies_on = .false. - call timer_start(time_inactive) - - ! Allocate particle - allocate(p) - - ! Display column titles - if (entropy_on) then - message = " Batch k(batch) Entropy Average k" - call write_message(1) - message = " ===== ======== ======= ===================" - call write_message(1) - else - message = " Batch k(batch) Average k" - call write_message(1) - message = " ===== ======== ===================" - call write_message(1) - end if - - ! ========================================================================== - ! LOOP OVER BATCHES - BATCH_LOOP: do current_batch = 1, n_batches - - message = "Simulating batch " // trim(to_str(current_batch)) // "..." - call write_message(8) - - ! Reset total starting weight - total_weight = ZERO - - ! ======================================================================= - ! LOOP OVER GENERATIONS - GENERATION_LOOP: do current_gen = 1, gen_per_batch - - ! Set all tallies to zero - n_bank = 0 - - ! Count source sites if using uniform fission source weighting - if (ufs) call count_source_for_ufs() - - ! ==================================================================== - ! LOOP OVER HISTORIES - - ! Start timer for transport - call timer_start(time_transport) - - HISTORY_LOOP: do i = 1, work - - ! grab source particle from bank - call get_source_particle(i) - - ! transport particle - call transport() - - end do HISTORY_LOOP - - ! 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() - - ! Stop timer for inter-cycle synchronization - call timer_stop(time_intercycle) - - end do GENERATION_LOOP - - ! Collect tallies - if (tallies_on) then - call timer_start(time_ic_tallies) - call synchronize_tallies() - call timer_stop(time_ic_tallies) - end if - - ! Calculate shannon entropy - if (entropy_on) call shannon_entropy() - - ! Collect results and statistics - call calculate_keff() - - ! Turn tallies on once inactive cycles are complete - if (current_batch == n_inactive) then - tallies_on = .true. - call timer_stop(time_inactive) - call timer_start(time_active) - end if - - end do BATCH_LOOP - - call timer_stop(time_active) - - ! ========================================================================== - ! END OF RUN WRAPUP - - if (master) call header("SIMULATION FINISHED", level=1) - - end subroutine run_criticality - end program main From 66b06a3bda96bc66d0ec24f1463c3afb40c74195 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Apr 2012 11:40:14 -0400 Subject: [PATCH 5/8] Added initialize/finalize routines in criticality module. --- src/criticality.F90 | 108 +++++++++++++++++++++++++++----------------- 1 file changed, 66 insertions(+), 42 deletions(-) diff --git a/src/criticality.F90 b/src/criticality.F90 index 92f3c95b0..e2492cbf8 100644 --- a/src/criticality.F90 +++ b/src/criticality.F90 @@ -47,29 +47,20 @@ contains ! LOOP OVER BATCHES BATCH_LOOP: do current_batch = 1, n_batches - message = "Simulating batch " // trim(to_str(current_batch)) // "..." - call write_message(8) - - ! Reset total starting weight - total_weight = ZERO + call initialize_batch() ! ======================================================================= ! LOOP OVER GENERATIONS GENERATION_LOOP: do current_gen = 1, gen_per_batch - ! Set all tallies to zero - n_bank = 0 - - ! Count source sites if using uniform fission source weighting - if (ufs) call count_source_for_ufs() - - ! ==================================================================== - ! LOOP OVER HISTORIES + call initialize_generation() ! Start timer for transport call timer_start(time_transport) - HISTORY_LOOP: do i = 1, work + ! ==================================================================== + ! LOOP OVER PARTICLES + PARTICLE_LOOP: do i = 1, work ! grab source particle from bank call get_source_particle(i) @@ -77,44 +68,19 @@ contains ! transport particle call transport() - end do HISTORY_LOOP + end do PARTICLE_LOOP ! 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 timer_start(time_intercycle) call synchronize_bank() - - ! Stop timer for inter-cycle synchronization call timer_stop(time_intercycle) end do GENERATION_LOOP - ! Collect tallies - if (tallies_on) then - call timer_start(time_ic_tallies) - call synchronize_tallies() - call timer_stop(time_ic_tallies) - end if - - ! Calculate shannon entropy - if (entropy_on) call shannon_entropy() - - ! Collect results and statistics - call calculate_keff() - - ! Turn tallies on once inactive cycles are complete - if (current_batch == n_inactive) then - tallies_on = .true. - call timer_stop(time_inactive) - call timer_start(time_active) - end if + call finalize_batch() end do BATCH_LOOP @@ -127,4 +93,62 @@ contains end subroutine run_criticality +!=============================================================================== +! INITIALIZE_BATCH +!=============================================================================== + + subroutine initialize_batch() + + message = "Simulating batch " // trim(to_str(current_batch)) // "..." + call write_message(8) + + ! Reset total starting particle weight used for normalizing tallies + total_weight = ZERO + + end subroutine initialize_batch + +!=============================================================================== +! INITIALIZE_GENERATION +!=============================================================================== + + subroutine initialize_generation() + + ! 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() + + end subroutine initialize_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 + if (tallies_on) then + call timer_start(time_ic_tallies) + call synchronize_tallies() + call timer_stop(time_ic_tallies) + end if + + ! Calculate shannon entropy + if (entropy_on) call shannon_entropy() + + ! Collect results and statistics + call calculate_keff() + + ! Turn tallies on once inactive cycles are complete + if (current_batch == n_inactive) then + tallies_on = .true. + call timer_stop(time_inactive) + call timer_start(time_active) + end if + + end subroutine finalize_batch + end module criticality From c74a7eacecbdd67f99f931f00d9f9fc8b5937003 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Sat, 7 Apr 2012 13:05:18 -0400 Subject: [PATCH 6/8] Updated documentation. --- docs/source/releasenotes/notes_0.4.2.rst | 12 ++++++---- docs/source/usersguide/input.rst | 30 +++++++++++++++++++++--- docs/source/usersguide/troubleshoot.rst | 2 +- 3 files changed, 36 insertions(+), 8 deletions(-) diff --git a/docs/source/releasenotes/notes_0.4.2.rst b/docs/source/releasenotes/notes_0.4.2.rst index 825a37d3a..f213581c4 100644 --- a/docs/source/releasenotes/notes_0.4.2.rst +++ b/docs/source/releasenotes/notes_0.4.2.rst @@ -22,6 +22,8 @@ hand (mostly on the number of nuclides in the problem). New Features ------------ +- Option to not reduce tallies across processors at end of each batch. +- Uniform fission site method for reducing variance on local tallies. - Reading/writing binary source files. - Added more messages for or high verbosity. - Estimator for diffusion coefficient. @@ -35,11 +37,13 @@ New Features Bug Fixes --------- -- `b2c40e`_: Fixed bug in incoming energy filter for track-length tallies. -- `5524fd`_: Mesh filter now works with track-length tallies. -- `d050c7`_: Added Bessel's correction to make estimate of variance unbiased. -- `2a5b9c`_: Fixed regression in plotting. +- 671f30_: Fixed surface currents on mesh not encompassing geometry. +- b2c40e_: Fixed bug in incoming energy filter for track-length tallies. +- 5524fd_: Mesh filter now works with track-length tallies. +- d050c7_: Added Bessel's correction to make estimate of variance unbiased. +- 2a5b9c_: Fixed regression in plotting. +.. _671f30: https://github.com/mit-crpg/openmc/commit/671f30 .. _b2c40e: https://github.com/mit-crpg/openmc/commit/b2c40e .. _5524fd: https://github.com/mit-crpg/openmc/commit/5524fd .. _d050c7: https://github.com/mit-crpg/openmc/commit/d050c7 diff --git a/docs/source/usersguide/input.rst b/docs/source/usersguide/input.rst index f6104bba4..4aab4684a 100644 --- a/docs/source/usersguide/input.rst +++ b/docs/source/usersguide/input.rst @@ -151,6 +151,19 @@ problem. It has the following attributes/sub-elements: *Default*: None +```` Element +----------------------- + +The ```` element has no attributes and has an accepted value of "on" +or "off". If set to "on", all user-defined tallies and global tallies will not +be reduced across processors in a parallel calculation. This means that the +accumulate score in one batch on a single processor is considered as an +independent realization for the tally random variable. For a problem with large +tally data, this option can significantly improve the parallel efficiency. + + *Default*: off + + ```` Element --------------------- @@ -195,12 +208,14 @@ criticality calculations. It takes the following attributes: ```` Element ------------------------------ -The ```` element has no attributes and assumes wither the -value ``on`` or ``off``. If turned on, this option will enable the use of -survival biasing, otherwise known as implicit capture or absorption. +The ```` element has no attributes and has an accepted value +of "on" or "off". If set to "on", this option will enable the use of survival +biasing, otherwise known as implicit capture or absorption. *Default*: off +.. _trace: + ```` Element ------------------- @@ -247,6 +262,15 @@ displayed. This element takes the following attributes: *Default*: 5 +```` Element +------------------------------ + +The ```` element has no attributes and has an accepted value of +"on" or "off". If set to "on", a binary source file will be written to diskat +the end of the run that can be used as a starting source for another run. + + *Default*: off + -------------------------------------- Geometry Specification -- geometry.xml -------------------------------------- diff --git a/docs/source/usersguide/troubleshoot.rst b/docs/source/usersguide/troubleshoot.rst index 40aa6a208..5943d553f 100644 --- a/docs/source/usersguide/troubleshoot.rst +++ b/docs/source/usersguide/troubleshoot.rst @@ -94,7 +94,7 @@ geometry has been specified incorrectly or is missing. The best way to debug this error is to turn on a trace for the particle getting lost. After the error message, the code will display what batch, generation, and -particle number caused the error. In your settings.xml, add a tag +particle number caused the error. In your settings.xml, add a :ref:`trace` followed by the batch, generation, and particle number. This will give you detailed output every time that particle enters a cell, crosses a boundary, or has a collision. For example, if you received this error at cycle 5, generation From 33ebeedd78096205beca1dacdac472fde19f605d Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 9 Apr 2012 10:17:46 -0400 Subject: [PATCH 7/8] Added check for number of active batches > 0. --- src/input_xml.F90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/input_xml.F90 b/src/input_xml.F90 index c2262f4ba..d1fa801fe 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -109,6 +109,12 @@ contains n_inactive = criticality % inactive n_active = n_batches - n_inactive gen_per_batch = criticality % generations_per_batch + + ! Check number of active batches + if (n_active <= 0) then + message = "Number of active batches must be greater than 0." + call fatal_error() + end if else message = "Need to specify number of batches with tag." call fatal_error() From 298db8f2b5e30d1a11d90459c2e66e3cf58f0d42 Mon Sep 17 00:00:00 2001 From: Bryan Herman Date: Tue, 10 Apr 2012 10:13:02 -0700 Subject: [PATCH 8/8] fixed bug in score surface current when using energy in filter --- src/tally.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/tally.F90 b/src/tally.F90 index dba6f741a..0acaaf254 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -1032,11 +1032,11 @@ contains n = t % n_filter_bins(FILTER_ENERGYIN) if (n > 0) then ! check if energy of the particle is within energy bins - if (p % last_E < t % energy_in(1) .or. & - p % last_E > t % energy_in(n + 1)) cycle + if (p % E < t % energy_in(1) .or. & + p % E > t % energy_in(n + 1)) cycle ! search to find incoming energy bin - bins(SURF_FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % last_E) + bins(SURF_FILTER_ENERGYIN) = binary_search(t % energy_in, n + 1, p % E) else bins(SURF_FILTER_ENERGYIN) = 1 end if