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 // --------------------------------------------------------------------------
31 // TODO: make this portable
33 // needed for stat & mkdir
34 #include <sys/types.h>
38 #include <string.h> // strerror
40 // TODO: make byteswap portable
46 // macros for portability
47 // #define BS32(a) bswap_32(*((uint32_t *)(&a)))
48 #define BS64(a) bswap_64(*((uint64_t *)(&a)))
51 void VtkWrite(LatticeDesc * ld, KernelData * kd, CaseData * cd, int iteration)
55 Assert(ld->Dims[0] > 0);
56 Assert(ld->Dims[1] > 0);
57 Assert(ld->Dims[2] > 0);
59 // TODO: this should be made portable...
60 // Check if subdirectory vtk exists, if not, create it.
63 struct stat fileStatus;
65 err = stat("vtk", &fileStatus);
68 // printf("ERROR: stat %d - %s\n", errno, strerror(errno));
70 // Set default mask and hope mkdir applies umask...
71 err = mkdir("vtk", 0700);
74 printf("ERROR: cannot create directory vtk - %d: %s\n", errno, strerror(errno));
78 printf("# created directory vtk.\n");
82 if (!S_ISDIR(fileStatus.st_mode)) {
83 printf("ERROR: cannot create subdirectory vtk as already a file with the same name exists.\n");
93 snprintf(fileName, sizeof(fileName), "vtk/file-%04d.vtk", iteration);
95 printf("# VTK: writing file %s\n", fileName);
99 fh = fopen(fileName, "w");
102 printf("ERROR: opening file %s failed.\n", fileName);
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;
112 // Temporaries for endian conversion.
113 uint64_t uDensity, uUx, uUy, uUz;
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");
122 fprintf(fh, "BINARY\n");
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);
130 // ----------------------------------------------------------------------
131 // Flag field: obstacle = 0, fluid = 1, inlet = 2, outlet = 4
133 fprintf(fh, "SCALARS NodesTypes unsigned_char 1\n");
134 fprintf(fh, "LOOKUP_TABLE default\n");
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)]);
144 c = (unsigned char)ld->Lattice[L_INDEX_4(ld->Dims, x, y, z)];
145 fwrite(&c, sizeof(unsigned char), 1, fh);
151 // ----------------------------------------------------------------------
154 fprintf(fh, "SCALARS Density double\n");
155 fprintf(fh, "LOOKUP_TABLE default\n");
159 for(int z = 0; z < nZ; ++z) {
160 for(int y = 0; y < nY; ++y) {
161 for(int x = 0; x < nX; ++x) {
164 if (ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) {
165 kd->GetNode(kd, x, y, z, pdfs);
167 for (int d = 0; d < N_D3Q19; ++d) {
172 #ifdef VTK_OUTPUT_ASCII
173 fprintf(fh, "%e\n", density);
175 uDensity = BS64(density);
176 fwrite(&uDensity, sizeof(double), 1, fh);
182 // ----------------------------------------------------------------------
183 // Velocity vectors: velocity in x, y, and z direction
185 fprintf(fh, "VECTORS VelocityVectors double\n");
187 // Declare pdf_N, pdf_E, pdf_S, pdf_W, ...
188 #define X(name, idx, idxinv, x, y, z) PdfT JOIN(pdf_,name);
194 for(int z = 0; z < nZ; ++z) {
195 for(int y = 0; y < nY; ++y) {
196 for(int x = 0; x < nX; ++x) {
198 if (ld->Lattice[L_INDEX_4(lDims, x, y, z)] != LAT_CELL_OBSTACLE) {
199 kd->GetNode(kd, x, y, z, pdfs);
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]);
211 #define X(name, idx, idxinv, _x, _y, _z) JOIN(pdf_,name) = pdfs[idx];
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;
224 ux += 0.5 * cd->XForce;
228 ux = 0.0; uy = 0.0; uz = 0.0;
231 #ifdef VTK_OUTPUT_ASCII
232 fprintf(fh, "%f %f %f\n", ux, uy, uz);
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);