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 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 }