First test implementation

apps/showcases/FreeSurface/CMakeLists.txt

@@ -29,14 +29,17 @@ waLBerla_add_executable(NAME    GravityWave
                         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_generate_target_from_python( NAME      GravityWaveGeneration
+                                         FILE      GravityWave.py
+                                         OUT_FILES GravityWaveStorageSpecification.h GravityWaveStorageSpecification.${CODEGEN_FILE_SUFFIX}
+                                                   GravityWaveSweepCollection.h GravityWaveSweepCollection.${CODEGEN_FILE_SUFFIX}
+                                                   NoSlip.h NoSlip.${CODEGEN_FILE_SUFFIX}
+                                                   GravityWaveBoundaryCollection.h
+                                                   GravityWaveInfoHeader.h)
 
    waLBerla_add_executable(NAME    GravityWaveCodegen
                            FILES   GravityWaveCodegen.cpp
-                           DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk
-                                   GravityWaveLatticeModelGeneration)
+                           DEPENDS blockforest boundary core domain_decomposition field geometry lbm_generated postprocessing timeloop vtk GravityWaveGeneration)
 endif()
 
 waLBerla_add_executable(NAME    MovingDrop

apps/showcases/FreeSurface/GravityWave.prm

@@ -43,6 +43,11 @@ EvaluationParameters
    filename                gravity-wave.txt;
 }
 
+Logging
+{
+    logLevel info;  // info progress detail tracing
+}
+
 BoundaryParameters
 {
    // X

apps/showcases/FreeSurface/GravityWaveLatticeModelGeneration.py → apps/showcases/FreeSurface/GravityWave.py

 import sympy as sp
 import pystencils as ps
 from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, create_lb_collision_rule
+from lbmpy.boundaries import NoSlip
 from lbmpy.enums import ForceModel, Method, Stencil
 from lbmpy.stencils import LBStencil
 
-from pystencils_walberla import CodeGeneration
-from lbmpy_walberla import generate_lattice_model
+from pystencils_walberla import CodeGeneration, generate_info_header, generate_sweep
+from lbmpy_walberla import generate_lbm_package, lbm_boundary_generator
+
+
+info_header = """
+const bool infoCseGlobal = {cse_global};
+const bool infoCsePdfs = {cse_pdfs};
+"""
 
 
 with CodeGeneration() as ctx:
@@ -15,6 +22,11 @@ with CodeGeneration() as ctx:
 
     stencil = LBStencil(Stencil.D3Q19)
     omega = sp.Symbol('omega')
+
+    assert stencil.D == 3, "This application supports only three-dimensional stencils"
+    pdfs, pdfs_tmp = ps.fields(f"pdfs({stencil.Q}), pdfs_tmp({stencil.Q}): {data_type}[3D]", layout='fzyx')
+    density_field, velocity_field = ps.fields(f"density, velocity(3) : {data_type}[3D]", layout='fzyx')
+    macroscopic_fields = {'density': density_field, 'velocity': velocity_field}
     force_field = ps.fields(f"force(3): {data_type}[3D]", layout='fzyx')
 
     # method definition
@@ -29,9 +41,30 @@ with CodeGeneration() as ctx:
 
     # optimizations to be used by the code generator
     lbm_opt = LBMOptimisation(cse_global=True,
-                              field_layout=layout)
+                              field_layout=layout, symbolic_field=pdfs, symbolic_temporary_field=pdfs_tmp)
 
     collision_rule = create_lb_collision_rule(lbm_config=lbm_config,
                                               lbm_optimisation=lbm_opt)
 
-    generate_lattice_model(ctx, "GravityWaveLatticeModel", collision_rule, field_layout=layout)
+    no_slip = lbm_boundary_generator(class_name='NoSlip', flag_uid='NoSlip',
+                                     boundary_object=NoSlip())
+
+    generate_lbm_package(ctx, name="GravityWave",
+                         collision_rule=collision_rule,
+                         lbm_config=lbm_config, lbm_optimisation=lbm_opt,
+                         nonuniform=False, boundaries=[no_slip, ],
+                         macroscopic_fields=macroscopic_fields,
+                         target=ps.Target.CPU)
+
+    infoHeaderParams = {
+        'cse_global': int(lbm_opt.cse_global),
+        'cse_pdfs': int(lbm_opt.cse_pdfs),
+    }
+
+    field_typedefs = {'VectorField_T': velocity_field,
+                      'ScalarField_T': density_field}
+
+    # Info header containing correct template definitions for stencil and field
+    generate_info_header(ctx, 'GravityWaveInfoHeader',
+                         field_typedefs=field_typedefs,
+                         additional_code=info_header.format(**infoHeaderParams))

src/lbm_generated/CMakeLists.txt

@@ -15,11 +15,14 @@ target_link_libraries( lbm_generated
         vtk
         )
 
+add_subdirectory( blockforest )
 add_subdirectory( boundary )
 add_subdirectory( communication )
 add_subdirectory( gpu )
+add_subdirectory( macroscopics )
 add_subdirectory( evaluation )
 add_subdirectory( field )
+add_subdirectory( free_surface )
 add_subdirectory( refinement )
 add_subdirectory( storage_specification )
 add_subdirectory( sweep_collection )
\ No newline at end of file

src/lbm_generated/blockforest/CMakeLists.txt

+target_sources( lbm_generated
+        PRIVATE
+        SimpleCommunication.h
+        UpdateSecondGhostLayer.h
+        )

