src/Geometry.c

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