Commit | Line | Data |
---|---|---|
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 "BenchKernelD3Q19ListCommon.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 | ||
37 | void FNAME(D3Q19ListKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd) | |
38 | { | |
39 | Assert(ld != NULL); | |
40 | Assert(kernelData != NULL); | |
41 | Assert(cd != NULL); | |
42 | ||
43 | Assert(cd->Omega > 0.0); | |
44 | Assert(cd->Omega < 2.0); | |
45 | ||
46 | KernelData * kd = (KernelData *)kernelData; | |
47 | KernelDataList * kdl = (KernelDataList *)kernelData; | |
48 | ||
49 | PdfT omega = cd->Omega; | |
50 | PdfT omegaEven = omega; | |
51 | // PdfT omegaOdd = 8.0*((2.0-omegaEven)/(8.0-omegaEven)); //"standard" trt odd relaxation parameter | |
52 | PdfT magicParam = 1.0/12.0; // 1/4: best stability; 1/12: removes third-order advection error (best advection); 1/6: removes fourth-order diffusion error (best diffusion); 3/16: exact location of bounce back for poiseuille flow | |
53 | PdfT omegaOdd = 1.0/( 0.5 + magicParam/(1.0/omega - 0.5) ); | |
54 | ||
55 | PdfT evenPart = 0.0; | |
56 | PdfT oddPart = 0.0; | |
57 | PdfT dir_indep_trm = 0.0; | |
58 | ||
59 | PdfT w_0 = 1.0 / 3.0; | |
60 | PdfT w_1 = 1.0 / 18.0; | |
61 | PdfT w_2 = 1.0 / 36.0; | |
62 | ||
63 | PdfT w_1_x3 = w_1 * 3.0; PdfT w_1_nine_half = w_1 * 9.0/2.0; PdfT w_1_indep = 0.0; | |
64 | PdfT w_2_x3 = w_2 * 3.0; PdfT w_2_nine_half = w_2 * 9.0/2.0; PdfT w_2_indep = 0.0; | |
65 | ||
66 | PdfT ux, uy, uz, ui; | |
67 | PdfT dens; | |
68 | ||
69 | // Declare pdf_N, pdf_E, pdf_S, pdf_W, ... | |
70 | #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name); | |
71 | D3Q19_LIST | |
72 | #undef X | |
73 | ||
74 | PdfT * src = kd->Pdfs[0]; | |
75 | PdfT * dst = kd->Pdfs[1]; | |
76 | PdfT * tmp; | |
77 | ||
78 | int maxIterations = cd->MaxIterations; | |
79 | int nFluid = kdl->nFluid; | |
80 | int nCells = kdl->nCells; | |
81 | ||
82 | uint32_t adjListIndex; | |
83 | uint32_t * adjList = kdl->AdjList; | |
84 | ||
85 | #ifdef VTK_OUTPUT | |
86 | if (cd->VtkOutput) { | |
87 | kd->PdfsActive = src; | |
88 | VtkWrite(ld, kd, cd, 0); | |
89 | } | |
90 | #endif | |
91 | ||
92 | #ifdef STATISTICS | |
93 | kd->PdfsActive = src; | |
94 | KernelStatistics(kd, ld, cd, 0); | |
95 | #endif | |
96 | ||
97 | // TODO: outer openmp parallel | |
98 | for(int iter = 0; iter < maxIterations; ++iter) { | |
99 | ||
100 | ||
e3f82424 | 101 | X_LIKWID_START("list-os"); |
10988083 MW |
102 | |
103 | #ifdef _OPENMP | |
104 | #pragma omp parallel for default(none) \ | |
105 | shared(nFluid, nCells, kd, kdl, adjList, src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \ | |
106 | w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd) \ | |
107 | private(ux, uy, uz, ui, dens, dir_indep_trm, adjListIndex, \ | |
108 | pdf_C, \ | |
109 | pdf_N, pdf_E, pdf_S, pdf_W, \ | |
110 | pdf_NE, pdf_SE, pdf_SW, pdf_NW, \ | |
111 | pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \ | |
112 | pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \ | |
113 | evenPart, oddPart, w_1_indep, w_2_indep) | |
114 | #endif | |
115 | for (int index = 0; index < nFluid; ++index) { | |
116 | ||
117 | #define I(index, dir) P_INDEX_3((nCells), (index), (dir)) | |
118 | ||
119 | #ifdef PROP_MODEL_PUSH | |
120 | ||
121 | // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ... | |
122 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(index, idx)]; | |
123 | D3Q19_LIST | |
124 | #undef X | |
125 | ||
126 | #elif PROP_MODEL_PULL | |
127 | ||
128 | adjListIndex = index * N_D3Q19_IDX; | |
129 | ||
130 | pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)]; | |
131 | ||
132 | // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_N]]; ... | |
133 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[adjList[adjListIndex + idx]]; | |
134 | D3Q19_LIST_WO_C | |
135 | #undef X | |
136 | ||
137 | #else | |
138 | #error No implementation for PROP_MODEL_NAME. | |
139 | #endif | |
140 | ||
141 | // #define LID_DRIVEN_CAVITY | |
142 | ||
143 | #ifdef LID_DRIVEN_CAVITY | |
144 | int nX = kd->Dims[0]; | |
145 | int nY = kd->Dims[1]; | |
146 | int nZ = kd->Dims[2]; | |
147 | ||
148 | int x = kdl->Coords[C_INDEX_X(index)]; | |
149 | int y = kdl->Coords[C_INDEX_Y(index)]; | |
150 | int z = kdl->Coords[C_INDEX_Z(index)]; | |
151 | ||
152 | if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) { | |
153 | ux = 0.1 * 0.577; | |
154 | uy = 0.0; | |
155 | uz = 0.0; | |
156 | } else { | |
157 | #endif | |
158 | ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - | |
159 | pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; | |
160 | uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN - | |
161 | pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS; | |
162 | uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS - | |
163 | pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS; | |
164 | #ifdef LID_DRIVEN_CAVITY | |
165 | } | |
166 | #endif | |
167 | ||
168 | dens = pdf_C + | |
169 | pdf_N + pdf_E + pdf_S + pdf_W + | |
170 | pdf_NE + pdf_SE + pdf_SW + pdf_NW + | |
171 | pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW + | |
172 | pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW; | |
173 | ||
174 | dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*3.0/2.0; | |
175 | ||
176 | #ifdef PROP_MODEL_PUSH | |
177 | ||
178 | adjListIndex = index * N_D3Q19_IDX; | |
179 | ||
180 | // direction: w_0 | |
181 | dst[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
182 | ||
183 | // direction: w_1 | |
184 | w_1_indep = w_1*dir_indep_trm; | |
185 | ||
186 | ui = uy; | |
187 | evenPart = omegaEven*( 0.5*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); | |
188 | oddPart = omegaOdd*(0.5*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
189 | dst[adjList[adjListIndex + D3Q19_N]] = pdf_N - evenPart - oddPart; | |
190 | dst[adjList[adjListIndex + D3Q19_S]] = pdf_S - evenPart + oddPart; | |
191 | ||
192 | ui = ux; | |
193 | evenPart = omegaEven*( 0.5*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); | |
194 | oddPart = omegaOdd*(0.5*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
195 | dst[adjList[adjListIndex + D3Q19_E]] = pdf_E - evenPart - oddPart; | |
196 | dst[adjList[adjListIndex + D3Q19_W]] = pdf_W - evenPart + oddPart; | |
197 | ||
198 | ui = uz; | |
199 | evenPart = omegaEven*( 0.5*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); | |
200 | oddPart = omegaOdd*(0.5*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
201 | dst[adjList[adjListIndex + D3Q19_T]] = pdf_T - evenPart - oddPart; | |
202 | dst[adjList[adjListIndex + D3Q19_B]] = pdf_B - evenPart + oddPart; | |
203 | ||
204 | // direction: w_2 | |
205 | w_2_indep = w_2*dir_indep_trm; | |
206 | ||
207 | ui = -ux + uy; | |
208 | evenPart = omegaEven*( 0.5*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); | |
209 | oddPart = omegaOdd*(0.5*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
210 | dst[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart; | |
211 | dst[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart; | |
212 | ||
213 | ui = ux + uy; | |
214 | evenPart = omegaEven*( 0.5*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); | |
215 | oddPart = omegaOdd*(0.5*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
216 | dst[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart; | |
217 | dst[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart; | |
218 | ||
219 | ui = -ux + uz; | |
220 | evenPart = omegaEven*( 0.5*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); | |
221 | oddPart = omegaOdd*(0.5*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
222 | dst[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart; | |
223 | dst[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart; | |
224 | ||
225 | ui = ux + uz; | |
226 | evenPart = omegaEven*( 0.5*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); | |
227 | oddPart = omegaOdd*(0.5*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
228 | dst[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart; | |
229 | dst[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart; | |
230 | ||
231 | ui = -uy + uz; | |
232 | evenPart = omegaEven*( 0.5*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); | |
233 | oddPart = omegaOdd*(0.5*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
234 | dst[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart; | |
235 | dst[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart; | |
236 | ||
237 | ui = uy + uz; | |
238 | evenPart = omegaEven*( 0.5*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); | |
239 | oddPart = omegaOdd*(0.5*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
240 | dst[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart; | |
241 | dst[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart; | |
242 | ||
243 | #elif PROP_MODEL_PULL | |
244 | ||
245 | // direction: w_0 | |
246 | dst[I(index, D3Q19_C )] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
247 | ||
248 | // direction: w_1 | |
249 | w_1_indep = w_1*dir_indep_trm; | |
250 | ||
251 | ui = uy; | |
252 | evenPart = omegaEven*( 0.5*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); | |
253 | oddPart = omegaOdd*(0.5*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
254 | dst[I(index, D3Q19_N )] = pdf_N - evenPart - oddPart; | |
255 | dst[I(index, D3Q19_S )] = pdf_S - evenPart + oddPart; | |
256 | ||
257 | ui = ux; | |
258 | evenPart = omegaEven*( 0.5*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); | |
259 | oddPart = omegaOdd*(0.5*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
260 | dst[I(index, D3Q19_E )] = pdf_E - evenPart - oddPart; | |
261 | dst[I(index, D3Q19_W )] = pdf_W - evenPart + oddPart; | |
262 | ||
263 | ui = uz; | |
264 | evenPart = omegaEven*( 0.5*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); | |
265 | oddPart = omegaOdd*(0.5*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
266 | dst[I(index, D3Q19_T )] = pdf_T - evenPart - oddPart; | |
267 | dst[I(index, D3Q19_B )] = pdf_B - evenPart + oddPart; | |
268 | ||
269 | // direction: w_2 | |
270 | w_2_indep = w_2*dir_indep_trm; | |
271 | ||
272 | ui = -ux + uy; | |
273 | evenPart = omegaEven*( 0.5*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); | |
274 | oddPart = omegaOdd*(0.5*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
275 | dst[I(index, D3Q19_NW)] = pdf_NW - evenPart - oddPart; | |
276 | dst[I(index, D3Q19_SE)] = pdf_SE - evenPart + oddPart; | |
277 | ||
278 | ui = ux + uy; | |
279 | evenPart = omegaEven*( 0.5*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); | |
280 | oddPart = omegaOdd*(0.5*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
281 | dst[I(index, D3Q19_NE)] = pdf_NE - evenPart - oddPart; | |
282 | dst[I(index, D3Q19_SW)] = pdf_SW - evenPart + oddPart; | |
283 | ||
284 | ui = -ux + uz; | |
285 | evenPart = omegaEven*( 0.5*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); | |
286 | oddPart = omegaOdd*(0.5*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
287 | dst[I(index, D3Q19_TW)] = pdf_TW - evenPart - oddPart; | |
288 | dst[I(index, D3Q19_BE)] = pdf_BE - evenPart + oddPart; | |
289 | ||
290 | ui = ux + uz; | |
291 | evenPart = omegaEven*( 0.5*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); | |
292 | oddPart = omegaOdd*(0.5*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
293 | dst[I(index, D3Q19_TE)] = pdf_TE - evenPart - oddPart; | |
294 | dst[I(index, D3Q19_BW)] = pdf_BW - evenPart + oddPart; | |
295 | ||
296 | ui = -uy + uz; | |
297 | evenPart = omegaEven*( 0.5*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); | |
298 | oddPart = omegaOdd*(0.5*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
299 | dst[I(index, D3Q19_TS)] = pdf_TS - evenPart - oddPart; | |
300 | dst[I(index, D3Q19_BN)] = pdf_BN - evenPart + oddPart; | |
301 | ||
302 | ui = uy + uz; | |
303 | evenPart = omegaEven*( 0.5*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); | |
304 | oddPart = omegaOdd*(0.5*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
305 | dst[I(index, D3Q19_TN)] = pdf_TN - evenPart - oddPart; | |
306 | dst[I(index, D3Q19_BS)] = pdf_BS - evenPart + oddPart; | |
307 | ||
308 | #endif | |
309 | #undef I | |
310 | } // loop over fluid nodes | |
311 | ||
e3f82424 MW |
312 | X_LIKWID_STOP("list-os"); |
313 | ||
10988083 MW |
314 | #ifdef VERIFICATION |
315 | kd->PdfsActive = dst; | |
316 | KernelAddBodyForce(kd, ld, cd); | |
317 | #endif | |
318 | ||
319 | #ifdef VTK_OUTPUT | |
320 | if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) { | |
321 | kd->PdfsActive = dst; | |
322 | VtkWrite(ld, kd, cd, iter); | |
323 | } | |
324 | #endif | |
325 | ||
326 | #ifdef STATISTICS | |
327 | kd->PdfsActive = dst; | |
328 | KernelStatistics(kd, ld, cd, iter); | |
329 | #endif | |
330 | ||
331 | // swap grids | |
332 | tmp = src; | |
333 | src = dst; | |
334 | dst = tmp; | |
335 | ||
336 | } // for (int iter = 0; ... | |
337 | ||
338 | #ifdef VTK_OUTPUT | |
339 | if (cd->VtkOutput) { | |
340 | kd->PdfsActive = src; | |
341 | VtkWrite(ld, kd, cd, maxIterations); | |
342 | } | |
343 | #endif | |
344 | ||
345 | #ifdef STATISTICS | |
346 | kd->PdfsActive = src; | |
347 | KernelStatistics(kd, ld, cd, maxIterations); | |
348 | #endif | |
349 | ||
350 | return; | |
351 | } |