Changeset 4545


Ignore:
Timestamp:
May 22, 2020 1:17:57 PM (11 months ago)
Author:
schwenkel
Message:

Add gaussian random number generator to parallel random generator and using parallel random number in lpm

Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/lagrangian_particle_model_mod.f90

    r4535 r4545  
    2525! -----------------
    2626! $Id$
     27! Using parallel random generator, thus preventing dependency of PE number
     28!
     29! 4535 2020-05-15 12:07:23Z raasch
    2730! bugfix for restart data format query
    2831!
     
    212215               netcdf_handle_error
    213216
    214     USE random_function_mod,                                                   &
    215         ONLY:  random_function
     217    USE random_generator_parallel,                                             &
     218        ONLY:  init_parallel_random_generator,                                 &
     219               random_dummy,                                                   &
     220               random_number_parallel,                                         &
     221               random_number_parallel_gauss,                                   &
     222               random_seed_parallel,                                           &
     223               id_random_array
    216224
    217225    USE restart_data_mpi_io_mod,                                                                   &
     
    273281    INTEGER(iwp) ::  trnp_count_sum                               !< parameter for particle exchange of PEs
    274282    INTEGER(iwp) ::  trnp_count_recv_sum                          !< parameter for particle exchange of PEs
     283
     284    INTEGER(isp), DIMENSION(:,:,:), ALLOCATABLE ::  seq_random_array_particle   !< sequence of random array for particle
    275285
    276286    LOGICAL ::  lagrangian_particle_model = .FALSE.       !< namelist parameter (see documentation)
     
    12441254!--    different on the different PEs.
    12451255       iran_part = iran_part + myid
     1256
     1257!
     1258!--    Initialize parallel random number sequence seed for particles
     1259!--    This is done individually, as thus particle random numbers does
     1260!--    not affect random numbers used for the flow field.
     1261       ALLOCATE ( seq_random_array_particle(5,nys:nyn,nxl:nxr) )
     1262       seq_random_array_particle = 0
     1263
     1264!--    Initializing with random_seed_parallel for every vertical
     1265!--    gridpoint column.
     1266       random_dummy = 0
     1267       DO  i = nxl, nxr
     1268          DO  j = nys, nyn
     1269             CALL random_seed_parallel (random_sequence=id_random_array(j, i))
     1270             CALL random_number_parallel (random_dummy)
     1271             CALL random_seed_parallel (get=seq_random_array_particle(:, j, i))
     1272          ENDDO
     1273       ENDDO
    12461274!
    12471275!--    Create the particle set, and set the initial particles
     
    15071535       DO  ip = nxl, nxr
    15081536          DO  jp = nys, nyn
     1537!
     1538!--          Put the random seeds at grid point jp, ip
     1539             CALL random_seed_parallel( put=seq_random_array_particle(:,jp,ip) )
    15091540             DO  kp = nzb+1, nzt
    15101541                number_of_particles = prt_count(kp,jp,ip)
     
    15191550                DO  n = local_start(kp,jp,ip), number_of_particles
    15201551                   IF ( psl(particles(n)%group) /= psr(particles(n)%group) )  THEN
    1521                       rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
     1552                      CALL random_number_parallel( random_dummy )
     1553                      rand_contr = ( random_dummy - 0.5_wp ) *                 &
    15221554                                     pdx(particles(n)%group)
    15231555                      particles(n)%x = particles(n)%x +                        &
     
    15271559                   ENDIF
    15281560                   IF ( pss(particles(n)%group) /= psn(particles(n)%group) )  THEN
    1529                       rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
     1561                      CALL random_number_parallel( random_dummy )
     1562                      rand_contr = ( random_dummy - 0.5_wp ) *                 &
    15301563                                     pdy(particles(n)%group)
    15311564                      particles(n)%y = particles(n)%y +                        &
     
    15351568                   ENDIF
    15361569                   IF ( psb(particles(n)%group) /= pst(particles(n)%group) )  THEN
    1537                       rand_contr = ( random_function( iran_part ) - 0.5_wp ) * &
     1570                      CALL random_number_parallel( random_dummy )
     1571                      rand_contr = ( random_dummy - 0.5_wp ) *                 &
    15381572                                     pdz(particles(n)%group)
    15391573                      particles(n)%z = particles(n)%z +                        &
     
    15721606                ENDDO
    15731607             ENDDO
     1608!
     1609!--       Get the new random seeds from last call at grid point jp, ip
     1610          CALL random_seed_parallel( get=seq_random_array_particle(:,jp,ip) )
    15741611          ENDDO
    15751612       ENDDO
     
    16991736    DO  ip = nxl, nxr
    17001737       DO  jp = nys, nyn
     1738!
     1739!--       Put the random seeds at grid point jp, ip
     1740          CALL random_seed_parallel( put=seq_random_array_particle(:,jp,ip) )
    17011741          DO  kp = nzb+1, nzt
    17021742
     
    17301770!--             to increase or decrease the number of simulated aerosols
    17311771                particles(n)%weight_factor = particles(n)%weight_factor * aero_weight
    1732 
     1772!
     1773!--             create random numver with parallel number generator
     1774                CALL random_number_parallel( random_dummy )
    17331775                IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) &
    1734                      > random_function( iran_part ) )  THEN
     1776                     > random_dummy )  THEN
    17351777                   particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0_wp
    17361778                ELSE
     
    17761818
    17771819          ENDDO
     1820!
     1821!--    Get the new random seeds from last call at grid point j
     1822       CALL random_seed_parallel( get=seq_random_array_particle(:,jp,ip) )
    17781823       ENDDO
    17791824    ENDDO
     
    21682213                DO  i = nxl, nxr
    21692214                   DO  j = nys, nyn
     2215!
     2216!--                   Put the random seeds at grid point j, i
     2217                      CALL random_seed_parallel( put=seq_random_array_particle(:,j,i) )
     2218
    21702219                      DO  k = nzb+1, nzt
    21712220
     
    22502299
    22512300                      ENDDO
     2301!
     2302!--                   Get the new random seeds from last call at grid point jp, ip
     2303                      CALL random_seed_parallel( get=seq_random_array_particle(:,j,i) )
     2304
    22522305                   ENDDO
    22532306                ENDDO
     
    33943447    REAL(wp) ::  height_p           !< dummy argument for logarithmic interpolation
    33953448    REAL(wp) ::  log_z_z0_int       !< logarithmus used for surface_layer interpolation
    3396     REAL(wp) ::  random_gauss       !< Gaussian-distributed random number used for SGS particle advection
    33973449    REAL(wp) ::  RL                 !< Lagrangian autocorrelation coefficient
    33983450    REAL(wp) ::  rg1                !< Gaussian distributed random number
     
    40104062       DO  nb = 0, 7
    40114063          DO  n = start_index(nb), end_index(nb)
    4012              rg(n,1) = random_gauss( iran_part, 5.0_wp )
    4013              rg(n,2) = random_gauss( iran_part, 5.0_wp )
    4014              rg(n,3) = random_gauss( iran_part, 5.0_wp )
     4064             CALL random_number_parallel_gauss( random_dummy )
     4065             rg(n,1) = random_dummy
     4066             CALL random_number_parallel_gauss( random_dummy )
     4067             rg(n,2) = random_dummy
     4068             CALL random_number_parallel_gauss( random_dummy )
     4069             rg(n,3) = random_dummy
    40154070          ENDDO
    40164071       ENDDO
     
    40544109!--             from becoming unrealistically large.
    40554110                rvar1_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
    4056                                           + 1E-20_wp ) * ( rg(n,1) - 1.0_wp )
     4111                                          + 1E-20_wp ) * rg(n,1)
    40574112                rvar2_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
    4058                                           + 1E-20_wp ) * ( rg(n,2) - 1.0_wp )
     4113                                          + 1E-20_wp ) * rg(n,2)
    40594114                rvar3_temp(n) = SQRT( 2.0_wp * sgs_wf_part * e_int(n)          &
    4060                                           + 1E-20_wp ) * ( rg(n,3) - 1.0_wp )
     4115                                          + 1E-20_wp ) * rg(n,3)
    40614116             ELSE
    40624117!
     
    42444299                                             1.0E-20_wp ) )
    42454300                      sigma          = SQRT( e(kp,jp,ip) )
    4246 
    4247                       rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    4248                       rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    4249                       rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     4301!
     4302!--                   Calculate random component of particle sgs velocity using parallel
     4303!--                   random generator
     4304                      CALL random_number_parallel_gauss( random_dummy )
     4305                      rg1 = random_dummy
     4306                      CALL random_number_parallel_gauss( random_dummy )
     4307                      rg2 = random_dummy
     4308                      CALL random_number_parallel_gauss( random_dummy )
     4309                      rg3 = random_dummy
    42504310
    42514311                      particles(n)%rvar1 = RL * particles(n)%rvar1 +              &
     
    43404400                    sigma          = SQRT( e(kp,jp,ip) )
    43414401
    4342                     rg1 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    4343                     rg2 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
    4344                     rg3 = random_gauss( iran_part, 5.0_wp ) - 1.0_wp
     4402!
     4403!--                 Calculate random component of particle sgs velocity using parallel
     4404!--                 random generator
     4405                    CALL random_number_parallel_gauss( random_dummy )
     4406                    rg1 = random_dummy
     4407                    CALL random_number_parallel_gauss( random_dummy )
     4408                    rg2 = random_dummy
     4409                    CALL random_number_parallel_gauss( random_dummy )
     4410                    rg3 = random_dummy
    43454411
    43464412                    particles(n)%rvar1 = RL * particles(n)%rvar1 +                &
     
    46034669    REAL(wp) ::  prt_y          !< current particle position (y)
    46044670    REAL(wp) ::  prt_z          !< current particle position (z)
    4605     REAL(wp) ::  ran_val        !< location of wall in z
    46064671    REAL(wp) ::  reset_top      !< location of wall in z
    46074672    REAL(wp) ::  t_old          !< previous reflection time
     
    46824747                index_reset = MINLOC( prt_count(nzb+1:particles_top,j,i), DIM = 1 )
    46834748                reset_top = zu(index_reset)
    4684                 iran_part = iran_part + myid
    4685                 ran_val = random_function( iran_part )
    4686                 particles(n)%z       = reset_top *  ( 1.0  + ( ran_val / 10.0_wp) )
     4749                CALL random_number_parallel( random_dummy )
     4750                particles(n)%z       = reset_top *  ( 1.0  + ( random_dummy / 10.0_wp) )
    46874751                particles(n)%speed_z = 0.0_wp
    46884752                IF ( curvature_solution_effects )  THEN
     
    55825646!--          Calculate the number of collections and consider multiple collections.
    55835647!--          (Accordingly, p_crit will be 0.0, 1.0, 2.0, ...)
     5648             CALL random_number_parallel( random_dummy )
    55845649             IF ( collection_probability - FLOOR(collection_probability)    &
    5585                   > random_function( iran_part ) )  THEN
     5650                  > random_dummy )  THEN
    55865651                collection_probability = FLOOR(collection_probability) + 1.0_wp
    55875652             ELSE
  • palm/trunk/SOURCE/random_generator_parallel_mod.f90

    r4438 r4545  
    2525! -----------------
    2626! $Id$
     27! Add generator (using parallel mode) returning gaussian distributed random
     28! number
     29!
     30! 4438 2020-03-03 20:49:28Z raasch
    2731! - Rename variables to avoid confusion with subdomain grid indices
    2832! - Some formatting adjustments
     
    4549!> ran_parallel returns a uniform random deviate between 0.0 and 1.0
    4650!> (exclusive of the end point values).
     51!> Moreover, it contains a routine returning a normally distributed random number
     52!> with mean zero and unity standard deviation.
    4753!> Additionally it provides the generator with five integer for use as initial state space.
    4854!> The first tree integers (iran, jran, kran) are maintained as non negative values,
     
    95101   END INTERFACE
    96102
     103   INTERFACE random_number_parallel_gauss
     104      MODULE PROCEDURE gasdev_s
     105   END INTERFACE
     106
    97107   INTERFACE random_seed_parallel
    98108      MODULE PROCEDURE random_seed_parallel
     
    114124
    115125   PUBLIC random_number_parallel, random_seed_parallel, random_dummy,          &
    116           id_random_array, seq_random_array, init_parallel_random_generator
     126          id_random_array, seq_random_array, init_parallel_random_generator,   &
     127          random_number_parallel_gauss
    117128
    118129 CONTAINS
     
    248259! Description:
    249260! ------------
     261!> Returns in harvest a normally distributed deviate with zero mean and unit
     262!> variance, using ran0_s as the source of uniform deviates. Following
     263!> Numerical Recipes in Fortran90 (Press et al., 2nd edition, 1996, pp 1152ff).
     264!> Note, instead of ran1_s, ran0_s is used.
     265!------------------------------------------------------------------------------!
     266   SUBROUTINE gasdev_s(harvest)
     267
     268      REAL(wp), INTENT(OUT) ::  harvest   !<
     269
     270      REAL(wp) ::  rsq      !<
     271      REAL(wp) ::  v1       !<
     272      REAL(wp) ::  v2       !<
     273      REAL(wp), SAVE ::  g  !<
     274
     275      LOGICAL, SAVE ::  gaus_stored = .FALSE. !<
     276!
     277!--   We have an extra deviate handy, so return it, and unset the flag.
     278      IF (gaus_stored)  THEN
     279         harvest = g
     280         gaus_stored = .FALSE.
     281!
     282!--   We don’t have an extra deviate handy, so pick two uniform numbers in the
     283!--   square extending from -1 to +1 in each direction
     284      ELSE
     285         DO
     286            CALL ran0_s(v1)
     287            CALL ran0_s(v2)
     288            v1 = 2.0_wp * v1 - 1.0_wp
     289            v2 = 2.0_wp * v2 - 1.0_wp
     290!
     291!--         see if they are in the unit circle
     292            rsq = v1**2 + v2**2
     293!
     294!--         otherwise try again.
     295            IF (rsq > 0.0  .AND.  rsq < 1.0) EXIT
     296         ENDDO
     297!
     298!--      Now make the Box-Muller transformation to get two normal deviates.
     299!--      Return one and save the other for next time. Set flag.
     300         rsq = SQRT(-2.0_sp * LOG(rsq)/rsq)
     301         harvest = v1 * rsq
     302         g = v2 * rsq
     303         gaus_stored = .TRUE.
     304      ENDIF
     305
     306   END SUBROUTINE gasdev_s
     307
     308!------------------------------------------------------------------------------!
     309! Description:
     310! ------------
    250311!> Initialize or reinitialize the random generator state space to vectors of size length.
    251312!> The saved variable seq is hashed (via calls to the module routine ran_hash)
Note: See TracChangeset for help on using the changeset viewer.