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

Last change on this file since 4439 was 4438, checked in by suehring, 4 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 Giersch $
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!------------------------------------------------------------------------------!
57MODULE 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  (amm*hgng >= 1.0 .or. amm*hgng <= 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
586END MODULE random_generator_parallel
Note: See TracBrowser for help on using the repository browser.