src/lbm_generated/blockforest/SimpleCommunication.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 SimpleCommunication.h
+//! \ingroup blockforest
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "core/Abort.h"
+#include "core/math/Vector3.h"
+#include "core/mpi/BufferDataTypeExtensions.h"
+
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+
+namespace walberla::lbm_generated
+{
+using blockforest::communication::UniformBufferedScheme;
+using field::communication::PackInfo;
+
+template< typename Stencil_T >
+class SimpleCommunication : public blockforest::communication::UniformBufferedScheme< Stencil_T >
+{
+   using RealScalarField_T      = GhostLayerField< real_t, 1 >;
+   using VectorField_T          = GhostLayerField< Vector3< real_t >, 1 >;
+   using VectorFieldFlattened_T = GhostLayerField< real_t, 3 >;
+   using UintScalarField_T      = GhostLayerField< uint_t, 1 >;
+   using IDScalarField_T        = walberla::GhostLayerField< walberla::id_t, 1 >;
+
+   using FlagField16_T = FlagField< uint16_t >;
+   using FlagField32_T = FlagField< uint32_t >;
+   using FlagField64_T = FlagField< uint64_t >;
+
+ public:
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1;
+   }
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2;
+   }
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2,
+                       BlockDataID f3)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2 << f3;
+   }
+
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2,
+                       BlockDataID f3, BlockDataID f4)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2 << f3 << f4;
+   }
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2,
+                       BlockDataID f3, BlockDataID f4, BlockDataID f5)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2 << f3 << f4 << f5;
+   }
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2,
+                       BlockDataID f3, BlockDataID f4, BlockDataID f5, BlockDataID f6)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2 << f3 << f4 << f5 << f6;
+   }
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2,
+                       BlockDataID f3, BlockDataID f4, BlockDataID f5, BlockDataID f6, BlockDataID f7)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2 << f3 << f4 << f5 << f6 << f7;
+   }
+   SimpleCommunication(const std::weak_ptr< StructuredBlockForest >& blockForest, BlockDataID f1, BlockDataID f2,
+                       BlockDataID f3, BlockDataID f4, BlockDataID f5, BlockDataID f6, BlockDataID f7, BlockDataID f8)
+      : UniformBufferedScheme< Stencil_T >(blockForest), blockForest_(blockForest)
+   {
+      (*this) << f1 << f2 << f3 << f4 << f5 << f6 << f7 << f8;
+   }
+
+   SimpleCommunication& operator<<(BlockDataID fieldId)
+   {
+      auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      if (blockForest->begin() == blockForest->end()) { return *this; }
+      IBlock& firstBlock = *(blockForest->begin());
+
+      if (firstBlock.isDataClassOrSubclassOf< RealScalarField_T >(fieldId))
+      {
+         this->addPackInfo(make_shared< PackInfo< RealScalarField_T > >(fieldId));
+      }
+      else
+      {
+         if (firstBlock.isDataClassOrSubclassOf< VectorField_T >(fieldId))
+         {
+            this->addPackInfo(make_shared< PackInfo< VectorField_T > >(fieldId));
+         }
+         else
+         {
+            if (firstBlock.isDataClassOrSubclassOf< FlagField16_T >(fieldId))
+            {
+               this->addPackInfo(make_shared< PackInfo< FlagField16_T > >(fieldId));
+            }
+            else
+            {
+               if (firstBlock.isDataClassOrSubclassOf< FlagField32_T >(fieldId))
+               {
+                  this->addPackInfo(make_shared< PackInfo< FlagField32_T > >(fieldId));
+               }
+               else
+               {
+                  if (firstBlock.isDataClassOrSubclassOf< FlagField64_T >(fieldId))
+                  {
+                     this->addPackInfo(make_shared< PackInfo< FlagField64_T > >(fieldId));
+                  }
+                  else
+                  {
+                     if (firstBlock.isDataClassOrSubclassOf< IDScalarField_T >(fieldId))
+                     {
+                        this->addPackInfo(make_shared< PackInfo< IDScalarField_T > >(fieldId));
+                     }
+                     else
+                     {
+                        if (firstBlock.isDataClassOrSubclassOf< UintScalarField_T >(fieldId))
+                        {
+                           this->addPackInfo(make_shared< PackInfo< UintScalarField_T > >(fieldId));
+                        }
+                        else
+                        {
+                           if (firstBlock.isDataClassOrSubclassOf< VectorFieldFlattened_T >(fieldId))
+                           {
+                              this->addPackInfo(make_shared< PackInfo< VectorFieldFlattened_T > >(fieldId));
+                           }
+                           else { WALBERLA_ABORT("Problem with UID"); }
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+
+      return *this;
+   }
+
+ protected:
+   std::weak_ptr< StructuredBlockForest > blockForest_;
+}; // class SimpleCommunication
+
+} // namespace walberla

src/lbm_generated/blockforest/UpdateSecondGhostLayer.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 UpdateSecondGhostLayer.h
+//! \ingroup blockforest
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+
+namespace walberla::lbm_generated
+{
+/********************************************************************************************************************
+ * Manual update of a field's second and further ghost layers in setups with domain size = 1 and periodicity in at
+ * least one direction.
+ *
+ * If a field has two ghost layers and if the inner field has only a size of one (in one or more directions),
+ * regular communication does not update the second (and further) ghost layers correctly. Here, the content of the
+ * first ghost layers is manually copied to the second (and further) ghost layers when the dimension of the field is
+ * one in this direction.
+ *******************************************************************************************************************/
+template< typename Field_T >
+class UpdateSecondGhostLayer
+{
+ public:
+   UpdateSecondGhostLayer(const std::weak_ptr< StructuredBlockForest >& blockForestPtr, BlockDataID fieldID)
+      : blockForest_(blockForestPtr), fieldID_(fieldID)
+   {
+      const auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         Field_T* const field = blockIt->template getData< Field_T >(fieldID_);
+
+         // this function is only necessary if the flag field has at least two ghost layers
+         if (field->nrOfGhostLayers() < uint_c(2)) { isNecessary_ = false; }
+         else
+         {
+            // running this function is only necessary when the field is of size 1 in at least one direction
+            isNecessary_ = (field->xSize() == uint_c(1) && blockForest->isXPeriodic()) ||
+                           (field->ySize() == uint_c(1) && blockForest->isYPeriodic()) ||
+                           (field->zSize() == uint_c(1) && blockForest->isZPeriodic());
+         }
+      }
+   }
+
+   void operator()()
+   {
+      if (!isNecessary_) { return; }
+
+      const auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         Field_T* const field = blockIt->template getData< Field_T >(fieldID_);
+
+         if (field->xSize() == uint_c(1) && blockForest->isXPeriodic())
+         {
+            // iterate ghost layer at x == -1
+            for (auto fieldIt = field->beginGhostLayerOnly(uint_c(1), stencil::W, true); fieldIt != field->end();
+                 ++fieldIt)
+            {
+               // copy data from ghost layer at x == -1 to x == -2
+               fieldIt.neighbor(stencil::W) = *fieldIt;
+            }
+
+            // iterate ghost layer at x == xSize()+1
+            for (auto fieldIt = field->beginGhostLayerOnly(uint_c(1), stencil::E, true); fieldIt != field->end();
+                 ++fieldIt)
+            {
+               // copy data from ghost layer at x == xSize()+1 to x == xSize()+2
+               fieldIt.neighbor(stencil::E) = *fieldIt;
+            }
+         }
+
+         if (field->ySize() == uint_c(1) && blockForest->isYPeriodic())
+         {
+            // iterate ghost layer at y == -1
+            for (auto fieldIt = field->beginGhostLayerOnly(uint_c(1), stencil::S, true); fieldIt != field->end();
+                 ++fieldIt)
+            {
+               // copy data from ghost layer at y == -1 to y == -2
+               fieldIt.neighbor(stencil::S) = *fieldIt;
+            }
+
+            // iterate ghost layer at y == ySize()+1
+            for (auto fieldIt = field->beginGhostLayerOnly(uint_c(1), stencil::N, true); fieldIt != field->end();
+                 ++fieldIt)
+            {
+               // copy data from ghost layer at y == ySize()+1 to y == ySize()+2
+               fieldIt.neighbor(stencil::N) = *fieldIt;
+            }
+         }
+
+         if (field->zSize() == uint_c(1) && blockForest->isZPeriodic())
+         {
+            // iterate ghost layer at z == -1
+            for (auto fieldIt = field->beginGhostLayerOnly(uint_c(1), stencil::B, true); fieldIt != field->end();
+                 ++fieldIt)
+            {
+               // copy data from ghost layer at z == -1 to z == -2
+               fieldIt.neighbor(stencil::B) = *fieldIt;
+            }
+
+            // iterate ghost layer at y == zSize()+1
+            for (auto fieldIt = field->beginGhostLayerOnly(uint_c(1), stencil::T, true); fieldIt != field->end();
+                 ++fieldIt)
+            {
+               // copy data from ghost layer at z == zSize()+1 to z == zSize()+2
+               fieldIt.neighbor(stencil::T) = *fieldIt;
+            }
+         }
+      }
+   }
+
+ private:
+   std::weak_ptr< StructuredBlockForest > blockForest_;
+
+   BlockDataID fieldID_;
+
+   bool isNecessary_;
+}; // class UpdateSecondGhostLayer
+
+} // namespace walberla

src/lbm_generated/free_surface/BlockStateDetectorSweep.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 BlockStateDetectorSweep.h
+//! \ingroup free_surface
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Set block states according to their content, i.e., free surface flags.
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/uid/SUID.h"
+
+#include "FlagInfo.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+/***********************************************************************************************************************
+ * Set block states according to free surface flag field:
+ *  - blocks that contain at least one interface cell are marked as "fullFreeSurface" (computationally most expensive
+ *    type of blocks)
+ *  - cells without interface cell and with at least one liquid cell are marked "onlyLBM"
+ *  - all other blocks are marked as "onlyGasAndBoundary"
+ **********************************************************************************************************************/
+template< typename FlagField_T >
+class BlockStateDetectorSweep
+{
+ public:
+   static const SUID fullFreeSurface;
+   static const SUID onlyGasAndBoundary;
+   static const SUID onlyLBM;
+
+   BlockStateDetectorSweep(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+                           FlagInfo< FlagField_T > flagInfo, ConstBlockDataID flagFieldID)
+      : flagFieldID_(flagFieldID), flagInfo_(flagInfo)
+   {
+      const auto blockForest = blockForestPtr.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         (*this)(&(*blockIt)); // initialize block states
+      }
+   }
+
+   void operator()(IBlock* const block)
+   {
+      bool liquidFound    = false;
+      bool interfaceFound = false;
+
+      const FlagField_T* const flagField = block->getData< const FlagField_T >(flagFieldID_);
+
+      // search the flag field for interface and liquid cells
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(flagField, uint_c(1), {
+         const typename FlagField_T::ConstPtr flagFieldPtr(*flagField, x, y, z);
+
+         // at inflow boundaries, interface cells are generated; therefore, the block must be fullFreeSurface even if it
+         // does not yet contain an interface cell
+         if (flagInfo_.isInterface(flagFieldPtr) || flagInfo_.isInflow(flagFieldPtr))
+         {
+            interfaceFound = true;
+            break; // stop search, as this block belongs already to the computationally most expensive block type
+         }
+
+         if (flagInfo_.isLiquid(flagFieldPtr)) { liquidFound = true; }
+      }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+
+      block->clearState();
+      if (interfaceFound)
+      {
+         block->addState(fullFreeSurface);
+         return;
+      }
+      if (liquidFound)
+      {
+         block->addState(onlyLBM);
+         return;
+      }
+      block->addState(onlyGasAndBoundary);
+   }
+
+ protected:
+   ConstBlockDataID flagFieldID_;
+   FlagInfo< FlagField_T > flagInfo_;
+}; // class BlockStateDetectorSweep
+
+template< typename FlagField_T >
+const SUID BlockStateDetectorSweep< FlagField_T >::fullFreeSurface = SUID("fullFreeSurface");
+template< typename FlagField_T >
+const SUID BlockStateDetectorSweep< FlagField_T >::onlyGasAndBoundary = SUID("onlyGasAndBoundary");
+template< typename FlagField_T >
+const SUID BlockStateDetectorSweep< FlagField_T >::onlyLBM = SUID("onlyLBM");
+
+} // namespace free_surface
+} // namespace walberla

src/lbm_generated/free_surface/CMakeLists.txt

+target_sources( lbm_generated
+        PRIVATE
+        BlockStateDetectorSweep.h
+        FlagDefinitions.h
+        FlagInfo.h
+        FlagInfo.impl.h
+        InitFunctions.h
+        InterfaceFromFillLevel.h
+        LoadBalancing.h
+        MaxVelocityComputer.h
+        SurfaceMeshWriter.h
+        TotalMassComputer.h
+        VtkWriter.h
+        )
+
+add_subdirectory( dynamics )
+add_subdirectory( surface_geometry )

src/lbm_generated/free_surface/FlagDefinitions.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 FlagInfo.h
+//! \ingroup free_surface
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Define free surface flags (e.g. conversion flags).
+//
+//======================================================================================================================
+
+#include "field/FlagUID.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+namespace flagIDs
+{
+/***********************************************************************************************************************
+ * Definition of free surface flag IDs.
+ **********************************************************************************************************************/
+const field::FlagUID interfaceFlagID                   = field::FlagUID("interface");
+const field::FlagUID liquidFlagID                      = field::FlagUID("liquid");
+const field::FlagUID gasFlagID                         = field::FlagUID("gas");
+const field::FlagUID convertedFlagID                   = field::FlagUID("converted");
+const field::FlagUID convertToGasFlagID                = field::FlagUID("convert to gas");
+const field::FlagUID convertToLiquidFlagID             = field::FlagUID("convert to liquid");
+const field::FlagUID convertedFromGasToInterfaceFlagID = field::FlagUID("convert from gas to interface");
+const field::FlagUID keepInterfaceForWettingFlagID     = field::FlagUID("convert to and keep interface for wetting");
+const field::FlagUID convertToInterfaceForInflowFlagID = field::FlagUID("convert to interface for inflow");
+
+const Set< field::FlagUID > liquidInterfaceFlagIDs = setUnion< field::FlagUID >(liquidFlagID, interfaceFlagID);
+
+const Set< field::FlagUID > liquidInterfaceGasFlagIDs =
+   setUnion(liquidInterfaceFlagIDs, Set< field::FlagUID >(gasFlagID));
+
+} // namespace flagIDs
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file

src/lbm_generated/free_surface/FlagInfo.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 FlagInfo.h
+//! \ingroup free_surface
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Manage free surface flags (e.g. conversion flags).
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/mpi/Broadcast.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/FlagField.h"
+
+#include "FlagDefinitions.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+/***********************************************************************************************************************
+ * Manage free surface flags (e.g. conversion flags).
+ **********************************************************************************************************************/
+template< typename FlagField_T >
+class FlagInfo
+{
+ public:
+   using flag_t = typename FlagField_T::flag_t;
+
+   FlagInfo();
+   FlagInfo(const Set< field::FlagUID >& obstacleIDs,
+            const Set< field::FlagUID >& outflowIDs  = Set< field::FlagUID >::emptySet(),
+            const Set< field::FlagUID >& inflowIDs   = Set< field::FlagUID >::emptySet(),
+            const Set< field::FlagUID >& freeSlipIDs = Set< field::FlagUID >::emptySet());
+
+   bool operator==(const FlagInfo& o) const;
+   bool operator!=(const FlagInfo& o) const;
+
+   flag_t interfaceFlag;
+   flag_t liquidFlag;
+   flag_t gasFlag;
+
+   flag_t convertedFlag;                   // interface cell that is already converted
+   flag_t convertToGasFlag;                // interface cell with too low fill level, should be converted
+   flag_t convertToLiquidFlag;             // interface cell with too high fill level, should be converted
+   flag_t convertFromGasToInterfaceFlag;   // interface cell that was a gas cell and needs refilling of its pdfs
+   flag_t keepInterfaceForWettingFlag;     // gas/liquid cell that needs to be converted to interface cell for a smooth
+                                           // continuation of the wetting surface (see dissertation of S. Donath, 2011,
+                                           // section 6.3.5.3)
+   flag_t convertToInterfaceForInflowFlag; // gas cell that needs to be converted to interface cell to enable inflow
+
+   flag_t obstacleFlagMask;
+   flag_t outflowFlagMask;
+   flag_t inflowFlagMask;
+   flag_t freeSlipFlagMask; // free slip obstacle cell (needs to be treated separately since PDFs going from gas into
+                            // boundary are not available and must be reconstructed)
+
+   template< typename FieldItOrPtr_T >
+   inline bool isInterface(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, interfaceFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isLiquid(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, liquidFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isGas(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, gasFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isObstacle(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isPartOfMaskSet(flagItOrPtr, obstacleFlagMask);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isOutflow(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isPartOfMaskSet(flagItOrPtr, outflowFlagMask);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isInflow(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isPartOfMaskSet(flagItOrPtr, inflowFlagMask);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isFreeSlip(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isPartOfMaskSet(flagItOrPtr, freeSlipFlagMask);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool hasConverted(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, convertedFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool hasConvertedToGas(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, convertToGasFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool hasConvertedToLiquid(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, convertToLiquidFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool hasConvertedFromGasToInterface(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, convertFromGasToInterfaceFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isKeepInterfaceForWetting(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, keepInterfaceForWettingFlag);
+   }
+   template< typename FieldItOrPtr_T >
+   inline bool isConvertToInterfaceForInflow(const FieldItOrPtr_T& flagItOrPtr) const
+   {
+      return isFlagSet(flagItOrPtr, convertToInterfaceForInflowFlag);
+   }
+
+   inline bool isInterface(const flag_t val) const { return isFlagSet(val, interfaceFlag); }
+   inline bool isLiquid(const flag_t val) const { return isFlagSet(val, liquidFlag); }
+   inline bool isGas(const flag_t val) const { return isFlagSet(val, gasFlag); }
+   inline bool isObstacle(const flag_t val) const { return isPartOfMaskSet(val, obstacleFlagMask); }
+   inline bool isOutflow(const flag_t val) const { return isPartOfMaskSet(val, outflowFlagMask); }
+   inline bool isInflow(const flag_t val) const { return isPartOfMaskSet(val, inflowFlagMask); }
+   inline bool isFreeSlip(const flag_t val) const { return isPartOfMaskSet(val, freeSlipFlagMask); }
+   inline bool isKeepInterfaceForWetting(const flag_t val) const { return isFlagSet(val, keepInterfaceForWettingFlag); }
+   inline bool hasConvertedToGas(const flag_t val) const { return isFlagSet(val, convertToGasFlag); }
+   inline bool hasConvertedToLiquid(const flag_t val) const { return isFlagSet(val, convertToLiquidFlag); }
+   inline bool hasConvertedFromGasToInterface(const flag_t val) const
+   {
+      return isFlagSet(val, convertFromGasToInterfaceFlag);
+   }
+   inline bool isConvertToInterfaceForInflow(const flag_t val) const
+   {
+      return isFlagSet(val, convertToInterfaceForInflowFlag);
+   }
+
+   // check whether FlagInfo is identical on all blocks and all processes
+   bool isConsistentAcrossBlocksAndProcesses(const std::weak_ptr< StructuredBlockStorage >& blockStorage,
+                                             ConstBlockDataID flagFieldID) const;
+
+   // register flags in flag field
+   static void registerFlags(FlagField_T* field, const Set< field::FlagUID >& obstacleIDs,
+                             const Set< field::FlagUID >& outflowIDs  = Set< field::FlagUID >::emptySet(),
+                             const Set< field::FlagUID >& inflowIDs   = Set< field::FlagUID >::emptySet(),
+                             const Set< field::FlagUID >& freeSlipIDs = Set< field::FlagUID >::emptySet());
+
+   void registerFlags(FlagField_T* field) const;
+
+   Set< field::FlagUID > getObstacleIDSet() const { return obstacleIDs_; }
+   Set< field::FlagUID > getOutflowIDs() const { return outflowIDs_; }
+   Set< field::FlagUID > getInflowIDs() const { return inflowIDs_; }
+   Set< field::FlagUID > getFreeSlipIDs() const { return freeSlipIDs_; }
+
+ protected:
+   FlagInfo(const FlagField_T* field, const Set< field::FlagUID >& obstacleIDs, const Set< field::FlagUID >& outflowIDs,
+            const Set< field::FlagUID >& inflowIDs, const Set< field::FlagUID >& freeSlipIDs);
+
+   // create sets of flag IDs with flags from free_surface/boundary/FreeSurfaceBoundaryHandling.impl.h
+   Set< field::FlagUID > obstacleIDs_;
+   Set< field::FlagUID > outflowIDs_;
+   Set< field::FlagUID > inflowIDs_;
+   Set< field::FlagUID > freeSlipIDs_;
+};
+
+template< typename FlagField_T >
+mpi::SendBuffer& operator<<(mpi::SendBuffer& buf, const FlagInfo< FlagField_T >& flagInfo)
+{
+   buf << flagInfo.interfaceFlag << flagInfo.liquidFlag << flagInfo.gasFlag << flagInfo.convertedFlag
+       << flagInfo.convertToGasFlag << flagInfo.convertToLiquidFlag << flagInfo.convertFromGasToInterfaceFlag
+       << flagInfo.keepInterfaceForWettingFlag << flagInfo.keepInterfaceForWettingFlag
+       << flagInfo.convertToInterfaceForInflowFlag << flagInfo.obstacleFlagMask << flagInfo.outflowFlagMask
+       << flagInfo.inflowFlagMask << flagInfo.freeSlipFlagMask;
+
+   return buf;
+}
+
+template< typename FlagField_T >
+mpi::RecvBuffer& operator>>(mpi::RecvBuffer& buf, FlagInfo< FlagField_T >& flagInfo)
+{
+   buf >> flagInfo.interfaceFlag >> flagInfo.liquidFlag >> flagInfo.gasFlag >> flagInfo.convertedFlag >>
+      flagInfo.convertToGasFlag >> flagInfo.convertToLiquidFlag >> flagInfo.convertFromGasToInterfaceFlag >>
+      flagInfo.keepInterfaceForWettingFlag >> flagInfo.keepInterfaceForWettingFlag >>
+      flagInfo.convertToInterfaceForInflowFlag >> flagInfo.obstacleFlagMask >> flagInfo.outflowFlagMask >>
+      flagInfo.inflowFlagMask >> flagInfo.freeSlipFlagMask;
+
+   return buf;
+}
+
+} // namespace free_surface
+} // namespace walberla
+
+#include "FlagInfo.impl.h"
\ No newline at end of file

src/lbm_generated/free_surface/FlagInfo.impl.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 FlagInfo.impl.h
+//! \ingroup free_surface
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Define and manage free surface flags (e.g. conversion flags).
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "core/mpi/Broadcast.h"
+
+#include "field/FlagField.h"
+
+#include "FlagInfo.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+template< typename FlagField_T >
+FlagInfo< FlagField_T >::FlagInfo()
+   : interfaceFlag(0), liquidFlag(0), gasFlag(0), convertedFlag(0), convertToGasFlag(0), convertToLiquidFlag(0),
+     convertFromGasToInterfaceFlag(0), keepInterfaceForWettingFlag(0), convertToInterfaceForInflowFlag(0),
+     obstacleFlagMask(0), outflowFlagMask(0), inflowFlagMask(0), freeSlipFlagMask(0)
+{}
+
+template< typename FlagField_T >
+FlagInfo< FlagField_T >::FlagInfo(const FlagField_T* field, const Set< field::FlagUID >& obstacleIDs,
+                                  const Set< field::FlagUID >& outflowIDs, const Set< field::FlagUID >& inflowIDs,
+                                  const Set< field::FlagUID >& freeSlipIDs)
+   : interfaceFlag(field->getFlag(flagIDs::interfaceFlagID)), liquidFlag(field->getFlag(flagIDs::liquidFlagID)),
+     gasFlag(field->getFlag(flagIDs::gasFlagID)), convertedFlag(field->getFlag(flagIDs::convertedFlagID)),
+     convertToGasFlag(field->getFlag(flagIDs::convertToGasFlagID)),
+     convertToLiquidFlag(field->getFlag(flagIDs::convertToLiquidFlagID)),
+     convertFromGasToInterfaceFlag(field->getFlag(flagIDs::convertedFromGasToInterfaceFlagID)),
+     keepInterfaceForWettingFlag(field->getFlag(flagIDs::keepInterfaceForWettingFlagID)),
+     convertToInterfaceForInflowFlag(field->getFlag(flagIDs::convertToInterfaceForInflowFlagID)),
+     obstacleIDs_(obstacleIDs), outflowIDs_(outflowIDs), inflowIDs_(inflowIDs), freeSlipIDs_(freeSlipIDs)
+{
+   // create obstacleFlagMask from obstacleIDs using bitwise OR
+   obstacleFlagMask = 0;
+   for (auto obstacleID = obstacleIDs.begin(); obstacleID != obstacleIDs.end(); ++obstacleID)
+   {
+      obstacleFlagMask = flag_t(obstacleFlagMask | field->getFlag(*obstacleID));
+   }
+
+   // create outflowFlagMask from outflowIDs using bitwise OR
+   outflowFlagMask = 0;
+   for (auto outflowID = outflowIDs.begin(); outflowID != outflowIDs.end(); ++outflowID)
+   {
+      outflowFlagMask = flag_t(outflowFlagMask | field->getFlag(*outflowID));
+   }
+
+   // create inflowFlagMask from inflowIDs using bitwise OR
+   inflowFlagMask = 0;
+   for (auto inflowID = inflowIDs.begin(); inflowID != inflowIDs.end(); ++inflowID)
+   {
+      inflowFlagMask = flag_t(inflowFlagMask | field->getFlag(*inflowID));
+   }
+
+   // create freeSlipFlagMask from freeSlipIDs using bitwise OR
+   freeSlipFlagMask = 0;
+   for (auto freeSlipID = freeSlipIDs.begin(); freeSlipID != freeSlipIDs.end(); ++freeSlipID)
+   {
+      freeSlipFlagMask = flag_t(freeSlipFlagMask | field->getFlag(*freeSlipID));
+   }
+}
+
+template< typename FlagField_T >
+FlagInfo< FlagField_T >::FlagInfo(const Set< field::FlagUID >& obstacleIDs, const Set< field::FlagUID >& outflowIDs,
+                                  const Set< field::FlagUID >& inflowIDs, const Set< field::FlagUID >& freeSlipIDs)
+   : obstacleIDs_(obstacleIDs), outflowIDs_(outflowIDs), inflowIDs_(inflowIDs), freeSlipIDs_(freeSlipIDs)
+{
+   // define flags
+   flag_t nextFreeBit = flag_t(0);
+   interfaceFlag = flag_t(flag_t(1) << nextFreeBit++); // outermost flag_t is necessary to avoid warning C4334 on MSVC
+   liquidFlag    = flag_t(flag_t(1) << nextFreeBit++);
+   gasFlag       = flag_t(flag_t(1) << nextFreeBit++);
+   convertedFlag = flag_t(flag_t(1) << nextFreeBit++);
+   convertToGasFlag                = flag_t(flag_t(1) << nextFreeBit++);
+   convertToLiquidFlag             = flag_t(flag_t(1) << nextFreeBit++);
+   convertFromGasToInterfaceFlag   = flag_t(flag_t(1) << nextFreeBit++);
+   keepInterfaceForWettingFlag     = flag_t(flag_t(1) << nextFreeBit++);
+   convertToInterfaceForInflowFlag = flag_t(flag_t(1) << nextFreeBit++);
+
+   obstacleFlagMask = flag_t(0);
+   outflowFlagMask  = flag_t(0);
+   inflowFlagMask   = flag_t(0);
+   freeSlipFlagMask = flag_t(0);
+
+   // define flags for obstacles, outflow, inflow and freeSlip
+   auto setUnion = obstacleIDs + outflowIDs + inflowIDs + freeSlipIDs;
+   for (auto id = setUnion.begin(); id != setUnion.end(); ++id)
+   {
+      if (obstacleIDs.contains(*id)) { obstacleFlagMask = obstacleFlagMask | flag_t((flag_t(1) << nextFreeBit)); }
+
+      if (outflowIDs.contains(*id)) { outflowFlagMask = outflowFlagMask | flag_t((flag_t(1) << nextFreeBit)); }
+
+      if (inflowIDs.contains(*id)) { inflowFlagMask = inflowFlagMask | flag_t((flag_t(1) << nextFreeBit)); }
+
+      if (freeSlipIDs.contains(*id)) { freeSlipFlagMask = freeSlipFlagMask | flag_t((flag_t(1) << nextFreeBit)); }
+
+      nextFreeBit++;
+   }
+}
+
+template< typename FlagField_T >
+void FlagInfo< FlagField_T >::registerFlags(FlagField_T* field, const Set< field::FlagUID >& obstacleIDs,
+                                            const Set< field::FlagUID >& outflowIDs,
+                                            const Set< field::FlagUID >& inflowIDs,
+                                            const Set< field::FlagUID >& freeSlipIDs)
+{
+   flag_t nextFreeBit = flag_t(0);
+   field->registerFlag(flagIDs::interfaceFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::liquidFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::gasFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::convertedFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::convertToGasFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::convertToLiquidFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::convertedFromGasToInterfaceFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::keepInterfaceForWettingFlagID, nextFreeBit++);
+   field->registerFlag(flagIDs::convertToInterfaceForInflowFlagID, nextFreeBit++);
+
+   // extract flags
+   auto setUnion = obstacleIDs + outflowIDs + inflowIDs + freeSlipIDs;
+   for (auto id = setUnion.begin(); id != setUnion.end(); ++id)
+   {
+      field->registerFlag(*id, nextFreeBit++);
+   }
+}
+
+template< typename FlagField_T >
+void FlagInfo< FlagField_T >::registerFlags(FlagField_T* field) const
+{
+   registerFlags(field, obstacleIDs_, outflowIDs_, inflowIDs_, freeSlipIDs_);
+}
+
+template< typename FlagField_T >
+bool FlagInfo< FlagField_T >::isConsistentAcrossBlocksAndProcesses(
+   const std::weak_ptr< StructuredBlockStorage >& blockStoragePtr, ConstBlockDataID flagFieldID) const
+{
+   // check consistency across processes
+   FlagInfo rootFlagInfo = *this;
+
+   // root broadcasts its FlagInfo to all other processes
+   mpi::broadcastObject(rootFlagInfo);
+
+   // this process' FlagInfo is not identical to the one of root
+   if (rootFlagInfo != *this) { return false; }
+
+   auto blockStorage = blockStoragePtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockStorage);
+
+   // check consistency across blocks
+   for (auto blockIt = blockStorage->begin(); blockIt != blockStorage->end(); ++blockIt)
+   {
+      const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
+      FlagInfo fi(flagField, obstacleIDs_, outflowIDs_, inflowIDs_, freeSlipIDs_);
+      if (fi != *this) { return false; }
+   }
+
+   return true;
+}
+
+template< typename FlagField_T >
+bool FlagInfo< FlagField_T >::operator==(const FlagInfo& o) const
+{
+   return interfaceFlag == o.interfaceFlag && gasFlag == o.gasFlag && liquidFlag == o.liquidFlag &&
+          convertToGasFlag == o.convertToGasFlag && convertToLiquidFlag == o.convertToLiquidFlag &&
+          convertFromGasToInterfaceFlag == o.convertFromGasToInterfaceFlag &&
+          keepInterfaceForWettingFlag == o.keepInterfaceForWettingFlag &&
+          convertToInterfaceForInflowFlag == o.convertToInterfaceForInflowFlag &&
+          obstacleFlagMask == o.obstacleFlagMask && outflowFlagMask == o.outflowFlagMask &&
+          inflowFlagMask == o.inflowFlagMask && freeSlipFlagMask == o.freeSlipFlagMask;
+}
+
+template< typename FlagField_T >
+bool FlagInfo< FlagField_T >::operator!=(const FlagInfo& o) const
+{
+   return !(*this == o);
+}
+
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file

src/lbm_generated/free_surface/InitFunctions.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 InitFunctions.h
+//! \ingroup free_surface
+//! \author Matthias Markl <matthias.markl@fau.de>
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Initialization functions.
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include <functional>
+
+#include "FlagInfo.h"
+#include "InterfaceFromFillLevel.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+/***********************************************************************************************************************
+ * Initialize fill level with "value" in cells belonging to boundary and obstacles such that the bubble model does not
+ * detect obstacles as gas cells.
+ **********************************************************************************************************************/
+template< typename BoundaryHandling_T, typename Stencil_T, typename ScalarField_T >
+void initFillLevelsInBoundaries(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+                                const ConstBlockDataID& handlingID, const BlockDataID& fillFieldID,
+                                real_t value = real_c(1))
+{
+   const auto blockForest = blockForestPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      ScalarField_T* const fillField           = blockIt->getData< ScalarField_T >(fillFieldID);
+      const BoundaryHandling_T* const handling = blockIt->getData< const BoundaryHandling_T >(handlingID);
+
+      // set fill level to "value" in every cell belonging to boundary
+      WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(fillField, {
+         if (handling->isBoundary(x, y, z)) { fillField->get(x, y, z) = value; }
+      }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+   }
+}
+
+/***********************************************************************************************************************
+ * Clear and initialize flags in every cell according to the fill level.
+ **********************************************************************************************************************/
+template< typename BoundaryHandling_T, typename Stencil_T, typename FlagField_T, typename ScalarField_T >
+void initFlagsFromFillLevels(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+                             const FlagInfo< FlagField_T >& flagInfo, const BlockDataID& handlingID,
+                             const ConstBlockDataID& fillFieldID)
+{
+   const auto blockForest = blockForestPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      const ScalarField_T* const fillField = blockIt->getData< const ScalarField_T >(fillFieldID);
+      BoundaryHandling_T* const handling   = blockIt->getData< BoundaryHandling_T >(handlingID);
+
+      // clear all flags in the boundary handling
+      handling->removeFlag(flagInfo.gasFlag);
+      handling->removeFlag(flagInfo.liquidFlag);
+      handling->removeFlag(flagInfo.interfaceFlag);
+
+      WALBERLA_FOR_ALL_CELLS(fillFieldIt, fillField, {
+         // set flags only in non-boundary and non-obstacle cells
+         if (!handling->isBoundary(fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z()))
+         {
+            if (floatIsEqual(*fillFieldIt, real_c(0), real_c(1e-14)))
+            {
+               // set gas flag
+               handling->forceFlag(flagInfo.gasFlag, fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z());
+            }
+            else
+            {
+               if (*fillFieldIt < real_c(1.0))
+               {
+                  // set interface flag
+                  handling->forceFlag(flagInfo.interfaceFlag, fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z());
+               }
+               else
+               {
+                  // check if the cell is an interface cell due to direct neighboring gas cells
+                  if (isInterfaceFromFillLevel< Stencil_T >(fillFieldIt))
+                  {
+                     // set interface flag
+                     handling->forceFlag(flagInfo.interfaceFlag, fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z());
+                  }
+                  else
+                  {
+                     // set liquid flag
+                     handling->forceFlag(flagInfo.liquidFlag, fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z());
+                  }
+               }
+            }
+         }
+      }) // WALBERLA_FOR_ALL_CELLS
+   }
+}
+
+/***********************************************************************************************************************
+ * Initialize the hydrostatic pressure in the direction in which a force is acting in ALL cells (regardless of a cell's
+ * flag). The velocity remains unchanged.
+ *
+ * The force vector must have only one component, i.e., the direction of the force can only be in x-, y- or z-axis.
+ * The variable fluidHeight determines the height at which the density is equal to reference density (=1).
+ **********************************************************************************************************************/
+template< typename ScalarField_T >
+void initHydrostaticPressure(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+                             const BlockDataID& densityFieldID, const Vector3< real_t >& force, real_t fluidHeight)
+{
+   // count number of non-zero components of the force vector
+   uint_t numForceComponents = uint_c(0);
+   if (!realIsEqual(force[0], real_c(0), real_c(1e-14))) { ++numForceComponents; }
+   if (!realIsEqual(force[1], real_c(0), real_c(1e-14))) { ++numForceComponents; }
+   if (!realIsEqual(force[2], real_c(0), real_c(1e-14))) { ++numForceComponents; }
+
+   WALBERLA_CHECK_EQUAL(numForceComponents, uint_c(1),
+                        "The current implementation of the hydrostatic pressure initialization does not support "
+                        "forces that have none or multiple components, i. e., a force that points in a direction other "
+                        "than the x-, y- or z-axis.");
+
+   const auto blockForest = blockForestPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      ScalarField_T* const densityField = blockIt->getData< ScalarField_T >(densityFieldID);
+      CellInterval local = densityField->xyzSizeWithGhostLayer(); // block-, i.e., process-local cell interval
+
+      for (auto cellIt = local.begin(); cellIt != local.end(); ++cellIt)
+      {
+         // transform the block-local coordinate to global coordinates
+         Cell global;
+         blockForest->transformBlockLocalToGlobalCell(global, *blockIt, *cellIt);
+
+         // get the current global coordinate, i.e., height of the fluid in the relevant direction
+         cell_idx_t coordinate = cell_idx_c(0);
+         real_t forceComponent = real_c(0);
+         if (!realIsEqual(force[0], real_c(0), real_c(1e-14)))
+         {
+            coordinate     = global.x();
+            forceComponent = force[0];
+         }
+         else
+         {
+            if (!realIsEqual(force[1], real_c(0), real_c(1e-14)))
+            {
+               coordinate     = global.y();
+               forceComponent = force[1];
+            }
+            else
+            {
+               if (!realIsEqual(force[2], real_c(0), real_c(1e-14)))
+               {
+                  coordinate     = global.z();
+                  forceComponent = force[2];
+               }
+               else
+               {
+                  WALBERLA_ABORT(
+                     "The current implementation of the hydrostatic pressure initialization does not support "
+                     "forces that have none or multiple components, i. e., a force that points in a direction other "
+                     "than the x-, y- or z-axis.")
+               }
+            }
+         }
+
+         // initialize the (hydrostatic) pressure, i.e., LBM density
+         // Bernoulli: p = p0 + density * gravity * height
+         // => LBM (density=1): rho = rho0 + gravity * height = 1 + 1/cs^2 * g * h = 1 + 3 * g * h
+         // shift global cell by 0.5 since density is set for cell center
+         densityField->get(*cellIt) = real_c(1) + real_c(3) * forceComponent * (real_c(coordinate) + real_c(0.5) - std::ceil(fluidHeight));
+      }
+   }
+}
+
+/***********************************************************************************************************************
+ * Initialize the force density field with a given acceleration and density of each cell.
+ * Differs from the version above by using a flattened vector field (GhostLayerField< real_t, 3 >). This is necessary
+ * because Pystencils does not support VectorField_T (GhostLayerField< Vector3<real_t>, 1 >).
+ **********************************************************************************************************************/
+template< typename FlagField_T, typename VectorField_T, typename ScalarField_T >
+void initForceDensityField(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+                                  const BlockDataID& forceDensityFieldID, const ConstBlockDataID& fillFieldID,
+                                  const ConstBlockDataID& densityFieldID, const ConstBlockDataID& flagFieldID,
+                                  const FlagInfo< FlagField_T >& flagInfo, const Vector3< real_t >& acceleration)
+{
+   const auto blockForest = blockForestPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      VectorField_T* const forceDensityField = blockIt->getData< VectorField_T >(forceDensityFieldID);
+      const ScalarField_T* const fillField            = blockIt->getData< const ScalarField_T >(fillFieldID);
+      const ScalarField_T* const densityField                = blockIt->getData< const ScalarField_T >(densityFieldID);
+      const FlagField_T* const flagField              = blockIt->getData< const FlagField_T >(flagFieldID);
+
+      WALBERLA_FOR_ALL_CELLS(forceDensityFieldIt, forceDensityField, fillFieldIt, fillField, densityFieldIt, densityField,
+                             flagFieldIt, flagField, {
+                                // set force density in cells to acceleration * density * fillLevel (see equation 15
+                                // in Koerner et al., 2005);
+
+                                if (flagInfo.isInterface(*flagFieldIt))
+                                {
+                                   const real_t density   = densityField->get(densityFieldIt.cell());
+                                   forceDensityFieldIt[0] = acceleration[0] * *fillFieldIt * density;
+                                   forceDensityFieldIt[1] = acceleration[1] * *fillFieldIt * density;
+                                   forceDensityFieldIt[2] = acceleration[2] * *fillFieldIt * density;
+                                }
+
+                                else
+                                {
+                                   if (flagInfo.isLiquid(*flagFieldIt))
+                                   {
+                                      const real_t density   = densityField->get(densityFieldIt.cell());
+                                      forceDensityFieldIt[0] = acceleration[0] * density;
+                                      forceDensityFieldIt[1] = acceleration[1] * density;
+                                      forceDensityFieldIt[2] = acceleration[2] * density;
+                                   }
+                                }
+                             }) // WALBERLA_FOR_ALL_CELLS
+   }
+}
+
+/***********************************************************************************************************************
+ * Set density in non-liquid and non-interface cells to 1.
+ **********************************************************************************************************************/
+template< typename FlagField_T, typename VectorField_T, typename ScalarField_T>
+void setDensityInNonFluidCellsToOne(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+                                    const FlagInfo< FlagField_T >& flagInfo, const ConstBlockDataID& flagFieldID,
+                                    const BlockDataID& velocityFieldID, const BlockDataID& densityFieldID)
+{
+   const auto blockForest = blockForestPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+   for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+   {
+      VectorField_T* const velocityField         = blockIt->getData< VectorField_T >(velocityFieldID);
+      ScalarField_T* const densityField         = blockIt->getData< ScalarField_T >(densityFieldID);
+
+      const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
+
+      WALBERLA_FOR_ALL_CELLS(velocityFieldIt, velocityField, densityFieldIt, densityField, flagFieldIt, flagField, {
+         if (!flagInfo.isLiquid(*flagFieldIt) && !flagInfo.isInterface(*flagFieldIt))
+         {
+            // set density in gas cells to 1 and velocity to zero
+            const cell_idx_t x = densityFieldIt.x();
+            const cell_idx_t y = densityFieldIt.y();
+            const cell_idx_t z = densityFieldIt.z();
+            densityField->get(x, y, z) = real_c(1.0);
+            velocityField->get(x, y, z, 0) = real_c(0.0);
+            velocityField->get(x, y, z, 1) = real_c(0.0);
+            velocityField->get(x, y, z, 2) = real_c(0.0);
+         }
+      }) // WALBERLA_FOR_ALL_CELLS
+   }
+}
+} // namespace free_surface
+} // namespace walberla

src/lbm_generated/free_surface/InterfaceFromFillLevel.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 InterfaceFromFillLevel.h
+//! \ingroup free_surface
+//! \author Matthias Markl <matthias.markl@fau.de>
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Check whether a cell should be an interface cell according to its properties.
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/cell/Cell.h"
+
+#include <type_traits>
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+/***********************************************************************************************************************
+ * Check whether a cell is an interface cell with respect to its fill level and its direct neighborhood (liquid cells
+ * can not have gas cell in their direct neighborhood; therefore, the liquid cell is forced to be an interface cell).
+ **********************************************************************************************************************/
+template< typename Stencil_T, typename ScalarField_T >
+inline bool isInterfaceFromFillLevel(const ScalarField_T& fillField, cell_idx_t x, cell_idx_t y, cell_idx_t z)
+{
+   WALBERLA_ASSERT(std::is_floating_point< typename ScalarField_T::value_type >::value,
+                   "Fill level field must contain floating point values.")
+
+   real_t fillLevel = fillField.get(x, y, z);
+
+   // this cell is regular gas cell
+   if (floatIsEqual(fillLevel, real_c(0.0), real_c(1e-14))) { return false; }
+
+   // this cell is regular interface cell
+   if (fillLevel < real_c(1.0)) { return true; }
+
+   // check this cell's direct neighborhood for gas cells (a liquid cell can not be direct neighbor to a gas cell)
+   for (auto d = Stencil_T::beginNoCenter(); d != Stencil_T::end(); ++d)
+   {
+      // this cell has a gas cell in its direct neighborhood; it therefore must not be a liquid cell but an interface
+      // cell although its fill level is 1.0
+      if (fillField.get(x + d.cx(), y + d.cy(), z + d.cz()) <= real_c(0.0)) { return true; }
+   }
+
+   // this cell is a regular fluid cell
+   return false;
+}
+
+template< typename Stencil_T, typename ScalarFieldIt_T >
+inline bool isInterfaceFromFillLevel(const ScalarFieldIt_T& fillFieldIt)
+{
+   return isInterfaceFromFillLevel< Stencil_T, typename ScalarFieldIt_T::FieldType >(
+      *(fillFieldIt.getField()), fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z());
+}
+
+template< typename Stencil_T, typename ScalarField >
+inline bool isInterfaceFromFillLevel(const ScalarField& fillField, const Cell& cell)
+{
+   return isInterfaceFromFillLevel< Stencil_T >(fillField, cell.x(), cell.y(), cell.z());
+}
+
+} // namespace free_surface
+} // namespace walberla

src/lbm_generated/free_surface/LoadBalancing.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 LoadBalancing.h
+//! \ingroup free_surface
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Free surface-specific functionality for load balancing.
+//
+//======================================================================================================================
+
+#include "blockforest/communication/UniformBufferedScheme.h"
+#include "blockforest/BlockForest.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/StructuredBlockForest.h"
+#include "blockforest/loadbalancing/DynamicCurve.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
+#include "core/logging/Logging.h"
+#include "core/math/AABB.h"
+#include "core/math/DistributedSample.h"
+#include "core/math/Vector3.h"
+#include "core/mpi/MPIManager.h"
+
+#include "lbm_generated/blockforest/SimpleCommunication.h"
+#include "lbm_generated/free_surface/BlockStateDetectorSweep.h"
+
+#include <algorithm>
+#include <numeric>
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+
+template< typename FlagField_T >
+class ProcessLoadEvaluator;
+
+/***********************************************************************************************************************
+ * Create non-uniform block forest to be used for load balancing.
+ **********************************************************************************************************************/
+std::shared_ptr< StructuredBlockForest > createNonUniformBlockForest(const Vector3< uint_t >& domainSize,
+                                                                     const Vector3< uint_t >& cellsPerBlock,
+                                                                     const Vector3< uint_t >& numBlocks,
+                                                                     const Vector3< bool >& periodicity)
+{
+   WALBERLA_CHECK_EQUAL(domainSize[0], cellsPerBlock[0] * numBlocks[0],
+                        "The domain size is not divisible by the specified \"cellsPerBlock\" in x-direction.");
+   WALBERLA_CHECK_EQUAL(domainSize[1], cellsPerBlock[1] * numBlocks[1],
+                        "The domain size is not divisible by the specified \"cellsPerBlock\" in y-direction.");
+   WALBERLA_CHECK_EQUAL(domainSize[2], cellsPerBlock[2] * numBlocks[2],
+                        "The domain size is not divisible by the specified \"cellsPerBlock\" in z-direction.");
+
+   // create SetupBlockForest for allowing load balancing
+   SetupBlockForest setupBlockForest;
+
+   AABB domainAABB(real_c(0), real_c(0), real_c(0), real_c(domainSize[0]), real_c(domainSize[1]),
+                   real_c(domainSize[2]));
+
+   setupBlockForest.init(domainAABB, numBlocks[0], numBlocks[1], numBlocks[2], periodicity[0], periodicity[1],
+                         periodicity[2]);
+
+   // compute initial process distribution
+   setupBlockForest.balanceLoad(blockforest::StaticLevelwiseCurveBalance(true),
+                                uint_c(MPIManager::instance()->numProcesses()));
+
+   WALBERLA_LOG_INFO_ON_ROOT(setupBlockForest);
+
+   // define MPI communicator
+   if (!MPIManager::instance()->rankValid()) { MPIManager::instance()->useWorldComm(); }
+
+   // create BlockForest (will be encapsulated in StructuredBlockForest)
+   const std::shared_ptr< BlockForest > blockForest =
+      std::make_shared< BlockForest >(uint_c(MPIManager::instance()->rank()), setupBlockForest, false);
+
+   // create StructuredBlockForest
+   std::shared_ptr< StructuredBlockForest > structuredBlockForest =
+      std::make_shared< StructuredBlockForest >(blockForest, cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2]);
+   structuredBlockForest->createCellBoundingBoxes();
+
+   return structuredBlockForest;
+}
+
+/***********************************************************************************************************************
+ * Example for a load balancing implementation.
+ *
+ * IMPORTANT REMARK: The following implementation should be considered a demonstrator based on best-practices. That is,
+ * it was not thoroughly benchmarked and does not guarantee a performance-optimal load balancing.
+ **********************************************************************************************************************/
+template< typename FlagField_T, typename CommunicationStencil_T, typename LBMStencil_T >
+class LoadBalancer
+{
+ public:
+   LoadBalancer(const std::shared_ptr< StructuredBlockForest >& blockForestPtr,
+                const lbm_generated::SimpleCommunication< CommunicationStencil_T >& communication,
+                const blockforest::communication::UniformBufferedScheme< LBMStencil_T >& pdfCommunication,
+                uint_t blockWeightFullFreeSurface,
+                uint_t blockWeightOnlyLBM, uint_t blockWeightOnlyGasAndBoundary, uint_t frequency,
+                bool printStatistics = false)
+      : blockForest_(blockForestPtr), communication_(communication), pdfCommunication_(pdfCommunication),
+        blockWeightFullFreeSurface_(blockWeightFullFreeSurface),
+        blockWeightOnlyLBM_(blockWeightOnlyLBM), blockWeightOnlyGasAndBoundary_(blockWeightOnlyGasAndBoundary),
+        frequency_(frequency), printStatistics_(printStatistics), executionCounter_(uint_c(0)),
+        evaluator_(ProcessLoadEvaluator< FlagField_T >(blockForest_, blockWeightFullFreeSurface_, blockWeightOnlyLBM_,
+                                                       blockWeightOnlyGasAndBoundary_, uint_c(1)))
+   {
+      BlockForest& blockForest = blockForest_->getBlockForest();
+
+      // refinement is not implemented in FSLBM such that this can be set to false
+      blockForest.recalculateBlockLevelsInRefresh(false);
+
+      // rebalancing of blocks must be forced here, as it would normally be done when refinement levels change
+      blockForest.alwaysRebalanceInRefresh(true);
+
+      // depth of levels is not changing and must therefore not be communicated
+      blockForest.allowRefreshChangingDepth(false);
+
+      // refinement is not implemented in FSLBM such that this can be set to false
+      blockForest.allowMultipleRefreshCycles(false);
+
+      // leave algorithm when load balancing has been performed
+      blockForest.checkForEarlyOutAfterLoadBalancing(true);
+
+      // use PhantomWeight as defined below
+      blockForest.setRefreshPhantomBlockDataAssignmentFunction(blockWeightAssignment);
+      blockForest.setRefreshPhantomBlockDataPackFunction(phantomWeightsPack);
+      blockForest.setRefreshPhantomBlockDataUnpackFunction(phantomWeightsUnpack);
+
+      // assign load balancing function
+      blockForest.setRefreshPhantomBlockMigrationPreparationFunction(
+         blockforest::DynamicCurveBalance< PhantomWeight >(true, true));
+   }
+
+   void operator()()
+   {
+      if (frequency_ == uint_c(0)) { return; }
+
+      // only balance load in given frequencies
+      if (executionCounter_ % frequency_ == uint_c(0))
+      {
+         BlockForest& blockForest = blockForest_->getBlockForest();
+
+         // balance load by updating the blockForest
+         const uint_t modificationStamp = blockForest.getModificationStamp();
+         blockForest.refresh();
+
+         const uint_t newModificationStamp = blockForest.getModificationStamp();
+
+         if (newModificationStamp != modificationStamp)
+         {
+            // communicate all fields
+            communication_();
+            pdfCommunication_();
+
+            if (printStatistics_) { evaluator_(); }
+
+            if (blockForest.getNumberOfBlocks() == uint_c(0))
+            {
+               WALBERLA_ABORT(
+                  "Load balancing lead to a situation where there is a process with no blocks. This is "
+                  "not supported yet. This can be avoided by either using smaller blocks or, equivalently, more blocks "
+                  "per process.");
+            }
+         }
+      }
+      ++executionCounter_;
+   }
+
+ private:
+   std::shared_ptr< StructuredBlockForest > blockForest_;
+   lbm_generated::SimpleCommunication< CommunicationStencil_T > communication_;
+   blockforest::communication::UniformBufferedScheme< LBMStencil_T > pdfCommunication_;
+
+   uint_t blockWeightFullFreeSurface_;    // empirical choice, not thoroughly benchmarked: 50
+   uint_t blockWeightOnlyLBM_;            // empirical choice, not thoroughly benchmarked: 10
+   uint_t blockWeightOnlyGasAndBoundary_; // empirical choice, not thoroughly benchmarked: 5
+
+   uint_t frequency_;
+   bool printStatistics_;
+
+   uint_t executionCounter_;
+
+   ProcessLoadEvaluator< FlagField_T > evaluator_;
+
+   class PhantomWeight // used as a 'PhantomBlockForest::PhantomBlockDataAssignmentFunction'
+   {
+    public:
+      using weight_t = uint_t;
+      PhantomWeight(const weight_t _weight) : weight_(_weight) {}
+      weight_t weight() const { return weight_; }
+
+    private:
+      weight_t weight_;
+   }; // class PhantomWeight
+
+   std::function< void(mpi::SendBuffer& buffer, const PhantomBlock& block) > phantomWeightsPack =
+      [](mpi::SendBuffer& buffer, const PhantomBlock& block) { buffer << block.getData< PhantomWeight >().weight(); };
+
+   std::function< void(mpi::RecvBuffer& buffer, const PhantomBlock&, walberla::any& data) > phantomWeightsUnpack =
+      [](mpi::RecvBuffer& buffer, const PhantomBlock&, walberla::any& data) {
+         typename PhantomWeight::weight_t w;
+         buffer >> w;
+         data = PhantomWeight(w);
+      };
+
+   std::function< void(std::vector< std::pair< const PhantomBlock*, walberla::any > >& blockData,
+                       const PhantomBlockForest&) >
+      blockWeightAssignment =
+         [this](std::vector< std::pair< const PhantomBlock*, walberla::any > >& blockData, const PhantomBlockForest&) {
+            for (auto it = blockData.begin(); it != blockData.end(); ++it)
+            {
+               if (it->first->getState().contains(BlockStateDetectorSweep< FlagField_T >::fullFreeSurface))
+               {
+                  it->second = PhantomWeight(blockWeightFullFreeSurface_);
+               }
+               else
+               {
+                  if (it->first->getState().contains(BlockStateDetectorSweep< FlagField_T >::onlyLBM))
+                  {
+                     it->second = PhantomWeight(blockWeightOnlyLBM_);
+                  }
+                  else
+                  {
+                     if (it->first->getState().contains(BlockStateDetectorSweep< FlagField_T >::onlyGasAndBoundary))
+                     {
+                        it->second = PhantomWeight(blockWeightOnlyGasAndBoundary_);
+                     }
+                     else { WALBERLA_ABORT("Unknown block state"); }
+                  }
+               }
+            }
+         };
+
+}; // class LoadBalancer
+
+/***********************************************************************************************************************
+ * Evaluates and prints statistics about the current load distribution situation:
+ * - Average weight per process
+ * - Maximum weight per process
+ * - Minimum weight per process
+ **********************************************************************************************************************/
+template< typename FlagField_T >
+class ProcessLoadEvaluator
+{
+ public:
+   ProcessLoadEvaluator(const std::weak_ptr< const StructuredBlockForest >& blockForest,
+                        uint_t blockWeightFullFreeSurface, uint_t blockWeightOnlyLBM,
+                        uint_t blockWeightOnlyGasAndBoundary, uint_t frequency)
+      : blockForest_(blockForest), blockWeightFullFreeSurface_(blockWeightFullFreeSurface),
+        blockWeightOnlyLBM_(blockWeightOnlyLBM), blockWeightOnlyGasAndBoundary_(blockWeightOnlyGasAndBoundary),
+        frequency_(frequency), executionCounter_(uint_c(0))
+   {}
+
+   void operator()()
+   {
+      if (frequency_ == uint_c(0)) { return; }
+
+      ++executionCounter_;
+
+      // only evaluate in given intervals
+      if (executionCounter_ % frequency_ != uint_c(0) && executionCounter_ != uint_c(1)) { return; }
+
+      std::vector< real_t > weightSum = computeWeightSumPerProcess();
+
+      print(weightSum);
+   }
+
+   std::vector< real_t > computeWeightSumPerProcess()
+   {
+      auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      std::vector< real_t > weightSum(uint_c(MPIManager::instance()->numProcesses()), real_c(0));
+
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         if (blockForest->blockExistsLocally(blockIt->getId()))
+         {
+            if (blockIt->getState().contains(BlockStateDetectorSweep< FlagField_T >::fullFreeSurface))
+            {
+               weightSum[blockForest->getProcessRank(blockIt->getId())] += real_c(blockWeightFullFreeSurface_);
+            }
+            else
+            {
+               if (blockIt->getState().contains(BlockStateDetectorSweep< FlagField_T >::onlyLBM))
+               {
+                  weightSum[blockForest->getProcessRank(blockIt->getId())] += real_c(blockWeightOnlyLBM_);
+               }
+               else
+               {
+                  if (blockIt->getState().contains(BlockStateDetectorSweep< FlagField_T >::onlyGasAndBoundary))
+                  {
+                     weightSum[blockForest->getProcessRank(blockIt->getId())] += real_c(blockWeightOnlyGasAndBoundary_);
+                  }
+               }
+            }
+         }
+      }
+
+      mpi::reduceInplace< real_t >(weightSum, mpi::SUM, 0);
+
+      return weightSum;
+   }
+
+   void print(const std::vector< real_t >& weightSum)
+   {
+      WALBERLA_ROOT_SECTION()
+      {
+         const std::vector< real_t >::const_iterator max = std::max_element(weightSum.cbegin(), weightSum.end());
+         const std::vector< real_t >::const_iterator min = std::min_element(weightSum.cbegin(), weightSum.end());
+         const real_t sum = std::accumulate(weightSum.cbegin(), weightSum.end(), real_c(0));
+         const real_t avg = sum / real_c(MPIManager::instance()->numProcesses());
+
+         WALBERLA_LOG_INFO("Load balancing:");
+         WALBERLA_LOG_INFO("\t Average weight per process " << avg);
+         WALBERLA_LOG_INFO("\t Maximum weight per process " << *max);
+         WALBERLA_LOG_INFO("\t Minimum weight per process " << *min);
+      }
+   }
+
+ private:
+   std::weak_ptr< const StructuredBlockForest > blockForest_;
+   uint_t blockWeightFullFreeSurface_;
+   uint_t blockWeightOnlyLBM_;
+   uint_t blockWeightOnlyGasAndBoundary_;
+
+   uint_t frequency_;
+   uint_t executionCounter_;
+}; // class ProcessLoadEvaluator
+
+} // namespace free_surface
+} // namespace walberla

