Commit | Line | Data |
---|---|---|
10988083 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 "BenchKernelD3Q19ListCommon.h" | |
28 | ||
29 | #include "Memory.h" | |
30 | #include "Vtk.h" | |
e3f82424 | 31 | #include "Padding.h" |
10988083 MW |
32 | |
33 | #include <math.h> | |
34 | ||
35 | ||
36 | // Forward definition. | |
37 | void FNAME(D3Q19ListKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd); | |
38 | ||
39 | ||
40 | ||
41 | ||
42 | // ----------------------------------------------------------------------- | |
43 | // Functions which are used as callback by the kernel to read or write | |
44 | // PDFs and nodes. | |
45 | ||
46 | static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf) | |
47 | { | |
48 | Assert(kd != NULL); | |
49 | Assert(kd->PdfsActive != NULL); | |
50 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
51 | Assert(pdf != NULL); | |
52 | ||
53 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
54 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
55 | Assert(dir >= 0); Assert(dir < N_D3Q19); | |
56 | ||
57 | #if 0 | |
58 | *pdf = kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)]; | |
59 | #else | |
60 | #ifdef PROP_MODEL_PUSH | |
61 | *pdf = kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)]; | |
62 | #elif PROP_MODEL_PULL | |
63 | ||
64 | ||
65 | // The relevant PDFs here are the ones, which will get streamed in later | |
66 | // during propagation. So we must return the *remote* PDFs. | |
67 | uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)]; | |
68 | ||
69 | if (dir != D3Q19_C) { | |
70 | ||
71 | uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX; | |
72 | ||
73 | *pdf = kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + dir]]; | |
74 | } | |
75 | else { | |
76 | *pdf = kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)]; | |
77 | ||
78 | } | |
79 | #endif | |
80 | #endif | |
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 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
91 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
92 | Assert(dir >= 0); Assert(dir < N_D3Q19); | |
93 | ||
94 | #if 0 | |
95 | if (isnan(pdf)) { | |
96 | printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]); | |
97 | DEBUG_BREAK_POINT(); | |
98 | exit(1); | |
99 | } | |
100 | #endif | |
101 | ||
102 | #if 0 | |
103 | kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)] = pdf; | |
104 | #else | |
105 | #ifdef PROP_MODEL_PUSH | |
106 | kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)] = pdf; | |
107 | #elif PROP_MODEL_PULL | |
108 | ||
109 | // The relevant PDFs here are the ones, which will get streamed in later | |
110 | // during propagation. So we must set this *remote* PDFs. | |
111 | uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)]; | |
112 | ||
113 | if (dir != D3Q19_C) { | |
114 | ||
115 | uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX; | |
116 | ||
117 | kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + dir]] = pdf; | |
118 | } | |
119 | else { | |
120 | kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)] = pdf; | |
121 | ||
122 | } | |
123 | #endif | |
124 | #endif | |
125 | ||
126 | return; | |
127 | } | |
128 | ||
129 | ||
130 | static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs) | |
131 | { | |
132 | Assert(kd != NULL); | |
133 | Assert(kd->PdfsActive != NULL); | |
134 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
135 | Assert(pdfs != NULL); | |
136 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
137 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
138 | ||
139 | PdfT sum = 0.0; | |
140 | ||
141 | // TODO: pull scheme? | |
142 | #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir)) | |
143 | #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)]; sum += pdfs[idx]; | |
144 | D3Q19_LIST | |
145 | #undef X | |
146 | #undef I | |
147 | ||
148 | // if (sum < 0.0) { | |
149 | // printf("%d %d %d negative density \n", x, y, z); | |
150 | // exit(1); | |
151 | // } | |
152 | ||
153 | #if 0 | |
154 | for (int d = 0; d < 19; ++d) { | |
155 | if(isnan(pdfs[d]) || isinf(pdfs[d])) { | |
156 | printf("%d %d %d %d nan! get node\n", x, y, z, d); | |
157 | for (int d2 = 0; d2 < 19; ++d2) { | |
158 | printf("%d: %e\n", d2, pdfs[d2]); | |
159 | } | |
160 | exit(1); | |
161 | } | |
162 | } | |
163 | #endif | |
164 | return; | |
165 | } | |
166 | ||
167 | ||
168 | static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs) | |
169 | { | |
170 | Assert(kd != NULL); | |
171 | Assert(kd->PdfsActive != NULL); | |
172 | Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); | |
173 | Assert(pdfs != NULL); | |
174 | ||
175 | Assert(x >= 0); Assert(y >= 0); Assert(z >= 0); | |
176 | Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]); | |
177 | ||
178 | #if 0 | |
179 | for (int d = 0; d < 19; ++d) { | |
180 | if(isnan(pdfs[d])) { | |
181 | printf("%d %d %d %d nan! get node\n", x, y, z, d); | |
182 | for (int d2 = 0; d2 < 19; ++d2) { | |
183 | printf("%d: %e\n", d2, pdfs[d2]); | |
184 | } | |
185 | exit(1); | |
186 | } | |
187 | } | |
188 | #endif | |
189 | ||
190 | // TODO: pull scheme? | |
191 | #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir)) | |
192 | #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx]; | |
193 | D3Q19_LIST | |
194 | #undef X | |
195 | #undef I | |
196 | ||
197 | return; | |
198 | } | |
199 | ||
200 | static void ParameterUsage() | |
201 | { | |
202 | printf("Kernel parameters:\n"); | |
203 | printf(" [-blk <n>] [-blk-[xyz] <n>]\n"); | |
e3f82424 MW |
204 | #ifdef DATA_LAYOUT_SOA |
205 | printf(" [-pad auto|modulus_1+offset_1(,modulus_n+offset_n)*]\n"); | |
206 | #endif | |
10988083 MW |
207 | |
208 | return; | |
209 | } | |
210 | ||
e3f82424 | 211 | static void ParseParameters(Parameters * params, int * blk, PadInfo ** padInfo) |
10988083 MW |
212 | { |
213 | Assert(blk != NULL); | |
214 | ||
215 | blk[0] = 0; blk[1] = 0; blk[2] = 0; | |
e3f82424 | 216 | *padInfo = NULL; |
10988083 MW |
217 | |
218 | #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param)) | |
219 | #define NEXT_ARG_PRESENT() \ | |
220 | do { \ | |
221 | if (i + 1 >= params->nKernelArgs) { \ | |
222 | printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \ | |
223 | exit(1); \ | |
224 | } \ | |
225 | } while (0) | |
226 | ||
227 | ||
228 | for (int i = 0; i < params->nKernelArgs; ++i) { | |
229 | if (ARG_IS("-blk") || ARG_IS("--blk")) { | |
230 | NEXT_ARG_PRESENT(); | |
231 | ||
232 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
233 | ||
e3f82424 MW |
234 | if (tmp < 0) { |
235 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
10988083 MW |
236 | exit(1); |
237 | } | |
238 | ||
239 | blk[0] = blk[1] = blk[2] = tmp; | |
240 | } | |
241 | else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) { | |
242 | NEXT_ARG_PRESENT(); | |
243 | ||
244 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
245 | ||
e3f82424 MW |
246 | if (tmp < 0) { |
247 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
10988083 MW |
248 | exit(1); |
249 | } | |
250 | ||
251 | blk[0] = tmp; | |
252 | } | |
253 | else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) { | |
254 | NEXT_ARG_PRESENT(); | |
255 | ||
256 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
257 | ||
e3f82424 MW |
258 | if (tmp < 0) { |
259 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
10988083 MW |
260 | exit(1); |
261 | } | |
262 | ||
263 | blk[1] = tmp; | |
264 | } | |
265 | else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) { | |
266 | NEXT_ARG_PRESENT(); | |
267 | ||
268 | int tmp = strtol(params->KernelArgs[++i], NULL, 0); | |
269 | ||
e3f82424 MW |
270 | if (tmp < 0) { |
271 | printf("ERROR: blocking parameter must be >= 0.\n"); | |
10988083 MW |
272 | exit(1); |
273 | } | |
274 | ||
275 | blk[2] = tmp; | |
276 | } | |
277 | else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) { | |
278 | ParameterUsage(); | |
279 | exit(1); | |
280 | } | |
e3f82424 MW |
281 | #ifdef DATA_LAYOUT_SOA |
282 | else if (ARG_IS("-pad") || ARG_IS("--pad")) { | |
283 | NEXT_ARG_PRESENT(); | |
284 | ||
285 | *padInfo = PadInfoFromStr(params->KernelArgs[++i]); | |
286 | } | |
287 | #endif | |
10988083 MW |
288 | else { |
289 | printf("ERROR: unknown kernel parameter.\n"); | |
290 | ParameterUsage(); | |
291 | exit(1); | |
292 | } | |
293 | } | |
294 | ||
295 | #undef ARG_IS | |
296 | #undef NEXT_ARG_PRESENT | |
297 | ||
298 | return; | |
299 | } | |
300 | ||
301 | void FNAME(D3Q19ListInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params) | |
302 | { | |
303 | KernelData * kd; | |
304 | KernelDataList * kdl; | |
305 | MemAlloc((void **)&kdl, sizeof(KernelDataList)); | |
306 | ||
307 | kd = (KernelData *)kdl; | |
308 | *kernelData = kd; | |
309 | ||
310 | #ifdef DEBUG | |
311 | kd->Pdfs[0] = NULL; | |
312 | kd->Pdfs[1] = NULL; | |
313 | kd->PdfsActive = NULL; | |
314 | kd->DstPdfs = NULL; | |
315 | kd->SrcPdfs = NULL; | |
316 | kd->Dims[0] = -1; | |
317 | kd->Dims[1] = -1; | |
318 | kd->Dims[2] = -1; | |
319 | kd->GlobalDims[0] = -1; | |
320 | kd->GlobalDims[1] = -1; | |
321 | kd->GlobalDims[2] = -1; | |
322 | kd->Offsets[0] = -1; | |
323 | kd->Offsets[1] = -1; | |
324 | kd->Offsets[2] = -1; | |
325 | ||
326 | kd->ObstIndices = NULL; | |
327 | kd->nObstIndices = -1; | |
328 | kd->BounceBackPdfsSrc = NULL; | |
329 | kd->BounceBackPdfsDst = NULL; | |
330 | kd->nBounceBackPdfs = -1; | |
331 | ||
332 | kdl->AdjList = NULL; | |
333 | kdl->Coords = NULL; | |
334 | kdl->Grid = NULL; | |
335 | kdl->nCells = -1; | |
336 | kdl->nFluid = -1; | |
337 | #endif | |
338 | ||
e3f82424 MW |
339 | int blk[3] = { 0 }; |
340 | PadInfo * padInfo = NULL; | |
341 | ||
342 | ParseParameters(params, blk, &padInfo); | |
343 | ||
344 | ||
10988083 MW |
345 | // Ajust the dimensions according to padding, if used. |
346 | kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0]; | |
347 | kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1]; | |
348 | kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2]; | |
349 | ||
350 | int * lDims = ld->Dims; | |
351 | ||
352 | int lX = lDims[0]; | |
353 | int lY = lDims[1]; | |
354 | int lZ = lDims[2]; | |
355 | ||
356 | int nTotalCells = lX * lY * lZ; | |
e3f82424 | 357 | int nCells = ld->nFluid; |
10988083 MW |
358 | int nFluid = ld->nFluid; |
359 | ||
e3f82424 MW |
360 | #ifdef DATA_LAYOUT_SOA |
361 | { | |
362 | nCells = PadCellsAndReport(nCells, sizeof(PdfT), &padInfo); | |
363 | PadInfoFree(padInfo); padInfo = NULL; | |
364 | } | |
365 | #endif | |
366 | ||
10988083 MW |
367 | kdl->nCells = nCells; |
368 | kdl->nFluid = nFluid; | |
369 | ||
370 | PdfT * pdfs[2]; | |
371 | ||
10988083 MW |
372 | if (blk[0] == 0) blk[0] = lX; |
373 | if (blk[1] == 0) blk[1] = lY; | |
374 | if (blk[2] == 0) blk[2] = lZ; | |
375 | ||
376 | printf("# blocking x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]); | |
377 | ||
378 | printf("# allocating data for %d fluid LB nodes with padding (%lu bytes = %f MiB for both lattices)\n", | |
379 | nCells, 2 * sizeof(PdfT) * nCells * N_D3Q19, | |
380 | 2 * sizeof(PdfT) * nCells * N_D3Q19 / 1024.0 / 1024.0); | |
381 | ||
382 | MemAlloc((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19); | |
383 | MemAlloc((void **)&pdfs[1], sizeof(PdfT) * nCells * N_D3Q19); | |
384 | ||
385 | kd->Pdfs[0] = pdfs[0]; | |
386 | kd->Pdfs[1] = pdfs[1]; | |
387 | ||
388 | // Initialize PDFs with some (arbitrary) data for correct NUMA placement. | |
389 | // Here we touch only the fluid nodes as this loop is OpenMP parallel and | |
390 | // we want the same scheduling as in the kernel. | |
391 | #ifdef _OPENMP | |
392 | #pragma omp parallel for | |
393 | #endif | |
394 | for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) { | |
395 | pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0; | |
396 | pdfs[1][P_INDEX_3(nCells, i, d)] = 1.0; | |
397 | } } | |
398 | ||
399 | // Initialize all PDFs to some standard value. | |
400 | for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) { | |
401 | pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0; | |
402 | pdfs[1][P_INDEX_3(nCells, i, d)] = 0.0; | |
403 | } } | |
404 | ||
405 | // ---------------------------------------------------------------------- | |
406 | // create grid which will hold the index numbers of the fluid nodes | |
407 | ||
408 | uint32_t * grid; | |
409 | ||
410 | if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) { | |
411 | printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t)); | |
412 | exit(1); | |
413 | } | |
414 | kdl->Grid = grid; | |
415 | ||
416 | int latticeIndex; | |
417 | ||
418 | #ifdef DEBUG | |
419 | for(int z = 0; z < lZ; ++z) { | |
420 | for(int y = 0; y < lY; ++y) { | |
421 | for(int x = 0; x < lX; ++x) { | |
422 | ||
423 | latticeIndex = L_INDEX_4(ld->Dims, x, y, z); | |
424 | ||
425 | grid[latticeIndex] = ~0; | |
426 | } | |
427 | } | |
428 | } | |
429 | #endif | |
430 | ||
431 | // ---------------------------------------------------------------------- | |
432 | // generate numbering over grid | |
433 | ||
434 | uint32_t * coords; | |
435 | ||
436 | if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) { | |
437 | printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3); | |
438 | exit(1); | |
439 | } | |
440 | ||
441 | kdl->Coords = coords; | |
442 | ||
443 | // Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes | |
444 | // through the ld->Lattice array. | |
445 | int counter = 0; | |
446 | ||
447 | // Blocking is implemented via setup of the adjacency list. The kernel later will | |
448 | // walk through the lattice blocked automatically. | |
10988083 | 449 | for (int bX = 0; bX < lX; bX += blk[0]) { |
e3f82424 MW |
450 | for (int bY = 0; bY < lY; bY += blk[1]) { |
451 | for (int bZ = 0; bZ < lZ; bZ += blk[2]) { | |
10988083 MW |
452 | |
453 | int eX = MIN(bX + blk[0], lX); | |
454 | int eY = MIN(bY + blk[1], lY); | |
455 | int eZ = MIN(bZ + blk[2], lZ); | |
456 | ||
10988083 | 457 | for (int x = bX; x < eX; ++x) { |
e3f82424 MW |
458 | for (int y = bY; y < eY; ++y) { |
459 | for (int z = bZ; z < eZ; ++z) { | |
10988083 MW |
460 | |
461 | latticeIndex = L_INDEX_4(lDims, x, y, z); | |
462 | ||
463 | if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) { | |
464 | grid[latticeIndex] = counter; | |
465 | ||
466 | coords[C_INDEX_X(counter)] = x; | |
467 | coords[C_INDEX_Y(counter)] = y; | |
468 | coords[C_INDEX_Z(counter)] = z; | |
469 | ||
470 | ++counter; | |
471 | } | |
472 | } } } | |
473 | } } } | |
474 | ||
475 | Verify(counter == nFluid); | |
476 | ||
477 | uint32_t * adjList; | |
478 | ||
479 | // AdjList only requires 18 instead of 19 entries per node, as | |
480 | // the center PDF needs no addressing. | |
481 | if (MemAlloc((void **)&adjList, nFluid * sizeof(uint32_t) * N_D3Q19_IDX)) { | |
482 | printf("ERROR: allocating adjList array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * N_D3Q19_IDX); | |
483 | exit(1); | |
484 | } | |
485 | ||
486 | kdl->AdjList = adjList; | |
487 | ||
488 | int x, y, z; | |
489 | ||
490 | uint32_t neighborIndex; | |
491 | uint32_t dstIndex; | |
492 | ||
493 | int nx, ny, nz, px, py, pz; | |
494 | ||
495 | // Loop over all fluid nodes and compute the indices to the neighboring | |
496 | // PDFs for configured data layout (AoS/SoA). | |
e3f82424 MW |
497 | #ifdef _OPENMP |
498 | #pragma omp parallel for | |
499 | #endif | |
500 | for (int index = 0; index < nFluid; ++index) { | |
501 | for (int d = 0; d < N_D3Q19_IDX; ++d) { | |
502 | adjList[index * N_D3Q19_IDX + d] = -1; | |
503 | } | |
504 | } | |
505 | ||
10988083 MW |
506 | // #ifdef _OPENMP --> add line continuation |
507 | // #pragma omp parallel for default(none) | |
508 | // shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z, | |
509 | // stderr, | |
510 | // lDims, grid, ld, lX, lY, lZ, adjList) | |
511 | // private(x, y, z, nx, ny, nz, neighborIndex, dstIndex) | |
512 | // #endif | |
513 | for (int index = 0; index < nFluid; ++index) { | |
514 | x = coords[C_INDEX_X(index)]; | |
515 | y = coords[C_INDEX_Y(index)]; | |
516 | z = coords[C_INDEX_Z(index)]; | |
517 | ||
518 | Assert(x >= 0 && x < lX); | |
519 | Assert(y >= 0 && y < lY); | |
520 | Assert(z >= 0 && z < lZ); | |
521 | ||
522 | Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE); | |
523 | ||
524 | // Loop over all directions except the center one. | |
525 | for(int d = 0; d < N_D3Q19 - 1; ++d) { | |
526 | Assert(d != D3Q19_C); | |
527 | #ifdef PROP_MODEL_PUSH | |
528 | nx = x + D3Q19_X[d]; | |
529 | ny = y + D3Q19_Y[d]; | |
530 | nz = z + D3Q19_Z[d]; | |
531 | #elif PROP_MODEL_PULL | |
532 | nx = x - D3Q19_X[d]; | |
533 | ny = y - D3Q19_Y[d]; | |
534 | nz = z - D3Q19_Z[d]; | |
535 | #else | |
536 | #error No implementation for this PROP_MODEL_NAME. | |
537 | #endif | |
538 | // If the neighbor is outside the latcie in X direction and we have a | |
539 | // periodic boundary then we need to wrap around. | |
540 | if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) || | |
541 | ((ny < 0 || ny >= lY) && ld->PeriodicY) || | |
542 | ((nz < 0 || nz >= lZ) && ld->PeriodicZ) | |
543 | ){ | |
544 | // x periodic | |
545 | ||
546 | if (nx < 0) { | |
547 | px = lX - 1; | |
548 | } | |
549 | else if (nx >= lX) { | |
550 | px = 0; | |
551 | } else { | |
552 | px = nx; | |
553 | } | |
554 | // y periodic | |
555 | if (ny < 0) { | |
556 | py = lY - 1; | |
557 | } | |
558 | else if (ny >= lY) { | |
559 | py = 0; | |
560 | } else { | |
561 | py = ny; | |
562 | } | |
563 | ||
564 | // z periodic | |
565 | if (nz < 0) { | |
566 | pz = lZ - 1; | |
567 | } | |
568 | else if (nz >= lZ) { | |
569 | pz = 0; | |
570 | } else { | |
571 | pz = nz; | |
572 | } | |
573 | ||
574 | if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) { | |
575 | dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]); | |
576 | } | |
577 | else { | |
578 | neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)]; | |
579 | ||
580 | AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d); | |
581 | ||
582 | dstIndex = P_INDEX_3(nCells, neighborIndex, d); | |
583 | } | |
584 | } | |
585 | else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { | |
586 | dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]); | |
587 | } | |
588 | else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) { | |
589 | dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]); | |
590 | } | |
591 | else { | |
592 | neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)]; | |
593 | ||
594 | Assert(neighborIndex != ~0); | |
595 | ||
596 | dstIndex = P_INDEX_3(nCells, neighborIndex, d); | |
597 | } | |
598 | ||
599 | Assert(dstIndex >= 0); | |
600 | Assert(dstIndex < nCells * N_D3Q19); | |
601 | ||
602 | adjList[index * N_D3Q19_IDX + d] = dstIndex; | |
603 | } | |
604 | } | |
605 | ||
606 | ||
607 | // Fill remaining KernelData structures | |
608 | kd->GetNode = GetNode; | |
609 | kd->SetNode = SetNode; | |
610 | ||
611 | kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf); | |
612 | kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf); | |
613 | ||
614 | kd->Kernel = FNAME(D3Q19ListKernel); | |
615 | ||
616 | kd->DstPdfs = NULL; | |
617 | kd->PdfsActive = kd->Pdfs[0]; | |
618 | ||
619 | return; | |
620 | } | |
621 | ||
622 | void FNAME(D3Q19ListDeinit)(LatticeDesc * ld, KernelData ** kernelData) | |
623 | { | |
624 | KernelDataList ** kdl = (KernelDataList **)kernelData; | |
625 | ||
626 | MemFree((void **)&((*kernelData)->Pdfs[0])); | |
627 | MemFree((void **)&((*kernelData)->Pdfs[1])); | |
628 | ||
629 | MemFree((void **)&((*kdl)->AdjList)); | |
630 | MemFree((void **)&((*kdl)->Coords)); | |
631 | MemFree((void **)&((*kdl)->Grid)); | |
632 | ||
633 | MemFree((void **)kernelData); | |
634 | ||
635 | return; | |
636 | } | |
637 |