apps/showcases/Channel/ConvergenceLin.txt

+/Users/holzer/walberla/cmake-build-release/apps/showcases/Channel/Channel Channel.prm
+[0][INFO    ]------(0.000 sec) Creating the block structure ...
+[0][INFO    ]------(0.000 sec) SetupBlockForest created successfully:
+[0]                            - AABB: [ <0,0,0>, <8,8,4> ]
+[0]                            - initial decomposition: 2 x 2 x 1 (= forest size)
+[0]                            - periodicity: true x false x true
+[0]                            - number of blocks discarded from the initial grid: 0 (= 0 %)
+[0]                            - number of levels: 2
+[0]                            - tree ID digits: 3 (-> block ID bytes = 1)
+[0]                            - total number of blocks: 18
+[0]                            - number of processes: 1 (1 worker process(es) / 0 empty buffer process(es))
+[0]                            - buffer processes are inserted into the process network: no
+[0]                            - process ID bytes: 0
+[0]                            - blocks/memory/workload per process:
+[0]                               + blocks:
+[0]                                  - min       = 18
+[0]                                  - max       = 18
+[0]                                  - avg       = 18
+[0]                                  - stdDev    = 0
+[0]                                  - relStdDev = 0
+[0]                               + memory:
+[0]                                  - min       = 4.35917e+06
+[0]                                  - max       = 4.35917e+06
+[0]                                  - avg       = 4.35917e+06
+[0]                                  - stdDev    = 0
+[0]                                  - relStdDev = 0
+[0]                               + workload:
+[0]                                  - min       = 34
+[0]                                  - max       = 34
+[0]                                  - avg       = 34
+[0]                                  - stdDev    = 0
+[0]                                  - relStdDev = 0
+[0]                            - distribution of space/memory/work to different grid levels:
+[0]                               + level 0:
+[0]                                  - 2 blocks ...
+[0]                                  - ... cover 50 % of the total simulation space
+[0]                                  - ... account for 11.1111 % of the total memory foot print
+[0]                                  - ... generate 5.88235 % of the total work load
+[0]                                  - blocks per process:
+[0]                                     + min       = 2
+[0]                                     + max       = 2
+[0]                                     + avg       = 2
+[0]                                     + stdDev    = 0
+[0]                                     + relStdDev = 0
+[0]                                  - memory per process:
+[0]                                     + min       = 484352
+[0]                                     + max       = 484352
+[0]                                     + avg       = 484352
+[0]                                     + stdDev    = 0
+[0]                                     + relStdDev = 0
+[0]                                  - workload per process:
+[0]                                     + min       = 2
+[0]                                     + max       = 2
+[0]                                     + avg       = 2
+[0]                                     + stdDev    = 0
+[0]                                     + relStdDev = 0
+[0]                               + level 1:
+[0]                                  - 16 blocks ...
+[0]                                  - ... cover 50 % of the total simulation space
+[0]                                  - ... account for 88.8889 % of the total memory foot print
+[0]                                  - ... generate 94.1176 % of the total work load
+[0]                                  - blocks per process:
+[0]                                     + min       = 16
+[0]                                     + max       = 16
+[0]                                     + avg       = 16
+[0]                                     + stdDev    = 0
+[0]                                     + relStdDev = 0
+[0]                                  - memory per process:
+[0]                                     + min       = 3.87482e+06
+[0]                                     + max       = 3.87482e+06
+[0]                                     + avg       = 3.87482e+06
+[0]                                     + stdDev    = 0
+[0]                                     + relStdDev = 0
+[0]                                  - workload per process:
+[0]                                     + min       = 32
+[0]                                     + max       = 32
+[0]                                     + avg       = 32
+[0]                                     + stdDev    = 0
+[0]                                     + relStdDev = 0
+[0]                            - using a uniform decomposition with a resolution equal to the finest level, one would ...
+[0]                               + ... need 1.77778 times the memory
+[0]                               + ... generate 1.88235 times the workload
+[0][INFO    ]------(0.003 sec) Setting up communication...
+[0][INFO    ]------(0.008 sec) Blocks created: 18
+[0][INFO    ]------(0.008 sec) Level 0 Blocks: 2
+[0][INFO    ]------(0.008 sec) Level 1 Blocks: 16
+[0][INFO    ]------(0.008 sec) Starting Simulation
+[0][INFO    ]------(0.013 sec)[   0] Evaluation of accuracy:
+[0]                                   - L1:   0.00479167
+[0]                                   - L2:   0.005469
+[0]                                   - Lmax: 0.0090625
+[0][INFO    ]------(0.435 sec)[ 100] Evaluation of accuracy:
+[0]                                   - L1:   0.00134704
+[0]                                   - L2:   0.00149049
+[0]                                   - Lmax: 0.002099
+[0][INFO    ]------(0.851 sec)[ 200] Evaluation of accuracy:
+[0]                                   - L1:   0.000448633
+[0]                                   - L2:   0.000496303
+[0]                                   - Lmax: 0.000696362
+[0][INFO    ]------(1.262 sec)[ 300] Evaluation of accuracy:
+[0]                                   - L1:   0.000149393
+[0]                                   - L2:   0.000165266
+[0]                                   - Lmax: 0.00023185
+[0][INFO    ]------(1.673 sec)[ 400] Evaluation of accuracy:
+[0]                                   - L1:   4.97466e-05
+[0]                                   - L2:   5.50323e-05
+[0]                                   - Lmax: 7.7204e-05
+[0][INFO    ]------(2.085 sec)[ 500] Evaluation of accuracy:
+[0]                                   - L1:   1.65652e-05
+[0]                                   - L2:   1.83253e-05
+[0]                                   - Lmax: 2.57083e-05
+[0][INFO    ]------(2.507 sec)[ 600] Evaluation of accuracy:
+[0]                                   - L1:   5.5161e-06
+[0]                                   - L2:   6.1022e-06
+[0]                                   - Lmax: 8.56068e-06
+[0][INFO    ]------(2.919 sec)[ 700] Evaluation of accuracy:
+[0]                                   - L1:   1.83682e-06
+[0]                                   - L2:   2.03199e-06
+[0]                                   - Lmax: 2.85064e-06
+[0][INFO    ]------(3.330 sec)[ 800] Evaluation of accuracy:
+[0]                                   - L1:   6.11647e-07
+[0]                                   - L2:   6.76636e-07
+[0]                                   - Lmax: 9.49242e-07
+[0][INFO    ]------(3.741 sec)[ 900] Evaluation of accuracy:
+[0]                                   - L1:   2.03674e-07
+[0]                                   - L2:   2.25315e-07
+[0]                                   - Lmax: 3.1609e-07
+[0][INFO    ]------(4.153 sec)[1000] Evaluation of accuracy:
+[0]                                   - L1:   6.78218e-08
+[0]                                   - L2:   7.5028e-08
+[0]                                   - Lmax: 1.05256e-07
+[0][INFO    ]------(4.162 sec) Simulation finished
+[0][RESULT  ]------(4.162 sec) Simulation performance:
+[0]                            - processes:   1
+[0]                            - threads:     1 (threads per process = 1, threads per core = n/a *))
+[0]                            - cores:       n/a *)
+[0]                            - time steps:  1001 (on the coarsest grid, 2002 on the finest grid)
+[0]                            - time:        4.1544 sec
+[0]                            - cells:       1152 (2048 if everything were fine -> data reduction by factor of 1.77778)
+[0]                            - fluid cells: 1152 (100 % of all cells)
+[0]                            - distribution of cells to different grid levels:
+[0]                               + level 0: 128 cells (128 fluid cells = 100 % of all cells on this level)
+[0]                               + level 1: 1024 cells (1024 fluid cells = 100 % of all cells on this level)
+[0]                            - performance: 0.524305 MLUPS (million lattice cell updates per second)
+[0]                                           0.524305 MLUPS / process
+[0]                                           n/a *) MLUPS / core
+[0]                                           0.524305 MFLUPS (million fluid lattice cell updates per second)
+[0]                                           0.524305 MFLUPS / process
+[0]                                           n/a *) MFLUPS / core
+[0]                                           240.949 time steps / second
+[0]                            - 'virtual' performance (if everything were fine): 0.986928 MLUPS (million lattice cell updates per second)
+[0]                                                                               0.986928 MLUPS / process
+[0]                                                                               n/a *) MLUPS / core
+[0]                                                                               0.986928 MFLUPS (million fluid lattice cell updates per second)
+[0]                                                                               0.986928 MFLUPS / process
+[0]                                                                               n/a *) MFLUPS / core
+[0]                                                                               481.898 fine time steps / second
+[0]                            - build / run information:
+[0]                               + host machine:   gazleu.cerfacs.fr
+[0]                               + build machine:  gazleu.cerfacs.fr
+[0]                               + git SHA1:       d5dba67210e6cc8141af6ce4b87c70822dca19e9
+[0]                               + build type:     Release
+[0]                               + compiler flags: -O3 -DNDEBUG  -Wall -Wconversion -Wshadow -Wno-c++11-extensions -Qunused-arguments -ffast-math -pthread
+[0]
+[0]                              *) only available if environment variable 'THREADS_PER_CORE' is set
+[0][RESULT  ]------(4.162 sec) Time loop timing:
+[0]                            Timer                   |       % |          Total|        Average|          Count|            Min|            Max|       Variance|
+[0]                            -----------------------------------------------------------------------------------------------------------------------------------
+[0]                            Refinement Cycle        |  99.24% |           4.12|          0.004|           1001|          0.004|          0.005|          0.000|
+[0]                            VTK Output              |   0.55% |           0.02|          0.000|           1001|          0.000|          0.009|          0.000|
+[0]                            evaluation              |   0.20% |           0.01|          0.000|           1001|          0.000|          0.001|          0.000|
+[0]                            remaining time logger   |   0.00% |           0.00|          0.000|           1001|          0.000|          0.000|          0.000|
+[0]
+[0][WARNING ]------(4.162 sec) Getting memory statistics is currently not supported on non-Linux systems!
+[0][INFO    ]------(4.162 sec) resident memory: Sample has 1 values in [0.000000, 0.000000], sum = 0.000000, mean = 0.000000, med = 0.000000, stddev = 0.000000 (relative: nan), mad = 0.000000
+
+Process finished with exit code 0
\ No newline at end of file

