projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add citation information
[LbmBenchmarkKernelsPublic.git] / src / BenchKernelD3Q19ListAaPvGather.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 "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
67static void KernelEven(LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
68static void KernelOdd( LatticeDesc * ld, KernelData * kernelData, CaseData * cd, int * threadIndices);
69
70void 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
200static 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
433static 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}
LBM benchmark kernels
This page took 0.106683 seconds and 5 git commands to generate.