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 "Geometry.h" | |
28 | #include "Memory.h" | |
29 | ||
30 | #include <strings.h> | |
31 | #include <math.h> | |
32 | ||
33 | #include <errno.h> | |
34 | ||
0fde6e45 MW |
35 | static const char * g_geoTypeStr[] = { "box", "channel", "pipe", "blocks", "fluid" }; |
36 | ||
10988083 MW |
37 | void GeoCreateByStr(const char * geometryType, int dims[3], int periodic[3], LatticeDesc * ld) |
38 | { | |
39 | int type = -1; | |
40 | void * typeDetails = NULL; | |
41 | int tmp; | |
42 | ||
43 | if (strncasecmp("channel", geometryType, 7) == 0) { | |
44 | type = GEO_TYPE_CHANNEL; | |
45 | } | |
46 | else if (strncasecmp("box", geometryType, 3) == 0) { | |
47 | type = GEO_TYPE_BOX; | |
48 | } | |
49 | else if (strncasecmp("pipe", geometryType, 4) == 0) { | |
50 | type = GEO_TYPE_PIPE; | |
51 | } | |
52 | else if (strncasecmp("blocks", geometryType, 6) == 0) { | |
53 | type = GEO_TYPE_BLOCKS; | |
54 | ||
55 | // Default block size | |
56 | tmp = 8; | |
57 | ||
58 | if (strlen(geometryType) > 7) { | |
59 | int blockSize = atoi(&geometryType[7]); | |
60 | ||
61 | int dimMin = dims[0]; | |
62 | ||
63 | if (dims[1] < dimMin) dimMin = dims[1]; | |
64 | if (dims[2] < dimMin) dimMin = dims[2]; | |
65 | ||
66 | if (blockSize < 0 || blockSize > dimMin / 2) { | |
67 | printf("ERROR: block size for geometry must be > 0 and smaller than half of the smalest dimension.\n"); | |
68 | // TODO: find a better solution for handling errors in here. | |
69 | Verify(0); | |
70 | } | |
71 | ||
72 | tmp = blockSize; | |
73 | } | |
74 | ||
75 | typeDetails = &tmp; | |
76 | } | |
e3f82424 MW |
77 | else if (strncasecmp("fluid", geometryType, 5) == 0) { |
78 | type = GEO_TYPE_FLUID; | |
79 | } | |
10988083 MW |
80 | else { |
81 | printf("ERROR: unknown geometry specified.\n"); | |
82 | Verify(0); | |
83 | } | |
84 | ||
85 | GeoCreateByType(type, typeDetails, dims, periodic, ld); | |
86 | ||
87 | return; | |
88 | } | |
89 | ||
90 | void GeoCreateByType(GEO_TYPES type, void * typeDetails, int dims[3], int periodic[3], LatticeDesc * ld) | |
91 | { | |
92 | Assert(dims != NULL); | |
93 | Assert(dims[0] > 0); | |
94 | Assert(dims[1] > 0); | |
95 | Assert(dims[2] > 0); | |
96 | ||
97 | Assert(periodic != NULL); | |
98 | Assert(periodic[0] >= 0); | |
99 | Assert(periodic[1] >= 0); | |
100 | Assert(periodic[2] >= 0); | |
101 | ||
102 | Assert(ld != NULL); | |
103 | ||
104 | Assert(type >= GEO_TYPE_MIN); | |
105 | Assert(type <= GEO_TYPE_MAX); | |
106 | ||
0fde6e45 | 107 | // const char * geoTypeStr[] = { "box", "channel", "pipe", "blocks", "fluid" }; |
10988083 | 108 | |
0fde6e45 | 109 | // printf("# geometry: %d x %d x %d nodes, type %d %s\n", dims[0], dims[1], dims[2], type, geoTypeStr[type]); |
10988083 MW |
110 | |
111 | ld->Dims[0] = dims[0]; | |
112 | ld->Dims[1] = dims[1]; | |
113 | ld->Dims[2] = dims[2]; | |
114 | ld->nCells = dims[0] * dims[1] * dims[2]; | |
115 | ld->PeriodicX = periodic[0]; | |
116 | ld->PeriodicY = periodic[1]; | |
117 | ld->PeriodicZ = periodic[2]; | |
0fde6e45 | 118 | ld->Name = g_geoTypeStr[type]; |
10988083 MW |
119 | |
120 | LatticeT * lattice; | |
121 | MemAlloc((void **)&lattice, sizeof(LatticeT) * dims[0] * dims[1] * dims[2]); | |
122 | ||
123 | ld->Lattice = lattice; | |
124 | ||
125 | for (int z = 0; z < dims[2]; ++z) { | |
126 | for (int y = 0; y < dims[1]; ++y) { | |
127 | for (int x = 0; x < dims[0]; ++x) { | |
128 | lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_FLUID; | |
129 | } | |
130 | } | |
131 | } | |
132 | ||
133 | if (type == GEO_TYPE_CHANNEL || type == GEO_TYPE_BLOCKS || type == GEO_TYPE_PIPE) { | |
134 | periodic[0] = 1; | |
135 | } | |
e3f82424 MW |
136 | else if (type == GEO_TYPE_FLUID) { |
137 | // Actually this geometry is non-periodic, but we need fluid at the boundary. | |
138 | periodic[0] = 1; periodic[1] = 1; periodic[2] = 1; | |
139 | } | |
10988083 MW |
140 | |
141 | // Walls or periodic on first and last x plane. | |
142 | for (int z = 0; z < dims[2]; ++z) { | |
143 | for (int y = 0; y < dims[1]; ++y) { | |
144 | if(periodic[0]){ | |
145 | lattice[L_INDEX_4(dims, 0, y, z)] = LAT_CELL_FLUID; | |
146 | lattice[L_INDEX_4(dims, dims[0] - 1, y, z)] = LAT_CELL_FLUID; | |
147 | } else { | |
148 | lattice[L_INDEX_4(dims, 0, y, z)] = LAT_CELL_OBSTACLE; | |
149 | lattice[L_INDEX_4(dims, dims[0] - 1, y, z)] = LAT_CELL_OBSTACLE; | |
150 | } | |
151 | } | |
152 | } | |
153 | ||
154 | // Walls or periodic on first and last y plane. | |
155 | for (int z = 0; z < dims[2]; ++z) { | |
156 | for (int x = 0; x < dims[0]; ++x) { | |
157 | if(periodic[1]){ | |
158 | lattice[L_INDEX_4(dims, x, 0, z)] = LAT_CELL_FLUID; | |
159 | lattice[L_INDEX_4(dims, x, dims[1] - 1, z)] = LAT_CELL_FLUID; | |
160 | } else { | |
161 | lattice[L_INDEX_4(dims, x, 0, z)] = LAT_CELL_OBSTACLE; | |
162 | lattice[L_INDEX_4(dims, x, dims[1] - 1, z)] = LAT_CELL_OBSTACLE; | |
163 | } | |
164 | } | |
165 | } | |
166 | ||
167 | // Walls or periodic on first and last z plane. | |
168 | for (int y = 0; y < dims[1]; ++y) { | |
169 | for (int x = 0; x < dims[0]; ++x) { | |
170 | if(periodic[2]){ | |
171 | lattice[L_INDEX_4(dims, x, y, 0)] = LAT_CELL_FLUID; | |
172 | lattice[L_INDEX_4(dims, x, y, dims[2] - 1)] = LAT_CELL_FLUID; | |
173 | } else { | |
174 | lattice[L_INDEX_4(dims, x, y, 0)] = LAT_CELL_OBSTACLE; | |
175 | lattice[L_INDEX_4(dims, x, y, dims[2] - 1)] = LAT_CELL_OBSTACLE; | |
176 | } | |
177 | } | |
178 | } | |
179 | ||
180 | if (type == GEO_TYPE_CHANNEL) { | |
e3f82424 | 181 | // Nothing todo here. |
10988083 MW |
182 | } |
183 | else if (type == GEO_TYPE_PIPE) { | |
184 | #define SQR(a) ((a)*(a)) | |
185 | double centerZ = dims[2] / 2.0 - 0.5; | |
186 | double centerY = dims[1] / 2.0 - 0.5; | |
187 | double minDiameter = MIN(dims[1], dims[2]); | |
188 | double minRadiusSquared = SQR(minDiameter / 2 - 1); | |
189 | ||
190 | for (int z = 0; z < dims[2]; ++z) { | |
191 | for (int y = 0; y < dims[1]; ++y) { | |
192 | if((SQR(z - centerZ) + SQR(y - centerY)) >= minRadiusSquared) { | |
193 | for (int x = 0; x < dims[0]; ++x) { | |
194 | lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_OBSTACLE; | |
195 | } | |
196 | } | |
197 | } | |
198 | } | |
199 | #undef SQR | |
200 | } | |
201 | else if (type == GEO_TYPE_BLOCKS) { | |
202 | ||
203 | int blockSize = *((int *)typeDetails); | |
204 | ||
e3f82424 | 205 | #if 0 |
10988083 MW |
206 | if (blockSize == 0) { |
207 | blockSize = 8; | |
208 | } | |
e3f82424 MW |
209 | #endif |
210 | if (blockSize > 0) { | |
10988083 | 211 | |
e3f82424 | 212 | int dimMin = dims[0]; |
10988083 | 213 | |
e3f82424 MW |
214 | if (dims[1] < dimMin) dimMin = dims[1]; |
215 | if (dims[2] < dimMin) dimMin = dims[2]; | |
10988083 | 216 | |
e3f82424 MW |
217 | if (blockSize < 0 || blockSize > dimMin / 2) { |
218 | printf("ERROR: block size for geometry must be > 0 and smaller than half of the smalest dimension.\n"); | |
219 | // TODO: find a better solution for handling errors in here. | |
220 | Verify(0); | |
221 | } | |
10988083 | 222 | |
e3f82424 MW |
223 | // Number of blocks in x, y, and z direction. |
224 | int nbx = blockSize, nby = blockSize, nbz = blockSize; | |
10988083 | 225 | |
e3f82424 MW |
226 | for (int z = 0; z < dims[2]; ++z) { |
227 | if ((z % (2 * nbz)) < nbz) continue; | |
10988083 | 228 | |
e3f82424 MW |
229 | for (int y = 0; y < dims[1]; ++y) { |
230 | if ((y % (2 * nby)) < nby) continue; | |
10988083 | 231 | |
e3f82424 | 232 | for (int x = 0; x < dims[0]; ++x) { |
10988083 | 233 | |
e3f82424 MW |
234 | if ((x % (2 * nbx)) >= nbx) { |
235 | lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_OBSTACLE; | |
236 | } | |
10988083 MW |
237 | } |
238 | } | |
239 | } | |
240 | } | |
241 | } | |
242 | ||
243 | // if (latticeDumpAscii) { | |
244 | // const char strLatCellType[] = "X.IxO"; // X = Obstacle, . = Fluid, I = inlet, O = outlet | |
245 | // for (int z = dims[2] - 1; z >= 0; --z) { | |
246 | // printf("plane % 2d\n", z); | |
247 | // | |
248 | // for (int y = dims[1] - 1; y >= 0; --y) { | |
249 | // printf(" %2d ", y); | |
250 | // for (int x = 0; x < dims[0]; ++x) { | |
251 | // printf("%c", strLatCellType[lattice[L_INDEX_4(dims, x, y, z)]]); | |
252 | // } | |
253 | // printf("\n"); | |
254 | // } | |
255 | // } | |
256 | // } | |
257 | ||
258 | // Lattice Helper Function | |
259 | ||
260 | ld->nObst = 0; | |
261 | ld->nFluid = 0; | |
262 | ld->nInlet = 0; | |
263 | ld->nOutlet = 0; | |
264 | ||
265 | for (int z = 0; z < dims[2]; ++z) { | |
266 | for (int y = 0; y < dims[1]; ++y) { | |
267 | for (int x = 0; x < dims[0]; ++x) { | |
268 | switch (lattice[L_INDEX_4(dims, x, y, z)]) { | |
269 | case LAT_CELL_OBSTACLE: ld->nObst++; break; | |
270 | case LAT_CELL_FLUID: ld->nFluid++; break; | |
271 | case LAT_CELL_INLET: ld->nInlet++; ld->nFluid++; break; | |
272 | case LAT_CELL_OUTLET: ld->nOutlet++; ld->nFluid++; break; | |
273 | default: | |
274 | Verify(0); | |
275 | } | |
276 | } | |
277 | } | |
278 | } | |
279 | ||
280 | return; | |
281 | } |