Merge branch 'lbm-mesapd-add-virtualmass-test' into 'master'

Add virtual mass functionality to lbm-mesapd coupling

See merge request walberla/walberla!472

apps/CMakeLists.txt

@@ -32,4 +32,4 @@ endif()
 # Python module
 if ( WALBERLA_BUILD_WITH_PYTHON )
     add_subdirectory( pythonmodule )
-endif()
\ No newline at end of file
+endif()

apps/showcases/CMakeLists.txt

@@ -2,6 +2,7 @@ add_subdirectory( BidisperseFluidizedBed )
 add_subdirectory( CombinedResolvedUnresolved )
 add_subdirectory( FluidizedBed )
 add_subdirectory( HeatConduction )
+add_subdirectory( LightRisingParticleInFluidAMR )
 add_subdirectory( Mixer )
 add_subdirectory( PegIntoSphereBed )
 if ( WALBERLA_BUILD_WITH_CODEGEN AND WALBERLA_BUILD_WITH_PYTHON )

apps/showcases/LightRisingParticleInFluidAMR/CMakeLists.txt

+waLBerla_add_executable ( NAME LightRisingParticleInFluidAMR
+        FILES LightRisingParticleInFluidAMR.cpp
+        DEPENDS core mesa_pd lbm lbm_mesapd_coupling domain_decomposition field vtk geometry postprocessing )