src/lbm_generated/free_surface/MaxVelocityComputer.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 MaxVelocityComputer.h
+//! \ingroup free_surface
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Compute the maximum velocity of all liquid and interface cells (in each direction) in the system.
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+template< typename FreeSurfaceBoundaryHandling_T, typename PdfField_T, typename FlagField_T >
+class MaxVelocityComputer
+{
+ public:
+   MaxVelocityComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest,
+                       const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
+                       const ConstBlockDataID& pdfFieldID, uint_t frequency,
+                       const std::shared_ptr< Vector3< real_t > >& maxVelocity)
+      : blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling), pdfFieldID_(pdfFieldID),
+        maxVelocity_(maxVelocity), frequency_(frequency), executionCounter_(uint_c(0))
+   {}
+
+   void operator()()
+   {
+      if (frequency_ == uint_c(0)) { return; }
+
+      auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      auto freeSurfaceBoundaryHandling = freeSurfaceBoundaryHandling_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(freeSurfaceBoundaryHandling);
+
+      if (executionCounter_ == uint_c(0)) { getMaxVelocity(blockForest, freeSurfaceBoundaryHandling); }
+      else
+      {
+         // only evaluate in given frequencies
+         if (executionCounter_ % frequency_ == uint_c(0)) { getMaxVelocity(blockForest, freeSurfaceBoundaryHandling); }
+      }
+
+      ++executionCounter_;
+   }
+
+   void getMaxVelocity(const std::shared_ptr< const StructuredBlockForest >& blockForest,
+                       const std::shared_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling)
+   {
+      const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
+      const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
+
+      real_t maxVelocityX = real_c(0);
+      real_t maxVelocityY = real_c(0);
+      real_t maxVelocityZ = real_c(0);
+
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         const FlagField_T* const flagField = blockIt->template getData< const FlagField_T >(flagFieldID);
+         const PdfField_T* const pdfField   = blockIt->template getData< const PdfField_T >(pdfFieldID_);
+
+            WALBERLA_FOR_ALL_CELLS_OMP(flagFieldIt, flagField, pdfFieldIt, pdfField,
+                                       omp parallel for schedule(static) reduction(max:maxVelocityX)
+                                                                         reduction(max:maxVelocityY)
+                                                                         reduction(max:maxVelocityZ),
+            {
+            if (flagInfo.isLiquid(flagFieldIt) || flagInfo.isInterface(flagFieldIt))
+            {
+               const Vector3< real_t > velocity = pdfField->getVelocity(pdfFieldIt.cell());
+
+               if (velocity[0] > maxVelocityX) { maxVelocityX = velocity[0]; }
+               if (velocity[1] > maxVelocityY) { maxVelocityY = velocity[1]; }
+               if (velocity[2] > maxVelocityZ) { maxVelocityZ = velocity[2]; }
+            }
+            }) // WALBERLA_FOR_ALL_CELLS_OMP
+      }
+
+      Vector3< real_t > maxVelocity(maxVelocityX, maxVelocityY, maxVelocityZ);
+      mpi::allReduceInplace< real_t >(maxVelocity, mpi::MAX);
+
+      *maxVelocity_ = maxVelocity;
+   };
+
+ private:
+   std::weak_ptr< const StructuredBlockForest > blockForest_;
+   std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
+
+   const ConstBlockDataID pdfFieldID_;
+
+   std::shared_ptr< Vector3< real_t > > maxVelocity_;
+
+   uint_t frequency_;
+   uint_t executionCounter_;
+}; // class MaxVelocityComputer
+
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file

