1 // --------------------------------------------------------------------------
4 // Markus Wittmann, 2016-2017
5 // RRZE, University of Erlangen-Nuremberg, Germany
6 // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
9 // LSS, University of Erlangen-Nuremberg, Germany
11 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
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.
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.
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/>.
26 // --------------------------------------------------------------------------
27 #include "BenchKernelD3Q19AaCommon.h"
40 void FNAME(D3Q19AaKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
43 Assert(kernelData != NULL);
46 Assert(cd->Omega > F(0.0));
47 Assert(cd->Omega < F(2.0));
49 KernelData * kd = (KernelData *)kernelData;
56 int * gDims = kd->GlobalDims;
58 int oX = kd->Offsets[0];
59 int oY = kd->Offsets[1];
60 int oZ = kd->Offsets[2];
62 KernelDataAa * kda = KDA(kd);
69 PdfT omega = cd->Omega;
70 PdfT omegaEven = omega;
71 PdfT magicParam = F(1.0) / F(12.0);
72 // 1/4: best stability;
73 // 1/12: removes third-order advection error (best advection);
74 // 1/6: removes fourth-order diffusion error (best diffusion);
75 // 3/16: exact location of bounce back for poiseuille flow
77 PdfT omegaOdd = F(1.0)/( F(0.5) + magicParam/(F(1.0)/omega - F(0.5)) );
79 PdfT evenPart = F(0.0);
80 PdfT oddPart = F(0.0);
81 PdfT dir_indep_trm = F(0.0);
83 PdfT w_0 = F(1.0) / F(3.0);
84 PdfT w_1 = F(1.0) / F(18.0);
85 PdfT w_2 = F(1.0) / F(36.0);
87 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);
88 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);
93 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
94 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
98 PdfT * src = kd->Pdfs[0];
100 int maxIterations = cd->MaxIterations;
104 kd->PdfsActive = src;
105 VtkWrite(ld, kd, cd, -1);
110 kd->PdfsActive = src;
111 KernelStatistics(kd, ld, cd, 0);
118 nThreads = omp_get_max_threads();
121 for (int iter = 0; iter < maxIterations; iter += 2) {
123 // --------------------------------------------------------------------
126 X_LIKWID_START("aa-even");
131 #pragma omp parallel for default(none) \
132 shared(gDims,src, w_0, w_1, w_2, omegaEven, omegaOdd, \
133 w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd, \
134 oX, oY, oZ, nX, nY, nZ, blk, nThreads, ld) \
135 private(ux, uy, uz, ui, dens, dir_indep_trm, \
137 pdf_N, pdf_E, pdf_S, pdf_W, \
138 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
139 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
140 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
141 evenPart, oddPart, w_1_indep, w_2_indep)
144 for (int i = 0; i < nThreads; ++i) {
146 int threadStartX = nX / nThreads * i;
147 int threadEndX = nX / nThreads * (i + 1);
149 if (nX % nThreads > 0) {
150 if (nX % nThreads > i) {
155 threadStartX += nX % nThreads;
156 threadEndX += nX % nThreads;
160 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
161 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
162 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
164 int eX = MIN(bX + blk[0], threadEndX + oX);
165 int eY = MIN(bY + blk[1], nY + oY);
166 int eZ = MIN(bZ + blk[2], nZ + oZ);
168 // printf("%d: %d-%d %d-%d %d-%d %d - %d\n", omp_get_thread_num(), bZ, eZ, bY, eY, bX, eX, threadStartX, threadEndX);
170 for (int x = bX; x < eX; ++x) {
171 for (int y = bY; y < eY; ++y) {
172 #ifdef INTEL_OPT_DIRECTIVES
174 #pragma vector always
177 for (int z = bZ; z < eZ; ++z) {
180 if (ld->Lattice[L_INDEX_4(ld->Dims, x - oX, y - oY, z - oZ)] == LAT_CELL_OBSTACLE) {
184 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
187 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
188 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x, y, z, idx)];
192 // #define LID_DRIVEN_CAVITY
194 #ifdef LID_DRIVEN_CAVITY
196 if (z == nZ - 4 + oZ && x > 3 + oX && x < (nX - 4 + oX) && y > 3 + oY && y < (nY - 4 + oY)) {
197 ux = F(0.1) * F(0.5)77;
203 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
204 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
205 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
206 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
207 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
208 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
209 #ifdef LID_DRIVEN_CAVITY
214 pdf_N + pdf_E + pdf_S + pdf_W +
215 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
216 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
217 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
219 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0);
222 src[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
225 w_1_indep = w_1*dir_indep_trm;
228 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
229 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
230 src[I(x, y, z, D3Q19_S)] = pdf_N - evenPart - oddPart;
231 src[I(x, y, z, D3Q19_N)] = pdf_S - evenPart + oddPart;
234 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
235 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
236 src[I(x, y, z, D3Q19_W)] = pdf_E - evenPart - oddPart;
237 src[I(x, y, z, D3Q19_E)] = pdf_W - evenPart + oddPart;
240 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
241 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
242 src[I(x, y, z, D3Q19_B)] = pdf_T - evenPart - oddPart;
243 src[I(x, y, z, D3Q19_T)] = pdf_B - evenPart + oddPart;
246 w_2_indep = w_2*dir_indep_trm;
249 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
250 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
251 src[I(x, y, z, D3Q19_SE)] = pdf_NW - evenPart - oddPart;
252 src[I(x, y, z, D3Q19_NW)] = pdf_SE - evenPart + oddPart;
255 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
256 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
257 src[I(x, y, z, D3Q19_SW)] = pdf_NE - evenPart - oddPart;
258 src[I(x, y, z, D3Q19_NE)] = pdf_SW - evenPart + oddPart;
261 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
262 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
263 src[I(x, y, z, D3Q19_BE)] = pdf_TW - evenPart - oddPart;
264 src[I(x, y, z, D3Q19_TW)] = pdf_BE - evenPart + oddPart;
267 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
268 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
269 src[I(x, y, z, D3Q19_BW)] = pdf_TE - evenPart - oddPart;
270 src[I(x, y, z, D3Q19_TE)] = pdf_BW - evenPart + oddPart;
273 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
274 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
275 src[I(x, y, z, D3Q19_BN)] = pdf_TS - evenPart - oddPart;
276 src[I(x, y, z, D3Q19_TS)] = pdf_BN - evenPart + oddPart;
279 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
280 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
281 src[I(x, y, z, D3Q19_BS)] = pdf_TN - evenPart - oddPart;
282 src[I(x, y, z, D3Q19_TN)] = pdf_BS - evenPart + oddPart;
285 } } } // z, y, x (from inner to outer)
286 } } } // z, y, x (from inner to outer)
288 } // loop over threads
292 X_LIKWID_STOP("aa-even");
295 kd->PdfsActive = src;
296 KernelStatistics(kd, ld, cd, iter);
299 // Fixup bounce back PDFs.
301 #pragma omp parallel for default(none) \
304 #ifdef INTEL_OPT_DIRECTIVES
307 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
308 src[kd->BounceBackPdfsSrc[i]] = src[kd->BounceBackPdfsDst[i]];
311 // save current iteration
312 kda->Iteration = iter;
315 kd->PdfsActive = src;
316 KernelAddBodyForce(kd, ld, cd);
320 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
321 kd->PdfsActive = src;
322 VtkWrite(ld, kd, cd, iter);
327 kd->PdfsActive = src;
328 KernelStatistics(kd, ld, cd, iter);
331 // --------------------------------------------------------------------
335 X_LIKWID_START("aa-odd");
340 #pragma omp parallel for default(none) \
341 shared(gDims,src, w_0, w_1, w_2, omegaEven, omegaOdd, \
342 w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd, \
343 oX, oY, oZ, nX, nY, nZ, blk, nThreads) \
344 private(ux, uy, uz, ui, dens, dir_indep_trm, \
346 pdf_N, pdf_E, pdf_S, pdf_W, \
347 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
348 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
349 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
350 evenPart, oddPart, w_1_indep, w_2_indep)
353 for (int i = 0; i < nThreads; ++i) {
355 int threadStartX = nX / nThreads * i;
356 int threadEndX = nX / nThreads * (i + 1);
358 if (nX % nThreads > 0) {
359 if (nX % nThreads > i) {
364 threadStartX += nX % nThreads;
365 threadEndX += nX % nThreads;
369 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
370 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
371 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
373 // Must do everything here, else it would break collapse.
374 int eZ = MIN(bZ + blk[2], nZ + oZ);
375 int eY = MIN(bY + blk[1], nY + oY);
376 int eX = MIN(bX + blk[0], threadEndX + oX);
378 for (int x = bX; x < eX; ++x) {
379 for (int y = bY; y < eY; ++y) {
380 #ifdef INTEL_OPT_DIRECTIVES
382 #pragma vector always
385 for (int z = bZ; z < eZ; ++z) {
387 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
389 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
390 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x - _x, y - _y, z - _z, idxinv)];
395 // #define LID_DRIVEN_CAVITY
397 #ifdef LID_DRIVEN_CAVITY
399 if (z == nZ - 4 + oZ && x > 3 + oX && x < (nX - 4 + oX) && y > 3 + oY && y < (nY - 4 + oY)) {
400 ux = F(0.1) * F(0.5)77;
406 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
407 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
408 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
409 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
410 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
411 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
412 #ifdef LID_DRIVEN_CAVITY
417 pdf_N + pdf_E + pdf_S + pdf_W +
418 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
419 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
420 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
422 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0);
425 src[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
428 w_1_indep = w_1*dir_indep_trm;
431 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
432 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
433 src[I(x, y + 1, z, D3Q19_N)] = pdf_N - evenPart - oddPart;
434 src[I(x, y - 1, z, D3Q19_S)] = pdf_S - evenPart + oddPart;
437 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
438 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
439 src[I(x + 1, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart;
440 src[I(x - 1, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart;
443 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
444 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
445 src[I(x, y, z + 1, D3Q19_T)] = pdf_T - evenPart - oddPart;
446 src[I(x, y, z - 1, D3Q19_B)] = pdf_B - evenPart + oddPart;
449 w_2_indep = w_2*dir_indep_trm;
452 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
453 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
454 src[I(x - 1, y + 1, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
455 src[I(x + 1, y - 1, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
458 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
459 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
460 src[I(x + 1, y + 1, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
461 src[I(x - 1, y - 1, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
464 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
465 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
466 src[I(x - 1, y, z + 1, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
467 src[I(x + 1, y, z - 1, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
470 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
471 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
472 src[I(x + 1, y, z + 1, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
473 src[I(x - 1, y, z - 1, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
476 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
477 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
478 src[I(x, y - 1, z + 1, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
479 src[I(x, y + 1, z - 1, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
482 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
483 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
484 src[I(x, y + 1, z + 1, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
485 src[I(x, y - 1, z - 1, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
489 } } } // z, y, x (from inner to outer)
490 } } } // z, y, x (from inner to outer)
491 } // loop over threads
495 // Stop counters before bounce back. Else computing loop balance will be incorrect.
497 X_LIKWID_STOP("aa-odd");
499 // Fixup bounce back PDFs.
501 #pragma omp parallel for default(none) \
504 #ifdef INTEL_OPT_DIRECTIVES
507 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
508 src[kd->BounceBackPdfsDst[i]] = src[kd->BounceBackPdfsSrc[i]];
511 // save current iteration
512 kda->Iteration = iter + 1;
515 kd->PdfsActive = src;
516 KernelAddBodyForce(kd, ld, cd);
520 if (cd->VtkOutput && (iter + 1 % cd->VtkModulus) == 0) {
521 kd->PdfsActive = src;
522 VtkWrite(ld, kd, cd, iter + 1);
527 kd->PdfsActive = src;
528 KernelStatistics(kd, ld, cd, iter + 1);
532 } // for (int iter = 0; ...
537 kd->PdfsActive = src;
538 VtkWrite(ld, kd, cd, maxIterations);