Add falling drop showcase

apps/showcases/FreeSurface/CMakeLists.txt

 waLBerla_link_files_to_builddir( *.prm )
 
+waLBerla_add_executable(NAME DropInPool
+                        FILES   DropInPool.cpp
+                        DEPENDS blockforest boundary core domain_decomposition field lbm timeloop vtk)
+
 waLBerla_add_executable(NAME    RisingBubble
                         FILES   RisingBubble.cpp
                         DEPENDS blockforest boundary core domain_decomposition field lbm timeloop vtk)

apps/showcases/FreeSurface/DropInPool.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file DropInPool.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+
+#include "field/adaptors/AdaptorCreators.h"
+#include "field/vtk/FlagFieldCellFilter.h"
+#include "field/vtk/FlagFieldMapping.h"
+#include "field/vtk/VTKWriter.h"
+
+#include "lbm/PerformanceLogger.h"
+#include "lbm/field/Adaptors.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/bubble_model/Geometry.h"
+#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
+#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
+#include "lbm/lattice_model/D3Q19.h"
+
+#include "vtk/Initialization.h"
+
+using namespace walberla;
+
+using CollisionModel_T = lbm::collision_model::SRT;
+using LatticeModel_T   = lbm::D3Q19< CollisionModel_T, true, lbm::force_model::GuoConstant, 2 >;
+using Stencil_T        = LatticeModel_T::Stencil;
+using PdfField_T       = lbm::PdfField< LatticeModel_T >;
+
+using Comm = blockforest::SimpleCommunication< LatticeModel_T::CommunicationStencil >;
+
+using flag_t      = uint16_t;
+using FlagField_T = FlagField< flag_t >;
+using FsBoundary  = free_surface::FreeSurfaceBoundary< LatticeModel_T, FlagField_T >;
+
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+
+int main(int argc, char** argv)
+{
+   Environment walberlaEnv(argc, argv);
+
+   if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
+
+   // read domain parameters from parameter file
+   auto domParam                            = walberlaEnv.config()->getOneBlock("DomainParameters");
+   const Vector3< uint_t > numBlocks        = domParam.getParameter< Vector3< uint_t > >("blocks");
+   const uint_t dropSize                    = domParam.getParameter< uint_t >("dropSize");
+   const Vector3< uint_t > domainSizeFactor = domParam.getParameter< Vector3< uint_t > >("domainSizeFactor");
+   const Vector3< uint_t > domainSize       = domainSizeFactor * dropSize;
+   const Vector3< bool > periodic           = domParam.getParameter< Vector3< bool > >("periodic");
+   const Vector3< uint_t > cellsPerBlock{ uint_c(ceil(real_c(domainSize[0]) / real_c(numBlocks[0]))),
+                                          uint_c(ceil(real_c(domainSize[1]) / real_c(numBlocks[1]))),
+                                          uint_c(ceil(real_c(domainSize[2]) / real_c(numBlocks[2]))) };
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(dropSize);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodic);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
+
+   // create uniform block grid
+   auto blocks = blockforest::createUniformBlockGrid(numBlocks[0], numBlocks[1], numBlocks[2], cellsPerBlock[0],
+                                                     cellsPerBlock[1], cellsPerBlock[2], real_c(1), true, periodic[0],
+                                                     periodic[1], periodic[2], true);
+
+   // read simulation parameters from parameter file
+   auto simParam                 = walberlaEnv.config()->getOneBlock("SimulationParameters");
+   const real_t omega            = simParam.getParameter< real_t >("omega");
+   const Vector3< real_t > force = simParam.getParameter< Vector3< real_t > >("force");
+
+   const uint_t timesteps       = simParam.getParameter< uint_t >("timesteps");
+   const uint_t perfLogInterval = simParam.getParameter< uint_t >("perfLogInterval");
+
+   const CollisionModel_T collisionModel = lbm::collision_model::SRT(omega);
+   const real_t viscosity                = collisionModel.viscosity();
+
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(omega);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
+   WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
+
+   // create lattice model
+   LatticeModel_T latticeModel = LatticeModel_T(collisionModel, lbm::force_model::GuoConstant(force));
+
+   // add pdf field
+   const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blocks, "PDF field", latticeModel);
+
+   // add fill level field (initialized with 0, i.e., gas everywhere)
+   const BlockDataID fillFieldID = field::addToStorage< ScalarField_T >(blocks, "Fill level field", real_c(0.0));
+
+   // add boundary handling
+   FsBoundary boundary(blocks, pdfFieldID, fillFieldID);
+   const BlockDataID flagFieldID           = boundary.getFlagFieldID();
+   const typename FsBoundary::Flags& flags = boundary.getFlagInfo();
+
+   // add various adaptors (for simplified access to macroscopic quantities)
+   const BlockDataID densityAdaptorID =
+      field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::Density >(blocks, pdfFieldID, "DensityAdaptor");
+   const BlockDataID velocityAdaptorID =
+      field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::VelocityVector >(blocks, pdfFieldID,
+                                                                                        "VelocityVectorAdaptor");
+   field::addFieldAdaptor< typename lbm::Adaptor< LatticeModel_T >::Velocity >(blocks, pdfFieldID, "VelocityAdaptor");
+
+   // initialize a pool of liquid at the bottom of the domain in z-direction
+   const AABB poolAABB(real_c(0), real_c(0), real_c(0), real_c(domainSize[0]), real_c(domainSize[1]),
+                       real_c(domainSize[2]) * real_c(0.3));
+
+   for (auto blockIterator = blocks->begin(); blockIterator != blocks->end(); ++blockIterator)
+   {
+      ScalarField_T* fillField = blockIterator->getData< ScalarField_T >(fillFieldID);
+
+      // determine the cells that are relevant for a block (still in global coordinates)
+      const CellInterval relevantCellBB = blocks->getCellBBFromAABB(poolAABB.getIntersection(blockIterator->getAABB()));
+
+      // transform the global coordinates of relevant cells to block local coordinates
+      CellInterval blockLocalCellBB;
+      blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, *blockIterator, relevantCellBB);
+
+      // clang-format off
+         WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(blockLocalCellBB,
+                                                fillField->get(x,y,z) = real_c(1);
+         ) // WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ
+      // clang-format on
+   }
+
+   // initialize a cubical drop of liquid in the center of the domain
+   const real_t dropSizeHalf = real_c(dropSize) * real_c(0.5);
+   const AABB dropAABB(
+      real_c(domainSize[0]) * real_c(0.5) - dropSizeHalf, real_c(domainSize[1]) * real_c(0.5) - dropSizeHalf,
+      real_c(domainSize[2]) * real_c(0.75) - dropSizeHalf, real_c(domainSize[0]) * real_c(0.5) + dropSizeHalf,
+      real_c(domainSize[1]) * real_c(0.5) + dropSizeHalf, real_c(domainSize[2]) * real_c(0.75) + dropSizeHalf);
+
+   for (auto blockIterator = blocks->begin(); blockIterator != blocks->end(); ++blockIterator)
+   {
+      ScalarField_T* fillField = blockIterator->getData< ScalarField_T >(fillFieldID);
+
+      // determine the cells that are relevant for a block (still in global coordinates)
+      CellInterval relevantCellBB = blocks->getCellBBFromAABB(dropAABB.getIntersection(blockIterator->getAABB()));
+
+      // transform the global coordinates of relevant cells to block local coordinates
+      CellInterval blockLocalCellBB;
+      blocks->transformGlobalToBlockLocalCellInterval(blockLocalCellBB, *blockIterator, relevantCellBB);
+
+      // clang-format off
+            WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ(blockLocalCellBB,
+                                                   fillField->get(x,y,z) = real_c(1);
+            ) // WALBERLA_FOR_ALL_CELLS_IN_INTERVAL_XYZ
+      // clang-format on
+   }
+
+   // initialize boundary conditions from config file
+   auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
+   boundary.initFromConfig(boundaryParameters);
+   boundary.initFlagsFromFillLevel();
+
+   // initial communication
+   Comm(blocks, pdfFieldID, fillFieldID, flagFieldID)();
+
+   // add bubble model
+   const bool disableSplits = true; // necessary if a gas bubble could split
+   auto bubbleModel         = make_shared< bubble_model::BubbleModel >(blocks, disableSplits);
+   bubbleModel->initFromFillLevelField(fillFieldID);
+
+   // create timeloop
+   SweepTimeloop tl(blocks, simParam.getParameter< uint_t >("timesteps"));
+
+   // add surface geometry handler
+   free_surface::SurfaceGeometryHandler< FlagField_T, ScalarField_T, VectorField_T > sgh(blocks, flagFieldID,
+                                                                                         fillFieldID);
+   sgh.setBoundaryFlagUIDs(FsBoundary::noSlipFlagID + FsBoundary::pressureFlag1ID + FsBoundary::pressureFlag2ID +
+                           FsBoundary::outletFlagID + FsBoundary::ubbFlagID + FsBoundary::freeSlipFlagID);
+   sgh.setContactAngle(simParam.getParameter< real_t >("contactAngle", 90));
+
+   ConstBlockDataID curvatureField = sgh.getCurvatureFieldID();
+   ConstBlockDataID normalField    = sgh.getNormalFieldID();
+
+   sgh.addSweeps(tl);
+
+   // add boundary handling for standard boundaries and free surface boundaries
+   free_surface::SurfaceDynamicsHandler< LatticeModel_T, FlagField_T > dynamicsHandler(
+      blocks, pdfFieldID, flagFieldID, fillFieldID, normalField, curvatureField, boundary, bubbleModel, real_c(0));
+
+   dynamicsHandler.addSweeps(tl);
+
+   // add VTK output
+   auto vtkConfigFunc = [&](std::vector< 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(
+         make_shared< VTKWriter< typename lbm::Adaptor< LatticeModel_T >::Density > >(densityAdaptorID, "density"));
+      writers.push_back(make_shared< VTKWriter< typename lbm::Adaptor< LatticeModel_T >::VelocityVector > >(
+         velocityAdaptorID, "velocity"));
+      writers.push_back(make_shared< VTKWriter< PdfField_T, float > >(pdfFieldID, "PDFs"));
+      writers.push_back(make_shared< VTKWriter< ScalarField_T, float > >(fillFieldID, "fill_levels"));
+      writers.push_back(make_shared< VTKWriter< ScalarField_T, float > >(curvatureField, "curvature"));
+      writers.push_back(make_shared< VTKWriter< VectorField_T, float > >(normalField, "normals"));
+      writers.push_back(make_shared< VTKWriter< FlagField_T, float > >(flagFieldID, "flags"));
+
+      auto flagMapper =
+         make_shared< field::FlagFieldMapping< FlagField_T, walberla::uint8_t > >(flagFieldID, "mapped_flag_field");
+      flagMapper->addMapping(FsBoundary::Flags::interfaceID, 1);
+      flagMapper->addMapping(FsBoundary::Flags::liquidID, 2);
+      flagMapper->addMapping(FsBoundary::Flags::gasID, 3);
+      flagMapper->addMapping(FsBoundary::noSlipFlagID, 4);
+      flagMapper->addMapping(FsBoundary::pressureFlag1ID, 5);
+      flagMapper->addMapping(FsBoundary::pressureFlag2ID, 6);
+      flagMapper->addMapping(FsBoundary::ubbFlagID, 7);
+
+      writers.push_back(flagMapper);
+
+      auto liquidInterfaceFilter = field::FlagFieldCellFilter< FlagField_T >(flagFieldID);
+      liquidInterfaceFilter.addFlag(FsBoundary::Flags::liquidID);
+      liquidInterfaceFilter.addFlag(FsBoundary::Flags::interfaceID);
+      filters["liquidInterfaceFilter"] = liquidInterfaceFilter;
+
+      auto preVTKComm = blockforest::communication::UniformBufferedScheme< stencil::D3Q27 >(blocks);
+      preVTKComm.addPackInfo(make_shared< field::communication::PackInfo< PdfField_T > >(pdfFieldID));
+      preVTKComm.addPackInfo(make_shared< field::communication::PackInfo< FlagField_T > >(flagFieldID));
+      preVTKComm.addPackInfo(make_shared< field::communication::PackInfo< ScalarField_T > >(fillFieldID));
+
+      beforeFuncs["sync"] = preVTKComm;
+   };
+   std::map< std::string, vtk::SelectableOutputFunction > vtkOutputFunctions;
+
+   vtk::initializeVTKOutput(vtkOutputFunctions, vtkConfigFunc, blocks, walberlaEnv.config());
+
+   for (auto output = vtkOutputFunctions.begin(); output != vtkOutputFunctions.end(); ++output)
+   {
+      tl.addFuncBeforeTimeStep(output->second.outputFunction, std::string("VTK: ") + output->first,
+                               output->second.requiredGlobalStates, output->second.incompatibleGlobalStates);
+   }
+
+   // set velocity and density to zero in obstacle and gas cells (only for visualization purposes)
+   // the PDF values in these fields are not important and thus not automatically set
+   bool enableZeroSetter = true;
+   if (enableZeroSetter)
+   {
+      auto function = [&](IBlock* block) {
+         using namespace free_surface;
+         PdfField_T* pdfFieldPtr   = block->getData< PdfField_T >(pdfFieldID);
+         FlagField_T* flagFieldPtr = block->getData< FlagField_T >(flagFieldID);
+         WALBERLA_FOR_ALL_CELLS(
+            pdfIter, pdfFieldPtr, flagIter, flagFieldPtr,
+            if (flags.isGas(flagIter) || flags.isObstacle(flagIter))
+               pdfFieldPtr->setDensityAndVelocity(pdfIter.cell(), Vector3< real_t >(0), real_c(1.0));)
+      };
+      tl.add() << Sweep(function, "Zerosetter");
+   }
+
+   // add logging for computational performance
+   lbm::PerformanceLogger< FlagField_T > perfLogger(blocks, flagFieldID, FsBoundary::Flags::liquidInterfaceIDs,
+                                                    perfLogInterval);
+   tl.addFuncAfterTimeStep(perfLogger, "PerformanceLogger");
+
+   WcTimingPool timingPool;
+
+   for (uint_t t = 0; t != timesteps; ++t)
+   {
+      if (t % uint_c(100) == uint_c(0)) { WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep=" << t); }
+      tl.singleStep(timingPool, true);
+
+      if (t % perfLogInterval == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
+   }
+
+   return EXIT_SUCCESS;
+}
\ No newline at end of file

apps/showcases/FreeSurface/DropInPool.prm

+DomainParameters
+{
+   blocks       < 1, 1, 4 >;
+   dropSize 10;
+   domainSizeFactor < 4, 4, 4 >; // values multiplied with dropDiameter
+   periodic     < 0, 0, 0 >;
+}
+
+SimulationParameters
+{
+   omega 1.8;
+   force <0, 0, -1e-4>;
+	initialVelocity < 0, 0, 0 >;
+	timesteps       4000;
+	perfLogInterval 1000;
+}
+
+BoundaryParameters
+{
+	Border { direction all;  walldistance -1;  NoSlip{} }
+}
+
+VTK {
+   fluid_field {
+      writeFrequency 50;
+      samplingResolution 1;
+      writers {
+         velocity;
+         density;
+         fill_levels;
+         curvature;
+         normals;
+         flags;
+         mapped_flag_field;
+      }
+   }
+}

src/lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h

@@ -93,26 +93,27 @@ class SurfaceGeometryHandler
                << Sweep(emptySweep, "Deactivated normals sweep")
                << AfterFunction(Comm(bf_, normalFieldID_), "Communication after normals sweep");
 
-      // curvature computation using local triangulation
-      CurvatureSweepTR< FlagField_T, ScalarField_T, VectorField_T > curvSweep(
-         curvatureFieldID_, normalFieldID_, fillFieldID_, flagFieldID_, interfaceUIDs_, boundaryUIDs_, contactAngle_);
-      tl.add() << Sweep(curvSweep, "Curvature sweep", StateSweep::fullFreeSurface)
-               << Sweep(emptySweep, "Deactivated curvature sweep")
-               << AfterFunction(Comm(bf_, curvatureFieldID_), "Communication after curvature sweep");
-
-      //      // extrapolation of normals to interface neighboring cells (required for computing the curvature)
-      //      ExtrapolateNormalsSweep< FlagField_T, VectorField_T > extNormalsSweep(normalFieldID_, flagFieldID_,
-      //                                                                            interfaceUIDs_);
-      //      tl.add() << Sweep(extNormalsSweep, "Extrapolate normals sweep", StateSweep::fullFreeSurface)
-      //               << Sweep(emptySweep, "Deactivated extrapolate normals sweep")
-      //               << AfterFunction(Comm(bf_, normalFieldID_), "Communication after extrapolate normals sweep");
-
-      //       // curvature computation using finite differences
-      //      CurvatureSweepFD< FlagField_T, ScalarField_T, VectorField_T > curvSweep(
-      //         curvatureFieldID_, normalFieldID_, flagFieldID_, interfaceUIDs_, boundaryUIDs_, contactAngle_);
+      //      // curvature computation using local triangulation
+      //      CurvatureSweepTR< FlagField_T, ScalarField_T, VectorField_T > curvSweep(
+      //         curvatureFieldID_, normalFieldID_, fillFieldID_, flagFieldID_, interfaceUIDs_, boundaryUIDs_,
+      //         contactAngle_);
       //      tl.add() << Sweep(curvSweep, "Curvature sweep", StateSweep::fullFreeSurface)
       //               << Sweep(emptySweep, "Deactivated curvature sweep")
       //               << AfterFunction(Comm(bf_, curvatureFieldID_), "Communication after curvature sweep");
+
+      // extrapolation of normals to interface neighboring cells (required for computing the curvature)
+      ExtrapolateNormalsSweep< FlagField_T, VectorField_T > extNormalsSweep(normalFieldID_, flagFieldID_,
+                                                                            interfaceUIDs_);
+      tl.add() << Sweep(extNormalsSweep, "Extrapolate normals sweep", StateSweep::fullFreeSurface)
+               << Sweep(emptySweep, "Deactivated extrapolate normals sweep")
+               << AfterFunction(Comm(bf_, normalFieldID_), "Communication after extrapolate normals sweep");
+
+      // curvature computation using finite differences
+      CurvatureSweepFD< FlagField_T, ScalarField_T, VectorField_T > curvSweep(
+         curvatureFieldID_, normalFieldID_, flagFieldID_, interfaceUIDs_, boundaryUIDs_, contactAngle_);
+      tl.add() << Sweep(curvSweep, "Curvature sweep", StateSweep::fullFreeSurface)
+               << Sweep(emptySweep, "Deactivated curvature sweep")
+               << AfterFunction(Comm(bf_, curvatureFieldID_), "Communication after curvature sweep");
    }
 
  private: