Changeset 4545
- Timestamp:
- May 22, 2020 1:17:57 PM (4 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
lagrangian_particle_model_mod.f90 (modified) (22 diffs)
-
random_generator_parallel_mod.f90 (modified) (5 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/lagrangian_particle_model_mod.f90
r4535 r4545 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Using parallel random generator, thus preventing dependency of PE number 28 ! 29 ! 4535 2020-05-15 12:07:23Z raasch 27 30 ! bugfix for restart data format query 28 31 ! … … 212 215 netcdf_handle_error 213 216 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 216 224 217 225 USE restart_data_mpi_io_mod, & … … 273 281 INTEGER(iwp) :: trnp_count_sum !< parameter for particle exchange of PEs 274 282 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 275 285 276 286 LOGICAL :: lagrangian_particle_model = .FALSE. !< namelist parameter (see documentation) … … 1244 1254 !-- different on the different PEs. 1245 1255 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 1246 1274 ! 1247 1275 !-- Create the particle set, and set the initial particles … … 1507 1535 DO ip = nxl, nxr 1508 1536 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) ) 1509 1540 DO kp = nzb+1, nzt 1510 1541 number_of_particles = prt_count(kp,jp,ip) … … 1519 1550 DO n = local_start(kp,jp,ip), number_of_particles 1520 1551 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 ) * & 1522 1554 pdx(particles(n)%group) 1523 1555 particles(n)%x = particles(n)%x + & … … 1527 1559 ENDIF 1528 1560 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 ) * & 1530 1563 pdy(particles(n)%group) 1531 1564 particles(n)%y = particles(n)%y + & … … 1535 1568 ENDIF 1536 1569 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 ) * & 1538 1572 pdz(particles(n)%group) 1539 1573 particles(n)%z = particles(n)%z + & … … 1572 1606 ENDDO 1573 1607 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) ) 1574 1611 ENDDO 1575 1612 ENDDO … … 1699 1736 DO ip = nxl, nxr 1700 1737 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) ) 1701 1741 DO kp = nzb+1, nzt 1702 1742 … … 1730 1770 !-- to increase or decrease the number of simulated aerosols 1731 1771 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 ) 1733 1775 IF ( particles(n)%weight_factor - FLOOR(particles(n)%weight_factor,KIND=wp) & 1734 > random_ function( iran_part )) THEN1776 > random_dummy ) THEN 1735 1777 particles(n)%weight_factor = FLOOR(particles(n)%weight_factor,KIND=wp) + 1.0_wp 1736 1778 ELSE … … 1776 1818 1777 1819 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) ) 1778 1823 ENDDO 1779 1824 ENDDO … … 2168 2213 DO i = nxl, nxr 2169 2214 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 2170 2219 DO k = nzb+1, nzt 2171 2220 … … 2250 2299 2251 2300 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 2252 2305 ENDDO 2253 2306 ENDDO … … 3394 3447 REAL(wp) :: height_p !< dummy argument for logarithmic interpolation 3395 3448 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 advection3397 3449 REAL(wp) :: RL !< Lagrangian autocorrelation coefficient 3398 3450 REAL(wp) :: rg1 !< Gaussian distributed random number … … 4010 4062 DO nb = 0, 7 4011 4063 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 4015 4070 ENDDO 4016 4071 ENDDO … … 4054 4109 !-- from becoming unrealistically large. 4055 4110 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) 4057 4112 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) 4059 4114 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) 4061 4116 ELSE 4062 4117 ! … … 4244 4299 1.0E-20_wp ) ) 4245 4300 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 4250 4310 4251 4311 particles(n)%rvar1 = RL * particles(n)%rvar1 + & … … 4340 4400 sigma = SQRT( e(kp,jp,ip) ) 4341 4401 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 4345 4411 4346 4412 particles(n)%rvar1 = RL * particles(n)%rvar1 + & … … 4603 4669 REAL(wp) :: prt_y !< current particle position (y) 4604 4670 REAL(wp) :: prt_z !< current particle position (z) 4605 REAL(wp) :: ran_val !< location of wall in z4606 4671 REAL(wp) :: reset_top !< location of wall in z 4607 4672 REAL(wp) :: t_old !< previous reflection time … … 4682 4747 index_reset = MINLOC( prt_count(nzb+1:particles_top,j,i), DIM = 1 ) 4683 4748 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) ) 4687 4751 particles(n)%speed_z = 0.0_wp 4688 4752 IF ( curvature_solution_effects ) THEN … … 5582 5646 !-- Calculate the number of collections and consider multiple collections. 5583 5647 !-- (Accordingly, p_crit will be 0.0, 1.0, 2.0, ...) 5648 CALL random_number_parallel( random_dummy ) 5584 5649 IF ( collection_probability - FLOOR(collection_probability) & 5585 > random_ function( iran_part )) THEN5650 > random_dummy ) THEN 5586 5651 collection_probability = FLOOR(collection_probability) + 1.0_wp 5587 5652 ELSE -
palm/trunk/SOURCE/random_generator_parallel_mod.f90
r4438 r4545 25 25 ! ----------------- 26 26 ! $Id$ 27 ! Add generator (using parallel mode) returning gaussian distributed random 28 ! number 29 ! 30 ! 4438 2020-03-03 20:49:28Z raasch 27 31 ! - Rename variables to avoid confusion with subdomain grid indices 28 32 ! - Some formatting adjustments … … 45 49 !> ran_parallel returns a uniform random deviate between 0.0 and 1.0 46 50 !> (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. 47 53 !> Additionally it provides the generator with five integer for use as initial state space. 48 54 !> The first tree integers (iran, jran, kran) are maintained as non negative values, … … 95 101 END INTERFACE 96 102 103 INTERFACE random_number_parallel_gauss 104 MODULE PROCEDURE gasdev_s 105 END INTERFACE 106 97 107 INTERFACE random_seed_parallel 98 108 MODULE PROCEDURE random_seed_parallel … … 114 124 115 125 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 117 128 118 129 CONTAINS … … 248 259 ! Description: 249 260 ! ------------ 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 ! ------------ 250 311 !> Initialize or reinitialize the random generator state space to vectors of size length. 251 312 !> The saved variable seq is hashed (via calls to the module routine ran_hash)
Note: See TracChangeset
for help on using the changeset viewer.