apps/showcases/FlowAroundCylinder/CMakeLists.txt

+waLBerla_link_files_to_builddir( *.py )
+waLBerla_link_files_to_builddir( *.prm )
+waLBerla_link_files_to_builddir( *.png )
+waLBerla_link_files_to_builddir( *.obj )
+waLBerla_link_files_to_builddir( *.stl )
+waLBerla_link_files_to_builddir( *.mtl )
+
+waLBerla_generate_target_from_python(NAME FlowAroundCylinderCodeGen
+        FILE FlowAroundCylinder.py
+        OUT_FILES FlowAroundCylinderSweepCollection.h FlowAroundCylinderSweepCollection.${CODEGEN_FILE_SUFFIX}
+        FlowAroundCylinderStorageSpecification.h FlowAroundCylinderStorageSpecification.${CODEGEN_FILE_SUFFIX}
+        FreeSlip.h FreeSlip.${CODEGEN_FILE_SUFFIX}
+        NoSlip.h NoSlip.${CODEGEN_FILE_SUFFIX}
+        Obstacle.h Obstacle.${CODEGEN_FILE_SUFFIX}
+        UBB.h UBB.${CODEGEN_FILE_SUFFIX}
+        Outflow.h Outflow.${CODEGEN_FILE_SUFFIX}
+        FlowAroundCylinderBoundaryCollection.h
+        FlowAroundCylinderInfoHeader.h
+        FlowAroundCylinderStaticDefines.h)
+
+if (WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)
+    waLBerla_add_executable ( NAME FlowAroundCylinder
+            FILES FlowAroundCylinder.cpp Cylinder.cpp Evaluation.cpp
+            DEPENDS blockforest boundary core field gpu lbm_generated stencil timeloop vtk FlowAroundCylinderCodeGen )
+else()
+    waLBerla_add_executable ( NAME FlowAroundCylinder
+            FILES FlowAroundCylinder.cpp Cylinder.cpp Evaluation.cpp
+            DEPENDS blockforest boundary core field lbm_generated stencil timeloop vtk FlowAroundCylinderCodeGen )
+
+endif(WALBERLA_BUILD_WITH_CUDA OR WALBERLA_BUILD_WITH_HIP)

