projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
version 0.1
[LbmBenchmarkKernelsPublic.git] / src / Kernel.c
CommitLineData
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
35int D3Q19_X[] = {
36 EXPAND(D3Q19_LIST)
37};
38#undef X
39
40#define X(name, idx, idx_inv, x, y, z) , y
41int D3Q19_Y[] = {
42 EXPAND(D3Q19_LIST)
43};
44#undef X
45
46#define X(name, idx, idx_inv, x, y, z) , z
47int D3Q19_Z[] = {
48 EXPAND(D3Q19_LIST)
49};
50#undef X
51
52#define X(name, idx, idxinv, x, y, z) , idxinv
53int D3Q19_INV[] = {
54 EXPAND(D3Q19_LIST)
55};
56#undef X
57
58
59#define X(name, idx, idxinv, x, y, z) , STRINGIFY(name)
60const char * D3Q19_NAMES[N_D3Q19] = {
61 EXPAND(D3Q19_LIST)
62};
63#undef X
64
65void 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
209PdfT 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
255void 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
341void 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//
448static PdfT CalcXVelForPipeProfile(PdfT maxRadiusSquared, PdfT curRadiusSquared, PdfT xForce, PdfT viscosity)
449{
450 return xForce*(maxRadiusSquared - curRadiusSquared) / (2.0*viscosity);
451}
452
453static 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//
508void 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
625void KernelStatistics(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration)
626{
627 KernelStatisticsAdv(kd, ld, cd, iteration, 0);
628}
629
630void 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
757void 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}
LBM benchmark kernels
This page took 0.121938 seconds and 5 git commands to generate.