src/BenchKernelD3Q19Common.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 "BenchKernelD3Q19Common.h"
  28
  29 #include "Memory.h"
  30 #include "Vtk.h"
  31
  32 #include <inttypes.h>
  33 #include <math.h>
  34
  35
  36 // Forward definition.
  37 void FNAME(D3Q19Kernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
  38
  39 void FNAME(D3Q19BlkKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd);
  40
  41
  42
  43 static void FNAME(BcGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf)
  44 {
  45         Assert(kd != NULL);
  46         Assert(kd->PdfsActive != NULL);
  47         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
  48         Assert(pdf != NULL);
  49
  50         Assert(x >= 0);
  51         Assert(y >= 0);
  52         Assert(z >= 0);
  53         Assert(x < kd->Dims[0]);
  54         Assert(y < kd->Dims[1]);
  55         Assert(z < kd->Dims[2]);
  56         Assert(dir >= 0);
  57         Assert(dir < N_D3Q19);
  58
  59         int oX = kd->Offsets[0];
  60         int oY = kd->Offsets[1];
  61         int oZ = kd->Offsets[2];
  62
  63 #ifdef PROP_MODEL_PUSH
  64         int nx = x;
  65         int ny = y;
  66         int nz = z;
  67 #elif PROP_MODEL_PULL
  68         int nx = x - D3Q19_X[dir];
  69         int ny = y - D3Q19_Y[dir];
  70         int nz = z - D3Q19_Z[dir];
  71 #endif
  72
  73         #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir))
  74         *pdf = kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, dir)];
  75         #undef I
  76
  77         return;
  78 }
  79
  80 static void FNAME(BcSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf)
  81 {
  82         Assert(kd != NULL);
  83         Assert(kd->PdfsActive != NULL);
  84         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
  85         Assert(x >= 0);
  86         Assert(y >= 0);
  87         Assert(z >= 0);
  88         Assert(x < kd->Dims[0]);
  89         Assert(y < kd->Dims[1]);
  90         Assert(z < kd->Dims[2]);
  91         Assert(dir >= 0);
  92         Assert(dir < N_D3Q19);
  93
  94         int oX = kd->Offsets[0];
  95         int oY = kd->Offsets[1];
  96         int oZ = kd->Offsets[2];
  97
  98 #ifdef PROP_MODEL_PUSH
  99         int nx = x;
 100         int ny = y;
 101         int nz = z;
 102 #elif PROP_MODEL_PULL
 103         int nx = x - D3Q19_X[dir];
 104         int ny = y - D3Q19_Y[dir];
 105         int nz = z - D3Q19_Z[dir];
 106 #endif
 107
 108         #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir))
 109         kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, dir)] = pdf;
 110         #undef I
 111
 112
 113         return;
 114 }
 115
 116
 117 static void FNAME(GetNode)(KernelData * kd, int x, int y, int z, PdfT * pdfs)
 118 {
 119         Assert(kd != NULL);
 120         Assert(kd->PdfsActive != NULL);
 121         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
 122         Assert(pdfs != NULL);
 123         Assert(x >= 0);
 124         Assert(y >= 0);
 125         Assert(z >= 0);
 126         Assert(x < kd->Dims[0]);
 127         Assert(y < kd->Dims[1]);
 128         Assert(z < kd->Dims[2]);
 129
 130         int oX = kd->Offsets[0];
 131         int oY = kd->Offsets[1];
 132         int oZ = kd->Offsets[2];
 133
 134
 135         #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir))
 136 #ifdef PROP_MODEL_PUSH
 137         #define X(name, idx, idxinv, _x, _y, _z)        pdfs[idx] = kd->PdfsActive[I(x + oX, y + oY, z + oZ, idx)];
 138 #elif PROP_MODEL_PULL
 139         #define X(name, idx, idxinv, _x, _y, _z)        pdfs[idx] = kd->PdfsActive[I(x + oX - (_x), y + oY - (_y), z + oZ - (_z), idx)];
 140 #endif
 141         D3Q19_LIST
 142         #undef X
 143         #undef I
 144
 145 #if 0           // DETECT NANs
 146
 147         for (int d = 0; d < 19; ++d) {
 148                 if (isnan(pdfs[d])) {
 149                         printf("%d %d %d %d nan! get node\n", x, y, z, d);
 150
 151                         for (int d2 = 0; d2 < 19; ++d2) {
 152                                 printf("%d: %e\n", d2, pdfs[d2]);
 153                         }
 154
 155                         exit(1);
 156                 }
 157         }
 158
 159 #endif
 160
 161         return;
 162 }
 163
 164
 165 static void FNAME(SetNode)(KernelData * kd, int x, int y, int z, PdfT * pdfs)
 166 {
 167         Assert(kd != NULL);
 168         Assert(kd->PdfsActive != NULL);
 169         Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]);
 170         Assert(pdfs != NULL);
 171
 172         Assert(x >= 0);
 173         Assert(y >= 0);
 174         Assert(z >= 0);
 175         Assert(x < kd->Dims[0]);
 176         Assert(y < kd->Dims[1]);
 177         Assert(z < kd->Dims[2]);
 178
 179         int oX = kd->Offsets[0];
 180         int oY = kd->Offsets[1];
 181         int oZ = kd->Offsets[2];
 182
 183         #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir))
 184 #ifdef PROP_MODEL_PUSH
 185         #define X(name, idx, idxinv, _x, _y, _z)        kd->PdfsActive[I(x + oX, y + oY, z + oZ, idx)] = pdfs[idx];
 186 #elif PROP_MODEL_PULL
 187         #define X(name, idx, idxinv, _x, _y, _z)        kd->PdfsActive[I(x + oX - (_x), y + oY - (_y), z + oZ - (_z), idx)] = pdfs[idx];
 188 #endif
 189         D3Q19_LIST
 190         #undef X
 191         #undef I
 192
 193         return;
 194 }
 195
 196
 197 static void ParameterUsage()
 198 {
 199         printf("Kernel parameters:\n");
 200         printf("  [-blk <n>] [-blk-[xyz] <n>]\n");
 201
 202         return;
 203 }
 204
 205 static void ParseParameters(Parameters * params, int * blk)
 206 {
 207         Assert(blk != NULL);
 208
 209         blk[0] = 0; blk[1] = 0; blk[2] = 0;
 210
 211         #define ARG_IS(param)                   (!strcmp(params->KernelArgs[i], param))
 212         #define NEXT_ARG_PRESENT() \
 213                 do { \
 214                         if (i + 1 >= params->nKernelArgs) { \
 215                                 printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \
 216                                 exit(1); \
 217                         } \
 218                 } while (0)
 219
 220
 221         for (int i = 0; i < params->nKernelArgs; ++i) {
 222                 if (ARG_IS("-blk") || ARG_IS("--blk")) {
 223                         NEXT_ARG_PRESENT();
 224
 225                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 226
 227                         if (tmp <= 0) {
 228                                 printf("ERROR: blocking parameter must be > 0.\n");
 229                                 exit(1);
 230                         }
 231
 232                         blk[0] = blk[1] = blk[2] = tmp;
 233                 }
 234                 else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) {
 235                         NEXT_ARG_PRESENT();
 236
 237                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 238
 239                         if (tmp <= 0) {
 240                                 printf("ERROR: blocking parameter must be > 0.\n");
 241                                 exit(1);
 242                         }
 243
 244                         blk[0] = tmp;
 245                 }
 246                 else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) {
 247                         NEXT_ARG_PRESENT();
 248
 249                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 250
 251                         if (tmp <= 0) {
 252                                 printf("ERROR: blocking parameter must be > 0.\n");
 253                                 exit(1);
 254                         }
 255
 256                         blk[1] = tmp;
 257                 }
 258                 else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) {
 259                         NEXT_ARG_PRESENT();
 260
 261                         int tmp = strtol(params->KernelArgs[++i], NULL, 0);
 262
 263                         if (tmp <= 0) {
 264                                 printf("ERROR: blocking parameter must be > 0.\n");
 265                                 exit(1);
 266                         }
 267
 268                         blk[2] = tmp;
 269                 }
 270                 else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) {
 271                         ParameterUsage();
 272                         exit(1);
 273                 }
 274                 else {
 275                         printf("ERROR: unknown kernel parameter.\n");
 276                         ParameterUsage();
 277                         exit(1);
 278                 }
 279         }
 280
 281         #undef ARG_IS
 282         #undef NEXT_ARG_PRESENT
 283
 284         return;
 285 }
 286
 287
 288 void FNAME(D3Q19BlkInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
 289 {
 290         KernelDataEx * kdex = NULL;
 291         MemAlloc((void **)&kdex, sizeof(KernelDataEx));
 292
 293         kdex->Blk[0] = 0; kdex->Blk[1] = 0; kdex->Blk[2] = 0;
 294
 295         KernelData * kd = &kdex->kd;
 296         *kernelData = kd;
 297
 298         kd->nObstIndices = ld->nObst;
 299
 300         // Ajust the dimensions according to padding, if used.
 301         kd->Dims[0] = ld->Dims[0];
 302         kd->Dims[1] = ld->Dims[1];
 303         kd->Dims[2] = ld->Dims[2];
 304
 305
 306         int * lDims = ld->Dims;
 307         int * gDims = kd->GlobalDims;
 308
 309         gDims[0] = lDims[0] + 2;
 310         gDims[1] = lDims[1] + 2;
 311         gDims[2] = lDims[2] + 2;
 312
 313         kd->Offsets[0] = 1;
 314         kd->Offsets[1] = 1;
 315         kd->Offsets[2] = 1;
 316
 317         int lX = lDims[0];
 318         int lY = lDims[1];
 319         int lZ = lDims[2];
 320
 321         int gX = gDims[0];
 322         int gY = gDims[1];
 323         int gZ = gDims[2];
 324
 325         int oX = kd->Offsets[0];
 326         int oY = kd->Offsets[1];
 327         int oZ = kd->Offsets[2];
 328
 329         int blk[3] = { 0 };
 330
 331         int nCells = gX * gY * gZ;
 332
 333         PdfT * pdfs[2];
 334
 335         ParseParameters(params, blk);
 336
 337         if (blk[0] == 0) blk[0] = gX;
 338         if (blk[1] == 0) blk[1] = gY;
 339         if (blk[2] == 0) blk[2] = gZ;
 340
 341         printf("# blocking x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]);
 342
 343
 344         kdex->Blk[0] = blk[0]; kdex->Blk[1] = blk[1]; kdex->Blk[2] = blk[2];
 345
 346
 347         printf("# allocating data for %d LB nodes with padding (%lu bytes = %f MiB for both lattices)\n",
 348                nCells, 2 * sizeof(PdfT) * nCells * N_D3Q19,
 349                2 * sizeof(PdfT) * nCells * N_D3Q19 / 1024.0 / 1024.0);
 350
 351         MemAlloc((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19);
 352         MemAlloc((void **)&pdfs[1], sizeof(PdfT) * nCells * N_D3Q19);
 353
 354         kd->Pdfs[0] = pdfs[0];
 355         kd->Pdfs[1] = pdfs[1];
 356
 357         // Initialize PDFs with some (arbitrary) data for correct NUMA placement.
 358         // This depends on the chosen data layout.
 359         // The structure of the loop should resemble the same "execution layout"
 360         // as in the kernel!
 361 #ifdef _OPENMP
 362         #pragma omp parallel for collapse(3)
 363 #endif
 364
 365         for (int bZ = 0; bZ < gZ; bZ += blk[2]) {
 366                 for (int bY = 0; bY < gY; bY += blk[1]) {
 367                         for (int bX = 0; bX < gX; bX += blk[0]) {
 368
 369                                 // Must do everything here, else it would break collapse.
 370                                 int eZ = MIN(bZ + blk[2], gZ);
 371                                 int eY = MIN(bY + blk[1], gY);
 372                                 int eX = MIN(bX + blk[0], gX);
 373
 374                                 for (int z = bZ; z < eZ; ++z) {
 375                                         for (int y = bY; y < eY; ++y) {
 376                                                 for (int x = bX; x < eX; ++x) {
 377
 378                                                         for (int d = 0; d < N_D3Q19; ++d) {
 379                                                                 pdfs[0][P_INDEX_5(gDims, x, y, z, d)] = 1.0;
 380                                                                 pdfs[1][P_INDEX_5(gDims, x, y, z, d)] = 1.0;
 381                                                         }
 382
 383                                                 }
 384                                         }
 385                                 }
 386                         }
 387                 }
 388         }
 389
 390         // Initialize all PDFs to some standard value.
 391         for (int z = 0; z < gZ; ++z) {
 392                 for (int y = 0; y < gY; ++y) {
 393                         for (int x = 0; x < gX; ++x) {
 394                                 for (int d = 0; d < N_D3Q19; ++d) {
 395                                         pdfs[0][P_INDEX_5(gDims, x, y, z, d)] = 0.0;
 396                                         pdfs[1][P_INDEX_5(gDims, x, y, z, d)] = 0.0;
 397                                 }
 398                         }
 399                 }
 400         }
 401
 402
 403         // Count how many *PDFs* need bounce back treatment.
 404
 405         uint64_t nPdfs = ((uint64_t)19) * gX * gY * gZ;
 406
 407         if (nPdfs > ((2LU << 31) - 1)) {
 408                 printf("ERROR: number of PDFs exceed 2^31.\n");
 409                 exit(1);
 410         }
 411
 412         // Compiler bug? Incorrect computation of nBounceBackPdfs when using icc 15.0.2.
 413         // Works when declaring nBounceBackPdfs as int64_t or using volatile.
 414         volatile int nBounceBackPdfs = 0;
 415         // int64_t nBounceBackPdfs = 0;
 416         int nx, ny, nz, px, py, pz;
 417
 418         // TODO: apply blocking?
 419
 420         for (int z = 0; z < lZ; ++z) {
 421                 for (int y = 0; y < lY; ++y) {
 422                         for (int x = 0; x < lX; ++x) {
 423
 424                                 if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
 425                                         for (int d = 0; d < N_D3Q19; ++d) {
 426 #ifdef PROP_MODEL_PUSH
 427                                                 nx = x + D3Q19_X[d];
 428                                                 ny = y + D3Q19_Y[d];
 429                                                 nz = z + D3Q19_Z[d];
 430 #elif PROP_MODEL_PULL
 431                                                 nx = x - D3Q19_X[d];
 432                                                 ny = y - D3Q19_Y[d];
 433                                                 nz = z - D3Q19_Z[d];
 434 #else
 435         #error PROP_MODEL_NAME unknown.
 436 #endif
 437                                                 // Check if neighbor is inside the lattice.
 438                                                 // if(nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
 439                                                 //      continue;
 440                                                 // }
 441                                                 if ((nx < 0 || nx >= lX) && ld->PeriodicX) {
 442                                                         ++nBounceBackPdfs; // Compiler bug --> see above
 443                                                 }
 444                                                 else if ((ny < 0 || ny >= lY) && ld->PeriodicY) {
 445                                                         ++nBounceBackPdfs; // Compiler bug --> see above
 446                                                 }
 447                                                 else if ((nz < 0 || nz >= lZ) && ld->PeriodicZ) {
 448                                                         ++nBounceBackPdfs; // Compiler bug --> see above
 449                                                 }
 450                                                 else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
 451                                                         continue;
 452                                                 }
 453                                                 else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
 454                                                         ++nBounceBackPdfs; // Compiler bug --> see above
 455                                                 }
 456                                         }
 457                                 }
 458                         }
 459                 }
 460         }
 461
 462
 463         printf("# allocating %d indices for bounce back pdfs (%s for source and destination array)\n", nBounceBackPdfs, ByteToHuman(sizeof(int) * nBounceBackPdfs * 2));
 464
 465         MemAlloc((void **) & (kd->BounceBackPdfsSrc), sizeof(int) * nBounceBackPdfs + 100);
 466         MemAlloc((void **) & (kd->BounceBackPdfsDst), sizeof(int) * nBounceBackPdfs + 100);
 467
 468         kd->nBounceBackPdfs = nBounceBackPdfs;
 469         nBounceBackPdfs = 0;
 470
 471         int srcIndex;
 472         int dstIndex;
 473
 474         for (int z = 0; z < lZ; ++z) {
 475                 for (int y = 0; y < lY; ++y) {
 476                         for (int x = 0; x < lX; ++x) {
 477
 478                                 if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
 479                                         for (int d = 0; d < N_D3Q19; ++d) {
 480 #ifdef PROP_MODEL_PUSH
 481                                                 nx = x + D3Q19_X[d];
 482                                                 ny = y + D3Q19_Y[d];
 483                                                 nz = z + D3Q19_Z[d];
 484 #elif PROP_MODEL_PULL
 485                                                 nx = x - D3Q19_X[d];
 486                                                 ny = y - D3Q19_Y[d];
 487                                                 nz = z - D3Q19_Z[d];
 488 #else
 489         #error PROP_MODEL_NAME unknown.
 490 #endif
 491
 492                                                 if (    ((nx < 0 || nx >= lX) && ld->PeriodicX) ||
 493                                                                 ((ny < 0 || ny >= lY) && ld->PeriodicY) ||
 494                                                                 ((nz < 0 || nz >= lZ) && ld->PeriodicZ)
 495                                                 ){
 496                                                         // Implement periodic boundary in X direction.
 497
 498                                                         // If the target node reached through propagation is outside the lattice
 499                                                         // the kernel stores it in some buffer around the domain.
 500                                                         // From this position the PDF must be transported to the other side of the
 501                                                         // geometry.
 502
 503                                                         // Take PDF from outside the domain.
 504
 505                                                         // x periodic
 506                                                         if (nx < 0) {
 507                                                                 px = lX - 1;
 508                                                         }
 509                                                         else if (nx >= lX) {
 510                                                                 px = 0;
 511                                                         } else {
 512                                                                 px = nx;
 513                                                         }
 514
 515                                                         // y periodic
 516                                                         if (ny < 0) {
 517                                                                 py = lY - 1;
 518                                                         }
 519                                                         else if (ny >= lY) {
 520                                                                 py = 0;
 521                                                         } else {
 522                                                                 py = ny;
 523                                                         }
 524
 525                                                         // z periodic
 526                                                         if (nz < 0) {
 527                                                                 pz = lZ - 1;
 528                                                         }
 529                                                         else if (nz >= lZ) {
 530                                                                 pz = 0;
 531                                                         } else {
 532                                                                 pz = nz;
 533                                                         }
 534
 535                                                         if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) {
 536 #ifdef PROP_MODEL_PUSH
 537                                                                 srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d);
 538                                                                 dstIndex = P_INDEX_5(gDims,  x + oX,  y + oY,  z + oZ, D3Q19_INV[d]);
 539 #elif PROP_MODEL_PULL
 540                                                                 srcIndex = P_INDEX_5(gDims,  x + oX,  y + oY,  z + oZ, D3Q19_INV[d]);
 541                                                                 dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d);
 542 #endif
 543                                                         }
 544                                                         else {
 545
 546 #ifdef PROP_MODEL_PUSH
 547                                                                 srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d);
 548                                                                 // Put it on the other side back into the domain.
 549                                                                 dstIndex = P_INDEX_5(gDims, px + oX, py + oY, pz + oZ, d);
 550 #elif PROP_MODEL_PULL
 551                                                                 srcIndex = P_INDEX_5(gDims, px + oX, py + oY, pz + oZ, d);
 552                                                                 // Put it on the other side back into the ghost layer.
 553                                                                 dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d);
 554 #endif
 555
 556                                                                 VerifyMsg(nBounceBackPdfs < kd->nBounceBackPdfs, "nBBPdfs %d < kd->nBBPdfs %d  xyz: %d %d %d d: %d\n", nBounceBackPdfs, kd->nBounceBackPdfs, x, y, z, d);
 557
 558                                                         }
 559
 560                                                         kd->BounceBackPdfsSrc[nBounceBackPdfs] = srcIndex;
 561                                                         kd->BounceBackPdfsDst[nBounceBackPdfs] = dstIndex;
 562
 563                                                         ++nBounceBackPdfs;
 564
 565                                                 }
 566                                                 else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) {
 567                                                         continue;
 568                                                 }
 569                                                 else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) {
 570 #ifdef PROP_MODEL_PUSH
 571                                                         srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d);
 572                                                         dstIndex = P_INDEX_5(gDims,  x + oX,  y + oY,  z + oZ, D3Q19_INV[d]);
 573 #elif PROP_MODEL_PULL
 574                                                         srcIndex = P_INDEX_5(gDims,  x + oX,  y + oY,  z + oZ, D3Q19_INV[d]);
 575                                                         dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d);
 576                                                         // srcIndex = P_INDEX_5(gDims,  x + oX,  y + oY,  z + oZ, d);
 577                                                         // dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, D3Q19_INV[d]);
 578 #endif
 579
 580                                                         VerifyMsg(nBounceBackPdfs < kd->nBounceBackPdfs, "nBBPdfs %d < kd->nBBPdfs %d  xyz: %d %d %d d: %d\n", nBounceBackPdfs, kd->nBounceBackPdfs, x, y, z, d);
 581
 582                                                         kd->BounceBackPdfsSrc[nBounceBackPdfs] = srcIndex;
 583                                                         kd->BounceBackPdfsDst[nBounceBackPdfs] = dstIndex;
 584
 585                                                         ++nBounceBackPdfs;
 586                                                 }
 587                                         }
 588                                 }
 589                         }
 590                 }
 591         }
 592
 593
 594         // Fill remaining KernelData structures
 595         kd->GetNode = FNAME(GetNode);
 596         kd->SetNode = FNAME(SetNode);
 597
 598         kd->BoundaryConditionsGetPdf = FNAME(BcGetPdf);
 599         kd->BoundaryConditionsSetPdf = FNAME(BcSetPdf);
 600
 601         kd->Kernel = FNAME(D3Q19BlkKernel);
 602
 603         kd->DstPdfs = NULL;
 604         kd->PdfsActive = kd->Pdfs[0];
 605
 606         return;
 607 }
 608
 609 void FNAME(D3Q19BlkDeinit)(LatticeDesc * ld, KernelData ** kernelData)
 610 {
 611         MemFree((void **) & ((*kernelData)->Pdfs[0]));
 612         MemFree((void **) & ((*kernelData)->Pdfs[1]));
 613
 614         MemFree((void **) & ((*kernelData)->BounceBackPdfsSrc));
 615         MemFree((void **) & ((*kernelData)->BounceBackPdfsDst));
 616
 617         MemFree((void **)kernelData);
 618
 619         return;
 620 }
 621
 622 // Kernels without blocking perform the same initialization/deinitialization as with
 623 // blocking, except that a different kernel is called. Hence, no arguments are allowed.
 624
 625 void FNAME(D3Q19Init)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params)
 626 {
 627         Parameters p;
 628
 629         if (params->nKernelArgs != 0) {
 630                 printf("ERROR: unknown kernel parameter.\n");
 631                 printf("This kernels accepts no parameters.\n");
 632                 exit(1);
 633         }
 634
 635         // Setup an empty parameters structure.
 636         p.nArgs        = params->nArgs;
 637         p.Args         = params->Args;
 638         p.nKernelArgs  = 0;
 639         p.KernelArgs   = NULL;
 640
 641         // Call init routine for blocking kernel and override the
 642         // kernel function to be called later on.
 643         FNAME(D3Q19BlkInit)(ld, kernelData, &p);
 644
 645         (*kernelData)->Kernel = FNAME(D3Q19Kernel);
 646
 647         return;
 648
 649 }
 650
 651 void FNAME(D3Q19Deinit)(LatticeDesc * ld, KernelData ** kernelData)
 652 {
 653         FNAME(D3Q19BlkDeinit)(ld, kernelData);
 654         return;
 655 }