Add free-surface module

apps/benchmarks/ComplexGeometry/ComplexGeometry.cpp

@@ -239,7 +239,7 @@ int main( int argc, char * argv[] )
                                                                                                              flagFieldId, fluidFlagUID ) ),
                                   "LBM stability check" );
 
-   timeloop.addFuncAfterTimeStep( perfLogger, "PerformanceLogger" );
+   timeloop.addFuncAfterTimeStep( perfLogger, "Evaluator: performance logging" );
 
 
    // add VTK output to time loop

apps/showcases/CMakeLists.txt

 add_subdirectory( BidisperseFluidizedBed )
 add_subdirectory( CombinedResolvedUnresolved )
 add_subdirectory( FluidizedBed )
+add_subdirectory( FreeSurface )
 add_subdirectory( HeatConduction )
 add_subdirectory( LightRisingParticleInFluidAMR )
 add_subdirectory( Mixer )
@@ -12,3 +13,4 @@ endif()
 if ( WALBERLA_BUILD_WITH_CODEGEN AND NOT WALBERLA_BUILD_WITH_OPENMP)
    add_subdirectory( PorousMedia )
 endif()
+

apps/showcases/FreeSurface/BubblyPoiseuille.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 BubblyPoiseuille.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+// This showcase simulates a plane Poiseuille flow with randomly distributed bubbles in the flow.
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+#include "core/math/Random.h"
+
+#include "lbm/PerformanceLogger.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/LoadBalancing.h"
+#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/VtkWriter.h"
+#include "lbm/free_surface/bubble_model/Geometry.h"
+#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
+#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
+#include "lbm/free_surface/surface_geometry/Utility.h"
+#include "lbm/lattice_model/D3Q19.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+namespace DropInPool
+{
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+
+using CollisionModel_T      = lbm::collision_model::SRT;
+using ForceModel_T          = lbm::force_model::GuoField< VectorField_T >;
+using LatticeModel_T        = lbm::D3Q19< CollisionModel_T, true, ForceModel_T, 2 >;
+using LatticeModelStencil_T = LatticeModel_T::Stencil;
+using PdfField_T            = lbm::PdfField< LatticeModel_T >;
+using PdfCommunication_T    = blockforest::SimpleCommunication< LatticeModelStencil_T >;
+
+// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
+// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
+using CommunicationStencil_T =
+   typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
+using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
+
+using flag_t                        = uint32_t;
+using FlagField_T                   = FlagField< flag_t >;
+using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
+
+int main(int argc, char** argv)
+{
+   Environment walberlaEnv(argc, argv);
+
+   if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
+
+   // print content of parameter file
+   WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
+
+   // get block forest parameters from parameter file
+   auto blockForestParameters              = walberlaEnv.config()->getOneBlock("BlockForestParameters");
+   const Vector3< uint_t > cellsPerBlock   = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+   const Vector3< bool > periodicity       = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
+   const uint_t loadBalancingFrequency     = blockForestParameters.getParameter< uint_t >("loadBalancingFrequency");
+   const bool printLoadBalancingStatistics = blockForestParameters.getParameter< bool >("printLoadBalancingStatistics");
+
+   // read domain parameters from parameter file
+   const auto domainParameters              = walberlaEnv.config()->getOneBlock("DomainParameters");
+   const real_t channelWidth                = domainParameters.getParameter< real_t >("channelWidth");
+   const Vector3< real_t > domainSizeFactor = domainParameters.getParameter< Vector3< real_t > >("domainSizeFactor");
+   const real_t bubbleDiameter    = domainParameters.getParameter< real_t >("bubbleDiameterFactor") * channelWidth;
+   const real_t gasVolumeFraction = domainParameters.getParameter< real_t >("gasVolumeFraction");
+
+   // define domain size
+   Vector3< uint_t > domainSize = domainSizeFactor * channelWidth;
+   domainSize[0]                = uint_c(domainSizeFactor[0] * channelWidth);
+   domainSize[1]                = uint_c(domainSizeFactor[1] * channelWidth);
+   domainSize[2]                = uint_c(domainSizeFactor[2] * channelWidth);
+
+   // compute number of blocks as defined by domainSize and cellsPerBlock
+   Vector3< uint_t > numBlocks;
+   numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
+   numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
+   numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
+
+   // get number of (MPI) processes
+   uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
+   WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
+                             "The number of MPI processes is greater than the number of blocks as defined by "
+                             "\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
+                             "the number of MPI processes or increase \"cellsPerBlock\".")
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(channelWidth);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSizeFactor);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(gasVolumeFraction);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleDiameter);
+
+   Vector3< uint_t > realDomainSize;
+   realDomainSize[0] = cellsPerBlock[0] * numBlocks[0];
+   realDomainSize[1] = cellsPerBlock[1] * numBlocks[1];
+   realDomainSize[2] = cellsPerBlock[2] * numBlocks[2];
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(realDomainSize);
+
+   if (domainSize[0] != realDomainSize[0] && periodicity[0])
+   {
+      WALBERLA_ABORT(
+         "The specified domain size in x-direction can not be obtained with the number of blocks you specified.")
+   }
+   if (domainSize[1] != realDomainSize[1] && periodicity[1])
+   {
+      WALBERLA_ABORT(
+         "The specified domain size in y-direction can not be obtained with the number of blocks you specified.")
+   }
+   if (domainSize[2] != realDomainSize[2] && periodicity[2])
+   {
+      WALBERLA_ABORT(
+         "The specified domain size in z-direction can not be obtained with the number of blocks you specified.")
+   }
+
+   // read physics parameters from parameter file
+   const auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
+   const real_t reynoldsNumber  = physicsParameters.getParameter< real_t >("reynoldsNumber");
+   const real_t mortonNumber    = physicsParameters.getParameter< real_t >("mortonNumber");
+   const real_t relaxationRate  = physicsParameters.getParameter< real_t >("relaxationRate");
+
+   const CollisionModel_T collisionModel = CollisionModel_T(relaxationRate);
+   const real_t viscosity                = collisionModel.viscosity();
+
+   const real_t forceX = real_c(8) * viscosity * viscosity * reynoldsNumber / real_c(std::pow(channelWidth, 3));
+   const Vector3< real_t > force = Vector3< real_t >(forceX, real_c(0), real_c(0));
+
+   const real_t analMaxVelocity = viscosity * reynoldsNumber / channelWidth;
+
+   const real_t surfaceTension = real_c(std::pow(forceX * real_c(std::pow(viscosity, 4)) / mortonNumber,
+                                                 real_c(1) / real_c(3))); // formula only valid for rho=1
+
+   const bool enableWetting  = physicsParameters.getParameter< bool >("enableWetting");
+   const real_t contactAngle = physicsParameters.getParameter< real_t >("contactAngle");
+
+   const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(reynoldsNumber);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(mortonNumber);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(relaxationRate);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(surfaceTension);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableWetting);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(contactAngle);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(analMaxVelocity);
+
+   // read model parameters from parameter file
+   const auto modelParameters               = walberlaEnv.config()->getOneBlock("ModelParameters");
+   const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
+   const std::string pdfRefillingModel      = modelParameters.getParameter< std::string >("pdfRefillingModel");
+   const std::string excessMassDistributionModel =
+      modelParameters.getParameter< std::string >("excessMassDistributionModel");
+   const std::string curvatureModel          = modelParameters.getParameter< std::string >("curvatureModel");
+   const bool enableForceWeighting           = modelParameters.getParameter< bool >("enableForceWeighting");
+   const bool useSimpleMassExchange          = modelParameters.getParameter< bool >("useSimpleMassExchange");
+   const bool enableBubbleModel              = modelParameters.getParameter< bool >("enableBubbleModel");
+   const bool enableBubbleSplits             = modelParameters.getParameter< bool >("enableBubbleSplits");
+   const real_t cellConversionThreshold      = modelParameters.getParameter< real_t >("cellConversionThreshold");
+   const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
+
+   // read evaluation parameters from parameter file
+   const auto evaluationParameters      = walberlaEnv.config()->getOneBlock("EvaluationParameters");
+   const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+
+   // create non-uniform block forest (non-uniformity required for load balancing)
+   const std::shared_ptr< StructuredBlockForest > blockForest =
+      createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
+
+   // add force field
+   const BlockDataID forceFieldID =
+      field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
+
+   // create lattice model
+   const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+
+   // add pdf field
+   const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel,
+                                                            Vector3< real_t >(real_c(0)), real_c(1), field::fzyx);
+
+   // add fill level field (initialized with 0, i.e., gas everywhere)
+   const BlockDataID fillFieldID =
+      field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
+
+   // add boundary handling
+   const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
+      std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
+   const BlockDataID flagFieldID                                      = freeSurfaceBoundaryHandling->getFlagFieldID();
+   const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
+
+   // initialize parabolic Poiseuille profile
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+
+      WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, {
+         // get global coordinate (with respect to whole simulation domain) of the currently processed cell
+         Vector3< real_t > globalCellCoordinate = blockForest->getBlockLocalCellCenter(*blockIt, pdfFieldIt.cell());
+
+         const real_t height = real_c(realDomainSize[2]);
+
+         const real_t velocityX =
+            forceX * real_c(0.5) / viscosity * globalCellCoordinate[2] * (height - globalCellCoordinate[2]);
+
+         pdfField->setDensityAndVelocity(pdfFieldIt.cell(), Vector3< real_t >(velocityX, real_c(0), real_c(0)),
+                                         real_c(1));
+      }); // WALBERLA_FOR_ALL_CELLS
+   }
+
+   const real_t bubbleRadius          = bubbleDiameter * real_c(0.5);
+   const real_t domainVolume          = real_c(realDomainSize[0] * realDomainSize[1] * realDomainSize[2]);
+   const real_t gasVolume             = gasVolumeFraction * domainVolume;
+   const real_t bubbleVolume          = real_c(4) / real_c(3) * real_c(math::pi) * real_c(std::pow(bubbleRadius, 3));
+   const uint_t numBubbles            = uint_c(gasVolume / bubbleVolume);
+   const real_t realGasVolumeFraction = real_c(numBubbles) * real_c(bubbleVolume) / domainVolume;
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBubbles);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(realGasVolumeFraction);
+
+   // randomly place bubbles
+   for (uint_t bubble = uint_c(0); bubble != numBubbles; ++bubble)
+   {
+      walberla::math::seedRandomGenerator(std::mt19937::result_type(bubble));
+      const Vector3< real_t > bubbleCenter =
+         Vector3< real_t >(math::realRandom(bubbleRadius, real_c(realDomainSize[0]) - bubbleRadius),
+                           math::realRandom(bubbleRadius, real_c(realDomainSize[1]) - bubbleRadius),
+                           math::realRandom(bubbleRadius, real_c(realDomainSize[2]) - bubbleRadius));
+
+      const geometry::Sphere sphereBubble(bubbleCenter, bubbleRadius);
+      bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
+   }
+
+   // initialize boundary conditions from config file
+   const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
+   freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
+
+   // IMPORTANT REMARK: this must be only called after every solid flag has been set; otherwise, the boundary handling
+   // might not detect solid flags correctly
+   freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
+
+   // communication after initialization
+   Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+   communication();
+
+   PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
+   pdfCommunication();
+
+   // add bubble model
+   std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel = nullptr;
+   if (enableBubbleModel)
+   {
+      const std::shared_ptr< bubble_model::BubbleModel< CommunicationStencil_T > > bubbleModelDerived =
+         std::make_shared< bubble_model::BubbleModel< CommunicationStencil_T > >(blockForest, enableBubbleSplits);
+      bubbleModelDerived->initFromFillLevelField(fillFieldID);
+
+      bubbleModel = std::static_pointer_cast< bubble_model::BubbleModelBase >(bubbleModelDerived);
+   }
+   else { bubbleModel = std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1)); }
+
+   // set density in non-liquid or non-interface cells to 1 (after initializing with hydrostatic pressure)
+   setDensityInNonFluidCellsToOne< FlagField_T, PdfField_T >(blockForest, flagInfo, flagFieldID, pdfFieldID);
+
+   // create timeloop
+   SweepTimeloop timeloop(blockForest, timesteps);
+
+   // Laplace pressure = 2 * surface tension * curvature; curvature computation is not necessary with 0 surface tension
+   bool computeCurvature = false;
+   if (!realIsEqual(surfaceTension, real_c(0), real_c(1e-14))) { computeCurvature = true; }
+
+   // add surface geometry handler
+   const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
+      blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, computeCurvature, enableWetting,
+      contactAngle);
+
+   geometryHandler.addSweeps(timeloop);
+
+   const ConstBlockDataID curvatureFieldID = geometryHandler.getConstCurvatureFieldID();
+   const ConstBlockDataID normalFieldID    = geometryHandler.getConstNormalFieldID();
+
+   // add boundary handling for standard boundaries and free surface boundaries
+   const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
+      blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+      freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
+      relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+      cellConversionForceThreshold);
+
+   dynamicsHandler.addSweeps(timeloop);
+
+   // add load balancing
+   const LoadBalancer< FlagField_T, CommunicationStencil_T, LatticeModelStencil_T > loadBalancer(
+      blockForest, communication, pdfCommunication, bubbleModel, uint_c(50), uint_c(10), uint_c(5),
+      loadBalancingFrequency, printLoadBalancingStatistics);
+   timeloop.addFuncAfterTimeStep(loadBalancer, "Sweep: load balancing");
+
+   // add VTK output
+   addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
+      blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+      geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
+      geometryHandler.getObstNormalFieldID());
+
+   // add triangle mesh output of free surface
+   SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
+      blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
+   surfaceMeshWriter(); // write initial mesh
+   timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
+
+   // add logging for computational performance
+   const lbm::PerformanceLogger< FlagField_T > perfLogger(blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs,
+                                                          performanceLogFrequency);
+   timeloop.addFuncAfterTimeStep(perfLogger, "Evaluator: performance logging");
+
+   WcTimingPool timingPool;
+
+   for (uint_t t = uint_c(0); t != timesteps; ++t)
+   {
+      if (t % uint_c(real_c(timesteps / 100)) == uint_c(0))
+      {
+         WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep = " << t);
+      }
+      timeloop.singleStep(timingPool, true);
+
+      if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
+   }
+
+   return EXIT_SUCCESS;
+}
+
+} // namespace DropInPool
+} // namespace free_surface
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::free_surface::DropInPool::main(argc, argv); }
\ No newline at end of file

