1 // --------------------------------------------------------------------------
4 // Markus Wittmann, 2016-2017
5 // RRZE, University of Erlangen-Nuremberg, Germany
6 // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
9 // LSS, University of Erlangen-Nuremberg, Germany
11 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
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.
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.
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/>.
26 // --------------------------------------------------------------------------
34 #define X(name, idx, idx_inv, x, y, z) , x
40 #define X(name, idx, idx_inv, x, y, z) , y
46 #define X(name, idx, idx_inv, x, y, z) , z
52 #define X(name, idx, idxinv, x, y, z) , idxinv
59 #define X(name, idx, idxinv, x, y, z) , STRINGIFY(name)
60 const char * D3Q19_NAMES[N_D3Q19] = {
65 void KernelComputeBoundaryConditions(KernelData * kd, LatticeDesc * ld, CaseData * cd)
71 Assert(cd->RhoIn > F(0.0));
72 Assert(cd->RhoOut > F(0.0));
74 PdfT rho_in = cd->RhoIn;
75 PdfT rho_out = cd->RhoOut;
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);
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);
97 double density_in = 0.0;
98 double density_out = 0.0;
100 // update inlet / outlet boundary conditions
101 for (int z = 1; z < nZ - 1; ++z) {
102 for (int y = 1; y < nY - 1; ++y) {
105 // -----------------------------------------------------------------------------
106 // update inlet conditions
108 if (ld->Lattice[L_INDEX_4(ld->Dims, x_in, y, z)] == LAT_CELL_INLET) {
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);
129 ux = F(1.0) - (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 F(2.0) * (pdfs[D3Q19_SW] + pdfs[D3Q19_TW] + pdfs[D3Q19_W] + pdfs[D3Q19_BW] + pdfs[D3Q19_NW])) * rho_in_inv;
134 indep_ux = one_sixth * dens * ux;
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;
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]);
149 for(int d = 0; d < N_D3Q19; ++d) {
150 density_in += pdfs[d];
154 // -----------------------------------------------------------------------------
155 // update outlet conditions
157 if (ld->Lattice[L_INDEX_4(ld->Dims, x_out, y, z)] == LAT_CELL_OUTLET) {
158 // update outlet conditions
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);
179 ux = F(-1.0) + (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 F(2.0) * (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;
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;
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]);
197 for(int d = 0; d < N_D3Q19; ++d) {
198 density_out += pdfs[d];
204 // DEBUG: printf("# density inlet: %e density outlet: %e\n", density_in, density_out);
209 PdfT KernelDensity(KernelData * kd, LatticeDesc * ld)
214 Assert(ld->Lattice != NULL);
215 Assert(ld->Dims != NULL);
217 Assert(ld->Dims[0] > 0);
218 Assert(ld->Dims[1] > 0);
219 Assert(ld->Dims[2] > 0);
221 int * lDims = ld->Dims;
226 PdfT pdfs[N_D3Q19] = { -1.0 };
229 for(int z = 0; z < nZ; ++z) {
230 for(int y = 0; y < nY; ++y) {
231 for(int x = 0; x < nX; ++x) {
233 if(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) {
235 kd->GetNode(kd, x, y, z, pdfs);
237 PdfT localDensity = F(0.0);
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]);
244 localDensity += pdfs[d];
246 density += localDensity;
253 return density / ld->nFluid;
257 // prescribes a given density
258 void KernelSetInitialDensity(LatticeDesc * ld, KernelData * kd, CaseData * cd)
260 int * lDims = ld->Dims;
262 PdfT rho_in = cd->RhoIn;
263 PdfT rho_out = cd->RhoOut;
270 PdfT omega = cd->Omega;
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);
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);
289 #pragma omp parallel for collapse(3)
291 for(int z = 0; z < nZ; ++z) { for(int y = 0; y < nY; ++y) { for(int x = 0; x < nX; ++x) {
293 if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
295 // if((caseData->geoType == GEO_TYPE_CHANNEL) || (caseData->geoType == GEO_TYPE_RCHANNEL))
296 dens = rho_in + (rho_out - rho_in) * (x) / (nX - F(1.0));
298 #define SQR(a) ((a)*(a))
300 dir_indep_trm = one_third * dens - F(0.5) * (ux * ux + uy * uy + uz * uz);
302 pdfs[D3Q19_C] = omega_w0 * (dir_indep_trm);
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));
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));
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));
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));
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));
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));
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));
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));
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));
335 kd->SetNode(kd, x, y, z, pdfs);
343 // prescribes a given velocity
344 void KernelSetInitialVelocity(LatticeDesc * ld, KernelData * kd, CaseData * cd)
347 int * lDims = ld->Dims;
349 // TODO: fix ux is overriden below
355 PdfT omega = cd->Omega;
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);
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);
376 #pragma omp parallel for collapse(3)
378 for(int z = 0; z < nZ; ++z) { for(int y = 0; y < nY; ++y) { for(int x = 0; x < nX; ++x) {
380 if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] == LAT_CELL_FLUID) {
386 kd->GetNode(kd, x, y, z, pdfs);
391 #define X(name, idx, idxinv, _x, _y, _z) density += pdfs[idx];
396 #define SQR(a) ((a)*(a))
397 dir_indep_trm = one_third * dens - F(0.5) * (ux * ux + uy * uy + uz * uz);
399 pdfs[D3Q19_C] = omega_w0 * (dir_indep_trm);
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));
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));
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));
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));
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));
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));
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));
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));
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));
434 kd->SetNode(kd, x, y, z, pdfs);
440 // Compute analytical x velocity for channel flow.
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.
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.
449 // calculate velocity in a pipe for a given radius
451 static PdfT CalcXVelForPipeProfile(PdfT maxRadiusSquared, PdfT curRadiusSquared, PdfT xForce, PdfT viscosity)
453 return xForce * (maxRadiusSquared - curRadiusSquared) / (F(2.0) * viscosity);
456 static void KernelGetXSlice(LatticeDesc * ld, KernelData * kd, CaseData * cd, PdfT * outputArray, int xPos)
461 int nY = ld->Dims[1];
462 int nZ = ld->Dims[2];
465 Assert(xPos < ld->Dims[0]);
470 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
471 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
476 for(int z = 0; z < nZ; ++z) {
477 for(int y = 0; y < nY; ++y) {
479 if (ld->Lattice[L_INDEX_4(ld->Dims, xPos, y, z)] != LAT_CELL_OBSTACLE) {
480 kd->GetNode(kd, xPos, y, z, pdfs);
482 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = pdfs[idx];
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);
488 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
489 pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
492 ux += F(0.5) * cd->XForce;
495 outputArray[y * nZ + z] = ux;
498 outputArray[y * nZ + z] = F(0.0);
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
511 void KernelVerifiy(LatticeDesc * ld, KernelData * kd, CaseData * cd, PdfT * errorNorm)
516 Assert(errorNorm != NULL);
518 int nX = ld->Dims[0];
519 int nY = ld->Dims[1];
520 int nZ = ld->Dims[2];
522 PdfT omega = cd->Omega;
523 PdfT viscosity = (F(1.0) / omega - F(0.5)) / F(3.0);
525 // ux averaged across cross sections in x direction
526 PdfT * outputArray = (PdfT *)malloc(nZ * nY * sizeof(PdfT));
527 Verify(outputArray != NULL);
529 memset(outputArray, -10, nZ*nY*sizeof(PdfT));
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));
538 PdfT tmpAvgUx = F(0.0);
539 PdfT tmpAnalyUx = F(0.0);
546 y = (nY - flagEvenNy - 1) / 2;
548 snprintf(fileName, sizeof(fileName), "flow-profile.dat");
550 printf("# Kernel validation: writing profile to %s\n", fileName);
552 fh = fopen(fileName, "w");
555 printf("ERROR: opening file %s failed.\n", fileName);
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");
565 PdfT deviation = F(0.0);
566 PdfT curRadiusSquared;
567 PdfT center = nY / F(2.0);
568 PdfT minDiameter = (PdfT)nY;
569 #define SQR(a) ((a)*(a))
570 PdfT minRadiusSquared = SQR(minDiameter / F(2.0) - F(1.0));
572 PdfT u_max = cd->XForce*minRadiusSquared / (F(2.0) * viscosity);
574 for(int z = 0; z < nZ; ++z) {
576 fprintf(fh, "%d\t", z);
578 #define SQR(a) ((a)*(a))
579 curRadiusSquared = SQR(z - center + F(0.5));
582 // dimensionless radius
583 fprintf(fh, "%e\t", (z - center + F(0.5)) / center);
586 if(curRadiusSquared >= minRadiusSquared)
589 tmpAnalyUx = CalcXVelForPipeProfile(minRadiusSquared, curRadiusSquared, cd->XForce, viscosity);
593 tmpAvgUx = (outputArray[y * nZ + z] + outputArray[(y + 1) * nZ + z]) / F(2.0);
595 tmpAvgUx = outputArray[y * nZ + z];
597 fprintf(fh, "%e\t", tmpAnalyUx);
598 fprintf(fh, "%e\t", tmpAvgUx);
600 fprintf(fh, "%e\t", fabs(tmpAnalyUx-tmpAvgUx));
601 if (tmpAnalyUx != 0.0) {
602 fprintf(fh, "%e\t", fabs(tmpAnalyUx - tmpAvgUx) / tmpAnalyUx);
603 deviation += SQR((PdfT)fabs(tmpAnalyUx - tmpAvgUx) / tmpAnalyUx);
606 fprintf(fh, "0.0\t");
609 fprintf(fh, "%e\t", tmpAnalyUx / u_max);
610 fprintf(fh, "%e\t", tmpAvgUx / u_max);
616 *errorNorm = (PdfT)sqrt(deviation);
618 printf("# Kernel validation: L2 error norm of relative error: %e\n", *errorNorm);
628 void KernelStatistics(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration)
630 KernelStatisticsAdv(kd, ld, cd, iteration, 0);
633 void KernelStatisticsAdv(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration, int forceOutput)
635 if (iteration % cd->StatisticsModulus == 0 || forceOutput) {
636 printf("# iter: %4d avg density: %e\n", iteration, KernelDensity(kd, ld));
639 if (iteration % 10 != 0 && !forceOutput) {
643 int nX = ld->Dims[0];
644 int nY = ld->Dims[1];
645 int nZ = ld->Dims[2];
651 // ----------------------------------------------------------------------
652 // velocity in x-direction in cross section appended for each iteration
660 for (int y = 0; y < nY; ++y) {
661 for (int z = 0; z < nZ; ++z) {
663 if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
664 kd->GetNode(kd, x, y, z, pdfs);
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];
675 const char * mode = "w";
681 const char * fileName = "ux-progress.dat";
684 fh = fopen(fileName, mode);
687 printf("ERROR: opening file %s failed.\n", fileName);
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");
697 fprintf(fh, "%d %e\n", iteration, uxSum / nFluidNodes);
701 // ----------------------------------------------------------------------
702 // average velocity/density for each in cross section in x direction
704 fileName = "density-ux.dat";
706 fh = fopen(fileName, "w");
709 printf("ERROR: opening file %s failed.\n", fileName);
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);
715 fprintf(fh, "# x, avg density, avg ux\n");
717 for (x = 0; x < nX; ++x) {
723 for (int y = 0; y < nY; ++y) {
724 for (int z = 0; z < nZ; ++z) {
726 if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] == LAT_CELL_OBSTACLE) {
730 kd->GetNode(kd, x, y, z, pdfs);
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];
739 densitySum += density;
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];
751 fprintf(fh, "%d %e %e\n", x, densitySum / nFluidNodes, uxSum / nFluidNodes);
759 void KernelAddBodyForce(KernelData * kd, LatticeDesc * ld, CaseData * cd)
765 int nX = kd->Dims[0];
766 int nY = kd->Dims[1];
767 int nZ = kd->Dims[2];
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
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};
774 PdfT xForce = cd->XForce;
780 #pragma omp parallel for collapse(3) default(none) \
781 shared(nX,nY,nZ,ld,kd,w,xForce,D3Q19_X,cd) \
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)
790 // load pdfs into temp array
791 kd->GetNode(kd, x, y, z, pdfs);
793 // add body force in x direction (method by Luo)
794 for (int d = 0; d < N_D3Q19; ++d) {
795 pdfs[d] = pdfs[d] + F(3.0) * w[d] * D3Q19_X[d] * xForce;
798 kd->SetNode(kd, x, y, z, pdfs);