Home

Context Navigation

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

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