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 / BenchKernelD3Q19ListAa.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 "BenchKernelD3Q19ListAaCommon.h"
28
29#include "Memory.h"
30#include "Vtk.h"
31#include "LikwidIf.h"
32
33#include <inttypes.h>
34#include <math.h>
35
36void FNAME(D3Q19ListAaKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
37{
38 Assert(ld != NULL);
39 Assert(kernelData != NULL);
40 Assert(cd != NULL);
41
0fde6e45
MW		 42 Assert(cd->Omega > F(0.0));
43 Assert(cd->Omega < F(2.0));
10988083
MW		 44
45 KernelData * kd = (KernelData *)kernelData;
46 KernelDataList * kdl = (KernelDataList *)kernelData;
47
48 PdfT omega = cd->Omega;
49 PdfT omegaEven = omega;
50 // 1/4: best stability;
51 // 1/12: removes third-order advection error (best advection);
52 // 1/6: removes fourth-order diffusion error (best diffusion);
53 // 3/16: exact location of bounce back for poiseuille flow
0fde6e45
MW		 54 PdfT magicParam = F(1.0)/F(12.0);
55 PdfT omegaOdd = F(1.0)/( F(0.5) + magicParam/(F(1.0)/omega - F(0.5)) );
10988083 56
0fde6e45
MW		 57 PdfT evenPart = F(0.0);
58 PdfT oddPart = F(0.0);
59 PdfT dir_indep_trm = F(0.0);
10988083 60
0fde6e45
MW		 61 const PdfT w_0 = F(1.0) / F(3.0);
62 const PdfT w_1 = F(1.0) / F(18.0);
63 const PdfT w_2 = F(1.0) / F(36.0);
10988083 64
0fde6e45
MW		 65 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0)/F(2.0); PdfT w_1_indep = F(0.0);
66 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0)/F(2.0); PdfT w_2_indep = F(0.0);
10988083
MW		 67
68 PdfT ui;
69
70 PdfT ux, uy, uz;
71 PdfT dens;
72
73 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
74 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
75 D3Q19_LIST
76 #undef X
77
78 PdfT * src = kd->Pdfs[0];
79
80 int maxIterations = cd->MaxIterations;
81 int nFluid = kdl->nFluid;
82 int nCells = kdl->nCells;
83
84 uint32_t adjListIndex;
85 uint32_t * adjList = kdl->AdjList;
86
87
88
89 #ifdef VTK_OUTPUT
90 if (cd->VtkOutput) {
91 kd->PdfsActive = src;
92 VtkWrite(ld, kd, cd, -1);
93 }
94 #endif
95
96 #ifdef STATISTICS
97 kd->PdfsActive = src;
98 KernelStatistics(kd, ld, cd, 0);
99 #endif
100
8cafd9ea
MW		 101 X_KERNEL_START(kernelData);
102
10988083
MW		 103 // TODO: outer openmp parallel
104 for(int iter = 0; iter < maxIterations; iter += 2) {
105
106
107 // --------------------------------------------------------------------
108 // even time step
109
110 X_LIKWID_START("list-aa-even");
111
0fde6e45 112#ifdef _OPENMP
10988083
MW		 113 #pragma omp parallel for default(none) \
114 shared(nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \
115 private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\
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)
0fde6e45
MW		 121#endif
122#ifdef INTEL_OPT_DIRECTIVES
123 #pragma ivdep
124 #pragma vector always
125 #pragma simd
126#endif
8cafd9ea 127 for (int index = 0; index < nFluid; ++index) { // LOOP list-aa-even
10988083
MW		 128
129
130 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
131
132 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
133 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(index, idx)];
134 D3Q19_LIST
135 #undef X
136
137// #define LID_DRIVEN_CAVITY
138
139#ifdef LID_DRIVEN_CAVITY
140 int nX = kd->Dims[0];
141 int nY = kd->Dims[1];
142 int nZ = kd->Dims[2];
143
144 int x = kdl->Coords[C_INDEX_X(index)];
145 int y = kdl->Coords[C_INDEX_Y(index)];
146 int z = kdl->Coords[C_INDEX_Z(index)];
147
148 if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
149 ux = 0.1 * 0.577;
150 uy = 0.0;
151 uz = 0.0;
152 } else {
153#endif
154 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
155 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
156 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
157 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
158 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
159 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
160#ifdef LID_DRIVEN_CAVITY
161 }
162#endif
163
164 dens = pdf_C +
165 pdf_N + pdf_E + pdf_S + pdf_W +
166 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
167 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
168 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
169
0fde6e45 170 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0);
10988083
MW		 171
172 // direction: w_0
173 src[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
174
175 // direction: w_1
176 w_1_indep = w_1*dir_indep_trm;
177
178 ui = uy;
0fde6e45
MW		 179 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
180 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
10988083
MW		 181 src[I(index, D3Q19_S)] = pdf_N - evenPart - oddPart;
182 src[I(index, D3Q19_N)] = pdf_S - evenPart + oddPart;
183
184 ui = ux;
0fde6e45
MW		 185 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
186 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
10988083
MW		 187 src[I(index, D3Q19_W)] = pdf_E - evenPart - oddPart;
188 src[I(index, D3Q19_E)] = pdf_W - evenPart + oddPart;
189
190 ui = uz;
0fde6e45
MW		 191 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
192 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
10988083
MW		 193 src[I(index, D3Q19_B)] = pdf_T - evenPart - oddPart;
194 src[I(index, D3Q19_T)] = pdf_B - evenPart + oddPart;
195
196 // direction: w_2
197 w_2_indep = w_2*dir_indep_trm;
198
199 ui = -ux + uy;
0fde6e45
MW		 200 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
201 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083
MW		 202 src[I(index, D3Q19_SE)] = pdf_NW - evenPart - oddPart;
203 src[I(index, D3Q19_NW)] = pdf_SE - evenPart + oddPart;
204
205 ui = ux + uy;
0fde6e45
MW		 206 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
207 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083
MW		 208 src[I(index, D3Q19_SW)] = pdf_NE - evenPart - oddPart;
209 src[I(index, D3Q19_NE)] = pdf_SW - evenPart + oddPart;
210
211 ui = -ux + uz;
0fde6e45
MW		 212 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
213 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083
MW		 214 src[I(index, D3Q19_BE)] = pdf_TW - evenPart - oddPart;
215 src[I(index, D3Q19_TW)] = pdf_BE - evenPart + oddPart;
216
217 ui = ux + uz;
0fde6e45
MW		 218 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
219 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083
MW		 220 src[I(index, D3Q19_BW)] = pdf_TE - evenPart - oddPart;
221 src[I(index, D3Q19_TE)] = pdf_BW - evenPart + oddPart;
222
223 ui = -uy + uz;
0fde6e45
MW		 224 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
225 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083
MW		 226 src[I(index, D3Q19_BN)] = pdf_TS - evenPart - oddPart;
227 src[I(index, D3Q19_TS)] = pdf_BN - evenPart + oddPart;
228
229 ui = uy + uz;
0fde6e45
MW		 230 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
231 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083
MW		 232 src[I(index, D3Q19_BS)] = pdf_TN - evenPart - oddPart;
233 src[I(index, D3Q19_TN)] = pdf_BS - evenPart + oddPart;
234
235 } // loop over fluid nodes
236
237 X_LIKWID_STOP("list-aa-even");
238
239 // save current iteration
240 kdl->Iteration = iter;
241 #ifdef VERIFICATION
242 kd->PdfsActive = src;
243 KernelAddBodyForce(kd, ld, cd);
244 #endif
245
246 // --------------------------------------------------------------------
247 // odd time step
248
249 X_LIKWID_START("list-aa-odd");
250
251#ifdef _OPENMP
252 #pragma omp parallel for default(none) \
253 shared(nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \
254 private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\
255 pdf_C, \
256 pdf_N, pdf_E, pdf_S, pdf_W, \
257 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
258 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
259 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW)
0fde6e45
MW		 260#endif
261#ifdef INTEL_OPT_DIRECTIVES
262 #pragma ivdep
10988083 263#endif
8cafd9ea 264 for (int index = 0; index < nFluid; ++index) { // LOOP list-aa-odd
10988083
MW		 265
266
267 adjListIndex = index * N_D3Q19_IDX;
268
269 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
270 pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)];
271
272 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[adjList[adjListIndex + idxinv]];
273 D3Q19_LIST_WO_C
274 #undef X
275
276#ifdef LID_DRIVEN_CAVITY
277 int nX = kd->Dims[0];
278 int nY = kd->Dims[1];
279 int nZ = kd->Dims[2];
280
281 int x = kdl->Coords[C_INDEX_X(index)];
282 int y = kdl->Coords[C_INDEX_Y(index)];
283 int z = kdl->Coords[C_INDEX_Z(index)];
284
285 if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
0fde6e45
MW		 286 ux = F(0.1) * F(0.5)77;
287 uy = F(0.0);
288 uz = F(0.0);
10988083
MW		 289 } else {
290#endif
291 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
292 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
293 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
294 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
295 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
296 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
297#ifdef LID_DRIVEN_CAVITY
298 }
299#endif
300
301 dens = pdf_C +
302 pdf_N + pdf_E + pdf_S + pdf_W +
303 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
304 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
305 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
306
0fde6e45 307 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0);
10988083
MW		 308
309 adjListIndex = index * N_D3Q19_IDX;
310
311 // direction: w_0
312 src[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
313
314 // direction: w_1
315 w_1_indep = w_1*dir_indep_trm;
316
317 ui = uy;
0fde6e45
MW		 318 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
319 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
10988083
MW		 320 src[adjList[adjListIndex + D3Q19_N]] = pdf_N - evenPart - oddPart;
321 src[adjList[adjListIndex + D3Q19_S]] = pdf_S - evenPart + oddPart;
322
323 ui = ux;
0fde6e45
MW		 324 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
325 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
10988083
MW		 326 src[adjList[adjListIndex + D3Q19_E]] = pdf_E - evenPart - oddPart;
327 src[adjList[adjListIndex + D3Q19_W]] = pdf_W - evenPart + oddPart;
328
329 ui = uz;
0fde6e45
MW		 330 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
331 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
10988083
MW		 332 src[adjList[adjListIndex + D3Q19_T]] = pdf_T - evenPart - oddPart;
333 src[adjList[adjListIndex + D3Q19_B]] = pdf_B - evenPart + oddPart;
334
335 // direction: w_2
336 w_2_indep = w_2*dir_indep_trm;
337
338 ui = -ux + uy;
0fde6e45
MW		 339 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
340 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
10988083
MW		 341 src[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart;
342 src[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart;
343
344 ui = ux + uy;
0fde6e45
MW		 345 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
346 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
10988083
MW		 347 src[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart;
348 src[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart;
349
350 ui = -ux + uz;
0fde6e45
MW		 351 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
352 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
10988083
MW		 353 src[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart;
354 src[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart;
355
356 ui = ux + uz;
0fde6e45
MW		 357 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
358 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
10988083
MW		 359 src[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart;
360 src[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart;
361
362 ui = -uy + uz;
0fde6e45
MW		 363 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
364 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
10988083
MW		 365 src[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart;
366 src[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart;
367
368 ui = uy + uz;
0fde6e45
MW		 369 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
370 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
10988083
MW		 371 src[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart;
372 src[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart;
373
374 #undef I
375 } // loop over fluid nodes
376
377 X_LIKWID_STOP("list-aa-odd");
378
379 // save current iteration
380 kdl->Iteration = iter + 1;
381
382 #ifdef VERIFICATION
383 kd->PdfsActive = src;
384 KernelAddBodyForce(kd, ld, cd);
385 #endif
386
387 #ifdef VTK_OUTPUT
388 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
389 kd->PdfsActive = src;
390 VtkWrite(ld, kd, cd, iter);
391 }
392 #endif
393
394 #ifdef STATISTICS
395 kd->PdfsActive = src;
396 KernelStatistics(kd, ld, cd, iter);
397 #endif
398
399
400 } // for (int iter = 0; ...
401
8cafd9ea
MW		 402 X_KERNEL_END(kernelData);
403
10988083
MW		 404
405#ifdef VTK_OUTPUT
406 if (cd->VtkOutput) {
407 kd->PdfsActive = src;
408 VtkWrite(ld, kd, cd, maxIterations);
409 }
410#endif
411
412#ifdef STATISTICS
413 kd->PdfsActive = src;
414 KernelStatistics(kd, ld, cd, maxIterations);
415#endif
416
417 return;
418}
LBM benchmark kernels
This page took 0.069637 seconds and 5 git commands to generate.