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