Added lbm-mesa_pd coupling to have a fluid-particle coupling to the new rigid...

Added lbm-mesa_pd coupling to have a fluid-particle coupling to the new rigid particle dynamics simulations

apps/benchmarks/CMakeLists.txt

@@ -4,6 +4,7 @@ add_subdirectory( DEM )
 add_subdirectory( FieldCommunication )
 add_subdirectory( MeshDistance )
 add_subdirectory( CouetteFlow )
+add_subdirectory( FluidParticleCoupling )
 add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
 add_subdirectory( GranularGas )
 add_subdirectory( LennardJones )

apps/benchmarks/FluidParticleCoupling/CMakeLists.txt

+waLBerla_link_files_to_builddir( "*.dat" )
+
+waLBerla_python_file_generates(GeneratedLBM.py
+      GeneratedLBM.cpp GeneratedLBM.h
+      )
+
+waLBerla_python_file_generates(GeneratedLBMWithForce.py
+      GeneratedLBMWithForce.cpp GeneratedLBMWithForce.h
+      )
+
+waLBerla_add_executable ( NAME SphereWallCollision FILES SphereWallCollision.cpp GeneratedLBM.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME SettlingSphereTenCate FILES SettlingSphereInBox.cpp GeneratedLBM.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME SphereMovingWithPrescribedVelocity FILES SphereMovingWithPrescribedVelocity.cpp GeneratedLBM.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm mesa_pd lbm_mesapd_coupling postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME LubricationForceEvaluation FILES LubricationForceEvaluation.cpp GeneratedLBM.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME DragForceSphere FILES DragForceSphere.cpp GeneratedLBMWithForce.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME ForcesOnSphereNearPlane FILES ForcesOnSphereNearPlane.cpp GeneratedLBM.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
+
+waLBerla_add_executable ( NAME ObliqueDryCollision FILES ObliqueDryCollision.cpp
+                          DEPENDS blockforest core mesa_pd postprocessing )
+
+waLBerla_add_executable ( NAME ObliqueWetCollision FILES ObliqueWetCollision.cpp GeneratedLBM.py
+                          DEPENDS blockforest boundary core domain_decomposition field lbm lbm_mesapd_coupling mesa_pd postprocessing timeloop vtk )
\ No newline at end of file

apps/benchmarks/FluidParticleCoupling/GeneratedLBM.py

+import sympy as sp
+import pystencils as ps
+from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_ast, create_lb_method
+from lbmpy_walberla import generate_lattice_model
+from pystencils_walberla import CodeGeneration
+
+from lbmpy.creationfunctions import create_lb_collision_rule
+from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
+from lbmpy.stencils import get_stencil
+from lbmpy.maxwellian_equilibrium import get_moments_of_discrete_maxwellian_equilibrium
+from lbmpy.forcemodels import *
+
+from collections import OrderedDict
+from lbmpy.methods import create_with_discrete_maxwellian_eq_moments
+
+with CodeGeneration() as ctx:
+
+    # methods: TRTlike, KBC, SRT, cumulant
+    generatedMethod = "TRTlike"
+
+    print("Generating " + generatedMethod + " LBM method")
+
+    #
+    # x,y,z = MOMENT_SYMBOLS
+    #
+    # e_sq = x**2 + y**2 + z**2
+    #
+    # moments = [0] * 19
+    #
+    # moments[0] = sp.sympify(1)
+    # moments[1] = 19 * e_sq - 30
+    # moments[2] = (21 * e_sq**2 - 53 * e_sq + 24) / 2
+    #
+    # moments[3] = x
+    # moments[5] = y
+    # moments[7] = z
+    #
+    # moments[4] = (5*e_sq - 9) * x
+    # moments[6] = (5*e_sq - 9) * y
+    # moments[8] = (5*e_sq - 9) * z
+    #
+    # moments[9] = 3 * x**2 - e_sq
+    # moments[11] = y**2 - z**2
+    #
+    # moments[13] = x * y
+    # moments[14] = y * z
+    # moments[15] = x * z
+    #
+    # moments[10] = (3* e_sq - 5)*(3*x**2 - e_sq)
+    # moments[12] = (3* e_sq - 5)*(y**2 - z**2)
+    # moments[16] = (y**2 - z**2) * x
+    # moments[17] = (z**2 - x**2) * y
+    # moments[18] = (x**2 - y**2) * z
+    #
+    # m_eq = get_moments_of_discrete_maxwellian_equilibrium(stencil, moments, order=2, c_s_sq=sp.Rational(1,3))
+    #
+    # rr_dict = OrderedDict(zip(moments, omega))
+    #
+    # forcing = sp.symbols("forcing_:%d" % 3)
+    # forcing=(sp.symbols("fx"),0,0)
+    # forcemodel=Luo(forcing) #None
+    # method = create_with_discrete_maxwellian_eq_moments(stencil, rr_dict, compressible=False, force_model=forcemodel)
+    # at this stage, we have the MRT model from the LBM book!
+
+     #this is the schiller / walberla MRT
+
+    if generatedMethod == "TRTlike":
+         stencil = get_stencil("D3Q19", 'walberla')
+         omega = sp.symbols("omega_:%d" % len(stencil))
+         method = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False)
+
+         def modification_func(moment, eq, rate):
+             omegaVisc = sp.Symbol("omega_visc")
+             omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# sp.Symbol("omega_bulk")
+             omegaMagic = sp.Symbol("omega_magic")
+             if get_order(moment) <= 1:
+                 return moment, eq, 0
+             elif rate == omega[1]:
+                 return moment, eq, omegaBulk.center_vector
+             elif rate == omega[2]:
+                 return moment, eq, omegaBulk.center_vector
+             elif is_even(moment):
+                 return moment, eq, omegaVisc
+             else:
+                 return moment, eq, omegaMagic
+
+         my_method = create_lb_method_from_existing(method, modification_func)
+
+         # optimizations
+
+         collision_rule = create_lb_collision_rule(lb_method=my_method,
+                             optimization={"cse_global": True,
+                                           "cse_pdfs": False,
+                                           #"split": True
+                                           } )
+
+         generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
+
+    elif generatedMethod == "KBC":
+        methodName = 'trt-kbc-n4'
+        #omega_field = ps.fields("omegaField(1): [3D]", layout='fzyx')  omega_output_field=omega_field.center,
+        collision_rule = create_lb_collision_rule(method=methodName,entropic=True, stencil=get_stencil("D3Q27", 'walberla'), compressible=True,
+                        optimization={"cse_global": True,
+                                      "cse_pdfs": False,
+                                      "split": True})
+        generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
+
+    elif generatedMethod == "SRT":
+        methodName = 'srt'
+        collision_rule = create_lb_collision_rule(method=methodName,stencil=get_stencil("D3Q19", 'walberla'),
+                     optimization={"cse_global": True,
+                                   "cse_pdfs": False,
+                                   "split": True})
+        generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
+
+    elif generatedMethod == "cumulant":
+
+        x,y,z = MOMENT_SYMBOLS
+
+        cumulants = [0] * 27
+
+        cumulants[0] = sp.sympify(1) #000
+
+        cumulants[1] = x #100
+        cumulants[2] = y #010
+        cumulants[3] = z #001
+
+        cumulants[4] = x*y #110
+        cumulants[5] = x*z #101
+        cumulants[6] = y*z #011
+
+        cumulants[7] = x**2 - y**2 #200 - 020
+        cumulants[8] = x**2 - z**2 #200 - 002
+        cumulants[9] = x**2 + y**2 + z**2 #200 + 020 + 002
+
+        cumulants[10] = x*y**2 + x*z**2 #120 + 102
+        cumulants[11] = x**2*y + y*z**2 #210 + 012
+        cumulants[12] = x**2*z + y**2*z #201 + 021
+        cumulants[13] = x*y**2 - x*z**2 #120 - 102
+        cumulants[14] = x**2*y - y*z**2 #210 - 012
+        cumulants[15] = x**2*z - y**2*z #201 - 021
+
+        cumulants[16] = x*y*z #111
+
+        cumulants[17] = x**2*y**2 - 2*x**2*z**2 + y**2*z**2 # 220- 2*202 +022
+        cumulants[18] = x**2*y**2 + x**2*z**2 - 2*y**2*z**2 # 220 + 202 - 2*002
+        cumulants[19] = x**2*y**2 + x**2*z**2 + y**2*z**2 # 220 + 202 + 022
+
+        cumulants[20] = x**2*y*z # 211
+        cumulants[21] = x*y**2*z # 121
+        cumulants[22] = x*y*z**2 # 112
+
+        cumulants[23] = x**2*y**2*z # 221
+        cumulants[24] = x**2*y*z**2 # 212
+        cumulants[25] = x*y**2*z**2 # 122
+
+        cumulants[26] = x**2*y**2*z**2 # 222
+
+
+        from lbmpy.moments import get_order
+        def get_relaxation_rate(cumulant, omega):
+            if get_order(cumulant) <= 1:
+                return 0
+            elif get_order(cumulant) == 2 and cumulant != x**2 + y**2 + z**2:
+                return omega
+            else:
+                return 1
+
+        stencil = get_stencil("D3Q27", 'walberla')
+
+        omega = sp.Symbol("omega")
+        rr_dict = OrderedDict((c, get_relaxation_rate(c, omega))
+                              for c in cumulants)
+
+        from lbmpy.methods import create_with_continuous_maxwellian_eq_moments
+        my_method = create_with_continuous_maxwellian_eq_moments(stencil, rr_dict, cumulant=True, compressible=True)
+
+        collision_rule = create_lb_collision_rule(lb_method=my_method,
+                            optimization={"cse_global": True,
+                                          "cse_pdfs": False})
+
+        generate_lattice_model(ctx, 'GeneratedLBM', collision_rule)
+
+    else:
+        print("Invalid generated method string! " + generatedMethod)
+
+
+

