Merge branch 'lbm_mesapd_coupling_amr' into 'master'

Adaptive grid refinement for LBM and MESA-PD

See merge request walberla/walberla!349

src/lbm/refinement/CurlBasedLevelDetermination.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 VorticityBasedLevelDetermination.h
+//! \ingroup lbm
+//! \author Florian Schornbaum <florian.schornbaum@fau.de>
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/BlockForest.h"
+#include "core/math/Vector3.h"
+#include "domain_decomposition/BlockDataID.h"
+#include "field/GhostLayerField.h"
+
+#include <vector>
+
+namespace walberla {
+namespace lbm {
+namespace refinement {
+
+
+/*!\brief Level determination for refinement check based on local curl
+ *
+ * If (scaled) vorticity magnitude is below lowerLimit in all cells of a block, that block could be coarsened.
+ * If the (scaled) vorticity value is above the upperLimit for at least one cell, that block gets marked for refinement.
+ * Else, the block remains on the current level.
+ *
+ * The scaling originates from neglecting the actual mesh size on the block to obtain different vorticity values for
+ * different mesh sizes.
+ *
+ * Parametes upperLimit corresponds to sigma_c, coarsenFactor to c, lowerLimit to c*sigma_c, lengthScaleWeight to r.
+ */
+template< typename Filter_T >
+class CurlBasedLevelDetermination // used as a 'BlockForest::RefreshMinTargetLevelDeterminationFunction'
+{
+
+public:
+
+   typedef GhostLayerField< Vector3<real_t>, 1 >  VectorField_T;
+
+   CurlBasedLevelDetermination(const ConstBlockDataID & fieldID, const StructuredBlockForest & structuredBlockForest,
+         const Filter_T & filter, const uint_t maxLevel,
+         const real_t upperLimit, const real_t lowerLimit, const real_t lengthScaleWeight = real_t(2)) :
+         fieldID_(fieldID), structuredBlockForest_(structuredBlockForest), filter_(filter), maxLevel_(maxLevel),
+         upperLimitSqr_(upperLimit*upperLimit), lowerLimitSqr_(lowerLimit*lowerLimit),
+         lengthScaleWeight_(lengthScaleWeight)
+   {
+      WALBERLA_CHECK_FLOAT_UNEQUAL(lengthScaleWeight_, real_t(0)); // else std::pow(x, y/0) is calculated further below
+   }
+
+   void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                    std::vector< const Block * > & blocksAlreadyMarkedForRefinement,
+                    const BlockForest & forest );
+
+private:
+
+   ConstBlockDataID fieldID_;
+   const StructuredBlockForest & structuredBlockForest_;
+
+   Filter_T filter_;
+
+   uint_t maxLevel_;
+
+   real_t upperLimitSqr_;
+   real_t lowerLimitSqr_;
+
+   real_t lengthScaleWeight_;
+};
+
+template< typename Filter_T >
+void CurlBasedLevelDetermination< Filter_T >::operator()( std::vector< std::pair< const Block *,
+      uint_t > > & minTargetLevels, std::vector< const Block * > &, const BlockForest & forest) {
+
+   for(auto & minTargetLevel : minTargetLevels) {
+      const Block * const block = minTargetLevel.first;
+      const VectorField_T * u = block->template getData< VectorField_T >(fieldID_);
+
+      if(u == nullptr) {
+         minTargetLevel.second = uint_t(0);
+         continue;
+      }
+
+      WALBERLA_ASSERT_GREATER_EQUAL(u->nrOfGhostLayers(), uint_t(1));
+
+      CellInterval interval = u->xyzSize();
+      Cell expand(cell_idx_c(-1), cell_idx_c(-1), cell_idx_c(-1));
+      interval.expand(expand);
+
+      const auto one = cell_idx_t(1);
+
+      const real_t dx = structuredBlockForest_.dx(forest.getLevel(*block));
+      const real_t dy = structuredBlockForest_.dy(forest.getLevel(*block));
+      const real_t dz = structuredBlockForest_.dz(forest.getLevel(*block));
+
+      const auto halfInvDx = real_t(0.5) * real_t(1) / dx;
+      const auto halfInvDy = real_t(0.5) * real_t(1) / dy;
+      const auto halfInvDz = real_t(0.5) * real_t(1) / dz;
+
+      bool refine = false;
+      bool coarsen = true;
+
+      filter_(*block);
+
+      const cell_idx_t xSize = cell_idx_c(interval.xSize());
+      const cell_idx_t ySize = cell_idx_c(interval.ySize());
+      const cell_idx_t zSize = cell_idx_c(interval.zSize());
+
+      const real_t lengthScale = std::cbrt(dx*dy*dz);
+      const real_t weightedLengthScale = std::pow(lengthScale, (lengthScaleWeight_+1)/lengthScaleWeight_);
+      const real_t weightedLengthScaleSqr = weightedLengthScale*weightedLengthScale;
+
+      for (auto z = cell_idx_t(0); z < zSize; ++z) {
+         for (auto y = cell_idx_t(0); y < ySize; ++y) {
+            for (auto x = cell_idx_t(0); x < xSize; ++x) {
+               if (filter_(x,y,z) && filter_(x+one,y,z) && filter_(x-one,y,z) && filter_(x,y+one,z) && filter_(x,y-one,z)
+                  && filter_(x,y,z+one) && filter_(x,y,z-one)) {
+                  const Vector3< real_t > xa = u->get(x+one,y,z);
+                  const Vector3< real_t > xb = u->get(x-one,y,z);
+                  const Vector3< real_t > ya = u->get(x,y+one,z);
+                  const Vector3< real_t > yb = u->get(x,y-one,z);
+                  const Vector3< real_t > za = u->get(x,y,z+one);
+                  const Vector3< real_t > zb = u->get(x,y,z-one);
+
+                  const real_t duzdy = halfInvDy * (ya[2] - yb[2]);
+                  const real_t duydz = halfInvDz * (za[1] - zb[1]);
+                  const real_t duxdz = halfInvDz * (za[0] - zb[0]);
+                  const real_t duzdx = halfInvDx * (xa[2] - xb[2]);
+                  const real_t duydx = halfInvDx * (xa[1] - xb[1]);
+                  const real_t duxdy = halfInvDy * (ya[0] - yb[0]);
+
+                  const Vector3< real_t > curl( duzdy - duydz, duxdz - duzdx, duydx - duxdy );
+                  const auto curlSqr = curl.sqrLength();
+
+                  const auto curlSensorSqr = curlSqr * weightedLengthScaleSqr;
+
+                  if (curlSensorSqr > lowerLimitSqr_) {
+                     // curl is not small enough to coarsen, i.e. stay at least on the current level
+                     coarsen = false;
+                     if (curlSensorSqr > upperLimitSqr_) {
+                        // curl is big enough for refinement
+                        refine = true;
+                     }
+                  }
+               }
+            }
+         }
+      }
+
+      if (refine && block->getLevel() < maxLevel_) {
+         WALBERLA_ASSERT(!coarsen);
+         minTargetLevel.second = block->getLevel() + uint_t(1);
+      }
+      if (coarsen && block->getLevel() > uint_t(0)) {
+         WALBERLA_ASSERT(!refine);
+         minTargetLevel.second = block->getLevel() - uint_t(1);
+      }
+   }
+}
+
+} // namespace refinement
+} // namespace lbm
+} // namespace walberla

