src/BenchKernelD3Q19ListAaPvGather.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 "BenchKernelD3Q19ListAaPvGatherCommon.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(D3Q19ListAaPvGatherKernel)(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("#   Load/Store prefetch lookahead L2: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L2);
  92         printf("#   Load/Store 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         // TODO: outer openmp parallel
 125
 126         X_KERNEL_START(kernelData);
 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         int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
 270         int nFluidVec    = nFluidThread - (nFluidThread % VSIZE);
 271
 272         int indexStartVec = threadIndices[threadId];
 273         int indexStopVec  = threadIndices[threadId] + nFluidVec;
 274         int indexStop     = threadIndices[threadId] + nFluidThread;
 275
 276         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 277
 278         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[I(indexStartVec, idx)]);
 279                         D3Q19_LIST
 280         #undef X
 281
 282         for (int index = indexStartVec; index < indexStopVec; index += VSIZE) {
 283
 284                 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 285                         #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE), _MM_HINT_T1);
 286                                 D3Q19_LIST
 287                         #undef X
 288                 #endif
 289
 290                 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 291                         #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE), _MM_HINT_T0);
 292                                 D3Q19_LIST
 293                         #undef X
 294                 #endif
 295
 296                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
 297                         D3Q19_LIST
 298                 #undef X
 299
 300                 //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 301                 //           vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 302                 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);
 303                 //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 304                 //           vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 305                 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);
 306                 //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 307                 //           vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 308                 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);
 309
 310                 //vdens = vpdf_C +
 311                 //          vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 312                 //          vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 313                 //          vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 314                 //          vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 315                 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));
 316
 317                 //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 318                 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 319
 320                 //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 321                 VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 322
 323                 //vw_1_indep = vw_1 * vdir_indep_trm;
 324                 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 325
 326                 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 327                         vui = tmpVui; \
 328                         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));\
 329                         voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 330                         VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 331                         VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 332
 333                 COLLIDE_AA_V(vuy, N, S)
 334                 COLLIDE_AA_V(vux, E, W)
 335                 COLLIDE_AA_V(vuz, T, B)
 336
 337                 #undef COLLIDE_AA_V
 338
 339                 //vw_2_indep = vw_2 * vdir_indep_trm;
 340                 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 341
 342                 // collide axis unaligned pdfs vectorized
 343                 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 344                         vui = tmpVui; \
 345                         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));\
 346                         voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 347                         VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 348                         VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 349
 350                 COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 351                 COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 352                 COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 353                 COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 354                 COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 355                 COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 356
 357                 #undef COLLIDE_UA_V
 358
 359                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) +=VSIZE;
 360                         D3Q19_LIST
 361                 #undef X
 362         } // loop over fluid nodes
 363
 364         for (int index = indexStopVec; index < indexStop; ++index) {
 365
 366                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 367                         D3Q19_LIST
 368                 #undef X
 369
 370                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 371                         pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 372                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 373                         pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 374                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 375                         pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 376
 377                 dens = pdf_C +
 378                         pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 379                         pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 380                         pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 381                         pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 382
 383                 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 384
 385                 // direction: w_0
 386                 *ppdf_C  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
 387
 388                 // direction: w_1
 389                 w_1_indep = w_1*dir_indep_trm;
 390
 391                 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 392                         ui = tmpUi; \
 393                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 394                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 395                         *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 396                         *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 397
 398                 COLLIDE_AA_S(uy, N, S)
 399                 COLLIDE_AA_S(ux, E, W)
 400                 COLLIDE_AA_S(uz, T, B)
 401
 402                 #undef COLLIDE_AA_S
 403
 404                 // direction: w_2
 405                 w_2_indep = w_2*dir_indep_trm;
 406
 407                 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 408                         ui = tmpUi; \
 409                         evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 410                         oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 411                         *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 412                         *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 413
 414                 COLLIDE_UA_S((-ux + uy), NW, SE)
 415                 COLLIDE_UA_S(( ux + uy), NE, SW)
 416                 COLLIDE_UA_S((-ux + uz), TW, BE)
 417                 COLLIDE_UA_S(( ux + uz), TE, BW)
 418                 COLLIDE_UA_S((-uy + uz), TS, BN)
 419                 COLLIDE_UA_S(( uy + uz), TN, BS)
 420
 421                 #undef COLLIDE_UA_S
 422
 423                 #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name)++;
 424                         D3Q19_LIST
 425                 #undef X
 426         } // loop over fluid nodes
 427
 428         #undef I
 429
 430         return;
 431 }
 432
 433 static void KernelOdd(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
 434 {
 435
 436         Assert(ld != NULL);
 437         Assert(kernelData != NULL);
 438         Assert(cd != NULL);
 439
 440         Assert(cd->Omega > 0.0);
 441         Assert(cd->Omega < 2.0);
 442
 443         KernelData * kd = (KernelData *)kernelData;
 444         KernelDataList * kdl = KDL(kernelData);
 445         KernelDataListRia * kdlr = KDLR(kernelData);
 446         PdfT omega = cd->Omega;
 447         PdfT omegaEven = omega;
 448
 449         PdfT magicParam = 1.0 / 12.0;
 450         PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
 451
 452         PdfT evenPart = 0.0;
 453         PdfT oddPart = 0.0;
 454         PdfT dir_indep_trm = 0.0;
 455
 456         const PdfT w_0 = 1.0 /  3.0;
 457         const PdfT w_1 = 1.0 / 18.0;
 458         const PdfT w_2 = 1.0 / 36.0;
 459
 460         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;
 461         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;
 462
 463         PdfT ux, uy, uz, ui;
 464         PdfT dens;
 465
 466         VPDFT VONE_HALF = VSET(0.5);
 467         VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
 468
 469         VPDFT vw_1_indep, vw_2_indep;
 470         VPDFT vw_0 = VSET(w_0);
 471         VPDFT vw_1 = VSET(w_1);
 472         VPDFT vw_2 = VSET(w_2);
 473
 474         VPDFT vw_1_x3 = VSET(w_1_x3);
 475         VPDFT vw_2_x3 = VSET(w_2_x3);
 476         VPDFT vw_1_nine_half = VSET(w_1_nine_half);
 477         VPDFT vw_2_nine_half = VSET(w_2_nine_half);
 478
 479         VPDFT vux, vuy, vuz, vui;
 480         VPDFT vdens;
 481
 482         VPDFT vevenPart, voddPart, vdir_indep_trm;
 483
 484         VPDFT vomegaEven = VSET(omegaEven);
 485         VPDFT vomegaOdd  = VSET(omegaOdd);
 486
 487         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
 488         #define X(name, idx, idxinv, x, y, z) \
 489                 PdfT JOIN(pdf_,name) = 0; \
 490                 PdfT * JOIN(ppdf_,name) = NULL; \
 491                 VPDFT JOIN(vpdf_,name);
 492                 D3Q19_LIST
 493         #undef X
 494
 495         #define X(name, idx, idxinv, x, y, z) \
 496                 __m256i JOIN(vgatheridx_,name) = _mm256_set1_epi32(0);
 497                 D3Q19_LIST_WO_C
 498         #undef X
 499
 500         __m256i vgatherinc = _mm256_set1_epi32(VSIZE);
 501
 502         uint32_t * consecNodes = kdlr->ConsecNodes;
 503         uint32_t consecIndex = 0;
 504         uint32_t consecValue = 0;
 505
 506         PdfT * src = kd->Pdfs[0];
 507
 508         int nCells = kdl->nCells;
 509
 510         int adjListIndex;
 511         uint32_t * adjList = kdl->AdjList;
 512
 513         int threadId = 0;
 514
 515 #ifdef _OPENMP
 516         threadId = omp_get_thread_num();
 517 #endif
 518         consecIndex = kdlr->ConsecThreadIndices[threadId];
 519         consecValue = 0;
 520
 521         int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
 522
 523         int indexStart = threadIndices[threadId];
 524         int indexStop  = threadIndices[threadId] + nFluidThread;
 525
 526         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 527         #define ADJ_LIST(dir) adjList[adjListIndex + (dir * VSIZE)]
 528
 529         int offset_ppdf_C = -1; //dummy init to detect wrong usage.
 530
 531         for (int index = indexStart; index < indexStop; index += 1) {
 532
 533                 if (consecValue > 0) {
 534                         --consecValue;
 535                         // Increment all pdf pointers by an offset. If the previous iteration was
 536                         // scalar, increment only by one. If the previous iteration was vectorized,
 537                         // increment by the vector width. These offsets are set in the corresponding
 538                         // if branches.
 539
 540                         //increment offsets
 541
 542                         #define X(name, idx, idxinv, _x, _y, _z) JOIN(vgatheridx_,name) = VADDI32(JOIN(vgatheridx_,name), vgatherinc);
 543                                 D3Q19_LIST_WO_C
 544                         #undef X
 545
 546                         ppdf_C += offset_ppdf_C;
 547
 548                 }
 549                 else {
 550                         // Load new pointers to PDFs of local cell:
 551                         Assert(consecIndex < nConsecNodes);
 552
 553                         consecValue = consecNodes[consecIndex] - 1;
 554
 555                         adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 556                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vgatheridx_,name) = VLIU(&(ADJ_LIST(idxinv)));
 557                                 D3Q19_LIST_WO_C
 558                         #undef X
 559
 560                         ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
 561                         ++consecIndex;
 562                 }
 563
 564                 if (consecValue >= (VSIZE - 1)) {
 565                         // Vectorized part.
 566                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
 567                                 int const indexPrefetchL2 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L2;
 568                                 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
 569                                 if (indexPrefetchL2 < indexStop){
 570                                         // update pointers from adjacency list only if necessary
 571                                         if (consecValue >= (SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE + VSIZE - 1)) {
 572                                                 #define INCR_PTR(name)          (VADDI32(JOIN(vgatheridx_,name), VMULI32(vgatherinc, VSETI32(SOFTWARE_PREFETCH_LOOKAHEAD_L2))))
 573                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(INCR_PTR(name), (char const *) src, 8, _MM_HINT_T1);
 574                                                         D3Q19_LIST_WO_C
 575                                                 #undef X
 576                                                 #undef INCR_PTR
 577                                         }
 578                                         else {
 579                                                 adjListIndex = (indexPrefetchL2 - (indexPrefetchL2 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL2 % VSIZE);
 580                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T1);
 581                                                         D3Q19_LIST_WO_C
 582                                                 #undef X
 583                                         }
 584
 585                                         _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL2, D3Q19_C)]), _MM_HINT_T1);
 586                                 }
 587                         #endif
 588
 589                         #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
 590                                 int const indexPrefetchL1 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L1;
 591                                 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
 592                                 if (indexPrefetchL1 < indexStop){
 593                                         // update pointers from adjacency list only if necessary
 594                                         if (consecValue > (SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE + VSIZE - 1)) {
 595                                                 #define INCR_PTR(name)          (VADDI32(JOIN(vgatheridx_,name), VMULI32(vgatherinc, VSETI32(SOFTWARE_PREFETCH_LOOKAHEAD_L1))))
 596                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(INCR_PTR(name), (char const *) src, 8, _MM_HINT_T0);
 597                                                         D3Q19_LIST_WO_C
 598                                                 #undef X
 599                                                 #undef INCR_PTR
 600                                         }
 601                                         else {
 602                                                 adjListIndex = (indexPrefetchL1 - (indexPrefetchL1 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL1 % VSIZE);
 603                                                 #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T0);
 604                                                         D3Q19_LIST_WO_C
 605                                                 #undef X
 606                                         }
 607
 608                                         _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL1, D3Q19_C)]), _MM_HINT_T0);
 609                                 }
 610                         #endif
 611
 612                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(vpdf_,name) = VG32(JOIN(vgatheridx_,name), src, 8);
 613                                 D3Q19_LIST_WO_C
 614                         #undef X
 615
 616                         vpdf_C = VLDU(ppdf_C);
 617
 618                         //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
 619                         //      vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
 620                         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);
 621                         //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
 622                         //      vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
 623                         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);
 624                         //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
 625                         //      vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
 626                         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);
 627
 628                         //vdens = vpdf_C +
 629                         //        vpdf_N  + vpdf_E  + vpdf_S  + vpdf_W  +
 630                         //        vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
 631                         //        vpdf_T  + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
 632                         //        vpdf_B  + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
 633                         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)),
 634                                                                                 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));
 635
 636                         //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
 637                         vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
 638
 639                         //src[I(index, D3Q19_C)]  =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
 640                         VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
 641
 642                         // collide axis aligend pdfs vectorized
 643                         #define SCAT(offsets, vsrc) VS32(src, offsets, vsrc, 8)
 644
 645                         //vw_1_indep = vw_1 * vdir_indep_trm;
 646                         vw_1_indep = VMUL(vw_1,vdir_indep_trm);
 647
 648                         // collide axis aligend pdfs vectorized
 649                         #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
 650                                 vui = tmpVui; \
 651                                 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));\
 652                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
 653                                 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
 654                                 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 655
 656                         COLLIDE_AA_V(vuy, N, S)
 657                         COLLIDE_AA_V(vux, E, W)
 658                         COLLIDE_AA_V(vuz, T, B)
 659
 660                         #undef COLLIDE_AA_V
 661
 662                         //vw_2_indep = vw_2 * vdir_indep_trm;
 663                         vw_2_indep = VMUL(vw_2,vdir_indep_trm);
 664
 665                         // collide axis unaligned pdfs vectorized
 666                         #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
 667                                 vui = tmpVui; \
 668                                 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));\
 669                                 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
 670                                 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
 671                                 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
 672
 673                         COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
 674                         COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
 675                         COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
 676                         COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
 677                         COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
 678                         COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
 679
 680                         #undef COLLIDE_UA_V
 681                         #undef SCAT
 682
 683                         consecValue   -= (VSIZE - 1);
 684                         index         += (VSIZE - 1);
 685                         offset_ppdf_C  = VSIZE;
 686
 687                 }
 688                 else {
 689                         // Scalar part.
 690
 691                         adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
 692                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
 693                                 D3Q19_LIST_WO_C
 694                         #undef X
 695                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = *(JOIN(ppdf_,name));
 696                                 D3Q19_LIST_WO_C
 697                         #undef X
 698
 699                         pdf_C = *ppdf_C;
 700
 701                         ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 702                                 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 703                         uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 704                                 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 705                         uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 706                                 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 707
 708                         dens = pdf_C +
 709                                 pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 710                                 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 711                                 pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 712                                 pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 713
 714                         dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
 715
 716                         // direction: w_0
 717                         *ppdf_C = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
 718
 719                         // direction: w_1
 720                         w_1_indep = w_1 * dir_indep_trm;
 721
 722                         #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
 723                                 ui = tmpUi; \
 724                                 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
 725                                 oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
 726                                 *(JOIN(ppdf_,dir2))  = JOIN(pdf_,dir1) - evenPart - oddPart; \
 727                                 *(JOIN(ppdf_,dir1))  = JOIN(pdf_,dir2) - evenPart + oddPart;
 728
 729                         COLLIDE_AA_S(uy, N, S)
 730                         COLLIDE_AA_S(ux, E, W)
 731                         COLLIDE_AA_S(uz, T, B)
 732
 733                         #undef COLLIDE_AA_S
 734
 735                         // direction: w_2
 736                         w_2_indep = w_2 * dir_indep_trm;
 737
 738                         #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
 739                                 ui = tmpUi; \
 740                                 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
 741                                 oddPart  = omegaOdd  * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
 742                                 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
 743                                 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
 744
 745                         COLLIDE_UA_S((-ux + uy), NW, SE)
 746                         COLLIDE_UA_S(( ux + uy), NE, SW)
 747                         COLLIDE_UA_S((-ux + uz), TW, BE)
 748                         COLLIDE_UA_S(( ux + uz), TE, BW)
 749                         COLLIDE_UA_S((-uy + uz), TS, BN)
 750                         COLLIDE_UA_S(( uy + uz), TN, BS)
 751
 752                         #undef COLLIDE_UA_S
 753
 754                         offset_ppdf_C = 1;
 755                 }
 756
 757         } // loop over fluid nodes
 758
 759 #undef ADJ_LIST
 760 #undef I
 761 }