apps/showcases/FlowAroundCylinder/Cylinder.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Cylinder.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "Cylinder.h"
+
+namespace walberla
+{
+
+bool Cylinder::contains(const Vector3< real_t >& point) const
+{
+   const real_t px = setup_.cylinderXPosition;
+   const real_t py = setup_.cylinderYPosition;
+   const real_t H  = setup_.H;
+   const real_t r  = setup_.cylinderRadius;
+
+   if (setup_.circularCrossSection)
+   {
+      const real_t xd = point[0] - px;
+      const real_t yd = point[1] - py;
+      const real_t d  = xd * xd + yd * yd;
+      return point[2] >= real_t(0.0) && point[2] <= H && d <= (r * r);
+   }
+   else
+   {
+      const AABB cylinder(px - r, py - r, real_t(0), px + r, py + r, H);
+      return cylinder.contains(point);
+   }
+}
+
+bool Cylinder::contains(const AABB& aabb) const
+{
+   if (setup_.circularCrossSection)
+   {
+      Vector3< real_t > p[8];
+      p[0].set(aabb.xMin(), aabb.yMin(), aabb.zMin());
+      p[1].set(aabb.xMax(), aabb.yMin(), aabb.zMin());
+      p[2].set(aabb.xMin(), aabb.yMax(), aabb.zMin());
+      p[3].set(aabb.xMax(), aabb.yMax(), aabb.zMin());
+      p[4].set(aabb.xMin(), aabb.yMin(), aabb.zMax());
+      p[5].set(aabb.xMax(), aabb.yMin(), aabb.zMax());
+      p[6].set(aabb.xMin(), aabb.yMax(), aabb.zMax());
+      p[7].set(aabb.xMax(), aabb.yMax(), aabb.zMax());
+      return contains(p[0]) && contains(p[1]) && contains(p[2]) && contains(p[3]) && contains(p[4]) && contains(p[5]) &&
+             contains(p[6]) && contains(p[7]);
+   }
+   else { return contains(aabb.min()) && contains(aabb.max()); }
+}
+
+bool Cylinder::intersects(const AABB& aabb, const real_t bufferDistance) const
+{
+   const real_t px = setup_.cylinderXPosition;
+   const real_t py = setup_.cylinderYPosition;
+   const real_t r  = setup_.cylinderRadius;
+
+   if (setup_.circularCrossSection)
+   {
+      Vector3< real_t > p(px, py, real_t(0));
+
+      if (p[0] < aabb.xMin())
+         p[0] = aabb.xMin();
+      else if (p[0] > aabb.xMax())
+         p[0] = aabb.xMax();
+      if (p[1] < aabb.yMin())
+         p[1] = aabb.yMin();
+      else if (p[1] > aabb.yMax())
+         p[1] = aabb.yMax();
+
+      const real_t xd = p[0] - px;
+      const real_t yd = p[1] - py;
+      const real_t d  = xd * xd + yd * yd;
+      const real_t rr = (r + bufferDistance) * (r + bufferDistance);
+      return d <= rr;
+   }
+   else
+   {
+      const AABB cylinder(px - r, py - r, real_t(0), px + r, py + r, setup_.H);
+      return cylinder.intersects(aabb, bufferDistance);
+   }
+}
+
+real_t Cylinder::delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const
+{
+   WALBERLA_CHECK(!contains(fluid), "Point ( " << fluid << " ) is not outside the cylinder!")
+   WALBERLA_CHECK(contains(boundary), "Point ( " << boundary << " ) is not inside the cylinder!")
+
+   const real_t px = setup_.cylinderXPosition;
+   const real_t py = setup_.cylinderYPosition;
+   const real_t r  = setup_.cylinderRadius;
+
+   if (setup_.circularCrossSection)
+   {
+      // http://devmag.org.za/2009/04/17/basic-collision-detection-in-2d-part-2/
+
+      const Vector3< real_t > circle(px, py, real_t(0));
+
+      const Vector3< real_t > f = fluid - circle;
+      const Vector3< real_t > d = (boundary - circle) - f;
+
+      const real_t a = d[0] * d[0] + d[1] * d[1];
+      const real_t b = real_t(2) * (d[0] * f[0] + d[1] * f[1]);
+      const real_t c = f[0] * f[0] + f[1] * f[1] - r * r;
+
+      const real_t bb4ac = b * b - (real_t(4) * a * c);
+      WALBERLA_CHECK_GREATER_EQUAL(bb4ac, real_t(0))
+
+      const real_t sqrtbb4ac = std::sqrt(bb4ac);
+
+      const real_t alpha = std::min((-b + sqrtbb4ac) / (real_t(2) * a), (-b - sqrtbb4ac) / (real_t(2) * a));
+
+      WALBERLA_CHECK_GREATER_EQUAL(alpha, real_t(0))
+      WALBERLA_CHECK_LESS_EQUAL(alpha, real_t(1))
+
+      return alpha;
+   }
+
+   const AABB cylinder(px - r, py - r, real_t(0), px + r, py + r, setup_.H);
+
+   if (fluid[0] <= cylinder.xMin())
+   {
+      const real_t xdiff = cylinder.xMin() - fluid[0];
+
+      if (fluid[1] <= cylinder.yMin())
+      {
+         const real_t ydiff = cylinder.yMin() - fluid[1];
+         if (xdiff >= ydiff)
+         {
+            WALBERLA_CHECK_LESS_EQUAL(fluid[0], boundary[0])
+            return xdiff / (boundary[0] - fluid[0]);
+         }
+         WALBERLA_CHECK_LESS_EQUAL(fluid[1], boundary[1])
+         return ydiff / (boundary[1] - fluid[1]);
+      }
+      else if (fluid[1] >= cylinder.yMax())
+      {
+         const real_t ydiff = fluid[1] - cylinder.yMax();
+         if (xdiff >= ydiff)
+         {
+            WALBERLA_CHECK_LESS_EQUAL(fluid[0], boundary[0])
+            return xdiff / (boundary[0] - fluid[0]);
+         }
+         WALBERLA_CHECK_GREATER_EQUAL(fluid[1], boundary[1])
+         return ydiff / (fluid[1] - boundary[1]);
+      }
+
+      WALBERLA_CHECK_LESS_EQUAL(fluid[0], boundary[0])
+      return xdiff / (boundary[0] - fluid[0]);
+   }
+   else if (fluid[0] >= cylinder.xMax())
+   {
+      const real_t xdiff = fluid[0] - cylinder.xMax();
+
+      if (fluid[1] <= cylinder.yMin())
+      {
+         const real_t ydiff = cylinder.yMin() - fluid[1];
+         if (xdiff >= ydiff)
+         {
+            WALBERLA_CHECK_GREATER_EQUAL(fluid[0], boundary[0])
+            return xdiff / (fluid[0] - boundary[0]);
+         }
+         WALBERLA_CHECK_LESS_EQUAL(fluid[1], boundary[1])
+         return ydiff / (boundary[1] - fluid[1]);
+      }
+      else if (fluid[1] >= cylinder.yMax())
+      {
+         const real_t ydiff = fluid[1] - cylinder.yMax();
+         if (xdiff >= ydiff)
+         {
+            WALBERLA_CHECK_GREATER_EQUAL(fluid[0], boundary[0])
+            return xdiff / (fluid[0] - boundary[0]);
+         }
+         WALBERLA_CHECK_GREATER_EQUAL(fluid[1], boundary[1])
+         return ydiff / (fluid[1] - boundary[1]);
+      }
+
+      WALBERLA_CHECK_GREATER_EQUAL(fluid[0], boundary[0])
+      return xdiff / (fluid[0] - boundary[0]);
+   }
+
+   if (fluid[1] <= cylinder.yMin())
+   {
+      WALBERLA_CHECK_LESS_EQUAL(fluid[1], boundary[1])
+      return (cylinder.yMin() - fluid[1]) / (boundary[1] - fluid[1]);
+   }
+
+   WALBERLA_CHECK_GREATER_EQUAL(fluid[1], cylinder.yMax())
+   WALBERLA_CHECK_GREATER_EQUAL(fluid[1], boundary[1])
+   return (fluid[1] - cylinder.yMax()) / (fluid[1] - boundary[1]);
+}
+
+real_t wallDistance::operator()(const Cell& fluidCell, const Cell& boundaryCell,
+                                const shared_ptr< StructuredBlockForest >& SbF, IBlock& block) const
+{
+   Vector3< real_t > boundary = SbF->getBlockLocalCellCenter( block, boundaryCell );
+   Vector3< real_t > fluid = SbF->getBlockLocalCellCenter( block, fluidCell );
+   SbF->mapToPeriodicDomain(boundary);
+   SbF->mapToPeriodicDomain(fluid);
+
+   WALBERLA_ASSERT(cylinder_.contains(boundary), "Boundary cell must be inside the cylinder!")
+   WALBERLA_ASSERT(!cylinder_.contains(fluid), "Fluid cell must not be inside the cylinder!")
+
+   return cylinder_.delta( fluid, boundary );
+}
+} // namespace walberla
\ No newline at end of file