src/lbm/refinement/all.h

@@ -30,3 +30,4 @@
 #include "TimeStepPdfPackInfo.h"
 #include "TimeTracker.h"
 #include "VorticityBasedLevelDetermination.h"
+#include "CurlBasedLevelDetermination.h"
\ No newline at end of file

src/lbm_mesapd_coupling/amr/BlockInfo.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 BlockInfo.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//! \author Lukas Werner <lks.werner@fau.de>
+//! //
+//======================================================================================================================
+
+#pragma once
+
+#include <core/mpi/RecvBuffer.h>
+#include <core/mpi/SendBuffer.h>
+
+#include <ostream>
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+struct BlockInfo {
+   // lbm quantities
+   uint_t numberOfCells;
+   uint_t numberOfFluidCells;
+   uint_t numberOfNearBoundaryCells;
+   // pe quantities
+   uint_t numberOfLocalParticles;
+   uint_t numberOfShadowParticles;
+   uint_t numberOfContacts;
+   // coupling quantities
+   uint_t numberOfMESAPDSubCycles;
+
+   BlockInfo()
+         : numberOfCells(0), numberOfFluidCells(0), numberOfNearBoundaryCells(0),
+           numberOfLocalParticles(0), numberOfShadowParticles(0), numberOfContacts(0), numberOfMESAPDSubCycles(0) {}
+
+   BlockInfo(const uint_t numCells, const uint_t numFluidCells, const uint_t numNearBoundaryCells,
+             const uint_t numLocalBodies, const uint_t numShadowParticles, const uint_t numContacts,
+             const uint_t numPeSubCycles)
+         : numberOfCells(numCells), numberOfFluidCells(numFluidCells), numberOfNearBoundaryCells(numNearBoundaryCells),
+           numberOfLocalParticles(numLocalBodies), numberOfShadowParticles(numShadowParticles), numberOfContacts(numContacts), numberOfMESAPDSubCycles(numPeSubCycles) {}
+};
+
+
+inline
+std::ostream& operator<<( std::ostream& os, const BlockInfo& bi )
+{
+   os << bi.numberOfCells << " / " << bi.numberOfFluidCells << " / " << bi.numberOfNearBoundaryCells << " / "
+      << bi.numberOfLocalParticles << " / "<< bi.numberOfShadowParticles << " / " << bi.numberOfContacts << " / "
+      << bi.numberOfMESAPDSubCycles;
+   return os;
+}
+
+template< typename T,    // Element type of SendBuffer
+          typename G>    // Growth policy of SendBuffer
+mpi::GenericSendBuffer<T,G>& operator<<( mpi::GenericSendBuffer<T,G> & buf, const BlockInfo& info )
+{
+   buf.addDebugMarker( "pca" );
+   buf << info.numberOfCells << info.numberOfFluidCells << info.numberOfNearBoundaryCells
+       << info.numberOfLocalParticles << info.numberOfShadowParticles << info.numberOfContacts
+       << info.numberOfMESAPDSubCycles;
+   return buf;
+}
+
+template< typename T>    // Element type of SendBuffer
+mpi::GenericRecvBuffer<T>& operator>>( mpi::GenericRecvBuffer<T> & buf, BlockInfo& info )
+{
+   buf.readDebugMarker( "pca" );
+   buf >> info.numberOfCells >> info.numberOfFluidCells >> info.numberOfNearBoundaryCells
+       >> info.numberOfLocalParticles >> info.numberOfShadowParticles >> info.numberOfContacts
+       >> info.numberOfMESAPDSubCycles;
+   return buf;
+}
+
+} // namespace lbm_mesapd_coupling
+} // namespace walberla

