apps/showcases/FlowAroundSphere/GridGeneration.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 GridGeneration.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "blockforest/Initialization.h"
+#include "blockforest/SetupBlock.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/timing/TimingPool.h"
+
+#include <string>
+#include <fstream>
+
+#include "Types.h"
+#include "Setup.h"
+#include "Sphere.h"
+
+using namespace walberla;
+
+uint_t blockCells(const Setup& setup, const bool withGhostLayers){
+   uint_t cells;
+   if(!withGhostLayers){
+      cells = setup.cellsPerBlock[0] * setup.cellsPerBlock[1] * setup.cellsPerBlock[2];
+   }
+   else{
+      cells = (setup.cellsPerBlock[0] + 2 * setup.numGhostLayers) *
+              (setup.cellsPerBlock[1] + 2 * setup.numGhostLayers) *
+              (setup.cellsPerBlock[2] + 2 * setup.numGhostLayers);
+   }
+   return cells;
+}
+
+
+void createSetupBlockForest(SetupBlockForest& setupBfs, const Setup& setup, const bool useMPIManager=false)
+{
+   WALBERLA_LOG_INFO_ON_ROOT("Generating SetupBlockForest...")
+
+   uint_t numProcesses = setup.numProcesses;
+   const std::string blockForestFilestem = setup.blockForestFilestem;
+
+   if(useMPIManager) {numProcesses = uint_c(mpi::MPIManager::instance()->numProcesses());}
+
+   Sphere Sphere( setup );
+   SphereRefinementSelection SphereRefinementSelection( Sphere, setup.refinementLevels );
+   SphereBlockExclusion SphereBlockExclusion( Sphere );
+
+   setupBfs.addRefinementSelectionFunction(std::function<void(SetupBlockForest &)>(SphereRefinementSelection));
+   setupBfs.addBlockExclusionFunction(SphereBlockExclusion);
+
+   const AABB domain(real_t(0.0), real_t(0.0), real_t(0.0), setup.domainSize[0], setup.domainSize[1], setup.domainSize[2]);
+   setupBfs.addWorkloadMemorySUIDAssignmentFunction(blockforest::uniformWorkloadAndMemoryAssignment);
+   setupBfs.init(domain, setup.rootBlocks[0], setup.rootBlocks[1], setup.rootBlocks[2],
+                         setup.periodic[0], setup.periodic[1], setup.periodic[2]);
+   setupBfs.balanceLoad(blockforest::StaticLevelwiseCurveBalanceWeighted(), numProcesses);
+
+   WALBERLA_LOG_INFO_ON_ROOT("===========================  BLOCK FOREST STATISTICS ============================");
+   WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks())
+   for (uint_t level = 0; level <= setup.refinementLevels; level++){
+      const uint_t numberOfBlocks = setupBfs.getNumberOfBlocks(level);
+      WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << numberOfBlocks)
+   }
+
+   const real_t avgBlocksPerProc = real_c(setupBfs.getNumberOfBlocks()) / real_c(setupBfs.getNumberOfProcesses());
+   WALBERLA_LOG_INFO_ON_ROOT("Average blocks per process: " << avgBlocksPerProc);
+
+   const uint_t cellsPerBlock        = blockCells(setup, false);
+   const uint_t totalNumberCells     = setupBfs.getNumberOfBlocks() * cellsPerBlock;
+   const real_t averageCellsPerGPU   = avgBlocksPerProc * real_c(cellsPerBlock);
+   const uint_t cellsPerBlockGL      = blockCells(setup, true);
+   const uint_t totalNumberCellsGL   = setupBfs.getNumberOfBlocks() * cellsPerBlockGL;
+   const real_t averageCellsPerGPUGL = avgBlocksPerProc * real_c(cellsPerBlockGL);
+
+   const uint_t PDFsPerCell            = StorageSpecification_T::inplace ? Stencil_T::Q : 2 * Stencil_T::Q;
+   const uint_t valuesPerCell          = (PDFsPerCell + VelocityField_T::F_SIZE + ScalarField_T::F_SIZE);
+   const uint_t sizePerValue           = sizeof(StorageSpecification_T::value_type);
+   const double expectedMemory         = double_c(totalNumberCells * valuesPerCell * sizePerValue) * 1e-9;
+   const double expectedMemoryPerGPU   = double_c(averageCellsPerGPU * valuesPerCell * sizePerValue) * 1e-9;
+   const double expectedMemoryGL       = double_c(totalNumberCellsGL * valuesPerCell * sizePerValue) * 1e-9;
+   const double expectedMemoryPerGPUGL = double_c(averageCellsPerGPUGL * valuesPerCell * sizePerValue) * 1e-9;
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Total number of cells will be " << totalNumberCells << " fluid cells (in total on all levels)")
+   WALBERLA_LOG_INFO_ON_ROOT( "Expected total memory demand will be " << expectedMemory << " GB")
+   WALBERLA_LOG_INFO_ON_ROOT( "Average memory demand per GPU will be " << expectedMemoryPerGPU << " GB")
+   WALBERLA_LOG_INFO_ON_ROOT( "Expected total memory demand (with GL) will be " << expectedMemoryGL << " GB")
+   WALBERLA_LOG_INFO_ON_ROOT( "Average memory demand per GPU (with GL) will be " << expectedMemoryPerGPUGL << " GB")
+
+   WALBERLA_LOG_INFO_ON_ROOT("=================================================================================")
+
+   if(mpi::MPIManager::instance()->numProcesses() > 1)
+      return;
+
+   if(setup.writeSetupForestAndReturn){
+      std::ostringstream oss;
+      oss << blockForestFilestem << ".bfs";
+      setupBfs.saveToFile(oss.str().c_str());
+
+      setupBfs.writeVTKOutput(blockForestFilestem);
+   }
+}
+
+void createBlockForest(shared_ptr< BlockForest >& bfs, const Setup& setup){
+   if (mpi::MPIManager::instance()->numProcesses() > 1){
+      // Load structured block forest from file
+      std::ostringstream oss;
+      oss << setup.blockForestFilestem << ".bfs";
+      const std::string setupBlockForestFilepath = oss.str();
+      std::ifstream infile(setupBlockForestFilepath.c_str());
+      if(!infile.good()){
+         WALBERLA_LOG_WARNING_ON_ROOT("Blockforest was not created beforehand and thus needs to be created on the fly. For large simulation runs this can be a severe problem!")
+         SetupBlockForest setupBfs;
+         createSetupBlockForest(setupBfs, setup, true);
+         bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
+      }
+      else{
+         bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()),
+                                               setupBlockForestFilepath.c_str(), false);
+      }
+   }
+   else{
+      SetupBlockForest setupBfs;
+      createSetupBlockForest(setupBfs, setup);
+      bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
+   }
+}
\ No newline at end of file

apps/showcases/FlowAroundSphere/Setup.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 Setup.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "core/config/Config.h"
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+
+#include <string>
+
+namespace walberla
+{
+struct Setup
+{
+
+   Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & domainParameters,
+         const Config::BlockHandle & logging, const Config::BlockHandle & boundaries,
+         std::map<std::string, std::string>& infoMap)
+   {
+      blockForestFilestem         = domainParameters.getParameter< std::string >("blockForestFilestem", "blockforest");
+      numProcesses                = domainParameters.getParameter< uint_t >("numberProcesses");
+      const Vector3< real_t > ds  = domainParameters.getParameter< Vector3< real_t > >("domainSize");
+      const real_t coarseMeshSize = parameters.getParameter< real_t >("coarseMeshSize");
+      cellsPerBlock               = domainParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+
+      rootBlocks[0] = uint_c(std::ceil( (ds[0] / coarseMeshSize) / real_c(cellsPerBlock[0])));
+      rootBlocks[1] = uint_c(std::ceil( (ds[1] / coarseMeshSize) / real_c(cellsPerBlock[1])));
+      rootBlocks[2] = uint_c(std::ceil( (ds[2] / coarseMeshSize) / real_c(cellsPerBlock[2])));
+
+      cells      = Vector3<uint_t>(rootBlocks[0] * cellsPerBlock[0], rootBlocks[1] * cellsPerBlock[1], rootBlocks[2] * cellsPerBlock[2]);
+      domainSize = Vector3<real_t>(cells) * coarseMeshSize;
+
+      periodic = domainParameters.getParameter< Vector3< bool > >("periodic");
+      refinementLevels = domainParameters.getParameter< uint_t >("refinementLevels");
+
+      sphereDiameter = parameters.getParameter< real_t >("diameterSphere") / coarseMeshSize;
+      sphereRadius   = sphereDiameter / real_c(2.0);
+
+      sphereXPosition = parameters.getParameter< real_t >("SphereXPosition") / coarseMeshSize;
+      sphereYPosition = real_c(cells[1]) / real_c(2.0);
+      sphereZPosition = real_c(cells[2]) / real_c(2.0);
+
+      reynoldsNumber    = parameters.getParameter< real_t >("reynoldsNumber");
+      referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
+      inflowVelocity    = parameters.getParameter< real_t >("latticeVelocity");
+
+      const real_t speedOfSound = real_c(real_c(1.0) / std::sqrt( real_c(3.0) ));
+      machNumber = inflowVelocity / speedOfSound;
+      viscosity  = real_c((inflowVelocity * sphereDiameter) / reynoldsNumber);
+      omega      = real_c(real_c(1.0) / (real_c(3.0) * viscosity + real_c(0.5)));
+
+      rho = real_c(1.0);
+      dxC = coarseMeshSize;
+      dxF = real_c(coarseMeshSize) / real_c( 1 << refinementLevels );
+      dt = inflowVelocity / referenceVelocity * coarseMeshSize;
+
+      resolutionSphere = parameters.getParameter< real_t >("diameterSphere") / dxF;
+
+      evaluatePressure = parameters.getParameter< bool >("evaluatePressure");
+      pAlpha = parameters.getParameter< Vector3<real_t> >( "pAlpha" );
+      pOmega = parameters.getParameter< Vector3<real_t> >( "pOmega" );
+
+      evaluateStrouhal = parameters.getParameter< bool >("evaluateStrouhal");;
+      pStrouhal = parameters.getParameter< Vector3<real_t> >( "pStrouhal");
+
+      timesteps = parameters.getParameter< uint_t >("timesteps");
+      numGhostLayers = uint_c(2);
+      nbrOfEvaluationPointsForCoefficientExtremas = 100;
+
+      valuesPerCell = (Stencil_T::Q + VelocityField_T::F_SIZE + uint_c(2) * ScalarField_T::F_SIZE);
+      memoryPerCell = memory_t(valuesPerCell * sizeof(PdfField_T::value_type) + uint_c(1));
+      processMemoryLimit = parameters.getParameter< memory_t >( "processMemoryLimit", memory_t( 512 ) ) * memory_t( 1024 * 1024 );
+
+      stencil          = infoMap["stencil"];
+      streamingPattern = infoMap["streamingPattern"];
+      collisionModel   = infoMap["collisionOperator"];
+
+      fluidUID    = FlagUID("Fluid");
+      obstacleUID = FlagUID(boundaries.getParameter< std::string >("sphere"));
+      inflowUID   = FlagUID(boundaries.getParameter< std::string >("inflow"));
+      outflowUID  = FlagUID(boundaries.getParameter< std::string >("outflow"));
+      wallUID     = FlagUID(boundaries.getParameter< std::string >("walls"));
+
+      writeSetupForestAndReturn = logging.getParameter< bool >("writeSetupForestAndReturn", false);
+
+   }
+
+   void writeParameterHeader(std::ofstream& file)
+   {
+      file << "ReynoldsNumber" << "," << "machNumber" << "," << "coarseMeshSize" << "," << "fineMeshSize" << "," << "resolutionSphere" << ",";
+      file << "refinementLevels" << "," << "stencil" << "," << "streamingPattern" << "," << "collisionOperator" << "," << "omega-coarse" << "\n";
+
+      file << reynoldsNumber << "," << machNumber << "," << dxC << "," << dxF << "," << resolutionSphere << ",";
+      file << refinementLevels << "," << stencil << "," << streamingPattern << "," << collisionModel << "," << omega << "\n";
+   }
+
+   std::string blockForestFilestem;
+   uint_t numProcesses;
+
+   Vector3<uint_t> rootBlocks;
+   Vector3<uint_t> cells;
+   Vector3<real_t> domainSize;
+
+   Vector3< bool > periodic;
+   uint_t refinementLevels;
+
+   Vector3< uint_t > cellsPerBlock;
+   uint_t numGhostLayers;
+
+   real_t sphereXPosition;
+   real_t sphereYPosition;
+   real_t sphereZPosition;
+   real_t sphereRadius;
+   real_t sphereDiameter;
+   real_t resolutionSphere;
+
+   uint_t nbrOfEvaluationPointsForCoefficientExtremas;
+
+   bool evaluatePressure;
+   Vector3< real_t > pAlpha;
+   Vector3< real_t > pOmega;
+
+   bool evaluateStrouhal;
+   Vector3< real_t > pStrouhal;
+
+   real_t machNumber;
+   real_t reynoldsNumber;
+
+   real_t viscosity;
+   real_t omega;
+   real_t rho;
+   real_t inflowVelocity;
+   real_t referenceVelocity;
+   real_t dxC;
+   real_t dxF;
+   real_t dt;
+
+   uint_t timesteps;
+
+   uint_t valuesPerCell;
+   memory_t memoryPerCell;
+   memory_t processMemoryLimit;
+
+   std::string stencil;
+   std::string streamingPattern;
+   std::string collisionModel;
+
+   FlagUID fluidUID;
+   FlagUID obstacleUID;
+   FlagUID inflowUID;
+   FlagUID outflowUID;
+   FlagUID wallUID;
+
+   bool writeSetupForestAndReturn;
+};
+
+}
\ No newline at end of file

