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 diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index a00e8a331..eede92f66 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -68,6 +68,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 @@ -190,16 +200,11 @@ logging.o: global.o main.o: cmfd_execute.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 734b6da54..c6c4fb47d 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -14,6 +14,7 @@ cmfd_prod_operator.o \ cmfd_slepc_solver.o \ cmfd_snes_solver.o \ cross_section.o \ +criticality.o \ datatypes.o \ datatypes_header.o \ doppler.o \ 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/criticality.F90 b/src/criticality.F90 new file mode 100644 index 000000000..cd470f8b2 --- /dev/null +++ b/src/criticality.F90 @@ -0,0 +1,175 @@ +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 + if (cmfd_on) then + message = " Cycle k(batch) Entropy Average k CMFD k CMFD Ent." + call write_message(1) + message = " ===== ======== ======= =================== ========== =========" + call write_message(1) + else + message = " Cycle k(batch) Entropy Average k" + call write_message(1) + message = " ===== ======== ======= ===================" + call write_message(1) + end if + else + if (cmfd_on) then + message = " Cycle k(batch) Average k CMFD k" + call write_message(1) + message = " ===== ======== =================== ==========" + call write_message(1) + else + message = " Cycle k(batch) Average k" + call write_message(1) + message = " ===== ======== ===================" + call write_message(1) + end if + end if + + ! ========================================================================== + ! LOOP OVER BATCHES + BATCH_LOOP: do current_batch = 1, n_batches + + call initialize_batch() + + ! ======================================================================= + ! LOOP OVER GENERATIONS + GENERATION_LOOP: do current_gen = 1, gen_per_batch + + call initialize_generation() + + ! Start timer for transport + call timer_start(time_transport) + + ! ==================================================================== + ! LOOP OVER PARTICLES + PARTICLE_LOOP: do i = 1, work + + ! grab source particle from bank + call get_source_particle(i) + + ! transport particle + call transport() + + end do PARTICLE_LOOP + + ! Accumulate time for transport + call timer_stop(time_transport) + + ! Distribute fission bank across processors evenly + call timer_start(time_intercycle) + call synchronize_bank() + call timer_stop(time_intercycle) + + end do GENERATION_LOOP + + call finalize_batch() + + 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 + +!=============================================================================== +! 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() + + use cmfd_execute, only: execute_cmfd + + ! 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() + + ! run cmfd + if (current_batch > n_inactive .and. cmfd_on) then + call execute_cmfd() + end if + + ! 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 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 0c0ab3e89..dda59c2a5 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -118,6 +118,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 @@ -135,9 +139,6 @@ module global integer :: current_batch = 0 ! current batch integer :: current_gen = 0 ! current generation - ! Random Number seed - integer(8) :: seed = 1_8 - ! External source type(ExtSource), target :: external_source @@ -150,9 +151,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. @@ -178,6 +180,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 @@ -201,11 +206,6 @@ module global real(8) :: weight_cutoff = 0.25 real(8) :: weight_survive = 1.0 - ! ============================================================================ - ! PLOTTING VARIABLES - - logical :: plotting = .false. - ! ============================================================================ ! HDF5 VARIABLES @@ -217,14 +217,17 @@ 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 ! Message used in message/warning/fatal_error character(MAX_LINE_LEN) :: message - ! Problem type - integer :: problem_type = PROB_CRITICALITY + ! Random number seed + integer(8) :: seed = 1_8 ! The verbosity controls how much information will be printed to the ! screen and in logs diff --git a/src/hdf5_interface.F90 b/src/hdf5_interface.F90 index 56e921822..dcf403aea 100644 --- a/src/hdf5_interface.F90 +++ b/src/hdf5_interface.F90 @@ -95,7 +95,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 06da79afb..4e30d8e52 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 58bae1f36..e1db03db0 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -70,6 +70,7 @@ contains energy_grid_ = "union" seed_ = 0_8 write_source_ = "" + no_reduce_ = "" source_ % type = "" ! Parse settings.xml file @@ -96,7 +97,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 @@ -110,6 +111,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() @@ -295,6 +302,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/intercycle.F90 b/src/intercycle.F90 index 98745023d..25ec88abe 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 if (cmfd_on) then write(UNIT=OUTPUT_UNIT, FMT=104) current_batch, k_batch, & @@ -427,20 +457,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/main.F90 b/src/main.F90 index b88f0cf31..e0b145f32 100644 --- a/src/main.F90 +++ b/src/main.F90 @@ -1,169 +1,26 @@ program main - use constants - use cmfd_execute, only: execute_cmfd - 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 ! 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() - -contains -!=============================================================================== -! RUN_PROBLEM encompasses all the main logic where iterations are performed over -! the cycles and histories. -!=============================================================================== - - subroutine run_problem() - - 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 - if (cmfd_on) then - message = " Cycle k(batch) Entropy Average k CMFD k CMFD Ent." - call write_message(1) - message = " ===== ======== ======= =================== ========== =========" - call write_message(1) - else - message = " Cycle k(batch) Entropy Average k" - call write_message(1) - message = " ===== ======== ======= ===================" - call write_message(1) - end if - else - if (cmfd_on) then - message = " Cycle k(batch) Average k CMFD k" - call write_message(1) - message = " ===== ======== =================== ==========" - call write_message(1) - else - message = " Cycle k(batch) Average k" - call write_message(1) - message = " ===== ======== ===================" - call write_message(1) - end if - 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) - - ! ======================================================================= - ! 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() - - ! run cmfd - if (current_batch > n_inactive .and. cmfd_on) then - call execute_cmfd() - end if - - ! 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_problem - -end program main +end program main diff --git a/src/output.F90 b/src/output.F90 index f677fa8c0..09fd64537 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 diff --git a/src/source.F90 b/src/source.F90 index c48744b97..684656e8f 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 7c15f38cc..617a1d6e8 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 @@ -1307,34 +1307,46 @@ contains subroutine synchronize_tallies() integer :: i ! index in tallies array - type(TallyObject), pointer :: t => null() #ifdef MPI - call reduce_tallies() - if (.not. master) return + 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 @@ -1349,27 +1361,127 @@ 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 + 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_tally_values +#endif + +!=============================================================================== +! 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 - end subroutine reduce_tallies + 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 !=============================================================================== @@ -1808,7 +1920,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 +1942,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 eeb26cbbf..339f393df 100644 --- a/src/xml-fortran/templates/settings_t.xml +++ b/src/xml-fortran/templates/settings_t.xml @@ -32,6 +32,7 @@ +