src/lbm_mesapd_coupling/amr/InfoCollection.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 InfoCollection.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "BlockInfo.h"
+
+#include "blockforest/BlockForest.h"
+#include "blockforest/BlockID.h"
+#include "core/mpi/BufferSystem.h"
+
+#include <mesa_pd/data/ParticleStorage.h>
+
+#include <map>
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+typedef std::map<blockforest::BlockID, BlockInfo>  InfoCollection;
+typedef std::pair<blockforest::BlockID, BlockInfo> InfoCollectionPair;
+
+template <typename BoundaryHandling_T, typename ParticleAccessor_T>
+void createWithNeighborhood(BlockForest& bf, const BlockDataID boundaryHandlingID,
+                            const ParticleAccessor_T& accessor, const uint_t numberOfMESAPDSubCycles,
+                            InfoCollection& ic ) {
+   ic.clear();
+
+   mpi::BufferSystem bs( MPIManager::instance()->comm(), 856 );
+
+   for (auto blockIt = bf.begin(); blockIt != bf.end(); ++blockIt)
+   {
+      auto * block = static_cast<blockforest::Block*> (&(*blockIt));
+
+      // evaluate LBM quantities
+      BoundaryHandling_T * boundaryHandling = blockIt->getData< BoundaryHandling_T >( boundaryHandlingID );
+      auto xyzSize = boundaryHandling->getFlagField()->xyzSize();
+      const uint_t numCells = xyzSize.numCells();
+      uint_t numFluidCells(0), numNearBoundaryCells(0);
+      for( auto cellIt = xyzSize.begin(); cellIt != xyzSize.end(); ++cellIt)
+      {
+         if( boundaryHandling->isDomain(*cellIt) )
+         {
+            ++numFluidCells;
+         }
+         if( boundaryHandling->isNearBoundary(*cellIt))
+         {
+            ++numNearBoundaryCells;
+         }
+      }
+
+      // evaluate MESAPD quantities
+      // count block local and (possible) shadow particles here
+      uint_t numLocalParticles = 0, numShadowParticles = 0;
+      for (size_t idx = 0; idx < accessor.size(); ++idx) {
+         if (block->getAABB().contains(accessor.getPosition(idx))) {
+            // a particle within a block on the current process cannot be a ghost particle
+            WALBERLA_ASSERT(!mesa_pd::data::particle_flags::isSet(accessor.getFlags(idx), mesa_pd::data::particle_flags::GHOST));
+            numLocalParticles++;
+         } else if (block->getAABB().sqDistance(accessor.getPosition(idx)) < accessor.getInteractionRadius(idx)*accessor.getInteractionRadius(idx)) {
+            numShadowParticles++;
+         }
+      }
+
+      // count contacts here
+      const uint_t numContacts = 0;
+
+      BlockInfo blockInfo(numCells, numFluidCells, numNearBoundaryCells, numLocalParticles, numShadowParticles, numContacts, numberOfMESAPDSubCycles);
+      InfoCollectionPair infoCollectionEntry(block->getId(), blockInfo);
+
+      ic.insert( infoCollectionEntry );
+
+      for( auto nb = uint_t(0); nb < block->getNeighborhoodSize(); ++nb )
+      {
+         bs.sendBuffer( block->getNeighborProcess(nb) ) << infoCollectionEntry;
+      }
+
+      //note: is it necessary to add child blocks already into the info collection?
+      // here, we still have full geometrical information and can probably determine number of fluid and near boundary cells more easily
+      // however, the interesting (and most costly) blocks are never refined and thus their child infos is never needed
+      // see pe/amr/InfoCollection.cpp for an example
+
+   }
+
+   // size of buffer is unknown and changes with each send
+   bs.setReceiverInfoFromSendBufferState(false, true);
+   bs.sendAll();
+
+   for( auto recvIt = bs.begin(); recvIt != bs.end(); ++recvIt )
+   {
+      while( !recvIt.buffer().isEmpty() )
+      {
+         InfoCollectionPair val;
+         recvIt.buffer() >> val;
+         ic.insert(val);
+      }
+   }
+}
+
+void getBlockInfoFromInfoCollection( const PhantomBlock * block, const shared_ptr<InfoCollection>& ic,
+                                     BlockInfo & blockInfo )
+{
+   WALBERLA_ASSERT_NOT_NULLPTR(block);
+
+   if (block->sourceBlockIsLarger())
+   {
+      // block is a result of refinement -> BlockInfo object only available for the father block
+      // there should be no particles on the block (otherwise it would not have been refined)
+      // and refinement in LBM does not change the number of cells
+      // we assume that the number of fluid and near boundary cells also stays the same
+      // (ATTENTION: not true for blocks intersecting with a boundary!)
+      // -> we can use the information of the father block for weight assignment
+
+      auto infoIt = ic->find( block->getId().getFatherId() );
+      WALBERLA_CHECK_UNEQUAL( infoIt, ic->end(), "Father block with ID " << block->getId().getFatherId() << " not found in info collection!" );
+
+      // check the above mentioned assumptions
+      WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfLocalParticles, uint_t(0));
+      WALBERLA_ASSERT_EQUAL(infoIt->second.numberOfContacts, uint_t(0));
+
+      blockInfo = infoIt->second;
+   }
+   else if (block->sourceBlockHasTheSameSize())
+   {
+      auto infoIt = ic->find( block->getId() );
+      WALBERLA_CHECK_UNEQUAL( infoIt, ic->end(), "Block with ID " << block->getId() << " not found in info collection!" );
+      blockInfo = infoIt->second;
+   }
+   else
+   {
+      // source block of block is smaller
+
+      // block is a result of coarsening -> BlockInfo object is available on all 8 child blocks
+      // there should be no particles on the block (otherwise it would not have been coarsened)
+      // and refinement in LBM does not change the number of cells
+      // we assume that the number of fluid and near boundary cells will be the average of all 8 child blocks
+      // -> we can use the information of the child blocks for weight assignment
+
+      blockforest::BlockID childIdForInit(block->getId(), 0);
+      auto childForInitIt = ic->find( childIdForInit );
+      WALBERLA_CHECK_UNEQUAL( childForInitIt, ic->end(), "Child block with ID " << childIdForInit << " not found in info collection!" );
+      BlockInfo combinedInfo = childForInitIt->second;
+      uint_t numFluidCells(0);
+      uint_t numNearBoundaryCells(0);
+      for (uint_t child = 0; child < 8; ++child)
+      {
+         blockforest::BlockID childId(block->getId(), child);
+         auto childIt = ic->find( childId );
+         WALBERLA_CHECK_UNEQUAL( childIt, ic->end(), "Child block with ID " << childId << " not found in info collection!" );
+         numFluidCells += childIt->second.numberOfFluidCells;
+         numNearBoundaryCells += childIt->second.numberOfNearBoundaryCells;
+
+         // check above mentioned assumptions
+         WALBERLA_ASSERT_EQUAL(childIt->second.numberOfLocalParticles, uint_t(0));
+         WALBERLA_ASSERT_EQUAL(childIt->second.numberOfContacts, uint_t(0));
+      }
+      // total number of cells remains unchanged
+      combinedInfo.numberOfFluidCells = uint_c(numFluidCells / uint_t(8)); //average
+      combinedInfo.numberOfNearBoundaryCells = uint_c( numNearBoundaryCells / uint_t(8) ); //average
+      combinedInfo.numberOfLocalParticles = uint_t(0);
+      combinedInfo.numberOfContacts = uint_t(0); //sum
+      // number of pe sub cycles stays the same
+
+      blockInfo = combinedInfo;
+   }
+}
+
+} // namespace lbm_mesapd_coupling
+} // namespace walberla

