projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add citation information
[LbmBenchmarkKernelsPublic.git] / src / BenchKernelD3Q19ListAaPvGatherHybrid.c
CommitLineData
8cafd9ea
MW		 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 "BenchKernelD3Q19ListAaPvGatherHybridCommon.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 loop iterations ahead
61#endif
62
63#ifndef SOFTWARE_PREFETCH_LOOKAHEAD_L1
64#define SOFTWARE_PREFETCH_LOOKAHEAD_L1 (0) //prefetchting X loop iterations ahead
65#endif
66
67static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
68static void KernelOdd( LatticeDesc * ld, KernelData * kernelData, CaseData * cd);
69
70void FNAME(D3Q19ListAaPvGatherHybridKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
71{
72 Assert(ld != NULL);
73 Assert(kernelData != NULL);
74 Assert(cd != NULL);
75
76 Assert(cd->Omega > 0.0);
77 Assert(cd->Omega < 2.0);
78
79#if defined(VTK_OUTPUT) || defined(STATISTICS) || defined(VERIFICATION)
80 KernelData * kd = (KernelData *)kernelData;
81#endif
82 KernelDataList * kdl = KDL(kernelData);
83
84 int maxIterations = cd->MaxIterations;
85 int nFluid = kdl->nFluid;
86
87 printf("\n");
88#if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0) || (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
89 printf("# Software prefetching enabled:\n");
90 printf("# Gather/Scatter + Load/Store prefetch lookahead L2: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L2);
91 printf("# Gather/Scatter + Load/Store prefetch lookahead L1: \t%d\n", SOFTWARE_PREFETCH_LOOKAHEAD_L1);
92#else
93 printf("# Software prefetching disabled.\n");
94#endif
95 printf("\n");
96
97 int nThreads = 1;
98#ifdef _OPENMP
99 nThreads = omp_get_max_threads();
100#endif
101
102 int * threadStartIndices = (int *)malloc(sizeof(int) * (nThreads + 1));
103
104 for (int i = 0; i < nThreads; ++i) {
105 threadStartIndices[i] = i * (nFluid / nThreads) + MinI(i, nFluid % nThreads);
106 }
107
108 threadStartIndices[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
130 for(int iter = 0; iter < maxIterations; iter += 2) {
131
132 // even time step
133
134#ifdef _OPENMP
135#pragma omp parallel default(none) shared(ld, kernelData, cd, threadStartIndices)
136#endif
137 {
138 //LIKWID_MARKER_START("KernelEven");
139 KernelEven(ld, kernelData, cd, threadStartIndices);
140 //LIKWID_MARKER_STOP("KernelEven");
141 }
142
143
144#ifdef VERIFICATION
145 kdl->Iteration = iter;
146 kd->PdfsActive = kd->Pdfs[0];
147 KernelAddBodyForce(kd, ld, cd);
148#endif
149
150 // odd time step
151
152#ifdef _OPENMP
153#pragma omp parallel default(none) shared(ld, kernelData, cd, threadStartIndices)
154#endif
155 {
156 //LIKWID_MARKER_START("KernelOdd");
157 KernelOdd(ld, kernelData, cd);
158 //LIKWID_MARKER_STOP("KernelOdd");
159 }
160
161
162#ifdef VERIFICATION
163 kdl->Iteration = iter + 1;
164 kd->PdfsActive = kd->Pdfs[0];
165 KernelAddBodyForce(kd, ld, cd);
166#endif
167
168#ifdef VTK_OUTPUT
169 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
170 kdl->Iteration = iter + 1;
171 kd->PdfsActive = kd->Pdfs[0];
172 VtkWrite(ld, kd, cd, iter);
173 }
174#endif
175
176#ifdef STATISTICS
177 kdl->Iteration = iter + 1;
178 kd->PdfsActive = kd->Pdfs[0];
179 KernelStatistics(kd, ld, cd, iter);
180#endif
181
182 } // for (int iter = 0; ...
183 LIKWID_MARKER_STOP("OuterLoop");
184
185 X_KERNEL_END(kernelData);
186
187#ifdef VTK_OUTPUT
188 if (cd->VtkOutput) {
189 kd->PdfsActive = kd->Pdfs[0];
190 VtkWrite(ld, kd, cd, maxIterations);
191 }
192#endif
193
194#ifdef STATISTICS
195 kd->PdfsActive = kd->Pdfs[0];
196 KernelStatistics(kd, ld, cd, maxIterations);
197#endif
198
199 LIKWID_MARKER_CLOSE;
200 free(threadStartIndices);
201
202 return;
203}
204
205static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices)
206{
207 Assert(ld != NULL);
208 Assert(kernelData != NULL);
209 Assert(cd != NULL);
210
211 Assert(cd->Omega > 0.0);
212 Assert(cd->Omega < 2.0);
213
214 KernelData * kd = (KernelData *)kernelData;
215 KernelDataList * kdl = KDL(kernelData);
216
217 PdfT omega = cd->Omega;
218 PdfT omegaEven = omega;
219
220 PdfT magicParam = 1.0 / 12.0;
221 PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
222
223 PdfT evenPart = 0.0;
224 PdfT oddPart = 0.0;
225 PdfT dir_indep_trm = 0.0;
226
227 const PdfT w_0 = 1.0 / 3.0;
228 const PdfT w_1 = 1.0 / 18.0;
229 const PdfT w_2 = 1.0 / 36.0;
230
231 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;
232 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;
233
234 PdfT ux, uy, uz, ui;
235 PdfT dens;
236
237 VPDFT VONE_HALF = VSET(0.5);
238 VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
239
240 VPDFT vw_1_indep, vw_2_indep;
241 VPDFT vw_0 = VSET(w_0);
242 VPDFT vw_1 = VSET(w_1);
243 VPDFT vw_2 = VSET(w_2);
244
245 VPDFT vw_1_x3 = VSET(w_1_x3);
246 VPDFT vw_2_x3 = VSET(w_2_x3);
247 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
248 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
249
250 VPDFT vui, vux, vuy, vuz, vdens;
251
252 VPDFT vevenPart, voddPart, vdir_indep_trm;
253
254 VPDFT vomegaEven = VSET(omegaEven);
255 VPDFT vomegaOdd = VSET(omegaOdd);
256
257 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
258 #define X(name, idx, idxinv, x, y, z) \
259 PdfT JOIN(pdf_,name); \
260 PdfT * JOIN(ppdf_,name); \
261 VPDFT JOIN(vpdf_,name);
262 D3Q19_LIST
263 #undef X
264
265 PdfT * src = kd->Pdfs[0];
266
267 int nCells = kdl->nCells;
268
269 int threadId = 0;
270#ifdef _OPENMP
271 threadId = omp_get_thread_num();
272#endif
273
274 int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId];
275 int nFluidVec = nFluidThread - (nFluidThread % VSIZE);
276
277 int indexStartVec = threadIndices[threadId];
278 int indexStopVec = threadIndices[threadId] + nFluidVec;
279 int indexStop = threadIndices[threadId] + nFluidThread;
280
281 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
282
283 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[I(indexStartVec, idx)]);
284 D3Q19_LIST
285 #undef X
286
287 for (int index = indexStartVec; index < indexStopVec; index += VSIZE) {
288
289 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
290 #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE), _MM_HINT_T1);
291 D3Q19_LIST
292 #undef X
293 #endif
294
295 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
296 #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE), _MM_HINT_T0);
297 D3Q19_LIST
298 #undef X
299 #endif
300
301 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
302 D3Q19_LIST
303 #undef X
304
305 //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
306 // vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
307 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);
308 //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
309 // vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
310 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);
311 //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
312 // vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
313 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);
314
315 //vdens = vpdf_C +
316 // vpdf_N + vpdf_E + vpdf_S + vpdf_W +
317 // vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
318 // vpdf_T + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
319 // vpdf_B + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
320 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));
321
322 //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
323 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
324
325 //src[I(index, D3Q19_C)] =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
326 VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
327
328 //vw_1_indep = vw_1 * vdir_indep_trm;
329 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
330
331 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
332 vui = tmpVui; \
333 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));\
334 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
335 VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
336 VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
337
338 COLLIDE_AA_V(vuy, N, S)
339 COLLIDE_AA_V(vux, E, W)
340 COLLIDE_AA_V(vuz, T, B)
341
342 #undef COLLIDE_AA_V
343
344 //vw_2_indep = vw_2 * vdir_indep_trm;
345 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
346
347 // collide axis unaligned pdfs vectorized
348 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
349 vui = tmpVui; \
350 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));\
351 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
352 VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
353 VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
354
355 COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
356 COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
357 COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
358 COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
359 COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
360 COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
361
362 #undef COLLIDE_UA_V
363
364 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) +=VSIZE;
365 D3Q19_LIST
366 #undef X
367 } // loop over fluid nodes
368
369 for (int index = indexStopVec; index < indexStop; ++index) {
370
371 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *(JOIN(ppdf_,name));
372 D3Q19_LIST
373 #undef X
374
375 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
376 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
377 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
378 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
379 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
380 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
381
382 dens = pdf_C +
383 pdf_N + pdf_E + pdf_S + pdf_W +
384 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
385 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
386 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
387
388 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
389
390 // direction: w_0
391 *ppdf_C = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
392
393 // direction: w_1
394 w_1_indep = w_1*dir_indep_trm;
395
396 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
397 ui = tmpUi; \
398 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
399 oddPart = omegaOdd * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
400 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
401 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
402
403 COLLIDE_AA_S(uy, N, S)
404 COLLIDE_AA_S(ux, E, W)
405 COLLIDE_AA_S(uz, T, B)
406
407 #undef COLLIDE_AA_S
408
409 // direction: w_2
410 w_2_indep = w_2*dir_indep_trm;
411
412 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
413 ui = tmpUi; \
414 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
415 oddPart = omegaOdd * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
416 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
417 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
418
419 COLLIDE_UA_S((-ux + uy), NW, SE)
420 COLLIDE_UA_S(( ux + uy), NE, SW)
421 COLLIDE_UA_S((-ux + uz), TW, BE)
422 COLLIDE_UA_S(( ux + uz), TE, BW)
423 COLLIDE_UA_S((-uy + uz), TS, BN)
424 COLLIDE_UA_S(( uy + uz), TN, BS)
425
426 #undef COLLIDE_UA_S
427
428 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name)++;
429 D3Q19_LIST
430 #undef X
431 } // loop over fluid nodes
432
433 #undef I
434
435 return;
436}
437
438
439static void KernelOdd(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
440{
441
442 Assert(ld != NULL);
443 Assert(kernelData != NULL);
444 Assert(cd != NULL);
445
446 Assert(cd->Omega > 0.0);
447 Assert(cd->Omega < 2.0);
448
449 KernelData * kd = (KernelData *)kernelData;
450 KernelDataList * kdl = KDL(kernelData);
451 KernelDataListRia * kdlr = KDLR(kernelData);
452
453 PdfT omega = cd->Omega;
454 PdfT omegaEven = omega;
455
456 PdfT magicParam = 1.0 / 12.0;
457 PdfT omegaOdd = 1.0 / (0.5 + magicParam / (1.0 / omega - 0.5));
458
459 PdfT evenPart = 0.0;
460 PdfT oddPart = 0.0;
461 PdfT dir_indep_trm = 0.0;
462
463 const PdfT w_0 = 1.0 / 3.0;
464 const PdfT w_1 = 1.0 / 18.0;
465 const PdfT w_2 = 1.0 / 36.0;
466
467 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;
468 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;
469
470 PdfT ux, uy, uz, ui;
471 PdfT dens;
472
473 VPDFT VONE_HALF = VSET(0.5);
474 VPDFT VTHREE_HALF = VSET(3.0 / 2.0);
475
476 VPDFT vw_1_indep, vw_2_indep;
477 VPDFT vw_0 = VSET(w_0);
478 VPDFT vw_1 = VSET(w_1);
479 VPDFT vw_2 = VSET(w_2);
480
481 VPDFT vw_1_x3 = VSET(w_1_x3);
482 VPDFT vw_2_x3 = VSET(w_2_x3);
483 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
484 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
485
486 VPDFT vux, vuy, vuz, vui;
487 VPDFT vdens;
488
489 VPDFT vevenPart, voddPart, vdir_indep_trm;
490
491 VPDFT vomegaEven = VSET(omegaEven);
492 VPDFT vomegaOdd = VSET(omegaOdd);
493
494
495 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
496 #define X(name, idx, idxinv, x, y, z) \
497 PdfT JOIN(pdf_,name); \
498 PdfT * JOIN(ppdf_,name) = NULL; \
499 VPDFT JOIN(vpdf_,name);
500 D3Q19_LIST
501 #undef X
502 #define X(name, idx, idxinv, x, y, z) __m256i JOIN(vgatheridx_,name);
503 D3Q19_LIST_WO_C
504 #undef X
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 * oddKernelThreadStartIndices = kdlr->oddKernelThreadStartIndices;
514 int * loopStartIndices = kdlr->loopStartIndices;
515
516 int threadId = 0;
517
518#ifdef _OPENMP
519 threadId = omp_get_thread_num();
520#endif
521
522 int threadStartIndex = oddKernelThreadStartIndices[threadId + 0];
523 int threadStopIndex = oddKernelThreadStartIndices[threadId + 1];
524
525 int indexStart = loopStartIndices[threadStartIndex];
526 int indexStop = loopStartIndices[threadStopIndex ];
527
528 int consecStartIndex = threadStartIndex + 1;
529 int consecStopIndex = threadStopIndex - 1;
530
531 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
532 #define ADJ_LIST(dir) adjList[adjListIndex + (dir * VSIZE)]
533
534 // scalar peel loop
535 #pragma novector
536 #pragma loop_count max=7
537 for (int index = indexStart; index < loopStartIndices[consecStartIndex]; ++index) {
538
539 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
540 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
541 D3Q19_LIST_WO_C
542 #undef X
543 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *(JOIN(ppdf_,name));
544 D3Q19_LIST_WO_C
545 #undef X
546
547 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
548 pdf_C = *ppdf_C;
549
550 // macroscopic values
551 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
552 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
553 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
554 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
555 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
556 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
557
558 dens = pdf_C +
559 pdf_N + pdf_E + pdf_S + pdf_W +
560 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
561 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
562 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
563
564 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
565
566 // collide direction: w_0
567 *ppdf_C = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
568
569 // collide direction: w_1
570 w_1_indep = w_1 * dir_indep_trm;
571
572 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
573 ui = tmpUi; \
574 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
575 oddPart = omegaOdd * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
576 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
577 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
578
579 COLLIDE_AA_S(uy, N, S)
580 COLLIDE_AA_S(ux, E, W)
581 COLLIDE_AA_S(uz, T, B)
582
583 #undef COLLIDE_AA_S
584
585 // collide direction: w_2
586 w_2_indep = w_2 * dir_indep_trm;
587
588 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
589 ui = tmpUi; \
590 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
591 oddPart = omegaOdd * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
592 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
593 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
594
595 COLLIDE_UA_S((-ux + uy), NW, SE)
596 COLLIDE_UA_S(( ux + uy), NE, SW)
597 COLLIDE_UA_S((-ux + uz), TW, BE)
598 COLLIDE_UA_S(( ux + uz), TE, BW)
599 COLLIDE_UA_S((-uy + uz), TS, BN)
600 COLLIDE_UA_S(( uy + uz), TN, BS)
601
602 #undef COLLIDE_UA_S
603
604 } // scalar peel loop
605
606 // vectorized loop
607 for (int consecIndex = consecStartIndex; consecIndex < consecStopIndex; /* increment in loop body */) {
608
609 int index = loopStartIndices[consecIndex];
610 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
611 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
612 D3Q19_LIST_WO_C
613 #undef X
614 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
615
616 // loop for consecutive access pattern -> load & store
617 ++consecIndex;
618
619 for (; index < loopStartIndices[consecIndex]; index+=VSIZE){
620
621 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
622 #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L2 * VSIZE), _MM_HINT_T1);
623 D3Q19_LIST
624 #undef X
625 #endif
626 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
627 #define X(name, idx, idxinv, _x, _y, _z) _mm_prefetch((char const *)(JOIN(ppdf_,name) + SOFTWARE_PREFETCH_LOOKAHEAD_L1 * VSIZE), _MM_HINT_T0);
628 D3Q19_LIST
629 #undef X
630 #endif
631
632 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(JOIN(ppdf_,name));
633 D3Q19_LIST
634 #undef X
635
636 // macroscopic values
637 //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
638 // vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
639 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);
640 //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
641 // vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
642 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);
643 //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
644 // vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
645 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);
646
647 //vdens = vpdf_C +
648 // vpdf_N + vpdf_E + vpdf_S + vpdf_W +
649 // vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
650 // vpdf_T + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
651 // vpdf_B + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
652 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)),
653 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));
654
655 //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
656 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
657
658 // collide direction: w_0
659 //src[I(index, D3Q19_C)] =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
660 VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
661
662 // collide direction: w_1
663 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
664
665 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
666 vui = tmpVui; \
667 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));\
668 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
669 VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
670 VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
671
672 COLLIDE_AA_V(vuy, N, S)
673 COLLIDE_AA_V(vux, E, W)
674 COLLIDE_AA_V(vuz, T, B)
675
676 #undef COLLIDE_AA_V
677
678 // collide direction: w_2
679 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
680
681 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
682 vui = tmpVui; \
683 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));\
684 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
685 VSTU(JOIN(ppdf_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
686 VSTU(JOIN(ppdf_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
687
688 COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
689 COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
690 COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
691 COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
692 COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
693 COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
694
695 #undef COLLIDE_UA_V
696
697 //increment pointer
698 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) += VSIZE;
699 D3Q19_LIST
700 #undef X
701 }
702
703 // loop for irregular access pattern -> gather & scatter
704 index = loopStartIndices[consecIndex];
705 ++consecIndex;
706 for (; index < loopStartIndices[consecIndex]; index+=VSIZE){
707
708
709 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L2 > 0)
710 int const indexPrefetchL2 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L2;
711 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
712 if (indexPrefetchL2 < indexStop){
713 adjListIndex = (indexPrefetchL2 - (indexPrefetchL2 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL2 % VSIZE);
714 #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T1);
715 D3Q19_LIST_WO_C
716 #undef X
717
718 _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL2, D3Q19_C)]), _MM_HINT_T1);
719 }
720 #endif
721
722 #if (SOFTWARE_PREFETCH_LOOKAHEAD_L1 > 0)
723 int const indexPrefetchL1 = index + VSIZE * SOFTWARE_PREFETCH_LOOKAHEAD_L1;
724 // make sure that adjList access is never out of bounds since it is an actual memory access and no prefetch
725 if (indexPrefetchL1 < indexStop){
726 adjListIndex = (indexPrefetchL1 - (indexPrefetchL1 % VSIZE)) * N_D3Q19_IDX + (indexPrefetchL1 % VSIZE);
727 #define X(name, idx, idxinv, _x, _y, _z) VPG32(VLIU(&ADJ_LIST(idxinv)), (char const *) src, 8, _MM_HINT_T0);
728 D3Q19_LIST_WO_C
729 #undef X
730
731 _mm_prefetch((char const *) &(src[P_INDEX_3(nCells, indexPrefetchL1, D3Q19_C)]), _MM_HINT_T0);
732 }
733 #endif
734
735 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
736 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vgatheridx_,name) = VLIU(&(ADJ_LIST(idxinv)));
737 D3Q19_LIST_WO_C
738 #undef X
739
740 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VG32(JOIN(vgatheridx_,name), src, 8);
741 D3Q19_LIST_WO_C
742 #undef X
743
744 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
745 vpdf_C = VLDU(ppdf_C);
746
747 // macroscopic values
748 //vux = vpdf_E + vpdf_NE + vpdf_SE + vpdf_TE + vpdf_BE -
749 // vpdf_W - vpdf_NW - vpdf_SW - vpdf_TW - vpdf_BW;
750 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);
751 //vuy = vpdf_N + vpdf_NE + vpdf_NW + vpdf_TN + vpdf_BN -
752 // vpdf_S - vpdf_SE - vpdf_SW - vpdf_TS - vpdf_BS;
753 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);
754 //vuz = vpdf_T + vpdf_TE + vpdf_TW + vpdf_TN + vpdf_TS -
755 // vpdf_B - vpdf_BE - vpdf_BW - vpdf_BN - vpdf_BS;
756 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);
757
758 //vdens = vpdf_C +
759 // vpdf_N + vpdf_E + vpdf_S + vpdf_W +
760 // vpdf_NE + vpdf_SE + vpdf_SW + vpdf_NW +
761 // vpdf_T + vpdf_TN + vpdf_TE + vpdf_TS + vpdf_TW +
762 // vpdf_B + vpdf_BN + vpdf_BE + vpdf_BS + vpdf_BW;
763 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)),
764 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));
765
766 //vdir_indep_trm = vdens - (vux * vux + vuy * vuy + vuz * vuz) * VTHREE_HALF;
767 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
768
769 // collide direction: w_0
770 //src[I(index, D3Q19_C)] =[UA] vpdf_C - vomegaEven * (vpdf_C - vw_0 * vdir_indep_trm);
771 VSTU(ppdf_C,VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
772
773
774 // collide axis aligend pdfs vectorized
775 #define SCAT(offsets, vsrc) VS32(src, offsets, vsrc, 8)
776
777 // collide direction: w_1
778 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
779
780 #define COLLIDE_AA_V(tmpVui, dir1, dir2) \
781 vui = tmpVui; \
782 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));\
783 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_1_x3)));\
784 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart));\
785 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
786
787 COLLIDE_AA_V(vuy, N, S)
788 COLLIDE_AA_V(vux, E, W)
789 COLLIDE_AA_V(vuz, T, B)
790
791 #undef COLLIDE_AA_V
792
793 // collide direction: w_2
794 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
795
796 #define COLLIDE_UA_V(tmpVui, dir1, dir2) \
797 vui = tmpVui; \
798 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));\
799 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(JOIN(vpdf_,dir1),JOIN(vpdf_,dir2))),VMUL(vui,vw_2_x3)));\
800 SCAT(JOIN(vgatheridx_,dir2),VSUB(VSUB(JOIN(vpdf_,dir1),vevenPart),voddPart)); \
801 SCAT(JOIN(vgatheridx_,dir1),VADD(VSUB(JOIN(vpdf_,dir2),vevenPart),voddPart));
802
803 COLLIDE_UA_V(VSUB(vuy,vux), NW, SE)
804 COLLIDE_UA_V(VADD(vux,vuy), NE, SW)
805 COLLIDE_UA_V(VSUB(vuz,vux), TW, BE)
806 COLLIDE_UA_V(VADD(vux,vuz), TE, BW)
807 COLLIDE_UA_V(VSUB(vuz,vuy), TS, BN)
808 COLLIDE_UA_V(VADD(vuy,vuz), TN, BS)
809
810 #undef COLLIDE_UA_V
811 #undef SCAT
812 } // gather & scatter loop
813
814 } // Vectorized Loop
815
816 // Scalar remainder loop
817 #pragma novector
818 #pragma loop_count max=7
819 for (int index = loopStartIndices[consecStopIndex]; index < indexStop; ++index) {
820
821 adjListIndex = (index - (index % VSIZE)) * N_D3Q19_IDX + (index % VSIZE);
822 #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[ADJ_LIST(idxinv)]);
823 D3Q19_LIST_WO_C
824 #undef X
825 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *(JOIN(ppdf_,name));
826 D3Q19_LIST_WO_C
827 #undef X
828
829 ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]);
830 pdf_C = *ppdf_C;
831
832 // macroscopic values
833 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
834 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
835 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
836 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
837 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
838 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
839
840 dens = pdf_C +
841 pdf_N + pdf_E + pdf_S + pdf_W +
842 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
843 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
844 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
845
846 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
847
848 // collide direction: w_0
849 *ppdf_C = pdf_C - omegaEven * (pdf_C - w_0 * dir_indep_trm);
850
851 // collide direction: w_1
852 w_1_indep = w_1 * dir_indep_trm;
853
854 #define COLLIDE_AA_S(tmpUi, dir1, dir2) \
855 ui = tmpUi; \
856 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_1_nine_half - w_1_indep); \
857 oddPart = omegaOdd * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_1_x3); \
858 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
859 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
860
861 COLLIDE_AA_S(uy, N, S)
862 COLLIDE_AA_S(ux, E, W)
863 COLLIDE_AA_S(uz, T, B)
864
865 #undef COLLIDE_AA_S
866
867 // collide direction: w_2
868 w_2_indep = w_2 * dir_indep_trm;
869
870 #define COLLIDE_UA_S(tmpUi, dir1, dir2) \
871 ui = tmpUi; \
872 evenPart = omegaEven * (0.5 * (JOIN(pdf_,dir1) + JOIN(pdf_,dir2)) - ui * ui * w_2_nine_half - w_2_indep); \
873 oddPart = omegaOdd * (0.5 * (JOIN(pdf_,dir1) - JOIN(pdf_,dir2)) - ui * w_2_x3); \
874 *(JOIN(ppdf_,dir2)) = JOIN(pdf_,dir1) - evenPart - oddPart; \
875 *(JOIN(ppdf_,dir1)) = JOIN(pdf_,dir2) - evenPart + oddPart;
876
877 COLLIDE_UA_S((-ux + uy), NW, SE)
878 COLLIDE_UA_S(( ux + uy), NE, SW)
879 COLLIDE_UA_S((-ux + uz), TW, BE)
880 COLLIDE_UA_S(( ux + uz), TE, BW)
881 COLLIDE_UA_S((-uy + uz), TS, BN)
882 COLLIDE_UA_S(( uy + uz), TN, BS)
883
884 #undef COLLIDE_UA_S
885
886 } // remainder loop over fluid nodes
887
888#undef SOFTWARE_PREFETCH_LOAD_STORE
889#undef SOFTWARE_PREFETCH_LOAD_STORE_HINT
890#undef SOFTWARE_PREFETCH_LOAD_STORE_LOOKAHEAD
891#undef SOFTWARE_PREFETCH_GATHER_SCATTER
892#undef SOFTWARE_PREFETCH_GATHER_SCATTER_HINT
893#undef SOFTWARE_PREFETCH_GATHER_SCATTER_LOOKAHEAD
894#undef ADJ_LIST
895#undef I
896}
LBM benchmark kernels
This page took 0.111911 seconds and 5 git commands to generate.