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