Commit | Line | Data |
---|---|---|
e3f82424 MW |
1 | // -------------------------------------------------------------------------- |
2 | // | |
3 | // Copyright | |
4 | // Markus Wittmann, 2016-2017 | |
5 | // RRZE, University of Erlangen-Nuremberg, Germany | |
6 | // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de | |
7 | // | |
8 | // Viktor Haag, 2016 | |
9 | // LSS, University of Erlangen-Nuremberg, Germany | |
10 | // | |
11 | // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels). | |
12 | // | |
13 | // LbmBenchKernels is free software: you can redistribute it and/or modify | |
14 | // it under the terms of the GNU General Public License as published by | |
15 | // the Free Software Foundation, either version 3 of the License, or | |
16 | // (at your option) any later version. | |
17 | // | |
18 | // LbmBenchKernels is distributed in the hope that it will be useful, | |
19 | // but WITHOUT ANY WARRANTY; without even the implied warranty of | |
20 | // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the | |
21 | // GNU General Public License for more details. | |
22 | // | |
23 | // You should have received a copy of the GNU General Public License | |
24 | // along with LbmBenchKernels. If not, see <http://www.gnu.org/licenses/>. | |
25 | // | |
26 | // -------------------------------------------------------------------------- | |
27 | #include "BenchKernelD3Q19AaVecCommon.h" | |
28 | ||
29 | #include "Memory.h" | |
30 | #include "Vtk.h" | |
31 | #include "Vector.h" | |
32 | ||
33 | #include <inttypes.h> | |
34 | #include <math.h> | |
35 | ||
36 | #ifdef _OPENMP | |
37 | #include <omp.h> | |
38 | #endif | |
39 | ||
40 | // Forward definition. | |
41 | void FNAME(D3Q19AaVecKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd); | |
42 | ||
43 | ||
44 | static void FNAME(BcGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf) | |
45 | { | |
46 | Assert(kd != NULL); | |
47 | Assert(kd->PdfsActive != NULL); | |
48 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
49 | Assert(pdf != NULL); | |
50 | ||
51 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
52 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
53 | Assert(dir >= 0); Assert(dir < N_D3Q19); | |
54 | ||
55 | KernelDataAa * kda = KDA(kd); | |
56 | ||
57 | int oX = kd->Offsets[0]; | |
58 | int oY = kd->Offsets[1]; | |
59 | int oZ = kd->Offsets[2]; | |
60 | ||
61 | if (kda->Iteration % 2 == 0) { | |
62 | // Pdfs are stored inverse, local PDFs are located in remote nodes | |
63 | int nx = x - D3Q19_X[dir]; | |
64 | int ny = y - D3Q19_Y[dir]; | |
65 | int nz = z - D3Q19_Z[dir]; | |
66 | ||
67 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
68 | *pdf = kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, D3Q19_INV[dir])]; | |
69 | #undef I | |
70 | } | |
71 | else { | |
72 | int nx = x; | |
73 | int ny = y; | |
74 | int nz = z; | |
75 | ||
76 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
77 | *pdf = kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, dir)]; | |
78 | #undef I | |
79 | } | |
80 | ||
81 | ||
82 | return; | |
83 | } | |
84 | ||
85 | static void FNAME(BcSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf) | |
86 | { | |
87 | Assert(kd != NULL); | |
88 | Assert(kd->PdfsActive != NULL); | |
89 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
90 | ||
91 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
92 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
93 | Assert(dir >= 0); Assert(dir < N_D3Q19); | |
94 | ||
95 | KernelDataAa * kda = KDA(kd); | |
96 | ||
97 | int oX = kd->Offsets[0]; | |
98 | int oY = kd->Offsets[1]; | |
99 | int oZ = kd->Offsets[2]; | |
100 | ||
101 | if (kda->Iteration % 2 == 0) { | |
102 | // Pdfs are stored inverse, local PDFs are located in remote nodes | |
103 | int nx = x - D3Q19_X[dir]; | |
104 | int ny = y - D3Q19_Y[dir]; | |
105 | int nz = z - D3Q19_Z[dir]; | |
106 | ||
107 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
108 | pdf = kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, D3Q19_INV[dir])] = pdf; | |
109 | #undef I | |
110 | } | |
111 | else { | |
112 | int nx = x; | |
113 | int ny = y; | |
114 | int nz = z; | |
115 | ||
116 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
117 | kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, dir)] = pdf; | |
118 | #undef I | |
119 | } | |
120 | ||
121 | return; | |
122 | } | |
123 | ||
124 | ||
125 | static void FNAME(GetNode)(KernelData * kd, int x, int y, int z, PdfT * pdfs) | |
126 | { | |
127 | Assert(kd != NULL); | |
128 | Assert(kd->PdfsActive != NULL); | |
129 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
130 | Assert(pdfs != NULL); | |
131 | ||
132 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
133 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
134 | ||
135 | KernelDataAa * kda = KDA(kd); | |
136 | ||
137 | int oX = kd->Offsets[0]; | |
138 | int oY = kd->Offsets[1]; | |
139 | int oZ = kd->Offsets[2]; | |
140 | ||
141 | ||
142 | if (kda->Iteration % 2 == 0) { | |
143 | // Pdfs are stored inverse, local PDFs are located in remote nodes | |
144 | ||
145 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
146 | #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x + oX - _x, y + oY - _y, z + oZ - _z, D3Q19_INV[idx])]; | |
147 | D3Q19_LIST | |
148 | #undef X | |
149 | #undef I | |
150 | } | |
151 | else { | |
152 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
153 | #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x + oX, y + oY, z + oZ, idx)]; | |
154 | D3Q19_LIST | |
155 | #undef X | |
156 | #undef I | |
157 | ||
158 | } | |
159 | ||
160 | #if 0 // DETECT NANs | |
161 | ||
162 | for (int d = 0; d < 19; ++d) { | |
163 | if (isnan(pdfs[d])) { | |
164 | printf("%d %d %d %d nan! get node\n", x, y, z, d); | |
165 | ||
166 | for (int d2 = 0; d2 < 19; ++d2) { | |
167 | printf("%d: %e\n", d2, pdfs[d2]); | |
168 | } | |
169 | ||
170 | exit(1); | |
171 | } | |
172 | } | |
173 | ||
174 | #endif | |
175 | ||
176 | return; | |
177 | } | |
178 | ||
179 | ||
180 | static void FNAME(SetNode)(KernelData * kd, int x, int y, int z, PdfT * pdfs) | |
181 | { | |
182 | Assert(kd != NULL); | |
183 | Assert(kd->PdfsActive != NULL); | |
184 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
185 | Assert(pdfs != NULL); | |
186 | ||
187 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
188 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
189 | ||
190 | KernelDataAa * kda = KDA(kd); | |
191 | ||
192 | int oX = kd->Offsets[0]; | |
193 | int oY = kd->Offsets[1]; | |
194 | int oZ = kd->Offsets[2]; | |
195 | ||
196 | if (kda->Iteration % 2 == 0) { | |
197 | // Pdfs are stored inverse, local PDFs are located in remote nodes | |
198 | ||
199 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
200 | #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x + oX - _x, y + oY - _y, z + oZ - _z, D3Q19_INV[idx])] = pdfs[idx]; | |
201 | D3Q19_LIST | |
202 | #undef X | |
203 | #undef I | |
204 | } | |
205 | else { | |
206 | #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) | |
207 | #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x + oX, y + oY, z + oZ, idx)] = pdfs[idx]; | |
208 | D3Q19_LIST | |
209 | #undef X | |
210 | #undef I | |
211 | } | |
212 | return; | |
213 | } | |
214 | ||
215 | ||
216 | static void ParameterUsage() | |
217 | { | |
218 | printf("Kernel parameters:\n"); | |
219 | printf(" [-blk <n>] [-blk-[xyz] <n>]\n"); | |
220 | ||
221 | return; | |
222 | } | |
223 | ||
224 | static void ParseParameters(Parameters * params, int * blk) | |
225 | { | |
226 | Assert(blk != NULL); | |
227 | ||
228 | blk[0] = 0; blk[1] = 0; blk[2] = 0; | |
229 | ||
230 | #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param)) | |
231 | #define NEXT_ARG_PRESENT() \ | |
232 | do { \ | |
233 | if (i + 1 >= params->nKernelArgs) { \ | |
234 | printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \ | |
235 | exit(1); \ | |
236 | } \ | |
237 | } while (0) | |
238 | ||
239 | ||
240 | for (int i = 0; i < params->nKernelArgs; ++i) { | |
241 | if (ARG_IS("-blk") || ARG_IS("--blk")) { | |
242 | NEXT_ARG_PRESENT(); | |
243 | ||
244 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
245 | ||
246 | if (tmp < 0) { | |
247 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
248 | exit(1); | |
249 | } | |
250 | ||
251 | blk[0] = blk[1] = blk[2] = tmp; | |
252 | } | |
253 | else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) { | |
254 | NEXT_ARG_PRESENT(); | |
255 | ||
256 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
257 | ||
258 | if (tmp < 0) { | |
259 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
260 | exit(1); | |
261 | } | |
262 | ||
263 | blk[0] = tmp; | |
264 | } | |
265 | else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) { | |
266 | NEXT_ARG_PRESENT(); | |
267 | ||
268 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
269 | ||
270 | if (tmp < 0) { | |
271 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
272 | exit(1); | |
273 | } | |
274 | ||
275 | blk[1] = tmp; | |
276 | } | |
277 | else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) { | |
278 | NEXT_ARG_PRESENT(); | |
279 | ||
280 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
281 | ||
282 | if (tmp < 0) { | |
283 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
284 | exit(1); | |
285 | } | |
286 | ||
287 | blk[2] = tmp; | |
288 | } | |
289 | else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) { | |
290 | ParameterUsage(); | |
291 | exit(1); | |
292 | } | |
293 | else { | |
294 | printf("ERROR: unknown kernel parameter.\n"); | |
295 | ParameterUsage(); | |
296 | exit(1); | |
297 | } | |
298 | } | |
299 | ||
300 | #undef ARG_IS | |
301 | #undef NEXT_ARG_PRESENT | |
302 | ||
303 | return; | |
304 | } | |
305 | ||
306 | ||
307 | void FNAME(D3Q19AaVecInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params) | |
308 | { | |
309 | KernelDataAa * kda = NULL; | |
310 | MemAlloc((void **)&kda, sizeof(KernelDataAa)); | |
311 | ||
312 | kda->Blk[0] = 0; kda->Blk[1] = 0; kda->Blk[2] = 0; | |
313 | kda->Iteration = -1; | |
314 | ||
315 | KernelData * kd = &kda->kd; | |
316 | *kernelData = kd; | |
317 | ||
318 | kd->nObstIndices = ld->nObst; | |
319 | ||
320 | // Ajust the dimensions according to padding, if used. | |
321 | kd->Dims[0] = ld->Dims[0]; | |
322 | kd->Dims[1] = ld->Dims[1]; | |
323 | kd->Dims[2] = ld->Dims[2]; | |
324 | ||
325 | ||
326 | int * lDims = ld->Dims; | |
327 | int * gDims = kd->GlobalDims; | |
328 | ||
329 | Assert(VSIZE <= 4); | |
330 | ||
331 | // TODO: only add enough ghost cells so we can compute everything vectorized. | |
332 | gDims[0] = lDims[0] + 2; | |
333 | gDims[1] = lDims[1] + 2; | |
334 | // TODO: fix this for aa-vec2-soa | |
0fde6e45 | 335 | gDims[2] = lDims[2] + 2 + VSIZE - 2; // one ghost cell in front, one in the back, plus at most two at the back for VSIZE = 4 |
e3f82424 MW |
336 | |
337 | kd->Offsets[0] = 1; | |
338 | kd->Offsets[1] = 1; | |
339 | kd->Offsets[2] = 1; | |
340 | ||
341 | int lX = lDims[0]; | |
342 | int lY = lDims[1]; | |
343 | int lZ = lDims[2]; | |
344 | ||
345 | int gX = gDims[0]; | |
346 | int gY = gDims[1]; | |
347 | int gZ = gDims[2]; | |
348 | ||
349 | int oX = kd->Offsets[0]; | |
350 | int oY = kd->Offsets[1]; | |
351 | int oZ = kd->Offsets[2]; | |
352 | ||
353 | int blk[3] = { 0 }; | |
354 | ||
355 | int nCells = gX * gY * gZ; | |
356 | ||
357 | PdfT * pdfs[2] = { NULL, NULL }; | |
358 | ||
359 | ParseParameters(params, blk); | |
360 | ||
361 | if (blk[2] % VSIZE != 0) { | |
362 | blk[2] -= blk[2] % VSIZE; | |
363 | printf("WARNING: blocking factor for z direction must be a multiple of VSIZE = %d, adjusting it to %d.\n", VSIZE, blk[2]); | |
364 | } | |
365 | ||
366 | if (blk[0] == 0) blk[0] = gX; | |
367 | if (blk[1] == 0) blk[1] = gY; | |
368 | if (blk[2] == 0) blk[2] = gZ; | |
369 | ||
370 | printf("# blocking x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]); | |
371 | ||
372 | kda->Blk[0] = blk[0]; kda->Blk[1] = blk[1]; kda->Blk[2] = blk[2]; | |
373 | ||
374 | ||
375 | printf("# allocating data for %d LB nodes with padding (%lu bytes = %f MiB for the single lattice)\n", | |
376 | nCells, | |
377 | sizeof(PdfT) * nCells * N_D3Q19, | |
378 | sizeof(PdfT) * nCells * N_D3Q19 / 1024.0 / 1024.0); | |
379 | ||
380 | #define PAGE_4K 4096 | |
381 | ||
382 | MemAllocAligned((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19, PAGE_4K); | |
383 | ||
384 | kd->Pdfs[0] = pdfs[0]; | |
385 | kd->Pdfs[1] = NULL; | |
386 | ||
387 | ||
388 | // Initialize PDFs with some (arbitrary) data for correct NUMA placement. | |
389 | // This depends on the chosen data layout. | |
390 | // The structure of the loop should resemble the same "execution layout" | |
391 | // as in the kernel! | |
392 | ||
393 | int nThreads; | |
394 | #ifdef _OPENMP | |
395 | nThreads = omp_get_max_threads(); | |
396 | #endif | |
397 | ||
398 | #ifdef _OPENMP | |
399 | #pragma omp parallel for \ | |
400 | shared(gDims, pdfs, \ | |
401 | oX, oY, oZ, lX, lY, lZ, blk, nThreads, ld) | |
402 | #endif | |
403 | for (int i = 0; i < nThreads; ++i) { | |
404 | ||
405 | int threadStartX = lX / nThreads * i; | |
406 | int threadEndX = lX / nThreads * (i + 1); | |
407 | ||
408 | if (lX % nThreads > 0) { | |
409 | if (lX % nThreads > i) { | |
410 | threadStartX += i; | |
411 | threadEndX += i + 1; | |
412 | } | |
413 | else { | |
414 | threadStartX += lX % nThreads; | |
415 | threadEndX += lX % nThreads; | |
416 | } | |
417 | } | |
418 | ||
419 | for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) { | |
420 | for (int bY = oY; bY < lY + oY; bY += blk[1]) { | |
421 | for (int bZ = oZ; bZ < lZ + oZ; bZ += blk[2]) { | |
422 | ||
423 | int eX = MIN(bX + blk[0], threadEndX + oX); | |
424 | int eY = MIN(bY + blk[1], lY + oY); | |
425 | int eZ = MIN(bZ + blk[2], lZ + oZ); | |
426 | ||
427 | // printf("%d: %d-%d %d-%d %d-%d %d - %d\n", omp_get_thread_num(), bZ, eZ, bY, eY, bX, eX, threadStartX, threadEndX); | |
428 | ||
429 | for (int x = bX; x < eX; ++x) { | |
430 | for (int y = bY; y < eY; ++y) { | |
431 | for (int z = bZ; z < eZ; ++z) { | |
432 | ||
433 | for (int d = 0; d < N_D3Q19; ++d) { | |
434 | pdfs[0][P_INDEX_5(gDims, x, y, z, d)] = -50.0; | |
435 | } | |
436 | ||
437 | } } } // x, y, z | |
438 | } } } // bx, by, bz | |
439 | } | |
440 | ||
441 | ||
442 | // Initialize all PDFs to some standard value. | |
443 | for (int x = oX; x < lX + oX; ++x) { | |
444 | for (int y = oY; y < lY + oY; ++y) { | |
445 | for (int z = oZ; z < lZ + oZ; ++z) { | |
446 | for (int d = 0; d < N_D3Q19; ++d) { | |
447 | pdfs[0][P_INDEX_5(gDims, x, y, z, d)] = 0.0; | |
448 | } | |
449 | } } } // x, y, z | |
450 | ||
451 | ||
452 | // Count how many *PDFs* need bounce back treatment. | |
453 | ||
454 | uint64_t nPdfs = ((uint64_t)19) * gX * gY * gZ; | |
455 | ||
456 | if (nPdfs > ((2LU << 31) - 1)) { | |
457 | printf("ERROR: number of PDFs exceed 2^31.\n"); | |
458 | exit(1); | |
459 | } | |
460 | ||
461 | // Compiler bug? Incorrect computation of nBounceBackPdfs when using icc 15.0.2. | |
462 | // Works when declaring nBounceBackPdfs as int64_t or using volatile. | |
463 | volatile int nBounceBackPdfs = 0; | |
464 | // int64_t nBounceBackPdfs = 0; | |
465 | int nx, ny, nz, px, py, pz; | |
466 | ||
467 | ||
468 | for (int x = 0; x < lX; ++x) { | |
469 | for (int y = 0; y < lY; ++y) { | |
470 | for (int z = 0; z < lZ; ++z) { | |
471 | ||
472 | if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) { | |
473 | for (int d = 0; d < N_D3Q19; ++d) { | |
474 | nx = x - D3Q19_X[d]; | |
475 | ny = y - D3Q19_Y[d]; | |
476 | nz = z - D3Q19_Z[d]; | |
477 | ||
478 | // Check if neighbor is inside the lattice. | |
479 | // if(nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { | |
480 | // continue; | |
481 | // } | |
482 | if ((nx < 0 || nx >= lX) && ld->PeriodicX) { | |
483 | ++nBounceBackPdfs; // Compiler bug --> see above | |
484 | } | |
485 | else if ((ny < 0 || ny >= lY) && ld->PeriodicY) { | |
486 | ++nBounceBackPdfs; // Compiler bug --> see above | |
487 | } | |
488 | else if ((nz < 0 || nz >= lZ) && ld->PeriodicZ) { | |
489 | ++nBounceBackPdfs; // Compiler bug --> see above | |
490 | } | |
491 | else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { | |
492 | continue; | |
493 | } | |
494 | else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) { | |
495 | ++nBounceBackPdfs; // Compiler bug --> see above | |
496 | } | |
497 | } | |
498 | } | |
499 | } | |
500 | } | |
501 | } | |
502 | ||
503 | printf("# allocating %d indices for bounce back pdfs (%s for source and destination array)\n", nBounceBackPdfs, ByteToHuman(sizeof(int) * nBounceBackPdfs * 2)); | |
504 | ||
505 | MemAlloc((void **) & (kd->BounceBackPdfsSrc), sizeof(int) * nBounceBackPdfs + 100); | |
506 | MemAlloc((void **) & (kd->BounceBackPdfsDst), sizeof(int) * nBounceBackPdfs + 100); | |
507 | ||
508 | kd->nBounceBackPdfs = nBounceBackPdfs; | |
509 | nBounceBackPdfs = 0; | |
510 | ||
511 | int srcIndex; | |
512 | int dstIndex; | |
513 | ||
514 | // TODO: currently this is not NUMA-aware | |
515 | // - maybe use the same blocking as for lattice initialization? | |
516 | // - do place the bounce back index vector parallel? | |
517 | ||
518 | for (int x = 0; x < lX; ++x) { | |
519 | for (int y = 0; y < lY; ++y) { | |
520 | for (int z = 0; z < lZ; ++z) { | |
521 | ||
522 | if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) { | |
523 | for (int d = 0; d < N_D3Q19; ++d) { | |
524 | nx = x + D3Q19_X[d]; | |
525 | ny = y + D3Q19_Y[d]; | |
526 | nz = z + D3Q19_Z[d]; | |
527 | ||
528 | if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) || | |
529 | ((ny < 0 || ny >= lY) && ld->PeriodicY) || | |
530 | ((nz < 0 || nz >= lZ) && ld->PeriodicZ) | |
531 | ){ | |
532 | // For periodicity: | |
533 | ||
534 | // We assume we have finished odd time step (accessing neighbor PDFs) and are | |
535 | // before executing the even time step (accessing local PDFs only). | |
536 | ||
537 | // Assuming we are at the most east position of the lattice. Through the odd | |
538 | // time step propagation has put a PDF in the east slot of the ghost cell east | |
539 | // of us, i.e. nx, ny, nz. We copy it to the east slot of the most west node. | |
540 | ||
541 | // In case of transition from even to odd time step , src and dest must be | |
542 | // exchanged. | |
543 | ||
544 | ||
545 | // x periodic | |
546 | if (nx < 0) { | |
547 | px = lX - 1; | |
548 | } | |
549 | else if (nx >= lX) { | |
550 | px = 0; | |
551 | } else { | |
552 | px = nx; | |
553 | } | |
554 | ||
555 | // y periodic | |
556 | if (ny < 0) { | |
557 | py = lY - 1; | |
558 | } | |
559 | else if (ny >= lY) { | |
560 | py = 0; | |
561 | } else { | |
562 | py = ny; | |
563 | } | |
564 | ||
565 | // z periodic | |
566 | if (nz < 0) { | |
567 | pz = lZ - 1; | |
568 | } | |
569 | else if (nz >= lZ) { | |
570 | pz = 0; | |
571 | } else { | |
572 | pz = nz; | |
573 | } | |
574 | ||
575 | if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) { | |
576 | // See description of bounce back handling below. | |
577 | srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); | |
578 | dstIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, D3Q19_INV[d]); | |
579 | } | |
580 | else { | |
581 | ||
582 | srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); | |
583 | // Put it on the other side back into the domain. | |
584 | dstIndex = P_INDEX_5(gDims, px + oX, py + oY, pz + oZ, d); | |
585 | ||
586 | VerifyMsg(nBounceBackPdfs < kd->nBounceBackPdfs, "nBBPdfs %d < kd->nBBPdfs %d xyz: %d %d %d d: %d\n", nBounceBackPdfs, kd->nBounceBackPdfs, x, y, z, d); | |
587 | ||
588 | } | |
589 | ||
590 | kd->BounceBackPdfsSrc[nBounceBackPdfs] = srcIndex; | |
591 | kd->BounceBackPdfsDst[nBounceBackPdfs] = dstIndex; | |
592 | ||
593 | ++nBounceBackPdfs; | |
594 | ||
595 | } | |
596 | else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { | |
597 | continue; | |
598 | } | |
599 | else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) { | |
600 | // Depending on the time step we are in we have to exchange src and dst index. | |
601 | ||
602 | // We build the list for the case, when we have finished odd time step | |
603 | // (accessing neighbor PDFs) and before we start with the even time step | |
604 | // (accessing local PDFs only). | |
605 | ||
606 | // Assume our neighbor east of us, i.e. nx, ny, nz, is an obstacle cell. | |
607 | // Then we have to move the east PDF from the obstacle to our west position, | |
608 | // i.e. the inverse of east. | |
609 | ||
610 | // In case of transition from even to odd time step src and dest just | |
611 | // have to be exchanged. | |
612 | ||
613 | srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); | |
614 | dstIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, D3Q19_INV[d]); | |
615 | ||
616 | VerifyMsg(nBounceBackPdfs < kd->nBounceBackPdfs, "nBBPdfs %d < kd->nBBPdfs %d xyz: %d %d %d d: %d\n", nBounceBackPdfs, kd->nBounceBackPdfs, x, y, z, d); | |
617 | ||
618 | kd->BounceBackPdfsSrc[nBounceBackPdfs] = srcIndex; | |
619 | kd->BounceBackPdfsDst[nBounceBackPdfs] = dstIndex; | |
620 | ||
621 | ++nBounceBackPdfs; | |
622 | } | |
623 | } | |
624 | } | |
625 | } | |
626 | } | |
627 | } | |
628 | ||
629 | ||
630 | // Fill remaining KernelData structures | |
631 | kd->GetNode = FNAME(GetNode); | |
632 | kd->SetNode = FNAME(SetNode); | |
633 | ||
634 | kd->BoundaryConditionsGetPdf = FNAME(BcGetPdf); | |
635 | kd->BoundaryConditionsSetPdf = FNAME(BcSetPdf); | |
636 | ||
637 | kd->Kernel = FNAME(D3Q19AaVecKernel); | |
638 | ||
639 | kd->DstPdfs = NULL; | |
640 | kd->PdfsActive = kd->Pdfs[0]; | |
641 | ||
642 | return; | |
643 | } | |
644 | ||
645 | void FNAME(D3Q19AaVecDeinit)(LatticeDesc * ld, KernelData ** kernelData) | |
646 | { | |
647 | MemFree((void **) & ((*kernelData)->Pdfs[0])); | |
648 | ||
649 | MemFree((void **) & ((*kernelData)->BounceBackPdfsSrc)); | |
650 | MemFree((void **) & ((*kernelData)->BounceBackPdfsDst)); | |
651 | ||
652 | MemFree((void **)kernelData); | |
653 | ||
654 | return; | |
655 | } | |
656 |