projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add single precision, add aa-vec-sl-soa kernel, updated doc
[LbmBenchmarkKernelsPublic.git] / src / BenchKernelD3Q19AaVec.c
CommitLineData
e3f82424
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// This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
12//
13// LbmBenchKernels is free software: you can redistribute it and/or modify
14// it under the terms of the GNU General Public License as published by
15// the Free Software Foundation, either version 3 of the License, or
16// (at your option) any later version.
17//
18// LbmBenchKernels is distributed in the hope that it will be useful,
19// but WITHOUT ANY WARRANTY; without even the implied warranty of
20// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
21// GNU General Public License for more details.
22//
23// You should have received a copy of the GNU General Public License
24// along with LbmBenchKernels. If not, see <http://www.gnu.org/licenses/>.
25//
26// --------------------------------------------------------------------------
27#include "BenchKernelD3Q19AaVecCommon.h"
28
29#include "Memory.h"
30#include "Vtk.h"
31#include "LikwidIf.h"
32#include "Vector.h"
33#include "Vector.h"
34
35#include <inttypes.h>
36#include <math.h>
37
38#ifdef _OPENMP
39 #include <omp.h>
40#endif
41
42static void KernelEven(LatticeDesc * ld, KernelData * kd, CaseData * cd);
43static void KernelOdd( LatticeDesc * ld, KernelData * kd, CaseData * cd);
44
45#if 0 // {{{
46void DumpPdfs(LatticeDesc * ld, KernelData * kd, int zStart, int zStop, int iter, const char * prefix)
47{
48 return;
49
50 int * lDims = ld->Dims;
51 int * gDims = kd->GlobalDims;
52
53 int nX = gDims[0];
54 int nY = gDims[1];
55 int nZ = gDims[2];
56
57 PdfT pdfs[N_D3Q19];
58
59 int localZStart = zStart;
60 int localZStop = zStop;
61
62 if (localZStart == -1) localZStart = 0;
63 if (localZStop == -1) localZStop = gDims[2] - 1;
64
65 printf("D iter: %d\n", iter);
66
67 for (int dir = 0; dir < 19; ++dir) {
68 for (int z = localZStop; z >= localZStart; --z) {
69 printf("D [%2d][%2d][%s] plane % 2d\n", iter, dir, prefix, z);
70
71 for(int y = 0; y < nY; ++y) {
72 printf("D [%2d][%2d][%s] %2d ", iter, dir, prefix, y);
73
74 for(int x = 0; x < nX; ++x) {
75
76 if (1) { // ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
77
78 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
79 pdfs[dir] = kd->PdfsActive[I(x, y, z, dir)];
80 #undef I
81// kd->GetNode(kd, x, y, z, pdfs);
82 }
83 else {
0fde6e45 84 pdfs[dir] = -F(1.0);
e3f82424
MW		 85 }
86
87 printf("%.16e ", pdfs[dir]);
88 }
89
90 printf("\n");
91 }
92 }
93 }
94}
95#endif // }}}
96
97void FNAME(D3Q19AaVecKernel)(LatticeDesc * ld, KernelData * kd, CaseData * cd)
98{
99 Assert(ld != NULL);
100 Assert(kd != NULL);
101 Assert(cd != NULL);
102
0fde6e45
MW		 103 Assert(cd->Omega > F(0.0));
104 Assert(cd->Omega < F(2.0));
e3f82424
MW		 105
106 KernelDataAa * kda = KDA(kd);
107
108 PdfT * src = kd->PdfsActive;
109
110 int maxIterations = cd->MaxIterations;
111
112 #ifdef VTK_OUTPUT
113 if (cd->VtkOutput) {
114 kd->PdfsActive = src;
115 VtkWrite(ld, kd, cd, -1);
116 }
117 #endif
118
119 #ifdef STATISTICS
120 kd->PdfsActive = src;
121 KernelStatistics(kd, ld, cd, 0);
122 #endif
123
124 Assert((maxIterations % 2) == 0);
125
126 for (int iter = 0; iter < maxIterations; iter += 2) {
127
128 // --------------------------------------------------------------------
129 // even time step
130 // --------------------------------------------------------------------
131
132 X_LIKWID_START("aa-vec-even");
133
134 #pragma omp parallel
135 {
136 KernelEven(ld, kd, cd);
137 }
138
139 X_LIKWID_STOP("aa-vec-even");
140
141 // Fixup bounce back PDFs.
142 #ifdef _OPENMP
143 #pragma omp parallel for default(none) \
144 shared(kd, src)
145 #endif
146 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
147 src[kd->BounceBackPdfsSrc[i]] = src[kd->BounceBackPdfsDst[i]];
148 }
149
150 // save current iteration
151 kda->Iteration = iter;
152
153 #ifdef VERIFICATION
154 kd->PdfsActive = src;
155 KernelAddBodyForce(kd, ld, cd);
156 #endif
157
158 #ifdef VTK_OUTPUT
159 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
160 kd->PdfsActive = src;
161 VtkWrite(ld, kd, cd, iter);
162 }
163 #endif
164
165 #ifdef STATISTICS
166 kd->PdfsActive = src;
167 KernelStatistics(kd, ld, cd, iter);
168 #endif
169
170 // --------------------------------------------------------------------
171 // odd time step
172 // --------------------------------------------------------------------
173
174 X_LIKWID_START("aa-vec-odd");
175
176 #pragma omp parallel
177 {
178 KernelOdd(ld, kd, cd);
179 }
180
181 // Stop counters before bounce back. Else computing loop balance will
182 // be incorrect.
183
184 X_LIKWID_STOP("aa-vec-odd");
185
186 // Fixup bounce back PDFs.
187 #ifdef _OPENMP
188 #pragma omp parallel for default(none) \
189 shared(kd, src)
190 #endif
191 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
192 src[kd->BounceBackPdfsDst[i]] = src[kd->BounceBackPdfsSrc[i]];
193 }
194
195 // save current iteration
196 kda->Iteration = iter + 1;
197
198 #ifdef VERIFICATION
199 kd->PdfsActive = src;
200 KernelAddBodyForce(kd, ld, cd);
201 #endif
202
203 #ifdef VTK_OUTPUT
204 if (cd->VtkOutput && ((iter + 1) % cd->VtkModulus) == 0) {
205 kd->PdfsActive = src;
206 VtkWrite(ld, kd, cd, iter + 1);
207 }
208 #endif
209
210 #ifdef STATISTICS
211 kd->PdfsActive = src;
212 KernelStatistics(kd, ld, cd, iter + 1);
213 #endif // }}}
214
215
216 } // for (int iter = 0; ...
217
218 #ifdef VTK_OUTPUT
219
220 if (cd->VtkOutput) {
221 kd->PdfsActive = src;
222 VtkWrite(ld, kd, cd, maxIterations);
223 }
224
225 #endif
226
227 return;
228}
229
230static void KernelEven(LatticeDesc * ld, KernelData * kd, CaseData * cd) // {{{
231{
232 Assert(ld != NULL);
233 Assert(kd != NULL);
234 Assert(cd != NULL);
235
0fde6e45
MW		 236 Assert(cd->Omega > F(0.0));
237 Assert(cd->Omega < F(2.0));
e3f82424
MW		 238
239 KernelDataAa * kda = KDA(kd);
240
241 int nX = ld->Dims[0];
242 int nY = ld->Dims[1];
243 int nZ = ld->Dims[2];
244
245 int * gDims = kd->GlobalDims;
246
247 int oX = kd->Offsets[0];
248 int oY = kd->Offsets[1];
249 int oZ = kd->Offsets[2];
250
251 int blk[3];
252 blk[0] = kda->Blk[0];
253 blk[1] = kda->Blk[1];
254 blk[2] = kda->Blk[2];
255
256 PdfT omega = cd->Omega;
257 PdfT omegaEven = omega;
258
0fde6e45
MW		 259 PdfT magicParam = F(1.0) / F(12.0);
260 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
e3f82424 261
0fde6e45
MW		 262 const PdfT w_0 = F(1.0) / F(3.0);
263 const PdfT w_1 = F(1.0) / F(18.0);
264 const PdfT w_2 = F(1.0) / F(36.0);
e3f82424 265
0fde6e45
MW		 266 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
267 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
e3f82424
MW		 268
269
0fde6e45
MW		 270 VPDFT VONE_HALF = VSET(F(0.5));
271 VPDFT VTHREE_HALF = VSET(F(3.0) / F(2.0));
e3f82424
MW		 272
273 VPDFT vw_1_indep, vw_2_indep;
274 VPDFT vw_0 = VSET(w_0);
275 VPDFT vw_1 = VSET(w_1);
276 VPDFT vw_2 = VSET(w_2);
277
278 VPDFT vw_1_x3 = VSET(w_1_x3);
279 VPDFT vw_2_x3 = VSET(w_2_x3);
280 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
281 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
282
283 VPDFT vui, vux, vuy, vuz, vdens;
284
285 VPDFT vevenPart, voddPart, vdir_indep_trm;
286
287 VPDFT vomegaEven = VSET(omegaEven);
288 VPDFT vomegaOdd = VSET(omegaOdd);
289
290 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
291 #define X(name, idx, idxinv, x, y, z) VPDFT JOIN(vpdf_,name);
292 D3Q19_LIST
293 #undef X
294
295 PdfT * src = kd->Pdfs[0];
296
297 int nThreads = 1;
298 int threadId = 0;
299
300 #ifdef _OPENMP
301 nThreads = omp_get_max_threads();
302 threadId = omp_get_thread_num();
303 #endif
304
305 // TODO: Currently only a 1-D decomposition is applied. For achritectures
306 // with a lot of cores we want at least 2-D.
307
308 int threadStartX = nX / nThreads * threadId;
309 int threadEndX = nX / nThreads * (threadId + 1);
310
311 if (nX % nThreads > 0) {
312 if (nX % nThreads > threadId) {
313 threadStartX += threadId;
314 threadEndX += threadId + 1;
315 }
316 else {
317 threadStartX += nX % nThreads;
318 threadEndX += nX % nThreads;
319 }
320 }
321
322 AssertMsg((blk[2] % VSIZE == 0) || blk[2] >= nZ, "Blocking in z direction must be a multiple of VSIZE = %d or larger than z dimension.", VSIZE);
323
324 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
325 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
326 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
327
328 int eX = MIN(bX + blk[0], threadEndX + oX);
329 int eY = MIN(bY + blk[1], nY + oY);
330 int eZ = MIN(bZ + blk[2], nZ + oZ);
331
332 for (int x = bX; x < eX; x += 1) {
333 for (int y = bY; y < eY; y += 1) {
334 for (int z = bZ; z < eZ; z += VSIZE) {
335
336 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
337
338 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
339 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(&src[I(x, y, z, idx)]);
340 D3Q19_LIST
341 #undef X
342
343
344 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);
345 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);
346 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);
347
348 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)),
349 VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),
350 VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
351
352 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
353
354 VSTU(&src[I(x, y, z, D3Q19_C)],VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
355
356 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
357
358 vui = vuy;
359 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_N,vpdf_S)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
360 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_N,vpdf_S)),VMUL(vui,vw_1_x3)));
361 VSTU(&src[I(x, y, z, D3Q19_S)],VSUB(VSUB(vpdf_N,vevenPart),voddPart));
362 VSTU(&src[I(x, y, z, D3Q19_N)],VADD(VSUB(vpdf_S,vevenPart),voddPart));
363
364 vui = vux;
365 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_E,vpdf_W)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
366 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_E,vpdf_W)),VMUL(vui,vw_1_x3)));
367 VSTU(&src[I(x, y, z, D3Q19_W)],VSUB(VSUB(vpdf_E,vevenPart),voddPart));
368 VSTU(&src[I(x, y, z, D3Q19_E)],VADD(VSUB(vpdf_W,vevenPart),voddPart));
369
370 vui = vuz;
371 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_T,vpdf_B)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
372 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_T,vpdf_B)),VMUL(vui,vw_1_x3)));
373 VSTU(&src[I(x, y, z, D3Q19_B)],VSUB(VSUB(vpdf_T,vevenPart),voddPart));
374 VSTU(&src[I(x, y, z, D3Q19_T)],VADD(VSUB(vpdf_B,vevenPart),voddPart));
375
376 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
377
378 vui = VSUB(vuy,vux);
379 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NW,vpdf_SE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
380 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NW,vpdf_SE)),VMUL(vui,vw_2_x3)));
381 VSTU(&src[I(x, y, z, D3Q19_SE)],VSUB(VSUB(vpdf_NW,vevenPart),voddPart));
382 VSTU(&src[I(x, y, z, D3Q19_NW)],VADD(VSUB(vpdf_SE,vevenPart),voddPart));
383
384 vui = VADD(vux,vuy);
385 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NE,vpdf_SW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
386 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NE,vpdf_SW)),VMUL(vui,vw_2_x3)));
387 VSTU(&src[I(x, y, z, D3Q19_SW)],VSUB(VSUB(vpdf_NE,vevenPart),voddPart));
388 VSTU(&src[I(x, y, z, D3Q19_NE)],VADD(VSUB(vpdf_SW,vevenPart),voddPart));
389
390 vui = VSUB(vuz,vux);
391 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TW,vpdf_BE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
392 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TW,vpdf_BE)),VMUL(vui,vw_2_x3)));
393 VSTU(&src[I(x, y, z, D3Q19_BE)],VSUB(VSUB(vpdf_TW,vevenPart),voddPart));
394 VSTU(&src[I(x, y, z, D3Q19_TW)],VADD(VSUB(vpdf_BE,vevenPart),voddPart));
395
396 vui = VADD(vux,vuz);
397 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TE,vpdf_BW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
398 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TE,vpdf_BW)),VMUL(vui,vw_2_x3)));
399 VSTU(&src[I(x, y, z, D3Q19_BW)],VSUB(VSUB(vpdf_TE,vevenPart),voddPart));
400 VSTU(&src[I(x, y, z, D3Q19_TE)],VADD(VSUB(vpdf_BW,vevenPart),voddPart));
401
402 vui = VSUB(vuz,vuy);
403 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TS,vpdf_BN)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
404 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TS,vpdf_BN)),VMUL(vui,vw_2_x3)));
405 VSTU(&src[I(x, y, z, D3Q19_BN)],VSUB(VSUB(vpdf_TS,vevenPart),voddPart));
406 VSTU(&src[I(x, y, z, D3Q19_TS)],VADD(VSUB(vpdf_BN,vevenPart),voddPart));
407
408 vui = VADD(vuy,vuz);
409 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TN,vpdf_BS)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
410 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TN,vpdf_BS)),VMUL(vui,vw_2_x3)));
411 VSTU(&src[I(x, y, z, D3Q19_BS)],VSUB(VSUB(vpdf_TN,vevenPart),voddPart));
412 VSTU(&src[I(x, y, z, D3Q19_TN)],VADD(VSUB(vpdf_BS,vevenPart),voddPart));
413
414 #undef I
415 } } } // x, y, z
416 } } } // blocked x, y, z
417
418
419
420 return;
421} // }}}
422
423
424static void KernelOdd(LatticeDesc * ld, KernelData * kd, CaseData * cd) // {{{
425{
426 Assert(ld != NULL);
427 Assert(kd != NULL);
428 Assert(cd != NULL);
429
0fde6e45
MW		 430 Assert(cd->Omega > F(0.0));
431 Assert(cd->Omega < F(2.0));
e3f82424
MW		 432
433 KernelDataAa * kda = KDA(kd);
434
435 int nX = ld->Dims[0];
436 int nY = ld->Dims[1];
437 int nZ = ld->Dims[2];
438
439 int * gDims = kd->GlobalDims;
440
441 int oX = kd->Offsets[0];
442 int oY = kd->Offsets[1];
443 int oZ = kd->Offsets[2];
444
445 int blk[3];
446 blk[0] = kda->Blk[0];
447 blk[1] = kda->Blk[1];
448 blk[2] = kda->Blk[2];
449
450 PdfT omega = cd->Omega;
451 PdfT omegaEven = omega;
452
0fde6e45
MW		 453 PdfT magicParam = F(1.0) / F(12.0);
454 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
e3f82424 455
0fde6e45
MW		 456 const PdfT w_0 = F(1.0) / F(3.0);
457 const PdfT w_1 = F(1.0) / F(18.0);
458 const PdfT w_2 = F(1.0) / F(36.0);
e3f82424 459
0fde6e45
MW		 460 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
461 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
e3f82424 462
0fde6e45
MW		 463 VPDFT VONE_HALF = VSET(F(0.5));
464 VPDFT VTHREE_HALF = VSET(F(3.0) / F(2.0));
e3f82424
MW		 465
466 VPDFT vw_1_indep, vw_2_indep;
467 VPDFT vw_0 = VSET(w_0);
468 VPDFT vw_1 = VSET(w_1);
469 VPDFT vw_2 = VSET(w_2);
470
471 VPDFT vw_1_x3 = VSET(w_1_x3);
472 VPDFT vw_2_x3 = VSET(w_2_x3);
473 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
474 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
475
476 VPDFT vui, vux, vuy, vuz, vdens;
477
478 VPDFT vevenPart, voddPart, vdir_indep_trm;
479
480 VPDFT vomegaEven = VSET(omegaEven);
481 VPDFT vomegaOdd = VSET(omegaOdd);
482
483 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
484 #define X(name, idx, idxinv, x, y, z) VPDFT JOIN(vpdf_,name);
485 D3Q19_LIST
486 #undef X
487
488 PdfT * src = kd->Pdfs[0];
489
490 int threadId = 0;
491 int nThreads = 1;
492
493 #ifdef _OPENMP
494 threadId = omp_get_thread_num();
495 nThreads = omp_get_max_threads();
496 #endif
497
498 // TODO: Currently only a 1-D decomposition is applied. For achritectures
499 // with a lot of cores we want at least 2-D.
500 int threadStartX = nX / nThreads * threadId;
501 int threadEndX = nX / nThreads * (threadId + 1);
502
503 if (nX % nThreads > 0) {
504 if (nX % nThreads > threadId) {
505 threadStartX += threadId;
506 threadEndX += threadId + 1;
507 }
508 else {
509 threadStartX += nX % nThreads;
510 threadEndX += nX % nThreads;
511 }
512 }
513
514 AssertMsg((blk[2] % VSIZE == 0) || blk[2] >= nZ, "Blocking in z direction must be a multiple of VSIZE = %d or larger than z dimension.", VSIZE);
515
516 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
517 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
518 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
519
520 int eX = MIN(bX + blk[0], threadEndX + oX);
521 int eY = MIN(bY + blk[1], nY + oY);
522 int eZ = MIN(bZ + blk[2], nZ + oZ);
523
524 for (int x = bX; x < eX; ++x) {
525 for (int y = bY; y < eY; ++y) {
526 for (int z = bZ; z < eZ; z += VSIZE) {
527
528 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
529
530
531 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(&src[I(x - _x, y - _y, z - _z, idxinv)]);
532 D3Q19_LIST
533 #undef X
534
535 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);
536 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);
537 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);
538
539 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)),
540 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));
541
542 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
543
544 VSTU(&src[I(x, y, z, D3Q19_C)],VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
545
546 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
547
548 vui = vuy;
549 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_N,vpdf_S)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
550 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_N,vpdf_S)),VMUL(vui,vw_1_x3)));
551 VSTU(&src[I(x, y + 1, z, D3Q19_N)], VSUB(VSUB(vpdf_N,vevenPart),voddPart));
552 VSTU(&src[I(x, y - 1, z, D3Q19_S)], VADD(VSUB(vpdf_S,vevenPart),voddPart));
553
554 vui = vux;
555 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_E,vpdf_W)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
556 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_E,vpdf_W)),VMUL(vui,vw_1_x3)));
557 VSTU(&src[I(x + 1, y, z, D3Q19_E)], VSUB(VSUB(vpdf_E,vevenPart),voddPart));
558 VSTU(&src[I(x - 1, y, z, D3Q19_W)], VADD(VSUB(vpdf_W,vevenPart),voddPart));
559
560 vui = vuz;
561 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_T,vpdf_B)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
562 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_T,vpdf_B)),VMUL(vui,vw_1_x3)));
563 VSTU(&src[I(x, y, z + 1, D3Q19_T)], VSUB(VSUB(vpdf_T,vevenPart),voddPart));
564 VSTU(&src[I(x, y, z - 1, D3Q19_B)], VADD(VSUB(vpdf_B,vevenPart),voddPart));
565
566 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
567
568 vui = VSUB(vuy,vux);
569 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NW,vpdf_SE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
570 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NW,vpdf_SE)),VMUL(vui,vw_2_x3)));
571 VSTU(&src[I(x - 1, y + 1, z, D3Q19_NW)], VSUB(VSUB(vpdf_NW,vevenPart),voddPart));
572 VSTU(&src[I(x + 1, y - 1, z, D3Q19_SE)], VADD(VSUB(vpdf_SE,vevenPart),voddPart));
573
574 vui = VADD(vux,vuy);
575 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NE,vpdf_SW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
576 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NE,vpdf_SW)),VMUL(vui,vw_2_x3)));
577 VSTU(&src[I(x + 1, y + 1, z, D3Q19_NE)], VSUB(VSUB(vpdf_NE,vevenPart),voddPart));
578 VSTU(&src[I(x - 1, y - 1, z, D3Q19_SW)], VADD(VSUB(vpdf_SW,vevenPart),voddPart));
579
580 vui = VSUB(vuz,vux);
581 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TW,vpdf_BE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
582 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TW,vpdf_BE)),VMUL(vui,vw_2_x3)));
583 VSTU(&src[I(x - 1, y, z + 1, D3Q19_TW)], VSUB(VSUB(vpdf_TW,vevenPart),voddPart));
584 VSTU(&src[I(x + 1, y, z - 1, D3Q19_BE)], VADD(VSUB(vpdf_BE,vevenPart),voddPart));
585
586 vui = VADD(vux,vuz);
587 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TE,vpdf_BW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
588 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TE,vpdf_BW)),VMUL(vui,vw_2_x3)));
589 VSTU(&src[I(x + 1, y, z + 1, D3Q19_TE)], VSUB(VSUB(vpdf_TE,vevenPart),voddPart));
590 VSTU(&src[I(x - 1, y, z - 1, D3Q19_BW)], VADD(VSUB(vpdf_BW,vevenPart),voddPart));
591
592 vui = VSUB(vuz,vuy);
593 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TS,vpdf_BN)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
594 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TS,vpdf_BN)),VMUL(vui,vw_2_x3)));
595 VSTU(&src[I(x, y - 1, z + 1, D3Q19_TS)], VSUB(VSUB(vpdf_TS,vevenPart),voddPart));
596 VSTU(&src[I(x, y + 1, z - 1, D3Q19_BN)], VADD(VSUB(vpdf_BN,vevenPart),voddPart));
597
598 vui = VADD(vuy,vuz);
599 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TN,vpdf_BS)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
600 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TN,vpdf_BS)),VMUL(vui,vw_2_x3)));
601 VSTU(&src[I(x, y + 1, z + 1, D3Q19_TN)], VSUB(VSUB(vpdf_TN,vevenPart),voddPart));
602 VSTU(&src[I(x, y - 1, z - 1, D3Q19_BS)], VADD(VSUB(vpdf_BS,vevenPart),voddPart));
603
604 #undef I
605 } } } // x, y, z
606 } } } // blocked x, y, z
607
608 return;
609
610} // }}}
LBM benchmark kernels
This page took 0.066636 seconds and 5 git commands to generate.