src/BenchKernelD3Q19ListAaPvGatherAoSoA.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 //   Michael Hussnaetter, 2017-2018
  12 //   University of Erlangen-Nuremberg, Germany
  13 //   michael.hussnaetter -at- fau.de
  14 //
  15 //  This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
  16 //
  17 //  LbmBenchKernels is free software: you can redistribute it and/or modify
  18 //  it under the terms of the GNU General Public License as published by
  19 //  the Free Software Foundation, either version 3 of the License, or
  20 //  (at your option) any later version.
  21 //
  22 //  LbmBenchKernels is distributed in the hope that it will be useful,
  23 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  24 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  25 //  GNU General Public License for more details.
  26 //
  27 //  You should have received a copy of the GNU General Public License
  28 //  along with LbmBenchKernels.  If not, see <http://www.gnu.org/licenses/>.
  29 //
  30 // --------------------------------------------------------------------------
  31 #include "BenchKernelD3Q19ListAaPvGatherAoSoACommon.h"
  32
  33 #include "Base.h"
  34 #include "Memory.h"
  35 #include "Vtk.h"
  36 #include "Vector.h"
  37
  38 #include <inttypes.h>
  39 #include <math.h>
  40
  41 #ifdef _OPENMP
  42 #include <omp.h>
  43 #endif
  44
  45 #ifdef LIKWID_PERFMON
  46 #include <likwid.h>
  47 #else
  48 #define LIKWID_MARKER_INIT
  49 #define LIKWID_MARKER_THREADINIT
  50 #define LIKWID_MARKER_SWITCH
  51 #define LIKWID_MARKER_REGISTER(regionTag)
  52 #define LIKWID_MARKER_START(regionTag)
  53 #define LIKWID_MARKER_STOP(regionTag)
  54 #define LIKWID_MARKER_CLOSE
  55 #define LIKWID_MARKER_GET(regionTag, nevents, events, time, count)
  56 #endif
  57
  58 //enable software prefetchting for vectorized gather/scatter loop in odd kernel
  59 #ifndef SOFTWARE_PREFETCH_LOOKAHEAD_L2
  60 #define SOFTWARE_PREFETCH_LOOKAHEAD_L2 (0) //prefetchting X SIMD widths ahead
  61 #endif
  62
  63 #ifndef SOFTWARE_PREFETCH_LOOKAHEAD_L1
  64 #define SOFTWARE_PREFETCH_LOOKAHEAD_L1 (0) //prefetchting X SIMD widths ahead
  65 #endif
  66
  67 static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
  68 static void KernelOdd( LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
  69
  70 void FNAME(D3Q19ListAaPvGatherAoSoAKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
  71 {
  72
  73         Assert(ld != NULL);
  74         Assert(kernelData != NULL);
  75         Assert(cd != NULL);
  76
  77         Assert(cd->Omega > 0.0);
  78         Assert(cd->Omega < 2.0);
  79
  80 #if defined(VTK_OUTPUT) || defined(STATISTICS) || defined(VERIFICATION)
  81         KernelData * kd = (KernelData *)kernelData;
  82 #endif
  83         KernelDataList * kdl = KDL(kernelData);
  84
  85         int maxIterations = cd->MaxIterations;
  86         int nFluid = kdl->nFluid;
  87
  88         printf("\n");
  89 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0) || (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
  90         printf("# Software prefetching enabled:\n");
  91         printf("#   Gather/Scatter prefetch lookahead L2: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L2);
  92         printf("#   Gather/Scatter prefetch lookahead L1: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L1);
  93 #else
  94         printf("# Software prefetching disabled.\n");
  95 #endif
  96         printf("\n");
  97
  98         int nThreads = 1;
  99 #ifdef _OPENMP
 100         nThreads = omp_get_max_threads();
 101 #endif
 102
 103         int * threadIndices = (int *)malloc(sizeof(int) * (nThreads + 1));
 104         for (int i = 0; i < nThreads; ++i) {
 105                 threadIndices[i] = i * (nFluid / nThreads) + MinI(i, nFluid % nThreads);
 106         }
 107
 108         threadIndices[nThreads] = nFluid;
 109
 110 #ifdef VTK_OUTPUT
 111         if (cd->VtkOutput) {
 112                 kd->PdfsActive = kd->Pdfs[0];
 113                 VtkWrite(ld, kd, cd, -1);
 114         }
 115 #endif
 116
 117 #ifdef STATISTICS
 118         kd->PdfsActive = kd->Pdfs[0];
 119         KernelStatistics(kd, ld, cd, 0);
 120 #endif
 121
 122         LIKWID_MARKER_INIT;
 123
 124         X_KERNEL_START(kernelData);
 125
 126         // TODO: outer openmp parallel
 127
 128         LIKWID_MARKER_START("OuterLoop");
 129         for(int iter = 0; iter < maxIterations; iter += 2) {
 130
 131                 // even time step
 132
 133 #ifdef _OPENMP
 134 #pragma omp parallel default(none) shared(ld, kernelData, cd, threadIndices)
 135 #endif
 136                 {
 137                         KernelEven(ld, kernelData, cd, threadIndices);
 138                 }
 139
 140
 141 #ifdef VERIFICATION
 142                 kdl->Iteration = iter;
 143                 kd->PdfsActive = kd->Pdfs[0];
 144                 KernelAddBodyForce(kd, ld, cd);
 145 #endif
 146
 147                 // odd time step
 148
 149 #ifdef _OPENMP
 150 #pragma omp parallel default(none) shared(ld, kernelData, cd, threadIndices)
 151 #endif
 152                 {
 153                         KernelOdd(ld, kernelData, cd, threadIndices);
 154                 }
 155
 156
 157 #ifdef VERIFICATION
 158                 kdl->Iteration = iter + 1;
 159                 kd->PdfsActive = kd->Pdfs[0];
 160                 KernelAddBodyForce(kd, ld, cd);
 161 #endif
 162
 163 #ifdef VTK_OUTPUT
 164                 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
 165                         kdl->Iteration = iter + 1;
 166                         kd->PdfsActive = kd->Pdfs[0];
 167                         VtkWrite(ld, kd, cd, iter);
 168                 }
 169 #endif
 170
 171 #ifdef STATISTICS
 172                 kdl->Iteration = iter + 1;
 173                 kd->PdfsActive = kd->Pdfs[0];
 174                 KernelStatistics(kd, ld, cd, iter);
 175 #endif
 176
 177         } // for (int iter = 0; ...
 178         LIKWID_MARKER_STOP("OuterLoop");
 179
 180         X_KERNEL_END(kernelData);
 181
 182 #ifdef VTK_OUTPUT
 183         if (cd->VtkOutput) {
 184                 kd->PdfsActive = kd->Pdfs[0];
 185                 VtkWrite(ld, kd, cd, maxIterations);
 186         }
 187 #endif
 188
 189 #ifdef STATISTICS
 190         kd->PdfsActive = kd->Pdfs[0];
 191         KernelStatistics(kd, ld, cd, maxIterations);
 192 #endif
 193
 194         LIKWID_MARKER_CLOSE;
 195         free(threadIndices);
 196
 197         return;
 198 }
 199
 200 static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
 201 {
 202         Assert(ld != NULL);
 203         Assert(kernelData != NULL);
 204         Assert(cd != NULL);
 205
 206         Assert(cd->Omega > 0.0);
 207         Assert(cd->Omega < 2.0);
 208
 209         KernelData * kd = (KernelData *)kernelData;
 210         KernelDataList * kdl = KDL(kernelData);
 211
 212         PdfT omega = cd->Omega;
 213         PdfT omegaEven = omega;
 214
 215         PdfT magicParam = 1.0 / 12.0;
 216         PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
 217
 218         PdfT evenPart = 0.0;
 219         PdfT oddPart = 0.0;
 220         PdfT dir_indep_trm = 0.0;
 221
 222         const PdfT w_0 = 1.0 /  3.0;
 223         const PdfT w_1 = 1.0 / 18.0;
 224         const PdfT w_2 = 1.0 / 36.0;
 225
 226         const PdfT w_1_x3 = w_1 * 3.0;  const PdfT w_1_nine_half = w_1 * 9.0 / 2.0;     PdfT w_1_indep = 0.0;
 227         const PdfT w_2_x3 = w_2 * 3.0;  const PdfT w_2_nine_half = w_2 * 9.0 / 2.0;     PdfT w_2_indep = 0.0;
 228
 229         PdfT ux, uy, uz, ui;
 230         PdfT dens;
 231
 232         VPDFT VONE_HALF = VSET(0.5);
 233         VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
 234
 235         VPDFT vw_1_indep, vw_2_indep;
 236         VPDFT vw_0 = VSET(w_0);
 237         VPDFT vw_1 = VSET(w_1);
 238         VPDFT vw_2 = VSET(w_2);
 239
 240         VPDFT vw_1_x3 = VSET(w_1_x3);
 241         VPDFT vw_2_x3 = VSET(w_2_x3);
 242         VPDFT vw_1_nine_half = VSET(w_1_nine_half);
 243         VPDFT vw_2_nine_half = VSET(w_2_nine_half);
 244
 245         VPDFT vui, vux, vuy, vuz, vdens;
 246
 247         VPDFT vevenPart, voddPart, vdir_indep_trm;
 248
 249         VPDFT vomegaEven = VSET(omegaEven);
 250         VPDFT vomegaOdd  = VSET(omegaOdd);
 251
 252         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
 253         #define X(name, idx, idxinv, x, y, z) \
 254                 PdfT JOIN(pdf_,name); \
 255                 PdfT * JOIN(ppdf_,name); \
 256                 VPDFT JOIN(vpdf_,name);
 257                 D3Q19_LIST
 258         #undef X
 259
 260         PdfT * src = kd->Pdfs[0];
 261
 262         int nCells = kdl->nCells;
 263
 264         int threadId = 0;
 265 #ifdef _OPENMP
 266         threadId =  omp_get_thread_num();
 267 #endif
 268
 269
 270         int indexStart    = threadIndices[threadId];
 271         int nFluidThread  = threadIndices[threadId + 1] - threadIndices[threadId];
 272         int indexStop     = indexStart + nFluidThread;
 273
 274         int indexStartVec = ((indexStart + VSIZE - 1) / VSIZE) * VSIZE;
 275         int nFluidVec     = (indexStop / VSIZE) * VSIZE - indexStartVec;
 276         int indexStopVec  = indexStartVec + nFluidVec;
 277
 278         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 279
 280         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[I(indexStart, idx)]);
 281                         D3Q19_LIST
 282         #undef X
 283
 284         for (int index = indexStart; index < indexStartVec; ++index) {
 285
 286                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 287                         D3Q19_LIST
 288                 #undef X
 289
 290                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 291                         pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 292                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 293                         pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 294                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 295                         pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 296
 297                 dens = pdf_C +
 298                         pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 299                         pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 300                         pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 301                         pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 302
 303                 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 304
 305                 // direction: w_0
 306                 *ppdf_C  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
 307
 308                 // direction: w_1
 309                 w_1_indep = w_1*dir_indep_trm;
 310
 311                 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 312                         ui = tmpUi; \
 313                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 314                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 315                         *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 316                         *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 317
 318                 COLLIDE_AA_S(uy, N, S)
 319                 COLLIDE_AA_S(ux, E, W)
 320                 COLLIDE_AA_S(uz, T, B)
 321
 322                 #undef COLLIDE_AA_S
 323
 324                 // direction: w_2
 325                 w_2_indep = w_2*dir_indep_trm;
 326
 327                 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 328                         ui = tmpUi; \
 329                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 330                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 331                         *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 332                         *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 333
 334                 COLLIDE_UA_S((-ux + uy), NW, SE)
 335                 COLLIDE_UA_S(( ux + uy), NE, SW)
 336                 COLLIDE_UA_S((-ux + uz), TW, BE)
 337                 COLLIDE_UA_S(( ux + uz), TE, BW)
 338                 COLLIDE_UA_S((-uy + uz), TS, BN)
 339                 COLLIDE_UA_S(( uy + uz), TN, BS)
 340
 341                 #undef COLLIDE_UA_S
 342
 343                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name)++;
 344                         D3Q19_LIST
 345                 #undef X
 346         }
 347
 348         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[I(indexStartVec, idx)]);
 349                         D3Q19_LIST
 350         #undef X
 351
 352         for (int index = indexStartVec; index < indexStopVec; index += VSIZE) {
 353
 354                 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 355                         #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE * N_D3Q19), _MM_HINT_T1);
 356                                 D3Q19_LIST
 357                         #undef X
 358                 #endif
 359
 360                 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 361                         #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE * N_D3Q19), _MM_HINT_T0);
 362                                 D3Q19_LIST
 363                         #undef X
 364                 #endif
 365
 366                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
 367                         D3Q19_LIST
 368                 #undef X
 369
 370                 //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 371                 //           vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 372                 vux = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_E,VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_TE,vpdf_BE)),vpdf_W),vpdf_NW),vpdf_SW),vpdf_TW),vpdf_BW);
 373                 //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 374                 //           vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 375                 vuy = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_N,VADD(vpdf_NE,vpdf_NW)),VADD(vpdf_TN,vpdf_BN)),vpdf_S),vpdf_SE),vpdf_SW),vpdf_TS),vpdf_BS);
 376                 //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 377                 //           vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 378                 vuz = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_T,VADD(vpdf_TE,vpdf_TW)),VADD(vpdf_TN,vpdf_TS)),vpdf_B),vpdf_BE),vpdf_BW),vpdf_BN),vpdf_BS);
 379
 380                 //vdens = vpdf_C +
 381                 //          vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 382                 //          vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 383                 //          vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 384                 //          vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 385                 vdens = VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(vpdf_C,VADD(vpdf_N,vpdf_E)),VADD(vpdf_S,vpdf_W)),VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
 386
 387                 //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 388                 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 389
 390                 //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 391                 VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 392
 393                 //vw_1_indep = vw_1 * vdir_indep_trm;
 394                 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 395
 396                 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 397                         vui = tmpVui; \
 398                         vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));\
 399                         voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 400                         VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 401                         VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 402
 403                 COLLIDE_AA_V(vuy, N, S)
 404                 COLLIDE_AA_V(vux, E, W)
 405                 COLLIDE_AA_V(vuz, T, B)
 406
 407                 #undef COLLIDE_AA_V
 408
 409                 //vw_2_indep = vw_2 * vdir_indep_trm;
 410                 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 411
 412                 // collide axis unaligned pdfs vectorized
 413                 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 414                         vui = tmpVui; \
 415                         vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));\
 416                         voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 417                         VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 418                         VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 419
 420                 COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 421                 COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 422                 COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 423                 COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 424                 COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 425                 COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 426
 427                 #undef COLLIDE_UA_V
 428
 429                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) +=(VSIZE * N_D3Q19);
 430                         D3Q19_LIST
 431                 #undef X
 432         } // loop over fluid nodes
 433
 434         for (int index = indexStopVec; index < indexStop; ++index) {
 435
 436                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 437                         D3Q19_LIST
 438                 #undef X
 439
 440                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 441                         pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 442                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 443                         pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 444                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 445                         pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 446
 447                 dens = pdf_C +
 448                         pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 449                         pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 450                         pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 451                         pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 452
 453                 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 454
 455                 // direction: w_0
 456                 *ppdf_C  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
 457
 458                 // direction: w_1
 459                 w_1_indep = w_1*dir_indep_trm;
 460
 461                 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 462                         ui = tmpUi; \
 463                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 464                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 465                         *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 466                         *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 467
 468                 COLLIDE_AA_S(uy, N, S)
 469                 COLLIDE_AA_S(ux, E, W)
 470                 COLLIDE_AA_S(uz, T, B)
 471
 472                 #undef COLLIDE_AA_S
 473
 474                 // direction: w_2
 475                 w_2_indep = w_2*dir_indep_trm;
 476
 477                 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 478                         ui = tmpUi; \
 479                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 480                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 481                         *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 482                         *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 483
 484                 COLLIDE_UA_S((-ux + uy), NW, SE)
 485                 COLLIDE_UA_S(( ux + uy), NE, SW)
 486                 COLLIDE_UA_S((-ux + uz), TW, BE)
 487                 COLLIDE_UA_S(( ux + uz), TE, BW)
 488                 COLLIDE_UA_S((-uy + uz), TS, BN)
 489                 COLLIDE_UA_S(( uy + uz), TN, BS)
 490
 491                 #undef COLLIDE_UA_S
 492
 493                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name)++;
 494                         D3Q19_LIST
 495                 #undef X
 496         } // loop over fluid nodes
 497
 498         #undef I
 499
 500         return;
 501 }
 502
 503 static void KernelOdd(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
 504 {
 505
 506         Assert(ld != NULL);
 507         Assert(kernelData != NULL);
 508         Assert(cd != NULL);
 509
 510         Assert(cd->Omega > 0.0);
 511         Assert(cd->Omega < 2.0);
 512
 513         KernelData * kd = (KernelData *)kernelData;
 514         KernelDataList * kdl = KDL(kernelData);
 515         KernelDataListRia * kdlr = KDLR(kernelData);
 516         PdfT omega = cd->Omega;
 517         PdfT omegaEven = omega;
 518
 519         PdfT magicParam = 1.0 / 12.0;
 520         PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
 521
 522         PdfT evenPart = 0.0;
 523         PdfT oddPart = 0.0;
 524         PdfT dir_indep_trm = 0.0;
 525
 526         const PdfT w_0 = 1.0 /  3.0;
 527         const PdfT w_1 = 1.0 / 18.0;
 528         const PdfT w_2 = 1.0 / 36.0;
 529
 530         const PdfT w_1_x3 = w_1 * 3.0;  const PdfT w_1_nine_half = w_1 * 9.0 / 2.0;     PdfT w_1_indep = 0.0;
 531         const PdfT w_2_x3 = w_2 * 3.0;  const PdfT w_2_nine_half = w_2 * 9.0 / 2.0;     PdfT w_2_indep = 0.0;
 532
 533         PdfT ux, uy, uz, ui;
 534         PdfT dens;
 535
 536         VPDFT VONE_HALF = VSET(0.5);
 537         VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
 538
 539         VPDFT vw_1_indep, vw_2_indep;
 540         VPDFT vw_0 = VSET(w_0);
 541         VPDFT vw_1 = VSET(w_1);
 542         VPDFT vw_2 = VSET(w_2);
 543
 544         VPDFT vw_1_x3 = VSET(w_1_x3);
 545         VPDFT vw_2_x3 = VSET(w_2_x3);
 546         VPDFT vw_1_nine_half = VSET(w_1_nine_half);
 547         VPDFT vw_2_nine_half = VSET(w_2_nine_half);
 548
 549         VPDFT vux, vuy, vuz, vui;
 550         VPDFT vdens;
 551
 552         VPDFT vevenPart, voddPart, vdir_indep_trm;
 553
 554         VPDFT vomegaEven = VSET(omegaEven);
 555         VPDFT vomegaOdd  = VSET(omegaOdd);
 556
 557         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
 558         #define X(name, idx, idxinv, x, y, z) \
 559                 PdfT JOIN(pdf_,name) = 0; \
 560                 PdfT * JOIN(ppdf_,name) = NULL; \
 561                 VPDFT JOIN(vpdf_,name);
 562                 D3Q19_LIST
 563         #undef X
 564
 565         #define X(name, idx, idxinv, x, y, z) \
 566                 __m256i JOIN(vgatheridx_,name) = _mm256_set1_epi32(0);
 567                 D3Q19_LIST_WO_C
 568         #undef X
 569
 570         __m256i vgatherinc = VSETI32(VSIZE * N_D3Q19);
 571
 572         uint32_t * consecNodes = kdlr->ConsecNodes;
 573         uint32_t consecIndex = 0;
 574         uint32_t consecValue = 0;
 575
 576         PdfT * src = kd->Pdfs[0];
 577
 578         int nCells = kdl->nCells;
 579
 580         int adjListIndex;
 581         uint32_t * adjList = kdl->AdjList;
 582
 583         int threadId = 0;
 584
 585 #ifdef _OPENMP
 586         threadId = omp_get_thread_num();
 587 #endif
 588         consecIndex = kdlr->ConsecThreadIndices[threadId];
 589         consecValue = 0;
 590
 591         int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
 592
 593         int indexStart = threadIndices[threadId];
 594         int indexStop  = threadIndices[threadId] + nFluidThread;
 595
 596         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 597         #define ADJ_LIST(dir) adjList[adjListIndex + (dir * VSIZE)]
 598
 599         int offset_ppdf_C = -1; //dummy init to detect errors
 600
 601         for (int index = indexStart; index < indexStop; index += 1) {
 602
 603                 if (consecValue > 0) {
 604                         --consecValue;
 605                         // Increment all pdf pointers by an offset. If the previous iteration was
 606                         // scalar, increment only by one. If the previous iteration was vectorized,
 607                         // increment by the vector width. These offsets are set in the corresponding
 608                         // if branches.
 609
 610                         //increment offsets
 611
 612                         #define X(name, idx, idxinv, _x, _y, _z) JOIN(vgatheridx_,name) = VADDI32(JOIN(vgatheridx_,name), vgatherinc);
 613                                 D3Q19_LIST_WO_C
 614                         #undef X
 615
 616                         //printVector(vgatheridx_N);
 617
 618                         ppdf_C += offset_ppdf_C;
 619
 620                 }
 621                 else {
 622                         // Load new pointers to PDFs of local cell:
 623                         Assert(consecIndex < nConsecNodes);
 624
 625                         consecValue = consecNodes[consecIndex] - 1;
 626
 627                         adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 628                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vgatheridx_,name) = VLIU(&(ADJ_LIST(idxinv)));
 629                                 D3Q19_LIST_WO_C
 630                         #undef X
 631
 632                         ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
 633                         ++consecIndex;
 634                 }
 635
 636                 if (consecValue >= (VSIZE - 1)) {
 637                         // Vectorized part.
 638
 639                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 640                                 int const indexPrefetchL2 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L2;
 641                                 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
 642                                 if (indexPrefetchL2 < indexStop){
 643                                         // update pointers from adjacency list only if necessary
 644                                         if (consecValue >= (SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE + VSIZE - 1)) {
 645                                                 #define INCR_PTR(name)          (VADDI32(JOIN(vgatheridx_,name), VMULI32(vgatherinc, VSETI32(SOFTWARE_PREFETCH_LOOKAHEAD_L2))))
 646                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(INCR_PTR(name), (char const *) src, 8, _MM_HINT_T1);
 647                                                         D3Q19_LIST_WO_C
 648                                                 #undef X
 649                                                 #undef INCR_PTR
 650                                         }
 651                                         else {
 652                                                 adjListIndex = (indexPrefetchL2 - (indexPrefetchL2 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL2 % VSIZE);
 653                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T1);
 654                                                         D3Q19_LIST_WO_C
 655                                                 #undef X
 656                                         }
 657
 658                                         _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL2, D3Q19_C)]), _MM_HINT_T1);
 659                                 }
 660                         #endif
 661
 662                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 663                                 int const indexPrefetchL1 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L1;
 664                                 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
 665                                 if (indexPrefetchL1 < indexStop){
 666                                         // update pointers from adjacency list only if necessary
 667                                         if (consecValue > (SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE + VSIZE - 1)) {
 668                                                 #define INCR_PTR(name)          (VADDI32(JOIN(vgatheridx_,name), VMULI32(vgatherinc, VSETI32(SOFTWARE_PREFETCH_LOOKAHEAD_L1))))
 669                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(INCR_PTR(name), (char const *) src, 8, _MM_HINT_T0);
 670                                                         D3Q19_LIST_WO_C
 671                                                 #undef X
 672                                                 #undef INCR_PTR
 673                                         }
 674                                         else {
 675                                                 adjListIndex = (indexPrefetchL1 - (indexPrefetchL1 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL1 % VSIZE);
 676                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T0);
 677                                                         D3Q19_LIST_WO_C
 678                                                 #undef X
 679                                         }
 680
 681                                         _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL1, D3Q19_C)]), _MM_HINT_T0);
 682                                 }
 683                         #endif
 684
 685                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VG32(JOIN(vgatheridx_,name), src, 8);
 686                                 D3Q19_LIST_WO_C
 687                         #undef X
 688
 689                         vpdf_C = VLDU(ppdf_C);
 690
 691                         //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 692                         //      vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 693                         vux = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_E,VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_TE,vpdf_BE)),vpdf_W),vpdf_NW),vpdf_SW),vpdf_TW),vpdf_BW);
 694                         //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 695                         //      vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 696                         vuy = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_N,VADD(vpdf_NE,vpdf_NW)),VADD(vpdf_TN,vpdf_BN)),vpdf_S),vpdf_SE),vpdf_SW),vpdf_TS),vpdf_BS);
 697                         //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 698                         //      vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 699                         vuz = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_T,VADD(vpdf_TE,vpdf_TW)),VADD(vpdf_TN,vpdf_TS)),vpdf_B),vpdf_BE),vpdf_BW),vpdf_BN),vpdf_BS);
 700
 701                         //vdens = vpdf_C +
 702                         //        vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 703                         //        vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 704                         //        vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 705                         //        vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 706                         vdens = VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(vpdf_C,VADD(vpdf_N,vpdf_E)),VADD(vpdf_S,vpdf_W)),VADD(vpdf_NE,vpdf_SE)),
 707                                                                                 VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
 708
 709                         //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 710                         vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 711
 712                         //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 713                         VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 714
 715                         // collide axis aligend pdfs vectorized
 716                         #define SCAT(offsets, vsrc) VS32(src, offsets, vsrc, 8)
 717
 718                         //vw_1_indep = vw_1 * vdir_indep_trm;
 719                         vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 720
 721                         // collide axis aligend pdfs vectorized
 722                         #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 723                                 vui = tmpVui; \
 724                                 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));\
 725                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 726                                 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 727                                 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 728
 729                         COLLIDE_AA_V(vuy, N, S)
 730                         COLLIDE_AA_V(vux, E, W)
 731                         COLLIDE_AA_V(vuz, T, B)
 732
 733                         #undef COLLIDE_AA_V
 734
 735                         //vw_2_indep = vw_2 * vdir_indep_trm;
 736                         vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 737
 738                         // collide axis unaligned pdfs vectorized
 739                         #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 740                                 vui = tmpVui; \
 741                                 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));\
 742                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 743                                 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 744                                 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 745
 746                         COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 747                         COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 748                         COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 749                         COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 750                         COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 751                         COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 752
 753                         #undef COLLIDE_UA_V
 754                         #undef SCAT
 755
 756                         consecValue   -= (VSIZE - 1);
 757                         index         += (VSIZE - 1);
 758                         offset_ppdf_C  = VSIZE * N_D3Q19;
 759
 760                 }
 761                 else {
 762                         // Scalar part.
 763
 764                         adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 765                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
 766                                 D3Q19_LIST_WO_C
 767                         #undef X
 768                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 769                                 D3Q19_LIST_WO_C
 770                         #undef X
 771
 772                         pdf_C = *ppdf_C;
 773
 774                         ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 775                                 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 776                         uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 777                                 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 778                         uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 779                                 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 780
 781                         dens = pdf_C +
 782                                 pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 783                                 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 784                                 pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 785                                 pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 786
 787                         dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 788
 789                         // direction: w_0
 790                         *ppdf_C = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
 791
 792                         // direction: w_1
 793                         w_1_indep = w_1 * dir_indep_trm;
 794
 795                         #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 796                                 ui = tmpUi; \
 797                                 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 798                                 oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 799                                 *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 800                                 *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 801
 802                         COLLIDE_AA_S(uy, N, S)
 803                         COLLIDE_AA_S(ux, E, W)
 804                         COLLIDE_AA_S(uz, T, B)
 805
 806                         #undef COLLIDE_AA_S
 807
 808                         // direction: w_2
 809                         w_2_indep = w_2 * dir_indep_trm;
 810
 811                         #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 812                                 ui = tmpUi; \
 813                                 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 814                                 oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 815                                 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 816                                 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 817
 818                         COLLIDE_UA_S((-ux + uy), NW, SE)
 819                         COLLIDE_UA_S(( ux + uy), NE, SW)
 820                         COLLIDE_UA_S((-ux + uz), TW, BE)
 821                         COLLIDE_UA_S(( ux + uz), TE, BW)
 822                         COLLIDE_UA_S((-uy + uz), TS, BN)
 823                         COLLIDE_UA_S(( uy + uz), TN, BS)
 824
 825                         #undef COLLIDE_UA_S
 826
 827                         offset_ppdf_C = 1;
 828                 }
 829
 830         } // loop over fluid nodes
 831
 832 #undef ADJ_LIST
 833 #undef I
 834 }