diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index 2adab475060437d4b2934f215c45d470b3b34b0a..3f5e6a95a0f0c6a8199de013da440e37e7f78afb 100644
--- a/apps/benchmarks/CMakeLists.txt
+++ b/apps/benchmarks/CMakeLists.txt
@@ -4,6 +4,7 @@ add_subdirectory( ComplexGeometry )
add_subdirectory( DEM )
add_subdirectory( MeshDistance )
add_subdirectory( CouetteFlow )
+add_subdirectory( FreeSurfaceAdvection )
add_subdirectory( FluidParticleCoupling )
add_subdirectory( FluidParticleCouplingWithLoadBalancing )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
diff --git a/apps/benchmarks/FreeSurfaceAdvection/CMakeLists.txt b/apps/benchmarks/FreeSurfaceAdvection/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..e14b5fe37fda409f4e80790593e378c8f20b98f5
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/CMakeLists.txt
@@ -0,0 +1,13 @@
+waLBerla_link_files_to_builddir( *.prm )
+
+waLBerla_add_executable(NAME DeformationField
+ FILES DeformationField.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME SingleVortex
+ FILES SingleVortex.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
+waLBerla_add_executable(NAME ZalesakDisk
+ FILES ZalesakDisk.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
\ No newline at end of file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp b/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..fcbc4ca54fe71364c93c275000f9b210276a56c7
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp
@@ -0,0 +1,351 @@
+//======================================================================================================================
+//
+// 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 DeformationField.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+// This benchmark simulates the advection of a spherical bubble in a vortex field with very high deformation. The vortex
+// field changes periodically so that the bubble returns to its initial position, where it should take its initial form.
+// The relative geometrical error of the bubble's shape after one period is evaluated. There is no LBM flow simulation
+// performed here. It is a test case for the FSLBM's mass advection only. This benchmark is based on Viktor Haag's
+// master thesis (https://www10.cs.fau.de/publications/theses/2017/Haag_MT_2017.pdf).
+//======================================================================================================================
+
+#include "core/Environment.h"
+
+#include "field/Gather.h"
+
+#include "lbm/PerformanceLogger.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/LoadBalancing.h"
+#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
+#include "lbm/free_surface/VtkWriter.h"
+#include "lbm/free_surface/bubble_model/Geometry.h"
+#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
+#include "lbm/free_surface/surface_geometry/Utility.h"
+#include "lbm/lattice_model/D3Q19.h"
+
+#include "functionality/AdvectionDynamicsHandler.h"
+#include "functionality/GeometricalErrorEvaluator.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+namespace DeformationField
+{
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+
+// Lattice model is only created for dummy purposes; no LBM simulation is performed
+using CollisionModel_T = lbm::collision_model::SRT;
+using LatticeModel_T = lbm::D3Q19< CollisionModel_T, true >;
+using LatticeModelStencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+using PdfCommunication_T = blockforest::SimpleCommunication< LatticeModelStencil_T >;
+
+// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
+// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
+using CommunicationStencil_T =
+ typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
+using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
+
+using flag_t = uint32_t;
+using FlagField_T = FlagField< flag_t >;
+using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
+
+// function describing the initialization velocity profile (in global cell coordinates)
+inline Vector3< real_t > velocityProfile(Cell globalCell, real_t timePeriod, uint_t timestep,
+ const Vector3< real_t >& domainSize)
+{
+ // add 0.5 to get the cell's center
+ const real_t x = (real_c(globalCell.x()) + real_c(0.5)) / domainSize[0];
+ const real_t y = (real_c(globalCell.y()) + real_c(0.5)) / domainSize[1];
+ const real_t z = (real_c(globalCell.z()) + real_c(0.5)) / domainSize[2];
+
+ const real_t timeTerm = real_c(std::cos(math::pi * real_t(timestep) / timePeriod));
+
+ const real_t sinpix = real_c(std::sin(math::pi * x));
+ const real_t sinpiy = real_c(std::sin(math::pi * y));
+ const real_t sinpiz = real_c(std::sin(math::pi * z));
+
+ const real_t sin2pix = real_c(std::sin(real_c(2) * math::pi * x));
+ const real_t sin2piy = real_c(std::sin(real_c(2) * math::pi * y));
+ const real_t sin2piz = real_c(std::sin(real_c(2) * math::pi * z));
+
+ const real_t velocityX = real_c(2) * sinpix * sinpix * sin2piy * sin2piz * timeTerm;
+ const real_t velocityY = -sin2pix * sinpiy * sinpiy * sin2piz * timeTerm;
+ const real_t velocityZ = -sin2pix * sin2piy * sinpiz * sinpiz * timeTerm;
+
+ return Vector3< real_t >(velocityX, velocityY, velocityZ);
+}
+
+int main(int argc, char** argv)
+{
+ Environment walberlaEnv(argc, argv);
+
+ if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
+
+ // print content of parameter file
+ WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
+
+ // get block forest parameters from parameter file
+ auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
+ const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+ const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
+
+ // get domain parameters from parameter file
+ auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
+ const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
+
+ const real_t bubbleDiameter = real_c(domainWidth) * real_c(0.3);
+ const Vector3< real_t > bubbleCenter = domainWidth * Vector3< real_t >(real_c(0.35), real_c(0.35), real_c(0.35));
+
+ // define domain size
+ Vector3< uint_t > domainSize;
+ domainSize[0] = domainWidth;
+ domainSize[1] = domainWidth;
+ domainSize[2] = domainWidth;
+
+ // compute number of blocks as defined by domainSize and cellsPerBlock
+ Vector3< uint_t > numBlocks;
+ numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
+ numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
+ numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
+
+ // get number of (MPI) processes
+ const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
+ WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
+ "The number of MPI processes is greater than the number of blocks as defined by "
+ "\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
+ "the number of MPI processes or increase \"cellsPerBlock\".")
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainWidth);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleDiameter);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleCenter);
+
+ // get physics parameters from parameter file
+ auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
+ const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+ const uint_t timestepsToInitialPosition = physicsParameters.getParameter< uint_t >("timestepsToInitialPosition");
+
+ // compute CFL number
+ const real_t dx_SI = real_c(1) / real_c(domainWidth);
+ const real_t dt_SI = real_c(3) / real_c(timestepsToInitialPosition);
+ const real_t CFL = dt_SI / dx_SI; // with velocity_SI = 1
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(CFL);
+
+ // dummy collision model (LBM not simulated in this benchmark)
+ const CollisionModel_T collisionModel = CollisionModel_T(real_c(1));
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timestepsToInitialPosition);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
+
+ // read model parameters from parameter file
+ const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
+ const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
+ const std::string excessMassDistributionModel =
+ modelParameters.getParameter< std::string >("excessMassDistributionModel");
+ const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
+ const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
+ const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
+ const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
+
+ // read evaluation parameters from parameter file
+ const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
+ const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
+
+ // create non-uniform block forest (non-uniformity required for load balancing)
+ const std::shared_ptr< StructuredBlockForest > blockForest =
+ createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
+
+ // create lattice model
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel);
+
+ // add pdf field
+ const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
+
+ // add fill level field (initialized with 1, i.e., liquid everywhere)
+ const BlockDataID fillFieldID =
+ field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
+
+ // add boundary handling
+ const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
+ std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
+ const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
+ const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
+
+ // initialize the velocity profile
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+ FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
+ // cell in block-local coordinates
+ const Cell localCell = pdfFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = pdfFieldIt.cell();
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
+ const Vector3< real_t > initialVelocity =
+ CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), uint_c(0), domainSize);
+ pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ // create the spherical bubble
+ const geometry::Sphere sphereBubble(real_c(domainWidth) * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25)),
+ real_c(domainWidth) * real_c(0.15));
+ bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
+
+ // initialize domain boundary conditions from config file
+ const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
+ freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
+
+ // IMPORTANT REMARK: this must be called only after every solid flag has been set; otherwise, the boundary handling
+ // might not detect solid flags correctly
+ freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
+
+ // communication after initialization
+ Communication_T communication(blockForest, flagFieldID, fillFieldID);
+ communication();
+
+ PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
+ pdfCommunication();
+
+ const ConstBlockDataID initialFillFieldID =
+ field::addCloneToStorage< ScalarField_T >(blockForest, fillFieldID, "Initial fill level field");
+
+ // add bubble model
+ const std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel =
+ std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1));
+
+ // create timeloop
+ SweepTimeloop timeloop(blockForest, timesteps);
+
+ // add surface geometry handler
+ const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
+ blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, false, false, real_c(0));
+
+ geometryHandler.addSweeps(timeloop);
+
+ // get fields created by surface geometry handler
+ const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
+
+ // add boundary handling for standard boundaries and free surface boundaries
+ const AdvectionDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, normalFieldID, freeSurfaceBoundaryHandling, bubbleModel,
+ pdfReconstructionModel, excessMassDistributionModel, useSimpleMassExchange, cellConversionThreshold,
+ cellConversionForceThreshold);
+
+ dynamicsHandler.addSweeps(timeloop);
+
+ // add evaluator for geometrical
+ const std::shared_ptr< real_t > geometricalError = std::make_shared< real_t >(real_c(0));
+ const GeometricalErrorEvaluator< FreeSurfaceBoundaryHandling_T, FlagField_T, ScalarField_T >
+ geometricalErrorEvaluator(blockForest, freeSurfaceBoundaryHandling, initialFillFieldID, fillFieldID,
+ evaluationFrequency, geometricalError);
+ timeloop.addFuncAfterTimeStep(geometricalErrorEvaluator, "Evaluator: geometrical error");
+
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
+ // add VTK output
+ addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(),
+ geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
+ geometryHandler.getObstNormalFieldID());
+
+ // add triangle mesh output of free surface
+ SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
+ blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
+ surfaceMeshWriter(); // write initial mesh
+ timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
+
+ // add logging for computational performance
+ const lbm::PerformanceLogger< FlagField_T > performanceLogger(
+ blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs, performanceLogFrequency);
+ timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluator: performance logging");
+
+ WcTimingPool timingPool;
+
+ for (uint_t t = uint_c(0); t != timesteps; ++t)
+ {
+ timeloop.singleStep(timingPool, true);
+
+ if (t % evaluationFrequency == uint_c(0))
+ {
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass << "\n\t\texcess mass = "
+ << *excessMass << "\n\t\tgeometrical error = " << *geometricalError);
+ }
+
+ // set the constant velocity profile
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+ FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
+ const Cell localCell = pdfFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = pdfFieldIt.cell();
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
+ const Vector3< real_t > velocity =
+ CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), t, domainSize);
+ pdfField->setDensityAndVelocity(localCell, velocity, real_c(1));
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ pdfCommunication();
+
+ if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace DeformationField
+} // namespace free_surface
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::free_surface::DeformationField::main(argc, argv); }
\ No newline at end of file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/DeformationField.prm b/apps/benchmarks/FreeSurfaceAdvection/DeformationField.prm
new file mode 100644
index 0000000000000000000000000000000000000000..359fd4a98093e95ace99e755819f022cca40e5c9
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/DeformationField.prm
@@ -0,0 +1,96 @@
+BlockForestParameters
+{
+ cellsPerBlock < 25, 25, 25 >;
+ periodicity < 1, 1, 1 >;
+}
+
+DomainParameters
+{
+ domainWidth 50;
+}
+
+PhysicsParameters
+{
+ timestepsToInitialPosition 3000;
+ timesteps 3001;
+
+}
+
+ModelParameters
+{
+ pdfReconstructionModel OnlyMissing;
+ excessMassDistributionModel EvenlyAllInterface;
+ curvatureModel FiniteDifferenceMethod;
+ useSimpleMassExchange false;
+ cellConversionThreshold 1e-2;
+ cellConversionForceThreshold 1e-1;
+}
+
+EvaluationParameters
+{
+ evaluationFrequency 300;
+ performanceLogFrequency 10000;
+}
+
+BoundaryParameters
+{
+ // X
+ //Border { direction W; walldistance -1; FreeSlip{} }
+ //Border { direction E; walldistance -1; FreeSlip{} }
+
+ // Y
+ //Border { direction N; walldistance -1; FreeSlip{} }
+ //Border { direction S; walldistance -1; FreeSlip{} }
+
+ // Z
+ //Border { direction T; walldistance -1; FreeSlip{} }
+ //Border { direction B; walldistance -1; FreeSlip{} }
+}
+
+MeshOutputParameters
+{
+ writeFrequency 300;
+ baseFolder mesh-out;
+}
+
+VTK
+{
+ fluid_field
+ {
+ writeFrequency 300;
+ ghostLayers 0;
+ baseFolder vtk-out;
+ samplingResolution 1;
+
+ writers
+ {
+ fill_level;
+ mapped_flag;
+ velocity;
+ density;
+ //curvature;
+ //normal;
+ //obstacle_normal;
+ //pdf;
+ //flag;
+ }
+
+ inclusion_filters
+ {
+ //liquidInterfaceFilter; // only include liquid and interface cells in VTK output
+ }
+
+ before_functions
+ {
+ //ghost_layer_synchronization; // only needed if writing the ghost layer
+ gas_cell_zero_setter; // sets velocity=0 and density=1 all gas cells before writing VTK output
+ }
+ }
+
+ domain_decomposition
+ {
+ writeFrequency 0;
+ baseFolder vtk-out;
+ outputDomainDecomposition true;
+ }
+}
\ No newline at end of file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp b/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d3c8e01ed776fac44d0b33f7cc7e94a6508c08ed
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp
@@ -0,0 +1,351 @@
+//======================================================================================================================
+//
+// 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 SingleVortex.cpp
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+// This benchmark simulates the advection of a spherical bubble in a vortex field.
+// The vortex field changes periodically so that the bubble returns to its initial position, where it should take its
+// initial form. The relative geometrical error of the bubble's shape after one period is evaluated. There is no LBM
+// flow simulation performed here, it is a test case for the FSLBM's mass advection. This benchmark is based on the
+// two-dimensional test case from Rider and Kothe (doi: 10.1006/jcph.1998.5906). The extension to three-dimensional
+// space is based on the Viktor Haag's master thesis
+// (https://www10.cs.fau.de/publications/theses/2017/Haag_MT_2017.pdf).
+//======================================================================================================================
+
+#include "core/Environment.h"
+
+#include "field/Gather.h"
+
+#include "lbm/PerformanceLogger.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/LoadBalancing.h"
+#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
+#include "lbm/free_surface/VtkWriter.h"
+#include "lbm/free_surface/bubble_model/Geometry.h"
+#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
+#include "lbm/free_surface/surface_geometry/Utility.h"
+#include "lbm/lattice_model/D3Q19.h"
+
+#include "functionality/AdvectionDynamicsHandler.h"
+#include "functionality/GeometricalErrorEvaluator.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+namespace SingleVortex
+{
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+
+// Lattice model is only created for dummy purposes; no LBM simulation is performed
+using CollisionModel_T = lbm::collision_model::SRT;
+using LatticeModel_T = lbm::D3Q19< CollisionModel_T, true >;
+using LatticeModelStencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField< LatticeModel_T >;
+using PdfCommunication_T = blockforest::SimpleCommunication< LatticeModelStencil_T >;
+
+// the geometry computations in SurfaceGeometryHandler require meaningful values in the ghost layers in corner
+// directions (flag field and fill level field); this holds, even if the lattice model uses a D3Q19 stencil
+using CommunicationStencil_T =
+ typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
+using Communication_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
+
+using flag_t = uint32_t;
+using FlagField_T = FlagField< flag_t >;
+using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
+
+// function describing the initialization velocity profile (in global cell coordinates)
+inline Vector3< real_t > velocityProfile(Cell globalCell, real_t timePeriod, uint_t timestep,
+ const Vector3< real_t >& domainSize)
+{
+ // add 0.5 to get the cell's center
+ const real_t x = (real_c(globalCell.x()) + real_c(0.5)) / domainSize[0];
+ const real_t y = (real_c(globalCell.y()) + real_c(0.5)) / domainSize[1];
+
+ const real_t xToDomainCenter = x - real_c(0.5);
+ const real_t yToDomainCenter = y - real_c(0.5);
+ const real_t r = real_c(std::sqrt(xToDomainCenter * xToDomainCenter + yToDomainCenter * yToDomainCenter));
+ const real_t rTerm = (real_c(1) - real_c(2) * r) * (real_c(1) - real_c(2) * r);
+
+ const real_t timeTerm = real_c(std::cos(math::pi * real_t(timestep) / timePeriod));
+
+ const real_t velocityX = real_c(std::sin(real_c(2) * math::pi * y)) * real_c(std::sin(math::pi * x)) *
+ real_c(std::sin(math::pi * x)) * timeTerm;
+ const real_t velocityY = -real_c(std::sin(real_c(2) * math::pi * x)) * real_c(std::sin(math::pi * y)) *
+ real_c(std::sin(math::pi * y)) * timeTerm;
+ const real_t velocityZ = rTerm * timeTerm;
+
+ return Vector3< real_t >(velocityX, velocityY, velocityZ);
+}
+
+int main(int argc, char** argv)
+{
+ Environment walberlaEnv(argc, argv);
+
+ if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
+
+ // print content of parameter file
+ WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
+
+ // get block forest parameters from parameter file
+ auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
+ const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+ const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
+
+ // get domain parameters from parameter file
+ auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
+ const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
+
+ const real_t bubbleDiameter = real_c(domainWidth) * real_c(0.075);
+ const Vector3< real_t > bubbleCenter = domainWidth * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25));
+
+ // define domain size
+ Vector3< uint_t > domainSize;
+ domainSize[0] = domainWidth;
+ domainSize[1] = domainWidth;
+ domainSize[2] = domainWidth;
+
+ // compute number of blocks as defined by domainSize and cellsPerBlock
+ Vector3< uint_t > numBlocks;
+ numBlocks[0] = uint_c(std::ceil(real_c(domainSize[0]) / real_c(cellsPerBlock[0])));
+ numBlocks[1] = uint_c(std::ceil(real_c(domainSize[1]) / real_c(cellsPerBlock[1])));
+ numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
+
+ // get number of (MPI) processes
+ const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
+ WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
+ "The number of MPI processes is greater than the number of blocks as defined by "
+ "\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
+ "the number of MPI processes or increase \"cellsPerBlock\".")
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellsPerBlock);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainSize);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainWidth);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleDiameter);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(bubbleCenter);
+
+ // get physics parameters from parameter file
+ auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
+ const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+ const uint_t timestepsToInitialPosition = physicsParameters.getParameter< uint_t >("timestepsToInitialPosition");
+
+ // compute CFL number
+ const real_t dx_SI = real_c(1) / real_c(domainWidth);
+ const real_t dt_SI = real_c(3) / real_c(timestepsToInitialPosition);
+ const real_t CFL = dt_SI / dx_SI; // with velocity_SI = 1
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(CFL);
+
+ // dummy collision model (LBM not simulated in this benchmark)
+ const CollisionModel_T collisionModel = CollisionModel_T(real_c(1));
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timestepsToInitialPosition);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
+
+ // read model parameters from parameter file
+ const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
+ const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
+ const std::string excessMassDistributionModel =
+ modelParameters.getParameter< std::string >("excessMassDistributionModel");
+ const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
+ const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
+ const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
+ const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
+
+ // read evaluation parameters from parameter file
+ const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
+ const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
+
+ // create non-uniform block forest (non-uniformity required for load balancing)
+ const std::shared_ptr< StructuredBlockForest > blockForest =
+ createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
+
+ // create lattice model
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel);
+
+ // add pdf field
+ const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
+
+ // add fill level field (initialized with 1, i.e., liquid everywhere)
+ const BlockDataID fillFieldID =
+ field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
+
+ // add boundary handling
+ const std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling =
+ std::make_shared< FreeSurfaceBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID);
+ const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
+ const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
+
+ // initialize the velocity profile
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+ FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
+ // cell in block-local coordinates
+ const Cell localCell = pdfFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = pdfFieldIt.cell();
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
+ const Vector3< real_t > initialVelocity =
+ CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), uint_c(0), domainSize);
+ pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ // create the spherical bubble
+ const geometry::Sphere sphereBubble(real_c(domainWidth) * Vector3< real_t >(real_c(0.5), real_c(0.75), real_c(0.25)),
+ real_c(domainWidth) * real_c(0.15));
+ bubble_model::addBodyToFillLevelField< geometry::Sphere >(*blockForest, fillFieldID, sphereBubble, true);
+
+ // initialize domain boundary conditions from config file
+ const auto boundaryParameters = walberlaEnv.config()->getOneBlock("BoundaryParameters");
+ freeSurfaceBoundaryHandling->initFromConfig(boundaryParameters);
+
+ // IMPORTANT REMARK: this must be called only after every solid flag has been set; otherwise, the boundary handling
+ // might not detect solid flags correctly
+ freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
+
+ // communication after initialization
+ Communication_T communication(blockForest, flagFieldID, fillFieldID);
+ communication();
+
+ PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
+ pdfCommunication();
+
+ const ConstBlockDataID initialFillFieldID =
+ field::addCloneToStorage< ScalarField_T >(blockForest, fillFieldID, "Initial fill level field");
+
+ // add bubble model
+ const std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel =
+ std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1));
+
+ // create timeloop
+ SweepTimeloop timeloop(blockForest, timesteps);
+
+ // add surface geometry handler
+ const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
+ blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, false, false, real_c(0));
+
+ geometryHandler.addSweeps(timeloop);
+
+ // get fields created by surface geometry handler
+ const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
+
+ // add boundary handling for standard boundaries and free surface boundaries
+ const AdvectionDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, normalFieldID, freeSurfaceBoundaryHandling, bubbleModel,
+ pdfReconstructionModel, excessMassDistributionModel, useSimpleMassExchange, cellConversionThreshold,
+ cellConversionForceThreshold);
+
+ dynamicsHandler.addSweeps(timeloop);
+
+ // add evaluator for geometrical
+ const std::shared_ptr< real_t > geometricalError = std::make_shared< real_t >(real_c(0));
+ const GeometricalErrorEvaluator< FreeSurfaceBoundaryHandling_T, FlagField_T, ScalarField_T >
+ geometricalErrorEvaluator(blockForest, freeSurfaceBoundaryHandling, initialFillFieldID, fillFieldID,
+ evaluationFrequency, geometricalError);
+ timeloop.addFuncAfterTimeStep(geometricalErrorEvaluator, "Evaluator: geometrical error");
+
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
+ // add VTK output
+ addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(),
+ geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
+ geometryHandler.getObstNormalFieldID());
+
+ // add triangle mesh output of free surface
+ SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
+ blockForest, fillFieldID, flagFieldID, flagIDs::liquidInterfaceGasFlagIDs, real_c(0), walberlaEnv.config());
+ surfaceMeshWriter(); // write initial mesh
+ timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
+
+ // add logging for computational performance
+ const lbm::PerformanceLogger< FlagField_T > performanceLogger(
+ blockForest, flagFieldID, flagIDs::liquidInterfaceFlagIDs, performanceLogFrequency);
+ timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluator: performance logging");
+
+ WcTimingPool timingPool;
+
+ for (uint_t t = uint_c(0); t != timesteps; ++t)
+ {
+ timeloop.singleStep(timingPool, true);
+
+ if (t % evaluationFrequency == uint_c(0))
+ {
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass << "\n\t\texcess mass = "
+ << *excessMass << "\n\t\tgeometrical error = " << *geometricalError);
+ }
+
+ // set the constant velocity profile
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+ FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
+ const Cell localCell = pdfFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = pdfFieldIt.cell();
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // set velocity profile (CFL_SI = CFL_LBM = velocity_LBM * dt_LBM / dx_LBM = velocity_LBM)
+ const Vector3< real_t > velocity =
+ CFL * velocityProfile(globalCell, real_c(timestepsToInitialPosition), t, domainSize);
+ pdfField->setDensityAndVelocity(localCell, velocity, real_c(1));
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ pdfCommunication();
+
+ if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace SingleVortex
+} // namespace free_surface
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::free_surface::SingleVortex::main(argc, argv); }
\ No newline at end of file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.prm b/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.prm
new file mode 100644
index 0000000000000000000000000000000000000000..359fd4a98093e95ace99e755819f022cca40e5c9
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/SingleVortex.prm
@@ -0,0 +1,96 @@
+BlockForestParameters
+{
+ cellsPerBlock < 25, 25, 25 >;
+ periodicity < 1, 1, 1 >;
+}
+
+DomainParameters
+{
+ domainWidth 50;
+}
+
+PhysicsParameters
+{
+ timestepsToInitialPosition 3000;
+ timesteps 3001;
+
+}
+
+ModelParameters
+{
+ pdfReconstructionModel OnlyMissing;
+ excessMassDistributionModel EvenlyAllInterface;
+ curvatureModel FiniteDifferenceMethod;
+ useSimpleMassExchange false;
+ cellConversionThreshold 1e-2;
+ cellConversionForceThreshold 1e-1;
+}
+
+EvaluationParameters
+{
+ evaluationFrequency 300;
+ performanceLogFrequency 10000;
+}
+
+BoundaryParameters
+{
+ // X
+ //Border { direction W; walldistance -1; FreeSlip{} }
+ //Border { direction E; walldistance -1; FreeSlip{} }
+
+ // Y
+ //Border { direction N; walldistance -1; FreeSlip{} }
+ //Border { direction S; walldistance -1; FreeSlip{} }
+
+ // Z
+ //Border { direction T; walldistance -1; FreeSlip{} }
+ //Border { direction B; walldistance -1; FreeSlip{} }
+}
+
+MeshOutputParameters
+{
+ writeFrequency 300;
+ baseFolder mesh-out;
+}
+
+VTK
+{
+ fluid_field
+ {
+ writeFrequency 300;
+ ghostLayers 0;
+ baseFolder vtk-out;
+ samplingResolution 1;
+
+ writers
+ {
+ fill_level;
+ mapped_flag;
+ velocity;
+ density;
+ //curvature;
+ //normal;
+ //obstacle_normal;
+ //pdf;
+ //flag;
+ }
+
+ inclusion_filters
+ {
+ //liquidInterfaceFilter; // only include liquid and interface cells in VTK output
+ }
+
+ before_functions
+ {
+ //ghost_layer_synchronization; // only needed if writing the ghost layer
+ gas_cell_zero_setter; // sets velocity=0 and density=1 all gas cells before writing VTK output
+ }
+ }
+
+ domain_decomposition
+ {
+ writeFrequency 0;
+ baseFolder vtk-out;
+ outputDomainDecomposition true;
+ }
+}
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/ZalesakDisk.cpp b/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp
similarity index 65%
rename from apps/showcases/FreeSurface/ZalesakDisk.cpp
rename to apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp
index edbaab294c474804a66bfc56a0066b7e73ce301e..f2177d55950a4a2a57e113b7dcc12344aa4b459b 100644
--- a/apps/showcases/FreeSurface/ZalesakDisk.cpp
+++ b/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.cpp
@@ -16,12 +16,14 @@
//! \file ZalesakDisk.cpp
//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
//
-// This showcase simulates a slotted disk of gas in a constant rotating velocity field. This benchmark is commonly
-// referred to as Zalesak's rotating disk (see doi: 10.1016/0021-9991(79)90051-2).
+// This benchmark simulates the advection of a slotted disk of gas in a constant rotating velocity field. The disk
+// returns to its initial position, where it should take its initial form. The relative geometrical error of the
+// bubble's shape after one rotation is evaluated. There is no LBM flow simulation performed here, it is a test case for
+// the FSLBM's mass advection. This benchmark is commonly referred to as Zalesak's rotating disk (see
+// doi: 10.1016/0021-9991(79)90051-2). The setup chosen here is identical to the one used by Janssen (see
+// doi: 10.1016/j.camwa.2009.08.064).
//======================================================================================================================
-#include "blockforest/Initialization.h"
-
#include "core/Environment.h"
#include "field/Gather.h"
@@ -33,11 +35,13 @@
#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
-#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
#include "lbm/free_surface/surface_geometry/SurfaceGeometryHandler.h"
#include "lbm/free_surface/surface_geometry/Utility.h"
#include "lbm/lattice_model/D2Q9.h"
+#include "functionality/AdvectionDynamicsHandler.h"
+#include "functionality/GeometricalErrorEvaluator.h"
+
namespace walberla
{
namespace free_surface
@@ -47,9 +51,9 @@ namespace ZalesakDisk
using ScalarField_T = GhostLayerField< real_t, 1 >;
using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+// Lattice model is only created for dummy purposes; no LBM simulation is performed
using CollisionModel_T = lbm::collision_model::SRT;
-using ForceModel_T = lbm::force_model::GuoField< VectorField_T >;
-using LatticeModel_T = lbm::D2Q9< CollisionModel_T, true, ForceModel_T, 2 >;
+using LatticeModel_T = lbm::D2Q9< CollisionModel_T, true >;
using LatticeModelStencil_T = LatticeModel_T::Stencil;
using PdfField_T = lbm::PdfField< LatticeModel_T >;
using PdfCommunication_T = blockforest::SimpleCommunication< LatticeModelStencil_T >;
@@ -65,7 +69,7 @@ using FlagField_T = FlagField< flag_t >;
using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
// function describing the initialization velocity profile (in global cell coordinates)
-inline Vector3< real_t > velocityProfile(real_t angularVelocity, Cell globalCell, Vector3< real_t > domainCenter)
+inline Vector3< real_t > velocityProfile(real_t angularVelocity, Cell globalCell, const Vector3< real_t >& domainCenter)
{
// add 0.5 to get Cell's center
const real_t velocityX = -angularVelocity * ((real_c(globalCell.y()) + real_c(0.5)) - domainCenter[0]);
@@ -84,21 +88,19 @@ int main(int argc, char** argv)
WALBERLA_LOG_INFO_ON_ROOT(*walberlaEnv.config());
// get block forest parameters from parameter file
- auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
- const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
- const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
- const uint_t loadBalancingFrequency = blockForestParameters.getParameter< uint_t >("loadBalancingFrequency");
- const bool printLoadBalancingStatistics = blockForestParameters.getParameter< bool >("printLoadBalancingStatistics");
+ auto blockForestParameters = walberlaEnv.config()->getOneBlock("BlockForestParameters");
+ const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+ const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
// get domain parameters from parameter file
auto domainParameters = walberlaEnv.config()->getOneBlock("DomainParameters");
const uint_t domainWidth = domainParameters.getParameter< uint_t >("domainWidth");
- const real_t diskRadius = real_c(domainWidth) * real_c(0.15);
+ const real_t diskRadius = real_c(domainWidth) * real_c(0.125);
const Vector3< real_t > diskCenter =
Vector3< real_t >(real_c(domainWidth) * real_c(0.5), real_c(domainWidth) * real_c(0.75), real_c(0.5));
- const real_t diskSlotLength = real_c(0.25) * real_c(domainWidth);
- const real_t diskSlotWidth = real_c(0.05) * real_c(domainWidth);
+ const real_t diskSlotLength = real_c(2) * diskRadius - real_c(0.1) * real_c(domainWidth);
+ const real_t diskSlotWidth = real_c(0.06) * real_c(domainWidth);
// define domain size
Vector3< uint_t > domainSize;
@@ -115,7 +117,7 @@ int main(int argc, char** argv)
numBlocks[2] = uint_c(std::ceil(real_c(domainSize[2]) / real_c(cellsPerBlock[2])));
// get number of (MPI) processes
- uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
+ const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
WALBERLA_CHECK_LESS_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
"The number of MPI processes is greater than the number of blocks as defined by "
"\"domainSize/cellsPerBlock\". This would result in unused MPI processes. Either decrease "
@@ -127,79 +129,68 @@ int main(int argc, char** argv)
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(domainWidth);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(periodicity);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(loadBalancingFrequency);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(printLoadBalancingStatistics);
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskRadius);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskCenter);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskSlotLength);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(diskSlotWidth);
// get physics parameters from parameter file
- auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
- const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+ auto physicsParameters = walberlaEnv.config()->getOneBlock("PhysicsParameters");
+ const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+ const uint_t timestepsFullRotation = physicsParameters.getParameter< uint_t >("timestepsFullRotation");
- const real_t relaxationRate = physicsParameters.getParameter< real_t >("relaxationRate");
- const CollisionModel_T collisionModel = CollisionModel_T(relaxationRate);
- const real_t viscosity = collisionModel.viscosity();
+ // compute CFL number
+ const real_t dx_SI = real_c(4) / real_c(domainWidth);
+ const real_t dt_SI = real_c(12.59652) / real_c(timestepsFullRotation);
+ const real_t CFL = dt_SI / dx_SI; // with velocity_SI = 1
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(CFL);
- const real_t reynoldsNumber = physicsParameters.getParameter< real_t >("reynoldsNumber");
- const real_t angularVelocity =
- reynoldsNumber * viscosity / (real_c(0.5) * real_c(domainWidth) * real_c(domainWidth));
- const Vector3< real_t > force(real_c(0), real_c(0), real_c(0));
+ // dummy collision model (LBM not simulated in this benchmark)
+ const CollisionModel_T collisionModel = CollisionModel_T(real_c(1));
- const bool enableWetting = physicsParameters.getParameter< bool >("enableWetting");
- const real_t contactAngle = physicsParameters.getParameter< real_t >("contactAngle");
+ const real_t angularVelocity = real_c(2) * math::pi / real_c(timestepsFullRotation);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(reynoldsNumber);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(relaxationRate);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableWetting);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(contactAngle);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timestepsFullRotation);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(timesteps);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(angularVelocity);
- WALBERLA_LOG_DEVEL_ON_ROOT("Timesteps for full rotation " << real_c(2) * math::pi / angularVelocity);
// read model parameters from parameter file
const auto modelParameters = walberlaEnv.config()->getOneBlock("ModelParameters");
const std::string pdfReconstructionModel = modelParameters.getParameter< std::string >("pdfReconstructionModel");
- const std::string pdfRefillingModel = modelParameters.getParameter< std::string >("pdfRefillingModel");
const std::string excessMassDistributionModel =
modelParameters.getParameter< std::string >("excessMassDistributionModel");
const std::string curvatureModel = modelParameters.getParameter< std::string >("curvatureModel");
const bool useSimpleMassExchange = modelParameters.getParameter< bool >("useSimpleMassExchange");
const real_t cellConversionThreshold = modelParameters.getParameter< real_t >("cellConversionThreshold");
const real_t cellConversionForceThreshold = modelParameters.getParameter< real_t >("cellConversionForceThreshold");
- const bool enableBubbleModel = modelParameters.getParameter< bool >("enableBubbleModel");
- const bool enableBubbleSplits = modelParameters.getParameter< bool >("enableBubbleSplits");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfReconstructionModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(pdfRefillingModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(excessMassDistributionModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(curvatureModel);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(useSimpleMassExchange);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionThreshold);
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(cellConversionForceThreshold);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleModel);
- WALBERLA_LOG_DEVEL_VAR_ON_ROOT(enableBubbleSplits);
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
createNonUniformBlockForest(domainSize, cellsPerBlock, numBlocks, periodicity);
- // add force field
- const BlockDataID forceDensityFieldID =
- field::addToStorage< VectorField_T >(blockForest, "Force field", force, field::fzyx, uint_c(1));
-
// create lattice model
- const LatticeModel_T latticeModel = LatticeModel_T(collisionModel, ForceModel_T(forceDensityFieldID));
+ const LatticeModel_T latticeModel = LatticeModel_T(collisionModel);
// add pdf field
const BlockDataID pdfFieldID = lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel, field::fzyx);
- // add fill level field (initialized with 0, i.e., gas everywhere)
+ // add fill level field (initialized with 1, i.e., liquid everywhere)
const BlockDataID fillFieldID =
field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(1.0), field::fzyx, uint_c(2));
@@ -212,10 +203,9 @@ int main(int argc, char** argv)
// initialize the velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
- PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
- FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
+ PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
- WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
+ WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, {
// cell in block-local coordinates
const Cell localCell = pdfFieldIt.cell();
@@ -252,66 +242,57 @@ int main(int argc, char** argv)
freeSurfaceBoundaryHandling->initFlagsFromFillLevel();
// communication after initialization
- Communication_T communication(blockForest, flagFieldID, fillFieldID, forceDensityFieldID);
+ Communication_T communication(blockForest, flagFieldID, fillFieldID);
communication();
PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
pdfCommunication();
- // add bubble model
- std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel = nullptr;
- if (enableBubbleModel)
- {
- const std::shared_ptr< bubble_model::BubbleModel< LatticeModelStencil_T > > bubbleModelDerived =
- std::make_shared< bubble_model::BubbleModel< LatticeModelStencil_T > >(blockForest, enableBubbleSplits);
- bubbleModelDerived->initFromFillLevelField(fillFieldID);
-
- bubbleModel = std::static_pointer_cast< bubble_model::BubbleModelBase >(bubbleModelDerived);
- }
- else { bubbleModel = std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1)); }
+ const ConstBlockDataID initialFillFieldID =
+ field::addCloneToStorage< ScalarField_T >(blockForest, fillFieldID, "Initial fill level field");
- // set density in non-liquid or non-interface cells to one (after initializing with hydrostatic pressure)
- setDensityInNonFluidCellsToOne< FlagField_T, PdfField_T >(blockForest, flagInfo, flagFieldID, pdfFieldID);
+ // add bubble model
+ const std::shared_ptr< bubble_model::BubbleModelBase > bubbleModel =
+ std::make_shared< bubble_model::BubbleModelConstantPressure >(real_c(1));
// create timeloop
SweepTimeloop timeloop(blockForest, timesteps);
- const real_t surfaceTension = real_c(0);
-
- // Laplace pressure = 2 * surface tension * curvature; curvature computation is not necessary with zero surface
- // tension
- bool computeCurvature = false;
- if (!realIsEqual(surfaceTension, real_c(0), real_c(1e-14))) { computeCurvature = true; }
-
// add surface geometry handler
const SurfaceGeometryHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > geometryHandler(
- blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, computeCurvature, enableWetting,
- contactAngle);
+ blockForest, freeSurfaceBoundaryHandling, fillFieldID, curvatureModel, false, false, real_c(0));
geometryHandler.addSweeps(timeloop);
// get fields created by surface geometry handler
- const ConstBlockDataID curvatureFieldID = geometryHandler.getConstCurvatureFieldID();
- const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
+ const ConstBlockDataID normalFieldID = geometryHandler.getConstNormalFieldID();
// add boundary handling for standard boundaries and free surface boundaries
- const SurfaceDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
- blockForest, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID, normalFieldID, curvatureFieldID,
- freeSurfaceBoundaryHandling, bubbleModel, pdfReconstructionModel, pdfRefillingModel, excessMassDistributionModel,
- relaxationRate, force, surfaceTension, useSimpleMassExchange, cellConversionThreshold,
+ const AdvectionDynamicsHandler< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T > dynamicsHandler(
+ blockForest, pdfFieldID, flagFieldID, fillFieldID, normalFieldID, freeSurfaceBoundaryHandling, bubbleModel,
+ pdfReconstructionModel, excessMassDistributionModel, useSimpleMassExchange, cellConversionThreshold,
cellConversionForceThreshold);
dynamicsHandler.addSweeps(timeloop);
- // add load balancing
- const LoadBalancer< FlagField_T, CommunicationStencil_T, LatticeModelStencil_T > loadBalancer(
- blockForest, communication, pdfCommunication, bubbleModel, uint_c(50), uint_c(10), uint_c(5),
- loadBalancingFrequency, printLoadBalancingStatistics);
- timeloop.addFuncAfterTimeStep(loadBalancer, "Sweep: load balancing");
+ // add evaluator for geometrical
+ const std::shared_ptr< real_t > geometricalError = std::make_shared< real_t >(real_c(0));
+ const GeometricalErrorEvaluator< FreeSurfaceBoundaryHandling_T, FlagField_T, ScalarField_T >
+ geometricalErrorEvaluator(blockForest, freeSurfaceBoundaryHandling, initialFillFieldID, fillFieldID,
+ evaluationFrequency, geometricalError);
+ timeloop.addFuncAfterTimeStep(geometricalErrorEvaluator, "Evaluator: geometrical errors");
+
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
- blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
+ blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(),
geometryHandler.getCurvatureFieldID(), geometryHandler.getNormalFieldID(),
geometryHandler.getObstNormalFieldID());
@@ -332,6 +313,12 @@ int main(int argc, char** argv)
{
timeloop.singleStep(timingPool, true);
+ if (t % evaluationFrequency == uint_c(0))
+ {
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass << "\n\t\texcess mass = "
+ << *excessMass << "\n\t\tgeometrical error = " << *geometricalError);
+ }
+
// set the constant velocity profile
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
@@ -339,22 +326,20 @@ int main(int argc, char** argv)
FlagField_T* const flagField = blockIt->getData< FlagField_T >(flagFieldID);
WALBERLA_FOR_ALL_CELLS(pdfFieldIt, pdfField, flagFieldIt, flagField, {
- if (flagInfo.isInterface(flagFieldIt) || flagInfo.isLiquid(flagFieldIt))
- {
- // cell in block-local coordinates
- const Cell localCell = pdfFieldIt.cell();
-
- // get cell in global coordinates
- Cell globalCell = pdfFieldIt.cell();
- blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
-
- // set velocity profile
- const Vector3< real_t > initialVelocity = velocityProfile(angularVelocity, globalCell, domainCenter);
- pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
- }
+ const Cell localCell = pdfFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = pdfFieldIt.cell();
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // set velocity profile
+ const Vector3< real_t > initialVelocity = velocityProfile(angularVelocity, globalCell, domainCenter);
+ pdfField->setDensityAndVelocity(localCell, initialVelocity, real_c(1));
}) // WALBERLA_FOR_ALL_CELLS
}
+ pdfCommunication();
+
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
diff --git a/apps/showcases/FreeSurface/ZalesakDisk.prm b/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.prm
similarity index 63%
rename from apps/showcases/FreeSurface/ZalesakDisk.prm
rename to apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.prm
index b9af1280c96a12bba2f0ae5725675ef6f34e7429..857f1510bcd1d0a9166b43b8e1dc24d227dfb181 100644
--- a/apps/showcases/FreeSurface/ZalesakDisk.prm
+++ b/apps/benchmarks/FreeSurfaceAdvection/ZalesakDisk.prm
@@ -1,53 +1,46 @@
BlockForestParameters
{
- cellsPerBlock < 50, 50, 1 >;
- periodicity < 0, 0, 1 >;
- loadBalancingFrequency 0;
- printLoadBalancingStatistics false;
+ cellsPerBlock < 100, 100, 1 >;
+ periodicity < 0, 0, 0 >;
}
DomainParameters
{
- domainWidth 100;
+ domainWidth 200;
}
PhysicsParameters
{
- reynoldsNumber 100; // Re = angularVelocity * domainWidth * 0.5 * domainWidth / kin. viscosity
- relaxationRate 1.8;
- enableWetting false;
- contactAngle 0; // only used if enableWetting=true
- timesteps 1000000;
+ timestepsFullRotation 12570; // angularVelocity = 2 * pi / timestepsFullRotation
+ timesteps 12571;
+
}
ModelParameters
{
pdfReconstructionModel OnlyMissing;
- pdfRefillingModel EquilibriumRefilling;
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
useSimpleMassExchange false;
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
-
- enableBubbleModel True;
- enableBubbleSplits false; // only used if enableBubbleModel=true
}
EvaluationParameters
{
+ evaluationFrequency 12570;
performanceLogFrequency 10000;
}
BoundaryParameters
{
// X
- Border { direction W; walldistance -1; FreeSlip{} }
- Border { direction E; walldistance -1; FreeSlip{} }
+ //Border { direction W; walldistance -1; FreeSlip{} }
+ //Border { direction E; walldistance -1; FreeSlip{} }
// Y
- Border { direction N; walldistance -1; FreeSlip{} }
- Border { direction S; walldistance -1; FreeSlip{} }
+ //Border { direction N; walldistance -1; FreeSlip{} }
+ //Border { direction S; walldistance -1; FreeSlip{} }
// Z
//Border { direction T; walldistance -1; FreeSlip{} }
@@ -64,7 +57,7 @@ VTK
{
fluid_field
{
- writeFrequency 4241;
+ writeFrequency 12570;
ghostLayers 0;
baseFolder vtk-out;
samplingResolution 1;
@@ -80,7 +73,6 @@ VTK
//obstacle_normal;
//pdf;
//flag;
- //force_density;
}
inclusion_filters
diff --git a/apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectSweep.h b/apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectSweep.h
new file mode 100644
index 0000000000000000000000000000000000000000..80db158e6b6f341e3f8d637efa511e0a709dfa2b
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectSweep.h
@@ -0,0 +1,152 @@
+//======================================================================================================================
+//
+// 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 AdvectSweep.h
+//! \ingroup free_surface
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Sweep for mass advection and interface cell conversion marking (simplified StreamReconstructAdvectSweep).
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+#include "core/math/Vector3.h"
+#include "core/mpi/Reduce.h"
+#include "core/timing/TimingPool.h"
+
+#include "field/FieldClone.h"
+#include "field/FlagField.h"
+
+#include "lbm/free_surface/FlagInfo.h"
+#include "lbm/free_surface/bubble_model/BubbleModel.h"
+#include "lbm/free_surface/dynamics/PdfReconstructionModel.h"
+#include "lbm/free_surface/dynamics/functionality/AdvectMass.h"
+#include "lbm/free_surface/dynamics/functionality/FindInterfaceCellConversion.h"
+#include "lbm/free_surface/dynamics/functionality/GetOredNeighborhood.h"
+#include "lbm/free_surface/dynamics/functionality/ReconstructInterfaceCellABB.h"
+#include "lbm/sweeps/StreamPull.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename FlagField_T, typename FlagInfo_T,
+ typename ScalarField_T, typename VectorField_T >
+class AdvectSweep
+{
+ public:
+ using flag_t = typename FlagInfo_T::flag_t;
+ using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+ AdvectSweep(BlockDataID handlingID, BlockDataID fillFieldID, BlockDataID flagFieldID, BlockDataID pdfField,
+ const FlagInfo_T& flagInfo, const PdfReconstructionModel& pdfReconstructionModel,
+ bool useSimpleMassExchange, real_t cellConversionThreshold, real_t cellConversionForceThreshold)
+ : handlingID_(handlingID), fillFieldID_(fillFieldID), flagFieldID_(flagFieldID), pdfFieldID_(pdfField),
+ flagInfo_(flagInfo), neighborhoodFlagFieldClone_(flagFieldID), fillFieldClone_(fillFieldID),
+ pdfFieldClone_(pdfField), pdfReconstructionModel_(pdfReconstructionModel),
+ useSimpleMassExchange_(useSimpleMassExchange), cellConversionThreshold_(cellConversionThreshold),
+ cellConversionForceThreshold_(cellConversionForceThreshold)
+ {}
+
+ void operator()(IBlock* const block);
+
+ protected:
+ BlockDataID handlingID_;
+ BlockDataID fillFieldID_;
+ BlockDataID flagFieldID_;
+ BlockDataID pdfFieldID_;
+
+ FlagInfo_T flagInfo_;
+
+ // efficient clones of fields to provide temporary fields (for writing to)
+ field::FieldClone< FlagField_T, true > neighborhoodFlagFieldClone_;
+ field::FieldClone< ScalarField_T, true > fillFieldClone_;
+ field::FieldClone< PdfField_T, true > pdfFieldClone_;
+
+ PdfReconstructionModel pdfReconstructionModel_;
+ bool useSimpleMassExchange_;
+ real_t cellConversionThreshold_;
+ real_t cellConversionForceThreshold_;
+}; // class AdvectSweep
+
+template< typename LatticeModel_T, typename BoundaryHandling_T, typename FlagField_T, typename FlagInfo_T,
+ typename ScalarField_T, typename VectorField_T >
+void AdvectSweep< LatticeModel_T, BoundaryHandling_T, FlagField_T, FlagInfo_T, ScalarField_T,
+ VectorField_T >::operator()(IBlock* const block)
+{
+ // fetch data
+ ScalarField_T* const fillSrcField = block->getData< ScalarField_T >(fillFieldID_);
+ PdfField_T* const pdfSrcField = block->getData< PdfField_T >(pdfFieldID_);
+
+ FlagField_T* const flagField = block->getData< FlagField_T >(flagFieldID_);
+
+ // temporary fields that act as destination fields
+ FlagField_T* const neighborhoodFlagField = neighborhoodFlagFieldClone_.get(block);
+ ScalarField_T* const fillDstField = fillFieldClone_.get(block);
+
+ // combine all neighbor flags using bitwise OR and write them to the neighborhood field
+ // this is simply a pre-computation of often required values
+ // IMPORTANT REMARK: the "OredNeighborhood" is also required for each cell in the first ghost layer; this requires
+ // access to all first ghost layer cell's neighbors (i.e., to the second ghost layer)
+ WALBERLA_CHECK_GREATER_EQUAL(flagField->nrOfGhostLayers(), uint_c(2));
+ getOredNeighborhood< typename LatticeModel_T::Stencil >(flagField, neighborhoodFlagField);
+
+ // explicitly avoid OpenMP due to bubble model update (reportFillLevelChange)
+ WALBERLA_FOR_ALL_CELLS_OMP(
+ pdfSrcFieldIt, pdfSrcField, fillSrcFieldIt, fillSrcField, fillDstFieldIt, fillDstField, flagFieldIt, flagField,
+ neighborhoodFlagFieldIt, neighborhoodFlagField, omp critical, {
+ if (flagInfo_.isInterface(flagFieldIt))
+ {
+ const real_t rho = lbm::getDensity(pdfSrcField->latticeModel(), pdfSrcFieldIt);
+
+ // compute mass advection using post-collision PDFs (explicitly not PDFs updated by stream above)
+ const real_t deltaMass =
+ (advectMass< LatticeModel_T, FlagField_T, typename ScalarField_T::iterator,
+ typename PdfField_T::iterator, typename FlagField_T::iterator,
+ typename FlagField_T::iterator, FlagInfo_T >) (flagField, fillSrcFieldIt, pdfSrcFieldIt,
+ flagFieldIt, neighborhoodFlagFieldIt,
+ flagInfo_, useSimpleMassExchange_);
+
+ // update fill level after LBM stream and mass exchange
+ *fillDstFieldIt = *fillSrcFieldIt + deltaMass / rho;
+ }
+ else // treat non-interface cells
+ {
+ // manually adjust the fill level to avoid outdated fill levels being copied from fillSrcField
+ if (flagInfo_.isGas(flagFieldIt)) { *fillDstFieldIt = real_c(0); }
+ else
+ {
+ if (flagInfo_.isLiquid(flagFieldIt)) { *fillDstFieldIt = real_c(1); }
+ else // flag is e.g. obstacle or outflow
+ {
+ *fillDstFieldIt = *fillSrcFieldIt;
+ }
+ }
+ }
+ }) // WALBERLA_FOR_ALL_CELLS_XYZ_OMP
+
+ fillSrcField->swapDataPointers(fillDstField);
+
+ BoundaryHandling_T* const handling = block->getData< BoundaryHandling_T >(handlingID_);
+
+ // mark interface cell for conversion
+ findInterfaceCellConversions< LatticeModel_T >(handling, fillSrcField, flagField, neighborhoodFlagField, flagInfo_,
+ cellConversionThreshold_, cellConversionForceThreshold_);
+}
+
+} // namespace free_surface
+} // namespace walberla
diff --git a/apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectionDynamicsHandler.h b/apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectionDynamicsHandler.h
new file mode 100644
index 0000000000000000000000000000000000000000..1e48a1683a0bde1ed7392cf3b0e295c5189abf72
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectionDynamicsHandler.h
@@ -0,0 +1,254 @@
+//======================================================================================================================
+//
+// 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 AdvectionDynamicsHandler.h
+//! \ingroup surface_dynamics
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Handles free surface advection (without LBM flow simulation; this is a simplified SurfaceDynamicsHandler).
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "core/DataTypes.h"
+
+#include "domain_decomposition/StructuredBlockStorage.h"
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+
+#include "lbm/blockforest/communication/SimpleCommunication.h"
+#include "lbm/blockforest/communication/UpdateSecondGhostLayer.h"
+#include "lbm/free_surface/BlockStateDetectorSweep.h"
+#include "lbm/free_surface/FlagInfo.h"
+#include "lbm/free_surface/boundary/FreeSurfaceBoundaryHandling.h"
+#include "lbm/free_surface/bubble_model/BubbleModel.h"
+#include "lbm/free_surface/dynamics/CellConversionSweep.h"
+#include "lbm/free_surface/dynamics/ConversionFlagsResetSweep.h"
+#include "lbm/free_surface/dynamics/ExcessMassDistributionModel.h"
+#include "lbm/free_surface/dynamics/ExcessMassDistributionSweep.h"
+#include "lbm/free_surface/dynamics/PdfReconstructionModel.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/sweeps/SweepWrappers.h"
+
+#include "stencil/D3Q27.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "AdvectSweep.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename VectorField_T >
+class AdvectionDynamicsHandler
+{
+ protected:
+ using Communication_T = blockforest::SimpleCommunication< typename LatticeModel_T::Stencil >;
+
+ // communication in corner directions (D2Q9/D3Q27) is required for all fields but the PDF field
+ using CommunicationStencil_T =
+ typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
+ using CommunicationCorner_T = blockforest::SimpleCommunication< CommunicationStencil_T >;
+
+ using FreeSurfaceBoundaryHandling_T = FreeSurfaceBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T >;
+
+ public:
+ AdvectionDynamicsHandler(const std::shared_ptr< StructuredBlockForest >& blockForest, BlockDataID pdfFieldID,
+ BlockDataID flagFieldID, BlockDataID fillFieldID, ConstBlockDataID normalFieldID,
+ const std::shared_ptr< FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
+ const std::shared_ptr< BubbleModelBase >& bubbleModel,
+ const std::string& pdfReconstructionModel, const std::string& excessMassDistributionModel,
+ bool useSimpleMassExchange, real_t cellConversionThreshold,
+ real_t cellConversionForceThreshold)
+ : blockForest_(blockForest), pdfFieldID_(pdfFieldID), flagFieldID_(flagFieldID), fillFieldID_(fillFieldID),
+ normalFieldID_(normalFieldID), bubbleModel_(bubbleModel),
+ freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling), pdfReconstructionModel_(pdfReconstructionModel),
+ excessMassDistributionModel_({ excessMassDistributionModel }), useSimpleMassExchange_(useSimpleMassExchange),
+ cellConversionThreshold_(cellConversionThreshold), cellConversionForceThreshold_(cellConversionForceThreshold)
+ {
+ WALBERLA_CHECK(LatticeModel_T::compressible,
+ "The free surface lattice Boltzmann extension works only with compressible LBM models.");
+
+ if (excessMassDistributionModel_.isEvenlyAllInterfaceFallbackLiquidType())
+ {
+ // add additional field for storing excess mass in liquid cells
+ excessMassFieldID_ =
+ field::addToStorage< ScalarField_T >(blockForest_, "Excess mass", real_c(0), field::fzyx, uint_c(1));
+ }
+ }
+
+ ConstBlockDataID getConstExcessMassFieldID() const { return excessMassFieldID_; }
+
+ /********************************************************************************************************************
+ * The order of these sweeps is similar to page 40 in the dissertation of T. Pohl, 2008.
+ *******************************************************************************************************************/
+ void addSweeps(SweepTimeloop& timeloop) const
+ {
+ using StateSweep = BlockStateDetectorSweep< FlagField_T >;
+
+ const auto& flagInfo = freeSurfaceBoundaryHandling_->getFlagInfo();
+
+ const auto blockStateUpdate = StateSweep(blockForest_, flagInfo, flagFieldID_);
+
+ // empty sweeps required for using selectors (e.g. StateSweep::onlyGasAndBoundary)
+ const auto emptySweep = [](IBlock*) {};
+
+ // sweep for
+ // - reconstruction of PDFs in interface cells
+ // - streaming of PDFs in interface cells (and liquid cells on the same block)
+ // - advection of mass
+ // - update bubble volumes
+ // - marking interface cells for conversion
+ const AdvectSweep< LatticeModel_T, typename FreeSurfaceBoundaryHandling_T::BoundaryHandling_T, FlagField_T,
+ typename FreeSurfaceBoundaryHandling_T::FlagInfo_T, ScalarField_T, VectorField_T >
+ advectSweep(freeSurfaceBoundaryHandling_->getHandlingID(), fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
+ pdfReconstructionModel_, useSimpleMassExchange_, cellConversionThreshold_,
+ cellConversionForceThreshold_);
+ // sweep acts only on blocks with at least one interface cell (due to StateSweep::fullFreeSurface)
+ timeloop.add()
+ << Sweep(advectSweep, "Sweep: Advect", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep, "Empty sweep: Advect")
+ // do not communicate PDFs here:
+ // - stream on blocks with "StateSweep::fullFreeSurface" was performed here using post-collision PDFs
+ // - stream on other blocks is performed below and should also use post-collision PDFs
+ // => if PDFs were communicated here, the ghost layer of other blocks would have post-stream PDFs
+ << AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_, flagFieldID_),
+ "Communication: after Advect sweep")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
+ "Second ghost layer update: after Advect sweep (fill level field)")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest_, flagFieldID_),
+ "Second ghost layer update: after Advect sweep (flag field)");
+
+ // convert cells
+ // - according to the flags from StreamReconstructAdvectSweep (interface -> gas/liquid)
+ // - to ensure a closed layer of interface cells (gas/liquid -> interface)
+ // - detect and register bubble merges/splits (bubble volumes are already updated in StreamReconstructAdvectSweep)
+ // - convert cells and initialize PDFs near inflow boundaries
+ const CellConversionSweep< LatticeModel_T, typename FreeSurfaceBoundaryHandling_T::BoundaryHandling_T,
+ ScalarField_T >
+ cellConvSweep(freeSurfaceBoundaryHandling_->getHandlingID(), pdfFieldID_, flagInfo, bubbleModel_.get());
+ timeloop.add() << Sweep(cellConvSweep, "Sweep: cell conversion", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep, "Empty sweep: cell conversion")
+ << AfterFunction(Communication_T(blockForest_, pdfFieldID_),
+ "Communication: after cell conversion sweep (PDF field)")
+ // communicate the flag field also in corner directions
+ << AfterFunction(CommunicationCorner_T(blockForest_, flagFieldID_),
+ "Communication: after cell conversion sweep (flag field)")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest_, flagFieldID_),
+ "Second ghost layer update: after cell conversion sweep (flag field)");
+
+ // distribute excess mass:
+ // - excess mass: mass that is free after conversion from interface to gas/liquid cells
+ // - update the bubble model
+ if (excessMassDistributionModel_.isEvenlyType())
+ {
+ const ExcessMassDistributionSweepInterfaceEvenly< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T >
+ distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo);
+ timeloop.add() << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep, "Empty sweep: distribute excess mass")
+ << AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_),
+ "Communication: after excess mass distribution sweep")
+ << AfterFunction(
+ blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
+ "Second ghost layer update: after excess mass distribution sweep (fill level field)")
+ // update bubble model, i.e., perform registered bubble merges/splits; bubble merges/splits are
+ // already detected and registered by CellConversionSweep
+ << AfterFunction(std::bind(&bubble_model::BubbleModelBase::update, bubbleModel_),
+ "Sweep: bubble model update");
+ }
+ else
+ {
+ if (excessMassDistributionModel_.isWeightedType())
+ {
+ const ExcessMassDistributionSweepInterfaceWeighted< LatticeModel_T, FlagField_T, ScalarField_T,
+ VectorField_T >
+ distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
+ normalFieldID_);
+ timeloop.add() << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep, "Empty sweep: distribute excess mass")
+ << AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_),
+ "Communication: after excess mass distribution sweep")
+ << AfterFunction(
+ blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
+ "Second ghost layer update: after excess mass distribution sweep (fill level field)")
+ // update bubble model, i.e., perform registered bubble merges/splits; bubble merges/splits
+ // are already detected and registered by CellConversionSweep
+ << AfterFunction(std::bind(&bubble_model::BubbleModelBase::update, bubbleModel_),
+ "Sweep: bubble model update");
+ }
+ else
+ {
+ if (excessMassDistributionModel_.isEvenlyAllInterfaceFallbackLiquidType())
+ {
+ const ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T, ScalarField_T,
+ VectorField_T >
+ distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
+ excessMassFieldID_);
+ timeloop.add()
+ // perform this sweep also on "onlyLBM" blocks because liquid cells also exchange excess mass here
+ << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
+ << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::onlyLBM)
+ << Sweep(emptySweep, "Empty sweep: distribute excess mass")
+ << AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_, excessMassFieldID_),
+ "Communication: after excess mass distribution sweep")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest_, fillFieldID_),
+ "Second ghost layer update: after excess mass distribution sweep (fill level field)")
+ // update bubble model, i.e., perform registered bubble merges/splits; bubble
+ // merges/splits are already detected and registered by CellConversionSweep
+ << AfterFunction(std::bind(&bubble_model::BubbleModelBase::update, bubbleModel_),
+ "Sweep: bubble model update");
+ }
+ }
+ }
+
+ // reset all flags that signal cell conversions (except "keepInterfaceForWettingFlag")
+ ConversionFlagsResetSweep< FlagField_T > resetConversionFlagsSweep(flagFieldID_, flagInfo);
+ timeloop.add() << Sweep(resetConversionFlagsSweep, "Sweep: conversion flag reset", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep, "Empty sweep: conversion flag reset")
+ << AfterFunction(CommunicationCorner_T(blockForest_, flagFieldID_),
+ "Communication: after excess mass distribution sweep")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest_, flagFieldID_),
+ "Second ghost layer update: after excess mass distribution sweep (flag field)");
+
+ // update block states
+ timeloop.add() << Sweep(blockStateUpdate, "Sweep: block state update");
+ }
+
+ private:
+ std::shared_ptr< StructuredBlockForest > blockForest_;
+
+ BlockDataID pdfFieldID_;
+ BlockDataID flagFieldID_;
+ BlockDataID fillFieldID_;
+
+ ConstBlockDataID normalFieldID_;
+
+ std::shared_ptr< BubbleModelBase > bubbleModel_;
+ std::shared_ptr< FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
+
+ PdfReconstructionModel pdfReconstructionModel_;
+ ExcessMassDistributionModel excessMassDistributionModel_;
+
+ bool useSimpleMassExchange_;
+ real_t cellConversionThreshold_;
+ real_t cellConversionForceThreshold_;
+
+ BlockDataID excessMassFieldID_ = BlockDataID();
+}; // class AdvectionDynamicsHandler
+
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file
diff --git a/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h b/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h
new file mode 100644
index 0000000000000000000000000000000000000000..b0a5d7bbae749142f0a25c108a35b274ef146038
--- /dev/null
+++ b/apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h
@@ -0,0 +1,139 @@
+//======================================================================================================================
+//
+// 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 GeometricalErrorEvaluator.h
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Compute the geometrical error in free-surface advection test cases.
+//======================================================================================================================
+
+#include "blockforest/StructuredBlockForest.h"
+
+#include "core/DataTypes.h"
+
+#include "domain_decomposition/BlockDataID.h"
+
+#include "field/iterators/IteratorMacros.h"
+
+namespace walberla
+{
+namespace free_surface
+{
+template< typename FreeSurfaceBoundaryHandling_T, typename FlagField_T, typename ScalarField_T >
+class GeometricalErrorEvaluator
+{
+ public:
+ GeometricalErrorEvaluator(const std::weak_ptr< StructuredBlockForest >& blockForest,
+ const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
+ const ConstBlockDataID& initialfillFieldID, const ConstBlockDataID& fillFieldID,
+ uint_t frequency, const std::shared_ptr< real_t >& geometricalError)
+ : blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling),
+ initialFillFieldID_(initialfillFieldID), fillFieldID_(fillFieldID), frequency_(frequency),
+ geometricalError_(geometricalError), 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))
+ {
+ computeInitialFillLevelSum(blockForest, freeSurfaceBoundaryHandling);
+ computeError(blockForest, freeSurfaceBoundaryHandling);
+ }
+ else
+ {
+ // only evaluate in given frequencies
+ if (executionCounter_ % frequency_ == uint_c(0)) { computeError(blockForest, freeSurfaceBoundaryHandling); }
+ }
+
+ ++executionCounter_;
+ }
+
+ void computeInitialFillLevelSum(
+ 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();
+
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ const ScalarField_T* const initialfillField = blockIt->getData< const ScalarField_T >(initialFillFieldID_);
+ const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
+
+ // avoid OpenMP here because initialFillLevelSum_ is a class member and not a regular variable
+ WALBERLA_FOR_ALL_CELLS_OMP(initialfillFieldIt, initialfillField, flagFieldIt, flagField, omp critical, {
+ if (flagInfo.isInterface(flagFieldIt) || flagInfo.isLiquid(flagFieldIt))
+ {
+ initialFillLevelSum_ += *initialfillFieldIt;
+ }
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ mpi::allReduceInplace< real_t >(initialFillLevelSum_, mpi::SUM);
+ }
+
+ void computeError(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 geometricalError = real_c(0);
+ real_t fillLevelSum = real_c(0);
+
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ const ScalarField_T* const initialfillField = blockIt->getData< const ScalarField_T >(initialFillFieldID_);
+ const ScalarField_T* const fillField = blockIt->getData< const ScalarField_T >(fillFieldID_);
+ const FlagField_T* const flagField = blockIt->getData< const FlagField_T >(flagFieldID);
+
+ WALBERLA_FOR_ALL_CELLS_OMP(initialfillFieldIt, initialfillField, fillFieldIt, fillField, flagFieldIt,
+ flagField, omp parallel for schedule(static) reduction(+:geometricalError)
+ reduction(+:fillLevelSum), {
+ if (flagInfo.isInterface(flagFieldIt) || flagInfo.isLiquid(flagFieldIt))
+ {
+ geometricalError += real_c(std::abs(*initialfillFieldIt - *fillFieldIt));
+ fillLevelSum += *fillFieldIt;
+ }
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ mpi::allReduceInplace< real_t >(geometricalError, mpi::SUM);
+
+ // compute L1 norms
+ *geometricalError_ = geometricalError / initialFillLevelSum_;
+ }
+
+ private:
+ std::weak_ptr< StructuredBlockForest > blockForest_;
+ std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
+ ConstBlockDataID initialFillFieldID_;
+ ConstBlockDataID fillFieldID_;
+ uint_t frequency_;
+ std::shared_ptr< real_t > geometricalError_;
+
+ uint_t executionCounter_;
+ real_t initialFillLevelSum_ = real_c(0);
+}; // class GeometricalErrorEvaluator
+
+} // namespace free_surface
+} // namespace walberla
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/BubblyPoiseuille.cpp b/apps/showcases/FreeSurface/BubblyPoiseuille.cpp
index e27edf893862b3a6569b631484d507e662c22328..e54117c159c3ac0dc4de3f507e642ead83f474af 100644
--- a/apps/showcases/FreeSurface/BubblyPoiseuille.cpp
+++ b/apps/showcases/FreeSurface/BubblyPoiseuille.cpp
@@ -28,6 +28,7 @@
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
@@ -195,8 +196,10 @@ int main(int argc, char** argv)
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
@@ -327,6 +330,14 @@ int main(int argc, char** argv)
loadBalancingFrequency, printLoadBalancingStatistics);
timeloop.addFuncAfterTimeStep(loadBalancer, "Sweep: load balancing");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -348,11 +359,13 @@ int main(int argc, char** argv)
for (uint_t t = uint_c(0); t != timesteps; ++t)
{
- if (t % uint_c(real_c(timesteps / 100)) == uint_c(0))
+ timeloop.singleStep(timingPool, true);
+
+ if (t % evaluationFrequency == uint_c(0))
{
- WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep = " << t);
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass
+ << "\n\t\texcess mass = " << *excessMass);
}
- timeloop.singleStep(timingPool, true);
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
diff --git a/apps/showcases/FreeSurface/BubblyPoiseuille.prm b/apps/showcases/FreeSurface/BubblyPoiseuille.prm
index e4fb01b3df536afc8512ee9bb5b59711a8aee4fa..5b3c7728fe3a3b21f29fcf66edc19cfad9a4a4c3 100644
--- a/apps/showcases/FreeSurface/BubblyPoiseuille.prm
+++ b/apps/showcases/FreeSurface/BubblyPoiseuille.prm
@@ -39,6 +39,7 @@ ModelParameters
EvaluationParameters
{
+ evaluationFrequency 5000;
performanceLogFrequency 5000;
}
diff --git a/apps/showcases/FreeSurface/CMakeLists.txt b/apps/showcases/FreeSurface/CMakeLists.txt
index b0dd586089fea55b83e95e1163772371a6b030cf..34bed9ecd86323ffd5f914a757ca1d4453d4fbe7 100644
--- a/apps/showcases/FreeSurface/CMakeLists.txt
+++ b/apps/showcases/FreeSurface/CMakeLists.txt
@@ -39,14 +39,14 @@ if( WALBERLA_BUILD_WITH_CODEGEN )
GravityWaveLatticeModelGeneration)
endif()
+waLBerla_add_executable(NAME MovingDrop
+ FILES MovingDrop.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
+
waLBerla_add_executable(NAME RisingBubble
FILES RisingBubble.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
waLBerla_add_executable(NAME TaylorBubble
FILES TaylorBubble.cpp
- DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
-
-waLBerla_add_executable(NAME ZalesakDisk
- FILES ZalesakDisk.cpp
DEPENDS blockforest boundary core domain_decomposition field lbm postprocessing timeloop vtk)
\ No newline at end of file
diff --git a/apps/showcases/FreeSurface/CapillaryWave.cpp b/apps/showcases/FreeSurface/CapillaryWave.cpp
index 64148c5b313fbb4f1c9e31eed5b58b004d49d6ac..376a8d9cb0810f76a58cc1485bf7cb4bedf67c55 100644
--- a/apps/showcases/FreeSurface/CapillaryWave.cpp
+++ b/apps/showcases/FreeSurface/CapillaryWave.cpp
@@ -391,6 +391,14 @@ int main(int argc, char** argv)
dynamicsHandler.addSweeps(timeloop);
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add load balancing
LoadBalancer< FlagField_T, CommunicationStencil_T, LatticeModelStencil_T > loadBalancer(
blockForest, communication, pdfCommunication, bubbleModel, uint_c(50), uint_c(10), uint_c(5),
@@ -436,9 +444,10 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ tNonDimensional, positionNonDimensional };
if (t % evaluationFrequency == uint_c(0))
{
- WALBERLA_LOG_DEVEL("time step = " << t);
- WALBERLA_LOG_DEVEL("\t\ttNonDimensional = " << tNonDimensional
- << "\n\t\tpositionNonDimensional = " << positionNonDimensional);
+ WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\ttNonDimensional = " << tNonDimensional
+ << "\n\t\tpositionNonDimensional = " << positionNonDimensional
+ << "\n\t\ttotal mass = " << *totalMass
+ << "\n\t\texcess mass = " << *excessMass);
writeVectorToFile(resultVector, filename);
}
}
diff --git a/apps/showcases/FreeSurface/DamBreakCylindrical.cpp b/apps/showcases/FreeSurface/DamBreakCylindrical.cpp
index f48e23232d26c404a9ffb70350d307ee591e1b19..0cf75946683d43aec338c23545e4239558257398 100644
--- a/apps/showcases/FreeSurface/DamBreakCylindrical.cpp
+++ b/apps/showcases/FreeSurface/DamBreakCylindrical.cpp
@@ -32,6 +32,7 @@
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
@@ -527,6 +528,14 @@ int main(int argc, char** argv)
evaluationFrequency, columnRadiusSample);
timeloop.addFuncAfterTimeStep(columnRadiusEvaluator, "Evaluator: radius");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -571,13 +580,14 @@ int main(int argc, char** argv)
H = real_c(*currentColumnHeight) / columnHeight;
}
- WALBERLA_LOG_DEVEL_ON_ROOT("time step =" << t);
- WALBERLA_LOG_DEVEL_ON_ROOT("\t\tT = " << T << "\n\t\tZ_mean = " << Z_mean << "\n\t\tZ_max = " << Z_max
- << "\n\t\tZ_min = " << Z_min
- << "\n\t\tZ_stdDeviation = " << Z_stdDeviation << "\n\t\tH = " << H);
-
WALBERLA_ROOT_SECTION()
{
+ WALBERLA_LOG_DEVEL("time step =" << t);
+ WALBERLA_LOG_DEVEL("\t\tT = " << T << "\n\t\tZ_mean = " << Z_mean << "\n\t\tZ_max = " << Z_max
+ << "\n\t\tZ_min = " << Z_min << "\n\t\tZ_stdDeviation = " << Z_stdDeviation
+ << "\n\t\tH = " << H << "\n\t\ttotal mass = " << *totalMass
+ << "\n\t\texcess mass = " << *excessMass);
+
const std::vector< real_t > resultVector{ T, Z_mean, Z_max, Z_min, Z_stdDeviation, H };
writeNumberVector(resultVector, t, filename);
}
diff --git a/apps/showcases/FreeSurface/DamBreakRectangular.cpp b/apps/showcases/FreeSurface/DamBreakRectangular.cpp
index e7abd62e6b9d605db803f8f680173c599cce61d8..fba4e95798adb5973e0f0aaa842ecdc45cf60b47 100644
--- a/apps/showcases/FreeSurface/DamBreakRectangular.cpp
+++ b/apps/showcases/FreeSurface/DamBreakRectangular.cpp
@@ -512,6 +512,14 @@ int main(int argc, char** argv)
blockForest, freeSurfaceBoundaryHandling, domainSize, evaluationFrequency, currentColumnWidth);
timeloop.addFuncAfterTimeStep(widthEvaluator, "Evaluator: column width");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -543,9 +551,14 @@ int main(int argc, char** argv)
const real_t H = real_c(*currentColumnHeight) / columnHeight;
const std::vector< real_t > resultVector{ T, Z, H };
- WALBERLA_LOG_DEVEL_ON_ROOT("time step =" << t);
- WALBERLA_LOG_DEVEL_ON_ROOT("\t\tT = " << T << "\n\t\tZ = " << Z << "\n\t\tH = " << H);
- WALBERLA_ROOT_SECTION() { writeNumberVector(resultVector, t, filename); }
+ WALBERLA_ROOT_SECTION()
+ {
+ WALBERLA_LOG_DEVEL("time step =" << t);
+ WALBERLA_LOG_DEVEL("\t\tT = " << T << "\n\t\tZ = " << Z << "\n\t\tH = " << H << "\n\t\ttotal mass = "
+ << *totalMass << "\n\t\texcess mass = " << *excessMass);
+
+ writeNumberVector(resultVector, t, filename);
+ }
// simulation is considered converged
if (Z >= real_c(domainSize[0]) / columnWidth - real_c(0.5))
diff --git a/apps/showcases/FreeSurface/DropImpact.cpp b/apps/showcases/FreeSurface/DropImpact.cpp
index 896c5b098b6a0e44edf6e6983b650aa25cdb4c90..61da9474b18682bf4e26e075d29bc917172df934 100644
--- a/apps/showcases/FreeSurface/DropImpact.cpp
+++ b/apps/showcases/FreeSurface/DropImpact.cpp
@@ -29,6 +29,7 @@
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
@@ -202,8 +203,10 @@ int main(int argc, char** argv)
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
@@ -325,6 +328,14 @@ int main(int argc, char** argv)
loadBalancingFrequency, printLoadBalancingStatistics);
timeloop.addFuncAfterTimeStep(loadBalancer, "Sweep: load balancing");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -346,11 +357,13 @@ int main(int argc, char** argv)
for (uint_t t = uint_c(0); t != timesteps; ++t)
{
- if (t % uint_c(real_c(timesteps / 100)) == uint_c(0))
+ timeloop.singleStep(timingPool, true);
+
+ if (t % evaluationFrequency == uint_c(0))
{
- WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep = " << t);
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass
+ << "\n\t\texcess mass = " << *excessMass);
}
- timeloop.singleStep(timingPool, true);
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
diff --git a/apps/showcases/FreeSurface/DropImpact.prm b/apps/showcases/FreeSurface/DropImpact.prm
index 8437db98ebca2c88b4c2e22d52d3e72d057f7d29..1790dffb48d4d8d5ff3d762e95f4ecbacaad89a1 100644
--- a/apps/showcases/FreeSurface/DropImpact.prm
+++ b/apps/showcases/FreeSurface/DropImpact.prm
@@ -41,6 +41,7 @@ ModelParameters
EvaluationParameters
{
+ evaluationFrequency 1000;
performanceLogFrequency 3000;
}
diff --git a/apps/showcases/FreeSurface/DropWetting.cpp b/apps/showcases/FreeSurface/DropWetting.cpp
index bdba5908629d95c20904d5ab157a19b90b37ebca..a55ad88083021fa3cc3f86caa3b7d8248a7e199d 100644
--- a/apps/showcases/FreeSurface/DropWetting.cpp
+++ b/apps/showcases/FreeSurface/DropWetting.cpp
@@ -29,6 +29,7 @@
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
@@ -368,6 +369,14 @@ int main(int argc, char** argv)
blockForest, freeSurfaceBoundaryHandling, fillFieldID, dropHeight, evaluationFrequency);
timeloop.addFuncAfterTimeStep(dropHeightEvaluator, "Evaluator: drop height");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -395,8 +404,7 @@ int main(int argc, char** argv)
// check convergence
if (t % evaluationFrequency == uint_c(0))
{
- WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t)
- WALBERLA_LOG_DEVEL_ON_ROOT("\t\tdrop height = " << *dropHeight)
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\tdrop height = " << *dropHeight)
if (std::abs(formerDropHeight - *dropHeight) / *dropHeight < convergenceThreshold)
{
WALBERLA_LOG_DEVEL_ON_ROOT("Final converged drop height=" << *dropHeight);
diff --git a/apps/showcases/FreeSurface/GravityWave.cpp b/apps/showcases/FreeSurface/GravityWave.cpp
index acca70e45fd4bb0fa75e8808cd639ece5353a59e..702682bcf0aa37a77ec8e54dea1ba2155810c66c 100644
--- a/apps/showcases/FreeSurface/GravityWave.cpp
+++ b/apps/showcases/FreeSurface/GravityWave.cpp
@@ -554,6 +554,14 @@ int main(int argc, char** argv)
evaluationFrequency, symmetryNorm);
timeloop.addFuncAfterTimeStep(symmetryEvaluator, "Evaluator: symmetry norm");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -586,10 +594,11 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ tNonDimensional, positionNonDimensional, *symmetryNorm };
if (t % evaluationFrequency == uint_c(0))
{
- WALBERLA_LOG_DEVEL("time step = " << t);
- WALBERLA_LOG_DEVEL("\t\ttNonDimensional = " << tNonDimensional
- << "\n\t\tpositionNonDimensional = " << positionNonDimensional
- << "\n\t\tsymmetryNorm = " << *symmetryNorm);
+ WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\ttNonDimensional = " << tNonDimensional
+ << "\n\t\tpositionNonDimensional = " << positionNonDimensional
+ << "\n\t\tsymmetryNorm = " << *symmetryNorm << "\n\t\ttotal mass = "
+ << *totalMass << "\n\t\texcess mass = " << *excessMass);
+
writeVectorToFile(resultVector, filename);
}
}
diff --git a/apps/showcases/FreeSurface/GravityWaveCodegen.cpp b/apps/showcases/FreeSurface/GravityWaveCodegen.cpp
index dd603d62c4a1552ad525d2fdceaf35b29ab98880..65d18594b5b2c9f024ef5347b871ff2b9b79501a 100644
--- a/apps/showcases/FreeSurface/GravityWaveCodegen.cpp
+++ b/apps/showcases/FreeSurface/GravityWaveCodegen.cpp
@@ -516,6 +516,14 @@ int main(int argc, char** argv)
evaluationFrequency, symmetryNorm);
timeloop.addFuncAfterTimeStep(symmetryEvaluator, "Evaluator: symmetry norm");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T,
true, VectorFieldFlattened_T >(
@@ -549,10 +557,11 @@ int main(int argc, char** argv)
const std::vector< real_t > resultVector{ tNonDimensional, positionNonDimensional, *symmetryNorm };
if (t % evaluationFrequency == uint_c(0))
{
- WALBERLA_LOG_DEVEL("time step = " << t);
- WALBERLA_LOG_DEVEL("\t\ttNonDimensional = " << tNonDimensional
- << "\n\t\tpositionNonDimensional = " << positionNonDimensional
- << "\n\t\tsymmetryNorm = " << *symmetryNorm);
+ WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\ttNonDimensional = " << tNonDimensional
+ << "\n\t\tpositionNonDimensional = " << positionNonDimensional
+ << "\n\t\tsymmetryNorm = " << *symmetryNorm << "\n\t\ttotal mass = "
+ << *totalMass << "\n\t\texcess mass = " << *excessMass);
+
writeVectorToFile(resultVector, filename);
}
}
diff --git a/apps/showcases/FreeSurface/MovingDrop.cpp b/apps/showcases/FreeSurface/MovingDrop.cpp
index 352a85c72cc7239d2608c62567677f7db3d5c324..4fb4c49d34f3031f74c480795d9873c90a7cc91b 100644
--- a/apps/showcases/FreeSurface/MovingDrop.cpp
+++ b/apps/showcases/FreeSurface/MovingDrop.cpp
@@ -27,6 +27,7 @@
#include "lbm/field/AddToStorage.h"
#include "lbm/free_surface/LoadBalancing.h"
#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
#include "lbm/free_surface/VtkWriter.h"
#include "lbm/free_surface/bubble_model/Geometry.h"
#include "lbm/free_surface/dynamics/SurfaceDynamicsHandler.h"
@@ -189,8 +190,10 @@ int main(int argc, char** argv)
// read evaluation parameters from parameter file
const auto evaluationParameters = walberlaEnv.config()->getOneBlock("EvaluationParameters");
const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
WALBERLA_LOG_DEVEL_VAR_ON_ROOT(performanceLogFrequency);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(evaluationFrequency);
// create non-uniform block forest (non-uniformity required for load balancing)
const std::shared_ptr< StructuredBlockForest > blockForest =
@@ -286,6 +289,14 @@ int main(int argc, char** argv)
loadBalancingFrequency, printLoadBalancingStatistics);
timeloop.addFuncAfterTimeStep(loadBalancer, "Sweep: load balancing");
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
+ const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
+ timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
+
// add VTK output
addVTKOutput< LatticeModel_T, FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T, VectorField_T >(
blockForest, timeloop, walberlaEnv.config(), flagInfo, pdfFieldID, flagFieldID, fillFieldID, forceDensityFieldID,
@@ -307,9 +318,14 @@ int main(int argc, char** argv)
for (uint_t t = uint_c(0); t != timesteps; ++t)
{
- if (t % uint_c(100) == uint_c(0)) { WALBERLA_LOG_DEVEL_ON_ROOT("Performing timestep=" << t); }
timeloop.singleStep(timingPool, true);
+ if (t % evaluationFrequency == uint_c(0))
+ {
+ WALBERLA_LOG_DEVEL_ON_ROOT("time step = " << t << "\n\t\ttotal mass = " << *totalMass
+ << "\n\t\texcess mass = " << *excessMass);
+ }
+
if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
}
diff --git a/apps/showcases/FreeSurface/MovingDrop.prm b/apps/showcases/FreeSurface/MovingDrop.prm
index dceed84e71501ad90e5d8020059792531f9a8d13..79718c8c060feceff66ce13e972a605283a07c7a 100644
--- a/apps/showcases/FreeSurface/MovingDrop.prm
+++ b/apps/showcases/FreeSurface/MovingDrop.prm
@@ -1,17 +1,17 @@
BlockForestParameters
{
- cellsPerBlock < 10, 10, 10 >;
+ cellsPerBlock < 20, 20, 20 >;
periodicity < 1, 1, 1 >;
loadBalancingFrequency 0;
- printLoadBalancingStatistics true;
+ printLoadBalancingStatistics false;
}
DomainParameters
{
- dropDiameter 50;
+ dropDiameter 20;
dropCenterFactor < 1, 1, 1 >; // values multiplied with dropDiameter
- poolHeightFactor 0; // value multiplied with dropDiameter
- domainSizeFactor < 2, 2, 4 >; // values multiplied with dropDiameter
+ poolHeightFactor 0; // value multiplied with dropDiameter
+ domainSizeFactor < 2, 2, 4 >; // values multiplied with dropDiameter
}
PhysicsParameters
@@ -31,14 +31,16 @@ ModelParameters
excessMassDistributionModel EvenlyAllInterface;
curvatureModel FiniteDifferenceMethod;
useSimpleMassExchange false;
- enableBubbleModel false;
- enableBubbleSplits false; // only used if enableBubbleModel=true
cellConversionThreshold 1e-2;
cellConversionForceThreshold 1e-1;
+
+ enableBubbleModel false;
+ enableBubbleSplits false; // only used if enableBubbleModel=true
}
EvaluationParameters
{
+ evaluationFrequency 100;
performanceLogFrequency 10000;
}
@@ -59,7 +61,7 @@ BoundaryParameters
MeshOutputParameters
{
- writeFrequency 10000;
+ writeFrequency 1000;
baseFolder mesh-out;
}
@@ -67,7 +69,7 @@ VTK
{
fluid_field
{
- writeFrequency 10000;
+ writeFrequency 1000;
ghostLayers 0;
baseFolder vtk-out;
samplingResolution 1;
diff --git a/apps/showcases/FreeSurface/RisingBubble.cpp b/apps/showcases/FreeSurface/RisingBubble.cpp
index 8fc7738761357e34d9cea02a73399d4572dddcec..6699fdf1eac647bc87f87bf3e1e9bbf8150679ee 100644
--- a/apps/showcases/FreeSurface/RisingBubble.cpp
+++ b/apps/showcases/FreeSurface/RisingBubble.cpp
@@ -376,10 +376,12 @@ int main(int argc, char** argv)
blockForest, freeSurfaceBoundaryHandling, evaluationFrequency, centerOfMass);
timeloop.addFuncAfterTimeStep(centerOfMassComputer, "Evaluator: center of mass");
- // add computation of total mass
- const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
- blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, evaluationFrequency, totalMass);
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
// add VTK output
@@ -419,10 +421,10 @@ int main(int argc, char** argv)
const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) /
(riseVelocity * riseVelocity);
- WALBERLA_LOG_DEVEL("time step = " << t);
- WALBERLA_LOG_DEVEL("\t\tcenterOfMass = " << *centerOfMass << "\n\t\triseVelocity = " << riseVelocity
- << "\n\t\tdragForce = " << dragForce);
- WALBERLA_LOG_DEVEL("\t\ttotalMass = " << *totalMass);
+ WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass
+ << "\n\t\triseVelocity = " << riseVelocity
+ << "\n\t\tdragForce = " << dragForce << "\n\t\ttotalMass = " << *totalMass
+ << "\n\t\texcessMass = " << *excessMass);
const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce };
diff --git a/apps/showcases/FreeSurface/TaylorBubble.cpp b/apps/showcases/FreeSurface/TaylorBubble.cpp
index a045cd6acca3b405e60168dea64489812984259a..ce2a45302d06494fb588819d1ecb43cbedec9afb 100644
--- a/apps/showcases/FreeSurface/TaylorBubble.cpp
+++ b/apps/showcases/FreeSurface/TaylorBubble.cpp
@@ -405,9 +405,12 @@ int main(int argc, char** argv)
blockForest, freeSurfaceBoundaryHandling, evaluationFrequency, centerOfMass);
timeloop.addFuncAfterTimeStep(centerOfMassComputer, "Evaluator: center of mass");
- const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ // add evaluator for total and excessive mass (mass that is currently undistributed)
+ const std::shared_ptr< real_t > totalMass = std::make_shared< real_t >(real_c(0));
+ const std::shared_ptr< real_t > excessMass = std::make_shared< real_t >(real_c(0));
const TotalMassComputer< FreeSurfaceBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T > totalMassComputer(
- blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, evaluationFrequency, totalMass);
+ blockForest, freeSurfaceBoundaryHandling, pdfFieldID, fillFieldID, dynamicsHandler.getConstExcessMassFieldID(),
+ evaluationFrequency, totalMass, excessMass);
timeloop.addFuncAfterTimeStep(totalMassComputer, "Evaluator: total mass");
// add VTK output
@@ -447,20 +450,20 @@ int main(int argc, char** argv)
const real_t dragForce = real_c(4) / real_c(3) * gravitationalAccelerationZ * real_c(bubbleDiameter) /
(riseVelocity * riseVelocity);
- WALBERLA_LOG_DEVEL("time step = " << t);
- WALBERLA_LOG_DEVEL("\t\tcenterOfMass = " << *centerOfMass << "\n\t\triseVelocity = " << riseVelocity
- << "\n\t\tdragForce = " << dragForce);
- WALBERLA_LOG_DEVEL("\t\ttotalMass = " << *totalMass);
+ WALBERLA_LOG_DEVEL("time step = " << t << "\n\t\tcenterOfMass = " << *centerOfMass
+ << "\n\t\triseVelocity = " << riseVelocity
+ << "\n\t\tdragForce = " << dragForce << "\n\t\ttotalMass = " << *totalMass
+ << "\n\t\texcessMass = " << *excessMass);
const std::vector< real_t > resultVector{ (*centerOfMass)[2], riseVelocity, dragForce };
writeVectorToFile(resultVector, t, filename);
}
-
- timestepOld = t;
- centerOfMassOld = *centerOfMass;
}
+ timestepOld = t;
+ centerOfMassOld = *centerOfMass;
+
// stop simulation before bubble hits the top wall
if ((*centerOfMass)[2] > stoppingHeight) { break; }
diff --git a/src/lbm/free_surface/TotalMassComputer.h b/src/lbm/free_surface/TotalMassComputer.h
index 192ec875594b98dc87bff0e1418c200544965590..7b3fda3e83162cc3657501ecca5fa833649a7ffc 100644
--- a/src/lbm/free_surface/TotalMassComputer.h
+++ b/src/lbm/free_surface/TotalMassComputer.h
@@ -50,10 +50,10 @@ class TotalMassComputer
const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
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 >& totalMass, const std::shared_ptr< real_t >& excessMass)
: blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling), pdfFieldID_(pdfFieldID),
- fillFieldID_(fillFieldID), excessMassFieldID_(excessMassFieldID), totalMass_(totalMass), frequency_(frequency),
- executionCounter_(uint_c(0))
+ fillFieldID_(fillFieldID), excessMassFieldID_(excessMassFieldID), totalMass_(totalMass),
+ excessMass_(excessMass), frequency_(frequency), executionCounter_(uint_c(0))
{}
void operator()()
@@ -66,11 +66,10 @@ class TotalMassComputer
auto freeSurfaceBoundaryHandling = freeSurfaceBoundaryHandling_.lock();
WALBERLA_CHECK_NOT_NULLPTR(freeSurfaceBoundaryHandling);
- if (executionCounter_ == uint_c(0)) { computeMass(blockForest, freeSurfaceBoundaryHandling); }
- else
+ // only evaluate in given frequencies
+ if (executionCounter_ % frequency_ == uint_c(0) || executionCounter_ == uint_c(0))
{
- // only evaluate in given frequencies
- if (executionCounter_ % frequency_ == uint_c(0)) { computeMass(blockForest, freeSurfaceBoundaryHandling); }
+ computeMass(blockForest, freeSurfaceBoundaryHandling);
}
++executionCounter_;
@@ -82,7 +81,8 @@ class TotalMassComputer
const BlockDataID flagFieldID = freeSurfaceBoundaryHandling->getFlagFieldID();
const typename FreeSurfaceBoundaryHandling_T::FlagInfo_T& flagInfo = freeSurfaceBoundaryHandling->getFlagInfo();
- real_t mass = real_c(0);
+ real_t mass = real_c(0);
+ real_t excessMass = real_c(0);
for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
{
@@ -104,6 +104,8 @@ class TotalMassComputer
{
const real_t density = pdfField->getDensity(pdfFieldIt.cell());
mass += *fillFieldIt * density + *excessMassFieldIt;
+
+ if (excessMass_ != nullptr) { excessMass += *excessMassFieldIt; }
}
}) // WALBERLA_FOR_ALL_CELLS_OMP
}
@@ -122,8 +124,13 @@ class TotalMassComputer
}
mpi::allReduceInplace< real_t >(mass, mpi::SUM);
-
*totalMass_ = mass;
+
+ if (excessMass_ != nullptr)
+ {
+ mpi::allReduceInplace< real_t >(excessMass, mpi::SUM);
+ *excessMass_ = excessMass;
+ }
};
private:
@@ -135,6 +142,7 @@ class TotalMassComputer
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_;
diff --git a/src/lbm/free_surface/VtkWriter.h b/src/lbm/free_surface/VtkWriter.h
index a27103510f6831a96bf7612391518626c6f45985..726e0ac735ef4495d234532a7e3eb3e52c0372a8 100644
--- a/src/lbm/free_surface/VtkWriter.h
+++ b/src/lbm/free_surface/VtkWriter.h
@@ -56,7 +56,7 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
const BlockDataID& forceDensityFieldID, const BlockDataID& curvatureFieldID,
const BlockDataID& normalFieldID, const BlockDataID& obstacleNormalFieldID)
{
- using value_type = typename FlagField_T::value_type;
+ using value_type = typename FlagField_T::value_type;
const auto blockForest = blockForestPtr.lock();
WALBERLA_CHECK_NOT_NULLPTR(blockForest);
@@ -89,27 +89,34 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
std::make_shared< VTKWriter< VectorField_T, float > >(obstacleNormalFieldID, "obstacle_normal"));
if constexpr (useCodegen)
{
- writers.push_back(
- std::make_shared< VTKWriter< VectorFieldFlattened_T, float > >(forceDensityFieldID, "force_density"));
+ if (forceDensityFieldID != BlockDataID())
+ {
+ writers.push_back(
+ std::make_shared< VTKWriter< VectorFieldFlattened_T, float > >(forceDensityFieldID, "force_density"));
+ }
}
else
{
- writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(forceDensityFieldID, "force_density"));
+ if (forceDensityFieldID != BlockDataID())
+ {
+ writers.push_back(
+ std::make_shared< VTKWriter< VectorField_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));
+ 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);
diff --git a/src/lbm/free_surface/dynamics/ExcessMassDistributionModel.h b/src/lbm/free_surface/dynamics/ExcessMassDistributionModel.h
index 2e4d0b798ca3fb8bf6e325f38abfaa30119291fc..a4b86c43c64ca2c22b7326b12a779be846da5b86 100644
--- a/src/lbm/free_surface/dynamics/ExcessMassDistributionModel.h
+++ b/src/lbm/free_surface/dynamics/ExcessMassDistributionModel.h
@@ -63,15 +63,21 @@ namespace free_surface
* cells, i.e., cells that are non-newly converted to interface. Falls back to WeightedAllInterface if not
* applicable.
*
- * - EvenlyLiquidAndAllInterface:
+ * - EvenlyAllInterfaceAndLiquid:
* Excess mass is distributed evenly among all neighboring interface and liquid cells (see p.47 in master thesis of
* M. Lehmann, 2019). The excess mass distributed to liquid cells does neither modify the cell's density nor fill
* level. Instead, it is stored in an additional excess mass field. Therefore, not only the converted interface
* cells' excess mass is distributed, but also the excess mass of liquid cells stored in this additional field.
*
- * - EvenlyLiquidAndAllInterfacePreferInterface:
- * Similar to EvenlyLiquidAndAllInterface, however, excess mass is preferably distributed to interface cells. It is
+ * - EvenlyAllInterfaceFallbackLiquid:
+ * Similar to EvenlyAllInterfaceAndLiquid, however, excess mass is preferably distributed to interface cells. It is
* distributed to liquid cells only if there are no neighboring interface cells available.
+ *
+ * - EvenlyNewInterfaceFallbackLiquid:
+ * Similar to EvenlyAllInterfaceFallbackLiquid, however, excess mass is preferably distributed to newly
+ * converted interface cells. If there are no newly converted interface cells, the excess mass is distributed to
+ * old interface cells. The excess mass is distributed to neighboring liquid cells only if there are no neighboring
+ * interface cells available.
* ********************************************************************************************************************/
class ExcessMassDistributionModel
{
@@ -83,8 +89,9 @@ class ExcessMassDistributionModel
WeightedAllInterface,
WeightedNewInterface,
WeightedOldInterface,
- EvenlyLiquidAndAllInterface,
- EvenlyLiquidAndAllInterfacePreferInterface
+ EvenlyAllInterfaceAndLiquid,
+ EvenlyAllInterfaceFallbackLiquid,
+ EvenlyNewInterfaceFallbackLiquid
};
ExcessMassDistributionModel(const std::string& modelName) : modelName_(modelName), modelType_(chooseType(modelName))
@@ -107,9 +114,11 @@ class ExcessMassDistributionModel
break;
case ExcessMassModel::WeightedOldInterface:
break;
- case ExcessMassModel::EvenlyLiquidAndAllInterface:
+ case ExcessMassModel::EvenlyAllInterfaceAndLiquid:
+ break;
+ case ExcessMassModel::EvenlyAllInterfaceFallbackLiquid:
break;
- case ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface:
+ case ExcessMassModel::EvenlyNewInterfaceFallbackLiquid:
break;
}
}
@@ -130,10 +139,12 @@ class ExcessMassDistributionModel
modelType_ == ExcessMassModel::WeightedNewInterface || modelType_ == ExcessMassModel::WeightedOldInterface;
}
- inline bool isEvenlyLiquidAndAllInterfacePreferInterfaceType() const
+ inline bool isEvenlyAllInterfaceFallbackLiquidType() const
{
- return modelType_ == ExcessMassModel::EvenlyLiquidAndAllInterface ||
- modelType_ == ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface;
+ return modelType_ == ExcessMassModel::EvenlyAllInterfaceAndLiquid ||
+ modelType_ == ExcessMassModel::EvenlyAllInterfaceFallbackLiquid ||
+ modelType_ == ExcessMassModel::EvenlyNewInterfaceFallbackLiquid;
+ ;
}
static inline std::initializer_list< const ExcessMassModel > getTypeIterator() { return listOfAllEnums; }
@@ -153,14 +164,19 @@ class ExcessMassDistributionModel
if (!string_icompare(modelName, "WeightedOldInterface")) { return ExcessMassModel::WeightedOldInterface; }
- if (!string_icompare(modelName, "EvenlyLiquidAndAllInterface"))
+ if (!string_icompare(modelName, "EvenlyAllInterfaceAndLiquid"))
+ {
+ return ExcessMassModel::EvenlyAllInterfaceAndLiquid;
+ }
+
+ if (!string_icompare(modelName, "EvenlyAllInterfaceFallbackLiquid"))
{
- return ExcessMassModel::EvenlyLiquidAndAllInterface;
+ return ExcessMassModel::EvenlyAllInterfaceFallbackLiquid;
}
- if (!string_icompare(modelName, "EvenlyLiquidAndAllInterfacePreferInterface"))
+ if (!string_icompare(modelName, "EvenlyNewInterfaceFallbackLiquid"))
{
- return ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface;
+ return ExcessMassModel::EvenlyNewInterfaceFallbackLiquid;
}
WALBERLA_ABORT("The specified PDF reinitialization model " << modelName << " is not available.");
@@ -190,11 +206,14 @@ class ExcessMassDistributionModel
modelName = "WeightedOldInterface";
break;
- case ExcessMassModel::EvenlyLiquidAndAllInterface:
- modelName = "EvenlyLiquidAndAllInterface";
+ case ExcessMassModel::EvenlyAllInterfaceAndLiquid:
+ modelName = "EvenlyAllInterfaceAndLiquid";
+ break;
+ case ExcessMassModel::EvenlyAllInterfaceFallbackLiquid:
+ modelName = "EvenlyAllInterfaceFallbackLiquid";
break;
- case ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface:
- modelName = "EvenlyLiquidAndAllInterfacePreferInterface";
+ case ExcessMassModel::EvenlyNewInterfaceFallbackLiquid:
+ modelName = "EvenlyNewInterfaceFallbackLiquid";
break;
}
return modelName;
@@ -203,10 +222,15 @@ class ExcessMassDistributionModel
std::string modelName_;
ExcessMassModel modelType_;
static constexpr std::initializer_list< const ExcessMassModel > listOfAllEnums = {
- ExcessMassModel::EvenlyAllInterface, ExcessMassModel::EvenlyNewInterface,
- ExcessMassModel::EvenlyOldInterface, ExcessMassModel::WeightedAllInterface,
- ExcessMassModel::WeightedNewInterface, ExcessMassModel::WeightedOldInterface,
- ExcessMassModel::EvenlyLiquidAndAllInterface, ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface
+ ExcessMassModel::EvenlyAllInterface,
+ ExcessMassModel::EvenlyNewInterface,
+ ExcessMassModel::EvenlyOldInterface,
+ ExcessMassModel::WeightedAllInterface,
+ ExcessMassModel::WeightedNewInterface,
+ ExcessMassModel::WeightedOldInterface,
+ ExcessMassModel::EvenlyAllInterfaceAndLiquid,
+ ExcessMassModel::EvenlyAllInterfaceFallbackLiquid,
+ ExcessMassModel::EvenlyNewInterfaceFallbackLiquid
};
}; // class ExcessMassDistributionModel
diff --git a/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.h b/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.h
index ab9efe397121b33391eceedefc407811aec16ac8..34ebe1b2055bbe04a09700c0a7b983b1fe9485e8 100644
--- a/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.h
+++ b/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.h
@@ -71,9 +71,8 @@ class ExcessMassDistributionSweepBase
/********************************************************************************************************************
* Determines the number of a cell's neighboring liquid and interface cells.
*******************************************************************************************************************/
- void getNumberOfEvenlyLiquidAndAllInterfacePreferInterfaceNeighbors(const FlagField_T* flagField, const Cell& cell,
- uint_t& liquidNeighbors,
- uint_t& interfaceNeighbors);
+ void getNumberOfLiquidAndInterfaceNeighbors(const FlagField_T* flagField, const Cell& cell, uint_t& liquidNeighbors,
+ uint_t& interfaceNeighbors, uint_t& newInterfaceNeighbors);
ExcessMassDistributionModel excessMassDistributionModel_;
BlockDataID fillFieldID_;
@@ -171,6 +170,8 @@ class ExcessMassDistributionSweepInterfaceWeighted
* Distribute the excess mass evenly among
* - all neighboring liquid and interface cells (see p. 47 in master thesis of M. Lehmann, 2019)
* - all neighboring interface cells and only to liquid cells if there exists no neighboring interface cell
+ * - new neighboring interface cells, if not available to old interface cells and only to liquid cells if there exists
+ * no neighboring interface cell
*
* Neither the fill level, nor the density of liquid cells is modified. Instead, the excess mass is stored in an
* additional field.
diff --git a/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.impl.h b/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.impl.h
index ba30c8445d19276136859641d97673fc8faff98c..12ffe03179912214a89b1f1110004d4e3cd331ce 100644
--- a/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.impl.h
+++ b/src/lbm/free_surface/dynamics/ExcessMassDistributionSweep.impl.h
@@ -79,18 +79,23 @@ void ExcessMassDistributionSweepInterfaceEvenly< LatticeModel_T, FlagField_T, Sc
template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename VectorField_T >
void ExcessMassDistributionSweepBase< LatticeModel_T, FlagField_T, ScalarField_T, VectorField_T >::
- getNumberOfEvenlyLiquidAndAllInterfacePreferInterfaceNeighbors(const FlagField_T* flagField, const Cell& cell,
- uint_t& liquidNeighbors, uint_t& interfaceNeighbors)
+ getNumberOfLiquidAndInterfaceNeighbors(const FlagField_T* flagField, const Cell& cell, uint_t& liquidNeighbors,
+ uint_t& interfaceNeighbors, uint_t& newInterfaceNeighbors)
{
- interfaceNeighbors = uint_c(0);
- liquidNeighbors = uint_c(0);
+ newInterfaceNeighbors = uint_c(0);
+ interfaceNeighbors = uint_c(0);
+ liquidNeighbors = uint_c(0);
for (auto d = LatticeModel_T::Stencil::beginNoCenter(); d != LatticeModel_T::Stencil::end(); ++d)
{
const Cell neighborCell = Cell(cell.x() + d.cx(), cell.y() + d.cy(), cell.z() + d.cz());
auto neighborFlags = flagField->get(neighborCell);
- if (isFlagSet(neighborFlags, flagInfo_.interfaceFlag)) { ++interfaceNeighbors; }
+ if (isFlagSet(neighborFlags, flagInfo_.interfaceFlag))
+ {
+ ++interfaceNeighbors;
+ if (isFlagSet(neighborFlags, flagInfo_.convertedFlag)) { ++newInterfaceNeighbors; }
+ }
else
{
if (isFlagSet(neighborFlags, flagInfo_.liquidFlag)) { ++liquidNeighbors; }
@@ -147,7 +152,8 @@ void ExcessMassDistributionSweepInterfaceEvenly< LatticeModel_T, FlagField_T, Sc
if (interfaceNeighbors == uint_c(0))
{
- WALBERLA_LOG_WARNING("No interface cell is in the neighborhood to distribute excess mass to. Mass is lost.");
+ WALBERLA_LOG_WARNING(
+ "No interface cell is in the neighborhood to distribute excess mass to. Mass is lost/gained.");
return;
}
@@ -285,7 +291,8 @@ void ExcessMassDistributionSweepInterfaceWeighted< LatticeModel_T, FlagField_T,
if (interfaceNeighbors == uint_c(0))
{
- WALBERLA_LOG_WARNING("No interface cell is in the neighborhood to distribute excess mass to. Mass is lost.");
+ WALBERLA_LOG_WARNING(
+ "No interface cell is in the neighborhood to distribute excess mass to. Mass is lost/gained.");
return;
}
@@ -505,7 +512,9 @@ void ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T,
// calculate excess fill level
const real_t excessFill = newGas ? fillField->get(cell) : (fillField->get(cell) - real_c(1.0));
- // store excess mass such that it can be distributed below
+ // store excess mass such that it can be distributed below (no += here because cell was an interface cell
+ // that can not have an excess mass stored in the field; any excess mass is added to the interface cell's
+ // fill level)
srcExcessMassField->get(cell) = excessFill * pdfField->getDensity(cell);
if (newGas) { fillField->get(cell) = real_c(0.0); }
@@ -533,28 +542,41 @@ void ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T,
using Base_T = ExcessMassDistributionSweepBase_T;
// get number of liquid and interface neighbors
- uint_t liquidNeighbors = uint_c(0);
- uint_t interfaceNeighbors = uint_c(0);
- Base_T::getNumberOfEvenlyLiquidAndAllInterfacePreferInterfaceNeighbors(flagField, cell, liquidNeighbors,
- interfaceNeighbors);
- const uint_t EvenlyLiquidAndAllInterfacePreferInterfaceNeighbors = liquidNeighbors + interfaceNeighbors;
+ uint_t liquidNeighbors = uint_c(0);
+ uint_t interfaceNeighbors = uint_c(0);
+ uint_t newInterfaceNeighbors = uint_c(0);
+ Base_T::getNumberOfLiquidAndInterfaceNeighbors(flagField, cell, liquidNeighbors, interfaceNeighbors,
+ newInterfaceNeighbors);
+ const uint_t liquidAndInterfaceNeighbors = liquidNeighbors + interfaceNeighbors;
- if (EvenlyLiquidAndAllInterfacePreferInterfaceNeighbors == uint_c(0))
+ if (liquidAndInterfaceNeighbors == uint_c(0))
{
WALBERLA_LOG_WARNING(
- "No liquid or interface cell is in the neighborhood to distribute excess mass to. Mass is lost.");
+ "No liquid or interface cell is in the neighborhood to distribute excess mass to. Mass is lost/gained.");
return;
}
- const bool preferInterface =
- Base_T::excessMassDistributionModel_.getModelType() ==
- ExcessMassDistributionModel::ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface &&
- interfaceNeighbors > uint_c(0);
+ // check if there are neighboring new interface cells
+ const bool preferNewInterface = Base_T::excessMassDistributionModel_.getModelType() ==
+ ExcessMassDistributionModel::ExcessMassModel::EvenlyNewInterfaceFallbackLiquid &&
+ newInterfaceNeighbors > uint_c(0);
+
+ // check if there are neighboring interface cells
+ const bool preferInterface = (Base_T::excessMassDistributionModel_.getModelType() ==
+ ExcessMassDistributionModel::ExcessMassModel::EvenlyAllInterfaceFallbackLiquid ||
+ !preferNewInterface) &&
+ interfaceNeighbors > uint_c(0) &&
+ Base_T::excessMassDistributionModel_.getModelType() !=
+ ExcessMassDistributionModel::ExcessMassModel::EvenlyAllInterfaceAndLiquid;
// compute mass to be distributed to neighboring cells
real_t deltaMass;
- if (preferInterface) { deltaMass = excessMass / real_c(interfaceNeighbors); }
- else { deltaMass = excessMass / real_c(EvenlyLiquidAndAllInterfacePreferInterfaceNeighbors); }
+ if (preferNewInterface) { deltaMass = excessMass / real_c(newInterfaceNeighbors); }
+ else
+ {
+ if (preferInterface) { deltaMass = excessMass / real_c(interfaceNeighbors); }
+ else { deltaMass = excessMass / real_c(liquidAndInterfaceNeighbors); }
+ }
// distribute the excess mass
for (auto pushDir = LatticeModel_T::Stencil::beginNoCenter(); pushDir != LatticeModel_T::Stencil::end(); ++pushDir)
@@ -571,20 +593,35 @@ void ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T,
continue;
}
- if (flagField->isFlagSet(neighborCell, Base_T::flagInfo_.interfaceFlag))
+ // distribute excess mass to newly converted interface cell
+ if (flagField->isMaskSet(neighborCell, Base_T::flagInfo_.convertedFlag | Base_T::flagInfo_.interfaceFlag))
{
// get density of neighboring interface cell
const real_t neighborDensity = pdfField->getDensity(neighborCell);
- // add excess mass directly to fill level for neighboring interface cells
+ // add excess mass directly to fill level for newly converted neighboring interface cells
fillField->get(neighborCell) += deltaMass / neighborDensity;
}
else
{
- if (flagField->isFlagSet(neighborCell, Base_T::flagInfo_.liquidFlag) && !preferInterface)
+ // distribute excess mass to old interface cell
+ if (flagField->isFlagSet(neighborCell, Base_T::flagInfo_.interfaceFlag) && !preferNewInterface)
{
- // add excess mass to excessMassField for neighboring liquid cells
- dstExcessMassField->get(neighborCell) += deltaMass;
+ // get density of neighboring interface cell
+ const real_t neighborDensity = pdfField->getDensity(neighborCell);
+
+ // add excess mass directly to fill level for neighboring interface cells
+ fillField->get(neighborCell) += deltaMass / neighborDensity;
+ }
+ else
+ {
+ // distribute excess mass to liquid cell
+ if (flagField->isFlagSet(neighborCell, Base_T::flagInfo_.liquidFlag) && !preferInterface &&
+ !preferNewInterface)
+ {
+ // add excess mass to excessMassField for neighboring liquid cells
+ dstExcessMassField->get(neighborCell) += deltaMass;
+ }
}
}
}
diff --git a/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h b/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h
index d654699409faadba97dc6f654035b8f0d22b0d17..7180ec57c255c30fb605f3cc9bd5fdb950b4ff9a 100644
--- a/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h
+++ b/src/lbm/free_surface/dynamics/SurfaceDynamicsHandler.h
@@ -103,7 +103,7 @@ class SurfaceDynamicsHandler
"generate the Smagorinsky model directly into the kernel.");
}
- if (excessMassDistributionModel_.isEvenlyLiquidAndAllInterfacePreferInterfaceType())
+ if (excessMassDistributionModel_.isEvenlyAllInterfaceFallbackLiquidType())
{
// add additional field for storing excess mass in liquid cells
excessMassFieldID_ =
@@ -389,14 +389,16 @@ class SurfaceDynamicsHandler
}
else
{
- if (excessMassDistributionModel_.isEvenlyLiquidAndAllInterfacePreferInterfaceType())
+ if (excessMassDistributionModel_.isEvenlyAllInterfaceFallbackLiquidType())
{
const ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T, ScalarField_T,
VectorField_T >
distributeMassSweep(excessMassDistributionModel_, fillFieldID_, flagFieldID_, pdfFieldID_, flagInfo,
excessMassFieldID_);
timeloop.add()
+ // perform this sweep also on "onlyLBM" blocks because liquid cells also exchange excess mass here
<< Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::fullFreeSurface)
+ << Sweep(distributeMassSweep, "Sweep: excess mass distribution", StateSweep::onlyLBM)
<< Sweep(emptySweep, "Empty sweep: distribute excess mass")
<< AfterFunction(CommunicationCorner_T(blockForest_, fillFieldID_, excessMassFieldID_),
"Communication: after excess mass distribution sweep")
diff --git a/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.cpp b/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.cpp
index 8c99cf8726ecd46fff29da30d9706edebb45ffcc..3bda9fd5762f118b3ab3e4ee97f0438bc309f2d7 100644
--- a/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.cpp
+++ b/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.cpp
@@ -275,7 +275,7 @@ void runSimulation(const ExcessMassDistributionModel& excessMassDistributionMode
}
else
{
- if (excessMassDistributionModel.isEvenlyLiquidAndAllInterfacePreferInterfaceType())
+ if (excessMassDistributionModel.isEvenlyAllInterfaceFallbackLiquidType())
{
const ExcessMassDistributionSweepInterfaceAndLiquid< LatticeModel_T, FlagField_T, ScalarField_T,
VectorField_T >
@@ -876,7 +876,7 @@ void runSimulation(const ExcessMassDistributionModel& excessMassDistributionMode
}
if (excessMassDistributionModel.getModelType() ==
- ExcessMassDistributionModel::ExcessMassModel::EvenlyLiquidAndAllInterface)
+ ExcessMassDistributionModel::ExcessMassModel::EvenlyAllInterfaceAndLiquid)
{
// left block
if (globalCell == Cell(cell_idx_c(0), cell_idx_c(0), cell_idx_c(0)))
@@ -990,7 +990,7 @@ void runSimulation(const ExcessMassDistributionModel& excessMassDistributionMode
}
if (excessMassDistributionModel.getModelType() ==
- ExcessMassDistributionModel::ExcessMassModel::EvenlyLiquidAndAllInterfacePreferInterface)
+ ExcessMassDistributionModel::ExcessMassModel::EvenlyAllInterfaceFallbackLiquid)
{
// left block
if (globalCell == Cell(cell_idx_c(0), cell_idx_c(0), cell_idx_c(0)))
@@ -1102,6 +1102,120 @@ void runSimulation(const ExcessMassDistributionModel& excessMassDistributionMode
WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
}
}
+
+ if (excessMassDistributionModel.getModelType() ==
+ ExcessMassDistributionModel::ExcessMassModel::EvenlyNewInterfaceFallbackLiquid)
+ {
+ // left block
+ if (globalCell == Cell(cell_idx_c(0), cell_idx_c(0), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(1), cell_idx_c(0), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.7), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(2), cell_idx_c(0), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.5), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(0), cell_idx_c(1), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(1), cell_idx_c(1), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(2), cell_idx_c(1), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.5), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(0), cell_idx_c(2), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(1), cell_idx_c(2), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.590909090909091), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(2), cell_idx_c(2), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.5), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ // right block
+ if (globalCell == Cell(cell_idx_c(3), cell_idx_c(0), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.595), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(4), cell_idx_c(0), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.615277777777778), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(5), cell_idx_c(0), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(3), cell_idx_c(1), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(4), cell_idx_c(1), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.605952380952381), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(5), cell_idx_c(1), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(3), cell_idx_c(2), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.5), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(4), cell_idx_c(2), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(0.573958333333333), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+
+ if (globalCell == Cell(cell_idx_c(5), cell_idx_c(2), cell_idx_c(0)))
+ {
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*fillFieldIt, real_c(1), real_c(1e-4));
+ WALBERLA_CHECK_FLOAT_EQUAL_EPSILON(*excessMassFieldIt, real_c(0), real_c(1e-4));
+ }
+ }
}) // WALBERLA_FOR_ALL_CELLS
}
@@ -1137,11 +1251,15 @@ int main(int argc, char** argv)
WALBERLA_LOG_INFO_ON_ROOT("Testing model " << model.getFullModelSpecification());
runSimulation(model);
- model = ExcessMassDistributionModel("EvenlyLiquidAndAllInterface");
+ model = ExcessMassDistributionModel("EvenlyAllInterfaceAndLiquid");
+ WALBERLA_LOG_INFO_ON_ROOT("Testing model " << model.getFullModelSpecification());
+ runSimulation(model);
+
+ model = ExcessMassDistributionModel("EvenlyAllInterfaceFallbackLiquid");
WALBERLA_LOG_INFO_ON_ROOT("Testing model " << model.getFullModelSpecification());
runSimulation(model);
- model = ExcessMassDistributionModel("EvenlyLiquidAndAllInterfacePreferInterface");
+ model = ExcessMassDistributionModel("EvenlyNewInterfaceFallbackLiquid");
WALBERLA_LOG_INFO_ON_ROOT("Testing model " << model.getFullModelSpecification());
runSimulation(model);
diff --git a/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.ods b/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.ods
index fa159c7b1b01c41933c16c886af01b6344a26bf9..8f0210d1c68a269eb7401a58947b64ebbc5bab05 100644
Binary files a/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.ods and b/tests/lbm/free_surface/dynamics/ExcessMassDistributionParallelTest.ods differ