| 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 | |