2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
10 c---------------------------------------------------------------------
11 c compute the right hand side based on exact solution
12 c---------------------------------------------------------------------
16 double precision dtemp(5), xi, eta, zeta, dtpp
17 integer c, m, i, j, k, ip1, im1, jp1,
20 c---------------------------------------------------------------------
21 c loop over all cells owned by this node
22 c---------------------------------------------------------------------
25 c---------------------------------------------------------------------
27 c---------------------------------------------------------------------
29 do k= 0, cell_size(3,c)-1
30 do j = 0, cell_size(2,c)-1
31 do i = 0, cell_size(1,c)-1
32 forcing(i,j,k,m,c) = 0.0d0
38 c---------------------------------------------------------------------
39 c xi-direction flux differences
40 c---------------------------------------------------------------------
41 do k = start(3,c), cell_size(3,c)-end(3,c)-1
42 zeta = dble(k+cell_low(3,c)) * dnzm1
43 do j = start(2,c), cell_size(2,c)-end(2,c)-1
44 eta = dble(j+cell_low(2,c)) * dnym1
46 do i=-2*(1-start(1,c)), cell_size(1,c)+1-2*end(1,c)
47 xi = dble(i+cell_low(1,c)) * dnxm1
49 call exact_solution(xi, eta, zeta, dtemp)
54 dtpp = 1.0d0 / dtemp(1)
57 buf(i,m) = dtpp * dtemp(m)
60 cuf(i) = buf(i,2) * buf(i,2)
61 buf(i,1) = cuf(i) + buf(i,3) * buf(i,3) +
63 q(i) = 0.5d0*(buf(i,2)*ue(i,2) + buf(i,3)*ue(i,3) +
68 do i = start(1,c), cell_size(1,c)-end(1,c)-1
72 forcing(i,j,k,1,c) = forcing(i,j,k,1,c) -
73 > tx2*( ue(ip1,2)-ue(im1,2) )+
74 > dx1tx1*(ue(ip1,1)-2.0d0*ue(i,1)+ue(im1,1))
76 forcing(i,j,k,2,c) = forcing(i,j,k,2,c) - tx2 * (
77 > (ue(ip1,2)*buf(ip1,2)+c2*(ue(ip1,5)-q(ip1)))-
78 > (ue(im1,2)*buf(im1,2)+c2*(ue(im1,5)-q(im1))))+
79 > xxcon1*(buf(ip1,2)-2.0d0*buf(i,2)+buf(im1,2))+
80 > dx2tx1*( ue(ip1,2)-2.0d0* ue(i,2)+ue(im1,2))
82 forcing(i,j,k,3,c) = forcing(i,j,k,3,c) - tx2 * (
83 > ue(ip1,3)*buf(ip1,2)-ue(im1,3)*buf(im1,2))+
84 > xxcon2*(buf(ip1,3)-2.0d0*buf(i,3)+buf(im1,3))+
85 > dx3tx1*( ue(ip1,3)-2.0d0*ue(i,3) +ue(im1,3))
87 forcing(i,j,k,4,c) = forcing(i,j,k,4,c) - tx2*(
88 > ue(ip1,4)*buf(ip1,2)-ue(im1,4)*buf(im1,2))+
89 > xxcon2*(buf(ip1,4)-2.0d0*buf(i,4)+buf(im1,4))+
90 > dx4tx1*( ue(ip1,4)-2.0d0* ue(i,4)+ ue(im1,4))
92 forcing(i,j,k,5,c) = forcing(i,j,k,5,c) - tx2*(
93 > buf(ip1,2)*(c1*ue(ip1,5)-c2*q(ip1))-
94 > buf(im1,2)*(c1*ue(im1,5)-c2*q(im1)))+
95 > 0.5d0*xxcon3*(buf(ip1,1)-2.0d0*buf(i,1)+
97 > xxcon4*(cuf(ip1)-2.0d0*cuf(i)+cuf(im1))+
98 > xxcon5*(buf(ip1,5)-2.0d0*buf(i,5)+buf(im1,5))+
99 > dx5tx1*( ue(ip1,5)-2.0d0* ue(i,5)+ ue(im1,5))
102 c---------------------------------------------------------------------
103 c Fourth-order dissipation
104 c---------------------------------------------------------------------
105 if (start(1,c) .gt. 0) then
108 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
109 > (5.0d0*ue(i,m) - 4.0d0*ue(i+1,m) +ue(i+2,m))
111 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
112 > (-4.0d0*ue(i-1,m) + 6.0d0*ue(i,m) -
113 > 4.0d0*ue(i+1,m) + ue(i+2,m))
118 do i = start(1,c)*3, cell_size(1,c)-3*end(1,c)-1
119 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp*
120 > (ue(i-2,m) - 4.0d0*ue(i-1,m) +
121 > 6.0d0*ue(i,m) - 4.0d0*ue(i+1,m) + ue(i+2,m))
125 if (end(1,c) .gt. 0) then
128 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
129 > (ue(i-2,m) - 4.0d0*ue(i-1,m) +
130 > 6.0d0*ue(i,m) - 4.0d0*ue(i+1,m))
132 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
133 > (ue(i-2,m) - 4.0d0*ue(i-1,m) + 5.0d0*ue(i,m))
139 c---------------------------------------------------------------------
140 c eta-direction flux differences
141 c---------------------------------------------------------------------
142 do k = start(3,c), cell_size(3,c)-end(3,c)-1
143 zeta = dble(k+cell_low(3,c)) * dnzm1
144 do i=start(1,c), cell_size(1,c)-end(1,c)-1
145 xi = dble(i+cell_low(1,c)) * dnxm1
147 do j=-2*(1-start(2,c)), cell_size(2,c)+1-2*end(2,c)
148 eta = dble(j+cell_low(2,c)) * dnym1
150 call exact_solution(xi, eta, zeta, dtemp)
154 dtpp = 1.0d0/dtemp(1)
157 buf(j,m) = dtpp * dtemp(m)
160 cuf(j) = buf(j,3) * buf(j,3)
161 buf(j,1) = cuf(j) + buf(j,2) * buf(j,2) +
162 > buf(j,4) * buf(j,4)
163 q(j) = 0.5d0*(buf(j,2)*ue(j,2) + buf(j,3)*ue(j,3) +
167 do j = start(2,c), cell_size(2,c)-end(2,c)-1
171 forcing(i,j,k,1,c) = forcing(i,j,k,1,c) -
172 > ty2*( ue(jp1,3)-ue(jm1,3) )+
173 > dy1ty1*(ue(jp1,1)-2.0d0*ue(j,1)+ue(jm1,1))
175 forcing(i,j,k,2,c) = forcing(i,j,k,2,c) - ty2*(
176 > ue(jp1,2)*buf(jp1,3)-ue(jm1,2)*buf(jm1,3))+
177 > yycon2*(buf(jp1,2)-2.0d0*buf(j,2)+buf(jm1,2))+
178 > dy2ty1*( ue(jp1,2)-2.0* ue(j,2)+ ue(jm1,2))
180 forcing(i,j,k,3,c) = forcing(i,j,k,3,c) - ty2*(
181 > (ue(jp1,3)*buf(jp1,3)+c2*(ue(jp1,5)-q(jp1)))-
182 > (ue(jm1,3)*buf(jm1,3)+c2*(ue(jm1,5)-q(jm1))))+
183 > yycon1*(buf(jp1,3)-2.0d0*buf(j,3)+buf(jm1,3))+
184 > dy3ty1*( ue(jp1,3)-2.0d0*ue(j,3) +ue(jm1,3))
186 forcing(i,j,k,4,c) = forcing(i,j,k,4,c) - ty2*(
187 > ue(jp1,4)*buf(jp1,3)-ue(jm1,4)*buf(jm1,3))+
188 > yycon2*(buf(jp1,4)-2.0d0*buf(j,4)+buf(jm1,4))+
189 > dy4ty1*( ue(jp1,4)-2.0d0*ue(j,4)+ ue(jm1,4))
191 forcing(i,j,k,5,c) = forcing(i,j,k,5,c) - ty2*(
192 > buf(jp1,3)*(c1*ue(jp1,5)-c2*q(jp1))-
193 > buf(jm1,3)*(c1*ue(jm1,5)-c2*q(jm1)))+
194 > 0.5d0*yycon3*(buf(jp1,1)-2.0d0*buf(j,1)+
196 > yycon4*(cuf(jp1)-2.0d0*cuf(j)+cuf(jm1))+
197 > yycon5*(buf(jp1,5)-2.0d0*buf(j,5)+buf(jm1,5))+
198 > dy5ty1*(ue(jp1,5)-2.0d0*ue(j,5)+ue(jm1,5))
201 c---------------------------------------------------------------------
202 c Fourth-order dissipation
203 c---------------------------------------------------------------------
204 if (start(2,c) .gt. 0) then
207 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
208 > (5.0d0*ue(j,m) - 4.0d0*ue(j+1,m) +ue(j+2,m))
210 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
211 > (-4.0d0*ue(j-1,m) + 6.0d0*ue(j,m) -
212 > 4.0d0*ue(j+1,m) + ue(j+2,m))
217 do j = start(2,c)*3, cell_size(2,c)-3*end(2,c)-1
218 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp*
219 > (ue(j-2,m) - 4.0d0*ue(j-1,m) +
220 > 6.0d0*ue(j,m) - 4.0d0*ue(j+1,m) + ue(j+2,m))
223 if (end(2,c) .gt. 0) then
226 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
227 > (ue(j-2,m) - 4.0d0*ue(j-1,m) +
228 > 6.0d0*ue(j,m) - 4.0d0*ue(j+1,m))
230 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
231 > (ue(j-2,m) - 4.0d0*ue(j-1,m) + 5.0d0*ue(j,m))
239 c---------------------------------------------------------------------
240 c zeta-direction flux differences
241 c---------------------------------------------------------------------
242 do j=start(2,c), cell_size(2,c)-end(2,c)-1
243 eta = dble(j+cell_low(2,c)) * dnym1
244 do i = start(1,c), cell_size(1,c)-end(1,c)-1
245 xi = dble(i+cell_low(1,c)) * dnxm1
247 do k=-2*(1-start(3,c)), cell_size(3,c)+1-2*end(3,c)
248 zeta = dble(k+cell_low(3,c)) * dnzm1
250 call exact_solution(xi, eta, zeta, dtemp)
255 dtpp = 1.0d0/dtemp(1)
258 buf(k,m) = dtpp * dtemp(m)
261 cuf(k) = buf(k,4) * buf(k,4)
262 buf(k,1) = cuf(k) + buf(k,2) * buf(k,2) +
263 > buf(k,3) * buf(k,3)
264 q(k) = 0.5d0*(buf(k,2)*ue(k,2) + buf(k,3)*ue(k,3) +
268 do k=start(3,c), cell_size(3,c)-end(3,c)-1
272 forcing(i,j,k,1,c) = forcing(i,j,k,1,c) -
273 > tz2*( ue(kp1,4)-ue(km1,4) )+
274 > dz1tz1*(ue(kp1,1)-2.0d0*ue(k,1)+ue(km1,1))
276 forcing(i,j,k,2,c) = forcing(i,j,k,2,c) - tz2 * (
277 > ue(kp1,2)*buf(kp1,4)-ue(km1,2)*buf(km1,4))+
278 > zzcon2*(buf(kp1,2)-2.0d0*buf(k,2)+buf(km1,2))+
279 > dz2tz1*( ue(kp1,2)-2.0d0* ue(k,2)+ ue(km1,2))
281 forcing(i,j,k,3,c) = forcing(i,j,k,3,c) - tz2 * (
282 > ue(kp1,3)*buf(kp1,4)-ue(km1,3)*buf(km1,4))+
283 > zzcon2*(buf(kp1,3)-2.0d0*buf(k,3)+buf(km1,3))+
284 > dz3tz1*(ue(kp1,3)-2.0d0*ue(k,3)+ue(km1,3))
286 forcing(i,j,k,4,c) = forcing(i,j,k,4,c) - tz2 * (
287 > (ue(kp1,4)*buf(kp1,4)+c2*(ue(kp1,5)-q(kp1)))-
288 > (ue(km1,4)*buf(km1,4)+c2*(ue(km1,5)-q(km1))))+
289 > zzcon1*(buf(kp1,4)-2.0d0*buf(k,4)+buf(km1,4))+
290 > dz4tz1*( ue(kp1,4)-2.0d0*ue(k,4) +ue(km1,4))
292 forcing(i,j,k,5,c) = forcing(i,j,k,5,c) - tz2 * (
293 > buf(kp1,4)*(c1*ue(kp1,5)-c2*q(kp1))-
294 > buf(km1,4)*(c1*ue(km1,5)-c2*q(km1)))+
295 > 0.5d0*zzcon3*(buf(kp1,1)-2.0d0*buf(k,1)
297 > zzcon4*(cuf(kp1)-2.0d0*cuf(k)+cuf(km1))+
298 > zzcon5*(buf(kp1,5)-2.0d0*buf(k,5)+buf(km1,5))+
299 > dz5tz1*( ue(kp1,5)-2.0d0*ue(k,5)+ ue(km1,5))
302 c---------------------------------------------------------------------
303 c Fourth-order dissipation
304 c---------------------------------------------------------------------
305 if (start(3,c) .gt. 0) then
308 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
309 > (5.0d0*ue(k,m) - 4.0d0*ue(k+1,m) +ue(k+2,m))
311 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
312 > (-4.0d0*ue(k-1,m) + 6.0d0*ue(k,m) -
313 > 4.0d0*ue(k+1,m) + ue(k+2,m))
318 do k = start(3,c)*3, cell_size(3,c)-3*end(3,c)-1
319 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp*
320 > (ue(k-2,m) - 4.0d0*ue(k-1,m) +
321 > 6.0d0*ue(k,m) - 4.0d0*ue(k+1,m) + ue(k+2,m))
325 if (end(3,c) .gt. 0) then
328 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
329 > (ue(k-2,m) - 4.0d0*ue(k-1,m) +
330 > 6.0d0*ue(k,m) - 4.0d0*ue(k+1,m))
332 forcing(i,j,k,m,c) = forcing(i,j,k,m,c) - dssp *
333 > (ue(k-2,m) - 4.0d0*ue(k-1,m) + 5.0d0*ue(k,m))
339 c---------------------------------------------------------------------
340 c now change the sign of the forcing function,
341 c---------------------------------------------------------------------
343 do k = start(3,c), cell_size(3,c)-end(3,c)-1
344 do j = start(2,c), cell_size(2,c)-end(2,c)-1
345 do i = start(1,c), cell_size(1,c)-end(1,c)-1
346 forcing(i,j,k,m,c) = -1.d0 * forcing(i,j,k,m,c)
352 c---------------------------------------------------------------------
354 c---------------------------------------------------------------------