crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 22:26:08 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

merged master

Browse source
This commit is contained in:
Bryan Herman 2012-04-01 20:31:38 -07:00
commit 8ddf58a0bc
19 changed files with 445 additions and 133 deletions
Download patch file Download diff file Expand all files Collapse all files

1
docs/source/index.rst
View file

@ -29,4 +29,5 @@ Contents
usersguide/index
devguide/index
publications
license
developers

24
docs/source/license.rst Normal file
View file

@ -0,0 +1,24 @@
.. _license:
=================
License Agreement
=================
Copyright © 2011-2012 Massachusetts Institute of Technology
Permission is hereby granted, free of charge, to any person obtaining a copy of
this software and associated documentation files (the "Software"), to deal in
the Software without restriction, including without limitation the rights to
use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
the Software, and to permit persons to whom the Software is furnished to do so,
subject to the following conditions:
The above copyright notice and this permission notice shall be included in all
copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

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

@ -210,6 +210,31 @@ integers: the batch number, generation number, and particle number.
*Default*: None
``<uniform_fs>`` Element
------------------------
The ``<uniform_fs>`` element describes a mesh that is used for re-weighting
source sites at every generation based on the uniform fission site methodology
described in Kelly et al., "MC21 Analysis of the Nuclear Energy Agency Monte
Carlo Performance Benchmark Problem," Proceedings of *Physor 2012*, Knoxville,
TN (2012). This mesh should cover all possible fissionable materials in the
problem. It has the following attributes/sub-elements:
:dimension:
The number of mesh cells in the x, y, and z directions, respectively.
*Default*: None
:lower_left:
The Cartersian coordinates of the lower-left corner of the mesh.
*Default*: None
:upper_right:
The Cartersian coordinates of the upper-right corner of the mesh.
*Default*: None
``<verbosity>`` Element
-----------------------

62
docs/source/usersguide/setup.rst
View file

@ -41,9 +41,11 @@ been compiled with the same compiler you intend to compile OpenMC with.
Obtaining the Source
--------------------
All OpenMC source code is hosted on GitHub_. This means that you will need to
have git_ installed on your computer in order to get source code and updates
directly from the repository. GitHub has a good set of `instructions
All OpenMC source code is hosted on GitHub_. You can download the source code
directly from GitHub or, if you have the git_ version control software installed
on your computer, you can use git to obtain the source code. The latter method
has the benefit that it is easy to receive updates directly from the GitHub
repository. GitHub has a good set of `instructions
<http://help.github.com/set-up-git-redirect>`_ for how to set up git to work
with GitHub since this involves setting up ssh_ keys. With git installed and
setup, the following command will download the full source code from the GitHub
@ -96,12 +98,12 @@ Makefile, you can enter the following from a terminal::
make DEBUG=yes
---------
Compiling
---------
-------------------------------
Compiling on Linux and Mac OS X
-------------------------------
To compile the code, run the following commands from within the root directory
for OpenMC:
To compile OpenMC on Linux or Max OS X, run the following commands from within
the root directory of the source code:
.. code-block:: sh
@ -110,6 +112,25 @@ for OpenMC:
This will build an executable named ``openmc``.
--------------------
Compiling on Windows
--------------------
To compile OpenMC on a Windows operating system, you will need to first install
Cygwin_, a Linux-like environment for Windows. When configuring Cygwin, make
sure you install both the gfortran compiler as well as git. Once you have
obtained the source code, run the following commands from within the source code
root directory:
.. code-block:: sh
cd src
make
This will build an executable named ``openmc``.
.. _Cygwin: http://www.cygwin.com/
---------------------------
Cross-Section Configuration
---------------------------
@ -130,3 +151,28 @@ directory of the OpenMC distribution to the location of the Serpent
cross-sections. Then, either set the :ref:`cross_sections` in a settings.xml
file or the :envvar:`CROSS_SECTIONS` environment variable to the absolute path
of the ``cross_sections_serpent.xml`` file.
--------------
Running OpenMC
--------------
Once you have a model built (see :ref:`usersguide_input`), you can either run
the openmc executable directly from the directory containing your XML input
files, or you can specify as a command-line argument the directory containing
the XML input files. For example, if the path of your OpenMC exectuable is
``/home/username/openmc/src/openmc`` and your XML input files are in the
directory ``/home/username/somemodel/``, one way to run the simulation would be:
.. code-block:: sh
cd /home/username/somemodel
openmc
Alternatively, you could run from any directory:
.. code-block:: sh
openmc /home/username/somemodel
Note that in the latter case, any output files will be placed in the present
working directory which may be different from ``/home/username/somemodel``.

