Added tests for the force and torque storage functionalities

tests/pe_coupling/CMakeLists.txt

@@ -156,3 +156,15 @@ waLBerla_execute_test( NAME SphereWallCollisionBehaviorDPMFuncTest COMMAND $<TAR
 
 waLBerla_compile_test( FILES discrete_particle_methods/HinderedSettlingDynamicsDPM.cpp DEPENDS blockforest pe timeloop )
 waLBerla_execute_test( NAME HinderedSettlingDynamicsDPMFuncTest COMMAND $<TARGET_FILE:HinderedSettlingDynamicsDPM> --funcTest PROCESSES 4 LABELS longrun CONFIGURATIONS RelWithDbgInfo )
+
+###################################################################################################
+# Utility tests
+###################################################################################################
+
+waLBerla_compile_test( FILES utility/BodiesForceTorqueContainerTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME BodiesForceTorqueContainerTest COMMAND $<TARGET_FILE:BodiesForceTorqueContainerTest> PROCESSES 1 )
+waLBerla_execute_test( NAME BodiesForceTorqueContainerParallelTest COMMAND $<TARGET_FILE:BodiesForceTorqueContainerTest> PROCESSES 3 )
+
+waLBerla_compile_test( FILES utility/PeSubCyclingTest.cpp DEPENDS blockforest pe timeloop )
+waLBerla_execute_test( NAME PeSubCyclingTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 1 )
+waLBerla_execute_test( NAME PeSubCyclingParallelTest COMMAND $<TARGET_FILE:PeSubCyclingTest> PROCESSES 3 )
\ No newline at end of file

tests/pe_coupling/utility/BodiesForceTorqueContainerTest.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 BodiesForceTorqueContainerTest.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+
+#include "pe/basic.h"
+#include "pe/utility/DestroyBody.h"
+
+#include <pe_coupling/utility/all.h>
+
+namespace force_torque_container_test
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+
+typedef boost::tuple<pe::Sphere> BodyTypeTuple ;
+
+
+/*!\brief Test cases for the force torque container provided by the coupling module
+ *
+ * Spheres at different positions are created and force(s) and torque(s) are applied.
+ * Then, they are store in the container and reset on the body.
+ * Then, in some cases, the sphere is moved to cross block borders, and the stored forces/torques are reapplied by the container again.
+ * The obtained force/torque values are compared against the originally set (and thus expected) values.
+ */
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   // uncomment to have logging
+   //logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::DETAIL);
+
+   const real_t dx     = real_t(1);
+   const real_t radius = real_t(5);
+
+   ///////////////////////////
+   // DATA STRUCTURES SETUP //
+   ///////////////////////////
+
+   Vector3<uint_t> blocksPerDirection(uint_t(3), uint_t(1), uint_t(1));
+   Vector3<uint_t> cellsPerBlock(uint_t(20), uint_t(20), uint_t(20));
+   Vector3<bool> periodicity(true, false, false);
+
+   // create fully periodic domain with refined blocks
+   auto blocks = blockforest::createUniformBlockGrid( blocksPerDirection[0], blocksPerDirection[1], blocksPerDirection[2],
+                                                      cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2],
+                                                      dx,
+                                                      0, false, false,
+                                                      periodicity[0], periodicity[1], periodicity[2],
+                                                      false );
+
+
+   // pe body storage
+   pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
+   auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
+   auto sphereMaterialID = pe::createMaterial( "sphereMat", real_t(1) , real_t(0.3), real_t(0.2), real_t(0.2), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
+
+   // pe coupling
+   const real_t overlap = real_t( 1.5 ) * dx;
+   std::function<void(void)> syncCall = std::bind( pe::syncNextNeighbors<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+
+   // sphere positions for test scenarios
+   Vector3<real_t> positionInsideBlock(real_t(10), real_t(10), real_t(10));
+   Vector3<real_t> positionAtBlockBorder(real_t(19.5), real_t(10), real_t(10));
+   Vector3<real_t> positionAtBlockBorderUpdated(real_t(20.5), real_t(10), real_t(10));
+
+   Vector3<real_t> positionAtBlockBorder2(real_t(20) + radius + overlap - real_t(0.5), real_t(10), real_t(10));
+   Vector3<real_t> positionAtBlockBorderUpdated2(real_t(20) + radius + overlap + real_t(0.5), real_t(10), real_t(10));
+
+
+   Vector3<real_t> testForce(real_t(2), real_t(1), real_t(0));
+   Vector3<real_t> torqueOffset = Vector3<real_t>(real_t(1), real_t(0), real_t(0));
+
+   pe_coupling::ForceTorqueOnBodiesResetter resetter(blocks, bodyStorageID);
+   shared_ptr<pe_coupling::BodiesForceTorqueContainer> container1 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
+   shared_ptr<pe_coupling::BodiesForceTorqueContainer> container2 = make_shared<pe_coupling::BodiesForceTorqueContainer>(blocks, bodyStorageID);
+   pe_coupling::BodyContainerSwapper containerSwapper(container1, container2);
+
+   //////////////////
+   // Inside block //
+   //////////////////
+   {
+      std::string testIdentifier("Test: sphere inside block");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       positionInsideBlock, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      uint_t applicationCounter( 0 );
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            ++applicationCounter;
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+      mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+      container1->store();
+
+      resetter();
+
+      containerSwapper();
+
+      container2->setOnBodies();
+
+      Vector3<real_t> expectedForce = applicationCounter * testForce;
+      Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+      Vector3<real_t> actingForce(real_t(0));
+      Vector3<real_t> actingTorque(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            actingForce += bodyIt->getForce();
+            actingTorque += bodyIt->getTorque();
+
+         }
+      }
+
+      mpi::allReduceInplace(actingForce[0], mpi::SUM);
+      mpi::allReduceInplace(actingForce[1], mpi::SUM);
+      mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+      mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+      WALBERLA_ROOT_SECTION()
+      {
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      container1->clear();
+      container2->clear();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   /////////////////////
+   // At block border //
+   /////////////////////
+   {
+      std::string testIdentifier("Test: sphere at block border");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       positionAtBlockBorder, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      uint_t applicationCounter( 0 );
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            ++applicationCounter;
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+      mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+      container1->store();
+
+      resetter();
+
+      containerSwapper();
+
+      container2->setOnBodies();
+
+      Vector3<real_t> expectedForce = applicationCounter * testForce;
+      Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+      Vector3<real_t> actingForce(real_t(0));
+      Vector3<real_t> actingTorque(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            actingForce += bodyIt->getForce();
+            actingTorque += bodyIt->getTorque();
+         }
+      }
+
+      mpi::allReduceInplace(actingForce[0], mpi::SUM);
+      mpi::allReduceInplace(actingForce[1], mpi::SUM);
+      mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+      mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+      WALBERLA_ROOT_SECTION()
+      {
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      container1->clear();
+      container2->clear();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   ///////////////////////////////////////////
+   // At block border with updated position //
+   ///////////////////////////////////////////
+   {
+      std::string testIdentifier("Test: sphere at block border with updated position");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       positionAtBlockBorder, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      uint_t applicationCounter( 0 );
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            ++applicationCounter;
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+      mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+      container1->store();
+
+      resetter();
+
+      containerSwapper();
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) {
+         for (auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt) {
+            bodyIt->setPosition(positionAtBlockBorderUpdated);
+         }
+      }
+      syncCall();
+
+      container2->setOnBodies();
+
+      Vector3<real_t> expectedForce = applicationCounter * testForce;
+      Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+      Vector3<real_t> actingForce(real_t(0));
+      Vector3<real_t> actingTorque(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            actingForce += bodyIt->getForce();
+            actingTorque += bodyIt->getTorque();
+         }
+      }
+
+      mpi::allReduceInplace(actingForce[0], mpi::SUM);
+      mpi::allReduceInplace(actingForce[1], mpi::SUM);
+      mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+      mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+      WALBERLA_ROOT_SECTION()
+      {
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      container1->clear();
+      container2->clear();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   /////////////////////////////////////////////
+   // At block border with updated position 2 //
+   /////////////////////////////////////////////
+   {
+      std::string testIdentifier("Test: sphere at block border with updated position 2");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       positionAtBlockBorder2, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      uint_t applicationCounter( 0 );
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            ++applicationCounter;
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+      mpi::allReduceInplace(applicationCounter, mpi::SUM);
+
+      container1->store();
+
+      resetter();
+
+      containerSwapper();
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt ) {
+         for (auto bodyIt = pe::LocalBodyIterator::begin(*blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt) {
+            bodyIt->setPosition(positionAtBlockBorderUpdated2);
+         }
+      }
+      syncCall();
+
+      container2->setOnBodies();
+
+      --applicationCounter; // sphere has vanished from one block
+
+      // in this case, the complete force can no longer be applied onto the body since the part from the block
+      // from which the body vanished is lost
+
+      // HOWEVER: this case will never appear in a coupled simulation since NO FORCE will be acting a body
+      // that is about to vanish from a block since it will not be mapped to the block from which it will vanish
+
+      // If you are expecting a different behavior, you need to change the container mechanism by first all-reducing the
+      // forces/torques to all processes/blocks that know this body such that every process/block has
+      // the full information and then the process/block that owns this body sets the complete force
+
+      Vector3<real_t> expectedForce = applicationCounter * testForce;
+      Vector3<real_t> expectedTorque = applicationCounter * ( torqueOffset % testForce );
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting force: " << expectedForce);
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting torque: " << expectedTorque);
+
+      Vector3<real_t> actingForce(real_t(0));
+      Vector3<real_t> actingTorque(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            actingForce += bodyIt->getForce();
+            actingTorque += bodyIt->getTorque();
+         }
+      }
+
+      mpi::allReduceInplace(actingForce[0], mpi::SUM);
+      mpi::allReduceInplace(actingForce[1], mpi::SUM);
+      mpi::allReduceInplace(actingForce[2], mpi::SUM);
+
+      mpi::allReduceInplace(actingTorque[0], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[1], mpi::SUM);
+      mpi::allReduceInplace(actingTorque[2], mpi::SUM);
+
+      WALBERLA_ROOT_SECTION()
+      {
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[0], expectedForce[0], "Mismatch in force0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[1], expectedForce[1], "Mismatch in force1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingForce[2], expectedForce[2], "Mismatch in force2");
+
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[0], expectedTorque[0], "Mismatch in torque0");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[1], expectedTorque[1], "Mismatch in torque1");
+         WALBERLA_CHECK_FLOAT_EQUAL(actingTorque[2], expectedTorque[2], "Mismatch in torque2");
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      container1->clear();
+      container2->clear();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   return 0;
+
+}
+
+} //namespace force_torque_container_test
+
+int main( int argc, char **argv ){
+   force_torque_container_test::main(argc, argv);
+}

