Home

Context Navigation

source: palm/trunk/SOURCE/synthetic_turbulence_generator_mod.f90 @ 2670

Last change on this file since 2670 was 2669, checked in by raasch, 7 years ago
file attributes and activation strings in .palm.iofiles revised, file extensions for nesting, masked output, wind turbine data, etc. changed
Property svn:keywords set to `Id`
File size: 35.4 KB

Line
1	!> @synthetic_turbulence_generator_mod.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM-4U.
4	!
5	! PALM-4U 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-4U 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 2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: synthetic_turbulence_generator_mod.f90 2669 2017-12-06 16:03:27Z raasch $
27	! unit number for file containing turbulence generator data changed to 90
28	! bugfix: preprocessor directives added for MPI specific code
29	!
30	! 2576 2017-10-24 13:49:46Z Giersch
31	! Definition of a new function called stg_skip_var_list to skip module
32	! parameters during reading restart data
33	!
34	! 2563 2017-10-19 15:36:10Z Giersch
35	! stg_read_restart_data is called in stg_parin in the case of a restart run
36	!
37	! 2259 2017-06-08 09:09:11Z gronemeier
38	! Initial revision
39	!
40	!
41	!
42	! Authors:
43	! --------
44	! @author Tobias Gronemeier, Atsushi Inagaki, Micha Gryschka, Christoph Knigge
45	!
46	!
47	! Description:
48	! ------------
49	!> The module generates turbulence at the inflow boundary based on a method by
50	!> Xie and Castro (2008) utilizing a Lund rotation (Lund, 1998) and a mass-flux
51	!> correction by Kim et al. (2013).
52	!> The turbulence is correlated based on length scales in y- and z-direction and
53	!> a time scale for each velocity component. The profiles of length and time
54	!> scales, mean u, v, w, e and pt, and all components of the Reynolds stress
55	!> tensor are read from file STG_PROFILES.
56	!>
57	!> @todo test restart
58	!> enable cyclic_fill
59	!> implement turbulence generation for e and pt
60	!> @note <Enter notes on the module>
61	!> @bug Height information from input file is not used. Profiles from input
62	!> must match with current PALM grid.
63	!> Transformation of length scales to number of gridpoints does not
64	!> consider grid stretching.
65	!> In case of restart, velocity seeds differ from precursor run if a11,
66	!> a22, or a33 are zero.
67	!------------------------------------------------------------------------------!
68	MODULE synthetic_turbulence_generator_mod
69
70
71	USE arrays_3d, &
72	ONLY: mean_inflow_profiles, u, v, w
73
74	USE constants, &
75	ONLY: pi
76
77	USE control_parameters, &
78	ONLY: initializing_actions, message_string, &
79	synthetic_turbulence_generator
80
81	USE cpulog, &
82	ONLY: cpu_log, log_point, log_point_s
83
84	USE indices, &
85	ONLY: nbgp, nzb, nzt, nyng, nysg
86
87	USE kinds
88
89	USE MPI
90
91	USE pegrid, &
92	ONLY: comm1dx, comm1dy, comm2d, id_inflow, ierr, myidx, pdims
93
94	USE transpose_indices, &
95	ONLY: nzb_x, nzt_x
96
97
98	IMPLICIT NONE
99
100	LOGICAL :: velocity_seed_initialized = .FALSE. !< true after first call of stg_main
101	LOGICAL :: use_synthetic_turbulence_generator = .FALSE. !< switch to use synthetic turbulence generator
102
103	INTEGER(iwp) :: stg_type_yz !< MPI type for full z range
104	INTEGER(iwp) :: stg_type_yz_small !< MPI type for small z range
105	INTEGER(iwp) :: merg !< maximum length scale (in gp)
106	INTEGER(iwp) :: mergp !< merg + nbgp
107
108	REAL(wp) :: mc_factor !< mass flux correction factor
109
110	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: displs !< displacement for MPI_GATHERV
111	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: recv_count !< receive count for MPI_GATHERV
112	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuy !< length scale of u in y direction (in gp)
113	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nuz !< length scale of u in z direction (in gp)
114	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvy !< length scale of v in y direction (in gp)
115	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nvz !< length scale of v in z direction (in gp)
116	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwy !< length scale of w in y direction (in gp)
117	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: nwz !< length scale of w in z direction (in gp)
118
119	INTEGER(iwp), DIMENSION(:), ALLOCATABLE :: seed !< seed of random number for rn-generator
120
121	REAL(wp), DIMENSION(:), ALLOCATABLE :: a11 !< coefficient for Lund rotation
122	REAL(wp), DIMENSION(:), ALLOCATABLE :: a21 !< coefficient for Lund rotation
123	REAL(wp), DIMENSION(:), ALLOCATABLE :: a22 !< coefficient for Lund rotation
124	REAL(wp), DIMENSION(:), ALLOCATABLE :: a31 !< coefficient for Lund rotation
125	REAL(wp), DIMENSION(:), ALLOCATABLE :: a32 !< coefficient for Lund rotation
126	REAL(wp), DIMENSION(:), ALLOCATABLE :: a33 !< coefficient for Lund rotation
127	REAL(wp), DIMENSION(:), ALLOCATABLE :: tu !< Lagrangian time scale of u
128	REAL(wp), DIMENSION(:), ALLOCATABLE :: tv !< Lagrangian time scale of v
129	REAL(wp), DIMENSION(:), ALLOCATABLE :: tw !< Lagrangian time scale of w
130
131	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buy !< filter function for u in y direction
132	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: buz !< filter function for u in z direction
133	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvy !< filter function for v in y direction
134	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bvz !< filter function for v in z direction
135	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwy !< filter function for w in y direction
136	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: bwz !< filter function for w in z direction
137	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fu !< velocity seed for u
138	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fuo !< velocity seed for u with new random number
139	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fv !< velocity seed for v
140	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fvo !< velocity seed for v with new random number
141	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fw !< velocity seed for w
142	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: fwo !< velocity seed for w with new random number
143
144
145	!
146	!-- PALM interfaces:
147	!-- Input parameter checks to be done in check_parameters
148	INTERFACE stg_check_parameters
149	MODULE PROCEDURE stg_check_parameters
150	END INTERFACE stg_check_parameters
151
152	!
153	!-- Calculate filter functions
154	INTERFACE stg_filter_func
155	MODULE PROCEDURE stg_filter_func
156	END INTERFACE stg_filter_func
157
158	!
159	!-- Generate velocity seeds
160	INTERFACE stg_generate_seed_yz
161	MODULE PROCEDURE stg_generate_seed_yz
162	END INTERFACE stg_generate_seed_yz
163
164	!
165	!-- Output of information to the header file
166	INTERFACE stg_header
167	MODULE PROCEDURE stg_header
168	END INTERFACE stg_header
169
170	!
171	!-- Initialization actions
172	INTERFACE stg_init
173	MODULE PROCEDURE stg_init
174	END INTERFACE stg_init
175
176	!
177	!-- Main procedure of synth. turb. gen.
178	INTERFACE stg_main
179	MODULE PROCEDURE stg_main
180	END INTERFACE stg_main
181
182	!
183	!-- Reading of NAMELIST parameters
184	INTERFACE stg_parin
185	MODULE PROCEDURE stg_parin
186	END INTERFACE stg_parin
187
188	!
189	!-- Skipping of parameters for restart runs
190	INTERFACE stg_skip_var_list
191	MODULE PROCEDURE stg_skip_var_list
192	END INTERFACE stg_skip_var_list
193
194	!
195	!-- Reading of parameters for restart runs
196	INTERFACE stg_read_restart_data
197	MODULE PROCEDURE stg_read_restart_data
198	END INTERFACE stg_read_restart_data
199
200	!
201	!-- Writing of binary output for restart runs
202	INTERFACE stg_write_restart_data
203	MODULE PROCEDURE stg_write_restart_data
204	END INTERFACE stg_write_restart_data
205
206	SAVE
207
208	PRIVATE
209
210	!
211	!-- Public interfaces
212	PUBLIC stg_check_parameters, stg_header, stg_init, stg_main, stg_parin, &
213	stg_write_restart_data, stg_skip_var_list
214
215	!
216	!-- Public variables
217	PUBLIC use_synthetic_turbulence_generator
218
219
220	CONTAINS
221
222
223	!------------------------------------------------------------------------------!
224	! Description:
225	! ------------
226	!> Check parameters routine for synthetic turbulence generator
227	!------------------------------------------------------------------------------!
228	SUBROUTINE stg_check_parameters
229
230
231	USE control_parameters, &
232	ONLY: bc_lr, bc_ns, turbulent_inflow
233
234
235	IMPLICIT NONE
236
237	IF ( use_synthetic_turbulence_generator ) THEN
238
239	IF ( INDEX( initializing_actions, 'set_constant_profiles' ) == 0 .AND. &
240	INDEX( initializing_actions, 'read_restart_data' ) == 0 ) THEN
241	message_string = 'Using synthetic turbulence generator ' // &
242	'requires &initializing_actions = ' // &
243	'"set_constant_profiles" or "read_restart_data"'
244	CALL message( 'stg_check_parameters', 'PA0015', 1, 2, 0, 6, 0 )
245	ENDIF
246
247	IF ( bc_lr /= 'dirichlet/radiation' ) THEN
248	message_string = 'Using synthetic turbulence generator ' // &
249	'requires &bc_lr = "dirichlet/radiation"'
250	CALL message( 'stg_check_parameters', 'PA0035', 1, 2, 0, 6, 0 )
251	ENDIF
252	IF ( bc_ns /= 'cyclic' ) THEN
253	message_string = 'Using synthetic turbulence generator ' // &
254	'requires &bc_ns = "cyclic"'
255	CALL message( 'stg_check_parameters', 'PA0037', 1, 2, 0, 6, 0 )
256	ENDIF
257	IF ( turbulent_inflow ) THEN
258	message_string = 'Using synthetic turbulence generator ' // &
259	'in combination &with turbulent_inflow = .T. ' // &
260	'is not allowed'
261	CALL message( 'stg_check_parameters', 'PA0039', 1, 2, 0, 6, 0 )
262	ENDIF
263
264	ENDIF
265
266	END SUBROUTINE stg_check_parameters
267
268
269	!------------------------------------------------------------------------------!
270	! Description:
271	! ------------
272	!> Header output for synthetic turbulence generator
273	!------------------------------------------------------------------------------!
274	SUBROUTINE stg_header ( io )
275
276
277	IMPLICIT NONE
278
279	INTEGER(iwp), INTENT(IN) :: io !< Unit of the output file
280
281	!
282	!-- Write synthetic turbulence generator Header
283	WRITE( io, 1 )
284	IF ( use_synthetic_turbulence_generator ) THEN
285	WRITE( io, 2 )
286	ELSE
287	WRITE( io, 3 )
288	ENDIF
289
290	1 FORMAT (//' Synthetic turbulence generator information:'/ &
291	' ------------------------------------------'/)
292	2 FORMAT (' synthetic turbulence generator switched on')
293	3 FORMAT (' synthetic turbulence generator switched off')
294
295	END SUBROUTINE stg_header
296
297
298	!------------------------------------------------------------------------------!
299	! Description:
300	! ------------
301	!> Initialization of the synthetic turbulence generator
302	!------------------------------------------------------------------------------!
303	SUBROUTINE stg_init
304
305
306	USE arrays_3d, &
307	ONLY: u_init, v_init
308
309	USE control_parameters, &
310	ONLY: coupling_char, dz, e_init
311
312	USE grid_variables, &
313	ONLY: ddy
314
315
316	IMPLICIT NONE
317
318	#if defined( __parallel )
319	INTEGER(KIND=MPI_ADDRESS_KIND) :: extent !< extent of new MPI type
320	INTEGER(KIND=MPI_ADDRESS_KIND) :: tob=0 !< dummy variable
321	#endif
322
323	INTEGER(iwp) :: j !> loop index
324	INTEGER(iwp) :: k !< index
325	INTEGER(iwp) :: newtype !< dummy MPI type
326	INTEGER(iwp) :: realsize !< size of REAL variables
327	INTEGER(iwp) :: nseed !< dimension of random number seed
328	INTEGER(iwp) :: startseed = 1234567890 !< start seed for random number generator
329
330	!
331	!-- Dummy variables used for reading profiles
332	REAL(wp) :: d1 !< u profile
333	REAL(wp) :: d2 !< v profile
334	REAL(wp) :: d3 !< w profile
335	REAL(wp) :: d5 !< e profile
336	REAL(wp) :: d11 !< vertical interpolation for a11
337	REAL(wp) :: d21 !< vertical interpolation for a21
338	REAL(wp) :: d22 !< vertical interpolation for a22
339	REAL(wp) :: luy !< length scale for u in y direction
340	REAL(wp) :: luz !< length scale for u in z direction
341	REAL(wp) :: lvy !< length scale for v in y direction
342	REAL(wp) :: lvz !< length scale for v in z direction
343	REAL(wp) :: lwy !< length scale for w in y direction
344	REAL(wp) :: lwz !< length scale for w in z direction
345	REAL(wp) :: zz !< height
346
347	REAL(wp),DIMENSION(nzb:nzt+1) :: r11 !< Reynolds parameter
348	REAL(wp),DIMENSION(nzb:nzt+1) :: r21 !< Reynolds parameter
349	REAL(wp),DIMENSION(nzb:nzt+1) :: r22 !< Reynolds parameter
350	REAL(wp),DIMENSION(nzb:nzt+1) :: r31 !< Reynolds parameter
351	REAL(wp),DIMENSION(nzb:nzt+1) :: r32 !< Reynolds parameter
352	REAL(wp),DIMENSION(nzb:nzt+1) :: r33 !< Reynolds parameter
353
354
355	#if defined( __parallel )
356	CALL MPI_BARRIER( comm2d, ierr )
357	#endif
358
359	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
360
361	IF ( .NOT. ALLOCATED( mean_inflow_profiles ) ) &
362	ALLOCATE( mean_inflow_profiles(nzb:nzt+1,5) )
363
364	ALLOCATE ( a11(nzb:nzt+1), a21(nzb:nzt+1), a22(nzb:nzt+1), &
365	a31(nzb:nzt+1), a32(nzb:nzt+1), a33(nzb:nzt+1), &
366	nuy(nzb:nzt+1), nuz(nzb:nzt+1), nvy(nzb:nzt+1), &
367	nvz(nzb:nzt+1), nwy(nzb:nzt+1), nwz(nzb:nzt+1), &
368	tu(nzb:nzt+1), tv(nzb:nzt+1), tw(nzb:nzt+1) )
369
370	#if defined( __parallel )
371	!
372	!-- Define MPI type used in stg_generate_seed_yz to gather vertical splitted
373	!-- velocity seeds
374	CALL MPI_TYPE_SIZE( MPI_REAL, realsize, ierr )
375	extent = 1 * realsize
376
377	! stg_type_yz: yz-slice with vertical bounds nzb:nzt+1
378	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt-nzb+2,nyng-nysg+1], &
379	[1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
380	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz, ierr )
381	CALL MPI_TYPE_COMMIT( stg_type_yz, ierr )
382	CALL MPI_TYPE_FREE( newtype, ierr )
383
384	! stg_type_yz_small: yz-slice with vertical bounds nzb_x:nzt_x+1
385	CALL MPI_TYPE_CREATE_SUBARRAY( 2, [nzt_x-nzb_x+2,nyng-nysg+1], &
386	[1,nyng-nysg+1], [0,0], MPI_ORDER_FORTRAN, MPI_REAL, newtype, ierr )
387	CALL MPI_TYPE_CREATE_RESIZED( newtype, tob, extent, stg_type_yz_small, ierr )
388	CALL MPI_TYPE_COMMIT( stg_type_yz_small, ierr )
389	CALL MPI_TYPE_FREE( newtype, ierr )
390
391	! receive count and displacement for MPI_GATHERV in stg_generate_seed_yz
392	ALLOCATE( recv_count(pdims(1)), displs(pdims(1)) )
393
394	recv_count = nzt_x-nzb_x + 1
395	recv_count(pdims(1)) = recv_count(pdims(1)) + 1
396
397	DO j = 1, pdims(1)
398	displs(j) = 0 + (nzt_x-nzb_x+1) * (j-1)
399	ENDDO
400	#endif
401
402	!
403	!-- Define seed of random number
404	CALL RANDOM_SEED()
405	CALL RANDOM_SEED( size=nseed )
406	ALLOCATE( seed(1:nseed) )
407	DO j = 1, nseed
408	seed(j) = startseed + j
409	ENDDO
410	CALL RANDOM_SEED(put=seed)
411
412	!
413	!-- Read inflow profile
414	!-- Height levels of profiles in input profile is as follows:
415	!-- zu: luy, luz, tu, lvy, lvz, tv, r11, r21, r22, d1, d2, d5
416	!-- zw: lwy, lwz, tw, r31, r32, r33, d3
417	!-- WARNING: zz is not used at the moment
418	OPEN( 90, FILE='STG_PROFILES'//TRIM( coupling_char ), STATUS='OLD', &
419	FORM='FORMATTED')
420
421	! Skip header
422	READ( 90, * )
423
424	DO k = nzb, nzt+1
425	READ( 90, * ) zz, luy, luz, tu(k), lvy, lvz, tv(k), lwy, lwz, tw(k), &
426	r11(k), r21(k), r22(k), r31(k), r32(k), r33(k), &
427	d1, d2, d3, d5
428
429	!
430	!-- Convert length scales from meter to number of grid points
431	nuy(k) = INT( luy * ddy )
432	nuz(k) = INT( luz / dz )
433	nvy(k) = INT( lvy * ddy )
434	nvz(k) = INT( lvz / dz )
435	nwy(k) = INT( lwy * ddy )
436	nwz(k) = INT( lwz / dz )
437
438	!
439	!-- Save Mean inflow profiles
440	IF ( TRIM( initializing_actions ) /= 'read_restart_data' ) THEN
441	mean_inflow_profiles(k,1) = d1
442	mean_inflow_profiles(k,2) = d2
443	! mean_inflow_profiles(k,4) = d4
444	mean_inflow_profiles(k,5) = d5
445	ENDIF
446	ENDDO
447
448	CLOSE( 90 )
449
450	!
451	!-- Assign initial profiles
452	u_init = mean_inflow_profiles(:,1)
453	v_init = mean_inflow_profiles(:,2)
454	! pt_init = mean_inflow_profiles(:,4)
455	e_init = MAXVAL( mean_inflow_profiles(:,5) )
456
457	!
458	!-- Calculate coefficient matrix from Reynolds stress tensor (Lund rotation)
459	DO k = nzb, nzt+1
460	IF ( r11(k) > 0.0_wp ) THEN
461	a11(k) = SQRT( r11(k) )
462	a21(k) = r21(k) / a11(k)
463	ELSE
464	a11(k) = 0.0_wp
465	a21(k) = 0.0_wp
466	ENDIF
467
468	a22(k) = r22(k) - a21(k)**2
469	IF ( a22(k) > 0.0_wp ) THEN
470	a22(k) = SQRT( a22(k) )
471	a32(k) = ( r32(k) - a21(k) * a31(k) ) / a22(k)
472	ELSE
473	a22(k) = 0.0_wp
474	a32(k) = 0.0_wp
475	ENDIF
476
477	!
478	!-- a31, a32, a33 must be calculated with interpolated a11, a21, a22 (d11,
479	!-- d21, d22) because of different vertical grid
480	IF ( k .le. nzt ) THEN
481	d11 = 0.5_wp * ( r11(k) + r11(k+1) )
482	IF ( d11 > 0.0_wp ) THEN
483	d11 = SQRT( d11 )
484	d21 = ( 0.5_wp * ( r21(k) + r21(k+1) ) ) / d11
485	a31(k) = r31(k) / d11
486	ELSE
487	d21 = 0.0_wp
488	a31(k) = 0.0_wp
489	ENDIF
490
491	d22 = 0.5_wp * ( r22(k) + r22(k+1) ) - d21 ** 2
492	IF ( d22 > 0.0_wp ) THEN
493	a32(k) = ( r32(k) - d21 * a31(k) ) / SQRT( d22 )
494	ELSE
495	a32(k) = 0.0_wp
496	ENDIF
497
498	a33(k) = r33(k) - a31(k) 2 - a32(k) 2
499	IF ( a33(k) > 0.0_wp ) THEN
500	a33(k) = SQRT( a33(k) )
501	ELSE
502	a33(k) = 0.0_wp
503	ENDIF
504	ELSE
505	a31(k) = a31(k-1)
506	a32(k) = a32(k-1)
507	a33(k) = a33(k-1)
508	ENDIF
509
510	ENDDO
511
512	!
513	!-- Define the size of the filter functions and allocate them.
514	merg = 0
515
516	! arrays must be large enough to cover the largest length scale
517	DO k = nzb, nzt+1
518	j = MAX( ABS(nuy(k)), ABS(nuz(k)), &
519	ABS(nvy(k)), ABS(nvz(k)), &
520	ABS(nwy(k)), ABS(nwz(k)) )
521	IF ( j > merg ) merg = j
522	ENDDO
523
524	merg = 2 * merg
525	mergp = merg + nbgp
526
527	ALLOCATE ( buy(-merg:merg,nzb:nzt+1), buz(-merg:merg,nzb:nzt+1), &
528	bvy(-merg:merg,nzb:nzt+1), bvz(-merg:merg,nzb:nzt+1), &
529	bwy(-merg:merg,nzb:nzt+1), bwz(-merg:merg,nzb:nzt+1) )
530
531	!
532	!-- Allocate velocity seeds
533	ALLOCATE ( fu( nzb:nzt+1,nysg:nyng), fuo(nzb:nzt+1,nysg:nyng), &
534	fv( nzb:nzt+1,nysg:nyng), fvo(nzb:nzt+1,nysg:nyng), &
535	fw( nzb:nzt+1,nysg:nyng), fwo(nzb:nzt+1,nysg:nyng) )
536
537	fu = 0._wp
538	fuo = 0._wp
539	fv = 0._wp
540	fvo = 0._wp
541	fw = 0._wp
542	fwo = 0._wp
543
544	!
545	!-- Create filter functions
546	CALL stg_filter_func( nuy, buy ) !filter uy
547	CALL stg_filter_func( nuz, buz ) !filter uz
548	CALL stg_filter_func( nvy, bvy ) !filter vy
549	CALL stg_filter_func( nvz, bvz ) !filter vz
550	CALL stg_filter_func( nwy, bwy ) !filter wy
551	CALL stg_filter_func( nwz, bwz ) !filter wz
552
553	#if defined( __parallel )
554	CALL MPI_BARRIER( comm2d, ierr )
555	#endif
556
557	!
558	!-- In case of restart, calculate velocity seeds fu, fv, fw from former
559	! time step.
560	! Bug: fu, fv, fw are different in those heights where a11, a22, a33
561	! are 0 compared to the prerun. This is mostly for k=nzt+1.
562	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
563	IF ( myidx == id_inflow ) THEN
564	DO j = nysg, nyng
565	DO k = nzb, nzt+1
566
567	IF ( a11(k) .NE. 0._wp ) THEN
568	fu(k,j) = ( u(k,j,-1) / mc_factor - &
569	mean_inflow_profiles(k,1) ) / a11(k)
570	ELSE
571	fu(k,j) = 0._wp
572	ENDIF
573
574	IF ( a22(k) .NE. 0._wp ) THEN
575	fv(k,j) = ( v(k,j,-1) / mc_factor - a21(k) * fu(k,j) - &
576	mean_inflow_profiles(k,2) ) / a22(k)
577	ELSE
578	fv(k,j) = 0._wp
579	ENDIF
580
581	IF ( a33(k) .NE. 0._wp ) THEN
582	fw(k,j) = ( w(k,j,-1) / mc_factor - a31(k) * fu(k,j) - &
583	a32(k) * fv(k,j) ) / a33(k)
584	ELSE
585	fw = 0._wp
586	ENDIF
587
588	ENDDO
589	ENDDO
590	ENDIF
591	ENDIF
592
593	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
594
595	END SUBROUTINE stg_init
596
597
598	!------------------------------------------------------------------------------!
599	! Description:
600	! ------------
601	!> Calculate filter function bxx from length scale nxx following Eg.9 and 10
602	!> (Xie and Castro, 2008)
603	!------------------------------------------------------------------------------!
604	SUBROUTINE stg_filter_func( nxx, bxx )
605
606
607	IMPLICIT NONE
608
609	INTEGER(iwp) :: k !< loop index
610	INTEGER(iwp) :: n_k !< length scale nXX in height k
611	INTEGER(iwp) :: n_k2 !< n_k * 2
612	INTEGER(iwp) :: nf !< index for length scales
613
614	REAL(wp) :: bdenom !< denominator for filter functions bXX
615	REAL(wp) :: qsi = 1.0_wp !< minimization factor
616
617	INTEGER(iwp), DIMENSION(:) :: nxx(nzb:nzt+1) !< length scale (in gp)
618
619	REAL(wp), DIMENSION(:,:) :: bxx(-merg:merg,nzb:nzt+1) !< filter function
620
621
622	bxx = 0.0_wp
623
624	DO k = nzb, nzt+1
625	bdenom = 0.0_wp
626	n_k = nxx(k)
627	IF ( n_k /= 0 ) THEN
628	n_k2 = n_k * 2
629
630	!
631	!-- ( Eq.10 )^2
632	DO nf = -n_k2, n_k2
633	bdenom = bdenom + EXP( -qsi * pi * ABS(nf) / n_k )**2
634	ENDDO
635
636	!
637	!-- ( Eq.9 )
638	bdenom = SQRT( bdenom )
639	DO nf = -n_k2, n_k2
640	bxx(nf,k) = EXP( -qsi * pi * ABS(nf) / n_k ) / bdenom
641	ENDDO
642	ENDIF
643	ENDDO
644
645	END SUBROUTINE stg_filter_func
646
647
648	!------------------------------------------------------------------------------!
649	! Description:
650	! ------------
651	!> Parin for &stg_par for synthetic turbulence generator
652	!------------------------------------------------------------------------------!
653	SUBROUTINE stg_parin
654
655
656	IMPLICIT NONE
657
658	CHARACTER (LEN=80) :: line !< dummy string that contains the current line of the parameter file
659
660
661	NAMELIST /stg_par/ use_synthetic_turbulence_generator
662
663	line = ' '
664
665	!
666	!-- Try to find stg package
667	REWIND ( 11 )
668	line = ' '
669	DO WHILE ( INDEX( line, '&stg_par' ) == 0 )
670	READ ( 11, '(A)', END=10 ) line
671	ENDDO
672	BACKSPACE ( 11 )
673
674	!
675	!-- Read namelist
676	READ ( 11, stg_par )
677
678	!
679	!-- Set flag that indicates that the synthetic turbulence generator is switched
680	!-- on
681	synthetic_turbulence_generator = .TRUE.
682
683	IF ( TRIM( initializing_actions ) == 'read_restart_data' ) THEN
684	CALL stg_read_restart_data
685	ENDIF
686
687	10 CONTINUE
688
689	END SUBROUTINE stg_parin
690
691
692	!------------------------------------------------------------------------------!
693	! Description:
694	! ------------
695	!> Skipping the stg variables from restart-file (binary format).
696	!------------------------------------------------------------------------------!
697	SUBROUTINE stg_skip_var_list
698
699	IMPLICIT NONE
700
701	CHARACTER (LEN=1) :: cdum
702	CHARACTER (LEN=30) :: variable_chr
703
704	READ ( 13 ) variable_chr
705
706	DO WHILE ( TRIM( variable_chr ) /= '* end stg module *' )
707
708	READ ( 13 ) cdum
709	READ ( 13 ) variable_chr
710
711	ENDDO
712
713	END SUBROUTINE stg_skip_var_list
714
715
716	!------------------------------------------------------------------------------!
717	! Description:
718	! ------------
719	!> This routine reads the respective restart data.
720	!------------------------------------------------------------------------------!
721	SUBROUTINE stg_read_restart_data
722
723
724	IMPLICIT NONE
725
726	CHARACTER (LEN=30) :: variable_chr !< dummy variable to read string
727
728
729	READ ( 13 ) variable_chr
730	DO WHILE ( TRIM( variable_chr ) /= '* end stg module *' )
731
732	SELECT CASE ( TRIM( variable_chr ) )
733
734	CASE ( 'use_synthetic_turbulence_generator ' )
735	READ ( 13 ) use_synthetic_turbulence_generator
736	CASE ( 'mc_factor' )
737	READ ( 13 ) mc_factor
738
739	END SELECT
740
741	READ ( 13 ) variable_chr
742
743	ENDDO
744
745	END SUBROUTINE stg_read_restart_data
746
747
748	!------------------------------------------------------------------------------!
749	! Description:
750	! ------------
751	!> This routine writes the respective restart data.
752	!------------------------------------------------------------------------------!
753	SUBROUTINE stg_write_restart_data
754
755
756	IMPLICIT NONE
757
758	WRITE ( 14 ) 'use_synthetic_turbulence_generator '
759	WRITE ( 14 ) use_synthetic_turbulence_generator
760	WRITE ( 14 ) 'mc_factor '
761	WRITE ( 14 ) mc_factor
762
763	WRITE ( 14 ) '* end stg module * '
764
765	END SUBROUTINE stg_write_restart_data
766
767
768	!------------------------------------------------------------------------------!
769	! Description:
770	! ------------
771	!> Create turbulent inflow fields for u, v, w with prescribed length scales and
772	!> Reynolds stress tensor after a method of Xie and Castro (2008), modified
773	!> following suggestions of Kim et al. (2013), and using a Lund rotation
774	!> (Lund, 1998).
775	!------------------------------------------------------------------------------!
776	SUBROUTINE stg_main
777
778
779	USE arrays_3d, &
780	ONLY: dzw
781
782	USE control_parameters, &
783	ONLY: dt_3d, intermediate_timestep_count, simulated_time, &
784	volume_flow_initial
785
786	USE indices, &
787	ONLY: nyn, nys
788
789	USE statistics, &
790	ONLY: weight_substep
791
792
793	IMPLICIT NONE
794
795	INTEGER(iwp) :: j !< loop index in y-direction
796	INTEGER(iwp) :: k !< loop index in z-direction
797
798	REAL(wp) :: dt_stg !< wheighted subtimestep
799	REAL(wp) :: mc_factor_l !< local mass flux correction factor
800
801	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,5,nbgp) :: inflow_dist !< disturbance for inflow profiles
802
803
804	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'start' )
805
806	!
807	!-- Calculate time step which is needed for filter functions
808	dt_stg = dt_3d * weight_substep(intermediate_timestep_count)
809
810	!
811	!-- Initial value of fu, fv, fw
812	IF ( simulated_time == 0.0_wp .AND. .NOT. velocity_seed_initialized ) THEN
813	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fu )
814	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fv )
815	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fw )
816	velocity_seed_initialized = .TRUE.
817	ENDIF
818	!
819	!-- New set of fu, fv, fw
820	CALL stg_generate_seed_yz( nuy, nuz, buy, buz, fuo )
821	CALL stg_generate_seed_yz( nvy, nvz, bvy, bvz, fvo )
822	CALL stg_generate_seed_yz( nwy, nwz, bwy, bwz, fwo )
823
824	IF ( myidx == id_inflow ) THEN
825
826	DO j = nysg, nyng
827	DO k = nzb, nzt + 1
828	!
829	!-- Update fu, fv, fw following Eq. 14 of Xie and Castro (2008)
830	IF ( tu(k) == 0.0_wp ) THEN
831	fu(k,j) = fuo(k,j)
832	ELSE
833	fu(k,j) = fu(k,j) * EXP( -pi * dt_stg * 0.5_wp / tu(k) ) + &
834	fuo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tu(k) ) )
835	ENDIF
836
837	IF ( tv(k) == 0.0_wp ) THEN
838	fv(k,j) = fvo(k,j)
839	ELSE
840	fv(k,j) = fv(k,j) * EXP( -pi * dt_stg * 0.5_wp / tv(k) ) + &
841	fvo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tv(k) ) )
842	ENDIF
843
844	IF ( tw(k) == 0.0_wp ) THEN
845	fw(k,j) = fwo(k,j)
846	ELSE
847	fw(k,j) = fw(k,j) * EXP( -pi * dt_stg * 0.5_wp / tw(k) ) + &
848	fwo(k,j) * SQRT( 1.0_wp - EXP( -pi * dt_stg / tw(k) ) )
849	ENDIF
850	!
851	!-- Lund rotation following Eq. 17 in Xie and Castro (2008).
852	!-- Additional factors are added to improve the variance of v and w
853	IF( k == 0 ) THEN
854	inflow_dist(k,j,1,1) = 0.0_wp
855	inflow_dist(k,j,2,1) = 0.0_wp
856	inflow_dist(k,j,3,1) = 0.0_wp
857	inflow_dist(k,j,4,1) = 0.0_wp
858	inflow_dist(k,j,5,1) = 0.0_wp
859	ELSE
860	inflow_dist(k,j,1,1) = a11(k) * fu(k,j)
861	!experimental test of 1.2
862	inflow_dist(k,j,2,1) = ( SQRT( a22(k) / MAXVAL(a22) ) &
863	* 1.2_wp ) &
864	* ( a21(k) * fu(k,j) &
865	+ a22(k) * fv(k,j) )
866	inflow_dist(k,j,3,1) = ( SQRT(a33(k) / MAXVAL(a33) ) &
867	* 1.3_wp ) &
868	* ( a31(k) * fu(k,j) &
869	+ a32(k) * fv(k,j) &
870	+ a33(k) * fw(k,j) )
871	! Calculation for pt and e not yet implemented
872	inflow_dist(k,j,4,1) = 0.0_wp
873	inflow_dist(k,j,5,1) = 0.0_wp
874	ENDIF
875
876	ENDDO
877	ENDDO
878
879	!
880	!-- Mass flux correction following Kim et al. (2013)
881	!-- This correction factor insures that the mass flux is preserved at the
882	!-- inflow boundary
883	mc_factor_l = 0.0_wp
884	mc_factor = 0.0_wp
885	DO j = nys, nyn
886	DO k = nzb+1, nzt
887	mc_factor_l = mc_factor_l + dzw(k) * &
888	( mean_inflow_profiles(k,1) + inflow_dist(k,j,1,1) )
889	ENDDO
890	ENDDO
891
892	#if defined( __parallel )
893	CALL MPI_ALLREDUCE( mc_factor_l, mc_factor, &
894	1, MPI_REAL, MPI_SUM, comm1dy, ierr )
895	#else
896	mc_factor = mc_factor_l
897	#endif
898
899	mc_factor = volume_flow_initial(1) / mc_factor
900
901	!
902	!-- Add disturbance at the inflow
903	DO j = nysg, nyng
904	DO k = nzb, nzt+1
905	u(k,j,-nbgp+1:0) = ( mean_inflow_profiles(k,1) + &
906	inflow_dist(k,j,1,1) ) * mc_factor
907	v(k,j,-nbgp:-1) = ( mean_inflow_profiles(k,2) + &
908	inflow_dist(k,j,2,1) ) * mc_factor
909	w(k,j,-nbgp:-1) = inflow_dist(k,j,3,1) * mc_factor
910	ENDDO
911	ENDDO
912
913	ENDIF
914
915	CALL cpu_log( log_point(57), 'synthetic_turbulence_gen', 'stop' )
916
917	END SUBROUTINE stg_main
918
919
920	!------------------------------------------------------------------------------!
921	! Description:
922	! ------------
923	!> Generate a set of random number rand_it wich is equal on each PE
924	!> and calculate the velocity seed f_n.
925	!> f_n is splitted in vertical direction by the number of PEs in x-direction and
926	!> and each PE calculates a vertical subsection of f_n. At the the end, all
927	!> parts are collected to form the full array.
928	!------------------------------------------------------------------------------!
929	SUBROUTINE stg_generate_seed_yz( n_y, n_z, b_y, b_z, f_n )
930
931
932	USE indices, &
933	ONLY: ny
934
935
936	IMPLICIT NONE
937
938	INTEGER(iwp) :: j !< loop index in y-direction
939	INTEGER(iwp) :: jj !< loop index in y-direction
940	INTEGER(iwp) :: k !< loop index in z-direction
941	INTEGER(iwp) :: kk !< loop index in z-direction
942	INTEGER(iwp) :: send_count !< send count for MPI_GATHERV
943
944	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y !< length scale in y-direction
945	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z !< length scale in z-direction
946	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_y2 !< n_y*2
947	INTEGER(iwp), DIMENSION(nzb:nzt+1) :: n_z2 !< n_z*2
948
949	REAL(wp) :: nyz_inv !< inverse of number of grid points in yz-slice
950	REAL(wp) :: rand_av !< average of random number
951	REAL(wp) :: rand_sigma_inv !< inverse of stdev of random number
952
953	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_y !< filter func in y-dir
954	REAL(wp), DIMENSION(-merg:merg,nzb:nzt+1) :: b_z !< filter func in z-dir
955	REAL(wp), DIMENSION(nzb_x:nzt_x+1,nysg:nyng) :: f_n_l !< local velocity seed
956	REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng) :: f_n !< velocity seed
957	REAL(wp), DIMENSION(:,:), ALLOCATABLE :: rand_it !< random number
958
959
960	!
961	!-- Generate random numbers using a seed generated in stg_init.
962	!-- The set of random numbers are modified to have an average of 0 and
963	!-- unit variance.
964	ALLOCATE( rand_it(nzb-mergp:nzt+1+mergp,-mergp:ny+mergp) )
965
966	rand_av = 0.0_wp
967	rand_sigma_inv = 0.0_wp
968	nyz_inv = 1.0_wp / REAL( ( nzt+1 - nzb+1 ) * ( ny+1 ), KIND=wp )
969
970	DO j = 0, ny
971	DO k = nzb, nzt+1
972	CALL RANDOM_NUMBER( rand_it(k,j) )
973	rand_av = rand_av + rand_it(k,j)
974	ENDDO
975	ENDDO
976
977	rand_av = rand_av * nyz_inv
978
979	DO j = 0, ny
980	DO k = nzb, nzt+1
981	rand_it(k,j) = rand_it(k,j) - rand_av
982	rand_sigma_inv = rand_sigma_inv + rand_it(k,j) ** 2
983	ENDDO
984	ENDDO
985
986	rand_sigma_inv = 1.0_wp / SQRT(rand_sigma_inv * nyz_inv)
987
988	DO j = 0, ny
989	DO k = nzb, nzt+1
990	rand_it(k,j) = rand_it(k,j) * rand_sigma_inv
991	ENDDO
992	ENDDO
993
994	!
995	!-- Periodic fill of random number in space
996	DO j = 0, ny
997	DO k = 1, mergp
998	rand_it(nzb -k,j) = rand_it(nzt+2-k,j) ! bottom margin
999	rand_it(nzt+1+k,j) = rand_it(nzb+k-1,j) ! top margin
1000	ENDDO
1001	ENDDO
1002	DO j = 1, mergp
1003	DO k = nzb-mergp, nzt+1+mergp
1004	rand_it(k, -j) = rand_it(k,ny-j+1) ! south margin
1005	rand_it(k,ny+j) = rand_it(k, j-1) ! north margin
1006	ENDDO
1007	ENDDO
1008
1009	!
1010	!-- Generate velocity seed following Eq.6 of Xie and Castro (2008)
1011	n_y2 = n_y * 2
1012	n_z2 = n_z * 2
1013	f_n_l = 0.0_wp
1014
1015	DO j = nysg, nyng
1016	DO k = nzb_x, nzt_x+1
1017	DO jj = -n_y2(k), n_y2(k)
1018	DO kk = -n_z2(k), n_z2(k)
1019	f_n_l(k,j) = f_n_l(k,j) &
1020	+ b_y(jj,k) * b_z(kk,k) * rand_it(k+kk,j+jj)
1021	ENDDO
1022	ENDDO
1023	ENDDO
1024	ENDDO
1025
1026	DEALLOCATE( rand_it )
1027
1028	!
1029	!-- Gather velocity seeds of full subdomain
1030	send_count = nzt_x - nzb_x + 1
1031	IF ( nzt_x == nzt ) send_count = send_count + 1
1032
1033	#if defined( __parallel )
1034	CALL MPI_GATHERV( f_n_l(nzb_x,nysg), send_count, stg_type_yz_small, &
1035	f_n(nzb+1,nysg), recv_count, displs, stg_type_yz, &
1036	id_inflow, comm1dx, ierr )
1037	#else
1038	f_n(nzb+1:nzt+1,nysg:nyng) = f_n_l(nzb_x:nzt_x+1,nysg:nyng)
1039	#endif
1040
1041
1042	END SUBROUTINE stg_generate_seed_yz
1043
1044	END MODULE synthetic_turbulence_generator_mod

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |