Merge branch 'DynamicRefinement' into 'master'

Dynamic refinement for pe module

See merge request walberla/walberla!52

src/pe/amr/regrid/RegridMinMax.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 RegridMinMax.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "RegridMinMax.h"
+
+namespace walberla {
+namespace pe {
+namespace amr {
+
+void ReGridMinMax::operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                               std::vector< const Block * > &, const BlockForest & /*forest*/ )
+{
+   for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it )
+   {
+      const auto infoIt = ic_->find(it->first->getId());
+      WALBERLA_ASSERT_UNEQUAL( infoIt, ic_->end() );
+
+      it->second = it->first->getLevel(); //keep everything as it is
+
+      //check for refinement
+      if (infoIt->second.numberOfLocalBodies > maxBodies_)
+      {
+         it->second = it->first->getLevel() + uint_t(1);
+         continue;
+      }
+
+      //check for coarsening
+      if ((it->first->getLevel() > 0) && (infoIt->second.numberOfLocalBodies < minBodies_))
+      {
+         if (getOrCreateCoarseInfo(it->first->getId())->second.numberOfLocalBodies < maxBodies_)
+         {
+            it->second = it->first->getLevel() - uint_t(1);
+         }
+         continue;
+      }
+   }
+}
+
+InfoCollection::const_iterator ReGridMinMax::getOrCreateCoarseInfo( const blockforest::BlockID& id )
+{
+   auto fatherId = id.getFatherId();
+   auto infoIt   = ic_->find( fatherId );
+   if (infoIt != ic_->end()) return infoIt;
+
+   BlockInfo newWeight( 0, 0);
+   for (uint_t child = 0; child < 8; ++child)
+   {
+      blockforest::BlockID childId(fatherId, child);
+      auto childIt = ic_->find( childId );
+      //meight be not available if not all blocks are on the same level
+      //return giant number to prevent coarsening
+      if (childIt == ic_->end())
+      {
+         newWeight = BlockInfo( std::numeric_limits<uint_t>::max(),
+                                std::numeric_limits<uint_t>::max());
+         break;
+      } else
+      {
+         newWeight += childIt->second;
+      }
+   }
+   WALBERLA_LOG_DETAIL("creating coarse weights (" << newWeight << ")");
+   return ic_->insert( std::make_pair(fatherId, newWeight) ).first;
+}
+
+} // namespace amr
+} // namespace pe
+} // namespace walberla

src/pe/amr/regrid/RegridMinMax.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 RegridMinMax.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include <pe/Types.h>
+#include <pe/amr/InfoCollection.h>
+
+#include <blockforest/Block.h>
+#include <blockforest/BlockForest.h>
+#include <core/logging/Logging.h>
+#include <domain_decomposition/BlockDataID.h>
+
+namespace walberla {
+namespace pe {
+namespace amr {
+
+class ReGridMinMax
+{
+public:
+
+   ReGridMinMax( const shared_ptr<InfoCollection>& ic, const size_t minBodies, const size_t maxBodies) :
+      ic_( ic ), minBodies_(minBodies), maxBodies_(maxBodies)
+   {}
+
+   void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                    std::vector< const Block * > &, const BlockForest & forest );
+
+public:
+   const shared_ptr<InfoCollection> ic_;
+   size_t      minBodies_;
+   size_t      maxBodies_;
+
+   InfoCollection::const_iterator getOrCreateCoarseInfo( const blockforest::BlockID& id );
+};
+
+} // namespace amr
+} // namespace pe
+} // namespace walberla

