-
Notifications
You must be signed in to change notification settings - Fork 2
/
LES_Filtering_read.f90
148 lines (117 loc) · 5.83 KB
/
LES_Filtering_read.f90
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
!------------------------------------------------------------------------------!
! !
! PROGRAM : LES_Filtering_read.f90 !
! !
! PURPOSE : Reading datas from the DNS data for turbulent channel flow. !
! !
! 2016.12.07 K.Noh !
! !
!------------------------------------------------------------------------------!
SUBROUTINE READ_DNS
USE LES_FILTERING_module, &
ONLY : G, FIND_dU, FIND_dx
USE LES_FILTERING_module, &
ONLY : Nx, Ny, Nz, dx, dz, Del, FW, tol, Nx_fil, Nz_fil, &
file_name, dir_name, path_name, Y_order
USE LES_FILTERING_module, &
ONLY : X, Y, Z, U, V, W, dy, S_T, O_T
IMPLICIT NONE
INTEGER :: i,j,k,it,x_i,x_j
REAL(KIND=8) :: tmp_x, tmp_y, tmp_z, time_sta, time_end, dg, r, &
dU_i, dU_j, dx_i, dx_j
CHARACTER(20) :: header
WRITE(*,*) '----------------------------------------------------'
WRITE(*,*) ' READING PROCESS STARTED '
CALL CPU_TIME(time_sta)
dir_name = 'Data'
path_name = TRIM(dir_name)//'/'//TRIM(file_name)
OPEN(100,FILE=path_name,FORM='FORMATTED',STATUS='OLD')
READ(100,*) header
READ(100,*) header
!--------------------------------------------------------------------!
! Main loop of reading DNS datas !
!--------------------------------------------------------------------!
IF ( Y_ORDER == 0 ) THEN
DO k = 1,Nz
DO j = 1,Ny
DO i = 1,Nx
READ(100,*) tmp_x, tmp_y, tmp_z, U(i,j,k), V(i,j,k), W(i,j,k)
IF (j==1 .AND. k==1) X(i) = tmp_x
IF (i==1 .AND. k==1) Y(j) = tmp_y
IF (i==1 .AND. j==1) Z(k) = tmp_z
END DO
END DO
END DO
ELSE
DO k = 1,Nz
DO j = Ny,1,-1
DO i = 1,Nx
READ(100,*) tmp_x, tmp_y, tmp_z, U(i,j,k), V(i,j,k), W(i,j,k)
IF (j==1 .AND. k==1) X(i) = tmp_x
IF (i==1 .AND. k==1) Y(j) = tmp_y
IF (i==1 .AND. j==1) Z(k) = tmp_z
END DO
END DO
END DO
END IF
CLOSE(100)
!--------------------------------------------------------------------!
! Calculating dx,dy,dz !
!--------------------------------------------------------------------!
DO j = 1,Ny-1
dy(j) = Y(j+1) - Y(j)
END DO
dx = X(2) - X(1)
dz = Z(2) - Z(1)
Del = FW*sqrt(dx * dz)
!--------------------------------------------------------------------!
! Calculating Stress Tensor !
!--------------------------------------------------------------------!
!$OMP PARALLEL DO private(k,j,i,x_i,x_j,dU_i,dU_j,dx_i,dx_j)
DO k = 2,Nz-1
DO j = 2,Ny-1
DO i = 2,Nx-1
DO x_i = 1,3
DO x_j = 1,3
dU_i = FIND_dU(i,j,k,x_i,x_j)
dU_j = FIND_dU(i,j,k,x_j,x_i)
dx_i = FIND_dx(i,j,k,x_i)
dx_j = FIND_dx(i,j,k,x_j)
S_T(i,j,k,x_i,x_j) = ( dU_i / dx_j + dU_j / dx_i ) / 2.0
O_T(i,j,k,x_i,x_j) = ( dU_i / dx_j - dU_j / dx_i ) / 2.0
END DO
END DO
END DO
END DO
END DO
!OMP END PARALLEL
S_T(1,1:Ny,1:Nz,1:3,1:3) = S_T(2,1:Ny,1:Nz,1:3,1:3)
S_T(Nx,1:Ny,1:Nz,1:3,1:3) = S_T(Nx-1,1:Ny,1:Nz,1:3,1:3)
S_T(1:Nx,1,1:Nz,1:3,1:3) = S_T(1:Nx,2,1:Nz,1:3,1:3)
S_T(1:Nx,Ny,1:Nz,1:3,1:3) = S_T(1:Nx,Ny-1,1:Nz,1:3,1:3)
S_T(1:Nx,1:Ny,1,1:3,1:3) = S_T(1:Nx,1:Ny,2,1:3,1:3)
S_T(1:Nx,1:Ny,Nz,1:3,1:3) = S_T(1:Nx,1:Ny,Nz-1,1:3,1:3)
O_T(1,1:Ny,1:Nz,1:3,1:3) = O_T(2,1:Ny,1:Nz,1:3,1:3)
O_T(Nx,1:Ny,1:Nz,1:3,1:3) = O_T(Nx-1,1:Ny,1:Nz,1:3,1:3)
O_T(1:Nx,1,1:Nz,1:3,1:3) = O_T(1:Nx,2,1:Nz,1:3,1:3)
O_T(1:Nx,Ny,1:Nz,1:3,1:3) = O_T(1:Nx,Ny-1,1:Nz,1:3,1:3)
O_T(1:Nx,1:Ny,1,1:3,1:3) = O_T(1:Nx,1:Ny,2,1:3,1:3)
O_T(1:Nx,1:Ny,Nz,1:3,1:3) = O_T(1:Nx,1:Ny,Nz-1,1:3,1:3)
!--------------------------------------------------------------------!
! Determining the number of filtering node in x,z dirctions !
!--------------------------------------------------------------------!
DO it = 1,1000
r = REAL(it,KIND=8)
dg = G(Del,r) - G(Del,r+1)
IF (dg < tol) THEN
Nx_fil = CEILING(REAL(it,KIND=8)/(2*dx))
Nz_fil = CEILING(REAL(it,KIND=8)/(2*dz))
EXIT
END IF
END DO
CALL CPU_TIME(time_end)
WRITE(*,*) ' READING PROCESS IS COMPLETED '
WRITE(*,*) ' Total Reading time : ',time_end - time_sta,' s'
WRITE(*,*) '----------------------------------------------------'
WRITE(*,*) ''
END SUBROUTINE READ_DNS