projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
version 0.1
[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
42 Assert(cd->Omega > 0.0);
43 Assert(cd->Omega < 2.0);
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
54 PdfT magicParam = 1.0/12.0;
55 PdfT omegaOdd = 1.0/( 0.5 + magicParam/(1.0/omega - 0.5) );
56
57 PdfT evenPart = 0.0;
58 PdfT oddPart = 0.0;
59 PdfT dir_indep_trm = 0.0;
60
61 const PdfT w_0 = 1.0 / 3.0;
62 const PdfT w_1 = 1.0 / 18.0;
63 const PdfT w_2 = 1.0 / 36.0;
64
65 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;
66 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;
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
101 // TODO: outer openmp parallel
102 for(int iter = 0; iter < maxIterations; iter += 2) {
103
104
105 // --------------------------------------------------------------------
106 // even time step
107
108 X_LIKWID_START("list-aa-even");
109
110 #ifdef _OPENMP
111 #pragma omp parallel for default(none) \
112 shared(nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \
113 private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\
114 pdf_C, \
115 pdf_N, pdf_E, pdf_S, pdf_W, \
116 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
117 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
118 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW)
119 #endif
120 for (int index = 0; index < nFluid; ++index) {
121
122
123 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
124
125 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
126 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(index, idx)];
127 D3Q19_LIST
128 #undef X
129
130// #define LID_DRIVEN_CAVITY
131
132#ifdef LID_DRIVEN_CAVITY
133 int nX = kd->Dims[0];
134 int nY = kd->Dims[1];
135 int nZ = kd->Dims[2];
136
137 int x = kdl->Coords[C_INDEX_X(index)];
138 int y = kdl->Coords[C_INDEX_Y(index)];
139 int z = kdl->Coords[C_INDEX_Z(index)];
140
141 if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
142 ux = 0.1 * 0.577;
143 uy = 0.0;
144 uz = 0.0;
145 } else {
146#endif
147 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
148 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
149 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
150 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
151 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
152 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
153#ifdef LID_DRIVEN_CAVITY
154 }
155#endif
156
157 dens = pdf_C +
158 pdf_N + pdf_E + pdf_S + pdf_W +
159 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
160 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
161 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
162
163 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
164
165 // direction: w_0
166 src[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
167
168 // direction: w_1
169 w_1_indep = w_1*dir_indep_trm;
170
171 ui = uy;
172 evenPart = omegaEven*( 0.5*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
173 oddPart = omegaOdd*(0.5*(pdf_N - pdf_S) - ui*w_1_x3 );
174 src[I(index, D3Q19_S)] = pdf_N - evenPart - oddPart;
175 src[I(index, D3Q19_N)] = pdf_S - evenPart + oddPart;
176
177 ui = ux;
178 evenPart = omegaEven*( 0.5*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
179 oddPart = omegaOdd*(0.5*(pdf_E - pdf_W) - ui*w_1_x3 );
180 src[I(index, D3Q19_W)] = pdf_E - evenPart - oddPart;
181 src[I(index, D3Q19_E)] = pdf_W - evenPart + oddPart;
182
183 ui = uz;
184 evenPart = omegaEven*( 0.5*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
185 oddPart = omegaOdd*(0.5*(pdf_T - pdf_B) - ui*w_1_x3 );
186 src[I(index, D3Q19_B)] = pdf_T - evenPart - oddPart;
187 src[I(index, D3Q19_T)] = pdf_B - evenPart + oddPart;
188
189 // direction: w_2
190 w_2_indep = w_2*dir_indep_trm;
191
192 ui = -ux + uy;
193 evenPart = omegaEven*( 0.5*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
194 oddPart = omegaOdd*(0.5*(pdf_NW - pdf_SE) - ui*w_2_x3 );
195 src[I(index, D3Q19_SE)] = pdf_NW - evenPart - oddPart;
196 src[I(index, D3Q19_NW)] = pdf_SE - evenPart + oddPart;
197
198 ui = ux + uy;
199 evenPart = omegaEven*( 0.5*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
200 oddPart = omegaOdd*(0.5*(pdf_NE - pdf_SW) - ui*w_2_x3 );
201 src[I(index, D3Q19_SW)] = pdf_NE - evenPart - oddPart;
202 src[I(index, D3Q19_NE)] = pdf_SW - evenPart + oddPart;
203
204 ui = -ux + uz;
205 evenPart = omegaEven*( 0.5*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
206 oddPart = omegaOdd*(0.5*(pdf_TW - pdf_BE) - ui*w_2_x3 );
207 src[I(index, D3Q19_BE)] = pdf_TW - evenPart - oddPart;
208 src[I(index, D3Q19_TW)] = pdf_BE - evenPart + oddPart;
209
210 ui = ux + uz;
211 evenPart = omegaEven*( 0.5*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
212 oddPart = omegaOdd*(0.5*(pdf_TE - pdf_BW) - ui*w_2_x3 );
213 src[I(index, D3Q19_BW)] = pdf_TE - evenPart - oddPart;
214 src[I(index, D3Q19_TE)] = pdf_BW - evenPart + oddPart;
215
216 ui = -uy + uz;
217 evenPart = omegaEven*( 0.5*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
218 oddPart = omegaOdd*(0.5*(pdf_TS - pdf_BN) - ui*w_2_x3 );
219 src[I(index, D3Q19_BN)] = pdf_TS - evenPart - oddPart;
220 src[I(index, D3Q19_TS)] = pdf_BN - evenPart + oddPart;
221
222 ui = uy + uz;
223 evenPart = omegaEven*( 0.5*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
224 oddPart = omegaOdd*(0.5*(pdf_TN - pdf_BS) - ui*w_2_x3 );
225 src[I(index, D3Q19_BS)] = pdf_TN - evenPart - oddPart;
226 src[I(index, D3Q19_TN)] = pdf_BS - evenPart + oddPart;
227
228 } // loop over fluid nodes
229
230 X_LIKWID_STOP("list-aa-even");
231
232 // save current iteration
233 kdl->Iteration = iter;
234 #ifdef VERIFICATION
235 kd->PdfsActive = src;
236 KernelAddBodyForce(kd, ld, cd);
237 #endif
238
239 // --------------------------------------------------------------------
240 // odd time step
241
242 X_LIKWID_START("list-aa-odd");
243
244#ifdef _OPENMP
245 #pragma omp parallel for default(none) \
246 shared(nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \
247 private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\
248 pdf_C, \
249 pdf_N, pdf_E, pdf_S, pdf_W, \
250 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
251 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
252 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW)
253#endif
254 for (int index = 0; index < nFluid; ++index) {
255
256
257 adjListIndex = index * N_D3Q19_IDX;
258
259 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
260 pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)];
261
262 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[adjList[adjListIndex + idxinv]];
263 D3Q19_LIST_WO_C
264 #undef X
265
266#ifdef LID_DRIVEN_CAVITY
267 int nX = kd->Dims[0];
268 int nY = kd->Dims[1];
269 int nZ = kd->Dims[2];
270
271 int x = kdl->Coords[C_INDEX_X(index)];
272 int y = kdl->Coords[C_INDEX_Y(index)];
273 int z = kdl->Coords[C_INDEX_Z(index)];
274
275 if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
276 ux = 0.1 * 0.577;
277 uy = 0.0;
278 uz = 0.0;
279 } else {
280#endif
281 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
282 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
283 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
284 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
285 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
286 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
287#ifdef LID_DRIVEN_CAVITY
288 }
289#endif
290
291 dens = pdf_C +
292 pdf_N + pdf_E + pdf_S + pdf_W +
293 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
294 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
295 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
296
297 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0;
298
299 adjListIndex = index * N_D3Q19_IDX;
300
301 // direction: w_0
302 src[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
303
304 // direction: w_1
305 w_1_indep = w_1*dir_indep_trm;
306
307 ui = uy;
308 evenPart = omegaEven*( 0.5*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
309 oddPart = omegaOdd*(0.5*(pdf_N - pdf_S) - ui*w_1_x3 );
310 src[adjList[adjListIndex + D3Q19_N]] = pdf_N - evenPart - oddPart;
311 src[adjList[adjListIndex + D3Q19_S]] = pdf_S - evenPart + oddPart;
312
313 ui = ux;
314 evenPart = omegaEven*( 0.5*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
315 oddPart = omegaOdd*(0.5*(pdf_E - pdf_W) - ui*w_1_x3 );
316 src[adjList[adjListIndex + D3Q19_E]] = pdf_E - evenPart - oddPart;
317 src[adjList[adjListIndex + D3Q19_W]] = pdf_W - evenPart + oddPart;
318
319 ui = uz;
320 evenPart = omegaEven*( 0.5*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
321 oddPart = omegaOdd*(0.5*(pdf_T - pdf_B) - ui*w_1_x3 );
322 src[adjList[adjListIndex + D3Q19_T]] = pdf_T - evenPart - oddPart;
323 src[adjList[adjListIndex + D3Q19_B]] = pdf_B - evenPart + oddPart;
324
325 // direction: w_2
326 w_2_indep = w_2*dir_indep_trm;
327
328 ui = -ux + uy;
329 evenPart = omegaEven*( 0.5*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
330 oddPart = omegaOdd*(0.5*(pdf_NW - pdf_SE) - ui*w_2_x3 );
331 src[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart;
332 src[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart;
333
334 ui = ux + uy;
335 evenPart = omegaEven*( 0.5*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
336 oddPart = omegaOdd*(0.5*(pdf_NE - pdf_SW) - ui*w_2_x3 );
337 src[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart;
338 src[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart;
339
340 ui = -ux + uz;
341 evenPart = omegaEven*( 0.5*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
342 oddPart = omegaOdd*(0.5*(pdf_TW - pdf_BE) - ui*w_2_x3 );
343 src[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart;
344 src[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart;
345
346 ui = ux + uz;
347 evenPart = omegaEven*( 0.5*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
348 oddPart = omegaOdd*(0.5*(pdf_TE - pdf_BW) - ui*w_2_x3 );
349 src[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart;
350 src[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart;
351
352 ui = -uy + uz;
353 evenPart = omegaEven*( 0.5*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
354 oddPart = omegaOdd*(0.5*(pdf_TS - pdf_BN) - ui*w_2_x3 );
355 src[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart;
356 src[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart;
357
358 ui = uy + uz;
359 evenPart = omegaEven*( 0.5*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
360 oddPart = omegaOdd*(0.5*(pdf_TN - pdf_BS) - ui*w_2_x3 );
361 src[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart;
362 src[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart;
363
364 #undef I
365 } // loop over fluid nodes
366
367 X_LIKWID_STOP("list-aa-odd");
368
369 // save current iteration
370 kdl->Iteration = iter + 1;
371
372 #ifdef VERIFICATION
373 kd->PdfsActive = src;
374 KernelAddBodyForce(kd, ld, cd);
375 #endif
376
377 #ifdef VTK_OUTPUT
378 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
379 kd->PdfsActive = src;
380 VtkWrite(ld, kd, cd, iter);
381 }
382 #endif
383
384 #ifdef STATISTICS
385 kd->PdfsActive = src;
386 KernelStatistics(kd, ld, cd, iter);
387 #endif
388
389
390 } // for (int iter = 0; ...
391
392
393#ifdef VTK_OUTPUT
394 if (cd->VtkOutput) {
395 kd->PdfsActive = src;
396 VtkWrite(ld, kd, cd, maxIterations);
397 }
398#endif
399
400#ifdef STATISTICS
401 kd->PdfsActive = src;
402 KernelStatistics(kd, ld, cd, maxIterations);
403#endif
404
405 return;
406}
LBM benchmark kernels
This page took 0.102454 seconds and 5 git commands to generate.