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 "BenchKernelD3Q19ListCommon.h"
37 void FNAME(D3Q19ListKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
40 Assert(kernelData != NULL);
43 Assert(cd->Omega > F(0.0));
44 Assert(cd->Omega < F(2.0));
46 KernelData * kd = (KernelData *)kernelData;
47 KernelDataList * kdl = (KernelDataList *)kernelData;
49 PdfT omega = cd->Omega;
50 PdfT omegaEven = omega;
51 PdfT magicParam = F(1.0) / F(12.0);
52 PdfT omegaOdd = F(1.0) /(F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
54 PdfT evenPart = F(0.0);
55 PdfT oddPart = F(0.0);
56 PdfT dir_indep_trm = F(0.0);
58 PdfT w_0 = F(1.0) / F( 3.0);
59 PdfT w_1 = F(1.0) / F(18.0);
60 PdfT w_2 = F(1.0) / F(36.0);
62 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);
63 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);
68 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
69 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
73 PdfT * src = kd->Pdfs[0];
74 PdfT * dst = kd->Pdfs[1];
77 int maxIterations = cd->MaxIterations;
78 int nFluid = kdl->nFluid;
79 int nCells = kdl->nCells;
81 uint32_t adjListIndex;
82 uint32_t * adjList = kdl->AdjList;
87 VtkWrite(ld, kd, cd, 0);
93 KernelStatistics(kd, ld, cd, 0);
96 // TODO: outer openmp parallel
97 for(int iter = 0; iter < maxIterations; ++iter) {
100 X_LIKWID_START("list-os");
103 #pragma omp parallel for default(none) \
104 shared(nFluid, nCells, kd, kdl, adjList, src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \
105 w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd) \
106 private(ux, uy, uz, ui, dens, dir_indep_trm, adjListIndex, \
108 pdf_N, pdf_E, pdf_S, pdf_W, \
109 pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
110 pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
111 pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \
112 evenPart, oddPart, w_1_indep, w_2_indep)
114 #ifdef INTEL_OPT_DIRECTIVES
117 for (int index = 0; index < nFluid; ++index) {
119 #define I(index, dir) P_INDEX_3((nCells), (index), (dir))
121 #ifdef PROP_MODEL_PUSH
123 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
124 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(index, idx)];
128 #elif PROP_MODEL_PULL
130 adjListIndex = index * N_D3Q19_IDX;
132 pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)];
134 // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_N]]; ...
135 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[adjList[adjListIndex + idx]];
140 #error No implementation for PROP_MODEL_NAME.
143 // #define LID_DRIVEN_CAVITY
145 #ifdef LID_DRIVEN_CAVITY
146 int nX = kd->Dims[0];
147 int nY = kd->Dims[1];
148 int nZ = kd->Dims[2];
150 int x = kdl->Coords[C_INDEX_X(index)];
151 int y = kdl->Coords[C_INDEX_Y(index)];
152 int z = kdl->Coords[C_INDEX_Z(index)];
154 if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
160 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
161 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
162 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
163 pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
164 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
165 pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
166 #ifdef LID_DRIVEN_CAVITY
171 pdf_N + pdf_E + pdf_S + pdf_W +
172 pdf_NE + pdf_SE + pdf_SW + pdf_NW +
173 pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
174 pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
176 dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0);
178 #ifdef PROP_MODEL_PUSH
180 adjListIndex = index * N_D3Q19_IDX;
183 dst[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
186 w_1_indep = w_1*dir_indep_trm;
189 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
190 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
191 dst[adjList[adjListIndex + D3Q19_N]] = pdf_N - evenPart - oddPart;
192 dst[adjList[adjListIndex + D3Q19_S]] = pdf_S - evenPart + oddPart;
195 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
196 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
197 dst[adjList[adjListIndex + D3Q19_E]] = pdf_E - evenPart - oddPart;
198 dst[adjList[adjListIndex + D3Q19_W]] = pdf_W - evenPart + oddPart;
201 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
202 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
203 dst[adjList[adjListIndex + D3Q19_T]] = pdf_T - evenPart - oddPart;
204 dst[adjList[adjListIndex + D3Q19_B]] = pdf_B - evenPart + oddPart;
207 w_2_indep = w_2*dir_indep_trm;
210 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
211 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
212 dst[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart;
213 dst[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart;
216 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
217 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
218 dst[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart;
219 dst[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart;
222 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
223 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
224 dst[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart;
225 dst[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart;
228 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
229 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
230 dst[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart;
231 dst[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart;
234 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
235 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
236 dst[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart;
237 dst[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart;
240 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
241 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
242 dst[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart;
243 dst[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart;
245 #elif PROP_MODEL_PULL
248 dst[I(index, D3Q19_C )] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
251 w_1_indep = w_1*dir_indep_trm;
254 evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
255 oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
256 dst[I(index, D3Q19_N )] = pdf_N - evenPart - oddPart;
257 dst[I(index, D3Q19_S )] = pdf_S - evenPart + oddPart;
260 evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
261 oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
262 dst[I(index, D3Q19_E )] = pdf_E - evenPart - oddPart;
263 dst[I(index, D3Q19_W )] = pdf_W - evenPart + oddPart;
266 evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
267 oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
268 dst[I(index, D3Q19_T )] = pdf_T - evenPart - oddPart;
269 dst[I(index, D3Q19_B )] = pdf_B - evenPart + oddPart;
272 w_2_indep = w_2*dir_indep_trm;
275 evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
276 oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
277 dst[I(index, D3Q19_NW)] = pdf_NW - evenPart - oddPart;
278 dst[I(index, D3Q19_SE)] = pdf_SE - evenPart + oddPart;
281 evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
282 oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
283 dst[I(index, D3Q19_NE)] = pdf_NE - evenPart - oddPart;
284 dst[I(index, D3Q19_SW)] = pdf_SW - evenPart + oddPart;
287 evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
288 oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
289 dst[I(index, D3Q19_TW)] = pdf_TW - evenPart - oddPart;
290 dst[I(index, D3Q19_BE)] = pdf_BE - evenPart + oddPart;
293 evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
294 oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
295 dst[I(index, D3Q19_TE)] = pdf_TE - evenPart - oddPart;
296 dst[I(index, D3Q19_BW)] = pdf_BW - evenPart + oddPart;
299 evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
300 oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
301 dst[I(index, D3Q19_TS)] = pdf_TS - evenPart - oddPart;
302 dst[I(index, D3Q19_BN)] = pdf_BN - evenPart + oddPart;
305 evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
306 oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
307 dst[I(index, D3Q19_TN)] = pdf_TN - evenPart - oddPart;
308 dst[I(index, D3Q19_BS)] = pdf_BS - evenPart + oddPart;
312 } // loop over fluid nodes
314 X_LIKWID_STOP("list-os");
317 kd->PdfsActive = dst;
318 KernelAddBodyForce(kd, ld, cd);
322 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
323 kd->PdfsActive = dst;
324 VtkWrite(ld, kd, cd, iter);
329 kd->PdfsActive = dst;
330 KernelStatistics(kd, ld, cd, iter);
338 } // for (int iter = 0; ...
342 kd->PdfsActive = src;
343 VtkWrite(ld, kd, cd, maxIterations);
348 kd->PdfsActive = src;
349 KernelStatistics(kd, ld, cd, maxIterations);