src/BenchKernelD3Q19ListAaPvCommon.c

   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 "BenchKernelD3Q19ListAaPvCommon.h"
  28
  29 #include "Memory.h"
  30 #include "Vtk.h"
  31 #include "Padding.h"
  32
  33 #include <math.h>
  34
  35 #ifdef _OPENMP
  36         #include <omp.h>
  37 #endif
  38
  39 // Forward definition.
  40 void FNAME(D3Q19ListAaPvKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
  41
  42
  43
  44
  45 // -----------------------------------------------------------------------
  46 // Functions which are used as callback by the kernel to read or write
  47 // PDFs and nodes.
  48
  49 static void FNAME(BCGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
  50 {
  51         Assert(kd != NULL);
  52         Assert(kd->PdfsActive != NULL);
  53         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
  54         Assert(pdf != NULL);
  55
  56         Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
  57         Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
  58         Assert(dir >= 0); Assert(dir < N_D3Q19);
  59
  60         KernelDataList * kdl = (KernelDataList *)kd;
  61
  62         if (kdl->Iteration % 2 == 0) {
  63                 // Pdfs are stored inverse, local PDFs are located in remote nodes
  64
  65                 uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];
  66
  67                 if (dir != D3Q19_C) {
  68                         uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
  69
  70                         *pdf = kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + D3Q19_INV[dir]]];
  71                 }
  72                 else {
  73                         *pdf = kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)];
  74                 }
  75
  76         }
  77         else {
  78                 *pdf = kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)];
  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         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         KernelDataList * kdl = (KernelDataList *)kd;
 103
 104         if (kdl->Iteration % 2 == 0) {
 105                 // Pdfs are stored inverse, local PDFs are located in remote nodes
 106
 107                 uint32_t nodeIndex = KDL(kd)->Grid[L_INDEX_4(kd->Dims, x, y, z)];
 108
 109                 if (dir != D3Q19_C) {
 110                         uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
 111
 112                         kd->PdfsActive[KDL(kd)->AdjList[adjListIndex + D3Q19_INV[dir]]] = pdf;
 113                 }
 114                 else {
 115                         kd->PdfsActive[P_INDEX_3(KDL(kd)->nCells, nodeIndex, dir)] = pdf;
 116                 }
 117
 118         }
 119         else {
 120                 kd->PdfsActive[P_INDEX_5(KDL(kd), x, y, z, dir)] = pdf;
 121         }
 122
 123         return;
 124 }
 125
 126
 127 static void GetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
 128 {
 129         Assert(kd != NULL);
 130         Assert(kd->PdfsActive != NULL);
 131         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
 132         Assert(pdfs != NULL);
 133         Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
 134         Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
 135
 136         KernelDataList * kdl = (KernelDataList *)kd;
 137
 138         if(kdl->Iteration % 2 == 0){
 139
 140                 uint32_t nodeIndex = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
 141                 uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
 142
 143                 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
 144                 pdfs[D3Q19_C] = kd->PdfsActive[P_INDEX_3(kdl->nCells, nodeIndex, D3Q19_C)];
 145
 146                 #define X(name, idx, idxinv, _x, _y, _z)        pdfs[idx] = kd->PdfsActive[kdl->AdjList[adjListIndex + idxinv]];
 147                 D3Q19_LIST_WO_C
 148                 #undef X
 149
 150         } else {
 151
 152                 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
 153                 #define X(name, idx, idxinv, _x, _y, _z)        pdfs[idx] = kd->PdfsActive[I(x, y, z, idx)];
 154                 D3Q19_LIST
 155                 #undef X
 156                 #undef I
 157
 158         }
 159
 160 #if 0
 161         for (int d = 0; d < 19; ++d) {
 162                 if(isnan(pdfs[d]) || isinf(pdfs[d])) {
 163                         printf("%d %d %d %d nan! get node\n", x, y, z, d);
 164                                                 for (int d2 = 0; d2 < 19; ++d2) {
 165                                                         printf("%d: %e\n", d2, pdfs[d2]);
 166                                                 }
 167                         exit(1);
 168                 }
 169         }
 170 #endif
 171
 172         return;
 173 }
 174
 175
 176 static void SetNode(KernelData * kd, int x, int y, int z, PdfT * pdfs)
 177 {
 178         Assert(kd != NULL);
 179         Assert(kd->PdfsActive != NULL);
 180         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
 181         Assert(pdfs != NULL);
 182
 183         Assert(x >= 0); Assert(y >= 0); Assert(z >= 0);
 184         Assert(x < kd->Dims[0]); Assert(y < kd->Dims[1]); Assert(z < kd->Dims[2]);
 185
 186 #if 0
 187         for (int d = 0; d < 19; ++d) {
 188                 if(isnan(pdfs[d])) {
 189                         printf("%d %d %d %d nan! get node\n", x, y, z, d);
 190                                                 for (int d2 = 0; d2 < 19; ++d2) {
 191                                                         printf("%d: %e\n", d2, pdfs[d2]);
 192                                                 }
 193                         exit(1);
 194                 }
 195         }
 196 #endif
 197
 198         KernelDataList * kdl = (KernelDataList *)kd;
 199
 200         if(kdl->Iteration % 2 == 0){
 201
 202                 uint32_t nodeIndex = kdl->Grid[L_INDEX_4(kdl->kd.Dims, x, y, z)];
 203                 uint32_t adjListIndex = nodeIndex * N_D3Q19_IDX;
 204
 205                 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
 206                 kd->PdfsActive[P_INDEX_3(kdl->nCells, nodeIndex, D3Q19_C)] = pdfs[D3Q19_C];
 207
 208                 #define X(name, idx, idxinv, _x, _y, _z)        kd->PdfsActive[kdl->AdjList[adjListIndex + idxinv]] = pdfs[idx];
 209                 D3Q19_LIST_WO_C
 210                 #undef X
 211
 212         } else {
 213
 214                 #define I(x, y, z, dir) P_INDEX_5(KDL(kd), (x), (y), (z), (dir))
 215                 #define X(name, idx, idxinv, _x, _y, _z)        kd->PdfsActive[I(x, y, z, idx)] = pdfs[idx];
 216                 D3Q19_LIST
 217                 #undef X
 218                 #undef I
 219
 220         }
 221
 222         return;
 223 }
 224
 225 static void ParameterUsage()
 226 {
 227         printf("Kernel parameters:\n");
 228         printf("  [-blk <n>] [-blk-[xyz] <n>]\n");
 229 #ifdef DATA_LAYOUT_SOA
 230         printf("  [-pad auto|modulus_1+offset_1(,modulus_n+offset_n)*]\n");
 231 #endif
 232
 233         return;
 234 }
 235
 236 static void ParseParameters(Parameters * params, int * blk, PadInfo ** padInfo)
 237 {
 238         Assert(blk != NULL);
 239
 240         blk[0] = 0; blk[1] = 0; blk[2] = 0;
 241         *padInfo = NULL;
 242
 243         #define ARG_IS(param)                   (!strcmp(params->KernelArgs[i], param))
 244         #define NEXT_ARG_PRESENT() \
 245                 do { \
 246                         if (i + 1 >= params->nKernelArgs) { \
 247                                 printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
 248                                 exit(1); \
 249                         } \
 250                 } while (0)
 251
 252
 253         for (int i = 0; i < params->nKernelArgs; ++i) {
 254                 if (ARG_IS("-blk") || ARG_IS("--blk")) {
 255                         NEXT_ARG_PRESENT();
 256
 257                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 258
 259                         if (tmp < 0) {
 260                                 printf("ERROR: blocking parameter must be >= 0.\n");
 261                                 exit(1);
 262                         }
 263
 264                         blk[0] = blk[1] = blk[2] = tmp;
 265                 }
 266                 else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
 267                         NEXT_ARG_PRESENT();
 268
 269                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 270
 271                         if (tmp < 0) {
 272                                 printf("ERROR: blocking parameter must be >= 0.\n");
 273                                 exit(1);
 274                         }
 275
 276                         blk[0] = tmp;
 277                 }
 278                 else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
 279                         NEXT_ARG_PRESENT();
 280
 281                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 282
 283                         if (tmp < 0) {
 284                                 printf("ERROR: blocking parameter must be >= 0.\n");
 285                                 exit(1);
 286                         }
 287
 288                         blk[1] = tmp;
 289                 }
 290                 else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
 291                         NEXT_ARG_PRESENT();
 292
 293                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 294
 295                         if (tmp < 0) {
 296                                 printf("ERROR: blocking parameter must be >= 0.\n");
 297                                 exit(1);
 298                         }
 299
 300                         blk[2] = tmp;
 301                 }
 302 #ifdef DATA_LAYOUT_SOA
 303                 else if (ARG_IS("-pad") || ARG_IS("--pad")) {
 304                         NEXT_ARG_PRESENT();
 305
 306                         *padInfo = PadInfoFromStr(params->KernelArgs[++i]);
 307                 }
 308 #endif
 309                 else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
 310                         ParameterUsage();
 311                         exit(1);
 312                 }
 313                 else {
 314                         printf("ERROR: unknown kernel parameter.\n");
 315                         ParameterUsage();
 316                         exit(1);
 317                 }
 318         }
 319
 320         #undef ARG_IS
 321         #undef NEXT_ARG_PRESENT
 322
 323         return;
 324 }
 325
 326 static void SetupConsecNodes(LatticeDesc * ld, KernelDataListRia * kdlr, int nThreads)
 327 {
 328         Assert(ld != NULL);
 329         Assert(kdlr != NULL);
 330         Assert(nThreads > 0);
 331
 332         uint32_t * adjList = kdlr->kdl.AdjList;
 333
 334         uint32_t nConsecNodes = 0;
 335         uint32_t consecIndex = 0;
 336
 337         int nFluid = kdlr->kdl.nFluid;
 338
 339         uint32_t * consecThreadIndices = (uint32_t *)malloc(sizeof(uint32_t) * (nThreads + 1));
 340         int * fluidNodeThreadIndices = (int *)malloc(sizeof(int) * (nThreads + 1));
 341
 342         int nNodesPerThread = nFluid / nThreads;
 343
 344         for (int i = 0; i < nThreads; ++i) {
 345                 consecThreadIndices[i]     = i * nNodesPerThread + MinI(i, nFluid % nThreads);
 346                 fluidNodeThreadIndices[i] = consecThreadIndices[i];
 347         }
 348         consecThreadIndices[nThreads]     = -1;
 349         fluidNodeThreadIndices[nThreads] = nFluid;
 350
 351         int indexThread = 1;
 352
 353         // We execute following code two times.
 354         // - The first time to get the count of how many entries we need for the
 355         //   consecNodes array.
 356         // - The second time to fill the array.
 357
 358         // Loop over adjacency list of all nodes.
 359     // Compare if adjacent nodes share the same access pattern.
 360         for (int index = 1; index < nFluid; ++index) {
 361
 362                 int different = 0;
 363
 364                 // Loop over all directions except the center one.
 365                 for(int d = 0; d < N_D3Q19 - 1; ++d) {
 366                         Assert(d != D3Q19_C);
 367
 368                         if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1) {
 369                                 // Different access pattern.
 370                                 different = 1;
 371                                 break;
 372                         }
 373                 }
 374
 375                 if (consecThreadIndices[indexThread] == index) {
 376                         // We are at a thread boundary. Starting from this index the fluids
 377                         // belong to another thread. Force a break, if nodes are consecutive.
 378                         ++indexThread;
 379                         different = 1;
 380                 }
 381
 382                 if (different) {
 383                         ++consecIndex;
 384                 }
 385         }
 386
 387         if (nFluid > 0) {
 388                 nConsecNodes = consecIndex + 1;
 389         }
 390
 391         uint32_t * consecNodes;
 392         MemAlloc((void **)&consecNodes, sizeof(uint32_t) * nConsecNodes);
 393
 394         consecIndex = 0;
 395
 396         if (nFluid > 0) {
 397                 consecNodes[consecIndex] = 1;
 398         }
 399
 400         indexThread = 1;
 401         consecThreadIndices[0] = 0;
 402
 403         // Loop over adjacency list of all nodes.
 404     // Compare if adjacent nodes share the same access pattern.
 405         for (int index = 1; index < nFluid; ++index) {
 406
 407                 int different = 0;
 408
 409                 // Loop over all directions except the center one.
 410                 for(int d = 0; d < N_D3Q19 - 1; ++d) {
 411                         Assert(d != D3Q19_C);
 412
 413                         if (adjList[index * N_D3Q19_IDX + d] != adjList[(index - 1) * N_D3Q19_IDX + d] + 1) {
 414                                 // Different access pattern.
 415                                 different = 1;
 416                                 break;
 417                         }
 418                 }
 419
 420                 if (consecThreadIndices[indexThread] == index) {
 421                         // We are at a thread boundary. Starting from this index the fluids
 422                         // belong to another thread. Force a break, if nodes are consecutive.
 423                         consecThreadIndices[indexThread] = consecIndex + 1;
 424                         ++indexThread;
 425                         different = 1;
 426                 }
 427
 428                 if (different) {
 429                         ++consecIndex;
 430                         Assert(consecIndex < nConsecNodes);
 431                         consecNodes[consecIndex] = 1;
 432                 }
 433                 else {
 434                         Assert(consecIndex < nConsecNodes);
 435                         consecNodes[consecIndex] += 1;
 436                 }
 437         }
 438
 439
 440         kdlr->ConsecNodes = consecNodes;
 441         kdlr->nConsecNodes = nConsecNodes;
 442
 443         kdlr->ConsecThreadIndices  = consecThreadIndices;
 444         kdlr->nConsecThreadIndices = nThreads;
 445
 446         kdlr->FluidNodeThreadIndices = fluidNodeThreadIndices;
 447         kdlr->nFluidNodeThreadIndices = nThreads;
 448
 449         printf("# total fluid nodes: %d   consecutive blocks: %d\n", nFluid, nConsecNodes);
 450
 451         uint32_t vwidth[]       = {2, 4, 8, 16, 32};
 452         uint32_t vectorizable[] = {0, 0, 0,  0,  0};
 453
 454         for (int i = 0; i < nConsecNodes; ++i) {
 455                 for (int k = 0; k < N_ELEMS(vwidth); ++k) {
 456                         vectorizable[k] += consecNodes[i] / vwidth[k];
 457                 }
 458         }
 459
 460         printf("# vectorizable fraction of fluid node updates:\n");
 461         for (int i = 0; i < N_ELEMS(vwidth); ++i) {
 462
 463                 printf("#  vector width: %2d     %6.2f %% (%u/%u fluid nodes)\n",
 464                         vwidth[i], (double)vectorizable[i] * vwidth[i] / nFluid * 100.0,
 465                         vectorizable[i] * vwidth[i], nFluid);
 466         }
 467
 468         return;
 469 }
 470
 471 void FNAME(D3Q19ListAaPvInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
 472 {
 473         KernelData * kd;
 474         KernelDataList * kdl;
 475         KernelDataListRia * kdlr;
 476         MemAlloc((void **)&kdlr, sizeof(KernelDataListRia));
 477
 478         kd = (KernelData *)kdlr;
 479         kdl = KDL(kdlr);
 480
 481         *kernelData = kd;
 482
 483 #ifdef DEBUG
 484         kd->Pdfs[0] = NULL;
 485         kd->Pdfs[1] = NULL;
 486         kd->PdfsActive = NULL;
 487         kd->DstPdfs = NULL;
 488         kd->SrcPdfs = NULL;
 489         kd->Dims[0] = -1;
 490         kd->Dims[1] = -1;
 491         kd->Dims[2] = -1;
 492         kd->GlobalDims[0] = -1;
 493         kd->GlobalDims[1] = -1;
 494         kd->GlobalDims[2] = -1;
 495         kd->Offsets[0] = -1;
 496         kd->Offsets[1] = -1;
 497         kd->Offsets[2] = -1;
 498
 499         kd->ObstIndices = NULL;
 500         kd->nObstIndices = -1;
 501         kd->BounceBackPdfsSrc = NULL;
 502         kd->BounceBackPdfsDst = NULL;
 503         kd->nBounceBackPdfs = -1;
 504
 505         kdl->AdjList = NULL;
 506         kdl->Coords = NULL;
 507         kdl->Grid = NULL;
 508         kdl->nCells = -1;
 509         kdl->nFluid = -1;
 510
 511         kdlr->ConsecNodes = NULL;
 512         kdlr->nConsecNodes = 0;
 513         kdlr->ConsecThreadIndices = NULL;
 514         kdlr->nConsecThreadIndices = 0;
 515 #endif
 516
 517         int blk[3] = { 0 };
 518         PadInfo * padInfo = NULL;
 519
 520         ParseParameters(params, blk, &padInfo);
 521
 522         // Ajust the dimensions according to padding, if used.
 523         kd->Dims[0] = kd->GlobalDims[0] = ld->Dims[0];
 524         kd->Dims[1] = kd->GlobalDims[1] = ld->Dims[1];
 525         kd->Dims[2] = kd->GlobalDims[2] = ld->Dims[2];
 526
 527         int * lDims = ld->Dims;
 528
 529         int lX = lDims[0];
 530         int lY = lDims[1];
 531         int lZ = lDims[2];
 532
 533         int nTotalCells = lX * lY * lZ;
 534         int nCells = ld->nFluid;
 535         int nFluid = ld->nFluid;
 536
 537 #ifdef DATA_LAYOUT_SOA
 538         {
 539                 nCells = PadCellsAndReport(nCells, sizeof(PdfT), &padInfo);
 540                 PadInfoFree(padInfo); padInfo = NULL;
 541         }
 542 #endif
 543
 544         kdl->nCells = nCells;
 545         kdl->nFluid = nFluid;
 546
 547         PdfT * pdfs[2];
 548
 549
 550         if (blk[0] == 0) blk[0] = lX;
 551         if (blk[1] == 0) blk[1] = lY;
 552         if (blk[2] == 0) blk[2] = lZ;
 553
 554         printf("# blocking               x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
 555
 556         double latMiB      = nCells * sizeof(PdfT) * N_D3Q19 / 1024.0 / 1024.0;
 557         double latFluidMib = nFluid * sizeof(PdfT) * N_D3Q19 / 1024.0 / 1024.0;
 558         double latPadMib   = (nCells - nFluid) * sizeof(PdfT) * N_D3Q19 / 1024.0 / 1024.0;
 559
 560         printf("# lattice size:          %e MiB\n", latMiB);
 561         printf("# fluid lattice size:    %e MiB\n", latFluidMib);
 562         printf("# lattice padding:       %e MiB\n", latPadMib);
 563
 564 #define PAGE_4K         4096
 565
 566         printf("# aligning lattices to:  %d b\n", PAGE_4K);
 567
 568         MemAllocAligned((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19, PAGE_4K);
 569
 570         kd->Pdfs[0] = pdfs[0];
 571
 572         // Initialize PDFs with some (arbitrary) data for correct NUMA placement.
 573         // Here we touch only the fluid nodes as this loop is OpenMP parallel and
 574         // we want the same scheduling as in the kernel.
 575         #ifdef _OPENMP
 576                 #pragma omp parallel for
 577         #endif
 578         for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
 579                 pdfs[0][P_INDEX_3(nCells, i, d)] = 1.0;
 580         } }
 581
 582         // Initialize all PDFs to some standard value.
 583         for (int i = 0; i < nFluid; ++i) { for(int d = 0; d < N_D3Q19; ++d) {
 584                 pdfs[0][P_INDEX_3(nCells, i, d)] = 0.0;
 585         } }
 586
 587         // ----------------------------------------------------------------------
 588         // create grid which will hold the index numbers of the fluid nodes
 589
 590         uint32_t * grid;
 591
 592         if (MemAlloc((void **)&grid, nTotalCells * sizeof(uint32_t))) {
 593                 printf("ERROR: allocating grid for numbering failed: %lu bytes.\n", nTotalCells * sizeof(uint32_t));
 594                 exit(1);
 595         }
 596         kdl->Grid = grid;
 597
 598         int latticeIndex;
 599
 600 #ifdef DEBUG
 601         for(int z = 0; z < lZ; ++z) {
 602                 for(int y = 0; y < lY; ++y) {
 603                         for(int x = 0; x < lX; ++x) {
 604
 605                                 latticeIndex = L_INDEX_4(ld->Dims, x, y, z);
 606
 607                                 grid[latticeIndex] = ~0;
 608                         }
 609                 }
 610         }
 611 #endif
 612
 613         // ----------------------------------------------------------------------
 614         // generate numbering over grid
 615
 616         uint32_t * coords;
 617
 618         if (MemAlloc((void **)&coords, nFluid * sizeof(uint32_t) * 3)) {
 619                 printf("ERROR: allocating coords array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * 3);
 620                 exit(1);
 621         }
 622
 623         kdl->Coords = coords;
 624
 625         // Index for the PDF nodes can start at 0 as we distinguish solid and fluid nodes
 626         // through the ld->Lattice array.
 627         int counter = 0;
 628
 629         // Blocking is implemented via setup of the adjacency list. The kernel later will
 630         // walk through the lattice blocked automatically.
 631         for (int bX = 0; bX < lX; bX += blk[0]) {
 632         for (int bY = 0; bY < lY; bY += blk[1]) {
 633         for (int bZ = 0; bZ < lZ; bZ += blk[2]) {
 634
 635                 int eX = MIN(bX + blk[0], lX);
 636                 int eY = MIN(bY + blk[1], lY);
 637                 int eZ = MIN(bZ + blk[2], lZ);
 638
 639
 640                 for (int x = bX; x < eX; ++x) {
 641                 for (int y = bY; y < eY; ++y) {
 642                 for (int z = bZ; z < eZ; ++z) {
 643
 644                         latticeIndex = L_INDEX_4(lDims, x, y, z);
 645
 646                         if (ld->Lattice[latticeIndex] != LAT_CELL_OBSTACLE) {
 647                                 grid[latticeIndex] = counter;
 648
 649                                 coords[C_INDEX_X(counter)] = x;
 650                                 coords[C_INDEX_Y(counter)] = y;
 651                                 coords[C_INDEX_Z(counter)] = z;
 652
 653                                 ++counter;
 654                         }
 655                 } } }
 656         } } }
 657
 658         Verify(counter == nFluid);
 659
 660         uint32_t * adjList;
 661
 662         double indexMib = nFluid * sizeof(uint32_t) * N_D3Q19_IDX / 1024.0 / 1024.0;
 663
 664         printf("# index size:            %e MiB\n", indexMib);
 665
 666         // AdjList only requires 18 instead of 19 entries per node, as
 667         // the center PDF needs no addressing.
 668         if (MemAlloc((void **)&adjList, nFluid * sizeof(uint32_t) * N_D3Q19_IDX)) {
 669                 printf("ERROR: allocating adjList array failed: %lu bytes.\n", nFluid * sizeof(uint32_t) * N_D3Q19_IDX);
 670                 exit(1);
 671         }
 672
 673         kdl->AdjList = adjList;
 674
 675         int x, y, z;
 676
 677         uint32_t neighborIndex;
 678         uint32_t dstIndex;
 679
 680         int nx, ny, nz, px, py, pz;
 681
 682         // Loop over all fluid nodes and compute the indices to the neighboring
 683         // PDFs for configured data layout (AoS/SoA).
 684         #ifdef _OPENMP
 685                 #pragma omp parallel for
 686         #endif
 687         for (int index = 0; index < nFluid; ++index) {
 688                 for (int d = 0; d < N_D3Q19_IDX; ++d) {
 689                         adjList[index * N_D3Q19_IDX + d] = -1;
 690                 }
 691         }
 692
 693         // #ifdef _OPENMP  --> add line continuation
 694         //      #pragma omp parallel for default(none)
 695         //              shared(nFluid, nCells, coords, D3Q19_INV, D3Q19_X, D3Q19_Y, D3Q19_Z,
 696         //                              stderr,
 697         //                              lDims, grid, ld, lX, lY, lZ, adjList)
 698         //              private(x, y, z, nx, ny, nz, neighborIndex, dstIndex)
 699         // #endif
 700         for (int index = 0; index < nFluid; ++index) {
 701                 x = coords[C_INDEX_X(index)];
 702                 y = coords[C_INDEX_Y(index)];
 703                 z = coords[C_INDEX_Z(index)];
 704
 705                 Assert(x >= 0 && x < lX);
 706                 Assert(y >= 0 && y < lY);
 707                 Assert(z >= 0 && z < lZ);
 708
 709                 Assert(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE);
 710
 711                 // Loop over all directions except the center one.
 712                 for(int d = 0; d < N_D3Q19 - 1; ++d) {
 713                         Assert(d != D3Q19_C);
 714
 715 #ifdef PROP_MODEL_PUSH
 716                         nx = x + D3Q19_X[d];
 717                         ny = y + D3Q19_Y[d];
 718                         nz = z + D3Q19_Z[d];
 719
 720 #elif PROP_MODEL_PULL
 721                         nx = x - D3Q19_X[d];
 722                         ny = y - D3Q19_Y[d];
 723                         nz = z - D3Q19_Z[d];
 724 #else
 725                         #error No implementation for this PROP_MODEL_NAME.
 726 #endif
 727                         // If the neighbor is outside the latcie in X direction and we have a
 728                         // periodic boundary then we need to wrap around.
 729                         if (    ((nx < 0 || nx >= lX) && ld->PeriodicX) ||
 730                                         ((ny < 0 || ny >= lY) && ld->PeriodicY) ||
 731                                         ((nz < 0 || nz >= lZ) && ld->PeriodicZ)
 732                                                                                                                                 ){
 733                                 // x periodic
 734
 735                                 if (nx < 0) {
 736                                         px = lX - 1;
 737                                 }
 738                                 else if (nx >= lX) {
 739                                         px = 0;
 740                                 } else {
 741                                         px = nx;
 742                                 }
 743                                 // y periodic
 744                                 if (ny < 0) {
 745                                         py = lY - 1;
 746                                 }
 747                                 else if (ny >= lY) {
 748                                         py = 0;
 749                                 } else {
 750                                         py = ny;
 751                                 }
 752
 753                                 // z periodic
 754                                 if (nz < 0) {
 755                                         pz = lZ - 1;
 756                                 }
 757                                 else if (nz >= lZ) {
 758                                         pz = 0;
 759                                 } else {
 760                                         pz = nz;
 761                                 }
 762
 763                                 if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
 764                                         dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
 765                                 }
 766                                 else {
 767                                         neighborIndex = grid[L_INDEX_4(lDims, px, py, pz)];
 768
 769                                         AssertMsg(neighborIndex != ~0, "Neighbor has no Index. (%d %d %d) direction %s (%d)\n", px, py, pz, D3Q19_NAMES[d], d);
 770
 771                                         dstIndex = P_INDEX_3(nCells, neighborIndex, d);
 772                                 }
 773                         }
 774                         else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
 775                                 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
 776                         }
 777                         else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
 778                                 dstIndex = P_INDEX_3(nCells, index, D3Q19_INV[d]);
 779                         }
 780                         else {
 781                                 neighborIndex = grid[L_INDEX_4(lDims, nx, ny, nz)];
 782
 783                                 Assert(neighborIndex != ~0);
 784
 785                                 dstIndex = P_INDEX_3(nCells, neighborIndex, d);
 786                         }
 787
 788                         Assert(dstIndex >= 0);
 789                         Assert(dstIndex < nCells * N_D3Q19);
 790
 791                         adjList[index * N_D3Q19_IDX + d] = dstIndex;
 792                 }
 793         }
 794
 795         int nThreads = 1;
 796
 797 #ifdef _OPENMP
 798         nThreads = omp_get_max_threads();
 799 #endif
 800
 801         SetupConsecNodes(ld, KDLR(kd), nThreads);
 802
 803         double loopBalanceEven = 2.0 * 19 * sizeof(PdfT);
 804         double loopBalanceOdd  = 2.0 * 19 * sizeof(PdfT) + (double)kdlr->nConsecNodes / nFluid * (18 * 4.0 + 4.0);
 805         double loopBalance     = (loopBalanceEven + loopBalanceOdd) / 2.0;
 806
 807         printf("# loop balance:          %.2f B/FLUP  even: %.2f B/FLUP  odd %.2f B/FLUP\n",
 808                         loopBalance, loopBalanceEven, loopBalanceOdd);
 809
 810         // Fill remaining KernelData structures
 811         kd->GetNode = GetNode;
 812         kd->SetNode = SetNode;
 813
 814         kd->BoundaryConditionsGetPdf = FNAME(BCGetPdf);
 815         kd->BoundaryConditionsSetPdf = FNAME(BCSetPdf);
 816
 817         kd->Kernel = FNAME(D3Q19ListAaPvKernel);
 818
 819         kd->DstPdfs = NULL;
 820         kd->PdfsActive = kd->Pdfs[0];
 821
 822         return;
 823 }
 824
 825 void FNAME(D3Q19ListAaPvDeinit)(LatticeDesc * ld, KernelData ** kernelData)
 826 {
 827         KernelDataListRia ** kdlr = (KernelDataListRia **)kernelData;
 828
 829         MemFree((void **)&((*kdlr)->ConsecNodes));
 830
 831         if ((*kdlr)->ConsecThreadIndices != NULL) {
 832                 MemFree((void **)&((*kdlr)->ConsecThreadIndices));
 833         }
 834
 835         if ((*kdlr)->FluidNodeThreadIndices != NULL) {
 836                 MemFree((void **)&((*kdlr)->FluidNodeThreadIndices));
 837         }
 838
 839         KernelDataList ** kdl = (KernelDataList **)kernelData;
 840
 841         MemFree((void **)&((*kdl)->AdjList));
 842         MemFree((void **)&((*kdl)->Coords));
 843         MemFree((void **)&((*kdl)->Grid));
 844
 845         MemFree((void **)&((*kernelData)->Pdfs[0]));
 846
 847         MemFree((void **)kernelData);
 848         return;
 849 }
 850