apps/showcases/LightRisingParticleInFluidAMR/LightRisingParticleInFluidAMR.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file LightRisingParticleInFluidAMR.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/SetupBlockForest.h"
+#include "blockforest/communication/UniformBufferedScheme.h"
+#include "blockforest/communication/UniformToNonUniformPackInfoAdapter.h"
+#include "blockforest/loadbalancing/DynamicCurve.h"
+#include "blockforest/loadbalancing/StaticCurve.h"
+
+#include "boundary/all.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/SharedFunctor.h"
+#include "core/debug/Debug.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/logging/all.h"
+#include "core/math/Random.h"
+#include "core/math/all.h"
+#include "core/mpi/Broadcast.h"
+#include "core/timing/RemainingTimeLogger.h"
+
+#include "domain_decomposition/BlockSweepWrapper.h"
+
+#include "field/AddToStorage.h"
+#include "field/StabilityChecker.h"
+#include "field/adaptors/AdaptorCreators.h"
+#include "field/communication/PackInfo.h"
+#include "field/vtk/all.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/field/QCriterionFieldWriter.h"
+#include "lbm/field/VelocityFieldWriter.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/refinement/all.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/vtk/all.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/all.h"
+
+#include <functional>
+
+#include "lbm_mesapd_coupling/DataTypes.h"
+#include "lbm_mesapd_coupling/amr/InfoCollection.h"
+#include "lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.h"
+#include "lbm_mesapd_coupling/mapping/ParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/MovingParticleMapping.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
+#include "lbm_mesapd_coupling/momentum_exchange_method/reconstruction/ExtrapolationDirectionFinder.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/AddForceOnParticlesKernel.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/ParticleSelector.h"
+#include "lbm_mesapd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/SubCyclingManager.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/AddVirtualForceTorqueKernel.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/InitializeVirtualMassKernel.h"
+#include "lbm_mesapd_coupling/utility/virtualmass/ParticleAccessorWithShapeVirtualMassWrapper.h"
+
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDataHandling.h"
+#include "mesa_pd/domain/BlockForestDomain.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/ReduceProperty.h"
+#include "mesa_pd/mpi/SyncGhostOwners.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/notifications/HydrodynamicForceTorqueNotification.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+/**
+ * \brief Showcase for a stabilized fluid - very light particle simulation using virtual mass
+ * A light sphere, situated at the bottom of a cuboid domain, rises.
+ * Depending on the density of the sphere (--sphereDensity) and galileo number (--galileoNumber), the sphere
+ * then undergoes various movement patterns along its way to the top.
+ * Employs adaptive mesh refinement and checkpointing.
+ * Can output VTK files based on the q criterion to visualize vortices behind the rising sphere.
+ * Most default parameters should be sufficient and do not need explicit values.
+ * As used in the master thesis of Lukas Werner, 2020, and follow-up paper at http://doi.org/10.1002/fld.5034
+ */
+
+namespace light_rising_particle_amr {
+using namespace walberla;
+using walberla::uint_t;
+
+#ifdef BUILD_WITH_MRT
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::D3Q19MRT>;
+#else
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT, false >;
+#endif
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+using ParticleField_T = lbm_mesapd_coupling::ParticleField_T;
+using VelocityField_T = GhostLayerField<Vector3<real_t>, 1>;
+using QCriterionField_T = GhostLayerField<real_t, 1>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+using ScalarField_T = GhostLayerField<real_t, 1>;
+
+const uint_t FieldGhostLayers = 4;
+
+using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+
+using NoSlip_T = lbm::NoSlip<LatticeModel_T, flag_t> ;
+using MO_T = lbm_mesapd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
+using BoundaryHandling_T = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T >;
+
+
+//// flags
+
+const FlagUID Fluid_Flag("fluid");
+const FlagUID NoSlip_Flag("no slip");
+const FlagUID MO_Flag("moving obstacle");
+const FlagUID FormerMO_Flag("former moving obstacle");
+
+
+//// block structure
+
+static void refinementSelection(SetupBlockForest &forest, uint_t levels, AABB refinementBox) {
+   real_t dx = real_t(1); // dx on finest level
+   const uint_t finestLevel = levels - uint_t(1);
+   for (auto block = forest.begin(); block != forest.end(); ++block) {
+      uint_t blockLevel = block->getLevel();
+      uint_t levelScalingFactor = (uint_t(1) << (finestLevel - blockLevel));
+      real_t dxOnLevel = dx * real_c(levelScalingFactor);
+      AABB blockAABB = block->getAABB();
+      // extend block AABB by ghostlayers
+      AABB extendedBlockAABB = blockAABB.getExtended(dxOnLevel * real_c(FieldGhostLayers));
+      if (extendedBlockAABB.intersects(refinementBox))
+         if (blockLevel < finestLevel) // only if the block is not on the finest level
+            block->setMarker(true);
+   }
+}
+
+static void workloadAndMemoryAssignment(SetupBlockForest &forest) {
+   for (auto block = forest.begin(); block != forest.end(); ++block) {
+      block->setWorkload(numeric_cast<workload_t>(uint_t(1) << block->getLevel()));
+      block->setMemory(numeric_cast<memory_t>(1));
+   }
+}
+static shared_ptr<StructuredBlockForest>
+createBlockStructure(const AABB &domainAABB, Vector3<uint_t> blockSizeInCells, uint_t numberOfLevels, real_t diameter,
+                     Vector3<real_t> spherePosition, Vector3<bool> periodicity,
+                     bool readFromCheckPointFile, const std::string & forestCheckPointFileName,
+                     bool keepGlobalBlockInformation = false) {
+   if (readFromCheckPointFile) {
+         WALBERLA_LOG_INFO_ON_ROOT("Reading block forest from file!");
+         return blockforest::createUniformBlockGrid(forestCheckPointFileName,
+                 blockSizeInCells[0], blockSizeInCells[1], blockSizeInCells[2],
+                 false);
+   } else {
+      SetupBlockForest sforest;
+      Vector3<uint_t> numberOfFineBlocksPerDirection(uint_c(domainAABB.size(0)) / blockSizeInCells[0],
+                                                     uint_c(domainAABB.size(1)) / blockSizeInCells[1],
+                                                     uint_c(domainAABB.size(2)) / blockSizeInCells[2]);
+      for (uint_t i = 0; i < 3; ++i) {
+         WALBERLA_CHECK_EQUAL(numberOfFineBlocksPerDirection[i] * blockSizeInCells[i], uint_c(domainAABB.size(i)),
+                              "Domain can not be decomposed in direction " << i << " into fine blocks of size "
+                                                                           << blockSizeInCells[i]);
+      }
+
+      uint_t levelScalingFactor = (uint_t(1) << (numberOfLevels - uint_t(1)));
+      Vector3<uint_t> numberOfCoarseBlocksPerDirection(numberOfFineBlocksPerDirection / levelScalingFactor);
+      for (uint_t i = 0; i < 3; ++i) {
+         WALBERLA_CHECK_EQUAL(numberOfCoarseBlocksPerDirection[i] * levelScalingFactor,
+                              numberOfFineBlocksPerDirection[i],
+                              "Domain can not be refined in direction " << i
+                                                                        << " according to the specified number of levels!");
+      }
+      AABB refinementBox( std::floor(spherePosition[0] - real_t(0.5) * diameter),
+                          std::floor(spherePosition[1] - real_t(0.5) * diameter),
+                          std::floor(spherePosition[2] - real_t(0.5) * diameter),
+                          std::ceil( spherePosition[0] + real_t(0.5) * diameter),
+                          std::ceil( spherePosition[1] + real_t(0.5) * diameter),
+                          std::ceil( spherePosition[2] + real_t(0.5) * diameter) );
+      WALBERLA_LOG_INFO_ON_ROOT(" - refinement box: " << refinementBox);
+
+      sforest.addRefinementSelectionFunction(
+              std::bind(refinementSelection, std::placeholders::_1, numberOfLevels, refinementBox));
+      sforest.addWorkloadMemorySUIDAssignmentFunction(workloadAndMemoryAssignment);
+
+      sforest.init(domainAABB, numberOfCoarseBlocksPerDirection[0], numberOfCoarseBlocksPerDirection[1],
+                   numberOfCoarseBlocksPerDirection[2], periodicity[0], periodicity[1], periodicity[2]);
+
+      // calculate process distribution
+      const memory_t memoryLimit = math::Limits<memory_t>::inf();
+      sforest.balanceLoad(blockforest::StaticLevelwiseCurveBalance(true),
+                          uint_c(MPIManager::instance()->numProcesses()), real_t(0), memoryLimit, true);
+
+      WALBERLA_LOG_INFO_ON_ROOT(sforest);
+
+      MPIManager::instance()->useWorldComm();
+
+      auto blockForest = make_shared<BlockForest>(uint_c(MPIManager::instance()->rank()), sforest,
+                                                  keepGlobalBlockInformation);
+
+      // create StructuredBlockForest (encapsulates a newly created BlockForest)
+      shared_ptr<StructuredBlockForest> sbf =
+              make_shared<StructuredBlockForest>(
+                      blockForest,
+                      blockSizeInCells[0], blockSizeInCells[1], blockSizeInCells[2]);
+      sbf->createCellBoundingBoxes();
+
+      return sbf;
+   }
+}
+
+//// boundary handling
+
+class MyBoundaryHandling : public blockforest::AlwaysInitializeBlockDataHandling< BoundaryHandling_T > {
+public:
+   MyBoundaryHandling(const weak_ptr< StructuredBlockStorage > & storage, const BlockDataID& flagFieldID,
+         const BlockDataID& pdfFieldID, const BlockDataID& particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) :
+         storage_(storage), flagFieldID_(flagFieldID), pdfFieldID_(pdfFieldID), particleFieldID_(particleFieldID),
+         ac_(ac) {
+
+   }
+
+   BoundaryHandling_T * initialize( IBlock * const block ) override;
+
+private:
+   weak_ptr< StructuredBlockStorage > storage_;
+
+   const BlockDataID flagFieldID_;
+   const BlockDataID pdfFieldID_;
+   const BlockDataID particleFieldID_;
+
+   shared_ptr<ParticleAccessor_T> ac_;
+};
+
+BoundaryHandling_T * MyBoundaryHandling::initialize( IBlock * const block )
+{
+   WALBERLA_ASSERT_NOT_NULLPTR( block );
+
+   auto* flagField = block->getData<FlagField_T>(flagFieldID_);
+   auto* pdfField = block->getData<PdfField_T>(pdfFieldID_);
+   auto* particleField = block->getData<lbm_mesapd_coupling::ParticleField_T>(particleFieldID_);
+
+   const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+   auto storagePtr = storage_.lock();
+   WALBERLA_CHECK_NOT_NULLPTR( storagePtr );
+
+   BoundaryHandling_T * handling = new BoundaryHandling_T( "moving obstacle boundary handling", flagField, fluid,
+                                                           NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                                           MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storagePtr, *block ) );
+
+   handling->fillWithDomain( FieldGhostLayers );
+
+   return handling;
+}
+
+void clearBoundaryHandling( BlockForest & forest, const BlockDataID & boundaryHandlingID ) {
+   for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
+   {
+      BoundaryHandling_T * boundaryHandling = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID);
+      boundaryHandling->clear( FieldGhostLayers );
+   }
+}
+
+void recreateBoundaryHandling( BlockForest & forest, const BlockDataID & boundaryHandlingID ) {
+   for( auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt )
+   {
+      BoundaryHandling_T * boundaryHandling = blockIt->getData<BoundaryHandling_T>(boundaryHandlingID);
+      boundaryHandling->fillWithDomain( FieldGhostLayers );
+   }
+}
+
+void clearParticleField( BlockForest & forest, const BlockDataID & particleFieldID, const ParticleAccessor_T & accessor ) {
+   for (auto blockIt = forest.begin(); blockIt != forest.end(); ++blockIt) {
+      ParticleField_T *particleField = blockIt->getData<ParticleField_T>(particleFieldID);
+      particleField->setWithGhostLayer(accessor.getInvalidUid());
+   }
+}
+
+class SpherePropertyLogger {
+public:
+   SpherePropertyLogger(lbm_mesapd_coupling::SubCyclingManager* mesapdTimeStep, const shared_ptr<ParticleAccessor_T>& ac, walberla::id_t sphereUid,
+                        const std::string& fileName, bool fileIO, bool writeHeader, real_t t_g, real_t u_g, real_t diameter) :
+         mesapdTimeStep_(mesapdTimeStep), ac_(ac), sphereUid_(sphereUid), fileName_(fileName), fileIO_(fileIO), t_g_(t_g),
+         u_g_(u_g), diameter_(diameter), position_(real_t(0)) {
+      if (fileIO_ && writeHeader) {
+         WALBERLA_ROOT_SECTION() {
+            std::ofstream file;
+            file.open(fileName_.c_str());
+            file << "#\t posZ\t uZ\t u\t rotX\t rotY\t rotZ\t hydFZ\t hydT\t tN\t posXN\t posYN\t posZN\t uXN\t uYN\t uZN\t uN\t wX\t wY\t wZ\n";
+            file.close();
+         }
+      }
+   }
+
+   void operator()() {
+      const uint_t timestep (mesapdTimeStep_->getCurrentTimeStep());
+      Vector3<real_t> transVel;
+      Vector3<real_t> hydForce;
+      Vector3<real_t> hydTorque;
+      Vector3<real_t> angVel;
+      Vector3<real_t> eulerRotation;
+
+      size_t idx = ac_->uidToIdx(sphereUid_);
+      if (idx != ac_->getInvalidIdx()) {
+         if (!mesa_pd::data::particle_flags::isSet(ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST)) {
+            transVel = ac_->getLinearVelocity(idx);
+            hydForce = ac_->getHydrodynamicForce(idx);
+            hydTorque = ac_->getHydrodynamicTorque(idx);
+            angVel = ac_->getAngularVelocity(idx);
+            eulerRotation = ac_->getRotation(idx).getMatrix().getEulerAnglesXYZ();
+         }
+      }
+
+      WALBERLA_MPI_SECTION() {
+         mpi::allReduceInplace(transVel[0], mpi::SUM);
+         mpi::allReduceInplace(transVel[1], mpi::SUM);
+         mpi::allReduceInplace(transVel[2], mpi::SUM);
+
+         mpi::allReduceInplace(angVel[0], mpi::SUM);
+         mpi::allReduceInplace(angVel[1], mpi::SUM);
+         mpi::allReduceInplace(angVel[2], mpi::SUM);
+
+         mpi::allReduceInplace(hydForce[0], mpi::SUM);
+         mpi::allReduceInplace(hydForce[1], mpi::SUM);
+         mpi::allReduceInplace(hydForce[2], mpi::SUM);
+
+         mpi::allReduceInplace(hydTorque[0], mpi::SUM);
+         mpi::allReduceInplace(hydTorque[1], mpi::SUM);
+         mpi::allReduceInplace(hydTorque[2], mpi::SUM);
+
+         mpi::allReduceInplace(eulerRotation[0], mpi::SUM);
+         mpi::allReduceInplace(eulerRotation[1], mpi::SUM);
+         mpi::allReduceInplace(eulerRotation[2], mpi::SUM);
+      }
+
+      syncPosition();
+
+      if (fileIO_) {
+         writeToFile(timestep, position_, eulerRotation, transVel, angVel, hydForce, hydTorque);
+      }
+   }
+
+   Vector3<real_t> getPosition() const {
+      return position_;
+   }
+
+   void syncPosition() {
+      Vector3<real_t> pos(real_t(0));
+
+      size_t idx = ac_->uidToIdx(sphereUid_);
+      if (idx != ac_->getInvalidIdx()) {
+         if (!mesa_pd::data::particle_flags::isSet(ac_->getFlags(idx), mesa_pd::data::particle_flags::GHOST)) {
+            pos = ac_->getPosition(idx);
+         }
+      }
+
+      WALBERLA_MPI_SECTION() {
+         mpi::allReduceInplace(pos[0], mpi::SUM);
+         mpi::allReduceInplace(pos[1], mpi::SUM);
+         mpi::allReduceInplace(pos[2], mpi::SUM);
+      }
+
+      position_ = pos;
+   }
+
+private:
+   void writeToFile(const uint_t timestep, const Vector3<real_t>& position, const Vector3<real_t>& eulerRotation,
+         const Vector3<real_t>& velocity, const Vector3<real_t>& angularVelocity,
+         const Vector3<real_t>& hydForce, const Vector3<real_t>& hydTorque) {
+      WALBERLA_ROOT_SECTION() {
+         std::ofstream file;
+         file.open(fileName_.c_str(), std::ofstream::app);
+
+         file << timestep << "\t"
+              << position[2] << "\t"
+              << velocity[2] << "\t"
+              << velocity.length() << "\t"
+              << eulerRotation[0] << "\t"
+              << eulerRotation[1] << "\t"
+              << eulerRotation[2] << "\t"
+              << hydForce[2] << "\t"
+              << hydTorque.length() << "\t"
+
+              << real_t(timestep) / t_g_ << "\t"
+              << position[0] / diameter_ << "\t"
+              << position[1] / diameter_ << "\t"
+              << position[2] / diameter_ << "\t"
+              << velocity[0] / u_g_ << "\t"
+              << velocity[1] / u_g_ << "\t"
+              << velocity[2] / u_g_ << "\t"
+              << velocity.length() / u_g_ << "\t"
+              << angularVelocity[0] * diameter_ / u_g_ << "\t"
+              << angularVelocity[1] * diameter_ / u_g_ << "\t"
+              << angularVelocity[2] * diameter_ / u_g_ << "\n";
+
+         file.close();
+      }
+   }
+
+   lbm_mesapd_coupling::SubCyclingManager* mesapdTimeStep_;
+   shared_ptr<ParticleAccessor_T> ac_;
+   const walberla::id_t sphereUid_;
+   std::string fileName_;
+   bool fileIO_;
+
+   real_t t_g_;
+   real_t u_g_;
+   real_t diameter_;
+
+   Vector3<real_t> position_;
+};
+
+void createCheckpointFiles(const std::string& checkPointFileName, std::vector<std::string>& oldCheckpointFiles,
+        uint_t lbmTimeStep, uint_t mesaTimeStep, const shared_ptr<BlockForest>& forest, const BlockDataID pdfFieldID,
+        const BlockDataID particleStorageID) {
+   WALBERLA_LOG_INFO_ON_ROOT("Writing checkpointing files in lbm time step " << lbmTimeStep << " (mesapd: " << mesaTimeStep << ").");
+
+   const std::string timeStepTag = std::to_string(lbmTimeStep) + "_" + std::to_string(mesaTimeStep);
+
+   const std::string lbmFileName = checkPointFileName + "_" + timeStepTag + "_lbm.txt";
+   const std::string mesaFileName = checkPointFileName + "_" + timeStepTag + "_mesa.txt";
+   const std::string forestFileName = checkPointFileName + "_" + timeStepTag + "_forest.txt";
+
+   forest->saveBlockData(lbmFileName, pdfFieldID);
+   forest->saveBlockData(mesaFileName, particleStorageID);
+   forest->saveToFile(forestFileName);
+
+   WALBERLA_ROOT_SECTION() {
+      if (!oldCheckpointFiles.empty()) {
+         for (const std::string& file: oldCheckpointFiles) {
+            filesystem::path path(file);
+            if (filesystem::exists(path)) {
+               filesystem::remove(path);
+            }
+         }
+         oldCheckpointFiles.clear();
+         WALBERLA_LOG_INFO("Deleted old checkpoint files.");
+      }
+
+      oldCheckpointFiles.push_back(lbmFileName);
+      oldCheckpointFiles.push_back(mesaFileName);
+      oldCheckpointFiles.push_back(forestFileName);
+   }
+}
+
+// main
+
+int main(int argc, char** argv) {
+   debug::enterTestMode();
+
+   mpi::Environment env(argc, argv);
+
+   logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
+   logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
+
+   math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
+   
+   // general params
+   bool fileIO = true;
+   uint_t vtkIOFreq = 0;
+   uint_t fullVtkIOFreq = 0;
+   uint_t qCriterionVtkIOFreq = 0;
+   std::string vtkIOSelection = "sliced"; // "none", "full", "thicksliced" (half the sphere still in), "trackingsliced" (sphere is followed in y direction)
+   std::string baseFolder = "vtk_out_UnsettlingSphereDR";
+   auto vtkLowerVelocityLimit = real_t(0.003);
+   auto vtkLowerQCriterionLimit = real_t(1e-10);
+   auto propertyLoggingFreq = uint_t(1);
+
+   // numerical params
+   auto blockSize = uint_t(32);
+   auto XBlocks = uint_t(32);
+   auto YBlocks = uint_t(32);
+   auto ZBlocks = uint_t(64);
+   auto refinementCheckFrequency = uint_t(0);
+   auto initialZ = real_t(-1);
+   auto diameter = real_t(16);
+   bool offsetInitialPosition = true;
+   bool averageForceTorqueOverTwoTimeSteps = true;
+   auto numMESAPDSubCycles = uint_t(10);
+   auto u_ref = real_t(0.01);
+   auto densitySphere = real_t(0.1);
+   auto timesteps = uint_t(2000);
+   auto galileoNumber = real_t(300);
+   auto Lambda_Bulk = real_t(100); // if = 1, normal TRT, else not. Up to 100.
+   bool conserveMomentum = false;
+
+   // virtual mass
+   bool useVirtualMass = true;
+   auto C_v = real_t(0.5);
+   auto C_v_omega = real_t(0.5);
+
+   // dynamic refinement
+   auto numberOfLevels = uint_t(3);
+
+   bool useVorticityCriterion = false;
+   auto lowerFluidRefinementLimit = real_t(0);
+   auto upperFluidRefinementLimit = std::numeric_limits<real_t>::infinity();
+
+   bool useCurlCriterion = true;
+   auto upperCurlLimit = real_t(0.001);
+   auto curlCoarsenFactor = real_t(0.2);
+   auto curlLengthScaleWeight = real_t(2);
+
+   // checkpointing
+   bool writeCheckPointFile = false;
+   bool readFromCheckPointFile = false;
+   uint_t checkPointingFreq = 0;