Merge branch 'mr_hashgrids_for_mesapd' into 'master'

HashGrids for mesapd

See merge request walberla/walberla!457

python/mesa_pd.py

@@ -90,6 +90,7 @@ if __name__ == '__main__':
     cs.add_property("diag_n_inv", "real_t", defValue="real_t(0)")
     cs.add_property("p", "walberla::mesa_pd::Vec3", defValue="real_t(0)")
 
+    mpd.add(data.HashGrids())
     mpd.add(data.LinkedCells())
     mpd.add(data.SparseLinkedCells())
     mpd.add(data.ShapeStorage(ps))

python/mesa_pd/data/HashGrids.py

+# -*- coding: utf-8 -*-
+
+from ..utility import generate_file
+
+
+class HashGrids:
+    def generate(self, module):
+        ctx = {'module': module}
+        generate_file(module['module_path'], 'data/HashGrids.templ.h', ctx)

python/mesa_pd/data/__init__.py

@@ -2,6 +2,7 @@
 
 from .ContactHistory import ContactHistory
 from .ContactStorage import ContactStorage
+from .HashGrids import HashGrids
 from .LinkedCells import LinkedCells
 from .ParticleStorage import ParticleStorage
 from .ShapeStorage import ShapeStorage
@@ -9,6 +10,7 @@ from .SparseLinkedCells import SparseLinkedCells
 
 __all__ = ['ContactHistory',
            'ContactStorage',
+           'HashGrids',
            'LinkedCells',
            'ParticleStorage',
            'ShapeStorage',

python/mesa_pd/templates/kernel/DetectAndStoreContacts.templ.h

@@ -41,7 +41,7 @@ namespace mesa_pd {
 namespace kernel {
 
 /**
- * Kernel which performes collision detection on a pair of two particles
+ * Kernel which performs collision detection on a pair of two particles
  * and inserts the contact (if existent) into the contact storage passed in the constructor.
  * Call this kernel on each particle pair to perform contact detection and insert each contact in the contact
  * storage.

python/mesa_pd/templates/kernel/ForceLJ.templ.h

@@ -37,7 +37,7 @@ namespace mesa_pd {
 namespace kernel {
 
 /**
- * Kernel which calculates the Lennard Jones froce between two particles.
+ * Kernel which calculates the Lennard Jones force between two particles.
  *
  * This kernel uses the type property of a particle to decide on the material parameters.
  *

python/mesa_pd/templates/kernel/HCSITSRelaxationStep.templ.h

@@ -998,7 +998,7 @@ inline real_t HCSITSRelaxationStep::relaxInelasticContactsByProjectedGaussSeidel
       // No need to apply zero impulse.
    }
    else {
-      // Dissect the impuls in a tangetial and normal directions
+      // Dissect the impuls in a tangential and normal directions
       // Use the inverted contact frame with -n as in the publication
       Mat3 reversedContactFrame( -ca.getNormal(cid), ca.getT(cid), ca.getO(cid) );
       Vec3 p_rcf( reversedContactFrame.getTranspose() * p_wf );

python/mesa_pd/templates/kernel/HeatConduction.templ.h

@@ -40,7 +40,7 @@ namespace kernel {
 /**
  * Heat conduction interaction kernel
  *
- * This kernel implements a simple heat conduction: \frac{dQ}{dt} = - \alpha (T_2 - T_1)
+ * This kernel implements a simple heat conduction: \f$ \frac{dQ}{dt} = - \alpha (T_2 - T_1) \f$
  *
  * \code
    {%- for prop in interface %}

python/mesa_pd/templates/kernel/InitParticlesForHCSITS.templ.h

@@ -40,7 +40,7 @@ namespace kernel {
 /**
  * Init the datastructures for the particles for later use of the HCSITS-Solver.
  * Call this kernel on all particles that will be treated with HCSITS before performing any relaxation timesteps.
- * Use setGlobalAcceleration() to set an accelleration action uniformly across all particles (e.g. gravity)
+ * Use setGlobalAcceleration() to set an acceleration action uniformly across all particles (e.g. gravity)
  * \ingroup mesa_pd_kernel
  */
 class InitParticlesForHCSITS

python/mesa_pd/templates/kernel/InsertParticleIntoLinkedCells.templ.h

@@ -79,7 +79,7 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
    {
       WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
       {%- for dim in range(3) %}
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
       {%- endfor %}
       {%- for dim in range(3) %}

python/mesa_pd/templates/kernel/InsertParticleIntoSparseLinkedCells.templ.h

@@ -79,7 +79,7 @@ inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx,
    {
       WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
       {%- for dim in range(3) %}
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash{{dim}} = static_cast<int>(std::floor((ac.getPosition(p_idx)[{{dim}}] - minCorner[{{dim}}]) * lc.invCellDiameter_[{{dim}}]));
       {%- endfor %}
       {%- for dim in range(3) %}

python/mesa_pd/templates/kernel/VelocityVerlet.templ.h

@@ -36,7 +36,7 @@ namespace kernel {
 /**
  * Velocity verlet integration for all particles.
  *
- * Velocit verlet integration is a two part kernel. preForceUpdate has to be
+ * Velocity verlet integration is a two part kernel. preForceUpdate has to be
  * called before the force calculation and postFroceUpdate afterwards. The
  * integration is only complete when both functions are called. The integration
  * is symplectic.

python/mesa_pd/templates/mpi/notifications/ContactHistoryNotification.templ.h

@@ -41,7 +41,7 @@ namespace walberla {
 namespace mesa_pd {
 
 /**
- * Trasmits the contact history
+ * Transmits the contact history
  */
 class ContactHistoryNotification
 {

python/mesa_pd/templates/mpi/notifications/PropertyNotification.templ.h

@@ -39,7 +39,7 @@ namespace walberla {
 namespace mesa_pd {
 
 /**
- * Trasmits force and torque information.
+ * Transmits force and torque information.
  */
 class {{name}}
 {

src/mesa_pd/kernel/InsertParticleIntoLinkedCells.h

@@ -73,11 +73,11 @@ inline void InsertParticleIntoLinkedCells::operator()(const size_t p_idx, Access
    } else
    {
       WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
       if (hash0 < 0) hash0 = 0;
       if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;

src/mesa_pd/kernel/InsertParticleIntoSparseLinkedCells.h

@@ -73,11 +73,11 @@ inline void InsertParticleIntoSparseLinkedCells::operator()(const size_t p_idx,
    } else
    {
       WALBERLA_ASSERT_GREATER(ac.getInteractionRadius(p_idx), 0_r, "Did you forget to set the interaction radius?");
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash0 = static_cast<int>(std::floor((ac.getPosition(p_idx)[0] - minCorner[0]) * lc.invCellDiameter_[0]));
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash1 = static_cast<int>(std::floor((ac.getPosition(p_idx)[1] - minCorner[1]) * lc.invCellDiameter_[1]));
-      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is to large for this cell size. Contacts might get lost.");
+      WALBERLA_ASSERT_LESS(2_r * ac.getInteractionRadius(p_idx), lc.cellDiameter_[0], "Interaction radius is too large for this cell size. Contacts might get lost.");
       int hash2 = static_cast<int>(std::floor((ac.getPosition(p_idx)[2] - minCorner[2]) * lc.invCellDiameter_[2]));
       if (hash0 < 0) hash0 = 0;
       if (hash0 >= lc.numCellsPerDim_[0]) hash0 = lc.numCellsPerDim_[0] - 1;

src/pe/ccd/HashGrids.h

@@ -81,7 +81,7 @@ namespace ccd {
  * bodies.
  *
  * For further information and a much more detailed explanation of this algorithm see
- *     http://www10.informatik.uni-erlangen.de/Publications/Theses/2009/Schornbaum_SA09.pdf
+ *     https://www10.cs.fau.de/publications/theses/2009/Schornbaum_SA_2009.pdf
  */
 class HashGrids : public ICCD
 {

tests/mesa_pd/CMakeLists.txt

@@ -43,6 +43,9 @@ waLBerla_execute_test( NAME   MESA_PD_Data_ContactHistory )
 waLBerla_compile_test( NAME   MESA_PD_Data_Flags FILES data/Flags.cpp DEPENDS core )
 waLBerla_execute_test( NAME   MESA_PD_Data_Flags )
 
+waLBerla_compile_test( NAME   MESA_PD_Data_HashGridsVsBruteForce FILES data/HashGridsVsBruteForce.cpp DEPENDS core )
+waLBerla_execute_test( NAME   MESA_PD_Data_HashGridsVsBruteForce PROCESSES 27 )
+
 waLBerla_compile_test( NAME   MESA_PD_Data_LinkedCells FILES data/LinkedCells.cpp DEPENDS core )
 waLBerla_execute_test( NAME   MESA_PD_Data_LinkedCells )
 

tests/mesa_pd/data/HashGridsVsBruteForce.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 Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include <mesa_pd/data/HashGrids.h>
+#include <mesa_pd/data/ParticleAccessor.h>
+#include <mesa_pd/data/ParticleStorage.h>
+#include <mesa_pd/domain/BlockForestDomain.h>
+#include <mesa_pd/kernel/ParticleSelector.h>
+#include <mesa_pd/mpi/SyncNextNeighbors.h>
+
+#include <blockforest/BlockForest.h>
+#include <blockforest/Initialization.h>
+#include <core/Environment.h>
+#include <core/grid_generator/SCIterator.h>
+#include <core/logging/Logging.h>
+#include <core/math/Random.h>
+
+#include <algorithm>
+#include <iostream>
+#include <memory>
+#include <numeric>
+
+namespace walberla {
+namespace mesa_pd {
+
+using IdxPair_T = std::pair<size_t, size_t>;
+
+// unique sorting of index pairs
+bool compPair(const IdxPair_T& a, const IdxPair_T& b)
+{
+   if (a.first == b.first) return a.second < b.second;
+   return a.first < b.first;
+}
+
+// check contact between two particles with position and radius
+// mimics sphere-sphere contact detection
+bool areOverlapping(Vec3 pos1, real_t radius1, Vec3 pos2, real_t radius2)
+{
+   return (pos2-pos1).length() < radius1 + radius2;
+}
+
+void checkTestScenario( real_t radiusRatio )
+{
+   const real_t radius  = real_t(0.5);
+   const real_t radiusMax = std::sqrt(radiusRatio) * radius;
+   const real_t radiusMin = radius / std::sqrt(radiusRatio);
+
+   math::seedRandomGenerator( numeric_cast<std::mt19937::result_type>( 42 * walberla::mpi::MPIManager::instance()->rank() ) );
+
+   //logging::Logging::instance()->setStreamLogLevel(logging::Logging::DETAIL);
+   //logging::Logging::instance()->includeLoggingToFile("MESA_PD_Data_HashGridsVsBruteForce");
+   //logging::Logging::instance()->setFileLogLevel(logging::Logging::DETAIL);
+
+   //init domain partitioning
+   auto forest = blockforest::createBlockForest( AABB(0,0,0,30,30,30), // simulation domain
+                                                 Vector3<uint_t>(3,3,3), // blocks in each direction
+                                                 Vector3<bool>(true, true, true) // periodicity
+                                                 );
+   domain::BlockForestDomain domain(forest);
+
+   WALBERLA_CHECK_EQUAL(forest->size(), 1);
+   auto& blk = *forest->begin();
+
+   WALBERLA_CHECK(blk.getAABB().xSize() > radiusMax && blk.getAABB().ySize() > radiusMax &&  blk.getAABB().zSize() > radiusMax,
+                  "Condition for next neighbor sync violated!" );
+
+   //init data structures
+   auto ps = std::make_shared<data::ParticleStorage>(100);
+   data::HashGrids hg;
+   std::vector<IdxPair_T> csBF;
+   std::vector<IdxPair_T> csHG1;
+   std::vector<IdxPair_T> csHG2;
+   std::vector<IdxPair_T> csHG3;
+
+   data::ParticleAccessor accessor(ps);
+
+   //initialize particles
+   for (int i = 0; i < 1000; ++i)
+   {
+      data::Particle&& p = *ps->create();
+      p.getPositionRef() = Vec3( math::realRandom(blk.getAABB().xMin(), blk.getAABB().xMax()),
+                                 math::realRandom(blk.getAABB().yMin(), blk.getAABB().yMax()),
+                                 math::realRandom(blk.getAABB().zMin(), blk.getAABB().zMax()) );
+
+      p.getOwnerRef() = walberla::mpi::MPIManager::instance()->rank();
+      p.getInteractionRadiusRef() = math::realRandom(radiusMin, radiusMax);
+   }
+
+   //init kernels
+   mpi::SyncNextNeighbors SNN;
+
+   SNN(*ps, domain);
+
+   ps->forEachParticlePairHalf(false,
+                               kernel::SelectAll(),
+                               accessor,
+                               [&csBF](const size_t idx1, const size_t idx2, auto& ac)
+   {
+      if (areOverlapping(ac.getPosition(idx1), ac.getInteractionRadius(idx1),
+                         ac.getPosition(idx2), ac.getInteractionRadius(idx2) ))
+      {
+         if(ac.getUid(idx2) < ac.getUid(idx1))
+            csBF.push_back(IdxPair_T(idx2, idx1));
+         else
+            csBF.push_back(IdxPair_T(idx1, idx2));
+      }
+   },
+   accessor );
+
+   // insert into hash grids initially
+
+   ps->forEachParticle(true, kernel::SelectAll(), accessor, hg, accessor);
+   hg.forEachParticlePairHalf(false,
+                              kernel::SelectAll(),
+                              accessor,
+                              [&csHG1](const size_t idx1, const size_t idx2, auto& ac)
+   {
+      if (areOverlapping(ac.getPosition(idx1), ac.getInteractionRadius(idx1),
+                         ac.getPosition(idx2), ac.getInteractionRadius(idx2) ))
+      {
+         if(ac.getUid(idx2) < ac.getUid(idx1))
+            csHG1.push_back(IdxPair_T(idx2, idx1));
+         else
+            csHG1.push_back(IdxPair_T(idx1, idx2));
+      }
+   },
+   accessor );
+
+   WALBERLA_CHECK_EQUAL(csBF.size(), csHG1.size());
+   WALBERLA_LOG_DEVEL(csBF.size() << " contacts detected");
+
+   std::sort(csBF.begin(), csBF.end(), compPair);
+   std::sort(csHG1.begin(), csHG1.end(), compPair);
+
+   for (size_t i = 0; i < csBF.size(); ++i)
+   {
+      WALBERLA_CHECK_EQUAL(csBF[i].first, csHG1[i].first);
+      WALBERLA_CHECK_EQUAL(csBF[i].second, csHG1[i].second);
+   }
+
+   WALBERLA_LOG_DEVEL_ON_ROOT("Initial insertion checked");
+   WALBERLA_MPI_BARRIER();
+
+   // redo to check clear
+   hg.clear();
+   ps->forEachParticle(true, kernel::SelectAll(), accessor, hg, accessor);
+   hg.forEachParticlePairHalf(false,
+                              kernel::SelectAll(),
+                              accessor,
+                              [&csHG2](const size_t idx1, const size_t idx2, auto& ac)
+                              {
+                                 if (areOverlapping(ac.getPosition(idx1), ac.getInteractionRadius(idx1),
+                                                    ac.getPosition(idx2), ac.getInteractionRadius(idx2) ))
+                                 {
+                                    if(ac.getUid(idx2) < ac.getUid(idx1))
+                                       csHG2.push_back(IdxPair_T(idx2, idx1));
+                                    else
+                                       csHG2.push_back(IdxPair_T(idx1, idx2));
+                                 }
+                              },
+                              accessor );
+
+   WALBERLA_CHECK_EQUAL(csBF.size(), csHG2.size());
+
+   std::sort(csHG2.begin(), csHG2.end(), compPair);
+
+   for (size_t i = 0; i < csBF.size(); ++i)
+   {
+      WALBERLA_CHECK_EQUAL(csBF[i].first, csHG2[i].first);
+      WALBERLA_CHECK_EQUAL(csBF[i].second, csHG2[i].second);
+   }
+   WALBERLA_LOG_DEVEL_ON_ROOT("Insertion after clear checked");
+   WALBERLA_MPI_BARRIER();
+
+   // redo to check clearAll
+   hg.clearAll();
+   ps->forEachParticle(true, kernel::SelectAll(), accessor, hg, accessor);
+   hg.forEachParticlePairHalf(false,
+                              kernel::SelectAll(),
+                              accessor,
+                              [&csHG3](const size_t idx1, const size_t idx2, auto& ac)
+                              {
+                                 if (areOverlapping(ac.getPosition(idx1), ac.getInteractionRadius(idx1),
+                                                    ac.getPosition(idx2), ac.getInteractionRadius(idx2) ))
+                                 {
+                                    if(ac.getUid(idx2) < ac.getUid(idx1))
+                                       csHG3.push_back(IdxPair_T(idx2, idx1));
+                                    else
+                                       csHG3.push_back(IdxPair_T(idx1, idx2));
+                                 }
+                              },
+                              accessor );
+
+   WALBERLA_CHECK_EQUAL(csBF.size(), csHG3.size());
+
+   std::sort(csHG3.begin(), csHG3.end(), compPair);
+
+   for (size_t i = 0; i < csBF.size(); ++i)
+   {
+      WALBERLA_CHECK_EQUAL(csBF[i].first, csHG3[i].first);
+      WALBERLA_CHECK_EQUAL(csBF[i].second, csHG3[i].second);
+   }
+
+   WALBERLA_LOG_DEVEL_ON_ROOT("Insertion after clear all checked");
+   WALBERLA_MPI_BARRIER();
+
+
+}
+
+/*
+ * Generates particles randomly inside the domain.
+ * Then checks if the Hash Grids find the same interaction pairs as the naive all-against-all check.
+ * Similar to test in "LinkedCellsVsBruteForce.cpp"
+ */
+int main( int argc, char ** argv ) {
+   Environment env(argc, argv);
+   WALBERLA_UNUSED(env);
+   walberla::mpi::MPIManager::instance()->useWorldComm();
+
+   WALBERLA_LOG_DEVEL_ON_ROOT("Checking monodisperse case");
+   checkTestScenario(walberla::real_t(1.01)); // monodisperse
+
+   WALBERLA_LOG_DEVEL_ON_ROOT("Checking polydisperse case");
+   checkTestScenario(walberla::real_t(10)); // polydisperse
+
+   return EXIT_SUCCESS;
+}
+
+} //namespace mesa_pd
+} //namespace walberla
+
+int main( int argc, char ** argv )
+{
+   return walberla::mesa_pd::main(argc, argv);
+}