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 "BenchKernelD3Q19ListAaRiaCommon.h" | |
28 | ||
29 | #include "Memory.h" | |
30 | #include "Vtk.h" | |
31 | #include "LikwidIf.h" | |
32 | ||
33 | #include <inttypes.h> | |
34 | #include <math.h> | |
35 | ||
36 | #ifdef _OPENMP | |
37 | #include <omp.h> | |
38 | #endif | |
39 | ||
40 | void FNAME(D3Q19ListAaRiaKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd) | |
41 | { | |
42 | Assert(ld != NULL); | |
43 | Assert(kernelData != NULL); | |
44 | Assert(cd != NULL); | |
45 | ||
46 | Assert(cd->Omega > 0.0); | |
47 | Assert(cd->Omega < 2.0); | |
48 | ||
49 | KernelData * kd = (KernelData *)kernelData; | |
50 | KernelDataList * kdl = KDL(kernelData); | |
51 | KernelDataListRia * kdlr = KDLR(kernelData); | |
52 | PdfT omega = cd->Omega; | |
53 | PdfT omegaEven = omega; | |
54 | // 1/ 4: best stability; | |
55 | // 1/12: removes third-order advection error (best advection); | |
56 | // 1/ 6: removes fourth-order diffusion error (best diffusion); | |
57 | // 3/16: exact location of bounce back for poiseuille flow | |
0fde6e45 MW |
58 | PdfT magicParam = F(1.0) / F(12.0); |
59 | PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5))); | |
10988083 | 60 | |
0fde6e45 MW |
61 | PdfT evenPart = F(0.0); |
62 | PdfT oddPart = F(0.0); | |
63 | PdfT dir_indep_trm = F(0.0); | |
10988083 | 64 | |
0fde6e45 MW |
65 | const PdfT w_0 = F(1.0) / F(3.0); |
66 | const PdfT w_1 = F(1.0) / F(18.0); | |
67 | const PdfT w_2 = F(1.0) / F(36.0); | |
10988083 | 68 | |
0fde6e45 MW |
69 | const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0); PdfT w_1_indep = F(0.0); |
70 | const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0); PdfT w_2_indep = F(0.0); | |
10988083 MW |
71 | |
72 | PdfT ui; | |
73 | ||
74 | PdfT ux, uy, uz; | |
75 | PdfT dens; | |
76 | ||
77 | // Declare pdf_N, pdf_E, pdf_S, pdf_W, ... | |
78 | #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name); | |
79 | D3Q19_LIST | |
80 | #undef X | |
81 | ||
82 | // Declare pointers to pdfs ppdf_N, ppdf_E, ppdf_S, ppdf_W, ... | |
83 | #define X(name, idx, idxinv, x, y, z) PdfT * JOIN(ppdf_,name) = NULL; | |
84 | D3Q19_LIST | |
85 | #undef X | |
86 | ||
87 | uint32_t nConsecNodes = kdlr->nConsecNodes; | |
88 | uint32_t * consecNodes = kdlr->ConsecNodes; | |
89 | uint32_t consecIndex = 0; | |
90 | uint32_t consecValue = 0; | |
91 | ||
92 | #ifndef DEBUG | |
93 | UNUSED(nConsecNodes); | |
94 | #endif | |
95 | ||
96 | PdfT * src = kd->Pdfs[0]; | |
97 | ||
98 | int maxIterations = cd->MaxIterations; | |
99 | int nFluid = kdl->nFluid; | |
100 | int nCells = kdl->nCells; | |
101 | ||
102 | uint32_t adjListIndex; | |
103 | uint32_t * adjList = kdl->AdjList; | |
104 | ||
105 | ||
106 | ||
107 | #ifdef VTK_OUTPUT | |
108 | if (cd->VtkOutput) { | |
109 | kd->PdfsActive = src; | |
110 | VtkWrite(ld, kd, cd, -1); | |
111 | } | |
112 | #endif | |
113 | ||
114 | #ifdef STATISTICS | |
115 | kd->PdfsActive = src; | |
116 | KernelStatistics(kd, ld, cd, 0); | |
117 | #endif | |
118 | ||
8cafd9ea MW |
119 | X_KERNEL_START(kernelData); |
120 | ||
10988083 MW |
121 | // TODO: outer openmp parallel |
122 | for(int iter = 0; iter < maxIterations; iter += 2) { | |
123 | ||
124 | X_LIKWID_START("list-aa-ria-even"); | |
125 | ||
126 | // -------------------------------------------------------------------- | |
127 | // even time step | |
128 | // -------------------------------------------------------------------- | |
129 | #ifdef _OPENMP | |
130 | #pragma omp parallel for default(none) \ | |
131 | shared(stderr, nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \ | |
132 | private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\ | |
133 | pdf_C, \ | |
134 | pdf_N, pdf_E, pdf_S, pdf_W, \ | |
135 | pdf_NE, pdf_SE, pdf_SW, pdf_NW, \ | |
136 | pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \ | |
137 | pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW) | |
138 | #endif | |
0fde6e45 MW |
139 | #ifdef INTEL_OPT_DIRECTIVES |
140 | #pragma ivdep | |
141 | #pragma vector always | |
142 | #pragma simd | |
143 | #endif | |
8cafd9ea | 144 | for (int index = 0; index < nFluid; ++index) { // LOOP list-aa-ria-even |
10988083 MW |
145 | |
146 | #define I(index, dir) P_INDEX_3((nCells), (index), (dir)) | |
147 | ||
148 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[I(index, idx)]; | |
149 | D3Q19_LIST | |
150 | #undef X | |
151 | ||
152 | // #define LID_DRIVEN_CAVITY | |
153 | ||
154 | #ifdef LID_DRIVEN_CAVITY | |
155 | int nX = kd->Dims[0]; | |
156 | int nY = kd->Dims[1]; | |
157 | int nZ = kd->Dims[2]; | |
158 | ||
159 | int x = kdl->Coords[C_INDEX_X(index)]; | |
160 | int y = kdl->Coords[C_INDEX_Y(index)]; | |
161 | int z = kdl->Coords[C_INDEX_Z(index)]; | |
162 | ||
163 | if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) { | |
0fde6e45 MW |
164 | ux = F(0.1) * F(0.5)77; |
165 | uy = F(0.0); | |
166 | uz = F(0.0); | |
10988083 MW |
167 | } else { |
168 | #endif | |
169 | ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - | |
170 | pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; | |
171 | uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN - | |
172 | pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS; | |
173 | uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS - | |
174 | pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS; | |
175 | #ifdef LID_DRIVEN_CAVITY | |
176 | } | |
177 | #endif | |
178 | ||
179 | dens = pdf_C + | |
180 | pdf_N + pdf_E + pdf_S + pdf_W + | |
181 | pdf_NE + pdf_SE + pdf_SW + pdf_NW + | |
182 | pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW + | |
183 | pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW; | |
184 | ||
0fde6e45 | 185 | dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0); |
10988083 MW |
186 | |
187 | // direction: w_0 | |
188 | src[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
189 | ||
190 | // direction: w_1 | |
191 | w_1_indep = w_1*dir_indep_trm; | |
192 | ||
193 | ui = uy; | |
0fde6e45 MW |
194 | evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); |
195 | oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
10988083 MW |
196 | src[I(index, D3Q19_S)] = pdf_N - evenPart - oddPart; |
197 | src[I(index, D3Q19_N)] = pdf_S - evenPart + oddPart; | |
198 | ||
199 | ui = ux; | |
0fde6e45 MW |
200 | evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); |
201 | oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
10988083 MW |
202 | src[I(index, D3Q19_W)] = pdf_E - evenPart - oddPart; |
203 | src[I(index, D3Q19_E)] = pdf_W - evenPart + oddPart; | |
204 | ||
205 | ui = uz; | |
0fde6e45 MW |
206 | evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); |
207 | oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
10988083 MW |
208 | src[I(index, D3Q19_B)] = pdf_T - evenPart - oddPart; |
209 | src[I(index, D3Q19_T)] = pdf_B - evenPart + oddPart; | |
210 | ||
211 | // direction: w_2 | |
212 | w_2_indep = w_2*dir_indep_trm; | |
213 | ||
214 | ui = -ux + uy; | |
0fde6e45 MW |
215 | evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); |
216 | oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
10988083 MW |
217 | src[I(index, D3Q19_SE)] = pdf_NW - evenPart - oddPart; |
218 | src[I(index, D3Q19_NW)] = pdf_SE - evenPart + oddPart; | |
219 | ||
220 | ui = ux + uy; | |
0fde6e45 MW |
221 | evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); |
222 | oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
10988083 MW |
223 | src[I(index, D3Q19_SW)] = pdf_NE - evenPart - oddPart; |
224 | src[I(index, D3Q19_NE)] = pdf_SW - evenPart + oddPart; | |
225 | ||
226 | ui = -ux + uz; | |
0fde6e45 MW |
227 | evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); |
228 | oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
10988083 MW |
229 | src[I(index, D3Q19_BE)] = pdf_TW - evenPart - oddPart; |
230 | src[I(index, D3Q19_TW)] = pdf_BE - evenPart + oddPart; | |
231 | ||
232 | ui = ux + uz; | |
0fde6e45 MW |
233 | evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); |
234 | oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
10988083 MW |
235 | src[I(index, D3Q19_BW)] = pdf_TE - evenPart - oddPart; |
236 | src[I(index, D3Q19_TE)] = pdf_BW - evenPart + oddPart; | |
237 | ||
238 | ui = -uy + uz; | |
0fde6e45 MW |
239 | evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); |
240 | oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
10988083 MW |
241 | src[I(index, D3Q19_BN)] = pdf_TS - evenPart - oddPart; |
242 | src[I(index, D3Q19_TS)] = pdf_BN - evenPart + oddPart; | |
243 | ||
244 | ui = uy + uz; | |
0fde6e45 MW |
245 | evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); |
246 | oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
10988083 MW |
247 | src[I(index, D3Q19_BS)] = pdf_TN - evenPart - oddPart; |
248 | src[I(index, D3Q19_TN)] = pdf_BS - evenPart + oddPart; | |
249 | ||
250 | } // (parallel) loop over fluid nodes | |
251 | ||
252 | X_LIKWID_STOP("list-aa-ria-even"); | |
253 | ||
254 | // save current iteration | |
255 | kdl->Iteration = iter; | |
256 | #ifdef VERIFICATION | |
257 | kd->PdfsActive = src; | |
258 | KernelAddBodyForce(kd, ld, cd); | |
259 | #endif | |
260 | ||
261 | // -------------------------------------------------------------------- | |
262 | // odd time step | |
263 | // -------------------------------------------------------------------- | |
264 | ||
265 | X_LIKWID_START("list-aa-ria-odd"); | |
266 | ||
267 | #ifdef _OPENMP | |
268 | #pragma omp parallel default(none) \ | |
269 | shared(stderr, nFluid, nCells, kd, kdl, kdlr, adjList, omegaOdd, omegaEven, src, consecNodes, nConsecNodes) \ | |
270 | private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\ | |
271 | pdf_C, \ | |
272 | pdf_N, pdf_E, pdf_S, pdf_W, \ | |
273 | pdf_NE, pdf_SE, pdf_SW, pdf_NW, \ | |
274 | pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \ | |
275 | pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW, \ | |
276 | ppdf_C, \ | |
277 | ppdf_N, ppdf_E, ppdf_S, ppdf_W, \ | |
278 | ppdf_NE, ppdf_SE, ppdf_SW, ppdf_NW, \ | |
279 | ppdf_T, ppdf_TN, ppdf_TE, ppdf_TS, ppdf_TW, \ | |
280 | ppdf_B, ppdf_BN, ppdf_BE, ppdf_BS, ppdf_BW, \ | |
281 | consecValue, consecIndex) | |
282 | #endif | |
283 | { | |
284 | int threadId = 0; | |
285 | ||
286 | #ifdef _OPENMP | |
287 | threadId = omp_get_thread_num(); | |
288 | #endif | |
289 | ||
290 | consecIndex = kdlr->ConsecThreadIndices[threadId]; | |
291 | consecValue = 0; | |
292 | ||
293 | int * threadIndices = kdlr->FluidNodeThreadIndices; | |
294 | ||
295 | int nFluidThread = threadIndices[threadId + 1] - threadIndices[threadId]; | |
296 | ||
297 | int indexStart = threadIndices[threadId]; | |
298 | int indexStop = threadIndices[threadId] + nFluidThread; | |
299 | ||
0fde6e45 | 300 | // Because of runlength coding iterations are not independent. |
8cafd9ea | 301 | for (int index = indexStart; index < indexStop; ++index) { // LOOP list-aa-ria-odd |
10988083 MW |
302 | |
303 | #define I(index, dir) P_INDEX_3((nCells), (index), (dir)) | |
304 | ||
305 | #if 1 | |
306 | if (consecValue > 0) { | |
307 | --consecValue; | |
308 | // Increment all pdf pointers. | |
309 | #define X(name, idx, idxinv, _x, _y, _z) ++JOIN(ppdf_,name); | |
310 | D3Q19_LIST | |
311 | #undef X | |
312 | } | |
313 | else { | |
314 | Assert(consecIndex < nConsecNodes); | |
315 | ||
316 | consecValue = consecNodes[consecIndex] - 1; | |
317 | // Load new pointers to PDFs of local cell: | |
318 | ||
319 | adjListIndex = index * N_D3Q19_IDX; | |
320 | ||
321 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[adjList[adjListIndex + idxinv]]); | |
322 | D3Q19_LIST_WO_C | |
323 | #undef X | |
324 | ||
325 | ppdf_C = &(src[P_INDEX_3(nCells, index, D3Q19_C)]); | |
326 | ++consecIndex; | |
327 | } | |
328 | ||
329 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = *JOIN(ppdf_,name); | |
330 | D3Q19_LIST | |
331 | #undef X | |
332 | #else | |
333 | adjListIndex = index * N_D3Q19_IDX; | |
334 | ||
335 | // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ... | |
336 | pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)]; | |
337 | ||
338 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(ppdf_,name) = &(src[adjList[adjListIndex + idxinv]]); | |
339 | D3Q19_LIST_WO_C | |
340 | #undef X | |
341 | ||
342 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = src[adjList[adjListIndex + idxinv]]; | |
343 | D3Q19_LIST_WO_C | |
344 | #undef X | |
345 | #endif | |
346 | ||
347 | #ifdef LID_DRIVEN_CAVITY | |
348 | int nX = kd->Dims[0]; | |
349 | int nY = kd->Dims[1]; | |
350 | int nZ = kd->Dims[2]; | |
351 | ||
352 | int x = kdl->Coords[C_INDEX_X(index)]; | |
353 | int y = kdl->Coords[C_INDEX_Y(index)]; | |
354 | int z = kdl->Coords[C_INDEX_Z(index)]; | |
355 | ||
356 | if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) { | |
0fde6e45 MW |
357 | ux = F(0.1) * F(0.5)77; |
358 | uy = F(0.0); | |
359 | uz = F(0.0); | |
10988083 MW |
360 | } else { |
361 | #endif | |
362 | ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - | |
363 | pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; | |
364 | uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN - | |
365 | pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS; | |
366 | uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS - | |
367 | pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS; | |
368 | #ifdef LID_DRIVEN_CAVITY | |
369 | } | |
370 | #endif | |
371 | ||
372 | dens = pdf_C + | |
373 | pdf_N + pdf_E + pdf_S + pdf_W + | |
374 | pdf_NE + pdf_SE + pdf_SW + pdf_NW + | |
375 | pdf_T + pdf_TN + pdf_TE + pdf_TS + pdf_TW + | |
376 | pdf_B + pdf_BN + pdf_BE + pdf_BS + pdf_BW; | |
377 | ||
0fde6e45 | 378 | dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0); |
10988083 MW |
379 | |
380 | adjListIndex = index * N_D3Q19_IDX; | |
381 | ||
382 | // direction: w_0 | |
383 | src[I(index, D3Q19_C) ] = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm); | |
384 | ||
385 | // direction: w_1 | |
386 | w_1_indep = w_1*dir_indep_trm; | |
387 | ||
388 | ui = uy; | |
0fde6e45 MW |
389 | evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep ); |
390 | oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 ); | |
10988083 MW |
391 | *ppdf_S = pdf_N - evenPart - oddPart; |
392 | *ppdf_N = pdf_S - evenPart + oddPart; | |
393 | ||
394 | ui = ux; | |
0fde6e45 MW |
395 | evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep ); |
396 | oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 ); | |
10988083 MW |
397 | *ppdf_W = pdf_E - evenPart - oddPart; |
398 | *ppdf_E = pdf_W - evenPart + oddPart; | |
399 | ||
400 | ui = uz; | |
0fde6e45 MW |
401 | evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep ); |
402 | oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 ); | |
10988083 MW |
403 | *ppdf_B = pdf_T - evenPart - oddPart; |
404 | *ppdf_T = pdf_B - evenPart + oddPart; | |
405 | ||
406 | // direction: w_2 | |
407 | w_2_indep = w_2*dir_indep_trm; | |
408 | ||
409 | ui = -ux + uy; | |
0fde6e45 MW |
410 | evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep ); |
411 | oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 ); | |
10988083 MW |
412 | *ppdf_SE = pdf_NW - evenPart - oddPart; |
413 | *ppdf_NW = pdf_SE - evenPart + oddPart; | |
414 | ||
415 | ui = ux + uy; | |
0fde6e45 MW |
416 | evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep ); |
417 | oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 ); | |
10988083 MW |
418 | *ppdf_SW = pdf_NE - evenPart - oddPart; |
419 | *ppdf_NE = pdf_SW - evenPart + oddPart; | |
420 | ||
421 | ui = -ux + uz; | |
0fde6e45 MW |
422 | evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep ); |
423 | oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 ); | |
10988083 MW |
424 | *ppdf_BE = pdf_TW - evenPart - oddPart; |
425 | *ppdf_TW = pdf_BE - evenPart + oddPart; | |
426 | ||
427 | ui = ux + uz; | |
0fde6e45 MW |
428 | evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep ); |
429 | oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 ); | |
10988083 MW |
430 | *ppdf_BW = pdf_TE - evenPart - oddPart; |
431 | *ppdf_TE = pdf_BW - evenPart + oddPart; | |
432 | ||
433 | ui = -uy + uz; | |
0fde6e45 MW |
434 | evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep ); |
435 | oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 ); | |
10988083 MW |
436 | *ppdf_BN = pdf_TS - evenPart - oddPart; |
437 | *ppdf_TS = pdf_BN - evenPart + oddPart; | |
438 | ||
439 | ui = uy + uz; | |
0fde6e45 MW |
440 | evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep ); |
441 | oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 ); | |
10988083 MW |
442 | *ppdf_BS = pdf_TN - evenPart - oddPart; |
443 | *ppdf_TN = pdf_BS - evenPart + oddPart; | |
444 | ||
445 | #undef I | |
446 | } // loop over fluid nodes | |
447 | } // end pragma omp parallel | |
448 | ||
449 | X_LIKWID_STOP("list-aa-ria-odd"); | |
450 | ||
451 | // save current iteration | |
452 | kdl->Iteration = iter + 1; | |
453 | ||
454 | #ifdef VERIFICATION | |
455 | kd->PdfsActive = src; | |
456 | KernelAddBodyForce(kd, ld, cd); | |
457 | #endif | |
458 | ||
459 | #ifdef VTK_OUTPUT | |
460 | if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) { | |
461 | kd->PdfsActive = src; | |
462 | VtkWrite(ld, kd, cd, iter); | |
463 | } | |
464 | #endif | |
465 | ||
466 | #ifdef STATISTICS | |
467 | kd->PdfsActive = src; | |
468 | KernelStatistics(kd, ld, cd, iter); | |
469 | #endif | |
470 | ||
471 | ||
472 | } // for (int iter = 0; ... | |
473 | ||
8cafd9ea MW |
474 | X_KERNEL_END(kernelData); |
475 | ||
10988083 MW |
476 | #ifdef VTK_OUTPUT |
477 | if (cd->VtkOutput) { | |
478 | kd->PdfsActive = src; | |
479 | VtkWrite(ld, kd, cd, maxIterations); | |
480 | } | |
481 | #endif | |
482 | ||
483 | #ifdef STATISTICS | |
484 | kd->PdfsActive = src; | |
485 | KernelStatistics(kd, ld, cd, maxIterations); | |
486 | #endif | |
487 | ||
488 | return; | |
489 | } |