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