apps/showcases/FlowAroundCylinder/Cylinder.h

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file Cylinder.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "blockforest/SetupBlock.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/StructuredBlockForest.h"
+
+#include "domain_decomposition/IBlock.h"
+
+#include "core/DataTypes.h"
+#include "core/math/AABB.h"
+#include "core/math/Vector3.h"
+#include "core/cell/Cell.h"
+
+#include "Setup.h"
+
+namespace walberla
+{
+
+class Cylinder
+{
+ public:
+   Cylinder(const Setup& setup) : setup_(setup) {}
+
+   bool operator()(const Vector3< real_t >& point) const { return contains(point); }
+
+   bool contains(const Vector3< real_t >& point) const;
+   bool contains(const AABB& aabb) const;
+
+   bool intersects(const AABB& aabb, const real_t bufferDistance = real_t(0)) const;
+
+   real_t delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const;
+
+ private:
+   Setup setup_;
+
+}; // class Cylinder
+
+class CylinderRefinementSelection
+{
+ public:
+   CylinderRefinementSelection(const Cylinder& cylinder, const uint_t level, const real_t bufferDistance)
+      : cylinder_(cylinder), level_(level), bufferDistance_(bufferDistance)
+   {}
+
+   void operator()(SetupBlockForest& forest)
+   {
+      for (auto block = forest.begin(); block != forest.end(); ++block)
+      {
+         const AABB& aabb = block->getAABB();
+
+         if (block->getLevel() < level_ && cylinder_.intersects(aabb, bufferDistance_) && !cylinder_.contains(aabb))
+            block->setMarker(true);
+      }
+   }
+
+ private:
+   Cylinder cylinder_;
+   uint_t level_;
+   real_t bufferDistance_;
+
+}; // class CylinderRefinementSelection
+
+class CylinderBlockExclusion
+{
+ public:
+   CylinderBlockExclusion(const Cylinder& cylinder) : cylinder_(cylinder) {}
+
+   bool operator()(const blockforest::SetupBlock& block)
+   {
+      const AABB aabb = block.getAABB();
+      return static_cast< bool >(cylinder_.contains(aabb));
+   }
+
+ private:
+   Cylinder cylinder_;
+
+}; // class CylinderBlockExclusion
+
+
+class wallDistance
+{
+ public:
+   wallDistance(const Cylinder& cylinder) : cylinder_(cylinder) {}
+
+   real_t operator()(const Cell& fluidCell, const Cell& boundaryCell, const shared_ptr< StructuredBlockForest >& SbF,
+                     IBlock& block) const;
+
+ private:
+   Cylinder cylinder_;
+}; // class wallDistance
+
+} // namespace walberla
\ No newline at end of file