src/pe/amr/weight_assignment/MetisAssignmentFunctor.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 MetisAssignmentFunctor.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "pe/amr/InfoCollection.h"
+
+#include "blockforest/loadbalancing/DynamicParMetis.h"
+
+namespace walberla {
+namespace pe {
+namespace amr {
+
+class MetisAssignmentFunctor
+{
+public:
+
+   typedef blockforest::DynamicParMetisBlockInfo           PhantomBlockWeight;
+   typedef blockforest::DynamicParMetisBlockInfoPackUnpack PhantomBlockWeightPackUnpackFunctor;
+
+   MetisAssignmentFunctor( const shared_ptr<InfoCollection>& ic ) : ic_( ic )
+   {}
+
+   void operator()( std::vector< std::pair< const PhantomBlock *, boost::any > > & blockData, const PhantomBlockForest & )
+   {
+      for( auto it = blockData.begin(); it != blockData.end(); ++it )
+      {
+         const uint_t& weight = ic_->find( it->first->getId() )->second.numberOfLocalBodies;
+         blockforest::DynamicParMetisBlockInfo info( int64_c(weight) );
+         info.setVertexSize(int64_c( weight ));
+         for( uint_t nb = uint_t(0); nb < it->first->getNeighborhoodSize(); ++nb )
+         {
+            info.setEdgeWeight(it->first->getNeighborId(nb), int64_c(weight) );
+         }
+         it->second = info;
+      }
+   }
+
+private:
+   shared_ptr< InfoCollection > ic_;
+};
+
+}
+}
+}
+

src/pe/amr/weight_assignment/WeightAssignmentFunctor.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 WeightAssignmentFunctor.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "WeightAssignmentFunctor.h"
+
+namespace walberla {
+namespace pe {
+namespace amr {
+
+const double WeightAssignmentFunctor::baseWeight = real_t(10.0);
+
+}
+}
+}

src/pe/amr/weight_assignment/WeightAssignmentFunctor.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 WeightAssignmentFunctor.h
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#pragma once
+
+#include "pe/amr/InfoCollection.h"
+
+#include "blockforest/loadbalancing/PODPhantomData.h"
+
+namespace walberla {
+namespace pe {
+namespace amr {
+
+class WeightAssignmentFunctor
+{
+public:
+   typedef walberla::blockforest::PODPhantomWeight<double>           PhantomBlockWeight;
+   typedef walberla::blockforest::PODPhantomWeightPackUnpack<double> PhantomBlockWeightPackUnpackFunctor;
+
+   ///Base weight due to allocated data structures. A weight of zero for blocks is dangerous as empty blocks might accumulate on one process!
+   static const double baseWeight;
+
+   WeightAssignmentFunctor( shared_ptr<InfoCollection>& ic ) : ic_(ic) {}
+
+   void operator()( std::vector< std::pair< const PhantomBlock *, boost::any > > & blockData, const PhantomBlockForest & )
+   {
+      for( auto it = blockData.begin(); it != blockData.end(); ++it )
+      {
+         const PhantomBlock * block = it->first;
+         //only change of one level is supported!
+         WALBERLA_ASSERT_LESS( int_c(block->getLevel()) - int_c(block->getSourceLevel()), 2 );
+
+         if (block->sourceBlockIsLarger())
+         {
+            auto infoIt = ic_->find( block->getId().getFatherId() );
+            WALBERLA_ASSERT_UNEQUAL( infoIt, ic_->end() );
+            it->second = PhantomBlockWeight( double_c(infoIt->second.numberOfLocalBodies) / double_c(8) + baseWeight );
+            continue;
+         }
+
+         if (block->sourceBlockHasTheSameSize())
+         {
+            auto infoIt = ic_->find( block->getId() );
+            WALBERLA_ASSERT_UNEQUAL( infoIt, ic_->end() );
+            it->second = PhantomBlockWeight( double_c(infoIt->second.numberOfLocalBodies) + baseWeight );
+            continue;
+         }
+
+         if (block->sourceBlockIsSmaller())
+         {
+            double weight = 0;
+            for (uint_t child = 0; child < 8; ++child)
+            {
+               blockforest::BlockID childId(block->getId(), child);
+               auto childIt = ic_->find( childId );
+               WALBERLA_ASSERT_UNEQUAL( childIt, ic_->end() );
+               weight += double_c(childIt->second.numberOfLocalBodies);
+            }
+            it->second = PhantomBlockWeight( weight + baseWeight );
+            continue;
+         }
+      }
+   }
+
+private:
+   shared_ptr<InfoCollection> ic_;
+};
+
+}
+}
+}

tests/pe/CMakeLists.txt

