Home

Context Navigation

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

Last change on this file since 2636 was 2576, checked in by Giersch, 7 years ago
Bugfixes for restart runs
Property svn:keywords set to `Id`
File size: 35.2 KB

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