apps/showcases/FlowAroundCylinder/CylinderRe3900.prm

+Parameters
+{
+    kinViscosity       0.005148; // [m^2/s]
+    rho                1;     // [kg/m^3]
+    inflowVelocity     20.0772;  // [m/s]
+    maxLatticeVelocity 0.0337749907475931;
+
+    timesteps 150001;
+
+    diameterCylinder 1.0; // On the finest mesh
+    cylinderXPosition 16;
+    cylinderyPosition 15;
+    coarseMeshSize 0.05859375;
+    circularCrossSection true;
+    inflowProfile 2; // 1: Parabolic (SchaeferTurek), 2: uniform (Jacob et al)
+
+    evaluatePressure false;
+    pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
+    pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
+
+    evaluateStrouhal true;
+    pStrouhal < 30, 23, 3 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
+
+    cylinderRefinementBuffer 0;
+    processMemoryLimit 512; // MiB
+
+    innerOuterSplit <1, 1, 1>;
+
+    // GPU related Parameters, only used if GPU is enabled
+    sendDirectlyFromGPU true;
+    gpuBlockSize <128, 1, 1>;
+
+}
+
+DomainSetup
+{
+    cellsPerBlock < 64, 64, 64 >;
+    domainSize    < 60, 30, 7.5 >;
+    periodic      < false, false, true>;
+    refinementLevels 2;
+    numberProcesses 1; // This is for load balancing, overwritten if more than one proc is used
+}
+
+Boundaries
+{
+    Border { direction N;    walldistance -1;  flag FreeSlip;  }
+    Border { direction S;    walldistance -1;  flag FreeSlip;  }
+    Border { direction W;    walldistance -1;  flag UBB;       }
+    Border { direction E;    walldistance -1;  flag Outflow;   }
+}
+
+StabilityChecker
+{
+    checkFrequency 0;
+    streamOutput   false;
+    vtkOutput      true;
+}
+
+VTKWriter
+{
+    vtkWriteFrequency 25000;
+    velocity true;
+    density false;
+    averageFields true;
+    flag false;
+    writeOnlySlice false;
+    amrFileFormat false;
+    oneFilePerProcess false;
+}
+
+Logging
+{
+    logLevel info;  // info progress detail tracing
+    writeSetupForestAndReturn true;
+    readSetupFromFile false;
+    remainingTimeLoggerFrequency 60; // in seconds
+}
+
+Evaluation
+{
+    evaluationCheckFrequency 2000;
+    rampUpTime 50000;
+    logToStream true;
+    logToFile true;
+    filename CylinderRe3900.txt;
+}
+
+AABBRefinementSelection
+{
+    // coordinates of Regions are always in [0,1] -> the entire simulation space is mapped to [ <0,0,0>, <1,1,1> ]
+    Region
+    {
+        level 2;
+        region [ <0.25, 0.49, 0>, <0.35, 0.51, 1> ];
+    }
+    /*Region
+    {
+        level 2;
+        region [ <0.26, 0.41, 0>, <0.7, 0.59, 1> ];
+    }
+    Region
+    {
+        level 3;
+        region [ <0.26, 0.41, 0>, <0.5, 0.59, 1> ];
+    }
+    Region
+    {
+        level 4;
+        region [ <0.26, 0.41, 0>, <0.4, 0.59, 1> ];
+    }
+    Region
+    {
+        level 5;
+        region [ <0.26, 0.41, 0>, <0.33, 0.59, 1> ];
+    }*/
+}
\ No newline at end of file

