projects / LbmBenchmarkKernelsPublic.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
add citation information
[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
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
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
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
258void 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
344void 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//
451static 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
456static 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//
511void 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
628void KernelStatistics(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration)
629{
630 KernelStatisticsAdv(kd, ld, cd, iteration, 0);
631}
632
633void 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
759void 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}
LBM benchmark kernels
This page took 0.14033 seconds and 5 git commands to generate.