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 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
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.
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.
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/>.
26 // --------------------------------------------------------------------------
27 #include "BenchKernelD3Q19ListPullSplitNtCommon.h"
41 // Forward definition.
42 void FNAME(KernelPullSplitNt1S)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
43 void FNAME(KernelPullSplitNt2S)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
45 void FNAME(KernelPullSplitNtRia1S)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
46 void FNAME(KernelPullSplitNtRia2S)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
51 // -----------------------------------------------------------------------
52 // Functions which are used as callback by the kernel to read or write
55 static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
58 Assert(kd->PdfsActive != NULL);
59 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
62 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
63 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
64 Assert(dir >= 0); Assert(dir < N_D3Q19);
66 // The relevant PDFs here are the ones, which will get streamed in later
67 // during propagation. So we must return the *remote* PDFs.
68 uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];
72 uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
74 *pdf = kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + dir]];
77 *pdf = kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)];
84 static void FNAME(BCSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
87 Assert(kd->PdfsActive != NULL);
88 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
89 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
90 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
91 Assert(dir >= 0); Assert(dir < N_D3Q19);
94 printf("ERROR: setting nan %d %d %d %d %s\n", x, y, z, dir, D3Q19_NAMES[dir]);
99 // The relevant PDFs here are the ones, which will get streamed in later
100 // during propagation. So we must set this *remote* PDFs.
101 uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];
103 if (dir != D3Q19_C) {
105 uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
107 kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + dir]] = pdf;
110 kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)] = pdf;
118 static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
121 Assert(kd->PdfsActive != NULL);
122 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
123 Assert(pdfs != NULL);
124 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
125 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
129 // TODO: pull scheme?
131 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
132 #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)]; sum += pdfs[idx];
138 for (int d = 0; d < 19; ++d) {
139 if(isnan(pdfs[d]) || isinf(pdfs[d])) {
140 printf("%d %d %d %d nan! get node\n", x, y, z, d);
141 for (int d2 = 0; d2 < 19; ++d2) {
142 printf("%d: %e\n", d2, pdfs[d2]);
152 static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
155 Assert(kd->PdfsActive != NULL);
156 Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
157 Assert(pdfs != NULL);
159 Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
160 Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
163 for (int d = 0; d < 19; ++d) {
165 printf("%d %d %d %d nan! get node\n", x, y, z, d);
166 for (int d2 = 0; d2 < 19; ++d2) {
167 printf("%d: %e\n", d2, pdfs[d2]);
174 // TODO: pull scheme?
175 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
176 #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
184 static void ParameterUsage()
186 printf("Kernel parameters:\n");
187 printf(" [-blk <n>] [-blk-[xyz] <n>] [-n-tmp-array <n>]\n");
188 printf(" [-pad auto|modulus_1+offset_1(,modulus_n+offset_n)*]\n");
193 static void ParseParameters(Parameters * params, int * blk, int * nTmpArray, PadInfo ** padInfo)
197 blk[0] = 0; blk[1] = 0; blk[2] = 0;
198 *nTmpArray = 152 - (152 % VSIZE);
201 #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param))
202 #define NEXT_ARG_PRESENT() \
204 if (i + 1 >= params->nKernelArgs) { \
205 printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
211 for (int i = 0; i < params->nKernelArgs; ++i) {
212 if (ARG_IS("-blk") || ARG_IS("--blk")) {
215 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
218 printf("ERROR: blocking parameter must be >= 0.\n");
222 blk[0] = blk[1] = blk[2] = tmp;
224 else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
227 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
230 printf("ERROR: blocking parameter must be >= 0.\n");
236 else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
239 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
242 printf("ERROR: blocking parameter must be >= 0.\n");
248 else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
251 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
254 printf("ERROR: blocking parameter must be >= 0.\n");
260 else if (ARG_IS("-n-tmp-array") || ARG_IS("--n-tmp-array")) {
263 int tmp = strtol(params->KernelArgs[++i], NULL, 0);
266 printf("ERROR: -n-tmp-array parameter must be > 0.\n");
270 if (tmp % VSIZE != 0) {
271 printf("ERROR: value for -n-tmp-array must be a multiple of %d.\n", VSIZE);
277 else if (ARG_IS("-pad") || ARG_IS("--pad")) {
280 *padInfo = PadInfoFromStr(params->KernelArgs[++i]);
282 else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
287 printf("ERROR: unknown kernel parameter.\n");
294 #undef NEXT_ARG_PRESENT
299 static void SetupConsecNodes(LatticeDesc * ld, KernelDataListRia * kdlr, int nThreads)
302 Assert(kdlr != NULL);
303 Assert(nThreads > 0);
305 uint32_t * adjList = kdlr->kdl.AdjList;
307 uint32_t nConsecNodes = 0;
308 uint32_t consecIndex = 0;
310 int nFluid = kdlr->kdl.nFluid;
312 uint32_t * consecThreadIndices = (uint32_t *)malloc(sizeof(uint32_t) * (nThreads + 1));
314 int nNodesPerThread = nFluid / nThreads;
316 for (int i = 0; i < nThreads; ++i) {
317 consecThreadIndices[i] = i * nNodesPerThread + MinI(i, nFluid % nThreads);
319 consecThreadIndices[nThreads] = -1;
323 // We execute following code two times.
324 // - The first time to get the count of how many entries we need for the
325 // consecNodes array.
326 // - The second time to fill the array.
328 // Loop over adjacency list of all nodes.
329 // Compare if adjacent nodes share the same access pattern.
330 for (int index = 1; index < nFluid; ++index) {
334 // Loop over all directions except the center one.
335 for(int d = 0; d < N_D3Q19 - 1; ++d) {
336 Assert(d != D3Q19_C);
338 if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1) {
339 // Different access pattern.
345 if (consecThreadIndices[indexThread] == index) {
346 // We are at a thread boundary. Starting from this index the fluids
347 // belong to another thread. Force a break, if nodes are consecutive.
358 nConsecNodes = consecIndex + 1;
361 uint32_t * consecNodes;
362 MemAlloc((void **)&consecNodes, sizeof(uint32_t) * nConsecNodes);
367 consecNodes[consecIndex] = 1;
371 consecThreadIndices[0] = 0;
373 // Loop over adjacency list of all nodes.
374 // Compare if adjacent nodes share the same access pattern.
375 for (int index = 1; index < nFluid; ++index) {
379 // Loop over all directions except the center one.
380 for(int d = 0; d < N_D3Q19 - 1; ++d) {
381 Assert(d != D3Q19_C);
383 if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1) {
384 // Different access pattern.
390 if (consecThreadIndices[indexThread] == index) {
391 // We are at a thread boundary. Starting from this index the fluids
392 // belong to another thread. Force a break, if nodes are consecutive.
393 consecThreadIndices[indexThread] = consecIndex + 1;
400 Assert(consecIndex < nConsecNodes);
401 consecNodes[consecIndex] = 1;
404 Assert(consecIndex < nConsecNodes);
405 consecNodes[consecIndex] += 1;
410 kdlr->ConsecNodes = consecNodes;
411 kdlr->nConsecNodes = nConsecNodes;
413 kdlr->ConsecThreadIndices = consecThreadIndices;
414 kdlr->nConsecThreadIndices = nThreads;
416 // printf("# total fluid nodes: %d consecutive blocks: %d\n", nFluid, nConsecNodes);
422 static void FNAME(Init)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
425 KernelDataList * kdl;
426 KernelDataListRia * kdlr;
427 MemAlloc((void **)&kdlr, sizeof(KernelDataListRia));
429 kd = (KernelData *)kdlr;
437 kd->PdfsActive = NULL;
443 kd->GlobalDims[0] = -1;
444 kd->GlobalDims[1] = -1;
445 kd->GlobalDims[2] = -1;
450 kd->ObstIndices = NULL;
451 kd->nObstIndices = -1;
452 kd->BounceBackPdfsSrc = NULL;
453 kd->BounceBackPdfsDst = NULL;
454 kd->nBounceBackPdfs = -1;
462 kdlr->ConsecNodes = NULL;
463 kdlr->nConsecNodes = 0;
464 kdlr->ConsecThreadIndices = NULL;
465 kdlr->nConsecThreadIndices = 0;
469 PadInfo * padInfo = NULL;
471 ParseParameters(params, blk, &kdlr->nTmpArray, &padInfo);
473 // Ajust the dimensions according to padding, if used.
474 kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
475 kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
476 kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];
478 int * lDims = ld->Dims;
484 int nTotalCells = lX * lY * lZ;
485 int nCells = ld->nFluid;
486 int nFluid = ld->nFluid;
489 nCells = PadCellsAndReport(nCells, sizeof(PdfT), &padInfo);
490 PadInfoFree(padInfo); padInfo = NULL;
493 // We padd each stream of a PDF array for a complete cache line.
494 int nElementsPerCl = 64 / sizeof(PdfT);
496 nCells = nCells + (nElementsPerCl - nCells % nElementsPerCl);
498 Assert(nCells % VSIZE == 0);
500 kdl->nCells = nCells;
501 kdl->nFluid = nFluid;
505 if (blk[0] == 0) blk[0] = lX;
506 if (blk[1] == 0) blk[1] = lY;
507 if (blk[2] == 0) blk[2] = lZ;
509 printf("# blocking x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
510 printf("# temporary array size: %d PDFs, %lu b\n", kdlr->nTmpArray, kdlr->nTmpArray * sizeof(PdfT) * 23);
512 double latMiB = nCells * sizeof(PdfT) * N_D3Q19 / 1024.0 / 1024.0;
513 double latFluidMib = nFluid * sizeof(PdfT) * N_D3Q19 / 1024.0 / 1024.0;
514 double latPadMib = (nCells - nFluid) * sizeof(PdfT) * N_D3Q19 / 1024.0 / 1024.0;
516 printf("# lattice size: %e MiB total: %e MiB\n", latMiB, latMiB * 2);
517 printf("# fluid lattice size: %e MiB total: %e MiB\n", latFluidMib, latFluidMib * 2);
518 printf("# lattice padding: %e MiB total: %e MiB\n", latPadMib, latPadMib * 2);
522 printf("# aligning lattices to: %d b\n", PAGE_4K);
524 MemAllocAligned((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19, PAGE_4K);
525 MemAllocAligned((void **)&pdfs[1], sizeof(PdfT) * nCells * N_D3Q19, PAGE_4K);
527 kd->Pdfs[0] = pdfs[0];
528 kd->Pdfs[1] = pdfs[1];
530 // Initialize PDFs with some (arbitrary) data for correct NUMA placement.
531 // Here we touch only the fluid nodes as this loop is OpenMP parallel and
532 // we want the same scheduling as in the kernel.
534 #pragma omp parallel for
536 for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
537 pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
538 pdfs[1][P_INDEX_3(nCells, i, d)] = 1.0;
541 // Initialize all PDFs to some standard value.
542 for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
543 pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
544 pdfs[1][P_INDEX_3(nCells, i, d)] = 0.0;
547 // ----------------------------------------------------------------------
548 // create grid which will hold the index numbers of the fluid nodes
552 if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) {
553 printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
561 for(int z = 0; z < lZ; ++z) {
562 for(int y = 0; y < lY; ++y) {
563 for(int x = 0; x < lX; ++x) {
565 latticeIndex = L_INDEX_4(ld->Dims, x, y, z);
567 grid[latticeIndex] = ~0;
573 // ----------------------------------------------------------------------
574 // generate numbering over grid
578 if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) {
579 printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
583 kdl->Coords = coords;
585 // Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
586 // through the ld->Lattice array.
589 // Blocking is implemented via setup of the adjacency list. The kernel later will
590 // walk through the lattice blocked automatically.
591 for (int bX = 0; bX < lX; bX += blk[0]) {
592 for (int bY = 0; bY < lY; bY += blk[1]) {
593 for (int bZ = 0; bZ < lZ; bZ += blk[2]) {
595 int eX = MIN(bX + blk[0], lX);
596 int eY = MIN(bY + blk[1], lY);
597 int eZ = MIN(bZ + blk[2], lZ);
600 for (int x = bX; x < eX; ++x) {
601 for (int y = bY; y < eY; ++y) {
602 for (int z = bZ; z < eZ; ++z) {
604 latticeIndex = L_INDEX_4(lDims, x, y, z);
606 if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
607 grid[latticeIndex] = counter;
609 coords[C_INDEX_X(counter)] = x;
610 coords[C_INDEX_Y(counter)] = y;
611 coords[C_INDEX_Z(counter)] = z;
618 Verify(counter == nFluid);
622 double indexMib = nFluid * sizeof(uint32_t) * N_D3Q19_IDX / 1024.0 / 1024.0;
624 printf("# index size: %e MiB\n", indexMib);
627 // AdjList only requires 18 instead of 19 entries per node, as
628 // the center PDF needs no addressing.
629 if (MemAlloc((void **)&adjList, nFluid * sizeof(uint32_t) * N_D3Q19_IDX)) {
630 printf("ERROR: allocating adjList array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * N_D3Q19_IDX);
634 kdl->AdjList = adjList;
638 uint32_t neighborIndex;
641 int nx, ny, nz, px, py, pz;
643 // Loop over all fluid nodes and compute the indices to the neighboring
644 // PDFs for configured data layout (AoS/SoA).
645 // Parallelized loop to ensure correct NUMA placement.
647 #pragma omp parallel for
649 for (int index = 0; index < nFluid; ++index) {
650 for (int d = 0; d < N_D3Q19_IDX; ++d) {
651 adjList[index * N_D3Q19_IDX + d] = -1;
655 // #ifdef _OPENMP --> add line continuation
656 // #pragma omp parallel for default(none)
657 // shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
659 // lDims, grid, ld, lX, lY, lZ, adjList)
660 // private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
662 for (int index = 0; index < nFluid; ++index) {
663 x = coords[C_INDEX_X(index)];
664 y = coords[C_INDEX_Y(index)];
665 z = coords[C_INDEX_Z(index)];
667 Assert(x >= 0 && x < lX);
668 Assert(y >= 0 && y < lY);
669 Assert(z >= 0 && z < lZ);
671 Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);
673 // Loop over all directions except the center one.
674 for(int d = 0; d < N_D3Q19 - 1; ++d) {
675 Assert(d != D3Q19_C);
677 #ifdef PROP_MODEL_PUSH
682 #elif PROP_MODEL_PULL
687 #error No implementation for this PROP_MODEL_NAME.
689 // If the neighbor is outside the latcie in X direction and we have a
690 // periodic boundary then we need to wrap around.
691 if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) ||
692 ((ny < 0 || ny >= lY) && ld->PeriodicY) ||
693 ((nz < 0 || nz >= lZ) && ld->PeriodicZ)
725 if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
726 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
729 neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];
731 AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);
733 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
736 else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
737 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
739 else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
740 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
743 neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];
745 Assert(neighborIndex != ~0);
747 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
750 Assert(dstIndex >= 0);
751 Assert(dstIndex < nCells * N_D3Q19);
753 adjList[index * N_D3Q19_IDX + d] = dstIndex;
760 nThreads = omp_get_max_threads();
763 SetupConsecNodes(ld, KDLR(kd), nThreads);
766 // Fill remaining KernelData structures
767 kd->GetNode = GetNode;
768 kd->SetNode = SetNode;
770 kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
771 kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);
773 kd->Kernel = NULL; // FNAME(KernelPullSplitNt2S);
776 kd->PdfsActive = kd->Pdfs[0];
781 void FNAME(D3Q19ListPullSplitNt1SInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
783 FNAME(Init)(ld, kernelData, params);
784 (*kernelData)->Kernel = FNAME(KernelPullSplitNt1S);
786 double loopBalance = 2.0 * 19 * sizeof(PdfT) + (18 * 4.0);
787 printf("# loop balance: %.2f B/FLUP\n", loopBalance);
790 void FNAME(D3Q19ListPullSplitNt2SInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
792 FNAME(Init)(ld, kernelData, params);
793 (*kernelData)->Kernel = FNAME(KernelPullSplitNt2S);
795 double loopBalance = 2.0 * 19 * sizeof(PdfT) + (18 * 4.0);
796 printf("# loop balance: %.2f B/FLUP\n", loopBalance);
800 void FNAME(D3Q19ListPullSplitNtDeinit)(LatticeDesc * ld, KernelData ** kernelData)
802 KernelDataListRia ** kdlr = (KernelDataListRia **)kernelData;
804 MemFree((void **)&((*kdlr)->ConsecNodes));
806 if ((*kdlr)->ConsecThreadIndices != NULL) {
807 MemFree((void **)&((*kdlr)->ConsecThreadIndices));
810 KernelDataList ** kdl = (KernelDataList **)kernelData;
812 MemFree((void **)&((*kdl)->AdjList));
813 MemFree((void **)&((*kdl)->Coords));
814 MemFree((void **)&((*kdl)->Grid));
816 MemFree((void **)&((*kernelData)->Pdfs[0]));
817 MemFree((void **)&((*kernelData)->Pdfs[1]));
819 MemFree((void **)kernelData);