@@ -35,6 +35,9 @@ waLBerla_execute_test( NAME   PE_DESTROYBODY )
 waLBerla_compile_test( NAME   PE_DOCUMENTATIONSNIPPETS FILES PeDocumentationSnippets.cpp DEPENDS core  )
 waLBerla_execute_test( NAME   PE_DOCUMENTATIONSNIPPETS )
 
+waLBerla_compile_test( NAME   PE_DYNAMICREFINEMENT FILES DynamicRefinement.cpp DEPENDS core blockforest  )
+waLBerla_execute_test( NAME   PE_DYNAMICREFINEMENT )
+
 waLBerla_compile_test( NAME   PE_FORCESYNC FILES ForceSync.cpp DEPENDS core blockforest  )
 waLBerla_execute_test( NAME   PE_FORCESYNC )
 
@@ -54,6 +57,9 @@ waLBerla_execute_test( NAME   PE_MARSHALLING )
 waLBerla_compile_test( NAME   PE_MATERIAL FILES Material.cpp DEPENDS core  )
 waLBerla_execute_test( NAME   PE_MATERIAL )
 
+waLBerla_compile_test( NAME   PE_MINMAXREFINEMENT FILES MinMaxRefinement.cpp DEPENDS core blockforest  )
+waLBerla_execute_test( NAME   PE_MINMAXREFINEMENT PROCESSES 8 )
+
 waLBerla_compile_test( NAME   PE_OVERLAP FILES Overlap.cpp DEPENDS core  )
 waLBerla_execute_test( NAME   PE_OVERLAP )
 

