MPI-AMRVAC  3.1
The MPI - Adaptive Mesh Refinement - Versatile Advection Code (development version)
mod_random.t
Go to the documentation of this file.
1 !> Module for pseudo random number generation. The internal pseudo random
2 !> generator is the xoroshiro128plus method.
3 module mod_random
4 
5 #include "amrvac.h"
6 
7  implicit none
8  private
9 
10  ! A 64 bit floating point type
11  integer, parameter :: dp = kind(0.0d0)
12 
13  ! A 32 bit integer type
14  integer, parameter :: i4 = selected_int_kind(9)
15 
16  ! A 64 bit integer type
17  integer, parameter :: i8 = selected_int_kind(18)
18 
19  !> Random number generator type, which contains the state
20  type rng_t
21  !> The rng state (always use your own seed)
22  integer(i8), private :: s(2) = [123456789_i8, 987654321_i8]
23  contains
24  procedure, non_overridable :: set_seed ! Seed the generator
25  procedure, non_overridable :: jump ! Jump function (see below)
26  procedure, non_overridable :: int_4 ! 4-byte random integer
27  procedure, non_overridable :: int_8 ! 8-byte random integer
28  procedure, non_overridable :: unif_01 ! Uniform (0,1] real
29  procedure, non_overridable :: unif_01_vec ! Uniform (0,1] real vector
30  procedure, non_overridable :: normal ! One normal(0,1) number
31  procedure, non_overridable :: two_normals ! Two normal(0,1) samples
32  procedure, non_overridable :: poisson ! Sample from Poisson-dist.
33  procedure, non_overridable :: circle ! Sample on a circle
34  procedure, non_overridable :: sphere ! Sample on a sphere
35  procedure, non_overridable :: next ! Internal method
36  end type rng_t
37 
38  !> Parallel random number generator type
39  type prng_t
40  type(rng_t), allocatable :: rngs(:)
41  contains
42  procedure, non_overridable :: init_parallel
43  end type prng_t
44 
45  public :: rng_t
46  public :: prng_t
47 
48 contains
49 
50 
51 #if defined(__NVCOMPILER) || (defined(USE_INTRINSIC_SHIFT) && USE_INTRINSIC_SHIFT==0)
52 
53  !> added for nvidia compilers
54  function shiftl(val, shift) result(res_value)
55  integer(i8), intent(in) :: val
56  integer, intent(in) :: shift
57  integer(i8) :: res_value
58 ! integer(i8),parameter :: bit_mask1=9223372036854775807
59 ! !b'0111111111111111111111111111111111111111111111111111111111111111'
60 ! integer(i8),parameter :: bit_mask2=-9223372036854775808
61 ! !b'1000000000000000000000000000000000000000000000000000000000000000'
62  integer(i8) :: bit_mask1, bit_mask2
63 
64  ! cannot be initialized with b values in gnu, cannot have the big decimal numbers in nvidia
65  bit_mask1=huge(val)
66  bit_mask2=-bit_mask1-1
67 
68 
69  if(val<0) then
70  res_value = ior(lshift(iand(val, bit_mask1), shift),bit_mask2)
71  else
72  res_value = lshift(val, shift)
73  endif
74 
75  end function shiftl
76 
77 
78  function shiftr(val, shift) result(res_value)
79  integer(i8), intent(in) :: val
80  integer, intent(in) :: shift
81  integer(i8) :: res_value
82 ! integer(i8),parameter :: bit_mask1=9223372036854775807
83 ! !b'0111111111111111111111111111111111111111111111111111111111111111'
84 ! integer(i8),parameter :: bit_mask2=-9223372036854775808
85 ! !b'1000000000000000000000000000000000000000000000000000000000000000'
86  integer(i8) :: bit_mask1, bit_mask2
87 
88  ! cannot be initialized with b values in gnu, cannot have the big decimal numbers in nvidia
89  bit_mask1=huge(val)
90  bit_mask2=-bit_mask1-1
91 
92  if(val<0) then
93  res_value = ior(rshift(iand(val, bit_mask1), shift),bit_mask2)
94  else
95  res_value = rshift(val, shift)
96  endif
97 
98  end function shiftr
99 
100 #endif
101 
102 
103  !> Initialize a collection of rng's for parallel use
104  subroutine init_parallel(self, n_proc, rng)
105  class(prng_t), intent(inout) :: self
106  type(rng_t), intent(inout) :: rng
107  integer, intent(in) :: n_proc
108  integer :: n
109 
110  allocate(self%rngs(n_proc))
111  self%rngs(1) = rng
112 
113  do n = 2, n_proc
114  self%rngs(n) = self%rngs(n-1)
115  call self%rngs(n)%jump()
116  end do
117  end subroutine init_parallel
118 
119  !> Set a seed for the rng
120  subroutine set_seed(self, the_seed)
121  class(rng_t), intent(inout) :: self
122  integer(i8), intent(in) :: the_seed(2)
123 
124  self%s = the_seed
125 
126  ! Simulate calls to next() to improve randomness of first number
127  call self%jump()
128  end subroutine set_seed
129 
130  ! This is the jump function for the generator. It is equivalent
131  ! to 2^64 calls to next(); it can be used to generate 2^64
132  ! non-overlapping subsequences for parallel computations.
133  subroutine jump(self)
134  class(rng_t), intent(inout) :: self
135  integer :: i, b
136  integer(i8) :: global_time(2), dummy
137 
138  ! The signed equivalent of the unsigned constants
139  integer(i8), parameter :: jmp_c(2) = &
140  (/-4707382666127344949_i8, -2852180941702784734_i8/)
141 
142  global_time = 0
143  do i = 1, 2
144  do b = 0, 63
145  if (iand(jmp_c(i), shiftl(1_i8, b)) /= 0) then
146  global_time = ieor(global_time, self%s)
147  end if
148  dummy = self%next()
149  end do
150  end do
151 
152  self%s = global_time
153  end subroutine jump
154 
155  !> Return 4-byte integer
156  integer(i4) function int_4(self)
157  class(rng_t), intent(inout) :: self
158  int_4 = int(self%next(), i4)
159  end function int_4
160 
161  !> Return 8-byte integer
162  integer(i8) function int_8(self)
163  class(rng_t), intent(inout) :: self
164  int_8 = self%next()
165  end function int_8
166 
167  !> Get a uniform [0,1) random real (double precision)
168  real(dp) function unif_01(self)
169  class(rng_t), intent(inout) :: self
170  integer(i8) :: x
171  real(dp) :: tmp
172 
173  x = self%next()
174  x = ior(shiftl(1023_i8, 52), shiftr(x, 12))
175  unif_01 = transfer(x, tmp) - 1.0_dp
176  end function unif_01
177 
178  !> Fill array with uniform random numbers
179  subroutine unif_01_vec(self, rr)
180  class(rng_t), intent(inout) :: self
181  real(dp), intent(out) :: rr(:)
182  integer :: i
183 
184  do i = 1, size(rr)
185  rr(i) = self%unif_01()
186  end do
187  end subroutine unif_01_vec
188 
189  !> Return two normal random variates with mean 0 and variance 1.
190  !> http://en.wikipedia.org/wiki/Marsaglia_polar_method
191  function two_normals(self) result(rands)
192  class(rng_t), intent(inout) :: self
193  real(dp) :: rands(2), sum_sq
194 
195  do
196  rands(1) = 2 * self%unif_01() - 1
197  rands(2) = 2 * self%unif_01() - 1
198  sum_sq = sum(rands**2)
199  if (sum_sq < 1.0_dp .and. sum_sq > 0.0_dp) exit
200  end do
201  rands = rands * sqrt(-2 * log(sum_sq) / sum_sq)
202  end function two_normals
203 
204  !> Single normal random number
205  real(dp) function normal(self)
206  class(rng_t), intent(inout) :: self
207  real(dp) :: rands(2)
208 
209  rands = self%two_normals()
210  normal = rands(1)
211  end function normal
212 
213  !> Return Poisson random variate with rate lambda. Works well for lambda < 30
214  !> or so. For lambda >> 1 it can produce wrong results due to roundoff error.
215  function poisson(self, lambda) result(rr)
216  class(rng_t), intent(inout) :: self
217  real(dp), intent(in) :: lambda
218  integer(i4) :: rr
219  real(dp) :: expl, p
220 
221  expl = exp(-lambda)
222  rr = 0
223  p = self%unif_01()
224 
225  do while (p > expl)
226  rr = rr + 1
227  p = p * self%unif_01()
228  end do
229  end function poisson
230 
231  !> Sample point on a circle with given radius
232  function circle(self, radius) result(xy)
233  class(rng_t), intent(inout) :: self
234  real(dp), intent(in) :: radius
235  real(dp) :: rands(2), xy(2)
236  real(dp) :: sum_sq
237 
238  ! Method for uniform sampling on circle
239  do
240  rands(1) = 2 * self%unif_01() - 1
241  rands(2) = 2 * self%unif_01() - 1
242  sum_sq = sum(rands**2)
243  if (sum_sq <= 1) exit
244  end do
245 
246  xy(1) = (rands(1)**2 - rands(2)**2) / sum_sq
247  xy(2) = 2 * rands(1) * rands(2) / sum_sq
248  xy = xy * radius
249  end function circle
250 
251  !> Sample point on a sphere with given radius
252  function sphere(self, radius) result(xyz)
253  class(rng_t), intent(inout) :: self
254  real(dp), intent(in) :: radius
255  real(dp) :: rands(2), xyz(3)
256  real(dp) :: sum_sq, tmp_sqrt
257 
258  ! Marsaglia method for uniform sampling on sphere
259  do
260  rands(1) = 2 * self%unif_01() - 1
261  rands(2) = 2 * self%unif_01() - 1
262  sum_sq = sum(rands**2)
263  if (sum_sq <= 1) exit
264  end do
265 
266  tmp_sqrt = sqrt(1 - sum_sq)
267  xyz(1:2) = 2 * rands(1:2) * tmp_sqrt
268  xyz(3) = 1 - 2 * sum_sq
269  xyz = xyz * radius
270  end function sphere
271 
272  !> Interal routine: get the next value (returned as 64 bit signed integer)
273  function next(self) result(res)
274  class(rng_t), intent(inout) :: self
275  integer(i8) :: res
276  integer(i8) :: global_time(2)
277 
278  global_time = self%s
279  res = global_time(1) + global_time(2)
280  global_time(2) = ieor(global_time(1), global_time(2))
281  self%s(1) = ieor(ieor(rotl(global_time(1), 55), global_time(2)), shiftl(global_time(2), 14))
282  self%s(2) = rotl(global_time(2), 36)
283  end function next
284 
285  !> Helper function for next()
286  pure function rotl(x, k) result(res)
287  integer(i8), intent(in) :: x
288  integer, intent(in) :: k
289  integer(i8) :: res
290 
291  res = ior(shiftl(x, k), shiftr(x, 64 - k))
292  end function rotl
293 
294 end module mod_random
Module for pseudo random number generation. The internal pseudo random generator is the xoroshiro128p...
Definition: mod_random.t:3