apps/benchmarks/FluidParticleCoupling/GeneratedLBMWithForce.py

+import sympy as sp
+import pystencils as ps
+from lbmpy.creationfunctions import create_lb_method_from_existing, create_lb_ast, create_lb_method
+from lbmpy_walberla import generate_lattice_model
+from pystencils_walberla import CodeGeneration
+
+from lbmpy.creationfunctions import create_lb_collision_rule
+from lbmpy.moments import MOMENT_SYMBOLS, is_even, get_order
+from lbmpy.stencils import get_stencil
+from lbmpy.maxwellian_equilibrium import get_moments_of_discrete_maxwellian_equilibrium
+from lbmpy.forcemodels import *
+
+from collections import OrderedDict
+from lbmpy.methods import create_with_discrete_maxwellian_eq_moments
+
+with CodeGeneration() as ctx:
+
+    forcing=(sp.symbols("fx"),0,0)
+    forcemodel=Luo(forcing) #None
+
+    # methods: TRTlike, KBC, SRT
+    generatedMethod = "TRTlike"
+
+    if generatedMethod == "TRTlike":
+        stencil = get_stencil("D3Q19", 'walberla')
+        omega = sp.symbols("omega_:%d" % len(stencil))
+
+        methodWithForce = create_lb_method(stencil=stencil, method='mrt', maxwellian_moments=False,force_model=forcemodel)
+
+        def modification_func(moment, eq, rate):
+            omegaVisc = sp.Symbol("omega_visc")
+            omegaBulk = ps.fields("omega_bulk: [3D]", layout='fzyx')# = sp.Symbol("omega_bulk")
+            #omegaBulk = sp.Symbol("omega_bulk")
+            omegaMagic = sp.Symbol("omega_magic")
+            if get_order(moment) <= 1:
+                return moment, eq, 0
+            elif rate == omega[1]:
+                return moment, eq, omegaBulk.center_vector
+            elif rate == omega[2]:
+                return moment, eq, omegaBulk.center_vector
+            elif is_even(moment):
+                return moment, eq, omegaVisc
+            else:
+                return moment, eq, omegaMagic
+
+        my_methodWithForce = create_lb_method_from_existing(methodWithForce, modification_func)
+
+        collision_rule = create_lb_collision_rule(lb_method=my_methodWithForce)
+        #    ,
+        #                     optimization={"cse_global": True,
+        #                                   "cse_pdfs": False,
+        #                                   "split": True,
+        #                                   "vectorization":{'instruction_set': 'avx',
+        #                                                    'nontemporal': True,
+        #                                                    'assume_inner_stride_one': True}})
+
+        generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
+
+    elif generatedMethod == "SRT":
+        collision_rule = create_lb_collision_rule(method='srt',stencil=get_stencil("D3Q19", 'walberla'), force_model=forcemodel,
+                            optimization={"cse_global": True,
+                                          "cse_pdfs": False,
+                                          #"split": True
+                                          })
+        generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
+
+    elif generatedMethod == "TRT":
+        collision_rule = create_lb_collision_rule(method='trt',stencil=get_stencil("D3Q19", 'walberla'), force_model=forcemodel, maxwellian_moments=False)
+        generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
+
+
+    elif generatedMethod == "KBC":
+        collision_rule = create_lb_collision_rule(method='trt-kbc-n4',entropic=True, stencil=get_stencil("D3Q27", 'walberla'), compressible=True, force_model=forcemodel,
+                            optimization={"cse_global": True,
+                                          "cse_pdfs": False,
+                                          #"split": True
+                                          })
+        generate_lattice_model(ctx, 'GeneratedLBMWithForce', collision_rule)
+
+    else:
+        print("Invalid generated method string! " + generatedMethod)

