1//---------------------------------------------------------------------
2//---------------------------------------------------------------------
3#include "header.h"
4
5void  initialize() {
6
7//---------------------------------------------------------------------
8//---------------------------------------------------------------------
9
10//---------------------------------------------------------------------
11//     This subroutine initializes the field variable u using
12//     tri-linear transfinite interpolation of the boundary values
13//---------------------------------------------------------------------
14
15      int c, i, j, k, m, ii, jj, kk, ix, iy, iz;
16      double xi, eta, zeta, Pface[5*3*2], Pxi, Peta,
17           Pzeta, temp[5];
18#define Pface(m,n,i) Pface[(m-1)+5*((n-1)+3*(i-1))]
19#define temp(m) temp[m-1]
20
21//---------------------------------------------------------------------
22//  Later (in compute_rhs) we compute 1/u for every element. A few of
23//  the corner elements are not used, but it convenient (and faster)
24//  to compute the whole thing with a simple loop. Make sure those
25//  values are nonzero by initializing the whole thing here.
26//---------------------------------------------------------------------
27      for (c = 1; c <= ncells; c++) {
28         for (kk = -1; kk <= KMAX; kk++) {
29            for (jj = -1; jj <= JMAX; jj++) {
30               for (ii = -1; ii <= IMAX; ii++) {
31                  for (m = 1; m <= 5; m++) {
32                     u(m, ii, jj, kk, c) = 1.0;
33                  }
34               }
35            }
36         }
37      }
38//---------------------------------------------------------------------
39
40
41
42//---------------------------------------------------------------------
43//     first store the "interpolated" values everywhere on the grid
44//---------------------------------------------------------------------
45      for (c = 1; c <= ncells; c++) {
46         kk = 0;
47         for (k = cell_low(3,c); k <= cell_high(3,c); k++) {
48            zeta = (double)(k) * dnzm1;
49            jj = 0;
50            for (j = cell_low(2,c); j <= cell_high(2,c); j++) {
51               eta = (double)(j) * dnym1;
52               ii = 0;
53               for (i = cell_low(1,c); i <= cell_high(1,c); i++) {
54                  xi = (double)(i) * dnxm1;
55
56                  for (ix = 1; ix <= 2; ix++) {
57                     exact_solution((double)(ix-1), eta, zeta,
58                          &Pface(1,1,ix));
59                  }
60
61                  for (iy = 1; iy <= 2; iy++) {
62                     exact_solution(xi, (double)(iy-1) , zeta,
63                          &Pface(1,2,iy));
64                  }
65
66                  for (iz = 1; iz <= 2; iz++) {
67                     exact_solution(xi, eta, (double)(iz-1),
68                          &Pface(1,3,iz));
69                  }
70
71                  for (m = 1; m <= 5; m++) {
72                     Pxi   = xi   * Pface(m,1,2) +
73                          (1.0e0-xi)   * Pface(m,1,1);
74                     Peta  = eta  * Pface(m,2,2) +
75                          (1.0e0-eta)  * Pface(m,2,1);
76                     Pzeta = zeta * Pface(m,3,2) +
77                          (1.0e0-zeta) * Pface(m,3,1);
78
79                     u(m,ii,jj,kk,c) = Pxi + Peta + Pzeta -
80                          Pxi*Peta - Pxi*Pzeta - Peta*Pzeta +
81                          Pxi*Peta*Pzeta;
82
83                  }
84                  ii = ii + 1;
85               }
86               jj = jj + 1;
87            }
88            kk = kk+1;
89         }
90      }
91
92//---------------------------------------------------------------------
93//     now store the exact values on the boundaries
94//---------------------------------------------------------------------
95
96//---------------------------------------------------------------------
97//     west face
98//---------------------------------------------------------------------
99      c = slice(1,1);
100      ii = 0;
101      xi = 0.0e0;
102      kk = 0;
103      for (k = cell_low(3,c); k <= cell_high(3,c); k++) {
104         zeta = (double)(k) * dnzm1;
105         jj = 0;
106         for (j = cell_low(2,c); j <= cell_high(2,c); j++) {
107            eta = (double)(j) * dnym1;
108            exact_solution(xi, eta, zeta, temp);
109            for (m = 1; m <= 5; m++) {
110               u(m,ii,jj,kk,c) = temp(m);
111            }
112            jj = jj + 1;
113         }
114         kk = kk + 1;
115      }
116
117//---------------------------------------------------------------------
118//     east face
119//---------------------------------------------------------------------
120      c  = slice(1,ncells);
121      ii = cell_size(1,c)-1;
122      xi = 1.0e0;
123      kk = 0;
124      for (k = cell_low(3,c); k <= cell_high(3,c); k++) {
125         zeta = (double)(k) * dnzm1;
126         jj = 0;
127         for (j = cell_low(2,c); j <= cell_high(2,c); j++) {
128            eta = (double)(j) * dnym1;
129            exact_solution(xi, eta, zeta, temp);
130            for (m = 1; m <= 5; m++) {
131               u(m,ii,jj,kk,c) = temp(m);
132            }
133            jj = jj + 1;
134         }
135         kk = kk + 1;
136      }
137
138//---------------------------------------------------------------------
139//     south face
140//---------------------------------------------------------------------
141      c = slice(2,1);
142      jj = 0;
143      eta = 0.0e0;
144      kk = 0;
145      for (k = cell_low(3,c); k <= cell_high(3,c); k++) {
146         zeta = (double)(k) * dnzm1;
147         ii = 0;
148         for (i = cell_low(1,c); i <= cell_high(1,c); i++) {
149            xi = (double)(i) * dnxm1;
150            exact_solution(xi, eta, zeta, temp);
151            for (m = 1; m <= 5; m++) {
152               u(m,ii,jj,kk,c) = temp(m);
153            }
154            ii = ii + 1;
155         }
156         kk = kk + 1;
157      }
158
159
160//---------------------------------------------------------------------
161//     north face
162//---------------------------------------------------------------------
163      c = slice(2,ncells);
164      jj = cell_size(2,c)-1;
165      eta = 1.0e0;
166      kk = 0;
167      for (k = cell_low(3,c); k <= cell_high(3,c); k++) {
168         zeta = (double)(k) * dnzm1;
169         ii = 0;
170         for (i = cell_low(1,c); i <= cell_high(1,c); i++) {
171            xi = (double)(i) * dnxm1;
172            exact_solution(xi, eta, zeta, temp);
173            for (m = 1; m <= 5; m++) {
174               u(m,ii,jj,kk,c) = temp(m);
175            }
176            ii = ii + 1;
177         }
178         kk = kk + 1;
179      }
180
181//---------------------------------------------------------------------
182//     bottom face
183//---------------------------------------------------------------------
184      c = slice(3,1);
185      kk = 0;
186      zeta = 0.0e0;
187      jj = 0;
188      for (j = cell_low(2,c); j <= cell_high(2,c); j++) {
189         eta = (double)(j) * dnym1;
190         ii = 0;
191         for (i = cell_low(1,c); i <= cell_high(1,c); i++) {
192            xi = (double)(i) *dnxm1;
193            exact_solution(xi, eta, zeta, temp);
194            for (m = 1; m <= 5; m++) {
195               u(m,ii,jj,kk,c) = temp(m);
196            }
197            ii = ii + 1;
198         }
199         jj = jj + 1;
200      }
201
202//---------------------------------------------------------------------
203//     top face
204//---------------------------------------------------------------------
205      c = slice(3,ncells);
206      kk = cell_size(3,c)-1;
207      zeta = 1.0e0;
208      jj = 0;
209      for (j = cell_low(2,c); j <= cell_high(2,c); j++) {
210         eta = (double)(j) * dnym1;
211         ii = 0;
212         for (i = cell_low(1,c); i <= cell_high(1,c); i++) {
213            xi = (double)(i) * dnxm1;
214            exact_solution(xi, eta, zeta, temp);
215            for (m = 1; m <= 5; m++) {
216               u(m,ii,jj,kk,c) = temp(m);
217            }
218            ii = ii + 1;
219         }
220         jj = jj + 1;
221      }
222
223      return;
224}
225
226
227//---------------------------------------------------------------------
228//---------------------------------------------------------------------
229
230void lhsinit() {
231
232//---------------------------------------------------------------------
233//---------------------------------------------------------------------
234
235      int i, j, k, d, c, m, n;
236
237//---------------------------------------------------------------------
238//     loop over all cells
239//---------------------------------------------------------------------
240      for (c = 1; c <= ncells; c++) {
241
242//---------------------------------------------------------------------
243//     first, initialize the start and end arrays
244//---------------------------------------------------------------------
245         for (d = 1; d <= 3; d++) {
246            if (cell_coord(d,c) == 1) {
247               start(d,c) = 1;
248            } else {
249               start(d,c) = 0;
250            }
251            if (cell_coord(d,c) == ncells) {
252               end(d,c) = 1;
253            } else {
254               end(d,c) = 0;
255            }
256         }
257
258//---------------------------------------------------------------------
259//     zero the whole left hand side for starters
260//---------------------------------------------------------------------
261         for (k = 0; k <= cell_size(3,c)-1; k++) {
262            for (j = 0; j <= cell_size(2,c)-1; j++) {
263               for (i = 0; i <= cell_size(1,c)-1; i++) {
264                  for (m = 1; m <= 5; m++) {
265                     for (n = 1; n <= 5; n++) {
266                        lhsc(m,n,i,j,k,c) = 0.0e0;
267                     }
268                  }
269               }
270            }
271         }
272
273      }
274
275      return;
276}
277
278
279//---------------------------------------------------------------------
280//---------------------------------------------------------------------
281
282void lhsabinit(double lhsa[], double lhsb[], int size) {
283
284#define lhsa(m,n,i) lhsa[(m-1)+5*((n-1)+5*(i+1))]
285#define lhsb(m,n,i) lhsb[(m-1)+5*((n-1)+5*(i+1))]
286
287      int i, m, n;
288
289//---------------------------------------------------------------------
290//     next, set all diagonal values to 1. This is overkill, but convenient
291//---------------------------------------------------------------------
292      for (i = 0; i <= size; i++) {
293         for (m = 1; m <= 5; m++) {
294            for (n = 1; n <= 5; n++) {
295               lhsa(m,n,i) = 0.0e0;
296               lhsb(m,n,i) = 0.0e0;
297            }
298            lhsb(m,m,i) = 1.0e0;
299         }
300      }
301
302      return;
303}
304
305
306
307