1
src/DEPENDENCIES
View file

@ -231,6 +231,7 @@ physics.o: geometry_header.o
physics.o: global.o
physics.o: interpolation.o
physics.o: material_header.o
physics.o: mesh.o
physics.o: output.o
physics.o: particle_header.o
physics.o: random_lcg.o

1
src/bank_header.F90
View file

@ -15,6 +15,7 @@ module bank_header
sequence
integer(8) :: id ! Unique ID
real(8) :: wgt ! weight of bank site
real(8) :: xyz(3) ! location of bank particle
real(8) :: uvw(3) ! diretional cosines
real(8) :: E ! energy

9
src/cross_section.F90
View file

@ -187,8 +187,13 @@ contains
end if
end if
! Set last evaluated energy
if (index_sab == 0) micro_xs(index_nuclide) % last_E = p % E
! Set last evaluated energy -- if we're in S(a,b) region, force
! re-calculation of cross-section
if (index_sab == 0) then
micro_xs(index_nuclide) % last_E = p % E
else
micro_xs(index_nuclide) % last_E = ZERO
end if
end subroutine calculate_nuclide_xs

5
src/global.F90
View file

@ -160,6 +160,11 @@ module global
real(8), allocatable :: entropy_p(:,:,:,:) ! % of source sites in each cell
type(StructuredMesh), pointer :: entropy_mesh
! Uniform fission source weighting
logical :: ufs = .false.
type(StructuredMesh), pointer :: ufs_mesh => null()
real(8), allocatable :: source_frac(:,:,:,:)
! Write source at end of simulation
logical :: write_source = .false.

39
src/initialize.F90
View file