apps/benchmarks/FluidParticleCoupling/ObliqueDryCollision.cpp

+//======================================================================================================================
+//
+//  This file is part of waLBerla. waLBerla is free software: you can
+//  redistribute it and/or modify it under the terms of the GNU General Public
+//  License as published by the Free Software Foundation, either version 3 of
+//  the License, or (at your option) any later version.
+//
+//  waLBerla is distributed in the hope that it will be useful, but WITHOUT
+//  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+//  FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
+//  for more details.
+//
+//  You should have received a copy of the GNU General Public License along
+//  with waLBerla (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file   ObliqueDryCollision.cpp
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "mesa_pd/collision_detection/AnalyticContactDetection.h"
+#include "mesa_pd/common/ParticleFunctions.h"
+
+#include "mesa_pd/data/ParticleAccessorWithShape.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/data/ShapeStorage.h"
+
+#include "mesa_pd/kernel/DoubleCast.h"
+#include "mesa_pd/kernel/ExplicitEulerWithShape.h"
+#include "mesa_pd/kernel/VelocityVerletWithShape.h"
+#include "mesa_pd/kernel/LinearSpringDashpot.h"
+#include "mesa_pd/mpi/ReduceContactHistory.h"
+
+#include "core/Environment.h"
+#include "core/logging/Logging.h"
+
+#include <iostream>
+
+namespace walberla {
+namespace mesa_pd {
+
+/*
+ * Simulates oblique sphere-wall collision and checks rebound angle, i.e. the tangential part of the collision model.
+ *
+ */
+int main( int argc, char ** argv )
+{
+   walberla::mpi::Environment env(argc, argv);
+   WALBERLA_UNUSED(env);
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   real_t uNin_SI = real_t(1.7); // m/s
+   real_t diameter_SI = real_t(0.00318); // m
+   //real_t density_SI = real_t(2500); // kg/m**3, not used
+
+   real_t uNin       = real_t(0.02);
+   real_t diameter   = real_t(20);
+   real_t radius     = real_t(0.5) * diameter;
+   real_t density    = real_c(2.5);
+
+   // these values have actually no influence here and are just computed for completeness
+   real_t dx_SI = diameter_SI / diameter;
+   real_t dt_SI = uNin / uNin_SI * dx_SI;
+
+   real_t impactAngle       = real_t(0);
+   real_t dt                = real_t(0.1); // = (1 / #sub steps)
+   real_t frictionCoeff_s   = real_t(0.8); // no influence
+   real_t frictionCoeff_d   = real_t(0.12); // paper: 0.125+-0.007
+   std::string filename     = "TangentialCollision.txt";
+   real_t collisionTime     = real_t(80);
+   real_t nu                = real_t(0.22); //Poissons ratio
+   bool useVelocityVerlet   = false;
+   real_t restitutionCoeff  = real_t(0.83);
+
+   for( int i = 1; i < argc; ++i )
+   {
+      if( std::strcmp( argv[i], "--impactAngle" ) == 0 ) { impactAngle = real_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--dt" ) == 0 ) { dt = real_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--staticFriction" ) == 0 ) { frictionCoeff_s = real_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--dynamicFriction" ) == 0 ) { frictionCoeff_d = real_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--collisionTime" ) == 0 ) { collisionTime = real_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--nu" ) == 0 ) { nu = real_c( std::atof( argv[++i] ) ); continue; }
+      if( std::strcmp( argv[i], "--filename" ) == 0 ) { filename = argv[++i]; continue; }
+      if( std::strcmp( argv[i], "--useVV" )   == 0 ) { useVelocityVerlet = true; continue; }
+      if( std::strcmp( argv[i], "--coefficientOfRestitution" ) == 0 ) { restitutionCoeff = real_c( std::atof( argv[++i] ) ); continue; }
+   }
+
+   WALBERLA_LOG_INFO_ON_ROOT("******************************************************");
+   WALBERLA_LOG_INFO_ON_ROOT("**                  NEW SIMULATION                  **");
+   WALBERLA_LOG_INFO_ON_ROOT("******************************************************");
+   WALBERLA_LOG_INFO_ON_ROOT("dt_SI = " << dt_SI);
+   WALBERLA_LOG_INFO_ON_ROOT("dx_SI = " << dx_SI);
+   WALBERLA_LOG_INFO_ON_ROOT("impactAngle = " << impactAngle);
+   WALBERLA_LOG_INFO_ON_ROOT("dt = " << dt);
+   WALBERLA_LOG_INFO_ON_ROOT("frictionCoeff_s = " << frictionCoeff_s);
+   WALBERLA_LOG_INFO_ON_ROOT("frictionCoeff_d = " << frictionCoeff_d);
+   WALBERLA_LOG_INFO_ON_ROOT("collisionTime = " << collisionTime);
+   WALBERLA_LOG_INFO_ON_ROOT("nu = " << nu);
+   WALBERLA_LOG_INFO_ON_ROOT("integrator = " << (useVelocityVerlet ? "Velocity Verlet" : "Explicit Euler"));
+   WALBERLA_LOG_INFO_ON_ROOT("restitutionCoeff = " << restitutionCoeff);
+
+
+   //init data structures
+   auto ps = walberla::make_shared<data::ParticleStorage>(2);
+   auto ss = walberla::make_shared<data::ShapeStorage>();
+   using ParticleAccessor_T = mesa_pd::data::ParticleAccessorWithShape;
+   auto accessor = walberla::make_shared<ParticleAccessor_T >(ps, ss);
+
+   auto sphereShape = ss->create<data::Sphere>( radius );
+   ss->shapes[sphereShape]->updateMassAndInertia(density);
+
+   const real_t particleMass =  real_t(1) / ss->shapes[sphereShape]->getInvMass();
+   const real_t Mij = particleMass; // * particleMass / ( real_t(2) * particleMass ); // Mij = M for sphere-wall collision
+   const real_t lnDryResCoeff = std::log(restitutionCoeff);
+
+   // normal material parameters
+   const real_t stiffnessN = Mij * ( math::pi * math::pi + lnDryResCoeff * lnDryResCoeff ) / (collisionTime * collisionTime); // Costa et al., Eq. 18
+   const real_t dampingN = - real_t(2) * Mij * lnDryResCoeff / collisionTime; // Costa et al., Eq. 18
+
+   WALBERLA_LOG_INFO_ON_ROOT("normal: stiffness = " << stiffnessN << ", damping = " << dampingN);
+
+   // from Thornton et al
+   const real_t kappa = real_t(2) * ( real_t(1) - nu ) / ( real_t(2) - nu ) ;
+   const real_t stiffnessT = kappa * stiffnessN;
+   const real_t dampingT = std::sqrt(kappa) * dampingN;
+
+   WALBERLA_LOG_INFO_ON_ROOT("tangential: kappa = " << kappa << ", stiffness T = " << stiffnessT << ", damping T = " << dampingT);
+
+   real_t uTin = uNin * impactAngle;
+
+   // create sphere
+   data::Particle&& p = *ps->create();
+   p.setPosition(Vec3(0,0,2*radius));
+   p.setLinearVelocity(Vec3(uTin, 0., -uNin));
+   p.setType(0);
+
+   // create plane
+   data::Particle&& p0 = *ps->create(true);
+   p0.setPosition(Vec3(0,0,0));
+   p0.setShapeID(ss->create<data::HalfSpace>(Vector3<real_t>(0,0,1)));
+   p0.setType(0);
+   data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::INFINITE);
+   data::particle_flags::set(p0.getFlagsRef(), data::particle_flags::FIXED);
+
+   // velocity verlet
+   kernel::VelocityVerletWithShapePreForceUpdate  vvPreForce( dt );
+   kernel::VelocityVerletWithShapePostForceUpdate vvPostForce( dt );
+
+   // explicit euler
+   kernel::ExplicitEulerWithShape explEuler( dt );
+
+   // collision response
+   collision_detection::AnalyticContactDetection acd;
+   kernel::DoubleCast           double_cast;
+   kernel::LinearSpringDashpot  dem(1);
+   mpi::ReduceContactHistory    rch;
+   dem.setStiffnessN(0,0,stiffnessN);
+   dem.setStiffnessT(0,0,stiffnessT);
+   dem.setDampingN(0,0,dampingN);
+   dem.setDampingT(0,0,dampingT);
+   dem.setFrictionCoefficientStatic(0,0,frictionCoeff_s);
+   dem.setFrictionCoefficientDynamic(0,0,frictionCoeff_d);
+
+   WALBERLA_LOG_DEVEL("begin: vel = " << p.getLinearVelocity() << ", contact vel: " << getVelocityAtWFPoint(0,*accessor,p.getPosition() + Vec3(0,0,-radius)) );
+
+   uint_t steps = 0;
+   real_t maxPenetration = real_t(0);
+   do
+   {
+      if(useVelocityVerlet) vvPreForce(0,*accessor);
+
+      real_t penetration;
+
+      if (double_cast(0, 1, *accessor, acd, *accessor ))
+      {
+         penetration = acd.getPenetrationDepth();
+         maxPenetration = std::max( maxPenetration, std::abs(penetration));
+
+         dem(acd.getIdx1(), acd.getIdx2(), *accessor, acd.getContactPoint(), acd.getContactNormal(), acd.getPenetrationDepth(), dt);
+         auto force = accessor->getForce(0);
+         auto torque = accessor->getTorque(0);
+         WALBERLA_LOG_INFO(steps << ": penetration = " << penetration << " || vel = " << accessor->getLinearVelocity(0) << " || force = " << force << ", torque = " << torque);
+      }
+      rch(*ps);
+
+      if(useVelocityVerlet) vvPostForce(0,*accessor);
+      else explEuler(0, *accessor);
+
+      ++steps;
+   } while (double_cast(0, 1, *accessor, acd, *accessor ) || p.getLinearVelocity()[2] < 0);
+
+   real_t uTout = p.getLinearVelocity()[0] - radius * p.getAngularVelocity()[1];
+   WALBERLA_LOG_DEVEL("end: linear vel = " << p.getLinearVelocity() << ", angular vel = " << p.getAngularVelocity());
+
+   real_t reboundAngle = uTout / uNin;
+   WALBERLA_LOG_INFO_ON_ROOT("gamma_in = " << impactAngle);
+   WALBERLA_LOG_INFO_ON_ROOT("gamma_out = " << reboundAngle);
+
+   WALBERLA_LOG_INFO_ON_ROOT("Thornton: sliding should occur if " << real_t(2) * impactAngle / ( frictionCoeff_d * ( real_t(1) + restitutionCoeff)) << " >= " << real_t(7) - real_t(1) / kappa );
+   WALBERLA_LOG_INFO_ON_ROOT("Max penetration = " << maxPenetration << " -> " << maxPenetration / radius * 100. << "% of radius");
+
+   std::ofstream file;
+   file.open( filename.c_str(), std::ios::out | std::ios::app );
+   file << impactAngle << " " << reboundAngle << "\n";
+   file.close();
+
+   return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+   return walberla::mesa_pd::main(argc, argv);
+}

apps/benchmarks/FluidParticleCoupling/inputObliqueWetCollision.dat

+Setup
+{
+
+   case 3;
+
+   impactRatio 0;
+
+   variant Acceleration; // Gravity or Acceleration
+
+   // control
+   applyOutflowBCAtTop true;
+
+   // LBM
+   diameter 30;
+   magicNumber 0.1875;
+   bulkViscRateFactor 100.;
+   useOmegaBulkAdaption true;
+
+   // RPD
+   numRPDSubCycles 10;
+   useLubricationCorrection true;
+   lubricationCutOffDistanceNormal 0.666667;
+   lubricationCutOffDistanceTangentialTranslational 0.5;
+   lubricationCutOffDistanceTangentialRotational 0.5;
+   lubricationMinimalGapSizeNonDim 0.002;
+
+   useVelocityVerlet true;
+