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

Uniform fission source weighting. Closes #26 on github.

Browse source
This commit is contained in:
Paul Romano 2012-03-31 17:22:03 -04:00
parent 1fb88c1360
commit 7819bef5bc
8 changed files with 167 additions and 21 deletions
Download patch file Download diff file Expand all files Collapse all files

1
src/DEPENDENCIES
View file

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

5
src/global.F90
View file

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

57
src/input_xml.F90
View file

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

51
src/intercycle.F90
View file

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

29
src/main.F90
View file

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

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

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

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

@ -20,7 +20,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" />
@ -30,13 +30,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>