apps/showcases/FlowAroundCylinder/Decomp.prm

+Parameters
+{
+    kinViscosity       0.005148; // [m^2/s]
+    rho                1;     // [kg/m^3]
+    inflowVelocity     20.0772;  // [m/s]
+    maxLatticeVelocity 0.0337749907475931;
+
+    timesteps 150001;
+
+    diameterCylinder 25; // On the finest mesh
+    cylinderXPosition 30;
+    cylinderyPosition 30;
+    coarseMeshSize 1;
+    circularCrossSection true;
+    inflowProfile 2; // 1: Parabolic (SchaeferTurek), 2: uniform (Jacob et al)
+
+    evaluatePressure false;
+    pAlpha < 0.45, 0.2, 0.205 >; // points for evaluating
+    pOmega < 0.55, 0.2, 0.205 >; // the pressure difference
+
+    evaluateStrouhal true;
+    pStrouhal < 30, 23, 3 >; // point used for evaluating the vortex shedding frequency and the Strouhal number
+
+    cylinderRefinementBuffer 0;
+    processMemoryLimit 512; // MiB
+
+    innerOuterSplit <1, 1, 1>;
+
+    // GPU related Parameters, only used if GPU is enabled
+    sendDirectlyFromGPU true;
+    gpuBlockSize <128, 1, 1>;
+
+}
+
+DomainSetup
+{
+    cellsPerBlock < 30, 30, 30 >;
+    domainSize    < 90, 60, 30 >;
+    periodic      < false, false, true>;
+    refinementLevels 2;
+    numberProcesses 3; // This is for load balancing, overwritten if more than one proc is used
+}
+
+Boundaries
+{
+    Border { direction N;    walldistance -1;  flag FreeSlip;  }
+    Border { direction S;    walldistance -1;  flag FreeSlip;  }
+    Border { direction W;    walldistance -1;  flag UBB;       }
+    Border { direction E;    walldistance -1;  flag Outflow;   }
+}
+
+StabilityChecker
+{
+    checkFrequency 0;
+    streamOutput   false;
+    vtkOutput      true;
+}
+
+VTKWriter
+{
+    vtkWriteFrequency 25000;
+    velocity true;
+    density false;
+    averageFields true;
+    flag false;
+    writeOnlySlice false;
+    amrFileFormat false;
+    oneFilePerProcess false;
+}
+
+Logging
+{
+    logLevel info;  // info progress detail tracing
+    writeSetupForestAndReturn true;
+    readSetupFromFile false;
+    remainingTimeLoggerFrequency 60; // in seconds
+}
+
+Evaluation
+{
+    evaluationCheckFrequency 2000;
+    rampUpTime 50000;
+    logToStream true;
+    logToFile true;
+    filename CylinderRe3900.txt;
+}
+
+AABBRefinementSelection
+{
+    // coordinates of Regions are always in [0,1] -> the entire simulation space is mapped to [ <0,0,0>, <1,1,1> ]
+    /*Region
+    {
+        level 2;
+        region [ <0.25, 0.49, 0>, <0.35, 0.51, 1> ];
+    }
+    Region
+    {
+        level 2;
+        region [ <0.26, 0.41, 0>, <0.7, 0.59, 1> ];
+    }
+    Region
+    {
+        level 3;
+        region [ <0.26, 0.41, 0>, <0.5, 0.59, 1> ];
+    }
+    Region
+    {
+        level 4;
+        region [ <0.26, 0.41, 0>, <0.4, 0.59, 1> ];
+    }
+    Region
+    {
+        level 5;
+        region [ <0.26, 0.41, 0>, <0.33, 0.59, 1> ];
+    }*/
+}
\ No newline at end of file

