X-Git-Url: http://git.rrze.uni-erlangen.de/gitweb/?p=LbmBenchmarkKernelsPublic.git;a=blobdiff_plain;f=src%2FVtk.c;fp=src%2FVtk.c;h=e650722a695add7aaf76337e2aeeb5514944f47b;hp=0000000000000000000000000000000000000000;hb=109880839321408644c94a34eb31208460b9f46d;hpb=42cf91486fb5c1ad178b3d21935a1be563e5fa39 diff --git a/src/Vtk.c b/src/Vtk.c new file mode 100644 index 0000000..e650722 --- /dev/null +++ b/src/Vtk.c @@ -0,0 +1,245 @@ +// -------------------------------------------------------------------------- +// +// 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 "Vtk.h" + +#include + +// TODO: make this portable + +// needed for stat & mkdir +#include +#include +#include +#include +#include // strerror + +// TODO: make byteswap portable + +#include +// glibc +#include + +// macros for portability +// #define BS32(a) bswap_32(*((uint32_t *)(&a))) +#define BS64(a) bswap_64(*((uint64_t *)(&a))) + + +void VtkWrite(LatticeDesc * ld, KernelData * kd, CaseData * cd, int iteration) +{ + Assert(kd != NULL); + Assert(ld != NULL); + Assert(ld->Dims[0] > 0); + Assert(ld->Dims[1] > 0); + Assert(ld->Dims[2] > 0); + + // TODO: this should be made portable... + // Check if subdirectory vtk exists, if not, create it. + { + int err; + struct stat fileStatus; + + err = stat("vtk", &fileStatus); + + if (err) { + // printf("ERROR: stat %d - %s\n", errno, strerror(errno)); + + // Set default mask and hope mkdir applies umask... + err = mkdir("vtk", 0700); + + if (err) { + printf("ERROR: cannot create directory vtk - %d: %s\n", errno, strerror(errno)); + exit(1); + } + + printf("# created directory vtk.\n"); + } + else { + + if (!S_ISDIR(fileStatus.st_mode)) { + printf("ERROR: cannot create subdirectory vtk as already a file with the same name exists.\n"); + exit(1); + } + + } + } + + + char fileName[1024]; + + snprintf(fileName, sizeof(fileName), "vtk/file-%04d.vtk", iteration); + + printf("# VTK: writing file %s\n", fileName); + + FILE * fh; + + fh = fopen(fileName, "w"); + + if(fh == NULL) { + printf("ERROR: opening file %s failed.\n", fileName); + exit(1); + } + + // http://www.vtk.org/pdf/file-formats.pdf + int nX = ld->Dims[0]; + int nY = ld->Dims[1]; + int nZ = ld->Dims[2]; + int * lDims = ld->Dims; + + // Temporaries for endian conversion. + uint64_t uDensity, uUx, uUy, uUz; + + PdfT pdfs[N_D3Q19]; + + fprintf(fh, "# vtk DataFile Version 1.0\n"); + fprintf(fh, "Comment: lid driven cavity, iteration % 4d\n", iteration); +#ifdef VTK_OUTPUT_ASCII + fprintf(fh, "ASCII\n"); +#else + fprintf(fh, "BINARY\n"); +#endif + fprintf(fh, "DATASET STRUCTURED_POINTS\n"); + fprintf(fh, "DIMENSIONS %d %d %d\n", nX, nY, nZ); + fprintf(fh, "ORIGIN 0 0 0 \n"); + fprintf(fh, "SPACING 1 1 1\n"); + fprintf(fh, "POINT_DATA %d\n", nX * nY * nZ); + + // ---------------------------------------------------------------------- + // Flag field: obstacle = 0, fluid = 1, inlet = 2, outlet = 4 + + fprintf(fh, "SCALARS NodesTypes unsigned_char 1\n"); + fprintf(fh, "LOOKUP_TABLE default\n"); + + unsigned char c; + + for(int z = 0; z < nZ; ++z) { + for(int y = 0; y < nY; ++y) { + for(int x = 0; x < nX; ++x) { +#ifdef VTK_OUTPUT_ASCII + fprintf(fh, "%d\n", ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)]); +#else + c = (unsigned char)ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)]; + fwrite(&c, sizeof(unsigned char), 1, fh); +#endif + } + } + } + + // ---------------------------------------------------------------------- + // Density field + + fprintf(fh, "SCALARS Density double\n"); + fprintf(fh, "LOOKUP_TABLE default\n"); + + double density; + + for(int z = 0; z < nZ; ++z) { + for(int y = 0; y < nY; ++y) { + for(int x = 0; x < nX; ++x) { + + density = 0.0; + if (ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) { + kd->GetNode(kd, x, y, z, pdfs); + + for (int d = 0; d < N_D3Q19; ++d) { + density += pdfs[d]; + } + } + +#ifdef VTK_OUTPUT_ASCII + fprintf(fh, "%e\n", density); +#else + uDensity = BS64(density); + fwrite(&uDensity, sizeof(double), 1, fh); +#endif + } + } + } + + // ---------------------------------------------------------------------- + // Velocity vectors: velocity in x, y, and z direction + + fprintf(fh, "VECTORS VelocityVectors double\n"); + + // Declare pdf_N, pdf_E, pdf_S, pdf_W, ... + #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name); + D3Q19_LIST + #undef X + + double ux, uy, uz; + + for(int z = 0; z < nZ; ++z) { + for(int y = 0; y < nY; ++y) { + for(int x = 0; x < nX; ++x) { + + if (ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) { + kd->GetNode(kd, x, y, z, pdfs); + +// DETECT NANS +// for (int d = 0; d < 19; ++d) { +// if(isnan(pdfs[d])) { +// printf("%d %d %d %d nan!\n", x, y, z, d); +// for (int d2 = 0; d2 < 19; ++d2) { +// printf("%d: %e\n", d2, pdfs[d2]); +// } +// exit(1); +// } +// } + #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = pdfs[idx]; + D3Q19_LIST + #undef X + UNUSED(pdf_C); + + + ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE - + pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW; + uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN - + pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS; + uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS - + pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS; + #ifdef VERIFICATION + ux += 0.5 * cd->XForce; + #endif + } + else { + ux = 0.0; uy = 0.0; uz = 0.0; + } + +#ifdef VTK_OUTPUT_ASCII + fprintf(fh, "%f %f %f\n", ux, uy, uz); +#else + uUx = BS64(ux); uUy = BS64(uy); uUz = BS64(uz); + fwrite(&uUx, sizeof(double), 1, fh); + fwrite(&uUy, sizeof(double), 1, fh); + fwrite(&uUz, sizeof(double), 1, fh); +#endif + } + } + } + + fclose(fh); +} +