apps/showcases/FlowAroundSphere/Sphere.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 Sphere.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "Sphere.h"
+
+namespace walberla
+{
+
+bool Sphere::contains(const Vector3< real_t >& point) const
+{
+   return (point - center_).sqrLength() <= radius2_;
+}
+
+bool Sphere::contains(const AABB& aabb) const
+{
+   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]);
+}
+
+real_t Sphere::delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const
+{
+   WALBERLA_ASSERT(!contains(fluid))
+   WALBERLA_ASSERT(contains(boundary))
+
+   // http://devmag.org.za/2009/04/17/basic-collision-detection-in-2d-part-2/
+   const Vector3< real_t > f = fluid - center_;
+   const Vector3< real_t > d = (boundary - center_) - f;
+
+   const real_t dDotd = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
+   const real_t fDotf = f[0] * f[0] + f[1] * f[1] + f[2] * f[2];
+   const real_t dDotf = d[0] * f[0] + d[1] * f[1] + d[2] * f[2];
+
+   const real_t b = real_c(2.0) * dDotf;
+   const real_t c = fDotf - radius2_;
+
+   const real_t bb4ac = b * b - (real_c(4.0) * dDotd * c);
+   WALBERLA_CHECK_GREATER_EQUAL(bb4ac, real_c(0.0))
+
+   const real_t sqrtbb4ac = std::sqrt(bb4ac);
+   const real_t alpha = std::min((-b + sqrtbb4ac) / (real_c(2.0) * dDotd), (-b - sqrtbb4ac) / (real_c(2.0) * dDotd));
+
+   WALBERLA_CHECK_GREATER_EQUAL(alpha, real_c(0.0))
+   WALBERLA_CHECK_LESS_EQUAL(alpha, real_c(1.0))
+
+   return alpha;
+}
+
+void Sphere::setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID)
+{
+
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+      auto flagField       = bIt->getData< FlagField_T >(flagFieldID);
+      const FlagField_T::flag_t inflowFlag  = flagField->registerFlag(setup_.inflowUID);
+      const FlagField_T::flag_t outflowFlag = flagField->registerFlag(setup_.outflowUID);
+      const FlagField_T::flag_t wallFlag    = flagField->registerFlag(setup_.wallUID);
+
+      const cell_idx_t gls = cell_idx_c(flagField->nrOfGhostLayers()) - cell_idx_c(1);
+
+      CellInterval blockBB(-1, -1, -1,
+                           cell_idx_c(setup_.cellsPerBlock[0]), cell_idx_c(setup_.cellsPerBlock[1]), cell_idx_c(setup_.cellsPerBlock[2]));
+
+      // inflow WEST
+      if(sbfs->atDomainXMinBorder(*bIt)){
+         CellInterval west(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMin(),
+                           blockBB.yMax() + gls, blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(west, inflowFlag, flagField);
+      }
+
+      // outflow EAST
+      if(sbfs->atDomainXMaxBorder(*bIt)){
+         CellInterval east(blockBB.xMax(), blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMax() + gls,
+                           blockBB.yMax() + gls, blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(east, outflowFlag, flagField);
+      }
+
+      // SOUTH
+      if(sbfs->atDomainYMinBorder(*bIt))
+      {
+         CellInterval south(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMax() + gls,
+                            blockBB.yMin(), blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(south, wallFlag, flagField);
+      }
+
+      // NORTH
+      if(sbfs->atDomainYMaxBorder(*bIt)){
+         CellInterval north( blockBB.xMin() - gls, blockBB.yMax(), blockBB.zMin() - gls, blockBB.xMax() + gls,
+                             blockBB.yMax() + gls, blockBB.zMax() + gls );
+         setBoundaryFromCellInterval(north, wallFlag, flagField);
+      }
+
+      // BOTTOM
+      if(sbfs->atDomainZMinBorder(*bIt)){
+         CellInterval bottom(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMax() + gls,
+                             blockBB.yMax() + gls, blockBB.zMin());
+         setBoundaryFromCellInterval(bottom, wallFlag, flagField);
+      }
+
+      // TOP
+      if(sbfs->atDomainZMaxBorder(*bIt)){
+         CellInterval top(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMax(), blockBB.xMax() + gls,
+                          blockBB.yMax() + gls, blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(top, wallFlag, flagField);
+      }
+   }
+
+   checkConsistency(sbfs, flagFieldID);
+   setupSphereBoundary(sbfs, flagFieldID);
+}
+
+void Sphere::setBoundaryFromCellInterval(CellInterval& cells, const FlagField_T::flag_t flag, FlagField_T* flagField)
+{
+
+   for (auto cell = cells.begin(); cell != cells.end(); ++cell){
+      if(flagField->get(cell->x(), cell->y(), cell->z()) == FlagField_T::flag_t(0))
+         flagField->addFlag(cell->x(), cell->y(), cell->z(), flag);
+   }
+}
+
+void Sphere::checkConsistency(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID){
+   const uint_t depth = sbfs->getDepth();
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+      Block& b             = dynamic_cast< Block& >(*bIt);
+      auto flagField       = b.getData< FlagField_T >(flagFieldID);
+      if (sbfs->getLevel(b) < depth){
+         for( auto it = flagField->beginWithGhostLayer(1); it != flagField->end(); ++it ){
+            Vector3< real_t > cellCenter = sbfs->getBlockLocalCellCenter( b, it.cell() );
+            sbfs->mapToPeriodicDomain(cellCenter);
+            WALBERLA_CHECK(!contains(cellCenter), "The sphere must be completely located on the finest level")
+         }
+      }
+   }
+}
+
+void Sphere::setupSphereBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID){
+   const uint_t depth = sbfs->getDepth();
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+      Block& b             = dynamic_cast< Block& >(*bIt);
+      auto flagField       = b.getData< FlagField_T >(flagFieldID);
+      uint8_t obstacleFlag = flagField->registerFlag(setup_.obstacleUID);
+
+      if (sbfs->getLevel(b) == depth){
+         for( auto it = flagField->beginWithGhostLayer(1); it != flagField->end(); ++it ){
+            Vector3< real_t > cellCenter = sbfs->getBlockLocalCellCenter( b, it.cell() );
+            sbfs->mapToPeriodicDomain(cellCenter);
+            if (contains(cellCenter)) { flagField->addFlag(it.x(), it.y(), it.z(), obstacleFlag); }
+         }
+      }
+   }
+}
+
+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_CHECK(!sphere_.contains(fluid), "fluid cell is in contained in sphere (" << fluid[0] << ", " << fluid[1] << ", " << fluid[2] << "). The block local cell is " << fluidCell)
+   WALBERLA_CHECK(sphere_.contains(boundary), "boundary cell is not in contained in sphere (" << boundary[0] << ", " << boundary[1] << ", " << boundary[2] << "). The block local cell is " << boundaryCell)
+
+   return sphere_.delta( fluid, boundary );
+}
+} // namespace walberla
\ No newline at end of file

apps/showcases/FlowAroundSphere/Sphere.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 Sphere.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 "field/FlagUID.h"
+
+#include "core/DataTypes.h"
+#include "core/math/AABB.h"
+#include "core/math/Vector3.h"
+#include "core/cell/Cell.h"
+
+#include "stencil/D3Q7.h"
+#include "stencil/D3Q27.h"
+
+#include "Types.h"
+#include "Setup.h"
+
+namespace walberla
+{
+
+class Sphere
+{
+ public:
+   Sphere(const Setup& setup) : setup_(setup)
+   {
+      const real_t px = setup_.sphereXPosition * setup_.dxC;
+      const real_t py = setup_.sphereYPosition * setup_.dxC;
+      const real_t pz = setup_.sphereZPosition * setup_.dxC;
+
+      center_  = Vector3< real_t >(px, py, pz);
+      radius_  = setup_.sphereRadius * setup_.dxC;
+      radius2_ = radius_ * radius_;
+   }
+
+   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;
+
+   real_t delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const;
+   Setup getSetup(){ return setup_; }
+   void setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+   void checkConsistency(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+   void setupSphereBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+   void setBoundaryFromCellInterval(CellInterval& cells, const FlagField_T::flag_t flag, FlagField_T* flagField);
+
+ private:
+   Setup setup_;
+   Vector3< real_t > center_;
+   real_t radius_;
+   real_t radius2_;
+
+}; // class Sphere
+
+class SphereRefinementSelection
+{
+ public:
+   SphereRefinementSelection(const Sphere& sphere, const uint_t level)
+      : sphere_(sphere), level_(level)
+   {
+      auto setup = sphere_.getSetup();
+      const real_t px = setup.sphereXPosition * setup.dxC;
+      const real_t py = setup.sphereYPosition * setup.dxC;
+      const real_t pz = setup.sphereZPosition * setup.dxC;
+
+      center_  = Vector3< real_t >(px, py, pz);
+      const real_t sphereRadius = setup.sphereRadius * setup.dxC;
+      const real_t bufferDistance = setup.dxF;
+      const real_t d = sphereRadius + bufferDistance;
+      radius2_ = d * d;
+
+      sphereBoundingBox1_ = AABB(center_[0], center_[1] - d, center_[2] - d,
+                                center_[0] + (real_c(2.5) * sphereRadius), center_[1] + d, center_[2] + d);
+
+      sphereBoundingBox2_ = AABB(center_[0], center_[1] - d, center_[2] - d,
+                                 center_[0] + (real_c(5) * sphereRadius), center_[1] + d, center_[2] + d);
+   }
+
+   bool contains(const Vector3< real_t >& point) const
+   {
+      return (point - center_).sqrLength() <= radius2_;
+   }
+
+   bool intersects(const AABB& aabb) const
+   {
+      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]);
+   }
+
+   void operator()(SetupBlockForest& forest)
+   {
+      if(level_ == 0)
+         return;
+      for (auto block = forest.begin(); block != forest.end(); ++block)
+      {
+         const AABB& aabb = block->getAABB();
+
+         if (block->getLevel() < level_ && (intersects(aabb) || sphereBoundingBox1_.intersects(aabb)) )
+            block->setMarker(true);
+
+         if (block->getLevel() < (level_ - 1) && (intersects(aabb) || sphereBoundingBox2_.intersects(aabb)) )
+            block->setMarker(true);
+      }
+   }
+
+ private:
+   Sphere sphere_;
+   Vector3< real_t > center_;
+   uint_t level_;
+   AABB sphereBoundingBox1_;
+   AABB sphereBoundingBox2_;
+   real_t radius2_;
+
+}; // class SphereRefinementSelection
+
+class SphereBlockExclusion
+{
+ public:
+   SphereBlockExclusion(const Sphere& sphere) : sphere_(sphere) {}
+
+   bool operator()(const blockforest::SetupBlock& block)
+   {
+      const AABB aabb = block.getAABB();
+      return static_cast< bool >(sphere_.contains(aabb));
+   }
+
+ private:
+   Sphere sphere_;
+
+}; // class SphereBlockExclusion
+
+
+class wallDistance
+{
+ public:
+   wallDistance(const Sphere& sphere) : sphere_(sphere) {}
+
+   real_t operator()(const Cell& fluidCell, const Cell& boundaryCell, const shared_ptr< StructuredBlockForest >& SbF,
+                     IBlock& block) const;
+
+ private:
+   Sphere sphere_;
+}; // class wallDistance
+
+} // namespace walberla
\ No newline at end of file

apps/showcases/FlowAroundSphere/Types.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 Types.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+# pragma once
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include "lbm_generated/field/PdfField.h"
+
+#include "FlowAroundSphereInfoHeader.h"
+
+using namespace walberla;
+
+using StorageSpecification_T = lbm::FlowAroundSphereStorageSpecification;
+using Stencil_T              = StorageSpecification_T::Stencil;
+using CommunicationStencil_T = StorageSpecification_T::CommunicationStencil;
+
+using PdfField_T           = lbm_generated::PdfField< StorageSpecification_T >;
+using FlagField_T          = FlagField< uint8_t >;
+
+struct IDs
+{
+   BlockDataID pdfField;
+   BlockDataID velocityField;
+   BlockDataID densityField;
+   BlockDataID omegaField;
+   BlockDataID flagField;
+
+   BlockDataID pdfFieldGPU;
+   BlockDataID velocityFieldGPU;
+   BlockDataID densityFieldGPU;
+   BlockDataID omegaFieldGPU;
+};

apps/showcases/FlowAroundSphereCPU/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 FlowAroundSphereCPUCodeGen
+        FILE FlowAroundSphereCPU.py
+        OUT_FILES FlowAroundSphereCPUSweepCollection.h FlowAroundSphereCPUSweepCollection.cpp
+        FlowAroundSphereCPUStorageSpecification.h FlowAroundSphereCPUStorageSpecification.cpp
+        FreeSlip.h FreeSlip.cpp
+        NoSlip.h NoSlip.cpp
+        Obstacle.h Obstacle.cpp
+        UBB.h UBB.cpp
+        Outflow.h Outflow.cpp
+        FixedDensity.h FixedDensity.cpp
+        FlowAroundSphereCPUBoundaryCollection.h
+        FlowAroundSphereInfoHeader.h
+        FlowAroundSphereStaticDefines.h)
+
+waLBerla_add_executable ( NAME FlowAroundSphereCPU
+        FILES FlowAroundSphereCPU.cpp Sphere.cpp Evaluation.cpp GridGeneration.h Types.h
+        DEPENDS blockforest boundary core field lbm_generated stencil timeloop vtk FlowAroundSphereCPUCodeGen )