src/lbm_mesapd_coupling/amr/level_determination/ParticlePresenceLevelDetermination.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 ParticlePresenceLevelDetermination.h
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "blockforest/BlockForest.h"
+#include "lbm_mesapd_coupling/amr/InfoCollection.h"
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+/*
+ * Class to determine the minimum level a block can be.
+ * For coupled LBM-PE simulations the following rules apply:
+ *  - a moving particle will always remain on the finest block
+ *  - a moving particle is not allowed to extend into an area with a coarser block
+ *  - if no moving particle is present, the level can be as coarse as possible (restricted by the 2:1 rule)
+ * Therefore, if a particle, local or remote (due to particles that are larger than a block), is present on any of the
+ * neighboring blocks of a certain block, this block's target level is the finest level.
+ * This, together with a refinement checking frequency that depends on the maximum translational particle velocity,
+ * ensures the above given requirements.
+ */
+class ParticlePresenceLevelDetermination
+{
+public:
+
+   ParticlePresenceLevelDetermination( const shared_ptr<lbm_mesapd_coupling::InfoCollection> & infoCollection, uint_t finestLevel) :
+         infoCollection_( infoCollection ), finestLevel_( finestLevel)
+   {}
+
+   void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                    std::vector< const Block * > &, const BlockForest & /*forest*/ ) {
+      for (auto &minTargetLevel : minTargetLevels) {
+         uint_t currentLevelOfBlock = minTargetLevel.first->getLevel();
+
+         const uint_t numberOfParticlesInDirectNeighborhood = getNumberOfLocalAndShadowParticlesInNeighborhood(minTargetLevel.first);
+
+         uint_t targetLevelOfBlock = currentLevelOfBlock; //keep everything as it is
+         if ( numberOfParticlesInDirectNeighborhood > uint_t(0) )
+         {
+            // set block to finest level if there are particles nearby
+            targetLevelOfBlock = finestLevel_;
+         }
+         else
+         {
+            // block could coarsen since there are no particles nearby
+            if( currentLevelOfBlock > uint_t(0) )
+               targetLevelOfBlock = currentLevelOfBlock - uint_t(1);
+         }
+
+         WALBERLA_CHECK_LESS_EQUAL(std::abs(int_c(targetLevelOfBlock) - int_c(currentLevelOfBlock)), uint_t(1), "Only level difference of maximum 1 allowed!");
+         minTargetLevel.second = targetLevelOfBlock;
+      }
+   }
+
+private:
+
+   uint_t getNumberOfLocalAndShadowParticlesInNeighborhood(const Block * block) {
+      auto numParticles = uint_t(0);
+
+      // add particles of current block
+      const auto infoIt = infoCollection_->find(block->getId());
+      WALBERLA_CHECK_UNEQUAL(infoIt, infoCollection_->end(), "Block with ID " << block->getId() << " not found in info collection!");
+
+      numParticles += infoIt->second.numberOfLocalParticles + infoIt->second.numberOfShadowParticles;
+
+      // add particles of all neighboring blocks
+      for(uint_t i = 0; i < block->getNeighborhoodSize(); ++i)
+      {
+         const BlockID &neighborBlockID = block->getNeighborId(i);
+         const auto infoItNeighbor = infoCollection_->find(neighborBlockID);
+         WALBERLA_CHECK_UNEQUAL(infoItNeighbor, infoCollection_->end(), "Neighbor block with ID " << neighborBlockID << " not found in info collection!");
+
+         numParticles += infoItNeighbor->second.numberOfLocalParticles + infoItNeighbor->second.numberOfShadowParticles;
+      }
+      return numParticles;
+   }
+
+   shared_ptr<lbm_mesapd_coupling::InfoCollection> infoCollection_;
+   uint_t finestLevel_;
+};
+
+} // namespace lbm_mesapd_coupling
+} // namespace walberla