src/lbm_generated/free_surface/SurfaceMeshWriter.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 SurfaceMeshWriter.h
+//! \ingroup free_surface
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Free surface-specific class for writing the free surface as triangle mesh.
+//
+//======================================================================================================================
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "core/Filesystem.h"
+
+#include "field/AddToStorage.h"
+
+#include "geometry/mesh/TriangleMeshIO.h"
+
+#include "postprocessing/FieldToSurfaceMesh.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+namespace abortIfNullptr
+{
+
+// helper function to check validity of std::weak_ptr in constructors' initializer list; WALBERLA_CHECK_NOT_NULLPTR()
+// does not work there, because the macro terminates with ";"
+template< typename T >
+void abortIfNullptr(const std::weak_ptr< T >& weakPointer)
+{
+   if (weakPointer.lock() == nullptr) { WALBERLA_ABORT("Weak pointer has expired."); }
+}
+} // namespace abortIfNullptr
+
+/***********************************************************************************************************************
+ * Write free surface as triangle mesh.
+ *
+ * Internally, a clone of the fill level field is stored and all cells not marked as liquidInterfaceGasFlagIDSet are set
+ * to "obstacleFillLevel" in the cloned field. This is done to avoid writing e.g. obstacle cells that were possibly
+ * assigned a fill level of 1 to not make them detect as gas cells in the bubble model.
+ **********************************************************************************************************************/
+template< typename ScalarField_T, typename FlagField_T >
+class SurfaceMeshWriter
+{
+ public:
+   SurfaceMeshWriter(const std::weak_ptr< StructuredBlockForest >& blockForest, const ConstBlockDataID& fillFieldID,
+                     const ConstBlockDataID& flagFieldID, const Set< FlagUID >& liquidInterfaceGasFlagIDSet,
+                     real_t obstacleFillLevel, uint_t writeFrequency, const std::string& baseFolder)
+      : blockForest_(blockForest), fillFieldID_(fillFieldID), flagFieldID_(flagFieldID),
+        liquidInterfaceGasFlagIDSet_(liquidInterfaceGasFlagIDSet), obstacleFillLevel_(obstacleFillLevel),
+        writeFrequency_(writeFrequency), baseFolder_(baseFolder), executionCounter_(uint_c(0)),
+        fillFieldCloneID_(
+           (abortIfNullptr::abortIfNullptr(blockForest),
+            field::addCloneToStorage< ScalarField_T >(blockForest_.lock(), fillFieldID_, "Fill level field clone")))
+   {}
+
+   // config block must be named "MeshOutputParameters"
+   SurfaceMeshWriter(const std::weak_ptr< StructuredBlockForest >& blockForest, const ConstBlockDataID& fillFieldID,
+                     const ConstBlockDataID& flagFieldID, const Set< FlagUID >& liquidInterfaceGasFlagIDSet,
+                     real_t obstacleFillLevel, const std::weak_ptr< Config >& config)
+      : SurfaceMeshWriter(
+           blockForest, fillFieldID, flagFieldID, liquidInterfaceGasFlagIDSet, obstacleFillLevel,
+           (abortIfNullptr::abortIfNullptr(config),
+            config.lock()->getOneBlock("MeshOutputParameters").getParameter< uint_t >("writeFrequency")),
+           (abortIfNullptr::abortIfNullptr(config),
+            config.lock()->getOneBlock("MeshOutputParameters").getParameter< std::string >("baseFolder")))
+   {}
+
+   void operator()()
+   {
+      if (writeFrequency_ == uint_c(0)) { return; }
+
+      if (executionCounter_ == uint_c(0))
+      {
+         createBaseFolder();
+         writeMesh();
+      }
+      else { writeMesh(); }
+
+      ++executionCounter_;
+   }
+
+ private:
+   void createBaseFolder() const
+   {
+      WALBERLA_ROOT_SECTION()
+      {
+         const filesystem::path basePath(baseFolder_);
+         if (filesystem::exists(basePath)) { filesystem::remove_all(basePath); }
+         filesystem::create_directories(basePath);
+      }
+      WALBERLA_MPI_BARRIER();
+   }
+
+   void writeMesh()
+   {
+      // only write mesh in given frequency
+      if (executionCounter_ % writeFrequency_ != uint_c(0)) { return; }
+
+      // rank=0 is just an arbitrary choice here
+      const int targetRank = 0;
+
+      auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      // update clone of fill level field and set fill level of all non-liquid, -interface, or -gas cells to zero
+      updateFillFieldClone(blockForest);
+
+      const auto surfaceMesh = postprocessing::realFieldToSurfaceMesh< ScalarField_T >(
+         blockForest, fillFieldCloneID_, real_c(0.5), uint_c(0), true, targetRank, MPI_COMM_WORLD);
+
+      WALBERLA_EXCLUSIVE_WORLD_SECTION(targetRank)
+      {
+         geometry::writeMesh(baseFolder_ + "/" + "simulation_step_" + std::to_string(executionCounter_) + ".obj",
+                             *surfaceMesh);
+      }
+   }
+
+   // update clone of fill level field and set fill level of all non-liquid, -interface, or -gas cells to zero;
+   // explicitly use shared_ptr instead of weak_ptr to avoid checking the latter's validity (is done in writeMesh()
+   // already)
+   void updateFillFieldClone(const shared_ptr< StructuredBlockForest >& blockForest)
+   {
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         ScalarField_T* const fillFieldClone  = blockIt->template getData< ScalarField_T >(fillFieldCloneID_);
+         const ScalarField_T* const fillField = blockIt->template getData< const ScalarField_T >(fillFieldID_);
+         const FlagField_T* const flagField   = blockIt->template getData< const FlagField_T >(flagFieldID_);
+
+         const auto liquidInterfaceGasFlagMask = flagField->getMask(liquidInterfaceGasFlagIDSet_);
+
+         WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(fillFieldClone, fillField->nrOfGhostLayers(), {
+            const typename ScalarField_T::Ptr fillFieldClonePtr(*fillFieldClone, x, y, z);
+            const typename ScalarField_T::ConstPtr fillFieldPtr(*fillField, x, y, z);
+            const typename FlagField_T::ConstPtr flagFieldPtr(*flagField, x, y, z);
+
+            // set fill level to zero in every non-liquid, -interface, or -gas cell
+            if (!isPartOfMaskSet(flagFieldPtr, liquidInterfaceGasFlagMask)) { *fillFieldClonePtr = obstacleFillLevel_; }
+            else
+            {
+               // copy fill level from fill level field
+               *fillFieldClonePtr = *fillFieldPtr;
+            }
+         }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+      }
+   }
+
+   std::weak_ptr< StructuredBlockForest > blockForest_;
+   ConstBlockDataID fillFieldID_;
+   ConstBlockDataID flagFieldID_;
+   Set< FlagUID > liquidInterfaceGasFlagIDSet_;
+   real_t obstacleFillLevel_;
+   uint_t writeFrequency_;
+   std::string baseFolder_;
+   uint_t executionCounter_;
+
+   BlockDataID fillFieldCloneID_;
+}; // class SurfaceMeshWriter
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file

