1 // --------------------------------------------------------------------------
4 // Markus Wittmann, 2016-2017
5 // RRZE, University of Erlangen-Nuremberg, Germany
6 // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
9 // LSS, University of Erlangen-Nuremberg, Germany
11 // Michael Hussnaetter, 2017-2018
12 // University of Erlangen-Nuremberg, Germany
13 // michael.hussnaetter -at- fau.de
15 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
17 // LbmBenchKernels is free software: you can redistribute it and/or modify
18 // it under the terms of the GNU General Public License as published by
19 // the Free Software Foundation, either version 3 of the License, or
20 // (at your option) any later version.
22 // LbmBenchKernels is distributed in the hope that it will be useful,
23 // but WITHOUT ANY WARRANTY; without even the implied warranty of
24 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
25 // GNU General Public License for more details.
27 // You should have received a copy of the GNU General Public License
28 // along with LbmBenchKernels. If not, see <http://www.gnu.org/licenses/>.
30 // --------------------------------------------------------------------------
31 #include "BenchKernelD3Q19ListAaPvGatherCommon.h"
45 #if ALLOC_ADJ_LIST_IN_HBM == 1
46 #define ADJ_LIST_ALLOCATOR HbwAllocAligned
47 #define ADJ_LIST_FREE HbwFree
49 #define ADJ_LIST_ALLOCATOR MemAllocAligned
50 #define ADJ_LIST_FREE MemFree
53 #if ALLOC_PDF_IN_HBM == 1
54 #define PDF_ALLOCATOR HbwAllocAligned
55 #define PDF_FREE HbwFree
57 #define PDF_ALLOCATOR MemAllocAligned
58 #define PDF_FREE MemFree
61 // Forward definition.
62 void FNAME(D3Q19ListAaPvGatherKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
65 // -----------------------------------------------------------------------
66 // Functions which are used as callback by the kernel to read or write
70 static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
73 Assert(kd->PdfsActive != NULL);
74 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
77 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
78 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
79 Assert(dir >= 0); Assert(dir < N_D3Q19);
81 KernelDataList * kdl = KDL(kd);
82 uint32_t * adjList = kdl->AdjList;
84 if (kdl->Iteration % 2 == 0) {
85 // Pdfs are stored inverse, local PDFs are located in remote nodes
88 uint32_t index = kdl->Grid[L_INDEX_4(kd->Dims, x, y, z)];
91 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
92 *pdf = kd->PdfsActive[ADJ_LIST(D3Q19_INV[dir])];
96 *pdf = kd->PdfsActive[P_INDEX_3(kdl->nCells, index, dir)];
101 *pdf = kd->PdfsActive[P_INDEX_5(kdl, x, y, z, dir)];
107 static void FNAME(BCSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
110 Assert(kd->PdfsActive != NULL);
111 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
112 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
113 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
114 Assert(dir >= 0); Assert(dir < N_D3Q19);
117 printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]);
122 KernelDataList * kdl = KDL(kd);
123 uint32_t * adjList = kdl->AdjList;
125 if (kdl->Iteration % 2 == 0) {
126 // Pdfs are stored inverse, local PDFs are located in remote nodes
128 // getting node index
129 uint32_t index = kdl->Grid[L_INDEX_4(kd->Dims, x, y, z)];
131 if (dir != D3Q19_C) {
132 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
133 kd->PdfsActive[ADJ_LIST(D3Q19_INV[dir])] = pdf;
136 kd->PdfsActive[P_INDEX_3(kdl->nCells, index, dir)] = pdf;
140 kd->PdfsActive[P_INDEX_5(kdl, x, y, z, dir)] = pdf;
147 static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
150 Assert(kd->PdfsActive != NULL);
151 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
152 Assert(pdfs != NULL);
153 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
154 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
156 KernelDataList * kdl = KDL(kd);
157 uint32_t * adjList = kdl->AdjList;
159 if(kdl->Iteration % 2 == 0){
161 uint32_t index = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
163 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
164 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
165 #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[ADJ_LIST(idxinv)];
169 pdfs[D3Q19_C] = kd->PdfsActive[P_INDEX_3(kdl->nCells, index, D3Q19_C)];
173 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
174 #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)];
181 for (int d = 0; d < 19; ++d) {
182 if(isnan(pdfs[d]) || isinf(pdfs[d])) {
183 printf("%d %d %d %d nan! get node\n", x, y, z, d);
184 for (int d2 = 0; d2 < 19; ++d2) {
185 printf("%d: %e\n", d2, pdfs[d2]);
195 static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
198 Assert(kd->PdfsActive != NULL);
199 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
200 Assert(pdfs != NULL);
202 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
203 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
205 for (int d = 0; d < 19; ++d) {
207 printf("%d %d %d %d nan! get node\n", x, y, z, d);
208 for (int d2 = 0; d2 < 19; ++d2) {
209 printf("%d: %e\n", d2, pdfs[d2]);
215 KernelDataList * kdl = KDL(kd);
216 uint32_t * adjList = kdl->AdjList;
218 if(kdl->Iteration % 2 == 0){
220 uint32_t index = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
222 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
223 kd->PdfsActive[P_INDEX_3(kdl->nCells, index, D3Q19_C)] = pdfs[D3Q19_C];
225 #define ADJ_LIST(dir) adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (dir * VSIZE) + (index % VSIZE)]
226 #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[ADJ_LIST(idxinv)] = pdfs[idx];
233 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
234 #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
244 static void ParameterUsage()
246 printf("Kernel parameters:\n");
247 printf(" [-blk <n>] [-blk-[xyz] <n>]\n");
252 static void ParseParameters(Parameters * params, int * blk)
256 blk[0] = 0; blk[1] = 0; blk[2] = 0;
258 #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param))
259 #define NEXT_ARG_PRESENT() \
261 if (i + 1 >= params->nKernelArgs) { \
262 printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
268 for (int i = 0; i < params->nKernelArgs; ++i) {
269 if (ARG_IS("-blk") || ARG_IS("--blk")) {
272 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
275 printf("ERROR: blocking parameter must be > 0.\n");
279 blk[0] = blk[1] = blk[2] = tmp;
281 else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
284 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
287 printf("ERROR: blocking parameter must be > 0.\n");
293 else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
296 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
299 printf("ERROR: blocking parameter must be > 0.\n");
305 else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
308 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
311 printf("ERROR: blocking parameter must be > 0.\n");
317 else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
322 printf("ERROR: unknown kernel parameter.\n");
329 #undef NEXT_ARG_PRESENT
334 static void SetupConsecNodes(LatticeDesc * ld, KernelDataListRia * kdlr, int nThreads)
337 Assert(kdlr != NULL);
338 Assert(nThreads > 0);
340 uint32_t * adjList = kdlr->kdl.AdjList;
342 uint32_t nConsecNodes = 0;
343 uint32_t consecIndex = 0;
345 int nFluid = kdlr->kdl.nFluid;
347 uint32_t * consecThreadIndices = (uint32_t *)malloc(sizeof(uint32_t) * (nThreads + 1));
349 int nNodesPerThread = nFluid / nThreads;
351 for (int i = 0; i < nThreads; ++i) {
352 consecThreadIndices[i] = i * nNodesPerThread + MinI(i, nFluid % nThreads);
354 consecThreadIndices[nThreads] = nFluid;
357 int similarPatterns = 0;
358 int wasLastChunkThreadBoundary = 0;
359 // We execute following code two times.
360 // - The first time to get the count of how many entries we need for the
361 // consecNodes array.
362 // - The second time to fill the array.
364 // Loop over adjacency list of all nodes.
365 // Compare if adjacent nodes share the same access pattern.
366 for (int fluidBaseIndex = VSIZE; fluidBaseIndex < nFluid; fluidBaseIndex += VSIZE) {
368 int hasSimilarAccessPattern = 1;
370 // Loop over all directions except the center one.
371 for(int d = 0; d < N_D3Q19 - 1; ++d) {
372 Assert(d != D3Q19_C);
374 // check if cache line itself has consecutive memory access pattern
375 for(int inChunkIndex = 0; (inChunkIndex < VSIZE) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
376 int index = fluidBaseIndex + inChunkIndex;
378 Assert(index < nFluid);
380 #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
381 //if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1)
382 if (ADJ_LIST(index, d) != ADJ_LIST(index-VSIZE, d) + VSIZE) {
383 //printf("different @: ADJ_LST(%d,%d)=%d != %d=ADJ_LST(%d, %d) + VSIZE\n", index, d, ADJ_LIST(index,d), ADJ_LIST(index-VSIZE,d) + VSIZE, index-VSIZE, d);
384 // Different access pattern.
385 hasSimilarAccessPattern = 0;
391 if(!hasSimilarAccessPattern){
392 break; //exit from nested loop
396 long threadBoundaryIndex = consecThreadIndices[indexThread];
397 if (fluidBaseIndex <= threadBoundaryIndex &&
398 threadBoundaryIndex < fluidBaseIndex + VSIZE) {
399 // Current chunk contains thread boundary.
400 // These chunks are treated by scalar peel and reminder loops
401 // in kernel of every thread to ensure VSIZE aligned access to
404 // final cells of current thread
407 // first cells of next thread
411 wasLastChunkThreadBoundary = 1;
414 // We are not at a thread boundary
415 if (hasSimilarAccessPattern && !wasLastChunkThreadBoundary){
422 wasLastChunkThreadBoundary = 0;
425 if (!hasSimilarAccessPattern) {
436 nConsecNodes = consecIndex + 1;
439 uint32_t * consecNodes;
440 MemAlloc((void **)&consecNodes, sizeof(uint32_t) * nConsecNodes);
442 unsigned long consecNodesByte = (nConsecNodes) * sizeof(uint32_t);
444 printf("# Consec. Nodes Array Allocation:\n");
445 printf("# similar patterns\t\t%d\n", similarPatterns);
446 printf("# elements: \t\t%d\n", nConsecNodes);
447 printf("# size: \t\t%e MiB\n", consecNodesByte / 1024.0 / 1024.0);
448 printf("# alignment: \t\t%d b\n", PAGE_4K);
450 if (MemAllocAligned((void **)&consecNodesByte, consecNodesByte, PAGE_4K)) {
451 printf("ERROR: allocating consecNodes array with MemAllocAligned failed: %lu bytes.\n", consecNodesByte);
455 printf("# allocator: \t\t\tMemAllocAligned()\n");
461 consecNodes[consecIndex] = VSIZE;
465 consecThreadIndices[0] = 0;
467 //add first chunk manually to enable backward check for consecutive pattern
468 consecNodes[consecIndex] = VSIZE;
470 wasLastChunkThreadBoundary = 0;
472 // Loop over adjacency list of all nodes.
473 // Compare if access pattern does not change on chunk level
474 // Since gather instructions are used, access pattern may not be consecutive
475 for (int fluidBaseIndex = VSIZE; fluidBaseIndex < nFluid; fluidBaseIndex += VSIZE) {
477 int hasSimilarAccessPattern = 1;
479 // Loop over all directions except the center one.
480 for(int d = 0; d < N_D3Q19 - 1; ++d) {
481 Assert(d != D3Q19_C);
483 // check if cache line itself has consecutive memory access pattern
484 for(int inChunkIndex = 0; (inChunkIndex < VSIZE) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
485 int index = fluidBaseIndex + inChunkIndex;
487 Assert(index < nFluid);
489 #define ADJ_LIST(idx, dir) adjList[((idx) - ((idx) % VSIZE)) * N_D3Q19_IDX + ((dir) * VSIZE) + ((idx) % VSIZE)]
490 //if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1)
491 if (ADJ_LIST(index, d) != ADJ_LIST(index-VSIZE, d) + VSIZE) {
492 // Different access pattern.
493 hasSimilarAccessPattern = 0;
499 if(!hasSimilarAccessPattern){
500 break; //exit from nested loop
504 long threadBoundaryIndex = consecThreadIndices[indexThread];
505 if (fluidBaseIndex <= threadBoundaryIndex &&
506 threadBoundaryIndex < fluidBaseIndex + VSIZE) {
507 // Current chunk contains thread boundary.
508 // These chunks are treated by scalar peel and reminder loops
509 // in kernel of every thread to ensure VSIZE aligned access to
512 // final cells of current thread
514 //consecThreadIndices[indexThread] = consecIndex;
515 consecNodes[consecIndex] = threadBoundaryIndex - fluidBaseIndex;
518 // first cells of next thread
520 consecThreadIndices[indexThread] = consecIndex;
521 consecNodes[consecIndex] = (fluidBaseIndex + VSIZE) - threadBoundaryIndex;
524 wasLastChunkThreadBoundary = 1;
528 // We are not at a thread boundary
529 if (hasSimilarAccessPattern && !wasLastChunkThreadBoundary){
530 Assert(consecIndex < nConsecNodes);
531 consecNodes[consecIndex] += VSIZE;
535 Assert(consecIndex < nConsecNodes);
536 consecNodes[consecIndex] = VSIZE;
540 if (!hasSimilarAccessPattern) {
542 Assert(consecIndex < nConsecNodes);
543 consecNodes[consecIndex] = VSIZE;
546 Assert(consecIndex < nConsecNodes);
547 consecNodes[consecIndex] += VSIZE;
550 wasLastChunkThreadBoundary = 0;
556 printf("consecNodes:\n");
557 for(int i = 0; i < nConsecNodes + 5; ++i){
558 printf("%d ", consecNodes[i]);
563 printf("consecThreadIndices:\n");
564 for(int i = 0; i < nThreads + 5; ++i){
565 printf("%d ", consecThreadIndices[i]);
570 kdlr->ConsecNodes = consecNodes;
571 kdlr->nConsecNodes = nConsecNodes;
573 kdlr->ConsecThreadIndices = consecThreadIndices;
574 kdlr->nConsecThreadIndices = nThreads;
576 double loopBalanceEven = 2.0 * 19 * sizeof(PdfT);
577 //N_D3Q19 - 1: C lookup not required, +1: transfer of consecValue
578 double loopBalanceOdd = 2.0 * 19 * sizeof(PdfT) + ((double)nConsecNodes *((N_D3Q19 - 1) * VSIZE + 1)) / nFluid * sizeof(int);
579 double loopBalance = (loopBalanceEven + loopBalanceOdd) / 2.0;
581 kdlr->kdl.kd.LoopBalance = loopBalance;
583 printf("# loop balance:\n");
584 printf("# even timestep: \t\t%.2f B/FLUP\n", loopBalanceEven);
585 printf("# odd timestep: \t\t%.2f B/FLUP\n", loopBalanceOdd);
586 printf("# average: \t\t%.2f B/FLUP\n", loopBalance);
591 void FNAME(D3Q19ListAaPvGatherInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
594 KernelDataList * kdl;
595 KernelDataListRia * kdlr;
596 MemAlloc((void **)&kdlr, sizeof(KernelDataListRia));
598 kd = (KernelData *)kdlr;
606 kd->PdfsActive = NULL;
612 kd->GlobalDims[0] = -1;
613 kd->GlobalDims[1] = -1;
614 kd->GlobalDims[2] = -1;
619 kd->ObstIndices = NULL;
620 kd->nObstIndices = -1;
621 kd->BounceBackPdfsSrc = NULL;
622 kd->BounceBackPdfsDst = NULL;
623 kd->nBounceBackPdfs = -1;
631 kdlr->ConsecNodes = NULL;
632 kdlr->nConsecNodes = 0;
633 kdlr->ConsecThreadIndices = NULL;
634 kdlr->nConsecThreadIndices = 0;
637 // Ajust the dimensions according to padding, if used.
638 kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
639 kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
640 kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];
642 int * lDims = ld->Dims;
648 int nTotalCells = lX * lY * lZ;
649 int nCells = ld->nFluid; // TODO: + padding
650 int nFluid = ld->nFluid;
652 // TODO: check nCells/nFluid do not exceed 2^31. This actually has to be
653 // done during lattice setup.
654 kdl->nCells = nCells;
655 kdl->nFluid = nFluid;
661 ParseParameters(params, blk);
663 if (blk[0] == 0) blk[0] = lX;
664 if (blk[1] == 0) blk[1] = lY;
665 if (blk[2] == 0) blk[2] = lZ;
667 printf("# blocking: \t\tx: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
669 unsigned long latByte = nCells * sizeof(PdfT) * N_D3Q19;
670 unsigned long latFluidByte = nFluid * sizeof(PdfT) * N_D3Q19;
671 unsigned long latPadByte = (nCells - nFluid) * sizeof(PdfT) * N_D3Q19;
673 printf("# Lattice Array Allocation:\n");
674 printf("# lattice size: \t\t%e MiB\n", latByte / 1024.0 / 1024.0);
675 printf("# fluid lattice size:\t\t%e MiB\n", latFluidByte / 1024.0 / 1024.0);
676 printf("# lattice padding: \t\t%e MiB\n", latPadByte / 1024.0 / 1024.0);
679 printf("# alignment: \t\t%d b\n", PAGE_4K);
681 if (PDF_ALLOCATOR((void **)&pdfs[0], latFluidByte, PAGE_4K)) {
682 printf("ERROR: allocating PDF array with %s() failed: %lu bytes.\n", STRINGIFY(PDF_ALLOCATOR), latFluidByte);
686 printf("# allocator: \t\t\t%s()\n", STRINGIFY(PDF_ALLOCATOR));
689 kd->Pdfs[0] = pdfs[0];
691 // Initialize PDFs with some (arbitrary) data for correct NUMA placement.
692 // Here we touch only the fluid nodes as this loop is OpenMP parallel and
693 // we want the same scheduling as in the kernel.
695 #pragma omp parallel for
697 for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
698 pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
701 // Initialize all PDFs to some standard value.
702 for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
703 pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
706 // ----------------------------------------------------------------------
707 // create grid which will hold the index numbers of the fluid nodes
711 if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) {
712 printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
720 for(int z = 0; z < lZ; ++z) {
721 for(int y = 0; y < lY; ++y) {
722 for(int x = 0; x < lX; ++x) {
724 latticeIndex = L_INDEX_4(ld->Dims, x, y, z);
726 grid[latticeIndex] = ~0;
732 // ----------------------------------------------------------------------
733 // generate numbering over grid
737 if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) {
738 printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
742 kdl->Coords = coords;
744 // Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
745 // through the ld->Lattice array.
748 // Blocking is implemented via setup of the adjacency list. The kernel later will
749 // walk through the lattice blocked automatically.
750 for (int bZ = 0; bZ < lZ; bZ += blk[2]) {
751 for (int bY = 0; bY < lY; bY += blk[1]) {
752 for (int bX = 0; bX < lX; bX += blk[0]) {
754 int eX = MIN(bX + blk[0], lX);
755 int eY = MIN(bY + blk[1], lY);
756 int eZ = MIN(bZ + blk[2], lZ);
759 for (int z = bZ; z < eZ; ++z) {
760 for (int y = bY; y < eY; ++y) {
761 for (int x = bX; x < eX; ++x) {
763 latticeIndex = L_INDEX_4(lDims, x, y, z);
765 if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
766 grid[latticeIndex] = counter;
768 coords[C_INDEX_X(counter)] = x;
769 coords[C_INDEX_Y(counter)] = y;
770 coords[C_INDEX_Z(counter)] = z;
777 Verify(counter == nFluid);
780 // AoSoA addressing for adjList needs padding for (nFluid % VSIZE) != 0
781 unsigned long adjListBytes = nFluid * sizeof(int) * N_D3Q19_IDX;
783 printf("# Adjacency List Allocation:\n");
784 printf("# size: \t\t%e MiB\n", adjListBytes / 1024.0 / 1024.0);
785 printf("# alignment: \t\t%d b\n", PAGE_4K);
787 // AdjList only requires 18 instead of 19 entries per node, as
788 // the center PDF needs no addressing.
789 if (ADJ_LIST_ALLOCATOR((void **)&adjList, adjListBytes, PAGE_4K)) {
790 printf("ERROR: allocating adjList array with %s() failed: %lu bytes.\n", STRINGIFY(ADJ_LIST_ALLOCATOR), adjListBytes);
794 printf("# allocator: \t\t\t%s()\n", STRINGIFY(ADJ_LIST_ALLOCATOR));
797 kdl->AdjList = adjList;
801 uint32_t neighborIndex;
804 int nx, ny, nz, px, py, pz;
806 // Loop over all fluid nodes and compute the indices to the neighboring
807 // PDFs for configured data layout (AoS/SoA).
808 // Parallelized loop to ensure correct NUMA placement.
809 // #ifdef _OPENMP --> add line continuation
810 // #pragma omp parallel for default(none)
811 // shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
813 // lDims, grid, ld, lX, lY, lZ, adjList)
814 // private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
816 for (int fluidBaseIndex = 0; fluidBaseIndex < nFluid; fluidBaseIndex+=VSIZE) {
819 // Loop over all directions except the center one.
820 for(int d = 0; d < N_D3Q19 - 1; ++d) {
821 Assert(d != D3Q19_C);
823 for(int inChunkIndex = 0; (inChunkIndex < VSIZE) && ((fluidBaseIndex + inChunkIndex) < nFluid); ++inChunkIndex){
824 int index = fluidBaseIndex + inChunkIndex;
826 Assert(index < nFluid);
828 x = coords[C_INDEX_X(index)];
829 y = coords[C_INDEX_Y(index)];
830 z = coords[C_INDEX_Z(index)];
832 Assert(x >= 0 && x < lX);
833 Assert(y >= 0 && y < lY);
834 Assert(z >= 0 && z < lZ);
836 Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);
838 #ifdef PROP_MODEL_PUSH
843 #elif PROP_MODEL_PULL
848 #error No implementation for this PROP_MODEL_NAME.
850 // If the neighbor is outside the latice in X direction and we have a
851 // periodic boundary then we need to wrap around.
852 if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) ||
853 ((ny < 0 || ny >= lY) && ld->PeriodicY) ||
854 ((nz < 0 || nz >= lZ) && ld->PeriodicZ)
886 if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
887 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
890 neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];
892 AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);
894 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
897 else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
898 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
900 else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
901 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
904 neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];
906 Assert(neighborIndex != ~0);
908 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
911 Assert(dstIndex >= 0);
912 Assert(dstIndex < nCells * N_D3Q19);
914 adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (d * VSIZE) + (index % VSIZE)] = dstIndex;
920 printf("============\n");
921 for(int baseIndex = 0; baseIndex < nFluid; baseIndex+=VSIZE){
922 for(int i = 0; i < VSIZE; ++i){
923 int index = baseIndex + i;
925 printf("%d ", adjList[(index - (index % VSIZE)) * N_D3Q19_IDX + (0 * VSIZE) + (index % VSIZE)]);
929 printf("============\n");
935 nThreads = omp_get_max_threads();
938 SetupConsecNodes(ld, KDLR(kd), nThreads);
940 // Fill remaining KernelData structures
941 kd->GetNode = GetNode;
942 kd->SetNode = SetNode;
944 kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
945 kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);
947 kd->Kernel = FNAME(D3Q19ListAaPvGatherKernel);
950 kd->PdfsActive = kd->Pdfs[0];
955 void FNAME(D3Q19ListAaPvGatherDeinit)(LatticeDesc * ld, KernelData ** kernelData)
957 KernelDataListRia ** kdlr = (KernelDataListRia **)kernelData;
959 MemFree((void **)&((*kdlr)->ConsecNodes));
961 if ((*kdlr)->ConsecThreadIndices != NULL) {
962 MemFree((void **)&((*kdlr)->ConsecThreadIndices));
965 KernelDataList ** kdl = (KernelDataList **)kernelData;
967 ADJ_LIST_FREE((void **)&((*kdl)->AdjList));
969 MemFree((void **)&((*kdl)->Coords));
970 MemFree((void **)&((*kdl)->Grid));
972 PDF_FREE((void **)&((*kernelData)->Pdfs[0]));
974 MemFree((void **)kernelData);