From 88bb098cf3a089df712733131d506107f52f01fc Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Thu, 17 Dec 2020 15:40:02 +0100
Subject: [PATCH] [ADD] MESA-PD kernel for the VBondModel wall contact
---
src/mesa_pd/kernel/VBondModel/WallContact.h | 113 ++++++++++++++++++
tests/mesa_pd/CMakeLists.txt | 3 +
.../kernel/VBondModel/WallContact.test.cpp | 102 ++++++++++++++++
3 files changed, 218 insertions(+)
create mode 100644 src/mesa_pd/kernel/VBondModel/WallContact.h
create mode 100644 tests/mesa_pd/kernel/VBondModel/WallContact.test.cpp
diff --git a/src/mesa_pd/kernel/VBondModel/WallContact.h b/src/mesa_pd/kernel/VBondModel/WallContact.h
new file mode 100644
index 000000000..24d153ca5
--- /dev/null
+++ b/src/mesa_pd/kernel/VBondModel/WallContact.h
@@ -0,0 +1,113 @@
+//======================================================================================================================
+//
+// 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/Constants.h>
+#include <core/logging/Logging.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+namespace VBondModel {
+
+/**
+ * Repulsive wall interaction kernel.
+ *
+ * Implementation follows:
+ * Wittmaack, Volkov, Zhigilei, "Phase transformation as the mechanism of mechanical deformation of vertically aligned CNT arrays" - Carbon, 2019
+ * https://doi.org/10.1016/j.carbon.2018.11.066
+ *
+ * The force is divided into three areas.
+ * ========================== z=0, position of the wall
+ * close to wall
+ * -------------------------- z=z1+r
+ * spline interpolation
+ * -------------------------- z=z2+r
+ * far away from wall
+ */
+class WallContact
+{
+public:
+ explicit WallContact(real_t zPos) : zPos_(zPos) {}
+
+ template<typename Accessor>
+ void operator()(const size_t p_idx1,
+ Accessor &ac);
+
+ static constexpr real_t r = 6.78_r; ///< A
+ static constexpr real_t eps = 0.254e-3_r; ///< eV/amu
+ static constexpr real_t m = 2648.8_r; ///< amu
+ static constexpr real_t s = 3.6_r; ///< A
+ static constexpr real_t s12 = ((s * s) * (s * s) * (s * s)) * ((s * s) * (s * s) * (s * s));
+
+ static constexpr real_t z1 = 10_r; ///< A
+ static constexpr real_t z2 = 12_r; ///< A
+
+ auto getLastForce() const {return F_;}
+
+ void setZPos(const real_t& zPos) {zPos_ = zPos;}
+ auto getZPos() const {return zPos_;}
+
+private:
+ real_t zPos_ = 0_r;
+ real_t F_ = 0_r; ///< resulting force from the last interaction
+};
+
+template<typename Accessor>
+inline void WallContact::operator()(const size_t p_idx,
+ Accessor &ac)
+{
+ auto dz = ac.getPosition(p_idx)[2] - zPos_;
+ F_ = std::copysign(1_r, dz);
+ dz = std::abs(dz);
+
+ if (dz < r + z1)
+ {
+ //close to wall
+ const auto tmp = dz - r;
+ const auto pow = ((tmp * tmp) * (tmp * tmp) * (tmp * tmp)) * ((tmp * tmp) * (tmp * tmp) * (tmp * tmp)) * tmp;
+ F_ *= m * eps * s12 * 12_r / pow;
+ } else if (dz < r + z2)
+ {
+ //cubic spline interpolation
+ auto S = [](real_t x) { return (3_r * (12_r + r) * (14_r + r) - 6_r * (13_r + r) * x + 3_r * x * x) * 5e-14_r; };
+ F_ *= m * eps * s12 * S(dz);
+ } else
+ {
+ //far away from wall
+ F_ = 0_r;
+ }
+
+ addForceAtomic( p_idx, ac, Vec3(0_r, 0_r, F_) );
+}
+
+} //namespace VBondModel
+} //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 0c3e69f27..b68b61f25 100644
--- a/tests/mesa_pd/CMakeLists.txt
+++ b/tests/mesa_pd/CMakeLists.txt
@@ -164,6 +164,9 @@ waLBerla_execute_test( NAME MESA_PD_Kernel_SyncNextNeighborsBlockForest PROCES
waLBerla_compile_test( NAME MESA_PD_Kernel_TemperatureIntegration FILES kernel/TemperatureIntegration.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_TemperatureIntegration )
+waLBerla_compile_test( NAME MESA_PD_Kernel_VBondModel_WallContact FILES kernel/VBondModel/WallContact.test.cpp DEPENDS core )
+waLBerla_execute_test( NAME MESA_PD_Kernel_VBondModel_WallContact )
+
waLBerla_compile_test( NAME MESA_PD_Kernel_VelocityVerlet FILES kernel/VelocityVerlet.cpp DEPENDS core )
waLBerla_execute_test( NAME MESA_PD_Kernel_VelocityVerlet )
diff --git a/tests/mesa_pd/kernel/VBondModel/WallContact.test.cpp b/tests/mesa_pd/kernel/VBondModel/WallContact.test.cpp
new file mode 100644
index 000000000..f923ee41d
--- /dev/null
+++ b/tests/mesa_pd/kernel/VBondModel/WallContact.test.cpp
@@ -0,0 +1,102 @@
+//======================================================================================================================
+//
+// 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 <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <mesa_pd/kernel/VBondModel/WallContact.h>
+
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+
+#include <iostream>
+
+namespace walberla {
+namespace mesa_pd {
+
+int main( int argc, char ** argv )
+{
+ Environment env(argc, argv);
+ WALBERLA_UNUSED(env);
+ mpi::MPIManager::instance()->useWorldComm();
+
+ if (std::is_same<real_t, float>::value)
+ {
+ WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
+ return EXIT_SUCCESS;
+ }
+
+ //init data structures
+ data::SingleParticleAccessor p;
+ p.setPosition(0, Vec3(0,0,0));
+
+ //init kernels
+ kernel::VBondModel::WallContact wallContact(0_r);
+
+ auto calcForce = [&](const real_t zPos)
+ {
+ clear(p.getForceRef(0));
+ p.getPositionRef(0)[2] = zPos;
+ wallContact(0, p);
+ WALBERLA_CHECK_FLOAT_EQUAL(wallContact.getLastForce(), p.getForce(0)[2]);
+ return p.getForce(0)[2];
+ };
+
+ WALBERLA_LOG_INFO("checking that force is always repulsive");
+ WALBERLA_CHECK_GREATER_EQUAL( calcForce(wallContact.z1 * 0.5_r + wallContact.r),
+ 0_r );
+ WALBERLA_CHECK_GREATER_EQUAL( calcForce((wallContact.z1 + wallContact.z2) * 0.5_r + wallContact.r),
+ 0_r );
+ WALBERLA_CHECK_GREATER_EQUAL( calcForce(wallContact.z2 * 2_r + wallContact.r),
+ 0_r );
+
+ WALBERLA_CHECK_LESS_EQUAL( calcForce(-wallContact.z1 * 0.5_r - wallContact.r),
+ 0_r );
+ WALBERLA_CHECK_LESS_EQUAL( calcForce(-(wallContact.z1 + wallContact.z2) * 0.5_r - wallContact.r),
+ 0_r );
+ WALBERLA_CHECK_LESS_EQUAL( calcForce(-wallContact.z2 * 2_r - wallContact.r),
+ 0_r );
+
+ WALBERLA_LOG_INFO("checking that force is symmetric");
+ WALBERLA_CHECK_FLOAT_EQUAL(calcForce(wallContact.z1 * 0.5_r + wallContact.r),
+ -calcForce(-wallContact.z1 * 0.5_r - wallContact.r) );
+ WALBERLA_CHECK_FLOAT_EQUAL(calcForce((wallContact.z1 + wallContact.z2) * 0.5_r + wallContact.r),
+ -calcForce(-(wallContact.z1 + wallContact.z2) * 0.5_r - wallContact.r) );
+ WALBERLA_CHECK_FLOAT_EQUAL(calcForce(wallContact.z2 * 2_r + wallContact.r),
+ -calcForce(-wallContact.z2 * 2_r - wallContact.r) );
+
+ WALBERLA_LOG_INFO("checking smoothness of the force");
+ WALBERLA_CHECK_FLOAT_EQUAL( calcForce(std::nextafter(wallContact.z1 + wallContact.r, 0_r)),
+ calcForce(std::nextafter(wallContact.z1 + wallContact.r, std::numeric_limits<real_t>::infinity())) );
+ WALBERLA_CHECK_FLOAT_EQUAL( calcForce(std::nextafter(wallContact.z2 + wallContact.r, 0_r)),
+ calcForce(std::nextafter(wallContact.z2 + wallContact.r, std::numeric_limits<real_t>::infinity())) );
+
+ return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+ return walberla::mesa_pd::main(argc, argv);
+}
--
GitLab