Merge pull request #235 from mit-crpg/rng-streams

Multiple random number generator streams

src/random_lcg.F90

@@ -5,8 +5,13 @@ module random_lcg
   private
   save
 
+  ! Random number streams
+  integer, parameter :: N_STREAMS = 2
+  integer, parameter :: STREAM_TRACKING = 1
+  integer, parameter :: STREAM_TALLIES  = 2
+
   integer(8) :: prn_seed0  ! original seed
-  integer(8) :: prn_seed   ! current seed
+  integer(8) :: prn_seed(N_STREAMS) ! current seed
   integer(8) :: prn_mult   ! multiplication factor, g
   integer(8) :: prn_add    ! additive factor, c
   integer    :: prn_bits   ! number of bits, M
@@ -14,6 +19,7 @@ module random_lcg
   integer(8) :: prn_mask   ! 2^M - 1
   integer(8) :: prn_stride ! stride between particles
   real(8)    :: prn_norm   ! 2^(-M)
+  integer    :: stream     ! current RNG stream
 
 !$omp threadprivate(prn_seed)
 
@@ -21,6 +27,8 @@ module random_lcg
   public :: initialize_prng
   public :: set_particle_seed
   public :: prn_skip
+  public :: prn_set_stream
+  public :: STREAM_TRACKING, STREAM_TALLIES
 
 contains
 
@@ -35,12 +43,12 @@ contains
     ! 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)
+    prn_seed(stream) = iand(prn_mult*prn_seed(stream) + 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
+    pseudo_rn = prn_seed(stream) * prn_norm
 
   end function prn
 
@@ -53,8 +61,13 @@ contains
 
     use global, only: seed
 
+    integer :: i
+
+    stream     = STREAM_TRACKING
     prn_seed0  = seed
-    prn_seed   = prn_seed
+    do i = 1, N_STREAMS
+      prn_seed(i) = prn_seed0 + i - 1
+    end do
     prn_mult   = 2806196910506780709_8
     prn_add    = 1_8
     prn_bits   = 63
@@ -74,7 +87,11 @@ contains
 
     integer(8), intent(in) :: id
 
-    prn_seed = prn_skip_ahead(id*prn_stride, prn_seed0)
+    integer :: i
+
+    do i = 1, N_STREAMS
+      prn_seed(i) = prn_skip_ahead(id*prn_stride, prn_seed0 + i - 1)
+    end do
 
   end subroutine set_particle_seed
     
@@ -86,7 +103,7 @@ contains
 
     integer(8), intent(in) :: n ! number of seeds to skip
 
-    prn_seed = prn_skip_ahead(n, prn_seed)
+    prn_seed(stream) = prn_skip_ahead(n, prn_seed(stream))
 
   end subroutine prn_skip
 
@@ -151,4 +168,18 @@ contains
 
   end function prn_skip_ahead
 
+!===============================================================================
+! PRN_SET_STREAM changes the random number stream. If random numbers are needed
+! in routines not used directly for tracking (e.g. physics), this allows the
+! numbers to be generated without affecting reproducibility of the physics.
+!===============================================================================
+
+  subroutine prn_set_stream(i)
+
+    integer, intent(in) :: i
+
+    stream = i
+
+  end subroutine prn_set_stream
+
 end module random_lcg