-
Notifications
You must be signed in to change notification settings - Fork 0
/
derivs.f03
123 lines (111 loc) · 4.29 KB
/
derivs.f03
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
module derivs
use workspace
contains
complex*16 function laplacian(field_in, i, j, k)
use params
implicit none
class(fluid_field) :: field_in
integer :: i, j, k
laplacian = (-6.0d0*field_in%GRID(i, j, k) + BC(field_in, i + 1, j, k) + BC(field_in, i - 1, j, k) &
+ BC(field_in, i, j + 1, k) + BC(field_in, i, j - 1, k) + BC(field_in, i, j, k + 1) &
+ BC(field_in, i, j, k - 1))/(DSPACE**2.0d0)
end function
COMPLEX*16 function ddx(field_in, i, j, k)
use params
implicit none
integer :: i, j, k
class(fluid_field) :: field_in
ddx = (BC(field_in, i + 1, j, k) - BC(field_in, i - 1, j, k))/(2.0d0*DSPACE)
end function
COMPLEX*16 function ddy(field_in, i, j, k)
use params
implicit none
integer :: i, j, k
class(fluid_field) :: field_in
ddy = (BC(field_in, i, j + 1, k) - BC(field_in, i, j - 1, k))/(2.0d0*DSPACE)
end function
COMPLEX*16 function ddz(field_in, i, j, k)
use params
implicit none
class(fluid_field) :: field_in
integer :: i, j, k
ddz = (BC(field_in, i, j, k + 1) - BC(field_in, i, j, k - 1))/(2.0d0*DSPACE)
end function
COMPLEX*16 function d2dx2(field_in, i, j, k)
use params
implicit none
integer :: i, j, k
class(fluid_field) :: field_in
d2dx2 = (-2.0d0*field_in%GRID(i, j, k) + BC(field_in, i + 1, j, k) + BC(field_in, i - 1, j, k))/(DSPACE**2.0d0)
end function
COMPLEX*16 function d2dy2(field_in, i, j, k)
use params
implicit none
integer :: i, j, k
class(fluid_field) :: field_in
d2dy2 = (-2.0d0*field_in%GRID(i, j, k) + BC(field_in, i, j + 1, k) + BC(field_in, i, j - 1, k))/(DSPACE**2.0d0)
end function
COMPLEX*16 function d2dz2(field_in, i, j, k)
use params
implicit none
integer :: i, j, k
class(fluid_field) :: field_in
d2dz2 = (-2.0d0*field_in%GRID(i, j, k) + BC(field_in, i, j, k + 1) + BC(field_in, i, j, k - 1))/(DSPACE**2.0d0)
end function
complex*16 function BC(field_in, i, j, k)
use params
implicit none
class(fluid_field) :: field_in
integer :: i, j, k, ii, jj, kk
!Zero
if ((i >= NX .or. i <= 1) .and. BCX == 2) then
BC = 0.0d0
RETURN
else if ((j >= NY .or. j <= 1) .and. BCY == 2) then
BC = 0.0d0
RETURN
else if ((k >= NZ .or. k <= 1) .and. BCZ == 2) then
BC = 0.0d0
RETURN
end if
!Periodic for free due periodic ghost points from MPI
ii = i
jj = j
kk = k
!Reflective
if (i == NX + 1 .and. BCX == 0) ii = NX
if (i == 0 .and. BCX == 0) ii = 1
if (j == NY + 1 .and. BCY == 0) jj = NY
if (j == 0 .and. BCY == 0) jj = 1
if (k == NZ + 1 .and. BCZ == 0) kk = NZ
if (k == 0 .and. BCZ == 0) kk = 1
BC = field_in%GRID(ii, jj, kk)
!Quasi-periodic fluid field
if (i == NX + 1 .and. BCX == 3) then
BC = BC*exp(EYE*PI*NVORTALL(3, field_in%field_number)*GY(jj)/((NY - 1)*DSPACE) &
- EYE*PI*NVORTALL(2, field_in%field_number)*GZ(kk)/((NZ - 1)*DSPACE))
end if
if (i == 0 .and. BCX == 3) then
BC = BC/exp(EYE*PI*NVORTALL(3, field_in%field_number)*GY(jj)/((NY - 1)*DSPACE) &
- EYE*PI*NVORTALL(2, field_in%field_number)*GZ(kk)/((NZ - 1)*DSPACE))
end if
if (j == NY + 1 .and. BCY == 3) then
BC = BC*exp(-EYE*PI*NVORTALL(3, field_in%field_number)*GX(ii)/((NX - 1)*DSPACE) &
+ EYE*PI*NVORTALL(1, field_in%field_number)*GZ(kk)/((NZ - 1)*DSPACE))
end if
if (j == 0 .and. BCY == 3) then
BC = BC/exp(-EYE*PI*NVORTALL(3, field_in%field_number)*GX(ii)/((NX - 1)*DSPACE) &
+ EYE*PI*NVORTALL(1, field_in%field_number)*GZ(kk)/((NZ - 1)*DSPACE))
end if
if (k == NZ + 1 .and. BCZ == 3) then
BC = BC*exp(EYE*PI*NVORTALL(2, field_in%field_number)*GX(ii)/((NX - 1)*DSPACE) &
- EYE*PI*NVORTALL(1, field_in%field_number)*GY(jj)/((NY - 1)*DSPACE) &
+ EYE*PI*NVORTALL(1, field_in%field_number)*NVORTALL(2, field_in%field_number))
end if
if (k == 0 .and. BCZ == 3) then
BC = BC/exp(EYE*PI*NVORTALL(2, field_in%field_number)*GX(ii)/((NX - 1)*DSPACE) &
- EYE*PI*NVORTALL(1, field_in%field_number)*GY(jj)/((NY - 1)*DSPACE) &
+ EYE*PI*NVORTALL(1, field_in%field_number)*NVORTALL(2, field_in%field_number))
end if
end function
end module