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
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
|
! Calculate second order upwinded differences and store them in the dfdx, dfdy
! and dfdz grid functions.
! $Header$
! First figure out the range of indices excluding ghost cells and symmetry
! cells.
# include "physical_part.h"
! Set idx, idy and idz to 0.5*gridspacing
! idx, idy and idz must be defined in the including routine.
idx = half / cctk_delta_space(1)
idy = half / cctk_delta_space(2)
idz = half / cctk_delta_space(3)
! Calculate appropriate one sided derivatives for all active cells.
! Cells are active if eh_mask is greater or equal to zero. They are interior
! cells if eh_mask == 0 and boundary cells if eh_mask > 0.
do l = 1, eh_number_level_sets
do k = kzl, kzr
do j = jyl, jyr
do i = ixl, ixr
! If this is an active cell...
if ( eh_mask(i,j,k,l) .ge. 0 ) then
! If this is not a boundary cell...
if ( ( .not. btest(eh_mask(i,j,k,l),0) ) .and. &
( .not. btest(eh_mask(i,j,k,l),1) ) ) then
! If the cell to the left is not a boundary cell...
if ( ( .not. btest(eh_mask(i-1,j,k,l),0) ) .or. &
( eh_mask(i-1,j,k,l) .eq. -2 ) ) then
! Calculate left handed one sided derivative
al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
else
al = zero
end if
! If the cell to the right is not a boundary cell...
if ( ( .not. btest(eh_mask(i+1,j,k,l),1) ) .or. &
( eh_mask(i+1,j,k,l) .eq. -2 ) ) then
! Calculate right handed one sided derivative
ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
else
ar = zero
end if
! Else if the cell is a left boundary cell.
else if ( btest(eh_mask(i,j,k,l),0) ) then
! calculate only right handed one sided difference.
al = zero
ar = idx * ( - three * f(i,j,k,l) + four * f(i+1,j,k,l) - f(i+2,j,k,l) )
! Else it must be a right boundary cell.
else
! So calculate only the left handed one sided difference.
al = idx * ( three * f(i,j,k,l) - four * f(i-1,j,k,l) + f(i-2,j,k,l) )
ar = zero
end if
! Do the same for the y-derivative.
if ( ( .not. btest(eh_mask(i,j,k,l),2) ) .and. &
( .not. btest(eh_mask(i,j,k,l),3) ) ) then
if ( ( .not. btest(eh_mask(i,j-1,k,l),2) ) .or. &
( eh_mask(i,j-1,k,l) .eq. -2 ) ) then
bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
else
bl = zero
end if
if ( ( .not. btest(eh_mask(i,j+1,k,l),3) ) .or. &
( eh_mask(i,j+1,k,l) .eq. -2 ) ) then
br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
else
br = zero
end if
else if ( btest(eh_mask(i,j,k,l),2) ) then
bl = zero
br = idy * ( - three * f(i,j,k,l) + four * f(i,j+1,k,l) - f(i,j+2,k,l) )
else
bl = idy * ( three * f(i,j,k,l) - four * f(i,j-1,k,l) + f(i,j-2,k,l) )
br = zero
end if
! And the z-derivative
if ( ( .not. btest(eh_mask(i,j,k,l),4) ) .and. &
( .not. btest(eh_mask(i,j,k,l),5) ) ) then
if ( ( .not. btest(eh_mask(i,j,k-1,l),4) ) .or. &
( eh_mask(i,j,k-1,l) .eq. -2 ) ) then
cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
else
cl = zero
end if
if ( ( .not. btest(eh_mask(i,j,k+1,l),5) ) .or. &
( eh_mask(i,j,k+1,l) .eq. -2 ) ) then
cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
else
cr = zero
end if
else if ( btest(eh_mask(i,j,k,l),4) ) then
cl = zero
cr = idz * ( - three * f(i,j,k,l) + four * f(i,j,k+1,l) - f(i,j,k+2,l) )
else
cl = idz * ( three * f(i,j,k,l) - four * f(i,j,k-1,l) + f(i,j,k-2,l) )
cr = zero
end if
! Get the negative and positive parts of the one sided derivatives in
! the x-direction.
alminus = half*(al-abs(al))
alplus = half*(al+abs(al))
arminus = half*(ar-abs(ar))
arplus = half*(ar+abs(ar))
! Ditto for the y-direction.
blminus = half*(bl-abs(bl))
blplus = half*(bl+abs(bl))
brminus = half*(br-abs(br))
brplus = half*(br+abs(br))
! Ditto for the z-direction.
clminus = half*(cl-abs(cl))
clplus = half*(cl+abs(cl))
crminus = half*(cr-abs(cr))
crplus = half*(cr+abs(cr))
! If f>0 choose the correct one sided derivative depending on the
! magnitude of the negative and positive parts
if ( f(i,j,k,l) .gt. 0 ) then
if ( abs(alplus) .gt. abs(arminus) ) then
dfx(i,j,k,l) = alplus
else
dfx(i,j,k,l) = arminus
end if
if ( abs(blplus) .gt. abs(brminus) ) then
dfy(i,j,k,l) = blplus
else
dfy(i,j,k,l) = brminus
end if
if ( abs(clplus) .gt. abs(crminus) ) then
dfz(i,j,k,l) = clplus
else
dfz(i,j,k,l) = crminus
end if
endif
! Ditto if f<0.
if ( f(i,j,k,l) .lt. 0 ) then
if ( abs(alminus) .gt. abs(arplus) ) then
dfx(i,j,k,l) = alminus
else
dfx(i,j,k,l) = arplus
end if
if ( abs(blminus) .gt. abs(brplus) ) then
dfy(i,j,k,l) = blminus
else
dfy(i,j,k,l) = brplus
end if
if ( abs(clminus) .gt. abs(crplus) ) then
dfz(i,j,k,l) = clminus
else
dfz(i,j,k,l) = crplus
end if
end if
! If the cell is not active set the derivatives to zero.
else
dfx(i,j,k,l) = zero
dfy(i,j,k,l) = zero
dfz(i,j,k,l) = zero
end if
end do
end do
end do
end do
|
