1 // --------------------------------------------------------------------------
4 // Markus Wittmann, 2016-2017
5 // RRZE, University of Erlangen-Nuremberg, Germany
6 // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
9 // LSS, University of Erlangen-Nuremberg, Germany
11 // Michael Hussnaetter, 2017-2018
12 // University of Erlangen-Nuremberg, Germany
13 // michael.hussnaetter -at- fau.de
15 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
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.
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.
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/>.
30 // --------------------------------------------------------------------------
31 #include "BenchKernelD3Q19ListAaPvGatherAoSoACommon.h"
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)
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
63 #ifndef SOFTWARE_PREFETCH_LOOKAHEAD_L1
64 #define SOFTWARE_PREFETCH_LOOKAHEAD_L1 (0) //prefetchting X SIMD widths ahead
67 static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
68 static void KernelOdd( LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
70 void FNAME(D3Q19ListAaPvGatherAoSoAKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
74 Assert(kernelData != NULL);
77 Assert(cd->Omega > 0.0);
78 Assert(cd->Omega < 2.0);
80 #if defined(VTK_OUTPUT) || defined(STATISTICS) || defined(VERIFICATION)
81 KernelData * kd = (KernelData *)kernelData;
83 KernelDataList * kdl = KDL(kernelData);
85 int maxIterations = cd->MaxIterations;
86 int nFluid = kdl->nFluid;
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);
94 printf("# Software prefetching disabled.\n");
100 nThreads = omp_get_max_threads();
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);
108 threadIndices[nThreads] = nFluid;
112 kd->PdfsActive = kd->Pdfs[0];
113 VtkWrite(ld, kd, cd, -1);
118 kd->PdfsActive = kd->Pdfs[0];
119 KernelStatistics(kd, ld, cd, 0);
124 X_KERNEL_START(kernelData);
126 // TODO: outer openmp parallel
128 LIKWID_MARKER_START("OuterLoop");
129 for(int iter = 0; iter < maxIterations; iter += 2) {
134 #pragma omp parallel default(none) shared(ld, kernelData, cd, threadIndices)
137 KernelEven(ld, kernelData, cd, threadIndices);
142 kdl->Iteration = iter;
143 kd->PdfsActive = kd->Pdfs[0];
144 KernelAddBodyForce(kd, ld, cd);
150 #pragma omp parallel default(none) shared(ld, kernelData, cd, threadIndices)
153 KernelOdd(ld, kernelData, cd, threadIndices);
158 kdl->Iteration = iter + 1;
159 kd->PdfsActive = kd->Pdfs[0];
160 KernelAddBodyForce(kd, ld, cd);
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);
172 kdl->Iteration = iter + 1;
173 kd->PdfsActive = kd->Pdfs[0];
174 KernelStatistics(kd, ld, cd, iter);
177 } // for (int iter = 0; ...
178 LIKWID_MARKER_STOP("OuterLoop");
180 X_KERNEL_END(kernelData);
184 kd->PdfsActive = kd->Pdfs[0];
185 VtkWrite(ld, kd, cd, maxIterations);
190 kd->PdfsActive = kd->Pdfs[0];
191 KernelStatistics(kd, ld, cd, maxIterations);
200 static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
203 Assert(kernelData != NULL);
206 Assert(cd->Omega > 0.0);
207 Assert(cd->Omega < 2.0);
209 KernelData * kd = (KernelData *)kernelData;
210 KernelDataList * kdl = KDL(kernelData);
212 PdfT omega = cd->Omega;
213 PdfT omegaEven = omega;
215 PdfT magicParam = 1.0 / 12.0;
216 PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
220 PdfT dir_indep_trm = 0.0;
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;
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;
232 VPDFT VONE_HALF = VSET(0.5);
233 VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
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);
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);
245 VPDFT vui, vux, vuy, vuz, vdens;
247 VPDFT vevenPart, voddPart, vdir_indep_trm;
249 VPDFT vomegaEven = VSET(omegaEven);
250 VPDFT vomegaOdd = VSET(omegaOdd);
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);
260 PdfT * src = kd->Pdfs[0];
262 int nCells = kdl->nCells;
266 threadId = omp_get_thread_num();
270 int indexStart = threadIndices[threadId];
271 int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
272 int indexStop = indexStart + nFluidThread;
274 int indexStartVec = ((indexStart + VSIZE - 1) / VSIZE) * VSIZE;
275 int nFluidVec = (indexStop / VSIZE) * VSIZE - indexStartVec;
276 int indexStopVec = indexStartVec + nFluidVec;
278 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
280 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[I(indexStart, idx)]);
284 for (int index = indexStart; index < indexStartVec; ++index) {
286 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *(JOIN(ppdf_,name));
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;
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;
303 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
306 *ppdf_C = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
309 w_1_indep = w_1*dir_indep_trm;
311 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
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;
318 COLLIDE_AA_S(uy, N, S)
319 COLLIDE_AA_S(ux, E, W)
320 COLLIDE_AA_S(uz, T, B)
325 w_2_indep = w_2*dir_indep_trm;
327 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
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;
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)
343 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name)++;
348 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[I(indexStartVec, idx)]);
352 for (int index = indexStartVec; index < indexStopVec; index += VSIZE) {
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);
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);
366 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
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);
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));
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));
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)))));
393 //vw_1_indep = vw_1 * vdir_indep_trm;
394 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
396 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
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));
403 COLLIDE_AA_V(vuy, N, S)
404 COLLIDE_AA_V(vux, E, W)
405 COLLIDE_AA_V(vuz, T, B)
409 //vw_2_indep = vw_2 * vdir_indep_trm;
410 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
412 // collide axis unaligned pdfs vectorized
413 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
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));
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)
429 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) +=(VSIZE * N_D3Q19);
432 } // loop over fluid nodes
434 for (int index = indexStopVec; index < indexStop; ++index) {
436 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *(JOIN(ppdf_,name));
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;
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;
453 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
456 *ppdf_C = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
459 w_1_indep = w_1*dir_indep_trm;
461 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
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;
468 COLLIDE_AA_S(uy, N, S)
469 COLLIDE_AA_S(ux, E, W)
470 COLLIDE_AA_S(uz, T, B)
475 w_2_indep = w_2*dir_indep_trm;
477 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
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;
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)
493 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name)++;
496 } // loop over fluid nodes
503 static void KernelOdd(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
507 Assert(kernelData != NULL);
510 Assert(cd->Omega > 0.0);
511 Assert(cd->Omega < 2.0);
513 KernelData * kd = (KernelData *)kernelData;
514 KernelDataList * kdl = KDL(kernelData);
515 KernelDataListRia * kdlr = KDLR(kernelData);
516 PdfT omega = cd->Omega;
517 PdfT omegaEven = omega;
519 PdfT magicParam = 1.0 / 12.0;
520 PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
524 PdfT dir_indep_trm = 0.0;
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;
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;
536 VPDFT VONE_HALF = VSET(0.5);
537 VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
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);
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);
549 VPDFT vux, vuy, vuz, vui;
552 VPDFT vevenPart, voddPart, vdir_indep_trm;
554 VPDFT vomegaEven = VSET(omegaEven);
555 VPDFT vomegaOdd = VSET(omegaOdd);
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);
565 #define X(name, idx, idxinv, x, y, z) \
566 __m256i JOIN(vgatheridx_,name) = _mm256_set1_epi32(0);
570 __m256i vgatherinc = VSETI32(VSIZE * N_D3Q19);
572 uint32_t * consecNodes = kdlr->ConsecNodes;
573 uint32_t consecIndex = 0;
574 uint32_t consecValue = 0;
576 PdfT * src = kd->Pdfs[0];
578 int nCells = kdl->nCells;
581 uint32_t * adjList = kdl->AdjList;
586 threadId = omp_get_thread_num();
588 consecIndex = kdlr->ConsecThreadIndices[threadId];
591 int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
593 int indexStart = threadIndices[threadId];
594 int indexStop = threadIndices[threadId] + nFluidThread;
596 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
597 #define ADJ_LIST(dir) adjList[adjListIndex + (dir * VSIZE)]
599 int offset_ppdf_C = -1; //dummy init to detect errors
601 for (int index = indexStart; index < indexStop; index += 1) {
603 if (consecValue > 0) {
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
612 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vgatheridx_,name) = VADDI32(JOIN(vgatheridx_,name), vgatherinc);
616 //printVector(vgatheridx_N);
618 ppdf_C += offset_ppdf_C;
622 // Load new pointers to PDFs of local cell:
623 Assert(consecIndex < nConsecNodes);
625 consecValue = consecNodes[consecIndex] - 1;
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)));
632 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
636 if (consecValue >= (VSIZE - 1)) {
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);
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);
658 _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL2, D3Q19_C)]), _MM_HINT_T1);
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);
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);
681 _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL1, D3Q19_C)]), _MM_HINT_T0);
685 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VG32(JOIN(vgatheridx_,name), src, 8);
689 vpdf_C = VLDU(ppdf_C);
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);
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));
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));
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)))));
715 // collide axis aligend pdfs vectorized
716 #define SCAT(offsets, vsrc) VS32(src, offsets, vsrc, 8)
718 //vw_1_indep = vw_1 * vdir_indep_trm;
719 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
721 // collide axis aligend pdfs vectorized
722 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
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));
729 COLLIDE_AA_V(vuy, N, S)
730 COLLIDE_AA_V(vux, E, W)
731 COLLIDE_AA_V(vuz, T, B)
735 //vw_2_indep = vw_2 * vdir_indep_trm;
736 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
738 // collide axis unaligned pdfs vectorized
739 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
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));
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)
756 consecValue -= (VSIZE - 1);
757 index += (VSIZE - 1);
758 offset_ppdf_C = VSIZE * N_D3Q19;
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)]);
768 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *(JOIN(ppdf_,name));
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;
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;
787 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
790 *ppdf_C = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
793 w_1_indep = w_1 * dir_indep_trm;
795 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
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;
802 COLLIDE_AA_S(uy, N, S)
803 COLLIDE_AA_S(ux, E, W)
804 COLLIDE_AA_S(uz, T, B)
809 w_2_indep = w_2 * dir_indep_trm;
811 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
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;
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)
830 } // loop over fluid nodes