src/lbm_generated/free_surface/TotalMassComputer.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 TotalMassComputer.h
+//! \ingroup free_surface
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Compute the total mass of the system (including mass from the excessMassField).
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+#include "core/StringUtility.h"
+
+#include "lbm_generated/free_surface/dynamics/SurfaceDynamicsHandler.h"
+#include "lbm_generated/free_surface/surface_geometry/SurfaceGeometryHandler.h"
+#include "lbm_generated/free_surface/surface_geometry/Utility.h"
+#include "lbm_generated/free_surface/FlagInfo.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+template< typename PdfField_T, typename FlagField_T, typename ScalarField_T >
+class TotalMassComputer
+{
+ public:
+   TotalMassComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest,
+                     const ConstBlockDataID& flagFieldID, const FlagInfo<FlagField_T>& flagInfo,
+                     const ConstBlockDataID& pdfFieldID, const ConstBlockDataID& fillFieldID, uint_t frequency,
+                     const std::shared_ptr< real_t >& totalMass)
+      : blockForest_(blockForest), flagFieldID_(flagFieldID), flagInfo_(flagInfo), pdfFieldID_(pdfFieldID),
+        fillFieldID_(fillFieldID), totalMass_(totalMass), frequency_(frequency), executionCounter_(uint_c(0))
+   {}
+
+   TotalMassComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest,
+                     const ConstBlockDataID& flagFieldID, const FlagInfo<FlagField_T>& flagInfo,
+                     const ConstBlockDataID& pdfFieldID, const ConstBlockDataID& fillFieldID,
+                     const ConstBlockDataID& excessMassFieldID, uint_t frequency,
+                     const std::shared_ptr< real_t >& totalMass, const std::shared_ptr< real_t >& excessMass)
+      : blockForest_(blockForest), flagFieldID_(flagFieldID), flagInfo_(flagInfo), pdfFieldID_(pdfFieldID),
+        fillFieldID_(fillFieldID), excessMassFieldID_(excessMassFieldID), totalMass_(totalMass),
+        excessMass_(excessMass), frequency_(frequency), executionCounter_(uint_c(0))
+   {}
+
+   void operator()()
+   {
+      if (frequency_ == uint_c(0)) { return; }
+
+      auto blockForest = blockForest_.lock();
+      WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+      // only evaluate in given frequencies
+      if (executionCounter_ % frequency_ == uint_c(0) || executionCounter_ == uint_c(0))
+      {
+         computeMass(blockForest, flagFieldID_, flagInfo_);
+      }
+
+      ++executionCounter_;
+   }
+
+   void computeMass(const std::shared_ptr< const StructuredBlockForest >& blockForest,
+                    const ConstBlockDataID& flagFieldID, const FlagInfo<FlagField_T>& flagInfo)
+   {
+      real_t mass       = real_c(0);
+      real_t excessMass = real_c(0);
+
+      for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+      {
+         const FlagField_T* const flagField   = blockIt->template getData< const FlagField_T >(flagFieldID);
+         const PdfField_T* const pdfField     = blockIt->template getData< const PdfField_T >(pdfFieldID_);
+         const ScalarField_T* const fillField = blockIt->template getData< const ScalarField_T >(fillFieldID_);
+
+         // if provided, also consider mass stored in excessMassField
+         if (excessMassFieldID_ != ConstBlockDataID())
+         {
+            const ScalarField_T* const excessMassField =
+               blockIt->template getData< const ScalarField_T >(excessMassFieldID_);
+
+            WALBERLA_FOR_ALL_CELLS_OMP(flagFieldIt, flagField, pdfFieldIt, pdfField, fillFieldIt, fillField,
+                                       excessMassFieldIt, excessMassField,
+                                       omp parallel for schedule(static) reduction(+:mass),
+            {
+               if (flagInfo.isLiquid(flagFieldIt) || flagInfo.isInterface(flagFieldIt))
+               {
+                  //TODO Implement
+                  const real_t density = 1.0; // pdfField->getDensity(pdfFieldIt.cell());
+                  mass += *fillFieldIt * density + *excessMassFieldIt;
+
+                  if (excessMass_ != nullptr) { excessMass += *excessMassFieldIt; }
+               }
+            }) // WALBERLA_FOR_ALL_CELLS_OMP
+         }
+         else
+         {
+            WALBERLA_FOR_ALL_CELLS_OMP(flagFieldIt, flagField, pdfFieldIt, pdfField, fillFieldIt, fillField,
+                                       omp parallel for schedule(static) reduction(+:mass),
+            {
+               if (flagInfo.isLiquid(flagFieldIt) || flagInfo.isInterface(flagFieldIt))
+               {
+                  //TODO Implement
+                  const real_t density = 1.0; // pdfField->getDensity(pdfFieldIt.cell());
+                  mass += *fillFieldIt * density;
+               }
+            }) // WALBERLA_FOR_ALL_CELLS_OMP
+         }
+      }
+
+      mpi::allReduceInplace< real_t >(mass, mpi::SUM);
+      *totalMass_ = mass;
+
+      if (excessMass_ != nullptr)
+      {
+         mpi::allReduceInplace< real_t >(excessMass, mpi::SUM);
+         *excessMass_ = excessMass;
+      }
+   };
+
+ private:
+   std::weak_ptr< const StructuredBlockForest > blockForest_;
+   ConstBlockDataID flagFieldID_;
+   FlagInfo<FlagField_T> flagInfo_;
+
+   const ConstBlockDataID pdfFieldID_;
+   const ConstBlockDataID fillFieldID_;
+   const ConstBlockDataID excessMassFieldID_ = ConstBlockDataID();
+
+   std::shared_ptr< real_t > totalMass_;
+   std::shared_ptr< real_t > excessMass_ = nullptr; // mass stored in the excessMassField
+
+   uint_t frequency_;
+   uint_t executionCounter_;
+}; // class TotalMassComputer
+
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file

