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