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

Rewrote random number generator.

Browse source
This commit is contained in:
Paul Romano 2011-11-20 19:05:54 -05:00
parent 0fc7491730
commit ca1d53fbc8
9 changed files with 232 additions and 715 deletions
Download patch file Download diff file Expand all files Collapse all files

10
src/DEPENDENCIES
View file

@ -72,9 +72,9 @@ initialize.o: geometry_header.o
initialize.o: global.o
initialize.o: input_xml.o
initialize.o: logging.o
initialize.o: mcnp_random.o
initialize.o: mpi_routines.o
initialize.o: output.o
initialize.o: random_lcg.o
initialize.o: source.o
initialize.o: string.o
initialize.o: tally.o
@ -108,12 +108,12 @@ logging.o: global.o
main.o: constants.o
main.o: global.o
main.o: initialize.o
main.o: mcnp_random.o
main.o: mpi_routines.o
main.o: output.o
main.o: particle_header.o
main.o: physics.o
main.o: plot.o
main.o: random_lcg.o
main.o: source.o
main.o: string.o
main.o: tally.o
@ -124,9 +124,9 @@ mesh.o: mesh_header.o
mpi_routines.o: constants.o
mpi_routines.o: error.o
mpi_routines.o: global.o
mpi_routines.o: mcnp_random.o
mpi_routines.o: output.o
mpi_routines.o: particle_header.o
mpi_routines.o: random_lcg.o
mpi_routines.o: tally_header.o
output.o: constants.o
@ -149,9 +149,9 @@ physics.o: geometry.o
physics.o: geometry_header.o
physics.o: global.o
physics.o: interpolation.o
physics.o: mcnp_random.o
physics.o: output.o
physics.o: particle_header.o
physics.o: random_lcg.o
physics.o: search.o
physics.o: string.o
physics.o: tally.o
@ -171,8 +171,8 @@ source.o: bank_header.o
source.o: constants.o
source.o: cross_section_header.o
source.o: global.o
source.o: mcnp_random.o
source.o: output.o
source.o: random_lcg.o
source.o: particle_header.o
source.o: physics.o

2
src/OBJECTS
View file

@ -20,7 +20,6 @@ input_xml.o \
logging.o \
main.o \
material_header.o \
mcnp_random.o \
mesh_header.o \
mesh.o \
mpi_routines.o \
@ -28,6 +27,7 @@ output.o \
particle_header.o \
physics.o \
plot.o \
random_lcg.o \
search.o \
source.o \
source_header.o \

8
src/initialize.f90
View file

@ -13,10 +13,10 @@ module initialize
use input_xml, only: read_input_xml, read_cross_sections_xml, &
cells_in_univ_dict
use logging, only: create_log
use mcnp_random, only: RN_init_problem
use mpi_routines, only: setup_mpi
use output, only: title, echo_input, message, print_summary, &
print_particle, header, print_plot
use random_lcg, only: initialize_prng
use source, only: initialize_source
use string, only: int_to_str, starts_with, ends_with, lower_case
use tally, only: create_tally_map, TallyObject
@ -51,10 +51,8 @@ contains
! Print initialization header block
if (master) call header("INITIALIZATION", 1)
! Initialize random number generator. The first argument corresponds to
! which random number generator to use- in this case one of the L'Ecuyer
! 63-bit RNGs.
call RN_init_problem(3, 0_8, 0_8, 0_8, 0)
! Initialize random number generator
call initialize_prng()
! set up dictionaries
call create_dictionaries()

7
src/main.f90
View file

@ -3,17 +3,16 @@ program main
use constants
use global
use initialize, only: initialize_run
use mcnp_random, only: RN_init_particle
use mpi_routines, only: synchronize_bank
use output, only: write_message, header, print_runtime
use particle_header, only: Particle
use plot, only: run_plot
use physics, only: transport
use tally, only: calculate_keff
use random_lcg, only: set_particle_seed
use source, only: get_source_particle
use string, only: int_to_str
use tally, only: synchronize_tallies, write_tallies, &
tally_statistics
tally_statistics, calculate_keff
use timing, only: timer_start, timer_stop
#ifdef MPI
@ -87,7 +86,7 @@ contains
! set random number seed
i_particle = (i_cycle-1)*n_particles + p % id
call RN_init_particle(i_particle)
call set_particle_seed(i_particle)
! transport particle
call transport(p)

627
src/mcnp_random.f90
View file

