diff --git a/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h b/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h new file mode 100644 index 0000000000000000000000000000000000000000..a92ecda09637c90da12722611665836875016328 --- /dev/null +++ b/src/mesa_pd/kernel/cnt/IsotropicVDWContact.h @@ -0,0 +1,106 @@ +//====================================================================================================================== +// +// 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 diff --git a/tests/mesa_pd/CMakeLists.txt b/tests/mesa_pd/CMakeLists.txt index b68b61f25dd8009fb01968e61d37c063410f53bc..07a3b6b6b4927bf3408b4b6ae3008ef76f1069e4 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_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 ) diff --git a/tests/mesa_pd/kernel/cnt/IsotropicVDWContact.test.cpp b/tests/mesa_pd/kernel/cnt/IsotropicVDWContact.test.cpp new file mode 100644 index 0000000000000000000000000000000000000000..60e4e6eca4450b02359ab963e61064243cd0c339 --- /dev/null +++ b/tests/mesa_pd/kernel/cnt/IsotropicVDWContact.test.cpp @@ -0,0 +1,87 @@ +//====================================================================================================================== +// +// 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); +}