diff --git a/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h b/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h
new file mode 100644
index 0000000000000000000000000000000000000000..df369216e4d5e808ab9c5567897ae13ac7c35f85
--- /dev/null
+++ b/src/mesa_pd/kernel/cnt/AnisotropicVDWContact.h
@@ -0,0 +1,386 @@
+//======================================================================================================================
+//
+// 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
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <mesa_pd/common/ParticleFunctions.h>
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/IAccessor.h>
+
+#include <core/math/Angles.h>
+#include <core/math/Constants.h>
+#include <core/logging/Logging.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+namespace cnt {
+
+/**
+ * anisotropic vdW contact
+ */
+class AnisotropicVDWContact
+{
+public:
+ template<typename Accessor>
+ void operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor &ac);
+
+ static constexpr real_t eps_ = 0.03733_r;
+ static constexpr real_t A_ = 0.0223_r;
+ static constexpr real_t B_ = 1.31_r;
+ static constexpr real_t alf_ = 9.5_r;
+ static constexpr real_t bet_ = 4.0_r;
+ static constexpr real_t Cg_ = 90_r;
+ static constexpr real_t del_ = -7.5_r;
+
+ static constexpr size_t M = 5;
+ /// Set of fitted constants, that determine the level of "smoothness" of vdW potential -
+ /// magnitude of shear forces between sliding CNTs
+ static constexpr std::array<real_t, M> C = { 0.35819_r, 0.03263_r, -0.00138_r, -0.00017_r, 0.00024_r };
+
+ /// vdW interaction radius w/r to inertial segment radius
+ static constexpr real_t CutoffFactor = 4_r;
+ /// CNT radius
+ static constexpr real_t R_CNT = 6.78_r;
+ static constexpr real_t R_ = R_CNT;
+ static constexpr real_t r_cut_ = CutoffFactor * R_;
+
+ auto isParallel() const {return isParallel_;}
+ auto getLastEnergy() const {return energy_;}
+private:
+ real_t energy_; ///< total potential
+ bool isParallel_;
+};
+
+template<typename Accessor>
+inline
+void AnisotropicVDWContact::operator()(const size_t p_idx1,
+ const size_t p_idx2,
+ Accessor &ac)
+{
+ //===Adaptation of PFC5 vdW contact model implementation====
+
+ // Getting the orientations of segments
+ Vec3 b1 = ac.getRotation(p_idx1).getMatrix() * Vec3(1.0, 0.0, 0.0); ///< ball 1 axial direction
+ Vec3 b2 = ac.getRotation(p_idx2).getMatrix() * Vec3(1.0, 0.0, 0.0); ///< ball 2 axial direction
+
+
+ // Distance between segments
+
+ Vec3 n = ac.getPosition(p_idx2) - ac.getPosition(p_idx1); ///< contact normal
+ auto L = n.length();
+ n *= (1_r/L);
+
+
+ //WALBERLA_LOG_DEVEL( "Normal: n = " << n );
+ //WALBERLA_LOG_DEVEL( "Orientation of seg 2: b2 = " << b2 );
+ //WALBERLA_LOG_DEVEL( "Orientation of seg 1: b1 = " << b1 );
+ //WALBERLA_LOG_DEVEL( "Length of rad vect: L = " << L );
+
+
+ constexpr real_t TOL = 10e-8_r;
+ //---------------------
+ // NORMALS CALCULATION
+ //---------------------
+ // names convention:
+ // c1 - contact 1-2 normal
+ // b1 - ball 1 axial direction
+ // b2 - ball 2 axial direction
+ // b3 - neytral direction
+ // g - alighning torque direction
+ // d - neytral plane normal direction
+ // s - shear force direction
+
+ // angle gamma - angle between two axial directions
+ auto cos_gamma = b1 * b2;
+
+
+
+ // if the angle between two axal directions is blunt, then inverce b2
+ if (cos_gamma < 0_r)
+ {
+ b2 = -b2;
+ cos_gamma = -cos_gamma;
+ }
+ // check that cosine belongs [-1,1]
+ cos_gamma = std::min(1.0_r, cos_gamma);
+ cos_gamma = std::max(-1.0_r, cos_gamma);
+ isParallel_ = false;
+ if (L < 20_r && L > 16_r)
+ {
+ const auto gamma = acos(cos_gamma);
+ if (gamma < math::degToRad(10_r) || gamma > math::degToRad(170_r))
+ isParallel_ = true;
+ }
+ //WALBERLA_LOG_DEVEL( "cos_gamma: = " << cos_gamma );
+
+
+ // calculate functions of double argument
+ auto sin_gamma = std::sqrt(1.0_r - cos_gamma * cos_gamma);
+ auto cos_2gamma = cos_gamma * cos_gamma - sin_gamma * sin_gamma;
+ auto sin_2gamma = 2.0_r * sin_gamma * cos_gamma;
+
+ //WALBERLA_LOG_DEVEL( "sin_gamma: = " << sin_gamma );
+ //WALBERLA_LOG_DEVEL( "cos_2gamma: = " << cos_2gamma );
+ //WALBERLA_LOG_DEVEL( "sin_2gamma: = " << sin_2gamma );
+
+
+ // g - direction of the aligning torques - b1 X b2
+ Vec3 g(0.0, 0.0, 0.0);
+ if (sin_gamma > TOL)
+ {
+ g = b1 % b2;
+ g = g * (1.0_r / g.length());
+ }
+ //WALBERLA_LOG_DEVEL( "Aligning moment direction: g = " << g );
+
+ // b3 - vector defining the neutral plane ( plane of shear forces )
+ Vec3 b3 = b1 + b2;
+ b3 = b3 * (1.0_r / b3.length());
+ //WALBERLA_LOG_DEVEL( "Neutral plane defined by b3 = " << b3 );
+
+ // angle theta - angle between b3 and c1
+ auto cos_theta = b3 * n;
+ // check that cosine belongs [-1,1]
+ cos_theta = std::min(1.0_r, cos_theta);
+ cos_theta = std::max(-1.0_r, cos_theta);
+ //WALBERLA_LOG_DEVEL( "cos_theta: = " << cos_theta );
+
+ // calculation of shear force direction
+ Vec3 s(0.0, 0.0, 0.0);
+ Vec3 d(0.0, 0.0, 0.0);
+
+ if ((cos_theta > -1.0 + TOL) || (cos_theta < 1.0 - TOL))
+ d = n % b3;
+ s = n % d;
+ s = s * (1.0_r / s.length());
+
+ //WALBERLA_LOG_DEVEL( "Shear force direction: = " << s );
+ //--------------------------------
+ // NORMALS CALCULATION - END
+ //--------------------------------
+
+ // Fast calculation of trigonometric functions ( Chebyshev formulas )
+ real_t coss[M], sinn[M];
+ real_t sin_theta = std::sqrt(1.0_r - cos_theta * cos_theta);
+
+ coss[0] = cos_theta * cos_theta - sin_theta * sin_theta;
+ sinn[0] = 2.0_r * sin_theta * cos_theta;
+
+ for (size_t i = 0; i < M - 1; ++i)
+ {
+ coss[i + 1] = coss[i] * coss[0] - sinn[i] * sinn[0];
+ sinn[i + 1] = sinn[i] * coss[0] + sinn[0] * coss[i];
+ }
+
+ //WALBERLA_LOG_DEVEL( "coss: = " << coss[0] <<" "<< coss[1] <<" "<< coss[2] <<" "<< coss[3] <<" "<< coss[4]);
+ //WALBERLA_LOG_DEVEL( "sinn: = " << sinn[0] <<" "<< sinn[1] <<" "<< sinn[2] <<" "<< sinn[3] <<" "<< sinn[4]);
+
+ //WALBERLA_LOG_DEVEL( "C: = " << C[0] <<" "<< C[1] <<" "<< C[2] <<" "<< C[3] <<" "<< C[4]);
+
+ // Cutoff for theta adjustment
+ real_t W_th = 1_r;
+ real_t W_th_L = 0_r;
+ real_t W_th_LL = 0_r;
+
+ // Adjustment w/r to O
+ real_t TH = 1_r;
+ real_t TH_L = 0_r;
+ real_t TH_LL = 0_r;
+ real_t TH_O = 0_r;
+ real_t TH_OO = 0_r;
+
+ real_t sign = 1_r;
+ for (size_t i = 0; i < M; ++i)
+ {
+ TH = TH + W_th * C[i] * (sign + coss[i]);
+ TH_L = TH_L + W_th_L * C[i] * (sign + coss[i]);
+ TH_LL = TH_LL + W_th_LL * C[i] * (sign + coss[i]);
+ TH_O = TH_O - W_th * 2_r * real_t(i + 1) * C[i] * (sinn[i]);
+ TH_OO = TH_OO - W_th * 4_r * real_t(i + 1) * real_t(i + 1) * C[i] * (coss[i]);
+ sign *= -1_r;
+ }
+
+ //WALBERLA_LOG_DEVEL( "TH: = " << TH );
+ //WALBERLA_LOG_DEVEL( "TH_L: = " << TH_L );
+ //WALBERLA_LOG_DEVEL( "TH_LL: = " << TH_LL );
+ //WALBERLA_LOG_DEVEL( "TH_O: = " << TH_O );
+ //WALBERLA_LOG_DEVEL( "TH_OO: = " << TH_OO );
+
+
+ //------------------------------------------------------------------
+ // THIS BLOCK IMPLEMENTS IF THE DISTANCE L IS WITHIN WORKING RANGE
+ //------------------------------------------------------------------
+ if ((L < 2_r * r_cut_) && (L > 2_r * R_ * 1.2_r * TH))
+ {
+ //-----Constants that appear in the potential--------------------------
+ // This set of constants is described in our paper.
+ //---------------------------------------------------------------------
+
+ //-----Function D and its derivatives-----------------------
+ real_t D = L / (R_ * TH) - 2_r;
+ real_t D_L = 1_r / (R_ * TH) - (L * TH_L) / (R_ * TH * TH);
+ real_t D_O = -(L * TH_O) / (R_ * TH * TH);
+ real_t D_LL = -(TH_L) / (R_ * TH * TH);
+ D_LL = D_LL - ((R_ * TH * TH) * (TH_L + L * TH_LL)) / (R_ * R_ * TH * TH * TH * TH);
+ D_LL = D_LL + (2_r * R_ * L * TH * TH_L * TH_L) / (R_ * R_ * TH * TH * TH * TH);
+ real_t D_OO = -(R_ * L * TH_OO * TH * TH) / (R_ * R_ * TH * TH * TH * TH);
+ D_OO = D_OO + (2_r * R_ * L * TH * TH_O * TH_O) / (R_ * R_ * TH * TH * TH * TH);
+ //-----------------------------------------------------------
+
+ //WALBERLA_LOG_DEVEL( "D: = " << D );
+ //WALBERLA_LOG_DEVEL( "D_L: = " << D_L );
+ //WALBERLA_LOG_DEVEL( "D_LL: = " << D_LL );
+ //WALBERLA_LOG_DEVEL( "D_O: = " << D_O );
+ //WALBERLA_LOG_DEVEL( "D_OO: = " << D_OO );
+
+
+ //----Function Vc and its derivatives---------------------------------------
+ const real_t DpowAlpha0 = std::pow(D, -(alf_));
+ const real_t DpowAlpha1 = std::pow(D, -(alf_ + 1));
+ const real_t DpowAlpha2 = std::pow(D, -(alf_ + 2));
+ const real_t DpowBeta0 = std::pow(D, -(bet_));
+ const real_t DpowBeta1 = std::pow(D, -(bet_ + 1));
+ const real_t DpowBeta2 = std::pow(D, -(bet_ + 2));
+ real_t Vc = 4_r * eps_ * (A_ * DpowAlpha0 - B_ * DpowBeta0);
+ real_t Vc_L = 4_r * eps_ * (-alf_ * A_ * DpowAlpha1 + bet_ * B_ * DpowBeta1) * D_L;
+ real_t Vc_O = 4_r * eps_ * (-alf_ * A_ * DpowAlpha1 + bet_ * B_ * DpowBeta1) * D_O;
+ real_t Vc_LL = 4_r * eps_ * (alf_ * (alf_ + 1_r) * A_ * DpowAlpha2 - bet_ * (bet_ + 1_r) * B_ * DpowBeta2) * D_L;
+ Vc_LL = Vc_LL + 4_r * eps_ * (-alf_ * A_ * DpowAlpha1 + bet_ * B_ * DpowBeta1) * D_LL;
+ real_t Vc_OO = 4_r * eps_ * (alf_ * (alf_ + 1_r) * A_ * DpowAlpha2 - bet_ * (bet_ + 1_r) * B_ * DpowBeta2) * D_O;
+ Vc_OO = Vc_OO + 4_r * eps_ * (-alf_ * A_ * DpowAlpha1 + bet_ * B_ * DpowBeta1) * D_OO;
+ //--------------------------------------------------------------------------
+
+ //WALBERLA_LOG_DEVEL( "VC = " << Vc );
+ //WALBERLA_LOG_DEVEL( "VC_L = " << Vc_L );
+ //WALBERLA_LOG_DEVEL( "VC_LL = " << Vc_LL );
+
+ //WALBERLA_LOG_DEVEL( "VC_O = " << Vc_O );
+ //WALBERLA_LOG_DEVEL( "VC_OO = " << Vc_OO );
+
+
+ // Cutoff for u adjustment
+ real_t W_u = 1_r;
+ real_t W_u_L = 0_r;
+ real_t W_u_LL = 0_r;
+ WALBERLA_UNUSED(W_u_LL);
+
+ // Cubic cutoff function 3T->4T (hardcoded since we do not need to mess w these parameters)
+ constexpr auto Q1_ = -80.0_r;
+ constexpr auto Q2_ = 288.0_r;
+ constexpr auto Q3_ = -336.0_r;
+ constexpr auto Q4_ = 128.0_r;
+ const real_t rcut2inv = 1_r / (2.0_r * r_cut_);
+ real_t nd = L * rcut2inv;
+ if ((nd > 0.75_r) && (nd < 1.0_r))
+ {
+ W_u = Q1_ + Q2_ * nd + Q3_ * nd * nd + Q4_ * nd * nd * nd;
+ W_u_L = (Q2_ + 2.0_r * Q3_ * nd + 3.0_r * Q4_ * nd * nd) * rcut2inv;
+ W_u_LL = (2.0_r * Q3_ + 6.0_r * Q4_ * nd) * rcut2inv * rcut2inv;
+ }
+ //--------------------------------------------------------------------------
+
+ //WALBERLA_LOG_DEVEL( "W_u = " << W_u );
+ //WALBERLA_LOG_DEVEL( "W_u_L = " << W_u_L );
+ //WALBERLA_LOG_DEVEL( "W_u_LL = " << W_u_LL );
+
+
+ // Cutoff for gamma adjustment
+
+ real_t W_ga, W_ga_L, W_ga_LL;
+ if (L / R_ > 2.75_r)
+ {
+ W_ga = Cg_ * std::pow((L / R_), del_);
+ W_ga_L = ((del_ * Cg_) / R_) * std::pow((L / R_), del_ - 1_r);
+ W_ga_LL = ((del_ * (del_ - 1) * Cg_) / (R_ * R_)) * std::pow((L / R_), del_ - 2_r);
+ } else
+ {
+ W_ga = Cg_ * std::pow((2.75_r), del_);
+ W_ga_L = 0;
+ W_ga_LL = 0;
+ }
+
+ //WALBERLA_LOG_DEVEL( "W_ga = " << W_ga );
+ //WALBERLA_LOG_DEVEL( "W_ga_L = " << W_ga_L );
+ //WALBERLA_LOG_DEVEL( "W_ga_LL = " << W_ga_LL );
+
+
+ real_t GA = 1_r;
+ real_t GA_L = 0_r;
+ real_t GA_LL = 0_r;
+ WALBERLA_UNUSED(GA_LL);
+ real_t GA_G = 0_r;
+ real_t GA_GG = 0_r;
+ WALBERLA_UNUSED(GA_GG);
+ if (std::abs(sin_gamma) > TOL)
+ {
+ GA = 1_r + W_ga * (1_r - cos_2gamma);
+ GA_L = W_ga_L * (1_r - cos_2gamma);
+ GA_LL = W_ga_LL * (1_r - cos_2gamma);
+ GA_G = 2_r * W_ga * sin_2gamma;
+ GA_GG = 4_r * W_ga * cos_2gamma;
+ }
+
+ //----Forces and torque-----------------------
+ real_t FL = -GA_L * W_u * Vc - GA * W_u_L * Vc - GA * W_u * Vc_L;
+ real_t FO = -(1_r / L) * GA * W_u * Vc_O;
+ real_t MG = -GA_G * W_u * Vc;
+
+ //WALBERLA_LOG_DEVEL( "FL = " << FL );
+ //WALBERLA_LOG_DEVEL( "FO = " << FO );
+ //WALBERLA_LOG_DEVEL( "MG = " << MG );
+
+
+ Vec3 force = FL * n + FO * s;
+ Vec3 moment = MG * g;
+
+
+ //WALBERLA_LOG_DEVEL( "Contact force: = " << force );
+ //WALBERLA_LOG_DEVEL( "Contact moment: = " << moment );
+
+ addForceAtomic(p_idx1, ac, -force);
+ addForceAtomic(p_idx2, ac, force);
+
+ addTorqueAtomic(p_idx1, ac, -moment);
+ addTorqueAtomic(p_idx2, ac, moment);
+
+ // Potential energy
+ energy_ = GA * W_u * Vc;
+ // WALBERLA_LOG_DEVEL( "U_vdw = " << U );
+ } else if (L <= 2_r * R_ * 1.2_r * TH)
+ { // Small distance
+ //WALBERLA_LOG_DEVEL( "Small distance");
+ energy_ = 0_r;
+ real_t F = -1_r;
+ addForceAtomic(p_idx1, ac, F * n);
+ addForceAtomic(p_idx2, ac, -F * n);
+ }
+}
+
+} //namespace cnt
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file
diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt
index 4e7985268047a837470ea5ed0b2a4793b1785873..48c07d87c940f09f35235a6f1d082bb9e2fc1fea 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -84,6 +84,9 @@ waLBerla_execute_test( NAME MESA_PD_DropTestGeneral )
waLBerla_compile_test( NAME MESA_PD_Kernel_ClearNextNeighborSync FILES kernel/ClearNextNeighborSync.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_ClearNextNeighborSync PROCESSES 2 )
+waLBerla_compile_test( NAME MESA_PD_Kernel_CNT_AnisotropicVDWContact FILES kernel/cnt/AnisotropicVDWContact.test.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Kernel_CNT_AnisotropicVDWContact )
+
waLBerla_compile_test( NAME MESA_PD_Kernel_CNT_IntegratedVDWContact FILES kernel/cnt/IntegratedVDWContact.test.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_CNT_IntegratedVDWContact )
diff --git a/tests/mesa_pd/kernel/cnt/AnisotropicVDWContact.test.cpp b/tests/mesa_pd/kernel/cnt/AnisotropicVDWContact.test.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9fcde4d9a8389a5858e5cb99d9a193739a6b2d1d
--- /dev/null
+++ b/tests/mesa_pd/kernel/cnt/AnisotropicVDWContact.test.cpp
@@ -0,0 +1,135 @@
+//======================================================================================================================
+//
+// 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
+//! \author Igor Ostanin <i.ostanin@skoltech.ru>
+//! \author Grigorii Drozdov <drozd013@umn.edu>
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "SphericalSegmentAccessor.h"
+
+#include "mesa_pd/data/Flags.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/cnt/AnisotropicVDWContact.h"
+#include "mesa_pd/kernel/cnt/ViscousDamping.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+#include "mesa_pd/vtk/ParticleVtkOutput.h"
+
+#include "core/Environment.h"
+#include "core/math/Constants.h"
+#include "vtk/VTKOutput.h"
+
+namespace walberla {
+using namespace walberla::mesa_pd;
+
+int main(int argc, char **argv)
+{
+ Environment env(argc, argv);
+ walberla::mpi::MPIManager::instance()->useWorldComm();
+
+ if (std::is_same<walberla::real_t, float>::value)
+ {
+ WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
+ return EXIT_SUCCESS;
+ }
+
+ logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
+
+ WALBERLA_LOG_INFO_ON_ROOT("loading configuration parameters");
+ constexpr auto numSimulationSteps = 20000ll;
+ constexpr auto outputInterval = 100ll;
+
+ WALBERLA_LOG_INFO_ON_ROOT("creating initial particle setup");
+ auto ps = std::make_shared<data::ParticleStorage>(10);
+ auto ac = SphericalSegmentAccessor(ps);
+
+ using namespace kernel::cnt;
+ data::Particle &&sp1 = *ps->create();
+ sp1.setPosition(Vec3(0_r, 0_r, 0_r));
+ sp1.setSegmentID(1);
+ sp1.setClusterID(1);
+
+ data::Particle &&sp2 = *ps->create();
+ sp2.setPosition(Vec3(20_r, 20_r, 20_r));
+ sp2.setSegmentID(2);
+ sp2.setClusterID(2);
+
+ WALBERLA_LOG_INFO_ON_ROOT("setting up VTK output");
+ auto vtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(ps);
+ vtkOutput->addOutput<data::SelectParticlePosition>("position");
+ auto vtkWriter = walberla::vtk::createVTKOutput_PointData(vtkOutput,
+ "cnt",
+ 1,
+ "vtk_integrated",
+ "particles",
+ false,
+ false);
+
+ WALBERLA_LOG_INFO_ON_ROOT("setting up interaction models");
+ kernel::cnt::AnisotropicVDWContact vdW_anisotropic;
+ kernel::cnt::ViscousDamping viscous_damping(0.1_r * 1052.0_r, 0.1_r * 1052.0_r);
+ kernel::VelocityVerletPreForceUpdate vv_pre(kernel::cnt::dT);
+ kernel::VelocityVerletPostForceUpdate vv_post(kernel::cnt::dT);
+
+ WALBERLA_LOG_INFO_ON_ROOT("running simulation");
+
+ real_t U = 0_r;
+ for (auto i = 0; i < numSimulationSteps; ++i)
+ {
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ vv_pre,
+ ac);
+
+ U = 0_r;
+ ps->forEachParticlePairHalf(false,
+ kernel::SelectAll(),
+ ac,
+ [&](size_t p_idx1, size_t p_idx2)
+ {
+ vdW_anisotropic(p_idx1, p_idx2, ac);
+ U += vdW_anisotropic.getLastEnergy();
+ viscous_damping(p_idx1, p_idx2, ac);
+ });
+
+ ps->forEachParticle(false,
+ kernel::SelectAll(),
+ ac,
+ vv_post,
+ ac);
+
+ if( i % outputInterval == 0 )
+ {
+// vtkWriter->write();
+// WALBERLA_LOG_DEVEL(i << " : " << U);
+ }
+ }
+
+ WALBERLA_CHECK_LESS(U, -1.16_r)
+
+ return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main(int argc, char *argv[])
+{
+ return walberla::main(argc, argv);
+}