Home

Context Navigation

source: palm/trunk/SOURCE/random_generator_parallel_mod.f90 @ 4516

Last change on this file since 4516 was 4438, checked in by suehring, 5 years ago
Synthetic turbulence: performance optimizations - random numbers only defined and computed locally, option to compute velocity seeds locally without need of global communication; paralell random number generator: new routine to initialize 1D random number arrays; virtual measurements: CPU-log points added
Property svn:keywords set to `Id`
File size: 20.4 KB

Line
1	!> @file random_generator_parallel_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of the PALM model system.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2020 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: random_generator_parallel_mod.f90 4438 2020-03-03 20:49:28Z suehring $
27	! - Rename variables to avoid confusion with subdomain grid indices
28	! - Some formatting adjustments
29	! - New routine to initialize spatial 1D-arrays independent on the 2D random-
30	! number array
31	!
32	! 4360 2020-01-07 11:25:50Z suehring
33	! Corrected "Former revisions" section
34	!
35	! 3655 2019-01-07 16:51:22Z knoop
36	! unused variable removed
37	!
38	! 1400 2014-05-09 14:03:54Z knoop
39	! Initial revision
40	!
41	!
42	! Description:
43	! ------------
44	!> This module contains and supports the random number generating routine ran_parallel.
45	!> ran_parallel returns a uniform random deviate between 0.0 and 1.0
46	!> (exclusive of the end point values).
47	!> Additionally it provides the generator with five integer for use as initial state space.
48	!> The first tree integers (iran, jran, kran) are maintained as non negative values,
49	!> while the last two (mran, nran) have 32-bit nonzero values.
50	!> Also provided by this module is support for initializing or reinitializing
51	!> the state space to a desired standard sequence number, hashing the initial
52	!> values to random values, and allocating and deallocating the internal workspace
53	!> Random number generator, produces numbers equally distributed in interval
54	!>
55	!> This routine is taken from the "numerical recipies vol. 2"
56	!------------------------------------------------------------------------------!
57	MODULE random_generator_parallel
58
59	USE kinds
60
61	IMPLICIT NONE
62
63	INTEGER(isp), SAVE :: lenran=0 !<
64	INTEGER(isp), SAVE :: seq=0 !<
65	INTEGER(isp), SAVE :: iran0 !<
66	INTEGER(isp), SAVE :: jran0 !<
67	INTEGER(isp), SAVE :: kran0 !<
68	INTEGER(isp), SAVE :: mran0 !<
69	INTEGER(isp), SAVE :: nran0 !<
70	INTEGER(isp), SAVE :: rans !<
71
72	INTEGER(isp), DIMENSION(:, :), POINTER, SAVE :: ranseeds !<
73
74	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: iran !<
75	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: jran !<
76	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: kran !<
77	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: mran !<
78	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: nran !<
79	INTEGER(isp), DIMENSION(:), POINTER, SAVE :: ranv !<
80
81	INTEGER(isp), DIMENSION(:,:), ALLOCATABLE :: id_random_array !<
82	INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE :: seq_random_array !<
83
84	REAL(wp), SAVE :: amm !<
85
86	REAL(wp) :: random_dummy=0.0 !<
87
88	INTERFACE init_parallel_random_generator
89	MODULE PROCEDURE init_parallel_random_generator_1d
90	MODULE PROCEDURE init_parallel_random_generator_2d
91	END INTERFACE
92
93	INTERFACE random_number_parallel
94	MODULE PROCEDURE ran0_s
95	END INTERFACE
96
97	INTERFACE random_seed_parallel
98	MODULE PROCEDURE random_seed_parallel
99	END INTERFACE
100
101	INTERFACE ran_hash
102	MODULE PROCEDURE ran_hash_v
103	END INTERFACE
104
105	INTERFACE reallocate
106	MODULE PROCEDURE reallocate_iv,reallocate_im
107	END INTERFACE
108
109	INTERFACE arth
110	MODULE PROCEDURE arth_i
111	END INTERFACE
112
113	PRIVATE
114
115	PUBLIC random_number_parallel, random_seed_parallel, random_dummy, &
116	id_random_array, seq_random_array, init_parallel_random_generator
117
118	CONTAINS
119
120	!------------------------------------------------------------------------------!
121	! Description:
122	! ------------
123	!> Initialize the parallel random number generator for a 1-dimensional array.
124	!------------------------------------------------------------------------------!
125	SUBROUTINE init_parallel_random_generator_1d( nxy, ns, ne, id_rand, seq_rand )
126
127	USE control_parameters, &
128	ONLY: ensemble_member_nr
129
130	INTEGER(isp), INTENT(IN) :: nxy !< constant scaling with grid dimensions
131	INTEGER(isp), INTENT(IN) :: ns !< start index on subdomain
132	INTEGER(isp), INTENT(IN) :: ne !< end index on subdomain
133
134	INTEGER(iwp) :: j !< loop index
135
136	INTEGER(isp), DIMENSION(ns:ne) :: id_rand !< initial IDs
137	INTEGER(isp), DIMENSION(5,ns:ne) :: seq_rand !< initial random seeds
138
139	!
140	!-- Asigning an ID to every vertical gridpoint column
141	!-- dependig on the ensemble run number.
142	DO j = ns, ne
143	id_rand(j) = j * ( nxy + 1.0_wp ) + 1.0_wp + 1E6 * &
144	( ensemble_member_nr - 1000 )
145	ENDDO
146	!
147	!-- Initializing with random_seed_parallel for every vertical
148	!-- gridpoint column.
149	random_dummy=0
150	DO j = ns, ne
151	CALL random_seed_parallel( random_sequence=id_rand(j) )
152	CALL random_number_parallel( random_dummy )
153	CALL random_seed_parallel( get=seq_rand(:,j) )
154	ENDDO
155
156	END SUBROUTINE init_parallel_random_generator_1d
157
158	!------------------------------------------------------------------------------!
159	! Description:
160	! ------------
161	!> Initialize the parallel random number generator for a specific subdomain
162	!------------------------------------------------------------------------------!
163	SUBROUTINE init_parallel_random_generator_2d( nx_l, nys_l, nyn_l, nxl_l, nxr_l )
164
165	USE kinds
166
167	USE control_parameters, &
168	ONLY: ensemble_member_nr
169
170	IMPLICIT NONE
171
172	INTEGER(isp), INTENT(IN) :: nx_l !< constant
173	INTEGER(isp), INTENT(IN) :: nys_l !< local lower subdomain bound index in y-direction
174	INTEGER(isp), INTENT(IN) :: nyn_l !< local upper subdomain bound index in y-direction
175	INTEGER(isp), INTENT(IN) :: nxl_l !< local lower subdomain bound index in x-direction
176	INTEGER(isp), INTENT(IN) :: nxr_l !< local upper subdomain bound index in x-direction
177
178	INTEGER(iwp) :: i !< grid index x-direction
179	INTEGER(iwp) :: j !< grid index y-direction
180
181	!-- Allocate ID-array and state-space-array
182	ALLOCATE ( seq_random_array(5,nys_l:nyn_l,nxl_l:nxr_l) )
183	ALLOCATE ( id_random_array(nys_l:nyn_l,nxl_l:nxr_l) )
184	seq_random_array = 0
185	id_random_array = 0
186
187	!-- Asigning an ID to every vertical gridpoint column
188	!-- dependig on the ensemble run number.
189	DO i = nxl_l, nxr_l
190	DO j = nys_l, nyn_l
191	id_random_array(j,i) = j * ( nx_l + 1.0_wp ) + i + 1.0_wp + 1E6 * &
192	( ensemble_member_nr - 1000 )
193	ENDDO
194	ENDDO
195	!-- Initializing with random_seed_parallel for every vertical
196	!-- gridpoint column.
197	random_dummy=0
198	DO i = nxl_l, nxr_l
199	DO j = nys_l, nyn_l
200	CALL random_seed_parallel (random_sequence=id_random_array(j, i))
201	CALL random_number_parallel (random_dummy)
202	CALL random_seed_parallel (get=seq_random_array(:, j, i))
203	ENDDO
204	ENDDO
205
206	END SUBROUTINE init_parallel_random_generator_2d
207
208	!------------------------------------------------------------------------------!
209	! Description:
210	! ------------
211	!> Lagged Fibonacci generator combined with a Marsaglia shiftsequence.
212	!> Returns as harvest a uniform random deviate between 0.0 and 1.0 (exclusive of the end point values).
213	!> This generator has the same calling and initialization conventions as Fortran 90's random_number routine.
214	!> Use random_seed_parallel to initialize or reinitialize to a particular sequence.
215	!> The period of this generator is about 2.0 x 10^28, and it fully vectorizes.
216	!> Validity of the integer model assumed by this generator is tested at initialization.
217	!------------------------------------------------------------------------------!
218	SUBROUTINE ran0_s(harvest)
219
220	USE kinds
221
222	IMPLICIT NONE
223
224	REAL(wp), INTENT(OUT) :: harvest !<
225
226	IF (lenran < 1) CALL ran_init(1) !- Initialization routine in ran_state.
227
228	!- Update Fibonacci generator, which has period p^2 + p + 1, p = 2^31 - 69.
229	rans = iran0 - kran0
230
231	IF (rans < 0) rans = rans + 2147483579_isp
232
233	iran0 = jran0
234	jran0 = kran0
235	kran0 = rans
236
237	nran0 = ieor( nran0, ishft (nran0, 13) ) !- Update Marsaglia shift sequence with period 2^32 - 1.
238	nran0 = ieor( nran0, ishft (nran0, -17) )
239	nran0 = ieor( nran0, ishft (nran0, 5) )
240
241	rans = ieor( nran0, rans ) !- Combine the generators.
242
243	harvest = amm * merge( rans, not(rans), rans < 0 ) !- Make the result positive definite (note that amm is negative).
244
245	END SUBROUTINE ran0_s
246
247	!------------------------------------------------------------------------------!
248	! Description:
249	! ------------
250	!> Initialize or reinitialize the random generator state space to vectors of size length.
251	!> The saved variable seq is hashed (via calls to the module routine ran_hash)
252	!> to create unique starting seeds, different for each vector component.
253	!------------------------------------------------------------------------------!
254	SUBROUTINE ran_init( length )
255
256	USE kinds
257
258	IMPLICIT NONE
259
260	INTEGER(isp), INTENT(IN) :: length !<
261
262	INTEGER(isp) :: hg !<
263	INTEGER(isp) :: hgm !<
264	INTEGER(isp) :: hgng !<
265
266	INTEGER(isp) :: new !<
267	INTEGER(isp) :: j !<
268	INTEGER(isp) :: hgt !<
269
270	IF ( length < lenran ) RETURN !- Simply return if enough space is already allocated.
271
272	hg=huge(1_isp)
273	hgm=-hg
274	hgng=hgm-1
275	hgt = hg
276
277	!- The following lines check that kind value isp is in fact a 32-bit integer
278	!- with the usual properties that we expect it to have (under negation and wrap-around addition).
279	!- If all of these tests are satisfied, then the routines that use this module are portable,
280	!- even though they go beyond Fortran 90's integer model.
281
282	IF ( hg /= 2147483647 ) CALL ran_error('arithmetic assumption 1 failed')
283	IF ( hgng >= 0 ) CALL ran_error('arithmetic assumption 2 failed')
284	IF ( hgt+1 /= hgng ) CALL ran_error('arithmetic assumption 3 failed')
285	IF ( not(hg) >= 0 ) CALL ran_error('arithmetic assumption 4 failed')
286	IF ( not(hgng) < 0 ) CALL ran_error('arithmetic assumption 5 failed')
287	IF ( hg+hgng >= 0 ) CALL ran_error('arithmetic assumption 6 failed')
288	IF ( not(-1_isp) < 0 ) CALL ran_error('arithmetic assumption 7 failed')
289	IF ( not(0_isp) >= 0 ) CALL ran_error('arithmetic assumption 8 failed')
290	IF ( not(1_isp) >= 0 ) CALL ran_error('arithmetic assumption 9 failed')
291
292	IF ( lenran > 0) THEN !- Reallocate space, or ...
293
294	ranseeds => reallocate( ranseeds, length, 5)
295	ranv => reallocate( ranv, length-1)
296	new = lenran+1
297
298	ELSE !- allocate space.
299
300	ALLOCATE(ranseeds(length,5))
301	ALLOCATE(ranv(length-1))
302	new = 1 !- Index of first location not yet initialized.
303	amm = nearest(1.0_wp,-1.0_wp)/hgng
304	!- Use of nearest is to ensure that returned random deviates are strictly lessthan 1.0.
305	IF (ammhgng >= 1.0 .or. ammhgng <= 0.0) &
306	CALL ran_error('arithmetic assumption 10 failed')
307
308	END IF
309
310	!- Set starting values, unique by seq and vector component.
311	ranseeds(new:,1) = seq
312	ranseeds(new:,2:5)=spread(arth(new,1,size(ranseeds(new:,1))),2,4)
313
314	DO j=1,4 !- Hash them.
315	CALL ran_hash(ranseeds(new:,j), ranseeds(new:,j+1))
316	END DO
317
318	WHERE (ranseeds (new: ,1:3) < 0) &
319	ranseeds(new: ,1:3)=not(ranseeds(new: ,1:3)) !- Enforce nonnegativity.
320
321	WHERE (ranseeds(new: ,4:5) == 0) ranseeds(new: ,4:5)=1 !- Enforce nonzero.
322
323	IF (new == 1) THEN !- Set scalar seeds.
324
325	iran0 = ranseeds(1,1)
326	jran0 = ranseeds(1,2)
327	kran0 = ranseeds(1,3)
328	mran0 = ranseeds(1,4)
329	nran0 = ranseeds(1,5)
330	rans = nran0
331
332	END IF
333
334	IF (length > 1) THEN !- Point to vector seeds.
335
336	iran => ranseeds(2:,1)
337	jran => ranseeds(2:,2)
338	kran => ranseeds(2:,3)
339	mran => ranseeds(2:,4)
340	nran => ranseeds(2:,5)
341	ranv = nran
342
343	END IF
344
345	lenran = length
346
347	END SUBROUTINE ran_init
348
349	!------------------------------------------------------------------------------!
350	! Description:
351	! ------------
352	!> User interface to release the workspace used by the random number routines.
353	!------------------------------------------------------------------------------!
354	SUBROUTINE ran_deallocate
355
356	IF ( lenran > 0 ) THEN
357
358	DEALLOCATE(ranseeds, ranv)
359	NULLIFY(ranseeds, ranv, iran, jran, kran, mran, nran)
360	lenran = 0
361
362	END IF
363
364	END SUBROUTINE ran_deallocate
365
366	!------------------------------------------------------------------------------!
367	! Description:
368	! ------------
369	!> User interface for seeding the random number routines.
370	!> Syntax is exactly like Fortran 90's random_seed routine,
371	!> with one additional argument keyword: random_sequence, set to any integer
372	!> value, causes an immediate new initialization, seeded by that integer.
373	!------------------------------------------------------------------------------!
374	SUBROUTINE random_seed_parallel( random_sequence, state_size, put, get )
375
376	IMPLICIT NONE
377
378	INTEGER(isp), OPTIONAL, INTENT(IN) :: random_sequence !<
379	INTEGER(isp), OPTIONAL, INTENT(OUT) :: state_size !<
380
381	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(IN) :: put !<
382	INTEGER(isp), DIMENSION(:), OPTIONAL, INTENT(OUT) :: get !<
383
384	IF ( present(state_size) ) THEN
385
386	state_size = 5 * lenran
387
388	ELSE IF ( present(put) ) THEN
389
390	IF ( lenran == 0 ) RETURN
391
392	ranseeds = reshape( put,shape(ranseeds) )
393
394	WHERE (ranseeds(:,1:3) < 0) ranseeds(: ,1:3) = not(ranseeds(: ,1:3))
395	!- Enforce nonnegativity and nonzero conditions on any user-supplied seeds.
396
397	WHERE (ranseeds(:,4:5) == 0) ranseeds(:,4:5) = 1
398
399	iran0 = ranseeds(1,1)
400	jran0 = ranseeds(1,2)
401	kran0 = ranseeds(1,3)
402	mran0 = ranseeds(1,4)
403	nran0 = ranseeds(1,5)
404
405	ELSE IF ( present(get) ) THEN
406
407	IF (lenran == 0) RETURN
408
409	ranseeds(1,1:5) = (/ iran0,jran0,kran0,mran0,nran0 /)
410	get = reshape( ranseeds, shape(get) )
411
412	ELSE IF ( present(random_sequence) ) THEN
413
414	CALL ran_deallocate
415	seq = random_sequence
416
417	END IF
418
419	END SUBROUTINE random_seed_parallel
420
421	!------------------------------------------------------------------------------!
422	! Description:
423	! ------------
424	!> DES-like hashing of two 32-bit integers, using shifts,
425	!> xor's, and adds to make the internal nonlinear function.
426	!------------------------------------------------------------------------------!
427	SUBROUTINE ran_hash_v( il, ir )
428
429	IMPLICIT NONE
430
431	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: il !<
432	INTEGER(isp), DIMENSION(:), INTENT(INOUT) :: ir !<
433
434	INTEGER(isp), DIMENSION(size(il)) :: is !<
435
436	INTEGER(isp) :: j !<
437
438	DO j=1,4
439
440	is = ir
441	ir = ieor( ir, ishft(ir,5) ) + 1422217823
442	ir = ieor( ir, ishft(ir,-16) ) + 1842055030
443	ir = ieor( ir, ishft(ir,9) ) + 80567781
444	ir = ieor( il, ir )
445	il = is
446	END DO
447
448	END SUBROUTINE ran_hash_v
449
450	!------------------------------------------------------------------------------!
451	! Description:
452	! ------------
453	!> User interface to process error-messages
454	!> produced by the random_number_parallel module
455	!------------------------------------------------------------------------------!
456	SUBROUTINE ran_error(string)
457
458	USE control_parameters, &
459	ONLY: message_string
460
461	CHARACTER(LEN=*), INTENT(IN) :: string !< Error message string
462
463	message_string = 'incompatible integer arithmetic: ' // TRIM( string )
464	CALL message( 'ran_init', 'PA0453', 1, 2, 0, 6, 0 )
465
466	END SUBROUTINE ran_error
467
468	!------------------------------------------------------------------------------!
469	! Description:
470	! ------------
471	!> Reallocates the generators state space "ranseeds" to vectors of size length.
472	!------------------------------------------------------------------------------!
473	FUNCTION reallocate_iv( p, n )
474
475	USE control_parameters, &
476	ONLY: message_string
477
478	INTEGER(isp), DIMENSION(:), POINTER :: p !<
479	INTEGER(isp), DIMENSION(:), POINTER :: reallocate_iv !<
480
481	INTEGER(isp), INTENT(IN) :: n !<
482
483	INTEGER(isp) :: nold !<
484	INTEGER(isp) :: ierr !<
485
486	ALLOCATE(reallocate_iv(n),stat=ierr)
487
488	IF (ierr /= 0) THEN
489	message_string = 'problem in attempt to allocate memory'
490	CALL message( 'reallocate_iv', 'PA0454', 1, 2, 0, 6, 0 )
491	END IF
492
493	IF (.not. associated(p)) RETURN
494
495	nold = size(p)
496
497	reallocate_iv(1:min(nold,n)) = p(1:min(nold,n))
498
499	DEALLOCATE(p)
500
501	END FUNCTION reallocate_iv
502
503	FUNCTION reallocate_im( p, n, m )
504
505	USE control_parameters, &
506	ONLY: message_string
507
508	INTEGER(isp), DIMENSION(:,:), POINTER :: p !<
509	INTEGER(isp), DIMENSION(:,:), POINTER :: reallocate_im !<
510
511	INTEGER(isp), INTENT(IN) :: m !<
512	INTEGER(isp), INTENT(IN) :: n !<
513
514	INTEGER(isp) :: mold !<
515	INTEGER(isp) :: nold !<
516	INTEGER(isp) :: ierr !<
517
518	ALLOCATE(reallocate_im(n,m),stat=ierr)
519
520	IF (ierr /= 0) THEN
521	message_string = 'problem in attempt to allocate memory'
522	CALL message( 'reallocate_im', 'PA0454', 1, 2, 0, 6, 0 )
523	END IF
524
525	IF (.not. associated(p)) RETURN
526
527	nold = size(p,1)
528	mold = size(p,2)
529
530	reallocate_im(1:min(nold,n),1:min(mold,m)) = &
531	p(1:min(nold,n),1:min(mold,m))
532
533	DEALLOCATE(p)
534
535	END FUNCTION reallocate_im
536
537	!------------------------------------------------------------------------------!
538	! Description:
539	! ------------
540	!> Reallocates the generators state space "ranseeds" to vectors of size length.
541	!------------------------------------------------------------------------------!
542	FUNCTION arth_i(first,increment,n)
543
544	INTEGER(isp), INTENT(IN) :: first !<
545	INTEGER(isp), INTENT(IN) :: increment !<
546	INTEGER(isp), INTENT(IN) :: n !<
547
548	INTEGER(isp), DIMENSION(n) :: arth_i !<
549
550	INTEGER(isp) :: k !<
551	INTEGER(isp) :: k2 !<
552	INTEGER(isp) :: temp !<
553
554	INTEGER(isp), PARAMETER :: npar_arth=16 !<
555	INTEGER(isp), PARAMETER :: npar2_arth=8 !<
556
557	IF (n > 0) arth_i(1) = first
558
559	IF (n <= npar_arth) THEN
560
561	DO k=2,n
562	arth_i(k) = arth_i(k-1)+increment
563	END DO
564
565	ELSE
566
567	DO k=2,npar2_arth
568	arth_i(k) = arth_i(k-1) + increment
569	END DO
570
571	temp = increment * npar2_arth
572	k = npar2_arth
573
574	DO
575	IF (k >= n) EXIT
576	k2 = k + k
577	arth_i(k+1:min(k2,n)) = temp + arth_i(1:min(k,n-k))
578	temp = temp + temp
579	k = k2
580	END DO
581
582	END IF
583
584	END FUNCTION arth_i
585
586	END MODULE random_generator_parallel

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |