diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 77f6a476bd055fdb6db221ffed5571d80627d119..a924b0fb33ec4aed389d23f8039fc5fbbde170c0 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -2082,7 +2082,7 @@ msvc-14.1_Hybrid_dbg:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.1_Hybrid_sp_dbg:
@@ -2095,7 +2095,7 @@ msvc-14.1_Hybrid_sp_dbg:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.1_Hybrid:
@@ -2106,7 +2106,7 @@ msvc-14.1_Hybrid:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.1_serial_dbg:
@@ -2120,7 +2120,7 @@ msvc-14.1_serial_dbg:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.1_serial:
@@ -2133,7 +2133,7 @@ msvc-14.1_serial:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.1_mpionly_dbg:
@@ -2146,7 +2146,7 @@ msvc-14.1_mpionly_dbg:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.1_mpionly:
@@ -2158,7 +2158,7 @@ msvc-14.1_mpionly:
only:
variables:
- $ENABLE_NIGHTLY_BUILDS
-
+
msvc-14.2_Hybrid_dbg:
diff --git a/apps/showcases/Antidunes/Antidunes.cpp b/apps/showcases/Antidunes/Antidunes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..5ddc9cfa174373e58bb5c01d8146daf704f4a26e
--- /dev/null
+++ b/apps/showcases/Antidunes/Antidunes.cpp
@@ -0,0 +1,1413 @@
+//======================================================================================================================
+//
+// 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 AntiDunes.cpp
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Jonas Plewinski <jonas.plewinski@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//
+// This showcase simulates antidunes, i.e., particulate dunes that travel in opposite stream-wise direction.
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+#include "core/SharedFunctor.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/StabilityChecker.h"
+
+#include "lbm/PerformanceLogger.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/free_surface/BlockStateDetectorSweep.h"
+#include "lbm/free_surface/SurfaceMeshWriter.h"
+#include "lbm/free_surface/TotalMassComputer.h"
+#include "lbm/free_surface/VtkWriter.h"
+#include "lbm/free_surface/dynamics/CellConversionSweep.h"
+#include "lbm/free_surface/dynamics/ConversionFlagsResetSweep.h"
+#include "lbm/free_surface/dynamics/ExcessMassDistributionSweep.h"
+#include "lbm/free_surface/dynamics/PdfRefillingSweep.h"
+#include "lbm/free_surface/dynamics/StreamReconstructAdvectSweep.h"
+#include "lbm/free_surface/surface_geometry/CurvatureSweep.h"
+#include "lbm/free_surface/surface_geometry/NormalSweep.h"
+#include "lbm/free_surface/surface_geometry/SmoothingSweep.h"
+#include "lbm/free_surface/surface_geometry/Utility.h"
+
+#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
+#include "lbm_mesapd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_mesapd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/InitializeHydrodynamicForceTorqueForAveragingKernel.h"
+#include "lbm_mesapd_coupling/utility/LubricationCorrectionKernel.h"
+#include "lbm_mesapd_coupling/utility/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
+#include "mesa_pd/data/ParticleAccessorWithBaseShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDataHandling.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/LinearSpringDashpot.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/mpi/ClearGhostOwnerSync.h"
+#include "mesa_pd/mpi/ClearNextNeighborSync.h"
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/ReduceContactHistory.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/SyncGhostOwners.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "vtk/VTKOutput.h"
+
+#include <core/waLBerlaBuildInfo.h>
+
+#include "AntidunesBoundaryHandling.h"
+#include "AntidunesLatticeModel.h"
+#include "PIDController.h"
+#include "Utility.h"
+
+namespace walberla
+{
+namespace antidunes
+{
+using ScalarField_T = GhostLayerField< real_t, 1 >;
+using VectorField_T = GhostLayerField< Vector3< real_t >, 1 >;
+
+using LatticeModel_T = lbm::AntidunesLatticeModel;
+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 ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithBaseShape;
+
+using flag_t = uint32_t;
+using FlagField_T = FlagField< flag_t >;
+using AntidunesBoundaryHandling_T =
+ free_surface::AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >;
+
+using StateSweep = walberla::free_surface::BlockStateDetectorSweep< FlagField_T >;
+
+const FlagUID FormerMO_Flag("former moving obstacle");
+
+// empty sweep required for using selectors (e.g. StateSweep::fullFreeSurface)
+struct emptySweep
+{
+ void operator()(IBlock*) {}
+};
+
+// data handling for loading a field of type ScalarField_T from file
+template< typename ScalarField_T >
+class ScalarFieldHandling : public field::BlockDataHandling< ScalarField_T >
+{
+ public:
+ ScalarFieldHandling(const weak_ptr< StructuredBlockStorage >& blocks, uint_t numberGhostLayer)
+ : blocks_(blocks), numberGhostLayer_(numberGhostLayer)
+ {}
+
+ protected:
+ ScalarField_T* allocate(IBlock* const block) override { return allocateDispatch(block); }
+
+ ScalarField_T* reallocate(IBlock* const block) override { return allocateDispatch(block); }
+
+ private:
+ weak_ptr< StructuredBlockStorage > blocks_;
+ uint_t numberGhostLayer_;
+
+ ScalarField_T* allocateDispatch(IBlock* const block)
+ {
+ WALBERLA_ASSERT_NOT_NULLPTR(block);
+
+ auto blocks = blocks_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR(blocks);
+
+ return new ScalarField_T(blocks->getNumberOfXCells(*block), blocks->getNumberOfYCells(*block),
+ blocks->getNumberOfZCells(*block), numberGhostLayer_, real_c(0), field::fzyx);
+ }
+}; // class ScalarFieldHandling
+
+// function describing the global initialization profile
+inline real_t initializationProfile(real_t x, real_t amplitude, real_t offset, real_t wavelength)
+{
+ return amplitude * std::cos(x / wavelength * real_c(2) * math::pi + math::pi) + offset;
+}
+
+real_t getHydrostaticDensity(real_t height, real_t referenceHeight, real_t gravitationalAcceleration)
+{
+ return real_c(1) + real_c(3) * gravitationalAcceleration * (height - referenceHeight);
+}
+
+void initializePoiseuilleProfile(StructuredBlockForest& forest, const BlockDataID& pdfFieldID,
+ const ConstBlockDataID& fillFieldID, const real_t& averageBedHeight,
+ const real_t& averageFluidHeight, const real_t& forcingX, const real_t& viscosity,
+ real_t amplitude, real_t wavelength)
+{
+ WALBERLA_LOG_INFO_ON_ROOT("Initializing Poiseuille velocity profile");
+
+ const real_t rho = real_t(1);
+
+ for (auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt)
+ {
+ auto pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+ const ScalarField_T* const fillField = blockIt->getData< const ScalarField_T >(fillFieldID);
+
+ WALBERLA_FOR_ALL_CELLS_XYZ(
+ pdfField, const Vector3< real_t > coord = forest.getBlockLocalCellCenter(*blockIt, Cell(x, y, z));
+
+ Vector3< real_t > velocity(real_t(0));
+
+ auto localBedHeight = initializationProfile(coord[0], amplitude, averageBedHeight, wavelength);
+ auto heightAboveBed = coord[2] - localBedHeight;
+
+ const real_t fillLevel = fillField->get(x, y, z);
+
+ if (heightAboveBed >= 0_r && fillLevel > 0_r) {
+ velocity[0] =
+ forcingX / (real_t(2) * viscosity) * heightAboveBed * (2_r * averageFluidHeight - heightAboveBed);
+ } pdfField->setToEquilibrium(x, y, z, velocity, rho);)
+ }
+}
+
+/***********************************************************************************************************************
+ * Initialize the hydrostatic pressure in the direction in which a force is acting in ALL cells (regardless of a cell's
+ * flag). The velocity remains unchanged.
+ *
+ * The force vector must have only one component, i.e., the direction of the force can only be in x-, y- or z-axis.
+ * The variable fluidHeight determines the height at which the density is equal to reference density (=1).
+ **********************************************************************************************************************/
+template< typename PdfField_T >
+void initHydrostaticPressure(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
+ const BlockDataID& pdfFieldID,
+ std::function< real_t(const Vector3< real_t >&) > hydrostaticDensityFct)
+{
+ WALBERLA_LOG_INFO_ON_ROOT("Initializing hydrostatic pressure");
+
+ const auto blockForest = blockForestPtr.lock();
+ WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ PdfField_T* const pdfField = blockIt->getData< PdfField_T >(pdfFieldID);
+
+ CellInterval local = pdfField->xyzSizeWithGhostLayer(); // block-, i.e., process-local cell interval
+
+ for (auto cellIt = local.begin(); cellIt != local.end(); ++cellIt)
+ {
+ // initialize the (hydrostatic) pressure, i.e., LBM density
+ // Bernoulli: p = p0 + density * gravity * height
+ // => LBM (density=1): rho = rho0 + gravity * height = 1 + 1/cs^2 * g * h = 1 + 3 * g * h
+ // shift global cell by 0.5 since density is set for cell center
+
+ Vector3< real_t > cellCenter = blockForest->getBlockLocalCellCenter(*blockIt, *cellIt);
+ const real_t rho = hydrostaticDensityFct(cellCenter);
+
+ const Vector3< real_t > velocity = pdfField->getVelocity(*cellIt);
+
+ pdfField->setDensityAndVelocity(*cellIt, velocity, rho);
+ }
+ }
+}
+
+template< typename FreeSurfaceBoundaryHandling_T, typename PdfField_T, typename FlagField_T >
+class MeanVelocityComputer
+{
+ public:
+ MeanVelocityComputer(const std::weak_ptr< const StructuredBlockForest >& blockForest,
+ const std::weak_ptr< const FreeSurfaceBoundaryHandling_T >& freeSurfaceBoundaryHandling,
+ const ConstBlockDataID& pdfFieldID, const std::shared_ptr< Vector3< real_t > >& meanVelocity,
+ real_t averagingFactor)
+ : blockForest_(blockForest), freeSurfaceBoundaryHandling_(freeSurfaceBoundaryHandling), pdfFieldID_(pdfFieldID),
+ meanVelocity_(meanVelocity), averagingFactor_(averagingFactor)
+ {}
+
+ void operator()()
+ {
+ auto blockForest = blockForest_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR(blockForest);
+
+ auto freeSurfaceBoundaryHandling = freeSurfaceBoundaryHandling_.lock();
+ WALBERLA_CHECK_NOT_NULLPTR(freeSurfaceBoundaryHandling);
+
+ getMeanVelocity(blockForest, freeSurfaceBoundaryHandling);
+ }
+
+ void getMeanVelocity(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();
+
+ // use separate variables for the velocity in each direction; syntax meanVelocity[0] does not work in OMP-macro
+ real_t meanVelocityX = real_c(0);
+ real_t meanVelocityY = real_c(0);
+ real_t meanVelocityZ = real_c(0);
+
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ const FlagField_T* const flagField = blockIt->template getData< const FlagField_T >(flagFieldID);
+ const PdfField_T* const pdfField = blockIt->template getData< const PdfField_T >(pdfFieldID_);
+
+ WALBERLA_FOR_ALL_CELLS_OMP(flagFieldIt, flagField, pdfFieldIt, pdfField,
+ omp parallel for schedule(static) reduction(+:meanVelocityX)
+ reduction(+:meanVelocityY) reduction(+:meanVelocityZ),
+ {
+ if (flagInfo.isLiquid(flagFieldIt) || flagInfo.isInterface(flagFieldIt))
+ {
+ const Vector3< real_t > velocity = pdfField->getVelocity(pdfFieldIt.cell());
+
+ meanVelocityX += velocity[0];
+ meanVelocityY += velocity[1];
+ meanVelocityZ += velocity[2];
+ }
+ }) // WALBERLA_FOR_ALL_CELLS_OMP
+ }
+
+ Vector3< real_t > meanVelocity(meanVelocityX, meanVelocityY, meanVelocityZ);
+ mpi::allReduceInplace< real_t >(meanVelocity, mpi::SUM);
+
+ meanVelocity *= averagingFactor_;
+ *meanVelocity_ = meanVelocity;
+ };
+
+ private:
+ std::weak_ptr< const StructuredBlockForest > blockForest_;
+ std::weak_ptr< const FreeSurfaceBoundaryHandling_T > freeSurfaceBoundaryHandling_;
+
+ const ConstBlockDataID pdfFieldID_;
+
+ std::shared_ptr< Vector3< real_t > > meanVelocity_;
+ real_t averagingFactor_;
+}; // class MeanVelocityComputer
+
+class ForcingAdjuster
+{
+ public:
+ ForcingAdjuster(const shared_ptr< StructuredBlockStorage >& blocks, const BlockDataID& pdfFieldID,
+ real_t targetVelocity, real_t externalForcing, real_t proportionalGain, real_t derivativeGain,
+ real_t integralGain, real_t maxRamp, real_t minActuatingVariable, real_t maxActuatingVariable)
+ : blocks_(blocks), pdfFieldID_(pdfFieldID), currentExternalForcing_(externalForcing),
+ pid_(targetVelocity, externalForcing, proportionalGain, derivativeGain, integralGain, maxRamp,
+ minActuatingVariable, maxActuatingVariable)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Creating PID controller with pg = " << pid_.getProportionalGain()
+ << ", dg = " << pid_.getDerivateGain()
+ << ", ig = " << pid_.getIntegralGain());
+ }
+
+ void operator()(const real_t currentMeanVelocity)
+ {
+ // compute new forcing value on root (since flow rate only known on root)
+ WALBERLA_ROOT_SECTION()
+ {
+ real_t newExternalForcing = pid_.update(currentMeanVelocity);
+ currentExternalForcing_ = newExternalForcing;
+ }
+
+ // send updated external forcing to all other processes
+ mpi::broadcastObject(currentExternalForcing_);
+
+ for (auto block = blocks_->begin(); block != blocks_->end(); ++block)
+ {
+ // get the field data out of the block
+ auto pdf = block->getData< PdfField_T >(pdfFieldID_);
+ pdf->latticeModel().force_0_ = currentExternalForcing_;
+ }
+ }
+
+ real_t getExternalForcing() { return currentExternalForcing_; }
+ void storePIDSnapshot(std::string filename)
+ {
+ WALBERLA_ROOT_SECTION() { pid_.writeStateToFile(filename); }
+ }
+ void loadPIDSnapshot(std::string filename) { pid_.readStateFromFile(filename); }
+
+ private:
+ shared_ptr< StructuredBlockStorage > blocks_;
+ BlockDataID pdfFieldID_;
+
+ real_t currentExternalForcing_;
+ PIDController pid_;
+}; // ForcingAdjuster
+
+int main(int argc, char** argv)
+{
+ Environment walberlaEnv(argc, argv);
+
+ WALBERLA_LOG_INFO_ON_ROOT("waLBerla Revision: " << std::string(WALBERLA_GIT_SHA1).substr(0, 8))
+
+ if (argc < 2) { WALBERLA_ABORT("Please specify a parameter file as input argument.") }
+
+ WALBERLA_LOG_DEVEL_ON_ROOT("Using generated lattice model.");
+ auto configPtr = walberlaEnv.config();
+
+ // print content of parameter file
+ WALBERLA_LOG_INFO_ON_ROOT(*configPtr);
+
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream file;
+ file.open("parameterConfiguration.cfg");
+ file << *configPtr;
+ file.close();
+ }
+
+ // get block forest parameters from parameter file
+ auto blockForestParameters = configPtr->getOneBlock("BlockForestParameters");
+ const Vector3< uint_t > cellsPerBlock = blockForestParameters.getParameter< Vector3< uint_t > >("cellsPerBlock");
+ const Vector3< bool > periodicity = blockForestParameters.getParameter< Vector3< bool > >("periodicity");
+ const bool loadSnapshot = blockForestParameters.getParameter< bool >("loadSnapshot");
+ const bool storeSnapshot = blockForestParameters.getParameter< bool >("storeSnapshot");
+ const uint_t snapshotFrequency = blockForestParameters.getParameter< uint_t >("snapshotFrequency");
+ const std::string snapshotBaseFolder = blockForestParameters.getParameter< std::string >("snapshotBaseFolder");
+
+ // get domain parameters from parameter file
+ auto domainParameters = configPtr->getOneBlock("DomainParameters");
+ const Vector3< uint_t > domainSize = domainParameters.getParameter< Vector3< uint_t > >("domainSize");
+ const uint_t wavePeriods = domainParameters.getParameter< uint_t >("wavePeriods");
+ const real_t liquidHeightFactor = domainParameters.getParameter< real_t >("liquidHeightFactor");
+ const real_t floorHeightFactor = domainParameters.getParameter< real_t >("floorHeightFactor");
+ const real_t initialAmplitude = domainParameters.getParameter< real_t >("initialAmplitude");
+
+ // compute number of blocks as defined by domainSize and cellsPerBlock
+ Vector3< uint_t > numBlocks;
+ for (uint i = 0; i <= 2; ++i)
+ {
+ numBlocks[i] = domainSize[i] / cellsPerBlock[i];
+ WALBERLA_CHECK_EQUAL(numBlocks[i] * cellsPerBlock[i], domainSize[i],
+ "Domain size in direction " << i << " is not a multiple of cells per block.")
+ }
+
+ // get number of (MPI) processes
+ const uint_t numProcesses = uint_c(MPIManager::instance()->numProcesses());
+ WALBERLA_CHECK_EQUAL(numProcesses, numBlocks[0] * numBlocks[1] * numBlocks[2],
+ "The number of MPI processes is different from the number of blocks as defined by "
+ "\"domainSize/cellsPerBlock\".")
+
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numProcesses);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(numBlocks);
+
+ // get PID controller parameters
+ auto PIDParameters = configPtr->getOneBlock("PIDParameters");
+ const real_t targetMeanVelocityMagnitude = PIDParameters.getParameter< real_t >("targetMeanVelocityMagnitude");
+ const real_t proportionalGain = PIDParameters.getParameter< real_t >("proportionalGain");
+ const real_t derivativeGain = PIDParameters.getParameter< real_t >("derivativeGain");
+ const real_t integralGain = PIDParameters.getParameter< real_t >("integralGain");
+ const real_t maxRamp = PIDParameters.getParameter< real_t >("maxRamp");
+ const real_t minActuatingVariable = PIDParameters.getParameter< real_t >("minActuatingVariable");
+ const real_t maxActuatingVariable = PIDParameters.getParameter< real_t >("maxActuatingVariable");
+
+ // read particle infos
+ const auto particleParameters = configPtr->getOneBlock("ParticleParameters");
+ const std::string particleInFileName = particleParameters.getParameter< std::string >("inFileName");
+ const uint_t bedCopiesInX = particleParameters.getParameter< uint_t >("bedCopiesInX");
+ const uint_t bedCopiesInY = particleParameters.getParameter< uint_t >("bedCopiesInY");
+ const real_t particleDensityRatio = particleParameters.getParameter< real_t >("densityRatio");
+ const real_t particleFixingHeightFactor = particleParameters.getParameter< real_t >("fixingHeightFactor");
+ const real_t particleFrictionCoefficient = particleParameters.getParameter< real_t >("frictionCoefficient");
+ const real_t particleRestitutionCoefficient = particleParameters.getParameter< real_t >("restitutionCoefficient");
+ const uint_t particleNumSubCycles = particleParameters.getParameter< uint_t >("numSubCycles");
+ const bool useLubricationCorrection = particleParameters.getParameter< bool >("useLubricationCorrection");
+ const bool useNoSlipParticles = particleParameters.getParameter< bool >("useNoSlipParticles");
+ const real_t particlePoissonsRatio = 0.22_r;
+ const real_t particleKappa = real_t(2) * (real_t(1) - particlePoissonsRatio) / (real_t(2) - particlePoissonsRatio);
+ real_t particleCollisionTimeNonDim = 4_r;
+ bool useOpenMP = false;
+ const uint_t vtkSpacingParticles =
+ configPtr->getOneBlock("VTK").getOneBlock("fluid_field").getParameter< uint_t >("writeFrequency");
+ const std::string vtkFolder =
+ configPtr->getOneBlock("VTK").getOneBlock("fluid_field").getParameter< std::string >("baseFolder");
+
+ // get physics parameters from parameter file
+ auto physicsParameters = configPtr->getOneBlock("PhysicsParameters");
+ const bool enableWetting = physicsParameters.getParameter< bool >("enableWetting");
+ const uint_t timesteps = physicsParameters.getParameter< uint_t >("timesteps");
+ const real_t Re = physicsParameters.getParameter< real_t >("Re");
+ const real_t Fr = physicsParameters.getParameter< real_t >("Fr");
+ const real_t We = physicsParameters.getParameter< real_t >("We");
+
+ // get avgDiameter and scaling factor
+ real_t avgParticleDiameter = 0_r;
+ real_t particleScalingFactor = 0_r;
+ getAvgDiameterScalingFactor(particleInFileName, domainSize, bedCopiesInX, bedCopiesInY, avgParticleDiameter,
+ particleScalingFactor);
+ const real_t liquidHeight = avgParticleDiameter * liquidHeightFactor;
+ const real_t floorHeight = avgParticleDiameter * floorHeightFactor;
+ const real_t particleFixingHeight = particleFixingHeightFactor * avgParticleDiameter;
+
+ WALBERLA_CHECK_FLOAT_UNEQUAL(liquidHeight, 0.0)
+ WALBERLA_CHECK_FLOAT_UNEQUAL(floorHeight, 0.0)
+
+ const real_t absoluteLiquidHeight = liquidHeight + floorHeight;
+
+ const real_t viscosity = targetMeanVelocityMagnitude * liquidHeight / Re;
+ const real_t relaxationRate = real_t(1.0) / (real_t(3) * viscosity + real_t(0.5));
+
+ const real_t gravity = (targetMeanVelocityMagnitude / Fr) * (targetMeanVelocityMagnitude / Fr) / liquidHeight;
+ const real_t fx = real_t(3) * targetMeanVelocityMagnitude * viscosity / (liquidHeight * liquidHeight);
+ Vector3< real_t > force(fx, real_t(0.0), -gravity);
+
+ const real_t surfaceTension =
+ real_t(1.0) * targetMeanVelocityMagnitude * targetMeanVelocityMagnitude * liquidHeight / We;
+
+ // compute SI dx and dt
+ const real_t viscosity_SI = real_t(1.0016e-6); // kinemtic viscosity of water at 20°C at 1 bar
+ const real_t dx_SI = 1_r / particleScalingFactor;
+ const real_t dt_SI = viscosity / viscosity_SI * dx_SI * dx_SI;
+
+ WALBERLA_LOG_INFO_ON_ROOT("\nPhysical parameters:")
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(liquidHeight);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(floorHeight);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(dx_SI);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(dt_SI);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(force);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(relaxationRate);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(viscosity);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(absoluteLiquidHeight);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(avgParticleDiameter);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(particleScalingFactor);
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(particleFixingHeight);
+ WALBERLA_LOG_INFO_ON_ROOT("\nFree surface physical parameters")
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(surfaceTension);
+
+ if ((periodicity[0] && numBlocks[0] < uint_c(3)) || (periodicity[1] && numBlocks[1] < uint_c(3)) ||
+ (periodicity[2] && numBlocks[2] < uint_c(3)))
+ {
+ WALBERLA_ABORT("When using particles, use at least three blocks per periodic direction.");
+ }
+
+ // read model parameters from parameter file
+ const auto modelParameters = configPtr->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");
+
+ // read evaluation parameters from parameter file
+ const auto evaluationParameters = configPtr->getOneBlock("EvaluationParameters");
+ const uint_t performanceLogFrequency = evaluationParameters.getParameter< uint_t >("performanceLogFrequency");
+ const uint_t evaluationFrequency = evaluationParameters.getParameter< uint_t >("evaluationFrequency");
+ const std::string baseFolderName = evaluationParameters.getParameter< std::string >("baseFolderName");
+
+ uint_t beginTimeStep = 0;
+ const std::string checkpointConfigFile("antidunesCheckpointConfig.file");
+ if (loadSnapshot)
+ {
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ifstream file;
+ file.open(snapshotBaseFolder + "/" + checkpointConfigFile);
+ if (file.fail()) WALBERLA_ABORT("Error: " << checkpointConfigFile << " could not be opened!");
+ file >> beginTimeStep;
+ file >> force[0];
+ file.close();
+ }
+ mpi::broadcastObject(beginTimeStep);
+ mpi::broadcastObject(force);
+
+ WALBERLA_LOG_INFO_ON_ROOT("Successfully read config parameters from checkpoint config file:")
+ WALBERLA_LOG_INFO_ON_ROOT(" - beginTimeStep = " << beginTimeStep)
+ WALBERLA_LOG_INFO_ON_ROOT(" - force = < " << force[0] << ", " << force[1] << ", " << force[2] << " >")
+ }
+
+ if (loadSnapshot)
+ {
+ // modify config file to start VTK output from "loadFromTimestep" rather than from 0
+ std::vector< config::Config::Block* > configVTKBlock;
+ configPtr->getWritableGlobalBlock().getWritableBlocks("VTK", configVTKBlock, 1, 1);
+ std::vector< config::Config::Block* > configVTKFluidFieldBlock;
+ configVTKBlock[0]->getWritableBlocks("fluid_field", configVTKFluidFieldBlock, 1, 1);
+ configVTKFluidFieldBlock[0]->setOrAddParameter("initialExecutionCount", std::to_string(beginTimeStep));
+ }
+
+ WALBERLA_ROOT_SECTION()
+ {
+ // create base directories if they do not yet exist
+ filesystem::path tpath(baseFolderName);
+ if (!filesystem::exists(tpath)) filesystem::create_directory(tpath);
+
+ filesystem::path snapshotPath(snapshotBaseFolder);
+ if (!filesystem::exists(snapshotPath)) filesystem::create_directory(snapshotPath);
+ }
+
+ std::shared_ptr< StructuredBlockForest > blockForest(nullptr);
+ const std::string blockForestFile("blockForest.file");
+
+ if (loadSnapshot)
+ {
+ // load block forest from file
+ MPIManager::instance()->useWorldComm();
+
+ blockForest = make_shared< StructuredBlockForest >(
+ make_shared< BlockForest >(uint_c(MPIManager::instance()->rank()),
+ (std::string(snapshotBaseFolder + "/" + blockForestFile)).c_str(), true, false),
+ cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2]);
+ blockForest->createCellBoundingBoxes();
+ }
+ else
+ {
+ // create uniform block forest
+ blockForest = blockforest::createUniformBlockGrid(numBlocks[0], numBlocks[1], numBlocks[2], // blocks
+ cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2], // cells
+ real_c(1.0), // dx
+ true, // one block per process
+ periodicity[0], periodicity[1], periodicity[2]); // periodicity
+ }
+
+ // save block forest to file (but do not overwrite existing file if snapshot is loaded)
+ if (storeSnapshot && !loadSnapshot)
+ {
+ blockForest->getBlockForest().saveToFile(snapshotBaseFolder + "/" + blockForestFile);
+ }
+
+ const auto vtkParameters = configPtr->getOneBlock("VTK");
+ const auto vtkFluidParameters = vtkParameters.getOneBlock("fluid_field");
+
+ LatticeModel_T latticeModel = LatticeModel_T(force[0], force[1], force[2], relaxationRate);
+
+ BlockDataID pdfFieldID;
+ const std::string pdfFieldFile("pdfField.file");
+ BlockDataID fillFieldID;
+ const std::string fillFieldFile("fillField.file");
+
+ if (loadSnapshot)
+ {
+ // load PDF field from file
+ shared_ptr< lbm::internal::PdfFieldHandling< LatticeModel_T > > pdfFieldDataHandling =
+ make_shared< lbm::internal::PdfFieldHandling< LatticeModel_T > >(
+ blockForest, latticeModel, false, Vector3< real_t >(real_c(0)), real_c(1), uint_c(1), field::fzyx);
+ pdfFieldID = (blockForest->getBlockStorage())
+ .loadBlockData(snapshotBaseFolder + "/" + pdfFieldFile, pdfFieldDataHandling, "PDF field");
+
+ // load fill level field from file
+ std::shared_ptr< ScalarFieldHandling< ScalarField_T > > fillFieldDataHandling =
+ std::make_shared< ScalarFieldHandling< ScalarField_T > >(blockForest, uint_c(2));
+ fillFieldID =
+ (blockForest->getBlockStorage())
+ .loadBlockData(snapshotBaseFolder + "/" + fillFieldFile, fillFieldDataHandling, "Fill level field");
+ }
+ else
+ {
+ // add PDF field
+ pdfFieldID =
+ lbm::addPdfFieldToStorage(blockForest, "PDF field", latticeModel,
+ Vector3< real_t >(targetMeanVelocityMagnitude, 0, 0), real_t(1.0), field::fzyx);
+
+ // add fill level field (initialized with 0, i.e., gas everywhere)
+ fillFieldID =
+ field::addToStorage< ScalarField_T >(blockForest, "Fill level field", real_c(0.0), field::fzyx, uint_c(2));
+ }
+
+ // MesaPD data structures
+ auto particleStorage = std::make_shared< mesa_pd::data::ParticleStorage >(1);
+ auto particleAccessor = std::make_shared< mesa_pd::data::ParticleAccessorWithBaseShape >(particleStorage);
+ auto mesapdDomain = std::make_shared< mesa_pd::domain::BlockForestDomain >(blockForest->getBlockForestPointer());
+
+ BlockDataID particleStorageID;
+ const std::string particleStorageFile("particleStorageFile.file");
+ if (loadSnapshot)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Initializing particles from checkpointing file!");
+ particleStorageID = blockForest->loadBlockData(snapshotBaseFolder + "/" + particleStorageFile,
+ mesa_pd::domain::createBlockForestDataHandling(particleStorage),
+ "Particle Storage");
+ mesa_pd::mpi::ClearNextNeighborSync CNNS;
+ CNNS(*particleAccessor);
+
+ mesa_pd::mpi::ClearGhostOwnerSync CGOS;
+ CGOS(*particleAccessor);
+ }
+ else
+ {
+ particleStorageID =
+ blockForest->addBlockData(mesa_pd::domain::createBlockForestDataHandling(particleStorage), "Particle Storage");
+ }
+
+ BlockDataID particleFieldID = field::addToStorage< lbm_mesapd_coupling::ParticleField_T >(
+ blockForest, "Particle field", particleAccessor->getInvalidUid(), field::zyxf, uint_c(2));
+
+ auto densityReferenceHeight = absoluteLiquidHeight;
+ auto hydrostaticDensityFct = [force, densityReferenceHeight](const Vector3< real_t >& position) {
+ uint_t forceComponent = 2; // gravity is here strictly only acting in z direction!
+ return getHydrostaticDensity(position[forceComponent], densityReferenceHeight, force[forceComponent]);
+ };
+
+ // add boundary handling
+ const std::shared_ptr< AntidunesBoundaryHandling_T > antidunesBoundaryHandling =
+ std::make_shared< AntidunesBoundaryHandling_T >(blockForest, pdfFieldID, fillFieldID, particleFieldID,
+ particleAccessor, hydrostaticDensityFct);
+ const BlockDataID flagFieldID = antidunesBoundaryHandling->getFlagFieldID();
+ const typename AntidunesBoundaryHandling_T::FlagInfo_T& flagInfo = antidunesBoundaryHandling->getFlagInfo();
+
+ real_t sinusAmplitude = real_c(0.5) * initialAmplitude;
+ real_t sinusOffset = floorHeight;
+ real_t sinusWavelength = real_c(domainSize[0]) / real_c(wavePeriods);
+
+ if (!loadSnapshot)
+ {
+ // samples used in the Monte-Carlo like estimation of the fill level
+ const uint_t fillLevelInitSamples = uint_c(100); // actually there will be 101 since 0 is also included
+
+ const uint_t numTotalPoints = (fillLevelInitSamples + uint_c(1)) * (fillLevelInitSamples + uint_c(1));
+ const real_t stepsize = real_c(1) / real_c(fillLevelInitSamples);
+
+ // initialize free-surface sine profile
+ for (auto blockIt = blockForest->begin(); blockIt != blockForest->end(); ++blockIt)
+ {
+ ScalarField_T* const fillField = blockIt->getData< ScalarField_T >(fillFieldID);
+
+ WALBERLA_FOR_ALL_CELLS(fillFieldIt, fillField, {
+ // cell in block-local coordinates
+ const Cell localCell = fillFieldIt.cell();
+
+ // get cell in global coordinates
+ Cell globalCell = localCell;
+ blockForest->transformBlockLocalToGlobalCell(globalCell, *blockIt, localCell);
+
+ // Monte-Carlo like estimation of the fill level:
+ // create uniformly-distributed sample points in each cell and count the number of points below the sine
+ // profile; this fraction of points is used as the fill level to initialize the profile
+ uint_t numPointsBelow = uint_c(0);
+
+ for (uint_t xSample = uint_c(0); xSample <= fillLevelInitSamples; ++xSample)
+ {
+ // Pascal et al. (2021) defined the amplitude to span from minimum peak to maximum peak; in
+ // initializationProfile(), the amplitude is defined to range from the average to the maximum peak
+ const real_t functionValue =
+ initializationProfile(real_c(globalCell[0]) + real_c(xSample) * stepsize, sinusAmplitude,
+ absoluteLiquidHeight + real_c(0.5), sinusWavelength);
+
+ for (uint_t zSample = uint_c(0); zSample <= fillLevelInitSamples; ++zSample)
+ {
+ const real_t zPoint = real_c(globalCell[2]) + real_c(zSample) * stepsize;
+ // with operator <, a fill level of 1 can not be reached when the line is equal to the cell's top
+ // border; with operator <=, a fill level of 0 can not be reached when the line is equal to the cell's
+ // bottom border
+ if (zPoint < functionValue) { ++numPointsBelow; }
+ }
+ }
+
+ // fill level is fraction of points below sine profile
+ fillField->get(localCell) = real_c(numPointsBelow) / real_c(numTotalPoints);
+ }) // WALBERLA_FOR_ALL_CELLS
+ }
+
+ initializePoiseuilleProfile(*blockForest, pdfFieldID, fillFieldID, floorHeight, liquidHeight + real_c(0.5),
+ force[0], viscosity, sinusAmplitude, sinusWavelength);
+ }
+
+ // initialize domain boundary conditions from config file
+ const auto boundaryParameters = configPtr->getOneBlock("BoundaryParameters");
+ antidunesBoundaryHandling->initFromConfig(boundaryParameters);
+
+ std::function< void(void) > syncCall;
+ auto simulationDomainAABB = blockForest->getDomain();
+
+ lbm_mesapd_coupling::ParticleMappingKernel< AntidunesBoundaryHandling_T::BoundaryHandling_T > particleMappingKernel(
+ blockForest, antidunesBoundaryHandling->getHandlingID());
+ lbm_mesapd_coupling::MovingParticleMappingKernel< AntidunesBoundaryHandling_T::BoundaryHandling_T >
+ movingParticleMappingKernel(blockForest, antidunesBoundaryHandling->getHandlingID(), particleFieldID);
+
+ uint_t numParticles = 0;
+ // initialize bottom solid sine profile
+ if (!loadSnapshot)
+ {
+ auto createParticleFct = [sinusAmplitude, sinusOffset, sinusWavelength](Vector3< real_t > pos) {
+ return pos[2] < initializationProfile(pos[0], sinusAmplitude, sinusOffset, sinusWavelength);
+ };
+
+ real_t maxParticleHeight = 0_r;
+ initSpheresFromFile(particleInFileName, *particleStorage, *mesapdDomain, particleDensityRatio, domainSize,
+ createParticleFct, simulationDomainAABB, bedCopiesInX, bedCopiesInY, numParticles,
+ maxParticleHeight, particleScalingFactor);
+ WALBERLA_LOG_INFO_ON_ROOT("Max particle height " << maxParticleHeight);
+ if ((sinusOffset + sinusAmplitude) > maxParticleHeight)
+ WALBERLA_ABORT("Created particle bed is below desired sinus shape!");
+ if (2_r * sinusAmplitude > (maxParticleHeight - particleFixingHeight))
+ WALBERLA_ABORT("Created mobile particle bed is not high enough for desired sinus shape!");
+ if (useNoSlipParticles && (particleFixingHeight < maxParticleHeight))
+ WALBERLA_ABORT("You are using no-slip BCs on particles (which does not set hydrodynamic forces) but do not "
+ "fix all particles - this leads to wrong behavior and is not permitted!")
+
+ // fix lower particles
+ particleStorage->forEachParticle(
+ useOpenMP, mesa_pd::kernel::SelectAll(), *particleAccessor,
+ [particleFixingHeight](const size_t idx, auto& ac) {
+ if (ac.getPosition(idx)[2] < particleFixingHeight)
+ mesa_pd::data::particle_flags::set(ac.getFlagsRef(idx), mesa_pd::data::particle_flags::FIXED);
+ },
+ *particleAccessor);
+ }
+ else
+ {
+ real_t avgParticleDiameterTest = 0_r;
+ particleStorage->forEachParticle(
+ false, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ [&numParticles, &avgParticleDiameterTest](const size_t idx, auto& ac) {
+ auto sp = static_cast< mesa_pd::data::Sphere* >(ac.getBaseShape(idx).get());
+ ++numParticles;
+ avgParticleDiameterTest += 2_r * sp->getRadius();
+ },
+ *particleAccessor);
+ mpi::allReduceInplace(numParticles, mpi::SUM);
+ mpi::allReduceInplace(avgParticleDiameterTest, mpi::SUM);
+ avgParticleDiameterTest /= real_c(numParticles);
+ WALBERLA_LOG_INFO_ON_ROOT("Read particles from check pointing file with avg diameter of "
+ << avgParticleDiameterTest)
+ if (std::abs(avgParticleDiameterTest - avgParticleDiameter) / avgParticleDiameterTest > 0.05)
+ {
+ WALBERLA_ABORT("Particle diameters not correct.")
+ }
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("Created " << numParticles << " particles");
+
+ // create planes
+ createPlane(*particleStorage, simulationDomainAABB.minCorner(), Vector3< real_t >(0, 0, 1));
+ createPlane(*particleStorage, simulationDomainAABB.maxCorner(), Vector3< real_t >(0, 0, -1));
+
+ const real_t blockSyncExtension = real_t(2.5);
+ real_t maxPossibleParticleDiameter = avgParticleDiameter * 1.1_r;
+ if (maxPossibleParticleDiameter < 2_r * real_c(cellsPerBlock.min()) - blockSyncExtension)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Using next neighbor sync for particles");
+ syncCall = [particleStorage, mesapdDomain, blockSyncExtension]() {
+ mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+ syncNextNeighborFunc(*particleStorage, *mesapdDomain, blockSyncExtension);
+ };
+ syncCall();
+ }
+ else
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Using ghost owner sync for particles")
+ syncCall = [particleStorage, mesapdDomain, blockSyncExtension]() {
+ mesa_pd::mpi::SyncGhostOwners syncGhostOwnersFunc;
+ syncGhostOwnersFunc(*particleStorage, *mesapdDomain, blockSyncExtension);
+ };
+ for (uint_t i = 0; i < uint_c(std::ceil(maxPossibleParticleDiameter / real_c(cellsPerBlock.min()))); ++i)
+ syncCall();
+ }
+
+ if (useNoSlipParticles)
+ {
+ particleStorage->forEachParticle(useOpenMP, SphereSelector(), *particleAccessor, particleMappingKernel,
+ *particleAccessor, AntidunesBoundaryHandling_T::noSlipFlagID);
+ }
+ else
+ {
+ particleStorage->forEachParticle(useOpenMP, SphereSelector(), *particleAccessor, movingParticleMappingKernel,
+ *particleAccessor, AntidunesBoundaryHandling_T::movingObstacleFlagID);
+ }
+
+ // IMPORTANT REMARK: this must be only called after every solid flag has been set; otherwise, the boundary handling
+ // might not detect solid flags correctly
+ antidunesBoundaryHandling->initFlagsFromFillLevel();
+
+ // communication after initialization
+ Communication_T communication(blockForest, flagFieldID, fillFieldID);
+ communication();
+
+ PdfCommunication_T pdfCommunication(blockForest, pdfFieldID);
+ pdfCommunication();
+
+ // add bubble model
+ std::shared_ptr< walberla::free_surface::bubble_model::BubbleModelBase > bubbleModel =
+ std::make_shared< walberla::free_surface::bubble_model::BubbleModelConstantPressure >(real_c(1));
+
+ // initialize hydrostatic pressure
+ if (!loadSnapshot) { initHydrostaticPressure< PdfField_T >(blockForest, pdfFieldID, hydrostaticDensityFct); }
+
+ // set density in non-liquid or non-interface cells to 1 (after initializing with hydrostatic pressure)
+ // setDensityInNonFluidCellsToOne< FlagField_T, PdfField_T >(blockForest, flagInfo, flagFieldID, pdfFieldID);
+
+ // create timeloop
+ SweepTimeloop timeloop(blockForest, timesteps);
+ timeloop.setCurrentTimeStep(beginTimeStep);
+
+ timeloop.addFuncBeforeTimeStep(RemainingTimeLogger(timeloop.getNrOfTimeSteps()), "Remaining Time Logger");
+
+ // Laplace pressure = 2 * surface tension * curvature; curvature computation is not necessary with no surface
+ // tension
+ bool computeCurvature = false;
+ if (!realIsEqual(surfaceTension, real_c(0), real_c(1e-14))) { computeCurvature = true; }
+
+ auto blockStateUpdate = StateSweep(blockForest, flagInfo, flagFieldID);
+
+ // add surface geometry handler
+ BlockDataID curvatureFieldID =
+ field::addToStorage< ScalarField_T >(blockForest, "Curvature field", real_c(0), field::fzyx, uint_c(1));
+ BlockDataID normalFieldID = field::addToStorage< VectorField_T >(
+ blockForest, "Normal field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
+ BlockDataID obstacleNormalFieldID = field::addToStorage< VectorField_T >(
+ blockForest, "Obstacle normal field", Vector3< real_t >(real_c(0)), field::fzyx, uint_c(1));
+ // add field for smoothed fill levels
+ BlockDataID smoothFillFieldID =
+ field::addToStorage< ScalarField_T >(blockForest, "Smooth fill level field", real_c(0), field::fzyx, uint_c(1));
+
+ // smooth fill level field for decreasing error in finite difference normal and curvature computation (see
+ // dissertation of S. Bogner, 2017 (section 4.4.2.1))
+ walberla::free_surface::SmoothingSweep< CommunicationStencil_T, FlagField_T, ScalarField_T, VectorField_T >
+ smoothingSweep(smoothFillFieldID, fillFieldID, flagFieldID,
+ walberla::free_surface::flagIDs::liquidInterfaceGasFlagIDs, flagInfo.getObstacleIDSet(),
+ enableWetting);
+ // IMPORTANT REMARK: SmoothingSweep must be executed on all blocks, because the algorithm works on all liquid,
+ // interface and gas cells. This is necessary since the normals are not only computed in interface cells, but also
+ // in the neighborhood of interface cells. Therefore, meaningful values for the fill levels of the second
+ // neighbors of interface cells are also required in NormalSweep.
+ timeloop.add() << Sweep(smoothingSweep, "Sweep: fill level smoothing")
+ << AfterFunction(Communication_T(blockForest, smoothFillFieldID),
+ "Communication: after smoothing sweep");
+
+ // compute interface normals (using smoothed fill level field)
+ walberla::free_surface::NormalSweep< CommunicationStencil_T, FlagField_T, ScalarField_T, VectorField_T > normalSweep(
+ normalFieldID, smoothFillFieldID, flagFieldID, walberla::free_surface::flagIDs::interfaceFlagID,
+ walberla::free_surface::flagIDs::liquidInterfaceGasFlagIDs, flagInfo.getObstacleIDSet(), true, false, true,
+ false);
+ timeloop.add() << Sweep(normalSweep, "Sweep: normal computation", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep(), "Empty sweep: normal")
+ << AfterFunction(Communication_T(blockForest, normalFieldID), "Communication: after normal sweep");
+
+ if (computeCurvature)
+ {
+ // compute interface curvature using finite differences according to Brackbill et al.
+ walberla::free_surface::CurvatureSweepFiniteDifferences< CommunicationStencil_T, FlagField_T, ScalarField_T,
+ VectorField_T >
+ curvSweep(curvatureFieldID, normalFieldID, obstacleNormalFieldID, flagFieldID,
+ walberla::free_surface::flagIDs::interfaceFlagID,
+ walberla::free_surface::flagIDs::liquidInterfaceGasFlagIDs, flagInfo.getObstacleIDSet(), false,
+ real_c(0));
+ timeloop.add() << Sweep(curvSweep, "Sweep: curvature computation (finite difference method)",
+ StateSweep::fullFreeSurface)
+ << Sweep(emptySweep(), "Empty sweep: curvature")
+ << AfterFunction(Communication_T(blockForest, curvatureFieldID),
+ "Communication: after curvature sweep");
+ }
+
+ // add surface dynamics handler
+
+ // add standard waLBerla boundary handling
+ timeloop.add() << Sweep(antidunesBoundaryHandling->getBoundarySweep(), "Sweep: boundary handling",
+ Set< SUID >::emptySet(), StateSweep::onlyGasAndBoundary)
+ << Sweep(emptySweep(), "Empty sweep: boundary handling", StateSweep::onlyGasAndBoundary);
+
+ // 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 walberla::free_surface::StreamReconstructAdvectSweep<
+ LatticeModel_T, typename AntidunesBoundaryHandling_T::BoundaryHandling_T, FlagField_T,
+ typename AntidunesBoundaryHandling_T::FlagInfo_T, ScalarField_T, VectorField_T, true >
+ streamReconstructAdvectSweep(surfaceTension, antidunesBoundaryHandling->getHandlingID(), fillFieldID, flagFieldID,
+ pdfFieldID, normalFieldID, curvatureFieldID, flagInfo, bubbleModel.get(),
+ pdfReconstructionModel, useSimpleMassExchange, cellConversionThreshold,
+ cellConversionForceThreshold);
+ // sweep acts only on blocks with at least one interface cell (due to StateSweep::fullFreeSurface)
+ timeloop.add() << Sweep(streamReconstructAdvectSweep, "Sweep: StreamReconstructAdvect", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep(), "Empty sweep: StreamReconstructAdvect")
+ // 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(Communication_T(blockForest, fillFieldID, flagFieldID),
+ "Communication: after StreamReconstructAdvect sweep")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< ScalarField_T >(blockForest, fillFieldID),
+ "Second ghost layer update: after StreamReconstructAdvect sweep (fill level field)")
+ << AfterFunction(blockforest::UpdateSecondGhostLayer< FlagField_T >(blockForest, flagFieldID),
+ "Second ghost layer update: after StreamReconstructAdvect sweep (flag field)");
+
+ auto lbmSweepGenerated = typename LatticeModel_T::Sweep(pdfFieldID);
+
+ // temporary class for being able to call the LBM collision with operator()
+ class CollideSweep
+ {
+ public:
+ CollideSweep(const typename LatticeModel_T::Sweep& sweep) : sweep_(sweep){};
+
+ void operator()(IBlock* const block, const uint_t numberOfGhostLayersToInclude = uint_t(0))
+ {
+ sweep_.collide(block, numberOfGhostLayersToInclude);
+ }
+
+ private:
+ typename LatticeModel_T::Sweep sweep_;
+ };
+
+ timeloop.add() << Sweep(CollideSweep(lbmSweepGenerated), "Sweep: collision (generated)", StateSweep::fullFreeSurface)
+ << Sweep(lbmSweepGenerated, "Sweep: streamCollide (generated)", StateSweep::onlyLBM)
+ << Sweep(emptySweep(), "Empty sweep: streamCollide (generated)")
+ << AfterFunction(PdfCommunication_T(blockForest, pdfFieldID),
+ "Communication: after streamCollide (generated)");
+
+ // 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 walberla::free_surface::CellConversionSweep< LatticeModel_T, AntidunesBoundaryHandling_T::BoundaryHandling_T,
+ ScalarField_T >
+ cellConvSweep(antidunesBoundaryHandling->getHandlingID(), pdfFieldID, flagInfo, bubbleModel.get());
+ timeloop.add() << Sweep(cellConvSweep, "Sweep: cell conversion", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep(), "Empty sweep: cell conversion")
+ //<< AfterFunction(PdfCommunication_T(blockForest, pdfFieldID),
+ //
+ // "Communication: after cell conversion sweep (PDF field)")
+ // communicate the flag field also in corner directions
+ << AfterFunction(Communication_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)");
+
+ // reinitialize PDFs, i.e., refill cells that were converted from gas to interface
+ // - when the flag "convertedFromGasToInterface" has been set (by CellConversionSweep)
+ // - according to the method specified with pdfRefillingModel_
+ const walberla::free_surface::EquilibriumRefillingSweep< LatticeModel_T, FlagField_T > equilibriumRefillingSweep(
+ pdfFieldID, flagFieldID, flagInfo, true);
+ timeloop.add() << Sweep(equilibriumRefillingSweep, "Sweep: EquilibriumRefilling", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep(), "Empty sweep: EquilibriumRefilling")
+ << AfterFunction(PdfCommunication_T(blockForest, pdfFieldID),
+ "Communication: after EquilibriumRefilling sweep");
+
+ // distribute excess mass:
+ // - excess mass: mass that is free after conversion from interface to gas/liquid cells
+ // - update the bubble model
+ // IMPORTANT REMARK: this sweep computes the mass via the density, i.e., the PDF field must be up-to-date and the
+ // PdfRefillingSweep must have been performed
+ const walberla::free_surface::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(Communication_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(&walberla::free_surface::bubble_model::BubbleModelBase::update, bubbleModel),
+ "Sweep: bubble model update");
+
+ // reset all flags that signal cell conversions (except "keepInterfaceForWettingFlag")
+ walberla::free_surface::ConversionFlagsResetSweep< FlagField_T > resetConversionFlagsSweep(flagFieldID, flagInfo);
+ timeloop.add() << Sweep(resetConversionFlagsSweep, "Sweep: conversion flag reset", StateSweep::fullFreeSurface)
+ << Sweep(emptySweep(), "Empty sweep: conversion flag reset")
+ << AfterFunction(Communication_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");
+
+ // add VTK output
+ walberla::free_surface::addVTKOutput< LatticeModel_T, AntidunesBoundaryHandling_T, PdfField_T, FlagField_T,
+ ScalarField_T, VectorField_T >(
+ blockForest, timeloop, configPtr, flagInfo, pdfFieldID, flagFieldID, fillFieldID, BlockDataID(), curvatureFieldID,
+ normalFieldID, obstacleNormalFieldID);
+
+ // add triangle mesh output of free surface
+ walberla::free_surface::SurfaceMeshWriter< ScalarField_T, FlagField_T > surfaceMeshWriter(
+ blockForest, fillFieldID, flagFieldID, walberla::free_surface::flagIDs::liquidInterfaceGasFlagIDs, real_c(0),
+ configPtr);
+ surfaceMeshWriter(); // write initial mesh
+ timeloop.addFuncAfterTimeStep(surfaceMeshWriter, "Writer: surface mesh");
+
+ if (vtkSpacingParticles != uint_t(0))
+ {
+ // particle field
+ auto particleFieldVTK =
+ vtk::createVTKOutput_BlockData(blockForest, "particle_field", vtkSpacingParticles, 0, false, vtkFolder);
+ auto cellBB_filterParameters = vtkFluidParameters.getOneBlock("CellBB_filter");
+ const Vector3< uint_t > cellBB_filterMin = cellBB_filterParameters.getParameter< Vector3< uint_t > >("min");
+ const Vector3< uint_t > cellBB_filterMax = cellBB_filterParameters.getParameter< Vector3< uint_t > >("max");
+ AABB sliceAABB(real_c(cellBB_filterMin[0]), real_c(cellBB_filterMin[1]), real_c(cellBB_filterMin[2]),
+ real_c(cellBB_filterMax[0] + uint_t(1)), real_c(cellBB_filterMax[1] + uint_t(1)),
+ real_c(cellBB_filterMax[2] + uint_t(1)));
+
+ particleFieldVTK->addCellInclusionFilter(vtk::AABBCellFilter(sliceAABB));
+ particleFieldVTK->addCellDataWriter(
+ make_shared< field::VTKWriter< GhostLayerField< walberla::id_t, 1 > > >(particleFieldID, "particleField"));
+ particleFieldVTK->setSamplingResolution(vtkFluidParameters.getParameter< real_t >("samplingResolution"));
+ timeloop.addFuncBeforeTimeStep(vtk::writeFiles(particleFieldVTK), "VTK (particle field data");
+ }
+
+ if (vtkSpacingParticles != uint_t(0))
+ {
+ // sphere
+ auto particleVtkOutput = make_shared< mesa_pd::vtk::ParticleVtkOutput >(particleStorage);
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleUid >("uid");
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleLinearVelocity >("velocity");
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleInteractionRadius >("radius");
+ // limit output to process-local spheres
+ particleVtkOutput->setParticleSelector([](const mesa_pd::data::ParticleStorage::iterator& pIt) {
+ using namespace walberla::mesa_pd::data::particle_flags;
+ return (pIt->getBaseShape()->getShapeType() == mesa_pd::data::Sphere::SHAPE_TYPE) &&
+ !isSet(pIt->getFlags(), GHOST);
+ });
+ auto particleVtkWriter =
+ vtk::createVTKOutput_PointData(particleVtkOutput, "particles", vtkSpacingParticles, vtkFolder,
+ std::string("simulation_step"), false, true, true, true, beginTimeStep);
+ timeloop.addFuncAfterTimeStep(vtk::writeFiles(particleVtkWriter), "VTK (sphere data)");
+ }
+
+ // add logging for computational performance
+ const lbm::PerformanceLogger< FlagField_T > performanceLogger(
+ blockForest, flagFieldID, walberla::free_surface::flagIDs::liquidInterfaceFlagIDs, performanceLogFrequency);
+ timeloop.addFuncAfterTimeStep(performanceLogger, "Evaluator: performance logging");
+
+ // LBM stability check
+ timeloop.addFuncAfterTimeStep(makeSharedFunctor(field::makeStabilityChecker< PdfField_T, FlagField_T >(
+ walberlaEnv.config(), blockForest, pdfFieldID, flagFieldID,
+ walberla::free_surface::flagIDs::liquidInterfaceFlagIDs)),
+ "LBM stability check");
+
+ // add sweep for evaluating the fluid's mean velocity
+ const std::shared_ptr< Vector3< real_t > > meanVelocity = std::make_shared< Vector3< real_t > >(real_c(0));
+ const real_t velocityAveragingFactor = 1_r / (liquidHeight * real_c(domainSize[0]) * real_c(domainSize[1]));
+ MeanVelocityComputer< AntidunesBoundaryHandling_T, PdfField_T, FlagField_T > meanVelocityComputer(
+ blockForest, antidunesBoundaryHandling, pdfFieldID, meanVelocity, velocityAveragingFactor);
+
+ // PID Controller
+ shared_ptr< ForcingAdjuster > forcingAdjuster =
+ make_shared< ForcingAdjuster >(blockForest, pdfFieldID, targetMeanVelocityMagnitude, force[0], proportionalGain,
+ derivativeGain, integralGain, maxRamp, minActuatingVariable, maxActuatingVariable);
+
+ if (loadSnapshot) { forcingAdjuster->loadPIDSnapshot(snapshotBaseFolder + "/" + "pidState.file"); }
+
+ WcTimingPool timingPool;
+
+ // this is carried out after the particle integration, it corrects the flag field and restores missing PDF
+ // information then, the checkpointing file can be written, as otherwise some cells are invalid and can not be
+ // recovered
+ SweepTimeloop timeloopAfterParticles(blockForest, timesteps);
+ timeloopAfterParticles.setCurrentTimeStep(beginTimeStep);
+
+ // sweep for updating the particle mapping into the LBM simulation
+ bool strictlyConserveMomentum = false;
+ timeloopAfterParticles.add() << Sweep(
+ lbm_mesapd_coupling::makeMovingParticleMapping< PdfField_T, AntidunesBoundaryHandling_T::BoundaryHandling_T >(
+ blockForest, pdfFieldID, antidunesBoundaryHandling->getHandlingID(), particleFieldID, particleAccessor,
+ AntidunesBoundaryHandling_T::movingObstacleFlagID, FormerMO_Flag,
+ lbm_mesapd_coupling::RegularParticlesSelector(), strictlyConserveMomentum),
+ "Particle Mapping");
+
+ // sweep for restoring PDFs in cells previously occupied by particles
+ bool reconstruction_recomputeTargetDensity = false;
+ bool reconstruction_useCentralDifferences = true;
+ auto gradReconstructor =
+ lbm_mesapd_coupling::makeGradsMomentApproximationReconstructor< AntidunesBoundaryHandling_T::BoundaryHandling_T >(
+ blockForest, antidunesBoundaryHandling->getHandlingID(), relaxationRate, reconstruction_recomputeTargetDensity,
+ reconstruction_useCentralDifferences);
+
+ timeloopAfterParticles.add()
+ << Sweep(makeSharedSweep(
+ lbm_mesapd_coupling::makePdfReconstructionManager< PdfField_T,
+ AntidunesBoundaryHandling_T::BoundaryHandling_T >(
+ blockForest, pdfFieldID, antidunesBoundaryHandling->getHandlingID(), particleFieldID,
+ particleAccessor, FormerMO_Flag, walberla::free_surface::flagIDs::liquidFlagID, gradReconstructor,
+ strictlyConserveMomentum)),
+ "PDF Restore")
+ << AfterFunction(Communication_T(blockForest, flagFieldID, particleFieldID),
+ "Communication: after PDF reconstruction sweep") // unsure if necessary but added for consistency
+ << AfterFunction(pdfCommunication, "PDF Communication");
+
+ real_t timeStepSizeParticles = real_t(1) / real_t(particleNumSubCycles);
+ mesa_pd::kernel::VelocityVerletPreForceUpdate vvIntegratorPreForce(timeStepSizeParticles);
+ mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(timeStepSizeParticles);
+ mesa_pd::kernel::LinearSpringDashpot collisionResponse(1);
+ collisionResponse.setFrictionCoefficientDynamic(0, 0, particleFrictionCoefficient);
+ mesa_pd::mpi::ReduceProperty reduceProperty;
+ mesa_pd::mpi::ReduceContactHistory reduceAndSwapContactHistory;
+ lbm_mesapd_coupling::ResetHydrodynamicForceTorqueKernel resetHydrodynamicForceTorque;
+ lbm_mesapd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
+ real_t particleCollisionTime = particleCollisionTimeNonDim * avgParticleDiameter;
+ lbm_mesapd_coupling::LubricationCorrectionKernel lubricationCorrectionKernel(
+ viscosity, [](real_t r) { return (real_t(0.001 + real_t(0.00007) * r)) * r; });
+
+ WALBERLA_LOG_INFO_ON_ROOT("Will use particle time step size of "
+ << timeStepSizeParticles << " and collision time of " << particleCollisionTime);
+
+ AverageDataSliceEvaluator< PdfField_T, AntidunesBoundaryHandling_T, FlagField_T, ScalarField_T >
+ averageDataSliceEvaluator(blockForest, flagFieldID, fillFieldID, pdfFieldID);
+
+ std::shared_ptr< real_t > totalFluidMass = std::make_shared< real_t >(real_c(0));
+ walberla::free_surface::TotalMassComputer< AntidunesBoundaryHandling_T, PdfField_T, FlagField_T, ScalarField_T >
+ totalFluidMassEvaluator(blockForest, antidunesBoundaryHandling, pdfFieldID, fillFieldID, evaluationFrequency,
+ totalFluidMass);
+
+ BedloadTransportEvaluator< ParticleAccessor_T > bedloadTransportEvaluator(
+ particleAccessor, 1_r / real_c(domainSize[0] * domainSize[1]), numParticles);
+ auto bedLoadTransportFileName = baseFolderName + "/bedload.txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing bedload info to file " << bedLoadTransportFileName);
+
+ auto fluidInfoFileName = baseFolderName + "/fluidInfo.txt";
+ WALBERLA_LOG_INFO_ON_ROOT("Writing fluid info to file " << fluidInfoFileName);
+
+ // write info file
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ofstream evalInfoFile(baseFolderName + "/info.txt");
+ evalInfoFile << evaluationFrequency << "\n";
+ evalInfoFile << gravity << "\n";
+ evalInfoFile << viscosity << "\n";
+ evalInfoFile << particleDensityRatio << "\n";
+ evalInfoFile << avgParticleDiameter << "\n";
+ evalInfoFile << domainSize[0] << "\n";
+ evalInfoFile << domainSize[1] << "\n";
+ evalInfoFile << domainSize[2] << "\n";
+ evalInfoFile << numParticles << "\n";
+ evalInfoFile << dx_SI << "\n";
+ evalInfoFile << dt_SI << "\n";
+ evalInfoFile.close();
+ }
+
+ Vector3< real_t > totalHydrodynamicForceOnParticles(0_r); // only root will have valid values
+
+ for (uint_t t = beginTimeStep; t != timesteps; ++t)
+ {
+ timeloop.singleStep(timingPool, true);
+
+ timingPool["Mesa_pd"].start();
+
+ reduceProperty.operator()< mesa_pd::HydrodynamicForceTorqueNotification >(*particleStorage);
+
+ if (t == 0)
+ {
+ lbm_mesapd_coupling::InitializeHydrodynamicForceTorqueForAveragingKernel
+ initializeHydrodynamicForceTorqueForAveragingKernel;
+ particleStorage->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ initializeHydrodynamicForceTorqueForAveragingKernel, *particleAccessor);
+ }
+ particleStorage->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ averageHydrodynamicForceTorque, *particleAccessor);
+
+ for (auto subCycle = uint_t(0); subCycle < particleNumSubCycles; ++subCycle)
+ {
+ particleStorage->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ vvIntegratorPreForce, *particleAccessor);
+ syncCall();
+
+ if (useLubricationCorrection)
+ {
+ // lubrication correction
+ particleStorage->forEachParticlePairHalf(
+ useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *particleAccessor,
+ [&lubricationCorrectionKernel, &mesapdDomain](const size_t idx1, const size_t idx2, auto& ac) {
+ mesa_pd::collision_detection::AnalyticContactDetection acd;
+ acd.getContactThreshold() = lubricationCorrectionKernel.getNormalCutOffDistance();
+ mesa_pd::kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *mesapdDomain))
+ {
+ double_cast(acd.getIdx1(), acd.getIdx2(), ac, lubricationCorrectionKernel, ac,
+ acd.getContactNormal(), acd.getPenetrationDepth());
+ }
+ }
+ },
+ *particleAccessor);
+ }
+
+ // collision response
+ particleStorage->forEachParticlePairHalf(
+ useOpenMP, mesa_pd::kernel::ExcludeInfiniteInfinite(), *particleAccessor,
+ [&collisionResponse, &mesapdDomain, timeStepSizeParticles, particleRestitutionCoefficient,
+ particleCollisionTime, particleKappa](const size_t idx1, const size_t idx2, auto& ac) {
+ mesa_pd::collision_detection::AnalyticContactDetection acd;
+ mesa_pd::kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *mesapdDomain))
+ {
+ auto meff = real_t(1) / (ac.getInvMass(idx1) + ac.getInvMass(idx2));
+ collisionResponse.setStiffnessAndDamping(0, 0, particleRestitutionCoefficient,
+ particleCollisionTime, particleKappa, meff);
+ collisionResponse(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(),
+ acd.getPenetrationDepth(), timeStepSizeParticles);
+ }
+ }
+ },
+ *particleAccessor);
+
+ reduceAndSwapContactHistory(*particleStorage);
+
+ // add hydrodynamic force
+ lbm_mesapd_coupling::AddHydrodynamicInteractionKernel addHydrodynamicInteraction;
+ particleStorage->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ addHydrodynamicInteraction, *particleAccessor);
+
+ // add external forces
+ particleStorage->forEachParticle(
+ useOpenMP, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ [particleDensityRatio, force](const size_t idx, auto& ac) {
+ mesa_pd::addForceAtomic(
+ idx, ac,
+ ac.getVolume(idx) *
+ Vector3< real_t >(force[0], force[1], (particleDensityRatio - real_t(1)) * force[2]));
+ },
+ *particleAccessor);
+
+ reduceProperty.operator()< mesa_pd::ForceTorqueNotification >(*particleStorage);
+
+ particleStorage->forEachParticle(useOpenMP, mesa_pd::kernel::SelectLocal(), *particleAccessor,
+ vvIntegratorPostForce, *particleAccessor);
+ syncCall();
+ }
+
+ // has to be evaluated here before the force info is erased from particles
+ if (t % evaluationFrequency == uint_c(0))
+ totalHydrodynamicForceOnParticles = getTotalHydrodynamicForceOnParticles(particleAccessor);
+
+ particleStorage->forEachParticle(useOpenMP, mesa_pd::kernel::SelectAll(), *particleAccessor,
+ resetHydrodynamicForceTorque, *particleAccessor);
+ timingPool["Mesa_pd"].end();
+
+ // update particle mapping
+ timeloopAfterParticles.singleStep(timingPool, true);
+
+ timingPool["Evaluation"].start();
+
+ if (t % evaluationFrequency == uint_c(0))
+ {
+ averageDataSliceEvaluator();
+ totalFluidMassEvaluator.computeMass(blockForest, antidunesBoundaryHandling);
+ bedloadTransportEvaluator();
+ meanVelocityComputer();
+
+ WALBERLA_ROOT_SECTION()
+ {
+ write2DVectorToFile(averageDataSliceEvaluator.getSolidVolumeFractionVector(),
+ averageDataSliceEvaluator.getXLen(), averageDataSliceEvaluator.getZLen(),
+ baseFolderName + "/svfSlice_" + std::to_string(t) + ".txt");
+ write2DVectorToFile(averageDataSliceEvaluator.getFillLevelVector(), averageDataSliceEvaluator.getXLen(),
+ averageDataSliceEvaluator.getZLen(),
+ baseFolderName + "/fillSlice_" + std::to_string(t) + ".txt");
+ write2DVectorToFile(averageDataSliceEvaluator.getVelocityXVector(), averageDataSliceEvaluator.getXLen(),
+ averageDataSliceEvaluator.getZLen(),
+ baseFolderName + "/velXSlice_" + std::to_string(t) + ".txt");
+
+ std::ofstream bedloadFile(bedLoadTransportFileName, std::ofstream::app);
+ bedloadFile << t << " " << bedloadTransportEvaluator.getTransportRate() << " "
+ << bedloadTransportEvaluator.getAverageVelocity() << " " << totalHydrodynamicForceOnParticles[0]
+ << " " << totalHydrodynamicForceOnParticles[1] << " " << totalHydrodynamicForceOnParticles[2]
+ << "\n";
+ bedloadFile.close();
+
+ WALBERLA_LOG_DEVEL("____________________________________________________________________");
+ WALBERLA_LOG_DEVEL("time step = " << t);
+ const real_t froudeNumber = (*meanVelocity)[0] / real_c(std::sqrt(liquidHeight * std::abs(force[2])));
+
+ const real_t reynoldsNumber = (*meanVelocity)[0] * liquidHeight / viscosity;
+
+ const real_t weberNumber =
+ real_t(1.0) * (*meanVelocity)[0] * (*meanVelocity)[0] * liquidHeight / surfaceTension;
+
+ WALBERLA_LOG_DEVEL(" - Total fluid mass = " << std::setprecision(16) << (*totalFluidMass));
+ auto maxFluidZPos = averageDataSliceEvaluator.getMaxFluidZPos();
+ WALBERLA_LOG_DEVEL(" - Max fluid z-position = " << maxFluidZPos);
+ WALBERLA_LOG_DEVEL(" - Froude number = " << froudeNumber);
+ WALBERLA_LOG_DEVEL(" - Reynolds number = " << reynoldsNumber);
+ WALBERLA_LOG_DEVEL(" - We = " << weberNumber);
+
+ WALBERLA_LOG_DEVEL(" - meanVelocity = " << *meanVelocity);
+
+ std::ofstream fluidInfoFile(fluidInfoFileName, std::ofstream::app);
+ fluidInfoFile << t << " " << force[0] << " " << (*meanVelocity)[0] << " " << maxFluidZPos << " "
+ << std::setprecision(16) << (*totalFluidMass) << "\n";
+ fluidInfoFile.close();
+
+ if (std::isnan(reynoldsNumber)) WALBERLA_ABORT("reynoldsNumber is inf!")
+ }
+
+ WALBERLA_LOG_DEVEL_ON_ROOT(" -> CurrentExternalForce in x-direction before update = " << force[0]);
+ (*forcingAdjuster)(meanVelocity->length());
+ force[0] = forcingAdjuster->getExternalForcing();
+ WALBERLA_LOG_DEVEL_ON_ROOT(" -> CurrentExternalForce in x-direction after update = " << force[0]);
+ }
+ timingPool["Evaluation"].end();
+
+ if (storeSnapshot)
+ {
+ if (t % snapshotFrequency == uint_c(0) && t > uint_c(0))
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Writing checkpointing file in time step " << t)
+
+ blockForest->saveBlockData(snapshotBaseFolder + "/tmp_" + pdfFieldFile, pdfFieldID);
+ blockForest->saveBlockData(snapshotBaseFolder + "/tmp_" + fillFieldFile, fillFieldID);
+ blockForest->saveBlockData(snapshotBaseFolder + "/tmp_" + particleStorageFile, particleStorageID);
+
+ WALBERLA_ROOT_SECTION()
+ {
+ std::string tmpCheckpointConfigFile = snapshotBaseFolder + "/tmp_" + checkpointConfigFile;
+ std::ofstream file;
+ file.open(tmpCheckpointConfigFile.c_str());
+
+ file << std::setprecision(16);
+ file << t + 1 << "\n";
+ file << force[0] << "\n";
+ file.close();
+ }
+
+ forcingAdjuster->storePIDSnapshot(snapshotBaseFolder + "/" + "pidState.file");
+
+ WALBERLA_MPI_BARRIER();
+
+ // rename checkpoint files to "real" ones
+ // otherwise, the checkpointed state might be incomplete if the simulation stops due to over time during
+ // checkpointing
+ WALBERLA_ROOT_SECTION()
+ {
+ renameFile(snapshotBaseFolder + "/tmp_" + pdfFieldFile, snapshotBaseFolder + "/" + pdfFieldFile);
+ renameFile(snapshotBaseFolder + "/tmp_" + fillFieldFile, snapshotBaseFolder + "/" + fillFieldFile);
+ renameFile(snapshotBaseFolder + "/tmp_" + particleStorageFile,
+ snapshotBaseFolder + "/" + particleStorageFile);
+ renameFile(snapshotBaseFolder + "/tmp_" + checkpointConfigFile,
+ snapshotBaseFolder + "/" + checkpointConfigFile);
+ }
+ }
+ }
+
+ if (t % performanceLogFrequency == uint_c(0) && t > uint_c(0)) { timingPool.logResultOnRoot(); }
+ }
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace antidunes
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::antidunes::main(argc, argv); }
diff --git a/apps/showcases/Antidunes/Antidunes.prm b/apps/showcases/Antidunes/Antidunes.prm
new file mode 100644
index 0000000000000000000000000000000000000000..306cb548ec21423211c3586c41fc9ac8d0a738d1
--- /dev/null
+++ b/apps/showcases/Antidunes/Antidunes.prm
@@ -0,0 +1,147 @@
+BlockForestParameters
+{
+ cellsPerBlock < 50, 20, 40 >;
+ periodicity < 1, 1, 0 >;
+ loadSnapshot false;
+ storeSnapshot true;
+ snapshotFrequency 10000;
+ snapshotBaseFolder snapshots;
+}
+
+DomainParameters
+{
+ domainSize <3200, 60, 160>;
+ wavePeriods 1; // never set to 0 -> division by zero, even if you initialize a flat particle bed
+ liquidHeightFactor 2.9655; // h_0 / d (water height / avg particle diameter) -> from experiment [E1 = 2.9655, E4 = 3.5862]
+ floorHeightFactor 4.1393; // from domain bottom to sine's average
+ initialAmplitude 0; // defined from minimum peak to maximum peak as by Pascal et al. (2021)
+}
+
+PIDParameters
+{
+ targetMeanVelocityMagnitude 0.02;
+ proportionalGain 2e-4;
+ derivativeGain 1e-6;
+ integralGain 2e-4;
+ maxRamp 1e-4;
+ minActuatingVariable 0;
+ maxActuatingVariable 1e-3;
+}
+
+PhysicsParameters
+{
+ enableWetting false;
+ timesteps 2000000;
+ Re 3100; // [E1=3100, E4=4800]
+ Fr 1.31; // [E1=1.31, E4=1.45]
+ We 15.6188; // [E1=15.6188 , E4=30.2493]
+}
+
+ParticleParameters
+{
+ inFileName spheres_out.dat;
+ bedCopiesInX 1;
+ bedCopiesInY 1;
+ densityRatio 2.55;
+ fixingHeightFactor 1.5; // proportional to the mean particle diameter
+ frictionCoefficient 0.5;
+ restitutionCoefficient 0.97;
+ numSubCycles 10;
+ useLubricationCorrection true;
+ useNoSlipParticles false;
+}
+
+ModelParameters
+{
+ pdfReconstructionModel OnlyMissing;
+ pdfRefillingModel EquilibriumRefilling;
+ excessMassDistributionModel EvenlyAllInterface;
+ curvatureModel FiniteDifferenceMethod;
+ enableForceWeighting false;
+ useSimpleMassExchange false;
+ cellConversionThreshold 1e-2;
+ cellConversionForceThreshold 1e-1;
+}
+
+EvaluationParameters
+{
+ performanceLogFrequency 10000;
+ evaluationFrequency 1000;
+ baseFolderName eval;
+}
+
+StabilityChecker
+{
+ checkFrequency 0;
+ streamOutput false;
+ vtkOutput true;
+}
+
+BoundaryParameters
+{
+ // X
+ //Border { direction W; walldistance -1; NoSlip{} }
+ //Border { direction E; walldistance -1; NoSlip{} }
+
+ // Y
+ //Border { direction N; walldistance -1; NoSlip{} }
+ //Border { direction S; walldistance -1; NoSlip{} }
+
+ // Z
+ Border { direction T; walldistance -1; NoSlip{} }
+ Border { direction B; walldistance -1; NoSlip{} }
+}
+
+MeshOutputParameters
+{
+ writeFrequency 0;
+ baseFolder mesh-out;
+}
+
+VTK
+{
+ fluid_field
+ {
+ writeFrequency 1000;
+ ghostLayers 0;
+ baseFolder vtk-out;
+ samplingResolution 1;
+
+ writers
+ {
+ fill_level;
+ mapped_flag;
+ velocity;
+ density;
+ //curvature;
+ //normal;
+ //obstacle_normal;
+ //pdf;
+ //flag;
+ //force;
+ }
+
+ CellBB_filter {
+ min < 0, 29, 0 >;
+ max < 3200, 30, 160 >;
+ }
+
+ inclusion_filters
+ {
+ CellBB_filter;
+ //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 10000000000;
+ baseFolder vtk-out;
+ outputDomainDecomposition true;
+ }
+}
diff --git a/apps/showcases/Antidunes/AntidunesBoundaryHandling.h b/apps/showcases/Antidunes/AntidunesBoundaryHandling.h
new file mode 100644
index 0000000000000000000000000000000000000000..f210d8e74f3a638ddfa691f3c1cac05bba14acd8
--- /dev/null
+++ b/apps/showcases/Antidunes/AntidunesBoundaryHandling.h
@@ -0,0 +1,207 @@
+//======================================================================================================================
+//
+// 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 Boundary.h
+//! \ingroup free_surface
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Boundary handling for the free surface LBM module.
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/StructuredBlockForest.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+
+#include "boundary/BoundaryHandling.h"
+
+#include "field/GhostLayerField.h"
+
+#include "geometry/initializer/InitializationManager.h"
+#include "geometry/initializer/OverlapFieldFromBody.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/free_surface/FlagInfo.h"
+#include "lbm/free_surface/InitFunctions.h"
+#include "lbm/free_surface/InterfaceFromFillLevel.h"
+#include "lbm/free_surface/boundary/SimplePressureWithFreeSurface.h"
+
+#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h"
+
+namespace walberla
+{
+namespace antidunes
+{
+namespace free_surface
+{
+/***********************************************************************************************************************
+ * Boundary handling for the free surface LBM extension.
+ **********************************************************************************************************************/
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+class AntidunesBoundaryHandling
+{
+ public:
+ using flag_t = typename FlagField_T::value_type;
+ using Stencil_T = typename LatticeModel_T::Stencil;
+ using CommunicationStencil_T =
+ typename std::conditional< LatticeModel_T::Stencil::D == uint_t(2), stencil::D2Q9, stencil::D3Q27 >::type;
+ using PdfField_T = lbm::PdfField< LatticeModel_T >;
+
+ // boundary
+ using NoSlip_T = lbm::NoSlip< LatticeModel_T, flag_t >;
+ using FreeSlip_T = lbm::FreeSlip< LatticeModel_T, FlagField_T >;
+ using UBB_T = lbm::UBB< LatticeModel_T, flag_t >;
+ using Pressure_T = walberla::free_surface::SimplePressureWithFreeSurface< LatticeModel_T, FlagField_T >;
+ using Outlet_T = lbm::Outlet< LatticeModel_T, FlagField_T, 4, 3 >;
+ using UBB_Inflow_T =
+ lbm::UBB< LatticeModel_T, flag_t >; // creates interface cells in the direction of the prescribed velocity, i.e.,
+ // represents an inflow boundary condition
+ using MovingObstacle_T = lbm_mesapd_coupling::SimpleBB< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
+
+ // handling type
+ using BoundaryHandling_T =
+ BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, UBB_T, UBB_Inflow_T, Pressure_T, Pressure_T, Outlet_T,
+ FreeSlip_T,
+ MovingObstacle_T >; // 2 pressure boundaries with different densities, e.g., inflow-outflow
+
+ using FlagInfo_T = walberla::free_surface::FlagInfo< FlagField_T >;
+
+ // constructor
+ AntidunesBoundaryHandling(const std::shared_ptr< StructuredBlockForest >& blockForest, BlockDataID pdfFieldID,
+ BlockDataID fillLevelID, BlockDataID particleFieldID,
+ const shared_ptr< ParticleAccessor_T >& ac,
+ std::function< real_t(const Vector3< real_t >&) > hydrostaticDensityFct);
+
+ // initialize fluid field from config file using (quotes indicate the string to be used in the file):
+ // - "CellInterval" (see src/geometry/initializer/BoundaryFromCellInterval.h)
+ // - "Border" (see src/geometry/initializer/BoundaryFromDomainBorder.h)
+ // - "Image" (see src/geometry/initializer/BoundaryFromImage.h)
+ // - "Body" (see src/geometry/initializer/OverlapFieldFromBody.h)
+ inline void initFromConfig(const Config::BlockHandle& block);
+
+ // initialize free surface object from geometric body (see src/geometry/initializer/OverlapFieldFromBody.h)
+ template< typename Body_T >
+ inline void addFreeSurfaceObject(const Body_T& body, bool addOrSubtract = false);
+
+ // clear and initialize flags in every cell according to the fill level and initialize fill level in boundaries (with
+ // value 1) and obstacles such that the bubble model does not detect obstacles as gas cells
+ void initFlagsFromFillLevel();
+
+ inline void setNoSlipAtBorder(stencil::Direction d, cell_idx_t wallDistance = cell_idx_c(0));
+ inline void setNoSlipAtAllBorders(cell_idx_t wallDistance = cell_idx_c(0));
+ void setNoSlipInCell(const Cell& globalCell);
+
+ inline void setFreeSlipAtBorder(stencil::Direction d, cell_idx_t wallDistance = cell_idx_c(0));
+ inline void setFreeSlipAtAllBorders(cell_idx_t wallDistance = cell_idx_c(0));
+ void setFreeSlipInCell(const Cell& globalCell);
+
+ void setUBBInCell(const Cell& globalCell, const Vector3< real_t >& velocity);
+
+ // UBB that generates interface cells to resemble an inflow boundary
+ void setInflowInCell(const Cell& globalCell, const Vector3< real_t >& velocity);
+
+ inline void setPressure(real_t density);
+ void setPressureOutflow(real_t density);
+ void setBodyForce(const Vector3< real_t >& bodyForce);
+
+ void enableBubbleOutflow(BubbleModelBase* bubbleModel);
+
+ // checks if an obstacle cell is located in an outermost ghost layer (corners are explicitly ignored, as they do not
+ // influence periodic communication)
+ Vector3< bool > isObstacleInGlobalGhostLayer();
+
+ // flag management
+ const walberla::free_surface::FlagInfo< FlagField_T >& getFlagInfo() const { return flagInfo_; }
+
+ // flag IDs
+ static const field::FlagUID noSlipFlagID;
+ static const field::FlagUID ubbFlagID;
+ static const field::FlagUID ubbInflowFlagID;
+ static const field::FlagUID pressureFlagID;
+ static const field::FlagUID pressureOutflowFlagID;
+ static const field::FlagUID outletFlagID;
+ static const field::FlagUID freeSlipFlagID;
+ static const field::FlagUID movingObstacleFlagID;
+
+ // boundary IDs
+ static const BoundaryUID noSlipBoundaryID;
+ static const BoundaryUID ubbBoundaryID;
+ static const BoundaryUID ubbInflowBoundaryID;
+ static const BoundaryUID pressureBoundaryID;
+ static const BoundaryUID pressureOutflowBoundaryID;
+ static const BoundaryUID outletBoundaryID;
+ static const BoundaryUID freeSlipBoundaryID;
+ static const BoundaryUID movingObstacleBoundaryID;
+
+ inline BlockDataID getHandlingID() const { return handlingID_; }
+ inline BlockDataID getPdfFieldID() const { return pdfFieldID_; }
+ inline BlockDataID getFillFieldID() const { return fillFieldID_; }
+ inline BlockDataID getFlagFieldID() const { return flagFieldID_; }
+ inline BlockDataID getParticleFieldID() const { return particleFieldID_; }
+
+ inline std::shared_ptr< StructuredBlockForest > getBlockForest() const { return blockForest_; };
+ inline shared_ptr< ParticleAccessor_T > getParticleAccessor() const { return ac_; }
+ inline std::function< real_t(const Vector3< real_t >&) > getHydrostaticDensityFct() const
+ {
+ return hydrostaticDensityFct_;
+ }
+
+ // executes standard waLBerla boundary handling
+ class ExecuteBoundaryHandling
+ {
+ public:
+ ExecuteBoundaryHandling(const BlockDataID& collection) : handlingID_(collection) {}
+ void operator()(IBlock* const block) const
+ {
+ BoundaryHandling_T* const handling = block->getData< BoundaryHandling_T >(handlingID_);
+ // reset "near boundary" flags
+ handling->refresh();
+ (*handling)();
+ }
+
+ protected:
+ BlockDataID handlingID_;
+ }; // class ExecuteBoundaryHandling
+
+ ExecuteBoundaryHandling getBoundarySweep() const { return ExecuteBoundaryHandling(getHandlingID()); }
+
+ private:
+ walberla::free_surface::FlagInfo< FlagField_T > flagInfo_;
+
+ // register standard waLBerla initializers
+ geometry::initializer::InitializationManager getInitManager();
+
+ std::shared_ptr< StructuredBlockForest > blockForest_;
+
+ BlockDataID flagFieldID_;
+ BlockDataID pdfFieldID_;
+ BlockDataID fillFieldID_;
+ BlockDataID particleFieldID_;
+ shared_ptr< ParticleAccessor_T > ac_;
+ std::function< real_t(const Vector3< real_t >&) > hydrostaticDensityFct_;
+
+ BlockDataID handlingID_;
+
+ blockforest::communication::UniformBufferedScheme< CommunicationStencil_T > comm_;
+}; // class AntidunesBoundaryHandling
+
+} // namespace free_surface
+} // namespace antidunes
+} // namespace walberla
+
+#include "AntidunesBoundaryHandling.impl.h"
diff --git a/apps/showcases/Antidunes/AntidunesBoundaryHandling.impl.h b/apps/showcases/Antidunes/AntidunesBoundaryHandling.impl.h
new file mode 100644
index 0000000000000000000000000000000000000000..3aa86f16a5e3d93bcf030e4235b898242767ea90
--- /dev/null
+++ b/apps/showcases/Antidunes/AntidunesBoundaryHandling.impl.h
@@ -0,0 +1,583 @@
+//======================================================================================================================
+//
+// 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 AntidunesBoundaryHandling.impl.h
+//! \ingroup free_surface
+//! \author Martin Bauer
+//! \author Christoph Schwarzmeier <christoph.schwarzmeier@fau.de>
+//! \brief Boundary handling for the free surface LBM module.
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "field/AddToStorage.h"
+#include "field/FlagField.h"
+#include "field/communication/PackInfo.h"
+
+#include "geometry/initializer/BoundaryFromCellInterval.h"
+#include "geometry/initializer/BoundaryFromDomainBorder.h"
+#include "geometry/initializer/BoundaryFromImage.h"
+#include "geometry/structured/GrayScaleImage.h"
+
+#include "lbm/free_surface/FlagInfo.h"
+#include "lbm/free_surface/InterfaceFromFillLevel.h"
+#include "lbm/lattice_model/CollisionModel.h"
+
+#include "AntidunesBoundaryHandling.h"
+
+namespace walberla
+{
+namespace antidunes
+{
+namespace free_surface
+{
+namespace internal
+{
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+class BoundaryBlockDataHandling
+ : public domain_decomposition::BlockDataHandling< typename AntidunesBoundaryHandling<
+ LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::BoundaryHandling_T >
+{
+ public:
+ using BoundaryHandling_T =
+ typename AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T,
+ ParticleAccessor_T >::BoundaryHandling_T; // handling as defined in
+ // AntidunesBoundaryHandling.h
+
+ BoundaryBlockDataHandling(
+ const AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >* boundary)
+ : boundary_(boundary)
+ {}
+
+ // initialize standard waLBerla boundary handling
+ BoundaryHandling_T* initialize(IBlock* const block)
+ {
+ using B = AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >;
+ using flag_t = typename B::flag_t;
+
+ // get fields
+ FlagField_T* const flagField = block->getData< FlagField_T >(boundary_->getFlagFieldID());
+ typename B::PdfField_T* const pdfField = block->getData< typename B::PdfField_T >(boundary_->getPdfFieldID());
+ lbm_mesapd_coupling::ParticleField_T* const particleField =
+ block->getData< lbm_mesapd_coupling::ParticleField_T >(boundary_->getParticleFieldID());
+
+ auto interfaceFlag = flag_t(flagField->getFlag(walberla::free_surface::flagIDs::interfaceFlagID));
+ auto liquidFlag = flag_t(flagField->getFlag(walberla::free_surface::flagIDs::liquidFlagID));
+
+ // domainMask is used to identify liquid and interface cells
+ auto domainMask = flag_t(liquidFlag | interfaceFlag);
+ WALBERLA_ASSERT(domainMask != 0);
+
+ // initialize boundary conditions
+ typename B::UBB_T ubb(B::ubbBoundaryID, B::ubbFlagID, pdfField, flagField);
+ typename B::UBB_Inflow_T ubbInflow(B::ubbInflowBoundaryID, B::ubbInflowFlagID, pdfField, flagField);
+ typename B::NoSlip_T noslip(B::noSlipBoundaryID, B::noSlipFlagID, pdfField);
+ typename B::Pressure_T pressure(B::pressureBoundaryID, B::pressureFlagID, block, pdfField, flagField,
+ interfaceFlag, real_c(1.0));
+ typename B::Pressure_T pressureOutflow(B::pressureOutflowBoundaryID, B::pressureOutflowFlagID, block, pdfField,
+ flagField, interfaceFlag, real_c(1.0));
+ typename B::Outlet_T outlet(B::outletBoundaryID, B::outletFlagID, pdfField, flagField, domainMask);
+ typename B::FreeSlip_T freeSlip(B::freeSlipBoundaryID, B::freeSlipFlagID, pdfField, flagField, domainMask);
+ typename B::MovingObstacle_T curvedLinear(B::movingObstacleBoundaryID, B::movingObstacleFlagID, pdfField,
+ flagField, particleField, boundary_->getParticleAccessor(), domainMask,
+ *boundary_->getBlockForest(), *block,
+ boundary_->getHydrostaticDensityFct());
+
+ return new BoundaryHandling_T("Boundary Handling", flagField, domainMask, noslip, ubb, ubbInflow, pressure,
+ pressureOutflow, outlet, freeSlip, curvedLinear);
+ }
+
+ void serialize(IBlock* const block, const BlockDataID& id, mpi::SendBuffer& buffer)
+ {
+ BoundaryHandling_T* const boundaryHandlingPtr = block->getData< BoundaryHandling_T >(id);
+ CellInterval everyCell = boundaryHandlingPtr->getFlagField()->xyzSizeWithGhostLayer();
+ boundaryHandlingPtr->pack(buffer, everyCell, true);
+ }
+
+ BoundaryHandling_T* deserialize(IBlock* const block) { return initialize(block); }
+
+ void deserialize(IBlock* const block, const BlockDataID& id, mpi::RecvBuffer& buffer)
+ {
+ BoundaryHandling_T* const boundaryHandlingPtr = block->getData< BoundaryHandling_T >(id);
+ CellInterval everyCell = boundaryHandlingPtr->getFlagField()->xyzSizeWithGhostLayer();
+ boundaryHandlingPtr->unpack(buffer, everyCell, true);
+ }
+
+ private:
+ const AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >* boundary_;
+}; // class BoundaryBlockDataHandling
+
+// helper function wrapper for adding the flag field to the block storage; since the input parameter for an
+// initialization function in field::addFlagFieldToStorage() is a std::function<void(FlagField_T*,IBlock* const)>, we
+// need a function wrapper that has both these input parameters; as FlagInfo< FlagField_T >::registerFlags() does not
+// have both of these input parameters, a function wrapper with both input parameters is created and the second input
+// parameter is simply ignored inside the function wrapper
+template< typename FlagField_T >
+void flagFieldInitFunction(FlagField_T* flagField, IBlock* const, const Set< field::FlagUID >& obstacleIDs,
+ const Set< field::FlagUID >& outflowIDs, const Set< field::FlagUID >& inflowIDs,
+ const Set< field::FlagUID >& freeSlipIDs)
+{
+ // register flags in the flag field
+ walberla::free_surface::FlagInfo< FlagField_T >::registerFlags(flagField, obstacleIDs, outflowIDs, inflowIDs,
+ freeSlipIDs);
+}
+
+} // namespace internal
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::AntidunesBoundaryHandling(
+ const std::shared_ptr< StructuredBlockForest >& blockForest, BlockDataID pdfFieldID, BlockDataID fillLevelID,
+ BlockDataID particleFieldID, const shared_ptr< ParticleAccessor_T >& ac,
+ std::function< real_t(const Vector3< real_t >&) > hydrostaticDensityFct)
+ : blockForest_(blockForest), pdfFieldID_(pdfFieldID), fillFieldID_(fillLevelID), particleFieldID_(particleFieldID),
+ ac_(ac), hydrostaticDensityFct_(std::move(hydrostaticDensityFct)), comm_(blockForest)
+{
+ // initialize obstacleIDs
+ Set< FlagUID > obstacleIDs;
+ obstacleIDs += noSlipFlagID;
+ obstacleIDs += ubbFlagID;
+ obstacleIDs += ubbInflowFlagID;
+ obstacleIDs += pressureFlagID;
+ obstacleIDs += pressureOutflowFlagID;
+ obstacleIDs += freeSlipFlagID;
+ obstacleIDs += outletFlagID;
+ obstacleIDs += movingObstacleFlagID;
+
+ // initialize outflowIDs
+ Set< FlagUID > outflowIDs;
+ outflowIDs += pressureOutflowFlagID;
+ outflowIDs += outletFlagID;
+
+ // initialize outflowIDs
+ Set< FlagUID > inflowIDs;
+ inflowIDs += ubbInflowFlagID;
+
+ // initialize freeSlipIDs
+ Set< FlagUID > freeSlipIDs;
+ freeSlipIDs += freeSlipFlagID;
+
+ // create callable function wrapper with input arguments 1 and 2 unset, whereas arguments 3, 4 and 5 are set to be
+ // obstacleIDs, outflowIDs, and inflowIDs, respectively; this is necessary for field::addFlagFieldToStorage()
+ auto ffInitFunc = std::bind(internal::flagFieldInitFunction< FlagField_T >, std::placeholders::_1,
+ std::placeholders::_2, obstacleIDs, outflowIDs, inflowIDs, freeSlipIDs);
+
+ // IMPORTANT REMARK: The flag field needs two ghost layers because of function advectMass(). There, the flags of all
+ // D3Q* neighbors are determined for each cell, including cells in the first ghost layer. Therefore, all D3Q*
+ // neighbors of the first ghost layer must be accessible. This requires a second ghost layer.
+ flagFieldID_ = field::addFlagFieldToStorage< FlagField_T >(blockForest, "Flags", uint_c(2), true, ffInitFunc);
+
+ // create FlagInfo
+ flagInfo_ = walberla::free_surface::FlagInfo< FlagField_T >(obstacleIDs, outflowIDs, inflowIDs, freeSlipIDs);
+ WALBERLA_ASSERT(flagInfo_.isConsistentAcrossBlocksAndProcesses(blockForest, flagFieldID_));
+
+ // add boundary handling to blockForest
+ handlingID_ = blockForest_->addBlockData(
+ std::make_shared<
+ internal::BoundaryBlockDataHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T > >(this),
+ "Boundary Handling");
+
+ // create communication object with fill level field, since fill levels determine the flags during the simulation
+ comm_.addPackInfo(std::make_shared< field::communication::PackInfo< ScalarField_T > >(fillFieldID_));
+}
+
+// define IDs (static const variables)
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::noSlipFlagID =
+ field::FlagUID("NoSlip");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::ubbFlagID =
+ field::FlagUID("UBB");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::ubbInflowFlagID =
+ field::FlagUID("UBB_Inflow");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::pressureFlagID =
+ field::FlagUID("Pressure");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::pressureOutflowFlagID =
+ field::FlagUID("PressureOutflow");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::outletFlagID =
+ field::FlagUID("Outlet");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::freeSlipFlagID =
+ field::FlagUID("FreeSlip");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const field::FlagUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::movingObstacleFlagID =
+ field::FlagUID("MovingObstacle");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::noSlipBoundaryID =
+ BoundaryUID("NoSlip");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::ubbBoundaryID =
+ BoundaryUID("UBB");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::ubbInflowBoundaryID =
+ BoundaryUID("UBB_Inflow");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::pressureBoundaryID =
+ BoundaryUID("Pressure");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T,
+ ParticleAccessor_T >::pressureOutflowBoundaryID =
+ BoundaryUID("PressureOutflow");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::outletBoundaryID =
+ BoundaryUID("Outlet");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::freeSlipBoundaryID =
+ BoundaryUID("FreeSlip");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+const BoundaryUID AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T,
+ ParticleAccessor_T >::movingObstacleBoundaryID =
+ BoundaryUID("MovingObstacle");
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+geometry::initializer::InitializationManager
+ AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::getInitManager()
+{
+ using namespace geometry::initializer;
+
+ InitializationManager initManager(blockForest_->getBlockStorage());
+
+ // define initializers
+ auto cellIntvInit = std::make_shared< BoundaryFromCellInterval< BoundaryHandling_T > >(*blockForest_, handlingID_);
+ auto borderInit = std::make_shared< BoundaryFromDomainBorder< BoundaryHandling_T > >(*blockForest_, handlingID_);
+ auto imgInit =
+ std::make_shared< BoundaryFromImage< BoundaryHandling_T, geometry::GrayScaleImage > >(*blockForest_, handlingID_);
+ auto bodyInit = std::make_shared< OverlapFieldFromBody >(*blockForest_, fillFieldID_);
+
+ // register initializers
+ initManager.registerInitializer("CellInterval", cellIntvInit);
+ initManager.registerInitializer("Border", borderInit);
+ initManager.registerInitializer("Image", imgInit);
+ initManager.registerInitializer("Body", bodyInit);
+
+ return initManager;
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::initFromConfig(
+ const Config::BlockHandle& configBlock)
+{
+ // initialize from config file
+ getInitManager().init(configBlock);
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+template< typename Body_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::addFreeSurfaceObject(
+ const Body_T& body, bool addOrSubtract)
+{
+ geometry::initializer::OverlapFieldFromBody(*blockForest_, fillFieldID_).init(body, addOrSubtract);
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setNoSlipAtBorder(
+ stencil::Direction d, cell_idx_t wallDistance)
+{
+ geometry::initializer::BoundaryFromDomainBorder< BoundaryHandling_T > init(*blockForest_, handlingID_);
+ init.init(noSlipFlagID, d, wallDistance);
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setNoSlipAtAllBorders(
+ cell_idx_t wallDistance)
+{
+ geometry::initializer::BoundaryFromDomainBorder< BoundaryHandling_T > init(*blockForest_, handlingID_);
+ init.initAllBorders(noSlipFlagID, wallDistance);
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setNoSlipInCell(
+ const Cell& globalCell)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+
+ // transform cell in global coordinates to cell in block local coordinates
+ Cell blockLocalCell;
+ blockForest_->transformGlobalToBlockLocalCell(blockLocalCell, *blockIt, globalCell);
+
+ handling->forceBoundary(noSlipFlagID, blockLocalCell[0], blockLocalCell[1], blockLocalCell[2]);
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setFreeSlipAtBorder(
+ stencil::Direction d, cell_idx_t wallDistance)
+{
+ geometry::initializer::BoundaryFromDomainBorder< BoundaryHandling_T > init(*blockForest_, handlingID_);
+ init.init(freeSlipFlagID, d, wallDistance);
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T,
+ ParticleAccessor_T >::setFreeSlipAtAllBorders(cell_idx_t wallDistance)
+{
+ geometry::initializer::BoundaryFromDomainBorder< BoundaryHandling_T > init(*blockForest_, handlingID_);
+ init.initAllBorders(freeSlipFlagID, wallDistance);
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setFreeSlipInCell(
+ const Cell& globalCell)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+
+ // transform cell in global coordinates to cell in block local coordinates
+ Cell blockLocalCell;
+ blockForest_->transformGlobalToBlockLocalCell(blockLocalCell, *blockIt, globalCell);
+
+ handling->forceBoundary(freeSlipFlagID, blockLocalCell[0], blockLocalCell[1], blockLocalCell[2]);
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setPressure(
+ real_t density)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+ Pressure_T& pressure =
+ handling->template getBoundaryCondition< Pressure_T >(handling->getBoundaryUID(pressureFlagID));
+ pressure.setLatticeDensity(density);
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setUBBInCell(
+ const Cell& globalCell, const Vector3< real_t >& velocity)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+
+ typename UBB_Inflow_T::Velocity ubbVel(velocity);
+
+ // transform cell in global coordinates to cell in block-local coordinates
+ Cell blockLocalCell;
+ blockForest_->transformGlobalToBlockLocalCell(blockLocalCell, *blockIt, globalCell);
+
+ // get block cell bounding box to check if cell is contained in block
+ CellInterval blockCellBB = blockForest_->getBlockCellBB(*blockIt);
+
+ // flag field has two ghost layers so blockCellBB is actually larger than returned; this is relevant for setups
+ // where the UBB is set in a ghost layer cell
+ blockCellBB.expand(cell_idx_c(2));
+
+ if (blockCellBB.contains(globalCell))
+ {
+ handling->forceBoundary(ubbFlagID, blockLocalCell[0], blockLocalCell[1], blockLocalCell[2], ubbVel);
+ }
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setInflowInCell(
+ const Cell& globalCell, const Vector3< real_t >& velocity)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+
+ typename UBB_Inflow_T::Velocity ubbVel(velocity);
+
+ // transform cell in global coordinates to cell in block-local coordinates
+ Cell blockLocalCell;
+ blockForest_->transformGlobalToBlockLocalCell(blockLocalCell, *blockIt, globalCell);
+
+ // get block cell bounding box to check if cell is contained in block
+ CellInterval blockCellBB = blockForest_->getBlockCellBB(*blockIt);
+
+ // flag field has two ghost layers so blockCellBB is actually larger than returned; this is relevant for setups
+ // where the UBB is set in a ghost layer cell
+ blockCellBB.expand(cell_idx_c(2));
+
+ if (blockCellBB.contains(globalCell))
+ {
+ handling->forceBoundary(ubbInflowFlagID, blockLocalCell[0], blockLocalCell[1], blockLocalCell[2], ubbVel);
+ }
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::setPressureOutflow(
+ real_t density)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+ Pressure_T& pressure =
+ handling->template getBoundaryCondition< Pressure_T >(handling->getBoundaryUID(pressureOutflowFlagID));
+ pressure.setLatticeDensity(density);
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T, ParticleAccessor_T >::enableBubbleOutflow(
+ BubbleModelBase* bubbleModel)
+{
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ BoundaryHandling_T* const handling = blockIt->template getData< BoundaryHandling_T >(handlingID_);
+
+ // get pressure from boundary handling
+ Pressure_T& pressure =
+ handling->template getBoundaryCondition< Pressure_T >(handling->getBoundaryUID(pressureFlagID));
+ Pressure_T& pressureOutflow =
+ handling->template getBoundaryCondition< Pressure_T >(handling->getBoundaryUID(pressureOutflowFlagID));
+
+ // set pressure in bubble model
+ pressure.setBubbleModel(bubbleModel);
+ pressureOutflow.setBubbleModel(bubbleModel);
+ }
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+Vector3< bool > AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T,
+ ParticleAccessor_T >::isObstacleInGlobalGhostLayer()
+{
+ Vector3< bool > isObstacleInGlobalGhostLayer(false, false, false);
+
+ for (auto blockIt = blockForest_->begin(); blockIt != blockForest_->end(); ++blockIt)
+ {
+ const FlagField_T* const flagField = blockIt->template getData< const FlagField_T >(flagFieldID_);
+
+ const CellInterval domainCellBB = blockForest_->getDomainCellBB();
+
+ // disable OpenMP such that loop termination works correctly
+ WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ_OMP(flagField, uint_c(1), omp critical, {
+ // get cell in global coordinates
+ Cell globalCell = Cell(x, y, z);
+ blockForest_->transformBlockLocalToGlobalCell(globalCell, *blockIt);
+
+ // check if the current cell is located in a global ghost layer
+ const bool isCellInGlobalGhostLayerX =
+ globalCell[0] < domainCellBB.xMin() || globalCell[0] > domainCellBB.xMax();
+
+ const bool isCellInGlobalGhostLayerY =
+ globalCell[1] < domainCellBB.yMin() || globalCell[1] > domainCellBB.yMax();
+
+ const bool isCellInGlobalGhostLayerZ =
+ globalCell[2] < domainCellBB.zMin() || globalCell[2] > domainCellBB.zMax();
+
+ // skip corners, as they do not influence periodic communication
+ if ((isCellInGlobalGhostLayerX && (isCellInGlobalGhostLayerY || isCellInGlobalGhostLayerZ)) ||
+ (isCellInGlobalGhostLayerY && isCellInGlobalGhostLayerZ))
+ {
+ continue;
+ }
+
+ if (!isObstacleInGlobalGhostLayer[0] && isCellInGlobalGhostLayerX &&
+ isPartOfMaskSet(flagField->get(x, y, z), flagField->getMask(flagInfo_.getObstacleIDSet())))
+ {
+ isObstacleInGlobalGhostLayer[0] = true;
+ }
+
+ if (!isObstacleInGlobalGhostLayer[1] && isCellInGlobalGhostLayerY &&
+ isPartOfMaskSet(flagField->get(x, y, z), flagField->getMask(flagInfo_.getObstacleIDSet())))
+ {
+ isObstacleInGlobalGhostLayer[1] = true;
+ }
+
+ if (!isObstacleInGlobalGhostLayer[2] && isCellInGlobalGhostLayerZ &&
+ isPartOfMaskSet(flagField->get(x, y, z), flagField->getMask(flagInfo_.getObstacleIDSet())))
+ {
+ isObstacleInGlobalGhostLayer[2] = true;
+ }
+
+ if (isObstacleInGlobalGhostLayer[0] && isObstacleInGlobalGhostLayer[1] && isObstacleInGlobalGhostLayer[2])
+ {
+ break; // there is no need to check other cells on this block
+ }
+ }) // WALBERLA_FOR_ALL_CELLS_INCLUDING_GHOST_LAYER_XYZ_OMP
+ }
+
+ mpi::allReduceInplace(isObstacleInGlobalGhostLayer, mpi::LOGICAL_OR);
+
+ return isObstacleInGlobalGhostLayer;
+}
+
+template< typename LatticeModel_T, typename FlagField_T, typename ScalarField_T, typename ParticleAccessor_T >
+void AntidunesBoundaryHandling< LatticeModel_T, FlagField_T, ScalarField_T,
+ ParticleAccessor_T >::initFlagsFromFillLevel()
+{
+ const Vector3< bool > isObstacleInGlobalGhostLayer = this->isObstacleInGlobalGhostLayer();
+
+ WALBERLA_ROOT_SECTION()
+ {
+ if ((blockForest_->isXPeriodic() && isObstacleInGlobalGhostLayer[0]) ||
+ (blockForest_->isYPeriodic() && isObstacleInGlobalGhostLayer[1]) ||
+ (blockForest_->isZPeriodic() && isObstacleInGlobalGhostLayer[2]))
+ {
+ WALBERLA_LOG_WARNING_ON_ROOT(
+ "WARNING: An obstacle cell is located in a global outermost ghost layer in a periodic "
+ "direction. Be aware that due to periodicity, this obstacle cell will be "
+ "overwritten during communication.");
+ }
+ }
+
+ // communicate fill level (neighborhood is used in initialization)
+ comm_();
+
+ // initialize fill level in boundaries (with value 1), i.e., obstacles such that the bubble model does not detect
+ // obstacles as gas cells
+ walberla::free_surface::initFillLevelsInBoundaries< BoundaryHandling_T, typename LatticeModel_T::Stencil,
+ ScalarField_T >(blockForest_, handlingID_, fillFieldID_);
+
+ // clear and initialize flags in every cell according to the fill level
+ walberla::free_surface::initFlagsFromFillLevels< BoundaryHandling_T, typename LatticeModel_T::Stencil, FlagField_T,
+ const ScalarField_T >(blockForest_, flagInfo_, handlingID_,
+ fillFieldID_);
+}
+
+} // namespace free_surface
+} // namespace antidunes
+} // namespace walberla
diff --git a/apps/showcases/Antidunes/AntidunesLatticeModelGeneration.py b/apps/showcases/Antidunes/AntidunesLatticeModelGeneration.py
new file mode 100644
index 0000000000000000000000000000000000000000..92376f290c3d1d54e4dc76bd0598689151ff16c4
--- /dev/null
+++ b/apps/showcases/Antidunes/AntidunesLatticeModelGeneration.py
@@ -0,0 +1,70 @@
+import sympy as sp
+
+from lbmpy.creationfunctions import LBMConfig, LBMOptimisation, create_lb_collision_rule
+from lbmpy.enums import ForceModel, Method, Stencil
+from lbmpy.stencils import LBStencil
+
+from pystencils_walberla import CodeGeneration
+from lbmpy_walberla import generate_lattice_model
+
+# general parameters
+generatedMethod = "Cumulants" # "TRT", "SRT"
+stencilStr = "D3Q27"
+stencil = LBStencil(Stencil.D3Q19 if stencilStr == "D3Q19" else Stencil.D3Q27)
+force = sp.symbols("force_:3")
+layout = "fzyx"
+
+if generatedMethod == "Cumulants":
+ omega = sp.Symbol("omega")
+ # method definition
+ lbm_config = LBMConfig(
+ stencil=stencil,
+ method=Method.CUMULANT,
+ relaxation_rate=omega,
+ compressible=True,
+ force=force,
+ zero_centered=False,
+ streaming_pattern="pull",
+ galilean_correction=True if stencil == LBStencil(Stencil.D3Q27) else False,
+ ) # free surface implementation only works with pull pattern
+elif generatedMethod == "TRT":
+ omega_e = sp.Symbol("omega_e")
+ omega_o = sp.Symbol("omega_o")
+ # method definition
+ lbm_config = LBMConfig(
+ stencil=stencil,
+ method=Method.TRT,
+ smagorinsky=False,
+ relaxation_rates=[omega_e, omega_o],
+ compressible=True,
+ force=force,
+ force_model=ForceModel.GUO,
+ zero_centered=False,
+ streaming_pattern="pull",
+ ) # free surface implementation only works with pull pattern
+elif generatedMethod == "SRT":
+ omega = sp.Symbol("omega")
+ # method definition
+ lbm_config = LBMConfig(
+ stencil=stencil,
+ method=Method.SRT,
+ smagorinsky=True,
+ relaxation_rate=omega,
+ compressible=True,
+ force=force,
+ force_model=ForceModel.GUO,
+ zero_centered=False,
+ streaming_pattern="pull",
+ ) # free surface implementation only works with pull pattern
+
+# optimizations to be used by the code generator
+lbm_opt = LBMOptimisation(cse_global=True, field_layout=layout)
+
+collision_rule = create_lb_collision_rule(
+ lbm_config=lbm_config, lbm_optimisation=lbm_opt
+)
+
+with CodeGeneration() as ctx:
+ generate_lattice_model(
+ ctx, "AntidunesLatticeModel", collision_rule, field_layout=layout
+ )
diff --git a/apps/showcases/Antidunes/BedGeneration.cpp b/apps/showcases/Antidunes/BedGeneration.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..baac22f4d82bcd9ed161df5136c9919a13d83ac1
--- /dev/null
+++ b/apps/showcases/Antidunes/BedGeneration.cpp
@@ -0,0 +1,265 @@
+//======================================================================================================================
+//
+// 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 BedGeneration.cpp
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/Environment.h"
+#include "core/grid_generator/HCPIterator.h"
+#include "core/grid_generator/SCIterator.h"
+#include "core/math/Random.h"
+
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/LinkedCells.h"
+#include "mesa_pd/data/ParticleAccessorWithBaseShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/AssocToBlock.h"
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/ExplicitEuler.h"
+#include "mesa_pd/kernel/InsertParticleIntoLinkedCells.h"
+#include "mesa_pd/kernel/LinearSpringDashpot.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/SyncNextNeighborsBlockForest.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "vtk/VTKOutput.h"
+
+#include <mesa_pd/mpi/notifications/ForceTorqueNotification.h>
+
+#include "Utility.h"
+
+namespace walberla
+{
+namespace antidunes
+{
+
+using namespace mesa_pd;
+
+int main(int argc, char** argv)
+{
+ Environment env(argc, argv);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ // Config
+ auto cfg = env.config();
+ if (cfg == nullptr) WALBERLA_ABORT("No config specified!");
+ WALBERLA_LOG_INFO_ON_ROOT(*cfg);
+ const Config::BlockHandle bedGenerationConf = cfg->getBlock("BedGeneration");
+
+ Vec3 domainSize_SI = bedGenerationConf.getParameter< Vec3 >("domainSize_SI");
+ Vector3< int > blocks = bedGenerationConf.getParameter< Vector3< int > >("blocks");
+ real_t diameter_SI = bedGenerationConf.getParameter< real_t >("diameter_SI");
+ real_t gravity_SI = bedGenerationConf.getParameter< real_t >("gravity_SI");
+ real_t densityFluid_SI = bedGenerationConf.getParameter< real_t >("densityFluid_SI");
+ real_t densityParticle_SI = bedGenerationConf.getParameter< real_t >("densityParticle_SI");
+ real_t generationSpacing_SI = bedGenerationConf.getParameter< real_t >("generationSpacing_SI");
+ real_t initialVelocity_SI = bedGenerationConf.getParameter< real_t >("initialVelocity_SI");
+ real_t dt_SI = bedGenerationConf.getParameter< real_t >("dt_SI");
+ real_t frictionCoefficient = bedGenerationConf.getParameter< real_t >("frictionCoefficient");
+ real_t restitutionCoefficient = bedGenerationConf.getParameter< real_t >("restitutionCoefficient");
+ real_t collisionTime_SI = bedGenerationConf.getParameter< real_t >("collisionTime_SI");
+ real_t poissonsRatio = bedGenerationConf.getParameter< real_t >("poissonsRatio");
+ uint_t timeSteps = bedGenerationConf.getParameter< uint_t >("timeSteps");
+ uint_t visSpacing = bedGenerationConf.getParameter< uint_t >("visSpacing");
+ std::string outFileName = bedGenerationConf.getParameter< std::string >("outFileName");
+ bool denseBottomLayer = bedGenerationConf.getParameter< bool >("denseBottomLayer");
+ real_t bottomLayerOffsetFactor = bedGenerationConf.getParameter< real_t >("bottomLayerOffsetFactor");
+
+ // BlockForest
+ math::AABB simulationDomain_SI(real_t(0.0), real_t(0.0), real_t(0.0), domainSize_SI[0], domainSize_SI[1],
+ domainSize_SI[2]);
+ Vector3< bool > isPeriodic{ true, true, false };
+
+ shared_ptr< BlockForest > forest = blockforest::createBlockForest(simulationDomain_SI, blocks, isPeriodic);
+ auto domain = std::make_shared< mesa_pd::domain::BlockForestDomain >(forest);
+
+ // MesaPD data structures
+ auto ps = std::make_shared< data::ParticleStorage >(1);
+ data::ParticleAccessorWithBaseShape accessor(ps);
+
+ // Init spheres
+ real_t minDiameter_SI = diameter_SI * real_t(0.9);
+ real_t maxDiameter_SI = diameter_SI * real_t(1.1);
+
+ math::AABB generationDomain_SI(simulationDomain_SI.xMin(), simulationDomain_SI.yMin(),
+ simulationDomain_SI.zMin() + diameter_SI, simulationDomain_SI.xMax(),
+ simulationDomain_SI.yMax(), simulationDomain_SI.zMax());
+
+ for (auto pt :
+ grid_generator::SCGrid(generationDomain_SI, Vec3(generationSpacing_SI) * real_c(0.5), generationSpacing_SI))
+ {
+ auto diameter = math::realRandom< real_t >(minDiameter_SI, maxDiameter_SI);
+
+ if (!domain->isContainedInLocalSubdomain(pt, real_t(0))) continue;
+ auto p = ps->create();
+ p->getPositionRef() = pt;
+ p->getInteractionRadiusRef() = diameter * real_t(0.5);
+ p->getBaseShapeRef() = std::make_shared< data::Sphere >(p->getInteractionRadius());
+ p->getBaseShapeRef()->updateMassAndInertia(densityParticle_SI);
+
+ p->setLinearVelocity(Vec3(real_t(0.1) * math::realRandom(-initialVelocity_SI, initialVelocity_SI),
+ real_t(0.1) * math::realRandom(-initialVelocity_SI, initialVelocity_SI),
+ -initialVelocity_SI));
+ p->getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
+ p->getTypeRef() = 0;
+ }
+
+ math::AABB bottomLayerDomain_SI(simulationDomain_SI.xMin(), simulationDomain_SI.yMin(), simulationDomain_SI.zMin(),
+ simulationDomain_SI.xMax(), simulationDomain_SI.yMax(), diameter_SI);
+
+ real_t bottomLayerSpacing = bottomLayerDomain_SI.xSize() / std::floor(bottomLayerDomain_SI.xSize() / diameter_SI);
+ real_t bottomLayerYStretchFactor =
+ real_t((bottomLayerDomain_SI.ySize() / (sqrt(3_r) * bottomLayerSpacing)) /
+ std::floor(bottomLayerDomain_SI.ySize() / (sqrt(3_r) * bottomLayerSpacing)));
+ if (denseBottomLayer)
+ {
+ bottomLayerYStretchFactor = real_t((bottomLayerDomain_SI.ySize() / (sqrt(3_r) * bottomLayerSpacing)) /
+ std::ceil(bottomLayerDomain_SI.ySize() / (sqrt(3_r) * bottomLayerSpacing)));
+ }
+ WALBERLA_LOG_INFO_ON_ROOT(bottomLayerSpacing << " " << bottomLayerYStretchFactor);
+ for (auto pt : grid_generator::HCPGrid(bottomLayerDomain_SI, Vec3(diameter_SI) * real_c(0.5), bottomLayerSpacing))
+ {
+ auto diameter = math::realRandom< real_t >(minDiameter_SI, maxDiameter_SI);
+ auto zCoord = math::realRandom< real_t >(real_t(1e-10), diameter_SI * bottomLayerOffsetFactor);
+ Vec3 position{ pt[0], pt[1] * bottomLayerYStretchFactor, zCoord };
+
+ if (!domain->isContainedInLocalSubdomain(position, real_t(0))) continue;
+ auto p = ps->create();
+ p->getPositionRef() = position;
+ p->getInteractionRadiusRef() = diameter * real_t(0.5);
+ p->getBaseShapeRef() = std::make_shared< data::Sphere >(p->getInteractionRadius());
+ p->getBaseShapeRef()->updateMassAndInertia(densityParticle_SI);
+
+ p->getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
+ p->getTypeRef() = 0;
+ mesa_pd::data::particle_flags::set(p->getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+ }
+
+ createPlane(*ps, simulationDomain_SI.minCorner(), Vec3(real_t(0), real_t(0), real_t(1)));
+ createPlane(*ps, simulationDomain_SI.maxCorner(), Vec3(real_t(0), real_t(0), real_t(-1)));
+
+ // VTK
+ auto vtkDomainOutput =
+ walberla::vtk::createVTKOutput_DomainDecomposition(forest, "domain_decomposition", 1, "vtk", "simulation_step");
+ vtkDomainOutput->write();
+
+ auto particleVtkOutput = make_shared< mesa_pd::vtk::ParticleVtkOutput >(ps);
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleLinearVelocity >("velocity");
+ particleVtkOutput->addOutput< mesa_pd::data::SelectParticleInteractionRadius >("radius");
+ particleVtkOutput->setParticleSelector([](const data::ParticleStorage::iterator& pIt) {
+ using namespace walberla::mesa_pd::data::particle_flags;
+ return (pIt->getBaseShape()->getShapeType() == data::Sphere::SHAPE_TYPE) && !isSet(pIt->getFlags(), GHOST);
+ });
+ auto vtkWriter = walberla::vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", 1, "vtk",
+ "simulation_step", false, false);
+
+ // Init kernels
+ kernel::ExplicitEuler explicitEulerWithShape(dt_SI);
+ kernel::LinearSpringDashpot dem(1);
+ dem.setFrictionCoefficientDynamic(0, 0, frictionCoefficient);
+ real_t kappa = real_t(2) * (real_t(1) - poissonsRatio) / (real_t(2) - poissonsRatio); // from Thornton et al
+
+ kernel::AssocToBlock assoc(forest);
+ mesa_pd::mpi::ReduceProperty RP;
+ mesa_pd::mpi::SyncNextNeighborsBlockForest SNN;
+
+ ps->forEachParticle(false, kernel::SelectLocal(), accessor, assoc, accessor);
+
+ // initial sync
+ SNN(*ps, forest, domain);
+
+ real_t averageVelocity = real_t(0);
+ uint_t currentNumParticles = 0;
+ real_t maxVelocity = real_t(0);
+ real_t maxHeight = real_t(0);
+
+ real_t linkedCellWidth = 1.01_r * maxDiameter_SI;
+ data::LinkedCells linkedCells(domain->getUnionOfLocalAABBs().getExtended(linkedCellWidth), linkedCellWidth);
+ kernel::InsertParticleIntoLinkedCells ipilc;
+
+ for (uint_t i = 0; i < timeSteps; ++i)
+ {
+ if (i % visSpacing == 0) { vtkWriter->write(); }
+
+ ps->forEachParticle(false, kernel::SelectLocal(), accessor, assoc, accessor);
+
+ SNN(*ps, forest, domain);
+
+ // gravity - buoyancy
+ ps->forEachParticle(
+ false, kernel::SelectLocal(), accessor,
+ [densityParticle_SI, densityFluid_SI, gravity_SI](const size_t idx, auto& ac) {
+ mesa_pd::addForceAtomic(
+ idx, ac, Vec3(0, 0, -(densityParticle_SI - densityFluid_SI) * ac.getVolume(idx) * gravity_SI));
+ },
+ accessor);
+
+ linkedCells.clear();
+ ps->forEachParticle(false, kernel::SelectAll(), accessor, ipilc, accessor, linkedCells);
+ linkedCells.forEachParticlePairHalf(
+ false, kernel::ExcludeInfiniteInfinite(), accessor,
+ [restitutionCoefficient, collisionTime_SI, kappa, domain, &dem, dt_SI](const size_t idx1, const size_t idx2,
+ auto& ac) {
+ kernel::DoubleCast double_cast;
+ mesa_pd::mpi::ContactFilter contact_filter;
+ collision_detection::AnalyticContactDetection acd;
+
+ if (double_cast(idx1, idx2, ac, acd, ac))
+ {
+ if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), *domain))
+ {
+ auto meff = real_t(1) / (ac.getInvMass(idx1) + ac.getInvMass(idx2));
+ dem.setStiffnessAndDamping(0, 0, restitutionCoefficient, collisionTime_SI, kappa, meff);
+ dem(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(),
+ acd.getPenetrationDepth(), dt_SI);
+ }
+ }
+ },
+ accessor);
+
+ RP.operator()< ForceTorqueNotification >(*ps);
+
+ ps->forEachParticle(false, kernel::SelectLocal(), accessor, explicitEulerWithShape, accessor);
+
+ getAverageVelocity(accessor, averageVelocity, maxVelocity, currentNumParticles, maxHeight);
+
+ SNN(*ps, forest, domain);
+
+ if (i % 1000 == 0)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT("Timestep " << i << " / " << timeSteps << ", average velocity = " << averageVelocity
+ << ", max velocity = " << maxVelocity << ", #particles = "
+ << currentNumParticles << ", max height = " << maxHeight);
+ }
+ }
+
+ writeSphereInformationToFile(outFileName, *ps, domainSize_SI);
+
+ return EXIT_SUCCESS;
+}
+} // namespace antidunes
+} // namespace walberla
+
+int main(int argc, char** argv) { return walberla::antidunes::main(argc, argv); }
diff --git a/apps/showcases/Antidunes/BedGeneration.prm b/apps/showcases/Antidunes/BedGeneration.prm
new file mode 100644
index 0000000000000000000000000000000000000000..af2ecf9645a560468350281c27319a140e2d5bc4
--- /dev/null
+++ b/apps/showcases/Antidunes/BedGeneration.prm
@@ -0,0 +1,20 @@
+BedGeneration{
+ domainSize_SI < 0.8, 0.015, 0.2 >;
+ blocks < 3, 3, 1 >;
+ diameter_SI 0.0029;
+ gravity_SI 9.81;
+ densityFluid_SI 1000;
+ densityParticle_SI 2550;
+ generationSpacing_SI 0.005;
+ initialVelocity_SI 1;
+ dt_SI 5e-5;
+ frictionCoefficient 0.5;
+ restitutionCoefficient 0.1;
+ collisionTime_SI 5e-4;
+ poissonsRatio 0.22;
+ timeSteps 10000;
+ visSpacing 100;
+ outFileName spheres_out.dat;
+ denseBottomLayer False;
+ bottomLayerOffsetFactor 1.0;
+}
diff --git a/apps/showcases/Antidunes/CMakeLists.txt b/apps/showcases/Antidunes/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a0feea41cc2502dc568f9ee3aa2271e2932760e9
--- /dev/null
+++ b/apps/showcases/Antidunes/CMakeLists.txt
@@ -0,0 +1,16 @@
+waLBerla_link_files_to_builddir( *.prm )
+
+if( WALBERLA_BUILD_WITH_CODEGEN )
+ walberla_generate_target_from_python( NAME AntidunesLatticeModelGeneration
+ FILE AntidunesLatticeModelGeneration.py
+ OUT_FILES AntidunesLatticeModel.cpp AntidunesLatticeModel.h )
+
+ waLBerla_add_executable(NAME Antidunes
+ FILES Antidunes.cpp PIDController.cpp
+ DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling
+ postprocessing timeloop vtk AntidunesLatticeModelGeneration)
+endif()
+
+waLBerla_add_executable(NAME BedGeneration
+ FILES BedGeneration.cpp
+ DEPENDS blockforest core domain_decomposition mesa_pd vtk)
diff --git a/apps/showcases/Antidunes/PIDController.cpp b/apps/showcases/Antidunes/PIDController.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..83df9f470f30f259a4cdab8073ff4e67ff0f1f86
--- /dev/null
+++ b/apps/showcases/Antidunes/PIDController.cpp
@@ -0,0 +1,115 @@
+//======================================================================================================================
+/*!
+ * \file PIDController.cpp
+ */
+//======================================================================================================================
+
+#include "PIDController.h"
+
+#include <algorithm>
+#include <fstream>
+
+using namespace walberla;
+using walberla::real_t;
+
+PIDController::PIDController()
+ : commandVariable_(0), actuatingVariable_(0), proportionalGain_(0), derivateGain_(0), integralGain_(0), maxRamp_(0),
+ minActuatingVariable_(0), maxActuatingVariable_(0), errorIntegral_(0)
+{
+ std::fill(errorHistory_, errorHistory_ + sizeof(errorHistory_) / sizeof(real_t), real_t(0));
+}
+
+PIDController::PIDController(const real_t commandVariable, const real_t initialActuatingVariable,
+ const real_t proportionalGain, const real_t derivateGain, const real_t integralGain,
+ const real_t maxRamp, const real_t minActuatingVariable, const real_t maxActuatingVariable)
+ : commandVariable_(commandVariable), actuatingVariable_(initialActuatingVariable),
+ proportionalGain_(proportionalGain), derivateGain_(derivateGain), integralGain_(integralGain), maxRamp_(maxRamp),
+ minActuatingVariable_(minActuatingVariable), maxActuatingVariable_(maxActuatingVariable), errorIntegral_(0)
+{
+ std::fill(errorHistory_, errorHistory_ + sizeof(errorHistory_) / sizeof(real_t), real_t(0));
+
+ if (integralGain_ > real_t(0))
+ errorIntegral_ = initialActuatingVariable / integralGain_;
+ else
+ errorIntegral_ = real_t(0);
+}
+
+PIDController::PIDController(const real_t commandVariable, const real_t initialActuatingVariable,
+ const real_t proportionalGain, const real_t derivateGain, const real_t integralGain)
+ : commandVariable_(commandVariable), actuatingVariable_(initialActuatingVariable),
+ proportionalGain_(proportionalGain), derivateGain_(derivateGain), integralGain_(integralGain),
+ maxRamp_(std::numeric_limits< real_t >::max()), minActuatingVariable_(-std::numeric_limits< real_t >::max()),
+ maxActuatingVariable_(std::numeric_limits< real_t >::max()), errorIntegral_(0)
+{
+ std::fill(errorHistory_, errorHistory_ + sizeof(errorHistory_) / sizeof(real_t), real_t(0));
+
+ if (integralGain_ > real_t(0))
+ errorIntegral_ = initialActuatingVariable / integralGain_;
+ else
+ errorIntegral_ = real_t(0);
+}
+
+real_t PIDController::update(const real_t controlledVariable)
+{
+ static const real_t ONE_OVER_SIX = real_t(1) / real_t(6);
+ const real_t error = commandVariable_ - controlledVariable;
+
+ const real_t d =
+ (error + real_t(3) * errorHistory_[0] - real_t(3) * errorHistory_[1] - errorHistory_[2]) * ONE_OVER_SIX;
+ std::rotate(errorHistory_, errorHistory_ + 1, errorHistory_ + sizeof(errorHistory_) / sizeof(real_t));
+ errorHistory_[sizeof(errorHistory_) / sizeof(real_t) - size_t(1)] = error;
+
+ real_t newActuationVariable = proportionalGain_ * error + derivateGain_ * d + integralGain_ * errorIntegral_;
+
+ if (std::fabs(actuatingVariable_ - newActuationVariable) < maxRamp_)
+ {
+ errorIntegral_ += error;
+ newActuationVariable = proportionalGain_ * error + derivateGain_ * d + integralGain_ * errorIntegral_;
+ }
+
+ const real_t maxValue = std::min(actuatingVariable_ + maxRamp_, maxActuatingVariable_);
+ const real_t minValue = std::max(actuatingVariable_ - maxRamp_, minActuatingVariable_);
+
+ actuatingVariable_ = std::min(std::max(minValue, newActuationVariable), maxValue);
+
+ return actuatingVariable_;
+}
+
+void PIDController::writeStateToFile(std::string filename) const
+{
+ std::ofstream file;
+ file.open(filename);
+ file << std::setprecision(16);
+ file << commandVariable_ << std::endl;
+ file << actuatingVariable_ << std::endl;
+ file << proportionalGain_ << std::endl;
+ file << derivateGain_ << std::endl;
+ file << integralGain_ << std::endl;
+ file << maxRamp_ << std::endl;
+ file << minActuatingVariable_ << std::endl;
+ file << maxActuatingVariable_ << std::endl;
+ file << errorHistory_[0] << std::endl;
+ file << errorHistory_[1] << std::endl;
+ file << errorHistory_[2] << std::endl;
+ file << errorIntegral_ << std::endl;
+ file.close();
+}
+
+void PIDController::readStateFromFile(std::string filename)
+{
+ std::ifstream file;
+ file.open(filename);
+ file >> commandVariable_;
+ file >> actuatingVariable_;
+ file >> proportionalGain_;
+ file >> derivateGain_;
+ file >> integralGain_;
+ file >> maxRamp_;
+ file >> minActuatingVariable_;
+ file >> maxActuatingVariable_;
+ file >> errorHistory_[0];
+ file >> errorHistory_[1];
+ file >> errorHistory_[2];
+ file >> errorIntegral_;
+ file.close();
+}
diff --git a/apps/showcases/Antidunes/PIDController.h b/apps/showcases/Antidunes/PIDController.h
new file mode 100644
index 0000000000000000000000000000000000000000..aeae57260688bbeca012a240fac40209e94e2e55
--- /dev/null
+++ b/apps/showcases/Antidunes/PIDController.h
@@ -0,0 +1,51 @@
+//======================================================================================================================
+/*!
+ * \file PIDController.h
+ */
+//======================================================================================================================
+
+#pragma once
+
+#include <core/DataTypes.h>
+#include <core/logging/Logging.h>
+
+using namespace walberla;
+using walberla::real_t;
+
+class PIDController
+{
+ public:
+ PIDController();
+
+ PIDController(const real_t commandVariable, const real_t initialActuatingVariable, const real_t proportionalGain,
+ const real_t derivateGain, const real_t integralGain, const real_t maxRamp,
+ const real_t minActuatingVariable, const real_t maxActuatingVariable);
+
+ PIDController(const real_t commandVariable, const real_t initialActuatingVariable, const real_t proportionalGain,
+ const real_t derivateGain, const real_t integralGain);
+
+ real_t update(const real_t controlledVariable);
+
+ real_t getProportionalGain() const { return proportionalGain_; }
+ real_t getDerivateGain() const { return derivateGain_; }
+ real_t getIntegralGain() const { return integralGain_; }
+
+ void writeStateToFile(std::string filename) const;
+
+ void readStateFromFile(std::string filename);
+
+ private:
+ real_t commandVariable_;
+ real_t actuatingVariable_;
+
+ real_t proportionalGain_;
+ real_t derivateGain_;
+ real_t integralGain_;
+
+ real_t maxRamp_;
+ real_t minActuatingVariable_;
+ real_t maxActuatingVariable_;
+
+ real_t errorHistory_[3];
+ real_t errorIntegral_;
+};
diff --git a/apps/showcases/Antidunes/Utility.h b/apps/showcases/Antidunes/Utility.h
new file mode 100644
index 0000000000000000000000000000000000000000..e3e0e2149a01f26ee097287494d2d982ef9025c3
--- /dev/null
+++ b/apps/showcases/Antidunes/Utility.h
@@ -0,0 +1,472 @@
+//======================================================================================================================
+//
+// 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 Utility.h
+//! \author Samuel Kemmler <samuel.kemmler@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/shape/Sphere.h"
+
+#include <algorithm>
+#include <core/mpi/Broadcast.h>
+#include <core/mpi/MPITextFile.h>
+#include <core/mpi/Reduce.h>
+#include <functional>
+#include <iterator>
+
+namespace walberla
+{
+namespace antidunes
+{
+
+struct SphereSelector
+{
+ template< typename ParticleAccessor_T >
+ bool inline operator()(const size_t particleIdx, const ParticleAccessor_T& ac) const
+ {
+ static_assert(std::is_base_of< mesa_pd::data::IAccessor, ParticleAccessor_T >::value,
+ "Provide a valid accessor as template");
+ return ac.getBaseShape(particleIdx)->getShapeType() == mesa_pd::data::Sphere::SHAPE_TYPE;
+ }
+};
+
+void renameFile(const std::string& oldName, const std::string& newName)
+{
+ int result = std::rename(oldName.c_str(), newName.c_str());
+ if (result != 0)
+ WALBERLA_LOG_WARNING_ON_ROOT("Could not rename file " << oldName << " to " << newName << " with error code "
+ << result);
+}
+
+void write2DVectorToFile(const std::vector< real_t >& vec, uint_t len1, uint_t len2, std::string filename)
+{
+ std::ofstream file;
+ file.open(filename.c_str());
+ file.precision(5);
+
+ file << "# " << len1 << " " << len2 << "\n";
+
+ for (uint_t j = uint_t(0); j < len2; ++j)
+ {
+ for (uint_t i = uint_t(0); i < len1; ++i)
+ {
+ file << vec[i + j * len1] << "\n";
+ }
+ }
+ file.close();
+}
+
+template< typename ParticleAccessor_T >
+class BedloadTransportEvaluator
+{
+ public:
+ BedloadTransportEvaluator(const shared_ptr< ParticleAccessor_T >& ac, real_t normalizationFactor,
+ uint_t numParticles)
+ : ac_(ac), normalizationFactor_(normalizationFactor), numParticles_(numParticles)
+ {}
+
+ void operator()()
+ {
+ real_t transportRate(real_t(0));
+ real_t velocity(real_t(0));
+
+ for (uint_t i = uint_t(0); i < ac_->size(); ++i)
+ {
+ if (!isSet(ac_->getFlags(i), mesa_pd::data::particle_flags::GHOST) &&
+ !isSet(ac_->getFlags(i), mesa_pd::data::particle_flags::GLOBAL))
+ {
+ auto velX = ac_->getLinearVelocity(i)[0];
+ transportRate += velX * ac_->getVolume(i);
+ velocity += velX;
+ }
+ }
+
+ // only reduce to root
+ WALBERLA_MPI_SECTION()
+ {
+ mpi::reduceInplace(transportRate, mpi::SUM);
+ mpi::reduceInplace(velocity, mpi::SUM);
+ }
+
+ avgTransportRate_ = transportRate * normalizationFactor_;
+ averageVelocity_ = velocity / real_c(numParticles_);
+ }
+
+ // sum_i V_p,i * u_p,i / (L*W)
+ real_t getTransportRate() const { return avgTransportRate_; }
+
+ real_t getAverageVelocity() const { return averageVelocity_; }
+
+ private:
+ shared_ptr< ParticleAccessor_T > ac_;
+ real_t normalizationFactor_;
+ uint_t numParticles_;
+ real_t averageVelocity_;
+ real_t avgTransportRate_;
+};
+
+template< typename ParticleAccessor_T >
+Vector3< real_t > getTotalHydrodynamicForceOnParticles(const shared_ptr< ParticleAccessor_T >& ac)
+{
+ Vector3< real_t > totalHydrodynamicForce(0_r);
+ for (uint_t i = uint_t(0); i < ac->size(); ++i)
+ {
+ if (!isSet(ac->getFlags(i), mesa_pd::data::particle_flags::GHOST) &&
+ !isSet(ac->getFlags(i), mesa_pd::data::particle_flags::GLOBAL))
+ {
+ totalHydrodynamicForce += ac->getHydrodynamicForce(i);
+ }
+ }
+ // only reduce to root
+ mpi::reduceInplace(totalHydrodynamicForce, mpi::SUM);
+
+ return totalHydrodynamicForce;
+}
+
+// evaluates slices of solid volume fraction and fill level
+template< typename PdfField_T, typename AntidunesBoundaryHandling_T, typename FlagField_T, typename ScalarField_T >
+class AverageDataSliceEvaluator
+{
+ public:
+ AverageDataSliceEvaluator(const shared_ptr< StructuredBlockStorage >& blocks, const ConstBlockDataID& flagFieldID,
+ const ConstBlockDataID& fillFieldID, const ConstBlockDataID& pdfFieldID)
+ : blocks_(blocks), flagFieldID_(flagFieldID), fillFieldID_(fillFieldID), pdfFieldID_(pdfFieldID)
+ {
+ xlen_ = blocks_->getNumberOfXCells();
+ ylen_ = blocks_->getNumberOfYCells();
+ zlen_ = blocks_->getNumberOfZCells();
+
+ x_z_SolidVolumeFraction_ = std::vector< real_t >(xlen_ * zlen_, real_t(0));
+ x_z_FillLevel_ = std::vector< real_t >(xlen_ * zlen_, real_t(0));
+ x_z_VelocityX_ = std::vector< real_t >(xlen_ * zlen_, real_t(0));
+ x_z_FluidCellCount_ = std::vector< uint_t >(xlen_ * zlen_, 0);
+ maxFluidZPos_ = uint_t(0);
+ }
+
+ void operator()()
+ {
+ // erase data
+ std::fill(x_z_SolidVolumeFraction_.begin(), x_z_SolidVolumeFraction_.end(), real_t(0));
+ std::fill(x_z_FillLevel_.begin(), x_z_FillLevel_.end(), real_t(0));
+ std::fill(x_z_VelocityX_.begin(), x_z_VelocityX_.end(), real_t(0));
+ std::fill(x_z_FluidCellCount_.begin(), x_z_FluidCellCount_.end(), 0);
+
+ // fill contributions
+ for (auto block = blocks_->begin(); block != blocks_->end(); ++block)
+ {
+ const PdfField_T* const pdfField = block->getData< const PdfField_T >(pdfFieldID_);
+ const FlagField_T* const flagField = block->getData< const FlagField_T >(flagFieldID_);
+ const ScalarField_T* const fillField = block->getData< const ScalarField_T >(fillFieldID_);
+
+ const auto solidMO = flagField->getFlag(AntidunesBoundaryHandling_T::movingObstacleFlagID);
+ const auto solidNoSlip = flagField->getFlag(AntidunesBoundaryHandling_T::noSlipFlagID);
+
+ CellInterval xyz = flagField->xyzSize();
+ Cell globalCell;
+
+ maxFluidZPos_ = uint_t(0);
+
+ // iterate all (inner) cells in the field
+ for (auto cell = xyz.begin(); cell != xyz.end(); ++cell)
+ {
+ blocks_->transformBlockLocalToGlobalCell(globalCell, *block, *cell);
+ auto entryIdx = uint_c(globalCell.x()) + uint_c(globalCell.z()) * xlen_;
+ if (flagField->isFlagSet(*cell, solidMO) || flagField->isFlagSet(*cell, solidNoSlip))
+ {
+ x_z_SolidVolumeFraction_[entryIdx] += real_t(1);
+ x_z_FillLevel_[entryIdx] += real_t(1);
+ }
+ else
+ {
+ auto fillLevel = fillField->get(*cell);
+ x_z_FillLevel_[entryIdx] += fillLevel;
+ if (fillLevel > 0_r)
+ {
+ x_z_VelocityX_[entryIdx] += pdfField->getVelocity(*cell)[0];
+ ++x_z_FluidCellCount_[entryIdx];
+
+ maxFluidZPos_ = std::max(uint_t(globalCell.z()), maxFluidZPos_);
+ }
+ }
+ }
+ }
+
+ // reduce this information to the root process
+ mpi::reduceInplace(x_z_SolidVolumeFraction_, mpi::SUM);
+ mpi::reduceInplace(x_z_FillLevel_, mpi::SUM);
+ mpi::reduceInplace(x_z_VelocityX_, mpi::SUM);
+ mpi::reduceInplace(x_z_FluidCellCount_, mpi::SUM);
+ mpi::reduceInplace(maxFluidZPos_, mpi::MAX);
+
+ // normalize
+ for (uint_t i = 0; i < x_z_VelocityX_.size(); ++i)
+ {
+ if (x_z_FluidCellCount_[i] > 0) x_z_VelocityX_[i] /= real_c(x_z_FluidCellCount_[i]);
+ }
+ real_t invNumYCells = 1_r / real_c(ylen_);
+ std::for_each(x_z_SolidVolumeFraction_.begin(), x_z_SolidVolumeFraction_.end(),
+ [invNumYCells](real_t& n) { n *= invNumYCells; });
+ std::for_each(x_z_FillLevel_.begin(), x_z_FillLevel_.end(), [invNumYCells](real_t& n) { n *= invNumYCells; });
+
+ // note: only root process has the complete information!
+ }
+
+ std::vector< real_t >& getSolidVolumeFractionVector() { return x_z_SolidVolumeFraction_; }
+ std::vector< real_t >& getFillLevelVector() { return x_z_FillLevel_; }
+ std::vector< real_t >& getVelocityXVector() { return x_z_VelocityX_; }
+ std::vector< uint_t >& getFluidCellCountVector() { return x_z_FluidCellCount_; }
+ uint_t getMaxFluidZPos() { return maxFluidZPos_; }
+
+ uint_t getXLen() const { return xlen_; }
+ uint_t getYLen() const { return ylen_; }
+ uint_t getZLen() const { return zlen_; }
+
+ private:
+ shared_ptr< StructuredBlockStorage > blocks_;
+ const ConstBlockDataID flagFieldID_;
+ const ConstBlockDataID fillFieldID_;
+ const ConstBlockDataID pdfFieldID_;
+
+ uint_t xlen_;
+ uint_t ylen_;
+ uint_t zlen_;
+ uint_t maxFluidZPos_;
+ std::vector< real_t > x_z_SolidVolumeFraction_;
+ std::vector< real_t > x_z_FillLevel_;
+ std::vector< real_t > x_z_VelocityX_;
+ std::vector< uint_t > x_z_FluidCellCount_;
+};
+
+void writeSphereInformationToFile(const std::string& filename, walberla::mesa_pd::data::ParticleStorage& ps,
+ Vector3< real_t >& domainSize, int precision = 12)
+{
+ std::ostringstream ossData;
+ ossData << std::setprecision(precision);
+
+ WALBERLA_ROOT_SECTION() { ossData << domainSize[0] << " " << domainSize[1] << " " << domainSize[2] << "\n"; }
+
+ for (auto pIt : ps)
+ {
+ using namespace walberla::mesa_pd::data;
+ if (pIt->getBaseShape()->getShapeType() != Sphere::SHAPE_TYPE) continue;
+ using namespace walberla::mesa_pd::data::particle_flags;
+ if (isSet(pIt->getFlags(), GHOST)) continue;
+ auto sp = static_cast< Sphere* >(pIt->getBaseShape().get());
+
+ auto position = pIt->getPosition();
+
+ ossData << pIt->getUid() << " " << position[0] << " " << position[1] << " " << position[2] << " "
+ << sp->getRadius() << '\n';
+ }
+
+ walberla::mpi::writeMPITextFile(filename, ossData.str());
+}
+
+void getAvgDiameterScalingFactor(const std::string& filename, const Vector3< uint_t >& domainSize,
+ const uint_t bedCopiesInX, const uint_t bedCopiesInY, real_t& avgDiameter,
+ real_t& scalingFactor)
+{
+ using namespace walberla;
+
+ std::string textFile;
+
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ifstream t(filename.c_str());
+ if (!t) { WALBERLA_ABORT("Invalid input file " << filename << "\n"); }
+ std::stringstream buffer;
+ buffer << t.rdbuf();
+ textFile = buffer.str();
+
+ std::istringstream fileIss(textFile);
+ std::string line;
+
+ // first line contains generation domain sizes
+ std::getline(fileIss, line);
+ Vector3< real_t > generationDomainSize_SI(0_r);
+ std::istringstream firstLine(line);
+ firstLine >> generationDomainSize_SI[0] >> generationDomainSize_SI[1] >> generationDomainSize_SI[2];
+
+ real_t diameter_SI = 0.0;
+ uint_t numParticles = 0;
+ while (std::getline(fileIss, line))
+ {
+ std::istringstream iss(line);
+
+ mesa_pd::data::ParticleStorage::uid_type uID;
+ mesa_pd::data::ParticleStorage::position_type pos(0_r);
+ walberla::real_t radius = 0;
+ iss >> uID >> pos[0] >> pos[1] >> pos[2] >> radius;
+ WALBERLA_CHECK_GREATER(radius, 0_r, "Invalid radius of " << radius << " found in input file!")
+
+ diameter_SI += 2_r * radius;
+
+ numParticles++;
+ }
+ diameter_SI /= real_t(numParticles);
+
+ scalingFactor = real_c(domainSize[0]) / (generationDomainSize_SI[0] * real_c(bedCopiesInX));
+ avgDiameter = diameter_SI * scalingFactor;
+
+ WALBERLA_CHECK_EQUAL(uint_c(scalingFactor * generationDomainSize_SI[1] * real_c(bedCopiesInY)), domainSize[1],
+ "Error: Generated bed with copies and simulation domain do not match!")
+ }
+
+ walberla::mpi::broadcastObject(scalingFactor);
+ walberla::mpi::broadcastObject(avgDiameter);
+}
+
+void initSpheresFromFile(const std::string& filename, walberla::mesa_pd::data::ParticleStorage& ps,
+ const walberla::mesa_pd::domain::IDomain& domain, walberla::real_t density,
+ const Vector3< uint_t >& domainSize,
+ const std::function< bool(walberla::Vector3< real_t >) >& particleCreateFunction,
+ math::AABB simulationDomain, uint_t bedCopiesInX, uint_t bedCopiesInY, uint_t& numParticles,
+ real_t& maxParticleHeight, const real_t& scalingFactor)
+{
+ using namespace walberla;
+ using namespace walberla::mesa_pd;
+ using namespace walberla::mesa_pd::data;
+
+ auto rank = walberla::mpi::MPIManager::instance()->rank();
+
+ std::string textFile;
+
+ WALBERLA_ROOT_SECTION()
+ {
+ std::ifstream t(filename.c_str());
+ if (!t) { WALBERLA_ABORT("Invalid input file " << filename << "\n"); }
+ std::stringstream buffer;
+ buffer << t.rdbuf();
+ textFile = buffer.str();
+ }
+
+ walberla::mpi::broadcastObject(textFile);
+
+ std::istringstream fileIss(textFile);
+ std::string line;
+
+ // first line contains generation domain sizes
+ std::getline(fileIss, line);
+ Vector3< real_t > generationDomainSize_SI(0_r);
+ std::istringstream firstLine(line);
+ firstLine >> generationDomainSize_SI[0] >> generationDomainSize_SI[1] >> generationDomainSize_SI[2];
+ WALBERLA_LOG_DEVEL_VAR_ON_ROOT(generationDomainSize_SI)
+
+ WALBERLA_CHECK_EQUAL(uint_c(scalingFactor * generationDomainSize_SI[0] * real_c(bedCopiesInX)), domainSize[0],
+ "Error: Generated bed with copies and simulation domain do not match in x!")
+ WALBERLA_CHECK_EQUAL(uint_c(scalingFactor * generationDomainSize_SI[1] * real_c(bedCopiesInY)), domainSize[1],
+ "Error: Generated bed with copies and simulation domain do not match in y!")
+
+ numParticles = 0;
+ maxParticleHeight = 0_r;
+
+ while (std::getline(fileIss, line))
+ {
+ std::istringstream iss(line);
+
+ data::ParticleStorage::uid_type uID;
+ data::ParticleStorage::position_type pos;
+ walberla::real_t radius;
+ iss >> uID >> pos[0] >> pos[1] >> pos[2] >> radius;
+ radius *= scalingFactor;
+
+ for (uint_t copyInYDir = 0; copyInYDir < bedCopiesInY; ++copyInYDir)
+ {
+ for (uint_t copyInXDir = 0; copyInXDir < bedCopiesInX; ++copyInXDir)
+ {
+ auto particlePos = pos;
+
+ particlePos[0] += real_c(copyInXDir) * generationDomainSize_SI[0];
+ particlePos[1] += real_c(copyInYDir) * generationDomainSize_SI[1];
+
+ particlePos *= scalingFactor;
+
+ maxParticleHeight = std::max(maxParticleHeight, particlePos[2] + radius);
+
+ if (!particleCreateFunction(particlePos)) continue;
+
+ WALBERLA_CHECK(simulationDomain.contains(particlePos),
+ "Particle read from file is not contained in simulation domain");
+
+ if (!domain.isContainedInProcessSubdomain(uint_c(rank), particlePos)) continue;
+
+ auto pIt = ps.create();
+ pIt->setPosition(particlePos);
+ pIt->getBaseShapeRef() = std::make_shared< data::Sphere >(radius);
+ pIt->getBaseShapeRef()->updateMassAndInertia(density);
+ pIt->setInteractionRadius(radius);
+ pIt->setOwner(rank);
+ pIt->setType(0);
+
+ numParticles++;
+ }
+ }
+
+ WALBERLA_CHECK_EQUAL(iss.tellg(), -1);
+ }
+ walberla::mpi::allReduceInplace(maxParticleHeight, walberla::mpi::MAX);
+ walberla::mpi::allReduceInplace(numParticles, walberla::mpi::SUM);
+}
+
+void getAverageVelocity(const mesa_pd::data::ParticleAccessorWithBaseShape& ac, real_t& averageVelocity,
+ real_t& maxVelocity, uint_t& numParticles, real_t& maxHeight)
+{
+ averageVelocity = real_t(0);
+ maxVelocity = real_t(0);
+ numParticles = uint_t(0);
+ maxHeight = real_t(0);
+ for (uint_t i = 0; i < ac.size(); ++i)
+ {
+ if (isSet(ac.getFlags(i), walberla::mesa_pd::data::particle_flags::GHOST)) continue;
+ if (isSet(ac.getFlags(i), walberla::mesa_pd::data::particle_flags::GLOBAL)) continue;
+
+ ++numParticles;
+ real_t velMagnitude = ac.getLinearVelocity(i).length();
+ averageVelocity += velMagnitude;
+ maxVelocity = std::max(maxVelocity, velMagnitude);
+ maxHeight = std::max(maxHeight, ac.getPosition(i)[2]);
+ }
+
+ walberla::mpi::allReduceInplace(numParticles, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(averageVelocity, walberla::mpi::SUM);
+ walberla::mpi::allReduceInplace(maxVelocity, walberla::mpi::MAX);
+ walberla::mpi::allReduceInplace(maxHeight, walberla::mpi::MAX);
+
+ averageVelocity /= real_t(numParticles);
+}
+
+auto createPlane(mesa_pd::data::ParticleStorage& ps, const mesa_pd::Vec3& pos, const mesa_pd::Vec3& normal)
+{
+ auto p0 = ps.create(true);
+ p0->setPosition(pos);
+ p0->setBaseShape(std::make_shared< mesa_pd::data::HalfSpace >(normal));
+ p0->getBaseShapeRef()->updateMassAndInertia(real_t(1));
+ p0->setOwner(walberla::mpi::MPIManager::instance()->rank());
+ p0->setType(0);
+ p0->setInteractionRadius(std::numeric_limits< real_t >::infinity());
+ mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::GLOBAL);
+ mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+ mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::FIXED);
+ mesa_pd::data::particle_flags::set(p0->getFlagsRef(), mesa_pd::data::particle_flags::NON_COMMUNICATING);
+ return p0;
+}
+
+} // namespace antidunes
+} // namespace walberla
diff --git a/apps/showcases/Antidunes/pyvista.ipynb b/apps/showcases/Antidunes/pyvista.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..6cb28649aca82ce6aa152f0b56a3a49552f6c00a
--- /dev/null
+++ b/apps/showcases/Antidunes/pyvista.ipynb
@@ -0,0 +1,827 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0953882f",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pyvista as pv\n",
+ "import numpy as np\n",
+ "\n",
+ "import os\n",
+ "import sys"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "81395353",
+ "metadata": {},
+ "source": [
+ "## Plot initial setup"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ddf1b644",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fluid_file = \"/simdata/on74yces/2386423/vtk-out/fluid_field/simulation_step_0.vtu\"\n",
+ "particle_file = \"/simdata/on74yces/2386423/vtk-out/particles/simulation_step_0.vtu\"\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025 \n",
+ "dt = 1.38e-5\n",
+ "\n",
+ "fluid_reader = pv.get_reader(fluid_file)\n",
+ "fluid_mesh = fluid_reader.read()\n",
+ "print(\"fluid_mesh:\", fluid_mesh.array_names)\n",
+ "\n",
+ "particle_reader = pv.get_reader(particle_file)\n",
+ "particle_mesh = particle_reader.read()\n",
+ "print(\"particle_mesh:\", particle_mesh.array_names)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "2b7bb742",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# clip meshes\n",
+ "fluid_mesh = fluid_mesh.clip(normal='x', origin=(2700,59,41.5), invert=False)\n",
+ "particle_mesh = particle_mesh.clip(normal='x', origin=(2700,59,41.5), invert=False)\n",
+ "\n",
+ "# remove gas data from fluid_mesh\n",
+ "fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+ "\n",
+ "# convert velocity to SI units\n",
+ "velocity_unit_conversion = dx / dt\n",
+ "velocity_si = fluid_mesh.get_array(\"velocity\") * velocity_unit_conversion\n",
+ "fluid_mesh.add_field_data(velocity_si, \"velocity_si\")\n",
+ "\n",
+ "# convert particle radius to diameter in SI units\n",
+ "diameter_si = particle_mesh.point_data[\"radius\"] * 2 * dx\n",
+ "particle_mesh.point_data[\"diameter_si\"] = diameter_si\n",
+ "\n",
+ "# create glyphs for particles\n",
+ "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
+ "\n",
+ "# add box for outlining the domain\n",
+ "domain_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 160))\n",
+ "\n",
+ "# add box for outlining fixed particles\n",
+ "particles_fixed_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 17.3943))\n",
+ "\n",
+ "# add box for outlining the bed height\n",
+ "particles_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 48.0002))\n",
+ "\n",
+ "# add box for outlining the liquid height\n",
+ "liquid_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 82.3888))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b8a6aff5",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pl = pv.Plotter(lighting='three lights')\n",
+ "pl.add_mesh(fluid_mesh, scalars=\"velocity_si\", component=0)\n",
+ "pl.add_mesh(sphere_glyphs, scalars=\"diameter_si\", cmap=\"Greys\")\n",
+ "pl.add_mesh(domain_box.outline(), line_width=3, color=\"black\")\n",
+ "pl.add_mesh(particles_fixed_box.outline(), line_width=1, color=\"black\")\n",
+ "pl.add_mesh(particles_box.outline(), line_width=1, color=\"black\")\n",
+ "pl.add_mesh(liquid_box.outline(), line_width=1, color=\"black\")\n",
+ "\n",
+ "pl.view_isometric()\n",
+ "pl.enable_parallel_projection()\n",
+ "#pl.camera.roll += 0\n",
+ "pl.camera.elevation -= 15\n",
+ "pl.camera.azimuth -= 90\n",
+ "pl.remove_scalar_bar(\"velocity_si\")\n",
+ "pl.remove_scalar_bar(\"diameter_si\")\n",
+ "pl.set_background('white') \n",
+ "#pl.show_axes()\n",
+ "pl.camera.zoom(1.75)\n",
+ "#pl.screenshot(\"/home/rzlin/ca36xymo/setup.png\", transparent_background=False, window_size=(2560, 1440))\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "5d1248ea",
+ "metadata": {},
+ "source": [
+ "## Plot velocity field (side view)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "6ca340e6",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fluid_file = \"/simdata/on74yces/2519434/vtk-out/fluid_field/simulation_step_4488000.vtu\"\n",
+ "particle_file = \"/simdata/on74yces/2519434/vtk-out/particles/simulation_step_4488000.vtu\"\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025 \n",
+ "dt = 1.38e-5\n",
+ "\n",
+ "fluid_reader = pv.get_reader(fluid_file)\n",
+ "fluid_mesh = fluid_reader.read()\n",
+ "print(\"fluid_mesh:\", fluid_mesh.array_names)\n",
+ "\n",
+ "particle_reader = pv.get_reader(particle_file)\n",
+ "particle_mesh = particle_reader.read()\n",
+ "print(\"particle_mesh:\", particle_mesh.array_names)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f41650a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# clip meshes\n",
+ "fluid_mesh = fluid_mesh.clip(normal='x', origin=(2200,59,41.5), invert=False)\n",
+ "particle_mesh = particle_mesh.clip(normal='x', origin=(2200,59,41.5), invert=False)\n",
+ "fluid_mesh = fluid_mesh.clip(normal='x', origin=(3000,59,41.5), invert=True)\n",
+ "particle_mesh = particle_mesh.clip(normal='x', origin=(3000,59,41.5), invert=True)\n",
+ "\n",
+ "# remove gas data from fluid_mesh\n",
+ "fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+ "\n",
+ "# convert velocity to SI units\n",
+ "velocity_unit_conversion = dx / dt\n",
+ "velocity_si = fluid_mesh.get_array(\"velocity\") * velocity_unit_conversion\n",
+ "fluid_mesh.add_field_data(velocity_si, \"velocity_si\")\n",
+ "\n",
+ "# convert particle radius to diameter in SI units\n",
+ "diameter_si = particle_mesh.point_data[\"radius\"] * 2 * dx\n",
+ "particle_mesh.point_data[\"diameter_si\"] = diameter_si\n",
+ "\n",
+ "# create glyphs for particles\n",
+ "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
+ "\n",
+ "# # add box for outlining the domain\n",
+ "# domain_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 160))\n",
+ "\n",
+ "# # add box for outlining fixed particles\n",
+ "# particles_fixed_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 17.3943))\n",
+ "\n",
+ "# # add box for outlining the bed height\n",
+ "# particles_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 48.0002))\n",
+ "\n",
+ "# # add box for outlining the liquid height\n",
+ "# liquid_box = pv.Box(bounds=(2700, 3200, 0, 60, 0, 82.3888))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3647ceda",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pl = pv.Plotter(lighting='three lights')\n",
+ "pl.add_mesh(fluid_mesh, scalars=\"velocity_si\", component=0)\n",
+ "pl.add_mesh(sphere_glyphs, scalars=\"diameter_si\", cmap=\"Greys\")\n",
+ "# pl.add_mesh(domain_box.outline(), line_width=3, color=\"black\")\n",
+ "# pl.add_mesh(particles_fixed_box.outline(), line_width=1, color=\"black\")\n",
+ "# pl.add_mesh(particles_box.outline(), line_width=1, color=\"black\")\n",
+ "# pl.add_mesh(liquid_box.outline(), line_width=1, color=\"black\")\n",
+ "\n",
+ "#pl.view_isometric()\n",
+ "pl.view_xz()\n",
+ "pl.enable_parallel_projection()\n",
+ "#pl.camera.roll += 0\n",
+ "#pl.camera.elevation -= 15\n",
+ "#pl.camera.azimuth -= 90\n",
+ "#pl.remove_scalar_bar(\"velocity_si\")\n",
+ "#pl.remove_scalar_bar(\"diameter_si\")\n",
+ "#pl.set_background('white') \n",
+ "#pl.show_axes()\n",
+ "pl.camera.zoom(3.3)\n",
+ "#pl.screenshot(\"/simdata/ca36xymo/velocity-field.png\", transparent_background=True, window_size=(2560, 720))\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "8c1f1df5",
+ "metadata": {},
+ "source": [
+ "## Plot velocity field (3D view) for graphical abstract"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "070962c0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "fluid_file = \"/simdata/on74yces/2593014/vtk-out/fluid_field/simulation_step_7550001.vtu\"\n",
+ "particle_file = \"/simdata/on74yces/2593014/vtk-out/particles/simulation_step_7550001.vtu\"\n",
+ "surface_file = \"/simdata/on74yces/2593014/mesh-out/simulation_step_0.obj\"\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025 \n",
+ "dt = 1.38e-5\n",
+ "\n",
+ "fluid_reader = pv.get_reader(fluid_file)\n",
+ "fluid_mesh = fluid_reader.read()\n",
+ "print(\"fluid_mesh:\", fluid_mesh.array_names)\n",
+ "\n",
+ "particle_reader = pv.get_reader(particle_file)\n",
+ "particle_mesh = particle_reader.read()\n",
+ "print(\"particle_mesh:\", particle_mesh.array_names)\n",
+ "\n",
+ "surface_reader = pv.get_reader(surface_file)\n",
+ "surface_mesh = surface_reader.read()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "23597fc3",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# remove gas data from fluid_mesh\n",
+ "fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+ "\n",
+ "# remove 5 cells in y-direction so that we can see the particles more clearly \n",
+ "fluid_mesh = fluid_mesh.clip(normal='y', origin=(2200,5,41.5), invert=False)\n",
+ "\n",
+ "# remove front part of the so that the figure starts from a low point (looks better)\n",
+ "# 2958 in particle field so that particle holes are not visible in the fluid mesh\n",
+ "fluid_mesh_back = fluid_mesh.clip(normal='x', origin=(2950,5,41.5), invert=True)\n",
+ "particle_mesh_back = particle_mesh.clip(normal='x', origin=(2958,59,41.5), invert=True) \n",
+ "\n",
+ "# store front part for appending it at the back (to maintain periodicity) \n",
+ "fluid_mesh_front = fluid_mesh.clip(normal='x', origin=(2950,5,41.5), invert=False)\n",
+ "particle_mesh_front = particle_mesh.clip(normal='x', origin=(2950,59,41.5), invert=False) # 2958 not necessary here\n",
+ "\n",
+ "# move front part to the back\n",
+ "transform_matrix = np.array([[1, 0, 0, -3200],\n",
+ " [0, 1, 0, 0],\n",
+ " [0, 0, 1, 0],\n",
+ " [0, 0, 0, 1]])\n",
+ "fluid_mesh_front = fluid_mesh_front.transform(transform_matrix)\n",
+ "particle_mesh_front = particle_mesh_front.transform(transform_matrix)\n",
+ "\n",
+ "# create glyphs for particles\n",
+ "sphere_glyphs_back = particle_mesh_back.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
+ "sphere_glyphs_front = particle_mesh_front.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
+ "\n",
+ "# convert velocity to SI units\n",
+ "velocity_unit_conversion = dx / dt\n",
+ "# IMPORTANT: \n",
+ "# The new velocity field (velocity_si) is not clipped if it is part of fluid_mesh.\n",
+ "# Therefore, we have to do the unit conversion for every clipped mesh individually.\n",
+ "fluid_mesh_back.add_field_data(fluid_mesh_back.get_array(\"velocity\") * velocity_unit_conversion, \"velocity_si\")\n",
+ "fluid_mesh_front.add_field_data(fluid_mesh_front.get_array(\"velocity\") * velocity_unit_conversion, \"velocity_si\")\n",
+ "\n",
+ "# get data range of the x-velocity from both fields for using it in the color bar\n",
+ "velocity_field_back_x_si = (fluid_mesh_back.field_data[\"velocity_si\"])[:,0]\n",
+ "velocity_field_front_x_si = (fluid_mesh_front.field_data[\"velocity_si\"])[:,0]\n",
+ "range_velocity_x_si = [min(np.min(velocity_field_back_x_si), np.min(velocity_field_front_x_si)),\n",
+ " max(np.max(velocity_field_back_x_si), np.max(velocity_field_front_x_si))]\n",
+ "\n",
+ "# convert particle radius to diameter in SI units\n",
+ "sphere_glyphs_back.point_data[\"diameter_si\"] = sphere_glyphs_back.point_data[\"radius\"] * 2 * dx\n",
+ "sphere_glyphs_front.point_data[\"diameter_si\"] = sphere_glyphs_front.point_data[\"radius\"] * 2 * dx"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "d0a2b037",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "pl = pv.Plotter(lighting='three lights')\n",
+ "pl.add_mesh(fluid_mesh_back, scalars=\"velocity_si\", component=0, show_scalar_bar=False)\n",
+ "pl.add_mesh(fluid_mesh_front, scalars=\"velocity_si\", component=0, show_scalar_bar=False)\n",
+ "pl.add_scalar_bar(title=\"Downstream fluid velocity / m*s\", vertical=False, width=0.15, height=0.05,\n",
+ " position_x=0.7, position_y=0.6, color='black', title_font_size=40, label_font_size=25)\n",
+ "pl.update_scalar_bar_range(clim=[range_velocity_x_si[0],range_velocity_x_si[1]])\n",
+ "\n",
+ "pl.add_mesh(sphere_glyphs_back, scalars=\"diameter_si\", cmap=\"Greys\", show_scalar_bar=False)\n",
+ "pl.add_mesh(sphere_glyphs_front, scalars=\"diameter_si\", cmap=\"Greys\", show_scalar_bar=False)\n",
+ "pl.add_scalar_bar(title=\"Particle diameter / m\", vertical=False, width=0.15, height=0.05,\n",
+ " position_x=0.3, position_y=0.2, color='black', title_font_size=40, label_font_size=25, \n",
+ " n_labels=3, fmt=\"%.4f\")\n",
+ "pl.update_scalar_bar_range(clim=[0.0026,0.0032])\n",
+ "\n",
+ "pl.view_isometric()\n",
+ "#pl.enable_parallel_projection()\n",
+ "pl.camera.roll -= -5\n",
+ "pl.camera.elevation -= 35\n",
+ "pl.camera.azimuth -= 50\n",
+ "pl.set_background('white') \n",
+ "#pl.show_axes()\n",
+ "pl.camera.zoom(7)\n",
+ "#pl.screenshot(\"/home/rzlin/ca36xymo/velocity-3d.png\", transparent_background=False, window_size=(2560, 1440))\n",
+ "pl.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "adc126de",
+ "metadata": {},
+ "source": [
+ "## Plot full velocity field at each sampling point (for creating the animations)\n",
+ "\n",
+ "### Warning: there are severe bugs in pyvista that overcomplicate this script:\n",
+ "1. memory leak in pl.screenshot(); memory leak is supposed to be circumvented with pl.clear();\n",
+ "2. HOWEVER: when pl.clear() was called, pl.view_xz() somehow disrupts pl.camera.zoom() such that the latter command does not work anymore\n",
+ "\n",
+ "=> new strategy: call python from bash such that the kernel restarts for every file; this avoids the memory leak"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e4b4d571",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def get_image_array(pl, fluid_file_path, particle_file_path, dx, dt, color_bar_limits_velocity):\n",
+ " fluid_reader = pv.get_reader(fluid_file_path)\n",
+ " fluid_mesh = fluid_reader.read()\n",
+ " \n",
+ " # remove gas data from fluid_mesh\n",
+ " fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+ "\n",
+ " particle_reader = pv.get_reader(particle_file_path)\n",
+ " particle_mesh = particle_reader.read()\n",
+ " \n",
+ " # convert velocity to SI units\n",
+ " velocity_unit_conversion = dx / dt\n",
+ " velocity_si = fluid_mesh.get_array(\"velocity\") * velocity_unit_conversion\n",
+ " fluid_mesh.add_field_data(velocity_si, \"velocity_si\")\n",
+ "\n",
+ " # convert particle radius to diameter in SI units\n",
+ " diameter_si = particle_mesh.point_data[\"radius\"] * 2 * dx\n",
+ " particle_mesh.point_data[\"diameter_si\"] = diameter_si\n",
+ "\n",
+ " # create glyphs for particles\n",
+ " sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
+ "\n",
+ " \n",
+ " pl.add_mesh(fluid_mesh, scalars=\"velocity_si\", component=0, show_scalar_bar=False)\n",
+ " pl.add_scalar_bar(title=\"Downstream fluid velocity / m*s\", vertical=False, width=0.15, height=0.15,\n",
+ " position_x=0.325, position_y=0.8, color='black', title_font_size=70, label_font_size=50)\n",
+ " pl.update_scalar_bar_range(clim=color_bar_limits_velocity)\n",
+ " pl.add_mesh(sphere_glyphs, scalars=\"diameter_si\", cmap=\"Greys\", show_scalar_bar=False)\n",
+ " pl.add_scalar_bar(title=\"Particle diameter / m\", vertical=False, width=0.15, height=0.15,\n",
+ " position_x=0.525, position_y=0.8, color='black', title_font_size=70, label_font_size=50)\n",
+ " pl.update_scalar_bar_range(clim=[0.0026,0.0032])\n",
+ " \n",
+ " fluid_mesh.clear_field_data()\n",
+ " particle_mesh.clear_field_data()\n",
+ " sphere_glyphs.clear_field_data()\n",
+ "\n",
+ " pl.view_xz()\n",
+ " pl.enable_parallel_projection()\n",
+ " #pl.remove_scalar_bar(\"velocity\")\n",
+ " #pl.remove_scalar_bar(\"radius\")\n",
+ " pl.set_background('white')\n",
+ " pl.camera.zoom(10)\n",
+ " image_array = pl.screenshot(transparent_background=False, window_size=(10000, 1000))\n",
+ " \n",
+ " # returns image as numpy array\n",
+ " return image_array"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "17d29056",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import os\n",
+ "from PIL import Image\n",
+ "import gc\n",
+ "\n",
+ "fluid_file_basepath = \"/simdata/on74yces/2386423/vtk-out/fluid_field\"\n",
+ "particle_file_basepath = \"/simdata/on74yces/2386423/vtk-out/particles\"\n",
+ "output_file_basepath = \"/simdata/ca36xymo/antidunes/test\"\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025 \n",
+ "dt = 1.38e-5 # E1\n",
+ "#dt = 1.08125e-5 # E4\n",
+ "\n",
+ "# empirically set after viewing data in ParaView (in lattice units)\n",
+ "color_bar_limits_velocity = np.array([-0.3, 1]) # E1 with e_dry=0.97\n",
+ "color_bar_limits_velocity = np.array([-0.3, 1.4]) # E4 with e_dry=0.97\n",
+ "\n",
+ "fluid_file_directory = os.fsencode(fluid_file_basepath)\n",
+ "\n",
+ "# create plotter outside of loop to avoid memoryleak in pyvista.screenshot()\n",
+ "# (https://github.com/pyvista/pyvista/issues/2252#issuecomment-1241929793)\n",
+ "pl = pv.Plotter(lighting='three lights')\n",
+ "\n",
+ "for file in os.listdir(fluid_file_directory):\n",
+ " filename = os.fsdecode(file)\n",
+ " \n",
+ " fluid_file_path = fluid_file_basepath + '/' + filename\n",
+ " particle_file_path = particle_file_basepath + '/' + filename\n",
+ " \n",
+ " filenumber = filename.split('_')\n",
+ " filenumber = filenumber[2].split('.')\n",
+ " filenumber = filenumber[0]\n",
+ " \n",
+ " output_file_path = output_file_basepath + '/' + filenumber + '.jpg'\n",
+ " \n",
+ " image_array = get_image_array(pl, fluid_file_path, particle_file_path, dx, dt, color_bar_limits_velocity)\n",
+ "\n",
+ " pl.clear_actors()\n",
+ "\n",
+ " image = Image.fromarray(image_array)\n",
+ " \n",
+ " image = image.crop((0, 0, image.size[0], 700)) # coordinate systems starts from top left\n",
+ " image.save(output_file_path)\n",
+ " image.close()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b9a63e3b",
+ "metadata": {},
+ "source": [
+ "### Bash script solution to circumvent memory leak bug"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a5b66e66",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# this Python script (named create_image.py) must be called by the bash script\n",
+ "\n",
+ "#!/usr/bin/env python3\n",
+ "\n",
+ "import os\n",
+ "import sys\n",
+ "from PIL import Image\n",
+ "import pyvista as pv\n",
+ "import numpy as np\n",
+ "\n",
+ "pv.start_xvfb()\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025\n",
+ "dt = 1.38e-5 # E1\n",
+ "#dt = 1.08125e-5 # E4\n",
+ "\n",
+ "# empirically set after viewing data in ParaView (in lattice units)\n",
+ "color_bar_limits_velocity = np.array([-0.3, 1]) # E1 with e_dry=0.97\n",
+ "#color_bar_limits_velocity = np.array([-0.3, 1.4]) # E4 with e_dry=0.97\n",
+ "\n",
+ "# path to fluid file must be given as first command line argument\n",
+ "fluid_file_path = str(sys.argv[1])\n",
+ "particle_file_path = fluid_file_path.replace(\"fluid_field\", \"particles\")\n",
+ "\n",
+ "filenumber = fluid_file_path.split('_')[-1]\n",
+ "filenumber = filenumber.split('.')[0]\n",
+ "filenumber = str(filenumber).rjust(10, '0')\n",
+ "output_file_path = \"/simdata/ca36xymo/antidunes/animation-e1/\" + filenumber + \".jpg\"\n",
+ "\n",
+ "fluid_reader = pv.get_reader(fluid_file_path)\n",
+ "fluid_mesh = fluid_reader.read()\n",
+ "\n",
+ "particle_reader = pv.get_reader(particle_file_path)\n",
+ "particle_mesh = particle_reader.read()\n",
+ "\n",
+ "# remove gas data from fluid_mesh\n",
+ "fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+ "\n",
+ "# convert velocity to SI units\n",
+ "velocity_unit_conversion = dx / dt\n",
+ "velocity_si = fluid_mesh.get_array(\"velocity\") * velocity_unit_conversion\n",
+ "fluid_mesh.add_field_data(velocity_si, \"velocity_si\")\n",
+ "\n",
+ "# convert particle radius to diameter in SI units\n",
+ "diameter_si = particle_mesh.point_data[\"radius\"] * 2 * dx\n",
+ "particle_mesh.point_data[\"diameter_si\"] = diameter_si\n",
+ "\n",
+ "# create glyphs for particles\n",
+ "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1,\n",
+ " geom=pv.Sphere(radius=1, theta_resolution=150, phi_resolution=150))\n",
+ "\n",
+ "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+ "\n",
+ "pl.add_mesh(fluid_mesh, scalars=\"velocity_si\", component=0, show_scalar_bar=False)\n",
+ "pl.add_scalar_bar(title=\"Downstream fluid velocity / m*s\", vertical=False, width=0.15, height=0.15,\n",
+ " position_x=0.325, position_y=0.8, color='black', title_font_size=70, label_font_size=50)\n",
+ "pl.update_scalar_bar_range(clim=color_bar_limits_velocity)\n",
+ "pl.add_mesh(sphere_glyphs, scalars=\"diameter_si\", cmap=\"Greys\", show_scalar_bar=False)\n",
+ "pl.add_scalar_bar(title=\"Particle diameter / m\", vertical=False, width=0.15, height=0.15,\n",
+ " position_x=0.525, position_y=0.8, color='black', title_font_size=70, label_font_size=50)\n",
+ "pl.update_scalar_bar_range(clim=[0.0026,0.0032])\n",
+ "\n",
+ "pl.view_xz()\n",
+ "#pl.enable_parallel_projection()\n",
+ "pl.set_background('white')\n",
+ "#pl.camera.zoom(10)\n",
+ "\n",
+ "# specify camera so that image stays at a fixed position\n",
+ "pl.camera.clipping_range = (5812.79, 6720.05)\n",
+ "pl.camera.distance = 6207.82\n",
+ "pl.camera.focal_point = (1600.44, 29.9889, 38.3558)\n",
+ "pl.camera.parallel_projection = True\n",
+ "pl.camera.parallel_scale = 160.67\n",
+ "pl.camera.thickness = 907.259\n",
+ "pl.camera.view_angle = 30\n",
+ "\n",
+ "pl.window_size = [10000, 1000]\n",
+ "pl.screenshot(output_file_path,transparent_background=False)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "537dc005",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "#!/bin/bash -l\n",
+ "\n",
+ "conda activate\n",
+ "\n",
+ "fluid_file_basepath1=\"/simdata/on74yces/2502347/vtk-out/fluid_field/\"\n",
+ "# fluid_file_basepath2=\"/simdata/on74yces/2504475/vtk-out/fluid_field/\"\n",
+ "# fluid_file_basepath3=\"/simdata/on74yces/2509267/vtk-out/fluid_field/\"\n",
+ "# fluid_file_basepath4=\"/simdata/on74yces/2519434/vtk-out/fluid_field/\"\n",
+ "\n",
+ "declare -a PathList=($fluid_file_basepath1 $fluid_file_basepath2 $fluid_file_basepath3 $fluid_file_basepath4)\n",
+ "\n",
+ "for path in ${PathList[@]}; do\n",
+ " for file in $path*.vtu; do\n",
+ " python3 create-image.py $file\n",
+ " # echo $file\n",
+ " done\n",
+ "done"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "ebe75f76",
+ "metadata": {},
+ "source": [
+ "## Plot full bed height elevation at each sampling point (for creating the animations)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "9254f001",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# this Python script (named create_image.py) must be called by the bash script\n",
+ "\n",
+ "#!/usr/bin/env python3\n",
+ "\n",
+ "import os\n",
+ "import sys\n",
+ "from PIL import Image\n",
+ "import pyvista as pv\n",
+ "import numpy as np\n",
+ "\n",
+ "pv.start_xvfb()\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025\n",
+ "\n",
+ "# path to particle file must be given as first command line argument\n",
+ "particle_file_path = \"/simdata/on74yces/2509268/vtk-out/particles/simulation_step_3799000.vtu\" #str(sys.argv[1])\n",
+ "\n",
+ "filenumber = particle_file_path.split('_')[-1]\n",
+ "filenumber = filenumber.split('.')[0]\n",
+ "filenumber = str(filenumber).rjust(10, '0')\n",
+ "output_file_path = \"/simdata/ca36xymo/antidunes/animation-particles-e4/\" + filenumber + \".jpg\"\n",
+ "\n",
+ "particle_reader = pv.get_reader(particle_file_path)\n",
+ "particle_mesh = particle_reader.read()\n",
+ "\n",
+ "# add z-coordinate of particle center in SI units\n",
+ "particle_mesh.point_data[\"center_coordinate_z_si\"] = particle_mesh.points[:,2] * dx\n",
+ "\n",
+ "# create glyphs for particles\n",
+ "sphere_glyphs = particle_mesh.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=30, phi_resolution=30))\n",
+ "\n",
+ "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+ "\n",
+ "# add mesh at periodic sides to enlarge bed region\n",
+ "transform_matrix = np.array([[1, 0, 0, 0],\n",
+ " [0, 1, 0, -60],\n",
+ " [0, 0, 1, 0],\n",
+ " [0, 0, 0, 1]])\n",
+ "particle_mesh_transformed_once = particle_mesh.transform(transform_matrix,inplace=False)\n",
+ "sphere_glyphs_transformed_once = particle_mesh_transformed_once.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=30, phi_resolution=30))\n",
+ "particle_mesh_transformed_twice = particle_mesh_transformed_once.transform(transform_matrix,inplace=False)\n",
+ "sphere_glyphs_transformed_twice = particle_mesh_transformed_twice.glyph(scale=\"radius\", factor=1, \n",
+ " geom=pv.Sphere(radius=1, theta_resolution=30, phi_resolution=30))\n",
+ "\n",
+ "pl.add_mesh(sphere_glyphs_transformed_once, scalars=\"center_coordinate_z_si\", show_scalar_bar=False)\n",
+ "pl.update_scalar_bar_range(clim=[0.008,0.012])\n",
+ "pl.add_mesh(sphere_glyphs, scalars=\"center_coordinate_z_si\", show_scalar_bar=False)\n",
+ "pl.update_scalar_bar_range(clim=[0.008,0.012])\n",
+ "pl.add_mesh(sphere_glyphs_transformed_twice, scalars=\"center_coordinate_z_si\", show_scalar_bar=False)\n",
+ "pl.add_scalar_bar(title=\"Bed height elevation / m\", vertical=False, width=0.15, height=0.15,\n",
+ " position_x=0.425, position_y=0.7, color='black', title_font_size=70, label_font_size=50,\n",
+ " fmt=\"%.3f\")\n",
+ "pl.update_scalar_bar_range(clim=[0.008,0.012])\n",
+ "\n",
+ "# pl.add_mesh(pv.Line((0,0,100), (3200,0,100)), color=\"lightgray\", line_width=10)\n",
+ "# pl.add_mesh(pv.Line((0,-60,100), (3200,-60,100)), color=\"lightgray\", line_width=10)\n",
+ "\n",
+ "pl.view_xy()\n",
+ "#pl.enable_parallel_projection()\n",
+ "pl.set_background('white')\n",
+ "#pl.camera.zoom(10)\n",
+ "\n",
+ "# specify camera so that image stays at a fixed position\n",
+ "pl.camera.clipping_range = (5812.79, 6720.05)\n",
+ "pl.camera.distance = 6207.82\n",
+ "pl.camera.focal_point = (1600.44, 29.9889, 38.3558)\n",
+ "pl.camera.parallel_projection = True\n",
+ "pl.camera.parallel_scale = 160.67\n",
+ "pl.camera.thickness = 907.259\n",
+ "pl.camera.view_angle = 30\n",
+ "\n",
+ "pl.window_size = [10000, 1000]\n",
+ "pl.screenshot(output_file_path,transparent_background=False)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "73d07351",
+ "metadata": {},
+ "source": [
+ "## Plot full free-surface height elevation at each sampling point (for creating the animations)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "95013d8f",
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "# this Python script (named create_image.py) must be called by the bash script\n",
+ "\n",
+ "#!/usr/bin/env python3\n",
+ "\n",
+ "import os\n",
+ "import sys\n",
+ "from PIL import Image\n",
+ "import pyvista as pv\n",
+ "import numpy as np\n",
+ "\n",
+ "pv.start_xvfb()\n",
+ "\n",
+ "# taken from job output file\n",
+ "dx = 0.00025\n",
+ "\n",
+ "# path to fluid file must be given as first command line argument\n",
+ "fluid_file_path = str(sys.argv[1]) #\"/simdata/on74yces/2517354/vtk-out/fluid_field/simulation_step_7551000.vtu\" \n",
+ "\n",
+ "filenumber = fluid_file_path.split('_')[-1]\n",
+ "filenumber = filenumber.split('.')[0]\n",
+ "filenumber = str(filenumber).rjust(10, '0')\n",
+ "output_file_path = \"/simdata/ca36xymo/antidunes/animation-free-surface-e4/\" + filenumber + \".jpg\"\n",
+ "\n",
+ "fluid_reader = pv.get_reader(fluid_file_path)\n",
+ "fluid_mesh = fluid_reader.read()\n",
+ "\n",
+ "# remove gas data from fluid_mesh\n",
+ "fluid_mesh = fluid_mesh.threshold(value=2, scalars=\"mapped_flag\", invert=True)\n",
+ "\n",
+ "# slice fluid mesh to get real 2D data\n",
+ "fluid_mesh = fluid_mesh.slice(normal=[0, 1, 0])\n",
+ "\n",
+ "# extract surface so we apply the extrusion\n",
+ "fluid_mesh = fluid_mesh.extract_surface()\n",
+ "\n",
+ "# extrude fluid_mesh to 3D\n",
+ "fluid_mesh = fluid_mesh.extrude([0, -60, 0], capping=False)\n",
+ "\n",
+ "# add z-coordinate of cell center in SI units\n",
+ "fluid_mesh.point_data[\"coordinate_z_si\"] = fluid_mesh.points[:,2] * dx\n",
+ "\n",
+ "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+ "\n",
+ "pl.add_mesh(fluid_mesh, scalars=\"coordinate_z_si\", component=0, show_scalar_bar=False)\n",
+ "pl.add_scalar_bar(title=\"Free-surface elevation / m\", vertical=False, width=0.15, height=0.15,\n",
+ " position_x=0.425, position_y=0.7, color='black', title_font_size=70, label_font_size=50,\n",
+ " fmt=\"%.3f\")\n",
+ "#pl.update_scalar_bar_range(clim=[0.018,0.022]) # E1\n",
+ "pl.update_scalar_bar_range(clim=[0.020,0.024]) # E4\n",
+ "\n",
+ "pl.view_xy()\n",
+ "#pl.enable_parallel_projection()\n",
+ "pl.set_background('white')\n",
+ "#pl.camera.zoom(10)\n",
+ "\n",
+ "# specify camera so that image stays at a fixed position\n",
+ "pl.camera.clipping_range = (5812.79, 6720.05)\n",
+ "pl.camera.distance = 6207.82\n",
+ "pl.camera.focal_point = (1600.44, 29.9889, 38.3558)\n",
+ "pl.camera.parallel_projection = True\n",
+ "pl.camera.parallel_scale = 160.67\n",
+ "pl.camera.thickness = 907.259\n",
+ "pl.camera.view_angle = 30\n",
+ "\n",
+ "pl.window_size = [10000, 1000]\n",
+ "\n",
+ "pl.screenshot(output_file_path,transparent_background=False)"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "b516bfd5",
+ "metadata": {},
+ "source": [
+ "## Plot coordinate systems\n",
+ "This is very troublesome to do in the actual plots => add them afterwards using ffmpeg"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e46cfd29",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import pyvista as pv\n",
+ "\n",
+ "output_file_path = \"/simdata/ca36xymo/antidunes/animation-free-surface-e1/coordinate-system.jpg\"\n",
+ "\n",
+ "pl = pv.Plotter(lighting='three lights', off_screen=True)\n",
+ "pl.view_xz()\n",
+ "pl.add_axes(color=\"black\", x_color=\"black\", y_color=\"black\", z_color=\"black\", viewport=(0, 0, 1, 1))\n",
+ "pl.set_background('white')\n",
+ "\n",
+ "pl.screenshot(output_file_path,transparent_background=False)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "63082336",
+ "metadata": {},
+ "outputs": [],
+ "source": []
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/apps/showcases/Antidunes/slice_evaluation.ipynb b/apps/showcases/Antidunes/slice_evaluation.ipynb
new file mode 100644
index 0000000000000000000000000000000000000000..422762f92b8eceb90364660bb43a1c07f777a96c
--- /dev/null
+++ b/apps/showcases/Antidunes/slice_evaluation.ipynb
@@ -0,0 +1,1138 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4f094327",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "%matplotlib notebook\n",
+ "import matplotlib as mpl\n",
+ "import numpy as np\n",
+ "import seaborn as sns"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "904ccd9a",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "baseFolder = \"/simdata/on74yces/merge_E1_flat_new\""
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "a7782605",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def readSliceData(fileName):\n",
+ " with open(fileName) as f:\n",
+ " line = f.readline()[2:].split(\" \")\n",
+ " nx = int(line[0])\n",
+ " nz = int(line[1])\n",
+ " svf = np.loadtxt(fileName)\n",
+ " return svf.reshape((nz,nx)).transpose() # x: dim0, z: dim1\n",
+ "\n",
+ "def extractMaxAvailableTimeStep():\n",
+ " import glob\n",
+ "\n",
+ " baseFileName = \"/svfSlice_\"\n",
+ " fileEnding = \".txt\"\n",
+ " allFiles = glob.glob(baseFolder + baseFileName+\"*\"+fileEnding)\n",
+ " maxTimeStep = np.max([int(f[(f.find(baseFileName)+len(baseFileName)):f.find(fileEnding)]) for f in allFiles])\n",
+ " return maxTimeStep\n",
+ "\n",
+ "def getCurrentForcing(t):\n",
+ " fileName = baseFolder+\"/fluidInfo.txt\"\n",
+ " # columns: t, fx, uMean, ...\n",
+ " content = np.loadtxt(fileName).transpose()\n",
+ " idx = np.argwhere(content[0] == t)\n",
+ " return content[1][idx]\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "114b2285",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# source: \n",
+ "# https://stackoverflow.com/questions/26563858/matplotlib-imshow-fixed-aspect-and-vertical-colorbar-matching-master-axis-height\n",
+ "def match_colorbar(cb, ax=None):\n",
+ " \"\"\"\n",
+ " Match the size of the colorbar with the size of the axes.\n",
+ " \n",
+ " Args:\n",
+ " ax: Axes from which the colorbar \"stole\" space.\n",
+ " cb: Colorbar to match to `ax`.\n",
+ " \"\"\"\n",
+ " ax = ax or plt.gca()\n",
+ " bbox = ax.get_position()\n",
+ " cb_bbox = cb.ax.get_position()\n",
+ " if cb.orientation == \"vertical\":\n",
+ " # Update bottom and height.\n",
+ " left = cb_bbox.xmin\n",
+ " width = cb_bbox.width\n",
+ " bottom = bbox.ymin\n",
+ " height = bbox.height\n",
+ " else:\n",
+ " # Update left and width.\n",
+ " left = bbox.xmin\n",
+ " width = bbox.width\n",
+ " bottom = cb_bbox.ymin\n",
+ " height = cb_bbox.height\n",
+ " pos = [left, bottom, width, height]\n",
+ " cb.ax.set_position(pos)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4d487479",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def evalTimeStep(t):\n",
+ " \n",
+ " bedHeightThreshold = 0.3\n",
+ " \n",
+ " timestep = str(t)\n",
+ " fileName = baseFolder+\"/svfSlice_\"+timestep+\".txt\"\n",
+ " svf = readSliceData(fileName)\n",
+ " fileName = baseFolder+\"/fillSlice_\"+timestep+\".txt\"\n",
+ " fill = readSliceData(fileName)\n",
+ " fileName = baseFolder+\"/velXSlice_\"+timestep+\".txt\"\n",
+ " vel = readSliceData(fileName)\n",
+ " \n",
+ " fx = getCurrentForcing(t)\n",
+ " \n",
+ " #print(\"Total water volume: \" , np.sum(fill))\n",
+ " \n",
+ " #plt.figure()\n",
+ " #plt.imshow(svf<bedHeightThreshold, interpolation='none')\n",
+ " #plt.imshow(vel, interpolation='none')\n",
+ " #plt.colorbar()\n",
+ " #plt.show()\n",
+ " \n",
+ " xlen = svf.shape[0]\n",
+ " zlen = svf.shape[1]\n",
+ "\n",
+ "\n",
+ " xPos = np.linspace(0.5, xlen - 0.5, num = xlen)\n",
+ " zPos = np.linspace(0.5, zlen - 0.5, num = zlen)\n",
+ "\n",
+ " waterHeightOverX = np.zeros(xlen)\n",
+ " bedHeightOverX = np.zeros(xlen)\n",
+ " froudeOverX = np.zeros(xlen)\n",
+ " shieldsOverX = np.zeros(xlen)\n",
+ " ReFricOverX = np.zeros(xlen)\n",
+ " \n",
+ " #plt.figure()\n",
+ " \n",
+ " avgUMean = 0\n",
+ " avgWaterHeight = 0\n",
+ "\n",
+ " for x in range(xlen):\n",
+ " svfX = svf[x,:]\n",
+ " fillX = fill[x,:]\n",
+ " velX = vel[x,:]\n",
+ "\n",
+ " waterElevation = np.sum(fillX)\n",
+ " idxWaterElevation = int(waterElevation)\n",
+ " \n",
+ " idxBedHeight = zlen - np.argmax(svfX[::-1]>bedHeightThreshold) - 1\n",
+ " bedHeight = zPos[idxBedHeight] + (bedHeightThreshold - svfX[idxBedHeight]) / (svfX[idxBedHeight+1] - svfX[idxBedHeight]) * 1\n",
+ "\n",
+ " waterHeight = waterElevation - bedHeight\n",
+ " \n",
+ " #print(idxWaterElevation, idxBedHeight, svfX[idxBedHeight+1], svfX[idxBedHeight], svfX[idxBedHeight-1])\n",
+ "\n",
+ " uMean = np.average(velX[idxBedHeight:idxWaterElevation+1])\n",
+ " \n",
+ " #velDerivative = np.gradient(velX)\n",
+ " # velocity derivative is very sensitive to the exact location -> where to take it?\n",
+ " # here: take average over some cells right above the bed\n",
+ " #velDerivativeAtBed = np.average(velDerivative[idxBedHeight+1:idxBedHeight+5]) # TODO check this!!!!\n",
+ " #wallShearStress = dynViscosity * velDerivativeAtBed \n",
+ " # note: this does only get a small part of the total stress on the particle -> different approach necessary\n",
+ " \n",
+ " # integrate current forcing over water height to the position of the bed (see discussion with Bernhard)\n",
+ " wallShearStress = waterHeight * fx\n",
+ " \n",
+ " frictionVelocity = np.sqrt(wallShearStress / fluidDensity)\n",
+ "\n",
+ " waterHeightOverX[x] = waterHeight\n",
+ " bedHeightOverX[x] = bedHeight\n",
+ " froudeOverX[x] = uMean / np.sqrt(waterHeight * gravAcceleration)\n",
+ " shieldsOverX[x] = wallShearStress / ((densityRatio - 1) * gravAcceleration * avgDiameter)\n",
+ " ReFricOverX[x] = frictionVelocity * waterHeight / kinematicViscosity\n",
+ " \n",
+ " \n",
+ " \n",
+ " #plt.plot(velDerivative[idxBedHeight:idxWaterElevation-1])\n",
+ " #plt.plot(velX[idxBedHeight:idxWaterElevation+1])\n",
+ " \n",
+ " avgWaterHeight += waterHeight\n",
+ " avgUMean += uMean\n",
+ " \n",
+ " #plt.xlabel(\"height\")\n",
+ " #plt.ylabel(\"velX\")\n",
+ " #plt.show()\n",
+ " \n",
+ " avgWaterHeight /= float(xlen)\n",
+ " avgUMean /= float(xlen)\n",
+ " \n",
+ " ReBulk = avgUMean * avgWaterHeight / kinematicViscosity\n",
+ " \n",
+ " print(\"Re_b = \", ReBulk)\n",
+ " \n",
+ " return (xPos,waterHeightOverX, bedHeightOverX, froudeOverX, shieldsOverX, ReFricOverX)\n",
+ "\n",
+ "\n",
+ "def plotResult(result):\n",
+ " fig, axes = plt.subplots(nrows=3, ncols=2,sharex=True)\n",
+ " \n",
+ " labels = [\"$h_f$\", \"$h_b$\", \"$Fr$\", \"$Sh$\", \"$Re_\\\\tau$\"]\n",
+ " for i, ax in enumerate(axes.flatten()):\n",
+ " if i >= len(labels):\n",
+ " break\n",
+ " ax.plot(result[0], result[i+1],'-')\n",
+ " ax.set_ylabel(labels[i])\n",
+ " ax.set_xlabel(\"$x$\")\n",
+ " \n",
+ " fig.tight_layout()\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "eeb92b70",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "with open(baseFolder+'/info.txt') as f:\n",
+ " lines = f.readlines()\n",
+ " evalFrequency = int(lines[0])\n",
+ " gravAcceleration = float(lines[1])\n",
+ " dynViscosity = float(lines[2])\n",
+ " densityRatio = float(lines[3])\n",
+ " avgDiameter = float(lines[4])\n",
+ " xSize = int(lines[5])\n",
+ " ySize = int(lines[6])\n",
+ " zSize = int(lines[7])\n",
+ " numParticles = int(lines[8])\n",
+ "\n",
+ "print(\"Eval infos:\")\n",
+ "print(evalFrequency, gravAcceleration, dynViscosity, densityRatio, avgDiameter,numParticles)\n",
+ "maxTimeStep = extractMaxAvailableTimeStep()\n",
+ "\n",
+ "averageVelocity = 0.02\n",
+ "fluidDensity = 1\n",
+ "kinematicViscosity = dynViscosity / fluidDensity\n",
+ "#timeStepSpacing = 50*evalFrequency\n",
+ "timeStepSpacing = 10*evalFrequency\n",
+ "evalTimeSteps = np.arange(0,maxTimeStep,timeStepSpacing)\n",
+ "\n",
+ "lengthRef = avgDiameter\n",
+ "timeRef = avgDiameter / averageVelocity\n",
+ "\n",
+ "print(\"Will evaluate\", len(evalTimeSteps), \"slices, from 0 to\", maxTimeStep, \"with spacing\", timeStepSpacing)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5bdec6b4",
+ "metadata": {
+ "scrolled": false
+ },
+ "outputs": [],
+ "source": [
+ "hbOverTime = np.zeros((len(evalTimeSteps), xSize))\n",
+ "hfOverTime = np.zeros((len(evalTimeSteps), xSize))\n",
+ "for i,t in enumerate(evalTimeSteps):\n",
+ " result = evalTimeStep(t)\n",
+ " \n",
+ " # uncomment if plot per time step required\n",
+ " #plotResult(result)\n",
+ " \n",
+ " hfAvg = np.average(result[1])\n",
+ " hbAvg = np.average(result[2])\n",
+ " FrAvg = np.average(result[3])\n",
+ " ShAvg = np.average(result[4])\n",
+ " print(f\"{i}, t = {t}: hf = {hfAvg:.2f}, hb = {hbAvg:.2f}, Fr = {FrAvg:.3f}, Sh = {ShAvg:.3e}\")\n",
+ " hfOverTime[i] = result[1]\n",
+ " hbOverTime[i] = result[2]\n",
+ " "
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "3587551b",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot 2D x-t-Data evaluation in lattice units\n",
+ "\n",
+ "xPlotLimits = np.array([0,xSize]) / lengthRef\n",
+ "tPlotLimits = np.array([min(evalTimeSteps), max(evalTimeSteps)]) / timeRef\n",
+ "\n",
+ "fluidSurfaceOverTime = hbOverTime+hfOverTime\n",
+ "plt.figure()\n",
+ "plt.imshow(hbOverTime, interpolation='none', origin='lower',\n",
+ " aspect='auto',extent=(*xPlotLimits,*tPlotLimits))\n",
+ "#plt.imshow(fluidSurfaceOverTime[:,:], interpolation='none', origin='lower',aspect='auto')\n",
+ "plt.colorbar()\n",
+ "plt.xlabel(r\"$x / D$\")\n",
+ "plt.ylabel(r\"$t / t_{ref}$\")\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "ccfa2b70",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# values taken from job output file\n",
+ "dx = 0.00025\n",
+ "\n",
+ "# E1\n",
+ "dt = 1.38442e-05\n",
+ "\n",
+ "# E4\n",
+ "#dt = 1.08125e-5"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "c1b25e53",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot 2D x-t-Data evaluation in SI units\n",
+ "\n",
+ "xPlotLimits = np.array([0,xSize]) * dx\n",
+ "tPlotLimits = np.array([min(evalTimeSteps), max(evalTimeSteps)]) * dt\n",
+ "averageHb = np.average(hbOverTime) * dx\n",
+ "\n",
+ "#fluidSurfaceOverTime = (hbOverTime+hfOverTime) * 0.00025\n",
+ "plt.figure()\n",
+ "plt.imshow(hbOverTime[:] * dx - averageHb, interpolation='none', origin='lower',\n",
+ " aspect='auto', extent=(*xPlotLimits,*tPlotLimits), vmin=-5e-3, vmax=5e-3)\n",
+ "#plt.imshow(fluidSurfaceOverTime[:,:], interpolation='none', origin='lower',aspect='auto')\n",
+ "plt.colorbar()\n",
+ "plt.xlabel(r\"$x$ / m\")\n",
+ "plt.ylabel(r\"$t$ / s\")\n",
+ "plt.show()\n",
+ "\n",
+ "print(tPlotLimits)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4223c1cc",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot bed height evaluation in SI units with\n",
+ " # time restricted to 30 - 75 s\n",
+ " # x-domain restricted to 0 - 0.75 m\n",
+ "\n",
+ "timestep_length = evalTimeSteps[1] - evalTimeSteps[0]\n",
+ "\n",
+ "t_min = 30\n",
+ "timestep_index_min = int((t_min / dt) // timestep_length)\n",
+ "\n",
+ "t_max = 75\n",
+ "timestep_index_max = int((t_max / dt) // timestep_length)\n",
+ "\n",
+ "x_min = 0\n",
+ "x_index_min = int(x_min // dx)\n",
+ "\n",
+ "x_max = 0.75\n",
+ "x_index_max = int(x_max // dx)\n",
+ "\n",
+ "xPlotLimits = np.array([x_min,x_max])\n",
+ "tPlotLimits = np.array([t_min, t_max])\n",
+ "averageHb = np.average(hbOverTime[timestep_index_min:timestep_index_max,x_index_min:x_index_max]) * dx\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.imshow(hbOverTime[timestep_index_min:timestep_index_max,x_index_min:x_index_max] * dx - averageHb, interpolation='none', origin='lower',\n",
+ " aspect='auto', extent=(*xPlotLimits,*tPlotLimits), vmin=-5e-3, vmax=5e-3)\n",
+ "\n",
+ "cb = plt.colorbar(label=r\"$h$ / m\")\n",
+ "match_colorbar(cb)\n",
+ "#cb.ax.tick_params(labelsize=8) \n",
+ "plt.xlabel(r\"$x$ / m\")\n",
+ "plt.ylabel(r\"$t$ / s\")\n",
+ "\n",
+ "import tikzplotlib\n",
+ "tikzplotlib.clean_figure()\n",
+ "tikzplotlib.save(\n",
+ " '/home/rzlin/ca36xymo/tikz/bed-elevation-simulation.tex',\n",
+ " axis_height = '\\\\figureheight',\n",
+ " axis_width = '\\\\figurewidth'\n",
+ " )\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "bed63b06",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(np.average(hfOverTime[timestep_index_min:timestep_index_max,x_index_min:x_index_max]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b8c8fc0e",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import scipy.fft\n",
+ "\n",
+ "# plot PSD with\n",
+ " # time restricted to 30 - 75 s\n",
+ " # x-domain restricted to 0 - 0.75 m\n",
+ "\n",
+ "timestep_length = evalTimeSteps[1] - evalTimeSteps[0]\n",
+ "\n",
+ "t_min = 30\n",
+ "timestep_index_min = int((t_min / dt) // timestep_length)\n",
+ "\n",
+ "t_max = 75\n",
+ "timestep_index_max = int((t_max / dt) // timestep_length)\n",
+ "\n",
+ "x_min = 0\n",
+ "x_index_min = int(x_min // dx)\n",
+ "\n",
+ "x_max = 0.75\n",
+ "x_index_max = int(x_max // dx)\n",
+ "\n",
+ "averageHb = np.average(hbOverTime[timestep_index_min:timestep_index_max,x_index_min:x_index_max]) * dx\n",
+ "\n",
+ "hbOverTimeShort = hbOverTime[timestep_index_min:timestep_index_max,x_index_min:x_index_max] * dx - averageHb\n",
+ "\n",
+ "# power spectral density (PSD) on period-wavelength (PW) plane\n",
+ "psd_pw = scipy.fft.fft2(hbOverTimeShort)\n",
+ "psd_pw = np.abs(psd_pw)**2 / (hbOverTimeShort.shape[0]*hbOverTimeShort.shape[1])\n",
+ "\n",
+ "ft_period = np.abs(1 / scipy.fft.fftfreq(np.transpose(hbOverTimeShort).shape[-1],timestep_length)[1:] * dt)\n",
+ "ft_wavelength = np.abs(1 / scipy.fft.fftfreq(hbOverTimeShort.shape[-1],1)[1:] * dx)\n",
+ "\n",
+ "plt.figure()\n",
+ "cs = plt.contour(ft_period, ft_wavelength, np.transpose(np.abs(psd_pw[1:,1:])), levels=[0.01], colors='tab:orange',linewidths=6)\n",
+ "plt.xlabel(r\"$T$ / s\")\n",
+ "plt.ylabel(r\"$\\lambda$ / m\")\n",
+ "ax = plt.gca()\n",
+ "ax.set_xlim([5, 70])\n",
+ "ax.set_ylim([0.05, 0.2])\n",
+ "#cb = plt.colorbar(label=r\"PSD / m$^2$\")\n",
+ "\n",
+ "\n",
+ "## uncomment the following things to only store the png (no axes etc.)\n",
+ "## => tikz image must be created manually later\n",
+ "# plt.subplots_adjust(bottom = 0)\n",
+ "# plt.subplots_adjust(top = 1)\n",
+ "# plt.subplots_adjust(right = 1)\n",
+ "# plt.subplots_adjust(left = 0)\n",
+ "# plt.axis('off')\n",
+ "# plt.savefig('/home/rzlin/ca36xymo/tikz/spectral-density-simulation.png',bbox_inches='tight',transparent=True, pad_inches=0)\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "afd23ef9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "psd_pw_new = np.transpose(np.abs(psd_pw[1:,1:]))\n",
+ "\n",
+ "# compute celerity (for each point in wavelength-period PSD-array)\n",
+ "# => identify all possible y-axis values for new celerity-wavelength PSD-array\n",
+ "c = []\n",
+ "x = 0\n",
+ "for wavelength in ft_wavelength:\n",
+ " y = 0\n",
+ " for period in ft_period:\n",
+ " if psd_pw_new[x,y] > 0.005:\n",
+ " c.append(wavelength / period)\n",
+ " y += 1\n",
+ " x += 1\n",
+ " \n",
+ "c = list(set(c)) # remove duplicate entries from list\n",
+ "c.sort() # sort in ascending order\n",
+ "# list 'c' contains all possible celerities in ascending order without duplicates\n",
+ "# => this will be the x-axis of the celerity-wavelength PSD-array \n",
+ "\n",
+ "# create celerity-wavelength PSD-array \n",
+ "psd_cw = np.zeros((ft_wavelength.shape[0], len(c)))\n",
+ "\n",
+ "# fill this array with PSD at correct position\n",
+ "for i_wavelength in range(0, ft_wavelength.shape[0]):\n",
+ " for i_period in range(0, ft_period.shape[0]):\n",
+ " if psd_pw_new[i_wavelength,i_period] > 0.005:\n",
+ " celerity = ft_wavelength[i_wavelength] / ft_period[i_period] # compute celerity\n",
+ " i_celerity = c.index(celerity) # map celerity to it's index on the x-axis of the new array\n",
+ " psd_cw[i_wavelength,i_celerity] = abs(psd_pw_new[i_wavelength,i_period]) # store PSD at correct position in new array\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.contourf(ft_wavelength, np.asarray(c), np.abs(np.transpose(psd_cw)),colors='tab:orange')\n",
+ "#plt.colorbar()\n",
+ "plt.xlabel(r\"$\\lambda$ / m\")\n",
+ "plt.ylabel(r\"$c$ / m/s\")\n",
+ "ax = plt.gca()\n",
+ "ax.set_xlim([0.05, 0.2])\n",
+ "ax.set_ylim([0, 0.015])\n",
+ "\n",
+ "# # uncomment the following things to only store the png (no axes etc.)\n",
+ "# # => tikz image must be created manually later\n",
+ "plt.subplots_adjust(bottom = 0)\n",
+ "plt.subplots_adjust(top = 1)\n",
+ "plt.subplots_adjust(right = 1)\n",
+ "plt.subplots_adjust(left = 0)\n",
+ "plt.axis('off')\n",
+ "plt.savefig('/home/rzlin/ca36xymo/tikz/celerity-simulation-e1.png',bbox_inches='tight',transparent=True, pad_inches=0)\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "e8146b19",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# FFT along x per time frame -> find wavelength\n",
+ "\n",
+ "import scipy.fft\n",
+ "\n",
+ "N = hbOverTime.shape[-1] # = number of samples, here xSize\n",
+ "colors = sns.color_palette(\"crest\",hbOverTime.shape[0] )\n",
+ "plt.figure()\n",
+ "\n",
+ "for i,t in enumerate(evalTimeSteps):\n",
+ "\n",
+ " yf = scipy.fft.rfft(hbOverTime[i])\n",
+ " sampleSpacing = 1 # = dx\n",
+ " \n",
+ " # discard first element (is 0 frequency -> offset of sinus)\n",
+ " xf = scipy.fft.rfftfreq(N,sampleSpacing)[1:]\n",
+ " yfMod = (2.0 / N * np.abs(yf))[1:]\n",
+ " plt.plot(1/xf, yfMod, color=colors[i])\n",
+ "\n",
+ " print(t, \"Dominant wavelength / D = \", (1/xf[np.argmax(yfMod)])/avgDiameter )\n",
+ "\n",
+ "plt.xlabel(\"wave length (LU)\")\n",
+ "plt.ylabel(\"amplitude (LU)\")\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "4c441c91",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# FFT along y per dx -> find period length (same as above with transposed array)\n",
+ "\n",
+ "import scipy.fft\n",
+ "\n",
+ "averageHb = np.average(hbOverTime) * dx\n",
+ "# use only data from some later time step to avoid including the influence of the initial condition\n",
+ "hbOverTimeShort = hbOverTime[:250,:] * dx - averageHb\n",
+ "\n",
+ "hbOverTimeTransposed = np.transpose(hbOverTimeShort)\n",
+ "\n",
+ "N = hbOverTimeTransposed.shape[-1] # = number of samples, here xSize\n",
+ "colors = sns.color_palette(\"crest\", hbOverTimeTransposed.shape[0] )\n",
+ "plt.figure()\n",
+ "\n",
+ "for i,t in enumerate(np.arange(0,hbOverTime.shape[-1])):\n",
+ "\n",
+ " yf = scipy.fft.rfft(hbOverTimeTransposed[i])\n",
+ " sampleSpacing = evalTimeSteps[1] - evalTimeSteps[0]\n",
+ " \n",
+ " # discard first element (is 0 frequency -> offset of sinus)\n",
+ " xf = scipy.fft.rfftfreq(N,sampleSpacing)[1:]\n",
+ " yfMod = (2.0 / N * np.abs(yf))[1:]\n",
+ " plt.plot(1/xf, yfMod, color=colors[i])\n",
+ "\n",
+ " print(t, \"Dominant period length / s = \", (1/xf[np.argmax(yfMod)]) * dt )\n",
+ "\n",
+ "plt.xlabel(\"Period length T / s\")\n",
+ "plt.ylabel(\"amplitude (LU)\")\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "19469336",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# # transport rate\n",
+ "\n",
+ "# def readValueOverTimeFromFile(fileName, columnIndex):\n",
+ "# return np.transpose(np.loadtxt(fileName, usecols=(0,columnIndex)))\n",
+ "\n",
+ "# sedimentTransportRateData = readValueOverTimeFromFile(baseFolder+'/bedload.txt',1)\n",
+ "\n",
+ "\n",
+ "# transportRateRef = avgDiameter * averageVelocity # = l_ref^2 / t_ref, m^2 / s\n",
+ "\n",
+ "# plt.figure()\n",
+ "# plt.plot(sedimentTransportRateData[0] / timeRef, sedimentTransportRateData[1] / transportRateRef, '-',label=\"simulation\")\n",
+ "\n",
+ "# plt.axhline(2.0e-5 / (2.9e-3 * 0.37), linestyle='--', color='k', label=\"E1 (experiment)\")\n",
+ "# plt.axhline(6.1e-5 / (2.9e-3 * 0.46), linestyle='-.', color='k', label=\"E4 (experiment)\")\n",
+ "# plt.xlabel(r\"$t / t_{ref}$\")\n",
+ "# plt.ylabel(r\"$q_s / q_{ref}$\")\n",
+ "# plt.grid()\n",
+ "# plt.legend()\n",
+ "# plt.show()\n",
+ "\n",
+ "# # import tikzplotlib\n",
+ "# # tikzplotlib.clean_figure()\n",
+ "# # tikzplotlib.save(\n",
+ "# # 'bedload-transport-rate.tex',\n",
+ "# # axis_height = '\\\\figureheight',\n",
+ "# # axis_width = '\\\\figurewidth'\n",
+ "# # )\n",
+ "\n",
+ "# # in Pascal Paper:\n",
+ "# # q_s,out approx q_s,in (Fig. 5)\n",
+ "# # -> take values from Tab. 1 as reference\n",
+ "# # E1: 2.0e-5 , E4: 6.1e-5 m^2/s\n",
+ "# # normalize by averageDiameter ( = 2.9mm) * averageVelocity (U)\n",
+ "# # -> E1: q_normalized = 0.0186, E4: q_normalized = 0.0457"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "b43b04f7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# transport rate in SI units\n",
+ "\n",
+ "def readValueOverTimeFromFile(fileName, columnIndex):\n",
+ " return np.transpose(np.loadtxt(fileName, usecols=(0,columnIndex)))\n",
+ "\n",
+ "sedimentTransportRateData = readValueOverTimeFromFile(baseFolder+'/bedload.txt',1)\n",
+ "\n",
+ "timestep_length = sedimentTransportRateData[0,1] - sedimentTransportRateData[0,0]\n",
+ "\n",
+ "t_min = 30\n",
+ "timestep_index_min = int((t_min / dt) // timestep_length)\n",
+ "\n",
+ "t_max = 75\n",
+ "timestep_index_max = int((t_max / dt) // timestep_length)\n",
+ "\n",
+ "plt.figure()\n",
+ "\n",
+ "## E1\n",
+ "#filepath = \"/simdata/on74yces/experiment-bedload-rate/qs_E1.csv\"\n",
+ "#exp_mean = np.array([[t_min, t_min + 10, t_max], [2.0e-5, 2.0e-5, 2.0e-5]])\n",
+ "\n",
+ "# E4\n",
+ "filepath = \"/simdata/on74yces/experiment-bedload-rate/qs_E4.csv\"\n",
+ "exp_mean = np.array([[t_min, t_min + 10, t_max], [6.1e-5, 6.1e-5, 6.1e-5]])\n",
+ "\n",
+ "data = np.transpose(np.loadtxt(filepath, delimiter=';',skiprows=1))\n",
+ "plt.plot(data[0,int(t_min//0.2):int(t_max//0.2)], data[1,int(t_min//0.2):int(t_max//0.2)], '-',label=\"Experiment\")\n",
+ "#plt.plot(exp_mean[0], exp_mean[1], linestyle='--', color='black', label=\"Experiment (mean)\")\n",
+ "\n",
+ "plt.plot(sedimentTransportRateData[0,timestep_index_min:timestep_index_max] * dt, \n",
+ " sedimentTransportRateData[1,timestep_index_min:timestep_index_max] * dx**2 / dt, '-',label=\"Simulation\")\n",
+ "mean = np.mean(sedimentTransportRateData[1,timestep_index_min:timestep_index_max] * dx**2 / dt)\n",
+ "mean_array = np.array([[t_min, t_min + 10, t_max], [mean, mean, mean]])\n",
+ "#plt.plot(mean_array[0], mean_array[1], linestyle='-', color='black', label=\"Simulation (mean)\")\n",
+ "\n",
+ "plt.xlabel(r\"$t$ / s\")\n",
+ "plt.ylabel(r\"$q_s$ / (m$^{2}$/s)\")\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "ax.set_xlim([t_min, t_max])\n",
+ "\n",
+ "plt.grid()\n",
+ "#plt.legend()\n",
+ "plt.show()\n",
+ "\n",
+ "print(np.mean(sedimentTransportRateData[1,timestep_index_min:timestep_index_max] * dx**2 / dt))\n",
+ "\n",
+ "# For tikzplotlib export:\n",
+ "# - remove legend\n",
+ "# - remove mean plots\n",
+ "import tikzplotlib\n",
+ "tikzplotlib.clean_figure()\n",
+ "tikzplotlib.save(\n",
+ " '/home/rzlin/ca36xymo/tikz/bedload-rate.tex',\n",
+ " axis_height = '\\\\figureheight',\n",
+ " axis_width = '\\\\figurewidth'\n",
+ " )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "40c316a1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# compute median particle diameter\n",
+ "\n",
+ "import numpy as np\n",
+ "\n",
+ "filepath = \"/simdata/on74yces/2395774/spheres_out.dat\"\n",
+ "\n",
+ "particles = np.loadtxt(filepath, delimiter=' ', skiprows=1)\n",
+ "\n",
+ "# 50th percentile diameter\n",
+ "diameterMedian = np.median(particles[:,4]*2)\n",
+ "print(\"50th percentile, i.e., median particle diameter: \", diameterMedian, \"m\")\n",
+ "print(\"50th percentile, i.e., median particle diameter: \", diameterMedian / dx, \"cells\")\n",
+ "\n",
+ "# 16th percentile diameter\n",
+ "diameter16th = np.percentile(particles[:,4]*2, 16)\n",
+ "print(\"16th percentile particle diameter: \", diameter16th, \"m\")\n",
+ "print(\"16th percentile particle diameter: \", diameter16th / dx, \"cells\")\n",
+ "\n",
+ "# 84th percentile diameter\n",
+ "diameter84th = np.percentile(particles[:,4]*2, 84)\n",
+ "print(\"84th percentile particle diameter: \", diameter84th, \"m\")\n",
+ "print(\"84th percentile particle diameter: \", diameter84th / dx, \"cells\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "5f7f03b1",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# plot data from experiment of Pascal et al.\n",
+ "import math\n",
+ "\n",
+ "filepath = \"/simdata/on74yces/experiment-bed-elevation/bed_surface_in_time_E4.csv\"\n",
+ "data = np.loadtxt(filepath, delimiter=';')\n",
+ "\n",
+ "# crop data to relevant time steps\n",
+ "data = data[1200:1450,:]\n",
+ "\n",
+ "# remove inclination from data\n",
+ "sine3degree = math.sin(3.0*math.pi/180)\n",
+ "for x in range(0, data.shape[1]):\n",
+ " data[:,x] = data[:,x] - abs(x-1280)*0.75/data.shape[1]*sine3degree\n",
+ "\n",
+ "xPlotLimitsN = np.array([0, data.shape[1]]) * 0.75/data.shape[1]\n",
+ "tPlotLimitsN = np.array([1200, 1200+data.shape[0]])\n",
+ "averageHbN = np.average(data)\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.imshow(data - averageHbN, interpolation='none', origin='lower',\n",
+ " aspect='auto', extent=(*xPlotLimitsN,*tPlotLimitsN), vmin=-5e-3, vmax=5e-3)\n",
+ "plt.colorbar()\n",
+ "plt.xlabel(r\"$x$ / m\")\n",
+ "plt.ylabel(r\"$t$ / s\")\n",
+ "plt.show()\n",
+ "\n",
+ "print(data.shape)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "45779e36",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import math\n",
+ "\n",
+ "# plot bed height evaluation of experiment\n",
+ " # time restricted to 1200 - 1245 s\n",
+ " # x-domain restricted to 0 - 0.75 m\n",
+ " \n",
+ "filepath = \"/simdata/on74yces/experiment-bed-elevation/bed_surface_in_time_E1.csv\"\n",
+ "hbOverTimeExp = np.loadtxt(filepath, delimiter=';')\n",
+ "\n",
+ "dxExp = 0.75 / 1280\n",
+ "\n",
+ "# remove inclination from data\n",
+ "sine3degree = math.sin(3*math.pi/180)\n",
+ "for x in range(0, hbOverTimeExp.shape[1]):\n",
+ " hbOverTimeExp[:,x] = hbOverTimeExp[:,x] - abs(x-1280)*0.75/hbOverTimeExp.shape[1]*sine3degree\n",
+ "\n",
+ "t_min = 1200\n",
+ "timestep_index_min = t_min\n",
+ "\n",
+ "t_max = 1245\n",
+ "timestep_index_max = t_max\n",
+ "\n",
+ "x_min = 0\n",
+ "x_index_min = int(x_min // dxExp)\n",
+ "\n",
+ "x_max = 0.75\n",
+ "x_index_max = int(x_max // dxExp)\n",
+ "\n",
+ "xPlotLimits = np.array([x_min,x_max])\n",
+ "tPlotLimits = np.array([t_min, t_max])\n",
+ "averageHbExp = np.average(hbOverTimeExp)\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.imshow(hbOverTimeExp[timestep_index_min:timestep_index_max,x_index_min:x_index_max] - averageHbExp, \n",
+ " interpolation='none', origin='lower', aspect='auto', extent=(*xPlotLimits,*tPlotLimits), \n",
+ " vmin=-5e-3, vmax=5e-3)\n",
+ "\n",
+ "cb = plt.colorbar(label=r\"$h$ / m\")\n",
+ "#match_colorbar(cb)\n",
+ "#cb.ax.tick_params(labelsize=8) \n",
+ "plt.xlabel(r\"$x$ / m\")\n",
+ "plt.ylabel(r\"$t$ / s\")\n",
+ "\n",
+ "import tikzplotlib\n",
+ "tikzplotlib.clean_figure()\n",
+ "tikzplotlib.save(\n",
+ " '/home/rzlin/ca36xymo/tikz/bed-elevation-experiment.tex',\n",
+ " axis_height = '\\\\figureheight',\n",
+ " axis_width = '\\\\figurewidth'\n",
+ " )\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "dec9d4b2",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import scipy.fft\n",
+ "\n",
+ "# plot spectral density of bed height in experiment\n",
+ "\n",
+ "filepath = \"/simdata/on74yces/experiment-bed-elevation/bed_surface_in_time_E1.csv\"\n",
+ "hbOverTimeExp = np.loadtxt(filepath, delimiter=';')\n",
+ "\n",
+ "dxExp = 0.75 / 1280\n",
+ "\n",
+ "# remove inclination from data\n",
+ "sine3degree = math.sin(2.9*math.pi/180)\n",
+ "for x in range(0, hbOverTimeExp.shape[1]):\n",
+ " hbOverTimeExp[:,x] = hbOverTimeExp[:,x] - abs(x-1280)*0.75/hbOverTimeExp.shape[1]*sine3degree\n",
+ "\n",
+ "averageHbExp = np.average(hbOverTimeExp)\n",
+ "\n",
+ "# power spectral density (PSD) on period-wavelength (PW) plane\n",
+ "psd_pw_exp = scipy.fft.fft2(hbOverTimeExp - averageHbExp)\n",
+ "psd_pw_exp = np.abs(psd_pw_exp)**2 / (hbOverTimeExp.shape[0] * hbOverTimeExp.shape[1])\n",
+ "\n",
+ "ft_period_exp = np.abs(1 / scipy.fft.fftfreq(np.transpose(hbOverTimeExp).shape[-1],1)[1:])\n",
+ "ft_wavelength_exp = np.abs(1 / scipy.fft.fftfreq(hbOverTimeExp.shape[-1],0.75/1280)[1:])\n",
+ "\n",
+ "plt.figure()\n",
+ "\n",
+ "plt.contourf(ft_period_exp, ft_wavelength_exp, np.transpose(np.abs(psd_pw_exp[1:,1:])))\n",
+ "#cb = plt.colorbar(label=r\"PSD / m$^2$\")\n",
+ "#match_colorbar(cb)\n",
+ "plt.xlabel(r\"$T$ / s\")\n",
+ "plt.ylabel(r\"$\\lambda$ / m\")\n",
+ "ax = plt.gca()\n",
+ "ax.set_xlim([5, 70])\n",
+ "ax.set_ylim([0.05, 0.2])\n",
+ "\n",
+ "# uncomment the following things to only store the png (no axes etc.)\n",
+ "# => tikz image must be created manually later\n",
+ "plt.subplots_adjust(bottom = 0)\n",
+ "plt.subplots_adjust(top = 1)\n",
+ "plt.subplots_adjust(right = 1)\n",
+ "plt.subplots_adjust(left = 0)\n",
+ "plt.axis('off')\n",
+ "plt.savefig('/home/rzlin/ca36xymo/tikz/spectral-density-experiment.png',bbox_inches='tight',transparent=True, pad_inches=0)\n",
+ "\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "cf86063d",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "psd_pw_exp_new = np.transpose(np.abs(psd_pw_exp[1:,1:]))\n",
+ "\n",
+ "# compute celerity (for each point in wavelength-period PSD-array)\n",
+ "# => identify all possible y-axis values for new celerity-wavelength PSD-array\n",
+ "c_exp = []\n",
+ "x = 0\n",
+ "for wavelength in ft_wavelength_exp:\n",
+ " y = 0\n",
+ " for period in ft_period_exp:\n",
+ " if psd_pw_exp_new[x,y] > 0.0001:\n",
+ " c_exp.append(wavelength / period)\n",
+ " y += 1\n",
+ " x += 1\n",
+ "\n",
+ "c_exp = list(set(c_exp)) # remove duplicate entries from list\n",
+ "c_exp.sort() # sort in ascending order\n",
+ "# list 'c' contains all possible celerities in ascending order without duplicates\n",
+ "# => this will be the x-axis of the celerity-wavelength PSD-array\n",
+ "\n",
+ "# create celerity-wavelength PSD-array \n",
+ "psd_cw = np.zeros((ft_wavelength_exp.shape[0], len(c_exp)))\n",
+ "\n",
+ "# fill this array with PSD at correct position\n",
+ "for i_wavelength in range(0, ft_wavelength_exp.shape[0]):\n",
+ " for i_period in range(0, ft_period_exp.shape[0]):\n",
+ " if psd_pw_exp_new[i_wavelength,i_period] > 0.0001:\n",
+ " celerity = ft_wavelength_exp[i_wavelength] / ft_period_exp[i_period] # compute celerity\n",
+ " i_celerity = c_exp.index(celerity) # map celerity to it's index on the x-axis of the new array\n",
+ " psd_cw[i_wavelength,i_celerity] = abs(psd_pw_exp_new[i_wavelength,i_period]) # store PSD at correct position in new array\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.contourf(ft_wavelength_exp, np.asarray(c_exp), np.abs(np.transpose(psd_cw)),colors='tab:blue')\n",
+ "#plt.colorbar()\n",
+ "plt.xlabel(r\"$\\lambda$ / m\")\n",
+ "plt.ylabel(r\"$c$ / m/s\")\n",
+ "ax = plt.gca()\n",
+ "ax.set_xlim([0.05, 0.2])\n",
+ "ax.set_ylim([0, 0.015])\n",
+ "\n",
+ "# uncomment the following things to only store the png (no axes etc.)\n",
+ "# => tikz image must be created manually later\n",
+ "plt.subplots_adjust(bottom = 0)\n",
+ "plt.subplots_adjust(top = 1)\n",
+ "plt.subplots_adjust(right = 1)\n",
+ "plt.subplots_adjust(left = 0)\n",
+ "plt.axis('off')\n",
+ "plt.savefig('/home/rzlin/ca36xymo/tikz/celerity-experiment-e1.png',bbox_inches='tight',transparent=True, pad_inches=0)\n",
+ "\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "7168b78e",
+ "metadata": {
+ "scrolled": true
+ },
+ "outputs": [],
+ "source": [
+ "# grid refinement study: transport rate in SI units\n",
+ "\n",
+ "def readValueOverTimeFromFile(fileName, columnIndex):\n",
+ " return np.transpose(np.loadtxt(fileName, usecols=(0,columnIndex)))\n",
+ "\n",
+ "t_min = 0\n",
+ "t_max = 50\n",
+ "t_min_for_mean = 10\n",
+ "\n",
+ "dx_coarse = 0.0005\n",
+ "dt_coarse = 2.16294e-05\n",
+ "sedimentTransportRateData_coarse = readValueOverTimeFromFile(\"/simdata/on74yces/resolution5\"+'/bedload.txt',1)\n",
+ "timestep_length_coarse = sedimentTransportRateData_coarse[0,1] - sedimentTransportRateData_coarse[0,0]\n",
+ "timestep_index_min_coarse = int((t_min / dt_coarse) // timestep_length_coarse)\n",
+ "timestep_index_max_coarse = int((t_max / dt_coarse) // timestep_length_coarse)\n",
+ "timestep_index_min_mean_coarse = int((35 / dt_coarse) // timestep_length_coarse)\n",
+ "mean_coarse = np.mean(sedimentTransportRateData_coarse[1,timestep_index_min_mean_coarse:timestep_index_max_coarse] * dx_coarse**2 / dt_coarse)\n",
+ "mean_array_coarse = np.array([[35, 35+1, t_max], [mean_coarse, mean_coarse, mean_coarse]])\n",
+ "\n",
+ "print(sedimentTransportRateData_coarse[0,-1])\n",
+ "print(sedimentTransportRateData_coarse[0,-1] * dt_coarse)\n",
+ "print(mean_coarse)\n",
+ "\n",
+ "dx_medium = 0.00025\n",
+ "dt_medium = 1.08147e-05\n",
+ "sedimentTransportRateData_medium = readValueOverTimeFromFile(\"/simdata/on74yces/resolution10\"+'/bedload.txt',1)\n",
+ "timestep_length_medium = sedimentTransportRateData_medium[0,1] - sedimentTransportRateData_medium[0,0]\n",
+ "timestep_index_min_medium = int((t_min / dt_medium) // timestep_length_medium)\n",
+ "timestep_index_max_medium = int((t_max / dt_medium) // timestep_length_medium)\n",
+ "timestep_index_min_mean_medium = int((t_min_for_mean / dt_medium) // timestep_length_medium)\n",
+ "mean_medium = np.mean(sedimentTransportRateData_medium[1,timestep_index_min_mean_medium:timestep_index_max_medium] * dx_medium**2 / dt_medium)\n",
+ "mean_array_medium = np.array([[t_min_for_mean, t_min_for_mean+1, t_max], [mean_medium, mean_medium, mean_medium]])\n",
+ "\n",
+ "print(sedimentTransportRateData_medium[0,-1])\n",
+ "print(sedimentTransportRateData_medium[0,-1] * dt_medium)\n",
+ "print(mean_medium)\n",
+ "\n",
+ "dx_fine = 0.000125\n",
+ "dt_fine = 5.40734e-06\n",
+ "sedimentTransportRateData_fine = readValueOverTimeFromFile(\"/simdata/on74yces/resolution20\"+'/bedload.txt',1)\n",
+ "timestep_length_fine = sedimentTransportRateData_fine[0,1] - sedimentTransportRateData_fine[0,0]\n",
+ "timestep_index_min_fine = int((t_min / dt_fine) // timestep_length_fine)\n",
+ "timestep_index_max_fine = int((t_max / dt_fine) // timestep_length_fine)\n",
+ "timestep_index_min_mean_fine = int((t_min_for_mean / dt_fine) // timestep_length_fine)\n",
+ "mean_fine = np.mean(sedimentTransportRateData_fine[1,timestep_index_min_mean_fine:timestep_index_max_fine] * dx_fine**2 / dt_fine)\n",
+ "mean_array_fine = np.array([[t_min_for_mean, t_min_for_mean+1, t_max], [mean_fine, mean_fine, mean_fine]])\n",
+ "\n",
+ "print(sedimentTransportRateData_fine[0,-1])\n",
+ "print(sedimentTransportRateData_fine[0,-1] * dt_fine)\n",
+ "print(mean_fine)\n",
+ "\n",
+ "plt.figure()\n",
+ "plt.plot(sedimentTransportRateData_coarse[0,timestep_index_min_coarse:timestep_index_max_coarse:4] * dt_coarse,\n",
+ " sedimentTransportRateData_coarse[1,timestep_index_min_coarse:timestep_index_max_coarse:4] * dx_coarse**2 / dt_coarse,\n",
+ " '-',label=\"Coarse\",color=\"tab:purple\")\n",
+ "\n",
+ "plt.plot(sedimentTransportRateData_medium[0,timestep_index_min_medium:timestep_index_max_medium:8] * dt_medium,\n",
+ " sedimentTransportRateData_medium[1,timestep_index_min_medium:timestep_index_max_medium:8] * dx_medium**2 / dt_medium,\n",
+ " '-',label=\"Medium\")\n",
+ "\n",
+ "plt.plot(sedimentTransportRateData_fine[0,timestep_index_min_fine:timestep_index_max_fine:16] * dt_fine,\n",
+ " sedimentTransportRateData_fine[1,timestep_index_min_fine:timestep_index_max_fine:16] * dx_fine**2 / dt_fine,\n",
+ " '-',label=\"Fine\")\n",
+ "\n",
+ "plt.plot(mean_array_coarse[0], mean_array_coarse[1], linestyle=':', color='black', label=\"Coarse (mean)\")\n",
+ "plt.plot(mean_array_medium[0], mean_array_medium[1], linestyle='-', color='black', label=\"Medium (mean)\")\n",
+ "plt.plot(mean_array_fine[0], mean_array_fine[1], linestyle='--', color='black', label=\"Fine (mean)\")\n",
+ "\n",
+ "plt.xlabel(r\"$t$ / s\")\n",
+ "plt.ylabel(r\"$q_s$ / (m$^{2}$/s)\")\n",
+ "\n",
+ "ax = plt.gca()\n",
+ "ax.set_xlim([t_min, t_max])\n",
+ "\n",
+ "plt.grid()\n",
+ "plt.legend()\n",
+ "plt.show()\n",
+ "\n",
+ "# # For tikzplotlib export:\n",
+ "# # - remove legend\n",
+ "# # - remove mean plots\n",
+ "# import tikzplotlib\n",
+ "# tikzplotlib.clean_figure()\n",
+ "# tikzplotlib.save(\n",
+ "# '/home/rzlin/ca36xymo/tikz/resolution-study.tex',\n",
+ "# axis_height = '\\\\figureheight',\n",
+ "# axis_width = '\\\\figurewidth'\n",
+ "# )"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "75ff7fe7",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def readValueOverTimeFromFile(fileName, columnIndex):\n",
+ " return np.loadtxt(fileName, usecols=(0,columnIndex))\n",
+ "\n",
+ "dx = 0.00025\n",
+ "dt = 1.08147e-05\n",
+ "velocity = readValueOverTimeFromFile(\"/simdata/on74yces/merge_E4_flat_new\"+'/fluidInfo.txt',2)\n",
+ "t = velocity[:,0] * dt\n",
+ "velocity = velocity[:,1] * dx/dt\n",
+ "plt.plot(t, velocity)\n",
+ "\n",
+ "plt.ylabel(r\"$U_{x,\\text{l}}$ / (m/s)\")\n",
+ "plt.xlabel(r\"$t$ / s\")\n",
+ "\n",
+ "import tikzplotlib\n",
+ "tikzplotlib.clean_figure()\n",
+ "tikzplotlib.save(\n",
+ " '/home/rzlin/ca36xymo/tikz/bed-elevation-simulation.tex',\n",
+ " axis_height = '\\\\figureheight',\n",
+ " axis_width = '\\\\figurewidth'\n",
+ " )\n",
+ "\n",
+ "plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "0605ba06",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# old stuff below"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "40b9bb68",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "(3*0.37 / 0.0083 * 1e-6)/((2.55-1)*9.81*2.9e-3)"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "1f88f760",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "print(np.max(hbOverTime[0,:]),np.min(hbOverTime[0,:]))"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "id": "345263f0",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import seaborn as sns\n",
+ "plt.figure()\n",
+ "sns.distplot(x=fluidSurfaceOverTime[0,:])\n",
+ "plt.show()"
+ ]
+ }
+ ],
+ "metadata": {
+ "kernelspec": {
+ "display_name": "Python 3 (ipykernel)",
+ "language": "python",
+ "name": "python3"
+ },
+ "language_info": {
+ "codemirror_mode": {
+ "name": "ipython",
+ "version": 3
+ },
+ "file_extension": ".py",
+ "mimetype": "text/x-python",
+ "name": "python",
+ "nbconvert_exporter": "python",
+ "pygments_lexer": "ipython3",
+ "version": "3.9.5"
+ }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
diff --git a/apps/showcases/CMakeLists.txt b/apps/showcases/CMakeLists.txt
index c6f0534cc6679f996777742823e5901aca8233f4..0b655239106e3c076b00e1c071e89d9cffb7d4d7 100644
--- a/apps/showcases/CMakeLists.txt
+++ b/apps/showcases/CMakeLists.txt
@@ -8,9 +8,12 @@ add_subdirectory( Mixer )
add_subdirectory( ParticlePacking )
add_subdirectory( PegIntoSphereBed )
if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )
-add_subdirectory( PhaseFieldAllenCahn )
+ add_subdirectory( PhaseFieldAllenCahn )
endif()
if ( WALBERLA_BUILD_WITH_CODEGEN AND NOT WALBERLA_BUILD_WITH_OPENMP)
add_subdirectory( PorousMedia )
endif()
+if ( WALBERLA_BUILD_WITH_CODEGEN )
+ add_subdirectory( Antidunes )
+endif()
diff --git a/src/lbm/blockforest/communication/SimpleCommunication.h b/src/lbm/blockforest/communication/SimpleCommunication.h
index 764152da65f081b164b102a41f3d639db5ae9065..7e256f92977f04a2c7e3d7d1321b409b3eb883fb 100644
--- a/src/lbm/blockforest/communication/SimpleCommunication.h
+++ b/src/lbm/blockforest/communication/SimpleCommunication.h
@@ -47,6 +47,7 @@ class SimpleCommunication : public communication::UniformBufferedScheme< Stencil
using VectorFieldFlattened_T = GhostLayerField< real_t, 3 >;
using PdfField_T = GhostLayerField< real_t, Stencil_T::Size >;
using UintScalarField_T = GhostLayerField< uint_t, 1 >;
+ using IDScalarField_T = walberla::GhostLayerField< walberla::id_t, 1 >;
using FlagField16_T = FlagField< uint16_t >;
using FlagField32_T = FlagField< uint32_t >;
@@ -148,17 +149,24 @@ class SimpleCommunication : public communication::UniformBufferedScheme< Stencil
}
else
{
- if (firstBlock.isDataClassOrSubclassOf< UintScalarField_T >(fieldId))
+ if (firstBlock.isDataClassOrSubclassOf< IDScalarField_T >(fieldId))
{
- this->addPackInfo(make_shared< PackInfo< UintScalarField_T > >(fieldId));
+ this->addPackInfo(make_shared< PackInfo< IDScalarField_T > >(fieldId));
}
else
{
- if (firstBlock.isDataClassOrSubclassOf< VectorFieldFlattened_T >(fieldId))
+ if (firstBlock.isDataClassOrSubclassOf< UintScalarField_T >(fieldId))
{
- this->addPackInfo(make_shared< PackInfo< VectorFieldFlattened_T > >(fieldId));
+ this->addPackInfo(make_shared< PackInfo< UintScalarField_T > >(fieldId));
+ }
+ else
+ {
+ if (firstBlock.isDataClassOrSubclassOf< VectorFieldFlattened_T >(fieldId))
+ {
+ this->addPackInfo(make_shared< PackInfo< VectorFieldFlattened_T > >(fieldId));
+ }
+ else { WALBERLA_ABORT("Problem with UID"); }
}
- else { WALBERLA_ABORT("Problem with UID"); }
}
}
}
diff --git a/src/lbm/free_surface/InitFunctions.h b/src/lbm/free_surface/InitFunctions.h
index 22cd09c62c794260668a70c30b07f52495c3a780..0dd2a0bbb25938f82ebc4a0957238e7852a086f6 100644
--- a/src/lbm/free_surface/InitFunctions.h
+++ b/src/lbm/free_surface/InitFunctions.h
@@ -85,7 +85,7 @@ void initFlagsFromFillLevels(const std::weak_ptr< StructuredBlockForest >& block
// set flags only in non-boundary and non-obstacle cells
if (!handling->isBoundary(fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z()))
{
- if (*fillFieldIt <= real_c(0))
+ if (floatIsEqual(*fillFieldIt, real_c(0), real_c(1e-14)))
{
// set gas flag
handling->forceFlag(flagInfo.gasFlag, fillFieldIt.x(), fillFieldIt.y(), fillFieldIt.z());
diff --git a/src/lbm/free_surface/InterfaceFromFillLevel.h b/src/lbm/free_surface/InterfaceFromFillLevel.h
index ef113e327032695d73df64109f0581c3b03cc4bd..f50c16cdd046f6bea042e163bdec5880f5172945 100644
--- a/src/lbm/free_surface/InterfaceFromFillLevel.h
+++ b/src/lbm/free_surface/InterfaceFromFillLevel.h
@@ -44,7 +44,7 @@ inline bool isInterfaceFromFillLevel(const ScalarField_T& fillField, cell_idx_t
real_t fillLevel = fillField.get(x, y, z);
// this cell is regular gas cell
- if (fillLevel <= real_c(0.0)) { return false; }
+ if (floatIsEqual(fillLevel, real_c(0.0), real_c(1e-14))) { return false; }
// this cell is regular interface cell
if (fillLevel < real_c(1.0)) { return true; }
diff --git a/src/lbm/free_surface/VtkWriter.h b/src/lbm/free_surface/VtkWriter.h
index 726e0ac735ef4495d234532a7e3eb3e52c0372a8..f8d3b453d79edc572e8d697ebc921b615a2f2486 100644
--- a/src/lbm/free_surface/VtkWriter.h
+++ b/src/lbm/free_surface/VtkWriter.h
@@ -87,6 +87,7 @@ void addVTKOutput(const std::weak_ptr< StructuredBlockForest >& blockForestPtr,
writers.push_back(std::make_shared< VTKWriter< VectorField_T, float > >(normalFieldID, "normal"));
writers.push_back(
std::make_shared< VTKWriter< VectorField_T, float > >(obstacleNormalFieldID, "obstacle_normal"));
+
if constexpr (useCodegen)
{
if (forceDensityFieldID != BlockDataID())
diff --git a/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h b/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h
index eeba79e77a31c0989ea397c859babb05c6eca8a0..215c03f3a543ea81545b5e291187bfce15e66734 100644
--- a/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h
+++ b/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h
@@ -80,7 +80,8 @@ public:
inline CurvedLinear( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField_T * const pdfField, const FlagField_T * const flagField,
ParticleField_T * const particleField, const shared_ptr<ParticleAccessor_T>& ac,
- const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block );
+ const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block,
+ std::function<real_t(const Vector3<real_t>&)> hydrostaticDensityFct = nullptr );
void pushFlags( std::vector< FlagUID >& uids ) const { uids.push_back( uid_ ); }
@@ -119,6 +120,8 @@ private:
const StructuredBlockStorage & blockStorage_;
const IBlock & block_;
+ std::function<real_t(const Vector3<real_t>&)> hydrostaticDensityFct_;
+
real_t lengthScalingFactor_;
real_t forceScalingFactor_;
@@ -128,8 +131,10 @@ private:
template< typename LatticeModel_T, typename FlagField_T, typename ParticleAccessor_T >
inline CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >::CurvedLinear( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField_T * const pdfField, const FlagField_T * const flagField,
ParticleField_T * const particleField, const shared_ptr<ParticleAccessor_T>& ac,
- const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block ):
-Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), flagField_( flagField ), particleField_( particleField ), ac_( ac ), domainMask_(domain), blockStorage_( blockStorage ), block_( block )
+ const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block,
+ std::function<real_t(const Vector3<real_t>&)> hydrostaticDensityFct ):
+Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), flagField_( flagField ), particleField_( particleField ),
+ac_( ac ), domainMask_(domain), blockStorage_( blockStorage ), block_( block ), hydrostaticDensityFct_(hydrostaticDensityFct)
{
WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ );
WALBERLA_ASSERT_NOT_NULLPTR( flagField_ );
@@ -159,12 +164,12 @@ inline void CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >::tre
WALBERLA_ASSERT_UNEQUAL( mask & this->mask_, numeric_cast<flag_t>(0) );
WALBERLA_ASSERT_EQUAL ( mask & this->mask_, this->mask_ ); // only true if "this->mask_" only contains one single flag, which is the case for the
// current implementation of this boundary condition
- WALBERLA_ASSERT_UNEQUAL( particleField_->get(nx,ny,nz), ac_->getInvalidUid() );
// determine distance to real boundary, i.e. delta value
// cell center of the near-boundary fluid cell
Cell nearBoundaryCell(x,y,z);
Vector3< real_t > cellCenter = blockStorage_.getBlockLocalCellCenter(block_, nearBoundaryCell);
+ WALBERLA_ASSERT_UNEQUAL( particleField_->get(nx,ny,nz), ac_->getInvalidUid(), x << " " << y << " " << z << " -> " << nx << " " << ny << " " << nz << " " << cellCenter);
// direction of the ray (from the fluid cell center to the boundary cell)
Vector3< real_t > direction( lengthScalingFactor_ * real_c( stencil::cx[ dir ] ),
@@ -173,7 +178,7 @@ inline void CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >::tre
//get particle index
auto particleIdx = ac_->uidToIdx(particleField_->get( nx, ny, nz ));
- WALBERLA_ASSERT_UNEQUAL( particleIdx, ac_->getInvalidIdx(), "Index of particle is invalid!" );
+ WALBERLA_ASSERT_UNEQUAL( particleIdx, ac_->getInvalidIdx(), "Index of particle is invalid! " << x << " " << y << " " << z << " -> " << nx << " " << ny << " " << nz << " " << cellCenter << " UID:" << particleField_->get( nx, ny, nz ) );
real_t delta = real_t(0);
real_t pdf_new = real_t(0);
@@ -268,8 +273,17 @@ inline void CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >::tre
// as a consequence, some (non-zero) PDF contributions would be missing after summing up the force contributions
// those would need to be added artificially, see e.g. Ernst, Dietzel, Sommerfeld - A lattice Boltzmann method for simulating transport and agglomeration of resolved particles, Acta Mech, 2013
// instead, we use the trick there that we just require the deviations from the equilibrium to get the correct force as it is already used for the incompressible case
- pdf_old -= LatticeModel_T::w[ Stencil_T::idx[dir] ];
- pdf_new -= LatticeModel_T::w[ Stencil_T::idx[dir] ];
+ if (hydrostaticDensityFct_ == nullptr)
+ {
+ pdf_old -= LatticeModel_T::w[Stencil_T::idx[dir]];
+ pdf_new -= LatticeModel_T::w[Stencil_T::idx[dir]];
+ }
+ else
+ {
+ const real_t rhoHydStat = hydrostaticDensityFct_(cellCenter);
+ pdf_old -= rhoHydStat * LatticeModel_T::w[Stencil_T::idx[dir]];
+ pdf_new -= rhoHydStat * LatticeModel_T::w[Stencil_T::idx[dir]];
+ }
}
// MEM: F = pdf_old + pdf_new - common
diff --git a/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h b/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h
index aeafb2b255fd965b248c71ee71c56b3c27972348..7d5c528489cae9578b34586238f282ba2dd0df7a 100644
--- a/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h
+++ b/src/lbm_mesapd_coupling/momentum_exchange_method/boundary/SimpleBB.h
@@ -72,7 +72,8 @@ public:
inline SimpleBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField_T * const pdfField, const FlagField_T * const flagField,
ParticleField_T * const particleField, const shared_ptr<ParticleAccessor_T>& ac,
- const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block );
+ const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block,
+ std::function<real_t(const Vector3<real_t>&)> hydrostaticDensityFct = nullptr );
void pushFlags( std::vector< FlagUID >& uids ) const { uids.push_back( uid_ ); }
@@ -110,6 +111,8 @@ private:
const StructuredBlockStorage & blockStorage_;
const IBlock & block_;
+ std::function<real_t(const Vector3<real_t>&)> hydrostaticDensityFct_;
+
real_t lengthScalingFactor_;
real_t forceScalingFactor_;
@@ -119,8 +122,10 @@ private:
template< typename LatticeModel_T, typename FlagField_T, typename ParticleAccessor_T >
inline SimpleBB< LatticeModel_T, FlagField_T, ParticleAccessor_T >::SimpleBB( const BoundaryUID & boundaryUID, const FlagUID & uid, PDFField_T * const pdfField, const FlagField_T * const flagField,
ParticleField_T * const particleField, const shared_ptr<ParticleAccessor_T>& ac,
- const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block ):
-Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), flagField_( flagField ), particleField_( particleField ), ac_( ac ), domainMask_(domain), blockStorage_( blockStorage ), block_( block )
+ const flag_t domain, const StructuredBlockStorage & blockStorage, const IBlock & block,
+ std::function<real_t(const Vector3<real_t>&)> hydrostaticDensityFct):
+Boundary<flag_t>( boundaryUID ), uid_( uid ), pdfField_( pdfField ), flagField_( flagField ), particleField_( particleField ), ac_( ac ), domainMask_(domain),
+blockStorage_( blockStorage ), block_( block ), hydrostaticDensityFct_(hydrostaticDensityFct)
{
WALBERLA_ASSERT_NOT_NULLPTR( pdfField_ );
WALBERLA_ASSERT_NOT_NULLPTR( flagField_ );
@@ -213,8 +218,17 @@ inline void SimpleBB< LatticeModel_T, FlagField_T, ParticleAccessor_T >::treatDi
// as a consequence, some (non-zero) PDF contributions would be missing after summing up the force contributions
// those would need to be added artificially, see e.g. Ernst, Dietzel, Sommerfeld - A lattice Boltzmann method for simulating transport and agglomeration of resolved particles, Acta Mech, 2013
// instead, we use the trick there that we just require the deviations from the equilibrium to get the correct force as it is already used for the incompressible case
- pdf_old -= LatticeModel_T::w[ Stencil_T::idx[dir] ];
- pdf_new -= LatticeModel_T::w[ Stencil_T::idx[dir] ];
+ if (hydrostaticDensityFct_ == nullptr)
+ {
+ pdf_old -= LatticeModel_T::w[Stencil_T::idx[dir]];
+ pdf_new -= LatticeModel_T::w[Stencil_T::idx[dir]];
+ }
+ else
+ {
+ const real_t rhoHydStat = hydrostaticDensityFct_(cellCenter);
+ pdf_old -= rhoHydStat * LatticeModel_T::w[Stencil_T::idx[dir]];
+ pdf_new -= rhoHydStat * LatticeModel_T::w[Stencil_T::idx[dir]];
+ }
}
// MEM: F = pdf_old + pdf_new - common