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