apps/showcases/FlowAroundSphereCPU/Evaluation.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 Evaluation.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#include "Evaluation.h"
+
+
+namespace walberla
+{
+
+void Evaluation::operator()()
+{
+   if (checkFrequency_ == uint_t(0)) return;
+   ++coarseExecutionCounter_;
+   if (rampUpTime_ > coarseExecutionCounter_) return;
+
+   real_t pressureDifference_L(real_t(0));
+   real_t pressureDifference(real_t(0));
+
+   evaluate(pressureDifference_L, pressureDifference);
+
+   if ((coarseExecutionCounter_ - uint_c(1)) % checkFrequency_ != 0) return;
+
+   WALBERLA_CHECK_NOT_NULLPTR(blocks_)
+
+   // Strouhal number (needs vortex shedding frequency)
+
+   real_t vortexVelocity(real_t(0));
+
+   if (setup_.evaluateStrouhal){
+      auto block = blocks_->getBlock(setup_.pStrouhal);
+      if (block != nullptr){
+         const VelocityField_T* const velocityField = block->template getData< VelocityField_T >(ids_.velocityField);
+         const auto cell                            = blocks_->getBlockLocalCell(*block, setup_.pStrouhal);
+         WALBERLA_ASSERT(velocityField->xyzSize().contains(cell))
+         vortexVelocity += velocityField->get(cell.x(), cell.y(), cell.z(), cell_idx_c(0));
+      }
+      mpi::reduceInplace(vortexVelocity, mpi::SUM);
+   }
+
+   WALBERLA_ROOT_SECTION()
+   {
+      for(auto it = dragResults.begin(); it != dragResults.end(); ++it){
+         coefficients_[0].push_back(it->cDRealArea);
+         coefficients_[1].push_back(it->cLRealArea);
+         coefficients_[2].push_back(it->cDDiscreteArea);
+         coefficients_[3].push_back(it->cLDiscreteArea);
+      }
+
+      if (coefficients_[0].size() > setup_.nbrOfEvaluationPointsForCoefficientExtremas){
+         for (uint_t i = uint_t(0); i < uint_t(4); ++i)
+            coefficients_[i].pop_front();
+      }
+
+      for (uint_t i = uint_t(0); i < uint_t(4); ++i){
+         coefficientExtremas_[i] = std::make_pair(*(coefficients_[i].begin()), *(coefficients_[i].begin()));
+         for (auto v = coefficients_[i].begin(); v != coefficients_[i].end(); ++v){
+            coefficientExtremas_[i].first  = std::min(coefficientExtremas_[i].first, *v);
+            coefficientExtremas_[i].second = std::max(coefficientExtremas_[i].second, *v);
+         }
+      }
+
+      std::ostringstream oss;
+
+      if (logToStream_ and !dragResults.empty()){
+         WALBERLA_LOG_RESULT_ON_ROOT(
+            "force acting on sphere (in dimensionless lattice units of the coarsest grid - evaluated in time step "
+            << coarseExecutionCounter_ - uint_c(1) << "):\n   " << force_ << oss.str()
+            << "\ndrag and lift coefficients (including extrema of last " << (coefficients_[0].size() * checkFrequency_)
+            << " time steps):"
+               "\n   \"real\" area:"
+               "\n      c_D: "
+            << dragResults.back().cDRealArea << " (min = " << coefficientExtremas_[0].first << ", max = " << coefficientExtremas_[0].second
+            << ")"
+            << "\n      c_L: " << dragResults.back().cLRealArea << " (min = " << coefficientExtremas_[1].first
+            << ", max = " << coefficientExtremas_[1].second << ")"
+            << "\n   discrete area:"
+               "\n      c_D: "
+            << dragResults.back().cDDiscreteArea << " (min = " << coefficientExtremas_[2].first
+            << ", max = " << coefficientExtremas_[2].second << ")"
+            << "\n      c_L: " << dragResults.back().cLDiscreteArea << " (min = " << coefficientExtremas_[3].first
+            << ", max = " << coefficientExtremas_[3].second << ")")
+      }
+
+      if (setup_.evaluatePressure && logToStream_){
+         WALBERLA_LOG_RESULT_ON_ROOT("pressure:\n   difference: " << pressureDifference << " Pa ("<< pressureDifference_L << ")")
+      }
+
+      if (setup_.evaluateStrouhal)
+      {
+         // We evaluate the derivative (-> strouhalRising_) in order to find the local minima and maxima of the velocity
+         // over time. If we know the time between a local minimum and a local maximum, we can calculate the frequency.
+         // -> "Smooth noise-robust differentiators"
+         // (http://www.holoborodko.com/pavel/numerical-methods/numerical-derivative/smooth-low-noise-differentiators/)
+
+         if (strouhalVelocities_.size() < uint_t(11)) { strouhalVelocities_.push_back(vortexVelocity); }
+         else{
+            for (uint_t i = uint_t(0); i < uint_t(10); ++i)
+               strouhalVelocities_[i] = strouhalVelocities_[i + 1];
+            strouhalVelocities_[10] = vortexVelocity;
+
+            const real_t f1 = strouhalVelocities_[6] - strouhalVelocities_[4];
+            const real_t f2 = strouhalVelocities_[7] - strouhalVelocities_[3];
+            const real_t f3 = strouhalVelocities_[8] - strouhalVelocities_[2];
+            const real_t f4 = strouhalVelocities_[9] - strouhalVelocities_[1];
+            const real_t f5 = strouhalVelocities_[10] - strouhalVelocities_[0];
+
+            const real_t diff =
+               (real_c(322) * f1 + real_c(256) * f2 + real_c(39) * f3 - real_c(32) * f4 - real_c(11) * f5) /
+               real_c(1536);
+
+            if ((diff > real_t(0)) != strouhalRising_){
+               strouhalRising_ = (diff > real_t(0));
+
+               if (strouhalTimeStep_.size() < uint_t(3)) { strouhalTimeStep_.push_back(coarseExecutionCounter_); }
+               else{
+                  strouhalTimeStep_[0] = strouhalTimeStep_[1];
+                  strouhalTimeStep_[1] = strouhalTimeStep_[2];
+                  strouhalTimeStep_[2] = coarseExecutionCounter_;
+               }
+            }
+         }
+
+         if (strouhalTimeStep_.size() == uint_t(3)){
+            strouhalNumberRealD_     = diameterSphere / (meanVelocity * real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]));
+            strouhalEvaluationExecutionCount_ = coarseExecutionCounter_ - uint_t(1);
+
+            if (logToStream_){
+               WALBERLA_LOG_RESULT_ON_ROOT(
+                  "Strouhal number (evaluated in time step "
+                  << strouhalEvaluationExecutionCount_
+                  << "):"
+                     "\n   D/U (in lattice units): "
+                  << (diameterSphere / meanVelocity) << " (\"real\" D), "
+                  << "\n   T: "
+                  << (real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]) * setup_.dt) << " s ("
+                  << real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0])
+                  << ")"
+                     "\n   f: "
+                  << (real_t(1) / (real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]) * setup_.dt)) << " Hz ("
+                  << (real_t(1) / real_c(strouhalTimeStep_[2] - strouhalTimeStep_[0]))
+                  << ")"
+                     "\n   St (\"real\" D):   "
+                  << strouhalNumberRealD_ << "\n")
+            }
+         }
+      }
+
+      std::ofstream outfile( dragFilename_.c_str(), std::ios_base::app );
+
+      for(auto it = dragResults.begin(); it != dragResults.end(); ++it){
+         outfile << it->timestep << ",";
+         outfile << it->Fx << "," << it->Fy << "," << it->Fz << ",";
+         outfile << it->cDRealArea << ",";
+         outfile << it->cLRealArea << ",";
+         outfile << it->cDDiscreteArea << ",";
+         outfile << it->cLDiscreteArea;
+         outfile << "\n";
+      }
+      outfile.close();
+
+      if (logToFile_ and !dragResults.empty()){
+         std::ofstream file(filename_.c_str(), std::ios_base::app);
+         file << coarseExecutionCounter_ - uint_t(1) << " " << force_[0] << " " << force_[1] << " " << force_[2] << " "
+              << dragResults.back().cDRealArea << " " << dragResults.back().cLRealArea << " " << dragResults.back().cDDiscreteArea << " " << dragResults.back().cLDiscreteArea << " "
+              << pressureDifference_L << " " << pressureDifference << " " << vortexVelocity << " "
+              << strouhalNumberRealD_ << '\n';
+         file.close();
+      }
+      dragResults.clear();
+   }
+
+   // WALBERLA_MPI_WORLD_BARRIER();
+}
+
+void Evaluation::resetForce()
+{
+   if (!initialized_) refresh();
+}
+
+void Evaluation::forceCalculation(const uint_t level)
+{
+   if (rampUpTime_ > coarseExecutionCounter_) return;
+
+   if(level == maxLevel_){
+      for (auto b : finestBlocks_){
+         force_ += Vector3<double>(boundaryCollection_.ObstacleObject->getForce(b));
+      }
+
+      mpi::reduceInplace(force_, mpi::SUM);
+      WALBERLA_ROOT_SECTION(){
+         const double meanU2 = double_c(meanVelocity) * double_c(meanVelocity);
+
+         const double cDRealArea = (double_c(8.0) * force_[0]) / (meanU2 * double_c(surfaceAreaSphere));
+         const double cLRealArea = (double_c(8.0) * force_[1]) / (meanU2 * double_c(surfaceAreaSphere));
+         const double cDDiscreteArea = (double_c(8.0) * force_[0]) / (meanU2 * double_c(AD_));
+         const double cLDiscreteArea = (double_c(8.0) * force_[1]) / (meanU2 * double_c(AL_));
+
+         DragCoefficient DC(fineExecutionCounter_, force_, cDRealArea, cLRealArea, cDDiscreteArea, cLDiscreteArea);
+         dragResults.push_back(DC);
+
+         fineExecutionCounter_++;
+      }
+   }
+
+   force_[0] = double_c(0.0);
+   force_[1] = double_c(0.0);
+   force_[2] = double_c(0.0);
+}
+
+void Evaluation::refresh()
+{
+   WALBERLA_CHECK_NOT_NULLPTR(blocks_)
+   const uint_t finestLevel = blocks_->getDepth();
+
+   real_t AD(real_t(0));
+   real_t AL(real_t(0));
+   for (auto block = blocks_->begin(); block != blocks_->end(); ++block){
+      const uint_t blockLevel = blocks_->getLevel(*block);
+      const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
+
+      auto fluid    = flagField->getFlag(setup_.fluidUID);
+      auto obstacle = flagField->getFlag(setup_.obstacleUID);
+      const real_t area  = real_c(4.0);
+
+      auto xyzSize = flagField->xyzSize();
+      for (auto z = xyzSize.zMin(); z <= xyzSize.zMax(); ++z){
+         for (auto y = xyzSize.yMin(); y <= xyzSize.yMax(); ++y){
+            for (auto x = xyzSize.xMin(); x <= xyzSize.xMax(); ++x){
+               if (flagField->isFlagSet(x, y, z, fluid)){
+                  for (auto it = Stencil_T::beginNoCenter(); it != Stencil_T::end(); ++it){
+                     auto nx = x + cell_idx_c(it.cx());
+                     auto ny = y + cell_idx_c(it.cy());
+                     auto nz = z + cell_idx_c(it.cz());
+
+                     if (flagField->isFlagSet(nx, ny, nz, obstacle)){
+                        WALBERLA_CHECK(blockLevel == finestLevel, "The sphere must be completely located on the finest level")
+                        if (it.cx() == 1 && it.cy() == 0 && it.cz() == 0) { AD += area; }
+                        else if (it.cx() == 0 && it.cz() == 0){
+                           if (it.cy() == 1){
+                              AL += area;
+                           }
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+   mpi::reduceInplace(AD, mpi::SUM);
+   mpi::reduceInplace(AL, mpi::SUM);
+
+   WALBERLA_ROOT_SECTION(){
+      AD_ = AD;
+      AL_ = AL;
+   }
+
+   // Check if points alpha and omega (required for evaluating the pressure difference) are located in fluid cells
+   if (setup_.evaluatePressure){
+      auto block = blocks_->getBlock(setup_.pAlpha);
+      if (block != nullptr){
+         const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
+
+         const auto cell = blocks_->getBlockLocalCell(*block, setup_.pAlpha);
+         WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
+
+         const auto fluid = flagField->getFlag(setup_.fluidUID);
+         if (!flagField->isFlagSet(cell, fluid)){
+            WALBERLA_ABORT("Cell for evaluating pressure difference at point alpha " << setup_.pAlpha << " is not a fluid cell!")
+         }
+      }
+
+      block = blocks_->getBlock(setup_.pOmega);
+      if (block != nullptr){
+         const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
+
+         const auto cell = blocks_->getBlockLocalCell(*block, setup_.pOmega);
+         WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
+
+         const auto fluid = flagField->getFlag(setup_.fluidUID);
+         if (!flagField->isFlagSet(cell, fluid)){
+            WALBERLA_ABORT("Cell for evaluating pressure difference at point omega " << setup_.pOmega << " is not a fluid cell!")
+         }
+      }
+   }
+
+   // Check if point for evaluating the Strouhal number is located inside a fluid cell
+
+   if (setup_.evaluateStrouhal){
+      auto block = blocks_->getBlock(setup_.pStrouhal);
+      if (block != nullptr){
+         const FlagField_T* const flagField = block->template getData< FlagField_T >(ids_.flagField);
+
+         const auto cell = blocks_->getBlockLocalCell(*block, setup_.pStrouhal);
+         WALBERLA_ASSERT(flagField->xyzSize().contains(cell))
+
+         const auto fluid = flagField->getFlag(setup_.fluidUID);
+
+         if (!flagField->isFlagSet(cell, fluid))
+            WALBERLA_ABORT("Cell for evaluating the Strouhal number at point " << setup_.pStrouhal << " is not a fluid cell!")
+      }
+   }
+
+   initialized_ = true;
+}
+
+void Evaluation::evaluate(real_t& pressureDifference_L, real_t& pressureDifference)
+{
+   if (!initialized_) refresh();
+
+   real_t pAlpha(real_t(0));
+   real_t pOmega(real_t(0));
+
+   WALBERLA_CHECK_NOT_NULLPTR(blocks_)
+
+   if (setup_.evaluatePressure){
+      auto block = blocks_->getBlock(setup_.pAlpha);
+      if (block != nullptr){
+         const ScalarField_T* const densityField = block->template getData< ScalarField_T >(ids_.densityField);
+         const auto cell                         = blocks_->getBlockLocalCell(*block, setup_.pAlpha);
+         WALBERLA_ASSERT(densityField->xyzSize().contains(cell))
+         pAlpha += densityField->get(cell) / real_c(3);
+      }
+
+      block = blocks_->getBlock(setup_.pOmega);
+      if (block != nullptr){
+         const ScalarField_T* const densityField = block->template getData< ScalarField_T >(ids_.densityField);
+         const auto cell                         = blocks_->getBlockLocalCell(*block, setup_.pOmega);
+         WALBERLA_ASSERT(densityField->xyzSize().contains(cell))
+         pOmega += densityField->get(cell) / real_c(3);
+      }
+
+      mpi::reduceInplace(pAlpha, mpi::SUM);
+      mpi::reduceInplace(pOmega, mpi::SUM);
+   }
+
+   WALBERLA_ROOT_SECTION(){
+      pressureDifference_L = pAlpha - pOmega;
+      pressureDifference   = (pressureDifference_L * setup_.rho * setup_.dxC * setup_.dxC) / (setup_.dt * setup_.dt);
+   }
+}
+
+}
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/Evaluation.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 Evaluation.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+# pragma once
+
+#include "core/config/Config.h"
+#include "core/math/Sample.h"
+#include "core/math/Constants.h"
+
+#include "lbm_generated/field/PdfField.h"
+
+#include "sqlite/SQLite.h"
+
+#include <algorithm>
+#include <cmath>
+#include <cstdlib>
+#include <cstring>
+#include <functional>
+#include <iostream>
+#include <memory>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include "Types.h"
+#include "Setup.h"
+#include "FlowAroundSphereInfoHeader.h"
+using namespace walberla;
+
+using StorageSpecification_T = lbm::FlowAroundSphereCPUStorageSpecification;
+using Stencil_T              = StorageSpecification_T::Stencil;
+using PdfField_T             = lbm_generated::PdfField< StorageSpecification_T >;
+using FlagField_T            = FlagField< uint8_t >;
+using BoundaryCollection_T = lbm::FlowAroundSphereCPUBoundaryCollection< FlagField_T >;
+
+using VoidFunction = std::function<void ()>;
+
+namespace walberla
+{
+
+struct DragCoefficient {
+   uint_t timestep;
+   double Fx;
+   double Fy;
+   double Fz;
+   double cDRealArea;
+   double cLRealArea;
+   double cDDiscreteArea;
+   double cLDiscreteArea;
+   DragCoefficient(uint_t t, Vector3<double> f, double cdR, double clR, double cdD, double clD) : timestep(t), Fx(f[0]), Fy(f[1]), Fz(f[2]), cDRealArea(cdR), cLRealArea(clR), cDDiscreteArea(cdD), cLDiscreteArea(clD) {}
+};
+
+class Evaluation
+{
+ public:
+   Evaluation(std::shared_ptr< StructuredBlockForest >& blocks, const uint_t checkFrequency, const uint_t rampUpTime,
+              BoundaryCollection_T & boundaryCollection,
+              const IDs& ids, const Setup& setup,
+              const bool logToStream = true, const bool logToFile = true,
+              const std::string& filename = std::string("FlowAroundSphere.txt"))
+      : blocks_(blocks), checkFrequency_(checkFrequency), rampUpTime_(rampUpTime),
+        boundaryCollection_(boundaryCollection), ids_(ids), setup_(setup),
+        logToStream_(logToStream), logToFile_(logToFile), filename_(filename)
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR(blocks)
+      maxLevel_ = blocks->getDepth();
+      blocks->getBlocks(finestBlocks_, maxLevel_);
+
+      coefficients_.resize(uint_c(4));
+      coefficientExtremas_.resize(uint_c(4));
+
+      const real_t factor = setup_.dxC / setup_.dxF;
+      diameterSphere      = real_c(2.0) * (setup_.sphereRadius * factor);
+      surfaceAreaSphere   = math::pi * diameterSphere * diameterSphere;
+      meanVelocity        = setup_.inflowVelocity; // (real_c(4.0) * setup_.inflowVelocity) / real_c(9.0);
+
+      dragFilename_ = std::string("dragSphereRe_") + std::to_string(uint_t(setup_.reynoldsNumber)) + std::string("_meshLevels_") +  std::to_string(uint_t(setup_.refinementLevels + 1)) + std::string(".csv");
+
+      WALBERLA_ROOT_SECTION()
+      {
+         filesystem::path dataFile( dragFilename_.c_str() );
+         if( filesystem::exists( dataFile ) )
+            std::remove( dataFile.string().c_str() );
+
+         std::ofstream outfile( dragFilename_.c_str() );
+         outfile << "SEP=," << "\n";
+
+         setup_.writeParameterHeader(outfile);
+
+         outfile << "timestep," << "Fx," << "Fy," << "Fz," << "cDRealArea," << "cLRealArea," << "cDDiscreteArea," << "cLDiscreteArea" << "\n";
+         outfile.close();
+
+         if (logToFile_)
+         {
+            std::ofstream file(filename_.c_str());
+            file << "# time step [1], force (x) [2], force (y) [3], force (z) [4], "
+                    "cD [5], cL [6], "
+                    "pressure difference (in lattice units) [7], pressure difference (in Pa) [8], vortex velocity (in "
+                    "lattice units) [9], "
+                    "Strouhal number (real D) [10]"
+                 << '\n';
+            if (!setup_.evaluatePressure)
+               file << "# ATTENTION: pressure was not evaluated, pressure difference is set to zero!" << '\n';
+            if (!setup_.evaluateStrouhal)
+               file << "# ATTENTION: vortex velocities were not evaluated, Strouhal number is set to zero!"
+                    << '\n';
+            file.close();
+         }
+      }
+
+      refresh();
+
+      WALBERLA_LOG_INFO_ON_ROOT(
+         "Evaluation initialised:"
+            "\n   + Sphere real area:                      " << surfaceAreaSphere
+         << "\n   + Sphere real diameter:                  " << diameterSphere
+         << "\n   + Sphere discrete area drag coefficient: " << AD_
+         << "\n   + Sphere discrete area lift coefficient: " << AL_
+      )
+
+   }
+
+   void operator()();
+   void forceCalculation(const uint_t level); // for calculating the force
+   void resetForce();
+
+   std::function<void (const uint_t)> forceCalculationFunctor()
+   {
+      return [this](uint_t level) { forceCalculation(level); };
+   }
+
+   std::function<void()> resetForceFunctor()
+   {
+      return [this]() { resetForce(); };
+   }
+
+   void refresh();
+
+ protected:
+   void evaluate(real_t& pressureDifference_L, real_t& pressureDifference);
+
+   bool initialized_{false};
+
+   std::shared_ptr< StructuredBlockForest > blocks_;
+   uint_t maxLevel_;
+   std::vector<Block *> finestBlocks_;
+
+   uint_t coarseExecutionCounter_{ uint_c(0) };
+   uint_t fineExecutionCounter_{ uint_c(0) };
+   uint_t checkFrequency_;
+   uint_t rampUpTime_;
+
+   BoundaryCollection_T & boundaryCollection_;
+
+   IDs ids_;
+   Setup setup_;
+
+   real_t diameterSphere;
+   real_t surfaceAreaSphere;
+   real_t meanVelocity;
+
+   Vector3< double > force_;
+   real_t AD_;
+   real_t AL_;
+   std::vector<DragCoefficient> dragResults;
+
+   std::vector< std::deque< real_t > > coefficients_;
+   std::vector< std::pair< real_t, real_t > > coefficientExtremas_;
+
+   std::vector< real_t > strouhalVelocities_;
+   std::vector< uint_t > strouhalTimeStep_;
+   bool strouhalRising_{false};
+   real_t strouhalNumberRealD_ { real_c(0.0) };
+   uint_t strouhalEvaluationExecutionCount_ { uint_c(0) };
+
+   bool logToStream_;
+   bool logToFile_;
+   std::string filename_;
+   std::string dragFilename_;
+
+   std::map< std::string, double > sqlResults_;
+
+}; // class Evaluation
+
+}
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/FlowAroundSphereCPU.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 FlowAroundSphere.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/AABBRefinementSelection.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/StructuredBlockForest.h"
+#include "blockforest/communication/NonUniformBufferedScheme.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
+#include "core/Abort.h"
+#include "core/DataTypes.h"
+#include "core/SharedFunctor.h"
+#include "core/debug/CheckFunctions.h"
+#include "core/logging/Logging.h"
+#include "core/logging/Initialization.h"
+#include "core/math/Vector3.h"
+#include "core/mpi/Environment.h"
+#include "core/mpi/MPIManager.h"
+#include "core/MemoryUsage.h"
+#include "core/mpi/Reduce.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/timing/TimingPool.h"
+
+#include "field/AddToStorage.h"
+#include "field/CellCounter.h"
+#include "field/FlagField.h"
+#include "field/StabilityChecker.h"
+#include "field/adaptors/AdaptorCreators.h"
+#include "field/iterators/FieldIterator.h"
+#include "field/vtk/VTKWriter.h"
+#include "field/vtk/FlagFieldCellFilter.h"
+
+#include "geometry/InitBoundaryHandling.h"
+
+#include "lbm_generated/communication/NonuniformGeneratedPdfPackInfo.h"
+#include "lbm_generated/communication/UniformGeneratedPdfPackInfo.h"
+#include "lbm_generated/evaluation/PerformanceEvaluation.h"
+#include "lbm_generated/field/AddToStorage.h"
+#include "lbm_generated/field/PdfField.h"
+#include "lbm_generated/refinement/BasicRecursiveTimeStep.h"
+#include "lbm_generated/refinement/RefinementScaling.h"
+
+#include "timeloop/SweepTimeloop.h"
+#include "vtk/VTKOutput.h"
+
+#include <cstdlib>
+#include <iostream>
+#include <memory>
+#include <vector>
+
+#include "Evaluation.h"
+#include "GridGeneration.h"
+#include "Setup.h"
+#include "Sphere.h"
+#include "Types.h"
+
+#include "FlowAroundSphereStaticDefines.h"
+
+using namespace walberla;
+
+using BoundaryCollection_T = lbm::FlowAroundSphereCPUBoundaryCollection< FlagField_T >;
+using SweepCollection_T = lbm::FlowAroundSphereCPUSweepCollection;
+
+
+using PdfField_T = lbm_generated::PdfField< StorageSpecification_T >;
+using blockforest::communication::NonUniformBufferedScheme;
+using blockforest::communication::UniformBufferedScheme;
+
+using lbm_generated::NonuniformGeneratedPdfPackInfo;
+using lbm_generated::UniformGeneratedPdfPackInfo;
+
+int main(int argc, char** argv)
+{
+   mpi::Environment env( argc, argv );
+   mpi::MPIManager::instance()->useWorldComm();
+
+   shared_ptr< Config > config = make_shared< Config >();
+   config->readParameterFile( argv[1] );
+   logging::configureLogging(config);
+
+   ///////////////////////
+   /// PARAMETER INPUT ///
+   ///////////////////////
+
+   // read general simulation parameters
+   auto parameters      = config->getOneBlock("Parameters");
+   auto domainParameters= config->getOneBlock("DomainSetup");
+   auto boundariesConfig= config->getBlock("Boundaries");
+   Setup setup(parameters, domainParameters, boundariesConfig, infoMap);
+
+   auto loggingParameters         = config->getOneBlock("Logging");
+   bool writeSetupForestAndReturn              = loggingParameters.getParameter< bool >("writeSetupForestAndReturn", false);
+   if (uint_c(MPIManager::instance()->numProcesses()) > 1) writeSetupForestAndReturn = false;
+
+   WALBERLA_LOG_INFO_ON_ROOT("Diameter of the Sphere is resolved with " << setup.resolutionSphere << " lattice cells.")
+   Sphere Sphere( setup );
+
+   ///////////////////////////
+   /// CREATE BLOCK FOREST ///
+   ///////////////////////////
+
+   shared_ptr< BlockForest > bfs;
+   createBlockForest(bfs, setup);
+
+   if (writeSetupForestAndReturn && mpi::MPIManager::instance()->numProcesses() == 1)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("BlockForest has been created and writen to file. Returning program")
+
+      const uint_t nBlocks = bfs->getNumberOfBlocks();
+      const uint_t nCells = nBlocks * (setup.cellsPerBlock[0] * setup.cellsPerBlock[1] * setup.cellsPerBlock[2]);
+      const memory_t totalMemory = memory_t(nCells) * setup.memoryPerCell;
+
+      WALBERLA_LOG_INFO_ON_ROOT( "Benchmark run data:"
+                                "\n- simulation parameters:"
+                                "\n   + collision model:     " << infoMap["collisionOperator"] <<
+                                "\n   + stencil:             " << infoMap["stencil"] <<
+                                "\n   + streaming:           " << infoMap["streamingPattern"] <<
+                                "\n   + compressible:        " << ( StorageSpecification_T::compressible ? "yes" : "no" ) <<
+                                "\n   + mesh levels:         " << setup.refinementLevels + uint_c(1) <<
+                                "\n   + resolution:          " << setup.dxC << " - on the coarsest grid" <<
+                                "\n   + resolution:          " << setup.dxF << " - on the finest grid" <<
+                                "\n- domain properties:      "
+                                "\n   + # blocks:            " << nBlocks <<
+                                "\n   + # cells:             " << uint_c(real_c(nCells) / real_c(1e6)) << " [1e6]"
+                                "\n   + # cells sphere D:    " << setup.resolutionSphere <<
+                                "\n   + total memory:        " << totalMemory / 1e9 << " [GB]" <<
+                                "\n- simulation properties:  "
+                                "\n   + sphere pos.(x):      " << setup.sphereXPosition * setup.dxC << " [m]" <<
+                                "\n   + sphere pos.(y):      " << setup.sphereYPosition * setup.dxC << " [m]" <<
+                                "\n   + sphere pos.(z):      " << setup.sphereZPosition * setup.dxC << " [m]" <<
+                                "\n   + sphere radius:       " << setup.sphereRadius * setup.dxC << " [m]" <<
+                                "\n   + kin. viscosity:      " << setup.viscosity * setup.dxC * setup.dxC / setup.dt << " [m^2/s] (on the coarsest grid)" <<
+                                "\n   + viscosity LB:        " << setup.viscosity  << " [dx*dx/dt] (on the coarsest grid)" <<
+                                "\n   + omega:               " << setup.omega << " (on the coarsest grid)" <<
+                                "\n   + rho:                 " << setup.rho << " [kg/m^3]" <<
+                                "\n   + inflow velocity:     " << setup.referenceVelocity << " [m/s]" <<
+                                "\n   + lattice velocity:    " << setup.inflowVelocity << " [dx/dt]" <<
+                                "\n   + Reynolds number:     " << setup.reynoldsNumber <<
+                                "\n   + Mach number:         " << setup.machNumber <<
+                                "\n   + dx (coarsest grid):  " << setup.dxC << " [m]" <<
+                                "\n   + dt (coarsest grid):  " << setup.dt << " [s]"
+                                "\n   + #time steps:         " << setup.timesteps << " (on the coarsest grid, " << ( real_c(1.0) / setup.dt ) << " for 1s of real time)" <<
+                                "\n   + simulation time:     " << ( real_c( setup.timesteps ) * setup.dt ) << " [s]" )
+
+      logging::Logging::printFooterOnStream();
+      return EXIT_SUCCESS;
+   }
+
+   auto blocks = std::make_shared< StructuredBlockForest >(bfs, setup.cellsPerBlock[0], setup.cellsPerBlock[1], setup.cellsPerBlock[2]);
+   blocks->createCellBoundingBoxes();
+
+   ////////////////////////////////////
+   /// CREATE AND INITIALIZE FIELDS ///
+   ////////////////////////////////////
+
+   // create fields
+   const StorageSpecification_T StorageSpec = StorageSpecification_T();
+   IDs ids;
+
+   ids.pdfField = lbm_generated::addPdfFieldToStorage(blocks, "pdfs", StorageSpec, uint_c(2), field::fzyx);
+   ids.velocityField = field::addToStorage< VelocityField_T >(blocks, "vel", real_c(0.0), field::fzyx, uint_c(2));
+   ids.densityField = field::addToStorage< ScalarField_T >(blocks, "density", real_c(1.0), field::fzyx, uint_c(2));
+   ids.omegaField = field::addToStorage< ScalarField_T >(blocks, "omega", real_c(0.0), field::fzyx, uint_c(2));
+   ids.flagField = field::addFlagFieldToStorage< FlagField_T >(blocks, "Boundary Flag Field", uint_c(3));
+
+   auto spongeLayer     = config->getOneBlock("SpongeLayer");
+   const bool deactivateSpongeLayer = spongeLayer.getParameter< bool >("deactivateSpongeLayer");
+   const real_t targetOmega         = spongeLayer.getParameter< real_t >("targetOmega");
+   const real_t spongeStart         = spongeLayer.getParameter< real_t >("spongeStart");
+
+   const real_t omegaDiff = setup.omega - targetOmega;
+   const real_t spongeWidth = real_c(blocks->getDomain().xMax()) - spongeStart;
+   const real_t spongeMid = spongeStart + (spongeWidth / real_c(2.0));
+
+   if(!deactivateSpongeLayer)
+      WALBERLA_LOG_WARNING_ON_ROOT("Using a sponge layer at the Outlet. The sponge layer starts at " << spongeStart << " [m] and targets a relaxation rate of " << targetOmega << " at the outlet" )
+
+   for (auto& block : *blocks)
+   {
+      Block& b = dynamic_cast< Block& >(block);
+      const uint_t level = blocks->getLevel(b);
+
+      auto * omegaField = b.getData< ScalarField_T > ( ids.omegaField );
+      for( auto it = omegaField->beginWithGhostLayer(0); it != omegaField->end(); ++it ){
+         if(deactivateSpongeLayer){
+            omegaField->get(it.cell()) = lbm_generated::relaxationRateScaling(setup.omega, level);
+         }
+         else{
+            Cell globalCell;
+            blocks->transformBlockLocalToGlobalCell(globalCell, block, it.cell());
+            Vector3<real_t> cellCenter = blocks->getCellCenter(globalCell, level);
+            const real_t factor = real_c(0.5) + real_c(0.5) * std::tanh(real_c(2.0) * (cellCenter[0] - spongeMid) / spongeWidth);
+            omegaField->get(it.cell()) = lbm_generated::relaxationRateScaling(setup.omega - (factor * omegaDiff), level);
+         }
+      }
+   }
+
+   WALBERLA_MPI_BARRIER()
+
+   const Cell innerOuterSplit =
+      Cell(parameters.getParameter< Vector3< cell_idx_t > >("innerOuterSplit", Vector3< cell_idx_t >(1, 1, 1)));
+   SweepCollection_T sweepCollection(blocks, ids.pdfField, ids.omegaField, ids.densityField, ids.velocityField, innerOuterSplit);
+
+   WALBERLA_LOG_INFO_ON_ROOT("Setting up communication...")
+   auto nonUniformCommunication = std::make_shared< NonUniformBufferedScheme< CommunicationStencil_T > >(blocks);
+   auto nonUniformPackInfo      = lbm_generated::setupNonuniformPdfCommunication< PdfField_T >(blocks, ids.pdfField);
+   nonUniformCommunication->addPackInfo(nonUniformPackInfo);
+
+   WALBERLA_MPI_BARRIER()
+   WALBERLA_LOG_INFO_ON_ROOT("Setting up communication done")
+
+   /////////////////////////
+   /// BOUNDARY HANDLING ///
+   /////////////////////////
+   WALBERLA_LOG_INFO_ON_ROOT("Start BOUNDARY HANDLING")
+   // create and initialize boundary handling
+   Sphere.setupBoundary(blocks, ids.flagField);
+   geometry::setNonBoundaryCellsToDomain< FlagField_T >(*blocks, ids.flagField, setup.fluidUID, cell_idx_c(0));
+
+   if(parameters.getParameter< bool >("initialiseWithInletVelocity"))
+   {
+      for (auto& block : *blocks)
+      {
+         auto * flagField = block.getData< FlagField_T > ( ids.flagField );
+         auto * velField = block.getData< VelocityField_T > ( ids.velocityField );
+         // auto domainFlag = flagField->getFlag(fluidFlagUID);
+
+         for( auto it = flagField->beginWithGhostLayer(2); it != flagField->end(); ++it )
+         {
+            // if (!isFlagSet(it, domainFlag))
+            //  continue;
+
+            velField->get(it.cell(), 0) = setup.inflowVelocity;
+         }
+      }
+   }
+
+   for (auto& block : *blocks)
+   {
+      sweepCollection.initialise(&block, cell_idx_c(1));
+   }
+
+   std::function< real_t(const Cell&, const Cell&, const shared_ptr< StructuredBlockForest >&, IBlock&) >
+      wallDistanceFunctor = wallDistance(Sphere);
+
+   const real_t omegaFinestLevel = lbm_generated::relaxationRateScaling(setup.omega, setup.refinementLevels);
+   BoundaryCollection_T boundaryCollection(blocks, ids.flagField, ids.pdfField, setup.fluidUID, omegaFinestLevel, setup.inflowVelocity, wallDistanceFunctor);
+
+   WALBERLA_MPI_BARRIER()
+   WALBERLA_LOG_INFO_ON_ROOT("BOUNDARY HANDLING done")
+
+   //////////////////////////////////
+   /// SET UP SWEEPS AND TIMELOOP ///
+   //////////////////////////////////
+   WALBERLA_LOG_INFO_ON_ROOT("Start SWEEPS AND TIMELOOP")
+   // flow evaluation
+   auto EvaluationParameters      = config->getOneBlock("Evaluation");
+   const uint_t evaluationCheckFrequency       = EvaluationParameters.getParameter< uint_t >("evaluationCheckFrequency");
+   const uint_t rampUpTime                     = EvaluationParameters.getParameter< uint_t >("rampUpTime");
+   const bool evaluationLogToStream            = EvaluationParameters.getParameter< bool >("logToStream");
+   const bool evaluationLogToFile              = EvaluationParameters.getParameter< bool >("logToFile");
+   const std::string evaluationFilename        = EvaluationParameters.getParameter< std::string >("filename");
+
+   shared_ptr< Evaluation > evaluation( new Evaluation( blocks, evaluationCheckFrequency, rampUpTime, boundaryCollection,
+                                                        ids, setup, evaluationLogToStream, evaluationLogToFile, evaluationFilename));
+
+   // create time loop
+   SweepTimeloop timeloop(blocks->getBlockStorage(), setup.timesteps);
+   std::shared_ptr< lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > >
+      LBMRefinement;
+
+   LBMRefinement = std::make_shared<
+   lbm_generated::BasicRecursiveTimeStep< PdfField_T, SweepCollection_T, BoundaryCollection_T > >(
+   blocks, ids.pdfField, sweepCollection, boundaryCollection, nonUniformCommunication, nonUniformPackInfo);
+   LBMRefinement->addPostBoundaryHandlingBlockFunction(evaluation->forceCalculationFunctor());
+   LBMRefinement->addRefinementToTimeLoop(timeloop);
+   //////////////////
+   /// VTK OUTPUT ///
+   //////////////////
+   WALBERLA_LOG_INFO_ON_ROOT("SWEEPS AND TIMELOOP done")
+
+   auto VTKWriter           = config->getOneBlock("VTKWriter");
+   const uint_t vtkWriteFrequency       = VTKWriter.getParameter< uint_t >("vtkWriteFrequency", 0);
+   const uint_t vtkGhostLayers          = VTKWriter.getParameter< uint_t >("vtkGhostLayers", 0);
+   const bool amrFileFormat             = VTKWriter.getParameter< bool >("amrFileFormat", false);
+   const bool oneFilePerProcess         = VTKWriter.getParameter< bool >("oneFilePerProcess", false);
+   const real_t samplingResolution      = VTKWriter.getParameter< real_t >("samplingResolution", real_c(-1.0));
+   const uint_t initialWriteCallsToSkip = VTKWriter.getParameter< uint_t >("initialWriteCallsToSkip", uint_c(0.0));
+
+   auto finalDomain = blocks->getDomain();
+   if (vtkWriteFrequency > 0)
+   {
+      auto vtkOutput =
+         vtk::createVTKOutput_BlockData(*blocks, "vtk", vtkWriteFrequency, vtkGhostLayers, false, "vtk_FlowAroundSphere",
+                                        "simulation_step", false, true, true, false, 0, amrFileFormat, oneFilePerProcess);
+
+      vtkOutput->setInitialWriteCallsToSkip(initialWriteCallsToSkip);
+
+      vtkOutput->addBeforeFunction([&]() {
+         for (auto& block : *blocks){
+            sweepCollection.calculateMacroscopicParameters(&block);
+         }
+      });
+
+      vtkOutput->setSamplingResolution(samplingResolution);
+
+      field::FlagFieldCellFilter<FlagField_T> fluidFilter( ids.flagField );
+      fluidFilter.addFlag( setup.obstacleUID );
+      vtkOutput->addCellExclusionFilter(fluidFilter);
+
+
+      if (VTKWriter.getParameter< bool >("writeOnlySlice", true)){
+         const AABB sliceXY(finalDomain.xMin(), finalDomain.yMin(), finalDomain.center()[2] - setup.dxF,
+                            finalDomain.xMax(), finalDomain.yMax(), finalDomain.center()[2] + setup.dxF);
+         vtkOutput->addCellInclusionFilter(vtk::AABBCellFilter(sliceXY));
+
+         if (VTKWriter.getParameter< bool >("writeXZSlice", false)){
+            const AABB sliceXZ(finalDomain.xMin(), finalDomain.center()[1] - setup.dxF, finalDomain.zMin(),
+                               finalDomain.xMax(), finalDomain.center()[1] + setup.dxF, finalDomain.zMax());
+            vtkOutput->addCellInclusionFilter(vtk::AABBCellFilter(sliceXZ));
+         }
+      }
+
+      if (VTKWriter.getParameter< bool >("velocity"))
+      {
+         auto velWriter = make_shared< field::VTKWriter< VelocityField_T, float32 > >(ids.velocityField, "velocity");
+         vtkOutput->addCellDataWriter(velWriter);
+      }
+      if (VTKWriter.getParameter< bool >("density"))
+      {
+         auto densityWriter = make_shared< field::VTKWriter< ScalarField_T, float32 > >(ids.densityField, "density");
+         vtkOutput->addCellDataWriter(densityWriter);
+      }
+      if (VTKWriter.getParameter< bool >("omega"))
+      {
+         auto omegaWriter = make_shared< field::VTKWriter< ScalarField_T, float32 > >(ids.omegaField, "omega");
+         vtkOutput->addCellDataWriter(omegaWriter);
+      }
+      if (VTKWriter.getParameter< bool >("flag"))
+      {
+         auto flagWriter = make_shared< field::VTKWriter< FlagField_T > >(ids.flagField, "flag");
+         vtkOutput->addCellDataWriter(flagWriter);
+      }
+      timeloop.addFuncAfterTimeStep(vtk::writeFiles(vtkOutput), "VTK Output");
+   }
+
+   // log remaining time
+   const real_t remainingTimeLoggerFrequency =
+      loggingParameters.getParameter< real_t >("remainingTimeLoggerFrequency", 3.0); // in seconds
+   if (uint_c(remainingTimeLoggerFrequency) > 0)
+   {
+      timeloop.addFuncAfterTimeStep(
+         timing::RemainingTimeLogger(timeloop.getNrOfTimeSteps(), remainingTimeLoggerFrequency),
+         "remaining time logger");
+   }
+
+   // LBM stability check
+   auto CheckerParameters      = config->getOneBlock("StabilityChecker");
+   const uint_t checkFrequency = CheckerParameters.getParameter< uint_t >("checkFrequency", uint_t(0));
+   if (checkFrequency > 0)
+   {
+      auto checkFunction = [](PdfField_T::value_type value) {  return value < math::abs(PdfField_T::value_type(10)); };
+      timeloop.addFuncAfterTimeStep(
+         makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
+            config, blocks, ids.pdfField, ids.flagField, setup.fluidUID, checkFunction)),
+         "Stability check");
+   }
+
+   timeloop.addFuncBeforeTimeStep( SharedFunctor< Evaluation >(evaluation), "evaluation" );
+
+   WALBERLA_MPI_BARRIER()
+   // WALBERLA_LOG_INFO_ON_ROOT("Execute single timestep to fully complete the preprocessing done")
+
+   //////////////////////
+   /// RUN SIMULATION ///
+   //////////////////////
+   const lbm_generated::PerformanceEvaluation< FlagField_T > performance(blocks, ids.flagField, setup.fluidUID);
+   field::CellCounter< FlagField_T > fluidCells(blocks, ids.flagField, setup.fluidUID);
+   fluidCells();
+
+   WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << blocks->getNumberOfBlocks())
+   for (uint_t level = 0; level <= setup.refinementLevels; level++)
+   {
+      WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << blocks->getNumberOfBlocks(level))
+   }
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Benchmark run data:"
+                             "\n- simulation parameters:"
+                             "\n   + collision model:     " << infoMap["collisionOperator"] <<
+                             "\n   + stencil:             " << infoMap["stencil"] <<
+                             "\n   + streaming:           " << infoMap["streamingPattern"] <<
+                             "\n   + compressible:        " << ( StorageSpecification_T::compressible ? "yes" : "no" ) <<
+                             "\n   + mesh levels:         " << setup.refinementLevels + uint_c(1) <<
+                             "\n   + resolution:          " << setup.dxC << " - on the coarsest grid" <<
+                             "\n   + resolution:          " << setup.dxF << " - on the finest grid" <<
+                             "\n- simulation properties:  "
+                             "\n   + sphere pos.(x):      " << setup.sphereXPosition * setup.dxC << " [m]" <<
+                             "\n   + sphere pos.(y):      " << setup.sphereYPosition * setup.dxC << " [m]" <<
+                             "\n   + sphere pos.(z):      " << setup.sphereZPosition * setup.dxC << " [m]" <<
+                             "\n   + sphere radius:       " << setup.sphereRadius * setup.dxC << " [m]" <<
+                             "\n   + kin. viscosity:      " << setup.viscosity * setup.dxC * setup.dxC / setup.dt << " [m^2/s] (on the coarsest grid)" <<
+                             "\n   + viscosity LB:        " << setup.viscosity  << " [dx*dx/dt] (on the coarsest grid)" <<
+                             "\n   + omega:               " << setup.omega << " (on the coarsest grid)" <<
+                             "\n   + rho:                 " << setup.rho << " [kg/m^3]" <<
+                             "\n   + inflow velocity:     " << setup.referenceVelocity << " [m/s]" <<
+                             "\n   + lattice velocity:    " << setup.inflowVelocity << " [dx/dt]" <<
+                             "\n   + Reynolds number:     " << setup.reynoldsNumber <<
+                             "\n   + Mach number:         " << setup.machNumber <<
+                             "\n   + dt (coarsest grid):  " << setup.dt << " [s]"
+                             "\n   + #time steps:         " << timeloop.getNrOfTimeSteps() << " (on the coarsest grid, " << ( real_t(1) / setup.dt ) << " for 1s of real time)" <<
+                             "\n   + simulation time:     " << ( real_c( timeloop.getNrOfTimeSteps() ) * setup.dt ) << " [s]" )
+
+   WALBERLA_LOG_INFO_ON_ROOT("Starting Simulation")
+   WcTimingPool timeloopTiming;
+   WcTimer simTimer;
+
+   WALBERLA_MPI_BARRIER()
+
+   simTimer.start();
+   timeloop.run(timeloopTiming);
+   WALBERLA_MPI_BARRIER()
+   simTimer.end();
+   WALBERLA_LOG_INFO_ON_ROOT("Simulation finished")
+   real_t time = simTimer.max();
+   WALBERLA_MPI_SECTION() { walberla::mpi::reduceInplace(time, walberla::mpi::MAX); }
+   performance.logResultOnRoot(setup.timesteps, time);
+
+   const auto reducedTimeloopTiming = timeloopTiming.getReduced();
+   WALBERLA_LOG_RESULT_ON_ROOT("Time loop timing:\n" << *reducedTimeloopTiming)
+
+   printResidentMemoryStatistics();
+   return EXIT_SUCCESS;
+}
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/FlowAroundSphereCPU.py

