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