X-Git-Url: http://git.rrze.uni-erlangen.de/gitweb/?p=LbmBenchmarkKernelsPublic.git;a=blobdiff_plain;f=src%2FBenchKernelD3Q19Common.c;fp=src%2FBenchKernelD3Q19Common.c;h=e697bc8a573a2d5fb646a0686ed2157805c829f8;hp=0000000000000000000000000000000000000000;hb=109880839321408644c94a34eb31208460b9f46d;hpb=42cf91486fb5c1ad178b3d21935a1be563e5fa39 diff --git a/src/BenchKernelD3Q19Common.c b/src/BenchKernelD3Q19Common.c new file mode 100644 index 0000000..e697bc8 --- /dev/null +++ b/src/BenchKernelD3Q19Common.c @@ -0,0 +1,655 @@ +// -------------------------------------------------------------------------- +// +// Copyright +// Markus Wittmann, 2016-2017 +// RRZE, University of Erlangen-Nuremberg, Germany +// markus.wittmann -at- fau.de or hpc -at- rrze.fau.de +// +// Viktor Haag, 2016 +// LSS, University of Erlangen-Nuremberg, Germany +// +// This file is part of the Lattice Boltzmann Benchmark Kernels (LbmBenchKernels). +// +// LbmBenchKernels is free software: you can redistribute it and/or modify +// it under the terms of the GNU General Public License as published by +// the Free Software Foundation, either version 3 of the License, or +// (at your option) any later version. +// +// LbmBenchKernels is distributed in the hope that it will be useful, +// but WITHOUT ANY WARRANTY; without even the implied warranty of +// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +// GNU General Public License for more details. +// +// You should have received a copy of the GNU General Public License +// along with LbmBenchKernels. If not, see . +// +// -------------------------------------------------------------------------- +#include "BenchKernelD3Q19Common.h" + +#include "Memory.h" +#include "Vtk.h" + +#include +#include + + +// Forward definition. +void FNAME(D3Q19Kernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd); + +void FNAME(D3Q19BlkKernel)(LatticeDesc * ld, struct KernelData_ * kd, CaseData * cd); + + + +static void FNAME(BcGetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT * pdf) +{ + Assert(kd != NULL); + Assert(kd->PdfsActive != NULL); + Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); + Assert(pdf != NULL); + + Assert(x >= 0); + Assert(y >= 0); + Assert(z >= 0); + Assert(x < kd->Dims[0]); + Assert(y < kd->Dims[1]); + Assert(z < kd->Dims[2]); + Assert(dir >= 0); + Assert(dir < N_D3Q19); + + int oX = kd->Offsets[0]; + int oY = kd->Offsets[1]; + int oZ = kd->Offsets[2]; + +#ifdef PROP_MODEL_PUSH + int nx = x; + int ny = y; + int nz = z; +#elif PROP_MODEL_PULL + int nx = x - D3Q19_X[dir]; + int ny = y - D3Q19_Y[dir]; + int nz = z - D3Q19_Z[dir]; +#endif + + #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) + *pdf = kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, dir)]; + #undef I + + return; +} + +static void FNAME(BcSetPdf)(KernelData * kd, int x, int y, int z, int dir, PdfT pdf) +{ + Assert(kd != NULL); + Assert(kd->PdfsActive != NULL); + Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); + Assert(x >= 0); + Assert(y >= 0); + Assert(z >= 0); + Assert(x < kd->Dims[0]); + Assert(y < kd->Dims[1]); + Assert(z < kd->Dims[2]); + Assert(dir >= 0); + Assert(dir < N_D3Q19); + + int oX = kd->Offsets[0]; + int oY = kd->Offsets[1]; + int oZ = kd->Offsets[2]; + +#ifdef PROP_MODEL_PUSH + int nx = x; + int ny = y; + int nz = z; +#elif PROP_MODEL_PULL + int nx = x - D3Q19_X[dir]; + int ny = y - D3Q19_Y[dir]; + int nz = z - D3Q19_Z[dir]; +#endif + + #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) + kd->PdfsActive[I(nx + oX, ny + oY, nz + oZ, dir)] = pdf; + #undef I + + + return; +} + + +static void FNAME(GetNode)(KernelData * kd, int x, int y, int z, PdfT * pdfs) +{ + Assert(kd != NULL); + Assert(kd->PdfsActive != NULL); + Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); + Assert(pdfs != NULL); + Assert(x >= 0); + Assert(y >= 0); + Assert(z >= 0); + Assert(x < kd->Dims[0]); + Assert(y < kd->Dims[1]); + Assert(z < kd->Dims[2]); + + int oX = kd->Offsets[0]; + int oY = kd->Offsets[1]; + int oZ = kd->Offsets[2]; + + + #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) +#ifdef PROP_MODEL_PUSH + #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x + oX, y + oY, z + oZ, idx)]; +#elif PROP_MODEL_PULL + #define X(name, idx, idxinv, _x, _y, _z) pdfs[idx] = kd->PdfsActive[I(x + oX - (_x), y + oY - (_y), z + oZ - (_z), idx)]; +#endif + D3Q19_LIST + #undef X + #undef I + +#if 0 // DETECT NANs + + for (int d = 0; d < 19; ++d) { + if (isnan(pdfs[d])) { + printf("%d %d %d %d nan! get node\n", x, y, z, d); + + for (int d2 = 0; d2 < 19; ++d2) { + printf("%d: %e\n", d2, pdfs[d2]); + } + + exit(1); + } + } + +#endif + + return; +} + + +static void FNAME(SetNode)(KernelData * kd, int x, int y, int z, PdfT * pdfs) +{ + Assert(kd != NULL); + Assert(kd->PdfsActive != NULL); + Assert(kd->PdfsActive == kd->Pdfs[0] || kd->PdfsActive == kd->Pdfs[1]); + Assert(pdfs != NULL); + + Assert(x >= 0); + Assert(y >= 0); + Assert(z >= 0); + Assert(x < kd->Dims[0]); + Assert(y < kd->Dims[1]); + Assert(z < kd->Dims[2]); + + int oX = kd->Offsets[0]; + int oY = kd->Offsets[1]; + int oZ = kd->Offsets[2]; + + #define I(x, y, z, dir) P_INDEX_5(kd->GlobalDims, (x), (y), (z), (dir)) +#ifdef PROP_MODEL_PUSH + #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x + oX, y + oY, z + oZ, idx)] = pdfs[idx]; +#elif PROP_MODEL_PULL + #define X(name, idx, idxinv, _x, _y, _z) kd->PdfsActive[I(x + oX - (_x), y + oY - (_y), z + oZ - (_z), idx)] = pdfs[idx]; +#endif + D3Q19_LIST + #undef X + #undef I + + return; +} + + +static void ParameterUsage() +{ + printf("Kernel parameters:\n"); + printf(" [-blk ] [-blk-[xyz] ]\n"); + + return; +} + +static void ParseParameters(Parameters * params, int * blk) +{ + Assert(blk != NULL); + + blk[0] = 0; blk[1] = 0; blk[2] = 0; + + #define ARG_IS(param) (!strcmp(params->KernelArgs[i], param)) + #define NEXT_ARG_PRESENT() \ + do { \ + if (i + 1 >= params->nKernelArgs) { \ + printf("ERROR: argument %s requires a parameter.\n", params->KernelArgs[i]); \ + exit(1); \ + } \ + } while (0) + + + for (int i = 0; i < params->nKernelArgs; ++i) { + if (ARG_IS("-blk") || ARG_IS("--blk")) { + NEXT_ARG_PRESENT(); + + int tmp = strtol(params->KernelArgs[++i], NULL, 0); + + if (tmp <= 0) { + printf("ERROR: blocking parameter must be > 0.\n"); + exit(1); + } + + blk[0] = blk[1] = blk[2] = tmp; + } + else if (ARG_IS("-blk-x") || ARG_IS("--blk-x")) { + NEXT_ARG_PRESENT(); + + int tmp = strtol(params->KernelArgs[++i], NULL, 0); + + if (tmp <= 0) { + printf("ERROR: blocking parameter must be > 0.\n"); + exit(1); + } + + blk[0] = tmp; + } + else if (ARG_IS("-blk-y") || ARG_IS("--blk-y")) { + NEXT_ARG_PRESENT(); + + int tmp = strtol(params->KernelArgs[++i], NULL, 0); + + if (tmp <= 0) { + printf("ERROR: blocking parameter must be > 0.\n"); + exit(1); + } + + blk[1] = tmp; + } + else if (ARG_IS("-blk-z") || ARG_IS("--blk-z")) { + NEXT_ARG_PRESENT(); + + int tmp = strtol(params->KernelArgs[++i], NULL, 0); + + if (tmp <= 0) { + printf("ERROR: blocking parameter must be > 0.\n"); + exit(1); + } + + blk[2] = tmp; + } + else if (ARG_IS("-h") || ARG_IS("-help") || ARG_IS("--help")) { + ParameterUsage(); + exit(1); + } + else { + printf("ERROR: unknown kernel parameter.\n"); + ParameterUsage(); + exit(1); + } + } + + #undef ARG_IS + #undef NEXT_ARG_PRESENT + + return; +} + + +void FNAME(D3Q19BlkInit)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params) +{ + KernelDataEx * kdex = NULL; + MemAlloc((void **)&kdex, sizeof(KernelDataEx)); + + kdex->Blk[0] = 0; kdex->Blk[1] = 0; kdex->Blk[2] = 0; + + KernelData * kd = &kdex->kd; + *kernelData = kd; + + kd->nObstIndices = ld->nObst; + + // Ajust the dimensions according to padding, if used. + kd->Dims[0] = ld->Dims[0]; + kd->Dims[1] = ld->Dims[1]; + kd->Dims[2] = ld->Dims[2]; + + + int * lDims = ld->Dims; + int * gDims = kd->GlobalDims; + + gDims[0] = lDims[0] + 2; + gDims[1] = lDims[1] + 2; + gDims[2] = lDims[2] + 2; + + kd->Offsets[0] = 1; + kd->Offsets[1] = 1; + kd->Offsets[2] = 1; + + int lX = lDims[0]; + int lY = lDims[1]; + int lZ = lDims[2]; + + int gX = gDims[0]; + int gY = gDims[1]; + int gZ = gDims[2]; + + int oX = kd->Offsets[0]; + int oY = kd->Offsets[1]; + int oZ = kd->Offsets[2]; + + int blk[3] = { 0 }; + + int nCells = gX * gY * gZ; + + PdfT * pdfs[2]; + + ParseParameters(params, blk); + + if (blk[0] == 0) blk[0] = gX; + if (blk[1] == 0) blk[1] = gY; + if (blk[2] == 0) blk[2] = gZ; + + printf("# blocking x: %3d y: %3d z: %3d\n", blk[0], blk[1], blk[2]); + + + kdex->Blk[0] = blk[0]; kdex->Blk[1] = blk[1]; kdex->Blk[2] = blk[2]; + + + printf("# allocating data for %d LB nodes with padding (%lu bytes = %f MiB for both lattices)\n", + nCells, 2 * sizeof(PdfT) * nCells * N_D3Q19, + 2 * sizeof(PdfT) * nCells * N_D3Q19 / 1024.0 / 1024.0); + + MemAlloc((void **)&pdfs[0], sizeof(PdfT) * nCells * N_D3Q19); + MemAlloc((void **)&pdfs[1], sizeof(PdfT) * nCells * N_D3Q19); + + kd->Pdfs[0] = pdfs[0]; + kd->Pdfs[1] = pdfs[1]; + + // Initialize PDFs with some (arbitrary) data for correct NUMA placement. + // This depends on the chosen data layout. + // The structure of the loop should resemble the same "execution layout" + // as in the kernel! +#ifdef _OPENMP + #pragma omp parallel for collapse(3) +#endif + + for (int bZ = 0; bZ < gZ; bZ += blk[2]) { + for (int bY = 0; bY < gY; bY += blk[1]) { + for (int bX = 0; bX < gX; bX += blk[0]) { + + // Must do everything here, else it would break collapse. + int eZ = MIN(bZ + blk[2], gZ); + int eY = MIN(bY + blk[1], gY); + int eX = MIN(bX + blk[0], gX); + + for (int z = bZ; z < eZ; ++z) { + for (int y = bY; y < eY; ++y) { + for (int x = bX; x < eX; ++x) { + + for (int d = 0; d < N_D3Q19; ++d) { + pdfs[0][P_INDEX_5(gDims, x, y, z, d)] = 1.0; + pdfs[1][P_INDEX_5(gDims, x, y, z, d)] = 1.0; + } + + } + } + } + } + } + } + + // Initialize all PDFs to some standard value. + for (int z = 0; z < gZ; ++z) { + for (int y = 0; y < gY; ++y) { + for (int x = 0; x < gX; ++x) { + for (int d = 0; d < N_D3Q19; ++d) { + pdfs[0][P_INDEX_5(gDims, x, y, z, d)] = 0.0; + pdfs[1][P_INDEX_5(gDims, x, y, z, d)] = 0.0; + } + } + } + } + + + // Count how many *PDFs* need bounce back treatment. + + uint64_t nPdfs = ((uint64_t)19) * gX * gY * gZ; + + if (nPdfs > ((2LU << 31) - 1)) { + printf("ERROR: number of PDFs exceed 2^31.\n"); + exit(1); + } + + // Compiler bug? Incorrect computation of nBounceBackPdfs when using icc 15.0.2. + // Works when declaring nBounceBackPdfs as int64_t or using volatile. + volatile int nBounceBackPdfs = 0; + // int64_t nBounceBackPdfs = 0; + int nx, ny, nz, px, py, pz; + + // TODO: apply blocking? + + for (int z = 0; z < lZ; ++z) { + for (int y = 0; y < lY; ++y) { + for (int x = 0; x < lX; ++x) { + + if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) { + for (int d = 0; d < N_D3Q19; ++d) { +#ifdef PROP_MODEL_PUSH + nx = x + D3Q19_X[d]; + ny = y + D3Q19_Y[d]; + nz = z + D3Q19_Z[d]; +#elif PROP_MODEL_PULL + nx = x - D3Q19_X[d]; + ny = y - D3Q19_Y[d]; + nz = z - D3Q19_Z[d]; +#else + #error PROP_MODEL_NAME unknown. +#endif + // Check if neighbor is inside the lattice. + // if(nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { + // continue; + // } + if ((nx < 0 || nx >= lX) && ld->PeriodicX) { + ++nBounceBackPdfs; // Compiler bug --> see above + } + else if ((ny < 0 || ny >= lY) && ld->PeriodicY) { + ++nBounceBackPdfs; // Compiler bug --> see above + } + else if ((nz < 0 || nz >= lZ) && ld->PeriodicZ) { + ++nBounceBackPdfs; // Compiler bug --> see above + } + else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { + continue; + } + else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) { + ++nBounceBackPdfs; // Compiler bug --> see above + } + } + } + } + } + } + + + printf("# allocating %d indices for bounce back pdfs (%s for source and destination array)\n", nBounceBackPdfs, ByteToHuman(sizeof(int) * nBounceBackPdfs * 2)); + + MemAlloc((void **) & (kd->BounceBackPdfsSrc), sizeof(int) * nBounceBackPdfs + 100); + MemAlloc((void **) & (kd->BounceBackPdfsDst), sizeof(int) * nBounceBackPdfs + 100); + + kd->nBounceBackPdfs = nBounceBackPdfs; + nBounceBackPdfs = 0; + + int srcIndex; + int dstIndex; + + for (int z = 0; z < lZ; ++z) { + for (int y = 0; y < lY; ++y) { + for (int x = 0; x < lX; ++x) { + + if (ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)] != LAT_CELL_OBSTACLE) { + for (int d = 0; d < N_D3Q19; ++d) { +#ifdef PROP_MODEL_PUSH + nx = x + D3Q19_X[d]; + ny = y + D3Q19_Y[d]; + nz = z + D3Q19_Z[d]; +#elif PROP_MODEL_PULL + nx = x - D3Q19_X[d]; + ny = y - D3Q19_Y[d]; + nz = z - D3Q19_Z[d]; +#else + #error PROP_MODEL_NAME unknown. +#endif + + if ( ((nx < 0 || nx >= lX) && ld->PeriodicX) || + ((ny < 0 || ny >= lY) && ld->PeriodicY) || + ((nz < 0 || nz >= lZ) && ld->PeriodicZ) + ){ + // Implement periodic boundary in X direction. + + // If the target node reached through propagation is outside the lattice + // the kernel stores it in some buffer around the domain. + // From this position the PDF must be transported to the other side of the + // geometry. + + // Take PDF from outside the domain. + + // x periodic + if (nx < 0) { + px = lX - 1; + } + else if (nx >= lX) { + px = 0; + } else { + px = nx; + } + + // y periodic + if (ny < 0) { + py = lY - 1; + } + else if (ny >= lY) { + py = 0; + } else { + py = ny; + } + + // z periodic + if (nz < 0) { + pz = lZ - 1; + } + else if (nz >= lZ) { + pz = 0; + } else { + pz = nz; + } + + if (ld->Lattice[L_INDEX_4(lDims, px, py, pz)] == LAT_CELL_OBSTACLE) { +#ifdef PROP_MODEL_PUSH + srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); + dstIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, D3Q19_INV[d]); +#elif PROP_MODEL_PULL + srcIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, D3Q19_INV[d]); + dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); +#endif + } + else { + +#ifdef PROP_MODEL_PUSH + srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); + // Put it on the other side back into the domain. + dstIndex = P_INDEX_5(gDims, px + oX, py + oY, pz + oZ, d); +#elif PROP_MODEL_PULL + srcIndex = P_INDEX_5(gDims, px + oX, py + oY, pz + oZ, d); + // Put it on the other side back into the ghost layer. + dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); +#endif + + VerifyMsg(nBounceBackPdfs < kd->nBounceBackPdfs, "nBBPdfs %d < kd->nBBPdfs %d xyz: %d %d %d d: %d\n", nBounceBackPdfs, kd->nBounceBackPdfs, x, y, z, d); + + } + + kd->BounceBackPdfsSrc[nBounceBackPdfs] = srcIndex; + kd->BounceBackPdfsDst[nBounceBackPdfs] = dstIndex; + + ++nBounceBackPdfs; + + } + else if (nx < 0 || ny < 0 || nz < 0 || nx >= lX || ny >= lY || nz >= lZ) { + continue; + } + else if (ld->Lattice[L_INDEX_4(lDims, nx, ny, nz)] == LAT_CELL_OBSTACLE) { +#ifdef PROP_MODEL_PUSH + srcIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); + dstIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, D3Q19_INV[d]); +#elif PROP_MODEL_PULL + srcIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, D3Q19_INV[d]); + dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, d); + // srcIndex = P_INDEX_5(gDims, x + oX, y + oY, z + oZ, d); + // dstIndex = P_INDEX_5(gDims, nx + oX, ny + oY, nz + oZ, D3Q19_INV[d]); +#endif + + VerifyMsg(nBounceBackPdfs < kd->nBounceBackPdfs, "nBBPdfs %d < kd->nBBPdfs %d xyz: %d %d %d d: %d\n", nBounceBackPdfs, kd->nBounceBackPdfs, x, y, z, d); + + kd->BounceBackPdfsSrc[nBounceBackPdfs] = srcIndex; + kd->BounceBackPdfsDst[nBounceBackPdfs] = dstIndex; + + ++nBounceBackPdfs; + } + } + } + } + } + } + + + // Fill remaining KernelData structures + kd->GetNode = FNAME(GetNode); + kd->SetNode = FNAME(SetNode); + + kd->BoundaryConditionsGetPdf = FNAME(BcGetPdf); + kd->BoundaryConditionsSetPdf = FNAME(BcSetPdf); + + kd->Kernel = FNAME(D3Q19BlkKernel); + + kd->DstPdfs = NULL; + kd->PdfsActive = kd->Pdfs[0]; + + return; +} + +void FNAME(D3Q19BlkDeinit)(LatticeDesc * ld, KernelData ** kernelData) +{ + MemFree((void **) & ((*kernelData)->Pdfs[0])); + MemFree((void **) & ((*kernelData)->Pdfs[1])); + + MemFree((void **) & ((*kernelData)->BounceBackPdfsSrc)); + MemFree((void **) & ((*kernelData)->BounceBackPdfsDst)); + + MemFree((void **)kernelData); + + return; +} + +// Kernels without blocking perform the same initialization/deinitialization as with +// blocking, except that a different kernel is called. Hence, no arguments are allowed. + +void FNAME(D3Q19Init)(LatticeDesc * ld, KernelData ** kernelData, Parameters * params) +{ + Parameters p; + + if (params->nKernelArgs != 0) { + printf("ERROR: unknown kernel parameter.\n"); + printf("This kernels accepts no parameters.\n"); + exit(1); + } + + // Setup an empty parameters structure. + p.nArgs = params->nArgs; + p.Args = params->Args; + p.nKernelArgs = 0; + p.KernelArgs = NULL; + + // Call init routine for blocking kernel and override the + // kernel function to be called later on. + FNAME(D3Q19BlkInit)(ld, kernelData, &p); + + (*kernelData)->Kernel = FNAME(D3Q19Kernel); + + return; + +} + +void FNAME(D3Q19Deinit)(LatticeDesc * ld, KernelData ** kernelData) +{ + FNAME(D3Q19BlkDeinit)(ld, kernelData); + return; +}