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