From df4cc8e0f21458703ef6264c57d2dbe529de2c40 Mon Sep 17 00:00:00 2001 From: jingang Date: Fri, 4 Mar 2016 10:35:13 -0500 Subject: [PATCH] LCG approach Al.5(by Paul): skip f(ZAID) states starting from xs_seed 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. --- src/cross_section.F90 | 12 +++++------- src/global.F90 | 11 ++++++++--- src/input_xml.F90 | 41 ++++++++++++++++++++++++++--------------- src/random_lcg.F90 | 9 +++++---- src/tracking.F90 | 7 ++++--- 5 files changed, 48 insertions(+), 32 deletions(-) diff --git a/src/cross_section.F90 b/src/cross_section.F90 index f059c9c07..a66870714 100644 --- a/src/cross_section.F90 +++ b/src/cross_section.F90 @@ -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 diff --git a/src/global.F90 b/src/global.F90 index 9f5bd6567..93fe3c598 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -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) ! ============================================================================ diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 03ef8dcbc..47b0aaacf 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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) diff --git a/src/random_lcg.F90 b/src/random_lcg.F90 index 1e92fa6ce..755ee673f 100644 --- a/src/random_lcg.F90 +++ b/src/random_lcg.F90 @@ -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 diff --git a/src/tracking.F90 b/src/tracking.F90 index b17f8ba59..86f0eded5 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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