apps/showcases/FlowAroundCylinder/FlowAroundCylinder.py

+import sympy as sp
+import numpy as np
+from pystencils import TypedSymbol, Target
+
+from pystencils.field import fields
+from pystencils.simp.subexpression_insertion import insert_constants, insert_aliases
+
+from lbmpy import Stencil, LBStencil, Method, LBMConfig, LBMOptimisation
+from lbmpy.boundaries.boundaryconditions import ExtrapolationOutflow, UBB, QuadraticBounceBack, FreeSlip, NoSlip
+from lbmpy.creationfunctions import create_lb_collision_rule
+
+from pystencils_walberla import CodeGeneration, generate_info_header
+from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
+
+info_header = """
+#pragma once
+const char * infoStencil = "{stencil}";
+const char * infoStreamingPattern = "{streaming_pattern}";
+const char * infoCollisionOperator = "{collision_operator}";
+"""
+
+omega = sp.symbols("omega")
+
+with CodeGeneration() as ctx:
+    dtype = 'float64'
+    pdf_dtype = 'float64'
+
+    stencil = LBStencil(Stencil.D3Q27)
+    q = stencil.Q
+    dim = stencil.D
+
+    streaming_pattern = 'aa'
+
+    pdfs, pdfs_tmp = fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {pdf_dtype}[3D]", layout='fzyx')
+    velocity_field, density_field = fields(f"velocity({dim}), density(1) : {dtype}[{dim}D]", layout='fzyx')
+    macroscopic_fields = {'density': density_field, 'velocity': velocity_field}
+
+    method_enum = Method.CUMULANT
+    lbm_config = LBMConfig(
+        method=method_enum,
+        stencil=stencil,
+        relaxation_rate=omega,
+        compressible=True,
+        galilean_correction=False,
+        fourth_order_correction=1e-4 if method_enum == Method.CUMULANT and stencil.Q == 27 else False,
+        field_name='pdfs',
+        streaming_pattern=streaming_pattern,
+    )
+
+    lbm_opt = LBMOptimisation(cse_global=False, cse_pdfs=False)
+
+    collision_rule = create_lb_collision_rule(lbm_config=lbm_config, lbm_optimisation=lbm_opt)
+    collision_rule = insert_constants(collision_rule)
+    collision_rule = insert_aliases(collision_rule)
+    lb_method = collision_rule.method
+
+    if ctx.gpu:
+        target = Target.GPU
+        openmp = False
+        cpu_vec = None
+        vp = [('int64_t', 'cudaBlockSize0'),
+              ('int64_t', 'cudaBlockSize1'),
+              ('int64_t', 'cudaBlockSize2')]
+
+        sweep_block_size = (TypedSymbol("cudaBlockSize0", np.int64),
+                            TypedSymbol("cudaBlockSize1", np.int64),
+                            TypedSymbol("cudaBlockSize2", np.int64))
+        sweep_params = {'block_size': sweep_block_size}
+
+    else:
+        if ctx.optimize_for_localhost:
+            cpu_vec = {"instruction_set": None}
+        else:
+            cpu_vec = None
+
+        openmp = True if ctx.openmp else False
+
+        target = Target.CPU
+        sweep_params = {}
+        vp = ()
+
+    freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil))
+    no_slip = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
+                                     boundary_object=NoSlip(), field_data_type=pdf_dtype)
+
+    quadratic_bounce_back = QuadraticBounceBack(omega, calculate_force_on_boundary=True)
+    no_slip_interpolated = lbm_boundary_generator(class_name='Obstacle', flag_uid='Obstacle',
+                                                  boundary_object=quadratic_bounce_back, field_data_type=pdf_dtype)
+
+    ubb = lbm_boundary_generator(class_name='UBB', flag_uid='UBB',
+                                 boundary_object=UBB(lambda *args: None, dim=stencil.D, data_type=dtype),
+                                 field_data_type=pdf_dtype)
+
+    outflow = lbm_boundary_generator(class_name='Outflow', flag_uid='Outflow',
+                                     boundary_object=ExtrapolationOutflow(stencil[4], lb_method),
+                                     field_data_type=pdf_dtype)
+
+    generate_lbm_package(ctx, name="FlowAroundCylinder", collision_rule=collision_rule,
+                         lbm_config=lbm_config, lbm_optimisation=lbm_opt,
+                         nonuniform=True, boundaries=[freeslip, no_slip, no_slip_interpolated, ubb, outflow],
+                         macroscopic_fields=macroscopic_fields, set_pre_collision_pdfs=True,
+                         gpu_indexing_params=sweep_params,
+                         target=target, data_type=dtype, pdfs_data_type=pdf_dtype,
+                         cpu_vectorize_info=cpu_vec)
+
+    field_typedefs = {'VelocityField_T': velocity_field,
+                      'ScalarField_T': density_field}
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, 'FlowAroundCylinderInfoHeader',
+                         field_typedefs=field_typedefs)
+
+    infoHeaderParams = {
+        'stencil': stencil.name.lower(),
+        'streaming_pattern': streaming_pattern,
+        'collision_operator': lbm_config.method.name.lower(),
+    }
+
+    ctx.write_file("FlowAroundCylinderStaticDefines.h", info_header.format(**infoHeaderParams))