+import sympy as sp
+from pystencils import Target
+
+from pystencils.field import fields
+from pystencils.simp import insert_aliases, insert_constants
+
+from lbmpy import LBStencil, LBMConfig, LBMOptimisation
+from lbmpy.boundaries.boundaryconditions import (ExtrapolationOutflow, UBB, QuadraticBounceBack,
+                                                 FreeSlip, NoSlip, FixedDensity)
+from lbmpy.creationfunctions import create_lb_collision_rule
+from lbmpy.enums import Method, Stencil, SubgridScaleModel
+
+from pystencils_walberla import CodeGeneration, generate_info_header
+from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
+
+info_header = """
+#pragma once
+#include <map>
+std::map<std::string, std::string> infoMap{{{{"stencil", "{stencil}"}},
+                                           {{"streamingPattern", "{streaming_pattern}"}}, 
+                                           {{"collisionOperator", "{collision_operator}"}}}};
+"""
+
+omega = sp.symbols("omega")
+inlet_velocity = sp.symbols("u_x")
+max_threads = 256
+
+with CodeGeneration() as ctx:
+    target = Target.CPU
+    sweep_params = {'block_size': (128, 1, 1)} if ctx.gpu else {}
+
+    dtype = 'float64'
+    pdf_dtype = 'float64'
+
+    stencil = LBStencil(Stencil.D3Q27)
+    q = stencil.Q
+    dim = stencil.D
+
+    streaming_pattern = 'pull'
+
+    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')
+    omega_field = fields(f"rr(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_field.center,
+        compressible=True,
+        # subgrid_scale_model=SubgridScaleModel.QR,
+        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)
+    if lbm_config.method == Method.CUMULANT:
+        collision_rule = insert_constants(collision_rule)
+        collision_rule = insert_aliases(collision_rule)
+    lb_method = collision_rule.method
+
+    freeslip = lbm_boundary_generator("FreeSlip", flag_uid="FreeSlip", boundary_object=FreeSlip(stencil),
+                                      field_data_type=pdf_dtype)
+    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((inlet_velocity, 0.0, 0.0), density=1.0, data_type=dtype),
+                                 field_data_type=pdf_dtype)
+
+    outflow_boundary = ExtrapolationOutflow(stencil[4], lb_method)
+    outflow = lbm_boundary_generator(class_name='Outflow', flag_uid='Outflow',
+                                     boundary_object=ExtrapolationOutflow(stencil[4], lb_method),
+                                     field_data_type=pdf_dtype)
+
+    fixed_density = lbm_boundary_generator(class_name='FixedDensity', flag_uid='FixedDensity',
+                                           boundary_object=FixedDensity(1.0),
+                                           field_data_type=pdf_dtype)
+
+    generate_lbm_package(ctx, name="FlowAroundSphereCPU", collision_rule=collision_rule,
+                         lbm_config=lbm_config, lbm_optimisation=lbm_opt,
+                         nonuniform=True, boundaries=[freeslip, no_slip, no_slip_interpolated,
+                                                      ubb, outflow, fixed_density],
+                         macroscopic_fields=macroscopic_fields, gpu_indexing_params=sweep_params,
+                         target=target, data_type=dtype, pdfs_data_type=pdf_dtype,
+                         max_threads=max_threads)
+
+    field_typedefs = {'VelocityField_T': velocity_field,
+                      'ScalarField_T': density_field}
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, 'FlowAroundSphereInfoHeader',
+                         field_typedefs=field_typedefs)
+
+    infoHeaderParams = {
+        'stencil': stencil.name.lower(),
+        'streaming_pattern': streaming_pattern,
+        'collision_operator': lbm_config.method.name.lower(),
+    }
+
+    ctx.write_file("FlowAroundSphereStaticDefines.h", info_header.format(**infoHeaderParams))

