Commit | Line | Data |
---|---|---|
10988083 MW |
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 |