Merge branch 'vbondmodel_isotropic' into 'master'

[ADD] kernel for the VBondModel isotropic vdW contact

See merge request walberla/walberla!377

src/mesa_pd/kernel/cnt/IsotropicVDWContact.h

+//======================================================================================================================
+//
+//  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 cnt {
+
+/**
+ * Isotropic vdW contact model.
+ *
+ * implementation follows:
+ * I. Ostanin, R. Ballarini, D. Potyondy, T. Dumitrica, A distinct element method for large scale simulations of carbon nanotube assemblies, J. Mech. Phys. Solid. 61 (2013) 762-782.
+ * https://doi.org/10.1016/j.jmps.2012.10.016
+ */
+class IsotropicVDWContact
+{
+public:
+   template<typename Accessor>
+   void operator()(const size_t p_idx1,
+                   const size_t p_idx2,
+                   Accessor &ac);
+   static real_t equilibriumDistance();
+
+   static constexpr real_t eps = 0.07124_r;
+   static constexpr real_t A = 0.0223_r;
+   static constexpr real_t B = 1.31_r;
+   static constexpr real_t alpha = 9.5_r;
+   static constexpr real_t beta = 4.0_r;
+   static constexpr real_t r = 6.78_r;
+   static constexpr real_t Rinv = 1.0_r / r;
+   static constexpr real_t Dc = 0.4_r;
+
+   auto getLastForce() const {return F_;}
+   auto getLastEnergy() const {return U_;}
+private:
+   real_t F_ = 0_r; ///< resulting force from the last interaction
+   real_t U_ = 0_r; ///< resulting energy from the last interaction
+};
+
+template<typename Accessor>
+inline
+void IsotropicVDWContact::operator()(const size_t p_idx1,
+                                     const size_t p_idx2,
+                                     Accessor &ac)
+{
+   Vec3 n = ac.getPosition(p_idx2) - ac.getPosition(p_idx1); ///< contact normal
+   real_t L = n.length();
+   n *= (1_r/L);
+
+   real_t D = L * Rinv - 2_r;
+   F_ = 0.01_r; //default value
+   U_ = 0_r; //default value
+
+   if (D > Dc)
+   {
+      const auto pow_alpha = std::pow(D, alpha);
+      const auto pow_beta = std::pow(D, beta);
+      F_ = 4_r * eps * (-(alpha * A) / (pow_alpha * D) + (beta * B) / (pow_beta * D));
+      U_ = 4_r * eps * (A / pow_alpha - B / pow_beta);
+   }
+
+   Vec3 force = n * F_;
+   addForceAtomic(p_idx1, ac,  force);
+   addForceAtomic(p_idx2, ac, -force);
+}
+
+real_t IsotropicVDWContact::equilibriumDistance()
+{
+   return r * ( std::pow( (alpha*A)/(beta*B), 1_r/(alpha-beta)) + 2_r);
+}
+
+} //namespace cnt
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file

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_IsotropicVDWContact FILES kernel/cnt/IsotropicVDWContact.test.cpp DEPENDS core )
+waLBerla_execute_test( NAME   MESA_PD_Kernel_CNT_IsotropicVDWContact )
+
 waLBerla_compile_test( NAME   MESA_PD_Kernel_CoefficientOfRestitutionSD FILES kernel/CoefficientOfRestitutionSD.cpp DEPENDS core )
 waLBerla_execute_test( NAME   MESA_PD_Kernel_CoefficientOfRestitutionSDEuler COMMAND $<TARGET_FILE:MESA_PD_Kernel_CoefficientOfRestitutionSD> )
 waLBerla_execute_test( NAME   MESA_PD_Kernel_CoefficientOfRestitutionSDVelocityVerlet COMMAND $<TARGET_FILE:MESA_PD_Kernel_CoefficientOfRestitutionSD> --useVV )

tests/mesa_pd/kernel/cnt/IsotropicVDWContact.test.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
+//! \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/cnt/IsotropicVDWContact.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
+   auto ps = std::make_shared<data::ParticleStorage>(2);
+   data::Particle &&p1 = *ps->create();
+   data::Particle &&p2 = *ps->create();
+
+   data::ParticleAccessor accessor(ps);
+
+   //init kernels
+   kernel::cnt::IsotropicVDWContact vdWContact;
+
+   auto calcForce = [&](const Vec3 pos1, const Vec3 pos2) {
+      p1.setPosition(pos1);
+      p2.setPosition(pos2);
+      clear(p1.getForceRef());
+      clear(p2.getForceRef());
+      vdWContact(0, 1, accessor);
+      return p1.getForce();
+   };
+
+   const Vec3 randomNormal = Vec3(1_r, 2_r, 3_r).getNormalized();
+
+   WALBERLA_LOG_INFO("checking repulsion - equilibrium - attraction");
+   vdWContact(0, 1, accessor);
+   WALBERLA_CHECK_LESS       (dot(randomNormal, calcForce(Vec3(0), randomNormal * (vdWContact.equilibriumDistance() - 1_r))),
+                              0_r);
+   WALBERLA_CHECK_FLOAT_EQUAL(dot(randomNormal, calcForce(Vec3(0), randomNormal * vdWContact.equilibriumDistance())),
+                              0_r);
+   WALBERLA_CHECK_GREATER    (dot(randomNormal, calcForce(Vec3(0), randomNormal * (vdWContact.equilibriumDistance() + 1_r))),
+                              0_r);
+
+   return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main(int argc, char **argv)
+{
+   return walberla::mesa_pd::main(argc, argv);
+}