apps/showcases/FlowAroundSphereCPU/GridGeneration.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 GridGeneration.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "blockforest/Initialization.h"
+#include "blockforest/SetupBlock.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
+#include "core/Environment.h"
+#include "core/logging/Initialization.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/timing/TimingPool.h"
+
+#include <string>
+#include <fstream>
+
+#include "Types.h"
+#include "Setup.h"
+#include "Sphere.h"
+
+using namespace walberla;
+
+void createSetupBlockForest(SetupBlockForest& setupBfs, const Setup& setup, const bool useMPIManager=false)
+{
+   WALBERLA_LOG_INFO_ON_ROOT("Generating SetupBlockForest...")
+
+   uint_t numProcesses = setup.numProcesses;
+   const std::string blockForestFilestem = setup.blockForestFilestem;
+
+   if(useMPIManager) {numProcesses = uint_c(mpi::MPIManager::instance()->numProcesses());}
+
+   Sphere Sphere( setup );
+   SphereRefinementSelection SphereRefinementSelection( Sphere, setup.refinementLevels );
+   SphereBlockExclusion SphereBlockExclusion( Sphere );
+
+   setupBfs.addRefinementSelectionFunction(std::function<void(SetupBlockForest &)>(SphereRefinementSelection));
+   setupBfs.addBlockExclusionFunction(SphereBlockExclusion);
+
+   const AABB domain(real_t(0.0), real_t(0.0), real_t(0.0), setup.domainSize[0], setup.domainSize[1], setup.domainSize[2]);
+   setupBfs.addWorkloadMemorySUIDAssignmentFunction(blockforest::uniformWorkloadAndMemoryAssignment);
+   setupBfs.init(domain, setup.rootBlocks[0], setup.rootBlocks[1], setup.rootBlocks[2],
+                         setup.periodic[0], setup.periodic[1], setup.periodic[2]);
+   setupBfs.balanceLoad(blockforest::StaticLevelwiseCurveBalanceWeighted(), numProcesses);
+
+   if(mpi::MPIManager::instance()->numProcesses() > 1)
+      return;
+
+   {
+      std::ostringstream oss;
+      oss << blockForestFilestem << ".bfs";
+      setupBfs.saveToFile(oss.str().c_str());
+   }
+
+   setupBfs.writeVTKOutput(blockForestFilestem);
+
+   WALBERLA_LOG_INFO_ON_ROOT("===========================  BLOCK FOREST STATISTICS ============================");
+   WALBERLA_LOG_INFO_ON_ROOT("Blocks created: " << setupBfs.getNumberOfBlocks())
+   for (uint_t level = 0; level <= setup.refinementLevels; level++){
+      const uint_t numberOfBlocks = setupBfs.getNumberOfBlocks(level);
+      WALBERLA_LOG_INFO_ON_ROOT("Level " << level << " Blocks: " << numberOfBlocks)
+   }
+
+   const real_t avgBlocksPerProc = real_c(setupBfs.getNumberOfBlocks()) / real_c(setupBfs.getNumberOfProcesses());
+   WALBERLA_LOG_INFO_ON_ROOT("Average blocks per process: " << avgBlocksPerProc);
+
+   const uint_t totalNumberCellsPerBlock = setup.cellsPerBlock[0] * setup.cellsPerBlock[1] * setup.cellsPerBlock[2];
+   const uint_t totalNumberCells         = setupBfs.getNumberOfBlocks() * totalNumberCellsPerBlock;
+   const real_t averageCellsPerGPU       = avgBlocksPerProc * real_c(totalNumberCellsPerBlock);
+
+   const uint_t PDFsPerCell = StorageSpecification_T::inplace ? Stencil_T::Q : 2 * Stencil_T::Q;
+   const uint_t valuesPerCell = (PDFsPerCell + VelocityField_T::F_SIZE + ScalarField_T::F_SIZE);
+   const uint_t sizePerValue = sizeof(StorageSpecification_T::value_type);
+   const double expectedMemory = double_c(totalNumberCells * valuesPerCell * sizePerValue) * 1e-9;
+   const double expectedMemoryPerGPU = double_c(averageCellsPerGPU * valuesPerCell * sizePerValue) * 1e-9;
+
+   WALBERLA_LOG_INFO_ON_ROOT( "Total number of cells will be " << totalNumberCells << " fluid cells (in total on all levels)")
+   WALBERLA_LOG_INFO_ON_ROOT( "Expected total memory demand will be " << expectedMemory << " GB")
+   WALBERLA_LOG_INFO_ON_ROOT( "Average memory demand per GPU will be " << expectedMemoryPerGPU << " GB")
+
+   WALBERLA_LOG_INFO_ON_ROOT("=================================================================================");
+}
+
+void createBlockForest(shared_ptr< BlockForest >& bfs, const Setup& setup)
+{
+   if (mpi::MPIManager::instance()->numProcesses() > 1)
+   {
+      // Load structured block forest from file
+      std::ostringstream oss;
+      oss << setup.blockForestFilestem << ".bfs";
+      const std::string setupBlockForestFilepath = oss.str();
+      std::ifstream infile(setupBlockForestFilepath.c_str());
+      if(!infile.good()){
+         WALBERLA_LOG_WARNING_ON_ROOT("Blockforest was not created beforehand and thus needs to be created on the fly. For large simulation runs this can be a severe problem!")
+         SetupBlockForest setupBfs;
+         createSetupBlockForest(setupBfs, setup, true);
+         bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
+      }
+      else{
+         bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()),
+                                               setupBlockForestFilepath.c_str(), false);
+      }
+   }
+   else{
+      SetupBlockForest setupBfs;
+      createSetupBlockForest(setupBfs, setup);
+      bfs = std::make_shared< BlockForest >(uint_c(MPIManager::instance()->worldRank()), setupBfs);
+   }
+}
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/Setup.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 Setup.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+#pragma once
+
+#include "core/config/Config.h"
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+
+#include <string>
+
+namespace walberla
+{
+struct Setup
+{
+
+   Setup(const Config::BlockHandle & parameters, const Config::BlockHandle & domainParameters,
+         const Config::BlockHandle & boundaries, std::map<std::string, std::string>& infoMap)
+   {
+      blockForestFilestem         = domainParameters.getParameter< std::string >("blockForestFilestem", "blockforest");
+      numProcesses                = domainParameters.getParameter< uint_t >("numberProcesses");
+      const Vector3< real_t > ds  = domainParameters.getParameter< Vector3< real_t > >("domainSize");
+      const real_t coarseMeshSize = parameters.getParameter< real_t >("coarseMeshSize");
+      cellsPerBlock               = domainParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+
+      rootBlocks[0] = uint_c(std::ceil( (ds[0] / coarseMeshSize) / real_c(cellsPerBlock[0])));
+      rootBlocks[1] = uint_c(std::ceil( (ds[1] / coarseMeshSize) / real_c(cellsPerBlock[1])));
+      rootBlocks[2] = uint_c(std::ceil( (ds[2] / coarseMeshSize) / real_c(cellsPerBlock[2])));
+
+      cells      = Vector3<uint_t>(rootBlocks[0] * cellsPerBlock[0], rootBlocks[1] * cellsPerBlock[1], rootBlocks[2] * cellsPerBlock[2]);
+      domainSize = Vector3<real_t>(cells) * coarseMeshSize;
+
+      periodic = domainParameters.getParameter< Vector3< bool > >("periodic");
+      refinementLevels = domainParameters.getParameter< uint_t >("refinementLevels");
+
+      sphereDiameter = parameters.getParameter< real_t >("diameterSphere") / coarseMeshSize;
+      sphereRadius   = sphereDiameter / real_c(2.0);
+
+      sphereXPosition = parameters.getParameter< real_t >("SphereXPosition") / coarseMeshSize;
+      sphereYPosition = real_c(cells[1]) / real_c(2.0);
+      sphereZPosition = real_c(cells[2]) / real_c(2.0);
+
+      reynoldsNumber    = parameters.getParameter< real_t >("reynoldsNumber");
+      referenceVelocity = parameters.getParameter< real_t >("referenceVelocity");
+      inflowVelocity    = parameters.getParameter< real_t >("latticeVelocity");
+
+      const real_t speedOfSound = real_c(real_c(1.0) / std::sqrt( real_c(3.0) ));
+      machNumber = inflowVelocity / speedOfSound;
+      viscosity  = real_c((inflowVelocity * sphereDiameter) / reynoldsNumber);
+      omega      = real_c(real_c(1.0) / (real_c(3.0) * viscosity + real_c(0.5)));
+
+      rho = real_c(1.0);
+      dxC = coarseMeshSize;
+      dxF = real_c(coarseMeshSize) / real_c( 1 << refinementLevels );
+      dt = inflowVelocity / referenceVelocity * coarseMeshSize;
+
+      resolutionSphere = parameters.getParameter< real_t >("diameterSphere") / dxF;
+
+      evaluatePressure = parameters.getParameter< bool >("evaluatePressure");
+      pAlpha = parameters.getParameter< Vector3<real_t> >( "pAlpha" );
+      pOmega = parameters.getParameter< Vector3<real_t> >( "pOmega" );
+
+      evaluateStrouhal = parameters.getParameter< bool >("evaluateStrouhal");;
+      pStrouhal = parameters.getParameter< Vector3<real_t> >( "pStrouhal");
+
+      timesteps = parameters.getParameter< uint_t >("timesteps");
+      numGhostLayers = uint_c(2);
+      nbrOfEvaluationPointsForCoefficientExtremas = 100;
+
+      valuesPerCell = (Stencil_T::Q + VelocityField_T::F_SIZE + uint_c(2) * ScalarField_T::F_SIZE);
+      memoryPerCell = memory_t(valuesPerCell * sizeof(PdfField_T::value_type) + uint_c(1));
+      processMemoryLimit = parameters.getParameter< memory_t >( "processMemoryLimit", memory_t( 512 ) ) * memory_t( 1024 * 1024 );
+
+      stencil          = infoMap["stencil"];
+      streamingPattern = infoMap["streamingPattern"];
+      collisionModel   = infoMap["collisionOperator"];
+
+      fluidUID    = FlagUID("Fluid");
+      obstacleUID = FlagUID(boundaries.getParameter< std::string >("sphere"));
+      inflowUID   = FlagUID(boundaries.getParameter< std::string >("inflow"));
+      outflowUID  = FlagUID(boundaries.getParameter< std::string >("outflow"));
+      wallUID     = FlagUID(boundaries.getParameter< std::string >("walls"));
+
+   }
+
+   void writeParameterHeader(std::ofstream& file)
+   {
+      file << "ReynoldsNumber" << "," << "machNumber" << "," << "coarseMeshSize" << "," << "fineMeshSize" << "," << "resolutionSphere" << ",";
+      file << "refinementLevels" << "," << "stencil" << "," << "streamingPattern" << "," << "collisionOperator" << "," << "omega-coarse" << "\n";
+
+      file << reynoldsNumber << "," << machNumber << "," << dxC << "," << dxF << "," << resolutionSphere << ",";
+      file << refinementLevels << "," << stencil << "," << streamingPattern << "," << collisionModel << "," << omega << "\n";
+   }
+
+   std::string blockForestFilestem;
+   uint_t numProcesses;
+
+   Vector3<uint_t> rootBlocks;
+   Vector3<uint_t> cells;
+   Vector3<real_t> domainSize;
+
+   Vector3< bool > periodic;
+   uint_t refinementLevels;
+
+   Vector3< uint_t > cellsPerBlock;
+   uint_t numGhostLayers;
+
+   real_t sphereXPosition;
+   real_t sphereYPosition;
+   real_t sphereZPosition;
+   real_t sphereRadius;
+   real_t sphereDiameter;
+   real_t resolutionSphere;
+
+   uint_t nbrOfEvaluationPointsForCoefficientExtremas;
+
+   bool evaluatePressure;
+   Vector3< real_t > pAlpha;
+   Vector3< real_t > pOmega;
+
+   bool evaluateStrouhal;
+   Vector3< real_t > pStrouhal;
+
+   real_t machNumber;
+   real_t reynoldsNumber;
+
+   real_t viscosity;
+   real_t omega;
+   real_t rho;
+   real_t inflowVelocity;
+   real_t referenceVelocity;
+   real_t dxC;
+   real_t dxF;
+   real_t dt;
+
+   uint_t timesteps;
+
+   uint_t valuesPerCell;
+   memory_t memoryPerCell;
+   memory_t processMemoryLimit;
+
+   std::string stencil;
+   std::string streamingPattern;
+   std::string collisionModel;
+
+   FlagUID fluidUID;
+   FlagUID obstacleUID;
+   FlagUID inflowUID;
+   FlagUID outflowUID;
+   FlagUID wallUID;
+};
+
+}
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/Sphere.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 Sphere.cpp
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+
+#include "Sphere.h"
+
+namespace walberla
+{
+
+bool Sphere::contains(const Vector3< real_t >& point) const
+{
+   return (point - center_).sqrLength() <= radius2_;
+}
+
+bool Sphere::contains(const AABB& aabb) const
+{
+   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]);
+}
+
+real_t Sphere::delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const
+{
+   WALBERLA_ASSERT(!contains(fluid))
+   WALBERLA_ASSERT(contains(boundary))
+
+   // http://devmag.org.za/2009/04/17/basic-collision-detection-in-2d-part-2/
+   const Vector3< real_t > f = fluid - center_;
+   const Vector3< real_t > d = (boundary - center_) - f;
+
+   const real_t dDotd = d[0] * d[0] + d[1] * d[1] + d[2] * d[2];
+   const real_t fDotf = f[0] * f[0] + f[1] * f[1] + f[2] * f[2];
+   const real_t dDotf = d[0] * f[0] + d[1] * f[1] + d[2] * f[2];
+
+   const real_t b = real_c(2.0) * dDotf;
+   const real_t c = fDotf - radius2_;
+
+   const real_t bb4ac = b * b - (real_c(4.0) * dDotd * c);
+   WALBERLA_CHECK_GREATER_EQUAL(bb4ac, real_c(0.0))
+
+   const real_t sqrtbb4ac = std::sqrt(bb4ac);
+   const real_t alpha = std::min((-b + sqrtbb4ac) / (real_c(2.0) * dDotd), (-b - sqrtbb4ac) / (real_c(2.0) * dDotd));
+
+   WALBERLA_CHECK_GREATER_EQUAL(alpha, real_c(0.0))
+   WALBERLA_CHECK_LESS_EQUAL(alpha, real_c(1.0))
+
+   return alpha;
+}
+
+void Sphere::setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID)
+{
+
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+      auto flagField       = bIt->getData< FlagField_T >(flagFieldID);
+      const FlagField_T::flag_t inflowFlag  = flagField->registerFlag(setup_.inflowUID);
+      const FlagField_T::flag_t outflowFlag = flagField->registerFlag(setup_.outflowUID);
+      const FlagField_T::flag_t wallFlag    = flagField->registerFlag(setup_.wallUID);
+
+      const cell_idx_t gls = cell_idx_c(flagField->nrOfGhostLayers()) - cell_idx_c(1);
+
+      CellInterval blockBB(-1, -1, -1,
+                           cell_idx_c(setup_.cellsPerBlock[0]), cell_idx_c(setup_.cellsPerBlock[1]), cell_idx_c(setup_.cellsPerBlock[2]));
+
+      // inflow WEST
+      if(sbfs->atDomainXMinBorder(*bIt)){
+         CellInterval west(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMin(),
+                           blockBB.yMax() + gls, blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(west, inflowFlag, flagField);
+      }
+
+      // outflow EAST
+      if(sbfs->atDomainXMaxBorder(*bIt)){
+         CellInterval east(blockBB.xMax(), blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMax() + gls,
+                           blockBB.yMax() + gls, blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(east, outflowFlag, flagField);
+      }
+
+      // SOUTH
+      if(sbfs->atDomainYMinBorder(*bIt))
+      {
+         CellInterval south(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMax() + gls,
+                            blockBB.yMin(), blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(south, wallFlag, flagField);
+      }
+
+      // NORTH
+      if(sbfs->atDomainYMaxBorder(*bIt)){
+         CellInterval north( blockBB.xMin() - gls, blockBB.yMax(), blockBB.zMin() - gls, blockBB.xMax() + gls,
+                             blockBB.yMax() + gls, blockBB.zMax() + gls );
+         setBoundaryFromCellInterval(north, wallFlag, flagField);
+      }
+
+      // BOTTOM
+      if(sbfs->atDomainZMinBorder(*bIt)){
+         CellInterval bottom(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMin() - gls, blockBB.xMax() + gls,
+                             blockBB.yMax() + gls, blockBB.zMin());
+         setBoundaryFromCellInterval(bottom, wallFlag, flagField);
+      }
+
+      // TOP
+      if(sbfs->atDomainZMaxBorder(*bIt)){
+         CellInterval top(blockBB.xMin() - gls, blockBB.yMin() - gls, blockBB.zMax(), blockBB.xMax() + gls,
+                          blockBB.yMax() + gls, blockBB.zMax() + gls);
+         setBoundaryFromCellInterval(top, wallFlag, flagField);
+      }
+   }
+
+   checkConsistency(sbfs, flagFieldID);
+   setupSphereBoundary(sbfs, flagFieldID);
+}
+
+void Sphere::setBoundaryFromCellInterval(CellInterval& cells, const FlagField_T::flag_t flag, FlagField_T* flagField)
+{
+
+   for (auto cell = cells.begin(); cell != cells.end(); ++cell){
+      if(flagField->get(cell->x(), cell->y(), cell->z()) == FlagField_T::flag_t(0))
+         flagField->addFlag(cell->x(), cell->y(), cell->z(), flag);
+   }
+}
+
+void Sphere::checkConsistency(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID){
+   const uint_t depth = sbfs->getDepth();
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+      Block& b             = dynamic_cast< Block& >(*bIt);
+      auto flagField       = b.getData< FlagField_T >(flagFieldID);
+      if (sbfs->getLevel(b) < depth){
+         for( auto it = flagField->beginWithGhostLayer(1); it != flagField->end(); ++it ){
+            Vector3< real_t > cellCenter = sbfs->getBlockLocalCellCenter( b, it.cell() );
+            sbfs->mapToPeriodicDomain(cellCenter);
+            WALBERLA_CHECK(!contains(cellCenter), "The sphere must be completely located on the finest level")
+         }
+      }
+   }
+}
+
+void Sphere::setupSphereBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID){
+   const uint_t depth = sbfs->getDepth();
+   for (auto bIt = sbfs->begin(); bIt != sbfs->end(); ++bIt)
+   {
+      Block& b             = dynamic_cast< Block& >(*bIt);
+      auto flagField       = b.getData< FlagField_T >(flagFieldID);
+      uint8_t obstacleFlag = flagField->registerFlag(setup_.obstacleUID);
+
+      if (sbfs->getLevel(b) == depth){
+         for( auto it = flagField->beginWithGhostLayer(1); it != flagField->end(); ++it ){
+            Vector3< real_t > cellCenter = sbfs->getBlockLocalCellCenter( b, it.cell() );
+            sbfs->mapToPeriodicDomain(cellCenter);
+            if (contains(cellCenter)) { flagField->addFlag(it.x(), it.y(), it.z(), obstacleFlag); }
+         }
+      }
+   }
+}
+
+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_CHECK(!sphere_.contains(fluid), "fluid cell is in contained in sphere (" << fluid[0] << ", " << fluid[1] << ", " << fluid[2] << "). The block local cell is " << fluidCell)
+   WALBERLA_CHECK(sphere_.contains(boundary), "boundary cell is not in contained in sphere (" << boundary[0] << ", " << boundary[1] << ", " << boundary[2] << "). The block local cell is " << boundaryCell)
+
+   return sphere_.delta( fluid, boundary );
+}
+} // namespace walberla
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/Sphere.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 Sphere.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 "field/FlagUID.h"
+
+#include "core/DataTypes.h"
+#include "core/math/AABB.h"
+#include "core/math/Vector3.h"
+#include "core/cell/Cell.h"
+
+#include "stencil/D3Q7.h"
+#include "stencil/D3Q27.h"
+
+#include "Types.h"
+#include "Setup.h"
+
+namespace walberla
+{
+
+class Sphere
+{
+ public:
+   Sphere(const Setup& setup) : setup_(setup)
+   {
+      const real_t px = setup_.sphereXPosition * setup_.dxC;
+      const real_t py = setup_.sphereYPosition * setup_.dxC;
+      const real_t pz = setup_.sphereZPosition * setup_.dxC;
+
+      center_  = Vector3< real_t >(px, py, pz);
+      radius_  = setup_.sphereRadius * setup_.dxC;
+      radius2_ = radius_ * radius_;
+   }
+
+   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;
+
+   real_t delta(const Vector3< real_t >& fluid, const Vector3< real_t >& boundary) const;
+   Setup getSetup(){ return setup_; }
+   void setupBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+   void checkConsistency(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+   void setupSphereBoundary(const std::shared_ptr< StructuredBlockForest >& sbfs, const BlockDataID flagFieldID);
+   void setBoundaryFromCellInterval(CellInterval& cells, const FlagField_T::flag_t flag, FlagField_T* flagField);
+
+ private:
+   Setup setup_;
+   Vector3< real_t > center_;
+   real_t radius_;
+   real_t radius2_;
+
+}; // class Sphere
+
+class SphereRefinementSelection
+{
+ public:
+   SphereRefinementSelection(const Sphere& sphere, const uint_t level)
+      : sphere_(sphere), level_(level)
+   {
+      auto setup = sphere_.getSetup();
+      const real_t px = setup.sphereXPosition * setup.dxC;
+      const real_t py = setup.sphereYPosition * setup.dxC;
+      const real_t pz = setup.sphereZPosition * setup.dxC;
+
+      auto center  = Vector3< real_t >(px, py, pz);
+      const real_t sphereRadius = setup.sphereRadius * setup.dxC;
+      const real_t bufferDistance = real_c(5.0) * setup.dxF;
+      const real_t d = sphereRadius + bufferDistance;
+
+      sphereBoundingBox_ = AABB(center[0] - d, center[1] - d, center[2] - d,
+                          center[0] + (real_c(2.5) * sphereRadius), center[1] + d, center[2] + d);
+   }
+
+   void operator()(SetupBlockForest& forest)
+   {
+      for (auto block = forest.begin(); block != forest.end(); ++block)
+      {
+         const AABB& aabb = block->getAABB();
+
+         if (block->getLevel() < level_ && sphereBoundingBox_.intersects(aabb))
+            block->setMarker(true);
+      }
+   }
+
+ private:
+   Sphere sphere_;
+   uint_t level_;
+   AABB sphereBoundingBox_;
+
+}; // class SphereRefinementSelection
+
+class SphereBlockExclusion
+{
+ public:
+   SphereBlockExclusion(const Sphere& sphere) : sphere_(sphere) {}
+
+   bool operator()(const blockforest::SetupBlock& block)
+   {
+      const AABB aabb = block.getAABB();
+      return static_cast< bool >(sphere_.contains(aabb));
+   }
+
+ private:
+   Sphere sphere_;
+
+}; // class SphereBlockExclusion
+
+
+class wallDistance
+{
+ public:
+   wallDistance(const Sphere& sphere) : sphere_(sphere) {}
+
+   real_t operator()(const Cell& fluidCell, const Cell& boundaryCell, const shared_ptr< StructuredBlockForest >& SbF,
+                     IBlock& block) const;
+
+ private:
+   Sphere sphere_;
+}; // class wallDistance
+
+} // namespace walberla
\ No newline at end of file

apps/showcases/FlowAroundSphereCPU/Types.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 Types.h
+//! \author Markus Holzer <markus.holzer@fau.de>
+//
+//======================================================================================================================
+# pragma once
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include "lbm_generated/field/PdfField.h"
+
+#include "FlowAroundSphereInfoHeader.h"
+
+using namespace walberla;
+
+using StorageSpecification_T = lbm::FlowAroundSphereCPUStorageSpecification;
+using Stencil_T              = StorageSpecification_T::Stencil;
+using CommunicationStencil_T = StorageSpecification_T::CommunicationStencil;
+
+using PdfField_T           = lbm_generated::PdfField< StorageSpecification_T >;
+using FlagField_T          = FlagField< uint8_t >;
+
+struct IDs
+{
+   BlockDataID pdfField;
+   BlockDataID velocityField;
+   BlockDataID densityField;
+   BlockDataID omegaField;
+   BlockDataID flagField;
+
+   BlockDataID pdfFieldGPU;
+   BlockDataID velocityFieldGPU;
+   BlockDataID densityFieldGPU;
+   BlockDataID omegaFieldGPU;
+};

apps/showcases/PhaseFieldAllenCahn/GPU/multiphase.cpp

@@ -34,6 +34,7 @@
 #include "field/AddToStorage.h"
 #include "field/FlagField.h"
 #include "field/communication/PackInfo.h"
+#include "field/vtk/FlagFieldCellFilter.h"
 #include "field/vtk/VTKWriter.h"
 
 #include "geometry/InitBoundaryHandling.h"
@@ -97,14 +98,11 @@ int main(int argc, char** argv)
       //////////////////////////////////////
 
       auto parameters                    = config->getOneBlock("Parameters");
-      const std::string timeStepStrategy = parameters.getParameter< std::string >("timeStepStrategy", "normal");
-      const uint_t timesteps             = parameters.getParameter< uint_t >("timesteps", uint_c(50));
-      const real_t remainingTimeLoggerFrequency =
-         parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0));
-      const uint_t scenario = parameters.getParameter< uint_t >("scenario", uint_c(1));
-      Vector3< int > gpuBlockSize =
-         parameters.getParameter< Vector3< int > >("gpuBlockSize", Vector3< int >(128, 1, 1));
-      const bool cudaEnabledMpi = parameters.getParameter< bool >("cudaEnabledMpi", false);
+      const uint_t timesteps                    = parameters.getParameter< uint_t >("timesteps", uint_c(50));
+      const real_t remainingTimeLoggerFrequency = parameters.getParameter< real_t >("remainingTimeLoggerFrequency", real_c(3.0));
+      const uint_t scenario                     = parameters.getParameter< uint_t >("scenario", uint_c(1));
+      Vector3< int > gpuBlockSize               = parameters.getParameter< Vector3< int > >("gpuBlockSize", Vector3< int >(128, 1, 1));
+      const bool gpuEnabledMpi                 = parameters.getParameter< bool >("gpuEnabledMpi", true);
 
       /////////////////////////
       // ADD DATA TO BLOCKS //
@@ -115,13 +113,13 @@ int main(int argc, char** argv)
       BlockDataID phase_field = field::addToStorage< PhaseField_T >(blocks, "phase", real_c(0.0), field::fzyx);
       // GPU fields
       const BlockDataID lb_phase_field_gpu = gpu::addGPUFieldToStorage< gpu::GPUField< real_t > >(
-         blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, 1);
+         blocks, "lb phase field on GPU", Stencil_phase_T::Size, field::fzyx, uint_c(1));
       const BlockDataID lb_velocity_field_gpu = gpu::addGPUFieldToStorage< gpu::GPUField< real_t > >(
-         blocks, "lb velocity field on GPU", Stencil_hydro_T::Size, field::fzyx, 1);
+         blocks, "lb velocity field on GPU", Stencil_hydro_T::Size, field::fzyx, uint_c(1));
       BlockDataID vel_field_gpu =
          gpu::addGPUFieldToStorage< VelocityField_T >(blocks, vel_field, "velocity field on GPU", true);
-      BlockDataID phase_field_gpu =
-         gpu::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "phase field on GPU", true);
+      BlockDataID phase_field_gpu = gpu::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "phase field on GPU", true);
+      BlockDataID phase_field_tmp = gpu::addGPUFieldToStorage< PhaseField_T >(blocks, phase_field, "temporary phasefield", true);
       // Flag field
       const BlockDataID flagFieldID     = field::addFlagFieldToStorage< FlagField_T >(blocks, "flag field");
       const BlockDataID flagFieldID_gpu = gpu::addGPUFieldToStorage< FlagField_T >(blocks, flagFieldID, "flag on GPU", true);
