src/Kernel.c

   1 // --------------------------------------------------------------------------
   2 //
   3 // Copyright
   4 //   Markus Wittmann, 2016-2017
   5 //   RRZE, University of Erlangen-Nuremberg, Germany
   6 //   markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
   7 //
   8 //   Viktor Haag, 2016
   9 //   LSS, University of Erlangen-Nuremberg, Germany
  10 //
  11 //  This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
  12 //
  13 //  LbmBenchKernels is free software: you can redistribute it and/or modify
  14 //  it under the terms of the GNU General Public License as published by
  15 //  the Free Software Foundation, either version 3 of the License, or
  16 //  (at your option) any later version.
  17 //
  18 //  LbmBenchKernels is distributed in the hope that it will be useful,
  19 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
  20 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  21 //  GNU General Public License for more details.
  22 //
  23 //  You should have received a copy of the GNU General Public License
  24 //  along with LbmBenchKernels.  If not, see <http://www.gnu.org/licenses/>.
  25 //
  26 // --------------------------------------------------------------------------
  27 #include "Kernel.h"
  28 #include "Lattice.h"
  29
  30 #include <stdlib.h>
  31 #include <string.h>
  32 #include <math.h>
  33
  34 #define X(name, idx, idx_inv, x, y, z)  , x
  35 int D3Q19_X[] = {
  36         EXPAND(D3Q19_LIST)
  37 };
  38 #undef X
  39
  40 #define X(name, idx, idx_inv, x, y, z)  , y
  41 int D3Q19_Y[] = {
  42         EXPAND(D3Q19_LIST)
  43 };
  44 #undef X
  45
  46 #define X(name, idx, idx_inv, x, y, z)  , z
  47 int D3Q19_Z[] = {
  48         EXPAND(D3Q19_LIST)
  49 };
  50 #undef X
  51
  52 #define X(name, idx, idxinv, x, y, z)   , idxinv
  53 int D3Q19_INV[] = {
  54         EXPAND(D3Q19_LIST)
  55 };
  56 #undef X
  57
  58
  59 #define X(name, idx, idxinv, x, y, z)   , STRINGIFY(name)
  60 const char * D3Q19_NAMES[N_D3Q19] = {
  61         EXPAND(D3Q19_LIST)
  62 };
  63 #undef X
  64
  65 void KernelComputeBoundaryConditions(KernelData * kd, LatticeDesc * ld, CaseData * cd)
  66 {
  67         Assert(kd != NULL);
  68         Assert(ld != NULL);
  69         Assert(cd != NULL);
  70
  71         Assert(cd->RhoIn  > F(0.0));
  72         Assert(cd->RhoOut > F(0.0));
  73
  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);
  79
  80         PdfT dens;
  81         PdfT ux;
  82
  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);
  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 = 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;
 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 = 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;
 184
 185                                 pdfs[D3Q19_W ] = pdfs[D3Q19_E] - one_third * dens * ux;
 186                                 pdfs[D3Q19_SW] = pdfs[D3Q19_NE] + one_fourth * (pdfs[D3Q19_N] - pdfs[D3Q19_S]) - indep_ux;
 187                                 pdfs[D3Q19_NW] = pdfs[D3Q19_SE] - one_fourth * (pdfs[D3Q19_N] - pdfs[D3Q19_S]) - indep_ux;
 188                                 pdfs[D3Q19_BW] = pdfs[D3Q19_TE] + one_fourth * (pdfs[D3Q19_T] - pdfs[D3Q19_B]) - indep_ux;
 189                                 pdfs[D3Q19_TW] = pdfs[D3Q19_BE] - one_fourth * (pdfs[D3Q19_T] - pdfs[D3Q19_B]) - indep_ux;
 190
 191                                 kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_W , pdfs[D3Q19_W ]);
 192                                 kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_NW, pdfs[D3Q19_NW]);
 193                                 kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_SW, pdfs[D3Q19_SW]);
 194                                 kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_TW, pdfs[D3Q19_TW]);
 195                                 kd->BoundaryConditionsSetPdf(kd, x, y, z, D3Q19_BW, pdfs[D3Q19_BW]);
 196
 197                                 for(int d = 0; d < N_D3Q19; ++d) {
 198                                         density_out += pdfs[d];
 199                                 }
 200                         }
 201                 }
 202         }
 203
 204         // DEBUG: printf("# density inlet: %e  density outlet: %e\n", density_in, density_out);
 205
 206 }
 207
 208
 209 PdfT KernelDensity(KernelData * kd, LatticeDesc * ld)
 210 {
 211         Assert(kd != NULL);
 212         Assert(ld != NULL);
 213
 214         Assert(ld->Lattice != NULL);
 215         Assert(ld->Dims    != NULL);
 216
 217         Assert(ld->Dims[0] > 0);
 218         Assert(ld->Dims[1] > 0);
 219         Assert(ld->Dims[2] > 0);
 220
 221         int * lDims = ld->Dims;
 222         int nX = lDims[0];
 223         int nY = lDims[1];
 224         int nZ = lDims[2];
 225
 226         PdfT pdfs[N_D3Q19] = { -1.0 };
 227         PdfT density = 0.0;
 228
 229         for(int z = 0; z < nZ; ++z) {
 230                 for(int y = 0; y < nY; ++y) {
 231                         for(int x = 0; x < nX; ++x) {
 232
 233                                 if(ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) {
 234
 235                                         kd->GetNode(kd, x, y, z, pdfs);
 236
 237                                         PdfT localDensity = F(0.0);
 238
 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 //                                              }
 244                                                 localDensity += pdfs[d];
 245                                         }
 246                                         density += localDensity;
 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
 265         PdfT ux   = F(0.0);
 266         PdfT uy   = F(0.0);
 267         PdfT uz   = F(0.0);
 268         PdfT dens = F(1.0);
 269
 270         PdfT omega = cd->Omega;
 271
 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);
 275
 276         PdfT dir_indep_trm;
 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);
 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))
 296                         dens = rho_in + (rho_out - rho_in) * (x) / (nX - F(1.0));
 297
 298                         #define SQR(a) ((a)*(a))
 299
 300                         dir_indep_trm = one_third * dens - F(0.5) * (ux * ux + uy * uy + uz * uz);
 301
 302                         pdfs[D3Q19_C]  = omega_w0 * (dir_indep_trm);
 303
 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));
 306
 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));
 309
 310
 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));
 313
 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));
 316
 317
 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));
 320
 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));
 323
 324
 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));
 327
 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));
 330
 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));
 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
 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);
 354
 355         PdfT omega = cd->Omega;
 356
 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);
 360
 361         PdfT dir_indep_trm;
 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);
 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
 382                         ux = F(0.0);
 383                         uy = F(0.0);
 384                         uz = F(0.0);
 385
 386                         kd->GetNode(kd, x, y, z, pdfs);
 387
 388
 389                         density = F(0.0);
 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))
 397                         dir_indep_trm = one_third * dens - F(0.5) * (ux * ux + uy * uy + uz * uz);
 398
 399                         pdfs[D3Q19_C]  = omega_w0 * (dir_indep_trm);
 400
 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));
 403
 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));
 406
 407
 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));
 410
 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));
 413
 414
 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));
 417
 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));
 420
 421
 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));
 424
 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));
 427
 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));
 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 {
 453         return xForce * (maxRadiusSquared - curRadiusSquared) / (F(2.0) * viscosity);
 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
 468         PdfT ux = F(0.0);
 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
 492                                 ux += F(0.5) * cd->XForce;
 493                                 #endif
 494
 495                                 outputArray[y * nZ + z] = ux;
 496                         }
 497                         else {
 498                                 outputArray[y * nZ + z] = F(0.0);
 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;
 523         PdfT viscosity        = (F(1.0) / omega - F(0.5)) / F(3.0);
 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];
 538         PdfT tmpAvgUx   = F(0.0);
 539         PdfT tmpAnalyUx = F(0.0);
 540         int flagEvenNy = 0;
 541         int y = 0;
 542
 543         if (nY % 2 == 0)
 544                 flagEvenNy = 1;
 545
 546         y = (nY - flagEvenNy - 1) / 2;
 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
 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));
 571         #undef SQR
 572         PdfT u_max = cd->XForce*minRadiusSquared / (F(2.0) * viscosity);
 573
 574         for(int z = 0; z < nZ; ++z) {
 575
 576                 fprintf(fh, "%d\t", z);
 577
 578                 #define SQR(a) ((a)*(a))
 579                 curRadiusSquared = SQR(z - center + F(0.5));
 580
 581
 582                 // dimensionless radius
 583                 fprintf(fh, "%e\t", (z - center + F(0.5)) / center);
 584
 585                 // analytic profile
 586                 if(curRadiusSquared >= minRadiusSquared)
 587                         tmpAnalyUx = F(0.0);
 588                 else
 589                         tmpAnalyUx = CalcXVelForPipeProfile(minRadiusSquared, curRadiusSquared, cd->XForce, viscosity);
 590
 591                 //averaged profile
 592                 if(flagEvenNy == 1)
 593                         tmpAvgUx = (outputArray[y * nZ + z] + outputArray[(y + 1) * nZ + z]) / F(2.0);
 594                 else
 595                         tmpAvgUx =  outputArray[y * nZ + z];
 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);
 603                         deviation += SQR((PdfT)fabs(tmpAnalyUx - tmpAvgUx) / tmpAnalyUx);
 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
 616         *errorNorm = (PdfT)sqrt(deviation);
 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);
 715         fprintf(fh, "# x, avg density, avg ux\n");
 716
 717         for (x = 0; x < nX; ++x) {
 718
 719                 uxSum = F(0.0);
 720                 densitySum = F(0.0);
 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
 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};
 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
 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;
 796                                 }
 797
 798                                 kd->SetNode(kd, x, y, z, pdfs);
 799
 800                         }
 801                 }
 802         }
 803 }