tests/pe/DynamicRefinement.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 DynamicRefinement.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include "pe/basic.h"
+#include "pe/synchronization/ClearSynchronization.h"
+#include "pe/utility/GetBody.h"
+#include "pe/utility/DestroyBody.h"
+
+#include "blockforest/Initialization.h"
+#include "core/all.h"
+#include "domain_decomposition/all.h"
+
+#include "core/debug/TestSubsystem.h"
+
+using namespace walberla;
+using namespace walberla::pe;
+
+typedef boost::tuple<Sphere> BodyTuple ;
+
+class ReGrid
+{
+public:
+
+   ReGrid( const BlockDataID storageID, const size_t minParticles, const size_t maxParticles) :
+      storageID_( storageID ), minParticles_(minParticles), maxParticles_(maxParticles)
+   {}
+
+   void operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                    std::vector< const Block * > &, const BlockForest & forest );
+
+private:
+   const BlockDataID storageID_;
+   const size_t      minParticles_;
+   const size_t      maxParticles_;
+};
+
+void ReGrid::operator()( std::vector< std::pair< const Block *, uint_t > > & minTargetLevels,
+                         std::vector< const Block * > &, const BlockForest & /*forest*/ )
+{
+   for( auto it = minTargetLevels.begin(); it != minTargetLevels.end(); ++it )
+   {
+      const auto numberOfParticles = (*(it->first->getData< Storage >( storageID_ )))[0].size();
+      //WALBERLA_LOG_DEVEL("storage size: " << localBodyStorage.size());
+
+      it->second = it->first->getLevel(); //keep everything as it is
+      if (numberOfParticles < minParticles_)
+      {
+         WALBERLA_LOG_DEVEL(it->first->getLevel() << " -> " << it->first->getLevel() - uint_t(1) << " (" << numberOfParticles << ")" );
+         if (it->first->getLevel() > 0)
+            it->second = it->first->getLevel() - uint_t(1);
+      } else if (numberOfParticles > maxParticles_)
+      {
+         it->second = it->first->getLevel() + uint_t(1);
+         WALBERLA_LOG_DEVEL(it->first->getLevel() << " -> " << it->first->getLevel() + uint_t(1) << " (" << numberOfParticles << ")" );
+      }
+   }
+}
+
+int main( int argc, char** argv )
+{
+   walberla::debug::enterTestMode();
+   walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+
+   shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+
+   // create blocks
+   shared_ptr< StructuredBlockForest > forest = blockforest::createUniformBlockGrid(
+            math::AABB(0,0,0,20,20,20),
+            uint_c( 1), uint_c( 1), uint_c( 1), // number of blocks in x,y,z direction
+            uint_c( 1), uint_c( 1), uint_c( 1), // how many cells per block (x,y,z)
+            false,                              // max blocks per process
+            false, false, false,                // full periodicity
+            false);
+
+   SetBodyTypeIDs<BodyTuple>::execute();
+
+   auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+   forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "HCCD");
+   forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
+
+   auto & blockforest = forest->getBlockForest();
+   blockforest.recalculateBlockLevelsInRefresh( true );
+   blockforest.alwaysRebalanceInRefresh( false );
+   blockforest.reevaluateMinTargetLevelsAfterForcedRefinement( false );
+   blockforest.allowRefreshChangingDepth( true );
+
+   blockforest.allowMultipleRefreshCycles( false );
+   blockforest.checkForEarlyOutInRefresh( true );
+   blockforest.checkForLateOutInRefresh( true );
+
+   ReGrid regrid( storageID, 20, 20 );
+
+   blockforest.setRefreshMinTargetLevelDeterminationFunction( regrid );
+
+   blockforest.setRefreshPhantomBlockMigrationPreparationFunction(
+            blockforest::DynamicLevelwiseCurveBalance< blockforest::NoPhantomData >( true, true ) );
+
+   real_t spacing(2.5);
+   for (auto blkIt = forest->begin(); blkIt != forest->end(); ++blkIt)
+   {
+      IBlock & currentBlock = *blkIt;
+      for (auto it = grid_generator::SCIterator(currentBlock.getAABB(), Vector3<real_t>(spacing) * real_t(0.5), spacing); it != grid_generator::SCIterator(); ++it)
+      {
+         createSphere( *globalBodyStorage, forest->getBlockStorage(), storageID, 0, *it, 1 );
+      }
+   }
+   syncNextNeighbors<BodyTuple>(forest->getBlockForest(), storageID);
+   syncNextNeighbors<BodyTuple>(forest->getBlockForest(), storageID);
+
+   clearSynchronization( forest->getBlockForest(), storageID );
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(forest->getBlockForest(), storageID);
+
+   WALBERLA_ASSERT_EQUAL( forest->size(), 8 );
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+//      IBlock & currentBlock = *blockIt;
+//      Storage * storage = currentBlock.getData< Storage >( storageID );
+//      BodyStorage& localStorage = (*storage)[0];
+//      BodyStorage& shadowStorage = (*storage)[1];
+
+      for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID); bodyIt != LocalBodyIterator::end(); ++bodyIt)
+      {
+         WALBERLA_ASSERT( blockIt->getAABB().contains(bodyIt->getPosition()) );
+
+//         WALBERLA_LOG_DEVEL( blockIt->getAABB() );
+//         WALBERLA_LOG_DEVEL(*bodyIt );
+      }
+   }
+
+   WALBERLA_LOG_DEVEL("========================================================");
+
+   clearSynchronization( forest->getBlockForest(), storageID );
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(forest->getBlockForest(), storageID);
+
+   WALBERLA_ASSERT_EQUAL( forest->size(), 64 );
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+//      IBlock & currentBlock = *blockIt;
+//      Storage * storage = currentBlock.getData< Storage >( storageID );
+//      BodyStorage& localStorage = (*storage)[0];
+//      BodyStorage& shadowStorage = (*storage)[1];
+
+      for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID); bodyIt != LocalBodyIterator::end(); ++bodyIt)
+      {
+         WALBERLA_ASSERT( blockIt->getAABB().contains(bodyIt->getPosition()) );
+
+//         WALBERLA_LOG_DEVEL( blockIt->getAABB() );
+//         WALBERLA_LOG_DEVEL(*bodyIt );
+      }
+   }
+
+   WALBERLA_LOG_DEVEL("========================================================");
+
+   clearSynchronization( forest->getBlockForest(), storageID );
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(forest->getBlockForest(), storageID);
+
+   WALBERLA_ASSERT_EQUAL( forest->size(), 8 );
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+//      IBlock & currentBlock = *blockIt;
+//      Storage * storage = currentBlock.getData< Storage >( storageID );
+//      BodyStorage& localStorage = (*storage)[0];
+//      BodyStorage& shadowStorage = (*storage)[1];
+
+      for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID); bodyIt != LocalBodyIterator::end(); ++bodyIt)
+      {
+         WALBERLA_ASSERT( blockIt->getAABB().contains(bodyIt->getPosition()) );
+
+//         WALBERLA_LOG_DEVEL( blockIt->getAABB() );
+//         WALBERLA_LOG_DEVEL(*bodyIt );
+      }
+   }
+
+   return EXIT_SUCCESS;
+}