src/lbm_mesapd_coupling/utility/SubCyclingManager.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 SubCyclingManager.cpp
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+
+#include "SubCyclingManager.h"
+
+#include <utility>
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+SubCyclingManager::SubCyclingManager(size_t numberOfSubCycles, shared_ptr<WcTimingPool> timingPool)
+   : numberOfSubCycles_(numberOfSubCycles), timingPool_(std::move(timingPool)), currentTimeStep_(0) {}
+
+void SubCyclingManager::operator()() {
+   executeBeforeFunctions();
+   for (size_t s = 0; s < numberOfSubCycles_; ++s) {
+      executeDuringFunctions();
+   }
+   executeAfterFunctions();
+   ++currentTimeStep_;
+}
+
+SubCyclingManager::FuncHandle
+SubCyclingManager::addFuncBeforeSubCycles(const VoidVoidFunc &f, const std::string &identifier) {
+   IdentifiedFunc idFunc(identifier, f);
+   beforeFunctions_.emplace_back(idFunc);
+   return beforeFunctions_.size() - 1;
+}
+
+SubCyclingManager::FuncHandle
+SubCyclingManager::addFuncDuringSubCycles(const VoidVoidFunc &f, const std::string &identifier) {
+   IdentifiedFunc idFunc(identifier, f);
+   duringFunctions_.emplace_back(idFunc);
+   return duringFunctions_.size() - 1;
+}
+
+SubCyclingManager::FuncHandle
+SubCyclingManager::addFuncAfterSubCycles(const VoidVoidFunc &f, const std::string &identifier) {
+   IdentifiedFunc idFunc(identifier, f);
+   afterFunctions_.emplace_back(idFunc);
+   return afterFunctions_.size() - 1;
+}
+
+inline void SubCyclingManager::executeBeforeFunctions() {
+   executeFunctions(beforeFunctions_);
+}
+
+inline void SubCyclingManager::executeDuringFunctions() {
+   executeFunctions(duringFunctions_);
+}
+
+inline void SubCyclingManager::executeAfterFunctions() {
+   executeFunctions(afterFunctions_);
+}
+
+inline void SubCyclingManager::executeFunctions(std::vector<IdentifiedFunc>& functions) {
+   for (const auto &f: functions) {
+      startTiming(f.first);
+      f.second();
+      stopTiming(f.first);
+   }
+}
+
+inline void SubCyclingManager::startTiming(const std::string & name){
+   if (timingPool_ != nullptr) {
+      (*timingPool_)[name].start();
+   }
+}
+
+inline void SubCyclingManager::stopTiming(const std::string & name){
+   if (timingPool_ != nullptr) {
+      (*timingPool_)[name].end();
+   }
+}
+
+}
+}
\ No newline at end of file