@@ -137,6 +135,24 @@ int main(int argc, char** argv)
       const real_t relaxation_time_gas    = physical_parameters.getParameter< real_t >("relaxation_time_gas");
       const real_t interface_thickness    = physical_parameters.getParameter< real_t >("interface_thickness");
 
+      const real_t omega_l = real_c(1.0) / real_c(real_c(0.5) + relaxation_time_gas);
+      const real_t omega_h = real_c(1.0) / real_c(real_c(0.5) + relaxation_time_gas + (relaxation_time_liquid - relaxation_time_gas));
+      const real_t omega_phi = real_c(1.0) / real_c(real_c(0.5) + (real_c(3.0) * mobility));
+
+      WALBERLA_LOG_INFO_ON_ROOT(
+         "Benchmark run data:"
+         "\n- simulation parameters:"
+         << "\n   + density_liquid:             " << density_liquid
+         << "\n   + density_gas:                " << density_gas
+         << "\n   + omega liquid:               " << omega_h
+         << "\n   + omega gas:                  " << omega_l
+         << "\n   + surface_tension:            " << surface_tension
+         << "\n   + mobility:                   " << mobility
+         << "\n   + gravitational_acceleration: " << gravitational_acceleration
+         << "\n   + interface_thickness:        " << interface_thickness
+         << "\n   + omega phi:                  " << omega_phi)
+
+
       std::array< real_t, 3 > center_of_mass = { 0.0, 0.0, 0.0 };
 
       WALBERLA_LOG_INFO_ON_ROOT("initialization of the phase field")
@@ -177,59 +193,57 @@ int main(int argc, char** argv)
       // ADD SWEEPS //
       ///////////////
 
-      pystencils::initialize_phase_field_distributions init_h(lb_phase_field_gpu, phase_field_gpu, vel_field_gpu,
-                                                              interface_thickness);
+      pystencils::initialize_phase_field_distributions init_h(lb_phase_field_gpu, phase_field_gpu, vel_field_gpu, interface_thickness);
       pystencils::initialize_velocity_based_distributions init_g(lb_velocity_field_gpu, vel_field_gpu);
 
-      const pystencils::phase_field_LB_step phase_field_LB_step(flagFieldID_gpu, lb_phase_field_gpu, phase_field_gpu,
-                                                          vel_field_gpu, mobility, interface_thickness, gpuBlockSize[0],
-                                                          gpuBlockSize[1], gpuBlockSize[2]);
+      pystencils::phase_field_LB_step phase_field_LB_step(flagFieldID_gpu, lb_phase_field_gpu, phase_field_gpu, phase_field_tmp,
+                                                          vel_field_gpu, mobility, interface_thickness,
+                                                          gpuBlockSize[0], gpuBlockSize[1], gpuBlockSize[2]);
 
-      const pystencils::hydro_LB_step hydro_LB_step(flagFieldID_gpu, lb_velocity_field_gpu, phase_field_gpu, vel_field_gpu,
+      pystencils::hydro_LB_step hydro_LB_step(flagFieldID_gpu, lb_velocity_field_gpu, phase_field_gpu, vel_field_gpu,
                                               gravitational_acceleration, interface_thickness, density_liquid,
                                               density_gas, surface_tension, relaxation_time_liquid, relaxation_time_gas,
                                               gpuBlockSize[0], gpuBlockSize[1], gpuBlockSize[2]);
-
+      
+      auto swapPhaseField = std::function< void(IBlock *) >([&](IBlock * b)
+      {
+         auto phaseField    = b->getData< gpu::GPUField<real_t> >(phase_field_gpu);
+         auto phaseFieldTMP = b->getData< gpu::GPUField<real_t> >(phase_field_tmp);
+         phaseField->swapDataPointers(phaseFieldTMP);
+      });
       ////////////////////////
       // ADD COMMUNICATION //
       //////////////////////
+      if(gpuEnabledMpi)
+         WALBERLA_LOG_INFO_ON_ROOT("GPU-Direct communication is used for MPI")
+      else
+         WALBERLA_LOG_INFO_ON_ROOT("No GPU-Direct communication is used for MPI")
+
       const int streamLowPriority  = 0;
       const int streamHighPriority = 0;
       auto defaultStream     = gpu::StreamRAII::newPriorityStream(streamLowPriority);
       auto innerOuterStreams = gpu::ParallelStreams(streamHighPriority);
-
-      auto UniformGPUSchemeVelocityDistributions =
-         make_shared< gpu::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, cudaEnabledMpi);
-      auto generatedPackInfo_velocity_based_distributions =
-         make_shared< lbm::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
-      UniformGPUSchemeVelocityDistributions->addPackInfo(generatedPackInfo_velocity_based_distributions);
-      auto Comm_velocity_based_distributions =
-         std::function< void() >([&]() { UniformGPUSchemeVelocityDistributions->communicate(); });
-      auto Comm_velocity_based_distributions_start =
-         std::function< void() >([&]() { UniformGPUSchemeVelocityDistributions->startCommunication(); });
-      auto Comm_velocity_based_distributions_wait =
-         std::function< void() >([&]() { UniformGPUSchemeVelocityDistributions->wait(); });
-
-      auto UniformGPUSchemePhaseField =
-         make_shared< gpu::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, cudaEnabledMpi);
+      
+      auto generatedPackInfo_phase_field_distributions = make_shared< lbm::PackInfo_phase_field_distributions>(lb_phase_field_gpu);
+      auto generatedPackInfo_velocity_based_distributions = make_shared< lbm::PackInfo_velocity_based_distributions >(lb_velocity_field_gpu);
       auto generatedPackInfo_phase_field = make_shared< pystencils::PackInfo_phase_field >(phase_field_gpu);
+      
+      auto UniformGPUSchemeVelocityBasedDistributions = make_shared< gpu::communication::UniformGPUScheme< CommStencil_hydro_T > >(blocks, gpuEnabledMpi, false);
+      auto UniformGPUSchemePhaseFieldDistributions = make_shared< gpu::communication::UniformGPUScheme< CommStencil_phase_T > >(blocks, gpuEnabledMpi, false);
+      auto UniformGPUSchemePhaseField = make_shared< gpu::communication::UniformGPUScheme< Full_Stencil_T > >(blocks, gpuEnabledMpi, false, 65432);
+
+      UniformGPUSchemeVelocityBasedDistributions->addPackInfo(generatedPackInfo_velocity_based_distributions);
+      UniformGPUSchemePhaseFieldDistributions->addPackInfo(generatedPackInfo_phase_field_distributions);
       UniformGPUSchemePhaseField->addPackInfo(generatedPackInfo_phase_field);
+
+      auto Comm_velocity_based_distributions_start = std::function< void() >([&]() { UniformGPUSchemeVelocityBasedDistributions->startCommunication(); });
+      auto Comm_velocity_based_distributions_wait = std::function< void() >([&]() { UniformGPUSchemeVelocityBasedDistributions->wait(); });
+
+      auto Comm_phase_field_distributions = std::function< void() >([&]() { UniformGPUSchemePhaseFieldDistributions->communicate(); });
+      auto Comm_phase_field_distributions_start = std::function< void() >([&]() { UniformGPUSchemePhaseFieldDistributions->startCommunication(); });
+      auto Comm_phase_field_distributions_wait = std::function< void() >([&]() { UniformGPUSchemePhaseFieldDistributions->wait(); });
+
       auto Comm_phase_field = std::function< void() >([&]() { UniformGPUSchemePhaseField->communicate(); });
