diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index d16b4255d1c2d48cc25d1d9587f453a7c321ce28..420dc32d3c973993890656ce458d65cd32a2d835 100644
--- a/apps/benchmarks/CMakeLists.txt
+++ b/apps/benchmarks/CMakeLists.txt
@@ -1,4 +1,5 @@
add_subdirectory( AdaptiveMeshRefinementFluidParticleCoupling )
+add_subdirectory( CNT )
add_subdirectory( ComplexGeometry )
add_subdirectory( DEM )
add_subdirectory( MeshDistance )
diff --git a/apps/benchmarks/CNT/Accessor.h b/apps/benchmarks/CNT/Accessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..d23002562eee622245ae2865057b6ff89937d537
--- /dev/null
+++ b/apps/benchmarks/CNT/Accessor.h
@@ -0,0 +1,48 @@
+//======================================================================================================================
+//
+// 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
diff --git a/apps/benchmarks/CNT/CMakeLists.txt b/apps/benchmarks/CNT/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..320d0ce94b9ced0ffd74e9c09f8e74b3137f936f
--- /dev/null
+++ b/apps/benchmarks/CNT/CMakeLists.txt
@@ -0,0 +1,3 @@
+waLBerla_add_executable ( NAME VBondModel
+ FILES VBondModel.cpp
+ DEPENDS core mesa_pd )
diff --git a/apps/benchmarks/CNT/VBondModel.cpp b/apps/benchmarks/CNT/VBondModel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..091327c2f0290032abc36b85e3fe3d32b7022311
--- /dev/null
+++ b/apps/benchmarks/CNT/VBondModel.cpp
@@ -0,0 +1,128 @@
+//======================================================================================================================
+//
+// 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);
+}
diff --git a/src/core/math/Constants.h b/src/core/math/Constants.h
index 92e49678d67bce580b9ce0dd46fd004f96609e28..5487e14cf34f651a9b1065520d182fae3b4e9c19 100644
--- a/src/core/math/Constants.h
+++ b/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
diff --git a/src/mesa_pd/kernel/cnt/Parameters.h b/src/mesa_pd/kernel/cnt/Parameters.h
new file mode 100644
index 0000000000000000000000000000000000000000..63f256a04dca173007766af03f3b70256047824f
--- /dev/null
+++ b/src/mesa_pd/kernel/cnt/Parameters.h
@@ -0,0 +1,108 @@
+//======================================================================================================================
+//
+// 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
diff --git a/src/mesa_pd/kernel/VBondModel/VBondContact.h b/src/mesa_pd/kernel/cnt/VBondContact.h
similarity index 73%
rename from src/mesa_pd/kernel/VBondModel/VBondContact.h
rename to src/mesa_pd/kernel/cnt/VBondContact.h
index 716fdc195b5382e0ce5820211e61d709963083e9..f245aca1eda27cdb805a4d426abd3496795d54eb 100644
--- a/src/mesa_pd/kernel/VBondModel/VBondContact.h
+++ b/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
diff --git a/tests/mesa_pd/kernel/VBondModel/VBondContact.test.cpp b/tests/mesa_pd/kernel/VBondModel/VBondContact.test.cpp
index 941b4c2856e21a877eb4553a2aeb9ef4ec07b3ca..796eca1063d805b5a366f70323cc99212c3d1e2f 100644
--- a/tests/mesa_pd/kernel/VBondModel/VBondContact.test.cpp
+++ b/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);