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  > 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
 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                                         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
 255 void 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
 341 void 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 //
 448 static PdfT CalcXVelForPipeProfile(PdfT maxRadiusSquared, PdfT curRadiusSquared, PdfT xForce, PdfT viscosity)
 449 {
 450         return xForce*(maxRadiusSquared - curRadiusSquared) / (2.0*viscosity);
 451 }
 452
 453 static 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 //
 508 void 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
 625 void KernelStatistics(KernelData * kd, LatticeDesc * ld, CaseData * cd, int iteration)
 626 {
 627         KernelStatisticsAdv(kd, ld, cd, iteration, 0);
 628 }
 629
 630 void 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
 757 void 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 }