tests/pe/MinMaxRefinement.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 Refinement.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+
+#include "blockforest/all.h"
+#include <blockforest/loadbalancing/PODPhantomData.h>
+#include "core/all.h"
+#include "domain_decomposition/all.h"
+#include "timeloop/SweepTimeloop.h"
+#include "vtk/VTKOutput.h"
+
+
+#include "pe/basic.h"
+#include "pe/amr/InfoCollection.h"
+#include "pe/amr/regrid/RegridMinMax.h"
+#include "pe/amr/weight_assignment/WeightAssignmentFunctor.h"
+#include "pe/ccd/SimpleCCDDataHandling.h"
+#include "pe/synchronization/SyncNextNeighbors.h"
+#include "pe/synchronization/ClearSynchronization.h"
+#include "pe/vtk/BodyVtkOutput.h"
+#include "pe/vtk/SphereVtkOutput.h"
+
+#include "CheckVitalParameters.h"
+
+#include "core/debug/TestSubsystem.h"
+
+#include <boost/tuple/tuple.hpp>
+
+#include <algorithm>
+#include <limits>
+#include <vector>
+
+using namespace walberla;
+using namespace walberla::pe;
+
+typedef boost::tuple<Sphere, Plane> BodyTuple ;
+
+int main( int argc, char ** argv )
+{
+   using namespace walberla::pe;
+
+   debug::enterTestMode();
+
+   walberla::MPIManager::instance()->initializeMPI( &argc, &argv );
+
+   //      logging::Logging::instance()->setStreamLogLevel( logging::Logging::DETAIL );
+   //   logging::Logging::instance()->setFileLogLevel( logging::Logging::DETAIL );
+   //   logging::Logging::instance()->includeLoggingToFile("SyncLog");
+
+   shared_ptr<BodyStorage> globalStorage = make_shared<BodyStorage>();
+
+   // create forest
+   shared_ptr< blockforest::StructuredBlockForest > forest = blockforest::createUniformBlockGrid(
+            math::AABB(0,0,0,4,4,4),
+            1,1,1,                                                 // number of blocks in x,y,z direction
+            1,1,1,                                                  // how many cells per block (x,y,z)
+            0,                                                      // max blocks per process
+            false, false,                                           // include metis / force metis
+            false, false, false );                                    // full periodicity
+
+   SetBodyTypeIDs<BodyTuple>::execute();
+
+   auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+   auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalStorage, storageID ), "CCD");
+   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
+   WALBERLA_UNUSED(fcdID);
+
+   auto & blockforest = forest->getBlockForest();
+
+   //***** SETUP LOADBALACING & REFINEMENT
+   blockforest.recalculateBlockLevelsInRefresh( true );
+   blockforest.alwaysRebalanceInRefresh( true );
+   blockforest.reevaluateMinTargetLevelsAfterForcedRefinement( false );
+   blockforest.allowRefreshChangingDepth( true );
+
+   blockforest.allowMultipleRefreshCycles( false );
+   blockforest.checkForEarlyOutInRefresh( true );
+   blockforest.checkForLateOutInRefresh( true );
+
+   auto infoCollection = make_shared<InfoCollection>();
+
+   amr::ReGridMinMax regrid(infoCollection, 2, 5);
+   blockforest.setRefreshMinTargetLevelDeterminationFunction( regrid );
+
+   blockforest.setRefreshPhantomBlockDataAssignmentFunction( amr::WeightAssignmentFunctor( infoCollection ) );
+   blockforest.setRefreshPhantomBlockDataPackFunction( amr::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor() );
+   blockforest.setRefreshPhantomBlockDataUnpackFunction( amr::WeightAssignmentFunctor::PhantomBlockWeightPackUnpackFunctor() );
+
+   blockforest.setRefreshPhantomBlockMigrationPreparationFunction(
+            blockforest::DynamicLevelwiseCurveBalance< amr::WeightAssignmentFunctor::PhantomBlockWeight >( false, true, false ) );
+
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(1,1,1), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(1,1,3), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(1,3,1), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(1,3,3), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(3,1,1), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(3,1,3), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(3,3,1), 1);
+   createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(3,3,3), 1);
+
+   WALBERLA_MPI_BARRIER();
+   WALBERLA_LOG_DEVEL_ON_ROOT( "Refinement 1" );
+   createWithNeighborhood(blockforest, storageID, *infoCollection);
+   clearSynchronization( blockforest, storageID);
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(blockforest, storageID);
+
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+      ccd::ICCD* ccd = blockIt->getData< ccd::ICCD >( ccdID );
+      ccd->reloadBodies();
+   }
+
+   WALBERLA_CHECK_EQUAL( blockforest.size(), 1);
+
+   WALBERLA_MPI_BARRIER();
+   WALBERLA_LOG_DEVEL_ON_ROOT( "Refinement 2" );
+   blockforest.setRefreshMinTargetLevelDeterminationFunction( amr::ReGridMinMax(infoCollection, 9, 20) );
+   createWithNeighborhood(blockforest, storageID, *infoCollection);
+   clearSynchronization( blockforest, storageID);
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(blockforest, storageID);
+
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+      ccd::ICCD* ccd = blockIt->getData< ccd::ICCD >( ccdID );
+      ccd->reloadBodies();
+   }
+
+   WALBERLA_CHECK_EQUAL( blockforest.size(), mpi::MPIManager::instance()->worldRank() == 6 ? 1 : 0);
+   WALBERLA_LOG_DEVEL( infoCollection->size() );
+
+   for (unsigned int i = 0; i < 30; ++i)
+   {
+      createSphere(*globalStorage.get(), forest->getBlockStorage(), storageID, 0, Vec3(real_t(2.1), real_t(2.1), real_t(2.1)), 1);
+   }
+
+   WALBERLA_MPI_BARRIER();
+   WALBERLA_LOG_DEVEL_ON_ROOT( "Refinement 3" );
+   blockforest.setRefreshMinTargetLevelDeterminationFunction( amr::ReGridMinMax(infoCollection, 2, 3) );
+   createWithNeighborhood(blockforest, storageID, *infoCollection);
+   clearSynchronization( blockforest, storageID);
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(blockforest, storageID);
+
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+      ccd::ICCD* ccd = blockIt->getData< ccd::ICCD >( ccdID );
+      ccd->reloadBodies();
+   }
+
+   WALBERLA_LOG_DEVEL( infoCollection->size() );
+
+   WALBERLA_MPI_BARRIER();
+   WALBERLA_LOG_DEVEL_ON_ROOT( "Refinement 4" );
+   createWithNeighborhood(blockforest, storageID, *infoCollection);
+   clearSynchronization( blockforest, storageID);
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(blockforest, storageID);
+
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+      ccd::ICCD* ccd = blockIt->getData< ccd::ICCD >( ccdID );
+      ccd->reloadBodies();
+   }
+
+   WALBERLA_LOG_DEVEL( infoCollection->size() );
+
+   WALBERLA_MPI_BARRIER();
+   WALBERLA_LOG_DEVEL_ON_ROOT( "Refinement 5" );
+   WALBERLA_LOG_DEVEL( "SIZE: " << blockforest.size() );
+   createWithNeighborhood(blockforest, storageID, *infoCollection);
+   clearSynchronization( blockforest, storageID);
+   forest->refresh();
+   syncNextNeighbors<BodyTuple>(blockforest, storageID);
+
+   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
+   {
+      ccd::ICCD* ccd = blockIt->getData< ccd::ICCD >( ccdID );
+      ccd->reloadBodies();
+   }
+
+   WALBERLA_LOG_DEVEL( infoCollection->size() );
+
+   return EXIT_SUCCESS;
+}