added particle sorting to ParticleStorage

python/mesa_pd/templates/data/ParticleStorage.templ.h

@@ -133,7 +133,7 @@ public:
 
    iterator begin() { return iterator(this, 0); }
    iterator end()   { return iterator(this, size()); }
-   iterator operator[](const size_t n) { return iterator(this, n); }
+   Particle operator[](const size_t n) { return *iterator(this, n); }
 
    {% for prop in properties %}
    const {{prop.type}}& get{{prop.name | capFirst}}(const size_t idx) const {return {{prop.name}}_[idx];}
@@ -159,6 +159,8 @@ public:
    inline void reserve(const size_t size);
    inline void clear();
    inline size_t size() const;
+   template <class Compare>
+   void sort(Compare comp);
 
    /**
     * Calls the provided functor \p func for all Particles.
@@ -256,6 +258,15 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
    return *this;
 }
 
+inline
+void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
+{
+   if (lhs.i_ == rhs.i_) return;
+   {%- for prop in properties %}
+   std::swap(lhs.get{{prop.name | capFirst}}Ref(), rhs.get{{prop.name | capFirst}}Ref());
+   {%- endfor %}
+}
+
 inline
 std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
 {
@@ -392,7 +403,7 @@ inline ParticleStorage::iterator ParticleStorage::find(const id_t& uid)
    //use unordered_map for faster lookup
    auto it = uidToIdx_.find(uid);
    if (it == uidToIdx_.end()) return end();
-   WALBERLA_ASSERT_EQUAL(it->first, uid, "Lookup via uidToIdx map is not up to date!!!");
+   WALBERLA_ASSERT_EQUAL(getUid(it->second), uid, "Lookup via uidToIdx map is not up to date!!!");
    return iterator(this, it->second);
 }
 
@@ -419,6 +430,47 @@ inline size_t ParticleStorage::size() const
    return {{properties[0].name}}_.size();
 }
 
+template <class Compare>
+void ParticleStorage::sort(Compare comp)
+{
+   using WeightPair = std::pair<size_t, double>; //idx, weight
+
+   const size_t length = size();
+   std::vector<size_t>     newIdx(length); //where is old idx now?
+   std::vector<size_t>     oldIdx(length); //what old idx is at that idx?
+   std::vector<WeightPair> weight(length);
+
+   for (size_t idx = 0; idx < length; ++idx)
+   {
+      newIdx[idx] = idx;
+      oldIdx[idx] = idx;
+      weight[idx] = std::make_pair(idx, comp.getWeight(operator[](idx)));
+   }
+   std::sort(weight.begin(), weight.end(), [](const WeightPair& lhs, const WeightPair& rhs){return lhs.second < rhs.second;});
+   for (size_t idx = 0; idx < length; ++idx)
+   {
+      using std::swap;
+      WALBERLA_ASSERT_IDENTICAL(weight[idx].second, comp.getWeight(operator[](newIdx[weight[idx].first])));
+
+      WALBERLA_ASSERT_LESS_EQUAL(comp.getWeight(operator[](newIdx[weight[idx].first])), comp.getWeight(operator[](idx)));
+      swap( operator[](idx), operator[](newIdx[weight[idx].first]) );
+      const auto lhsIDX = idx;
+      const auto rhsIDX = newIdx[weight[idx].first];
+
+      newIdx[oldIdx[lhsIDX]] = rhsIDX;
+      newIdx[oldIdx[rhsIDX]] = lhsIDX;
+
+      swap(oldIdx[lhsIDX], oldIdx[rhsIDX]);
+   }
+
+   //rebuild lookup table
+   uidToIdx_.clear();
+   for (size_t idx = 0; idx < length; ++idx)
+   {
+      uidToIdx_[getUid(idx)] = idx;
+   }
+}
+
 {%- for const in ["", "const"] %}
 template <typename Selector, typename Accessor, typename Func, typename... Args>
 inline void ParticleStorage::forEachParticle(const bool openmp,

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

@@ -29,6 +29,7 @@
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/mpi/notifications/NotificationType.h>
+#include <mesa_pd/mpi/notifications/reset.h>
 
 #include <core/mpi/Datatype.h>
 #include <core/mpi/RecvBuffer.h>
@@ -54,6 +55,12 @@ public:
    const data::Particle& p_;
 };
 
+template <>
+void reset<HeatFluxNotification>(data::Particle& p)
+{
+   p.setHeatFlux( real_t(0) );
+}
+
 void reduce(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
 {
    p.getHeatFluxRef()  += objparam.heatFlux_;

src/mesa_pd/data/ParticleStorage.h

@@ -223,7 +223,7 @@ public:
 
    iterator begin() { return iterator(this, 0); }
    iterator end()   { return iterator(this, size()); }
-   iterator operator[](const size_t n) { return iterator(this, n); }
+   Particle operator[](const size_t n) { return *iterator(this, n); }
 
    
    const walberla::id_t& getUid(const size_t idx) const {return uid_[idx];}
@@ -341,6 +341,8 @@ public:
    inline void reserve(const size_t size);
    inline void clear();
    inline size_t size() const;
+   template <class Compare>
+   void sort(Compare comp);
 
    /**
     * Calls the provided functor \p func for all Particles.
@@ -501,6 +503,33 @@ ParticleStorage::Particle& ParticleStorage::Particle::operator=(ParticleStorage:
    return *this;
 }
 
+inline
+void swap(ParticleStorage::Particle lhs, ParticleStorage::Particle rhs)
+{
+   if (lhs.i_ == rhs.i_) return;
+   std::swap(lhs.getUidRef(), rhs.getUidRef());
+   std::swap(lhs.getPositionRef(), rhs.getPositionRef());
+   std::swap(lhs.getInteractionRadiusRef(), rhs.getInteractionRadiusRef());
+   std::swap(lhs.getFlagsRef(), rhs.getFlagsRef());
+   std::swap(lhs.getOwnerRef(), rhs.getOwnerRef());
+   std::swap(lhs.getGhostOwnersRef(), rhs.getGhostOwnersRef());
+   std::swap(lhs.getShapeIDRef(), rhs.getShapeIDRef());
+   std::swap(lhs.getRotationRef(), rhs.getRotationRef());
+   std::swap(lhs.getAngularVelocityRef(), rhs.getAngularVelocityRef());
+   std::swap(lhs.getTorqueRef(), rhs.getTorqueRef());
+   std::swap(lhs.getLinearVelocityRef(), rhs.getLinearVelocityRef());
+   std::swap(lhs.getInvMassRef(), rhs.getInvMassRef());
+   std::swap(lhs.getForceRef(), rhs.getForceRef());
+   std::swap(lhs.getOldForceRef(), rhs.getOldForceRef());
+   std::swap(lhs.getOldTorqueRef(), rhs.getOldTorqueRef());
+   std::swap(lhs.getTypeRef(), rhs.getTypeRef());
+   std::swap(lhs.getNextParticleRef(), rhs.getNextParticleRef());
+   std::swap(lhs.getOldContactHistoryRef(), rhs.getOldContactHistoryRef());
+   std::swap(lhs.getNewContactHistoryRef(), rhs.getNewContactHistoryRef());
+   std::swap(lhs.getTemperatureRef(), rhs.getTemperatureRef());
+   std::swap(lhs.getHeatFluxRef(), rhs.getHeatFluxRef());
+}
+
 inline
 std::ostream& operator<<( std::ostream& os, const ParticleStorage::Particle& p )
 {
@@ -700,7 +729,7 @@ inline ParticleStorage::iterator ParticleStorage::find(const id_t& uid)
    //use unordered_map for faster lookup
    auto it = uidToIdx_.find(uid);
    if (it == uidToIdx_.end()) return end();
-   WALBERLA_ASSERT_EQUAL(it->first, uid, "Lookup via uidToIdx map is not up to date!!!");
+   WALBERLA_ASSERT_EQUAL(getUid(it->second), uid, "Lookup via uidToIdx map is not up to date!!!");
    return iterator(this, it->second);
 }
 
@@ -789,6 +818,47 @@ inline size_t ParticleStorage::size() const
    //WALBERLA_ASSERT_EQUAL( uid_.size(), neighborState.size() );
    return uid_.size();
 }
+
+template <class Compare>
+void ParticleStorage::sort(Compare comp)
+{
+   using WeightPair = std::pair<size_t, double>; //idx, weight
+
+   const size_t length = size();
+   std::vector<size_t>     newIdx(length); //where is old idx now?
+   std::vector<size_t>     oldIdx(length); //what old idx is at that idx?
+   std::vector<WeightPair> weight(length);
+
+   for (size_t idx = 0; idx < length; ++idx)
+   {
+      newIdx[idx] = idx;
+      oldIdx[idx] = idx;
+      weight[idx] = std::make_pair(idx, comp.getWeight(operator[](idx)));
+   }
+   std::sort(weight.begin(), weight.end(), [](const WeightPair& lhs, const WeightPair& rhs){return lhs.second < rhs.second;});
+   for (size_t idx = 0; idx < length; ++idx)
+   {
+      using std::swap;
+      WALBERLA_ASSERT_IDENTICAL(weight[idx].second, comp.getWeight(operator[](newIdx[weight[idx].first])));
+
+      WALBERLA_ASSERT_LESS_EQUAL(comp.getWeight(operator[](newIdx[weight[idx].first])), comp.getWeight(operator[](idx)));
+      swap( operator[](idx), operator[](newIdx[weight[idx].first]) );
+      const auto lhsIDX = idx;
+      const auto rhsIDX = newIdx[weight[idx].first];
+
+      newIdx[oldIdx[lhsIDX]] = rhsIDX;
+      newIdx[oldIdx[rhsIDX]] = lhsIDX;
+
+      swap(oldIdx[lhsIDX], oldIdx[rhsIDX]);
+   }
+
+   //rebuild lookup table
+   uidToIdx_.clear();
+   for (size_t idx = 0; idx < length; ++idx)
+   {
+      uidToIdx_[getUid(idx)] = idx;
+   }
+}
 template <typename Selector, typename Accessor, typename Func, typename... Args>
 inline void ParticleStorage::forEachParticle(const bool openmp,
                                              const Selector& selector,

src/mesa_pd/mpi/notifications/HeatFluxNotification.h

@@ -29,6 +29,7 @@
 #include <mesa_pd/data/DataTypes.h>
 #include <mesa_pd/data/ParticleStorage.h>
 #include <mesa_pd/mpi/notifications/NotificationType.h>
+#include <mesa_pd/mpi/notifications/reset.h>
 
 #include <core/mpi/Datatype.h>
 #include <core/mpi/RecvBuffer.h>
@@ -54,6 +55,12 @@ public:
    const data::Particle& p_;
 };
 
+template <>
+void reset<HeatFluxNotification>(data::Particle& p)
+{
+   p.setHeatFlux( real_t(0) );
+}
+
 void reduce(data::Particle&& p, const HeatFluxNotification::Parameters& objparam)
 {
    p.getHeatFluxRef()  += objparam.heatFlux_;

src/mesa_pd/sorting/HilbertCompareFunctor.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 HilbertCompareFunctor.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "HilbertCompareFunctor.h"
+
+#include <blockforest/HilbertCurveConstruction.h>
+#include <core/DataTypes.h>
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+#include <core/math/Vector3.h>
+#include <core/mpi/MPIManager.h>
+
+#include <stack>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+math::Vector3<uint_t> getChildPosition(const uint_t currentSize, const math::Vector3<uint_t>& parent, const uint_t childId)
+{
+   auto size = currentSize >> 1;
+   auto child = parent;
+
+   if (!(childId & 0b001))
+   {
+      child[0] -= size;
+   }
+
+   if (!(childId & 0b010))
+   {
+      child[1] -= size;
+   }
+
+   if (!(childId & 0b100))
+   {
+      child[2] -= size;
+   }
+
+   return child;
+}
+
+std::vector< math::Vector3<uint_t> > generateHilbertCurve(const uint_t size)
+{
+   using Node = std::pair<uint_t, Vector3<uint_t>>; //level, coordinate
+
+   std::vector< math::Vector3<uint_t> > result;
+
+   std::stack< Node > stack;
+   std::stack< uint_t > orientation;
+
+   stack.push( std::make_pair(size, math::Vector3<uint_t>(size)) );
+   orientation.push( uint_t(0) );
+
+   while( !stack.empty() )
+   {
+      const Node node = stack.top();
+      uint_t hilbertIndex = orientation.top();
+
+      stack.pop();
+      orientation.pop();
+
+      if( node.first > 1 )
+      {
+         for( uint_t c = 8; c-- != 0; )
+         {
+            stack.push( std::make_pair( node.first >> 1,
+                                        getChildPosition(node.first, node.second, blockforest::hilbertOrder[hilbertIndex][c]) ));
+            orientation.push( blockforest::hilbertOrientation[hilbertIndex][c] );
+         }
+      }
+      else
+      {
+         result.push_back( node.second );
+      }
+   }
+   return result;
+}
+
+
+HilbertCompareFunctor::HilbertCompareFunctor(const math::AABB& domain, const uint_t cells)
+   : domain_(domain)
+   , cells_(cells)
+   , hilbertLookup_(cells_ * cells_ * cells_)
+{
+   WALBERLA_CHECK(math::uintIsPowerOfTwo(cells));
+   inverse_dx[0] = real_t(1.0) / (domain_.xSize() / real_c(cells_));
+   inverse_dx[1] = real_t(1.0) / (domain_.ySize() / real_c(cells_));
+   inverse_dx[2] = real_t(1.0) / (domain_.zSize() / real_c(cells_));
+   initializeLookup();
+}
+
+bool HilbertCompareFunctor::operator()(const data::Particle p1, const data::Particle p2) const
+{
+   const auto hash1 = discretize(p1.getPosition() - domain_.minCorner());
+   WALBERLA_ASSERT_LESS( hash1, cells_*cells_*cells_);
+   const auto hash2 = discretize(p2.getPosition() - domain_.minCorner());
+   WALBERLA_ASSERT_LESS( hash2, cells_*cells_*cells_);
+
+   return hilbertLookup_[hash1] < hilbertLookup_[hash2];
+}
+
+uint_t HilbertCompareFunctor::discretize(const Vec3& pos) const
+{
+   const real_t& x = pos[0];
+   const real_t& y = pos[1];
+   const real_t& z = pos[2];
+
+   const uint_t hashMask = (cells_ - 1);
+
+   math::Vector3<uint_t> hash;
+
+   if( x < 0 ) {
+      real_t i = ( -x ) * inverse_dx[0];
+      hash[0]  = cells_ - 1 - ( static_cast<uint_t>( i ) & hashMask );
+   }
+   else {
+      real_t i = x * inverse_dx[0];
+      hash[0]  = static_cast<uint_t>( i ) & hashMask;
+   }
+
+   if( y < 0 ) {
+      real_t i = ( -y ) * inverse_dx[1];
+      hash[1]  = cells_ - 1 - ( static_cast<uint_t>( i ) & hashMask );
+   }
+   else {
+      real_t i = y * inverse_dx[1];
+      hash[1]  = static_cast<uint_t>( i ) & hashMask;
+   }
+
+   if( z < 0 ) {
+      real_t i = ( -z ) * inverse_dx[2];
+      hash[2]  = cells_ - 1 - ( static_cast<uint_t>( i ) & hashMask );
+   }
+   else {
+      real_t i = z * inverse_dx[2];
+      hash[2]  = static_cast<uint_t>( i ) & hashMask;
+   }
+
+   WALBERLA_ASSERT_LESS(hash[0], cells_);
+   WALBERLA_ASSERT_LESS(hash[1], cells_);
+   WALBERLA_ASSERT_LESS(hash[2], cells_);
+
+   return hash[2] * cells_ * cells_ + hash[1] * cells_ + hash[0];
+}
+
+
+void HilbertCompareFunctor::initializeLookup()
+{
+   const auto sfc = generateHilbertCurve(cells_);
+   uint_t counter = 0;
+   for (auto p : sfc)
+   {
+      WALBERLA_ASSERT_GREATER(p[0], 0);
+      WALBERLA_ASSERT_GREATER(p[1], 0);
+      WALBERLA_ASSERT_GREATER(p[2], 0);
+      p[0] -= 1;
+      p[1] -= 1;
+      p[2] -= 1;
+
+      hilbertLookup_[p[2] * cells_ * cells_ + p[1] * cells_ + p[0]] = counter;
+      ++counter;
+   }
+   WALBERLA_ASSERT_EQUAL(counter, sfc.size());
+   WALBERLA_ASSERT_EQUAL(counter, hilbertLookup_.size());
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla

src/mesa_pd/sorting/HilbertCompareFunctor.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 HilbertCompareFunctor.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <core/math/AABB.h>
+#include <core/math/Vector3.h>
+
+#include <mesa_pd/data/DataTypes.h>
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+/**
+ * Compare functor which sorts particles according to the hilbert space filling curve.
+ */
+class HilbertCompareFunctor
+{
+public:
+   /**
+    * Subdivides a domain into cells and arranges the cells according to the sfc.
+    * The position within the sfc gives the particles their weight.
+    * \param domain domain which should be partitioned into cells
+    * \param cells number of cells in every direction.
+    * \attention cells has to be a power of 2!
+    */
+   HilbertCompareFunctor(const math::AABB& domain, const uint_t cells);
+   bool operator()(const data::Particle bd1, const data::Particle bd2) const;
+   double getWeight(const data::Particle p1) const;
+private:
+   uint_t discretize(const Vec3& pos) const;
+   void initializeLookup();
+
+   math::AABB  domain_;
+   const uint_t cells_;
+   Vec3 inverse_dx;
+   std::vector<uint_t> hilbertLookup_;
+};
+
+inline
+double HilbertCompareFunctor::getWeight(const data::Particle p1) const
+{
+   const auto hash1 = discretize(p1.getPosition() - domain_.minCorner());
+   WALBERLA_ASSERT_LESS( hash1, cells_*cells_*cells_);
+
+   return double_c(hilbertLookup_[hash1]);
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla

src/mesa_pd/sorting/LinearizedCompareFunctor.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 LinearizedCompareFunctor.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "LinearizedCompareFunctor.h"
+
+#include <core/DataTypes.h>
+#include <core/Environment.h>
+#include <core/logging/Logging.h>
+#include <core/math/Vector3.h>
+#include <core/mpi/MPIManager.h>
+
+#include <stack>
+#include <vector>
+
+namespace walberla {
+namespace mesa_pd {
+namespace sorting {
+
+LinearizedCompareFunctor::LinearizedCompareFunctor(const math::AABB& domain, const Vector3<uint_t> cells)
+   : domain_(domain)
+   , cells_(cells)
+{
+   inverse_dx[0] = real_c(cells_[0]) / domain_.xSize();
+   inverse_dx[1] = real_c(cells_[1]) / domain_.ySize();
+   inverse_dx[2] = real_c(cells_[2]) / domain_.zSize();
+}
+
+bool LinearizedCompareFunctor::operator()(const data::Particle p1, const data::Particle p2) const
+{
+   const auto hash1 = discretize(p1.getPosition() - domain_.minCorner());
+   WALBERLA_ASSERT_LESS( hash1, cells_[0]*cells_[1]*cells_[2]);
+   const auto hash2 = discretize(p2.getPosition() - domain_.minCorner());
+   WALBERLA_ASSERT_LESS( hash2, cells_[0]*cells_[1]*cells_[2]);
+
+   return hash1 < hash2;
+}
+
+uint_t LinearizedCompareFunctor::discretize(const Vec3& pos) const
+{
+   int x = int_c(pos[0] * inverse_dx[0]);
+   int y = int_c(pos[1] * inverse_dx[1]);
+   int z = int_c(pos[2] * inverse_dx[2]);
+   if (x<0) x=0;
+   if (y<0) y=0;
+   if (z<0) z=0;
+   if (x>=int_c(cells_[0])) x=int_c(cells_[0])-1;
+   if (y>=int_c(cells_[1])) y=int_c(cells_[1])-1;
+   if (z>=int_c(cells_[2])) z=int_c(cells_[2])-1;
+
+   return uint_c(z) * cells_[1] * cells_[0] + uint_c(y) * cells_[0] + uint_c(x);
+}
+
+} //namespace sorting
+} //namespace mesa_pd
+} //namespace walberla

src/mesa_pd/sorting/LinearizedCompareFunctor.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/>.
+//