@ -14,8 +14,8 @@ module initialize
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict, read_plots_xml
use output, only: title, header, print_summary, print_geometry, &
print_plot, create_summary_file, &
create_xs_summary_file
print_plot, create_summary_file, print_usage, &
create_xs_summary_file, print_version
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: to_str, starts_with, ends_with, lower_case
@ -149,9 +149,9 @@ contains
subroutine setup_mpi()
#ifdef MPI
integer :: bank_blocks(4) ! Count for each datatype
integer :: bank_types(4) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(4) ! Displacements
integer :: bank_blocks(5) ! Count for each datatype
integer :: bank_types(5) ! Datatypes
integer(MPI_ADDRESS_KIND) :: bank_disp(5) ! Displacements
type(Bank) :: b
mpi_enabled = .true.
@ -186,17 +186,18 @@ contains
! Determine displacements for MPI_BANK type
call MPI_GET_ADDRESS(b % id, bank_disp(1), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(4), mpi_err)
call MPI_GET_ADDRESS(b % wgt, bank_disp(2), mpi_err)
call MPI_GET_ADDRESS(b % xyz, bank_disp(3), mpi_err)
call MPI_GET_ADDRESS(b % uvw, bank_disp(4), mpi_err)
call MPI_GET_ADDRESS(b % E, bank_disp(5), mpi_err)
! Adjust displacements
bank_disp = bank_disp - bank_disp(1)
! Define MPI_BANK for fission sites
bank_blocks = (/ 1, 3, 3, 1 /)
bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(4, bank_blocks, bank_disp, &
bank_blocks = (/ 1, 1, 3, 3, 1 /)
bank_types = (/ MPI_INTEGER8, MPI_REAL8, MPI_REAL8, MPI_REAL8, MPI_REAL8 /)
call MPI_TYPE_CREATE_STRUCT(5, bank_blocks, bank_disp, &
bank_types, MPI_BANK, mpi_err)
call MPI_TYPE_COMMIT(MPI_BANK, mpi_err)
@ -235,9 +236,21 @@ contains
! Check for flags
if (starts_with(argv(i), "-")) then
if (argv(i)(2:5) == 'plot') then
select case (argv(i))
case ('-p', '-plot', '--plot')
plotting = .true.
end if
case ('-?', '-help', '--help')
call print_usage()
stop
case ('-v', '-version', '--version')
call print_version()
stop
case default
message = "Unknown command line option: " // argv(i)
call fatal_error()
end select
last_flag = i
end if
! Determine directory where XML input files are

57
src/input_xml.F90
View file

@ -238,6 +238,60 @@ contains
entropy_on = .true.
end if
! Uniform fission source weighting mesh
if (size(uniform_fs_) > 0) then
! Check to make sure enough values were supplied
if (size(uniform_fs_(1) % lower_left) /= 3) then
message = "Need to specify (x,y,z) coordinates of lower-left corner &
&of UFS mesh."
elseif (size(uniform_fs_(1) % upper_right) /= 3) then
message = "Need to specify (x,y,z) coordinates of upper-right corner &
&of UFS mesh."
elseif (size(uniform_fs_(1) % dimension) /= 3) then
message = "Dimension of UFS mesh must be given as three &
&integers."
call fatal_error()
end if
! Allocate mesh object and coordinates on mesh
allocate(ufs_mesh)
allocate(ufs_mesh % lower_left(3))
allocate(ufs_mesh % upper_right(3))
allocate(ufs_mesh % width(3))
! Allocate dimensions
ufs_mesh % n_dimension = 3
allocate(ufs_mesh % dimension(3))
! Copy dimensions
ufs_mesh % dimension = uniform_fs_(1) % dimension
! Copy values
ufs_mesh % lower_left = uniform_fs_(1) % lower_left
ufs_mesh % upper_right = uniform_fs_(1) % upper_right
! Check on values provided
if (.not. all(ufs_mesh % upper_right > ufs_mesh % lower_left)) then
message = "Upper-right coordinate must be greater than lower-left &
&coordinate for UFS mesh."
call fatal_error()
end if
! Calculate width
ufs_mesh % width = (ufs_mesh % upper_right - &
ufs_mesh % lower_left) / ufs_mesh % dimension
! Calculate volume fraction of each cell
ufs_mesh % volume_frac = ONE/real(product(ufs_mesh % dimension),8)
! Turn on uniform fission source weighting
ufs = .true.
! Allocate source_frac
allocate(source_frac(1, ufs_mesh % dimension(1), &
ufs_mesh % dimension(2), ufs_mesh % dimension(3)))
end if
! Check if the user has specified to write binary source file
if (trim(write_source_) == 'on') write_source = .true.
@ -814,6 +868,9 @@ contains
! Set upper right coordinate
m % upper_right = m % lower_left + m % dimension * m % width
! Set volume fraction
m % volume_frac = ONE/real(product(m % dimension),8)
! Add mesh to dictionary
call dict_add_key(mesh_dict, m % id, i)
end do

64
src/intercycle.F90
View file

@ -4,13 +4,13 @@ module intercycle
use error, only: fatal_error, warning
use global
use mesh, only: count_fission_sites
use mesh, only: count_bank_sites
use mesh_header, only: StructuredMesh
use output, only: write_message
use random_lcg, only: prn, set_particle_seed, prn_skip
use search, only: binary_search
use string, only: to_str
use tally, only: accumulate_batch_estimate
use tally, only: accumulate_score
use tally_header, only: TallyObject
use timing, only: timer_start, timer_stop
@ -285,8 +285,7 @@ contains
integer :: i, j, k ! index for bank sites
integer :: n ! # of boxes in each dimension
real(8) :: total ! total weight of fission bank sites
logical :: sites_outside ! were there sites outside entropy box?
logical :: sites_outside ! were there sites outside entropy box?
type(StructuredMesh), pointer :: m => null()
! Get pointer to entropy mesh
@ -320,7 +319,7 @@ contains
end if
! count number of fission sites over mesh
call count_fission_sites(m, fission_bank, entropy_p, total, &
call count_bank_sites(m, fission_bank, entropy_p, &
size_bank=n_bank, sites_outside=sites_outside)
! display warning message if there were sites outside entropy box
@ -332,7 +331,7 @@ contains
! sum values to obtain shannon entropy
if (master) then
! Normalize to total weight of bank sites
entropy_p = entropy_p / total
entropy_p = entropy_p / sum(entropy_p)
entropy = 0
do i = 1, m % dimension(1)
@ -393,7 +392,7 @@ contains
n = current_batch - n_inactive
! Accumulate single batch realizations of k
call accumulate_batch_estimate(global_tallies)
call accumulate_score(global_tallies)
! Determine sample mean of keff
keff = global_tallies(K_ANALOG) % sum/n
@ -458,4 +457,55 @@ contains
end subroutine calculate_keff
!===============================================================================
! COUNT_SOURCE_FOR_UFS determines the source fraction in each UFS mesh cell and
! reweights the source bank so that the sum of the weights is equal to
! n_particles. The 'source_frac' variable is used later to bias the production
! of fission sites
!===============================================================================
subroutine count_source_for_ufs()
real(8) :: total ! total weight in source bank
logical :: sites_outside ! were there sites outside the ufs mesh?
#ifdef MPI
integer :: n ! total number of ufs mesh cells
#endif
if (current_batch == 1 .and. current_gen == 1) then
! On the first cycle, just assume that the source is already evenly
! distributed so that effectively the production of fission sites is not
! biased
source_frac = ufs_mesh % volume_frac
else
! count number of source sites in each ufs mesh cell
call count_bank_sites(ufs_mesh, source_bank, source_frac, &
sites_outside=sites_outside)
! Check for sites outside of the mesh
if (master .and. sites_outside) then
message = "Source sites outside of the UFS mesh!"
call fatal_error()
end if
#ifdef MPI
! Send source fraction to all processors
n = product(ufs_mesh % dimension)
call MPI_BCAST(source_frac, n, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
#endif
! Normalize to total weight to get fraction of source in each cell
total = sum(source_frac)
source_frac = source_frac / total
! Since the total starting weight is not equal to n_particles, we need to
! renormalize the weight of the source sites
source_bank % wgt = source_bank % wgt * n_particles / total
end if
end subroutine count_source_for_ufs
end module intercycle

29
src/main.F90
View file

@ -3,16 +3,17 @@ 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
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, add_to_score
use timing, only: timer_start, timer_stop
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
@ -93,6 +94,9 @@ contains
! 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
@ -121,11 +125,6 @@ contains
! Distribute fission bank across processors evenly
call synchronize_bank()
! Add to analog estimate of keff -- since the creation of bank sites
! was originally weighted by the last cycle keff, we need to multiply
! by that keff to get the current cycle's value
call add_to_score(global_tallies(K_ANALOG), n_bank * keff)
! Stop timer for inter-cycle synchronization
call timer_stop(time_intercycle)

26
src/mesh.F90
View file

@ -145,20 +145,19 @@ contains
end subroutine bin_to_mesh_indices
!===============================================================================
! COUNT_FISSION_SITES determines the number of fission bank sites in each cell
! of a given mesh as well as an optional energy group structure. This can be
! used for a variety of purposes (Shannon entropy, CMFD, uniform fission source
! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a
! given mesh as well as an optional energy group structure. This can be used for
! a variety of purposes (Shannon entropy, CMFD, uniform fission source
! weighting)
!===============================================================================
subroutine count_fission_sites(m, bank_array, cnt, total, &
energies, size_bank, sites_outside)
subroutine count_bank_sites(m, bank_array, cnt, energies, size_bank, &
sites_outside)
type(StructuredMesh), pointer :: m ! mesh to count sites
type(Bank), intent(in) :: bank_array(:) ! fission or source bank
real(8), intent(out) :: cnt(:,:,:,:) ! weight of sites in each
! cell and energy group
real(8), intent(out) :: total ! total weight of sites
real(8), optional :: energies(:) ! energy grid to search
integer(8), optional :: size_bank ! # of bank sites (on each proc)
logical, optional :: sites_outside ! were there sites outside mesh?
@ -168,7 +167,6 @@ contains
integer :: ijk(3) ! indices on mesh
integer :: n_groups ! number of groups in energies
integer :: e_bin ! energy_bin
real(8) :: weight ! accumulated weight of sites
logical :: in_mesh ! was single site outside mesh?
logical :: outside ! was any site outside mesh?
#ifdef MPI
@ -177,7 +175,6 @@ contains
! initialize variables
cnt = ZERO
weight = ZERO
outside = .false.
! Set size of bank
@ -205,9 +202,6 @@ contains
cycle
end if
! add weight
weight = weight + ONE
! determine energy bin
if (present(energies)) then
if (bank_array(i) % E < energies(1)) then
@ -222,8 +216,8 @@ contains
end if
! add to appropriate mesh box
! TODO: if tracking weight through bank, add weight instead
cnt(e_bin,ijk(1),ijk(2),ijk(3)) = cnt(e_bin,ijk(1),ijk(2),ijk(3)) + 1
cnt(e_bin,ijk(1),ijk(2),ijk(3)) = cnt(e_bin,ijk(1),ijk(2),ijk(3)) + &
bank_array(i) % wgt
end do FISSION_SITES
#ifdef MPI
@ -245,15 +239,11 @@ contains
MPI_COMM_WORLD, mpi_err)
end if
! determine total weight of bank sites
call MPI_REDUCE(weight, total, 1, MPI_REAL8, MPI_SUM, 0, &
MPI_COMM_WORLD, mpi_err)
#else
total = weight
sites_outside = outside
#endif
end subroutine count_fission_sites
end subroutine count_bank_sites
!===============================================================================
! MESH_INTERSECT determines if a line between xyz0 and xyz1 intersects the outer

3
src/mesh_header.F90
View file

@ -3,7 +3,7 @@ module mesh_header
implicit none
!===============================================================================
! STRUCTUREDMESH represents a tesslation of n-dimensional Euclidean space by
! STRUCTUREDMESH represents a tessellation of n-dimensional Euclidean space by
! congruent squares or cubes
!===============================================================================
@ -11,6 +11,7 @@ module mesh_header
integer :: id ! user-specified id
integer :: type ! rectangular, hexagonal
integer :: n_dimension ! rank of mesh
real(8) :: volume_frac ! volume fraction of each cell
integer, allocatable :: dimension(:) ! number of cells in each direction
real(8), allocatable :: lower_left(:) ! lower-left corner of mesh
real(8), allocatable :: upper_right(:) ! upper-right corner of mesh

35
src/output.F90
View file

@ -138,6 +138,41 @@ contains
end subroutine header
!===============================================================================
! PRINT_VERSION shows the current version as well as copright and license
! information
!===============================================================================
subroutine print_version()
if (master) then
write(UNIT=OUTPUT_UNIT, FMT='(1X,A,1X,I1,".",I1,".",I1)') &
"OpenMC version", VERSION_MAJOR, VERSION_MINOR, VERSION_RELEASE
write(UNIT=OUTPUT_UNIT, FMT=*) "Copyright (c) 2011-2012 &
&Massachusetts Institute of Technology"
write(UNIT=OUTPUT_UNIT, FMT=*) "MIT/X license at &
&<http://mit-crpg.github.com/openmc/license.html>"
end if
end subroutine print_version
!===============================================================================
! PRINT_USAGE displays information about command line usage of OpenMC
!===============================================================================
subroutine print_usage()
if (master) then
write(OUTPUT_UNIT,*) 'Usage: openmc [options] [directory]'
write(OUTPUT_UNIT,*)
write(OUTPUT_UNIT,*) 'Options:'
write(OUTPUT_UNIT,*) ' -p, --plot Run in plotting mode'
write(OUTPUT_UNIT,*) ' -v, --version Show version information'
write(OUTPUT_UNIT,*) ' -?, --help Show this message'
end if
end subroutine print_usage
!===============================================================================
! WRITE_MESSAGE displays an informational message to the log file and the
! standard output stream.

38
src/physics.F90
View file

@ -12,6 +12,7 @@ module physics
use global
use interpolation, only: interpolate_tab1
use material_header, only: Material
use mesh, only: get_mesh_indices
use output, only: write_message
use particle_header, only: LocalCoord
use random_lcg, only: prn
@ -839,6 +840,7 @@ contains
integer :: nu ! actual number of neutrons produced
integer :: law ! energy distribution law
integer :: n_sample ! number of times resampling
integer :: ijk(3) ! indices in ufs mesh
real(8) :: E ! incoming energy of neutron
real(8) :: E_out ! outgoing energy of fission neutron
real(8) :: nu_t ! total nu
@ -850,7 +852,9 @@ contains
real(8) :: xi ! random number
real(8) :: yield ! delayed neutron precursor yield
real(8) :: prob ! cumulative probability
real(8) :: weight ! weight adjustment for ufs method
logical :: actual_event ! did fission actually occur? (no survival biasing)
logical :: in_mesh ! source site in ufs mesh?
type(Nuclide), pointer :: nuc
type(DistEnergy), pointer :: edist => null()
@ -886,12 +890,32 @@ contains
! TODO: Heat generation from fission
! If uniform fission source weighting is turned on, we increase of decrease
! the expected number of fission sites produced
if (ufs) then
! Determine indices on ufs mesh for current location
call get_mesh_indices(ufs_mesh, p % coord0 % xyz, ijk, in_mesh)
if (.not. in_mesh) then
message = "Source site outside UFS mesh!"
call fatal_error()
end if
if (source_frac(1,ijk(1),ijk(2),ijk(3)) /= ZERO) then
weight = ufs_mesh % volume_frac / source_frac(1,ijk(1),ijk(2),ijk(3))
else
weight = ONE
end if
else
weight = ONE
end if
! Sample number of neutrons produced
if (actual_event) then
nu_t = p % wgt / keff * nu_t
nu_t = p % wgt / keff * nu_t * weight
else
nu_t = p % last_wgt * micro_xs(index_nuclide) % fission / (keff * &
micro_xs(index_nuclide) % total) * nu_t
micro_xs(index_nuclide) % total) * nu_t * weight
end if
if (prn() > nu_t - int(nu_t)) then
nu = int(nu_t)
@ -899,6 +923,9 @@ contains
nu = int(nu_t) + 1
end if
! Add to analog estimate of keff
call add_to_score(global_tallies(K_ANALOG), nu/weight * keff)
! Bank source neutrons
if (nu == 0 .or. n_bank == 3*work) return
do i = int(n_bank,4) + 1, int(min(n_bank + nu, 3*work),4)
@ -906,6 +933,13 @@ contains
fission_bank(i) % id = p % id
fission_bank(i) % xyz = p % coord0 % xyz
! Set weight of fission bank site
if (ufs) then
fission_bank(i) % wgt = ONE/weight
else
fission_bank(i) % wgt = ONE
end if
! sample cosine of angle
mu = sample_angle(rxn, E)

7
src/source.F90
View file

@ -94,6 +94,9 @@ contains
id = bank_first + j - 1
source_bank(j) % id = id
! Set weight to one
source_bank(j) % wgt = ONE
! initialize random number seed
call set_particle_seed(id)
@ -148,6 +151,8 @@ contains
src => source_bank(index_source)
! copy attributes from source bank site
p % wgt = src % wgt
p % last_wgt = src % wgt
p % coord % xyz = src % xyz
p % coord % uvw = src % uvw
p % last_xyz = src % xyz
@ -315,7 +320,7 @@ contains
call fatal_error()
end if
! Read position, angle, and energy
! Read source sites
read(UNIT=UNIT_SOURCE) source_bank(1:work)
! Close binary source file

147
src/tally.F90
View file

@ -1014,8 +1014,11 @@ contains
call get_mesh_indices(m, xyz0, ijk0, start_in_mesh)
call get_mesh_indices(m, xyz1, ijk1, end_in_mesh)
! Check to make sure start or end is in mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) cycle
! Check to if start or end is in mesh -- if not, check if track still
! intersects with mesh
if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then
if (.not. mesh_intersects(m, xyz0, xyz1)) cycle
end if
! Calculate number of surface crossings
n_cross = sum(abs(ijk1 - ijk0))
@ -1295,46 +1298,6 @@ contains
end function get_next_bin
!===============================================================================
! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and
! specific response function. Note that we don't need to add the square of the
! contribution since that is done at the cycle level, not the history level
!===============================================================================
subroutine add_to_score(score, val)
type(TallyScore), intent(inout) :: score
real(8), intent(in) :: val
score % n_events = score % n_events + 1
score % value = score % value + val
end subroutine add_to_score
!===============================================================================
! ACCUMULATE_BATCH_ESTIMATE accumulates scores from many histories (or many
! generations) into a single realization of a random variable.
!===============================================================================
elemental subroutine accumulate_batch_estimate(score)
type(TallyScore), intent(inout) :: score
real(8) :: val
! Add the sum and square of the sum of contributions from each cycle
! within a cycle to the variables sum and sum_sq. This will later allow us
! to calculate a variance on the tallies
val = score % value/(n_particles*gen_per_batch)
score % sum = score % sum + val
score % sum_sq = score % sum_sq + val*val
! Reset the single batch estimate
score % value = ZERO
end subroutine accumulate_batch_estimate
!===============================================================================
! SYNCHRONIZE_TALLIES accumulates the sum of the contributions from each history
! within the batch to a new random variable
@ -1354,7 +1317,7 @@ contains
t => tallies(i)
! Loop over all filter and scoring bins
call accumulate_batch_estimate(t % scores)
call accumulate_score(t % scores)
end do
end subroutine synchronize_tallies
@ -1790,25 +1753,6 @@ contains
end function get_label
!===============================================================================
! CALCULATE_STATISTICS determines the sample mean and the standard deviation of
! the mean for a TallyScore.
!===============================================================================
elemental subroutine calculate_statistics(score)
type(TallyScore), intent(inout) :: score
! Calculate sample mean and standard deviation of the mean -- note that we
! 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))
end subroutine calculate_statistics
!===============================================================================
! TALLY_STATISTICS computes the mean and standard deviation of the mean of each
! tally and stores them in the val and val_sq attributes of the TallyScores
@ -1824,14 +1768,89 @@ contains
do i = 1, n_tallies
t => tallies(i)
call calculate_statistics(t % scores)
call statistics_score(t % scores)
end do
! Calculate statistics for global tallies
call calculate_statistics(global_tallies)
call statistics_score(global_tallies)
end subroutine tally_statistics
!===============================================================================
! ADD_TO_SCORE accumulates a scoring contribution to a specific tally bin and
! specific response function. Note that we don't need to add the square of the
! contribution since that is done at the cycle level, not the history level
!===============================================================================
subroutine add_to_score(score, val)
type(TallyScore), intent(inout) :: score
real(8), intent(in) :: val
score % n_events = score % n_events + 1
score % value = score % value + val
end subroutine add_to_score
!===============================================================================
! ACCUMULATE_SCORE accumulates scores from many histories (or many generations)
! into a single realization of a random variable.
!===============================================================================
elemental subroutine accumulate_score(score)
type(TallyScore), intent(inout) :: score
real(8) :: val
! Add the sum and square of the sum of contributions from each cycle
! within a cycle to the variables sum and sum_sq. This will later allow us
! to calculate a variance on the tallies
val = score % value/(n_particles*gen_per_batch)
score % sum = score % sum + val
score % sum_sq = score % sum_sq + val*val
! Reset the single batch estimate
score % value = ZERO
end subroutine accumulate_score
!===============================================================================
! STATISTICS_SCORE determines the sample mean and the standard deviation of the
! mean for a TallyScore.
!===============================================================================
elemental subroutine statistics_score(score)
type(TallyScore), intent(inout) :: score
! Calculate sample mean and standard deviation of the mean -- note that we
! 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))
end subroutine statistics_score
!===============================================================================
! RESET_SCORE zeroes out the value and accumulated sum and sum-squared for a
! single TallyScore.
!===============================================================================
elemental subroutine reset_score(score)
type(TallyScore), intent(inout) :: score
score % n_events = 0
score % value = ZERO
score % sum = ZERO
score % sum_sq = ZERO
end subroutine reset_score
!==============================================================================
! TALLY_RESET resets all of the user-defined tallies
!==============================================================================

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

@ -21,7 +21,7 @@
<component name="weight_avg" type="double" default="1.0" />
</typedef>
<typedef name="entropy_xml">
<typedef name="mesh_xml">
<component name="dimension" type="integer-array" />
<component name="lower_left" type="double-array" />
<component name="upper_right" type="double-array" />
@ -31,13 +31,14 @@
<variable name="cross_sections_" tag="cross_sections" type="word" length="255" />
<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="entropy_xml" dimension="1" />
<variable name="entropy_" tag="entropy" type="mesh_xml" dimension="1" />
<variable name="ptables_" tag="ptables" type="word" length="3" />
<variable name="seed_" tag="seed" type="integer" />
<variable name="source_" tag="source" type="source_xml" />
<variable name="survival_" tag="survival_biasing" type="word" length="3" />
<variable name="verbosity_" tag="verbosity" type="integer" />
<variable name="trace_" tag="trace" type="integer-array" />
<variable name="uniform_fs_" tag="uniform_fs" type="mesh_xml" dimension="1" />
<variable name="write_source_" tag="write_source" type="word" length="3" />
</template>