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