apps/showcases/FreeSurface/BubblyPoiseuille.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 20, 20, 20 >;
+   periodicity                   < 1, 1, 0 >;
+   loadBalancingFrequency        0;
+   printLoadBalancingStatistics  false;
+}
+
+DomainParameters
+{
+   channelWidth            100;
+   domainSizeFactor        < 2, 1, 1 >;   // values multiplied with channelWidth
+   bubbleDiameterFactor    0.1;           // value multiplied with channelWidth
+   gasVolumeFraction       0.1;
+}
+
+PhysicsParameters
+{
+   reynoldsNumber    10000;
+   mortonNumber      1e-5;
+   relaxationRate    1.989;
+   enableWetting     false;
+   contactAngle      0;       // only used if enableWetting=true
+   timesteps         10000;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   enableBubbleModel             true;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 5000;
+}
+
+BoundaryParameters
+{
+   // X
+   //Border { direction W;  walldistance -1; NoSlip{} }
+   //Border { direction E;  walldistance -1; NoSlip{} }
+
+   // Y
+   //Border { direction N;  walldistance -1; NoSlip{} }
+   //Border { direction S;  walldistance -1; NoSlip{} }
+
+   // Z
+   Border { direction T;  walldistance -1; NoSlip{} }
+   Border { direction B;  walldistance -1; NoSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 2170;
+   baseFolder mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency 2170;
+      ghostLayers 0;
+      baseFolder vtk-out;
+      samplingResolution 1;
+
+      writers
+      {
+         velocity;
+         density;
+         //pdf;
+         flag;
+         fill_level;
+         force;
+         curvature;
+         normal;
+         obstacle_normal;
+         mapped_flag;
+        }
+
+        inclusion_filters
+        {
+            // only include liquid and interface cells in VTK output
+            //liquidInterfaceFilter;
+        }
+
+        before_functions
+        {
+            //ghost_layer_synchronization; // only needed if writing the ghost layer
+        }
+
+   }
+   domain_decomposition
+   {
+      writeFrequency 0;
+      baseFolder vtk-out;
+      outputDomainDecomposition true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/CMakeLists.txt

+waLBerla_link_files_to_builddir( *.prm )
+
+waLBerla_add_executable(NAME    BubblyPoiseuille
+                        FILES   BubblyPoiseuille.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    CapillaryWave
+                        FILES   CapillaryWave.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    DamBreakCylindrical
+                        FILES   DamBreakCylindrical.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    DamBreakRectangular
+                        FILES   DamBreakRectangular.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    DropImpact
+                        FILES   DropImpact.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    DropWetting
+                        FILES   DropWetting.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    GravityWave
+                        FILES   GravityWave.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+if( WALBERLA_BUILD_WITH_CODEGEN )
+   walberla_generate_target_from_python( NAME      GravityWaveLatticeModelGeneration
+                                         FILE      GravityWaveLatticeModelGeneration.py
+                                         OUT_FILES GravityWaveLatticeModel.cpp GravityWaveLatticeModel.h )
+
+   waLBerla_add_executable(NAME    GravityWaveCodegen
+                           FILES   GravityWaveCodegen.cpp
+                           DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk
+                                   GravityWaveLatticeModelGeneration)
+endif()
+
+waLBerla_add_executable(NAME    RisingBubble
+                        FILES   RisingBubble.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME    TaylorBubble
+                        FILES   TaylorBubble.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)

apps/showcases/FreeSurface/CapillaryWave.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 25, 25, 1 >;
+   periodicity                   < 1, 0, 1 >;
+   loadBalancingFrequency        50;
+   printLoadBalancingStatistics  true;
+}
+
+DomainParameters
+{
+   domainWidth       50; // equivalent to wavelength
+   liquidDepth       25;
+   initialAmplitude  0.5;
+}
+
+PhysicsParameters
+{
+   reynoldsNumber    10;
+   relaxationRate    1.8;
+   enableWetting     false;
+   contactAngle      0; // only used if enableWetting=true
+   timesteps         10000;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+
+   enableBubbleModel             false;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 5000;
+   evaluationFrequency     100;
+   filename                capillary-wave.txt;
+}
+
+BoundaryParameters
+{
+   // X
+   //Border { direction W;  walldistance -1; NoSlip{} }
+   //Border { direction E;  walldistance -1; NoSlip{} }
+
+   // Y
+   Border { direction N;  walldistance -1; NoSlip{} }
+   Border { direction S;  walldistance -1; NoSlip{} }
+
+   // Z
+   //Border { direction T;  walldistance -1; NoSlip{} }
+   //Border { direction B;  walldistance -1; NoSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 100;
+   baseFolder     mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency       100;
+      ghostLayers          0;
+      baseFolder           vtk-out;
+      samplingResolution   1;
+
+      writers
+      {
+         fill_level;
+         mapped_flag;
+         velocity;
+         density;
+         //curvature;
+         //normal;
+         //obstacle_normal;
+         //pdf;
+         //flag;
+         //force;
+      }
+
+      inclusion_filters
+      {
+         //liquidInterfaceFilter; // only include liquid and interface cells in VTK output
+      }
+
+      before_functions
+      {
+         //ghost_layer_synchronization; // only needed if writing the ghost layer
+      }
+   }
+   domain_decomposition
+   {
+      writeFrequency             100;
+      baseFolder                 vtk-out;
+      outputDomainDecomposition  true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/DamBreakCylindrical.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 10, 5, 10 >;
+   periodicity                   < 1, 0, 1 >;
+   loadBalancingFrequency        248;
+   printLoadBalancingStatistics  true;
+}
+
+DomainParameters
+{
+   columnRadius 25;  // initial radius of the liquid column; "a" in Martin & Moyce's paper (10.1098/rsta.1952.0006)
+   columnRatio 1;    // ratio between columnHeight and columnRadius; "n^2" in Martin & Moyce's paper
+}
+
+PhysicsParameters
+{
+   galileiNumber     1831123817;
+   bondNumber        445;
+   relaxationRate    1.9995;
+   enableWetting     false;
+   contactAngle      0; // only used if enableWetting=true
+   timesteps         2480;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+
+   enableBubbleModel             false;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+
+   smagorinskyConstant           0.1;
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 2480;
+   evaluationFrequency     24;
+   filename                rot-breaking-dam.txt;
+}
+
+BoundaryParameters
+{
+   // X
+   //Border { direction W;  walldistance -1;  FreeSlip{} }
+   //Border { direction E;  walldistance -1;  FreeSlip{} }
+
+   // Y
+   Border { direction N;  walldistance -1;  FreeSlip{} }
+   Border { direction S;  walldistance -1;  FreeSlip{} }
+
+   // Z
+   //Border { direction T;  walldistance -1; FreeSlip{} }
+   //Border { direction B;  walldistance -1; FreeSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 248;
+   baseFolder mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency 248;
+      ghostLayers 0;
+      baseFolder vtk-out;
+      samplingResolution 1;
+
+      writers
+      {
+         velocity;
+         density;
+         //pdf;
+         flag;
+         fill_level;
+         force;
+         curvature;
+         normal;
+         obstacle_normal;
+         mapped_flag;
+        }
+
+        inclusion_filters
+        {
+            // only include liquid and interface cells in VTK output
+            //liquidInterfaceFilter;
+        }
+
+        before_functions
+        {
+            //ghost_layer_synchronization; // only needed if writing the ghost layer
+        }
+
+   }
+   domain_decomposition
+   {
+      writeFrequency 248;
+      baseFolder vtk-out;
+      outputDomainDecomposition true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/DamBreakRectangular.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 50, 50, 1 >;
+   periodicity                   < 0, 0, 1 >;
+   loadBalancingFrequency        991;
+   printLoadBalancingStatistics  true;
+}
+
+DomainParameters
+{
+   columnWidth 50;   // initial width of the liquid column; "a" in Martin & Moyce's paper (10.1098/rsta.1952.0006)
+   columnRatio 2;    // ratio between columnHeight and columnWidth; "n^2" in Martin & Moyce's paper
+}
+
+PhysicsParameters
+{
+   galileiNumber     1831123817;
+   bondNumber        445;
+   relaxationRate    1.9995;
+   enableWetting     false;
+   contactAngle      0; // only used if enableWetting=true
+   timesteps         9910;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+
+   enableBubbleModel             false;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+
+   smagorinskyConstant           0.1;
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 25000;
+   evaluationFrequency     99;
+   filename                breaking-dam.txt;
+}
+
+BoundaryParameters
+{
+   // X
+   Border { direction W;  walldistance -1;  FreeSlip{} }
+   Border { direction E;  walldistance -1;  FreeSlip{} }
+
+   // Y
+   Border { direction N;  walldistance -1;  FreeSlip{} }
+   Border { direction S;  walldistance -1;  FreeSlip{} }
+
+   // Z
+   //Border { direction T;  walldistance -1; FreeSlip{} }
+   //Border { direction B;  walldistance -1; FreeSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 0;
+   baseFolder mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency 991;
+      ghostLayers 0;
+      baseFolder vtk-out;
+      samplingResolution 1;
+
+      writers
+      {
+         velocity;
+         density;
+         //pdf;
+         flag;
+         fill_level;
+         force;
+         curvature;
+         normal;
+         obstacle_normal;
+         mapped_flag;
+        }
+
+        inclusion_filters
+        {
+            // only include liquid and interface cells in VTK output
+            //liquidInterfaceFilter;
+        }
+
+        before_functions
+        {
+            //ghost_layer_synchronization; // only needed if writing the ghost layer
+        }
+
+   }
+   domain_decomposition
+   {
+      writeFrequency 991;
+      baseFolder vtk-out;
+      outputDomainDecomposition true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/DropImpact.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 DropImpact.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+// This showcase simulates the impact of a droplet into a pool of liquid. Reference experiments are available from
+// - Wang, Chen (2000), doi:10.1063/1.1287511 (vertical impact)
+// - Reijers, Liu, Lohse, Gelderblom (2019), url: http://arxiv.org/abs/1903.08978 (oblique impact)
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+
+#include "lbm/PerformanceLogger.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/LoadBalancing.h"
+#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/VtkWriter.h"
+#include "lbm/free_surface/bubble_model/Geometry.h"
+#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
+#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
+#include "lbm/free_surface/surface_geometry/Utility.h"
+#include "lbm/lattice_model/D3Q19.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+namespace DropImpact
+{
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+
+using CollisionModel_T      = lbm::collision_model::SRT;
+using ForceModel_T          = lbm::force_model::GuoField< VectorField_T >;
+using LatticeModel_T        = lbm::D3Q19< CollisionModel_T, true, ForceModel_T, 2 >;
+using LatticeModelStencil_T = LatticeModel_T::Stencil;
+using PdfField_T            = lbm::PdfField< LatticeModel_T >;
+using PdfCommunication_T    = blockforest::SimpleCommunication< LatticeModelStencil_T >;
+
+// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
+// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
+using CommunicationStencil_T =
+   typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
+using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
+
+using flag_t                        = uint32_t;
+using FlagField_T                   = FlagField< flag_t >;
+using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
+
+int main(int argc, char** argv)
+{
+   Environment walberlaEnv(argc, argv);
+
+   if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
+
+   // print content of parameter file
+   WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
+
+   // get block forest parameters from parameter file
+   auto blockForestParameters              = walberlaEnv.config()->getOneBlock("BlockForestParameters");
+   const Vector3< uint_t > cellsPerBlock   = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+   const Vector3< bool > periodicity       = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
+   const uint_t loadBalancingFrequency     = blockForestParameters.getParameter< uint_t >("loadBalancingFrequency");
+   const bool printLoadBalancingStatistics = blockForestParameters.getParameter< bool >("printLoadBalancingStatistics");
+
+   // read domain parameters from parameter file
+   const auto domainParameters              = walberlaEnv.config()->getOneBlock("DomainParameters");
+   const real_t dropDiameter                = domainParameters.getParameter< real_t >("dropDiameter");
+   const Vector3< real_t > dropCenterFactor = domainParameters.getParameter< Vector3< real_t > >("dropCenterFactor");
+   const real_t poolHeightFactor            = domainParameters.getParameter< real_t >("poolHeightFactor");
+   const Vector3< real_t > domainSizeFactor = domainParameters.getParameter< Vector3< real_t > >("domainSizeFactor");
+
+   // define domain size
+   Vector3< uint_t > domainSize = domainSizeFactor * dropDiameter;
+   domainSize[0]                = uint_c(domainSizeFactor[0] * dropDiameter);
+   domainSize[1]                = uint_c(domainSizeFactor[1] * dropDiameter);
+   domainSize[2]                = uint_c(domainSizeFactor[2] * dropDiameter);
+
+   // compute number of blocks as defined by domainSize and cellsPerBlock
+   Vector3< uint_t > numBlocks;
+   numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
+   numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
+   numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
+
+   // get number of (MPI) processes
+   uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
+   WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
+                             "The number of MPI processes is greater than the number of blocks as defined by "
+                             "\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
+                             "the number of MPI processes or increase \"cellsPerBlock\".")
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(dropDiameter);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(dropCenterFactor);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(poolHeightFactor);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSizeFactor);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
+
+   Vector3< uint_t > realDomainSize;
+   realDomainSize[0] = cellsPerBlock[0] * numBlocks[0];
+   realDomainSize[1] = cellsPerBlock[1] * numBlocks[1];
+   realDomainSize[2] = cellsPerBlock[2] * numBlocks[2];
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(realDomainSize);
+
+   if (domainSize[0] != realDomainSize[0] && periodicity[0])
+   {
+      WALBERLA_ABORT(
+         "The specified domain size in x-direction can not be obtained with the number of blocks you specified.")
+   }
+   if (domainSize[1] != realDomainSize[1] && periodicity[1])
+   {
+      WALBERLA_ABORT(
+         "The specified domain size in y-direction can not be obtained with the number of blocks you specified.")
+   }
+   // the z-direction must be no slip in this setup and is simply extended (see below) to obtain the specified size
+   int boundaryThicknessZ = int_c(realDomainSize[2]) - int_c(domainSize[2]);
+   if (boundaryThicknessZ < 0)
+   {
+      WALBERLA_ABORT("Something went wrong: the resulting domain size in z-direction is less than specified.")
+   }
+
+   // read physics parameters from parameter file
+   const auto physicsParameters   = walberlaEnv.config()->getOneBlock("PhysicsParameters");
+   const real_t bondNumber        = physicsParameters.getParameter< real_t >("bondNumber");
+   const real_t weberNumber       = physicsParameters.getParameter< real_t >("weberNumber");
+   const real_t ohnesorgeNumber   = physicsParameters.getParameter< real_t >("ohnesorgeNumber");
+   const real_t relaxationRate    = physicsParameters.getParameter< real_t >("relaxationRate");
+   const real_t impactAngleDegree = physicsParameters.getParameter< real_t >("impactAngleDegree");
+
+   const CollisionModel_T collisionModel = CollisionModel_T(relaxationRate);
+   const real_t viscosity                = collisionModel.viscosity();
+
+   const real_t surfaceTension             = real_c(std::pow(viscosity / ohnesorgeNumber, 2)) / dropDiameter;
+   const real_t gravitationalAccelerationZ = bondNumber * surfaceTension / (dropDiameter * dropDiameter);
+   const real_t impactVelocityMagnitude    = real_c(std::pow(weberNumber * surfaceTension / dropDiameter, 0.5));
+
+   const Vector3< real_t > force(real_c(0), real_c(0), -gravitationalAccelerationZ);
+
+   const real_t impactVelocityY =
+      impactVelocityMagnitude * real_c(std::sin(impactAngleDegree * math::pi / real_c(180)));
+   const real_t impactVelocityZ =
+      impactVelocityMagnitude * real_c(std::cos(impactAngleDegree * math::pi / real_c(180)));
+
+   const Vector3< real_t > impactVelocity(real_c(0), impactVelocityY, -impactVelocityZ);
+
+   const bool enableWetting  = physicsParameters.getParameter< bool >("enableWetting");
+   const real_t contactAngle = physicsParameters.getParameter< real_t >("contactAngle");
+
+   const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(relaxationRate);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(surfaceTension);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(impactVelocity);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableWetting);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(contactAngle);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
+
+   // read model parameters from parameter file
+   const auto modelParameters               = walberlaEnv.config()->getOneBlock("ModelParameters");
+   const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
+   const std::string pdfRefillingModel      = modelParameters.getParameter< std::string >("pdfRefillingModel");
+   const std::string excessMassDistributionModel =
+      modelParameters.getParameter< std::string >("excessMassDistributionModel");
+   const std::string curvatureModel          = modelParameters.getParameter< std::string >("curvatureModel");
+   const bool enableForceWeighting           = modelParameters.getParameter< bool >("enableForceWeighting");
+   const bool useSimpleMassExchange          = modelParameters.getParameter< bool >("useSimpleMassExchange");
+   const bool enableBubbleModel              = modelParameters.getParameter< bool >("enableBubbleModel");
+   const bool enableBubbleSplits             = modelParameters.getParameter< bool >("enableBubbleSplits");
+   const real_t cellConversionThreshold      = modelParameters.getParameter< real_t >("cellConversionThreshold");
+   const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableForceWeighting);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
+
+   // read evaluation parameters from parameter file
+   const auto evaluationParameters      = walberlaEnv.config()->getOneBlock("EvaluationParameters");
+   const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+
+   // create non-uniform block forest (non-uniformity required for load balancing)
+   const std::shared_ptr< StructuredBlockForest > blockForest =
+      createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
+
+   // add force field
+   const BlockDataID forceFieldID =
+      field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
+
+   // create lattice model
+   const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceFieldID));
+
+   // add pdf field
+   const BlockDataID pdfFieldID =
+      lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, impactVelocity, real_c(1), field::fzyx);
+
+   // add fill level field (initialized with 0, i.e., gas everywhere)
+   const BlockDataID fillFieldID =
+      field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(0.0), field::fzyx, uint_c(2));
+
+   // add boundary handling
+   const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
+      std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
+   const BlockDataID flagFieldID                                      = freeSurfaceBoundaryHandling->getFlagFieldID();
+   const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
+
+   const real_t poolHeight = poolHeightFactor * dropDiameter;
+
+   // initialize a pool of liquid at the bottom of the domain in z-direction
+   const AABB poolAABB(real_c(0), real_c(0), real_c(0), real_c(domainSize[0]), real_c(domainSize[1]), poolHeight);
+
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      ScalarField_T* const fillField = blockIt->getData< ScalarField_T >(fillFieldID);
+      PdfField_T* const pdfField     = blockIt->getData< PdfField_T >(pdfFieldID);
+
+      // determine the cells that are relevant for a block (still in global coordinates)
+      const CellInterval relevantCellBB = blockForest->getCellBBFromAABB(poolAABB.getIntersection(blockIt->getAABB()));
+
+      // transform the global coordinates of relevant cells to block local coordinates
+      CellInterval blockLocalCellBB;
+      blockForest->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, *blockIt, relevantCellBB);
+
+      WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(blockLocalCellBB, {
+         fillField->get(x, y, z) = real_c(1);
+         pdfField->setDensityAndVelocity(x, y, z, Vector3< real_t >(real_c(0)), real_c(1));
+      }) // WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ
+   }
+
+   const Vector3< real_t > dropCenter = dropCenterFactor * dropDiameter;
+
+   const geometry::Sphere sphereDrop(dropCenter, dropDiameter * real_c(0.5));
+   bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereDrop, false);
+
+   // initialize boundary conditions from config file
+   const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
+   freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
+
+   // IMPORTANT REMARK: this must be only called after every solid flag has been set; otherwise, the boundary handling
+   // might not detect solid flags correctly
+   freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
+
+   // communication after initialization
+   Communication_T communication(blockForest, flagFieldID, fillFieldID, forceFieldID);
+   communication();
+
+   PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
+   pdfCommunication();
+
+   // add bubble model
+   std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel = nullptr;
+   if (enableBubbleModel)
+   {
+      const std::shared_ptr< bubble_model::BubbleModel< CommunicationStencil_T > > bubbleModelDerived =
+         std::make_shared< bubble_model::BubbleModel< CommunicationStencil_T > >(blockForest, enableBubbleSplits);
+      bubbleModelDerived->initFromFillLevelField(fillFieldID);
+      bubbleModelDerived->setAtmosphere(
+         Cell(domainSize[0] - uint_c(1), domainSize[1] - uint_c(1), domainSize[2] - uint_c(1)), real_c(1));
+
+      bubbleModel = std::static_pointer_cast< bubble_model::BubbleModelBase >(bubbleModelDerived);
+   }
+   else { bubbleModel = std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1)); }
+
+   // initialize hydrostatic pressure
+   initHydrostaticPressure< PdfField_T >(blockForest, pdfFieldID, force, poolHeight);
+
+   // set density in non-liquid or non-interface cells to 1 (after initializing with hydrostatic pressure)
+   setDensityInNonFluidCellsToOne< FlagField_T, PdfField_T >(blockForest, flagInfo, flagFieldID, pdfFieldID);
+
+   // create timeloop
+   SweepTimeloop timeloop(blockForest, timesteps);
+
+   // Laplace pressure = 2 * surface tension * curvature; curvature computation is not necessary with 0 surface tension
+   bool computeCurvature = false;
+   if (!realIsEqual(surfaceTension, real_c(0), real_c(1e-14))) { computeCurvature = true; }
+
+   // add surface geometry handler
+   const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
+      blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, computeCurvature, enableWetting,
+      contactAngle);
+
+   geometryHandler.addSweeps(timeloop);
+
+   const ConstBlockDataID curvatureFieldID = geometryHandler.getConstCurvatureFieldID();
+   const ConstBlockDataID normalFieldID    = geometryHandler.getConstNormalFieldID();
+
+   // add boundary handling for standard boundaries and free surface boundaries
+   const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
+      blockForest, pdfFieldID, flagFieldID, fillFieldID, forceFieldID, normalFieldID, curvatureFieldID,
+      freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
+      relaxationRate, force, surfaceTension, enableForceWeighting, useSimpleMassExchange, cellConversionThreshold,
+      cellConversionForceThreshold);
+
+   dynamicsHandler.addSweeps(timeloop);
+
+   // add load balancing
+   const LoadBalancer< FlagField_T, CommunicationStencil_T, LatticeModelStencil_T > loadBalancer(
+      blockForest, communication, pdfCommunication, bubbleModel, uint_c(50), uint_c(10), uint_c(5),
+      loadBalancingFrequency, printLoadBalancingStatistics);
+   timeloop.addFuncAfterTimeStep(loadBalancer, "Sweep: load balancing");
+
+   // add VTK output
+   addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
+      blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceFieldID,
+      geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
+      geometryHandler.getObstNormalFieldID());
+
+   // add triangle mesh output of free surface
+   SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
+      blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
+   surfaceMeshWriter(); // write initial mesh
+   timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
+
+   // add logging for computational performance
+   const lbm::PerformanceLogger< FlagField_T > perfLogger(blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs,
+                                                          performanceLogFrequency);
+   timeloop.addFuncAfterTimeStep(perfLogger, "Evaluator: performance logging");
+
+   WcTimingPool timingPool;
+
+   for (uint_t t = uint_c(0); t != timesteps; ++t)
+   {
+      if (t % uint_c(real_c(timesteps / 100)) == uint_c(0))
+      {
+         WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep = " << t);
+      }
+      timeloop.singleStep(timingPool, true);
+
+      if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
+   }
+
+   return EXIT_SUCCESS;
+}
+
+} // namespace DropImpact
+} // namespace free_surface
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::free_surface::DropImpact::main(argc, argv); }
\ No newline at end of file

apps/showcases/FreeSurface/DropImpact.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 10, 10, 10 >;
+   periodicity                   < 1, 1, 0 >;
+   loadBalancingFrequency        372;
+   printLoadBalancingStatistics  true;
+}
+
+DomainParameters
+{
+   dropDiameter      20;
+   dropCenterFactor  < 2, 2, 1 >;      // values multiplied with dropDiameter
+   poolHeightFactor  0.5;              // value multiplied with dropDiameter
+   domainSizeFactor  < 10, 10, 5 >;    // values multiplied with dropDiameter
+}
+
+PhysicsParameters
+{
+   bondNumber        3.18;
+   weberNumber       2010;
+   ohnesorgeNumber   0.0384;
+   relaxationRate    1.989;
+   impactAngleDegree 0;
+   enableWetting     false;
+   contactAngle      0; // only used if enableWetting=true
+   timesteps         6000;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   enableBubbleModel             false;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 3000;
+}
+
+BoundaryParameters
+{
+   // X
+   //Border { direction W;  walldistance -1; NoSlip{} }
+   //Border { direction E;  walldistance -1; NoSlip{} }
+
+   // Y
+   //Border { direction N;  walldistance -1; NoSlip{} }
+   //Border { direction S;  walldistance -1; NoSlip{} }
+
+   // Z
+   Border { direction T;  walldistance -1; NoSlip{} }
+   Border { direction B;  walldistance -1; NoSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 50;
+   baseFolder mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency 1000;
+      ghostLayers 0;
+      baseFolder vtk-out;
+      samplingResolution 1;
+
+      writers
+      {
+         velocity;
+         density;
+         //pdf;
+         flag;
+         fill_level;
+         force;
+         curvature;
+         normal;
+         obstacle_normal;
+         mapped_flag;
+        }
+
+        inclusion_filters
+        {
+            // only include liquid and interface cells in VTK output
+            //liquidInterfaceFilter;
+        }
+
+        before_functions
+        {
+            //ghost_layer_synchronization; // only needed if writing the ghost layer
+        }
+
+   }
+   domain_decomposition
+   {
+      writeFrequency 372;
+      baseFolder vtk-out;
+      outputDomainDecomposition true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/DropWetting.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 10, 10, 10 >;
+   periodicity                   < 1, 1, 0 >;
+   loadBalancingFrequency        100;
+   printLoadBalancingStatistics  true;
+}
+
+DomainParameters
+{
+   dropDiameter      20;
+   domainSizeFactor  < 3, 3, 2 >;    // values multiplied with dropDiameter
+}
+
+PhysicsParameters
+{
+   relaxationRate 1.96;
+   surfaceTension 8.37e-3;
+   force <0, 0, -1.14e-7>;
+   enableWetting true;
+   contactAngle 45; // only used if enableWetting=true
+   timesteps       10001;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   enableBubbleModel             false;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 1000;
+   evaluationFrequency 200;
+   convergenceThreshold 1e-5;
+}
+
+BoundaryParameters
+{
+   // X
+   //Border { direction W;  walldistance -1; NoSlip{} }
+   //Border { direction E;  walldistance -1; NoSlip{} }
+
+   // Y
+   //Border { direction N;  walldistance -1; NoSlip{} }
+   //Border { direction S;  walldistance -1; NoSlip{} }
+
+   // Z
+   Border { direction T;  walldistance -1; NoSlip{} }
+   Border { direction B;  walldistance -1; NoSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 100;
+   baseFolder mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency 100;
+      ghostLayers 0;
+      baseFolder vtk-out;
+      samplingResolution 1;
+
+      writers
+      {
+         velocity;
+         density;
+         //pdf;
+         flag;
+         fill_level;
+         force;
+         curvature;
+         normal;
+         obstacle_normal;
+         mapped_flag;
+        }
+
+        inclusion_filters
+        {
+            // only include liquid and interface cells in VTK output
+            //liquidInterfaceFilter;
+        }
+
+        before_functions
+        {
+            //ghost_layer_synchronization; // only needed if writing the ghost layer
+        }
+
+   }
+   domain_decomposition
+   {
+      writeFrequency 100;
+      baseFolder vtk-out;
+      outputDomainDecomposition true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/GravityWave.prm

+BlockForestParameters
+{
+   cellsPerBlock                 < 25, 25, 1 >;
+   periodicity                   < 1, 0, 1 >;
+   loadBalancingFrequency        0;
+   printLoadBalancingStatistics  false;
+}
+
+DomainParameters
+{
+   domainWidth       200; // equivalent to wavelength
+   liquidDepth       100;
+   initialAmplitude  2;
+}
+
+PhysicsParameters
+{
+   reynoldsNumber    10;
+   relaxationRate    1.8;
+   enableWetting     false;
+   contactAngle      0; // only used if enableWetting=true
+   timesteps         86400;
+}
+
+ModelParameters
+{
+   pdfReconstructionModel        OnlyMissing;
+   pdfRefillingModel             EquilibriumRefilling;
+   excessMassDistributionModel   EvenlyAllInterface;
+   curvatureModel                FiniteDifferenceMethod;
+   enableForceWeighting          false;
+   useSimpleMassExchange         false;
+   cellConversionThreshold       1e-2;
+   cellConversionForceThreshold  1e-1;
+
+   enableBubbleModel             false;
+   enableBubbleSplits            false; // only used if enableBubbleModel=true
+}
+
+EvaluationParameters
+{
+   performanceLogFrequency 43200;
+   evaluationFrequency     864;
+   filename                gravity-wave.txt;
+}
+
+BoundaryParameters
+{
+   // X
+   //Border { direction W;  walldistance -1; NoSlip{} }
+   //Border { direction E;  walldistance -1; NoSlip{} }
+
+   // Y
+   Border { direction N;  walldistance -1; NoSlip{} }
+   Border { direction S;  walldistance -1; NoSlip{} }
+
+   // Z
+   //Border { direction T;  walldistance -1; NoSlip{} }
+   //Border { direction B;  walldistance -1; NoSlip{} }
+}
+
+MeshOutputParameters
+{
+   writeFrequency 864;
+   baseFolder     mesh-out;
+}
+
+VTK
+{
+   fluid_field
+   {
+      writeFrequency       864;
+      ghostLayers          0;
+      baseFolder           vtk-out;
+      samplingResolution   1;
+
+      writers
+      {
+         fill_level;
+         mapped_flag;
+         velocity;
+         density;
+         //curvature;
+         //normal;
+         //obstacle_normal;
+         //pdf;
+         //flag;
+         //force;
+      }
+
+      inclusion_filters
+      {
+         //liquidInterfaceFilter; // only include liquid and interface cells in VTK output
+      }
+
+      before_functions
+      {
+         //ghost_layer_synchronization; // only needed if writing the ghost layer
+      }
+   }
+   domain_decomposition
+   {
+      writeFrequency             0;
+      baseFolder                 vtk-out;
+      outputDomainDecomposition  true;
+   }
+}
\ No newline at end of file

apps/showcases/FreeSurface/GravityWaveLatticeModelGeneration.py

+import sympy as sp
+
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, create_lb_collision_rule
+from lbmpy.enums import ForceModel, Method, Stencil
+from lbmpy.stencils import LBStencil
+
+from pystencils_walberla import CodeGeneration
+from lbmpy_walberla import generate_lattice_model
+
+# general parameters
+stencil = LBStencil(Stencil.D3Q19)
+omega = sp.Symbol('omega')
+force = sp.symbols('force_:3')
+layout = 'fzyx'
+
+# method definition
+lbm_config = LBMConfig(stencil=stencil,
+                       method=Method.SRT,
+                       relaxation_rate=omega,
+                       compressible=True,
+                       force=force,
+                       force_model=ForceModel.GUO,
+                       zero_centered=False,
+                       streaming_pattern='pull')  # free surface implementation only works with pull pattern
+
+# optimizations to be used by the code generator
+lbm_opt = LBMOptimisation(cse_global=True,
+                          field_layout=layout)
+
+collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
+                                          lbm_optimisation=lbm_opt)
+
+with CodeGeneration() as ctx:
+    generate_lattice_model(ctx, "GravityWaveLatticeModel", collision_rule, field_layout=layout)