@ -1,627 +0,0 @@
module mcnp_random
!=======================================================================
! Description:
! mcnp_random.F90 -- random number generation routines
!=======================================================================
! This module contains:
!
! * Constants for the RN generator, including initial RN seed for the
! problem & the current RN seed
!
! * MCNP interface routines:
! - random number function: rang()
! - RN initialization for problem: RN_init_problem
! - RN initialization for particle: RN_init_particle
! - get info on RN parameters: RN_query
! - get RN seed for n-th history: RN_query_first
! - set new RN parameters: RN_set
! - skip-ahead in the RN sequence: RN_skip_ahead
! - Unit tests: RN_test_basic, RN_test_skip, RN_test_mixed
!
! * For interfacing with the rest of MCNP, arguments to/from these
! routines will have types of I8 or I4.
! Any args which are to hold random seeds, multipliers,
! skip-distance will be type I8, so that 63 bits can be held without
! truncation.
!
! Revisions:
! * 10-04-2001 - F Brown, initial mcnp version
! * 06-06-2002 - F Brown, mods for extended generators
! * 12-21-2004 - F Brown, added 3 of LeCuyer's 63-bit mult. RNGs
! * 01-29-2005 - J Sweezy, Modify to use mcnp modules prior to automatic
! io unit numbers.
! * 12-02-2005 - F Brown, mods for consistency with C version
!=======================================================================
!-------------------
! MCNP output units
!-------------------
integer, parameter :: iuo = 6
integer, parameter :: jtty = 6
PRIVATE
!---------------------------------------------------
! Kinds for LONG INTEGERS (64-bit) & REAL*8 (64-bit)
!---------------------------------------------------
integer, parameter :: R8 = selected_real_kind(15,307)
integer, parameter :: I8 = selected_int_kind(18)
!-----------------------------------
! Public functions and subroutines for this module
!-----------------------------------
PUBLIC :: rang
PUBLIC :: RN_init_problem
PUBLIC :: RN_init_particle
PUBLIC :: RN_set
PUBLIC :: RN_query
PUBLIC :: RN_query_first
PUBLIC :: RN_update_stats
PUBLIC :: RN_test_basic
PUBLIC :: RN_test_skip
PUBLIC :: RN_test_mixed
PUBLIC :: RN_skip
!-------------------------------------
! Constants for standard RN generators
!-------------------------------------
type :: RN_GEN
integer :: index
integer(I8) :: mult ! generator (multiplier)
integer(I8) :: add ! additive constant
integer :: log2mod ! log2 of modulus, must be <64
integer(I8) :: stride ! stride for particle skip-ahead
integer(I8) :: initseed ! default seed for problem
character(len=8) :: name
end type RN_GEN
! parameters for standard generators
integer, parameter :: n_RN_GEN = 7
type(RN_GEN), SAVE :: standard_generator(n_RN_GEN)
data standard_generator / &
& RN_GEN( 1, 19073486328125_I8, 0_I8, 48, 152917_I8, 19073486328125_I8 , 'mcnp std' ), &
& RN_GEN( 2, 9219741426499971445_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer1' ), &
& RN_GEN( 3, 2806196910506780709_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer2' ), &
& RN_GEN( 4, 3249286849523012805_I8, 1_I8, 63, 152917_I8, 1_I8, 'LEcuyer3' ), &
& RN_GEN( 5, 3512401965023503517_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer4' ), &
& RN_GEN( 6, 2444805353187672469_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer5' ), &
& RN_GEN( 7, 1987591058829310733_I8, 0_I8, 63, 152917_I8, 1_I8, 'LEcuyer6' ) &
& /
!-----------------------------------------------------------------
! * Linear multiplicative congruential RN algorithm:
!
! RN_SEED = RN_SEED*RN_MULT + RN_ADD mod RN_MOD
!
! * Default values listed below will be used, unless overridden
!-----------------------------------------------------------------
integer, SAVE :: RN_INDEX = 1
integer(I8), SAVE :: RN_MULT = 19073486328125_I8
integer(I8), SAVE :: RN_ADD = 0_I8
integer, SAVE :: RN_BITS = 48
integer(I8), SAVE :: RN_STRIDE = 152917_I8
integer(I8), SAVE :: RN_SEED0 = 19073486328125_I8
integer(I8), SAVE :: RN_MOD = 281474976710656_I8
integer(I8), SAVE :: RN_MASK = 281474976710655_I8
integer(I8), SAVE :: RN_PERIOD = 70368744177664_I8
real(R8), SAVE :: RN_NORM = 1._R8 / 281474976710656._R8
!------------------------------------
! Private data for a single particle
!------------------------------------
integer(I8), save :: RN_SEED = 19073486328125_I8 ! current seed
integer(I8), save :: RN_COUNT = 0_I8 ! current counter
integer(I8), save :: RN_NPS = 0_I8 ! current particle number
!$OMP THREADPRIVATE (RN_SEED, RN_COUNT, RN_NPS)
!------------------------------------------
! Shared data, to collect info on RN usage
!------------------------------------------
integer(I8), SAVE :: RN_COUNT_TOTAL = 0 ! total RN count all particles
integer(I8), SAVE :: RN_COUNT_STRIDE = 0 ! count for stride exceeded
integer(I8), SAVE :: RN_COUNT_MAX = 0 ! max RN count all particles
integer(I8), SAVE :: RN_COUNT_MAX_NPS = 0 ! part index for max count
!---------------------------------------------------------------------
! Reference data: Seeds for case of init.seed = 1,
! Seed numbers for index 1-5, 123456-123460
!---------------------------------------------------------------------
integer(I8), dimension(10,n_RN_GEN) :: RN_CHECK
data RN_CHECK / &
! ***** 1 ***** mcnp standard gen *****
& 19073486328125_I8, 29763723208841_I8, 187205367447973_I8, &
& 131230026111313_I8, 264374031214925_I8, 260251000190209_I8, &
& 106001385730621_I8, 232883458246025_I8, 97934850615973_I8, &
& 163056893025873_I8, &
! ***** 2 *****
& 9219741426499971446_I8, 666764808255707375_I8, 4935109208453540924_I8, &
& 7076815037777023853_I8, 5594070487082964434_I8, 7069484152921594561_I8, &
& 8424485724631982902_I8, 19322398608391599_I8, 8639759691969673212_I8, &
& 8181315819375227437_I8, &
! ***** 3 *****
& 2806196910506780710_I8, 6924308458965941631_I8, 7093833571386932060_I8, &
& 4133560638274335821_I8, 678653069250352930_I8, 6431942287813238977_I8, &
& 4489310252323546086_I8, 2001863356968247359_I8, 966581798125502748_I8, &
& 1984113134431471885_I8, &
! ***** 4 *****
& 3249286849523012806_I8, 4366192626284999775_I8, 4334967208229239068_I8, &
& 6386614828577350285_I8, 6651454004113087106_I8, 2732760390316414145_I8, &
& 2067727651689204870_I8, 2707840203503213343_I8, 6009142246302485212_I8, &
& 6678916955629521741_I8, &
! ***** 5 *****
& 3512401965023503517_I8, 5461769869401032777_I8, 1468184805722937541_I8, &
& 5160872062372652241_I8, 6637647758174943277_I8, 794206257475890433_I8, &
& 4662153896835267997_I8, 6075201270501039433_I8, 889694366662031813_I8, &
& 7299299962545529297_I8, &
! ***** 6 *****
& 2444805353187672469_I8, 316616515307798713_I8, 4805819485453690029_I8, &
& 7073529708596135345_I8, 3727902566206144773_I8, 1142015043749161729_I8, &
& 8632479219692570773_I8, 2795453530630165433_I8, 5678973088636679085_I8, &
& 3491041423396061361_I8, &
! ***** 7 *****
& 1987591058829310733_I8, 5032889449041854121_I8, 4423612208294109589_I8, &
& 3020985922691845009_I8, 5159892747138367837_I8, 8387642107983542529_I8, &
& 8488178996095934477_I8, 708540881389133737_I8, 3643160883363532437_I8, &
& 4752976516470772881_I8 /
!---------------------------------------------------------------------
CONTAINS
!-------------------------------------------------------------------
function rang()
! MCNP random number generator
!
! ***************************************
! ***** modifies RN_SEED & RN_COUNT *****
! ***************************************
implicit none
real(R8) :: rang
RN_SEED = iand( iand( RN_MULT*RN_SEED, RN_MASK) + RN_ADD, RN_MASK)
rang = RN_SEED * RN_NORM
RN_COUNT = RN_COUNT + 1
return
end function rang
!-------------------------------------------------------------------
subroutine RN_skip( skip )
! initialize MCNP random number parameters for particle "nps"
!
! * generate a new particle seed from the base seed
! & particle index
! * set the RN count to zero
implicit none
integer(I8), intent(in) :: skip
RN_SEED = RN_skip_ahead( RN_SEED, skip )
end subroutine RN_skip
!-------------------------------------------------------------------
function RN_skip_ahead( seed, skip )
! advance the seed "skip" RNs: seed*RN_MULT^n mod RN_MOD
implicit none
integer(I8) :: RN_skip_ahead
integer(I8), intent(in) :: seed, skip
integer(I8) :: nskip, gen, g, inc, c, gp, rn, seed_old
seed_old = seed
! add period till nskip>0
nskip = skip
do while( nskip<0_I8 )
if( RN_PERIOD>0_I8 ) then
nskip = nskip + RN_PERIOD
else
nskip = nskip + RN_MASK
nskip = nskip + 1_I8
endif
enddo
! get gen=RN_MULT^n, in log2(n) ops, not n ops !
nskip = iand( nskip, RN_MASK )
gen = 1
g = RN_MULT
inc = 0
c = RN_ADD
do while( nskip>0_I8 )
if( btest(nskip,0) ) then
gen = iand( gen*g, RN_MASK )
inc = iand( inc*g, RN_MASK )
inc = iand( inc+c, RN_MASK )
endif
gp = iand( g+1, RN_MASK )
g = iand( g*g, RN_MASK )
c = iand( gp*c, RN_MASK )
nskip = ishft( nskip, -1 )
enddo
rn = iand( gen*seed_old, RN_MASK )
rn = iand( rn + inc, RN_MASK )
RN_skip_ahead = rn
return
end function RN_skip_ahead
!-------------------------------------------------------------------
subroutine RN_init_problem( new_standard_gen, new_seed, &
& new_stride, new_part1, print_info )
! * initialize MCNP random number parameters for problem,
! based on user input. This routine should be called
! only from the main thread, if OMP threading is being used.
!
! * for initial & continue runs, these args should be set:
! new_standard_gen - index of built-in standard RN generator,
! from RAND gen= (or dbcn(14)
! new_seed - from RAND seed= (or dbcn(1))
! output - logical, print RN seed & mult if true
!
! new_stride - from RAND stride= (or dbcn(13))
! new_part1 - from RAND hist= (or dbcn(8))
!
! * for continue runs only, these should also be set:
! new_count_total - from "rnr" at end of previous run
! new_count_stride - from nrnh(1) at end of previous run
! new_count_max - from nrnh(2) at end of previous run
! new_count_max_nps - from nrnh(3) at end of previous run
!
! * check on size of long-ints & long-int arithmetic
! * check the multiplier
! * advance the base seed for the problem
! * set the initial particle seed
! * initialize the counters for RN stats
implicit none
integer, intent(in) :: new_standard_gen
integer(I8), intent(in) :: new_seed
integer(I8), intent(in) :: new_stride
integer(I8), intent(in) :: new_part1
integer, intent(in) :: print_info
character(len=20) :: printseed
integer(I8) :: itemp1, itemp2, itemp3, itemp4
if( new_standard_gen<1 .or. new_standard_gen>n_RN_GEN ) then
call expire( 'RN_init_problem', &
& ' ***** ERROR: illegal index for built-in RN generator')
endif
! set defaults, override if input supplied: seed, mult, stride
RN_INDEX = new_standard_gen
RN_MULT = standard_generator(RN_INDEX)%mult
RN_ADD = standard_generator(RN_INDEX)%add
RN_STRIDE = standard_generator(RN_INDEX)%stride
RN_SEED0 = standard_generator(RN_INDEX)%initseed
RN_BITS = standard_generator(RN_INDEX)%log2mod
RN_MOD = ishft( 1_I8, RN_BITS )
RN_MASK = ishft( not(0_I8), RN_BITS-64 )
RN_NORM = 2._R8**(-RN_BITS)
if( RN_ADD==0_I8) then
RN_PERIOD = ishft( 1_I8, RN_BITS-2 )
else
RN_PERIOD = ishft( 1_I8, RN_BITS )
endif
if( new_seed>0_I8 ) then
RN_SEED0 = new_seed
endif
if( new_stride>0_I8 ) then
RN_STRIDE = new_stride
endif
RN_COUNT_TOTAL = 0
RN_COUNT_STRIDE = 0
RN_COUNT_MAX = 0
RN_COUNT_MAX_NPS = 0
if( print_info /= 0 ) then
write(printseed,'(i20)') RN_SEED0
write( iuo,1) RN_INDEX, RN_SEED0, RN_MULT, RN_ADD, RN_BITS, RN_STRIDE
write(jtty,2) RN_INDEX, adjustl(printseed)
1 format( &
& /,' ***************************************************', &
& /,' * Random Number Generator = ',i20, ' *', &
& /,' * Random Number Seed = ',i20, ' *', &
& /,' * Random Number Multiplier = ',i20, ' *', &
& /,' * Random Number Adder = ',i20, ' *', &
& /,' * Random Number Bits Used = ',i20, ' *', &
& /,' * Random Number Stride = ',i20, ' *', &
& /,' ***************************************************',/)
2 format(' comment. using random number generator ',i2, &
& ', initial seed = ',a20)
endif
! double-check on number of bits in a long int
if( bit_size(RN_SEED)<64 ) then
call expire( 'RN_init_problem', &
& ' ***** ERROR: <64 bits in long-int, can-t generate RN-s')
endif
itemp1 = 5_I8**25
itemp2 = 5_I8**19
itemp3 = ishft(2_I8**62-1_I8,1) + 1_I8
itemp4 = itemp1*itemp2
if( iand(itemp4,itemp3)/=8443747864978395601_I8 ) then
call expire( 'RN_init_problem', &
& ' ***** ERROR: can-t do 64-bit integer ops for RN-s')
endif
if( new_part1>1_I8 ) then
! advance the problem seed to that for part1
RN_SEED0 = RN_skip_ahead( RN_SEED0, (new_part1-1_I8)*RN_STRIDE )
itemp1 = RN_skip_ahead( RN_SEED0, RN_STRIDE )
if( print_info /= 0 ) then
write(printseed,'(i20)') itemp1
write( iuo,3) new_part1, RN_SEED0, itemp1
write(jtty,4) new_part1, adjustl(printseed)
3 format( &
& /,' ***************************************************', &
& /,' * Random Number Seed will be advanced to that for *', &
& /,' * previous particle number = ',i20, ' *', &
& /,' * New RN Seed for problem = ',i20, ' *', &
& /,' * Next Random Number Seed = ',i20, ' *', &
& /,' ***************************************************',/)
4 format(' comment. advancing random number to particle ',i12, &
& ', initial seed = ',a20)
endif
endif
! set the initial particle seed
RN_SEED = RN_SEED0
RN_COUNT = 0
RN_NPS = 0
return
end subroutine RN_init_problem
!-------------------------------------------------------------------
subroutine RN_init_particle( nps )
! initialize MCNP random number parameters for particle "nps"
!
! * generate a new particle seed from the base seed
! & particle index
! * set the RN count to zero
implicit none
integer(I8), intent(in) :: nps
RN_SEED = RN_skip_ahead( RN_SEED0, nps*RN_STRIDE )
RN_COUNT = 0
RN_NPS = nps
return
end subroutine RN_init_particle
!-------------------------------------------------------------------
subroutine RN_set( key, value )
implicit none
character(len=*), intent(in) :: key
integer(I8), intent(in) :: value
character(len=20) :: printseed
integer(I8) :: itemp1
if( key == "stride" ) then
if( value>0_I8 ) then
RN_STRIDE = value
endif
endif
if( key == "count_total" ) RN_COUNT_TOTAL = value
if( key == "count_stride" ) RN_COUNT_STRIDE = value
if( key == "count_max" ) RN_COUNT_MAX = value
if( key == "count_max_nps" ) RN_COUNT_MAX_NPS = value
if( key == "seed" ) then
if( value>0_I8 ) then
RN_SEED0 = value
RN_SEED = RN_SEED0
RN_COUNT = 0
RN_NPS = 0
endif
endif
if( key == "part1" ) then
if( value>1_I8 ) then
! advance the problem seed to that for part1
RN_SEED0 = RN_skip_ahead( RN_SEED0, (value-1_I8)*RN_STRIDE )
itemp1 = RN_skip_ahead( RN_SEED0, RN_STRIDE )
write(printseed,'(i20)') itemp1
write( iuo,3) value, RN_SEED0, itemp1
write(jtty,4) value, adjustl(printseed)
3 format( &
& /,' ***************************************************', &
& /,' * Random Number Seed will be advanced to that for *', &
& /,' * previous particle number = ',i20, ' *', &
& /,' * New RN Seed for problem = ',i20, ' *', &
& /,' * Next Random Number Seed = ',i20, ' *', &
& /,' ***************************************************',/)
4 format(' comment. advancing random number to particle ',i12, &
& ', initial seed = ',a20)
RN_SEED = RN_SEED0
RN_COUNT = 0
RN_NPS = 0
endif
endif
return
end subroutine RN_set
!-------------------------------------------------------------------
function RN_query( key )
implicit none
integer(I8) :: RN_query
character(len=*), intent(in) :: key
RN_query = 0_I8
if( key == "seed" ) RN_query = RN_SEED
if( key == "stride" ) RN_query = RN_STRIDE
if( key == "mult" ) RN_query = RN_MULT
if( key == "add" ) RN_query = RN_ADD
if( key == "count" ) RN_query = RN_COUNT
if( key == "period" ) RN_query = RN_PERIOD
if( key == "count_total" ) RN_query = RN_COUNT_TOTAL
if( key == "count_stride" ) RN_query = RN_COUNT_STRIDE
if( key == "count_max" ) RN_query = RN_COUNT_MAX
if( key == "count_max_nps" ) RN_query = RN_COUNT_MAX_NPS
if( key == "first" ) RN_query = RN_SEED0
return
end function RN_query
!-------------------------------------------------------------------
function RN_query_first( nps )
implicit none
integer(I8) :: RN_query_first
integer(I8), intent(in) :: nps
RN_query_first = RN_skip_ahead( RN_SEED0, nps*RN_STRIDE )
return
end function RN_query_first
!-------------------------------------------------------------------
subroutine RN_update_stats()
! update overall RN count info
implicit none
!$OMP CRITICAL (RN_STATS)
RN_COUNT_TOTAL = RN_COUNT_TOTAL + RN_COUNT
if( RN_COUNT>RN_COUNT_MAX ) then
RN_COUNT_MAX = RN_COUNT
RN_COUNT_MAX_NPS = RN_NPS
endif
if( RN_COUNT>RN_STRIDE ) then
RN_COUNT_STRIDE = RN_COUNT_STRIDE + 1
endif
!$OMP END CRITICAL (RN_STATS)
RN_COUNT = 0
RN_NPS = 0
return
end subroutine RN_update_stats
!-------------------------------------------------------------------
subroutine expire( c1, c2 )
character(len=*), intent(in) :: c1, c2
write(*,*) ' ********** error: ',c1
write(*,*) ' ********** error: ',c2
stop '**error**'
end subroutine expire
!-------------------------------------------------------------------
!###################################################################
!#
!# Unit tests
!#
!###################################################################
subroutine RN_test_basic( new_gen )
! test routine for basic random number generator
implicit none
integer, intent(in) :: new_gen
real(R8) :: s
integer(I8) :: seeds(10)
integer :: i, j
write(jtty,"(/,a)") " ***** random number - basic test *****"
! set the seed
call RN_init_problem( new_gen, 1_I8, 0_I8, 0_I8, 1 )
! get the first 5 seeds, then skip a few, get 5 more - directly
s = 0.0_R8
do i = 1,5
s = s + rang()
seeds(i) = RN_query( "seed" )
enddo
do i = 6,123455
s = s + rang()
enddo
do i = 6,10
s = s + rang()
seeds(i) = RN_query( "seed" )
enddo
! compare
do i = 1,10
j = i
if( i>5 ) j = i + 123450
write(jtty,"(1x,i6,a,i20,a,i20)") &
& j, " reference: ", RN_CHECK(i,new_gen), " computed: ", seeds(i)
if( seeds(i)/=RN_CHECK(i,new_gen) ) then
write(jtty,"(a)") " ***** basic_test of RN generator failed:"
endif
enddo
return
end subroutine RN_test_basic
!-------------------------------------------------------------------
subroutine RN_test_skip( new_gen )
! test routine for basic random number generation & skip-ahead
implicit none
integer, intent(in) :: new_gen
integer(I8) :: seeds(10)
integer :: i, j
! set the seed
call RN_init_problem( new_gen, 1_I8, 0_I8, 0_I8, 0 )
! use the skip-ahead function to get first 5 seeds, then 5 more
do i = 1,10
j = i
if( i>5 ) j = i + 123450
seeds(i) = RN_skip_ahead( 1_I8, int(j,I8) )
enddo
! compare
write(jtty,"(/,a)") " ***** random number - skip test *****"
do i = 1,10
j = i
if( i>5 ) j = i + 123450
write(jtty,"(1x,i6,a,i20,a,i20)") &
& j, " reference: ", RN_CHECK(i,new_gen), " computed: ", seeds(i)
if( seeds(i)/=RN_CHECK(i,new_gen) ) then
write(jtty,"(a)") " ***** skip_test of RN generator failed:"
endif
enddo
return
end subroutine RN_test_skip
!-------------------------------------------------------------------
subroutine RN_test_mixed( new_gen )
! test routine -- print RN's 1-5 & 123456-123460,
! with reference vals
implicit none
integer, intent(in) :: new_gen
integer(I8) :: r
integer :: i, j
write(jtty,"(/,a)") " ***** random number - mixed test *****"
! set the seed & set the stride to 1
call RN_init_problem( new_gen, 1_I8, 1_I8, 0_I8, 0 )
write(jtty,"(a,i20,z20)") " RN_MULT = ", RN_MULT, RN_MULT
write(jtty,"(a,i20,z20)") " RN_ADD = ", RN_ADD, RN_ADD
write(jtty,"(a,i20,z20)") " RN_MOD = ", RN_MOD, RN_MOD
write(jtty,"(a,i20,z20)") " RN_MASK = ", RN_MASK, RN_MASK
write(jtty,"(a,i20)") " RN_BITS = ", RN_BITS
write(jtty,"(a,i20)") " RN_PERIOD = ", RN_PERIOD
write(jtty,"(a,es20.13)") " RN_NORM = ", RN_NORM
write(jtty,"(a)") " "
do i = 1,10
j = i
if( i>5 ) j = i + 123450
call RN_init_particle( int(j,I8) )
r = RN_query( "seed" )
write(jtty,"(1x,i6,a,i20,a,i20)") &
& j, " reference: ", RN_CHECK(i,new_gen)," computed: ", r
if( r/=RN_CHECK(i,new_gen) ) then
write(jtty,"(a)") " ***** mixed test of RN generator failed:"
endif
enddo
return
end subroutine RN_test_mixed
!-------------------------------------------------------------------
end module mcnp_random

8
src/mpi_routines.f90
View file

@ -3,9 +3,9 @@ module mpi_routines
use constants, only: MAX_LINE_LEN
use error, only: fatal_error
use global
use mcnp_random, only: rang, RN_init_particle, RN_skip
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use random_lcg, only: prn, set_particle_seed, prn_skip
use tally_header, only: TallyObject
#ifdef MPI
@ -162,10 +162,10 @@ contains
end if
! Make sure all processors start at the same point for random sampling
call RN_init_particle(int(i_cycle,8))
call set_particle_seed(int(i_cycle,8))
! Skip ahead however many random numbers are needed
call RN_skip(start)
call prn_skip(start)
allocate(temp_sites(2*work))
count = 0_8 ! Index for local source_bank
@ -198,7 +198,7 @@ contains
end if
! Randomly sample sites needed
if (rang() < p_sample) then
if (prn() < p_sample) then
count = count + 1
temp_sites(count) = fission_bank(i)
end if

128
src/physics.f90
View file

@ -10,12 +10,12 @@ module physics
use geometry_header, only: Universe, BASE_UNIVERSE
use global
use interpolation, only: interpolate_tab1
use mcnp_random, only: rang
use output, only: write_message, print_particle
use particle_header, only: Particle
use tally, only: score_tally, score_surface_current
use random_lcg, only: prn
use search, only: binary_search
use string, only: int_to_str
use tally, only: score_tally, score_surface_current
implicit none
@ -73,7 +73,7 @@ contains
call dist_to_boundary(p, d_to_boundary, surf, in_lattice)
! Sample a distance to collision
d_to_collision = -log(rang()) / material_xs % total
d_to_collision = -log(prn()) / material_xs % total
! Select smaller of the two distances
distance = min(d_to_boundary, d_to_collision)
@ -454,7 +454,7 @@ contains
! ==========================================================================
! SAMPLE NUCLIDE WITHIN THE MATERIAL
cutoff = rang() * material_xs % total
cutoff = prn() * material_xs % total
prob = ZERO
i = 0
@ -491,7 +491,7 @@ contains
else
! set cutoff variable for analog cases
cutoff = rang() * micro_xs(index_nuclide) % total
cutoff = prn() * micro_xs(index_nuclide) % total
prob = ZERO
! Add disappearance cross-section to prob
@ -518,7 +518,7 @@ contains
! created.
if (nuc % has_partial_fission) then
cutoff = rang() * micro_xs(index_nuclide) % fission
cutoff = prn() * micro_xs(index_nuclide) % fission
prob = ZERO
i = 0
@ -586,7 +586,7 @@ contains
if (survival_biasing) then
if (p % wgt < weight_cutoff) then
if (rang() < p % wgt / weight_survive) then
if (prn() < p % wgt / weight_survive) then
p % wgt = weight_survive
else
p % wgt = ZERO
@ -598,7 +598,7 @@ contains
! survival biasing. The cutoff will be a random number times the
! scattering cross section
cutoff = rang() * (micro_xs(index_nuclide) % total - &
cutoff = prn() * (micro_xs(index_nuclide) % total - &
micro_xs(index_nuclide) % absorption)
prob = ZERO
end if
@ -772,7 +772,7 @@ contains
sab => sab_tables(index_sab)
! Determine whether inelastic or elastic scattering will occur
if (rang() < micro_xs(index_nuclide) % elastic_sab / &
if (prn() < micro_xs(index_nuclide) % elastic_sab / &
micro_xs(index_nuclide) % elastic) then
! elastic scattering
@ -793,7 +793,7 @@ contains
! data derived in the incoherent approximation
! Sample outgoing cosine bin
k = 1 + rang() * sab % n_elastic_mu
k = 1 + prn() * sab % n_elastic_mu
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % elastic_mu(k,i)
@ -808,7 +808,7 @@ contains
! edges.
! Sample a Bragg edge between 1 and i
prob = rang() * sab % elastic_P(i+1)
prob = prn() * sab % elastic_P(i+1)
k = binary_search(sab % elastic_P(1:i+1), i+1, prob)
! Characteristic scattering cosine for this Bragg egg
@ -843,9 +843,9 @@ contains
if (sab % secondary_mode == SAB_SECONDARY_EQUAL) then
! All bins equally likely
j = 1 + rang() * n_energy_out
j = 1 + prn() * n_energy_out
elseif (sab % secondary_mode == SAB_SECONDARY_SKEWED) then
r = rang() * (n_energy_out - 3)
r = prn() * (n_energy_out - 3)
if (r > ONE) then
! equally likely N-4 middle bins
j = r + 2
@ -875,7 +875,7 @@ contains
E = (1 - f)*E_ij + f*E_i1j
! Sample outgoing cosine bin
k = 1 + rang() * sab % n_inelastic_mu
k = 1 + prn() * sab % n_inelastic_mu
! Determine outgoing cosine corresponding to E_in(i) and E_in(i+1)
mu_ijk = sab % inelastic_mu(k,j,i)
@ -942,10 +942,10 @@ contains
do
! Sample two random numbers
r1 = rang()
r2 = rang()
r1 = prn()
r2 = prn()
if (rang() < alpha) then
if (prn() < alpha) then
! With probability alpha, we sample the distribution p(y) =
! y*e^(-y). This can be done with sampling scheme C45 frmo the Monte
! Carlo sampler
@ -957,7 +957,7 @@ contains
! e^(-y^2). This can be done with sampling scheme C61 from the Monte
! Carlo sampler
c = cos(PI/2.0 * rang())
c = cos(PI/2.0 * prn())
beta_vt_sq = -log(r1) - log(r2)*c*c
end if
@ -965,14 +965,14 @@ contains
beta_vt = sqrt(beta_vt_sq)
! Sample cosine of angle between neutron and target velocity
mu = 2.0*rang() - ONE
mu = 2.0*prn() - ONE
! Determine rejection probability
accept_prob = sqrt(beta_vn*beta_vn + beta_vt_sq - 2*beta_vn*beta_vt*mu) &
/(beta_vn + beta_vt)
! Perform rejection sampling on vt and mu
if (rang() < accept_prob) exit
if (prn() < accept_prob) exit
end do
! determine direction of target velocity based on the neutron's velocity
@ -1065,7 +1065,7 @@ contains
nu_t = p % last_wgt * micro_xs(index_nuclide) % fission / (keff * &
micro_xs(index_nuclide) % total) * nu_t
end if
if (rang() > nu_t - int(nu_t)) then
if (prn() > nu_t - int(nu_t)) then
nu = int(nu_t)
else
nu = int(nu_t) + 1
@ -1082,12 +1082,12 @@ contains
mu = sample_angle(rxn, E)
! sample between delayed and prompt neutrons
if (rang() < beta) then
if (prn() < beta) then
! ====================================================================
! DELAYED NEUTRON SAMPLED
! sampled delayed precursor group
xi = rang()
xi = prn()
loc = 1
prob = ZERO
do j = 1, nuc % n_precursor
@ -1139,7 +1139,7 @@ contains
end if
! Sample azimuthal angle uniformly in [0,2*pi)
phi = TWO*PI*rang()
phi = TWO*PI*prn()
fission_bank(i) % uvw(1) = mu
fission_bank(i) % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
fission_bank(i) % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)
@ -1262,7 +1262,7 @@ contains
! check if reaction has angular distribution -- if not, sample outgoing
! angle isotropically
if (.not. rxn % has_angle_dist) then
mu = TWO * rang() - ONE
mu = TWO * prn() - ONE
return
end if
@ -1284,16 +1284,16 @@ contains
end if
! Sample between the ith and (i+1)th bin
if (r > rang()) i = i + 1
if (r > prn()) i = i + 1
! check whether this is a 32-equiprobable bin or a tabular distribution
loc = rxn % adist % location(i)
type = rxn % adist % type(i)
if (type == ANGLE_ISOTROPIC) then
mu = TWO * rang() - ONE
mu = TWO * prn() - ONE
elseif (type == ANGLE_32_EQUI) then
! sample cosine bin
xi = rang()
xi = prn()
k = 1 + int(32.0_8*xi)
! calculate cosine
@ -1306,7 +1306,7 @@ contains
NP = rxn % adist % data(loc + 2)
! determine outgoing cosine bin
xi = rang()
xi = prn()
loc = loc + 2
c_k = rxn % adist % data(loc + 2*NP + 1)
do k = 1, NP-1
@ -1376,7 +1376,7 @@ contains
w0 = w
! Sample azimuthal angle in [0,2pi)
phi = TWO * PI * rang()
phi = TWO * PI * prn()
! Precompute factors to save flops
sinphi = sin(phi)
@ -1466,7 +1466,7 @@ contains
if (edist % p_valid % n_regions > 0) then
p_valid = interpolate_tab1(edist % p_valid, E_in)
if (rang() > p_valid) then
if (prn() > p_valid) then
if (edist % law == 44 .or. edist % law == 61) then
call sample_energy(edist%next, E_in, E_out, mu_out)
elseif (edist % law == 66) then
@ -1503,7 +1503,7 @@ contains
& (edist%data(loc+i+1) - edist%data(loc+i))
! Sample outgoing energy bin
r1 = rang()
r1 = prn()
k = 1 + int(NET * r1)
! Determine E_1 and E_K
@ -1520,7 +1520,7 @@ contains
! Randomly select between the outgoing table for incoming energy E_i and
! E_(i+1)
if (rang() < r) then
if (prn() < r) then
l = i + 1
else
l = i
@ -1532,7 +1532,7 @@ contains
E_l_k1 = edist % data(loc+k+1)
! Determine E' (denoted here as E_out)
r2 = rang()
r2 = prn()
E_out = E_l_k + r2*(E_l_k1 - E_l_k)
! Now interpolate between incident energy bins i and i + 1
@ -1580,7 +1580,7 @@ contains
end if
! Sample between the ith and (i+1)th bin
r2 = rang()
r2 = prn()
if (r > r2) then
l = i + 1
else
@ -1623,7 +1623,7 @@ contains
end if
! determine outgoing energy bin
r1 = rang()
r1 = prn()
loc = loc + 2 ! start of EOUT
c_k = edist % data(loc + 2*NP + 1)
do k = 1, NP-1
@ -1697,8 +1697,8 @@ contains
! sample outgoing energy based on evaporation spectrum probability
! density function
do
r1 = rang()
r2 = rang()
r1 = prn()
r2 = prn()
E_out = -T * log(r1*r2)
if (E_out <= E_in - U) exit
end do
@ -1759,7 +1759,7 @@ contains
end if
! Sample between the ith and (i+1)th bin
r2 = rang()
r2 = prn()
if (r > r2) then
l = i + 1
else
@ -1803,7 +1803,7 @@ contains
end if
! determine outgoing energy bin
r1 = rang()
r1 = prn()
loc = loc + 2 ! start of EOUT
c_k = edist % data(loc + 2*NP + 1)
do k = 1, NP-1
@ -1858,8 +1858,8 @@ contains
end if
! Sampled correlated angle from Kalbach-Mann parameters
r3 = rang()
r4 = rang()
r3 = prn()
r4 = prn()
T = (TWO*r4 - ONE) * sinh(KM_A)
if (r3 > KM_R) then
mu_out = log(T + sqrt(T*T + ONE))/KM_A
@ -1902,7 +1902,7 @@ contains
end if
! Sample between the ith and (i+1)th bin
r2 = rang()
r2 = prn()
if (r > r2) then
l = i + 1
else
@ -1946,7 +1946,7 @@ contains
end if
! determine outgoing energy bin
r1 = rang()
r1 = prn()
loc = loc + 2 ! start of EOUT
c_k = edist % data(loc + 2*NP + 1)
do k = 1, NP-1
@ -1992,7 +1992,7 @@ contains
! Check if angular distribution is isotropic
if (loc == 0) then
mu_out = TWO * rang() - ONE
mu_out = TWO * prn() - ONE
return
end if
@ -2001,7 +2001,7 @@ contains
NP = edist % data(loc + 2)
! determine outgoing cosine bin
r3 = rang()
r3 = prn()
loc = loc + 2
c_k = edist % data(loc + 2*NP + 1)
do k = 1, NP-1
@ -2051,17 +2051,17 @@ contains
case (3)
y = maxwell_spectrum(ONE)
case (4)
r1 = rang()
r2 = rang()
r3 = rang()
r1 = prn()
r2 = prn()
r3 = prn()
y = -log(r1*r2*r3)
case (5)
r1 = rang()
r2 = rang()
r3 = rang()
r4 = rang()
r5 = rang()
r6 = rang()
r1 = prn()
r2 = prn()
r3 = prn()
r4 = prn()
r5 = prn()
r6 = prn()
y = -log(r1*r2*r3*r4) - log(r5) * cos(PI/2.*r6)**2
end select
@ -2092,9 +2092,9 @@ contains
real(8) :: r1, r2, r3 ! random numbers
real(8) :: c ! cosine of pi/2*r3
r1 = rang()
r2 = rang()
r3 = rang()
r1 = prn()
r2 = prn()
r3 = prn()
! determine cosine of pi/2*r
c = cos(PI/2.*r3)
@ -2121,7 +2121,7 @@ contains
real(8) :: w ! sampled from Maxwellian
w = maxwell_spectrum(a)
E_out = w + a*a*b/4. + (2.*rang() - ONE)*sqrt(a*a*b*w)
E_out = w + a*a*b/4. + (2.*prn() - ONE)*sqrt(a*a*b*w)
end function watt_spectrum
@ -2137,7 +2137,7 @@ contains
real(8) :: c
c = -4.*D_avg*D_avg/PI * log(rang())
c = -4.*D_avg*D_avg/PI * log(prn())
D = sqrt(c)
end function wigner
@ -2166,7 +2166,7 @@ contains
! sample x as -2/n*log(product(r_i, i = 1 to n/2))
x = ONE
do i = 1, n/2
x = x * rang()
x = x * prn()
end do
x = -2./n * log(x)
@ -2178,11 +2178,11 @@ contains
! Note that we take advantage of integer division on n/2
y = ONE
do i = 1, n/2
y = y * rang()
y = y * prn()
end do
r1 = rang()
r2 = rang()
r1 = prn()
r2 = prn()
c = cos(PI/2.*r2)
x = -2./n * (log(y) + log(r1)*c*c)
end select

147
src/random_lcg.f90 Normal file
View file

@ -0,0 +1,147 @@
module random_lcg
implicit none
private
save
integer(8) :: prn_seed0 ! original seed
integer(8) :: prn_seed ! current seed
integer(8) :: prn_mult ! multiplication factor, g
integer(8) :: prn_add ! additive factor, c
integer :: prn_bits ! number of bits, M
integer(8) :: prn_mod ! 2^M
integer(8) :: prn_mask ! 2^M - 1
integer(8) :: prn_stride ! stride between particles
real(8) :: prn_norm ! 2^(-M)
public :: prn
public :: initialize_prng
public :: set_particle_seed
public :: prn_skip
contains
!===============================================================================
! PRN generates a pseudo-random number using a linean congruential generator
!===============================================================================
function prn() result(pseudo_rn)
real(8) :: pseudo_rn
! This algorithm uses bit-masking to find the next integer(8) value to be
! used to calculate the random number
prn_seed = iand(prn_mult*prn_seed + prn_add, prn_mask)
! Once the integer is calculated, we just need to divide by 2**m,
! represented here as multiplying by a pre-calculated factor
pseudo_rn = prn_seed * prn_norm
end function prn
!===============================================================================
! INITIALIZE_PRNG sets up the random number generator, determining the seed and
! values for g, c, and m.
!===============================================================================
subroutine initialize_prng()
prn_seed0 = 1_8
prn_seed = prn_seed
prn_mult = 2806196910506780709_8
prn_add = 1_8
prn_bits = 63
prn_mod = ibset(0_8, prn_bits) ! clever way of calculating 2**bits
prn_mask = prn_mod - 1_8
prn_stride = 152917_8
prn_norm = 2._8**(-prn_bits)
end subroutine initialize_prng
!===============================================================================
! SET_PARTICLE_SEED sets the seed to a unique value based on the ID of the
! particle
!===============================================================================
subroutine set_particle_seed(id)
integer(8), intent(in) :: id
prn_seed = prn_skip_ahead(id*prn_stride, prn_seed0)
end subroutine set_particle_seed
!===============================================================================
! PRN_SKIP advances the random number seed 'n' times from the current seed
!===============================================================================
subroutine prn_skip(n)
integer(8), intent(in) :: n ! number of seeds to skip
prn_seed = prn_skip_ahead(n, prn_seed)
end subroutine prn_skip
!===============================================================================
! PRN_SKIP_AHEAD advances the random number seed 'skip' times. This is usually
! used to skip a fixed number of random numbers (the stride) so that a given
! particle always has the same starting seed regardless of how many processors
! are used
!===============================================================================
function prn_skip_ahead(n, seed) result(new_seed)
integer(8), intent(in) :: n ! number of seeds to skip
integer(8), intent(in) :: seed ! original seed
integer(8) :: new_seed ! new seed
integer(8) :: nskip
integer(8) :: gen
integer(8) :: g
integer(8) :: inc
integer(8) :: c
integer(8) :: gp
! In cases where we want to skip backwards, we add the period of the random
! number generator until the number of PRNs to skip is positive since
! skipping ahead that much is the same as skipping backwards by the original
! amount
nskip = n
do while (nskip < 0_8)
nskip = nskip + prn_mod
enddo
! The algorithm here to determine the parameters used to skip ahead is
! described in F. Brown, "Random Number Generation with Arbitrary Stride,"
! Trans. Am. Nucl. Soc. (Nov. 1994). This algorithm is able to skip ahead in
! O(log2(N)) operations instead of O(N). Basically, it computes parameters G
! and C which can then be used to find x_N = G*x_0 + C mod 2^M.
nskip = iand(nskip, prn_mask)
gen = 1
g = prn_mult
inc = 0
c = prn_add
do while (nskip > 0_8)
if (btest(nskip,0)) then
gen = iand(gen*g, prn_mask)
inc = iand(inc*g, prn_mask)
inc = iand(inc+c, prn_mask)
endif
gp = iand(g+1, prn_mask)
g = iand(g*g, prn_mask)
c = iand(gp*c, prn_mask)
nskip = ishft(nskip, -1)
enddo
! With G and C, we can now find the new seed
new_seed = iand(gen*seed + inc, prn_mask)
end function prn_skip_ahead
end module random_lcg

10
src/source.f90
View file

@ -4,10 +4,10 @@ module source
use constants, only: ONE, MAX_LINE_LEN
use cross_section_header, only: Nuclide
use global
use mcnp_random, only: rang, RN_init_particle
use output, only: write_message
use particle_header, only: Particle, initialize_particle
use physics, only: watt_spectrum
use random_lcg, only: prn, set_particle_seed
implicit none
@ -61,18 +61,18 @@ contains
p => source_bank(j - bank_first + 1)
! initialize random number seed
call RN_init_particle(int(j,8))
call set_particle_seed(int(j,8))
! sample position
r = (/ (rang(), k = 1,3) /)
r = (/ (prn(), k = 1,3) /)
p % id = j
p % xyz = p_min + r*(p_max - p_min)
p % xyz_local = p % xyz
p % last_xyz = p % xyz
! sample angle
phi = TWO*PI*rang()
mu = TWO*rang() - ONE
phi = TWO*PI*prn()
mu = TWO*prn() - ONE
p % uvw(1) = mu
p % uvw(2) = sqrt(ONE - mu*mu) * cos(phi)
p % uvw(3) = sqrt(ONE - mu*mu) * sin(phi)