tests/pe_coupling/utility/PeSubCyclingTest.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 PeSubCyclingTest.cpp
+//! \ingroup pe_coupling
+//! \author Christoph Rettinger <christoph.rettinger@fau.de>
+//
+//======================================================================================================================
+
+#include "blockforest/Initialization.h"
+
+#include "core/DataTypes.h"
+#include "core/Environment.h"
+#include "core/debug/TestSubsystem.h"
+#include "core/math/all.h"
+
+#include "pe/basic.h"
+#include "pe/utility/DestroyBody.h"
+#include "pe/cr/DEM.h"
+
+#include <pe_coupling/utility/all.h>
+
+namespace pe_sub_cycling_test
+{
+
+///////////
+// USING //
+///////////
+
+using namespace walberla;
+
+typedef boost::tuple<pe::Sphere> BodyTypeTuple ;
+
+/*!\brief test case to check functionality of sub cycling in the pe time step provided by the coupling module
+ *
+ * During this time step, currently acting forces/torques on a body are kept constant.
+ *
+ * This test first computes the resulting position offset as well as the linear and rotational velocity after a
+ * certain force has been applied to a sphere when using several 'real' pe time steps and re-applying the force 10 times.
+ *
+ * It then checks the pe time step sub-cycling functionality of the coupling module by creating same-sized spheres
+ * at different locations and applying the same force. But this time, 10 sub-cycles are carried out.
+ * As a result, the same position offset and velocities must be obtained as in the regular case.
+ *
+ */
+//////////
+// MAIN //
+//////////
+int main( int argc, char **argv )
+{
+   debug::enterTestMode();
+
+   mpi::Environment env( argc, argv );
+
+   // uncomment to have logging
+   //logging::Logging::instance()->setLogLevel(logging::Logging::LogLevel::DETAIL);
+
+   const real_t dx     = real_t(1);
+   const real_t radius = real_t(5);
+
+   ///////////////////////////
+   // DATA STRUCTURES SETUP //
+   ///////////////////////////
+
+   Vector3<uint_t> blocksPerDirection(uint_t(3), uint_t(1), uint_t(1));
+   Vector3<uint_t> cellsPerBlock(uint_t(20), uint_t(20), uint_t(20));
+   Vector3<bool> periodicity(true, false, false);
+
+   // create fully periodic domain with refined blocks
+   auto blocks = blockforest::createUniformBlockGrid( blocksPerDirection[0], blocksPerDirection[1], blocksPerDirection[2],
+                                                      cellsPerBlock[0], cellsPerBlock[1], cellsPerBlock[2],
+                                                      dx,
+                                                      0, false, false,
+                                                      periodicity[0], periodicity[1], periodicity[2],
+                                                      false );
+
+
+   // pe body storage
+   pe::SetBodyTypeIDs<BodyTypeTuple>::execute();
+   shared_ptr<pe::BodyStorage> globalBodyStorage = make_shared<pe::BodyStorage>();
+   auto bodyStorageID = blocks->addBlockData(pe::createStorageDataHandling<BodyTypeTuple>(), "Storage");
+   auto sphereMaterialID = pe::createMaterial( "sphereMat", real_t(1) , real_t(0.3), real_t(0.2), real_t(0.2), real_t(0.24), real_t(200), real_t(200), real_t(0), real_t(0) );
+
+   auto ccdID = blocks->addBlockData(pe::ccd::createHashGridsDataHandling( globalBodyStorage, bodyStorageID ), "CCD");
+   auto fcdID = blocks->addBlockData(pe::fcd::createGenericFCDDataHandling<BodyTypeTuple, pe::fcd::AnalyticCollideFunctor>(), "FCD");
+   pe::cr::DEM cr(globalBodyStorage, blocks->getBlockStoragePointer(), bodyStorageID, ccdID, fcdID, NULL);
+
+   // set up synchronization procedure
+   const real_t overlap = real_t( 1.5 ) * dx;
+   std::function<void(void)> syncCall = std::bind( pe::syncNextNeighbors<BodyTypeTuple>, boost::ref(blocks->getBlockForest()), bodyStorageID, static_cast<WcTimingTree*>(NULL), overlap, false );
+
+
+   // sphere positions for test scenarios
+   Vector3<real_t> positionInsideBlock(real_t(10), real_t(10), real_t(10));
+   Vector3<real_t> positionAtBlockBorder(real_t(19.9), real_t(10), real_t(10));
+   Vector3<real_t> positionAtBlockBorder2(real_t(20) + radius + overlap - real_t(0.1), real_t(10), real_t(10));
+
+   Vector3<real_t> testForce(real_t(2), real_t(1), real_t(0));
+   Vector3<real_t> torqueOffset = Vector3<real_t>(real_t(1), real_t(0), real_t(0));
+
+   uint_t peSubCycles( 10 );
+   real_t dtPe( real_t(10) );
+   real_t dtPeSubCycle = dtPe / real_c(peSubCycles);
+
+   pe_coupling::TimeStep timestep(blocks, bodyStorageID, cr, syncCall, dtPe, peSubCycles);
+
+   // evaluate how far the sphere will travel with a specific force applied which is the reference distance for later
+   // (depends on the chosen time integrator in the DEM and thus can not generally be computed a priori here)
+
+   Vector3<real_t> expectedPosOffset(real_t(0));
+   Vector3<real_t> expectedLinearVel(real_t(0));
+   Vector3<real_t> expectedAngularVel(real_t(0));
+   {
+
+      Vector3<real_t> startPos = positionInsideBlock;
+
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       startPos, radius, sphereMaterialID, false, true, false);
+
+      for( uint_t t = 0; t < peSubCycles; ++t )
+      {
+         for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+         {
+            for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+            {
+               bodyIt->addForceAtPos(testForce, bodyIt->getPosition() + torqueOffset);
+            }
+         }
+         cr.timestep(dtPeSubCycle);
+         syncCall();
+      }
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            expectedPosOffset = bodyIt->getPosition() - startPos;
+            expectedLinearVel = bodyIt->getLinearVel();
+            expectedAngularVel = bodyIt->getAngularVel();
+         }
+      }
+
+      mpi::allReduceInplace(expectedPosOffset[0], mpi::SUM);
+      mpi::allReduceInplace(expectedPosOffset[1], mpi::SUM);
+      mpi::allReduceInplace(expectedPosOffset[2], mpi::SUM);
+      mpi::allReduceInplace(expectedLinearVel[0], mpi::SUM);
+      mpi::allReduceInplace(expectedLinearVel[1], mpi::SUM);
+      mpi::allReduceInplace(expectedLinearVel[2], mpi::SUM);
+      mpi::allReduceInplace(expectedAngularVel[0], mpi::SUM);
+      mpi::allReduceInplace(expectedAngularVel[1], mpi::SUM);
+      mpi::allReduceInplace(expectedAngularVel[2], mpi::SUM);
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting position offset: " << expectedPosOffset);
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting linear vel: " << expectedLinearVel);
+      WALBERLA_LOG_DEVEL_ON_ROOT(" - expecting angular vel: " << expectedAngularVel);
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+   }
+
+   //////////////////
+   // Inside block //
+   //////////////////
+   {
+      std::string testIdentifier("Test: sphere inside block");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+      Vector3<real_t> startPos = positionInsideBlock;
+
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       startPos, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+
+      timestep();
+
+
+      Vector3<real_t> curPosOffset(real_t(0));
+      Vector3<real_t> curLinearVel(real_t(0));
+      Vector3<real_t> curAngularVel(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            curPosOffset = bodyIt->getPosition() - startPos;
+            curLinearVel = bodyIt->getLinearVel();
+            curAngularVel = bodyIt->getAngularVel();
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+
+         }
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   ///////////////////////
+   // At block border 1 //
+   ///////////////////////
+   {
+      std::string testIdentifier("Test: sphere at block border 1");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+      Vector3<real_t> startPos = positionAtBlockBorder;
+
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       startPos, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+
+      timestep();
+
+
+      Vector3<real_t> curPosOffset(real_t(0));
+      Vector3<real_t> curLinearVel(real_t(0));
+      Vector3<real_t> curAngularVel(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            curPosOffset = bodyIt->getPosition() - startPos;
+            curLinearVel = bodyIt->getLinearVel();
+            curAngularVel = bodyIt->getAngularVel();
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+         }
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   ////////////////////////////
+   // At block border 1, mod //
+   ////////////////////////////
+   {
+      std::string testIdentifier("Test: sphere at block border 1 mod");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+      Vector3<real_t> startPos = positionAtBlockBorder;
+
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       startPos, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      // also add on shadow copy, but only half of it to have same total force/torque
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(real_t(0.5)*testForce, pos+torqueOffset);
+         }
+      }
+
+      timestep();
+
+
+      Vector3<real_t> curPosOffset(real_t(0));
+      Vector3<real_t> curLinearVel(real_t(0));
+      Vector3<real_t> curAngularVel(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            curPosOffset = bodyIt->getPosition() - startPos;
+            curLinearVel = bodyIt->getLinearVel();
+            curAngularVel = bodyIt->getAngularVel();
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+         }
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   ///////////////////////
+   // At block border 2 //
+   ///////////////////////
+   {
+      std::string testIdentifier("Test: sphere at block border 2");
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - started");
+
+      Vector3<real_t> startPos = positionAtBlockBorder2;
+
+      pe::createSphere(*globalBodyStorage, blocks->getBlockStorage(), bodyStorageID, 0,
+                       startPos, radius, sphereMaterialID, false, true, false);
+
+      syncCall();
+
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            auto pos = bodyIt->getPosition();
+            bodyIt->addForceAtPos(testForce, pos+torqueOffset);
+         }
+      }
+
+      timestep();
+
+
+      Vector3<real_t> curPosOffset(real_t(0));
+      Vector3<real_t> curLinearVel(real_t(0));
+      Vector3<real_t> curAngularVel(real_t(0));
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::LocalBodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::LocalBodyIterator::end(); ++bodyIt )
+         {
+            curPosOffset = bodyIt->getPosition() - startPos;
+            curLinearVel = bodyIt->getLinearVel();
+            curAngularVel = bodyIt->getAngularVel();
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[0], expectedPosOffset[0], "Mismatch in posOffset0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[1], expectedPosOffset[1], "Mismatch in posOffset1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curPosOffset[2], expectedPosOffset[2], "Mismatch in posOffset2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[0], expectedLinearVel[0], "Mismatch in linearVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[1], expectedLinearVel[1], "Mismatch in linearVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curLinearVel[2], expectedLinearVel[2], "Mismatch in linearVel2");
+
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[0], expectedAngularVel[0], "Mismatch in angularVel0");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[1], expectedAngularVel[1], "Mismatch in angularVel1");
+            WALBERLA_CHECK_FLOAT_EQUAL(curAngularVel[2], expectedAngularVel[2], "Mismatch in angularVel2");
+
+         }
+      }
+
+      // clean up
+      for( auto blockIt = blocks->begin(); blockIt != blocks->end(); ++blockIt )
+      {
+         for( auto bodyIt = pe::BodyIterator::begin( *blockIt, bodyStorageID); bodyIt != pe::BodyIterator::end(); ++bodyIt )
+         {
+            bodyIt->markForDeletion();
+         }
+      }
+      syncCall();
+
+      WALBERLA_LOG_DEVEL_ON_ROOT(testIdentifier << " - ended");
+
+   }
+
+   return 0;
+
+}
+
+} //namespace pe_sub_cycling_test
+
+int main( int argc, char **argv ){
+   pe_sub_cycling_test::main(argc, argv);
+}