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

LCG approach Al.5(by Paul): skip f(ZAID) states starting from xs_seed

Browse source 
To guarantee random numbers are not re-used, advance the seed N times from its original position after energy changed.
0.	Initialization of f and N
0.	At the beginning of a particle life(including secondary particle):
        xs_seed = tracking_seed
1.	When calculating xs:
        Xi_urr(ZAID) = prn(skipping f(ZAID) from xs_seed)
2.	If particle changes energy,
        xs_seed = prn_seed (skipping N from xs_seed)
where N is the number of nuclides which have different zaid, f is a map of zaid.
This commit is contained in:
jingang 2016-03-04 10:35:13 -05:00
parent 4686a4ce88
commit df4cc8e0f2
5 changed files with 48 additions and 32 deletions
Download patch file Download diff file Expand all files Collapse all files

12
src/cross_section.F90
View file

@ -10,7 +10,7 @@ module cross_section
use material_header, only: Material
use nuclide_header
use particle_header, only: Particle
use random_lcg, only: prn, prn_ahead
use random_lcg, only: prn, get_prn_ahead
use sab_header, only: SAlphaBeta
use search, only: binary_search
@ -387,15 +387,13 @@ contains
! sample probability table using the cumulative distribution
! random numbers for xs calculation are sampled in a way separate from
! tracking. 'xs_seed' is a copy of normal tracking prn seed but updated
! until the particle undergoes a scattering event. Random number is
! calculated by skipping ahead 'xs_seed + ZZAAA'(zaid) times from the seed
! 'xs_seed' + 'ZZAAA'.
! Random numbers for xs calculation are sampled by skipping ahead
! f(zaid) times from the seed 'xs_seed' + 'zaid'.
! This guarantees the randomness and, at the same time, makes sure we reuse
! random number for the same nuclide at different temperatures, therefore
! preserving correlation of temperature in probability tables.
r = prn_ahead(xs_seed + nuc % zaid, xs_seed + nuc % zaid)
r = get_prn_ahead(int(nuc_zaid_dict % get_key(nuc % zaid), 8), &
xs_seed + nuc % zaid)
i_low = 1
do

11
src/global.F90
View file

@ -105,11 +105,16 @@ module global
integer :: default_expand = ENDF_BVII1
! Random number seed for cross sections, specially for URR ptables
! This number is copied from normal tracking random number sequence but
! updated until the particle undergoes a scattering event. It is shared for
! all nuclides.
! This number is shared by all nuclides and updated after particle
! changed its energy.
integer(8) :: xs_seed = 1_8
! Dictionary to look up the skip distance to get prn when sampling URR
type(DictIntInt) :: nuc_zaid_dict
! Total amount of nuclide zaid instances
integer(8) :: n_nuc_zaid_total
!$omp threadprivate(xs_seed)
! ============================================================================

41
src/input_xml.F90
View file