src/lbm_mesapd_coupling/utility/SubCyclingManager.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 SubCyclingManager.h
+//! \author Lukas Werner <lks.werner@fau.de>
+//
+//======================================================================================================================
+
+#include "core/DataTypes.h"
+#include "core/logging/Logging.h"
+#include "core/timing/TimingPool.h"
+
+namespace walberla {
+namespace lbm_mesapd_coupling {
+
+/*! \brief Handles the execution of subcycles in a timeloop to allow for finer grained time steps than the LBM ones.
+*
+* Supports registration of functions that are run before, during or after a subcycle.
+* The SubCyclingManager itself has to be added to a parent (LBM) time loop, and will execute the functions registered
+* for execution during the subcycling procedure numberOfSubCycles times.
+*/
+
+class SubCyclingManager {
+public:
+   typedef std::function<void ()> VoidVoidFunc;
+
+   /*! \name Construction & Destruction */
+   //@{
+   explicit SubCyclingManager(size_t numberOfSubCycles, shared_ptr<WcTimingPool> timingPool = nullptr);
+
+   virtual ~SubCyclingManager() {}
+   //@}
+
+   /*! \name Registration Functions */
+   //@{
+   typedef size_t FuncHandle;
+
+   FuncHandle addFuncBeforeSubCycles(const VoidVoidFunc &f, const std::string &identifier = "Other");
+   FuncHandle addFuncDuringSubCycles(const VoidVoidFunc &f, const std::string &identifier = "Other");
+   FuncHandle addFuncAfterSubCycles(const VoidVoidFunc &f, const std::string &identifier = "Other");
+   //@}
+
+   /*! \name Execution Functions */
+   //@{
+   void operator()();
+   //@}
+
+   /*! \name Bookkeeping Functions */
+   //@{
+   uint_t getCurrentTimeStep() const        { return currentTimeStep_;   }
+   void setCurrentTimeStep(uint_t timestep) { currentTimeStep_ = timestep; }
+   //@}
+
+private:
+   typedef std::pair<std::string, VoidVoidFunc> IdentifiedFunc;
+
+   inline void executeBeforeFunctions();
+   inline void executeDuringFunctions();
+   inline void executeAfterFunctions();
+
+   inline void executeFunctions(std::vector<IdentifiedFunc>& functions);
+
+   inline void startTiming(const std::string & name);
+   inline void stopTiming(const std::string & name);
+
+   size_t numberOfSubCycles_;
+   shared_ptr<WcTimingPool> timingPool_;
+   uint_t currentTimeStep_;
+
+   std::vector<IdentifiedFunc> beforeFunctions_;
+   std::vector<IdentifiedFunc> duringFunctions_;
+   std::vector<IdentifiedFunc> afterFunctions_;
+};
+
+}
+}
\ No newline at end of file