From 473a6ecc8c98c8158fa83dabd931e1ec56c701de Mon Sep 17 00:00:00 2001
From: Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
Date: Wed, 22 Feb 2023 10:40:48 +0100
Subject: [PATCH] Free surface mass advection benchmark
---
apps/benchmarks/CMakeLists.txt | 1 +
.../FreeSurfaceAdvection/CMakeLists.txt | 13 +
.../FreeSurfaceAdvection/DeformationField.cpp | 351 ++++++++++++++++++
.../FreeSurfaceAdvection/DeformationField.prm | 96 +++++
.../FreeSurfaceAdvection/SingleVortex.cpp | 351 ++++++++++++++++++
.../FreeSurfaceAdvection/SingleVortex.prm | 96 +++++
.../FreeSurfaceAdvection}/ZalesakDisk.cpp | 185 +++++----
.../FreeSurfaceAdvection}/ZalesakDisk.prm | 32 +-
.../functionality/AdvectSweep.h | 152 ++++++++
.../functionality/AdvectionDynamicsHandler.h | 254 +++++++++++++
.../functionality/GeometricalErrorEvaluator.h | 139 +++++++
.../FreeSurface/BubblyPoiseuille.cpp | 19 +-
.../FreeSurface/BubblyPoiseuille.prm | 1 +
apps/showcases/FreeSurface/CMakeLists.txt | 8 +-
apps/showcases/FreeSurface/CapillaryWave.cpp | 15 +-
.../FreeSurface/DamBreakCylindrical.cpp | 20 +-
.../FreeSurface/DamBreakRectangular.cpp | 19 +-
apps/showcases/FreeSurface/DropImpact.cpp | 19 +-
apps/showcases/FreeSurface/DropImpact.prm | 1 +
apps/showcases/FreeSurface/DropWetting.cpp | 12 +-
apps/showcases/FreeSurface/GravityWave.cpp | 17 +-
.../FreeSurface/GravityWaveCodegen.cpp | 17 +-
apps/showcases/FreeSurface/MovingDrop.cpp | 18 +-
apps/showcases/FreeSurface/MovingDrop.prm | 20 +-
apps/showcases/FreeSurface/RisingBubble.cpp | 16 +-
apps/showcases/FreeSurface/TaylorBubble.cpp | 21 +-
src/lbm/free_surface/TotalMassComputer.h | 26 +-
src/lbm/free_surface/VtkWriter.h | 35 +-
.../dynamics/ExcessMassDistributionModel.h | 68 ++--
.../dynamics/ExcessMassDistributionSweep.h | 7 +-
.../ExcessMassDistributionSweep.impl.h | 89 +++--
.../dynamics/SurfaceDynamicsHandler.h | 6 +-
.../ExcessMassDistributionParallelTest.cpp | 128 ++++++-
.../ExcessMassDistributionParallelTest.ods | Bin 16650 -> 19371 bytes
34 files changed, 1994 insertions(+), 258 deletions(-)
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/CMakeLists.txt
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/DeformationField.cpp
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/DeformationField.prm
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/SingleVortex.cpp
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/SingleVortex.prm
rename apps/{showcases/FreeSurface => benchmarks/FreeSurfaceAdvection}/ZalesakDisk.cpp (65%)
rename apps/{showcases/FreeSurface => benchmarks/FreeSurfaceAdvection}/ZalesakDisk.prm (63%)
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectSweep.h
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/functionality/AdvectionDynamicsHandler.h
create mode 100644 apps/benchmarks/FreeSurfaceAdvection/functionality/GeometricalErrorEvaluator.h
diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index 2adab4750..3f5e6a95a 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 000000000..e14b5fe37
--- /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 000000000..fcbc4ca54
--- /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 000000000..359fd4a98
--- /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 000000000..d3c8e01ed
--- /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 000000000..359fd4a98
--- /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 edbaab294..f2177d559 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 b9af1280c..857f1510b 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 000000000..80db158e6
--- /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 000000000..1e48a1683
--- /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 000000000..b0a5d7bba
--- /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 e27edf893..e54117c15 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 e4fb01b3d..5b3c7728f 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 b0dd58608..34bed9ecd 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 64148c5b3..376a8d9cb 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 f48e23232..0cf759466 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 e7abd62e6..fba4e9579 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 896c5b098..61da9474b 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 8437db98e..1790dffb4 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 bdba59086..a55ad8808 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 acca70e45..702682bcf 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 dd603d62c..65d18594b 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 352a85c72..4fb4c49d3 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 dceed84e7..79718c8c0 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 8fc773876..6699fdf1e 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 a045cd6ac..ce2a45302 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 192ec8755..7b3fda3e8 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 a27103510..726e0ac73 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 2e4d0b798..a4b86c43c 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 ab9efe397..34ebe1b20 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 ba30c8445..12ffe0317 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 d65469940..7180ec57c 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 8c99cf872..3bda9fd57 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
GIT binary patch
delta 17171
zcmZs?1CS<B(=OOO?e1yYwr$(CZCmrUZQHipJ#E{zZDZ&A?*4c0ZfsUmWu6n6C#oW|
zp2$3v6%hbjnF<W2APoYF0t5sJ1mq$ur4|pT0Qz4(O%NY|^uNWn{}Gw~JC%3}Opr(d
zjPkz*HU9><d~m=2&vpK*k@Ej&u>Fr6%l{GxesTU+Xw?114g9~RvOq8a%CfNw3@BX>
z)I^Oduf&oIgHqi{5QBz@5-y~6s#bqt!B=W<%$2`+-!89(;3fE057OD~a}v+LotWmB
zU)8ufg?}rf_$SL*Zd}eDd>iZJdEeDI$;)f~QkQchuOs9WEGrbW(%$wco`m~?+?#|t
z`@nmlbn796DyFmys}RNo@LQkm6pAuNxR*VobapJ!wmv~=;gPtT*HwZ6_U_B8-+UO?
zzdMeG#~Vl%uB24&$GP3NWoSsDFfcA!gAB#gYMig7@H*2TY8aWAa%(To9{(bytjQ9l
zd{z3BOH?i3Hbz)U(eN%B8Vn9R3Xrf!8gW$_$;0r~pfI;cHac|z#;BZx6)DHmExDvH
zw-yq}G0nluB-Ly*;_av#=FQX<xPHUqS5C_}53cxrQ~q>)wRn6w|EQ4xq(Q1jAO4=g
z@Fil`F_dMa!@r;IO`v2Satvr^ZSq=d=)&zpe#V*4QTOj%Zpp{_n>=0%mqte<4jmS6
z>-@Zjw-}z*!)l2Eknzq7{WUtKLkur1z)3M0jVuMx7jeZaQA@+ZBDy?@oG#oE06+@T
z5Ri<kl&D3(KtM7eK>x=7aXlm?<bOj5=>H(-f3mZSr;VvIorkT>mX@vk1{=b6k3L~p
zZyj(V<(ee6Wq4pUy|;B>zGZzZ7(T5clcE%g1ee_Vg}72u)eD8(m5Xvf2~jThHm~~#
zd*Tcjr;r#$Qw;3HFKEtVL2%_k5JYemaZO4r7t>$7HgkARJqYN}T46ju2Ax~eEs>a}
zP_MUmWqjg!+Q3wgSWJwr60kL~iBjMUz=4V0;CN75+JY-Yqyq{|i3Q+7zwUYxkHPV2
z=p>pnW@l&VA{5bi&?#h4rK2Ry1i`Mb)u>#Ru6Jb!#2)g<wxY@ta`Du86B!}c0mfYZ
zM!X8g=DXRKi$7i7oz!kAa;7SV*vR^$4t$>LG_%nrC1(OxUIIioE``wBEm0%QLoRmo
z6UodaPD@dSVOfyDLG&|`Y3fsIcQrJfbvs7s0vVbYSY}GiI)9`rBF46Yzj!Apf7*ko
zv1}{zwf3?>O_`XO#9z^zw*?M05>FtQ0C;O1P@hl!e>0K|D5zz<>$|DD!xetB58Av-
zO$}@CDsm7K4bts#nlKnL=ySJc6(S;QT%$zNg{Np)T#KZ^q?L27XJ_x*jAX4hyBsxT
zgES*9`1R7P^@=<f$gl@FSs}?}yxgtZRXQA^05xvw%=Id1hJQ5Bc}Fcb92}{t16tRO
z5nML3$z`ulouW71D{J^`T<W|E>UikYFdJ7{gp3oWnzacgx{5g#U@VA^h_Vd&j&zD4
z&gV3Ti@AF%`|-Rv2z%>pZy+m(hr{YrNtIo<9Wp@Ri=$%r^PRx{9G{Dq2E`zb=obaN
zNzz%_4)t&(R)Ujs9@F-ZY1Lp+1EdPHP_7AwUY{(l|J+NVo79|uZgav09bFxHk8N!C
zc$>rDZZ+K;eEDpCydit?^XI+b7qP}loyL2`C!(%Px3X32S0-50w*>5+1eJefDI;M_
zkpFg941$^>a8w&alyOYVm&=XYQq7&0(kx=@zrGBv;)qlZAih{july5Y4*1kUKz~02
zo#Qn?*ioQrs1<&$(+>%t{A1XjI=rV&F~w~u?l{kvXIycN;HA)O#bx?#iLDquhsg}I
zwl9H8xGOdbJQ%0=q2xiEw|XnA`f8sf=6%LycBOn^SZ3`bMcHxvaKig=Mksv;Z40_l
zJX$Zoh`fWM<nwb{^a&lC10aov2@eV~ys7Y(5c>9~*>ERcBRk4;x@0v@!Uc%U{*eoJ
zt;!WwysztiqBo@A_01q2=-a~LWJK}fh5sGW<m=SO;8>*NVc~gcfz4rsvHr|$$rI25
zacQ5M)G|b28)>|Nc4po|e1nMp*O%kjM&WufwlTjNo0adzZ#g^16QFufq_8_n$t~4`
zbU5mghdsCZp_}9724Ci_S<06gN8OM(FENgI9JW-0wR{h`)1%E~ThHS)#y=%cHc`=h
zH<LZ}wE#v}^Qtt8v?i`l<&_aVxtm!{O%>EJ?rvTKa8uIZ3NWmU;*>Wm&-1KJmrWEP
zM?*H^_{klr`!mCk3&?)$^bXNVK4hTknI#z7MAxHSZv1jwk{4Bk^;ecuCpYA}{$^)7
z*CxIC{sa$fOg)nztW{BQ&ahO;7_O<iK%HYgLp`RSc|xaldPuEBq?(=hN6xc$EVQ1H
zL29q6B$<_W!MMWx*!Yx`b`^6!YCUP_rf!bre;>#1tZZ=s$Wp4X^Ww5VxOg{(&)sj^
z<B*LHaE(LQ!cZlmZv07?1;8`_if{5d0#}Q>Z{yIo;jm~jU0GUcA_^(}Cv^66?MUkO
zttx?OSQ6|f&p>W{j_b`kllV5e*}pZ`RSLaF@@ke2UVp8?s8l=u)ZGRIb8EMKJMZv6
z%`77Le^+M-TuTp>0OV6MO$P;Z1po*wk6#@Eb8dM$;9rdpr$6f(_xk^VZO2e7@*;2`
zAU5*<FKj1D{>B6BI2UlC{P_zI9a-Vl6|L;#s^+%%TB@%8HJ{Q7q^K<1Pa@a^_<}9V
zu7S~2&tJ%H_pKM0c84P;Kk?z*yZ;y;Fqa<XC^)Rh^PxnYCU7i^9KEbyVC*>Af!@jd
z(V0I(;jG(JBEKa-Ji&n`483Khl=B$KH;E_`=#>ck72^R21cJX0all#(H6dU2;O;=A
z-q_SM6hTuDZb$f2sPTK?R$AjtXsxPj@#4wo(}g`r>79$Ae_I0$o-S1ieQ4_G;`H@g
zVW7nx6T9W>C)#>E_;?z=bH!fIZpZY&IR^BRXeZ0oUVOX;dj4=2sqc2lI+&zA*HmUO
zK01G#fT0YazsieZt<NVSf|7x*derbrJ&7;4(_b#pPnR9_>8HaeA5fE2(Hxn$_kl2V
z9A2VNp<lTJH^l$q2Z#xSQJ>vf0UJ2<j4(Bsta~Cv8*tdYP9FnH6jZPfBXqn)8P+16
z<>iJY;JxP-2B6KIAV44aVc0-KDYF&eK(&M7C3XRP>q>Wo|40d`osS5<Pzuv0LHfDk
z6EhakZNB_uH6R<|F!^V7&AwFacv4KTzAt~{+^UM_kw(q23Ybw|Q{qGPG_oIbH3u%u
zCya-;QDT0IETW+Xcwzq1PVx#78Dl9RRzTi0z-MJ^`cnnd`&+F2N%wX=e)l~qgg4}Q
zp+y!z8RCbM5dd~O+AN53CVh(KKzdEum}Z2v-H++8Y}ru)^nhXt#lDkL(lp`-%jtqS
zvP4?!V>>=NW;tz$VvJp<(x0bdwe)lFHs}$9(nMmf6@?MyF@kvK7KR1Ggtal#sLzA=
zjf?6jc+$8}yvR0`L@MUiu11+!+5EeN6`~aI2%MS%fj4Kd8rqwTov)KL=qW;$finy7
z$%o+d0QuBxt~}}5roY9kh$QZ#?*xH!2a%4$O){QYzUR_Ugzb<Gw3COm`|BD)6%CrL
z_dc7=yLywwB5q<=YD=Bss8-yp)A#joYAF2H%b?RD+@s<R)1Mq0U~_7-yM(&`uEhjE
z6oc}TjGh0Rd2GY}GzP+f66$Wp#{LVuz~Vg-O>U$yU*<a}c28;^=P<ISh=;SNT_b)R
zk>IVeJ=yLRi1h`umK>+!=6sX^$~23VujuwLgVu?J4qu4I$t3d4DJ6bWS9Aedv~Gu`
z4mh+q5{<~O8fNOK>c_smO9S?!n$-hPmJ$;3Uj8#2!I{(b`aNqA91{xi(t%v^+gD*J
zUn{3u4}P(RHU2M*cU47)y0})Yw%7cUYNpNhmKV;V{~T2IceH)tomJVpbUDa|NM<K*
zS5t0sy>2@tr{0RY>d#tDlh1{t27b)Z&k-uXYw`B9FJt42d-FM|{h+<)%F-DiFp1t;
zDv{;ry9zh3VVw-Wq#?F%y?SQtV)9#RJ_fzhJf%9cjxIM;ZAMEYHswbp=h9{HKo_20
zTa41?U|Hj0leX9w!;J`=o4=_CQz|O*Be;a4JGo0=H^Wz%$CJ9s&0+U(YF=c&ZRK0T
z+f6SdJBI7bHQaRjZ8mc;Dn}FG<bAC5E~(D#PCQ`KJZaZ$oKc-`HeEVgv_02U-?3oK
zvxuG|SggZaCw1IUC11lwT9{}e>X(eyip#<wyFCT?s&Zk!v8_M<a4a`j5%7bjVA8zb
z;-R#BpOZ^v!fvM{x#BeFMO6Iae9AjPi_--9bLD35Z5oH;DJK*t0e}Xai}2q=0x^uC
zzF>}_HNSelp=A*|QLjzSLE`TBPy?vWey&^;cdZh<6eO9QH6E-8EU*p8dj*P)laMi;
z^d>bRq%9^vX2^i{{+2JbkjlB=6Pb1s6d#HQCK-1&VZ##d3NM>3N*`Ka%Jf<5MX<&G
zzNfe{+a5`|uxGr}#mobog{CFXtffZ@EpGaxI?Au7-SeEgU1zrd{{Ul#w_X+caS#DP
zk=JWdtWP?j)PX;K`Y6r&>D%)c&zasH&%diMrWm5q_F-GK0z2*^o>4Yqt9%|!<QT29
z%TN9MH*NPUIrmuv0|H_t{U7}QpUTh8FxDR$5gQ2TKPLuIqMiU&q8bu0;MbRejEX1-
z2nYrS20k?pIXMCu85t!#FD(;06D=JJ9Rmk5H3tU=C7U=0k0`qU5UrpZw~!Qvh>{R9
zg)k4h0H1)243M%6kft_}q#%!i2)C{wkbxYKku0aVC6K8Lkf$e*tu39kHM4^Qn}~>r
zn5?Rdg1VZtsIr`bwvv<vKtn@P)kIFuK~2k4O~*mjz(v!@O4H2EKwij1UD{Yx!B|Vh
zSVzlBN7+hG!$DupNl)28&&bF~$i_g$#z@=ANX^+;)5Sv1%T&q3Ow-Fs$K2f9)W*%q
z$<x)^#MRN-*U`k=$;!jSQ_nHT%r3yfE!5IC-q}0Y;ZK5FAdsy;AlAk&$ul6rH6+a=
zD&IH0#NQq$@eh!{t!9X$evrLUsGC`ck4wC-Rh+MLXaG=JFi=VwP(T2GY#~rmIZ#b4
zP<b~{TQN{;D^Ns$SF)dNYKU8Ylt)Q|Z)c}UcQ4RD4^ZzU(C{+Q(kRgK^1oCN=wut{
z{1E8k9_afUXlO_SaCm4I92^`JnHU(I5f+^k9UmJNpPCdMlAaKenwlDtS{|0&kd#%F
zlv9zBRq*eK$#04&Y|AJpPcCdoC~i+JZ%VG{N-u86s%Xxv>dK4tD@_e8%1Um?4$IHY
zD=WxsEr@O@NhvKYEvjm+ZLF+pXl<-6>1r%*X=zET?#l)=441TxmUYcFb`Cc6jCFJm
zRP`-&3`{f*Ej0|UHIHq!56^Xt&v#60_RQ|IPwjQj9rRYF^jD`3)@P5k77unc4EA)4
z_cY9RmrV?`FZEZij<)PiwI9rOjgF2EPcBT(O^z?DFH8>1Pmixn^{+3CFD)(3FRibv
z%`L5MZ7c)kw$~OmH#bLCPbao-7dG}6w@){=_jdNzruT2AckkBr&(}`xk5`5cHy1B<
zCr|d45BIh%_ZBY>Hy=({4-XIbkIoLxjt@^SPLB>Q&kipxF7_`U4z6D>uJ11IPEQ}+
zA1)8>ug_j@j^7^7A08g=pWhx|A0OX7-=FV3-yZ?r|NH<1^z-vmcCp<81cW{#DI%!i
zv3`{e&_OoO=)Ie^Zmn{+GqvC7z9?-)XMl?at2SV7UEJn`I&JE`PT#$91G&83`sTpu
z+GxDi?Z#FM_DiO{?ob7tj0>We4_jY%Qw5EJ<?k4bAegu!_rw3wZWRZcV32SU>%9Y|
zLICiWV_@vdbpI3Qc`{|v!!z^IHS^K)(1Q#dYZA6w??HnXIkNx%$Dl3Sb{}J`E3+Lr
zG#X#)kZv??+(5ST27uW^fH!f?cS7K(H9o|v)m&K4I~?7cjN=Wf8@rQQ-%`Un*}klN
z8Sg?zh5S-g0pWjFi-=e!RDUHqkd&sCCE80*$?Mw8ug@-gOcimbj+#~sx*NsSL2_tI
z_{Apy?M?Qm3Gz>kw&J>nqhY4ddKR(MT0R<|IS;7+V1(8g<yFNwut_Kr0;800c$ACu
zTm?yYwbe1*9B#Jjr5Rt)QWcl+mlUQcogEunMZB}X11Km?rQ#`zy_d{4QS5g;|BvVO
zdaRQH9|$(6A?%NUbb~}p>C0Ul2+~KNDThi%5?g!=bD@^J=cUsqUODb=tAgv7)EWmd
zC~Y_BUWQ`GE4ID_^Zw2Y9%xyv!7Q~D)-bGY<~(mg#;M2JzhaD7uvK!hID!(L0mzWD
z1SKPpJ-0mu1%ssgGPVUXf)y>Mo{+vLF)R#X<s+II$)tz$=8(EDB1BkQQn$uRiS{`-
zL70|;w5>p*PVHLa_{H!Fe2R-YL`(->VJ68jnZ%#~&KD$wAzTut)Ff$ELvm#O)ye<0
zzVnKh6GbQ-tZKuSxiLw5+>&0|4p>`xx7XjUHy**>m}@caG&la~>JyuOjpDcY$r!&x
zo+w^mC`tGjaT6P0vA!cG9<Pel?HkutMcXO=Iicn*4Q2l}KIz)aJsH&&o&#Sft?cu^
zuU*;$YJ2SJTXY1$PgNa=wj~cpZmrup_Ipai-~Woi<@a(vk^g!6>vsA10?4*^7>E!@
zLWYljQ)*^0XYC!)Rw)0rve)V9y40A@VI2uH+5H-Le{TKrKKb%<L_Xtl5zP;Hz_ind
zjW1F=RGQP5YJ%mBEeeidE!RAgQ91kO80AHzA$^3o_}JL_JlE&_xKz)1>w5uwR^iVp
zp|^vEyVuA|S)VR#eyL|n0Ho9gF8Y-Fng10TRqXXmEjORz{I&q+kJi8M`hYiMJB@Nk
z>2>Oov4A$8pw&$O{lh>?mH2lLv}6}5jd>+y6+-mHyU9ig9{%^WY5cCY?i*XcK`3{w
zZzA$?(|wl`ExfO<HQMq?@XUylb<=d7fvU?Z;`he{e<qsS=gEv4Adgp%jDOEn>vba{
zMgeuCzYR2j1ge%6nD+3}c~GgI<E$|mC2}oUrvHP)O&I^C^fiEi&ZMShL;RAKoIUsU
zd9Vw;)VhBnFn6N&qr!d)nQ34$@FXJmLR$<f*(+S8Jl~+W_;R4G#{P}C8k`d${3b)h
z1lyBk=vl&`-WIbA;MV0<TT+tbm0QrS19ec2G(9h*p_-5t*bQL2+ovBWnZH?B@}(u(
zr*D(G*jFRir2pe+WP{n|L1j~?GR4<VKGhSbcgnqe_+edCk!~LdZ{h&TdHbc!gBCkv
ze5b*qMi59q)kdhI7n1f&T*PxNZ^b1g6WTBjuz&XA2F#QJVl7U|dc(ivcCMS1Smo+Y
z6dDw3on2(pRAsjLTazOpVgs`4?CN?waIT4NttKFA<<tnW4(@R3>1=_Y6ftyCh@3-w
zkzxJ6zjW%Ec1IN*v<WYoR@xmdOnldz6sHlbbq-ZrxI(u5Z~h#*(oN#A`vX(?CK02l
zDB8oK_7hqGYGb=l&16}e!2r3Q(6}kbEy>;t$b3Y+=xt{Qzl-Zbps6Fm2)hE~ZZScn
zeTt**NM1T(94t}xqK}-odP76F%|2aaA~oxf3XnX1eQh^&`?{<*-`L}kx=em<ZRB((
zcE>z?55<k8RR<LtDOH?FHAFbDZ`mBc+#*pXpV{C69-tGA?SEN9ZK$C-PHTta{|&U~
zm>(vz@BSwIBSv@@OC&E&jZ<Q+%Am%XBEQqCt`}OP$T%mFig9GjIaCRvhqJ>qf?Ejy
z89vv$0hyTS{Yyt}2!~5v-&8|f^?*r7w5H@H<Jch|5b9mMHer$238j5)=ec~gz1pkY
z|FBsP@EZysk(;MtN|&l+myxWbVN0x9ms(#|NcyesncTP7(1#<b_;wG=n|QGPt2_v%
z_dSH(($g?okY;2olVaod5F0x-srqyKN{QgL|1y@(QS&L*7(0}+sH&J<JlE}f1%R_9
zgj%gHqU*SIvtpEI;zYXA7StAPf%mcS`TqVM5HFubXVK)RU5h#N;tw&Fk=2^MUXPYq
zGlt+mV-5I`o%XCO954+2{)RK>@WgbDElIu8P+3%-P^UR*5v4<XmATS18xf-okMctC
z71D#kk7!z7%!1{i0yxPJv=#X#2|MCgksj*V5x&U!0)-3?^~c|;3temA+xNxV6vTN0
zc(Rga&FeQ(GrfX|Q!;0f2UtOf1}bvu)ce_@OTsJ1Fcx;|*adU9`?dRbjYXjg)xGF{
z<@c(08HCT8-*cuE$M+)`DkiZEc_%9=MF@SYQOC&Kx%Qt#EIi-}{_mH-V4`|AB(Ay7
zI5gk_safPy{=PUq;I_hk>j5x8-o<{xF0|@?xv+GltvZS+J=M~u{%RFi+F#YPbR&-a
znVib7TpHdYs=x@Pg#Ed1IgtiWy=gU4ow;w-ZlubiyZRi{X*9I`pIiQJDKV_kxSFyh
z-ZHRb<L=jo0#QDL&qv~`zbrq>N-ZM|py;b&wR(dFwsS)F(6^{>0`$4_J7D|NCMn6<
zle^wH%d>5@vn#8FN6<pe<hh6n_F{Rz<?%~dpscsG@UWV2V2Q#-17~PABilv|*>Zhk
zZ9|a(b&l>iWbRfBH(U+sMspEEk`*==jN;Sy7DQWbTM4vw`j^qsYH4d`pj&Sn$5zTu
z@$wM>W5i-!x6xuD09cXROP5#7%t3{rus{1VXgb;`3<uQwQd5*Nl~*>R?rj?8@!Tc(
zQ1soPC8?d*b%nDcQ>tY+<@^{aG*15P;iuqBLDor~f=Z+bxNm=uLsEkbMsRCDDSjhm
zQ#*0BW+QOXg9}4}kkH~wty1T^ML!;m(j2{jH-Rx#Tkmj~0t}fGq=R80%g)B%r=?Xy
zT@iG)miO<q`0J?B67XWAP8asW#L-k8C+8*{qQHALr4-8}y3z0zFNMUX9s6I{gG?OD
zB8(&Y@_mlEK0K^cMW-Dj?%r?%jQS@g7~GuHgwO+X8@a2LJo)X2V3^fg`NBI%C(&2y
zT(3`W`I6hy0M;Y=;u^veK2z92^v~<{8<$F1P^BrQy#`w~{#D9*Nv%FqduyQkEC1k4
zU5d`6jam6HGB2ErVaO5k0ZsDexJM!WyQl=rq-@w;R1}>3_Wy`Ffk5Im7=CnFfqE4(
z{s6IEc`M{{Cf<k|x+e?1MJ&5UCosrJz)kdFGQ;8gelo+XhR}Shz<(f`NlD@xaL06s
zyvzH%71W|VFtBXK^HAAs%5-c<onK#^-R^H*Y92xBPA(nxrv!UOUtej?az)s9B5TNT
z`GJ8h9aqxhU4wdqQSF8u-{J7%7i~UojlK<0)~*#BX^R5mWx8X($SP2VEVebk(tk>n
z>!cxQP30-BD27NIX<qpr3iWN<FaoDf<8_{s*feRetvGrz)>ohk3}xgg_^&56$}2Bd
zx08b`Gy{l}{4Z|LHs@^QB>X8;sK!nt@B{U=n#4I|Yuomv)I-~hTJ^^Hg0k_$nA7V`
zm-Yf@iswD*y%uEg7(uk2J6T>pPC4>wxa``T<(e)Kq_nS#oaNWHM0dkDCklB0sQLBT
zHf=xJ5n~kUfrE?<6aCIe+0oDK$K7gZZIJ1b1e+wW2DTyeg%=gKZKZY;_g%I_k=MT!
z@btWUR9Y}4!yN3mmoi(_nhB$CT^@xSMC(K6M`xt&g$b&qokmG^uu=gBob(|i)J0#f
zLse2t1Icg)<qTrgaWQpN=Z%{U+<Q5wcG`W=?tObM+&BXc3lVa=S9I^;imf2u3|6rI
zq%t{c%r`Yz*>(i1i_^-Eo$X)HD@OsdKE+=RZvo#<AE*zijqP9eUq|Ohm1gn{6Y8#H
zjqL+AWYm+90}~n@8eM=WgZH34HrCDkm)$^3Zn)0Ao5&ulw=k-FfAxQGDqI~DYlhv`
zXHdrM0QZMvsNDMek??BA$yw)9g6oY5yJsWj|1J*zLbQ25vY$MJC&^x$xr<GHn^)*5
z7fx+$9mZap7_}>T{3YRFZ5`LkmW>MQ+Wmt|L(Wdo1T=sFDrvc-6^qih+`ySJa_+wD
zr(!#v8`j{BQ*h{kqIm!>>Uaz=x8SBo!2-YK*HfrDqF4!vuWhugcJqnVUVc#Op#$1q
zX_>V%>#WyL4XBoKLf7ZAt<G(Dn&m)B!kp||6d>-Ezb`0F4=b%mS-H_DDrHVw1Cm;C
z)bTN<y(0ijgPEu#jn9fl_`nISr}hKzOy%lWonn3Qg(fzKBQ)kqaSO6fNnRPo0muo~
zzB|3TDxE}M(n1M0cm`*$AirC3nmg*JgG-;=%+g`mmUkS1@?54oyG{MBvphKb10C~C
ztZ6BW$<r-y_3e@~YU^8uqyf&Qr{7O_RHK$17jyvG7xA9L9_H^0;EW#1U#H$ifM8kL
zN=7rC1WH6LV&zGl)g{~zTKREi*xI(R=<Jn)My`6{c|w{kt{5npr&<awZEK|L^py~j
zT$b{q*L}&B!wQy__l_u4Vp9pNuVXUi#}jG$ki;ep=``0XA3}sm8#lz^b7ba#04w%>
zszgBYF!^?$o1TANlIJ+B5xkma6HPId@F*NbEj8T^gfELJ6gTN4@0<c?(EaWjt)@6I
zUCFZOvC4r>Q+da3`(}SbVL0vhmuYR4w?5`gJep-fQ@Pf|@)8xoGHqjfAstQft=G8$
zdYX>c*C4G8?qjp97ff0fL&@24?6RuxMkl~cHu4OOtN%|Hl+3kTiFMUcOUWDX=UeD3
zqRM+G!eVQGF$Wxj+HA6e76-fSH}I!N;+!>TiPiC*!HiHnUStH4gI;0EYy4~Zk#vrR
z8-8UxV-hNwen^Ec`-dv)bMi<07X=Q@a!fabudeM7tsBUF_V4QoBa^3LnIYNPvU9+@
zoK_q4UhjEj{9A2Pc#O@LN>fXv5y@P@4rNtmmaed0Kgqe?PS+(<MKadHI<=xNd5f@t
zkE=#m6<jis+MYURyEBE7x9WZt#J}`-Rl}i|R>`HG3O1cD!)D+*Ua94|?7rx5-)W7j
z8uUZnEuzZo<VB8x0~{cq4ISS!EDCs8Wp~+`xTgBk2V<=@RK&jh?eeY2uIObDWk8lh
z&h}_Piqh8O{P?ac9=flK{?3qc+KL}lxbvz32X%4;mw)dORyPY<1`qTBiRaPnVK(e=
z9-#{rf@YvY)lnJ0ePSabbX{3pkO_dJwp9>Uk2$b{GN#Yt7bqzJ$0*IPI|n@6nYD-J
z1&KAd!2E&0sn~gbEqju04RT`=-+#c<jVK_%c!i*e!5y4}<>ja?q?+%@@tmk`uaHMZ
z0Ih;(4yu0aRETPVz~Mw$I*``}hc>U3<)_M7RuHeuf)7HS++XdaJ+OPQ_yKiU?c~Sw
z@0rNsS>u^=Tw42s`CdJRYzF8Y3*uu1Zu+v%I(eG|`<G#TLP>v0U9l=xX3r{%Z@WbT
zu&cw>O3m>-c$-;VWw!}E%XJ93Nzk{J>=Ra`spgj~i|c!KZb#qn-yKetR3$1#4x)m)
z%XGZ<VvtoCbt|DxVA#~7_SgM_{D9*hSS6<$NFt>{L4zFA!jG)#$pOzM;ob<#@3Q!?
z%T)YM;ZUoy%q>{zq(d>FF!R#pq5J>(6<&z?ov2McXWwyj05GS2q~@!yo#t2`6;-4y
zXJ+x%DN3GW6_Z+GH)(|rBf()pzh!S6-wF!M9SyhkRf65!2lo5DNPhh*hkU<h#6)QB
z%%8Yj!|Fb<mXGa=cL0Z^)$?sA_ffG$-Su(|!}ROm%~kMoHJh5Im_^cg@ijrv9nW@r
zwv1vq;G#4YB5tt9C(#&Z-FusP$GldA)G-BO(?^nq*`ixG-1%kR+ZFWG9i~hrrLmF+
z?^YZ8jg=pK8m8?I<V8yo{M2M|uH$<^FNsS69Qp2;$xaPKq5!+nF)L%k`#txAX=xZE
z`xBC`9jSQYNNFhe7X6-dZDjV2e6=+Z7iL@Uk6_j2>VL`7$zK*JGu6sM>Wqr}`%d33
zaXz1<S?o&HVx=Xlyy}VRMFBuP@VSaET}V6UqERi_fQDz}_9R0~<q)f(4b$t$2M9I}
zqV(N?qoR&FKH#Bmrduu({(IJhbwGCN^NVVttJVDTsK2fXu}Q3U!Mce{F4D|@{#=&X
zuY0U!fm_5qv-i-(n8pp)?%Jj~9jjT+Pw+UAPB~k(^ybDyKZ0j8&qVcx4F#*W+t2a#
z84`=}jMb~l>M2a2c8Z)tD_Sp+X2g*~>17QwdmBOx7Z5D&1E<upZOoM1>Lg#&ip!~`
zs;wpR_y;*b>tHgdK<fpVDJq-VI!PQj?07IxUfs${`eEO20CF`vCHh93ty@E{E&pV5
z<J&k|Y;^X^%0-KY;dj~6SwWgex?b*EfJ!WB&}zvv0hfK1W4!SKt3zejn_fqCTS<*v
zL-_*x7(fK$!Xe70IA!o_Q!_iLF(2X5NZ}^}h&4H@el82z;Y;;(N)NHCr?aldVIdog
zlCZ<w<?nCm-%*Nq)#}5<Q8M8v$#HH-QCn{w<drxvyXD`1`F1We8@NKlR6Goe9`Qw2
zNZ|y1Gd7Cwd2@kpy~zDb<z0t3P#~ZP#GIJV)d4vMKkDOHLB;FmWb^28t=D3+TqzLr
z81l88IMs?ZVu}$eu}ui|YV1%YFK=>cEr@@%R~mITp^Ir}DVS%y^9<DMg=rVQQ%V%^
zD$W%JZzLa6r}^sO&n-c0QFo6wy}W+9yB&H;(Db;WU=)QiRBNP~PN`irEG8ooTs=P?
zqX0xXQ(9bE?HY7V#$EN71#{6Wc7DY5=6oQ3-<+{@l01V($YLnXT8=A7>liF$p<od6
zd4;G6inNr3p^PbPIP01oJ17VC225C#*6)zq<sxFmg|59ZYi6t!|FSOiF<3_+Tm)q7
zixX3Hi))>3*N@I_HG0@offEwS`J~*9G65v?#E9p2-%EN?N7USs+R|i@A|;uiXipo{
zObU5ANoF*i?r&i#I0x`*(*i7%hG3~Lu-d|a<lJt;2tKWpVvSgSGmn%!M37KS>1Y(2
z>m=YrG?^+Zi8oxHSu6Z#ToFcGl2K?26p31EWbVyz(M}6w@E%Egu~ZnrT7Rzc3;|Lm
zy|VcX_zbTtQ#tC3KD=xm+*K`a=WlvRTSofds{5fJ5zwr=Fk3E-s<*J1tZLO-<!B^E
z6r;$sDk}%U1c2&?)i>x7T<5f=-dI+{k{Il4SXoREfR%*ON2rZp@vQNw3YypMyt<_-
zb*dl2ER?d+rvK_^pWJ<WG2&**kOMT7EXEiFYhQ75C`qMLW@(}0#_}{<$CNLP3QMVD
zmvObCUF1j&c;@f3%rvfh|9w1J%Ds|x<LDcDrs`bXHWl})V<L~#UXC-qQ=wU*)|+Dn
zOOqNn=3v`)_L^I4)yR#btt%l(!Dx&kWRV|6I9mmKZo^9y)j9%EdH*PdRs*QB-^-Sy
zZhbXQ4sM;<XaD8xU?BaJ`FhOx7&!*<>=+f1-#_QYBL-ijt>X}jAx<lKgjs64_3(Sp
z&)T0nH!AIQVon?hU*PT?M%(9&7K>}*@5`l>g)Ir<E&0$W<~@JJIi{ts*wB|Z9pA$a
zjS>6<Y&ljMoXr@QNMEk8tOO7pd#5-LcR#1J8k8vJrY_?hHKVJ7{U-eV=Z6D!sa*kA
zq@5a8&n~I9Av+O%Fm|YKJm*1$r)w4U04=8CMd}h1XUrQ8(!7)M1btA%*qVnzo`09R
zGGNWP->5lkchoM(RZxzD-c%<ka5fzXeT){8V_3G={@Gv_+>)xwIfjWBQ!(hVn;mt$
zf<Bmd`3G3^*32nc>{8%AU2c-1aw0WC1_A%?2DnCk5(i+wf4V6Dqw_xm^G`6gw{tPI
zbNNp%+!gOv>^?`s-@ku<)fX;IPu|7f%y#mYcBkhobcM7pOjT_ss#TY=V1$v8=pelt
z|Gb-%`;f#5aB_r0m5^cyB`Y5_RF12fo~-=srpa1Hz6aTGw%<@{-LWHdzm&CccLQvE
ze|?<9tbK2U{=6?`%3U#bq<b7J{c1NfzicP`-sxYpb7*{2^!ux6xHvlQZqGwq7l^&0
zmA6ZVv*RIsXy;%L=V({(F}yM&m>1M~<%vOO@6~;Ba&p}`|MaPO^YPdFuvIQhuYVid
z>sX#l4y|ePhPqqR-Y8h@WK3-4-wrK<-tit0vM|<f@b1$US_^{{7TTw8ca~4KoWkIx
z`YnQ{O~!R5hSpDnbr?j4gs)vo@S&kR4Y)|s>wa+vQtb5CW#%T-Uy;GtA}^`jzz9(q
zuh`qoIk71@_;z{WJbCQIb(D4k>Ja{Kq`i`tV9>M_O#ZU_Jsd3`!~l#qeRe>eGQ3w}
zYVVm$H%`8XvIX-t%E6enVN9zDoq7*4A}8icTp*;>Y+}c}8<$6pyNHt&{dW&G**GYl
zOp;_mX88%d-ty=McD1$uz4YM*=$aFQBu4&$Y>7uf(CAAFdndTk!?x3KXZ+zzONZHc
z^`>MM8B+$i1$#N#GkO0UUp(NQ_gCZ5xY#u4aPo>K=m!nrWz#4WpU&mn7Se0}6J#zP
zq65-~E-hX>JQ~At5(M+_#y1}q0p1U`Yk$z!R|*b9#fbI;5VzS#KQ@dQbpyXDZ`L^>
zsbGyg9&wnoZp3d!FKpPv0lf9X4`LpSk7uMRD&x~dgJ-n%<FHA_M-jkdpITr11_URh
z6aBZR`59R$t!`)Y^+Zr7K`XzFPhi?^UTZhQ5TB5_RkWgcz_8ry-in8ZEq&|$O%bx)
zhA?*Ql0;aweM%(Mz_Kk!M_u;H#<PA^G5Ti%1iHRWtCB48>-6t-4}{%444w_axB!jy
zJq+ouvHL3vk0-gf<0-&Z-4hhAW)15A$6M1CuuzyctWXn2Xh*Hysi^MU`|Hk2wLQLb
z`36K)8kfM2DKK+kLn=Gf%1X8BsiCbwEVm-}#vWF8sRx$#{Ii+$J27kURZpW<+ZC+Y
zy^(qN6zIJ$HduAJ`@<>mh4=~Lo@5v{N$Yk~8#thaT!c|lR}`?qx)F^bG3y+~(N4Ah
zuDgIUJ9~WMnsvi_gBVgvQXT;yLO&VJnG$w#9%Hzhibd{%Rv6I^dG2@OlDwSNC&cP-
zes4n9swlreMa0uu!|iEr>qWzV%PwXW&c--E&d{PJ-#nn|Q@k{lximP_sI~T5u13@c
zswLrqu@BI8`UJ?mN-)|fZWjr~G`^WyzM0%%`f^G!2y!KXX!^R)jO4Uu_iDclGH0&S
zXAa+0fNO3f(|(UR1<=~d7ZF}|51zJO1xK%BLtTYZHk$P`XWdK{Z_z?U4r&*jC3KbR
zIHwx)<S^FG0g>7Py(O$D=C{&xr((^@qDGS!59ScPPXe%h(%uGf$y0^~tYdS@<{~2`
znyrY0@jq}J6IFt&d}6qdiKBG%G`=U@xZyI8h{cZL^BTa;-;iyF@VXiPZFu}yQAjnU
zzQysLe=o=HUN@gUc%n7BF0OG~{e~Oe!FZ=(ax>6iG-Fo{4I0N8OQp@1$>S0!Ql~1I
z#6v;ivID?_<4hLVP=8V+@)PzUBdNSU$ukl5GJ{$;Fo0l2oT{AY7gQRm^?;`Dqf_@e
zuq;<Lvzc;((u8uc6D5o{k;7t7E2ZD@|8aeRqMl4XPZ|H+q}vkqnl5G<`r5_<lS*2k
zwQq<_omS`+bm&Dl-n#U`K8-ekr$Xzr$C8s-w+v8}Uzuk=3C>+cP-Hc`f2sWnao0|V
zQ7VBWfPbi+EH=EGt|+Fqf+9_^Ad*`TXr5YSC(TNSY+nZ=lnqNYO)Bg|6iN<NQ8#7W
z2B#BRvO+eTLL}f<u3M+@?bb{=hPwSQ%!Au8B1d!j3o=4o(J@v2IJ_vr5AmG)YF_^L
zw>w}ZXKtAO7Cj{~;Wrx2J?@6+9mU!!Tsnd^u(+&PuCfE71jmT6A=^z*t#s0@u+cah
zI|o6(b}%S6<@c!<Lo1od#UQY3x%W_Vf(hi>^|<zs@1DlBjp@`aJVgx>?g1m#ds{Bx
zr^V6=WTF1oMt{JUa}V`gv>!Xa&!?B4kQ;!U^Yo~NxAi5yh2M9rz!t34lCf~x2i{>r
zN^I{FE=j_B!HnkN3%~UW=MyJ30@h87M$*p~29u!PYRsPKo3}JrBFUmYOZzTQHiEPN
zIo`O?GUT(n6gkih#(~2OzQlsQ)y?;7ppe~j%XU`deKf-WA1fIDPxOVy*q^LqWLdzl
zoS%?VPbcNwRx@XFeS}Mu_#x<qBi=w9l_~y|!@pX#I}SJ7L|+1JvGGb~IE5d|jdb!Z
zcJ2d<VHti0(tjSD(-K&M(;3=c_iqj`*3xOMb+U0MqILvL;5<kxrTuzTt6H^SpqvfN
z>|slFis#BETVWc}1h_QOd>1GYgbV;8XfYK!7HY?Px{Z2Gv|gOR?gj;S9ax>4M55H_
z^c@?>m2ED;`(Js`I^Js+-NJ6uGNU;SUHuhB@K`@>+sfvx-uP#Ngiu4P1{;PA)~{Xq
z6YEX{OYhHX=Rq`}>AlT)|K?1O$i=?;SG}OcJmn&6`>437PuP%EodYG<84iF_!R?T8
zjS|c}!^M_a)_?pJrf0vEY}L16tD#g|*0<jN)#+m&-<cj#JLW0kQ%KiQ(<-mQZm$iV
z&N#!NXf4RvO|DgqxvaP_P-_kLd?tp~+LsK%A1H*y_D5X1JmRi(fk(x~6v;44b}%(1
zGG~Tdgo@@GloYSHY&j|Jff=9{<GS%&!rlsSxoHk9Y_c(8m|Kwpb=-)t(__l^<~3_-
zJi|vT#E-0}2@kx|z?05xGp5X(?wr8jI;q4H3u2M{`W>mD*m3jAX(0bWxw15IUL@@V
zxzA#9^fQg7+ix4Q|B1K~Z-P}*^z>xY<c=!AH&M0fqLNLsNGCGU5&+;Mt%_QHgFc0h
zN|_Fv!*XCYL-`Yx(w8}3%A&P72k4pGk156}RE?Uhh+clH2%W>QOL*T~FzLIYFd95@
zp|)o<!d7MzksmGRQQ~psp_ch^+PF4%#vQOL4^e^3%A@3VN|aSr4^&)Ql7cy06kP4;
zWFBqNbd3|!6zqjs8UmP$CcCg2nn@5e5ERodY41Z5<*2@3QUJxvs-R1#rGgo*-9h5H
z)Gn4P^vof13uKe&SKb$=qaQt+-L)%LKNqKA3w_|}@~>3-Hm1(w{I0u9W<?)YN6Mj>
zR#cXWlsdEOUu-kkg85q*Cdf(w{sh(TkK1#P=vOfbr%cmWsS5y|qiLLQjfTIh6EY2j
z5P;=vQ{$FF^`~lPp9vF47wC{d^bf%Ljp$zyOp$@=rVdP{7aYE?dWJ%;5(-=61|`Bs
zEp;>FOr!uF$0RJ0XEoP@j_04?=d4uWlfI?uA5<=b;xOIw&z_3xTuda-g+R9>OhC|;
z{xgQ?*HRwjB4FiIq*@6?K|unf$5M}qNa@rQ);tR9H6G#91txVGTV;e{Z2|Fzys;1!
zAho&6ZFXyZIVcUVjemJUA*{pSW^%lIc;mSZn*5c#0JXl;??oK07$8<=;+Fk)quzo3
z5LW5dqvomW_b@dB#o?Z%z4&ImwEEYlOXMvs{hB*X7a%w;`R24`yLt9!jr*x>F1rT3
zrH-+Acn|s`MR{Ml#qNg*24FemMo!m!nI{SqsP4}5!n^W~EIZD?;nJg)az!e(tNX+I
zCXPxFL<@mzH<S*Q!82d{a>rkR0ngH1-azp+p#p;{`}i6!?CHCbf73P2=2TWHMvc0(
z6thCi42b^O^D4fefkuB$ayVn6O{SWcd5NOd6`AwIb+~a|b-~~2|8f!^bFG<|Tg!wC
zp0Dw-kQyHY78eGXntlK3BBp8&tac5JZeOPNUd2my)rr+Yj7LBWuOLQpr#D6cC8mai
z3q+c;4NSReG3$9+zt=fwZib#M3O-S_w)inO1<3z#Cf6TrOfOl+Si+B>V7~%!UmVZV
zEUx%<A*|cUm)q`I@|_{0N(p?$DOOyJnR0%@-;rDNEiCT!5#v)(Sz49glcAH9<@~6J
zhdqnr6mBcZ=W}kx?a4p#QER+5x7?pH>9ThlB&GUmDfF!ZwC=4#JYOqb*gM8QR^P4U
zH^7c|V)VLPWkcw_?w9AU(W5X^`ef?oEz`Q2Ki}a+JKeKR4cRy5+*#h0b-_ds-OarM
zgUIegt(~G_&=^UYqibrH8#jf3uaDf^KR>Uy8;!Annb+AQyq^L5QqXTzyj{MzzuHF#
z3Na3I{GuQsMcIU|$9s`1e`yL5Hl!g6%>i`63+lJ(8#4?}y;P>~xBgB9LKMh_Rm4AG
zf0_-xgXv?8VJ(v#{z_)^SM5X=$d?YHvvg~{p6dS^$ud<WQWt{qDW}c-`JE~O+%k_<
zS4XBy=t~hP=~SOs^jzojt=O5__J!<}y|RN7Do6lX2dYq~+oLSAzx(&?HKgs;{u;1t
zQUwCJX#wkwKM=jmp!d7({Gi_&N<tP>lWlJBn~G<gBTx&;5syb?&g^K9H)X4oz2SRy
zVf%SsT7{c{l=%%}3qF0As=5o=*$ds#EA=?Z=M%@51!}>)Nrt+kB7OWA%#<s7$mpR~
z2;FSOM3=FZ-s060zLE*w>!%O0X%etZN)px|RzLxNv!58bEq?CIjb#bXTp?Q_pvxSE
zfs4n)Rb=LpmkI7M<N8}S!>(^a6)PJs=AAp?^X#sdwM)F!KCGGppCn1oMoWU+*35W^
znh|ACJOXuuAw4N%Lc3JBoIp2IP8m@{SyCZA$E|@_T9>9VktKfNz^q?hMF@apW`zc3
z+Zeki%;9~>TyvhexE!NBz!+{cZ_DAOlosbT1#c+FOu*_MC{7+UDx}~z{fTU`Ut4M(
z?3u&C)jfie;iMLmv2QY2giDP|5=W4fx(Cx^DyfsAc8{&fL(vwzDnl(g)QUH$yO1(S
z*XYV3M!dIB+1qZ`nrt>R{s2&{RTNwhHh!PH4!*kI)$2k`@Qnkj_SE@R7-BrI@Elr7
z@JfP*53&19<@6gBVsXq)79@n(f)MQpq}sodM?wS(-NRb_lz5kplRv<M+BJnBzlB7O
zv5TC6%f+e|N#js(lPl0Bam>Psrnb;iWLnao#pFvN>kOb+G}L?hHv)JptJXI~wahM%
z<}t%Tnicl*=*s((hHvMSLPjCOqIy+LK;fncOg)rSNsaR|#6IDc>nLf+%;FL;aIl%h
z-)K5If*fcu{r~-xR0%babGI6FS71uN9>Pj=jt}wWiQDV4_aU4pikZInTbS_!L^fu-
zug3{GGJGuIEgllPQ4}zypc~*P>rj&e36ps!{6}&+{SY9s0pC7bNF+AX0*qK#8aW3L
zok1NsRYfXXkstXl6OH<AFnZ+l{>i{cG_MzipW&TGc4&I;OC9-``stQ*v*o`yh;P
zwW<5W!Xr!3s1}W?egz!5tBT|o+MY{~Fd|Xbv@{2AOo=m~P#aDLUREUN!+=9PcBm`G
z&s`6yOREf2kKL{KlBx}bA)0jmFblg>Y;f$M;^0VV>8_l~-2=*ob?MI0B^!In%b=g4
z6>{KUB$+3cZzM^PN%gETp|MohEV)ywY`ohBF<bed2Cc%89Xz4O*Z<+#k_J%=J^x7>
zwtB^bEExe%H<}&y^AmWvhVvq6N={{I5*p~{6$;{yZD|xLOx!?z*h$0txIaKU#!8vc
ziG!Oh%{K~f4@0V9vbOr=zXEAl74ae$cat}rx2>E{hC|%OpFC_-|2kV^DGx^5=n7T-
znge~$%V=<O0L;Dldv(vPm(6NJ)Qu_J2jU6Lwnq|>x(mT2T1T`o=q?~{lIbzVc4V_1
zc3JnkVIT?E+me9VS6KgVVC<8GA!GeO*9d*tC5>rQGaaFhXACt{9Mh;k1WQ_7Mzl=e
zkY?SBN^Eb%-xMD60SOZ_)VTD+viaR!6gJyNNcrD=>_HCP<Yn`^@zKfxY%ztDgbjK`
zluiWzZBZqXWR&vDRq|*X$~ifEN<s%>SsTjsx*`eV7~gei-A!BHz5>k1MVKyrvgn~f
z-XiTtETwa>B5(LgzLj{VJ>Ar}G&7nAQwl5OCJv1A`EKcWZG`q<;*bO-h>46Esjnw)
z4T#eEG>maF+vt3w5-~y@X(_=hWY=yoS<r8Qb<~?<ce3b##R>#=;{c-YJ`H(HNj_wE
z<Sh4wUH9s&w`nP-<+1FG+)XZAw{hDXmR!&_J=r7NcC3QKPhVc{D@bl!?l2YO`tSp=
zJl0WT<?G7@htX*UoX2|Q7{1*EV?@e6d2T)!pN9dwbkOWu413tWRt`d!Y52f-k6c@T
z3R#49@Xlr#Ta^OefR$oN7`X4!Z1EHKCbBF+I05I0Sf!w#libp*X~?K+%q2N8_Ykr$
zENuLM)durq0)^$t%4_cp#s2w;OC`tXSqJDOV4vg<<gk6r!w_sdfAjDuyE{|UPh$i1
z`aoVH@9+0B-rt~hMdp`=9{8DAv0P(-QH-yvU~eZw5jurpt(awV#>u86we`&)5$PEZ
zeAv({%=kqpG}BgZZMHsE_kYpmE;%SO605-sV;iYJ;+toG3Q>Espdqu2+N4_1VPvOo
z(6ONLRLx0|UM(##)G9hy82_GB(oy6!Ow2EH@AxRTygb8=BdXpk1plkg<}wDb^9z49
zc^GMppa;a72zyYPjp9afBa|(9nGM!tDoJ&W;D$*fvbZj7b&>2UFYmwZ^Q@tq!K>lN
zdbP?=N;EqU2SxHQw}6g}CfgGbK(__Lj@KOb{*-V+?D%fOykCSTWqj`_fIED})wxX_
zrlMWUewn-xK9H<zg_UGS1#Bs(_~!Md0{sKE;Y6p4AZS2PuGXc5Dw|nSLUqTW<Fv5+
zq!2tWrM}+@&+{OUnKUKm($S!^Z12P*3mW!NzkdsQI8zI37|%GF>^+m1;4GHrYxSlG
z^x}h3A72J=NYn}BVu@06b!w`qfw3f~VF?oENtMwGc=+)9<A*pqfQ*y|F62-tAD=D;
z`eLpNsBZ=WCSr6|`KfUY-=d-6PbeM!bKQ6po(XMY8HNmB7VPy%(VWIB<48&m1Xf%k
z4>dOIC=~MtGj!F;yTfdk)j_ucI_8i{pa#|B=^Ss<apl-1n1R}Tt`yL?0G-3?2Od`B
z5`hPWlFXcUU`0nVz}a^IhWSGnnTnGN6DC%gCSSC=9Kj*|R}HV2Qr_}o;}J*f@~e0<
zSSdIq2df-<7E)WhW6|y+0a&S7!ETZcM*6=ivbLP6L&cSkT)RSFc@4>T9-9|@<#~Ye
z)n$OPt}Ko|P>rLsVg2Y0=8B4g)P#+$5E59a3Rvkq9c38+g^EK}b?r2>$w|&QJV3eU
z!d9e@@&>H*UkIq_Qg@W3RnUpv7`SHl8yP8Q{xctvi55ZtD@}iPK#xm>LCwN)%H{G-
z5u=ig4P%yD)fg62FCbGd-9WL9=Tg}U!jm>v)3f!m;IvF`AqlMfhNMBJUYi_O(rT$V
z^->hg&<2QlSrAxGEU^0K{fk|g^nNN?!>~jRk^3Om*83D5PF(}6AM&>9SLJ}fn>1wu
zw_ujr{B+cEaB|Ts8yW`k$wiuji3%I4N}XsUIERD?V`c~r!xB3}{bdr0@(^e~2pKu4
ziL^;zy|i-8;r=_jk^|cDc`Ul!b#=AghNjxncEEB2OJY)Sdg;VG<mM=kQi@qZLz7?u
zO|%XMtcWHSc}ZdgjVxa)jVzA^w9LIfwq&}=;kS!R4o1M#t{IO7)i&B#+4zJ4I9SEo
zyhg!z<c&^C+5LX;c;rK<OuB5vT!E1ok4eNB&Q#^Vc)f=Cvf-I;buS~hw?+X{r1Qnn
zDL{%QU6mB=F$_j5)cI(YEe4Oz#bAR%Iu6|69GsP|k0QO|16bAbNyRGPtD-^FMUFLj
zye2~`m1S+Urx!KVl15SX|C(yjKc#2!#imzJE^%mHz0${zYr>5iv3@fIuZD<+`lc35
za!f*m%8I?SvI5E(gF5T3I+|_Di;Mgz!Ni*8@#)mE*rN5G?0a=;-}Ky^CoOQMRG#g}
zB2&(19Id;%`|mv|{hs@%H=O%P+YG)p%{w2q|7Qnn)6DoH-ug<Mf#KL5M&Py`MkWyk
z#1RdX=i5L}VMw%<n*82I19E%<h+}FC=Qu$`PuodNK4_~2;WPp{oOU`0nP5A(j2leU
z1;P-v_k=JK?9F&Va~bdr|C1NnYs!N}koS^;Hu@t##^l%bl93=a$oq|W5DI`S7a+Gd
zwWK67FCBb10J=9p`;idak&v|pO^$U?R{)tl<KDdhL*O9<IxGwfqR2`ZX6jD%RS}=O
b+kpqH5IH<S6EO(Dt3Ua<gDjhYJ4grs^<6P|
delta 14467
zcmaKT1yo+U^6!ThcQ0OCN^ytcE=51w-QC%^m*VbH+}*9XYmwsa?i60nIsg0ay>Gqu
zX05$4nM|@JlV6gZnX+VvhyVx_MHwh)ED#6|1fp=2Qi(-Tg#J67AdL;c_!k?0B!PhR
zUuqL9f0xA{K#;~$Kw$losr-}altlS2&A&C^q>%nwSP8-Lulj8NFr<*&e=~JDkh~E8
zh%%t?fq1!Z>7UR-&VOSMwW)t1q4E<0mxPvvR%n)Jlr)7c)|G1wMqAevymju;w^=w9
z??AEgosYE}6_ZRz;I{|V2$wNR=i$NRr0fa1Xx-F5V<jErJ}MsXoAvlPjz4{%!6oi#
z@K2_b;q7brm4~}8kF{-pmq!*^e`rrIh=2CU4$yhuiGX)Dm50}_HTPc3#VWY4ls}*_
z`il=_Lt~B|L*_u*b{*Ey&k9>9RaJM-w!93TW@Wn_QHKzrUhw)Bh4jUUWd+j(SO1Ay
zgU|(`-U02C9^&X373S!5TX0daMiZ~6&CxJ3>lBYdkgwNNfkRqIjo|LbkB9EHG*VYl
zDu9^0rMJizC-JJn`H#aTg81QZv?E+hBlQFH9BwclGrAsP?oCVQy;pX#udig{0KL{1
z5pBdJJE|Uz^b;oU;Ys^<2<@8bdtV5VyfeqP<wWzyrO!!%k0zf<c1^p;dU*8#mM7{@
zp;F#&HMnu1hci1BYuayDv6}m4_43vFBS707DVcGr`UDg4xD>2XWlI#6j2g*+nw5Zf
zbIIjA<))UP==D2AnRjq7sl8tL5Fn5|6bST>Nx{Lv{cTd9{~&}vypywswTaV5cN^<f
zO)a|>4lLjGDl`@SMZe*MlY+1HKYA8MZEZfWG2VfvFo-EJ3q}14whOv^Da3|bB^)e^
z(v2w4A2a)$I-StP;-_G^g@tJ9qmf4;5xDn+V3mJBp~JFCXwVWkn?QD0GZQ;@qT$|Y
zhVg;)+g2wUBk@h(A1;WCxg~P6An0xgSeRX;kSnX=Wne!+M8<o=62tRr3H>Ce>{Dcq
z&xaL;bm)lRhb3X4m8=I(O-+4__=3xaOCyUd6D4^l1aqpTM(3h*wjoO@ew{<L8dVaP
zO|0G#&+?8F$j|&-LrincdLb`=q<6k{7Li$~jIX*VzKw`KLESVCxhiRzhB8<Aj-pqH
zvvk8|)LpwRb>h^-V?U~--D=>IsC<xcZ}Ps_SeZdtvs#kQqFs$F>3F#dGHZ$Eny<8Z
z#PDkHZ~igbx0YafSevq3&CN_`6BZU0i3c2~wa+`?_=9&WfG_)Tw_^LadyNBKN0=fK
z=MU>5>G1(YodEeTxEdLE6nt5{v7hoE+ei!{NBP0bl+#Lx<C$<y5UlgRe=X5YxJH%y
z{;V#)!X4wnj0`sxhb+6pFmm&mCQ|p4#oC}%TSaedKTO>kjk$hF`41nBWL9x2?#rvv
zGOg7P_}0_$Zs7YJLd$^7@A)|rbsjMP7jPEt@0I#XadX+Y7FR^$)G{oyr-cN`7C36K
z{SwTLFh9>|voo%{>Y@=UH8dX1E0Iux2MKy(sd1_6mAuntZo>2>waVH`R?AbNvhZNo
zHj`4%3k7aNjp1PGurJcEt$U<>ZPIUdF)5YEUj72SE<j*>RYCVKVC*)do43fny1m1j
z$zx^3?fe9yGyT%1J^P!NAV7Y;|L*2Buzv3gJiiiCtwwH#vsFG_CtWy>6+ND>e9AH{
zH>H1#ckt%Nh)=or6cZ6#0kpO(bE%k)t!c`dw2^Z-vzZ++*k%sW_Qw3YPV<xYgH<pN
zxq6wTt-w2^t@qMmcP>N3ccY4*`k5HqO1NGS3O%S#X|fDZ4^m3GgNU)J)8tcqm&H{K
zod=LUH*zVZJ2j&_gP4uyI<a&SX*gaKv0vX8PV2UQ)p(wB$2%g(R*I;6e(_Kc<Ts|(
z&Hbw1hQ*CT`~6amtpEO-{7BWsfkxgzilijo-5TJ6Dm5e+WAc1Pu~Oiz9Nl3WRlDs5
zeY&3j9(=l`ktS|zdZk1vEK;~0vbS#%&%#6K&O${BsrPa0`Q%Wb?QZUIV*cKJ?o<8g
zeZeK%GI4d6kJ2(uem!HfoK}I}ERn~zhlj(ynl_;pWx<J%M%y{}b%+X9ey5TJsod^d
zWqBaMn{p`j@H1y#&wU%;?Fq5seYJe%(<pCk{EXxX`d-+4CBecKdZ$mL@qppu0p%MT
zR4z){!hoO~$Grk>GuM_3namEMP)(Wz3QM3}bY(HxQOvK*me(aIhXc5f3YI;-;KG2t
z4pR;}$V{+o*yam>)Az*P=k2<6K7ntk8VB?hKa~1mIqA6<f+z0=RRyud*rDYqbuxnv
zo1fQ*itQ5XpReB|xiYcF3hR~>d#2&5qzyz@A7f9?-XL5j??2Hawmgphgu_1A{X^X5
zsA;k7V#4$<baCTtm%&6Lz3)jVX>CjA-_$x&5KUYiOuKGIZY_<7Fr`(IXC&l601@$d
z{$0nT_z+y(HA-DSwg~)pF;~Ha6GIRKK=X<Tj^U|qeytIaH0mFSsr;UVjgD1V^+b<;
zThX8R%T|qaA2#mmQ+rS5B$0ZjkJn-A=PsR};9S4a?DCn78{c^ZmC`$TsC$Zwm!*>9
z><ka4mwn6d&t3FI&okAB>$AJ206w!K?i5K7qfggI#V8`0WFg=l^jnWc>lN`YxRz!g
zq|AT?fnuose{hZT7hIb-J6qVA{SDZsnl?^(JXo*3fY@M_TPuQP8&3uAgRNB=d8aQT
z8oQBX3Psr6&vn1$!&adBIM3PD%(e?Ge3j8R93M~pj=0l?$h*!z>Vq^f>qxTp9xO)W
zau@>%>xX=gXZGugt4Sywh8UPpjt0JW7!ISkXtju4v3>ANS_up>gd5CC-~B%5Qg!>c
zGnBKqka45X11CBtUr)6f@nUa{Oeur3CQN<S-?sL<6dT9r3Y7rb_iOyFg$bQ*NY(dQ
z!I|^gl`*S&UH2$e6D}tg3Ask0Z-PPM<0{b(Tmrtkn->hv*;*7im_PRRrU>!8#w3$R
z-$VbtFE+^=@DVhP;$mT*!-Ix(VI532+#g8~3zIEWR3WYv?0vy;#{DRY<$7kg17F)u
zhrp9fzlQxW#7+dL*A{Ed>QYC1bh@-7+E^nu>|?A~un`jU%<5&1W3_gR9-KN~q_%bJ
zp@3&e^usVyFfeBa@n6n@Sjf8^Z=T8eu!DrT0|VK=IgC6dd>~8iW}{oStWA9(_P_~#
zVJB-_Nyz1thVB-@{ctI4<4*8~ZGvl0$l_mKKc&zvD1HE>S>5En(Qa5JeSVLZckIVt
z)EgheDFIYvG(}f?YAdy%aZ&~9g=rQFD^eqXrSR6c3!}*lDy0uW&Pdg65d@x^(gaY|
zK|+qdX&-ciFE%5d7(?W*8>MJN{EX=WJU(r1ktRh^)+0A#^V2foq&(>~Loa-sSb)?g
zHe$%T++G8)+<1S9PUsA7-eJ@Mr}wdu*BSjZ162Ao5u2~C0FR#Ze)nh9qwg-+J)hku
za#pxa_hP1}0N>vYm>XS6u{h2=hytH0zQ}4gr*)2&SzewY%%%w|F8j$$@`oPG*@@_n
z1}f5t_p*^+%OV^Qtz8H4@%3&_Fw^ddxVN_H^-BUn4~=suTvMnloE(xd97c!EJ><t1
zWDu1?F5NJEl=JLJC;q#vr;~pcl2c^KI}!ZqG)F$Y)@}$S9?i{Iz6Yrs%cR7|uZi9H
znF`XG{Oyfm-p%pCL0<yawltpH$};+L1zMp(w_ZlZJ9~^s^=Ho$><OabO|IOuk(!it
zN!0<i!YwJX@!$})Xpq#pta`{M=GYC~+E=c6z;Zi6>R?1;OPM1}rA{oOcUL(`EW(mm
zXM83(v`k+!_d-4rmTKdoQc(=LTIYK}!(G=>hAzpx<4P8Dr4{h7{uv)t+|Jf0$9$_F
zOb;d4cq}@r%O>F7b!I;PDEV3`x5VKrEeY7MJJsbjR>U5m$VRk0)+^7vuL`W3rs<5<
z%^H7+=S@DXD>+NE%w@Aobaj2UfcYM{pQ?G)WBjtfNzn85D*@~z;$Rwoy5hYNv2+tI
zId*HkMW?2V3n9Ng^&m5T-{xW?NcC2qrBKGLIC-ke{<CS0u`!c5`FB=D(%G1n5jG&J
zd)fB2`t4U{N&Zaa9*Q_}hi$3<dRfSggKFRWyf)E${iCflO~VEkYex2pn(~FO87z2F
zx4~<Ks8Iv%9OJY;9n38CA(h|?#t>u^R=clfp6(|Q4b3(4z1t)uV+99j`n*VG?|8ID
zvq#R~RF=~2F4de??MIDVNV}=6aee}G6&BK?59(}kBeKRTF>`8u5HiruLxoJ%t$qyr
z7e@0>idO*?o6jgmFn30Z5`961-rWUzys9|w=1!^x)i*I>D^+O2(FuOM9zMwQKQk3Y
zLTV*r1W31}^rrpf6S4wn5yO&m6loVb&nh*GRQnH%ILShyy{1c>zJs`-V=@F5WI7hQ
zRHgfu#U6Y%gcpqG!c{cYQ~d7Bu?QjXcZ6m-JlQoXJ8y>?DEYVb#Ik;Kyf__4LI@2$
zSJ8VGgwQG&{gm&}UD%wAiM_a3y$_@0L%S5XB?!5zBSc&9i%yppKDv0?3{1(Tu%m;K
zMK+Qu(}x%MQU@CbgGo`em?qD87kchg>!9~h^%d$bz;}cHs#pdC0zFgw3-JH7<(uo@
zM_fS%f&RMiK=E(V1o2@Q9{|v>qO6J-BqSsvA|fs>E&%}nH8q-`Ac<oPmLXi;J3r
zgM){MM^+Z3r3ErD09jgsoSZ;DJ|G*Lk5*Q!_VygY!oo5#GAb%6T3T9$hQdZhMiv$p
z&d$!Bo}T{xpr|NNObjS96BH02n4b?SEd|xqg4)|bZEY$)et<@R5zx#GXk`VowFNpn
z1YKW)USC1|{otJ))4;&M@bK{X`1rK6w1R?yva+(;+S<0Zw!Xf;(b3Vlxw(~<m7SfP
z<KyG&>+9D))j*)Px3`0#4QLPuZCpxJNX32mG!xNio5RtsLn^%@K|8LILXpognM#f+
z-MOj%`S-(H7jOg`x%0GHl39+UwBgR9i4CIkA=du{<=3tHn|#-@U$aK79NX&umK+ys
zjr2M7_&BDV>>BjMIhB_D&@?JJwegxS9r(4@B-8LiN$D(^36cJTJl2k=2>iA9Tm+>O
z*sv}!kPtC(ibaU4g@GXfa6rMr{ss#%zBu7m3R6zJeHTs9|86|U=A4NEGmnrL$Dcdt
zb$gJGvN&lAf&Q|Zyg#UojX@I}9<V?PQ&?!mE>%38AgU{c9j+)}aCO_uvZ-lT0j6Jw
zXfXUcVW;=3y(?Y`uda+{Eb8$zN-mJMpAKd9DwTC_z>PH60rvI7oxg+a#yYd}tqh#1
zCC$;hHY?#`T@w2TofdFt`qMPP$bGpi7`SjFfgc{V>rL)C)UlvuWNkFQK;$H_6mkh1
z*P48>H_9rj2<f-3cTJDwI<8d>RP>Oia%NOTm{*~6Qk^fjZVbY;+gOjn!{obvBze23
zlG)J%2_!Ze17Om(->2v2I73~-)wpHm6&)|`8P7FYjxKdVABjgtr@pDwz1^>dwp~4+
zlVl~lc~fBSO}D^VSx~Yj&qT0T54MoVXw_8@(rHml^giwW>;&!ytIh>q?pCq^?}yXA
z^SWw`qeLUSn}wi;ts{%=J_v{JNErFNNt7qI)7IrAz}u6a;LGj4Ugz88`P<7Y({+7A
z&td@i15VPfYGGcE6lQj~U^+L6*#+`)xAU&&yFrq!=T|%6`LyV9IY4>=BFVXou~H-k
z!S&grpjOO>(W%<kWBpFB!|N@XWUT|ZpS=fspI7*7ZbMFz4)HL;IDd4na3?V*F0({L
ze`k~f!ZQ4C@6Q4AqsLo7;N@js4|v+X&lY^S@WolkeA2eM^VA(Mazs6Ly7DNu9CJ8g
z#Ks+Z^L-w~cyz;7*fl{i_wTk#DiBYH$V^m78Rg6VlwG0FT+FkEYp0{AHDDy@Mv~N*
z;3>;cVidPozFo1jd}KgTmh<S0CQ;z<Xjl)JnDi+vBOK9kU)KhwW>i$*+`R9Vb=UIe
zIj<F~stkcy4#2jxT3dMj8c98ubb<H-ML`GoUP1p*p|FCf4!M(1Lq5pm=2!JU*4AK{
zxn~u6>aeRafoX^sWSdlA3T0@vdxS$~21BI-d3@@Lyj33X`rDx3ry-RZ1b(3=ph&w<
zaub*~EnD7Q%EE33$Nxq}F!lPIwZ^+|=2Y5lL%_{ubLZigEF66_<zYnMQgVf;d8!o)
z(iaUC`MJmBioqSn)*mjqq!0{Jnc)*fGYdEjF_*`bB|R&Q?qRByYd<~I1UwksCI}}m
z7L@kPdcDowQe|n-sHpOM`{cgg17Gbb;EcYHPP8u^d)c9AG(x7Fn-n+w4jmfLI7)dG
z5I5+R_VG?Il%FeS-N)AsEZ$1F;$wwDwSyn`E{7gHeXRq_(i+&+FqDC)EyBJv16MSO
zLhve7c(<Du%SSe8b*yA*ej6FGLj~7v)f5Gb2r3dI4+V_-kSJQx&?1jwK<dSKPxYRu
zu*X*z&u8-l?5UB4I!v`g&g&jh#%*a5T`e3DUR)V#jToEB+5H-CRtDylWr21YN+-pd
z);X^j8B{3O$pRe|RWG-T+C-zVx=}S*N`>e7QLfhHOIN(c>!zy4KR=AFUBL}5FYd=%
zHN=-JH`Y+u4xB46{dEXW00}-J+xqB5yQL0m0%>*^WriNnyX8=8jwM#4i(eyZJ8WZ&
zTHOK0nDw|h6*&#rk%=~Rx#IC3bF6jRDSt&H8b5eS9bKlMx|1_VN7GaH@_5>U{{@-y
zDt1H~^1FOCAr7xb&$n@NeVpltC{LrAzdhPo9WvHmN{ToD_l6d~kC)9k$O7^et>{_m
zQ*Wg<7s;46OqPLy4oi;4Z;W<c&rOZw2aU;p45eseR=XTcSbqMU-b9zRj`!4?E0}bp
z>Mk>*(Z{is<KgFm8mWqRoR&$4MeB!zQVUS(6=No!@X`pFsmIYMHPXbxAqP*rgLD2^
zB{RhIsG1(gNw_0$!T4x3{!JggBTlB8A>-Tpjeib0f_ab72Davi&AH@9qX>KoJW{tI
zW?yI<#?mI+vPzF-%aU~^Y5mKnMp+2HNpPdJdSixGu>=psyVC7J8m&EZu$<O5mj`9m
z&x^+1JILobzCwWGK#$=hvPM@+%O6db+(?RJQ27&}ES76QCn+WmE6Qq;Uo=8kXieJ$
zDsx-z%2{*#lx@&{yLHr@Afi0a2|h!gCVmdG*}=dM=Sxf$nirV%nxv~ym4Q}t(bWw>
zyk^~GOH)02&kfiTw~JYF;E+f+O)zHx?G43iM+^LRhcM<=9?pbx0!0oN?<_|p?=V8{
zt>A!pZ3dpLy+gMBO~{djeu=>_CbU)qnK;yxY&fg7Ra=hhrn`1v&*4yA;lIcEf!C-5
zNzoNc=hWjyi`&}x>eweio0S-@r7@x7i`Ph!HQ;fXrL1wGw!!k4W|;szixJK@5EpJa
ziaH(jwV@+wjsGfVz33j`DNMLxYs?QiX05V;*9r3(c*skT-sQoiUezjFNd=F|VZL_C
zN*XXeyh=468)XXQtsB(?aPzUFYC`a>uM=K)!_CT&_8oe^yiOMFSvf~^87vVRq2Ecu
zhhlh0Idon792pgAC~awcCT#DdtXC>qHLmnO^_4&M>6k^<4wnb=Tfy#F^gU5RQtrh$
z_A7X7C-+Lp$Yih*Xx#b6O`ZqzRx)6dBj-Te2~JE(a7!C+#dW=?v|#iu4J>;Db@{J%
z9I&Pvx7A35u$1h5hkt%T&}z+ib{;GGH5HP^WxjrAEj;tAzc*5ag!#~D1doe{mePAH
z8^((MgC+10zJ4Cap@|KipW|r%vpiIlUtSVG6sOXPlk!9=lO6c9pH_7ES(>(m<aEqP
z@3c?>Gfu|$cP;yFN86P?#TW9TI~L=xYZ6!MXS589=k$u~=1$*q#kNCGnLP^=RQ}^P
zE2jLeuGOh={2)6_T|J)M<QmV5vUgqLgMar+QiaUS0+25DZFg89&;&O^BjrvMEHR1v
zzZVY{ovuruhQgI$@!}5dqRSx4jL5Qd8@njZft{yzjTSx`iEhu5GQ=lq?S9Re#YToc
zL&AQ|!zt?F0^_u^#O&}!@>TW?&o#4Fy%6>%fvUB7tR?Vw6DLgW_)jx+rO$ZQ#CNyA
zD&JR_=LqAru_o+xRL*ZZvh|4$r-J535B?_13v%B-#P#=C9tWF5V#Ynzt1&EX)x*)x
ze`7pl*2#XF3Gf}q<NftiZl|J2O8f`^!{pTwsXk{rC%X}NZ;l<RLur$4{Fw?iS0FeT
z9OD~a*4BueimK-$8d_4hSHqob6`lA18dD4-hm*Q?8Pkls%BqzyLyO6;3oXS%;B%qq
z_mq)_3r>^ZK1XTV!CCwyMjQs(57Fqerl%2K^U*p`**#?>FK{?qOiu#0+B)ToErO#1
z$=|D8k<ATB`I^W)crXxD6r;^S@~Q<{nyWG$SYj6cgwKgfakyq*V)F-tFchtU6``){
zW$6+EoC&ic#?7az@3M!(Gu(zv*1f%IAhLpxSlZxRGrW^}Y{i%?5sB7DL%7M|pg|_B
z*n-7h`c^VW+<ndZGi=)d@*G-Pf%H5N1@VW1{S6XR%%SCdh=l9v#Z5nY=$4%5w2I(I
z^kzq!+s3>m4PJYQQtnGjZ?we#rc6{{sV#;5Jcki!G9z{t35(Gjui*58;mQF#W12>x
z71x$rV`?9l+Ijf>$_nUhwHQ!B7@Byh!4a5lM%Te3MSn4*lLPi^X;nUkgu=|wsZ+zP
z5(De1#5f!Wx?u2p02+fc424<_dS-w;HC|wv#|bO~V^+ZXh?5;&Ap<%f?@V-=?!t4^
zfnYb#@A#wI50edKnIj}s6O`R(cAWC<9L&;LZ#Hb`b)Vv?Y6kj^$7Ro4%646bino~1
z%}Z)35eH)tRk5psv=AL-N9Pc)Y2nwSl#joSRM35*v}#pV`@M~+r`4{ZD{ETQUke%$
z4Bg-UAcH&9+e|GxqbMFw^d)$==NEl8!NXv0QQxGEBS_{TwX9Yl?F(L^g@H@CVCvVc
z%7qm7)JEN(qeB92oSg+YN&|}wq32BF&oWCt6*D3kaV#j={gdU0oafCfGdpa$e0N4(
zXHm^8VM6Fs_3_4vY3HM)SsW)|J}|OC^kbOT<V8p;>asOHeV*M1YE!R|Fsu=8$m=C*
zl12CE4+6CiPndD7r@=Bkjk9#B)mb&(5vsy<^jN!UM5KgqRs-!YDKsXMN|xa9Tl0e0
zAT~97ZgAokU7uO$vKUOVFAZD$mIV|eD$utvBf=Kf>d=ul0oAQbJ-!npan1%{hXs-W
zPYtkMQ2C!-!hR-DRK_V|>~QpkD#zy1M0ea7(^PeMIJE>#B~$bLIp3I__6kRuWj5Jl
zJcGLFGBE&0bl``FtE_HPBTXLL`)K;HQj0Kf%GbEzN7uBR13p6z=$Mnyyy~B}uLYF}
znjFgX-it~n*n7brj-l0X2*~vfhr9wtOe{u&Vp_ytxc3YI6`}9bT;Dp_M1sc-zIgPx
zlfUMTB;hJ?X^Ak5aD4lzO249&=oCp-#>kSA;~YMttG47>;dvZ+qDiA}*OY#GVskwo
zx3qm;-?M(#zG|g1s}M>eFNjL7$6P0xRZ{U;hq2Rq_Qdmm$cgSY(3OmLT`&G_c&S_t
z;Z&Pft&?XPuqHv+MG5SOu=V=-z{eo-+{HFsRNWjZW?@wr<gd2##qY-rJD~+3qZNS4
zrxZ9XUr&I%tNFRxU^_RiSk>y!SVXP9$x-ATUdEm5p{u?T5`yGw*mx&3l9NZXmjiR%
zG3&+F+(cKaMya)=b@fcX!Vht_-LkAKqlllZw&92d*py7$&*{=mJ5}2fY-Vc@28AF=
zopz!h=pZFE5UL@a<8H9|(9skaG(i7Co|Sz~<OZjrVb6HN*&KcY*VSV$+KQ5SB8-O9
zX{P2yH*w~n(lGm9s)fkLv9qD{r6kR>j9#<d`Wux-SeBRlAQ~_}s0h#|kgLv0hG&ho
zP=I;?s9|57py=^JdD9o`?}lN5QBB}sC<v?&K$#{igbf<1A&2qwT|lMZPCE4c!!_%1
zg*x*^Hs;SZs^6QX+urPXUIEi*A?F>n=rYDX42ji9IM|4*MT)hBuBtaPWLE^1jwm5h
ze34`xBJjiteUy8wCFgAki?|>h@of0Vf;44;M&S`7Q{glEAUdjzaBFsr6R5iM3S1+8
zj3iidb5XSDn_@5w*1Hq)`rHQfgHLV@_BpZ-)jY&VR%`KdNJQwP!YIrl@J?|MB(05f
zMl=CN7RTRtk^*@uDH)oe97kQS6;bA6+x4_o?J)9v`Jn*wB{P;OE}CM_4_wwu(wq?h
zjm3R?;&bMQ$FuMZ)F2Zob{8n0Zg(s7-LuYFg{*=lZ)@IPiN)z~&}%MGtiEbgqEYL>
ziEGBG2a;_#>6w7xaJl)tnV4ddlA0EudP{dWXO*wDJlR6JK22nTb$ua)m#|g=`S*1G
zB+lvfkvaQLgILquW^<PWqZ}B@lAp?e5i^$hVr7(=)UT~rfldOBqAc2}UnMx{?$k$*
z!xA>JIGVI;N5vuY?@Fqd+d7C0k(;|_A<a}6$&l&oDk)N4Wi=-!(2hx9!`o~jV+Yy(
zI)s%HlNYTN)(`l9pU-0d(UKEW5n`08>tNLUpL6*4Ohkaqw=J%j^=v)#n*7;zZNK)=
zwNviX@tdqYf<1Z8y4HBL3JD8C(&ZR?se|O-w`SB76fu5889)3i!b22y6Ut6X9aWbe
z51E)`)ZHP44j$HQ6pm26ATPdmp1NPJtW~W8=;Ld80<Y@|uhx6AJzq~1e#&f!ZOz@2
zKj`@_>be7U2%q7)G`{>A;<>Q(rg!drKdBjTL6x@lK)+Gg%=yg4dFc~>Jo`x{Z)J@<
zLX9)9=la^|P2=#n$NS>99@I@bo`14=jPNAXitA|3?xrIC*t<q-jjed7EI^W37iYQ_
zON8L7c;LfNBsOQdv?SQK)Tk~wTz!l1Xk76mP8bc~9gZ5>rXv~bgqxsu+V+;PB2DNn
zLhZ#g%%k{8x95NjOpma^{&e8H8VK$Ky6>`^4Vqsfi6|{w+`F@JH0#`^=`zIX=SaWn
zMHw`W<{>zaZ_isiYv=QQpx<c$6W+hl+k+~(uIYoeChgNN&CnDkg1E<rv@4)Esb;cJ
z{1DLqX|9+jB@QC*bwN^Yml*GL6wLks5~_&I^YJm%B&PQ}C*tWi-J=N5D8h}}k{F_B
zK3=fY($J8AOPx*5r*lTPRwAh$zFBrw;B+ZGjmCmFcNI%>F%#~c2RZUPII81BDDC%E
z>YdF<Zf!x;k>nB&4~E~)lls1uOS~MSL+#f99k*BL(2I~EX^wUwJ+=ZQgo=~?P~0(4
z*~HKi<aEww7k)&@hB;V$w^m3o$m&jvSSn%IjCA#W+$r58+fBn<f+r5#t-QT`78~W&
zswR*+cw?`%GlaR^uvswzLs*2ZK45=R)VS=XqzM7rL(krjWNtU1Gap%>z@QEL7@ZVA
zh(stuLz#40|0sk6kPcCI)sXD`K=tcx0EMuu1})ew=^cFF0watL&)Kofxx&H-D)7r+
z=Hp$E@GP8e1{5j^Bo|?y88s3j8a;KIaEYLXKPGBb5BkhYFeZF<PyDWiQ<$>c3JNKT
zLzMw^Qjy~Q{I{1B?B1fTb;NZ6Ri!2XQSxi6P+~jh&kC2FK94=w{^m?;(jVu+ow7c(
zZA&3{V2IC=N0rICO#NeeSn7%pfuxR~X~tgC2SMPj$ZsGNyG7Dgm*p+km!+ENgyg4)
zN~>sGH&Tg7vaOxQu=jUzPNKqEiv7*qz#18n?yz#=OVG!Ueyn<;2YYb_z3C(Xl{*pS
zw4T_+jAK8aD^@uK{3<qIg;<m2k*7?Mgc^6>Y^032j1XyMIX4<+)TF*?#j1IcdJ8e-
z-HK~(qE(0O!BTxlIm69#f1zDp(yz}^#P9=*#c3+tF0*=SYJN#Bm&t+99}0Z9NxO9G
zwaw1LbGuMtCRv?yX9hMsA;3#O8v3hoi@`z@^&Hv_BF4?i+KfMaYsNxgg2r{XCZp@#
z%w1A(<Phoy1F3nHvlPl|y~*mTOOO~Wj!%i9%1#x*{tQfYN#SJy+YoKRWclIwHHf=)
zsP<v1?5!tSPgJzj1`{7<YF#pof?UPO>DLz$gM+IXguQ6$R{fADyjTeUG7YDgWN|8j
z|LIH60&!nwZY4R{jm{vBFWx}aTKKiQQ)I!tlHQ6PfQVHP(is)Sh^G>@E9XL%OC@pT
zAybT1jEH_%5mB5Fs4X0N+7srlp<)|MVsbM`=`I>DQvED70OKA2LWVUZ6+<U7dgX2<
zFg=!u!S|lVKoT+P(6It$>qS006|;UIFOl+cH0aZ?WO*H|{DbnrsY<I#@&1k1IOx2y
z#fNAl@ur*YqH|JwV~xEq=HilB=l+&j9bVgxH7B24$Z^HDTez)=oLVv0GH`sCKmTxv
zeq8TDVal+P`zTq$<Hs=8MU@<Kn>r6q&jU}%VNM>s9Po8~k)0BVN`-4#1(8XGC7UK}
z2Yd*kLoBIWFl@(dhErXU3(|=OEu{vWiXXO(*J@{hy3XGnZ7?<(3f-?bwQrTKL=ad?
z%Dp%pyg&d+A!X)Xr}zcFLe=%Z0*rGuK5y_4ONvV!@g@jcLlbi110TElNv9#BgCB*X
zbk#}sLn#spITwKsy?PRGyiBhTp3T+P^|S`xzV~*-NlUQ>5!#(rVd^#~aq4L~JwQp8
ziWb>Ni9_sMSFO){SHhb{AVNgA`_|tm$_XnEJlwdT!kj)k^o-irZJv;u2>J|9vGf!&
zCaE}{gG?(D5(?~uV(Vu)FjLQ8QyS{-oml9UPV&H4APoSKk~G}1#FqTABQF@7px>ei
z9?lt$<Q_orJNDBEfB9YXA&>)Vk)+N`*;!Get@Cknx=F7|i$o2g+WJ5w@L*d(iim0X
zu-0{Zm?nf&c(=_W{Je2I6-u?hiXW15Y`_CDurvYcQz#=exx0`T<#co;OR?EtLy?+(
z=nBT(&LHsM1C??RY8UpIR%}8Gq%xu1y)svo>^H#vl&-n^!<Hf`Z!8C?dL?0&BVJDC
z{gk4tHAJ<Ymu0+Lg*9t5?TmFb8Tq&qsiUaNXGxKgyA-&W{kXO~9V}kx-n;1O+|ygj
z$e=zIxK$r0@R+!!V`$`fjp~50ReJ~7YKm!#RW)$Q+Dqfd(?ng+COqm(Pf4b&k(CVB
z2nX^?er7+zv~<6{tg00Jt%E1LN3kBME$nQj2C+n{AK9>|81Ah&@C>epMT}R9k&lHW
z`|xd8zJz<Wz$)DWp)#tgo8Dhv<bL=vEz6eZ((vbn>eZGVqF@yCPfr*jMzxG&0`sqo
zNY}sy-*<SCw6ATZMb!<NLX1R&b0e;bZ(Ev9o4C$DU<YETc7gW9Qr)-4w9v3hCGx2F
z@CS1DjYntUg5|A=FI$)qmLl*xsfMSTP5gkX)3er_t{0c23_)OcR|NA98>UW=KIb$0
zv4g}}!Ac21EF)!%!A2uhou@JF<dMU6-fuvzQ|nQ#=n!Q|wgV9z!M#rPFZ*AQB`$ph
z<mh{cqKX8pEWA`QRk>~h4^Oc;gLHCgL{Ipf8s!njKaJ;&d{(Fs6Rnt8P8EuRlBRid
zAawkvXEZ)(elQ(~>W!Xe0V_WRiAQojDR63zxs>#p&*#_K<h%U2AzKKO)bDDG+6&-(
z$mmB=bz%cv4HWJBK<bS~V%n^hPHb$T<rI0;E3?~FWm^S>Ft>!yQCVJDuP?u}iXCU@
zFrZq;%}m9%u9fs_b$}mw>G5O)$6;$wDw#oSu=3MPC#5af))L2LH6+@u=q>Z`(JW~>
zA8}8>Sz-wE2&XR&EY9qXm1-x&bsnH$yO90<)KV{?rKjP%=hmrdxyhS$yc<NXr@@-*
zO?mb>gT>xcZ=P5wlR~gx?sfU+c@=T70-uoDr^0oY;668)hQcQ;oqi0EOj;Um9q@7m
zlhyWpkVWQiRpXV#^0!P#rxwkTm2H(J_U~oK<Aa=s|3X23(iu?Vj!hS#P6;TtH}hCK
zH^`6%XiRenDujQ=Pf5W4qZWS0e+`!q#7$keki6uB9semY0<v;3x#xr8t*8$=It@|$
z$4rSFdl#lxA`j@J0CGjWZOjG&<j<SJ1P|$Zg*@c_{8gJlUjhLNXkKfB8IQhsW|cc<
zinr+#v^a&9_uipI-iVp*hrlGH=ejI9^z$VV_F!jgwpUh+M!qAD#4|6S?G7`+6kNqv
z<*Xa;!k+LvX+QB~W>&J02TTZJ3$VNb;db&%Rgl+uvvZsKIl7?<ZRU!If3pMhy1rC8
zsN&}7OYzyrQZAqUhS6-O3Cs44?VB&w-U?<iG&L5Q{=<<RuKzhV0D#)y14Z_%5|?z1
zq1R9aa{y5Y;r3nqi3%=eA!V`F8?r|UtnD5-yHYmkH|cp{8BBhU9YySOQOCW8##j-4
z4$bglTeo)N6pItopY%Fr=*`uGK{81e7OJsHt!A=^<dA#ei@z`rtG`)q`gZ2;@Tgl1
zo_v+esn;S&i*t1D24bhy(^3VY&ch_n78f)Nj4E%hi5uUU63yl(Nl>_p=ADN9SRKU*
z+Wdjg_d7?_dqvTn&wv<Q6>#2C)yo6j2-;8l<oep?dXmlmjOfWRJC+u3CW1cx-RFX+
zfy!JhytF%{xpFmT{1=ZZAV*c2(7j7oI5UGGb9>kHYG(033P_)?<trH9CfLL`x_
zN3WF~!P>&Nh&GXv07pr@9(M7+e}<n(S{b}Qb3VNYl2(Ca|8cE~zOj2u9BeKYJ<5$B
z;5(bj;%&7g>Z(AwjKlJSK1RY%HX8h-)MB`<GjO-6AQ*eoMxtTca{c>66cHeTCv}|)
zri2#GYJT6g3P6(66)2jLC2Q5#b;-XzRsr6xx3g~MlLB9DpWNbZE<U`_&>qy?U7U$P
z><|aYV#Ihsvj~wwFTu4i86`=67{Wn_eT>?O*ptTkx-GZbfSlX*g}v+cjqV$943AK$
z_=4R=dh7wNy8*eCVtNju=rVNGG7*HMPOr<bNu7%s4sa&0S007Yf-QZKynJ;yX&;WA
zqd4c+z@8e>Cz7fVSC>hV8M<5cZCP&nPNCaL4{+5%@F$I;^n*OT$}}W@_cmU{In)n0
z-G4VG^4p+#-~0@1dlhnJ=!f&os*2Lya%8GDoZ{~1)|UawxxBFJ1F~oidQ_5}6M?SS
z+fe^ITOdZ8d>Mi^^-=KG?-;T2fe~H<h4ZbfHEgd4fmaNG?Mu<QByYt*sVUtdV$-Y`
z13@}~UHhqx>Np~`!@h5l7hLf;)rKlgdTOmM-1`O30c&^D$I)>*6zpS~AA9pDf1)6w
z#A?7nV(^)W2-D5EARYED>lnI~c@;9bjDqJa1GpF9d^|ZX<?5lXb^5x(M^c=iYp$tC
z^{I||7&|UXzibd*6;Ea{cu}+TT?M5^ma;OIj<T$3@|K&mn09xbAtOC;Mc<e=Wo-yW
z=Aj~-dmRxYN!*fC_C%JKnwf$L{vcN>|MvI^EtqGEmNiapJZA2JNKey)f#;8h`fGNB
z4?tg4DR%dD=iq=tLBABtEbO^WncvuDOu-i*sdUA3LTLn1R^xcTh#(z~J%^$x)T>J$
z7j@c74yMDY+L8>SBT)4|*R0sroUQQ*TD1(k<*IwNaWs5+-dU{dFd^CWv(yErlp18*
zShTMfF<(;P@8YaKSlvP55Y2b{neHWu0?zL@9${;}`vlGTXt1c)b5QJ_+)BOK*8l2E
zt}HvQbZ?OiTm2ibSnUu0S{KN9&2XqyMQ2*CP1m|Gu|8%MsQ{}ONI8y1rV|DQ^fFer
z*!nsh_`{$d{+jd-Q;8~<2?Gzm5@09cE~2_sjiK+xNe$eUlPQe~(!{T#D+~|;u*w+e
z2WBoF)94EgXC3%WO%_NgVLw$&HG%D1cq6gVtOn$^uwnq?#>@p2AAEi<hht!@#ER5U
z7|q7*J9w9f2v(6z8jBV^NMY<naiN53ho6|kG`OUCOwoy>r+T!zyc40iq4_pFYvVU~
zFl>0E*W)C^?1|lfaJ&6i0vxyy=L=Th`SfRjS{AZE>~EnuR%$;iKS{d{D(X8-c6^8`
zs6Jm~o+E3w@S2>u;L81>TcK)L^%d)-HV!W(ban|v?82Uhe)EPRZ6eBciItxM9epow
z6Nw+^S`%_jmGZ)}0>gcEF{!p}vCvP9ftcFt7nhD`MtdobZlx-b9zf$X$mDNR#3K-1
z$(pZmFN?Q{WAO6eKtGvZ^~oeC)ACkhcy~6penh82b$6{kOuocH-D<x1{RgQ@-9?Xk
zavao(>#3)fkOf}8_z1+5efg<30LrY>kI&*%loI*Wek%Ih{<A&i*K_;kuoL*Fk20u%
zJN7cpauK^92PO9NP=TwyHcfaO4UR@&(JF{sZPbf;%8kJk?Yu<MDH^>~Q1wS%&BIis
zg#wI)SsRAZZ+`3pFWc@lYADwXltl*)eJ0BlF%K4apTR8{i%zNJ3smO_1CPe=Q$ONK
zj+!`|(#3Z46$AWH6inEY!q{>oWS{G3pK=L^$Os}+eqSMDGXOrn;!m%SIBD1TC)6NO
zC<Ho`dBj>&Qy|{?AE1L*f5xd><AskG1ej3am$EIg1}8Ua1Qnue+QAvP%A*yw8NM&l
zU*Zb+4n664Sulj&cvT=tMmw`yqd>5%lt;9#6tF=>u&(52Dis^aeP==K;*k6p;KGU%
zjC=Q*LXLxd55!nAT1p1vQmyiWS6iZk5YH=vMA?nLFo{g3<~;T&?w{#g{eqliB>a>R
z!~<S%Rwb}qt8uYMkBrD|FQ{M^{cF(wj8od!7uBhf+V+1#G4hqQGnvL)z7-lg`7}r7
zn}jJc8xy+bagg~@Ym(z;ZP`LT&11G}Hgr(r9L=^C@ZZk8uY>O}I+rMJdDCYb`)!aN
zqxZ-DmItmN6bLr!$`<dB@s6h0i7x7u;{_*1n2(aS_Oo+?eJ;9Da8PsWi0n4VoYjO6
zqDWv1YXnxpWU=ZH-W>Cn-Gi?v)gnKMu@RJpe6(4+<ew(QK#gOyt^Y_PIZ;f&Gl|wR
z5%aqMFf0NeO$CPxU0SZljA1ZYnEia%SdSf?S1_O#96D@;9`oHYK6!CFT@>`X!)1c(
z`7^=p)Zm>Oe{rt&7ZU9Bhv<>Ji9<#r!SK44k~qXOO3`VSz@<8X{F1e+?|{?KCS&9>
zo*ZL>k4Lh#8AF67m7(*cQLD)Wvu&d?Z2k%d2rrl&zMdb`^diTSA~~e}y@>coj*}<`
zV-a3jnX-hKenzh5hKSNT;YwNX<V7AgE$J<)tgghr=|MvE#O3k9Dm|1)eydEgqR6k?
z)_QWSSD>OStY2BkT4LUvE-X<+ao(u$-iEHkhG002wk*zn&BU5f?GJ7%Gg&-Ce=*|$
zuzC}z(@IY^_wpgC7DhNNzGsaR8`vQ->Q7HoAi*7!BX5WmB$j6J$LR&x3VavIdI+Z&
zP*k6b&1g4Lf5xsMh;<zwlajk2|svnUNFY#8M@peUR~k#KGQ{YPJ*Z2Yj<N&-qk&
zOcScLbi%h2sIyWUvfO9dK?)vw;T}Uk62_fMvJ)zMkFF_B%>fM#Y(4*E+_a4!^%oAK
zi!&31Ud0sfMMVz-Lgh-KA+I?shcSzgJu=;w9#1D)*gvhtPSlfcM^K&&5X7crZC-37
z6=byd$H-b<|4#g3L#c^P<lfel7RJbro;|)P8J{@Ld$w%wu55TpXo(=aIPeh|tli{g
z{$nNbW4SKpS+K?w!_pE^EHkOA_`wxFLmqNeT?Tj^KC(rG!`7%A=jQX8Cs$F?s{<yf
z067HBhn4~6b6#_{q<}?{SpA5#_o6NpQOEDqc<OlB94I9yjLsY`Z&|BYb|Gz*x-;KN
zTB@WrD(hM1JSlCv`_KlY<GllpWRO+#EccFNrh~(Wk7e*{)<ka$mlsA8tmc!Cq|Z=)
z(w*x!=;)PZYU`Djwd5mI%Pb>+*=9Tj*yoeuTpQ`hE}xY0o|HV&p;RpAlS@mjk>r%P
zacllGzp~5T-kEUwG~<&lU0R9#uW8#w=kM+>=}ahQB|-rFe<mj__u0VxL1R73A(1U_
z4fGOb{@zK}TqAs@)(f^Sp7FR+ml1gv%zHSVy4$T8pHZpOdP@AZA=DnFnz&>QN$KSZ
z4)A28IDMq}Ju-y_d+Yaab0@(zWQgFev`D^32nxwlru0^aaXhyL$r@6YpOn0%2BZIc
z;b;?YU2H%~E-@8-3_={k(s>W`A@fjkXrCRPpwy=~8;a@!rGd&Gh(7ZWJiz|_^Q-jo
z$JS)VfMivW57pqka8S#p0UVR2Ycq&+S=amK^yew(N_J&(cO)s%RN__Mx3|CDBXbzL
zV{33gAd@7BKjBLdkXWF9$HWC&bH|HFF#Qvh)-FL1A0t8YuZ-}&7=n00$^VgBOA*8y
zO49vP5+f=85A&Dw=HF6+cs?nHf0c;-!_-Sj{lh>=b0Gf}h4;UMs^UeZX-WQ){HNdk
zU%^!WCbRe`X?%^pbpF%T^Dhtne-rPY?0@?6|DSLA->#tlt_K2<#=lBalm69u%G#PD
u^PfO3hCkkve;pV6o0uu$#biGGt@NLt6Y<{~f4D-*_%InF*l&D)N&gq+#CL4~
--
GitLab