crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-27 05:35:49 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

merged in master

Browse source
This commit is contained in:
Bryan Herman 2012-04-10 11:36:24 -07:00
commit eedd1184af
19 changed files with 504 additions and 262 deletions
Download patch file Download diff file Expand all files Collapse all files

12
docs/source/releasenotes/notes_0.4.2.rst
View file

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

30
docs/source/usersguide/input.rst
View file

@ -151,6 +151,19 @@ problem. It has the following attributes/sub-elements:
*Default*: None
``<no_reduce>`` Element
-----------------------
The ``<no_reduce>`` 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
``<ptables>`` Element
---------------------
@ -195,12 +208,14 @@ criticality calculations. It takes the following attributes:
``<survival_biasing>`` Element
------------------------------
The ``<survival_biasing>`` 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 ``<survival_biasing>`` 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:
``<trace>`` Element
-------------------
@ -247,6 +262,15 @@ displayed. This element takes the following attributes:
*Default*: 5
``<write_source>`` Element
------------------------------
The ``<write_source>`` 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
--------------------------------------

2
docs/source/usersguide/troubleshoot.rst
View file

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

19
src/DEPENDENCIES
View file

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

1
src/OBJECTS
View file

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

6
src/constants.F90
View file

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

175
src/criticality.F90 Normal file
View file

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

4
src/finalize.F90
View file

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

6
src/geometry.F90
View file

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

29
src/global.F90
View file

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

2
src/hdf5_interface.F90
View file

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

8
src/initialize.F90
View file

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

21
src/input_xml.F90
View file

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

113
src/intercycle.F90
View file

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

169
src/main.F90
View file

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

7
src/output.F90
View file

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

3
src/source.F90
View file

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

158
src/tally.F90
View file

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

1
src/xml-fortran/templates/settings_t.xml
View file

@ -32,6 +32,7 @@
<variable name="cutoff_" tag="cutoff" type="cutoff_xml" dimension="1" />
<variable name="energy_grid_" tag="energy_grid" type="word" length="7" />
<variable name="entropy_" tag="entropy" type="mesh_xml" dimension="1" />
<variable name="no_reduce_" tag="no_reduce" type="word" length="3" />
<variable name="ptables_" tag="ptables" type="word" length="3" />
<variable name="seed_" tag="seed" type="integer" />
<variable name="source_" tag="source" type="source_xml" />