diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES
index d6a8591c6c..7d70c783ef 100644
--- a/src/DEPENDENCIES
+++ b/src/DEPENDENCIES
@@ -179,6 +179,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
diff --git a/src/global.F90 b/src/global.F90
index dae61fa6ad..5da8782b5c 100644
--- a/src/global.F90
+++ b/src/global.F90
@@ -156,6 +156,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.
diff --git a/src/input_xml.F90 b/src/input_xml.F90
index 039f8cb653..1b7f423c44 100644
--- a/src/input_xml.F90
+++ b/src/input_xml.F90
@@ -236,6 +236,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.
@@ -812,6 +866,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
diff --git a/src/intercycle.F90 b/src/intercycle.F90
index be5fa749af..d44e00c7f3 100644
--- a/src/intercycle.F90
+++ b/src/intercycle.F90
@@ -446,4 +446,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, total, &
+ 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
+
+ ! Normalize to total weight to get fraction of source in each cell
+ if (master) source_frac = source_frac / total
+
+#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)
+ call MPI_BCAST(total, 1, MPI_REAL8, 0, MPI_COMM_WORLD, mpi_err)
+#endif
+
+ ! 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
diff --git a/src/main.F90 b/src/main.F90
index bf5b922bc4..9cd90d6d0e 100644
--- a/src/main.F90
+++ b/src/main.F90
@@ -2,16 +2,17 @@ program main
use constants
use global
- use finalize, only: finalize_run
- use initialize, only: initialize_run
- use intercycle, only: shannon_entropy, calculate_keff, synchronize_bank
- 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
@@ -78,6 +79,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
@@ -106,11 +110,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)
diff --git a/src/mesh_header.F90 b/src/mesh_header.F90
index 8088972ce0..a26178b6f3 100644
--- a/src/mesh_header.F90
+++ b/src/mesh_header.F90
@@ -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
diff --git a/src/physics.F90 b/src/physics.F90
index 00c2508ce4..382cce60b3 100644
--- a/src/physics.F90
+++ b/src/physics.F90
@@ -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)
@@ -907,7 +934,11 @@ contains
fission_bank(i) % xyz = p % coord0 % xyz
! Set weight of fission bank site
- fission_bank(i) % wgt = ONE
+ 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)
diff --git a/src/xml-fortran/templates/settings_t.xml b/src/xml-fortran/templates/settings_t.xml
index f2b387316d..321565955f 100644
--- a/src/xml-fortran/templates/settings_t.xml
+++ b/src/xml-fortran/templates/settings_t.xml
@@ -20,7 +20,7 @@
-
+
@@ -30,13 +30,14 @@
-
+
+