diff --git a/apps/showcases/FreeSurface/CMakeLists.txt b/apps/showcases/FreeSurface/CMakeLists.txt
index ad5828ac44a34cbe280895b2b88cb7b5b59bfbd8..4f5f1c389c010a50518d91b7ea467e10ea8e138f 100644
--- a/apps/showcases/FreeSurface/CMakeLists.txt
+++ b/apps/showcases/FreeSurface/CMakeLists.txt
@@ -1,5 +1,9 @@
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)
diff --git a/apps/showcases/FreeSurface/DropInPool.cpp b/apps/showcases/FreeSurface/DropInPool.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..f2dbd21c1ad3b48e590939dbcda4c07547beaf0c
--- /dev/null
+++ b/apps/showcases/FreeSurface/DropInPool.cpp
@@ -0,0 +1,284 @@
+//======================================================================================================================
+//
+// 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
diff --git a/apps/showcases/FreeSurface/DropInPool.prm b/apps/showcases/FreeSurface/DropInPool.prm
new file mode 100644
index 0000000000000000000000000000000000000000..1b28076e1a0b58e75a4b8034576394f669075de8
--- /dev/null
+++ b/apps/showcases/FreeSurface/DropInPool.prm
@@ -0,0 +1,37 @@
+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;
+ }
+ }
+}
diff --git a/src/lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h b/src/lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h
index e07f890fb362fa2b92471c4c4081bc43c194b80a..3e3bbdb48b88b04abb5df306d5b9a397dc092df2 100644
--- a/src/lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h
+++ b/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: