diff --git a/apps/benchmarks/CMakeLists.txt b/apps/benchmarks/CMakeLists.txt
index 47f77911653f7c8eec3d68550711810a8606e153..a81420b475e6378a3f453ebd4ae16c2c0c895e8d 100644
--- a/apps/benchmarks/CMakeLists.txt
+++ b/apps/benchmarks/CMakeLists.txt
@@ -4,6 +4,7 @@ add_subdirectory( CouetteFlow )
add_subdirectory( ForcesOnSphereNearPlaneInShearFlow )
add_subdirectory( NonUniformGrid )
add_subdirectory( MotionSingleHeavySphere )
+add_subdirectory( PeriodicGranularGas )
add_subdirectory( PoiseuilleChannel )
add_subdirectory( SchaeferTurek )
-add_subdirectory( UniformGrid )
\ No newline at end of file
+add_subdirectory( UniformGrid )
diff --git a/apps/benchmarks/PeriodicGranularGas/CMakeLists.txt b/apps/benchmarks/PeriodicGranularGas/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..d26f6314fd3d73bbac01b460da7240949a45350a
--- /dev/null
+++ b/apps/benchmarks/PeriodicGranularGas/CMakeLists.txt
@@ -0,0 +1,6 @@
+waLBerla_link_files_to_builddir( *.cfg )
+waLBerla_link_files_to_builddir( *.py )
+
+waLBerla_add_executable ( NAME PeriodicGranularGas
+ FILES PeriodicGranularGas.cpp
+ DEPENDS blockforest core pe )
diff --git a/apps/benchmarks/PeriodicGranularGas/PeriodicGranularGas.cfg b/apps/benchmarks/PeriodicGranularGas/PeriodicGranularGas.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..e281e9b0ecd0abe2ee48313040696d09d86d9db1
--- /dev/null
+++ b/apps/benchmarks/PeriodicGranularGas/PeriodicGranularGas.cfg
@@ -0,0 +1,22 @@
+
+PeriodicGranularGas
+{
+ simulationCorner < 0, 0, 0 >;
+ simulationDomain < 30, 30, 30 >;
+ blocks < 2, 2, 2 >;
+ isPeriodic < 1, 1, 1 >;
+
+ radius 0.5;
+ spacing 1.0;
+ vMax 0.0;
+
+ dt 0.1;
+ simulationSteps 1000;
+ visSpacing 100;
+
+ HCSITSmaxIterations 10;
+ HCSITSRelaxationParameter 0.7;
+ HCSITSErrorReductionParameter 0.8;
+ HCSITSRelaxationModelStr ApproximateInelasticCoulombContactByDecoupling;
+ globalLinearAcceleration < 0, 0, 0 >;
+}
diff --git a/apps/benchmarks/PeriodicGranularGas/PeriodicGranularGas.cpp b/apps/benchmarks/PeriodicGranularGas/PeriodicGranularGas.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7efaeb660c56c92b7467a83b17312125fa403d19
--- /dev/null
+++ b/apps/benchmarks/PeriodicGranularGas/PeriodicGranularGas.cpp
@@ -0,0 +1,250 @@
+//======================================================================================================================
+//
+// 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 PeriodicGranularGas.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+#include <pe/basic.h>
+#include <pe/vtk/SphereVtkOutput.h>
+
+#include <core/Abort.h>
+#include <core/Environment.h>
+#include <core/grid_generator/SCIterator.h>
+#include <core/logging/Logging.h>
+#include <core/waLBerlaBuildInfo.h>
+#include <vtk/VTKOutput.h>
+
+#include <functional>
+#include <memory>
+
+using namespace walberla;
+using namespace walberla::pe;
+using namespace walberla::timing;
+
+typedef boost::tuple<Sphere> BodyTuple ;
+
+int main( int argc, char ** argv )
+{
+ Environment env(argc, argv);
+
+ logging::Logging::instance()->setStreamLogLevel(logging::Logging::INFO);
+ logging::Logging::instance()->setFileLogLevel(logging::Logging::INFO);
+
+ WALBERLA_LOG_INFO_ON_ROOT( "config file: " << argv[1] )
+ WALBERLA_LOG_INFO_ON_ROOT( "waLBerla Revision: " << WALBERLA_GIT_SHA1 );
+
+ math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** READING COMMANDLINE ARGUMENTS ***");
+ bool bDEM = false;
+ bool bHCSITS = false;
+
+ bool bNN = false;
+ bool bSO = false;
+
+ bool bInelasticFrictionlessContact = false;
+ bool bApproximateInelasticCoulombContactByDecoupling = false;
+ bool bInelasticCoulombContactByDecoupling = false;
+ bool bInelasticGeneralizedMaximumDissipationContact = false;
+
+ for( int i = 1; i < argc; ++i )
+ {
+ if( std::strcmp( argv[i], "--DEM" ) == 0 ) bDEM = true;
+ if( std::strcmp( argv[i], "--HCSITS" ) == 0 ) bHCSITS = true;
+
+ if( std::strcmp( argv[i], "--syncNextNeighbor" ) == 0 ) bNN = true;
+ if( std::strcmp( argv[i], "--syncShadowOwners" ) == 0 ) bSO = true;
+
+ if( std::strcmp( argv[i], "--InelasticFrictionlessContact" ) == 0 ) bInelasticFrictionlessContact = true;
+ if( std::strcmp( argv[i], "--ApproximateInelasticCoulombContactByDecoupling" ) == 0 ) bApproximateInelasticCoulombContactByDecoupling = true;
+ if( std::strcmp( argv[i], "--InelasticCoulombContactByDecoupling" ) == 0 ) bInelasticCoulombContactByDecoupling = true;
+ if( std::strcmp( argv[i], "--InelasticGeneralizedMaximumDissipationContact" ) == 0 ) bInelasticGeneralizedMaximumDissipationContact = true;
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** READING CONFIG FILE ***");
+ auto cfg = env.config();
+ if (cfg == NULL) WALBERLA_ABORT("No config specified!");
+ const Config::BlockHandle mainConf = cfg->getBlock( "PeriodicGranularGas" );
+
+ int simulationSteps = mainConf.getParameter<int>("simulationSteps", 10 );
+ WALBERLA_LOG_INFO_ON_ROOT("simulationSteps: " << simulationSteps);
+
+ real_t dt = mainConf.getParameter<real_t>("dt", real_c(0.01) );
+ WALBERLA_LOG_INFO_ON_ROOT("dt: " << dt);
+
+ const int visSpacing = mainConf.getParameter<int>("visSpacing", 1000 );
+ WALBERLA_LOG_INFO_ON_ROOT("visSpacing: " << visSpacing);
+ const std::string path = mainConf.getParameter<std::string>("path", "vtk_out" );
+ WALBERLA_LOG_INFO_ON_ROOT("path: " << path);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** GLOBALBODYSTORAGE ***");
+ shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** BLOCKFOREST ***");
+ // create forest
+ shared_ptr< BlockForest > forest = createBlockForestFromConfig( mainConf );
+ if (!forest)
+ {
+ WALBERLA_LOG_INFO_ON_ROOT( "No BlockForest created ... exiting!");
+ return EXIT_SUCCESS;
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("simulationDomain: " << forest->getDomain());
+
+ WALBERLA_LOG_INFO_ON_ROOT("blocks: " << Vector3<uint_t>(forest->getXSize(), forest->getYSize(), forest->getZSize()) );
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** BODYTUPLE ***");
+ // initialize body type ids
+ SetBodyTypeIDs<BodyTuple>::execute();
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** STORAGEDATAHANDLING ***");
+ // add block data
+ auto storageID = forest->addBlockData(createStorageDataHandling<BodyTuple>(), "Storage");
+ auto ccdID = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
+ auto fcdID = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTuple, fcd::AnalyticCollideFunctor>(), "FCD");
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** INTEGRATOR ***");
+ std::unique_ptr<cr::ICR> cr;
+ if (bDEM)
+ {
+ cr = std::make_unique<cr::DEM>(globalBodyStorage, forest, storageID, ccdID, fcdID);
+ WALBERLA_LOG_INFO_ON_ROOT("Using DEM!");
+ } else if (bHCSITS)
+ {
+ cr = std::make_unique<cr::HCSITS>(globalBodyStorage, forest, storageID, ccdID, fcdID);
+ configure(mainConf, *static_cast<cr::HCSITS*>(cr.get()));
+ WALBERLA_LOG_INFO_ON_ROOT("Using HCSITS!");
+
+ cr::HCSITS* hcsits = static_cast<cr::HCSITS*>(cr.get());
+
+ if (bInelasticFrictionlessContact)
+ {
+ hcsits->setRelaxationModel(cr::HCSITS::InelasticFrictionlessContact);
+ WALBERLA_LOG_INFO_ON_ROOT("Using InelasticFrictionlessContact!");
+ } else if (bApproximateInelasticCoulombContactByDecoupling)
+ {
+ hcsits->setRelaxationModel(cr::HCSITS::ApproximateInelasticCoulombContactByDecoupling);
+ WALBERLA_LOG_INFO_ON_ROOT("Using ApproximateInelasticCoulombContactByDecoupling!");
+ } else if (bInelasticCoulombContactByDecoupling)
+ {
+ hcsits->setRelaxationModel(cr::HCSITS::InelasticCoulombContactByDecoupling);
+ WALBERLA_LOG_INFO_ON_ROOT("Using InelasticCoulombContactByDecoupling!");
+ } else if (bInelasticGeneralizedMaximumDissipationContact)
+ {
+ hcsits->setRelaxationModel(cr::HCSITS::InelasticGeneralizedMaximumDissipationContact);
+ WALBERLA_LOG_INFO_ON_ROOT("Using InelasticGeneralizedMaximumDissipationContact!");
+ } else
+ {
+ WALBERLA_ABORT("Friction model could not be determined!");
+ }
+ } else
+ {
+ WALBERLA_ABORT("Model could not be determined!");
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SYNCCALL ***");
+ std::function<void(void)> syncCallWithoutTT;
+ if (bNN)
+ {
+ syncCallWithoutTT = std::bind( pe::syncNextNeighbors<BodyTuple>, boost::ref(*forest), storageID, static_cast<WcTimingTree*>(NULL), real_c(0.0), false );
+ WALBERLA_LOG_INFO_ON_ROOT("Using NextNeighbor sync!");
+ } else if (bSO)
+ {
+ syncCallWithoutTT = std::bind( pe::syncShadowOwners<BodyTuple>, boost::ref(*forest), storageID, static_cast<WcTimingTree*>(NULL), real_c(0.0), false );
+ WALBERLA_LOG_INFO_ON_ROOT("Using ShadowOwner sync!");
+ } else
+ {
+ WALBERLA_ABORT("Synchronization method could not be determined!");
+ }
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** VTK ***");
+ auto vtkDomainOutput = vtk::createVTKOutput_DomainDecomposition( forest, "domain_decomposition", 1, "vtk_out", "simulation_step" );
+ auto vtkSphereHelper = make_shared<SphereVtkOutput>(storageID, *forest) ;
+ auto vtkSphereOutput = vtk::createVTKOutput_PointData(vtkSphereHelper, "Bodies", 1, "vtk_out", "simulation_step", false, false);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - START ***");
+ const real_t static_cof ( real_c(0.1) / 2 ); // Coefficient of static friction. Note: pe doubles the input coefficient of friction for material-material contacts.
+ const real_t dynamic_cof ( static_cof ); // Coefficient of dynamic friction. Similar to static friction for low speed friction.
+ MaterialID material = createMaterial( "granular", real_t( 1.0 ), 0, static_cof, dynamic_cof, real_t( 0.5 ), 1, 1, 0, 0 );
+
+ auto simulationDomain = forest->getDomain();
+ auto generationDomain = simulationDomain; // simulationDomain.getExtended(-real_c(0.5) * spacing);
+
+ const real_t spacing(1.0);
+ const real_t radius(0.5);
+ uint_t numParticles = uint_c(0);
+
+ for (auto& currentBlock : *forest)
+ {
+ for (auto it = grid_generator::SCIterator(currentBlock.getAABB().getIntersection(generationDomain), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5), spacing); it != grid_generator::SCIterator(); ++it)
+ {
+ SphereID sp = pe::createSphere( *globalBodyStorage, *forest, storageID, 0, *it, radius, material);
+ if (sp != NULL) ++numParticles;
+ }
+ }
+ mpi::reduceInplace(numParticles, mpi::SUM);
+ WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
+
+ // synchronize particles
+ syncCallWithoutTT();
+ syncCallWithoutTT();
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
+ WALBERLA_MPI_BARRIER();
+ WcTimer timer;
+ for (int i=0; i < simulationSteps; ++i)
+ {
+ if( i % 200 == 0 )
+ {
+ WALBERLA_LOG_DEVEL_ON_ROOT( "Timestep " << i << " / " << simulationSteps );
+ }
+
+ cr->timestep( real_c(dt) );
+ syncCallWithoutTT();
+
+// if( i % visSpacing == 0 )
+// {
+// vtkDomainOutput->write( );
+// vtkSphereOutput->write( );
+// }
+ }
+ WALBERLA_MPI_BARRIER();
+ timer.end();
+ WALBERLA_LOG_INFO_ON_ROOT("runtime: " << timer.average());
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - END ***");
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** CHECKING RESULT - START ***");
+ for (auto& currentBlock : *forest)
+ {
+ Storage * storage = currentBlock.getData< Storage >( storageID );
+ BodyStorage& localStorage = (*storage)[0];
+
+ auto bodyIt = localStorage.begin();
+ for (auto it = grid_generator::SCIterator(currentBlock.getAABB().getIntersection(generationDomain), Vector3<real_t>(spacing, spacing, spacing) * real_c(0.5), spacing);
+ it != grid_generator::SCIterator();
+ ++it, ++bodyIt)
+ {
+ WALBERLA_CHECK_UNEQUAL(bodyIt, localStorage.end());
+ WALBERLA_CHECK_FLOAT_EQUAL(bodyIt->getPosition(), *it);
+ }
+ }
+ WALBERLA_LOG_INFO_ON_ROOT("*** CHECKING RESULT - END ***");
+
+ return EXIT_SUCCESS;
+}