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