1 // --------------------------------------------------------------------------
4 // Markus Wittmann, 2016-2017
5 // RRZE, University of Erlangen-Nuremberg, Germany
6 // markus.wittmann -at- fau.de or hpc -at- rrze.fau.de
9 // LSS, University of Erlangen-Nuremberg, Germany
11 // This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels).
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.
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.
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/>.
26 // --------------------------------------------------------------------------
27 #include "BenchKernelD3Q19AaVecCommon.h"
42 static void KernelEven(LatticeDesc * ld, KernelData * kd, CaseData * cd);
43 static void KernelOdd( LatticeDesc * ld, KernelData * kd, CaseData * cd);
46 void DumpPdfs(LatticeDesc * ld, KernelData * kd, int zStart, int zStop, int iter, const char * prefix)
50 int * lDims = ld->Dims;
51 int * gDims = kd->GlobalDims;
59 int localZStart = zStart;
60 int localZStop = zStop;
62 if (localZStart == -1) localZStart = 0;
63 if (localZStop == -1) localZStop = gDims[2] - 1;
65 printf("D iter: %d\n", iter);
67 for (int dir = 0; dir < 19; ++dir) {
68 for (int z = localZStop; z >= localZStart; --z) {
69 printf("D [%2d][%2d][%s] plane % 2d\n", iter, dir, prefix, z);
71 for(int y = 0; y < nY; ++y) {
72 printf("D [%2d][%2d][%s] %2d ", iter, dir, prefix, y);
74 for(int x = 0; x < nX; ++x) {
76 if (1) { // ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) {
78 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
79 pdfs[dir] = kd->PdfsActive[I(x, y, z, dir)];
81 // kd->GetNode(kd, x, y, z, pdfs);
87 printf("%.16e ", pdfs[dir]);
97 void FNAME(D3Q19AaVecKernel)(LatticeDesc * ld, KernelData * kd, CaseData * cd)
103 Assert(cd->Omega > F(0.0));
104 Assert(cd->Omega < F(2.0));
106 KernelDataAa * kda = KDA(kd);
108 PdfT * src = kd->PdfsActive;
110 int maxIterations = cd->MaxIterations;
114 kd->PdfsActive = src;
115 VtkWrite(ld, kd, cd, -1);
120 kd->PdfsActive = src;
121 KernelStatistics(kd, ld, cd, 0);
124 Assert((maxIterations % 2) == 0);
128 for (int iter = 0; iter < maxIterations; iter += 2) {
130 // --------------------------------------------------------------------
132 // --------------------------------------------------------------------
134 X_LIKWID_START("aa-vec-even");
138 KernelEven(ld, kd, cd);
141 X_LIKWID_STOP("aa-vec-even");
143 // Fixup bounce back PDFs.
145 #pragma omp parallel for default(none) \
148 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
149 src[kd->BounceBackPdfsSrc[i]] = src[kd->BounceBackPdfsDst[i]];
152 // save current iteration
153 kda->Iteration = iter;
156 kd->PdfsActive = src;
157 KernelAddBodyForce(kd, ld, cd);
161 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
162 kd->PdfsActive = src;
163 VtkWrite(ld, kd, cd, iter);
168 kd->PdfsActive = src;
169 KernelStatistics(kd, ld, cd, iter);
172 // --------------------------------------------------------------------
174 // --------------------------------------------------------------------
176 X_LIKWID_START("aa-vec-odd");
180 KernelOdd(ld, kd, cd);
183 // Stop counters before bounce back. Else computing loop balance will
186 X_LIKWID_STOP("aa-vec-odd");
188 // Fixup bounce back PDFs.
190 #pragma omp parallel for default(none) \
193 for (int i = 0; i < kd->nBounceBackPdfs; ++i) {
194 src[kd->BounceBackPdfsDst[i]] = src[kd->BounceBackPdfsSrc[i]];
197 // save current iteration
198 kda->Iteration = iter + 1;
201 kd->PdfsActive = src;
202 KernelAddBodyForce(kd, ld, cd);
206 if (cd->VtkOutput && ((iter + 1) % cd->VtkModulus) == 0) {
207 kd->PdfsActive = src;
208 VtkWrite(ld, kd, cd, iter + 1);
213 kd->PdfsActive = src;
214 KernelStatistics(kd, ld, cd, iter + 1);
218 } // for (int iter = 0; ...
225 kd->PdfsActive = src;
226 VtkWrite(ld, kd, cd, maxIterations);
234 static void KernelEven(LatticeDesc * ld, KernelData * kd, CaseData * cd) // {{{
240 Assert(cd->Omega > F(0.0));
241 Assert(cd->Omega < F(2.0));
243 KernelDataAa * kda = KDA(kd);
245 int nX = ld->Dims[0];
246 int nY = ld->Dims[1];
247 int nZ = ld->Dims[2];
249 int * gDims = kd->GlobalDims;
251 int oX = kd->Offsets[0];
252 int oY = kd->Offsets[1];
253 int oZ = kd->Offsets[2];
256 blk[0] = kda->Blk[0];
257 blk[1] = kda->Blk[1];
258 blk[2] = kda->Blk[2];
260 PdfT omega = cd->Omega;
261 PdfT omegaEven = omega;
263 PdfT magicParam = F(1.0) / F(12.0);
264 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
266 const PdfT w_0 = F(1.0) / F(3.0);
267 const PdfT w_1 = F(1.0) / F(18.0);
268 const PdfT w_2 = F(1.0) / F(36.0);
270 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
271 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
274 VPDFT VONE_HALF = VSET(F(0.5));
275 VPDFT VTHREE_HALF = VSET(F(3.0) / F(2.0));
277 VPDFT vw_1_indep, vw_2_indep;
278 VPDFT vw_0 = VSET(w_0);
279 VPDFT vw_1 = VSET(w_1);
280 VPDFT vw_2 = VSET(w_2);
282 VPDFT vw_1_x3 = VSET(w_1_x3);
283 VPDFT vw_2_x3 = VSET(w_2_x3);
284 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
285 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
287 VPDFT vui, vux, vuy, vuz, vdens;
289 VPDFT vevenPart, voddPart, vdir_indep_trm;
291 VPDFT vomegaEven = VSET(omegaEven);
292 VPDFT vomegaOdd = VSET(omegaOdd);
294 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
295 #define X(name, idx, idxinv, x, y, z) VPDFT JOIN(vpdf_,name);
299 PdfT * src = kd->Pdfs[0];
305 nThreads = omp_get_max_threads();
306 threadId = omp_get_thread_num();
309 // TODO: Currently only a 1-D decomposition is applied. For achritectures
310 // with a lot of cores we want at least 2-D.
312 int threadStartX = nX / nThreads * threadId;
313 int threadEndX = nX / nThreads * (threadId + 1);
315 if (nX % nThreads > 0) {
316 if (nX % nThreads > threadId) {
317 threadStartX += threadId;
318 threadEndX += threadId + 1;
321 threadStartX += nX % nThreads;
322 threadEndX += nX % nThreads;
326 AssertMsg((blk[2] % VSIZE == 0) || blk[2] >= nZ, "Blocking in z direction must be a multiple of VSIZE = %d or larger than z dimension.", VSIZE);
328 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
329 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
330 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
332 int eX = MIN(bX + blk[0], threadEndX + oX);
333 int eY = MIN(bY + blk[1], nY + oY);
334 int eZ = MIN(bZ + blk[2], nZ + oZ);
336 for (int x = bX; x < eX; x += 1) {
337 for (int y = bY; y < eY; y += 1) {
338 for (int z = bZ; z < eZ; z += VSIZE) { // LOOP aa-vec-even
340 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
342 // Load PDFs of local cell: pdf_N = src[I(x, y, z, D3Q19_N)]; ...
343 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(&src[I(x, y, z, idx)]);
348 vux = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_E,VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_TE,vpdf_BE)),vpdf_W),vpdf_NW),vpdf_SW),vpdf_TW),vpdf_BW);
349 vuy = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_N,VADD(vpdf_NE,vpdf_NW)),VADD(vpdf_TN,vpdf_BN)),vpdf_S),vpdf_SE),vpdf_SW),vpdf_TS),vpdf_BS);
350 vuz = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_T,VADD(vpdf_TE,vpdf_TW)),VADD(vpdf_TN,vpdf_TS)),vpdf_B),vpdf_BE),vpdf_BW),vpdf_BN),vpdf_BS);
352 vdens = VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(vpdf_C,VADD(vpdf_N,vpdf_E)),VADD(vpdf_S,vpdf_W)),VADD(vpdf_NE,vpdf_SE)),
353 VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),
354 VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
356 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
358 VSTU(&src[I(x, y, z, D3Q19_C)],VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
360 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
363 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_N,vpdf_S)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
364 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_N,vpdf_S)),VMUL(vui,vw_1_x3)));
365 VSTU(&src[I(x, y, z, D3Q19_S)],VSUB(VSUB(vpdf_N,vevenPart),voddPart));
366 VSTU(&src[I(x, y, z, D3Q19_N)],VADD(VSUB(vpdf_S,vevenPart),voddPart));
369 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_E,vpdf_W)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
370 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_E,vpdf_W)),VMUL(vui,vw_1_x3)));
371 VSTU(&src[I(x, y, z, D3Q19_W)],VSUB(VSUB(vpdf_E,vevenPart),voddPart));
372 VSTU(&src[I(x, y, z, D3Q19_E)],VADD(VSUB(vpdf_W,vevenPart),voddPart));
375 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_T,vpdf_B)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
376 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_T,vpdf_B)),VMUL(vui,vw_1_x3)));
377 VSTU(&src[I(x, y, z, D3Q19_B)],VSUB(VSUB(vpdf_T,vevenPart),voddPart));
378 VSTU(&src[I(x, y, z, D3Q19_T)],VADD(VSUB(vpdf_B,vevenPart),voddPart));
380 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
383 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NW,vpdf_SE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
384 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NW,vpdf_SE)),VMUL(vui,vw_2_x3)));
385 VSTU(&src[I(x, y, z, D3Q19_SE)],VSUB(VSUB(vpdf_NW,vevenPart),voddPart));
386 VSTU(&src[I(x, y, z, D3Q19_NW)],VADD(VSUB(vpdf_SE,vevenPart),voddPart));
389 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NE,vpdf_SW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
390 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NE,vpdf_SW)),VMUL(vui,vw_2_x3)));
391 VSTU(&src[I(x, y, z, D3Q19_SW)],VSUB(VSUB(vpdf_NE,vevenPart),voddPart));
392 VSTU(&src[I(x, y, z, D3Q19_NE)],VADD(VSUB(vpdf_SW,vevenPart),voddPart));
395 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TW,vpdf_BE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
396 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TW,vpdf_BE)),VMUL(vui,vw_2_x3)));
397 VSTU(&src[I(x, y, z, D3Q19_BE)],VSUB(VSUB(vpdf_TW,vevenPart),voddPart));
398 VSTU(&src[I(x, y, z, D3Q19_TW)],VADD(VSUB(vpdf_BE,vevenPart),voddPart));
401 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TE,vpdf_BW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
402 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TE,vpdf_BW)),VMUL(vui,vw_2_x3)));
403 VSTU(&src[I(x, y, z, D3Q19_BW)],VSUB(VSUB(vpdf_TE,vevenPart),voddPart));
404 VSTU(&src[I(x, y, z, D3Q19_TE)],VADD(VSUB(vpdf_BW,vevenPart),voddPart));
407 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TS,vpdf_BN)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
408 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TS,vpdf_BN)),VMUL(vui,vw_2_x3)));
409 VSTU(&src[I(x, y, z, D3Q19_BN)],VSUB(VSUB(vpdf_TS,vevenPart),voddPart));
410 VSTU(&src[I(x, y, z, D3Q19_TS)],VADD(VSUB(vpdf_BN,vevenPart),voddPart));
413 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TN,vpdf_BS)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
414 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TN,vpdf_BS)),VMUL(vui,vw_2_x3)));
415 VSTU(&src[I(x, y, z, D3Q19_BS)],VSUB(VSUB(vpdf_TN,vevenPart),voddPart));
416 VSTU(&src[I(x, y, z, D3Q19_TN)],VADD(VSUB(vpdf_BS,vevenPart),voddPart));
420 } } } // blocked x, y, z
428 static void KernelOdd(LatticeDesc * ld, KernelData * kd, CaseData * cd) // {{{
434 Assert(cd->Omega > F(0.0));
435 Assert(cd->Omega < F(2.0));
437 KernelDataAa * kda = KDA(kd);
439 int nX = ld->Dims[0];
440 int nY = ld->Dims[1];
441 int nZ = ld->Dims[2];
443 int * gDims = kd->GlobalDims;
445 int oX = kd->Offsets[0];
446 int oY = kd->Offsets[1];
447 int oZ = kd->Offsets[2];
450 blk[0] = kda->Blk[0];
451 blk[1] = kda->Blk[1];
452 blk[2] = kda->Blk[2];
454 PdfT omega = cd->Omega;
455 PdfT omegaEven = omega;
457 PdfT magicParam = F(1.0) / F(12.0);
458 PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
460 const PdfT w_0 = F(1.0) / F(3.0);
461 const PdfT w_1 = F(1.0) / F(18.0);
462 const PdfT w_2 = F(1.0) / F(36.0);
464 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
465 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
467 VPDFT VONE_HALF = VSET(F(0.5));
468 VPDFT VTHREE_HALF = VSET(F(3.0) / F(2.0));
470 VPDFT vw_1_indep, vw_2_indep;
471 VPDFT vw_0 = VSET(w_0);
472 VPDFT vw_1 = VSET(w_1);
473 VPDFT vw_2 = VSET(w_2);
475 VPDFT vw_1_x3 = VSET(w_1_x3);
476 VPDFT vw_2_x3 = VSET(w_2_x3);
477 VPDFT vw_1_nine_half = VSET(w_1_nine_half);
478 VPDFT vw_2_nine_half = VSET(w_2_nine_half);
480 VPDFT vui, vux, vuy, vuz, vdens;
482 VPDFT vevenPart, voddPart, vdir_indep_trm;
484 VPDFT vomegaEven = VSET(omegaEven);
485 VPDFT vomegaOdd = VSET(omegaOdd);
487 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
488 #define X(name, idx, idxinv, x, y, z) VPDFT JOIN(vpdf_,name);
492 PdfT * src = kd->Pdfs[0];
498 threadId = omp_get_thread_num();
499 nThreads = omp_get_max_threads();
502 // TODO: Currently only a 1-D decomposition is applied. For achritectures
503 // with a lot of cores we want at least 2-D.
504 int threadStartX = nX / nThreads * threadId;
505 int threadEndX = nX / nThreads * (threadId + 1);
507 if (nX % nThreads > 0) {
508 if (nX % nThreads > threadId) {
509 threadStartX += threadId;
510 threadEndX += threadId + 1;
513 threadStartX += nX % nThreads;
514 threadEndX += nX % nThreads;
518 AssertMsg((blk[2] % VSIZE == 0) || blk[2] >= nZ, "Blocking in z direction must be a multiple of VSIZE = %d or larger than z dimension.", VSIZE);
520 for (int bX = oX + threadStartX; bX < threadEndX + oX; bX += blk[0]) {
521 for (int bY = oY; bY < nY + oY; bY += blk[1]) {
522 for (int bZ = oZ; bZ < nZ + oZ; bZ += blk[2]) {
524 int eX = MIN(bX + blk[0], threadEndX + oX);
525 int eY = MIN(bY + blk[1], nY + oY);
526 int eZ = MIN(bZ + blk[2], nZ + oZ);
528 for (int x = bX; x < eX; ++x) {
529 for (int y = bY; y < eY; ++y) {
530 for (int z = bZ; z < eZ; z += VSIZE) { // LOOP aa-vec-odd
532 #define I(x, y, z, dir) P_INDEX_5(gDims, (x), (y), (z), (dir))
535 #define X(name, idx, idxinv, _x, _y, _z) JOIN(vpdf_,name) = VLDU(&src[I(x - _x, y - _y, z - _z, idxinv)]);
539 vux = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_E,VADD(vpdf_NE,vpdf_SE)),VADD(vpdf_TE,vpdf_BE)),vpdf_W),vpdf_NW),vpdf_SW),vpdf_TW),vpdf_BW);
540 vuy = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_N,VADD(vpdf_NE,vpdf_NW)),VADD(vpdf_TN,vpdf_BN)),vpdf_S),vpdf_SE),vpdf_SW),vpdf_TS),vpdf_BS);
541 vuz = VSUB(VSUB(VSUB(VSUB(VSUB(VADD(VADD(vpdf_T,VADD(vpdf_TE,vpdf_TW)),VADD(vpdf_TN,vpdf_TS)),vpdf_B),vpdf_BE),vpdf_BW),vpdf_BN),vpdf_BS);
543 vdens = VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(VADD(vpdf_C,VADD(vpdf_N,vpdf_E)),VADD(vpdf_S,vpdf_W)),VADD(vpdf_NE,vpdf_SE)),
544 VADD(vpdf_SW,vpdf_NW)),VADD(vpdf_T,vpdf_TN)),VADD(vpdf_TE,vpdf_TS)),VADD(vpdf_TW,vpdf_B)),VADD(vpdf_BN,vpdf_BE)),VADD(vpdf_BS,vpdf_BW));
546 vdir_indep_trm = VSUB(vdens,VMUL(VADD(VADD(VMUL(vux,vux),VMUL(vuy,vuy)),VMUL(vuz,vuz)),VTHREE_HALF));
548 VSTU(&src[I(x, y, z, D3Q19_C)],VSUB(vpdf_C,VMUL(vomegaEven,VSUB(vpdf_C,VMUL(vw_0,vdir_indep_trm)))));
550 vw_1_indep = VMUL(vw_1,vdir_indep_trm);
553 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_N,vpdf_S)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
554 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_N,vpdf_S)),VMUL(vui,vw_1_x3)));
555 VSTU(&src[I(x, y + 1, z, D3Q19_N)], VSUB(VSUB(vpdf_N,vevenPart),voddPart));
556 VSTU(&src[I(x, y - 1, z, D3Q19_S)], VADD(VSUB(vpdf_S,vevenPart),voddPart));
559 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_E,vpdf_W)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
560 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_E,vpdf_W)),VMUL(vui,vw_1_x3)));
561 VSTU(&src[I(x + 1, y, z, D3Q19_E)], VSUB(VSUB(vpdf_E,vevenPart),voddPart));
562 VSTU(&src[I(x - 1, y, z, D3Q19_W)], VADD(VSUB(vpdf_W,vevenPart),voddPart));
565 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_T,vpdf_B)),VMUL(vui,VMUL(vui,vw_1_nine_half))),vw_1_indep));
566 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_T,vpdf_B)),VMUL(vui,vw_1_x3)));
567 VSTU(&src[I(x, y, z + 1, D3Q19_T)], VSUB(VSUB(vpdf_T,vevenPart),voddPart));
568 VSTU(&src[I(x, y, z - 1, D3Q19_B)], VADD(VSUB(vpdf_B,vevenPart),voddPart));
570 vw_2_indep = VMUL(vw_2,vdir_indep_trm);
573 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NW,vpdf_SE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
574 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NW,vpdf_SE)),VMUL(vui,vw_2_x3)));
575 VSTU(&src[I(x - 1, y + 1, z, D3Q19_NW)], VSUB(VSUB(vpdf_NW,vevenPart),voddPart));
576 VSTU(&src[I(x + 1, y - 1, z, D3Q19_SE)], VADD(VSUB(vpdf_SE,vevenPart),voddPart));
579 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_NE,vpdf_SW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
580 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_NE,vpdf_SW)),VMUL(vui,vw_2_x3)));
581 VSTU(&src[I(x + 1, y + 1, z, D3Q19_NE)], VSUB(VSUB(vpdf_NE,vevenPart),voddPart));
582 VSTU(&src[I(x - 1, y - 1, z, D3Q19_SW)], VADD(VSUB(vpdf_SW,vevenPart),voddPart));
585 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TW,vpdf_BE)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
586 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TW,vpdf_BE)),VMUL(vui,vw_2_x3)));
587 VSTU(&src[I(x - 1, y, z + 1, D3Q19_TW)], VSUB(VSUB(vpdf_TW,vevenPart),voddPart));
588 VSTU(&src[I(x + 1, y, z - 1, D3Q19_BE)], VADD(VSUB(vpdf_BE,vevenPart),voddPart));
591 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TE,vpdf_BW)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
592 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TE,vpdf_BW)),VMUL(vui,vw_2_x3)));
593 VSTU(&src[I(x + 1, y, z + 1, D3Q19_TE)], VSUB(VSUB(vpdf_TE,vevenPart),voddPart));
594 VSTU(&src[I(x - 1, y, z - 1, D3Q19_BW)], VADD(VSUB(vpdf_BW,vevenPart),voddPart));
597 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TS,vpdf_BN)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
598 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TS,vpdf_BN)),VMUL(vui,vw_2_x3)));
599 VSTU(&src[I(x, y - 1, z + 1, D3Q19_TS)], VSUB(VSUB(vpdf_TS,vevenPart),voddPart));
600 VSTU(&src[I(x, y + 1, z - 1, D3Q19_BN)], VADD(VSUB(vpdf_BN,vevenPart),voddPart));
603 vevenPart = VMUL(vomegaEven,VSUB(VSUB(VMUL(VONE_HALF,VADD(vpdf_TN,vpdf_BS)),VMUL(vui,VMUL(vui,vw_2_nine_half))),vw_2_indep));
604 voddPart = VMUL(vomegaOdd,VSUB(VMUL(VONE_HALF,VSUB(vpdf_TN,vpdf_BS)),VMUL(vui,vw_2_x3)));
605 VSTU(&src[I(x, y + 1, z + 1, D3Q19_TN)], VSUB(VSUB(vpdf_TN,vevenPart),voddPart));
606 VSTU(&src[I(x, y - 1, z - 1, D3Q19_BS)], VADD(VSUB(vpdf_BS,vevenPart),voddPart));
610 } } } // blocked x, y, z