| 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 | |
| 35 | static const char * g_geoTypeStr[] = { "box", "channel", "pipe", "blocks", "fluid" }; |
| 36 | |
| 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 | } |
| 77 | else if (strncasecmp("fluid", geometryType, 5) == 0) { |
| 78 | type = GEO_TYPE_FLUID; |
| 79 | } |
| 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 | |
| 107 | // const char * geoTypeStr[] = { "box", "channel", "pipe", "blocks", "fluid" }; |
| 108 | |
| 109 | // printf("# geometry: %d x %d x %d nodes, type %d %s\n", dims[0], dims[1], dims[2], type, geoTypeStr[type]); |
| 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]; |
| 118 | ld->Name = g_geoTypeStr[type]; |
| 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 | } |
| 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 | } |
| 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) { |
| 181 | // Nothing todo here. |
| 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 | |
| 205 | #if 0 |
| 206 | if (blockSize == 0) { |
| 207 | blockSize = 8; |
| 208 | } |
| 209 | #endif |
| 210 | if (blockSize > 0) { |
| 211 | |
| 212 | int dimMin = dims[0]; |
| 213 | |
| 214 | if (dims[1] < dimMin) dimMin = dims[1]; |
| 215 | if (dims[2] < dimMin) dimMin = dims[2]; |
| 216 | |
| 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 | } |
| 222 | |
| 223 | // Number of blocks in x, y, and z direction. |
| 224 | int nbx = blockSize, nby = blockSize, nbz = blockSize; |
| 225 | |
| 226 | for (int z = 0; z < dims[2]; ++z) { |
| 227 | if ((z % (2 * nbz)) < nbz) continue; |
| 228 | |
| 229 | for (int y = 0; y < dims[1]; ++y) { |
| 230 | if ((y % (2 * nby)) < nby) continue; |
| 231 | |
| 232 | for (int x = 0; x < dims[0]; ++x) { |
| 233 | |
| 234 | if ((x % (2 * nbx)) >= nbx) { |
| 235 | lattice[L_INDEX_4(dims, x, y, z)] = LAT_CELL_OBSTACLE; |
| 236 | } |
| 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 | } |