2 c---------------------------------------------------------------------
3 c---------------------------------------------------------------------
7 c---------------------------------------------------------------------
8 c---------------------------------------------------------------------
10 c---------------------------------------------------------------------
11 c This subroutine initializes the field variable u using
12 c tri-linear transfinite interpolation of the boundary values
13 c---------------------------------------------------------------------
17 integer c, i, j, k, m, ii, jj, kk, ix, iy, iz
18 double precision xi, eta, zeta, Pface(5,3,2), Pxi, Peta,
22 c---------------------------------------------------------------------
23 c Later (in compute_rhs) we compute 1/u for every element. A few of
24 c the corner elements are not used, but it convenient (and faster)
25 c to compute the whole thing with a simple loop. Make sure those
26 c values are nonzero by initializing the whole thing here.
27 c---------------------------------------------------------------------
32 u(ii, jj, kk, 1, c) = 1.0
33 u(ii, jj, kk, 2, c) = 0.0
34 u(ii, jj, kk, 3, c) = 0.0
35 u(ii, jj, kk, 4, c) = 0.0
36 u(ii, jj, kk, 5, c) = 1.0
42 c---------------------------------------------------------------------
43 c first store the "interpolated" values everywhere on the grid
44 c---------------------------------------------------------------------
47 do k = cell_low(3,c), cell_high(3,c)
48 zeta = dble(k) * dnzm1
50 do j = cell_low(2,c), cell_high(2,c)
53 do i = cell_low(1,c), cell_high(1,c)
57 call exact_solution(dble(ix-1), eta, zeta,
62 call exact_solution(xi, dble(iy-1) , zeta,
67 call exact_solution(xi, eta, dble(iz-1),
72 Pxi = xi * Pface(m,1,2) +
73 > (1.0d0-xi) * Pface(m,1,1)
74 Peta = eta * Pface(m,2,2) +
75 > (1.0d0-eta) * Pface(m,2,1)
76 Pzeta = zeta * Pface(m,3,2) +
77 > (1.0d0-zeta) * Pface(m,3,1)
79 u(ii,jj,kk,m,c) = Pxi + Peta + Pzeta -
80 > Pxi*Peta - Pxi*Pzeta - Peta*Pzeta +
92 c---------------------------------------------------------------------
93 c now store the exact values on the boundaries
94 c---------------------------------------------------------------------
96 c---------------------------------------------------------------------
98 c---------------------------------------------------------------------
103 do k = cell_low(3,c), cell_high(3,c)
104 zeta = dble(k) * dnzm1
106 do j = cell_low(2,c), cell_high(2,c)
107 eta = dble(j) * dnym1
108 call exact_solution(xi, eta, zeta, temp)
110 u(ii,jj,kk,m,c) = temp(m)
117 c---------------------------------------------------------------------
119 c---------------------------------------------------------------------
121 ii = cell_size(1,c)-1
124 do k = cell_low(3,c), cell_high(3,c)
125 zeta = dble(k) * dnzm1
127 do j = cell_low(2,c), cell_high(2,c)
128 eta = dble(j) * dnym1
129 call exact_solution(xi, eta, zeta, temp)
131 u(ii,jj,kk,m,c) = temp(m)
138 c---------------------------------------------------------------------
140 c---------------------------------------------------------------------
145 do k = cell_low(3,c), cell_high(3,c)
146 zeta = dble(k) * dnzm1
148 do i = cell_low(1,c), cell_high(1,c)
150 call exact_solution(xi, eta, zeta, temp)
152 u(ii,jj,kk,m,c) = temp(m)
160 c---------------------------------------------------------------------
162 c---------------------------------------------------------------------
164 jj = cell_size(2,c)-1
167 do k = cell_low(3,c), cell_high(3,c)
168 zeta = dble(k) * dnzm1
170 do i = cell_low(1,c), cell_high(1,c)
172 call exact_solution(xi, eta, zeta, temp)
174 u(ii,jj,kk,m,c) = temp(m)
181 c---------------------------------------------------------------------
183 c---------------------------------------------------------------------
188 do j = cell_low(2,c), cell_high(2,c)
189 eta = dble(j) * dnym1
191 do i =cell_low(1,c), cell_high(1,c)
193 call exact_solution(xi, eta, zeta, temp)
195 u(ii,jj,kk,m,c) = temp(m)
202 c---------------------------------------------------------------------
204 c---------------------------------------------------------------------
206 kk = cell_size(3,c)-1
209 do j = cell_low(2,c), cell_high(2,c)
210 eta = dble(j) * dnym1
212 do i =cell_low(1,c), cell_high(1,c)
214 call exact_solution(xi, eta, zeta, temp)
216 u(ii,jj,kk,m,c) = temp(m)
231 integer i, j, k, d, c, n
233 c---------------------------------------------------------------------
234 c loop over all cells
235 c---------------------------------------------------------------------
238 c---------------------------------------------------------------------
239 c first, initialize the start and end arrays
240 c---------------------------------------------------------------------
242 if (cell_coord(d,c) .eq. 1) then
247 if (cell_coord(d,c) .eq. ncells) then
254 c---------------------------------------------------------------------
255 c zap the whole left hand side for starters
256 c---------------------------------------------------------------------
258 do k = 0, cell_size(3,c)-1
259 do j = 0, cell_size(2,c)-1
260 do i = 0, cell_size(1,c)-1
261 lhs(i,j,k,n,c) = 0.0d0
267 c---------------------------------------------------------------------
268 c next, set all diagonal values to 1. This is overkill, but convenient
269 c---------------------------------------------------------------------
271 do k = 0, cell_size(3,c)-1
272 do j = 0, cell_size(2,c)-1
273 do i = 0, cell_size(1,c)-1
274 lhs(i,j,k,5*n-2,c) = 1.0d0