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 "BenchKernelD3Q19ListPullSplitNtCommon.h"
49 #define TMP_INDEX(tmp_index, tmp_dir) nTmpArray * (tmp_dir) + (tmp_index)
51 void FNAME(KernelPullSplitNt1S)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
55 Assert(kernelData != NULL);
58 Assert(cd->Omega > F(0.0));
59 Assert(cd->Omega < F(2.0));
61 KernelData * kd = (KernelData *)kernelData;
62 KernelDataList * kdl = KDL(kernelData);
63 KernelDataListRia * kdlr = KDLR(kernelData);
65 PdfT omega = cd->Omega;
66 const PdfT omegaEven = omega;
68 PdfT magicParam = F(1.0) / F(12.0);
69 const PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
72 const PdfT w_0 = F(1.0) / F( 3.0);
73 const PdfT w_1 = F(1.0) / F(18.0);
74 const PdfT w_2 = F(1.0) / F(36.0);
76 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
77 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
79 const VPDFT vw_1_x3 = VSET(w_1_x3);
80 const VPDFT vw_2_x3 = VSET(w_2_x3);
82 const VPDFT vw_1_nine_half = VSET(w_1_nine_half);
83 const VPDFT vw_2_nine_half = VSET(w_2_nine_half);
85 const VPDFT vomegaEven = VSET(omegaEven);
86 const VPDFT vomegaOdd = VSET(omegaOdd);
88 const VPDFT voneHalf = VSET(F(0.5));
90 // uint32_t nConsecNodes = kdlr->nConsecNodes;
91 // uint32_t * consecNodes = kdlr->ConsecNodes;
92 // uint32_t consecIndex = 0;
93 // uint32_t consecValue = 0;
95 PdfT * src = kd->Pdfs[0];
96 PdfT * dst = kd->Pdfs[1];
99 int maxIterations = cd->MaxIterations;
101 int nFluid = kdl->nFluid;
102 int nCells = kdl->nCells;
104 int nTmpArray = kdlr->nTmpArray;
106 Assert(nTmpArray % VSIZE == 0);
108 uint32_t * adjList = kdl->AdjList;
112 kd->PdfsActive = src;
113 VtkWrite(ld, kd, cd, -1);
118 kd->PdfsActive = src;
119 KernelStatistics(kd, ld, cd, 0);
122 X_KERNEL_START(kernelData);
124 X_LIKWID_START("list-pull-split-nt-1s");
127 #pragma omp parallel default(none) \
128 shared(nFluid, nCells, kd, kdl, adjList, src, dst, \
129 cd, maxIterations, ld, tmp, nTmpArray, \
133 uint32_t adjListIndex;
136 VPDFT vux, vuy, vuz, vui;
138 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
141 VPDFT vpdf_a, vpdf_b;
143 PdfT evenPart, oddPart, dir_indep_trm, dens;
144 PdfT w_1_indep, w_2_indep;
145 VPDFT vevenPart, voddPart;
146 VPDFT vw_1_indep, vw_2_indep;
151 MemAllocAligned((void **)&tmpArray, sizeof(PdfT) * nTmpArray * N_TMP, VSIZE * sizeof(PdfT));
157 nThreads = omp_get_max_threads();
158 threadId = omp_get_thread_num();
161 int nCellsThread = nFluid / nThreads;
162 int blIndexStart = threadId * nCellsThread;
164 if (threadId < nFluid % nThreads) {
165 blIndexStart += threadId;
169 blIndexStart += nFluid % nThreads;
172 int blIndexStop = blIndexStart + nCellsThread;
174 // We have three loops:
175 // 1. Peeling to ensure alignment for non-temporal stores in loop 2 is correct.
176 // 2. Vectorized handling of nodes.
177 // 3. Remaining nodes, less than vector size.
179 unsigned long addrStart = (unsigned long)&(src[P_INDEX_3(nCells, blIndexStart, 0)]);
180 int nCellsUnaligned = (VSIZE - (int)((addrStart / sizeof(PdfT)) % VSIZE)) % VSIZE;
182 int nCellsVectorized = nCellsThread - nCellsUnaligned;
183 nCellsVectorized = nCellsVectorized - (nCellsVectorized % VSIZE);
185 int blIndexVec = blIndexStart + nCellsUnaligned;
186 int blIndexRemaining = blIndexStart + nCellsUnaligned + nCellsVectorized;
188 // printf("%d [%d, %d, %d, %d[\n", threadId, blIndexStart, blIndexVec, blIndexRemaining, blIndexStop);
190 for(int iter = 0; iter < maxIterations; ++iter) {
194 #define INDEX_START blIndexStart
195 #define INDEX_STOP blIndexVec
196 #include "BenchKernelD3Q19ListPullSplitNt1SScalar.h"
198 #define INDEX_START blIndexVec
199 #define INDEX_STOP blIndexRemaining
200 #include "BenchKernelD3Q19ListPullSplitNt1SIntrinsics.h"
202 #define INDEX_START blIndexRemaining
203 #define INDEX_STOP blIndexStop
204 #include "BenchKernelD3Q19ListPullSplitNt1SScalar.h"
206 #define INDEX_START blIndexStart
207 #define INDEX_STOP blIndexStop
208 #include "BenchKernelD3Q19ListPullSplitNt1SScalar.h"
217 kd->PdfsActive = dst;
218 KernelAddBodyForce(kd, ld, cd);
222 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
223 kd->PdfsActive = dst;
224 VtkWrite(ld, kd, cd, iter);
229 kd->PdfsActive = dst;
230 KernelStatistics(kd, ld, cd, iter);
241 } // for (int iter = 0; ...
243 MemFree((void **)&tmpArray);
247 X_LIKWID_STOP("list-pull-split-nt-1s");
249 X_KERNEL_END(kernelData);
253 kd->PdfsActive = src;
254 VtkWrite(ld, kd, cd, maxIterations);
259 kd->PdfsActive = src;
260 KernelStatistics(kd, ld, cd, maxIterations);
266 void FNAME(KernelPullSplitNt2S)(LatticeDesc * ld, KernelData * kernelData, CaseData * cd)
270 Assert(kernelData != NULL);
273 Assert(cd->Omega > F(0.0));
274 Assert(cd->Omega < F(2.0));
276 KernelData * kd = (KernelData *)kernelData;
277 KernelDataList * kdl = KDL(kernelData);
278 KernelDataListRia * kdlr = KDLR(kernelData);
280 PdfT omega = cd->Omega;
281 const PdfT omegaEven = omega;
283 PdfT magicParam = F(1.0) / F(12.0);
284 const PdfT omegaOdd = F(1.0) / (F(0.5) + magicParam / (F(1.0) / omega - F(0.5)));
286 const PdfT w_0 = F(1.0) / F( 3.0);
287 const PdfT w_1 = F(1.0) / F(18.0);
288 const PdfT w_2 = F(1.0) / F(36.0);
290 const PdfT w_1_x3 = w_1 * F(3.0); const PdfT w_1_nine_half = w_1 * F(9.0) / F(2.0);
291 const PdfT w_2_x3 = w_2 * F(3.0); const PdfT w_2_nine_half = w_2 * F(9.0) / F(2.0);
293 const VPDFT vw_1_x3 = VSET(w_1_x3);
294 const VPDFT vw_2_x3 = VSET(w_2_x3);
296 const VPDFT vw_1_nine_half = VSET(w_1_nine_half);
297 const VPDFT vw_2_nine_half = VSET(w_2_nine_half);
299 const VPDFT vomegaEven = VSET(omegaEven);
300 const VPDFT vomegaOdd = VSET(omegaOdd);
302 const VPDFT voneHalf = VSET(F(0.5));
304 // uint32_t nConsecNodes = kdlr->nConsecNodes;
305 // uint32_t * consecNodes = kdlr->ConsecNodes;
306 // uint32_t consecIndex = 0;
307 // uint32_t consecValue = 0;
309 PdfT * src = kd->Pdfs[0];
310 PdfT * dst = kd->Pdfs[1];
313 int maxIterations = cd->MaxIterations;
315 int nFluid = kdl->nFluid;
316 int nCells = kdl->nCells;
318 int nTmpArray = kdlr->nTmpArray;
320 Assert(nTmpArray % VSIZE == 0);
322 uint32_t * adjList = kdl->AdjList;
326 kd->PdfsActive = src;
327 VtkWrite(ld, kd, cd, -1);
332 kd->PdfsActive = src;
333 KernelStatistics(kd, ld, cd, 0);
337 X_KERNEL_START(kernelData);
339 X_LIKWID_START("list-pull-split-nt-2s");
343 #pragma omp parallel default(none) \
344 shared(nFluid, nCells, kd, kdl, adjList, src, dst, \
345 cd, maxIterations, ld, tmp, nTmpArray, \
349 uint32_t adjListIndex;
352 VPDFT vux, vuy, vuz, vui;
354 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
357 VPDFT vpdf_a, vpdf_b;
359 PdfT evenPart, oddPart, dir_indep_trm, dens;
360 PdfT w_1_indep, w_2_indep;
361 VPDFT vevenPart, voddPart;
362 VPDFT vw_1_indep, vw_2_indep;
367 MemAlloc((void **)&tmpArray, sizeof(PdfT) * nTmpArray * N_TMP);
373 nThreads = omp_get_max_threads();
374 threadId = omp_get_thread_num();
377 int nCellsThread = nFluid / nThreads;
378 int blIndexStart = threadId * nCellsThread;
380 if (threadId < nFluid % nThreads) {
381 blIndexStart += threadId;
385 blIndexStart += nFluid % nThreads;
388 int blIndexStop = blIndexStart + nCellsThread;
390 // We have three loops:
391 // 1. Peeling to ensure alignment for non-temporal stores in loop 2 is correct.
392 // 2. Vectorized handling of nodes.
393 // 3. Remaining nodes, less than vector size.
395 unsigned long addrStart = (unsigned long)&(src[P_INDEX_3(nCells, blIndexStart, 0)]);
396 int nCellsUnaligned = (VSIZE - (int)((addrStart / sizeof(PdfT)) % VSIZE)) % VSIZE;
398 int nCellsVectorized = nCellsThread - nCellsUnaligned;
399 nCellsVectorized = nCellsVectorized - (nCellsVectorized % VSIZE);
401 int blIndexVec = blIndexStart + nCellsUnaligned;
402 int blIndexRemaining = blIndexStart + nCellsUnaligned + nCellsVectorized;
404 // printf("%d [%d, %d, %d, %d[\n", threadId, blIndexStart, blIndexVec, blIndexRemaining, blIndexStop);
406 for(int iter = 0; iter < maxIterations; ++iter) {
409 #define INDEX_START blIndexStart
410 #define INDEX_STOP blIndexVec
411 #include "BenchKernelD3Q19ListPullSplitNt2SScalar.h"
413 #define INDEX_START blIndexVec
414 #define INDEX_STOP blIndexRemaining
415 #include "BenchKernelD3Q19ListPullSplitNt2SIntrinsics.h"
417 #define INDEX_START blIndexRemaining
418 #define INDEX_STOP blIndexStop
419 #include "BenchKernelD3Q19ListPullSplitNt2SScalar.h"
421 #define INDEX_START blIndexStart
422 #define INDEX_STOP blIndexStop
423 #include "BenchKernelD3Q19ListPullSplitNt2SScalar.h"
431 kd->PdfsActive = dst;
432 KernelAddBodyForce(kd, ld, cd);
436 if (cd->VtkOutput && (iter % cd->VtkModulus) == 0) {
437 kd->PdfsActive = dst;
438 VtkWrite(ld, kd, cd, iter);
443 kd->PdfsActive = dst;
444 KernelStatistics(kd, ld, cd, iter);
455 } // for (int iter = 0; ...
457 MemFree((void **)&tmpArray);
460 X_LIKWID_STOP("list-pull-split-nt-2s");
462 X_KERNEL_END(kernelData);
466 kd->PdfsActive = src;
467 VtkWrite(ld, kd, cd, maxIterations);
472 kd->PdfsActive = src;
473 KernelStatistics(kd, ld, cd, maxIterations);