src/Vtk.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 "Vtk.h"
  28
  29 #include <math.h>
  30
  31 // TODO: make this portable
  32
  33 // needed for stat & mkdir
  34 #include <sys/types.h>
  35 #include <sys/stat.h>
  36 #include <unistd.h>
  37 #include <errno.h>
  38 #include <string.h> // strerror
  39
  40 // TODO: make byteswap portable
  41
  42 #include <inttypes.h>
  43 // glibc
  44 #include <byteswap.h>
  45
  46 // macros for portability
  47 // #define BS32(a)              bswap_32(*((uint32_t *)(&a)))
  48 #define BS64(a)         bswap_64(*((uint64_t *)(&a)))
  49
  50
  51 void VtkWrite(LatticeDesc * ld, KernelData * kd, CaseData * cd, int iteration)
  52 {
  53         Assert(kd != NULL);
  54         Assert(ld != NULL);
  55         Assert(ld->Dims[0] > 0);
  56         Assert(ld->Dims[1] > 0);
  57         Assert(ld->Dims[2] > 0);
  58
  59         // TODO: this should be made portable...
  60         // Check if subdirectory vtk exists, if not, create it.
  61         {
  62                 int err;
  63                 struct stat fileStatus;
  64
  65                 err = stat("vtk", &fileStatus);
  66
  67                 if (err) {
  68                         // printf("ERROR: stat %d - %s\n", errno, strerror(errno));
  69
  70                         // Set default mask and hope mkdir applies umask...
  71                         err = mkdir("vtk", 0700);
  72
  73                         if (err) {
  74                                 printf("ERROR: cannot create directory vtk - %d: %s\n", errno, strerror(errno));
  75                                 exit(1);
  76                         }
  77
  78                         printf("# created directory vtk.\n");
  79                 }
  80                 else {
  81
  82                         if (!S_ISDIR(fileStatus.st_mode)) {
  83                                 printf("ERROR: cannot create subdirectory vtk as already a file with the same name exists.\n");
  84                                 exit(1);
  85                         }
  86
  87                 }
  88         }
  89
  90
  91         char fileName[1024];
  92
  93         snprintf(fileName, sizeof(fileName), "vtk/file-%04d.vtk", iteration);
  94
  95         printf("# VTK: writing file %s\n", fileName);
  96
  97         FILE * fh;
  98
  99         fh = fopen(fileName, "w");
 100
 101         if(fh == NULL) {
 102                 printf("ERROR: opening file %s failed.\n", fileName);
 103                 exit(1);
 104         }
 105
 106         // http://www.vtk.org/pdf/file-formats.pdf
 107         int nX = ld->Dims[0];
 108         int nY = ld->Dims[1];
 109         int nZ = ld->Dims[2];
 110         int * lDims = ld->Dims;
 111
 112         // Temporaries for endian conversion.
 113         uint64_t uDensity, uUx, uUy, uUz;
 114
 115         PdfT pdfs[N_D3Q19];
 116
 117         fprintf(fh, "# vtk DataFile Version 1.0\n");
 118         fprintf(fh, "Comment: lid driven cavity, iteration % 4d\n", iteration);
 119 #ifdef VTK_OUTPUT_ASCII
 120         fprintf(fh, "ASCII\n");
 121 #else
 122         fprintf(fh, "BINARY\n");
 123 #endif
 124         fprintf(fh, "DATASET STRUCTURED_POINTS\n");
 125         fprintf(fh, "DIMENSIONS %d %d %d\n", nX, nY, nZ);
 126         fprintf(fh, "ORIGIN 0 0 0 \n");
 127         fprintf(fh, "SPACING 1 1 1\n");
 128         fprintf(fh, "POINT_DATA %d\n", nX * nY * nZ);
 129
 130         // ----------------------------------------------------------------------
 131         // Flag field: obstacle = 0, fluid = 1, inlet = 2, outlet = 4
 132
 133         fprintf(fh, "SCALARS NodesTypes unsigned_char 1\n");
 134         fprintf(fh, "LOOKUP_TABLE default\n");
 135
 136         unsigned char c;
 137
 138         for(int z = 0; z < nZ; ++z) {
 139                 for(int y = 0; y < nY; ++y) {
 140                         for(int x = 0; x < nX; ++x) {
 141 #ifdef VTK_OUTPUT_ASCII
 142                                 fprintf(fh, "%d\n", ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)]);
 143 #else
 144                                 c = (unsigned char)ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)];
 145                                 fwrite(&c, sizeof(unsigned char), 1, fh);
 146 #endif
 147                         }
 148                 }
 149         }
 150
 151         // ----------------------------------------------------------------------
 152         // Density field
 153
 154         fprintf(fh, "SCALARS Density double\n");
 155         fprintf(fh, "LOOKUP_TABLE default\n");
 156
 157         double density;
 158
 159         for(int z = 0; z < nZ; ++z) {
 160                 for(int y = 0; y < nY; ++y) {
 161                         for(int x = 0; x < nX; ++x) {
 162
 163                                 density = 0.0;
 164                                 if (ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) {
 165                                         kd->GetNode(kd, x, y, z, pdfs);
 166
 167                                         for (int d = 0; d < N_D3Q19; ++d) {
 168                                                 density += pdfs[d];
 169                                         }
 170                                 }
 171
 172 #ifdef VTK_OUTPUT_ASCII
 173                                 fprintf(fh, "%e\n", density);
 174 #else
 175                                 uDensity = BS64(density);
 176                                 fwrite(&uDensity, sizeof(double), 1, fh);
 177 #endif
 178                         }
 179                 }
 180         }
 181
 182         // ----------------------------------------------------------------------
 183         // Velocity vectors: velocity in x, y, and z direction
 184
 185         fprintf(fh, "VECTORS VelocityVectors double\n");
 186
 187         // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
 188         #define X(name, idx, idxinv, x, y, z)   PdfT JOIN(pdf_,name);
 189         D3Q19_LIST
 190         #undef X
 191
 192         double ux, uy, uz;
 193
 194         for(int z = 0; z < nZ; ++z) {
 195                 for(int y = 0; y < nY; ++y) {
 196                         for(int x = 0; x < nX; ++x) {
 197
 198                                 if (ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) {
 199                                         kd->GetNode(kd, x, y, z, pdfs);
 200
 201 // DETECT NANS
 202 //                                      for (int d = 0; d < 19; ++d) {
 203 //                                              if(isnan(pdfs[d])) {
 204 //                                                      printf("%d %d %d %d nan!\n", x, y, z, d);
 205 //                                                      for (int d2 = 0; d2 < 19; ++d2) {
 206 //                                                              printf("%d: %e\n", d2, pdfs[d2]);
 207 //                                                      }
 208 //                                                      exit(1);
 209 //                                              }
 210 //                                      }
 211                                         #define X(name, idx, idxinv, _x, _y, _z)        JOIN(pdf_,name) = pdfs[idx];
 212                                         D3Q19_LIST
 213                                         #undef X
 214                                         UNUSED(pdf_C);
 215
 216
 217                                         ux = pdf_E + pdf_NE + pdf_SE + pdf_TE + pdf_BE -
 218                                                  pdf_W - pdf_NW - pdf_SW - pdf_TW - pdf_BW;
 219                                         uy = pdf_N + pdf_NE + pdf_NW + pdf_TN + pdf_BN -
 220                                                  pdf_S - pdf_SE - pdf_SW - pdf_TS - pdf_BS;
 221                                         uz = pdf_T + pdf_TE + pdf_TW + pdf_TN + pdf_TS -
 222                                                  pdf_B - pdf_BE - pdf_BW - pdf_BN - pdf_BS;
 223                                         #ifdef VERIFICATION
 224                                         ux += 0.5 * cd->XForce;
 225                                         #endif
 226                                 }
 227                                 else {
 228                                         ux = 0.0; uy = 0.0; uz = 0.0;
 229                                 }
 230
 231 #ifdef VTK_OUTPUT_ASCII
 232                                 fprintf(fh, "%f %f %f\n", ux, uy, uz);
 233 #else
 234                                 uUx = BS64(ux); uUy = BS64(uy); uUz = BS64(uz);
 235                                 fwrite(&uUx, sizeof(double), 1, fh);
 236                                 fwrite(&uUy, sizeof(double), 1, fh);
 237                                 fwrite(&uUz, sizeof(double), 1, fh);
 238 #endif
 239                         }
 240                 }
 241         }
 242
 243         fclose(fh);
 244 }
 245