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 "BenchKernelD3Q19Common.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 | #ifdef _OPENMP | |
37 | #include <omp.h> | |
38 | #endif | |
39 | ||
40 | void FNAME(D3Q19Kernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd) | |
41 | { | |
42 | Assert(ld != NULL); | |
43 | Assert(kernelData != NULL); | |
44 | Assert(cd != NULL); | |
45 | ||
0fde6e45 MW |
46 | Assert(cd->Omega > F(0.0)); |
47 | Assert(cd->Omega < F(2.0)); | |
10988083 MW |
48 | |
49 | KernelData * kd = (KernelData *)kernelData; | |
50 | ||
51 | ||
52 | int nX = ld->Dims[0]; | |
53 | int nY = ld->Dims[1]; | |
54 | int nZ = ld->Dims[2]; | |
55 | ||
56 | int * gDims = kd->GlobalDims; | |
57 | ||
58 | int oX = kd->Offsets[0]; | |
59 | int oY = kd->Offsets[1]; | |
60 | int oZ = kd->Offsets[2]; | |
61 | ||
62 | PdfT omega = cd->Omega; | |
63 | PdfT omegaEven = omega; | |
0fde6e45 MW |
64 | // PdfT omegaOdd = 8.0*((F(2.0)-omegaEven)/(8.0-omegaEven)); //"standard" trt odd relaxation parameter |
65 | PdfT magicParam = F(1.0) / F(12.0); | |
66 | // 1/ 4: best stability; | |
67 | // 1/12: removes third-order advection error (best advection); | |
68 | // 1/ 6: removes fourth-order diffusion error (best diffusion); | |
69 | // 3/16: exact location of bounce back for poiseuille flow | |
70 | PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5))); | |
10988083 | 71 | |
0fde6e45 MW |
72 | PdfT evenPart = F(0.0); |
73 | PdfT oddPart = F(0.0); | |
74 | PdfT dir_indep_trm = F(0.0); | |
10988083 | 75 | |
0fde6e45 MW |
76 | PdfT w_0 = F(1.0) / F( 3.0); |
77 | PdfT w_1 = F(1.0) / F(18.0); | |
78 | PdfT w_2 = F(1.0) / F(36.0); | |
10988083 | 79 | |
0fde6e45 MW |
80 | 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); |
81 | 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); | |
10988083 MW |
82 | |
83 | PdfT ux, uy, uz, ui; | |
84 | PdfT dens; | |
85 | ||
86 | // Declare pdf_N, pdf_E, pdf_S, pdf_W, ... | |
87 | #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name); | |
88 | D3Q19_LIST | |
89 | #undef X | |
90 | ||
91 | PdfT * src = kd->Pdfs[0]; | |
92 | PdfT * dst = kd->Pdfs[1]; | |
93 | PdfT * tmp; | |
94 | ||
95 | int maxIterations = cd->MaxIterations; | |
96 | ||
97 | #ifdef VTK_OUTPUT | |
98 | if (cd->VtkOutput) { | |
99 | kd->PdfsActive = src; | |
100 | VtkWrite(ld, kd, cd, 0); | |
101 | } | |
102 | #endif | |
103 | ||
8cafd9ea MW |
104 | X_KERNEL_START(kernelData); |
105 | ||
10988083 MW |
106 | for (int iter = 0; iter < maxIterations; ++iter) { |
107 | ||
e3f82424 MW |
108 | X_LIKWID_START("os"); |
109 | ||
10988083 | 110 | #ifdef _OPENMP |
0fde6e45 | 111 | #pragma omp parallel for collapse(2) default(none) \ |
10988083 MW |
112 | shared(gDims,src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \ |
113 | w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd, \ | |
114 | oX, oY, oZ, nX, nY, nZ) \ | |
115 | private(ux, uy, uz, ui, dens, dir_indep_trm, \ | |
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, \ | |
121 | evenPart, oddPart, w_1_indep, w_2_indep) | |
122 | #endif | |
e3f82424 | 123 | for (int x = oX; x < nX + oX; ++x) { |
10988083 | 124 | for (int y = oY; y < nY + oY; ++y) { |
0fde6e45 MW |
125 | #ifdef INTEL_OPT_DIRECTIVES |
126 | #pragma ivdep | |
127 | #pragma vector always | |
128 | #pragma simd | |
129 | #endif | |
8cafd9ea | 130 | for (int z = oZ; z < nZ + oZ; ++z) { // LOOP os |
10988083 MW |
131 | #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir)) |
132 | ||
133 | #ifdef PROP_MODEL_PUSH | |
134 | ||
135 | // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ... | |
136 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x, y, z, idx)]; | |
137 | //if (isnan(JOIN(pdf_,name))) { printf("iter: %d %d %d %d %d %s nan\n", iter, x-oX, y-oY, z-oZ, idx, D3Q19_NAMES[idx]); exit(1);} | |
138 | D3Q19_LIST | |
139 | #undef X | |
140 | ||
141 | #elif PROP_MODEL_PULL | |
142 | ||
143 | // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ... | |
144 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x - _x, y - _y, z - _z, idx)]; | |
145 | //if (isnan(JOIN(pdf_,name))) { printf("iter: %d %d %d %d %d %s nan\n", iter, x-oX, y-oY, z-oZ, idx, D3Q19_NAMES[idx]); exit(1);} | |
146 | D3Q19_LIST | |
147 | #undef X | |
148 | ||
149 | #else | |
150 | #error No implementation for PROP_MODEL_NAME. | |
151 | #endif | |
152 | ||
153 | // #define LID_DRIVEN_CAVITY | |
154 | ||
155 | #ifdef LID_DRIVEN_CAVITY | |
156 | ||
157 | if (z == nZ - 4 + oZ && x > 3 + oX && x < (nX - 4 + oX) && y > 3 + oY && y < (nY - 4 + oY)) { | |
0fde6e45 MW |
158 | ux = F(0.1 * 0.577); |
159 | uy = F(0.0); | |
160 | uz = F(0.0); | |
10988083 MW |
161 | |
162 | } else { | |
163 | #endif | |
164 | ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - | |
165 | pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; | |
166 | uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN - | |
167 | pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS; | |
168 | uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS - | |
169 | pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS; | |
170 | #ifdef LID_DRIVEN_CAVITY | |
171 | } | |
172 | ||
173 | #endif | |
174 | ||
175 | dens = pdf_C + | |
176 | pdf_N + pdf_E + pdf_S + pdf_W + | |
177 | pdf_NE + pdf_SE + pdf_SW + pdf_NW + | |
178 | pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW + | |
179 | pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW; | |
180 | ||
0fde6e45 | 181 | dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0); |
10988083 MW |
182 | |
183 | #ifdef PROP_MODEL_PUSH | |
184 | ||
185 | // direction: w_0 | |
186 | dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
187 | ||
188 | // direction: w_1 | |
189 | w_1_indep = w_1*dir_indep_trm; | |
190 | ||
191 | ui = uy; | |
0fde6e45 MW |
192 | evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); |
193 | oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
10988083 MW |
194 | dst[I(x, y + 1, z, D3Q19_N)] = pdf_N - evenPart - oddPart; |
195 | dst[I(x, y - 1, z, D3Q19_S)] = pdf_S - evenPart + oddPart; | |
196 | ||
197 | ui = ux; | |
0fde6e45 MW |
198 | evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); |
199 | oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
10988083 MW |
200 | dst[I(x + 1, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart; |
201 | dst[I(x - 1, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart; | |
202 | ||
203 | ui = uz; | |
0fde6e45 MW |
204 | evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); |
205 | oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
10988083 MW |
206 | dst[I(x, y, z + 1, D3Q19_T)] = pdf_T - evenPart - oddPart; |
207 | dst[I(x, y, z - 1, D3Q19_B)] = pdf_B - evenPart + oddPart; | |
208 | ||
209 | // direction: w_2 | |
210 | w_2_indep = w_2*dir_indep_trm; | |
211 | ||
212 | ui = -ux + uy; | |
0fde6e45 MW |
213 | evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); |
214 | oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
10988083 MW |
215 | dst[I(x - 1, y + 1, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart; |
216 | dst[I(x + 1, y - 1, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart; | |
217 | ||
218 | ui = ux + uy; | |
0fde6e45 MW |
219 | evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); |
220 | oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
10988083 MW |
221 | dst[I(x + 1, y + 1, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart; |
222 | dst[I(x - 1, y - 1, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart; | |
223 | ||
224 | ui = -ux + uz; | |
0fde6e45 MW |
225 | evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); |
226 | oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
10988083 MW |
227 | dst[I(x - 1, y, z + 1, D3Q19_TW)] = pdf_TW - evenPart - oddPart; |
228 | dst[I(x + 1, y, z - 1, D3Q19_BE)] = pdf_BE - evenPart + oddPart; | |
229 | ||
230 | ui = ux + uz; | |
0fde6e45 MW |
231 | evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); |
232 | oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
10988083 MW |
233 | dst[I(x + 1, y, z + 1, D3Q19_TE)] = pdf_TE - evenPart - oddPart; |
234 | dst[I(x - 1, y, z - 1, D3Q19_BW)] = pdf_BW - evenPart + oddPart; | |
235 | ||
236 | ui = -uy + uz; | |
0fde6e45 MW |
237 | evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); |
238 | oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
10988083 MW |
239 | dst[I(x, y - 1, z + 1, D3Q19_TS)] = pdf_TS - evenPart - oddPart; |
240 | dst[I(x, y + 1, z - 1, D3Q19_BN)] = pdf_BN - evenPart + oddPart; | |
241 | ||
242 | ui = uy + uz; | |
0fde6e45 MW |
243 | evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); |
244 | oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
10988083 MW |
245 | dst[I(x, y + 1, z + 1, D3Q19_TN)] = pdf_TN - evenPart - oddPart; |
246 | dst[I(x, y - 1, z - 1, D3Q19_BS)] = pdf_BS - evenPart + oddPart; | |
247 | ||
248 | #elif PROP_MODEL_PULL | |
249 | ||
250 | // direction: w_0 | |
251 | dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
252 | ||
253 | // direction: w_1 | |
254 | w_1_indep = w_1*dir_indep_trm; | |
255 | ||
256 | ui = uy; | |
0fde6e45 MW |
257 | evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); |
258 | oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
10988083 MW |
259 | dst[I(x, y, z, D3Q19_N)] = pdf_N - evenPart - oddPart; |
260 | dst[I(x, y, z, D3Q19_S)] = pdf_S - evenPart + oddPart; | |
261 | ||
262 | ui = ux; | |
0fde6e45 MW |
263 | evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); |
264 | oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
10988083 MW |
265 | dst[I(x, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart; |
266 | dst[I(x, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart; | |
267 | ||
268 | ui = uz; | |
0fde6e45 MW |
269 | evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); |
270 | oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
10988083 MW |
271 | dst[I(x, y, z, D3Q19_T)] = pdf_T - evenPart - oddPart; |
272 | dst[I(x, y, z, D3Q19_B)] = pdf_B - evenPart + oddPart; | |
273 | ||
274 | // direction: w_2 | |
275 | w_2_indep = w_2*dir_indep_trm; | |
276 | ||
277 | ui = -ux + uy; | |
0fde6e45 MW |
278 | evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); |
279 | oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
10988083 MW |
280 | dst[I(x, y, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart; |
281 | dst[I(x, y, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart; | |
282 | ||
283 | ui = ux + uy; | |
0fde6e45 MW |
284 | evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); |
285 | oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
10988083 MW |
286 | dst[I(x, y, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart; |
287 | dst[I(x, y, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart; | |
288 | ||
289 | ui = -ux + uz; | |
0fde6e45 MW |
290 | evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); |
291 | oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
10988083 MW |
292 | dst[I(x, y, z, D3Q19_TW)] = pdf_TW - evenPart - oddPart; |
293 | dst[I(x, y, z, D3Q19_BE)] = pdf_BE - evenPart + oddPart; | |
294 | ||
295 | ui = ux + uz; | |
0fde6e45 MW |
296 | evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); |
297 | oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
10988083 MW |
298 | dst[I(x, y, z, D3Q19_TE)] = pdf_TE - evenPart - oddPart; |
299 | dst[I(x, y, z, D3Q19_BW)] = pdf_BW - evenPart + oddPart; | |
300 | ||
301 | ui = -uy + uz; | |
0fde6e45 MW |
302 | evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); |
303 | oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
10988083 MW |
304 | dst[I(x, y, z, D3Q19_TS)] = pdf_TS - evenPart - oddPart; |
305 | dst[I(x, y, z, D3Q19_BN)] = pdf_BN - evenPart + oddPart; | |
306 | ||
307 | ui = uy + uz; | |
0fde6e45 MW |
308 | evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); |
309 | oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
10988083 MW |
310 | dst[I(x, y, z, D3Q19_TN)] = pdf_TN - evenPart - oddPart; |
311 | dst[I(x, y, z, D3Q19_BS)] = pdf_BS - evenPart + oddPart; | |
312 | ||
313 | #else | |
314 | #error No implementation for PROP_MODEL_NAME. | |
315 | #endif | |
316 | ||
317 | #undef I | |
318 | } | |
319 | } | |
320 | } // z, y, x (from inner to outer) | |
321 | ||
e3f82424 MW |
322 | // Stop counters before bounce back. Else computing loop balance will be incorrect. |
323 | X_LIKWID_STOP("os"); | |
324 | ||
8cafd9ea | 325 | |
10988083 MW |
326 | // Fixup bounce back PDFs. |
327 | #ifdef _OPENMP | |
328 | #pragma omp parallel for default(none) \ | |
329 | shared(kd, dst) | |
330 | #endif | |
331 | for (int i = 0; i < kd->nBounceBackPdfs; ++i) { | |
332 | dst[kd->BounceBackPdfsDst[i]] = dst[kd->BounceBackPdfsSrc[i]]; | |
333 | } | |
334 | ||
335 | #ifdef VERIFICATION | |
336 | kd->PdfsActive = dst; | |
337 | KernelAddBodyForce(kd, ld, cd); | |
338 | #endif | |
339 | ||
340 | #ifdef VTK_OUTPUT | |
341 | ||
342 | if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) { | |
343 | kd->PdfsActive = dst; | |
344 | VtkWrite(ld, kd, cd, iter); | |
345 | } | |
346 | ||
347 | #endif | |
348 | ||
349 | #ifdef STATISTICS | |
350 | kd->PdfsActive = dst; | |
351 | KernelStatistics(kd, ld, cd, iter); | |
352 | #endif | |
353 | ||
354 | // swap grids | |
355 | tmp = src; | |
356 | src = dst; | |
357 | dst = tmp; | |
358 | ||
359 | } // for (int iter = 0; ... | |
360 | ||
8cafd9ea MW |
361 | X_KERNEL_END(kernelData); |
362 | ||
10988083 MW |
363 | #ifdef VTK_OUTPUT |
364 | ||
365 | if (cd->VtkOutput) { | |
366 | kd->PdfsActive = src; | |
367 | VtkWrite(ld, kd, cd, maxIterations); | |
368 | } | |
369 | ||
370 | #endif | |
371 | ||
372 | return; | |
373 | } | |
374 | ||
375 | ||
376 | void FNAME(D3Q19BlkKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd) | |
377 | { | |
378 | Assert(ld != NULL); | |
379 | Assert(kernelData != NULL); | |
380 | Assert(cd != NULL); | |
381 | ||
0fde6e45 MW |
382 | Assert(cd->Omega > F(0.0)); |
383 | Assert(cd->Omega < F(2.0)); | |
10988083 MW |
384 | |
385 | KernelData * kd = (KernelData *)kernelData; | |
386 | ||
387 | ||
388 | int nX = ld->Dims[0]; | |
389 | int nY = ld->Dims[1]; | |
390 | int nZ = ld->Dims[2]; | |
391 | ||
392 | int * gDims = kd->GlobalDims; | |
393 | ||
394 | int oX = kd->Offsets[0]; | |
395 | int oY = kd->Offsets[1]; | |
396 | int oZ = kd->Offsets[2]; | |
397 | ||
398 | KernelDataEx * kdex = (KernelDataEx *)kd; | |
399 | ||
400 | int blk[3]; | |
401 | blk[0] = kdex->Blk[0]; | |
402 | blk[1] = kdex->Blk[1]; | |
403 | blk[2] = kdex->Blk[2]; | |
404 | ||
405 | PdfT omega = cd->Omega; | |
406 | PdfT omegaEven = omega; | |
0fde6e45 MW |
407 | // PdfT omegaOdd = 8.0*((F(2.0)-omegaEven)/(8.0-omegaEven)); //"standard" trt odd relaxation parameter |
408 | PdfT magicParam = F(1.0)/F(12.0); | |
409 | // 1/ 4: best stability; | |
410 | // 1/12: removes third-order advection error (best advection); | |
411 | // 1/ 6: removes fourth-order diffusion error (best diffusion); | |
412 | // 3/16: exact location of bounce back for poiseuille flow | |
413 | PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5))); | |
10988083 | 414 | |
0fde6e45 MW |
415 | PdfT evenPart = F(0.0); |
416 | PdfT oddPart = F(0.0); | |
417 | PdfT dir_indep_trm = F(0.0); | |
10988083 | 418 | |
0fde6e45 MW |
419 | PdfT w_0 = F(1.0) / F( 3.0); |
420 | PdfT w_1 = F(1.0) / F(18.0); | |
421 | PdfT w_2 = F(1.0) / F(36.0); | |
10988083 | 422 | |
0fde6e45 MW |
423 | 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); |
424 | 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); | |
10988083 MW |
425 | |
426 | PdfT ux, uy, uz, ui; | |
427 | PdfT dens; | |
428 | ||
429 | // Declare pdf_N, pdf_E, pdf_S, pdf_W, ... | |
430 | #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name); | |
431 | D3Q19_LIST | |
432 | #undef X | |
433 | ||
434 | PdfT * src = kd->Pdfs[0]; | |
435 | PdfT * dst = kd->Pdfs[1]; | |
436 | PdfT * tmp; | |
437 | ||
438 | int maxIterations = cd->MaxIterations; | |
439 | ||
440 | #ifdef VTK_OUTPUT | |
441 | if (cd->VtkOutput) { | |
442 | kd->PdfsActive = src; | |
443 | VtkWrite(ld, kd, cd, 0); | |
444 | } | |
445 | #endif | |
446 | ||
447 | int nThreads = 1; | |
448 | ||
449 | #ifdef _OPENMP | |
450 | nThreads = omp_get_max_threads(); | |
451 | #endif | |
452 | ||
8cafd9ea MW |
453 | X_KERNEL_START(kernelData); |
454 | ||
10988083 MW |
455 | for (int iter = 0; iter < maxIterations; ++iter) { |
456 | ||
10988083 | 457 | #ifdef _OPENMP |
e3f82424 | 458 | #pragma omp parallel default(none) \ |
10988083 MW |
459 | shared(gDims,src, dst, w_0, w_1, w_2, omegaEven, omegaOdd, \ |
460 | w_1_x3, w_2_x3, w_1_nine_half, w_2_nine_half, cd, \ | |
461 | oX, oY, oZ, nX, nY, nZ, blk, nThreads) \ | |
462 | private(ux, uy, uz, ui, dens, dir_indep_trm, \ | |
463 | pdf_C, \ | |
464 | pdf_N, pdf_E, pdf_S, pdf_W, \ | |
465 | pdf_NE, pdf_SE, pdf_SW, pdf_NW, \ | |
466 | pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \ | |
467 | pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \ | |
468 | evenPart, oddPart, w_1_indep, w_2_indep) | |
469 | #endif | |
e3f82424 MW |
470 | { |
471 | X_LIKWID_START("blk-os"); | |
10988083 | 472 | |
e3f82424 | 473 | int threadId = omp_get_thread_num(); |
10988083 | 474 | |
e3f82424 MW |
475 | int threadStartX = nX / nThreads * threadId; |
476 | int threadEndX = nX / nThreads * (threadId + 1); | |
10988083 MW |
477 | |
478 | if (nX % nThreads > 0) { | |
e3f82424 MW |
479 | if (nX % nThreads > threadId) { |
480 | threadStartX += threadId; | |
481 | threadEndX += threadId + 1; | |
10988083 MW |
482 | } |
483 | else { | |
484 | threadStartX += nX % nThreads; | |
485 | threadEndX += nX % nThreads; | |
486 | } | |
487 | } | |
488 | ||
e3f82424 | 489 | for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) { |
10988083 | 490 | for (int bY = oY; bY < nY + oY; bY += blk[1]) { |
e3f82424 | 491 | for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) { |
10988083 MW |
492 | |
493 | // Must do everything here, else it would break collapse. | |
494 | int eZ = MIN(bZ + blk[2], nZ + oZ); | |
495 | int eY = MIN(bY + blk[1], nY + oY); | |
496 | int eX = MIN(bX + blk[0], threadEndX + oX); | |
497 | ||
e3f82424 | 498 | for (int x = bX; x < eX; ++x) { |
10988083 | 499 | for (int y = bY; y < eY; ++y) { |
0fde6e45 MW |
500 | #ifdef INTEL_OPT_DIRECTIVES |
501 | #pragma ivdep | |
502 | #pragma vector always | |
503 | #pragma simd | |
504 | #endif | |
e3f82424 | 505 | for (int z = bZ; z < eZ; ++z) { |
10988083 MW |
506 | |
507 | #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir)) | |
508 | ||
509 | #ifdef PROP_MODEL_PUSH | |
510 | ||
511 | // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ... | |
512 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x, y, z, idx)]; | |
10988083 MW |
513 | D3Q19_LIST |
514 | #undef X | |
515 | ||
516 | #elif PROP_MODEL_PULL | |
517 | ||
518 | // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ... | |
519 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(x - _x, y - _y, z - _z, idx)]; | |
10988083 MW |
520 | D3Q19_LIST |
521 | #undef X | |
522 | ||
523 | #else | |
524 | #error No implementation for PROP_MODEL_NAME. | |
525 | #endif | |
526 | ||
527 | // #define LID_DRIVEN_CAVITY | |
528 | ||
529 | #ifdef LID_DRIVEN_CAVITY | |
530 | ||
531 | if (z == nZ - 4 + oZ && x > 3 + oX && x < (nX - 4 + oX) && y > 3 + oY && y < (nY - 4 + oY)) { | |
0fde6e45 | 532 | ux = 0.1 * F(0.5)77; |
10988083 MW |
533 | uy = 0.0; |
534 | uz = 0.0; | |
535 | ||
536 | } else { | |
537 | #endif | |
538 | ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - | |
e3f82424 | 539 | pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; |
10988083 | 540 | uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN - |
e3f82424 | 541 | pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS; |
10988083 | 542 | uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS - |
e3f82424 | 543 | pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS; |
10988083 MW |
544 | #ifdef LID_DRIVEN_CAVITY |
545 | } | |
546 | ||
547 | #endif | |
548 | ||
549 | dens = pdf_C + | |
550 | pdf_N + pdf_E + pdf_S + pdf_W + | |
551 | pdf_NE + pdf_SE + pdf_SW + pdf_NW + | |
552 | pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW + | |
553 | pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW; | |
554 | ||
0fde6e45 | 555 | dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz) * F(3.0) / F(2.0); |
10988083 MW |
556 | |
557 | #ifdef PROP_MODEL_PUSH | |
558 | ||
559 | // direction: w_0 | |
560 | dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
561 | ||
562 | // direction: w_1 | |
563 | w_1_indep = w_1*dir_indep_trm; | |
564 | ||
565 | ui = uy; | |
0fde6e45 MW |
566 | evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); |
567 | oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
10988083 MW |
568 | dst[I(x, y + 1, z, D3Q19_N)] = pdf_N - evenPart - oddPart; |
569 | dst[I(x, y - 1, z, D3Q19_S)] = pdf_S - evenPart + oddPart; | |
570 | ||
571 | ui = ux; | |
0fde6e45 MW |
572 | evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); |
573 | oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
10988083 MW |
574 | dst[I(x + 1, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart; |
575 | dst[I(x - 1, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart; | |
576 | ||
577 | ui = uz; | |
0fde6e45 MW |
578 | evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); |
579 | oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
10988083 MW |
580 | dst[I(x, y, z + 1, D3Q19_T)] = pdf_T - evenPart - oddPart; |
581 | dst[I(x, y, z - 1, D3Q19_B)] = pdf_B - evenPart + oddPart; | |
582 | ||
583 | // direction: w_2 | |
584 | w_2_indep = w_2*dir_indep_trm; | |
585 | ||
586 | ui = -ux + uy; | |
0fde6e45 MW |
587 | evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); |
588 | oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
10988083 MW |
589 | dst[I(x - 1, y + 1, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart; |
590 | dst[I(x + 1, y - 1, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart; | |
591 | ||
592 | ui = ux + uy; | |
0fde6e45 MW |
593 | evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); |
594 | oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
10988083 MW |
595 | dst[I(x + 1, y + 1, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart; |
596 | dst[I(x - 1, y - 1, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart; | |
597 | ||
598 | ui = -ux + uz; | |
0fde6e45 MW |
599 | evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); |
600 | oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
10988083 MW |
601 | dst[I(x - 1, y, z + 1, D3Q19_TW)] = pdf_TW - evenPart - oddPart; |
602 | dst[I(x + 1, y, z - 1, D3Q19_BE)] = pdf_BE - evenPart + oddPart; | |
603 | ||
604 | ui = ux + uz; | |
0fde6e45 MW |
605 | evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); |
606 | oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
10988083 MW |
607 | dst[I(x + 1, y, z + 1, D3Q19_TE)] = pdf_TE - evenPart - oddPart; |
608 | dst[I(x - 1, y, z - 1, D3Q19_BW)] = pdf_BW - evenPart + oddPart; | |
609 | ||
610 | ui = -uy + uz; | |
0fde6e45 MW |
611 | evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); |
612 | oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
10988083 MW |
613 | dst[I(x, y - 1, z + 1, D3Q19_TS)] = pdf_TS - evenPart - oddPart; |
614 | dst[I(x, y + 1, z - 1, D3Q19_BN)] = pdf_BN - evenPart + oddPart; | |
615 | ||
616 | ui = uy + uz; | |
0fde6e45 MW |
617 | evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); |
618 | oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
10988083 MW |
619 | dst[I(x, y + 1, z + 1, D3Q19_TN)] = pdf_TN - evenPart - oddPart; |
620 | dst[I(x, y - 1, z - 1, D3Q19_BS)] = pdf_BS - evenPart + oddPart; | |
621 | ||
622 | #elif PROP_MODEL_PULL | |
623 | ||
624 | // direction: w_0 | |
625 | dst[I(x, y, z, D3Q19_C)] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
626 | ||
627 | // direction: w_1 | |
628 | w_1_indep = w_1*dir_indep_trm; | |
629 | ||
630 | ui = uy; | |
0fde6e45 MW |
631 | evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); |
632 | oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
10988083 MW |
633 | dst[I(x, y, z, D3Q19_N)] = pdf_N - evenPart - oddPart; |
634 | dst[I(x, y, z, D3Q19_S)] = pdf_S - evenPart + oddPart; | |
635 | ||
636 | ui = ux; | |
0fde6e45 MW |
637 | evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); |
638 | oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
10988083 MW |
639 | dst[I(x, y, z, D3Q19_E)] = pdf_E - evenPart - oddPart; |
640 | dst[I(x, y, z, D3Q19_W)] = pdf_W - evenPart + oddPart; | |
641 | ||
642 | ui = uz; | |
0fde6e45 MW |
643 | evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); |
644 | oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
10988083 MW |
645 | dst[I(x, y, z, D3Q19_T)] = pdf_T - evenPart - oddPart; |
646 | dst[I(x, y, z, D3Q19_B)] = pdf_B - evenPart + oddPart; | |
647 | ||
648 | // direction: w_2 | |
649 | w_2_indep = w_2*dir_indep_trm; | |
650 | ||
651 | ui = -ux + uy; | |
0fde6e45 MW |
652 | evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); |
653 | oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
10988083 MW |
654 | dst[I(x, y, z, D3Q19_NW)] = pdf_NW - evenPart - oddPart; |
655 | dst[I(x, y, z, D3Q19_SE)] = pdf_SE - evenPart + oddPart; | |
656 | ||
657 | ui = ux + uy; | |
0fde6e45 MW |
658 | evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); |
659 | oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
10988083 MW |
660 | dst[I(x, y, z, D3Q19_NE)] = pdf_NE - evenPart - oddPart; |
661 | dst[I(x, y, z, D3Q19_SW)] = pdf_SW - evenPart + oddPart; | |
662 | ||
663 | ui = -ux + uz; | |
0fde6e45 MW |
664 | evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); |
665 | oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
10988083 MW |
666 | dst[I(x, y, z, D3Q19_TW)] = pdf_TW - evenPart - oddPart; |
667 | dst[I(x, y, z, D3Q19_BE)] = pdf_BE - evenPart + oddPart; | |
668 | ||
669 | ui = ux + uz; | |
0fde6e45 MW |
670 | evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); |
671 | oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
10988083 MW |
672 | dst[I(x, y, z, D3Q19_TE)] = pdf_TE - evenPart - oddPart; |
673 | dst[I(x, y, z, D3Q19_BW)] = pdf_BW - evenPart + oddPart; | |
674 | ||
675 | ui = -uy + uz; | |
0fde6e45 MW |
676 | evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); |
677 | oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
10988083 MW |
678 | dst[I(x, y, z, D3Q19_TS)] = pdf_TS - evenPart - oddPart; |
679 | dst[I(x, y, z, D3Q19_BN)] = pdf_BN - evenPart + oddPart; | |
680 | ||
681 | ui = uy + uz; | |
0fde6e45 MW |
682 | evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); |
683 | oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
10988083 MW |
684 | dst[I(x, y, z, D3Q19_TN)] = pdf_TN - evenPart - oddPart; |
685 | dst[I(x, y, z, D3Q19_BS)] = pdf_BS - evenPart + oddPart; | |
686 | ||
687 | #else | |
688 | #error No implementation for PROP_MODEL_NAME. | |
689 | #endif | |
690 | ||
691 | #undef I | |
692 | } | |
693 | } | |
694 | } // z, y, x (from inner to outer) | |
695 | } | |
696 | } | |
697 | } // z, y, x (from inner to outer) | |
698 | ||
e3f82424 MW |
699 | X_LIKWID_STOP("blk-os"); |
700 | } // parallel region | |
701 | ||
702 | // Stop counters before bounce back. Else computing loop balance will be incorrect. | |
10988083 MW |
703 | |
704 | // Fixup bounce back PDFs. | |
705 | #ifdef _OPENMP | |
706 | #pragma omp parallel for default(none) \ | |
707 | shared(kd, dst) | |
708 | #endif | |
709 | for (int i = 0; i < kd->nBounceBackPdfs; ++i) { | |
710 | dst[kd->BounceBackPdfsDst[i]] = dst[kd->BounceBackPdfsSrc[i]]; | |
711 | } | |
712 | ||
713 | #ifdef VERIFICATION | |
714 | kd->PdfsActive = dst; | |
715 | KernelAddBodyForce(kd, ld, cd); | |
716 | #endif | |
717 | ||
718 | #ifdef VTK_OUTPUT | |
719 | ||
720 | if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) { | |
721 | kd->PdfsActive = dst; | |
722 | VtkWrite(ld, kd, cd, iter); | |
723 | } | |
724 | ||
725 | #endif | |
726 | ||
727 | #ifdef STATISTICS | |
728 | kd->PdfsActive = dst; | |
729 | KernelStatistics(kd, ld, cd, iter); | |
730 | #endif | |
731 | ||
732 | // swap grids | |
733 | tmp = src; | |
734 | src = dst; | |
735 | dst = tmp; | |
736 | ||
737 | } // for (int iter = 0; ... | |
738 | ||
8cafd9ea MW |
739 | X_KERNEL_END(kernelData); |
740 | ||
10988083 MW |
741 | #ifdef VTK_OUTPUT |
742 | ||
743 | if (cd->VtkOutput) { | |
744 | kd->PdfsActive = src; | |
745 | VtkWrite(ld, kd, cd, maxIterations); | |
746 | } | |
747 | ||
748 | #endif | |
749 | ||
750 | return; | |
751 | } | |
752 |