diff --git a/apps/CMakeLists.txt b/apps/CMakeLists.txt
index d834bd20ecc74452643579637bad66d3897e727c..60ca6fe6a1c8974d3708a5efe2ad432f4cb9ca40 100644
--- a/apps/CMakeLists.txt
+++ b/apps/CMakeLists.txt
@@ -36,3 +36,4 @@ if ( WALBERLA_BUILD_WITH_PYTHON )
     # waLBerla was build without -fPIC , so no linking into shared library is possible
 endif()
 
+add_subdirectory( virtualmass )
\ No newline at end of file
diff --git a/apps/virtualmass/CMakeLists.txt b/apps/virtualmass/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..617700c263b6a2ef6c183dcdfb5be7ce4566b2d3
--- /dev/null
+++ b/apps/virtualmass/CMakeLists.txt
@@ -0,0 +1,14 @@
+waLBerla_add_executable ( NAME MESA_PD_VIRTUALMASS
+                          FILES VirtualMass.cpp
+                          DEPENDS mesa_pd core )
+waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_VIRTUALMASS )
+
+waLBerla_add_executable ( NAME MESA_PD_APP_SETTLINGSPHERE_VM
+        FILES SettlingSphere.cpp
+        DEPENDS blockforest boundary core domain_decomposition field lbm pe pe_coupling postprocessing timeloop vtk )
+waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_SETTLINGSPHERE_VM )
+
+waLBerla_add_executable ( NAME MESA_PD_APP_FIXEDROTATINGSPHERE_VM
+        FILES FixedRotatingSphere.cpp
+        DEPENDS core mesa_pd lbm lbm_rpd_coupling domain_decomposition field vtk )
+waLBerla_execute_test( NO_MODULE_LABEL NAME MESA_PD_APP_FIXEDROTATINGSPHERE_VM )
diff --git a/apps/virtualmass/FixedRotatingSphere.cpp b/apps/virtualmass/FixedRotatingSphere.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..25662d3f9fba331477baca6bf2e8792036fee070
--- /dev/null
+++ b/apps/virtualmass/FixedRotatingSphere.cpp
@@ -0,0 +1,885 @@
+//======================================================================================================================
+//
+//  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 FixedRotatingSphere.cpp
+//! \ingroup lbm_rpd_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.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/math/all.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/logging/all.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/StabilityChecker.h"
+#include "field/communication/PackInfo.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/sweeps/SweepWrappers.h"
+
+#include "lbm_rpd_coupling/mapping/ParticleMapping.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/MovingParticleMapping.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/boundary/SimpleBB.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/reconstruction/ExtrapolationDirectionFinder.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
+#include "lbm_rpd_coupling/utility/AddForceOnParticlesKernel.h"
+#include "lbm_rpd_coupling/utility/ParticleSelectors.h"
+#include "lbm_rpd_coupling/DataTypes.h"
+#include "lbm_rpd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_rpd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_rpd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+#include "lbm_rpd_coupling/utility/LubricationCorrectionKernel.h"
+
+#include "mesa_pd/common/AnalyticContactDetection.h"
+#include "mesa_pd/data/ParticleAccessor.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/shape/HalfSpace.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/SpringDashpot.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include <mesa_pd/kernel/AddVirtualForceAndTorque.h>
+#include <mesa_pd/kernel/VirtualMass.h>
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/all.h"
+#include "field/vtk/all.h"
+#include "lbm/vtk/all.h"
+
+#include <functional>
+
+namespace settling_sphere
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT>;
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+const uint_t FieldGhostLayers = 1;
+
+///////////
+// 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" );
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+template <typename ParticleAccessor_T>
+class MyBoundaryHandling
+{
+public:
+
+   using NoSlip_T = lbm::NoSlip< LatticeModel_T, flag_t >;
+   using MO_T = lbm_rpd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
+   using Type = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T >;
+
+
+   MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
+                       const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) :
+         flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), particleFieldID_( particleFieldID ), ac_( ac ) {}
+
+   Type * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR( block );
+      WALBERLA_ASSERT_NOT_NULLPTR( storage );
+
+      auto * flagField     = block->getData< FlagField_T >( flagFieldID_ );
+      auto *  pdfField     = block->getData< PdfField_T > ( pdfFieldID_ );
+      auto * particleField = block->getData< lbm_rpd_coupling::ParticleField_T > ( particleFieldID_ );
+
+      const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+      Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid,
+                                  NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                  MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) );
+
+      handling->fillWithDomain( FieldGhostLayers );
+
+      return handling;
+   }
+
+private:
+
+   const BlockDataID flagFieldID_;
+   const BlockDataID pdfFieldID_;
+   const BlockDataID particleFieldID_;
+
+   shared_ptr<ParticleAccessor_T> ac_;
+};
+//*******************************************************************************************************************
+
+
+class ParticleAccessorWithShape : public mesa_pd::data::ParticleAccessor
+{
+public:
+   ParticleAccessorWithShape(const std::shared_ptr<mesa_pd::data::ParticleStorage>& ps, const std::shared_ptr<mesa_pd::data::ShapeStorage>& ss)
+         : ParticleAccessor(ps)
+         , ss_(ss)
+   {}
+
+   const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
+   walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
+   void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass() = v;}
+
+   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
+   auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
+   void setInvInertiaBF(const size_t p_idx, const mesa_pd::Mat3& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF() = v;}
+
+   mesa_pd::data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
+   walberla::real_t getVolume(const size_t p_idx) const { return getShape(p_idx)->getVolume(); }
+private:
+   std::shared_ptr<mesa_pd::data::ShapeStorage> ss_;
+};
+
+class ParticleAccessorWithShapeAndVirtualMass : public ParticleAccessorWithShape {
+public:
+    ParticleAccessorWithShapeAndVirtualMass(std::shared_ptr<mesa_pd::data::ParticleStorage>& ps, std::shared_ptr<mesa_pd::data::ShapeStorage>& ss)
+            : ParticleAccessorWithShape(ps, ss)
+    {}
+
+    /**
+     * Returns the inverse mass of the body by also calculating and including its virtual mass.
+     *
+     * \param p_idx index of particle
+     */
+    walberla::real_t getInvMass(const size_t p_idx) const {
+        return getInvMassIncludingVirtual(p_idx);
+    }
+    const auto getInvInertiaBF(const size_t p_idx) const {
+        return getInvInertiaBFIncludingVirtual(p_idx);
+    }
+
+    walberla::real_t& getInvMassRef(const size_t p_idx) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+        return ParticleAccessorWithShape::getInvMassRef(p_idx);
+    }
+    void setInvMass(const size_t p_idx, const walberla::real_t& v) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_UNUSED(v);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+    }
+    auto& getInvInertiaBFRef(const size_t p_idx) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+        return ParticleAccessorWithShape::getInvInertiaBFRef(p_idx);
+    }
+    void setInvInertiaBF(const size_t p_idx, const mesa_pd::Mat3& v) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_UNUSED(v);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+    }
+};
+
+//*******************************************************************************************************************
+/*!\brief Evaluating the position and velocity of the sphere
+ *
+ */
+//*******************************************************************************************************************
+template< typename ParticleAccessor_T>
+class SpherePropertyLogger
+{
+public:
+   SpherePropertyLogger( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid,
+         const std::string & fileName, bool fileIO, real_t diameter, real_t u_g, real_t t_g) :
+      ac_( ac ), sphereUid_( sphereUid ), fileName_( fileName ), fileIO_(fileIO), diameter_( diameter ), u_g_( u_g), t_g_(t_g),
+      position_( real_t(0) ), maxVelocity_( real_t(0) )
+   {
+      if ( fileIO_ )
+      {
+         WALBERLA_ROOT_SECTION()
+         {
+            std::ofstream file;
+            file.open( fileName_.c_str() );
+            //file << "#\t posX\t posY\t posZ\t velX\t velY\t velZ\t t_n\t velZ_n\n";
+            file << "#\t posZ\t velZ\t t_n\t velZ_n\n";
+            file.close();
+         }
+      }
+   }
+
+   void operator()(const uint_t timestep)
+   {
+      Vector3<real_t> pos(real_t(0));
+      Vector3<real_t> transVel(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);
+            transVel = ac_->getLinearVelocity(idx);
+         }
+      }
+
+      WALBERLA_MPI_SECTION()
+      {
+         //mpi::allReduceInplace( pos[0], mpi::SUM ); // not required
+         //mpi::allReduceInplace( pos[1], mpi::SUM ); // not required
+         mpi::allReduceInplace( pos[2], mpi::SUM );
+
+         //mpi::allReduceInplace( transVel[0], mpi::SUM ); // not required
+         //mpi::allReduceInplace( transVel[1], mpi::SUM ); // not required
+         mpi::reduceInplace( transVel[2], mpi::SUM );
+      }
+
+      position_ = pos[2];
+      maxVelocity_ = std::max(maxVelocity_, -transVel[2]);
+
+      if( fileIO_ )
+         writeToFile( timestep, pos, transVel);
+   }
+
+   real_t getPosition() const
+   {
+      return position_;
+   }
+
+   real_t getMaxVelocity() const
+   {
+      return maxVelocity_;
+   }
+
+private:
+   void writeToFile( const uint_t timestep, const Vector3<real_t> & position, const Vector3<real_t> & velocity )
+   {
+      WALBERLA_ROOT_SECTION()
+      {
+         std::ofstream file;
+         file.open( fileName_.c_str(), std::ofstream::app );
+
+         file << timestep
+              << "\t" << position[2]
+              << "\t" << velocity[2]
+              << "\t" << real_t(timestep)/t_g_ << "\t" << velocity[2]/u_g_
+              << "\n";
+         file.close();
+      }
+   }
+
+   shared_ptr< ParticleAccessor_T > ac_;
+   const walberla::id_t sphereUid_;
+   std::string fileName_;
+   bool fileIO_;
+   real_t diameter_;
+   real_t u_g_;
+   real_t t_g_;
+
+   real_t position_;
+   real_t maxVelocity_;
+};
+
+
+void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain)
+{
+   // create bounding planes
+   mesa_pd::data::Particle p0 = *ps->create(true);
+   p0.setPosition(simulationDomain.minCorner());
+   p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
+   p0.setOwner(mpi::MPIManager::instance()->rank());
+   p0.setType(0);
+   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::INFINITE_MASS);
+
+   mesa_pd::data::Particle p1 = *ps->create(true);
+   p1.setPosition(simulationDomain.maxCorner());
+   p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
+   p1.setOwner(mpi::MPIManager::instance()->rank());
+   p1.setType(0);
+   mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p2 = *ps->create(true);
+   p2.setPosition(simulationDomain.minCorner());
+   p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
+   p2.setOwner(mpi::MPIManager::instance()->rank());
+   p2.setType(0);
+   mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p3 = *ps->create(true);
+   p3.setPosition(simulationDomain.maxCorner());
+   p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
+   p3.setOwner(mpi::MPIManager::instance()->rank());
+   p3.setType(0);
+   mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p4 = *ps->create(true);
+   p4.setPosition(simulationDomain.minCorner());
+   p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
+   p4.setOwner(mpi::MPIManager::instance()->rank());
+   p4.setType(0);
+   mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p5 = *ps->create(true);
+   p5.setPosition(simulationDomain.maxCorner());
+   p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
+   p5.setOwner(mpi::MPIManager::instance()->rank());
+   p5.setType(0);
+   mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+}
+
+//////////
+// MAIN //
+//////////
+
+//*******************************************************************************************************************
+/*!\brief Testcase that simulates the settling of a sphere inside a rectangular column filled with viscous fluid
+ *
+ * see: ten Cate, Nieuwstad, Derksen, Van den Akker - "Particle imaging velocimetry experiments and lattice-Boltzmann
+ * simulations on a single sphere settling under gravity" (2002), Physics of Fluids, doi: 10.1063/1.1512918
+ */
+//*******************************************************************************************************************
+
+int main( int argc, char **argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   ///////////////////
+   // Customization //
+   ///////////////////
+
+   // simulation control
+   bool shortrun = false;
+   bool funcTest = false;
+   bool fileIO = false;
+   uint_t vtkIOFreq = 0;
+   std::string baseFolder = "vtk_out_SettlingSphere";
+
+   // physical setup
+   uint_t fluidType = 1;
+
+   //numerical parameters
+   uint_t numberOfCellsInHorizontalDirection = uint_t(135);
+   bool averageForceTorqueOverTwoTimSteps = true; //true
+   uint_t numRPDSubCycles = uint_t(1);
+   bool useVelocityVerlet = false;
+   bool useEANReconstructor = false;
+   bool useExtReconstructor = false;
+
+   bool useVirtualMass = false;
+   bool large = false;
+   real_t relaxationTime = 0.55;
+   real_t densitySphere = real_t(1.5); // 0.5 for rising, 1.5 for settling
+   auto timesteps = uint_t(2000);
+   auto diameterOverride = real_t(-1);
+
+   for( int i = 1; i < argc; ++i )
+   {
+      //if( std::strcmp( argv[i], "--shortrun" )            == 0 ) { shortrun = true; continue; }
+      //if( std::strcmp( argv[i], "--funcTest" )            == 0 ) { funcTest = true; continue; }
+      if( std::strcmp( argv[i], "--fileIO" )              == 0 ) { fileIO = true; continue; }
+      if( std::strcmp( argv[i], "--vtkIOFreq" )           == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--fluidType" )           == 0 ) { fluidType = uint_c( std::atof( argv[++i] ) ); continue; }
+      //if( std::strcmp( argv[i], "--numRPDSubCycles" )     == 0 ) { numRPDSubCycles = uint_c( std::atof( argv[++i] ) ); continue; }
+      //if( std::strcmp( argv[i], "--resolution" )          == 0 ) { numberOfCellsInHorizontalDirection = uint_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--noForceAveraging" )    == 0 ) { averageForceTorqueOverTwoTimSteps = false; continue; }
+      if( std::strcmp( argv[i], "--baseFolder" )          == 0 ) { baseFolder = argv[++i]; continue; }
+      if( std::strcmp( argv[i], "--useVV" )               == 0 ) { useVelocityVerlet = true; continue; }
+      if( std::strcmp( argv[i], "--useEANReconstructor" ) == 0 ) { useEANReconstructor = true; continue; }
+      if( std::strcmp( argv[i], "--useExtReconstructor" ) == 0 ) { useExtReconstructor = true; continue; }
+      if( std::strcmp( argv[i], "--rising" )              == 0 ) { densitySphere = real_t(0.5); continue; }
+      if( std::strcmp( argv[i], "--sphereDensity" )       == 0 ) { densitySphere = real_c( std::stod( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--relaxationTime")       == 0 ) { relaxationTime = real_c( std::stod( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--virtualMass" )         == 0 ) { useVirtualMass = true; continue; }
+      if( std::strcmp( argv[i], "--large" )               == 0 ) { large = true; continue; }
+      if( std::strcmp( argv[i], "--timesteps" )           == 0 ) { timesteps = uint_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--diameter" )       == 0 ) { diameterOverride = real_c( std::stod( argv[++i] ) ); continue; }
+      WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+   }
+
+   WALBERLA_CHECK(!(useEANReconstructor && useExtReconstructor), "Can not use two reconstructors!");
+
+   if( funcTest )
+   {
+      walberla::logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::WARNING);
+   }
+
+   if( fileIO )
+   {
+      // create base directory if it does not yet exist
+      filesystem::path tpath( baseFolder );
+      if( !filesystem::exists( tpath ) )
+         filesystem::create_directory( tpath );
+   }
+
+   //////////////////////////
+   // NUMERICAL PARAMETERS //
+   //////////////////////////
+
+
+   const auto densityFluid = real_t(1);
+   const auto galileiNumber = real_t(170);
+   // relaxation time via arguments
+
+   const auto diameter = (diameterOverride > 0) ? diameterOverride : real_t(20));
+   const auto dt = real_t(5e-4);
+
+   const auto Nx = uint_t(real_t(30)*diameter);
+   const auto Ny = uint_t(real_t(15)*diameter);
+   const auto Nz = uint_t(real_t(15)*diameter);
+
+   const auto omega = real_t(1)/relaxationTime;
+   const auto viscosity = lbm::collision_model::viscosityFromOmega(omega);
+
+   const auto d3 = diameter*diameter*diameter;
+   const auto G2 = galileiNumber*galileiNumber;
+   const auto v2 = viscosity*viscosity;
+   const auto gravitationalAcceleration = G2*v2 / (abs(densitySphere-real_t(1))*d3);
+
+   const Vector3<uint_t> domainSize( Nx, Ny, Nz);
+   Vector3<real_t> initialPosition = Vector3<real_t>(real_t(10) * diameter,
+                                                     real_t(7.5) * diameter,
+                                                     real_t(7.5) * diameter);
+
+   const real_t sphereVolume = math::M_PI / real_t(6) * diameter * diameter * diameter;
+
+   const auto dx = real_t(1);
+
+   // virtual mass parameters
+   const auto C_v = real_t(0.5);
+   const auto C_v_omega = C_v;
+
+   const auto u_g = std::sqrt(std::fabs(densitySphere - real_t(1)) * gravitationalAcceleration * diameter);
+   const auto t_g = std::sqrt(diameter / (std::fabs(densitySphere - real_t(1)) * gravitationalAcceleration));
+
+   const bool rising = densitySphere < real_t(1.0);
+
+   //WALBERLA_LOG_INFO_ON_ROOT(" - dx_SI = " << dx_SI << ", dt_SI = " << dt_SI);
+   WALBERLA_LOG_INFO_ON_ROOT("Setup (in simulation, i.e. lattice, units):");
+   WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize << (large ? " (large simulation)" : ""));
+   WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter << ", density = " << densitySphere );
+   WALBERLA_LOG_INFO_ON_ROOT(" - fluid: density = " << densityFluid << ", relaxation time (tau) = " << relaxationTime << ", kin. visc = " << viscosity );
+   WALBERLA_LOG_INFO_ON_ROOT(" - gravitational acceleration = " << gravitationalAcceleration );
+   //WALBERLA_LOG_INFO_ON_ROOT(" - expected settling velocity = " << expectedSettlingVelocity << " --> Re_p = " << expectedSettlingVelocity * diameter / viscosity );
+   WALBERLA_LOG_INFO_ON_ROOT(" - integrator = " << (useVelocityVerlet ? "Velocity Verlet" : "Explicit Euler") );
+   WALBERLA_LOG_INFO_ON_ROOT(" - reconstructor = " << (useEANReconstructor ? "EAN" : "") << (useExtReconstructor ? "Ext" : "") << (!useEANReconstructor && !useExtReconstructor ? "default" : ""))
+   WALBERLA_LOG_INFO_ON_ROOT(" - timesteps = " << timesteps);
+
+   if( vtkIOFreq > 0 )
+   {
+      WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files to folder \"" << baseFolder << "\" with frequency " << vtkIOFreq);
+   }
+
+   ///////////////////////////
+   // BLOCK STRUCTURE SETUP //
+   ///////////////////////////
+
+   Vector3<uint_t> numberOfBlocksPerDirection( uint_t(1), uint_t(1), uint_t(4) );
+   Vector3<uint_t> cellsPerBlockPerDirection( domainSize[0] / numberOfBlocksPerDirection[0],
+                                              domainSize[1] / numberOfBlocksPerDirection[1],
+                                              domainSize[2] / numberOfBlocksPerDirection[2] );
+   for( uint_t i = 0; i < 3; ++i ) {
+      WALBERLA_CHECK_EQUAL(cellsPerBlockPerDirection[i] * numberOfBlocksPerDirection[i], domainSize[i],
+                           "Unmatching domain decomposition in direction " << i << "!");
+   }
+
+   auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1], numberOfBlocksPerDirection[2],
+                                                      cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], dx,
+                                                      0, false, false,
+                                                      true, true, true, //periodicity
+                                                      false );
+
+   WALBERLA_LOG_INFO_ON_ROOT("Domain decomposition:");
+   WALBERLA_LOG_INFO_ON_ROOT(" - blocks per direction = " << numberOfBlocksPerDirection );
+   WALBERLA_LOG_INFO_ON_ROOT(" - cells per block = " << cellsPerBlockPerDirection );
+
+   //write domain decomposition to file
+   if( vtkIOFreq > 0 )
+   {
+      vtk::writeDomainDecomposition( blocks, "initial_domain_decomposition", baseFolder );
+   }
+
+   //////////////////
+   // RPD COUPLING //
+   //////////////////
+
+   auto rpdDomain = std::make_shared<mesa_pd::domain::BlockForestDomain>(blocks->getBlockForestPointer());
+
+   //init data structures
+   auto ps = walberla::make_shared<mesa_pd::data::ParticleStorage>(1);
+   auto ss = walberla::make_shared<mesa_pd::data::ShapeStorage>();
+   using ParticleAccessor_T = ParticleAccessorWithShape;
+   auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
+   auto virtualMassAccessor = walberla::make_shared<ParticleAccessorWithShapeAndVirtualMass>(ps, ss);
+
+   // bounding planes
+   // excluded for virtual mass test case (periodic boundary conditions)
+   //createPlaneSetup(ps,ss,blocks->getDomain());
+
+   // create sphere and store Uid
+   auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) );
+   ss->shapes[sphereShape]->updateMassAndInertia(densitySphere);
+
+   walberla::id_t sphereUid = 0;
+   if (rpdDomain->isContainedInProcessSubdomain( uint_c(mpi::MPIManager::instance()->rank()), initialPosition ))
+   {
+      mesa_pd::data::Particle&& p = *ps->create();
+      p.setPosition(initialPosition);
+      p.setInteractionRadius(diameter * real_t(0.5));
+      p.setOwner(mpi::MPIManager::instance()->rank());
+      p.setShapeID(sphereShape);
+      sphereUid = p.getUid();
+   }
+   mpi::allReduceInplace(sphereUid, mpi::SUM);
+
+   ///////////////////////
+   // ADD DATA TO BLOCKS //
+   ////////////////////////
+
+   // create the lattice model
+   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( real_t(1) / relaxationTime ) );
+
+   // add PDF field
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
+                                                                         Vector3< real_t >( real_t(0) ), real_t(1),
+                                                                         uint_t(1), field::zyxf );
+   // add flag field
+   BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" );
+
+   // add particle field
+   BlockDataID particleFieldID = field::addToStorage<lbm_rpd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::zyxf, FieldGhostLayers );
+
+   // add boundary handling
+   using BoundaryHandling_T = MyBoundaryHandling<ParticleAccessor_T>::Type;
+   BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(MyBoundaryHandling<ParticleAccessor_T>( flagFieldID, pdfFieldID, particleFieldID, accessor), "boundary handling" );
+
+   // set up RPD functionality
+   std::function<void(void)> syncCall = [ps,rpdDomain](){
+      const real_t overlap = real_t( 1.5 );
+      mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+      syncNextNeighborFunc(*ps, *rpdDomain, overlap);
+   };
+
+   syncCall();
+
+   mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
+   mesa_pd::kernel::VelocityVerletPreForceUpdate  vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
+   mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
+
+   mesa_pd::kernel::SpringDashpot collisionResponse(1);
+   mesa_pd::mpi::ReduceProperty reduceProperty;
+
+   // set up coupling functionality
+   Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume );
+   lbm_rpd_coupling::AddForceOnParticlesKernel<lbm_rpd_coupling::LocalParticlesSelector> addGravitationalForce(gravitationalForce);
+
+   lbm_rpd_coupling::RegularParticlesSelector sphereSelector;
+   lbm_rpd_coupling::AddHydrodynamicInteractionKernel<decltype(sphereSelector)> addHydrodynamicInteraction(sphereSelector);
+
+   lbm_rpd_coupling::ResetHydrodynamicForceTorqueKernel<> resetHydrodynamicForceTorque;
+
+   lbm_rpd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
+
+   // virtual mass approach
+   mesa_pd::kernel::VirtualMass virtualMass(C_v, C_v_omega, densityFluid);
+   mesa_pd::kernel::AddVirtualForceAndTorque addVirtualForceAndTorque;
+
+   if ( useVirtualMass ) ps->forEachParticle(false, virtualMass, *accessor);
+
+   ///////////////
+   // TIME LOOP //
+   ///////////////
+
+   // map planes into the LBM simulation -> act as no-slip boundaries
+   lbm_rpd_coupling::mapParticles<BoundaryHandling_T>(*blocks, boundaryHandlingID, *accessor, NoSlip_Flag, lbm_rpd_coupling::GlobalParticlesSelector());
+
+   // map particles into the LBM simulation
+   lbm_rpd_coupling::mapMovingParticles<BoundaryHandling_T>(*blocks, boundaryHandlingID, particleFieldID, *accessor, MO_Flag, sphereSelector);
+
+
+   // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+   std::function< void () > commFunction;
+   blockforest::communication::UniformBufferedScheme< Stencil_T > scheme( blocks );
+   scheme.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) );
+   commFunction = scheme;
+
+   // create the timeloop
+   SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
+
+   auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
+   auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+
+   timeloop.addFuncBeforeTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+   // vtk output
+   if( vtkIOFreq != uint_t(0) )
+   {
+      // spheres
+      auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+      particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner");
+      particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
+      auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
+      timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
+
+      // flag field (written only once in the first time step, ghost layers are also written)
+      //auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", timesteps, FieldGhostLayers, false, baseFolder );
+      //flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
+      //timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
+
+      // pdf field
+      auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", vtkIOFreq, 0, false, baseFolder );
+
+      blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > pdfGhostLayerSync( blocks );
+      pdfGhostLayerSync.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) );
+      pdfFieldVTK->addBeforeFunction( pdfGhostLayerSync );
+
+      field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldID );
+      fluidFilter.addFlag( Fluid_Flag );
+      pdfFieldVTK->addCellInclusionFilter( fluidFilter );
+
+      pdfFieldVTK->addCellDataWriter( make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldID, "VelocityFromPDF" ) );
+      pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
+
+      timeloop.addFuncBeforeTimeStep( vtk::writeFiles( pdfFieldVTK ), "VTK (fluid field data)" );
+   }
+
+   // sweep for updating the particle mapping into the LBM simulation
+   timeloop.add() << Sweep( lbm_rpd_coupling::makeMovingParticleMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, particleFieldID, accessor, MO_Flag, FormerMO_Flag, sphereSelector), "Particle Mapping" );
+
+   // sweep for restoring PDFs in cells previously occupied by particles
+   if( useEANReconstructor )
+   {
+      auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_rpd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
+      auto equilibriumAndNonEquilibriumSphereNormalReconstructor = lbm_rpd_coupling::makeEquilibriumAndNonEquilibriumReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder);
+      auto reconstructionManager = lbm_rpd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, equilibriumAndNonEquilibriumSphereNormalReconstructor);
+
+      timeloop.add() << Sweep( makeSharedSweep(reconstructionManager), "PDF Restore" );
+   } else if( useExtReconstructor ) {
+      auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_rpd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
+      auto extrapolationSphereNormalReconstructor = lbm_rpd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, true);
+      timeloop.add() << Sweep( makeSharedSweep(lbm_rpd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, extrapolationSphereNormalReconstructor) ), "PDF Restore" );
+   } else {
+      timeloop.add() << Sweep( makeSharedSweep(lbm_rpd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag) ), "PDF Restore" );
+   }
+
+   bool useStreamCollide = true;
+
+   if( useStreamCollide )
+   {
+      // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
+      timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
+                     << Sweep(bhSweep, "Boundary Handling" );
+
+      // stream + collide LBM step
+      timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+   }
+   else
+   {
+      // collide LBM step
+      timeloop.add() << Sweep( makeCollideSweep( lbmSweep ), "cell-wise collide LB sweep" );
+
+      // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
+      timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
+                     << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+
+      // stream LBM step
+      timeloop.add() << Sweep( makeStreamSweep( lbmSweep ), "cell-wise stream LB sweep" );
+   }
+
+   // evaluation functionality
+   std::string loggingFileName( baseFolder + "/LoggingSettlingSphere_");
+   loggingFileName += std::to_string(fluidType);
+   loggingFileName += ".txt";
+   if( fileIO  )
+   {
+      WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\"");
+   }
+   SpherePropertyLogger<ParticleAccessor_T> logger( accessor, sphereUid, loggingFileName, fileIO, diameter, u_g, t_g);
+
+
+   ////////////////////////
+   // EXECUTE SIMULATION //
+   ////////////////////////
+
+   WcTimingPool timeloopTiming;
+
+   real_t terminationPosition;
+   if (large) {
+       // when the sphere has covered a distance of L_z - diameter
+       if (rising) {
+           terminationPosition = initialPosition[2] + real_t(domainSize[2]) - diameter;
+       } else {
+           terminationPosition = initialPosition[2] + diameter;
+       }
+   } else {
+       // when the sphere is only a diameter away from the initial position again
+       if (rising) {
+           terminationPosition = real_t(initialPosition[2]) - diameter;
+       } else {
+           terminationPosition = real_t(initialPosition[2]) + diameter;
+       }
+   }
+   WALBERLA_LOG_INFO_ON_ROOT(" - termination position = " << terminationPosition );
+
+   //  required for determining termination
+   bool hasCrossedBoundary = false;
+   bool lockSphereLocation = true;
+   bool lockSphereRotation = true;
+
+   const bool useOpenMP = false;
+
+   // time loop
+   for (uint_t i = 0; i < timesteps; ++i )
+   {
+      // perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
+      timeloop.singleStep( timeloopTiming );
+
+      timeloopTiming["RPD"].start();
+
+      if( averageForceTorqueOverTwoTimSteps && i!= 0)
+      {
+         ps->forEachParticle(useOpenMP, averageHydrodynamicForceTorque, *accessor );
+      }
+
+      for(auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle )
+      {
+
+         if( useVelocityVerlet ) {
+            if ( useVirtualMass ) ps->forEachParticle(useOpenMP, vvIntegratorPreForce, *virtualMassAccessor);
+            else ps->forEachParticle(useOpenMP, vvIntegratorPreForce, *accessor);
+         }
+
+         ps->forEachParticle(useOpenMP, addHydrodynamicInteraction, *accessor );
+         ps->forEachParticle(useOpenMP, addGravitationalForce, *accessor );
+
+         reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
+
+         auto& integrationAccessor = useVirtualMass ? *virtualMassAccessor : *accessor;
+
+         if ( lockSphereLocation ) {
+            ps->forEachParticle(useOpenMP, [](size_t idx, auto &ac) {
+                ac.setForce(idx, Vector3<real_t>(0));
+            }, *accessor);
+         }
+         if ( lockSphereRotation ) {
+            ps->forEachParticle(useOpenMP, [](size_t idx, auto &ac) {
+                ac.setTorque(idx, Vector3<real_t>(0));
+            }, *accessor);
+         }
+
+         if ( useVirtualMass ) { //TODO: check ordering of virtual mass application and sphere fixation
+            ps->forEachParticle(useOpenMP, addVirtualForceAndTorque, *accessor);
+         }
+
+         if (useVelocityVerlet) {
+            ps->forEachParticle(useOpenMP, vvIntegratorPostForce, integrationAccessor);
+         } else {
+            ps->forEachParticle(useOpenMP, explEulerIntegrator, integrationAccessor);
+         }
+
+         syncCall();
+      }
+
+      timeloopTiming["RPD"].end();
+
+      // evaluation
+      timeloopTiming["Logging"].start();
+      logger(i);
+      timeloopTiming["Logging"].end();
+
+      // reset after logging here
+      ps->forEachParticle(useOpenMP, resetHydrodynamicForceTorque, *accessor );
+
+      // check for termination
+      // check if the sphere crossed a boundary
+      if (!hasCrossedBoundary && ((rising && logger.getPosition() < terminationPosition) || (!rising && logger.getPosition() > terminationPosition))) {
+          hasCrossedBoundary = true;
+      }
+      if ((!large && hasCrossedBoundary && ((rising && logger.getPosition() > terminationPosition) || (!rising && logger.getPosition() < terminationPosition)))
+        || (large && ((rising && logger.getPosition() > terminationPosition) || (!rising && hasCrossedBoundary && logger.getPosition() < terminationPosition))))
+      {
+          // end of simulation reached after crossing a boundary and then passing the termination position
+          WALBERLA_LOG_INFO_ON_ROOT("Sphere reached terminal position " << logger.getPosition() << " after " << i << " timesteps!");
+          break;
+      }
+   }
+
+   timeloopTiming.logResultOnRoot();
+
+   // check the result
+   /*if ( !funcTest && !shortrun )
+   {
+      real_t relErr = std::fabs( expectedSettlingVelocity - logger.getMaxVelocity()) / expectedSettlingVelocity;
+      WALBERLA_LOG_INFO_ON_ROOT( "Expected maximum settling velocity: " << expectedSettlingVelocity );
+      WALBERLA_LOG_INFO_ON_ROOT( "Simulated maximum settling velocity: " << logger.getMaxVelocity() );
+      WALBERLA_LOG_INFO_ON_ROOT( "Relative error: " << relErr );
+
+      // the relative error has to be below 10%
+      WALBERLA_CHECK_LESS( relErr, real_t(0.1) );
+   }*/
+
+   return EXIT_SUCCESS;
+}
+
+} // namespace settling_sphere_mem
+
+int main( int argc, char **argv ){
+   settling_sphere::main(argc, argv);
+}
diff --git a/apps/virtualmass/SettlingSphere.cpp b/apps/virtualmass/SettlingSphere.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..322e6687625692406623957ec0f30a300cfe8025
--- /dev/null
+++ b/apps/virtualmass/SettlingSphere.cpp
@@ -0,0 +1,973 @@
+//======================================================================================================================
+//
+//  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 SettlingSphere.cpp
+//! \ingroup lbm_rpd_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+#include "blockforest/communication/UniformBufferedScheme.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/math/all.h"
+#include "core/timing/RemainingTimeLogger.h"
+#include "core/logging/all.h"
+
+#include "domain_decomposition/SharedSweep.h"
+
+#include "field/AddToStorage.h"
+#include "field/StabilityChecker.h"
+#include "field/communication/PackInfo.h"
+
+#include "lbm/boundary/all.h"
+#include "lbm/communication/PdfFieldPackInfo.h"
+#include "lbm/field/AddToStorage.h"
+#include "lbm/field/PdfField.h"
+#include "lbm/lattice_model/D3Q19.h"
+#include "lbm/sweeps/CellwiseSweep.h"
+#include "lbm/sweeps/SweepWrappers.h"
+
+#include "lbm_rpd_coupling/mapping/ParticleMapping.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/MovingParticleMapping.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/boundary/SimpleBB.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/boundary/CurvedLinear.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/reconstruction/Reconstructor.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/reconstruction/ExtrapolationDirectionFinder.h"
+#include "lbm_rpd_coupling/momentum_exchange_method/reconstruction/PdfReconstructionManager.h"
+#include "lbm_rpd_coupling/utility/AddForceOnParticlesKernel.h"
+#include "lbm_rpd_coupling/utility/ParticleSelectors.h"
+#include "lbm_rpd_coupling/DataTypes.h"
+#include "lbm_rpd_coupling/utility/AverageHydrodynamicForceTorqueKernel.h"
+#include "lbm_rpd_coupling/utility/AddHydrodynamicInteractionKernel.h"
+#include "lbm_rpd_coupling/utility/ResetHydrodynamicForceTorqueKernel.h"
+#include "lbm_rpd_coupling/utility/LubricationCorrectionKernel.h"
+
+#include "mesa_pd/common/AnalyticContactDetection.h"
+#include "mesa_pd/data/ParticleAccessor.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+#include "mesa_pd/data/DataTypes.h"
+#include "mesa_pd/data/shape/HalfSpace.h"
+#include "mesa_pd/data/shape/Sphere.h"
+#include "mesa_pd/domain/BlockForestDomain.h"
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/SpringDashpot.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include <mesa_pd/kernel/AddVirtualForceAndTorque.h>
+#include <mesa_pd/kernel/VirtualMass.h>
+#include "mesa_pd/mpi/ContactFilter.h"
+#include "mesa_pd/mpi/SyncNextNeighbors.h"
+#include "mesa_pd/mpi/ReduceProperty.h"
+#include "mesa_pd/mpi/notifications/ForceTorqueNotification.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "timeloop/SweepTimeloop.h"
+
+#include "vtk/all.h"
+#include "field/vtk/all.h"
+#include "lbm/vtk/all.h"
+
+#include <functional>
+
+namespace settling_sphere
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+using walberla::uint_t;
+
+using LatticeModel_T = lbm::D3Q19< lbm::collision_model::TRT>;
+
+using Stencil_T = LatticeModel_T::Stencil;
+using PdfField_T = lbm::PdfField<LatticeModel_T>;
+
+using flag_t = walberla::uint8_t;
+using FlagField_T = FlagField<flag_t>;
+
+const uint_t FieldGhostLayers = 1;
+
+///////////
+// 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" );
+
+/////////////////////////////////////
+// BOUNDARY HANDLING CUSTOMIZATION //
+/////////////////////////////////////
+template <typename ParticleAccessor_T>
+class MyBoundaryHandling
+{
+public:
+
+   using NoSlip_T = lbm::NoSlip< LatticeModel_T, flag_t >;
+   using MO_T = lbm_rpd_coupling::CurvedLinear< LatticeModel_T, FlagField_T, ParticleAccessor_T >;
+   using Type = BoundaryHandling< FlagField_T, Stencil_T, NoSlip_T, MO_T >;
+
+
+   MyBoundaryHandling( const BlockDataID & flagFieldID, const BlockDataID & pdfFieldID,
+                       const BlockDataID & particleFieldID, const shared_ptr<ParticleAccessor_T>& ac) :
+         flagFieldID_( flagFieldID ), pdfFieldID_( pdfFieldID ), particleFieldID_( particleFieldID ), ac_( ac ) {}
+
+   Type * operator()( IBlock* const block, const StructuredBlockStorage* const storage ) const
+   {
+      WALBERLA_ASSERT_NOT_NULLPTR( block );
+      WALBERLA_ASSERT_NOT_NULLPTR( storage );
+
+      auto * flagField     = block->getData< FlagField_T >( flagFieldID_ );
+      auto *  pdfField     = block->getData< PdfField_T > ( pdfFieldID_ );
+      auto * particleField = block->getData< lbm_rpd_coupling::ParticleField_T > ( particleFieldID_ );
+
+      const auto fluid = flagField->flagExists( Fluid_Flag ) ? flagField->getFlag( Fluid_Flag ) : flagField->registerFlag( Fluid_Flag );
+
+      Type * handling = new Type( "moving obstacle boundary handling", flagField, fluid,
+                                  NoSlip_T( "NoSlip", NoSlip_Flag, pdfField ),
+                                  MO_T( "MO", MO_Flag, pdfField, flagField, particleField, ac_, fluid, *storage, *block ) );
+
+      handling->fillWithDomain( FieldGhostLayers );
+
+      return handling;
+   }
+
+private:
+
+   const BlockDataID flagFieldID_;
+   const BlockDataID pdfFieldID_;
+   const BlockDataID particleFieldID_;
+
+   shared_ptr<ParticleAccessor_T> ac_;
+};
+//*******************************************************************************************************************
+
+
+class ParticleAccessorWithShape : public mesa_pd::data::ParticleAccessor
+{
+public:
+   ParticleAccessorWithShape(const std::shared_ptr<mesa_pd::data::ParticleStorage>& ps, const std::shared_ptr<mesa_pd::data::ShapeStorage>& ss)
+         : ParticleAccessor(ps)
+         , ss_(ss)
+   {}
+
+   const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
+   walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
+   void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass() = v;}
+
+   const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
+   auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
+   void setInvInertiaBF(const size_t p_idx, const mesa_pd::Mat3& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF() = v;}
+
+   mesa_pd::data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
+   walberla::real_t getVolume(const size_t p_idx) const { return getShape(p_idx)->getVolume(); }
+private:
+   std::shared_ptr<mesa_pd::data::ShapeStorage> ss_;
+};
+
+class ParticleAccessorWithShapeAndVirtualMass : public ParticleAccessorWithShape {
+public:
+    ParticleAccessorWithShapeAndVirtualMass(std::shared_ptr<mesa_pd::data::ParticleStorage>& ps, std::shared_ptr<mesa_pd::data::ShapeStorage>& ss)
+            : ParticleAccessorWithShape(ps, ss)
+    {}
+
+    /**
+     * Returns the inverse mass of the body by also calculating and including its virtual mass.
+     *
+     * \param p_idx index of particle
+     */
+    walberla::real_t getInvMass(const size_t p_idx) const {
+        return getInvMassIncludingVirtual(p_idx);
+    }
+    const auto getInvInertiaBF(const size_t p_idx) const {
+        return getInvInertiaBFIncludingVirtual(p_idx);
+    }
+
+    walberla::real_t& getInvMassRef(const size_t p_idx) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+        return ParticleAccessorWithShape::getInvMassRef(p_idx);
+    }
+    void setInvMass(const size_t p_idx, const walberla::real_t& v) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_UNUSED(v);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+    }
+    auto& getInvInertiaBFRef(const size_t p_idx) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+        return ParticleAccessorWithShape::getInvInertiaBFRef(p_idx);
+    }
+    void setInvInertiaBF(const size_t p_idx, const mesa_pd::Mat3& v) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_UNUSED(v);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+    }
+};
+
+//*******************************************************************************************************************
+/*!\brief Evaluating the position and velocity of the sphere
+ *
+ */
+//*******************************************************************************************************************
+template< typename ParticleAccessor_T>
+class SpherePropertyLogger
+{
+public:
+   SpherePropertyLogger( const shared_ptr< ParticleAccessor_T > & ac, walberla::id_t sphereUid,
+         const std::string & fileName, bool fileIO, real_t diameter, real_t u_g, real_t t_g,
+         Vector3<real_t>& particleHydrodynamicForce, Vector3<real_t>& particleGravitationalForce,
+         Vector3<real_t>& particleVirtualForce, Vector3<real_t>& particleForce) :
+      ac_( ac ), sphereUid_( sphereUid ), fileName_( fileName ), fileIO_(fileIO), diameter_( diameter ), u_g_( u_g), t_g_(t_g),
+      position_( real_t(0) ), maxVelocity_( real_t(0) ),
+      particleHydrodynamicForce_(particleHydrodynamicForce), particleGravitationalForce_(particleGravitationalForce),
+      particleVirtualForce_(particleVirtualForce), particleForce_(particleForce)
+   {
+      if ( fileIO_ )
+      {
+         WALBERLA_ROOT_SECTION()
+         {
+            std::ofstream file;
+            file.open( fileName_.c_str() );
+            //file << "#\t posX\t posY\t posZ\t velX\t velY\t velZ\t t_n\t velZ_n\n";
+            file << "#\t posZ\t velZ\t FhZ\t FgZ\t FvZ\t FZ\t t_n\t velZ_n\n";
+            file.close();
+         }
+      }
+   }
+
+   void operator()(const uint_t timestep)
+   {
+      Vector3<real_t> pos(real_t(0));
+      Vector3<real_t> transVel(real_t(0));
+
+      auto hydForceZ = particleHydrodynamicForce_[2];
+      auto gravForceZ = particleGravitationalForce_[2];
+      auto virtForceZ = particleVirtualForce_[2];
+      auto forceZ = particleForce_[2];
+
+      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);
+            transVel = ac_->getLinearVelocity(idx);
+         }
+      }
+
+      WALBERLA_MPI_SECTION()
+      {
+         //mpi::allReduceInplace( pos[0], mpi::SUM ); // not required
+         //mpi::allReduceInplace( pos[1], mpi::SUM ); // not required
+         mpi::allReduceInplace( pos[2], mpi::SUM );
+
+         //mpi::allReduceInplace( transVel[0], mpi::SUM ); // not required
+         //mpi::allReduceInplace( transVel[1], mpi::SUM ); // not required
+         mpi::reduceInplace( transVel[2], mpi::SUM );
+
+         mpi::reduceInplace( hydForceZ, mpi::SUM );
+         mpi::reduceInplace( gravForceZ, mpi::SUM );
+         mpi::reduceInplace( virtForceZ, mpi::SUM );
+         mpi::reduceInplace( forceZ, mpi::SUM );
+      }
+
+      position_ = pos[2];
+      maxVelocity_ = std::max(maxVelocity_, -transVel[2]);
+
+      if( fileIO_ )
+         writeToFile( timestep, pos, transVel, hydForceZ, gravForceZ, virtForceZ, forceZ);
+   }
+
+   real_t getPosition() const
+   {
+      return position_;
+   }
+
+   real_t getMaxVelocity() const
+   {
+      return maxVelocity_;
+   }
+
+private:
+   void writeToFile( const uint_t timestep, const Vector3<real_t> & position, const Vector3<real_t> & velocity,
+           const real_t hydForceZ, const real_t gravForceZ, const real_t virtForceZ, const real_t forceZ)
+   {
+      WALBERLA_ROOT_SECTION()
+      {
+         std::ofstream file;
+         file.open( fileName_.c_str(), std::ofstream::app );
+
+         file << timestep
+              << "\t" << position[2]
+              << "\t" << velocity[2]
+              << "\t" << hydForceZ
+              << "\t" << gravForceZ
+              << "\t" << virtForceZ
+              << "\t" << forceZ
+              << "\t" << real_t(timestep)/t_g_ << "\t" << velocity[2]/u_g_
+              << "\n";
+         file.close();
+      }
+   }
+
+   shared_ptr< ParticleAccessor_T > ac_;
+   const walberla::id_t sphereUid_;
+   std::string fileName_;
+   bool fileIO_;
+   real_t diameter_;
+   real_t u_g_;
+   real_t t_g_;
+
+   real_t position_;
+   real_t maxVelocity_;
+
+   Vector3<real_t>& particleHydrodynamicForce_;
+   Vector3<real_t>& particleGravitationalForce_;
+   Vector3<real_t>& particleVirtualForce_;
+   Vector3<real_t>& particleForce_;
+};
+
+
+void createPlaneSetup(const shared_ptr<mesa_pd::data::ParticleStorage> & ps, const shared_ptr<mesa_pd::data::ShapeStorage> & ss, const math::AABB & simulationDomain)
+{
+   // create bounding planes
+   mesa_pd::data::Particle p0 = *ps->create(true);
+   p0.setPosition(simulationDomain.minCorner());
+   p0.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,1) ));
+   p0.setOwner(mpi::MPIManager::instance()->rank());
+   p0.setType(0);
+   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::INFINITE_MASS);
+
+   mesa_pd::data::Particle p1 = *ps->create(true);
+   p1.setPosition(simulationDomain.maxCorner());
+   p1.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,0,-1) ));
+   p1.setOwner(mpi::MPIManager::instance()->rank());
+   p1.setType(0);
+   mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p1.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p2 = *ps->create(true);
+   p2.setPosition(simulationDomain.minCorner());
+   p2.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(1,0,0) ));
+   p2.setOwner(mpi::MPIManager::instance()->rank());
+   p2.setType(0);
+   mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p2.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p3 = *ps->create(true);
+   p3.setPosition(simulationDomain.maxCorner());
+   p3.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(-1,0,0) ));
+   p3.setOwner(mpi::MPIManager::instance()->rank());
+   p3.setType(0);
+   mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p3.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p4 = *ps->create(true);
+   p4.setPosition(simulationDomain.minCorner());
+   p4.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,1,0) ));
+   p4.setOwner(mpi::MPIManager::instance()->rank());
+   p4.setType(0);
+   mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p4.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+   mesa_pd::data::Particle p5 = *ps->create(true);
+   p5.setPosition(simulationDomain.maxCorner());
+   p5.setShapeID(ss->create<mesa_pd::data::HalfSpace>( Vector3<real_t>(0,-1,0) ));
+   p5.setOwner(mpi::MPIManager::instance()->rank());
+   p5.setType(0);
+   mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE);
+   mesa_pd::data::particle_flags::set(p5.getFlagsRef(), mesa_pd::data::particle_flags::INFINITE_MASS);
+
+}
+
+//////////
+// MAIN //
+//////////
+
+//*******************************************************************************************************************
+/*!\brief Testcase that simulates the settling of a sphere inside a rectangular column filled with viscous fluid
+ *
+ * see: ten Cate, Nieuwstad, Derksen, Van den Akker - "Particle imaging velocimetry experiments and lattice-Boltzmann
+ * simulations on a single sphere settling under gravity" (2002), Physics of Fluids, doi: 10.1063/1.1512918
+ */
+//*******************************************************************************************************************
+
+int main( int argc, char **argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   ///////////////////
+   // Customization //
+   ///////////////////
+
+   // simulation control
+   bool shortrun = false;
+   bool funcTest = false;
+   bool fileIO = false;
+   uint_t vtkIOFreq = 0;
+   std::string baseFolder = "vtk_out_SettlingSphere";
+
+   // physical setup
+   uint_t fluidType = 1;
+
+   //numerical parameters
+   uint_t numberOfCellsInHorizontalDirection = uint_t(135);
+   bool averageForceTorqueOverTwoTimSteps = true;
+   uint_t numRPDSubCycles = uint_t(1);
+   bool useVelocityVerlet = false;
+   bool useEANReconstructor = false;
+   bool useExtReconstructor = false;
+
+   auto u_g_target = real_t(0);
+   bool useVirtualMass = false;
+   bool large = false;
+   auto relaxationTime = real_t(0);
+   auto densitySphere = real_t(1.5); // 0.5 for rising, 1.5 for settling
+   auto timesteps = uint_t(2000);
+   auto diameterOverride = real_t(-1);
+   auto galileoNumber = real_t(170);
+
+   for( int i = 1; i < argc; ++i )
+   {
+      //if( std::strcmp( argv[i], "--shortrun" )            == 0 ) { shortrun = true; continue; }
+      //if( std::strcmp( argv[i], "--funcTest" )            == 0 ) { funcTest = true; continue; }
+      if( std::strcmp( argv[i], "--fileIO" )              == 0 ) { fileIO = true; continue; }
+      if( std::strcmp( argv[i], "--vtkIOFreq" )           == 0 ) { vtkIOFreq = uint_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--fluidType" )           == 0 ) { fluidType = uint_c( std::atof( argv[++i] ) ); continue; }
+      //if( std::strcmp( argv[i], "--numRPDSubCycles" )     == 0 ) { numRPDSubCycles = uint_c( std::atof( argv[++i] ) ); continue; }
+      //if( std::strcmp( argv[i], "--resolution" )          == 0 ) { numberOfCellsInHorizontalDirection = uint_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--noForceAveraging" )    == 0 ) { averageForceTorqueOverTwoTimSteps = false; continue; }
+      if( std::strcmp( argv[i], "--baseFolder" )          == 0 ) { baseFolder = argv[++i]; continue; }
+      if( std::strcmp( argv[i], "--useVV" )               == 0 ) { useVelocityVerlet = true; continue; }
+      if( std::strcmp( argv[i], "--useEANReconstructor" ) == 0 ) { useEANReconstructor = true; continue; }
+      if( std::strcmp( argv[i], "--useExtReconstructor" ) == 0 ) { useExtReconstructor = true; continue; }
+      if( std::strcmp( argv[i], "--rising" )              == 0 ) { densitySphere = real_t(0.5); continue; }
+      if( std::strcmp( argv[i], "--sphereDensity" )       == 0 ) { densitySphere = real_c( std::stod( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--relaxationTime")       == 0 ) { relaxationTime = real_c( std::stod( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--u_g" )                 == 0 ) { u_g_target = real_c( std::stod( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--virtualMass" )         == 0 ) { useVirtualMass = true; continue; }
+      if( std::strcmp( argv[i], "--large" )               == 0 ) { large = true; continue; }
+      if( std::strcmp( argv[i], "--timesteps" )           == 0 ) { timesteps = uint_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--diameter" )            == 0 ) { diameterOverride = real_c( std::stod( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--galileoNumber" )       == 0 ) { galileoNumber = real_c( std::stod( argv[++i] ) ); continue; }
+      WALBERLA_ABORT("Unrecognized command line argument found: " << argv[i]);
+   }
+
+   WALBERLA_CHECK(!(useEANReconstructor && useExtReconstructor), "Can not use two reconstructors!");
+   WALBERLA_CHECK(!realIsIdentical(u_g_target, real_t(0)) || !realIsIdentical(relaxationTime, real_t(0)), "Either u_g or relaxationTime has to be set!");
+   WALBERLA_CHECK(!(realIsIdentical(u_g_target, real_t(0)) && realIsIdentical(relaxationTime, real_t(0))), "Can not set both u_g and relaxationTime!");
+
+   if( funcTest )
+   {
+      walberla::logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::WARNING);
+   }
+
+   if( fileIO )
+   {
+      // create base directory if it does not yet exist
+      filesystem::path tpath( baseFolder );
+      if( !filesystem::exists( tpath ) )
+         filesystem::create_directory( tpath );
+   }
+
+   //////////////////////////
+   // NUMERICAL PARAMETERS //
+   //////////////////////////
+
+
+   const auto densityFluid = real_t(1);
+   // relaxation time via arguments
+
+   const auto diameter = (diameterOverride > 0) ? diameterOverride : (large ? real_t(40) : real_t(20));
+   const auto dt = real_t(5e-4);
+
+   const auto Nx = uint_t(real_t(6.4)*diameter);
+   const auto Ny = uint_t(real_t(6.4)*diameter);
+   const auto Nz = uint_t(real_t(6.4)*diameter * (large ? real_t(2) : real_t(1)));
+
+   real_t viscosity, omega;
+   if (realIsIdentical(relaxationTime, real_t(0))) {
+      // calculate relaxation time from u_g_target
+      viscosity = u_g_target * diameter / galileoNumber;
+      omega = lbm::collision_model::omegaFromViscosity(viscosity);
+      relaxationTime = real_t(1)/omega;
+   } else {
+      omega = real_t(1)/relaxationTime;
+      viscosity = lbm::collision_model::viscosityFromOmega(omega);
+   }
+
+   if (relaxationTime <= real_t(0.505)) {
+      WALBERLA_ABORT("relaxationTime " << relaxationTime << " smaller than 0.505, simulation will be unstable. Aborting.");
+   } else if (relaxationTime < real_t(0.51) && relaxationTime > real_t(0.505)) {
+      WALBERLA_LOG_WARNING_ON_ROOT("relaxationTime " << relaxationTime << " between 0.505 and 0.51, simulation might behave unstable.");
+   }
+
+   //const auto omega = real_t(1)/relaxationTime;
+   //const auto viscosity = lbm::collision_model::viscosityFromOmega(omega);
+
+   const auto d3 = diameter*diameter*diameter;
+   const auto G2 = galileoNumber*galileoNumber;
+   const auto v2 = viscosity*viscosity;
+   const auto gravitationalAcceleration = G2*v2 / (abs(densitySphere-real_t(1))*d3);
+
+   const Vector3<uint_t> domainSize( Nx, Ny, Nz);
+   Vector3<real_t> initialPosition;
+   if (large) {
+       initialPosition = Vector3<real_t>(real_t(3.2) * diameter,
+                                         real_t(3.2) * diameter,
+                                         real_t(0.6) * diameter);
+   } else {
+       initialPosition = Vector3<real_t>(real_t(0.5) * real_c(domainSize[0]),
+                                         real_t(0.5) * real_c(domainSize[1]),
+                                         real_t(0.5) * real_c(domainSize[2]));
+   }
+
+   const real_t sphereVolume = math::M_PI / real_t(6) * diameter * diameter * diameter;
+
+   const auto dx = real_t(1);
+
+   // virtual mass parameters
+   const auto criticalDensity = real_t(1.1); // particle density, which just about works with native walberla
+   const auto C_v = (criticalDensity - densitySphere) / densityFluid;// real_t(0.5);
+   const auto C_v_omega = C_v;
+
+   const auto u_g = std::sqrt(std::fabs(densitySphere - real_t(1)) * gravitationalAcceleration * diameter);
+   const auto t_g = std::sqrt(diameter / (std::fabs(densitySphere - real_t(1)) * gravitationalAcceleration));
+   if (u_g > real_t(0.1)){
+      WALBERLA_LOG_WARNING_ON_ROOT("u_g " << u_g << " > 0.1, simulation will most likely behave unstable.");
+   } else if (u_g > real_t(0.05)) {
+      WALBERLA_LOG_WARNING_ON_ROOT("u_g " << u_g << " > 0.05, simulation might behave unstable.");
+   }
+
+   const bool rising = densitySphere < real_t(1.0);
+
+   //WALBERLA_LOG_INFO_ON_ROOT(" - dx_SI = " << dx_SI << ", dt_SI = " << dt_SI);
+   WALBERLA_LOG_INFO_ON_ROOT("Setup (in simulation, i.e. lattice, units):");
+   WALBERLA_LOG_INFO_ON_ROOT(" - domain size = " << domainSize << (large ? " (large simulation)" : ""));
+   WALBERLA_LOG_INFO_ON_ROOT(" - sphere: diameter = " << diameter << ", density = " << densitySphere );
+   WALBERLA_LOG_INFO_ON_ROOT(" - fluid: density = " << densityFluid << ", relaxation time (tau) = " << relaxationTime << ", kin. visc = " << viscosity );
+   WALBERLA_LOG_INFO_ON_ROOT(" - galileo number = " << galileoNumber );
+   WALBERLA_LOG_INFO_ON_ROOT(" - gravitational velocity (u_g) = " << u_g);
+   WALBERLA_LOG_INFO_ON_ROOT(" - gravitational acceleration = " << gravitationalAcceleration );
+   //WALBERLA_LOG_INFO_ON_ROOT(" - expected settling velocity = " << expectedSettlingVelocity << " --> Re_p = " << expectedSettlingVelocity * diameter / viscosity );
+   WALBERLA_LOG_INFO_ON_ROOT(" - integrator = " << (useVelocityVerlet ? "Velocity Verlet" : "Explicit Euler") );
+   WALBERLA_LOG_INFO_ON_ROOT(" - use virtual mass = " << useVirtualMass );
+   WALBERLA_LOG_INFO_ON_ROOT(" - C_v = " << C_v );
+   WALBERLA_LOG_INFO_ON_ROOT(" - relaxation time = " << relaxationTime );
+   WALBERLA_LOG_INFO_ON_ROOT(" - reconstructor = " << (useEANReconstructor ? "EAN" : "") << (useExtReconstructor ? "Ext" : "") << (!useEANReconstructor && !useExtReconstructor ? "default" : ""))
+   WALBERLA_LOG_INFO_ON_ROOT(" - timesteps = " << timesteps);
+
+   if( vtkIOFreq > 0 )
+   {
+      WALBERLA_LOG_INFO_ON_ROOT(" - writing vtk files to folder \"" << baseFolder << "\" with frequency " << vtkIOFreq);
+   }
+
+   ///////////////////////////
+   // BLOCK STRUCTURE SETUP //
+   ///////////////////////////
+
+   Vector3<uint_t> numberOfBlocksPerDirection( uint_t(1), uint_t(1), uint_t(4) );
+   Vector3<uint_t> cellsPerBlockPerDirection( domainSize[0] / numberOfBlocksPerDirection[0],
+                                              domainSize[1] / numberOfBlocksPerDirection[1],
+                                              domainSize[2] / numberOfBlocksPerDirection[2] );
+   for( uint_t i = 0; i < 3; ++i ) {
+      WALBERLA_CHECK_EQUAL(cellsPerBlockPerDirection[i] * numberOfBlocksPerDirection[i], domainSize[i],
+                           "Unmatching domain decomposition in direction " << i << "!");
+   }
+
+   auto blocks = blockforest::createUniformBlockGrid( numberOfBlocksPerDirection[0], numberOfBlocksPerDirection[1], numberOfBlocksPerDirection[2],
+                                                      cellsPerBlockPerDirection[0], cellsPerBlockPerDirection[1], cellsPerBlockPerDirection[2], dx,
+                                                      0, false, false,
+                                                      true, true, true, //periodicity
+                                                      false );
+
+   WALBERLA_LOG_INFO_ON_ROOT("Domain decomposition:");
+   WALBERLA_LOG_INFO_ON_ROOT(" - blocks per direction = " << numberOfBlocksPerDirection );
+   WALBERLA_LOG_INFO_ON_ROOT(" - cells per block = " << cellsPerBlockPerDirection );
+
+   //write domain decomposition to file
+   if( vtkIOFreq > 0 )
+   {
+      vtk::writeDomainDecomposition( blocks, "initial_domain_decomposition", baseFolder );
+   }
+
+   //////////////////
+   // RPD COUPLING //
+   //////////////////
+
+   auto rpdDomain = std::make_shared<mesa_pd::domain::BlockForestDomain>(blocks->getBlockForestPointer());
+
+   //init data structures
+   auto ps = walberla::make_shared<mesa_pd::data::ParticleStorage>(1);
+   auto ss = walberla::make_shared<mesa_pd::data::ShapeStorage>();
+   using ParticleAccessor_T = ParticleAccessorWithShape;
+   auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
+   auto virtualMassAccessor = walberla::make_shared<ParticleAccessorWithShapeAndVirtualMass>(ps, ss);
+
+   // bounding planes
+   // excluded for virtual mass test case (periodic boundary conditions)
+   //createPlaneSetup(ps,ss,blocks->getDomain());
+
+   // create sphere and store Uid
+   auto sphereShape = ss->create<mesa_pd::data::Sphere>( diameter * real_t(0.5) );
+   ss->shapes[sphereShape]->updateMassAndInertia(densitySphere);
+
+   walberla::id_t sphereUid = 0;
+   if (rpdDomain->isContainedInProcessSubdomain( uint_c(mpi::MPIManager::instance()->rank()), initialPosition ))
+   {
+      mesa_pd::data::Particle&& p = *ps->create();
+      p.setPosition(initialPosition);
+      p.setInteractionRadius(diameter * real_t(0.5));
+      p.setOwner(mpi::MPIManager::instance()->rank());
+      p.setShapeID(sphereShape);
+      sphereUid = p.getUid();
+   }
+   mpi::allReduceInplace(sphereUid, mpi::SUM);
+
+   ///////////////////////
+   // ADD DATA TO BLOCKS //
+   ////////////////////////
+
+   // create the lattice model
+   LatticeModel_T latticeModel = LatticeModel_T( lbm::collision_model::TRT::constructWithMagicNumber( real_t(1) / relaxationTime ) );
+
+   // add PDF field
+   BlockDataID pdfFieldID = lbm::addPdfFieldToStorage< LatticeModel_T >( blocks, "pdf field (zyxf)", latticeModel,
+                                                                         Vector3< real_t >( real_t(0) ), real_t(1),
+                                                                         uint_t(1), field::zyxf );
+   // add flag field
+   BlockDataID flagFieldID = field::addFlagFieldToStorage<FlagField_T>( blocks, "flag field" );
+
+   // add particle field
+   BlockDataID particleFieldID = field::addToStorage<lbm_rpd_coupling::ParticleField_T>( blocks, "particle field", accessor->getInvalidUid(), field::zyxf, FieldGhostLayers );
+
+   // add boundary handling
+   using BoundaryHandling_T = MyBoundaryHandling<ParticleAccessor_T>::Type;
+   BlockDataID boundaryHandlingID = blocks->addStructuredBlockData< BoundaryHandling_T >(MyBoundaryHandling<ParticleAccessor_T>( flagFieldID, pdfFieldID, particleFieldID, accessor), "boundary handling" );
+
+   // set up RPD functionality
+   std::function<void(void)> syncCall = [ps,rpdDomain](){
+      const real_t overlap = real_t( 1.5 );
+      mesa_pd::mpi::SyncNextNeighbors syncNextNeighborFunc;
+      syncNextNeighborFunc(*ps, *rpdDomain, overlap);
+   };
+
+   syncCall();
+
+   mesa_pd::kernel::ExplicitEulerWithShape explEulerIntegrator(real_t(1)/real_t(numRPDSubCycles));
+   mesa_pd::kernel::VelocityVerletPreForceUpdate  vvIntegratorPreForce(real_t(1)/real_t(numRPDSubCycles));
+   mesa_pd::kernel::VelocityVerletPostForceUpdate vvIntegratorPostForce(real_t(1)/real_t(numRPDSubCycles));
+
+   mesa_pd::kernel::SpringDashpot collisionResponse(1);
+   mesa_pd::mpi::ReduceProperty reduceProperty;
+
+   // set up coupling functionality
+   Vector3<real_t> gravitationalForce( real_t(0), real_t(0), -(densitySphere - densityFluid) * gravitationalAcceleration * sphereVolume );
+   lbm_rpd_coupling::AddForceOnParticlesKernel<lbm_rpd_coupling::LocalParticlesSelector> addGravitationalForce(gravitationalForce);
+
+   lbm_rpd_coupling::RegularParticlesSelector sphereSelector;
+   lbm_rpd_coupling::AddHydrodynamicInteractionKernel<decltype(sphereSelector)> addHydrodynamicInteraction(sphereSelector);
+
+   lbm_rpd_coupling::ResetHydrodynamicForceTorqueKernel<> resetHydrodynamicForceTorque;
+
+   lbm_rpd_coupling::AverageHydrodynamicForceTorqueKernel averageHydrodynamicForceTorque;
+
+   // virtual mass approach
+   mesa_pd::kernel::VirtualMass virtualMass(C_v, C_v_omega, densityFluid);
+   mesa_pd::kernel::AddVirtualForceAndTorque addVirtualForceAndTorque;
+
+   if ( useVirtualMass ) ps->forEachParticle(false, virtualMass, *accessor);
+
+   ///////////////
+   // TIME LOOP //
+   ///////////////
+
+   // map planes into the LBM simulation -> act as no-slip boundaries
+   lbm_rpd_coupling::mapParticles<BoundaryHandling_T>(*blocks, boundaryHandlingID, *accessor, NoSlip_Flag, lbm_rpd_coupling::GlobalParticlesSelector());
+
+   // map particles into the LBM simulation
+   lbm_rpd_coupling::mapMovingParticles<BoundaryHandling_T>(*blocks, boundaryHandlingID, particleFieldID, *accessor, MO_Flag, sphereSelector);
+
+
+   // setup of the LBM communication for synchronizing the pdf field between neighboring blocks
+   std::function< void () > commFunction;
+   blockforest::communication::UniformBufferedScheme< Stencil_T > scheme( blocks );
+   scheme.addPackInfo( make_shared< lbm::PdfFieldPackInfo< LatticeModel_T > >( pdfFieldID ) );
+   commFunction = scheme;
+
+   // create the timeloop
+   SweepTimeloop timeloop( blocks->getBlockStorage(), timesteps );
+
+   auto bhSweep = BoundaryHandling_T::getBlockSweep( boundaryHandlingID );
+   auto lbmSweep = lbm::makeCellwiseSweep< LatticeModel_T, FlagField_T >( pdfFieldID, flagFieldID, Fluid_Flag );
+
+   timeloop.addFuncBeforeTimeStep( RemainingTimeLogger( timeloop.getNrOfTimeSteps() ), "Remaining Time Logger" );
+
+   // vtk output
+   if( vtkIOFreq != uint_t(0) )
+   {
+      // spheres
+      auto particleVtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+      particleVtkOutput->addOutput<mesa_pd::data::SelectParticleOwner>("owner");
+      particleVtkOutput->addOutput<mesa_pd::data::SelectParticleLinearVelocity>("velocity");
+      auto particleVtkWriter = vtk::createVTKOutput_PointData(particleVtkOutput, "Particles", vtkIOFreq, baseFolder, "simulation_step");
+      timeloop.addFuncBeforeTimeStep( vtk::writeFiles( particleVtkWriter ), "VTK (sphere data)" );
+
+      // flag field (written only once in the first time step, ghost layers are also written)
+      //auto flagFieldVTK = vtk::createVTKOutput_BlockData( blocks, "flag_field", timesteps, FieldGhostLayers, false, baseFolder );
+      //flagFieldVTK->addCellDataWriter( make_shared< field::VTKWriter< FlagField_T > >( flagFieldID, "FlagField" ) );
+      //timeloop.addFuncBeforeTimeStep( vtk::writeFiles( flagFieldVTK ), "VTK (flag field data)" );
+
+      // pdf field
+      auto pdfFieldVTK = vtk::createVTKOutput_BlockData( blocks, "fluid_field", vtkIOFreq, 0, false, baseFolder );
+
+      blockforest::communication::UniformBufferedScheme< stencil::D3Q27 > pdfGhostLayerSync( blocks );
+      pdfGhostLayerSync.addPackInfo( make_shared< field::communication::PackInfo< PdfField_T > >( pdfFieldID ) );
+      pdfFieldVTK->addBeforeFunction( pdfGhostLayerSync );
+
+      field::FlagFieldCellFilter< FlagField_T > fluidFilter( flagFieldID );
+      fluidFilter.addFlag( Fluid_Flag );
+      pdfFieldVTK->addCellInclusionFilter( fluidFilter );
+
+      pdfFieldVTK->addCellDataWriter( make_shared< lbm::VelocityVTKWriter< LatticeModel_T, float > >( pdfFieldID, "VelocityFromPDF" ) );
+      pdfFieldVTK->addCellDataWriter( make_shared< lbm::DensityVTKWriter < LatticeModel_T, float > >( pdfFieldID, "DensityFromPDF" ) );
+
+      timeloop.addFuncBeforeTimeStep( vtk::writeFiles( pdfFieldVTK ), "VTK (fluid field data)" );
+   }
+
+   // sweep for updating the particle mapping into the LBM simulation
+   timeloop.add() << Sweep( lbm_rpd_coupling::makeMovingParticleMapping<BoundaryHandling_T>(blocks, boundaryHandlingID, particleFieldID, accessor, MO_Flag, FormerMO_Flag, sphereSelector), "Particle Mapping" );
+
+   // sweep for restoring PDFs in cells previously occupied by particles
+   if( useEANReconstructor )
+   {
+      auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_rpd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
+      auto equilibriumAndNonEquilibriumSphereNormalReconstructor = lbm_rpd_coupling::makeEquilibriumAndNonEquilibriumReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder);
+      auto reconstructionManager = lbm_rpd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, equilibriumAndNonEquilibriumSphereNormalReconstructor);
+
+      timeloop.add() << Sweep( makeSharedSweep(reconstructionManager), "PDF Restore" );
+   } else if( useExtReconstructor ) {
+      auto sphereNormalExtrapolationDirectionFinder = make_shared<lbm_rpd_coupling::SphereNormalExtrapolationDirectionFinder>(blocks);
+      auto extrapolationSphereNormalReconstructor = lbm_rpd_coupling::makeExtrapolationReconstructor<BoundaryHandling_T>(blocks, boundaryHandlingID, sphereNormalExtrapolationDirectionFinder, true);
+      timeloop.add() << Sweep( makeSharedSweep(lbm_rpd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag, extrapolationSphereNormalReconstructor) ), "PDF Restore" );
+   } else {
+      timeloop.add() << Sweep( makeSharedSweep(lbm_rpd_coupling::makePdfReconstructionManager<PdfField_T,BoundaryHandling_T>(blocks, pdfFieldID, boundaryHandlingID, particleFieldID, accessor, FormerMO_Flag, Fluid_Flag) ), "PDF Restore" );
+   }
+
+   bool useStreamCollide = true;
+
+   if( useStreamCollide )
+   {
+      // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
+      timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
+                     << Sweep(bhSweep, "Boundary Handling" );
+
+      // stream + collide LBM step
+      timeloop.add() << Sweep( makeSharedSweep( lbmSweep ), "cell-wise LB sweep" );
+   }
+   else
+   {
+      // collide LBM step
+      timeloop.add() << Sweep( makeCollideSweep( lbmSweep ), "cell-wise collide LB sweep" );
+
+      // add LBM communication function and boundary handling sweep (does the hydro force calculations and the no-slip treatment)
+      timeloop.add() << BeforeFunction( commFunction, "LBM Communication" )
+                     << Sweep( BoundaryHandling_T::getBlockSweep( boundaryHandlingID ), "Boundary Handling" );
+
+      // stream LBM step
+      timeloop.add() << Sweep( makeStreamSweep( lbmSweep ), "cell-wise stream LB sweep" );
+   }
+
+   // evaluation functionality
+   std::string loggingFileName( baseFolder + "/LoggingSettlingSphere_");
+   loggingFileName += std::to_string(fluidType);
+   loggingFileName += ".txt";
+   if( fileIO  )
+   {
+      WALBERLA_LOG_INFO_ON_ROOT(" - writing logging output to file \"" << loggingFileName << "\"");
+   }
+   Vector3<real_t> particleHydrodynamicForce, particleGravitationalForce, particleVirtualForce, particleForce;
+   SpherePropertyLogger<ParticleAccessor_T> logger( accessor, sphereUid, loggingFileName, fileIO, diameter, u_g, t_g,
+                                                    particleHydrodynamicForce, particleGravitationalForce, particleVirtualForce, particleForce);
+
+
+   ////////////////////////
+   // EXECUTE SIMULATION //
+   ////////////////////////
+
+   WcTimingPool timeloopTiming;
+
+   real_t terminationPosition;
+   if (large) {
+       // when the sphere has covered a distance of L_z - diameter
+       if (rising) {
+           terminationPosition = initialPosition[2] + real_t(domainSize[2]) - diameter;
+       } else {
+           terminationPosition = initialPosition[2] + diameter;
+       }
+   } else {
+       // when the sphere is only a diameter away from the initial position again
+       if (rising) {
+           terminationPosition = real_t(initialPosition[2]) - diameter;
+       } else {
+           terminationPosition = real_t(initialPosition[2]) + diameter;
+       }
+   }
+   WALBERLA_LOG_INFO_ON_ROOT(" - termination position = " << terminationPosition );
+
+   //  required for determining termination
+   bool hasCrossedBoundary = false;
+
+   const bool useOpenMP = false;
+
+   // time loop
+   for (uint_t i = 0; i < timesteps; ++i )
+   {
+      particleVirtualForce.set(real_t(0), real_t(0), real_t(0));
+      particleForce.set(real_t(0), real_t(0), real_t(0));
+
+      // perform a single simulation step -> this contains LBM and setting of the hydrodynamic interactions
+      timeloop.singleStep( timeloopTiming );
+
+      timeloopTiming["RPD"].start();
+
+      if( averageForceTorqueOverTwoTimSteps && i!= 0)
+      {
+         ps->forEachParticle(useOpenMP, averageHydrodynamicForceTorque, *accessor );
+      }
+
+      //for(auto subCycle = uint_t(0); subCycle < numRPDSubCycles; ++subCycle )
+      //{
+
+         if( useVelocityVerlet ) {
+            if ( useVirtualMass ) ps->forEachParticle(useOpenMP, vvIntegratorPreForce, *virtualMassAccessor);
+            else ps->forEachParticle(useOpenMP, vvIntegratorPreForce, *accessor);
+
+            syncCall();
+         }
+
+         ps->forEachParticle(useOpenMP, addHydrodynamicInteraction, *accessor );
+         ps->forEachParticle(useOpenMP, [&particleHydrodynamicForce](size_t idx, auto& ac) {
+             particleHydrodynamicForce = ac.getForce(idx);
+         }, *accessor);
+         ps->forEachParticle(useOpenMP, addGravitationalForce, *accessor );
+         ps->forEachParticle(useOpenMP, [&particleGravitationalForce, &particleHydrodynamicForce](size_t idx, auto& ac) {
+             particleGravitationalForce = ac.getForce(idx) - particleHydrodynamicForce;
+         }, *accessor);
+
+         reduceProperty.operator()<mesa_pd::ForceTorqueNotification>(*ps);
+
+
+         if ( useVirtualMass ) {
+            Vector3<real_t> particleOldForce;
+            ps->forEachParticle(useOpenMP, [&particleOldForce](size_t idx, auto& ac) {
+                lbm_rpd_coupling::LocalParticlesSelector localSelector;
+                if(localSelector(idx, ac)){
+                   particleOldForce = ac.getForce(idx);
+                }
+            }, *accessor);
+
+            ps->forEachParticle(useOpenMP, addVirtualForceAndTorque, *accessor);
+            ps->forEachParticle(useOpenMP, [&particleOldForce, &particleVirtualForce](size_t idx, auto& ac) {
+               lbm_rpd_coupling::LocalParticlesSelector localSelector;
+               if(localSelector(idx, ac)){
+                  particleVirtualForce = ac.getForce(idx) - particleOldForce;
+               }
+            }, *accessor);
+
+            ps->forEachParticle(useOpenMP, [&particleForce](size_t idx, auto& ac) {
+               lbm_rpd_coupling::LocalParticlesSelector localSelector;
+               if(localSelector(idx, ac)){
+                  particleForce = ac.getForce(idx);
+               }
+            }, *accessor);
+
+            if( useVelocityVerlet ) ps->forEachParticle(useOpenMP, vvIntegratorPostForce, *virtualMassAccessor);
+            else ps->forEachParticle(useOpenMP, explEulerIntegrator, *virtualMassAccessor);
+         } else {
+            ps->forEachParticle(useOpenMP, [&particleForce](size_t idx, auto& ac) {
+                particleForce = ac.getForce(idx);
+            }, *accessor);
+
+            if( useVelocityVerlet ) ps->forEachParticle(useOpenMP, vvIntegratorPostForce, *accessor);
+            else ps->forEachParticle(useOpenMP, explEulerIntegrator, *accessor);
+         }
+
+         syncCall();
+      //}
+
+      timeloopTiming["RPD"].end();
+
+      // evaluation
+      timeloopTiming["Logging"].start();
+      logger(i);
+      timeloopTiming["Logging"].end();
+
+      // reset after logging here
+      ps->forEachParticle(useOpenMP, resetHydrodynamicForceTorque, *accessor );
+
+      // check for termination
+      // check if the sphere crossed a boundary
+      if (!hasCrossedBoundary && ((rising && logger.getPosition() < terminationPosition) || (!rising && logger.getPosition() > terminationPosition))) {
+          hasCrossedBoundary = true;
+      }
+      if ((!large && hasCrossedBoundary && ((rising && logger.getPosition() > terminationPosition) || (!rising && logger.getPosition() < terminationPosition)))
+        || (large && ((rising && logger.getPosition() > terminationPosition) || (!rising && hasCrossedBoundary && logger.getPosition() < terminationPosition))))
+      {
+          // end of simulation reached after crossing a boundary and then passing the termination position
+          WALBERLA_LOG_INFO_ON_ROOT("Sphere reached terminal position " << logger.getPosition() << " after " << i << " timesteps!");
+          break;
+      }
+   }
+
+   timeloopTiming.logResultOnRoot();
+
+   // check the result
+   /*if ( !funcTest && !shortrun )
+   {
+      real_t relErr = std::fabs( expectedSettlingVelocity - logger.getMaxVelocity()) / expectedSettlingVelocity;
+      WALBERLA_LOG_INFO_ON_ROOT( "Expected maximum settling velocity: " << expectedSettlingVelocity );
+      WALBERLA_LOG_INFO_ON_ROOT( "Simulated maximum settling velocity: " << logger.getMaxVelocity() );
+      WALBERLA_LOG_INFO_ON_ROOT( "Relative error: " << relErr );
+
+      // the relative error has to be below 10%
+      WALBERLA_CHECK_LESS( relErr, real_t(0.1) );
+   }*/
+
+   return EXIT_SUCCESS;
+}
+
+} // namespace settling_sphere_mem
+
+int main( int argc, char **argv ){
+   settling_sphere::main(argc, argv);
+}
diff --git a/apps/virtualmass/VirtualMass.cpp b/apps/virtualmass/VirtualMass.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9311b546b4e1c5d4154e1333329947f7be901134
--- /dev/null
+++ b/apps/virtualmass/VirtualMass.cpp
@@ -0,0 +1,201 @@
+//======================================================================================================================
+//
+//  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   VirtualMass.cpp
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/data/ShapeStorage.h>
+#include <mesa_pd/data/shape/Sphere.h>
+
+#include <mesa_pd/kernel/ExplicitEuler.h>
+#include <mesa_pd/kernel/ExplicitEulerWithShape.h>
+#include <mesa_pd/kernel/AddVirtualForceAndTorque.h>
+#include <mesa_pd/kernel/VirtualMass.h>
+
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+
+#include <iostream>
+
+namespace walberla {
+
+using namespace walberla::mesa_pd;
+
+class ParticleAccessorWithShape : public data::ParticleAccessor
+{
+public:
+    ParticleAccessorWithShape(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
+            : ParticleAccessor(ps)
+            , ss_(ss)
+    {}
+
+    const walberla::real_t& getInvMass(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
+    walberla::real_t& getInvMassRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass();}
+    void setInvMass(const size_t p_idx, const walberla::real_t& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvMass() = v;}
+
+    const auto& getInvInertiaBF(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
+    auto& getInvInertiaBFRef(const size_t p_idx) {return ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF();}
+    void setInvInertiaBF(const size_t p_idx, const Mat3& v) { ss_->shapes[ps_->getShapeID(p_idx)]->getInvInertiaBF() = v;}
+
+    data::BaseShape* getShape(const size_t p_idx) const {return ss_->shapes[ps_->getShapeID(p_idx)].get();}
+    walberla::real_t getVolume(const size_t p_idx) const { return getShape(p_idx)->getVolume(); }
+private:
+    std::shared_ptr<data::ShapeStorage> ss_;
+};
+
+/**
+ * Accessor which takes the virtual mass for calls to invMass and getInvInertiaBF into account.
+ */
+class ParticleAccessorWithShapeAndVirtualMass : public ParticleAccessorWithShape {
+public:
+    ParticleAccessorWithShapeAndVirtualMass(std::shared_ptr<data::ParticleStorage>& ps, std::shared_ptr<data::ShapeStorage>& ss)
+            : ParticleAccessorWithShape(ps, ss)
+    {}
+
+    /**
+     * Returns the inverse mass of the body by also calculating and including its virtual mass.
+     *
+     * \param p_idx index of particle
+     */
+    walberla::real_t getInvMass(const size_t p_idx) const {
+        return getInvMassIncludingVirtual(p_idx);
+    }
+    const auto getInvInertiaBF(const size_t p_idx) const {
+        return getInvInertiaBFIncludingVirtual(p_idx);
+    }
+
+    walberla::real_t& getInvMassRef(const size_t p_idx) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+        return ParticleAccessorWithShape::getInvMassRef(p_idx);
+    }
+    void setInvMass(const size_t p_idx, const walberla::real_t& v) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_UNUSED(v);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+    }
+    auto& getInvInertiaBFRef(const size_t p_idx) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+        return ParticleAccessorWithShape::getInvInertiaBFRef(p_idx);
+    }
+    void setInvInertiaBF(const size_t p_idx, const Mat3& v) {
+        WALBERLA_UNUSED(p_idx);
+        WALBERLA_UNUSED(v);
+        WALBERLA_ABORT("This should probably not be called on this accessor. Not sure though, you can also comment this out and try.") //TODO
+    }
+};
+
+int main( int argc, char ** argv )
+{
+    Environment env(argc, argv);
+    WALBERLA_UNUSED(env);
+    mpi::MPIManager::instance()->useWorldComm();
+
+    auto ps = std::make_shared<data::ParticleStorage>(2);
+    auto ss = std::make_shared<data::ShapeStorage>();
+
+    ParticleAccessorWithShape accessor(ps, ss);
+    // accessor for the virtual mass approach
+    ParticleAccessorWithShapeAndVirtualMass virtualMassAccessor(ps, ss);
+
+    //initialize particles
+    const auto radius  = real_t(0.6);
+    //const auto spacing = real_t(1.0);
+    auto smallSphere = ss->create<data::Sphere>( radius );
+
+    ss->shapes[smallSphere]->updateMassAndInertia(real_t(2707));
+
+    ps->create();
+    ps->create();
+
+    auto force = Vec3(1,2,3);
+    auto torque = Vec3(1,2,3);
+
+    ps->forEachParticle(false, [force, torque](size_t idx, auto& ac) {
+        const auto linVel = Vec3(1,2,3);
+        const auto angVel = Vec3(1,2,3);
+        const auto invMass = real_t(1.23456);
+        const auto invInertiaBF = Mat3(
+                real_t(1.23456), real_t(0), real_t(0),
+                real_t(0), real_t(1.23456), real_t(0),
+                real_t(0), real_t(0), real_t(1.23456));
+
+        ac.setLinearVelocity(idx, linVel);
+        ac.setAngularVelocity(idx, angVel);
+        ac.setForce(idx, force);
+        ac.setTorque(idx, torque);
+
+        ac.setInvMass(idx, invMass);
+        ac.setInvInertiaBF(idx, invInertiaBF);
+    }, accessor);
+
+    // indices of test particles
+    const size_t plainIdx = 0; // particle using the default integration
+    const size_t virtualIdx = 1; // particle using the virtual mass approach
+
+    //init kernels
+    const auto dt = real_t(1);
+    // constants required for virtual mass approach
+    const auto C_v = real_t(0.5);
+    const auto C_v_omega = C_v;
+    const auto fluidDensity = real_t(1);
+
+    //kernel::ExplicitEuler integrator( dt );
+    kernel::ExplicitEulerWithShape integrator( dt );
+    kernel::VirtualMass virtualMass(C_v, C_v_omega, fluidDensity);
+    kernel::AddVirtualForceAndTorque addVirtualForceAndTorque;
+
+    ps->forEachParticle(false, [virtualMass](size_t idx, auto& ac) {
+        virtualMass(idx, ac);
+    }, accessor);
+
+    for (int i = 0; i < 3; ++i) {
+        // integration of particle with no virtual mass added
+        accessor.setForce(plainIdx, force);
+        accessor.setTorque(plainIdx, torque);
+        integrator(plainIdx, accessor);
+
+        // integration of particle with virtual mass added
+        accessor.setForce(virtualIdx, force);
+        accessor.setTorque(virtualIdx, torque);
+        addVirtualForceAndTorque(virtualIdx, accessor);
+        integrator(virtualIdx, virtualMassAccessor);
+
+        force += Vec3(0.1, 0.1, 0.1);
+        torque += Vec3(0.1, 0.1, 0.1);
+
+        WALBERLA_LOG_INFO("p vel: " << accessor.getLinearVelocity(plainIdx));
+        WALBERLA_LOG_INFO("p angvel: " << accessor.getAngularVelocity(plainIdx));
+        WALBERLA_LOG_INFO("v vel: " << accessor.getLinearVelocity(virtualIdx));
+        WALBERLA_LOG_INFO("v angvel: " << accessor.getAngularVelocity(virtualIdx));
+        WALBERLA_ASSERT_FLOAT_EQUAL(accessor.getLinearVelocity(plainIdx), accessor.getLinearVelocity(virtualIdx));
+        WALBERLA_ASSERT_FLOAT_EQUAL(accessor.getAngularVelocity(plainIdx), accessor.getAngularVelocity(virtualIdx));
+    }
+
+    return EXIT_SUCCESS;
+}
+
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+    return walberla::main(argc, argv);
+}
diff --git a/apps/virtualmass/mesa_pd.py b/apps/virtualmass/mesa_pd.py
new file mode 100755
index 0000000000000000000000000000000000000000..db4dbb5a0b4e95542a7c3f052a5a4bbdc89be1e7
--- /dev/null
+++ b/apps/virtualmass/mesa_pd.py
@@ -0,0 +1,159 @@
+#! /usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+from mesa_pd.accessor import Accessor
+import mesa_pd.data as data
+import mesa_pd.kernel as kernel
+import mesa_pd.mpi as mpi
+
+import argparse
+import numpy as np
+import os
+
+if __name__ == '__main__':
+   parser = argparse.ArgumentParser(description='Generate all necessary files for the waLBerla mesa_pd module.')
+   parser.add_argument('path', help='Where should the files be created?')
+   parser.add_argument("-f", "--force", help="Generate the files even if not inside a waLBerla directory.",
+                       action="store_true")
+   args = parser.parse_args()
+
+   if ((not os.path.isfile(args.path + "/src/walberla.h")) and (not args.force)):
+      raise RuntimeError(args.path + " is not the path to a waLBerla root directory! Specify -f to generate the files anyway.")
+
+   os.makedirs(args.path + "/src/mesa_pd/common", exist_ok = True)
+   os.makedirs(args.path + "/src/mesa_pd/data", exist_ok = True)
+   os.makedirs(args.path + "/src/mesa_pd/domain", exist_ok = True)
+   os.makedirs(args.path + "/src/mesa_pd/kernel", exist_ok = True)
+   os.makedirs(args.path + "/src/mesa_pd/mpi/notifications", exist_ok = True)
+   os.makedirs(args.path + "/src/mesa_pd/vtk", exist_ok = True)
+
+   shapes = ["Sphere", "HalfSpace", "CylindricalBoundary", "Box", "Ellipsoid"]
+
+   ps    = data.ParticleStorage()
+   ch    = data.ContactHistory()
+   lc    = data.LinkedCells()
+   slc   = data.SparseLinkedCells()
+   ss    = data.ShapeStorage(ps, shapes)
+   cs    = data.ContactStorage()
+
+   ps.addProperty("position",         "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("linearVelocity",   "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("invMass",          "walberla::real_t",        defValue="real_t(1)", syncMode="COPY")
+   ps.addProperty("force",            "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+   ps.addProperty("oldForce",         "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="MIGRATION")
+
+   ps.addProperty("shapeID",          "size_t",                  defValue="",          syncMode="COPY")
+   ps.addProperty("rotation",         "walberla::mesa_pd::Rot3", defValue="",          syncMode="ALWAYS")
+   ps.addProperty("angularVelocity",  "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("torque",           "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+   ps.addProperty("oldTorque",        "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="MIGRATION")
+
+   ps.addInclude("blockforest/BlockForest.h")
+   ps.addProperty("currentBlock",     "blockforest::Block*",     defValue="nullptr",   syncMode="NEVER")
+
+   ps.addProperty("type",             "uint_t",                  defValue="0",         syncMode="COPY")
+
+   ps.addProperty("flags",            "walberla::mesa_pd::data::particle_flags::FlagT", defValue="", syncMode="COPY")
+   ps.addProperty("nextParticle",     "int",                     defValue="-1",        syncMode="NEVER")
+
+   ps.addProperty("oldContactHistory", "std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>", defValue="", syncMode="ALWAYS")
+   ps.addProperty("newContactHistory", "std::map<walberla::id_t, walberla::mesa_pd::data::ContactHistory>", defValue="", syncMode="NEVER")
+
+   ps.addProperty("temperature",      "walberla::real_t",        defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("heatFlux",         "walberla::real_t",        defValue="real_t(0)", syncMode="NEVER")
+
+   # Properties for HCSITS
+   ps.addProperty("dv",               "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+   ps.addProperty("dw",               "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+
+   # for lbm_mesapd_coupling:
+   ps.addProperty("hydrodynamicForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+   ps.addProperty("hydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+   ps.addProperty("oldHydrodynamicForce", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+   ps.addProperty("oldHydrodynamicTorque", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="NEVER")
+
+   # for virtual mass:
+   ps.addProperty("virtualMass", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS") #TODO: fix syncModes -v
+   ps.addProperty("invMassIncludingVirtual", "walberla::real_t", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("oldAcceleration", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("oldWdot", "walberla::mesa_pd::Vec3", defValue="real_t(0)", syncMode="ALWAYS")
+
+   ps.addProperty("invInertiaBF", "walberla::mesa_pd::Mat3", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("virtualInertiaBF", "walberla::mesa_pd::Mat3", defValue="real_t(0)", syncMode="ALWAYS")
+   ps.addProperty("invInertiaBFIncludingVirtual", "walberla::mesa_pd::Mat3", defValue="real_t(0)", syncMode="ALWAYS")
+
+
+   ch.addProperty("tangentialSpringDisplacement", "walberla::mesa_pd::Vec3", defValue="real_t(0)")
+   ch.addProperty("isSticking",                   "bool",                    defValue="false")
+   ch.addProperty("impactVelocityMagnitude",      "real_t",                  defValue="real_t(0)")
+
+
+   cs.addProperty("id1",              "walberla::id_t",          defValue = "walberla::id_t(-1)", syncMode="NEVER")
+   cs.addProperty("id2",              "walberla::id_t",          defValue = "walberla::id_t(-1)", syncMode="NEVER")
+   cs.addProperty("distance",         "real_t",                  defValue = "real_t(1)",          syncMode="NEVER")
+   cs.addProperty("normal",           "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("position",         "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("t",                "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("o",                "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("r1",               "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("r2",               "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("mu",               "real_t",                  defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("p",                "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("diag_nto",         "walberla::mesa_pd::Mat3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("diag_nto_inv",     "walberla::mesa_pd::Mat3", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("diag_to_inv",      "walberla::mesa_pd::Mat2", defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("diag_n_inv",       "real_t",                  defValue = "real_t(0)",          syncMode="NEVER")
+   cs.addProperty("p",                "walberla::mesa_pd::Vec3", defValue = "real_t(0)",          syncMode="NEVER")
+
+
+   kernels = []
+   kernels.append( kernel.DetectAndStoreContacts() )
+   kernels.append( kernel.DoubleCast(shapes) )
+   kernels.append( kernel.ExplicitEuler() )
+   kernels.append( kernel.ExplicitEulerWithShape() )
+   kernels.append( kernel.ForceLJ() )
+   kernels.append( kernel.HCSITSRelaxationStep() )
+   kernels.append( kernel.HeatConduction() )
+   kernels.append( kernel.InitParticlesForHCSITS() )
+   kernels.append( kernel.InitContactsForHCSITS() )
+   kernels.append( kernel.IntegrateParticlesHCSITS() )
+   kernels.append( kernel.InsertParticleIntoLinkedCells() )
+   kernels.append( kernel.InsertParticleIntoSparseLinkedCells() )
+   kernels.append( kernel.LinearSpringDashpot() )
+   kernels.append( kernel.NonLinearSpringDashpot() )
+   kernels.append( kernel.SingleCast(shapes) )
+   kernels.append( kernel.SpringDashpot() )
+   kernels.append( kernel.TemperatureIntegration() )
+   kernels.append( kernel.VelocityVerlet() )
+   kernels.append( kernel.VelocityVerletWithShape() )
+   kernels.append( kernel.AddVirtualForceAndTorque() )
+   kernels.append( kernel.VirtualMass() )
+
+
+
+   ac = Accessor()
+   for k in kernels:
+      ac.mergeRequirements(k.getRequirements())
+   ac.printSummary()
+
+   comm = []
+   comm.append(mpi.BroadcastProperty())
+   comm.append(mpi.ClearNextNeighborSync())
+   comm.append(mpi.ReduceContactHistory())
+   comm.append(mpi.ReduceProperty())
+   comm.append(mpi.SyncGhostOwners(ps))
+   comm.append(mpi.SyncNextNeighbors(ps))
+
+
+   ps.generate(args.path + "/src/mesa_pd/")
+   ch.generate(args.path + "/src/mesa_pd/")
+   lc.generate(args.path + "/src/mesa_pd/")
+   slc.generate(args.path + "/src/mesa_pd/")
+   ss.generate(args.path + "/src/mesa_pd/")
+   cs.generate(args.path + "/src/mesa_pd/")
+
+   for k in kernels:
+      k.generate(args.path + "/src/mesa_pd/")
+
+   for c in comm:
+      c.generate(args.path + "/src/mesa_pd/")
diff --git a/python/mesa_pd/kernel/AddVirtualForceAndTorque.py b/python/mesa_pd/kernel/AddVirtualForceAndTorque.py
new file mode 100644
index 0000000000000000000000000000000000000000..16226658f46b41fc8b141b9cba0ba90b49179b67
--- /dev/null
+++ b/python/mesa_pd/kernel/AddVirtualForceAndTorque.py
@@ -0,0 +1,26 @@
+# -*- coding: utf-8 -*-
+
+from mesa_pd.accessor import Accessor
+from mesa_pd.utility import generateFile
+
+class AddVirtualForceAndTorque:
+   def __init__(self):
+      self.accessor = Accessor()
+      self.accessor.require("virtualMass", "walberla::real_t", access="g")
+      self.accessor.require("invMass", "walberla::real_t", access="g")
+      self.accessor.require("force", "walberla::mesa_pd::Vec3", access="gs")
+
+      self.accessor.require("oldWdot", "walberla::mesa_pd::Vec3", access="gs" )
+      self.accessor.require("oldAcceleration", "walberla::mesa_pd::Vec3", access="gs" )
+      self.accessor.require("rotation", "walberla::mesa_pd::Rot3", access="gs")
+      self.accessor.require("invInertiaBF", "walberla::mesa_pd::Mat3", access="g" )
+      self.accessor.require("torque", "walberla::mesa_pd::Vec3", access="gs" )
+      self.accessor.require("virtualInertiaBF", "walberla::mesa_pd::Mat3", access="g" )
+
+   def getRequirements(self):
+      return self.accessor
+
+   def generate(self, path):
+      context = dict()
+      context["interface"] = self.accessor.properties
+      generateFile(path, 'kernel/AddVirtualForceAndTorque.templ.h', context)
\ No newline at end of file
diff --git a/python/mesa_pd/kernel/VirtualMass.py b/python/mesa_pd/kernel/VirtualMass.py
new file mode 100644
index 0000000000000000000000000000000000000000..ac102bde1e4b1acccc08cdd20358566985621e81
--- /dev/null
+++ b/python/mesa_pd/kernel/VirtualMass.py
@@ -0,0 +1,25 @@
+# -*- coding: utf-8 -*-
+
+from mesa_pd.accessor import Accessor
+from mesa_pd.utility import generateFile
+
+class VirtualMass:
+   def __init__(self):
+      self.accessor = Accessor()
+      self.accessor.require("volume", "walberla::real_t", "g")
+
+      self.accessor.require("invMass", "walberla::real_t", "g")
+      self.accessor.require("virtualMass", "walberla::real_t", "s")
+      self.accessor.require("invMassIncludingVirtual", "walberla::real_t", "s")
+
+      self.accessor.require("invInertiaBF", "walberla::mesa_pd::Mat3", "g")
+      self.accessor.require("virtualInertiaBF", "walberla::mesa_pd::Mat3", "gs")
+      self.accessor.require("invInertiaBFIncludingVirtual", "walberla::mesa_pd::Mat3", "gs")
+
+   def getRequirements(self):
+      return self.accessor
+
+   def generate(self, path):
+      context = dict()
+      context["interface"] = self.accessor.properties
+      generateFile(path, 'kernel/VirtualMass.templ.h', context)
\ No newline at end of file
diff --git a/python/mesa_pd/kernel/__init__.py b/python/mesa_pd/kernel/__init__.py
index eef9907a547ca08e93d35a26eed6d89ef005580b..2972c4cdff5f91a3f04669ff249d980240aecb25 100644
--- a/python/mesa_pd/kernel/__init__.py
+++ b/python/mesa_pd/kernel/__init__.py
@@ -1,4 +1,5 @@
 # -*- coding: utf-8 -*-
+from .AddVirtualForceAndTorque import AddVirtualForceAndTorque
 from .DetectAndStoreContacts import DetectAndStoreContacts
 from .DoubleCast import DoubleCast
 from .ExplicitEuler import ExplicitEuler
@@ -18,8 +19,11 @@ from .SpringDashpot import SpringDashpot
 from .TemperatureIntegration import TemperatureIntegration
 from .VelocityVerlet import VelocityVerlet
 from .VelocityVerletWithShape import VelocityVerletWithShape
+from .VirtualMass import VirtualMass
 
-__all__ = ['DoubleCast',
+
+__all__ = ['AddVirtualForceAndTorque',
+           'DoubleCast',
            'DetectAndStoreContacts',
            'ExplicitEuler',
            'ExplicitEulerWithShape',
@@ -37,4 +41,5 @@ __all__ = ['DoubleCast',
            'SpringDashpot',
            'TemperatureIntegration',
            'VelocityVerlet',
-           'VelocityVerletWithShape']
+           'VelocityVerletWithShape',
+           'VirtualMass']
diff --git a/python/mesa_pd/templates/kernel/AddVirtualForceAndTorque.templ.h b/python/mesa_pd/templates/kernel/AddVirtualForceAndTorque.templ.h
new file mode 100644
index 0000000000000000000000000000000000000000..015dda606d28a1d80fc49369d5b8dce0272abd3a
--- /dev/null
+++ b/python/mesa_pd/templates/kernel/AddVirtualForceAndTorque.templ.h
@@ -0,0 +1,78 @@
+//======================================================================================================================
+//
+//  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 AddVirtualForce.h
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+#include "lbm_rpd_coupling/utility/ParticleSelectors.h"
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+/**
+ * Kernel which adds a virtual force to the particles.
+ *
+ * \remark
+ * \pre
+ * \post
+ * \ingroup mesa_pd_kernel
+ */
+
+struct AddVirtualForceAndTorque {
+   explicit AddVirtualForceAndTorque() {}
+
+   template <typename Accessor>
+   void operator()(const size_t i, Accessor& ac) const;
+};
+
+
+template <typename Accessor>
+inline void AddVirtualForceAndTorque::operator()(const size_t i,
+                                Accessor& ac) const {
+   static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+   lbm_rpd_coupling::LocalParticlesSelector localSelector;
+   if (!localSelector(i, ac)) {
+       // only set virtual mass on particles local on the process
+       return;
+   }
+
+   WALBERLA_ASSERT_FLOAT_UNEQUAL(ac.getVirtualMass(i), real_t(0.), "No virtual mass set on body. Was the virtualMass kernel not called before?");
+
+   const Vec3 virtualForce = ac.getVirtualMass(i)*ac.getOldAcceleration(i);
+   ac.setForce(i, ac.getForce(i)+virtualForce);
+
+   const Vec3 acceleration = real_t(1.)/(ac.getVirtualMass(i)+real_t(1.)/ac.getInvMass(i)) * ac.getForce(i);
+   ac.setOldAcceleration(i, acceleration);
+
+   if (!realIsEqual(ac.getInvInertiaBF(i).getDeterminant(), real_t(0))) {
+      const Vec3 wdot = math::transformMatrixRART(ac.getRotation(i).getMatrix(),
+                                                  ac.getInvInertiaBF(i)) * ac.getTorque(i);
+      const Vec3 virtualTorque = ac.getVirtualInertiaBF(i) * ac.getOldWdot(i);
+      ac.setTorque(i, ac.getTorque(i) + virtualTorque);
+      ac.setOldWdot(i, wdot);
+   }
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
diff --git a/python/mesa_pd/templates/kernel/VirtualMass.templ.h b/python/mesa_pd/templates/kernel/VirtualMass.templ.h
new file mode 100644
index 0000000000000000000000000000000000000000..cccce1c4f7a782e0fc2ceca1fbe204b182f5d55a
--- /dev/null
+++ b/python/mesa_pd/templates/kernel/VirtualMass.templ.h
@@ -0,0 +1,74 @@
+//======================================================================================================================
+//
+//  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 AddVirtualForce.h
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+
+/**
+ * Kernel which adds a virtual force to the particles.
+ *
+ * \remark
+ * \pre
+ * \post
+ * \ingroup mesa_pd_kernel
+ */
+
+struct VirtualMass {
+   explicit VirtualMass(const real_t C_v, const real_t C_v_omega, const real_t fluidDensity)
+      : C_v_(C_v), C_v_omega_(C_v_omega), fluidDensity_(fluidDensity) {}
+
+   template <typename Accessor>
+   void operator()(const size_t i, Accessor& ac) const;
+private:
+   real_t C_v_ = real_t(0.0);
+   real_t C_v_omega_ = real_t(0.0);
+   real_t fluidDensity_ = real_t(0.0);
+};
+
+
+template <typename Accessor>
+inline void VirtualMass::operator()(const size_t i,
+                                Accessor& ac) const {
+   static_assert(std::is_base_of<data::IAccessor, Accessor>::value, "please provide a valid accessor");
+
+   const real_t invMass = ac.getInvMass(i);
+
+   const real_t virtualMass = C_v_*fluidDensity_*ac.getVolume(i);
+   ac.setVirtualMass(i, virtualMass);
+   const real_t mass = real_t(1.)/invMass;
+   ac.setInvMassIncludingVirtual(i, real_t(1.) / (mass + virtualMass));
+
+   const Mat3 inertiaBF = ac.getInvInertiaBF(i).getInverse();
+   const real_t angularVirtualMass = C_v_omega_ * fluidDensity_ * ac.getVolume(i);
+   const Mat3 virtualInertiaBF = inertiaBF * invMass * angularVirtualMass; //TODO: can this be done universally?
+   const Mat3 inertiaBFIncludingVirtual = inertiaBF + virtualInertiaBF;
+   ac.setInvInertiaBFIncludingVirtual(i, inertiaBFIncludingVirtual.getInverse());
+   ac.setVirtualInertiaBF(i, virtualInertiaBF);
+}
+
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla