projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
merge with kernels from MH's master thesis
[LbmBenchmarkKernelsPublic.git] / src / BenchKernelD3Q19.c
CommitLineData
10988083
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 "BenchKernelD3Q19Common.h"
28
29#include "Memory.h"
30#include "Vtk.h"
e3f82424 31#include "LikwidIf.h"
10988083
MW		 32
33#include <inttypes.h>
34#include <math.h>
35
36#ifdef _OPENMP
37 #include <omp.h>
38#endif
39
40void FNAME(D3Q19Kernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
41{
42 Assert(ld != NULL);
43 Assert(kernelData != NULL);
44 Assert(cd != NULL);
45
0fde6e45
MW		 46 Assert(cd->Omega > F(0.0));
47 Assert(cd->Omega < F(2.0));
10988083
MW		 48
49 KernelData * kd = (KernelData *)kernelData;
50
51
52 int nX = ld->Dims[0];
53 int nY = ld->Dims[1];
54 int nZ = ld->Dims[2];
55
56 int * gDims = kd->GlobalDims;
57
58 int oX = kd->Offsets[0];
59 int oY = kd->Offsets[1];
60 int oZ = kd->Offsets[2];
61
62 PdfT omega = cd->Omega;
63 PdfT omegaEven = omega;
0fde6e45
MW		 64 // PdfT omegaOdd = 8.0*((F(2.0)-omegaEven)/(8.0-omegaEven)); //"standard" trt odd relaxation parameter
65 PdfT magicParam = F(1.0) / F(12.0);
66 // 1/ 4: best stability;
67 // 1/12: removes third-order advection error (best advection);
68 // 1/ 6: removes fourth-order diffusion error (best diffusion);
69 // 3/16: exact location of bounce back for poiseuille flow
70 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
10988083 71
0fde6e45
MW		 72 PdfT evenPart = F(0.0);
73 PdfT oddPart = F(0.0);
74 PdfT dir_indep_trm = F(0.0);
10988083 75
0fde6e45
MW		 76 PdfT w_0 = F(1.0) / F( 3.0);
77 PdfT w_1 = F(1.0) / F(18.0);
78 PdfT w_2 = F(1.0) / F(36.0);
10988083 79
0fde6e45
MW		 80 PdfT w_1_x3 = w_1 * F(3.0); PdfT w_1_nine_half = w_1 * F(9.0)/F(2.0); PdfT w_1_indep = F(0.0);
81 PdfT w_2_x3 = w_2 * F(3.0); PdfT w_2_nine_half = w_2 * F(9.0)/F(2.0); PdfT w_2_indep = F(0.0);
10988083
MW		 82
83 PdfT ux, uy, uz, ui;
84 PdfT dens;
85
86 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
87 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
88 D3Q19_LIST
89 #undef X
90
91 PdfT * src = kd->Pdfs[0];
92 PdfT * dst = kd->Pdfs[1];
93 PdfT * tmp;
94
95 int maxIterations = cd->MaxIterations;
96
97 #ifdef VTK_OUTPUT
98 if (cd->VtkOutput) {
99 kd->PdfsActive = src;
100 VtkWrite(ld, kd, cd, 0);
101 }
102 #endif
103
8cafd9ea
MW		 104 X_KERNEL_START(kernelData);
105
10988083
MW		 106 for (int iter = 0; iter < maxIterations; ++iter) {
107
e3f82424
MW		 108 X_LIKWID_START("os");
109
10988083 110 #ifdef _OPENMP
0fde6e45 111 #pragma omp parallel for collapse(2) default(none) \
10988083
MW		 112 shared(gDims,src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \
113 w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd, \
114 oX, oY, oZ, nX, nY, nZ) \
115 private(ux, uy, uz, ui, dens, dir_indep_trm, \
116 pdf_C, \
117 pdf_N, pdf_E, pdf_S, pdf_W, \
118 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
119 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
120 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
121 evenPart, oddPart, w_1_indep, w_2_indep)
122 #endif
e3f82424 123 for (int x = oX; x < nX + oX; ++x) {
10988083 124 for (int y = oY; y < nY + oY; ++y) {
0fde6e45
MW		 125 #ifdef INTEL_OPT_DIRECTIVES
126 #pragma ivdep
127 #pragma vector always
128 #pragma simd
129 #endif
8cafd9ea 130 for (int z = oZ; z < nZ + oZ; ++z) { // LOOP os
10988083
MW		 131 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
132
133#ifdef PROP_MODEL_PUSH
134
135 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
136 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x, y, z, idx)];
137 //if (isnan(JOIN(pdf_,name))) { printf("iter: %d %d %d %d %d %s nan\n", iter, x-oX, y-oY, z-oZ, idx, D3Q19_NAMES[idx]); exit(1);}
138 D3Q19_LIST
139 #undef X
140
141#elif PROP_MODEL_PULL
142
143 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
144 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x - _x, y - _y, z - _z, idx)];
145 //if (isnan(JOIN(pdf_,name))) { printf("iter: %d %d %d %d %d %s nan\n", iter, x-oX, y-oY, z-oZ, idx, D3Q19_NAMES[idx]); exit(1);}
146 D3Q19_LIST
147 #undef X
148
149#else
150 #error No implementation for PROP_MODEL_NAME.
151#endif
152
153 // #define LID_DRIVEN_CAVITY
154
155 #ifdef LID_DRIVEN_CAVITY
156
157 if (z == nZ - 4 + oZ && x > 3 + oX && x < (nX - 4 + oX) && y > 3 + oY && y < (nY - 4 + oY)) {
0fde6e45
MW		 158 ux = F(0.1 * 0.577);
159 uy = F(0.0);
160 uz = F(0.0);
10988083
MW		 161
162 } else {
163 #endif
164 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
165 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
166 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
167 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
168 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
169 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
170 #ifdef LID_DRIVEN_CAVITY
171 }
172
173 #endif
174
175 dens = pdf_C +
176 pdf_N + pdf_E + pdf_S + pdf_W +
177 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
178 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
179 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
180
0fde6e45 181 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0);
10988083
MW		 182
183#ifdef PROP_MODEL_PUSH
184
185 // direction: w_0
186 dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
187
188 // direction: w_1
189 w_1_indep = w_1*dir_indep_trm;
190
191 ui = uy;
0fde6e45
MW		 192 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
193 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
10988083
MW		 194 dst[I(x, y + 1, z, D3Q19_N)] = pdf_N - evenPart - oddPart;
195 dst[I(x, y - 1, z, D3Q19_S)] = pdf_S - evenPart + oddPart;
196
197 ui = ux;
0fde6e45
MW		 198 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
199 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
10988083
MW		 200 dst[I(x + 1, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart;
201 dst[I(x - 1, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart;
202
203 ui = uz;
0fde6e45
MW		 204 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
205 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
10988083
MW		 206 dst[I(x, y, z + 1, D3Q19_T)] = pdf_T - evenPart - oddPart;
207 dst[I(x, y, z - 1, D3Q19_B)] = pdf_B - evenPart + oddPart;
208
209 // direction: w_2
210 w_2_indep = w_2*dir_indep_trm;
211
212 ui = -ux + uy;
0fde6e45
MW		 213 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
214 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083
MW		 215 dst[I(x - 1, y + 1, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
216 dst[I(x + 1, y - 1, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
217
218 ui = ux + uy;
0fde6e45
MW		 219 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
220 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083
MW		 221 dst[I(x + 1, y + 1, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
222 dst[I(x - 1, y - 1, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
223
224 ui = -ux + uz;
0fde6e45
MW		 225 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
226 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083
MW		 227 dst[I(x - 1, y, z + 1, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
228 dst[I(x + 1, y, z - 1, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
229
230 ui = ux + uz;
0fde6e45
MW		 231 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
232 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083
MW		 233 dst[I(x + 1, y, z + 1, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
234 dst[I(x - 1, y, z - 1, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
235
236 ui = -uy + uz;
0fde6e45
MW		 237 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
238 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083
MW		 239 dst[I(x, y - 1, z + 1, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
240 dst[I(x, y + 1, z - 1, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
241
242 ui = uy + uz;
0fde6e45
MW		 243 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
244 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083
MW		 245 dst[I(x, y + 1, z + 1, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
246 dst[I(x, y - 1, z - 1, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
247
248#elif PROP_MODEL_PULL
249
250 // direction: w_0
251 dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
252
253 // direction: w_1
254 w_1_indep = w_1*dir_indep_trm;
255
256 ui = uy;
0fde6e45
MW		 257 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
258 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
10988083
MW		 259 dst[I(x, y, z, D3Q19_N)] = pdf_N - evenPart - oddPart;
260 dst[I(x, y, z, D3Q19_S)] = pdf_S - evenPart + oddPart;
261
262 ui = ux;
0fde6e45
MW		 263 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
264 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
10988083
MW		 265 dst[I(x, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart;
266 dst[I(x, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart;
267
268 ui = uz;
0fde6e45
MW		 269 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
270 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
10988083
MW		 271 dst[I(x, y, z, D3Q19_T)] = pdf_T - evenPart - oddPart;
272 dst[I(x, y, z, D3Q19_B)] = pdf_B - evenPart + oddPart;
273
274 // direction: w_2
275 w_2_indep = w_2*dir_indep_trm;
276
277 ui = -ux + uy;
0fde6e45
MW		 278 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
279 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083
MW		 280 dst[I(x, y, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
281 dst[I(x, y, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
282
283 ui = ux + uy;
0fde6e45
MW		 284 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
285 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083
MW		 286 dst[I(x, y, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
287 dst[I(x, y, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
288
289 ui = -ux + uz;
0fde6e45
MW		 290 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
291 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083
MW		 292 dst[I(x, y, z, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
293 dst[I(x, y, z, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
294
295 ui = ux + uz;
0fde6e45
MW		 296 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
297 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083
MW		 298 dst[I(x, y, z, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
299 dst[I(x, y, z, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
300
301 ui = -uy + uz;
0fde6e45
MW		 302 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
303 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083
MW		 304 dst[I(x, y, z, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
305 dst[I(x, y, z, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
306
307 ui = uy + uz;
0fde6e45
MW		 308 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
309 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083
MW		 310 dst[I(x, y, z, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
311 dst[I(x, y, z, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
312
313#else
314 #error No implementation for PROP_MODEL_NAME.
315#endif
316
317 #undef I
318 }
319 }
320 } // z, y, x (from inner to outer)
321
e3f82424
MW		 322 // Stop counters before bounce back. Else computing loop balance will be incorrect.
323 X_LIKWID_STOP("os");
324
8cafd9ea 325
10988083
MW		 326 // Fixup bounce back PDFs.
327 #ifdef _OPENMP
328 #pragma omp parallel for default(none) \
329 shared(kd, dst)
330 #endif
331 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
332 dst[kd->BounceBackPdfsDst[i]] = dst[kd->BounceBackPdfsSrc[i]];
333 }
334
335 #ifdef VERIFICATION
336 kd->PdfsActive = dst;
337 KernelAddBodyForce(kd, ld, cd);
338 #endif
339
340 #ifdef VTK_OUTPUT
341
342 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
343 kd->PdfsActive = dst;
344 VtkWrite(ld, kd, cd, iter);
345 }
346
347 #endif
348
349 #ifdef STATISTICS
350 kd->PdfsActive = dst;
351 KernelStatistics(kd, ld, cd, iter);
352 #endif
353
354 // swap grids
355 tmp = src;
356 src = dst;
357 dst = tmp;
358
359 } // for (int iter = 0; ...
360
8cafd9ea
MW		 361 X_KERNEL_END(kernelData);
362
10988083
MW		 363 #ifdef VTK_OUTPUT
364
365 if (cd->VtkOutput) {
366 kd->PdfsActive = src;
367 VtkWrite(ld, kd, cd, maxIterations);
368 }
369
370 #endif
371
372 return;
373}
374
375
376void FNAME(D3Q19BlkKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
377{
378 Assert(ld != NULL);
379 Assert(kernelData != NULL);
380 Assert(cd != NULL);
381
0fde6e45
MW		 382 Assert(cd->Omega > F(0.0));
383 Assert(cd->Omega < F(2.0));
10988083
MW		 384
385 KernelData * kd = (KernelData *)kernelData;
386
387
388 int nX = ld->Dims[0];
389 int nY = ld->Dims[1];
390 int nZ = ld->Dims[2];
391
392 int * gDims = kd->GlobalDims;
393
394 int oX = kd->Offsets[0];
395 int oY = kd->Offsets[1];
396 int oZ = kd->Offsets[2];
397
398 KernelDataEx * kdex = (KernelDataEx *)kd;
399
400 int blk[3];
401 blk[0] = kdex->Blk[0];
402 blk[1] = kdex->Blk[1];
403 blk[2] = kdex->Blk[2];
404
405 PdfT omega = cd->Omega;
406 PdfT omegaEven = omega;
0fde6e45
MW		 407 // PdfT omegaOdd = 8.0*((F(2.0)-omegaEven)/(8.0-omegaEven)); //"standard" trt odd relaxation parameter
408 PdfT magicParam = F(1.0)/F(12.0);
409 // 1/ 4: best stability;
410 // 1/12: removes third-order advection error (best advection);
411 // 1/ 6: removes fourth-order diffusion error (best diffusion);
412 // 3/16: exact location of bounce back for poiseuille flow
413 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
10988083 414
0fde6e45
MW		 415 PdfT evenPart = F(0.0);
416 PdfT oddPart = F(0.0);
417 PdfT dir_indep_trm = F(0.0);
10988083 418
0fde6e45
MW		 419 PdfT w_0 = F(1.0) / F( 3.0);
420 PdfT w_1 = F(1.0) / F(18.0);
421 PdfT w_2 = F(1.0) / F(36.0);
10988083 422
0fde6e45
MW		 423 PdfT w_1_x3 = w_1 * F(3.0); PdfT w_1_nine_half = w_1 * F(9.0)/F(2.0); PdfT w_1_indep = F(0.0);
424 PdfT w_2_x3 = w_2 * F(3.0); PdfT w_2_nine_half = w_2 * F(9.0)/F(2.0); PdfT w_2_indep = F(0.0);
10988083
MW		 425
426 PdfT ux, uy, uz, ui;
427 PdfT dens;
428
429 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
430 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
431 D3Q19_LIST
432 #undef X
433
434 PdfT * src = kd->Pdfs[0];
435 PdfT * dst = kd->Pdfs[1];
436 PdfT * tmp;
437
438 int maxIterations = cd->MaxIterations;
439
440 #ifdef VTK_OUTPUT
441 if (cd->VtkOutput) {
442 kd->PdfsActive = src;
443 VtkWrite(ld, kd, cd, 0);
444 }
445 #endif
446
447 int nThreads = 1;
448
449 #ifdef _OPENMP
450 nThreads = omp_get_max_threads();
451 #endif
452
8cafd9ea
MW		 453 X_KERNEL_START(kernelData);
454
10988083
MW		 455 for (int iter = 0; iter < maxIterations; ++iter) {
456
10988083 457 #ifdef _OPENMP
e3f82424 458 #pragma omp parallel default(none) \
10988083
MW		 459 shared(gDims,src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \
460 w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd, \
461 oX, oY, oZ, nX, nY, nZ, blk, nThreads) \
462 private(ux, uy, uz, ui, dens, dir_indep_trm, \
463 pdf_C, \
464 pdf_N, pdf_E, pdf_S, pdf_W, \
465 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
466 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
467 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
468 evenPart, oddPart, w_1_indep, w_2_indep)
469 #endif
e3f82424
MW		 470 {
471 X_LIKWID_START("blk-os");
10988083 472
e3f82424 473 int threadId = omp_get_thread_num();
10988083 474
e3f82424
MW		 475 int threadStartX = nX / nThreads * threadId;
476 int threadEndX = nX / nThreads * (threadId + 1);
10988083
MW		 477
478 if (nX % nThreads > 0) {
e3f82424
MW		 479 if (nX % nThreads > threadId) {
480 threadStartX += threadId;
481 threadEndX += threadId + 1;
10988083
MW		 482 }
483 else {
484 threadStartX += nX % nThreads;
485 threadEndX += nX % nThreads;
486 }
487 }
488
e3f82424 489 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
10988083 490 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
e3f82424 491 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
10988083
MW		 492
493 // Must do everything here, else it would break collapse.
494 int eZ = MIN(bZ + blk[2], nZ + oZ);
495 int eY = MIN(bY + blk[1], nY + oY);
496 int eX = MIN(bX + blk[0], threadEndX + oX);
497
e3f82424 498 for (int x = bX; x < eX; ++x) {
10988083 499 for (int y = bY; y < eY; ++y) {
0fde6e45
MW		 500 #ifdef INTEL_OPT_DIRECTIVES
501 #pragma ivdep
502 #pragma vector always
503 #pragma simd
504 #endif
e3f82424 505 for (int z = bZ; z < eZ; ++z) {
10988083
MW		 506
507 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
508
509#ifdef PROP_MODEL_PUSH
510
511 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
512 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x, y, z, idx)];
10988083
MW		 513 D3Q19_LIST
514 #undef X
515
516#elif PROP_MODEL_PULL
517
518 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
519 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x - _x, y - _y, z - _z, idx)];
10988083
MW		 520 D3Q19_LIST
521 #undef X
522
523#else
524 #error No implementation for PROP_MODEL_NAME.
525#endif
526
527 // #define LID_DRIVEN_CAVITY
528
529 #ifdef LID_DRIVEN_CAVITY
530
531 if (z == nZ - 4 + oZ && x > 3 + oX && x < (nX - 4 + oX) && y > 3 + oY && y < (nY - 4 + oY)) {
0fde6e45 532 ux = 0.1 * F(0.5)77;
10988083
MW		 533 uy = 0.0;
534 uz = 0.0;
535
536 } else {
537 #endif
538 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
e3f82424 539 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
10988083 540 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
e3f82424 541 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
10988083 542 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
e3f82424 543 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
10988083
MW		 544 #ifdef LID_DRIVEN_CAVITY
545 }
546
547 #endif
548
549 dens = pdf_C +
550 pdf_N + pdf_E + pdf_S + pdf_W +
551 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
552 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
553 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
554
0fde6e45 555 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0);
10988083
MW		 556
557#ifdef PROP_MODEL_PUSH
558
559 // direction: w_0
560 dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
561
562 // direction: w_1
563 w_1_indep = w_1*dir_indep_trm;
564
565 ui = uy;
0fde6e45
MW		 566 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
567 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
10988083
MW		 568 dst[I(x, y + 1, z, D3Q19_N)] = pdf_N - evenPart - oddPart;
569 dst[I(x, y - 1, z, D3Q19_S)] = pdf_S - evenPart + oddPart;
570
571 ui = ux;
0fde6e45
MW		 572 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
573 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
10988083
MW		 574 dst[I(x + 1, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart;
575 dst[I(x - 1, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart;
576
577 ui = uz;
0fde6e45
MW		 578 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
579 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
10988083
MW		 580 dst[I(x, y, z + 1, D3Q19_T)] = pdf_T - evenPart - oddPart;
581 dst[I(x, y, z - 1, D3Q19_B)] = pdf_B - evenPart + oddPart;
582
583 // direction: w_2
584 w_2_indep = w_2*dir_indep_trm;
585
586 ui = -ux + uy;
0fde6e45
MW		 587 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
588 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083
MW		 589 dst[I(x - 1, y + 1, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
590 dst[I(x + 1, y - 1, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
591
592 ui = ux + uy;
0fde6e45
MW		 593 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
594 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083
MW		 595 dst[I(x + 1, y + 1, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
596 dst[I(x - 1, y - 1, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
597
598 ui = -ux + uz;
0fde6e45
MW		 599 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
600 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083
MW		 601 dst[I(x - 1, y, z + 1, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
602 dst[I(x + 1, y, z - 1, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
603
604 ui = ux + uz;
0fde6e45
MW		 605 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
606 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083
MW		 607 dst[I(x + 1, y, z + 1, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
608 dst[I(x - 1, y, z - 1, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
609
610 ui = -uy + uz;
0fde6e45
MW		 611 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
612 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083
MW		 613 dst[I(x, y - 1, z + 1, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
614 dst[I(x, y + 1, z - 1, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
615
616 ui = uy + uz;
0fde6e45
MW		 617 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
618 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083
MW		 619 dst[I(x, y + 1, z + 1, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
620 dst[I(x, y - 1, z - 1, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
621
622#elif PROP_MODEL_PULL
623
624 // direction: w_0
625 dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
626
627 // direction: w_1
628 w_1_indep = w_1*dir_indep_trm;
629
630 ui = uy;
0fde6e45
MW		 631 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
632 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
10988083
MW		 633 dst[I(x, y, z, D3Q19_N)] = pdf_N - evenPart - oddPart;
634 dst[I(x, y, z, D3Q19_S)] = pdf_S - evenPart + oddPart;
635
636 ui = ux;
0fde6e45
MW		 637 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
638 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
10988083
MW		 639 dst[I(x, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart;
640 dst[I(x, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart;
641
642 ui = uz;
0fde6e45
MW		 643 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
644 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
10988083
MW		 645 dst[I(x, y, z, D3Q19_T)] = pdf_T - evenPart - oddPart;
646 dst[I(x, y, z, D3Q19_B)] = pdf_B - evenPart + oddPart;
647
648 // direction: w_2
649 w_2_indep = w_2*dir_indep_trm;
650
651 ui = -ux + uy;
0fde6e45
MW		 652 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
653 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083
MW		 654 dst[I(x, y, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
655 dst[I(x, y, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
656
657 ui = ux + uy;
0fde6e45
MW		 658 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
659 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083
MW		 660 dst[I(x, y, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
661 dst[I(x, y, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
662
663 ui = -ux + uz;
0fde6e45
MW		 664 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
665 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083
MW		 666 dst[I(x, y, z, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
667 dst[I(x, y, z, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
668
669 ui = ux + uz;
0fde6e45
MW		 670 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
671 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083
MW		 672 dst[I(x, y, z, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
673 dst[I(x, y, z, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
674
675 ui = -uy + uz;
0fde6e45
MW		 676 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
677 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083
MW		 678 dst[I(x, y, z, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
679 dst[I(x, y, z, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
680
681 ui = uy + uz;
0fde6e45
MW		 682 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
683 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083
MW		 684 dst[I(x, y, z, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
685 dst[I(x, y, z, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
686
687#else
688 #error No implementation for PROP_MODEL_NAME.
689#endif
690
691 #undef I
692 }
693 }
694 } // z, y, x (from inner to outer)
695 }
696 }
697 } // z, y, x (from inner to outer)
698
e3f82424
MW		 699 X_LIKWID_STOP("blk-os");
700 } // parallel region
701
702 // Stop counters before bounce back. Else computing loop balance will be incorrect.
10988083
MW		 703
704 // Fixup bounce back PDFs.
705 #ifdef _OPENMP
706 #pragma omp parallel for default(none) \
707 shared(kd, dst)
708 #endif
709 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
710 dst[kd->BounceBackPdfsDst[i]] = dst[kd->BounceBackPdfsSrc[i]];
711 }
712
713 #ifdef VERIFICATION
714 kd->PdfsActive = dst;
715 KernelAddBodyForce(kd, ld, cd);
716 #endif
717
718 #ifdef VTK_OUTPUT
719
720 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
721 kd->PdfsActive = dst;
722 VtkWrite(ld, kd, cd, iter);
723 }
724
725 #endif
726
727 #ifdef STATISTICS
728 kd->PdfsActive = dst;
729 KernelStatistics(kd, ld, cd, iter);
730 #endif
731
732 // swap grids
733 tmp = src;
734 src = dst;
735 dst = tmp;
736
737 } // for (int iter = 0; ...
738
8cafd9ea
MW		 739 X_KERNEL_END(kernelData);
740
10988083
MW		 741 #ifdef VTK_OUTPUT
742
743 if (cd->VtkOutput) {
744 kd->PdfsActive = src;
745 VtkWrite(ld, kd, cd, maxIterations);
746 }
747
748 #endif
749
750 return;
751}
752
LBM benchmark kernels
This page took 0.144362 seconds and 5 git commands to generate.