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
|
! Calculation of characteristic upwinded differences except at boundaries.
! $Header$
! If this is an active cell we have to compute its derivative.
if ( eh_mask(i,j,k) .ge. 0 ) then
! If it is not a boundary cell in the x-direction then calculate as a
! starting point the centered derivative in the x-direction.
if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
( .not. btest(eh_mask(i,j,k),1) ) ) then
dfx(i,j,k) = idx * ( f(i+1,j,k) - f(i-1,j,k) )
! Otherwise we will have to calculate the appropriate one sided derivative.
else if ( btest(eh_mask(i,j,k),0) ) then
dfx(i,j,k) = idx * ( - three * f(i,j,k) + four * f(i+1,j,k) - f(i+2,j,k) )
else
dfx(i,j,k) = idx * ( three * f(i,j,k) - four * f(i-1,j,k) + f(i-2,j,k) )
end if
! Do the same for the y direction.
if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
( .not. btest(eh_mask(i,j,k),3) ) ) then
dfy(i,j,k) = idy * ( f(i,j+1,k) - f(i,j-1,k) )
else if ( btest(eh_mask(i,j,k),2) ) then
dfy(i,j,k) = idy * ( - three * f(i,j,k) + four * f(i,j+1,k) - f(i,j+2,k) )
else
dfy(i,j,k) = idy * ( three * f(i,j,k) - four * f(i,j-1,k) + f(i,j-2,k) )
end if
! And for the z direction.
if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
( .not. btest(eh_mask(i,j,k),5) ) ) then
dfz(i,j,k) = idz * ( f(i,j,k+1) - f(i,j,k-1) )
else if ( btest(eh_mask(i,j,k),4) ) then
dfz(i,j,k) = idz * ( - three * f(i,j,k) + four * f(i,j,k+1) - f(i,j,k+2) )
else
dfz(i,j,k) = idz * ( three * f(i,j,k) - four * f(i,j,k-1) + f(i,j,k-2) )
end if
! Calculate the characteristic direction using these preliminary
! derivatives and multiply it with the timestep.
dfup(1) = g3xx * dfx(i,j,k) + g3xy * dfy(i,j,k) + g3xz * dfz(i,j,k)
dfup(2) = g3xy * dfx(i,j,k) + g3yy * dfy(i,j,k) + g3yz * dfz(i,j,k)
dfup(3) = g3xz * dfx(i,j,k) + g3yz * dfy(i,j,k) + g3zz * dfz(i,j,k)
cfactor = one / sqrt ( alp2 * ( dfup(1) * dfx(i,j,k) + &
dfup(2) * dfy(i,j,k) + &
dfup(3) * dfz(i,j,k) ) )
cdx(1) = ( -betax(i,j,k) + alp2 * dfup(1) * cfactor ) * cctk_delta_time
cdx(2) = ( -betay(i,j,k) + alp2 * dfup(2) * cfactor ) * cctk_delta_time
cdx(3) = ( -betaz(i,j,k) + alp2 * dfup(3) * cfactor ) * cctk_delta_time
if ( ( .not. btest(eh_mask(i,j,k),0) ) .and. &
( .not. btest(eh_mask(i,j,k),1) ) ) then
! Look at the characteristic direction to figure out in which direction
! the stencil should be chosen.
if ( cdx(1) .gt. zero ) then ! Choose the stencil to the left.
! If the neighbour to the left is also an active non-boundary cell its
! safe to assume that i-2, i-1,i and i+1 have good values.
! Otherwise the neigbour must be a boundary cell, so we can only use a
! centered derivative, which we have already computed.
if ( eh_mask(i-1,j,k) .eq. 0 ) then
! Choose the stencil, that gives the smalles second derivative.
! If it happens to be the centered, do not do anything since we have
! already calculated the centered derivative.
! if ( f(i-2,j,k) - two * f(i-1,j,k) + f(i,j,k) .le. &
! f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) ) then
dfx(i,j,k) = idx * ( three * f(i,j,k) - &
four * f(i-1,j,k) + f(i-2,j,k) )
! end if
end if
else if ( cdx(1) .lt. zero ) then ! Choose the stencil to the right.
! If the neighbour to the right is also an active non-boundary cell its
! safe to assume that i-1, i,i+1 and i+2 have good values.
! Otherwise the neigbour must be a boundary cell, so we can only use a
! centered derivative, which we have already computed.
if ( eh_mask(i+1,j,k) .eq. 0 ) then
! Choose the stencil, that gives the smalles second derivative.
! If it happens to be the centered, do not do anything since we have
! already calculated the centered derivative.
! if ( f(i-1,j,k) - two * f(i,j,k) + f(i+1,j,k) .ge. &
! f(i,j,k) - two * f(i+1,j,k) + f(i+2,j,k) ) then
dfx(i,j,k) = idx * ( -three * f(i,j,k) + &
four * f(i+1,j,k) - f(i+2,j,k) )
! end if
end if
end if
end if
! Do the same for the y direction.
if ( ( .not. btest(eh_mask(i,j,k),2) ) .and. &
( .not. btest(eh_mask(i,j,k),3) ) ) then
if ( cdx(2) .gt. zero ) then ! Choose the stencil to the left.
if ( eh_mask(i,j-1,k) .eq. 0 ) then
! if ( f(i,j-2,k) - two * f(i,j-1,k) + f(i,j,k) .le. &
! f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) ) then
dfy(i,j,k) = idy * ( three * f(i,j,k) - &
four * f(i,j-1,k) + f(i,j-2,k) )
! end if
end if
else if ( cdx(2) .lt. zero ) then ! Choose the stencil to the left.
if ( eh_mask(i,j+1,k) .eq. 0 ) then
! if ( f(i,j-1,k) - two * f(i,j,k) + f(i,j+1,k) .ge. &
! f(i,j,k) - two * f(i,j+1,k) + f(i,j+2,k) ) then
dfy(i,j,k) = idy * ( -three * f(i,j,k) + &
four * f(i,j+1,k) - f(i,j+2,k) )
! end if
end if
end if
end if
! And for the z direction.
if ( ( .not. btest(eh_mask(i,j,k),4) ) .and. &
( .not. btest(eh_mask(i,j,k),5) ) ) then
if ( cdx(3) .gt. zero ) then
if ( eh_mask(i,j,k-1) .eq. 0 ) then
! if ( f(i,j,k-2) - two * f(i,j,k-1) + f(i,j,k) .le. &
! f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) ) then
dfz(i,j,k) = idz * ( three * f(i,j,k) - &
four * f(i,j,k-1) + f(i,j,k-2) )
! end if
end if
else if ( cdx(3) .lt. zero ) then
if ( eh_mask(i,j,k+1) .eq. 0 ) then
! if ( f(i,j,k-1) - two * f(i,j,k) + f(i,j,k+1) .ge. &
! f(i,j,k) - two * f(i,j,k+1) + f(i,j,k+2) ) then
dfz(i,j,k) = idz * ( -three * f(i,j,k) + &
four * f(i,j,k+1) - f(i,j,k+2) )
! end if
end if
end if
end if
else
dfx(i,j,k) = zero
dfy(i,j,k) = zero
dfz(i,j,k) = zero
end if
|
