diff --git a/src/DEPENDENCIES b/src/DEPENDENCIES index 83a718fef7..9e526688a4 100644 --- a/src/DEPENDENCIES +++ b/src/DEPENDENCIES @@ -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 diff --git a/src/OBJECTS b/src/OBJECTS index a80d124184..6ce0dcafc7 100644 --- a/src/OBJECTS +++ b/src/OBJECTS @@ -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 \ diff --git a/src/initialize.f90 b/src/initialize.f90 index 15fa8bdc85..8c84c2d7b6 100644 --- a/src/initialize.f90 +++ b/src/initialize.f90 @@ -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() diff --git a/src/main.f90 b/src/main.f90 index 2dfdcce773..de1c473fc3 100644 --- a/src/main.f90 +++ b/src/main.f90 @@ -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) diff --git a/src/mcnp_random.f90 b/src/mcnp_random.f90 deleted file mode 100644 index 1f7e2550ef..0000000000 --- a/src/mcnp_random.f90 +++ /dev/null @@ -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 diff --git a/src/mpi_routines.f90 b/src/mpi_routines.f90 index dc13f04f7f..d1cdd82bc0 100644 --- a/src/mpi_routines.f90 +++ b/src/mpi_routines.f90 @@ -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 diff --git a/src/physics.f90 b/src/physics.f90 index 0cf2f41596..4e896bfb92 100644 --- a/src/physics.f90 +++ b/src/physics.f90 @@ -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 diff --git a/src/random_lcg.f90 b/src/random_lcg.f90 new file mode 100644 index 0000000000..8b7225c7ff --- /dev/null +++ b/src/random_lcg.f90 @@ -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 diff --git a/src/source.f90 b/src/source.f90 index 1da0ea2e68..30e60da351 100644 --- a/src/source.f90 +++ b/src/source.f90 @@ -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)