1 c---------------------------------------------------------------------
2 c---------------------------------------------------------------------
6 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
9 c---------------------------------------------------------------------
10 c This subroutine initializes the field variable u using
11 c tri-linear transfinite interpolation of the boundary values
12 c---------------------------------------------------------------------
16 integer c, i, j, k, m, ii, jj, kk, ix, iy, iz
17 double precision xi, eta, zeta, Pface(5,3,2), Pxi, Peta,
20 c---------------------------------------------------------------------
21 c Later (in compute_rhs) we compute 1/u for every element. A few of
22 c the corner elements are not used, but it convenient (and faster)
23 c to compute the whole thing with a simple loop. Make sure those
24 c values are nonzero by initializing the whole thing here.
25 c---------------------------------------------------------------------
31 u(m, ii, jj, kk, c) = 1.0
37 c---------------------------------------------------------------------
41 c---------------------------------------------------------------------
42 c first store the "interpolated" values everywhere on the grid
43 c---------------------------------------------------------------------
46 do k = cell_low(3,c), cell_high(3,c)
47 zeta = dble(k) * dnzm1
49 do j = cell_low(2,c), cell_high(2,c)
52 do i = cell_low(1,c), cell_high(1,c)
56 call exact_solution(dble(ix-1), eta, zeta,
61 call exact_solution(xi, dble(iy-1) , zeta,
66 call exact_solution(xi, eta, dble(iz-1),
71 Pxi = xi * Pface(m,1,2) +
72 > (1.0d0-xi) * Pface(m,1,1)
73 Peta = eta * Pface(m,2,2) +
74 > (1.0d0-eta) * Pface(m,2,1)
75 Pzeta = zeta * Pface(m,3,2) +
76 > (1.0d0-zeta) * Pface(m,3,1)
78 u(m,ii,jj,kk,c) = Pxi + Peta + Pzeta -
79 > Pxi*Peta - Pxi*Pzeta - Peta*Pzeta +
91 c---------------------------------------------------------------------
92 c now store the exact values on the boundaries
93 c---------------------------------------------------------------------
95 c---------------------------------------------------------------------
97 c---------------------------------------------------------------------
102 do k = cell_low(3,c), cell_high(3,c)
103 zeta = dble(k) * dnzm1
105 do j = cell_low(2,c), cell_high(2,c)
106 eta = dble(j) * dnym1
107 call exact_solution(xi, eta, zeta, temp)
109 u(m,ii,jj,kk,c) = temp(m)
116 c---------------------------------------------------------------------
118 c---------------------------------------------------------------------
120 ii = cell_size(1,c)-1
123 do k = cell_low(3,c), cell_high(3,c)
124 zeta = dble(k) * dnzm1
126 do j = cell_low(2,c), cell_high(2,c)
127 eta = dble(j) * dnym1
128 call exact_solution(xi, eta, zeta, temp)
130 u(m,ii,jj,kk,c) = temp(m)
137 c---------------------------------------------------------------------
139 c---------------------------------------------------------------------
144 do k = cell_low(3,c), cell_high(3,c)
145 zeta = dble(k) * dnzm1
147 do i = cell_low(1,c), cell_high(1,c)
149 call exact_solution(xi, eta, zeta, temp)
151 u(m,ii,jj,kk,c) = temp(m)
159 c---------------------------------------------------------------------
161 c---------------------------------------------------------------------
163 jj = cell_size(2,c)-1
166 do k = cell_low(3,c), cell_high(3,c)
167 zeta = dble(k) * dnzm1
169 do i = cell_low(1,c), cell_high(1,c)
171 call exact_solution(xi, eta, zeta, temp)
173 u(m,ii,jj,kk,c) = temp(m)
180 c---------------------------------------------------------------------
182 c---------------------------------------------------------------------
187 do j = cell_low(2,c), cell_high(2,c)
188 eta = dble(j) * dnym1
190 do i =cell_low(1,c), cell_high(1,c)
192 call exact_solution(xi, eta, zeta, temp)
194 u(m,ii,jj,kk,c) = temp(m)
201 c---------------------------------------------------------------------
203 c---------------------------------------------------------------------
205 kk = cell_size(3,c)-1
208 do j = cell_low(2,c), cell_high(2,c)
209 eta = dble(j) * dnym1
211 do i =cell_low(1,c), cell_high(1,c)
213 call exact_solution(xi, eta, zeta, temp)
215 u(m,ii,jj,kk,c) = temp(m)
226 c---------------------------------------------------------------------
227 c---------------------------------------------------------------------
231 c---------------------------------------------------------------------
232 c---------------------------------------------------------------------
236 integer i, j, k, d, c, m, n
238 c---------------------------------------------------------------------
239 c loop over all cells
240 c---------------------------------------------------------------------
243 c---------------------------------------------------------------------
244 c first, initialize the start and end arrays
245 c---------------------------------------------------------------------
247 if (cell_coord(d,c) .eq. 1) then
252 if (cell_coord(d,c) .eq. ncells) then
259 c---------------------------------------------------------------------
260 c zero the whole left hand side for starters
261 c---------------------------------------------------------------------
262 do k = 0, cell_size(3,c)-1
263 do j = 0, cell_size(2,c)-1
264 do i = 0, cell_size(1,c)-1
267 lhsc(m,n,i,j,k,c) = 0.0d0
280 c---------------------------------------------------------------------
281 c---------------------------------------------------------------------
283 subroutine lhsabinit(lhsa, lhsb, size)
287 double precision lhsa(5, 5, -1:size), lhsb(5, 5, -1:size)
291 c---------------------------------------------------------------------
292 c next, set all diagonal values to 1. This is overkill, but convenient
293 c---------------------------------------------------------------------