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 "Kernel.h" | |
28 | #include "Lattice.h" | |
29 | ||
30 | #include <stdlib.h> | |
31 | #include <string.h> | |
32 | #include <math.h> | |
33 | ||
34 | #define X(name, idx, idx_inv, x, y, z) , x | |
35 | int D3Q19_X[] = { | |
36 | EXPAND(D3Q19_LIST) | |
37 | }; | |
38 | #undef X | |
39 | ||
40 | #define X(name, idx, idx_inv, x, y, z) , y | |
41 | int D3Q19_Y[] = { | |
42 | EXPAND(D3Q19_LIST) | |
43 | }; | |
44 | #undef X | |
45 | ||
46 | #define X(name, idx, idx_inv, x, y, z) , z | |
47 | int D3Q19_Z[] = { | |
48 | EXPAND(D3Q19_LIST) | |
49 | }; | |
50 | #undef X | |
51 | ||
52 | #define X(name, idx, idxinv, x, y, z) , idxinv | |
53 | int D3Q19_INV[] = { | |
54 | EXPAND(D3Q19_LIST) | |
55 | }; | |
56 | #undef X | |
57 | ||
58 | ||
59 | #define X(name, idx, idxinv, x, y, z) , STRINGIFY(name) | |
60 | const char * D3Q19_NAMES[N_D3Q19] = { | |
61 | EXPAND(D3Q19_LIST) | |
62 | }; | |
63 | #undef X | |
64 | ||
65 | void KernelComputeBoundaryConditions(KernelData * kd, LatticeDesc * ld, CaseData * cd) | |
66 | { | |
67 | Assert(kd != NULL); | |
68 | Assert(ld != NULL); | |
69 | Assert(cd != NULL); | |
70 | ||
71 | Assert(cd->RhoIn > 0.0); | |
72 | Assert(cd->RhoOut > 0.0); | |
73 | ||
74 | PdfT rho_in = cd->RhoIn; | |
75 | PdfT rho_out = cd->RhoOut; | |
76 | PdfT rho_in_inv = 1.0 / rho_in; | |
77 | PdfT rho_out_inv = 1.0 / rho_out; | |
78 | PdfT indep_ux = 0.0; | |
79 | ||
80 | PdfT dens; | |
81 | PdfT ux; | |
82 | ||
83 | const PdfT one_third = 1.0 / 3.0; | |
84 | const PdfT one_fourth = 1.0 / 4.0; | |
85 | const PdfT one_sixth = 1.0 / 6.0; | |
86 | ||
87 | PdfT pdfs[N_D3Q19]; | |
88 | ||
89 | int nX = kd->Dims[0]; | |
90 | int nY = kd->Dims[1]; | |
91 | int nZ = kd->Dims[2]; | |
92 | ||
93 | int x; | |
94 | int x_in = 0; | |
95 | int x_out = nX - 1; | |
96 | ||
97 | double density_in = 0.0; | |
98 | double density_out = 0.0; | |
99 | ||
100 | // update inlet / outlet boundary conditions | |
101 | for (int z = 1; z < nZ - 1; ++z) { | |
102 | for (int y = 1; y < nY - 1; ++y) { | |
103 | ||
104 | ||
105 | // ----------------------------------------------------------------------------- | |
106 | // update inlet conditions | |
107 | ||
108 | if (ld->Lattice[L_INDEX_4(ld->Dims, x_in, y, z)] == LAT_CELL_INLET) { | |
109 | ||
110 | x = x_in; | |
111 | ||
112 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_C , pdfs + D3Q19_C); | |
113 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_T , pdfs + D3Q19_T); | |
114 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_B , pdfs + D3Q19_B); | |
115 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_S , pdfs + D3Q19_S); | |
116 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_N , pdfs + D3Q19_N); | |
117 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_TS, pdfs + D3Q19_TS); | |
118 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_BS, pdfs + D3Q19_BS); | |
119 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_TN, pdfs + D3Q19_TN); | |
120 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_BN, pdfs + D3Q19_BN); | |
121 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_SW, pdfs + D3Q19_SW); | |
122 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_TW, pdfs + D3Q19_TW); | |
123 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_W , pdfs + D3Q19_W); | |
124 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_BW, pdfs + D3Q19_BW); | |
125 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_NW, pdfs + D3Q19_NW); | |
126 | ||
127 | dens = rho_in; | |
128 | ||
129 | ux = 1 - (pdfs[D3Q19_C] + | |
130 | (pdfs[D3Q19_T] + pdfs[D3Q19_B] + pdfs[D3Q19_S] + pdfs[D3Q19_N]) + | |
131 | (pdfs[D3Q19_TS] + pdfs[D3Q19_BS] + pdfs[D3Q19_TN] + pdfs[D3Q19_BN]) + | |
132 | 2 * (pdfs[D3Q19_SW] + pdfs[D3Q19_TW] + pdfs[D3Q19_W] + pdfs[D3Q19_BW] + pdfs[D3Q19_NW])) * rho_in_inv; | |
133 | ||
134 | indep_ux = one_sixth * dens * ux; | |
135 | ||
136 | pdfs[D3Q19_E ] = pdfs[D3Q19_W] + one_third * dens * ux; | |
137 | pdfs[D3Q19_NE] = pdfs[D3Q19_SW] - one_fourth * (pdfs[D3Q19_N] - pdfs[D3Q19_S]) + indep_ux; | |
138 | pdfs[D3Q19_SE] = pdfs[D3Q19_NW] + one_fourth * (pdfs[D3Q19_N] - pdfs[D3Q19_S]) + indep_ux; | |
139 | pdfs[D3Q19_TE] = pdfs[D3Q19_BW] - one_fourth * (pdfs[D3Q19_T] - pdfs[D3Q19_B]) + indep_ux; | |
140 | pdfs[D3Q19_BE] = pdfs[D3Q19_TW] + one_fourth * (pdfs[D3Q19_T] - pdfs[D3Q19_B]) + indep_ux; | |
141 | ||
142 | ||
143 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_E , pdfs[D3Q19_E ]); | |
144 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_NE, pdfs[D3Q19_NE]); | |
145 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_SE, pdfs[D3Q19_SE]); | |
146 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_TE, pdfs[D3Q19_TE]); | |
147 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_BE, pdfs[D3Q19_BE]); | |
148 | ||
149 | for(int d = 0; d < N_D3Q19; ++d) { | |
150 | density_in += pdfs[d]; | |
151 | } | |
152 | } | |
153 | ||
154 | // ----------------------------------------------------------------------------- | |
155 | // update outlet conditions | |
156 | ||
157 | if (ld->Lattice[L_INDEX_4(ld->Dims, x_out, y, z)] == LAT_CELL_OUTLET) { | |
158 | // update outlet conditions | |
159 | ||
160 | x = x_out; | |
161 | ||
162 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_C , pdfs + D3Q19_C ); | |
163 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_T , pdfs + D3Q19_T ); | |
164 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_B , pdfs + D3Q19_B ); | |
165 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_S , pdfs + D3Q19_S ); | |
166 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_N , pdfs + D3Q19_N ); | |
167 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_TS, pdfs + D3Q19_TS); | |
168 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_BS, pdfs + D3Q19_BS); | |
169 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_TN, pdfs + D3Q19_TN); | |
170 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_BN, pdfs + D3Q19_BN); | |
171 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_NE, pdfs + D3Q19_NE); | |
172 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_BE, pdfs + D3Q19_BE); | |
173 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_E , pdfs + D3Q19_E ); | |
174 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_TE, pdfs + D3Q19_TE); | |
175 | kd->BoundaryConditionsGetPdf(kd, x, y, z, D3Q19_SE, pdfs + D3Q19_SE); | |
176 | ||
177 | dens = rho_out; | |
178 | ||
179 | ux = -1 + (pdfs[D3Q19_C] + | |
180 | (pdfs[D3Q19_T] + pdfs[D3Q19_B] + pdfs[D3Q19_S] + pdfs[D3Q19_N]) + | |
181 | (pdfs[D3Q19_TS] + pdfs[D3Q19_BS] + pdfs[D3Q19_TN] + pdfs[D3Q19_BN]) + | |
182 | 2 * (pdfs[D3Q19_NE] + pdfs[D3Q19_BE] + pdfs[D3Q19_E] + pdfs[D3Q19_TE] + pdfs[D3Q19_SE])) * rho_out_inv; | |
183 | indep_ux = one_sixth * dens * ux; | |
184 | ||
185 | pdfs[D3Q19_W ] = pdfs[D3Q19_E] - one_third * dens * ux; | |
186 | pdfs[D3Q19_SW] = pdfs[D3Q19_NE] + one_fourth * (pdfs[D3Q19_N] - pdfs[D3Q19_S]) - indep_ux; | |
187 | pdfs[D3Q19_NW] = pdfs[D3Q19_SE] - one_fourth * (pdfs[D3Q19_N] - pdfs[D3Q19_S]) - indep_ux; | |
188 | pdfs[D3Q19_BW] = pdfs[D3Q19_TE] + one_fourth * (pdfs[D3Q19_T] - pdfs[D3Q19_B]) - indep_ux; | |
189 | pdfs[D3Q19_TW] = pdfs[D3Q19_BE] - one_fourth * (pdfs[D3Q19_T] - pdfs[D3Q19_B]) - indep_ux; | |
190 | ||
191 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_W , pdfs[D3Q19_W ]); | |
192 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_NW, pdfs[D3Q19_NW]); | |
193 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_SW, pdfs[D3Q19_SW]); | |
194 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_TW, pdfs[D3Q19_TW]); | |
195 | kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_BW, pdfs[D3Q19_BW]); | |
196 | ||
197 | for(int d = 0; d < N_D3Q19; ++d) { | |
198 | density_out += pdfs[d]; | |
199 | } | |
200 | } | |
201 | } | |
202 | } | |
203 | ||
204 | // DEBUG: printf("# density inlet: %e density outlet: %e\n", density_in, density_out); | |
205 | ||
206 | } | |
207 | ||
208 | ||
209 | PdfT KernelDensity(KernelData * kd, LatticeDesc * ld) | |
210 | { | |
211 | Assert(kd != NULL); | |
212 | Assert(ld != NULL); | |
213 | ||
214 | Assert(ld->Lattice != NULL); | |
215 | Assert(ld->Dims != NULL); | |
216 | ||
217 | Assert(ld->Dims[0] > 0); | |
218 | Assert(ld->Dims[1] > 0); | |
219 | Assert(ld->Dims[2] > 0); | |
220 | ||
221 | int * lDims = ld->Dims; | |
222 | int nX = lDims[0]; | |
223 | int nY = lDims[1]; | |
224 | int nZ = lDims[2]; | |
225 | ||
226 | PdfT pdfs[N_D3Q19] = { -1.0 }; | |
227 | PdfT density = 0.0; | |
228 | ||
229 | for(int z = 0; z < nZ; ++z) { | |
230 | for(int y = 0; y < nY; ++y) { | |
231 | for(int x = 0; x < nX; ++x) { | |
232 | ||
233 | if(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) { | |
234 | ||
235 | kd->GetNode(kd, x, y, z, pdfs); | |
236 | ||
237 | for(int d = 0; d < N_D3Q19; ++d) { | |
238 | // if (pdfs[d] < 0.0) { | |
239 | // printf("# %d %d %d %d < 0 %e %s\n", x, y, z, d, pdfs[d], D3Q19_NAMES[d]); | |
240 | // exit(1); | |
241 | // } | |
242 | density += pdfs[d]; | |
243 | } | |
244 | } | |
245 | ||
246 | } | |
247 | } | |
248 | } | |
249 | ||
250 | return density / ld->nFluid; | |
251 | } | |
252 | ||
253 | ||
254 | // prescribes a given density | |
255 | void KernelSetInitialDensity(LatticeDesc * ld, KernelData * kd, CaseData * cd) | |
256 | { | |
257 | int * lDims = ld->Dims; | |
258 | ||
259 | PdfT rho_in = cd->RhoIn; | |
260 | PdfT rho_out = cd->RhoOut; | |
261 | ||
262 | PdfT ux = 0.0; | |
263 | PdfT uy = 0.0; | |
264 | PdfT uz = 0.0; | |
265 | PdfT dens = 1.0; | |
266 | ||
267 | PdfT omega = cd->Omega; | |
268 | ||
269 | PdfT w_0 = 1.0 / 3.0; | |
270 | PdfT w_1 = 1.0 / 18.0; | |
271 | PdfT w_2 = 1.0 / 36.0; | |
272 | ||
273 | PdfT dir_indep_trm; | |
274 | PdfT omega_w0 = 3.0 * w_0 * omega; | |
275 | PdfT omega_w1 = 3.0 * w_1 * omega; | |
276 | PdfT omega_w2 = 3.0 * w_2 * omega; | |
277 | PdfT one_third = 1.0 / 3.0; | |
278 | ||
279 | int nX = lDims[0]; | |
280 | int nY = lDims[1]; | |
281 | int nZ = lDims[2]; | |
282 | ||
283 | PdfT pdfs[N_D3Q19]; | |
284 | ||
285 | #ifdef _OPENMP | |
286 | #pragma omp parallel for collapse(3) | |
287 | #endif | |
288 | for(int z = 0; z < nZ; ++z) { for(int y = 0; y < nY; ++y) { for(int x = 0; x < nX; ++x) { | |
289 | ||
290 | if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) { | |
291 | // TODO: fix later. | |
292 | // if((caseData->geoType == GEO_TYPE_CHANNEL) || (caseData->geoType == GEO_TYPE_RCHANNEL)) | |
293 | dens = rho_in + (rho_out - rho_in)*(x)/(nX-1.0); | |
294 | ||
295 | #define SQR(a) ((a)*(a)) | |
296 | ||
297 | dir_indep_trm = one_third * dens - 0.5 * (ux * ux + uy * uy + uz * uz); | |
298 | ||
299 | pdfs[D3Q19_C] = omega_w0 * (dir_indep_trm); | |
300 | ||
301 | pdfs[D3Q19_NW] = omega_w2 * (dir_indep_trm - (ux - uy) + 1.5 * SQR(ux - uy)); | |
302 | pdfs[D3Q19_SE] = omega_w2 * (dir_indep_trm + (ux - uy) + 1.5 * SQR(ux - uy)); | |
303 | ||
304 | pdfs[D3Q19_NE] = omega_w2 * (dir_indep_trm + (ux + uy) + 1.5 * SQR(ux + uy)); | |
305 | pdfs[D3Q19_SW] = omega_w2 * (dir_indep_trm - (ux + uy) + 1.5 * SQR(ux + uy)); | |
306 | ||
307 | ||
308 | pdfs[D3Q19_TW] = omega_w2 * (dir_indep_trm - (ux - uz) + 1.5 * SQR(ux - uz)); | |
309 | pdfs[D3Q19_BE] = omega_w2 * (dir_indep_trm + (ux - uz) + 1.5 * SQR(ux - uz)); | |
310 | ||
311 | pdfs[D3Q19_TE] = omega_w2 * (dir_indep_trm + (ux + uz) + 1.5 * SQR(ux + uz)); | |
312 | pdfs[D3Q19_BW] = omega_w2 * (dir_indep_trm - (ux + uz) + 1.5 * SQR(ux + uz)); | |
313 | ||
314 | ||
315 | pdfs[D3Q19_TS] = omega_w2 * (dir_indep_trm - (uy - uz) + 1.5 * SQR(uy - uz)); | |
316 | pdfs[D3Q19_BN] = omega_w2 * (dir_indep_trm + (uy - uz) + 1.5 * SQR(uy - uz)); | |
317 | ||
318 | pdfs[D3Q19_TN] = omega_w2 * (dir_indep_trm + (uy + uz) + 1.5 * SQR(uy + uz)); | |
319 | pdfs[D3Q19_BS] = omega_w2 * (dir_indep_trm - (uy + uz) + 1.5 * SQR(uy + uz)); | |
320 | ||
321 | ||
322 | pdfs[D3Q19_N] = omega_w1 * (dir_indep_trm + uy + 1.5 * SQR(uy)); | |
323 | pdfs[D3Q19_S] = omega_w1 * (dir_indep_trm - uy + 1.5 * SQR(uy)); | |
324 | ||
325 | pdfs[D3Q19_E] = omega_w1 * (dir_indep_trm + ux + 1.5 * SQR(ux)); | |
326 | pdfs[D3Q19_W] = omega_w1 * (dir_indep_trm - ux + 1.5 * SQR(ux)); | |
327 | ||
328 | pdfs[D3Q19_T] = omega_w1 * (dir_indep_trm + uz + 1.5 * SQR(uz)); | |
329 | pdfs[D3Q19_B] = omega_w1 * (dir_indep_trm - uz + 1.5 * SQR(uz)); | |
330 | ||
331 | ||
332 | kd->SetNode(kd, x, y, z, pdfs); | |
333 | ||
334 | #undef SQR | |
335 | } | |
336 | } } } | |
337 | } | |
338 | ||
339 | ||
340 | // prescribes a given velocity | |
341 | void KernelSetInitialVelocity(LatticeDesc * ld, KernelData * kd, CaseData * cd) | |
342 | { | |
343 | ||
344 | int * lDims = ld->Dims; | |
345 | ||
346 | // TODO: ux is overriden below... | |
347 | PdfT ux = 0.09; // caseData->initUx; | |
348 | PdfT uy = 0.0; // caseData->initUy; | |
349 | PdfT uz = 0.0; // caseData->initUz; | |
350 | PdfT dens = 1.0; | |
351 | ||
352 | PdfT omega = cd->Omega; | |
353 | ||
354 | PdfT w_0 = 1.0 / 3.0; | |
355 | PdfT w_1 = 1.0 / 18.0; | |
356 | PdfT w_2 = 1.0 / 36.0; | |
357 | ||
358 | PdfT dir_indep_trm; | |
359 | PdfT omega_w0 = 3.0 * w_0 * omega; | |
360 | PdfT omega_w1 = 3.0 * w_1 * omega; | |
361 | PdfT omega_w2 = 3.0 * w_2 * omega; | |
362 | PdfT one_third = 1.0 / 3.0; | |
363 | ||
364 | int nX = lDims[0]; | |
365 | int nY = lDims[1]; | |
366 | int nZ = lDims[2]; | |
367 | ||
368 | PdfT pdfs[N_D3Q19]; | |
369 | ||
370 | PdfT density; | |
371 | ||
372 | #ifdef _OPENMP | |
373 | #pragma omp parallel for collapse(3) | |
374 | #endif | |
375 | for(int z = 0; z < nZ; ++z) { for(int y = 0; y < nY; ++y) { for(int x = 0; x < nX; ++x) { | |
376 | ||
377 | if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] == LAT_CELL_FLUID) { | |
378 | ||
379 | ux = 0.0; | |
380 | uy = 0.0; | |
381 | uz = 0.0; | |
382 | ||
383 | kd->GetNode(kd, x, y, z, pdfs); | |
384 | ||
385 | ||
386 | density = 0.0; | |
387 | ||
388 | #define X(name, idx, idxinv, _x, _y, _z) density += pdfs[idx]; | |
389 | D3Q19_LIST | |
390 | #undef X | |
391 | ||
392 | ||
393 | #define SQR(a) ((a)*(a)) | |
394 | dir_indep_trm = one_third * dens - 0.5 * (ux * ux + uy * uy + uz * uz); | |
395 | ||
396 | pdfs[D3Q19_C] = omega_w0 * (dir_indep_trm); | |
397 | ||
398 | pdfs[D3Q19_NW] = omega_w2 * (dir_indep_trm - (ux - uy) + 1.5 * SQR(ux - uy)); | |
399 | pdfs[D3Q19_SE] = omega_w2 * (dir_indep_trm + (ux - uy) + 1.5 * SQR(ux - uy)); | |
400 | ||
401 | pdfs[D3Q19_NE] = omega_w2 * (dir_indep_trm + (ux + uy) + 1.5 * SQR(ux + uy)); | |
402 | pdfs[D3Q19_SW] = omega_w2 * (dir_indep_trm - (ux + uy) + 1.5 * SQR(ux + uy)); | |
403 | ||
404 | ||
405 | pdfs[D3Q19_TW] = omega_w2 * (dir_indep_trm - (ux - uz) + 1.5 * SQR(ux - uz)); | |
406 | pdfs[D3Q19_BE] = omega_w2 * (dir_indep_trm + (ux - uz) + 1.5 * SQR(ux - uz)); | |
407 | ||
408 | pdfs[D3Q19_TE] = omega_w2 * (dir_indep_trm + (ux + uz) + 1.5 * SQR(ux + uz)); | |
409 | pdfs[D3Q19_BW] = omega_w2 * (dir_indep_trm - (ux + uz) + 1.5 * SQR(ux + uz)); | |
410 | ||
411 | ||
412 | pdfs[D3Q19_TS] = omega_w2 * (dir_indep_trm - (uy - uz) + 1.5 * SQR(uy - uz)); | |
413 | pdfs[D3Q19_BN] = omega_w2 * (dir_indep_trm + (uy - uz) + 1.5 * SQR(uy - uz)); | |
414 | ||
415 | pdfs[D3Q19_TN] = omega_w2 * (dir_indep_trm + (uy + uz) + 1.5 * SQR(uy + uz)); | |
416 | pdfs[D3Q19_BS] = omega_w2 * (dir_indep_trm - (uy + uz) + 1.5 * SQR(uy + uz)); | |
417 | ||
418 | ||
419 | pdfs[D3Q19_N] = omega_w1 * (dir_indep_trm + uy + 1.5 * SQR(uy)); | |
420 | pdfs[D3Q19_S] = omega_w1 * (dir_indep_trm - uy + 1.5 * SQR(uy)); | |
421 | ||
422 | pdfs[D3Q19_E] = omega_w1 * (dir_indep_trm + ux + 1.5 * SQR(ux)); | |
423 | pdfs[D3Q19_W] = omega_w1 * (dir_indep_trm - ux + 1.5 * SQR(ux)); | |
424 | ||
425 | pdfs[D3Q19_T] = omega_w1 * (dir_indep_trm + uz + 1.5 * SQR(uz)); | |
426 | pdfs[D3Q19_B] = omega_w1 * (dir_indep_trm - uz + 1.5 * SQR(uz)); | |
427 | ||
428 | #undef SQR | |
429 | ||
430 | ||
431 | kd->SetNode(kd, x, y, z, pdfs); | |
432 | } | |
433 | } } } | |
434 | ||
435 | } | |
436 | ||
437 | // Compute analytical x velocity for channel flow. | |
438 | // | |
439 | // Formula 7 from Kutay et al. "Laboratory validation of lattice Boltzmann method for modeling | |
440 | // pore-scale flow in granular materials", doi:10.1016/j.compgeo.2006.08.002. | |
441 | // | |
442 | // also formula 10 from | |
443 | // Pan et al. "An evaluation of lattice Boltzmann equation methods for simulating flow | |
444 | // through porous media", doi:10.1016/S0167-5648(04)80040-6. | |
445 | // | |
446 | // calculate velocity in a pipe for a given radius | |
447 | // | |
448 | static PdfT CalcXVelForPipeProfile(PdfT maxRadiusSquared, PdfT curRadiusSquared, PdfT xForce, PdfT viscosity) | |
449 | { | |
450 | return xForce*(maxRadiusSquared - curRadiusSquared) / (2.0*viscosity); | |
451 | } | |
452 | ||
453 | static void KernelGetXSlice(LatticeDesc * ld, KernelData * kd, CaseData * cd, PdfT * outputArray, int xPos) | |
454 | { | |
455 | Assert(ld != NULL); | |
456 | Assert(kd != NULL); | |
457 | ||
458 | int nY = ld->Dims[1]; | |
459 | int nZ = ld->Dims[2]; | |
460 | ||
461 | Assert(xPos >= 0); | |
462 | Assert(xPos < ld->Dims[0]); | |
463 | ||
464 | ||
465 | PdfT ux = 0.0; | |
466 | ||
467 | // Declare pdf_N, pdf_E, pdf_S, pdf_W, ... | |
468 | #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name); | |
469 | D3Q19_LIST | |
470 | #undef X | |
471 | PdfT pdfs[N_D3Q19]; | |
472 | ||
473 | for(int z = 0; z < nZ; ++z) { | |
474 | for(int y = 0; y < nY; ++y) { | |
475 | ||
476 | if (ld->Lattice[L_INDEX_4(ld->Dims, xPos, y, z)] != LAT_CELL_OBSTACLE) { | |
477 | kd->GetNode(kd, xPos, y, z, pdfs); | |
478 | ||
479 | #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = pdfs[idx]; | |
480 | D3Q19_LIST | |
481 | #undef X | |
482 | UNUSED(pdf_C); UNUSED(pdf_S); UNUSED(pdf_N); UNUSED(pdf_T); UNUSED(pdf_B); | |
483 | UNUSED(pdf_TN); UNUSED(pdf_BN); UNUSED(pdf_TS); UNUSED(pdf_BS); | |
484 | ||
485 | ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - | |
486 | pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; | |
487 | ||
488 | #ifdef VERIFICATION | |
489 | ux += 0.5 * cd->XForce; | |
490 | #endif | |
491 | ||
492 | outputArray[y * nZ + z] = ux; | |
493 | } | |
494 | else { | |
495 | outputArray[y * nZ + z] = 0.0; | |
496 | } | |
497 | } | |
498 | } | |
499 | ||
500 | } | |
501 | ||
502 | // Verification of channel profile with analytical solution. | |
503 | // Taken from Kutay et al. "Laboratory validation of lattice Boltzmann method for modeling | |
504 | // pore-scale flow in granular materials", doi:10.1016/j.compgeo.2006.08.002. and | |
505 | // Pan et al. "An evaluation of lattice Boltzmann equation methods for simulating flow | |
506 | // through porous media", doi:10.1016/S0167-5648(04)80040-6 | |
507 | // | |
508 | void KernelVerifiy(LatticeDesc * ld, KernelData * kd, CaseData * cd, PdfT * errorNorm) | |
509 | { | |
510 | Assert(ld != NULL); | |
511 | Assert(kd != NULL); | |
512 | Assert(cd != NULL); | |
513 | Assert(errorNorm != NULL); | |
514 | ||
515 | int nX = ld->Dims[0]; | |
516 | int nY = ld->Dims[1]; | |
517 | int nZ = ld->Dims[2]; | |
518 | ||
519 | PdfT omega = cd->Omega; | |
520 | PdfT viscosity = (1.0 / omega - 0.5) / 3.0; | |
521 | ||
522 | // ux averaged across cross sections in x direction | |
523 | PdfT * outputArray = (PdfT *)malloc(nZ * nY * sizeof(PdfT)); | |
524 | Verify(outputArray != NULL); | |
525 | ||
526 | memset(outputArray, -10, nZ*nY*sizeof(PdfT)); | |
527 | ||
528 | // uncomment this to get values averaged along the x-axis | |
529 | //AveragePipeCrossSections(ld, kd, outputArray); | |
530 | KernelGetXSlice(ld, kd, cd, outputArray, (int)(nX/2)); | |
531 | ||
532 | ||
533 | FILE * fh; | |
534 | char fileName[1024]; | |
535 | PdfT tmpAvgUx = 0.0; | |
536 | PdfT tmpAnalyUx = 0.0; | |
537 | int flagEvenNy = 0; | |
538 | int y = 0; | |
539 | ||
540 | if (nY % 2 == 0) | |
541 | flagEvenNy = 1; | |
542 | ||
543 | y = (nY-flagEvenNy-1)/2; | |
544 | ||
545 | snprintf(fileName, sizeof(fileName), "flow-profile.dat"); | |
546 | ||
547 | printf("# Kernel validation: writing profile to %s\n", fileName); | |
548 | ||
549 | fh = fopen(fileName, "w"); | |
550 | ||
551 | if(fh == NULL) { | |
552 | printf("ERROR: opening file %s failed.\n", fileName); | |
553 | exit(1); | |
554 | } | |
555 | ||
556 | fprintf(fh, "# Flow profile in Z direction. Taken at the middle of the X length (= %d) of total length %d.\n", nZ / 2, nZ); | |
557 | // fprintf(fh, "# Snapshot taken at iteration %d.\n", iteration); | |
558 | fprintf(fh, "# Plot on terminal: gnuplot -e \"set terminal dumb; plot \\\"%s\\\" u 1:3 t \\\"analytical\\\", \\\"\\\" u 1:4 t \\\"simulation\\\";\"\n", fileName); | |
559 | fprintf(fh, "# Plot graphically: gnuplot -e \"plot \\\"%s\\\" u 1:3 w linesp t \\\"analytical\\\", \\\"\\\" u 1:4 w linesp t \\\"simulation\\\"; pause -1;\"\n", fileName); | |
560 | fprintf(fh, "# z coord., radius, analytic, simulation, diff abs, diff rel, undim_analytic, undim_sim\n"); | |
561 | ||
562 | double deviation = 0.0; | |
563 | double curRadiusSquared; | |
564 | double center = nY / 2.0; | |
565 | double minDiameter = nY; | |
566 | #define SQR(a) ((a)*(a)) | |
567 | double minRadiusSquared = SQR(minDiameter / 2.0 - 1.0); | |
568 | #undef SQR | |
569 | double u_max = cd->XForce*minRadiusSquared/(2.0*viscosity); | |
570 | ||
571 | for(int z = 0; z < nZ; ++z) { | |
572 | ||
573 | fprintf(fh, "%d\t", z); | |
574 | ||
575 | #define SQR(a) ((a)*(a)) | |
576 | curRadiusSquared = SQR(z-center+0.5); | |
577 | ||
578 | ||
579 | // dimensionless radius | |
580 | fprintf(fh, "%e\t", (z-center+0.5)/center); | |
581 | ||
582 | // analytic profile | |
583 | if(curRadiusSquared >= minRadiusSquared) | |
584 | tmpAnalyUx = 0.0; | |
585 | else | |
586 | tmpAnalyUx = CalcXVelForPipeProfile(minRadiusSquared, curRadiusSquared, cd->XForce, viscosity); | |
587 | ||
588 | //averaged profile | |
589 | if(flagEvenNy == 1) | |
590 | tmpAvgUx = (outputArray[y*nZ + z] + outputArray[(y+1)*nZ + z])/2.0; | |
591 | else | |
592 | tmpAvgUx = outputArray[y*nZ + z]; | |
593 | ||
594 | fprintf(fh, "%e\t", tmpAnalyUx); | |
595 | fprintf(fh, "%e\t", tmpAvgUx); | |
596 | ||
597 | fprintf(fh, "%e\t", fabs(tmpAnalyUx-tmpAvgUx)); | |
598 | if (tmpAnalyUx != 0.0) { | |
599 | fprintf(fh, "%e\t", fabs(tmpAnalyUx - tmpAvgUx) / tmpAnalyUx); | |
600 | deviation += SQR(fabs(tmpAnalyUx - tmpAvgUx) / tmpAnalyUx); | |
601 | } | |
602 | else { | |
603 | fprintf(fh, "0.0\t"); | |
604 | } | |
605 | ||
606 | fprintf(fh, "%e\t", tmpAnalyUx / u_max); | |
607 | fprintf(fh, "%e\t", tmpAvgUx / u_max); | |
608 | fprintf(fh, "\n"); | |
609 | ||
610 | #undef SQR | |
611 | } | |
612 | ||
613 | *errorNorm = sqrt(deviation); | |
614 | ||
615 | printf("# Kernel validation: L2 error norm of relative error: %e\n", *errorNorm); | |
616 | ||
617 | ||
618 | fclose(fh); | |
619 | free(outputArray); | |
620 | ||
621 | ||
622 | } | |
623 | ||
624 | ||
625 | void KernelStatistics(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration) | |
626 | { | |
627 | KernelStatisticsAdv(kd, ld, cd, iteration, 0); | |
628 | } | |
629 | ||
630 | void KernelStatisticsAdv(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration, int forceOutput) | |
631 | { | |
632 | if (iteration % cd->StatisticsModulus == 0 || forceOutput) { | |
633 | printf("# iter: %4d avg density: %e\n", iteration, KernelDensity(kd, ld)); | |
634 | } | |
635 | ||
636 | if (iteration % 10 != 0 && !forceOutput) { | |
637 | return; | |
638 | } | |
639 | ||
640 | int nX = ld->Dims[0]; | |
641 | int nY = ld->Dims[1]; | |
642 | int nZ = ld->Dims[2]; | |
643 | ||
644 | int x = nX / 2; | |
645 | ||
646 | PdfT pdfs[N_D3Q19]; | |
647 | ||
648 | // ---------------------------------------------------------------------- | |
649 | // velocity in x-direction in cross section appended for each iteration | |
650 | ||
651 | double density; | |
652 | double densitySum; | |
653 | double ux; | |
654 | double uxSum = 0.0; | |
655 | int nFluidNodes = 0; | |
656 | ||
657 | for (int y = 0; y < nY; ++y) { | |
658 | for (int z = 0; z < nZ; ++z) { | |
659 | ||
660 | if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) { | |
661 | kd->GetNode(kd, x, y, z, pdfs); | |
662 | ||
663 | ux = pdfs[D3Q19_E] + pdfs[D3Q19_NE] + pdfs[D3Q19_SE] + pdfs[D3Q19_TE] + pdfs[D3Q19_BE] - | |
664 | pdfs[D3Q19_W] - pdfs[D3Q19_NW] - pdfs[D3Q19_SW] - pdfs[D3Q19_TW] - pdfs[D3Q19_BW]; | |
665 | ||
666 | uxSum += ux; | |
667 | ++nFluidNodes; | |
668 | } | |
669 | } | |
670 | } | |
671 | ||
672 | const char * mode = "w"; | |
673 | ||
674 | if (iteration > 0) { | |
675 | mode = "a"; | |
676 | } | |
677 | ||
678 | const char * fileName = "ux-progress.dat"; | |
679 | FILE * fh; | |
680 | ||
681 | fh = fopen(fileName, mode); | |
682 | ||
683 | if(fh == NULL) { | |
684 | printf("ERROR: opening file %s failed.\n", fileName); | |
685 | exit(1); | |
686 | } | |
687 | ||
688 | if (iteration == 0) { | |
689 | fprintf(fh, "# Average velocity in x direction of cross section in the middle (x = %d) of the geometry (NX = %d).\n", x, nX); | |
690 | fprintf(fh, "# Plot on terminal: gnuplot -e \"set terminal dumb; plot \\\"%s\\\";\"\n", fileName); | |
691 | fprintf(fh, "# iteration, avg ux\n"); | |
692 | } | |
693 | ||
694 | fprintf(fh, "%d %e\n", iteration, uxSum / nFluidNodes); | |
695 | ||
696 | fclose(fh); | |
697 | ||
698 | // ---------------------------------------------------------------------- | |
699 | // average velocity/density for each in cross section in x direction | |
700 | ||
701 | fileName = "density-ux.dat"; | |
702 | ||
703 | fh = fopen(fileName, "w"); | |
704 | ||
705 | if(fh == NULL) { | |
706 | printf("ERROR: opening file %s failed.\n", fileName); | |
707 | exit(1); | |
708 | } | |
709 | ||
710 | fprintf(fh, "# Average density and average x velocity over each cross section in x direction. Snapshot taken at iteration %d.\n", iteration); | |
711 | fprintf(fh, "# Plot on terminal: gnuplot -e \"set terminal dumb; plot \\\"%s\\\" u 1:2; plot \\\"%s\\\" u 1:3;\"\n", fileName, fileName); | |
712 | // fprintf(fh, "# Plot graphically: gnuplot -e \"plot \\\"%s\\\" u 1:3 w linesp t \\\"l\\\", \\\"\\\" u 1:4 w linesp t \\\"simulation\\\"; pause -1;" | |
713 | fprintf(fh, "# x, avg density, avg ux\n"); | |
714 | ||
715 | for (x = 0; x < nX; ++x) { | |
716 | ||
717 | uxSum = 0.0; | |
718 | densitySum = 0.0; | |
719 | nFluidNodes = 0; | |
720 | ||
721 | for (int y = 0; y < nY; ++y) { | |
722 | for (int z = 0; z < nZ; ++z) { | |
723 | ||
724 | if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] == LAT_CELL_OBSTACLE) { | |
725 | continue; | |
726 | } | |
727 | ||
728 | kd->GetNode(kd, x, y, z, pdfs); | |
729 | ||
730 | density = | |
731 | pdfs[D3Q19_C] + | |
732 | pdfs[D3Q19_N] + pdfs[D3Q19_E] + pdfs[D3Q19_S] + pdfs[D3Q19_W] + | |
733 | pdfs[D3Q19_NE] + pdfs[D3Q19_SE] + pdfs[D3Q19_SW] + pdfs[D3Q19_NW] + | |
734 | pdfs[D3Q19_T] + pdfs[D3Q19_TN] + pdfs[D3Q19_TE] + pdfs[D3Q19_TS] + pdfs[D3Q19_TW] + | |
735 | pdfs[D3Q19_B] + pdfs[D3Q19_BN] + pdfs[D3Q19_BE] + pdfs[D3Q19_BS] + pdfs[D3Q19_BW]; | |
736 | ||
737 | densitySum += density; | |
738 | ||
739 | ux = | |
740 | pdfs[D3Q19_E] + pdfs[D3Q19_NE] + pdfs[D3Q19_SE] + pdfs[D3Q19_TE] + pdfs[D3Q19_BE] - | |
741 | pdfs[D3Q19_W] - pdfs[D3Q19_NW] - pdfs[D3Q19_SW] - pdfs[D3Q19_TW] - pdfs[D3Q19_BW]; | |
742 | ||
743 | uxSum += ux; | |
744 | ||
745 | ++nFluidNodes; | |
746 | } | |
747 | } | |
748 | ||
749 | fprintf(fh, "%d %e %e\n", x, densitySum / nFluidNodes, uxSum / nFluidNodes); | |
750 | } | |
751 | ||
752 | fclose(fh); | |
753 | } | |
754 | ||
755 | ||
756 | ||
757 | void KernelAddBodyForce(KernelData * kd, LatticeDesc * ld, CaseData * cd) | |
758 | { | |
759 | Assert(kd != NULL); | |
760 | Assert(ld != NULL); | |
761 | Assert(cd != NULL); | |
762 | ||
763 | int nX = kd->Dims[0]; | |
764 | int nY = kd->Dims[1]; | |
765 | int nZ = kd->Dims[2]; | |
766 | ||
767 | PdfT w_0 = 1.0 / 3.0; // C | |
768 | PdfT w_1 = 1.0 / 18.0; // N,S,E,W,T,B | |
769 | PdfT w_2 = 1.0 / 36.0; // NE,NW,SE,SW,TE,TW,BE,BW,TN,TS,BN,BS | |
770 | PdfT w[] = {w_1,w_1,w_1,w_1,w_2,w_2,w_2,w_2,w_1,w_2,w_2,w_2,w_2,w_1,w_2,w_2,w_2,w_2,w_0}; | |
771 | ||
772 | PdfT xForce = cd->XForce; | |
773 | ||
774 | PdfT pdfs[N_D3Q19]; | |
775 | ||
776 | ||
777 | #ifdef _OPENMP | |
778 | #pragma omp parallel for collapse(3) default(none) \ | |
779 | shared(nX,nY,nZ,ld,kd,w,xForce,D3Q19_X,cd) \ | |
780 | private(pdfs) | |
781 | #endif | |
782 | for(int z = 0; z < nZ; ++z) { | |
783 | for(int y = 0; y < nY; ++y) { | |
784 | for(int x = 0; x < nX; ++x) { | |
785 | if(ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] == LAT_CELL_OBSTACLE) | |
786 | continue; | |
787 | ||
788 | // load pdfs into temp array | |
789 | kd->GetNode(kd, x, y, z, pdfs); | |
790 | ||
791 | // add body force in x direction ( method by Luo) | |
792 | for (int d = 0; d < N_D3Q19; ++d) { | |
793 | pdfs[d] = pdfs[d] + 3.0*w[d]*D3Q19_X[d]*xForce; | |
794 | } | |
795 | ||
796 | kd->SetNode(kd, x, y, z, pdfs); | |
797 | ||
798 | } | |
799 | } | |
800 | } | |
801 | } |