src/BenchKernelD3Q19ListAa.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 "BenchKernelD3Q19ListAaCommon.h"
  28
  29 #include "Memory.h"
  30 #include "Vtk.h"
  31 #include "LikwidIf.h"
  32
  33 #include <inttypes.h>
  34 #include <math.h>
  35
  36 void FNAME(D3Q19ListAaKernel)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
  37 {
  38         Assert(ld != NULL);
  39         Assert(kernelData != NULL);
  40         Assert(cd != NULL);
  41
  42         Assert(cd->Omega > F(0.0));
  43         Assert(cd->Omega < F(2.0));
  44
  45         KernelData * kd = (KernelData *)kernelData;
  46         KernelDataList * kdl = (KernelDataList *)kernelData;
  47
  48         PdfT omega = cd->Omega;
  49         PdfT omegaEven = omega;
  50         //  1/4: best stability;
  51         // 1/12: removes third-order advection error (best advection);
  52         //  1/6: removes fourth-order diffusion error (best diffusion);
  53         // 3/16: exact location of bounce back for poiseuille flow
  54         PdfT magicParam = F(1.0)/F(12.0);
  55         PdfT omegaOdd = F(1.0)/( F(0.5) + magicParam/(F(1.0)/omega - F(0.5)) );
  56
  57         PdfT evenPart = F(0.0);
  58         PdfT oddPart = F(0.0);
  59         PdfT dir_indep_trm = F(0.0);
  60
  61         const PdfT w_0 = F(1.0) /  F(3.0);
  62         const PdfT w_1 = F(1.0) / F(18.0);
  63         const PdfT w_2 = F(1.0) / F(36.0);
  64
  65         const PdfT w_1_x3 = w_1 * F(3.0);       const PdfT w_1_nine_half = w_1 * F(9.0)/F(2.0); PdfT w_1_indep = F(0.0);
  66         const PdfT w_2_x3 = w_2 * F(3.0);       const PdfT w_2_nine_half = w_2 * F(9.0)/F(2.0); PdfT w_2_indep = F(0.0);
  67
  68         PdfT ui;
  69
  70         PdfT ux, uy, uz;
  71         PdfT dens;
  72
  73         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
  74         #define X(name, idx, idxinv, x, y, z)   PdfT JOIN(pdf_,name);
  75         D3Q19_LIST
  76         #undef X
  77
  78         PdfT * src = kd->Pdfs[0];
  79
  80         int maxIterations = cd->MaxIterations;
  81         int nFluid = kdl->nFluid;
  82         int nCells = kdl->nCells;
  83
  84         uint32_t adjListIndex;
  85         uint32_t * adjList = kdl->AdjList;
  86
  87
  88
  89         #ifdef VTK_OUTPUT
  90                 if (cd->VtkOutput) {
  91                         kd->PdfsActive = src;
  92                         VtkWrite(ld, kd, cd, -1);
  93                 }
  94         #endif
  95
  96         #ifdef STATISTICS
  97                 kd->PdfsActive = src;
  98                 KernelStatistics(kd, ld, cd, 0);
  99         #endif
 100
 101         X_KERNEL_START(kernelData);
 102
 103         // TODO: outer openmp parallel
 104         for(int iter = 0; iter < maxIterations; iter += 2) {
 105
 106
 107                 // --------------------------------------------------------------------
 108                 // even time step
 109
 110                 X_LIKWID_START("list-aa-even");
 111
 112 #ifdef _OPENMP
 113                 #pragma omp parallel for default(none) \
 114                                 shared(nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \
 115                                 private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\
 116                                         pdf_C, \
 117                                         pdf_N, pdf_E, pdf_S, pdf_W, \
 118                                         pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
 119                                         pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
 120                                         pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW)
 121 #endif
 122 #ifdef INTEL_OPT_DIRECTIVES
 123                 #pragma ivdep
 124                 #pragma vector always
 125                 #pragma simd
 126 #endif
 127                 for (int index = 0; index < nFluid; ++index) { // LOOP list-aa-even
 128
 129
 130                         #define I(index, dir)   P_INDEX_3((nCells), (index), (dir))
 131
 132                         // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
 133                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = src[I(index, idx)];
 134                         D3Q19_LIST
 135                         #undef X
 136
 137 // #define LID_DRIVEN_CAVITY
 138
 139 #ifdef LID_DRIVEN_CAVITY
 140                         int nX = kd->Dims[0];
 141                         int nY = kd->Dims[1];
 142                         int nZ = kd->Dims[2];
 143
 144                         int x = kdl->Coords[C_INDEX_X(index)];
 145                         int y = kdl->Coords[C_INDEX_Y(index)];
 146                         int z = kdl->Coords[C_INDEX_Z(index)];
 147
 148                         if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
 149                                 ux = 0.1 * 0.577;
 150                                 uy = 0.0;
 151                                 uz = 0.0;
 152                         } else {
 153 #endif
 154                                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 155                                          pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 156                                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 157                                          pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 158                                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 159                                          pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 160 #ifdef LID_DRIVEN_CAVITY
 161                         }
 162 #endif
 163
 164                         dens = pdf_C +
 165                                    pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 166                                    pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 167                                    pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 168                                    pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 169
 170                         dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0);
 171
 172                         // direction: w_0
 173                         src[I(index, D3Q19_C)             ]  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
 174
 175                         // direction: w_1
 176                         w_1_indep = w_1*dir_indep_trm;
 177
 178                         ui = uy;
 179                         evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
 180                         oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
 181                         src[I(index, D3Q19_S)]  = pdf_N - evenPart - oddPart;
 182                         src[I(index, D3Q19_N)]  = pdf_S - evenPart + oddPart;
 183
 184                         ui = ux;
 185                         evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
 186                         oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
 187                         src[I(index, D3Q19_W)]  = pdf_E - evenPart - oddPart;
 188                         src[I(index, D3Q19_E)]  = pdf_W - evenPart + oddPart;
 189
 190                         ui = uz;
 191                         evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
 192                         oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
 193                         src[I(index, D3Q19_B)]  = pdf_T - evenPart - oddPart;
 194                         src[I(index, D3Q19_T)]  = pdf_B - evenPart + oddPart;
 195
 196                         // direction: w_2
 197                         w_2_indep = w_2*dir_indep_trm;
 198
 199                         ui = -ux + uy;
 200                         evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
 201                         oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
 202                         src[I(index, D3Q19_SE)] = pdf_NW - evenPart - oddPart;
 203                         src[I(index, D3Q19_NW)] = pdf_SE - evenPart + oddPart;
 204
 205                         ui = ux + uy;
 206                         evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
 207                         oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
 208                         src[I(index, D3Q19_SW)] = pdf_NE - evenPart - oddPart;
 209                         src[I(index, D3Q19_NE)] = pdf_SW - evenPart + oddPart;
 210
 211                         ui = -ux + uz;
 212                         evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
 213                         oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
 214                         src[I(index, D3Q19_BE)] = pdf_TW - evenPart - oddPart;
 215                         src[I(index, D3Q19_TW)] = pdf_BE - evenPart + oddPart;
 216
 217                         ui = ux + uz;
 218                         evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
 219                         oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
 220                         src[I(index, D3Q19_BW)] = pdf_TE - evenPart - oddPart;
 221                         src[I(index, D3Q19_TE)] = pdf_BW - evenPart + oddPart;
 222
 223                         ui = -uy + uz;
 224                         evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
 225                         oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
 226                         src[I(index, D3Q19_BN)] = pdf_TS - evenPart - oddPart;
 227                         src[I(index, D3Q19_TS)] = pdf_BN - evenPart + oddPart;
 228
 229                         ui = uy + uz;
 230                         evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
 231                         oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
 232                         src[I(index, D3Q19_BS)] = pdf_TN - evenPart - oddPart;
 233                         src[I(index, D3Q19_TN)] = pdf_BS - evenPart + oddPart;
 234
 235                 } // loop over fluid nodes
 236
 237                 X_LIKWID_STOP("list-aa-even");
 238
 239                 // save current iteration
 240                 kdl->Iteration = iter;
 241                 #ifdef VERIFICATION
 242                         kd->PdfsActive = src;
 243                         KernelAddBodyForce(kd, ld, cd);
 244                 #endif
 245
 246                 // --------------------------------------------------------------------
 247                 // odd time step
 248
 249                 X_LIKWID_START("list-aa-odd");
 250
 251 #ifdef _OPENMP
 252                 #pragma omp parallel for default(none) \
 253                                 shared(nFluid, nCells, kd, kdl, adjList, omegaOdd, omegaEven, src) \
 254                                 private(ux, uy, uz, dens, adjListIndex, evenPart, oddPart, dir_indep_trm, w_1_indep, w_2_indep, ui,\
 255                                         pdf_C, \
 256                                         pdf_N, pdf_E, pdf_S, pdf_W, \
 257                                         pdf_NE, pdf_SE, pdf_SW, pdf_NW, \
 258                                         pdf_T, pdf_TN, pdf_TE, pdf_TS, pdf_TW, \
 259                                         pdf_B, pdf_BN, pdf_BE, pdf_BS, pdf_BW)
 260 #endif
 261 #ifdef INTEL_OPT_DIRECTIVES
 262                 #pragma ivdep
 263 #endif
 264                 for (int index = 0; index < nFluid; ++index) { // LOOP list-aa-odd
 265
 266
 267                         adjListIndex = index * N_D3Q19_IDX;
 268
 269                         // Load PDFs of local cell: pdf_N = src[adjList[adjListIndex + D3Q19_S]]; ...
 270                         pdf_C = src[P_INDEX_3(nCells, index, D3Q19_C)];
 271
 272                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = src[adjList[adjListIndex + idxinv]];
 273                         D3Q19_LIST_WO_C
 274                         #undef X
 275
 276 #ifdef LID_DRIVEN_CAVITY
 277                         int nX = kd->Dims[0];
 278                         int nY = kd->Dims[1];
 279                         int nZ = kd->Dims[2];
 280
 281                         int x = kdl->Coords[C_INDEX_X(index)];
 282                         int y = kdl->Coords[C_INDEX_Y(index)];
 283                         int z = kdl->Coords[C_INDEX_Z(index)];
 284
 285                         if (z == nZ - 4 && x > 3 && x < (nX - 4) && y > 3 && y < (nY - 4)) {
 286                                 ux = F(0.1) * F(0.5)77;
 287                                 uy = F(0.0);
 288                                 uz = F(0.0);
 289                         } else {
 290 #endif
 291                                 ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 292                                          pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 293                                 uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 294                                          pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 295                                 uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 296                                          pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 297 #ifdef LID_DRIVEN_CAVITY
 298                         }
 299 #endif
 300
 301                         dens = pdf_C +
 302                                    pdf_N  + pdf_E  + pdf_S  + pdf_W  +
 303                                    pdf_NE + pdf_SE + pdf_SW + pdf_NW +
 304                                    pdf_T  + pdf_TN + pdf_TE + pdf_TS + pdf_TW +
 305                                    pdf_B  + pdf_BN + pdf_BE + pdf_BS + pdf_BW;
 306
 307                         dir_indep_trm = dens - (ux * ux + uy * uy + uz * uz)*F(3.0)/F(2.0);
 308
 309                         adjListIndex = index * N_D3Q19_IDX;
 310
 311                         // direction: w_0
 312                         src[I(index, D3Q19_C)             ]  = pdf_C - omegaEven*(pdf_C - w_0*dir_indep_trm);
 313
 314                         // direction: w_1
 315                         w_1_indep = w_1*dir_indep_trm;
 316
 317                         ui = uy;
 318                         evenPart = omegaEven*( F(0.5)*(pdf_N + pdf_S) - ui*ui*w_1_nine_half - w_1_indep );
 319                         oddPart = omegaOdd*(F(0.5)*(pdf_N - pdf_S) - ui*w_1_x3 );
 320                         src[adjList[adjListIndex + D3Q19_N]]  = pdf_N - evenPart - oddPart;
 321                         src[adjList[adjListIndex + D3Q19_S]]  = pdf_S - evenPart + oddPart;
 322
 323                         ui = ux;
 324                         evenPart = omegaEven*( F(0.5)*(pdf_E + pdf_W) - ui*ui*w_1_nine_half - w_1_indep );
 325                         oddPart = omegaOdd*(F(0.5)*(pdf_E - pdf_W) - ui*w_1_x3 );
 326                         src[adjList[adjListIndex + D3Q19_E]]  = pdf_E - evenPart - oddPart;
 327                         src[adjList[adjListIndex + D3Q19_W]]  = pdf_W - evenPart + oddPart;
 328
 329                         ui = uz;
 330                         evenPart = omegaEven*( F(0.5)*(pdf_T + pdf_B) - ui*ui*w_1_nine_half - w_1_indep );
 331                         oddPart = omegaOdd*(F(0.5)*(pdf_T - pdf_B) - ui*w_1_x3 );
 332                         src[adjList[adjListIndex + D3Q19_T]]  = pdf_T - evenPart - oddPart;
 333                         src[adjList[adjListIndex + D3Q19_B]]  = pdf_B - evenPart + oddPart;
 334
 335                         // direction: w_2
 336                         w_2_indep = w_2*dir_indep_trm;
 337
 338                         ui = -ux + uy;
 339                         evenPart = omegaEven*( F(0.5)*(pdf_NW + pdf_SE) - ui*ui*w_2_nine_half - w_2_indep );
 340                         oddPart = omegaOdd*(F(0.5)*(pdf_NW - pdf_SE) - ui*w_2_x3 );
 341                         src[adjList[adjListIndex + D3Q19_NW]] = pdf_NW - evenPart - oddPart;
 342                         src[adjList[adjListIndex + D3Q19_SE]] = pdf_SE - evenPart + oddPart;
 343
 344                         ui = ux + uy;
 345                         evenPart = omegaEven*( F(0.5)*(pdf_NE + pdf_SW) - ui*ui*w_2_nine_half - w_2_indep );
 346                         oddPart = omegaOdd*(F(0.5)*(pdf_NE - pdf_SW) - ui*w_2_x3 );
 347                         src[adjList[adjListIndex + D3Q19_NE]] = pdf_NE - evenPart - oddPart;
 348                         src[adjList[adjListIndex + D3Q19_SW]] = pdf_SW - evenPart + oddPart;
 349
 350                         ui = -ux + uz;
 351                         evenPart = omegaEven*( F(0.5)*(pdf_TW + pdf_BE) - ui*ui*w_2_nine_half - w_2_indep );
 352                         oddPart = omegaOdd*(F(0.5)*(pdf_TW - pdf_BE) - ui*w_2_x3 );
 353                         src[adjList[adjListIndex + D3Q19_TW]] = pdf_TW - evenPart - oddPart;
 354                         src[adjList[adjListIndex + D3Q19_BE]] = pdf_BE - evenPart + oddPart;
 355
 356                         ui = ux + uz;
 357                         evenPart = omegaEven*( F(0.5)*(pdf_TE + pdf_BW) - ui*ui*w_2_nine_half - w_2_indep );
 358                         oddPart = omegaOdd*(F(0.5)*(pdf_TE - pdf_BW) - ui*w_2_x3 );
 359                         src[adjList[adjListIndex + D3Q19_TE]] = pdf_TE - evenPart - oddPart;
 360                         src[adjList[adjListIndex + D3Q19_BW]] = pdf_BW - evenPart + oddPart;
 361
 362                         ui = -uy + uz;
 363                         evenPart = omegaEven*( F(0.5)*(pdf_TS + pdf_BN) - ui*ui*w_2_nine_half - w_2_indep );
 364                         oddPart = omegaOdd*(F(0.5)*(pdf_TS - pdf_BN) - ui*w_2_x3 );
 365                         src[adjList[adjListIndex + D3Q19_TS]] = pdf_TS - evenPart - oddPart;
 366                         src[adjList[adjListIndex + D3Q19_BN]] = pdf_BN - evenPart + oddPart;
 367
 368                         ui = uy + uz;
 369                         evenPart = omegaEven*( F(0.5)*(pdf_TN + pdf_BS) - ui*ui*w_2_nine_half - w_2_indep );
 370                         oddPart = omegaOdd*(F(0.5)*(pdf_TN - pdf_BS) - ui*w_2_x3 );
 371                         src[adjList[adjListIndex + D3Q19_TN]] = pdf_TN - evenPart - oddPart;
 372                         src[adjList[adjListIndex + D3Q19_BS]] = pdf_BS - evenPart + oddPart;
 373
 374                         #undef I
 375                 } // loop over fluid nodes
 376
 377                 X_LIKWID_STOP("list-aa-odd");
 378
 379                 // save current iteration
 380                 kdl->Iteration = iter + 1;
 381
 382                 #ifdef VERIFICATION
 383                         kd->PdfsActive = src;
 384                         KernelAddBodyForce(kd, ld, cd);
 385                 #endif
 386
 387                 #ifdef VTK_OUTPUT
 388                         if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
 389                                 kd->PdfsActive = src;
 390                                 VtkWrite(ld, kd, cd, iter);
 391                         }
 392                 #endif
 393
 394                 #ifdef STATISTICS
 395                         kd->PdfsActive = src;
 396                         KernelStatistics(kd, ld, cd, iter);
 397                 #endif
 398
 399
 400         } // for (int iter = 0; ...
 401
 402         X_KERNEL_END(kernelData);
 403
 404
 405 #ifdef VTK_OUTPUT
 406         if (cd->VtkOutput) {
 407                 kd->PdfsActive = src;
 408                 VtkWrite(ld, kd, cd, maxIterations);
 409         }
 410 #endif
 411
 412 #ifdef STATISTICS
 413         kd->PdfsActive = src;
 414         KernelStatistics(kd, ld, cd, maxIterations);
 415 #endif
 416
 417         return;
 418 }