added benchmark as integration test

apps/benchmarks/CMakeLists.txt

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

apps/benchmarks/CNT/CMakeLists.txt

-waLBerla_add_executable ( NAME VBondModel
-                          FILES VBondModel.cpp
-                          DEPENDS core mesa_pd )
-#waLBerla_execute_test( NAME MESA_PD_CNT_VBondModel
-#                       COMMAND $<TARGET_FILE:VBondModel> )

src/mesa_pd/kernel/cnt/VBondContact.h

@@ -55,6 +55,12 @@ namespace cnt {
 class VBondContact
 {
 public:
+   VBondContact( const Vector3<bool>& isPeriodic = Vector3<bool>{false, false, false},
+                 const Vector3<int64_t>& maxSegments = Vector3<real_t>{0, 0, 0} )
+   : isPeriodic_(isPeriodic)
+   , maxSegments_(maxSegments)
+   {}
+
    template<typename Accessor>
    void operator()(const size_t p_idx1,
                    const size_t p_idx2,
@@ -73,13 +79,24 @@ public:
    static constexpr real_t B3 = -2_r * E * J / a - G * Jp / (2_r * a);
    static constexpr real_t B4 = G * Jp / a;
 
-   Vector3<bool> isPeriodic_;
-   Vector3<int64_t> maxSegments_;
 
+
+   /// Get tensile energy from last contact.
+   auto getLastTensileEnergy() const {return tensileEnergy;}
+   /// Get shear energy from last contact.
+   auto getLastShearEnergy() const {return shearEnergy;}
+   /// Get bending energy from last contact.
+   auto getLastBendingEnergy() const {return bendingEnergy;}
+   /// Get twisting energy from last contact.
+   auto getLastTwistingEnergy() const {return twistingEnergy;}
+private:
    real_t tensileEnergy;
    real_t shearEnergy;
    real_t bendingEnergy;
    real_t twistingEnergy;
+
+   Vector3<bool> isPeriodic_;
+   Vector3<int64_t> maxSegments_;
 };
 
 

tests/mesa_pd/CMakeLists.txt

@@ -90,6 +90,9 @@ waLBerla_execute_test( NAME   MESA_PD_Kernel_CNT_IsotropicVDWContact )
 waLBerla_compile_test( NAME   MESA_PD_Kernel_CNT_VBondContact FILES kernel/cnt/VBondContact.test.cpp DEPENDS core )
 waLBerla_execute_test( NAME   MESA_PD_Kernel_CNT_VBondContact )
 
+waLBerla_compile_test( NAME   MESA_PD_Kernel_CNT_VBondContactIntegration FILES kernel/cnt/VBondContactIntegration.test.cpp DEPENDS core )
+waLBerla_execute_test( NAME   MESA_PD_Kernel_CNT_VBondContactIntegration )
+
 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 )

apps/benchmarks/CNT/Accessor.h → tests/mesa_pd/kernel/cnt/SphericalSegmentAccessor.h

@@ -28,10 +28,10 @@
 namespace walberla {
 namespace mesa_pd {
 
-class Accessor : public data::ParticleAccessor
+class SphericalSegmentAccessor : public data::ParticleAccessor
 {
 public:
-   Accessor(std::shared_ptr<data::ParticleStorage>& ps)
+   SphericalSegmentAccessor(std::shared_ptr<data::ParticleStorage>& ps)
    : ParticleAccessor(ps)
    {}
 

tests/mesa_pd/kernel/cnt/VBondContactIntegration.ipynb

+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "northern-indie",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import matplotlib.pyplot as pl\n",
+    "import numpy as np\n",
+    "import seaborn as sns\n",
+    "\n",
+    "%matplotlib inline\n",
+    "%config InlineBackend.figure_formats = ['svg']\n",
+    "\n",
+    "sns.set_context(\"talk\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "distant-tsunami",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "energies = np.loadtxt('VBondContactEnergies.txt').transpose()\n",
+    "ref = np.loadtxt('VBondContactEnergies.ref.txt')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "failing-studio",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pl.figure( figsize=(10,10) )\n",
+    "pl.subplot(221)\n",
+    "pl.title('tensile')\n",
+    "pl.plot(energies[1], label='MESA-PD')\n",
+    "pl.plot(ref[0], label='reference')\n",
+    "pl.legend()\n",
+    "\n",
+    "pl.subplot(222)\n",
+    "pl.title('shear')\n",
+    "pl.plot(energies[2])\n",
+    "pl.plot(ref[1])\n",
+    "\n",
+    "pl.subplot(223)\n",
+    "pl.title('bending')\n",
+    "pl.plot(energies[3])\n",
+    "pl.plot(ref[2])\n",
+    "\n",
+    "pl.subplot(224)\n",
+    "pl.title('twisting')\n",
+    "pl.plot(energies[4])\n",
+    "pl.plot(ref[3])\n",
+    "\n",
+    "pl.tight_layout()\n",
+    "pl.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "enormous-salvation",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.5"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

apps/benchmarks/CNT/VBondModel.cpp → tests/mesa_pd/kernel/cnt/VBondContactIntegration.test.cpp

@@ -20,7 +20,7 @@
 //
 //======================================================================================================================
 
-#include "Accessor.h"
+#include "SphericalSegmentAccessor.h"
 #include "mesa_pd/data/Flags.h"
 #include "mesa_pd/data/ParticleStorage.h"
 #include "mesa_pd/kernel/ParticleSelector.h"
@@ -41,6 +41,12 @@ int main(int argc, char **argv)
    Environment env(argc, argv);
    walberla::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;
+   }
+
    logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
    logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
 
@@ -49,7 +55,7 @@ int main(int argc, char **argv)
 
    WALBERLA_LOG_INFO_ON_ROOT("creating initial particle setup");
    auto ps = std::make_shared<data::ParticleStorage>(10);
-   auto ac = Accessor(ps);
+   auto ac = SphericalSegmentAccessor(ps);
 
    for (auto i = 0; i < 10; ++i)
    {
@@ -89,7 +95,7 @@ int main(int argc, char **argv)
    real_t bendingEnergy = 0_r;
    real_t twistingEnergy = 0_r;
 
-   std::ofstream fout("output.txt");
+   std::ofstream fout("VBondContactEnergies.txt");
    for (auto i = 0; i < numSimulationSteps; ++i)
    {
       vtkWriter->write();
@@ -119,10 +125,10 @@ int main(int argc, char **argv)
                                      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;
+                                        tensileEnergy += vbond.getLastTensileEnergy();
+                                        shearEnergy += vbond.getLastShearEnergy();
+                                        bendingEnergy += vbond.getLastBendingEnergy();
+                                        twistingEnergy += vbond.getLastTwistingEnergy();
                                      }
                                   });
 
@@ -151,10 +157,5 @@ int main(int argc, char **argv)
 
 int main(int argc, char *argv[])
 {
-   if (std::is_same<walberla::real_t, float>::value)
-   {
-      WALBERLA_LOG_WARNING("waLBerla build in sp mode: skipping test due to low precision");
-      return EXIT_SUCCESS;
-   }
    return walberla::main(argc, argv);
 }