-      auto Comm_phase_field_start =
-         std::function< void() >([&]() { UniformGPUSchemePhaseField->startCommunication(); });
-      auto Comm_phase_field_wait = std::function< void() >([&]() { UniformGPUSchemePhaseField->wait(); });
-
-      auto UniformGPUSchemePhaseFieldDistributions =
-         make_shared< gpu::communication::UniformGPUScheme< Stencil_hydro_T > >(blocks, cudaEnabledMpi);
-      auto generatedPackInfo_phase_field_distributions =
-         make_shared< lbm::PackInfo_phase_field_distributions >(lb_phase_field_gpu);
-      UniformGPUSchemePhaseFieldDistributions->addPackInfo(generatedPackInfo_phase_field_distributions);
-      auto Comm_phase_field_distributions =
-         std::function< void() >([&]() { UniformGPUSchemePhaseFieldDistributions->communicate(); });
-      auto Comm_phase_field_distributions_start =
-         std::function< void() >([&]() { UniformGPUSchemePhaseFieldDistributions->startCommunication(); });
-      auto Comm_phase_field_distributions_wait =
-         std::function< void() >([&]() { UniformGPUSchemePhaseFieldDistributions->wait(); });
 
       ////////////////////////
       // BOUNDARY HANDLING //
@@ -273,16 +287,17 @@ int main(int argc, char** argv)
       SweepTimeloop timeloop(blocks->getBlockStorage(), timesteps);
 
       timeloop.add() << BeforeFunction(Comm_velocity_based_distributions_start, "Start Hydro PDFs Communication")
-                     << Sweep(phase_field_LB_NoSlip, "NoSlip Phase");
-      timeloop.add() << Sweep(phase_field_LB_step, "Phase LB Step")
+                     << Sweep(phase_field_LB_NoSlip.getSweep(defaultStream), "NoSlip Phase");
+      timeloop.add() << Sweep(phase_field_LB_step.getSweep(defaultStream), "Phase LB Step")
                      << AfterFunction(Comm_velocity_based_distributions_wait, "Wait Hydro PDFs Communication");
-      timeloop.add() << Sweep(contact_angle, "Contact Angle")
-                     << AfterFunction(Comm_phase_field, "Communication Phase Field");
 
       timeloop.add() << BeforeFunction(Comm_phase_field_distributions_start, "Start Phase PDFs Communication")
-                     << Sweep(hydro_LB_NoSlip, "NoSlip Hydro");
-      timeloop.add() << Sweep(hydro_LB_step, "Hydro LB Step")
+                     << Sweep(hydro_LB_NoSlip.getSweep(defaultStream), "NoSlip Hydro");
+      timeloop.add() << Sweep(hydro_LB_step.getSweep(defaultStream), "Hydro LB Step");
+      timeloop.add() << Sweep(swapPhaseField, "Swap PhaseField")
                      << AfterFunction(Comm_phase_field_distributions_wait, "Wait Phase PDFs Communication");
+      timeloop.add() << Sweep(contact_angle.getSweep(defaultStream), "Contact Angle")
+                     << AfterFunction(Comm_phase_field, "Communication Phase Field");
 
       if (scenario == 4)
       {
@@ -421,7 +436,10 @@ int main(int argc, char** argv)
             gpu::fieldCpy< PhaseField_T, GPUField >(blocks, phase_field, phase_field_gpu);
             gpu::fieldCpy< VelocityField_T, GPUField >(blocks, vel_field, vel_field_gpu);
          });
-         WALBERLA_GPU_CHECK(cudaPeekAtLastError())
+
+//         field::FlagFieldCellFilter<FlagField_T> fluidFilter( flagFieldID );
+//         fluidFilter.addFlag( wallFlagUID );
+//         vtkOutput->addCellExclusionFilter(fluidFilter);
 
          auto phaseWriter = make_shared< field::VTKWriter< PhaseField_T, float > >(phase_field, "PhaseField");
          vtkOutput->addCellDataWriter(phaseWriter);

apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_RTI_3D.py

@@ -61,7 +61,7 @@ class Scenario:
         # everything else
         self.scenario = 2  # 1 rising bubble or droplet, 2 RTI, 3 bubble field, 4 taylor bubble set up
 
-        self.cudaEnabledMpi = cuda_enabled_mpi
+        self.gpuEnabledMpi = cuda_enabled_mpi
         self.cuda_blocks = (64, 2, 2)
 
         self.config_dict = self.config()
@@ -82,7 +82,7 @@ class Scenario:
                 'dbWriteFrequency': self.dbWriteFrequency,
                 'remainingTimeLoggerFrequency': 10.0,
                 'scenario': self.scenario,
-                'cudaEnabledMpi': self.cudaEnabledMpi,
+                'gpuEnabledMpi': self.gpuEnabledMpi,
                 'gpuBlockSize': self.cuda_blocks
             },
             'PhysicalParameters': {

apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_codegen.py

@@ -5,7 +5,7 @@ from pystencils import fields, Target
 from pystencils.astnodes import Conditional, Block, SympyAssignment
 from pystencils.node_collection import NodeCollection
 from pystencils.typing import TypedSymbol
-from pystencils.simp.simplifications import sympy_cse
+from pystencils import simp as simplification
 
 from lbmpy import LBMConfig, LBMOptimisation, LBStencil, Method, Stencil, create_lb_update_rule
 from lbmpy.creationfunctions import create_lb_method
@@ -21,6 +21,11 @@ import pystencils_walberla
 from pystencils_walberla import CodeGeneration, generate_sweep, generate_pack_info_for_field, generate_info_header
 from lbmpy_walberla import generate_boundary, generate_lb_pack_info
 
+one = sp.Number(1)
+two = sp.Number(2)
+three = sp.Number(3)
+half = sp.Rational(1, 2)
+
 with CodeGeneration() as ctx:
     field_type = "float64" if ctx.double_accuracy else "float32"
 
@@ -56,7 +61,7 @@ with CodeGeneration() as ctx:
     ########################################
 
     # relaxation rate for interface tracking LB step
-    relaxation_rate_allen_cahn = 1.0 / (0.5 + (3.0 * parameters.symbolic_mobility))
+    relaxation_rate_allen_cahn = one / (half + (three * parameters.symbolic_mobility))
     # calculate the relaxation rate for hydrodynamic LB step
     rho_L = parameters.symbolic_density_light
     rho_H = parameters.symbolic_density_heavy
@@ -71,26 +76,22 @@ with CodeGeneration() as ctx:
     # LBM METHODS #
     ###############
 
-    rates = [0.0]
-    rates += [relaxation_rate_allen_cahn] * stencil_phase.D
-    rates += [1.0] * (stencil_phase.Q - stencil_phase.D - 1)
-
-    lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.MRT, compressible=True,
-                                 delta_equilibrium=False,
+    lbm_config_phase = LBMConfig(stencil=stencil_phase, method=Method.CENTRAL_MOMENT, compressible=True,
+                                 zero_centered=False,
                                  force=sp.symbols(f"F_:{stencil_phase.D}"), velocity_input=u,
-                                 weighted=True, relaxation_rates=rates,
-                                 output={'density': C_tmp}, kernel_type='stream_pull_collide')
+                                 weighted=True, relaxation_rates=[relaxation_rate_allen_cahn, ] * stencil_phase.Q,
+                                 output={'density': C_tmp})
     method_phase = create_lb_method(lbm_config=lbm_config_phase)
 
-    lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.MRT, compressible=False,
+    lbm_config_hydro = LBMConfig(stencil=stencil_hydro, method=Method.CENTRAL_MOMENT, compressible=False,
                                  weighted=True, relaxation_rate=omega,
                                  force=sp.symbols(f"F_:{stencil_hydro.D}"),
-                                 output={'velocity': u}, kernel_type='collide_stream_push')
+                                 output={'velocity': u})
     method_hydro = create_lb_method(lbm_config=lbm_config_hydro)
 
     # create the kernels for the initialization of the g and h field
     h_updates = initializer_kernel_phase_field_lb(method_phase, C, u, h, parameters)
-    g_updates = initializer_kernel_hydro_lb(method_hydro, 1, u, g)
+    g_updates = initializer_kernel_hydro_lb(method_hydro, 1.0, u, g)
 
     force_h = interface_tracking_force(C, stencil_phase, parameters)
     hydro_force = hydrodynamic_force(C, method_hydro, parameters, body_force)
@@ -105,7 +106,8 @@ with CodeGeneration() as ctx:
 
     allen_cahn_update_rule = add_interface_tracking_force(allen_cahn_update_rule, force_h)
 
-    allen_cahn_update_rule = sympy_cse(allen_cahn_update_rule)
+    allen_cahn_update_rule = simplification.add_subexpressions_for_field_reads(allen_cahn_update_rule)
+    # allen_cahn_update_rule = sympy_cse(allen_cahn_update_rule)
 
     allen_cahn_update_rule = NodeCollection([Conditional(sp.Eq(flag.center(), 2),
                                              Block(allen_cahn_update_rule.all_assignments),
@@ -117,7 +119,15 @@ with CodeGeneration() as ctx:
     hydro_lb_update_rule = create_lb_update_rule(lbm_config=lbm_config_hydro,
                                                  lbm_optimisation=lbm_optimisation)
 
-    hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters)
+    hydro_lb_update_rule = add_hydrodynamic_force(hydro_lb_update_rule, force_Assignments, C, g, parameters,
+                                                  lbm_config=lbm_config_hydro)
+    hydro_lb_update_rule = simplification.sympy_cse(hydro_lb_update_rule)
+    hydro_lb_update_rule = simplification.insert_aliases(hydro_lb_update_rule)
+    hydro_lb_update_rule = simplification.insert_constants(hydro_lb_update_rule)
+    hydro_lb_update_rule = simplification.insert_symbol_times_minus_one(hydro_lb_update_rule)
+    hydro_lb_update_rule = simplification.insert_squares(hydro_lb_update_rule)
+    hydro_lb_update_rule = simplification.insert_constant_multiples(hydro_lb_update_rule)
+    hydro_lb_update_rule = simplification.insert_constant_additions(hydro_lb_update_rule)
 
     hydro_lb_update_rule = NodeCollection([Conditional(sp.Eq(flag.center(), 2),
                                            Block(hydro_lb_update_rule.all_assignments))])
@@ -139,8 +149,12 @@ with CodeGeneration() as ctx:
 
     sweep_params = {'block_size': sweep_block_size}
 
-    stencil_typedefs = {'Stencil_phase_T': stencil_phase,
-                        'Stencil_hydro_T': stencil_hydro}
+    comm_stencil_phase = LBStencil(Stencil.D3Q27) if stencil_phase == LBStencil(Stencil.D3Q15) else stencil_phase
+    comm_stencil_hydro = LBStencil(Stencil.D3Q27) if stencil_hydro == LBStencil(Stencil.D3Q15) else stencil_hydro
+
+    stencil_typedefs = {'Stencil_phase_T': stencil_phase, 'CommStencil_phase_T': comm_stencil_phase,
+                        'Stencil_hydro_T': stencil_hydro, 'CommStencil_hydro_T': comm_stencil_hydro,
+                        'Full_Stencil_T': LBStencil(Stencil.D3Q27)}
     field_typedefs = {'PdfField_phase_T': h,
                       'PdfField_hydro_T': g,
                       'VelocityField_T': u,
@@ -155,7 +169,7 @@ with CodeGeneration() as ctx:
     generate_sweep(ctx, 'initialize_velocity_based_distributions', g_updates, target=Target.GPU)
 
     generate_sweep(ctx, 'phase_field_LB_step', allen_cahn_update_rule,
-                   field_swaps=[(h, h_tmp), (C, C_tmp)],
+                   field_swaps=[(h, h_tmp)],
                    target=Target.GPU,
                    gpu_indexing_params=sweep_params,
                    varying_parameters=vp)
@@ -165,8 +179,8 @@ with CodeGeneration() as ctx:
                    field_swaps=[(g, g_tmp)],
                    target=Target.GPU,
                    gpu_indexing_params=sweep_params,
-                   varying_parameters=vp)
-    generate_boundary(ctx, 'hydro_LB_NoSlip', NoSlip(), method_hydro, target=Target.GPU, streaming_pattern='push')
+                   varying_parameters=vp, max_threads=256)
+    generate_boundary(ctx, 'hydro_LB_NoSlip', NoSlip(), method_hydro, target=Target.GPU, streaming_pattern='pull')
 
     # communication
 
@@ -174,7 +188,7 @@ with CodeGeneration() as ctx:
                           streaming_pattern='pull', target=Target.GPU)
 
     generate_lb_pack_info(ctx, 'PackInfo_velocity_based_distributions', stencil_hydro, g,
-                          streaming_pattern='push', target=Target.GPU)
+                          streaming_pattern='pull', target=Target.GPU)
 
     generate_pack_info_for_field(ctx, 'PackInfo_phase_field', C, target=Target.GPU)
 

apps/showcases/PhaseFieldAllenCahn/GPU/multiphase_rising_bubble.py

@@ -54,7 +54,7 @@ class Scenario:
         self.counter = 0
         self.yPositions = []
 
-        self.cudaEnabledMpi = cuda_enabled_mpi
+        self.gpuEnabledMpi = cuda_enabled_mpi
         self.cuda_blocks = (64, 2, 2)
 
     @wlb.member_callback
@@ -72,7 +72,7 @@ class Scenario:
                 'dbWriteFrequency': self.dbWriteFrequency,
                 'remainingTimeLoggerFrequency': 10.0,
                 'scenario': self.scenario,
-                'cudaEnabledMpi': self.cudaEnabledMpi,
+                'gpuEnabledMpi': self.gpuEnabledMpi,
                 'gpuBlockSize': self.cuda_blocks
             },
             'PhysicalParameters': {

apps/showcases/PhaseFieldAllenCahn/GPU/taylor_bubble.py

@@ -61,20 +61,20 @@ class Scenario:
         self.reference_length = 128
         d1 = 0.0254  # (0.0508, 0.0381, 0.0254)
         d2 = 0.0127  # (0.0254, 0.0127, 0.0127)
-        self.interface_width = 4
-        self.mobility = 0.05
+        self.interface_width = 5
+        self.mobility = 0.1
 
         num_processes = wlb.mpi.numProcesses()
 
         # output frequencies
-        self.vtkWriteFrequency = self.reference_time
-        self.dbWriteFrequency = self.reference_time // 25
+        self.vtkWriteFrequency = 2000 # self.reference_time // 25
+        self.dbWriteFrequency = 0 # self.reference_time // 25
         self.meshWriteFrequency = self.reference_time
         self.pngWriteFrequency = self.reference_time
 
         # simulation parameters
         self.diameter = self.reference_length
-        self.timesteps = self.reference_time * 15 + 1
+        self.timesteps = 30001 # self.reference_time #  * 15 + 1
         self.cells = (self.diameter, (self.reference_length * 15) // num_processes, self.diameter)
         self.blocks = (1, num_processes, 1)
         self.periodic = (0, 0, 0)
@@ -126,7 +126,7 @@ class Scenario:
         self.relaxation_time_liquid = parameters.relaxation_time_heavy
         self.relaxation_time_gas = parameters.relaxation_time_light
 
-        self.cudaEnabledMpi = cuda_enabled_mpi
+        self.gpuEnabledMpi = cuda_enabled_mpi
         self.cuda_blocks = (64, 2, 2)
 
         self.config_dict = self.config()
@@ -150,9 +150,9 @@ class Scenario:
                 'vtkWriteFrequency': self.vtkWriteFrequency,
                 'dbWriteFrequency': self.dbWriteFrequency,
                 'meshWriteFrequency': self.meshWriteFrequency,
-                'remainingTimeLoggerFrequency': 60.0,
+                'remainingTimeLoggerFrequency': 20.0,
                 'scenario': self.scenario,
-                'cudaEnabledMpi': self.cudaEnabledMpi,
+                'gpuEnabledMpi': self.gpuEnabledMpi,
                 'gpuBlockSize': self.cuda_blocks
             },
             'PhysicalParameters': {
@@ -180,6 +180,9 @@ class Scenario:
                 'Donut_h': self.Donut_h,
                 'Donut_D': self.Donut_D,
                 'DonutTime': self.DonutTime
+            },
+            'Logging': {
+                'logLevel': 'info',  # info progress detail tracing
             }
         }
 
@@ -325,4 +328,4 @@ class Scenario:
 
 
 scenarios = wlb.ScenarioManager()
-scenarios.add(Scenario())
+scenarios.add(Scenario(cuda_enabled_mpi=False))