[ADD] benchmark

apps/benchmarks/CMakeLists.txt

 add_subdirectory( AdaptiveMeshRefinementFluidParticleCoupling )
+add_subdirectory( CNT )
 add_subdirectory( ComplexGeometry )
 add_subdirectory( DEM )
 add_subdirectory( MeshDistance )

apps/benchmarks/CNT/Accessor.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/data/ParticleAccessor.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+
+namespace walberla {
+namespace mesa_pd {
+
+class Accessor : public data::ParticleAccessor
+{
+public:
+   Accessor(std::shared_ptr<data::ParticleStorage>& ps)
+   : ParticleAccessor(ps)
+   {}
+
+   constexpr auto getInvMass(const size_t /*p_idx*/) const {return 1_r / kernel::cnt::mass_T;}
+
+   constexpr auto& getInvInertiaBF(const size_t /*p_idx*/) const {return invI;}
+private:
+   //  - sphere   :  I = (2/5)*mass*radius^2
+   static constexpr auto Ia = 0.4_r * kernel::cnt::mass_T * kernel::cnt::inner_radius * kernel::cnt::inner_radius;
+   static constexpr auto invI = Mat3(1_r/Ia, 0_r, 0_r, 0_r, 1_r/Ia, 0_r, 0_r, 0_r, 1_r/Ia);
+};
+
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file

apps/benchmarks/CNT/CMakeLists.txt

+waLBerla_add_executable ( NAME VBondModel
+                          FILES VBondModel.cpp
+                          DEPENDS core mesa_pd )

apps/benchmarks/CNT/VBondModel.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 "Accessor.h"
+#include "mesa_pd/data/Flags.h"
+#include "mesa_pd/data/ParticleStorage.h"
+#include "mesa_pd/kernel/ParticleSelector.h"
+#include "mesa_pd/kernel/cnt/Parameters.h"
+#include "mesa_pd/kernel/cnt/VBondContact.h"
+#include "mesa_pd/kernel/VelocityVerlet.h"
+
+#include "core/Environment.h"
+#include "core/math/Constants.h"
+
+namespace walberla {
+using namespace walberla::mesa_pd;
+
+int main(int argc, char **argv)
+{
+   Environment env(argc, argv);
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   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 = 500ll;
+
+   WALBERLA_LOG_INFO_ON_ROOT("creating initial particle setup");
+   auto ps = std::make_shared<data::ParticleStorage>(10);
+   auto ac = Accessor(ps);
+
+   for (auto i = 0; i < 10; ++i)
+   {
+      data::Particle &&p = *ps->create();
+      p.setPosition(Vec3(0_r, 0_r, real_c(i) * 13.56_r));
+      p.setSegmentID(i);
+      p.setClusterID(1);
+      if (i == 0)
+         data::particle_flags::set(p.getFlagsRef(), data::particle_flags::FIXED);
+      p.getRotationRef().rotate(Vec3(0_r, 1_r, 0_r), -0.5_r * math::pi);
+      p.getRotationRef().rotate(Vec3(0_r, 0_r, 1_r), 0_r);
+   }
+   data::Particle &&last_segment = *(ps->end() - 1);
+
+   WALBERLA_LOG_INFO_ON_ROOT("setting up interaction models");
+   kernel::cnt::VBondContact vbond;
+   kernel::VelocityVerletPreForceUpdate vv_pre(kernel::cnt::dT);
+   kernel::VelocityVerletPostForceUpdate vv_post(kernel::cnt::dT);
+
+   WALBERLA_LOG_INFO_ON_ROOT("running simulation");
+   auto appliedForce = Vec3(1_r, 0_r, 0_r);
+   auto appliedTorque = Vec3(0_r);
+   std::ofstream fout("output.txt");
+   for (auto i = 0; i < numSimulationSteps; ++i)
+   {
+      ps->forEachParticle(false,
+                          kernel::SelectAll(),
+                          ac,
+                          vv_pre,
+                          ac);
+
+      last_segment.setForce(appliedForce);
+      last_segment.setTorque(appliedTorque);
+      constexpr auto cutoff2 = kernel::cnt::outer_radius * kernel::cnt::outer_radius;
+
+      real_t tensileEnergy = 0_r;
+      real_t shearEnergy = 0_r;
+      real_t bendingEnergy = 0_r;
+      real_t twistingEnergy = 0_r;
+
+      ps->forEachParticlePairHalf(false,
+                                  kernel::SelectAll(),
+                                  ac,
+                                  [&](size_t p_idx1, size_t p_idx2, Accessor &ac)
+                                  {
+                                     if ((ac.getPosition(p_idx1) - ac.getPosition(p_idx2)).sqrLength() < cutoff2)
+                                     {
+                                        vbond(p_idx1, p_idx2, ac);
+                                        tensileEnergy += vbond.tensileEnergy;
+                                        shearEnergy += vbond.shearEnergy;
+                                        bendingEnergy += vbond.bendingEnergy;
+                                        twistingEnergy += vbond.twistingEnergy;
+                                     }
+                                  },
+                                  ac);
+
+      ps->forEachParticle(false,
+                          kernel::SelectAll(),
+                          ac,
+                          vv_post,
+                          ac);
+
+      fout << last_segment.getPosition()[0] << " "
+           << tensileEnergy << " "
+           << shearEnergy << " "
+           << bendingEnergy << " "
+           << twistingEnergy
+           << std::endl;
+   }
+
+   return EXIT_SUCCESS;
+}
+} // namespace walberla
+
+int main(int argc, char *argv[])
+{
+   return walberla::main(argc, argv);
+}

src/core/math/Constants.h

@@ -34,19 +34,21 @@ namespace walberla
 {
 namespace math
 {
-constexpr real_t e                = real_t(2.7182818284590452354);  /* e */
-constexpr real_t log2_e           = real_t(1.4426950408889634074);  /* log_2 e */
-constexpr real_t log10_e          = real_t(0.43429448190325182765); /* log_10 e */
-constexpr real_t ln_two           = real_t(0.69314718055994530942); /* log_e 2 */
-constexpr real_t ln_ten           = real_t(2.30258509299404568402); /* log_e 10 */
-constexpr real_t pi               = real_t(3.14159265358979323846); /* pi */
-constexpr real_t half_pi          = real_t(1.57079632679489661923); /* pi/2 */
-constexpr real_t fourth_pi        = real_t(0.78539816339744830962); /* pi/4 */
-constexpr real_t one_div_pi       = real_t(0.31830988618379067154); /* 1/pi */
-constexpr real_t two_div_pi       = real_t(0.63661977236758134308); /* 2/pi */
-constexpr real_t two_div_root_pi  = real_t(1.12837916709551257390); /* 2/sqrt(pi) */
-constexpr real_t root_two         = real_t(1.41421356237309504880); /* sqrt(2) */
-constexpr real_t one_div_root_two = real_t(0.70710678118654752440); /* 1/sqrt(2) */
+constexpr real_t e                  = real_t(2.7182818284590452354);  /* e */
+constexpr real_t log2_e             = real_t(1.4426950408889634074);  /* log_2 e */
+constexpr real_t log10_e            = real_t(0.43429448190325182765); /* log_10 e */
+constexpr real_t ln_two             = real_t(0.69314718055994530942); /* log_e 2 */
+constexpr real_t ln_ten             = real_t(2.30258509299404568402); /* log_e 10 */
+constexpr real_t pi                 = real_t(3.14159265358979323846); /* pi */
+constexpr real_t half_pi            = real_t(1.57079632679489661923); /* pi/2 */
+constexpr real_t fourth_pi          = real_t(0.78539816339744830962); /* pi/4 */
+constexpr real_t one_div_pi         = real_t(0.31830988618379067154); /* 1/pi */
+constexpr real_t two_div_pi         = real_t(0.63661977236758134308); /* 2/pi */
+constexpr real_t two_div_root_pi    = real_t(1.12837916709551257390); /* 2/sqrt(pi) */
+constexpr real_t root_two           = real_t(1.41421356237309504880); /* sqrt(2) */
+constexpr real_t one_div_root_two   = real_t(0.70710678118654752440); /* 1/sqrt(2) */
+constexpr real_t root_three         = real_t(1.73205080757); /* sqrt(3) */
+constexpr real_t one_div_root_three = real_t(0.57735026919); /* 1/sqrt(3) */
 
 } // namespace math
 } // namespace walberla
\ No newline at end of file

src/mesa_pd/kernel/cnt/Parameters.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 "core/DataTypes.h"
+#include "core/math/Constants.h"
+
+namespace walberla {
+namespace mesa_pd {
+namespace kernel {
+namespace cnt {
+
+//====================CNT inertial and elastic parameters========================================================
+
+// Within our model, CNTs are represented with intersecting capsule primitives - cylinders with capped ends.
+// Capped ends are necessary for a continuous normal definition
+// Each capsule represents the inertial properties of a cylindrical segment.
+
+// For the units system used see PFC 5 implementation docs
+
+/// timestep is fixed to 20 fs.
+constexpr auto dT = 20_r;
+
+/// vdW interaction radius w/r to inertial segment radius
+constexpr auto CutoffFactor     = 4_r;
+
+// CNT lattice numbers (m,n) CNT
+constexpr auto mm = 10_r;
+constexpr auto nn = 10_r;
+
+/// CNT Young modulus in GPa based on atomistic simulations
+constexpr auto En = 1029_r;
+/// CNT shear modulus in GPa based on atomistic simulations
+constexpr auto Gs = 459_r;
+
+/// Equilibrium distance of a covalent C-C bond
+constexpr auto a_CC = 1.42_r;
+
+/// Equilibrium vdW separation of two CNT surfaces
+constexpr auto a_VD = 3.35_r;
+
+/// CNT radius
+constexpr auto R_CNT = 6.78_r; //( a_CC / (2 * PII)) * std::sqrt( 3.0 * ( mm * mm + nn * nn + mm * nn ) );
+
+/// external radius of an idealized hollow cylindrical shell
+constexpr auto R_1 = R_CNT + 0.5_r * a_VD;
+/// internal radius of an idealized hollow cylindrical shell
+constexpr auto R_2 = R_CNT - 0.5_r * a_VD;
+
+/// height of a cylindrical segment
+constexpr auto T = 2_r * R_CNT;
+
+/// linear density in atoms per A
+constexpr auto ro = 4_r * mm / ( math::root_three * a_CC );
+
+/// Atomic mass of Carbon in AMU
+constexpr auto M_C = 12.011_r;
+
+/// Mass of the repetative cell in AMU
+constexpr auto mass_T = ro * T * M_C * 104.397_r;
+
+// Volume of a capsule
+//double vol = (4./3.) * PII * (R_CNT * R_CNT * R_CNT) + PII * R_CNT * R_CNT * T;
+
+// Density of a capsule
+//double dens = mass_T / vol;
+
+/// V-bond parameter
+constexpr auto knorm = (En * 0.006242_r) / T;
+/// V-bond parameter
+constexpr auto kshear = (Gs * 0.006242_r) / T;
+
+constexpr auto margin = 10.e-10_r;
+
+constexpr auto outer_radius     = CutoffFactor * R_CNT;
+// real_t inner_radius     = R_CNT; // Capsule
+constexpr auto inner_radius     = 1.5811_r * R_CNT; // Sphere
+
+/// volume of a sphere
+constexpr auto vol = (4_r/3_r) * math::pi * (inner_radius * inner_radius * inner_radius);
+
+/// density of a sphere
+constexpr auto dens = mass_T / vol;
+
+} //namespace cnt
+} //namespace kernel
+} //namespace mesa_pd
+} //namespace walberla
\ No newline at end of file

src/mesa_pd/kernel/VBondModel/VBondContact.h → src/mesa_pd/kernel/cnt/VBondContact.h

@@ -40,10 +40,17 @@
 namespace walberla {
 namespace mesa_pd {
 namespace kernel {
-namespace VBondModel {
+namespace cnt {
 
 /**
- * Advanced DEM kernel
+ * VBond interaction model
+ *
+ * Implementation according to:
+ *
+ * Ostanin, I., Dumitrică, T., Eibl, S., and Rüde, U. (October 4, 2019).
+ * "Size-Independent Mechanical Response of Ultrathin Carbon Nanotube Films in Mesoscopic Distinct Element Method Simulations."
+ * ASME. J. Appl. Mech. December 2019; 86(12): 121006.
+ * https://doi.org/10.1115/1.4044413
  */
 class VBondContact
 {
@@ -61,7 +68,7 @@ public:
    static constexpr real_t Jp = 2_r * J;
 
    // Stiffnesses, equilibrium length etc
-   static constexpr real_t B1 = E * S / a; // Need to calibrate
+   static constexpr real_t B1 = E * S / a;
    static constexpr real_t B2 = 12_r * E * J / a;
    static constexpr real_t B3 = -2_r * E * J / a - G * Jp / (2_r * a);
    static constexpr real_t B4 = G * Jp / a;
@@ -69,28 +76,13 @@ public:
    Vector3<bool> isPeriodic_;
    Vector3<int64_t> maxSegments_;
 
-
-   // Calibration through atomistic simulation
-   //constexpr real_t B1 = 60.1_r;
-   //constexpr real_t B2 = 17100_r;
-   //constexpr real_t B3 = 3610_r;
-   //constexpr real_t B4 = 1107_r;
-
-   real_t U;
-   real_t U_1;
-   real_t U_2;
-   real_t U_3;
-   real_t U_4;
+   real_t tensileEnergy;
+   real_t shearEnergy;
+   real_t bendingEnergy;
+   real_t twistingEnergy;
 };
 
-/**
- *
- * @tparam Accessor
- * @param p_idx1
- * @param p_idx2
- * @param ac
- * @return vdW adhesion energy
- */
+
 template<typename Accessor>
 inline
 void VBondContact::operator()(const size_t p_idx1,
@@ -121,28 +113,10 @@ void VBondContact::operator()(const size_t p_idx1,
 
    real_t C = 1_r;
 
-   /*
-           //srand(SID_1);
-           double AA = 0.2 + 0.2 * ((rand() % 1000)/1000); //strain of initiation of stress corrosion
-           //double BB = 0.01; // width of stress corrosion band
-           //double DD = 1. + AA/BB;
-           //double KK = 1 / BB;
-           //if (((std::fabs((s-a)/s))>AA)&&((std::fabs((s-a)/s))<=AA+BB)) C = DD - std::fabs((s-a)/s) * KK;
-           //if ((std::fabs((s-a)/s))>AA+BB) C = 0;
-           if ((std::fabs((s-a)/s))>AA) C = 0;
-
-           // Potential energy
-           real_t U_1 = C * 0.5 * B1 * (s - a) * (s - a);
-           real_t U_2 = C * 0.5 * B2 * (nj1 - ni1) * dij + B2;
-           real_t U_3 = C * B3 * ni1 * nj1 + B3;
-           real_t U_4 = C * -0.5 * B4 * (ni2 * nj2 + ni3 * nj3) + B4;
-           real_t U = U_1 + U_2 + U_3 + U_4;
-   */
-   U_1 = 0.5_r * B1 * (s - a) * (s - a);
-   U_2 = B2 * (0.5_r * (nj1 - ni1) * dij - 0.25_r * ni1 * nj1 + 0.75_r);
-   U_3 = (0.25_r * B2 + B3 + 0.5_r * B4) * (ni1 * nj1 + 1_r);
-   U_4 = -0.5_r * B4 * (ni1 * nj1 + ni2 * nj2 + ni3 * nj3 - 1_r);
-   U = U_1 + U_2 + U_3 + U_4;
+   tensileEnergy  = 0.5_r * B1 * (s - a) * (s - a);
+   shearEnergy    = B2 * (0.5_r * (nj1 - ni1) * dij - 0.25_r * ni1 * nj1 + 0.75_r);
+   bendingEnergy  = (0.25_r * B2 + B3 + 0.5_r * B4) * (ni1 * nj1 + 1_r);
+   twistingEnergy = -0.5_r * B4 * (ni1 * nj1 + ni2 * nj2 + ni3 * nj3 - 1_r);
 
    Vec3 rij = dij;
 
@@ -158,7 +132,7 @@ void VBondContact::operator()(const size_t p_idx1,
    addTorqueAtomic(p_idx2, ac, Mji);
 }
 
-} //namespace VBondModel
+} //namespace cnt
 } //namespace kernel
 } //namespace mesa_pd
 } //namespace walberla
\ No newline at end of file

tests/mesa_pd/kernel/VBondModel/VBondContact.test.cpp

@@ -23,7 +23,7 @@
 #include <mesa_pd/data/ParticleAccessor.h>
 #include <mesa_pd/data/ParticleStorage.h>
 
-#include <mesa_pd/kernel/VBondModel/VBondContact.h>
+#include <mesa_pd/kernel/cnt/VBondContact.h>
 
 #include <core/Environment.h>
 #include <core/logging/Logging.h>
@@ -61,7 +61,7 @@ int main( int argc, char ** argv )
    data::ParticleAccessor accessor(ps);
 
    //init kernels
-   kernel::VBondModel::VBondContact vbond;
+   kernel::cnt::VBondContact vbond;
 
    //check equilibrium distance
    vbond(0, 1, accessor);