@ -1891,21 +1891,23 @@ contains
subroutine read_materials_xml()
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: k ! loop index for elements
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: n_nuc_ele ! number of nuclides in an element
integer :: index_list ! index in xs_listings array
integer :: index_nuclide ! index in nuclides
integer :: index_sab ! index in sab_tables
real(8) :: val ! value entered for density
real(8) :: temp_dble ! temporary double prec. real
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
integer :: i ! loop index for materials
integer :: j ! loop index for nuclides
integer :: k ! loop index for elements
integer :: n ! number of nuclides
integer :: n_sab ! number of sab tables for a material
integer :: n_nuc_ele ! number of nuclides in an element
integer :: index_list ! index in xs_listings array
integer :: index_nuclide ! index in nuclides
integer :: index_nuc_zaid ! index in nuclide ZAID
integer :: index_sab ! index in sab_tables
real(8) :: val ! value entered for density
real(8) :: temp_dble ! temporary double prec. real
logical :: file_exists ! does materials.xml exist?
logical :: sum_density ! density is taken to be sum of nuclide densities
integer :: zaid ! ZAID of nuclide
character(12) :: name ! name of isotope, e.g. 92235.03c
character(12) :: alias ! alias of nuclide, e.g. U-235.03c
character(MAX_WORD_LEN) :: units ! units on density
character(MAX_LINE_LEN) :: filename ! absolute path to materials.xml
character(MAX_LINE_LEN) :: temp_str ! temporary string when reading
@ -1955,6 +1957,7 @@ contains
! Initialize count for number of nuclides/S(a,b) tables
index_nuclide = 0
index_nuc_zaid = 0
index_sab = 0
do i = 1, n_materials
@ -2300,6 +2303,7 @@ contains
index_list = xs_listing_dict % get_key(to_lower(name))
name = xs_listings(index_list) % name
alias = xs_listings(index_list) % alias
zaid = xs_listings(index_list) % zaid
! If this nuclide hasn't been encountered yet, we need to add its name
! and alias to the nuclide_dict
@ -2313,6 +2317,12 @@ contains
mat % nuclide(j) = nuclide_dict % get_key(to_lower(name))
end if
! Construct dict of nuclide zaid
if (.not. nuc_zaid_dict % has_key(zaid)) then
index_nuc_zaid = index_nuc_zaid + 1
call nuc_zaid_dict % add_key(zaid, index_nuc_zaid)
end if
! Copy name and atom/weight percent
mat % names(j) = name
mat % atom_density(j) = list_density % get_item(j)
@ -2407,6 +2417,7 @@ contains
! Set total number of nuclides and S(a,b) tables
n_nuclides_total = index_nuclide
n_sab_tables = index_sab
n_nuc_zaid_total = index_nuc_zaid
! Close materials XML file
call close_xmldoc(doc)

9
src/random_lcg.F90
View file

@ -24,7 +24,8 @@ module random_lcg
!$omp threadprivate(prn_seed, stream)
public :: prn
public :: prn_ahead
public :: get_prn_ahead
public :: prn_skip_ahead
public :: initialize_prng
public :: set_particle_seed
public :: prn_skip
@ -54,11 +55,11 @@ contains
end function prn
!===============================================================================
! PRN_AHEAD generates a pseudo-random number which is 'n' times ahead from a
! GET_PRN_AHEAD generates a pseudo-random number which is 'n' times ahead from a
! specific seed. This function does not changed current LCG status.
!===============================================================================
function prn_ahead(n, seed) result(pseudo_rn)
function get_prn_ahead(n, seed) result(pseudo_rn)
integer(8), intent(in) :: n ! number of prns to skip
integer(8), intent(in) :: seed ! starting seed
@ -70,7 +71,7 @@ contains
pseudo_rn = prn_skip_ahead(n, seed) * prn_norm
end function prn_ahead
end function get_prn_ahead
!===============================================================================
! INITIALIZE_PRNG sets up the random number generator, determining the seed and

7
src/tracking.F90
View file

@ -12,7 +12,7 @@ module tracking
use particle_header, only: LocalCoord, Particle
use physics, only: collision
use physics_mg, only: collision_mg
use random_lcg, only: prn, prn_seed
use random_lcg, only: prn, prn_seed, prn_skip_ahead
use string, only: to_str
use tally, only: score_analog_tally, score_tracklength_tally, &
score_collision_tally, score_surface_current
@ -200,8 +200,9 @@ contains
! re-evaluated
p % last_material = NONE
! Update xs_seed to be current tracking seed after a collision
if (p % E /= p % last_E) xs_seed = prn_seed(STREAM_TRACKING)
! Advance xs_seed N times ahead to avoid re-using prn
if (p % E /= p % last_E) &
xs_seed = prn_skip_ahead(n_nuc_zaid_total, xs_seed)
! Set all uvws to base level -- right now, after a collision, only the
! base level uvws are changed