src/lbm_generated/free_surface/VtkWriter.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 VtkWriter.h
+//! \ingroup free_surface
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Free surface-specific VTK writer function.
+//
+//======================================================================================================================
+
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "field/adaptors/AdaptorCreators.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/FlagFieldCellFilter.h"
+#include "field/vtk/FlagFieldMapping.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "stencil/D3Q27.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/Initialization.h"
+
+#include "FlagInfo.h"
+
+namespace walberla
+{
+namespace free_surface_generated
+{
+/***********************************************************************************************************************
+ * Add VTK output to time loop that includes all relevant free surface information. It must be configured via
+ * config-file.
+ **********************************************************************************************************************/
+template< typename FlagInfo_T, typename PdfField_T, typename FlagField_T,
+          typename ScalarField_T, typename VectorField_T,
+          typename VectorFieldFlattened_T = GhostLayerField< real_t, 3 > >
+void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr, SweepTimeloop& timeloop,
+                  const std::weak_ptr< Config >& configPtr,
+                  const FlagInfo_T& flagInfo, const BlockDataID& pdfFieldID,
+                  const BlockDataID& flagFieldID, const BlockDataID& fillFieldID,
+                  const BlockDataID& forceDensityFieldID, const BlockDataID& curvatureFieldID,
+                  const BlockDataID& normalFieldID, const BlockDataID& obstacleNormalFieldID)
+{
+   // TODO density and velocity are not added to VTK Output.
+   using value_type       = typename FlagField_T::value_type;
+   const auto blockForest = blockForestPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+   const auto config = configPtr.lock();
+   WALBERLA_CHECK_NOT_NULLPTR(config);
+
+   // define VTK output (see src/vtk/Initialization.cpp, line 574 for usage)
+   const auto vtkConfigFunc = [&](std::vector< std::shared_ptr< vtk::BlockCellDataWriterInterface > >& writers,
+                                  std::map< std::string, vtk::VTKOutput::CellFilter >& filters,
+                                  std::map< std::string, vtk::VTKOutput::BeforeFunction >& beforeFuncs) {
+      using field::VTKWriter;
+      writers.push_back(std::make_shared< VTKWriter< PdfField_T, float > >(pdfFieldID, "pdf"));
+      writers.push_back(std::make_shared< VTKWriter< FlagField_T, float > >(flagFieldID, "flag"));
+      writers.push_back(std::make_shared< VTKWriter< ScalarField_T, float > >(fillFieldID, "fill_level"));
+      writers.push_back(std::make_shared< VTKWriter< ScalarField_T, float > >(curvatureFieldID, "curvature"));
+      writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(normalFieldID, "normal"));
+      writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(obstacleNormalFieldID, "obstacle_normal"));
+
+      if (forceDensityFieldID != BlockDataID())
+      {
+         writers.push_back(std::make_shared< VTKWriter< VectorFieldFlattened_T, float > >(forceDensityFieldID, "force_density"));
+      }
+
+      // map flagIDs to integer values
+      const auto flagMapper =
+         std::make_shared< field::FlagFieldMapping< FlagField_T, value_type > >(flagFieldID, "mapped_flag");
+      flagMapper->addMapping(flagIDs::liquidFlagID, numeric_cast< value_type >(1));
+      flagMapper->addMapping(flagIDs::interfaceFlagID, numeric_cast< value_type >(2));
+      flagMapper->addMapping(flagIDs::gasFlagID, numeric_cast< value_type >(3));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::noSlipFlagID, numeric_cast< value_type >(4));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::freeSlipFlagID, numeric_cast< value_type >(6));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::ubbFlagID, numeric_cast< value_type >(6));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::ubbInflowFlagID, numeric_cast< value_type >(7));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::pressureFlagID, numeric_cast< value_type >(8));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::pressureOutflowFlagID, numeric_cast< value_type >(9));
+//      flagMapper->addMapping(FreeSurfaceBoundaryHandling_T::outletFlagID, numeric_cast< value_type >(10));
+
+      writers.push_back(flagMapper);
+
+      // filter for writing only liquid and interface cells to VTK
+      auto liquidInterfaceFilter = field::FlagFieldCellFilter< FlagField_T >(flagFieldID);
+      liquidInterfaceFilter.addFlag(flagIDs::liquidFlagID);
+      liquidInterfaceFilter.addFlag(flagIDs::interfaceFlagID);
+      filters["liquidInterfaceFilter"] = liquidInterfaceFilter;
+
+      // communicate fields to update the ghost layer
+      auto preVTKComm = blockforest::communication::UniformBufferedScheme< stencil::D3Q27 >(blockForest);
+      preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< PdfField_T > >(pdfFieldID));
+      preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< FlagField_T > >(flagFieldID));
+      preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< ScalarField_T > >(fillFieldID));
+      preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< ScalarField_T > >(curvatureFieldID));
+      preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< VectorField_T > >(normalFieldID));
+      preVTKComm.addPackInfo(
+         std::make_shared< field::communication::PackInfo< VectorField_T > >(obstacleNormalFieldID));
+      preVTKComm.addPackInfo(std::make_shared< field::communication::PackInfo< VectorFieldFlattened_T > >(forceDensityFieldID));
+
+      beforeFuncs["ghost_layer_synchronization"] = preVTKComm;
+
+      // set velocity and density to zero in obstacle and gas cells (only for visualization purposes); the PDF values in
+      // these cells are not important and thus not set during the simulation;
+      // only enable this functionality if the non-liquid and non-interface cells are not excluded anyway
+      const auto vtkConfigBlock        = config->getOneBlock("VTK");
+      const auto fluidFieldConfigBlock = vtkConfigBlock.getBlock("fluid_field");
+      if (fluidFieldConfigBlock)
+      {
+         auto inclusionFiltersConfigBlock = fluidFieldConfigBlock.getBlock("inclusion_filters");
+
+         // liquidInterfaceFilter limits VTK-output to only liquid and interface cells
+         if (!inclusionFiltersConfigBlock.isDefined("liquidInterfaceFilter"))
+         {
+            class ZeroSetter
+            {
+             public:
+               ZeroSetter(const weak_ptr< StructuredBlockForest >& blockForest, const BlockDataID& pdfFieldID,
+                          const ConstBlockDataID& flagFieldID,
+                          const FlagInfo_T& flagInfo)
+                  : blockForest_(blockForest), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID), flagInfo_(flagInfo)
+               {}
+
+               void operator()()
+               {
+                  auto blockForest = blockForest_.lock();
+                  WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+                  for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+                  {
+                     PdfField_T* const pdfField         = blockIt->template getData< PdfField_T >(pdfFieldID_);
+                     const FlagField_T* const flagField = blockIt->template getData< const FlagField_T >(flagFieldID_);
+                     WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ(pdfField, uint_c(1), {
+                        const typename PdfField_T::Ptr pdfFieldPtr(*pdfField, x, y, z);
+                        const typename FlagField_T::ConstPtr flagFieldPtr(*flagField, x, y, z);
+
+//                        if (flagInfo_.isGas(*flagFieldPtr) || flagInfo_.isObstacle(*flagFieldPtr))
+//                        {
+//                           pdfField->setDensityAndVelocity(pdfFieldPtr.cell(), Vector3< real_t >(real_c(0)),
+//                                                           real_c(1.0));
+//                        }
+                     }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ
+                  }
+               }
+
+             private:
+               weak_ptr< StructuredBlockForest > blockForest_;
+               BlockDataID pdfFieldID_;
+               ConstBlockDataID flagFieldID_;
+               FlagInfo_T flagInfo_;
+            };
+
+            beforeFuncs["gas_cell_zero_setter"] = ZeroSetter(blockForest, pdfFieldID, flagFieldID, flagInfo);
+         }
+      }
+   };
+
+   // add VTK output to timeloop
+   std::map< std::string, vtk::SelectableOutputFunction > vtkOutputFunctions;
+   vtk::initializeVTKOutput(vtkOutputFunctions, vtkConfigFunc, blockForest, config);
+   for (auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output)
+   {
+      timeloop.addFuncBeforeTimeStep(output->second.outputFunction, std::string("VTK: ") + output->first,
+                                     output->second.requiredGlobalStates, output->second.incompatibleGlobalStates);
+   }
+}
+
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file