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