Tutorial - Confined Gas: Setting up a simple pe simulation

This tutorial will help you to set up a simple pe simulation. The scenario for this tutorial is a particle gas (spheres) confined by solid walls.

For this tutorial only a few includes are needed, namely:

#include <pe/basic.h>
 
#include <blockforest/Initialization.h>
#include <core/Environment.h>
#include <core/grid_generator/HCPIterator.h>
#include <core/grid_generator/SCIterator.h>
#include <core/logging/Logging.h>
#include <core/math/Random.h>
 
#include <tuple>

The first important step is to define a BodyTuple. This tuple stores all body types which will be used during the simulation. This tuple helps to generate all the functions which need to be specialized for a particular body type.

Attention

If you miss a body type the behaviour of the simulation is undefined. Most probably there will be an error.

using BodyTypeTuple = std::tuple<Sphere, Plane> ;

Next the waLBerla environment is initialized, the random number generator is seeded and some simulation parameters are set.

   Environment env(argc, argv);
   WALBERLA_UNUSED(env);
 
   math::seedRandomGenerator( static_cast<unsigned int>(1337 * mpi::MPIManager::instance()->worldRank()) );
 
   real_t spacing          = real_c(1.0);
   real_t radius           = real_c(0.4);
   real_t vMax             = real_c(1.0);
   int    simulationSteps  = 10;
   real_t dt               = real_c(0.01);

The BlockForest is the main datastructure in the waLBerla framework. It is responsible for the domain decomposition and holds all the blocks with their data. For more information about the general design of the waLBerla framework please refer to Tutorial - Basics 1: waLBerla Data Structures: Blocks and Fields and the documentation of domain_decomposition::BlockStorage. You can choose the number of blocks you want to have in each direction. In a parallel simulation these blocks get assigned to different processes. You should make sure that you always have at least as many blocks as processes. The number of processes you want your simulation to run with is specified when you start your program with mpiexec.

Attention

If you run a simulation with periodic boundaries you need at least three blocks in the direction of periodicity!

   auto forest = blockforest::createBlockForest( AABB(0,0,0,20,20,20), // simulation domain
                                                 Vector3<uint_t>(2,2,2), // blocks in each direction
                                                 Vector3<bool>(false, false, false) // periodicity
                                                 );

There are two types of storages to store particle information. One is the global body storage which is responsible for very large particles. These particles are stored on all processes.

   shared_ptr<BodyStorage> globalBodyStorage = make_shared<BodyStorage>();

The other one is a block local storage for all bodies which are located within the subdomain of a block. You can find more information about this concept of data on blocks in Tutorial - Basics 1: waLBerla Data Structures: Blocks and Fields.

   auto storageID           = forest->addBlockData(createStorageDataHandling<BodyTypeTuple>(), "Storage");

In addition to the block local storage also the coarse as well as the fine collision detection needs storage on each block. Therefore the corresponding data handling has to be registered.

   auto ccdID               = forest->addBlockData(ccd::createHashGridsDataHandling( globalBodyStorage, storageID ), "CCD");
   auto fcdID               = forest->addBlockData(fcd::createGenericFCDDataHandling<BodyTypeTuple, fcd::AnalyticCollideFunctor>(), "FCD");

Only one final component is missing for a successful simulation - the time integrator. Currently there exists two integrators cr::DEM for soft contacts and cr::HCSITS for hard contacts. These have to be setup as local objects.

   cr::HCSITS cr(globalBodyStorage, forest, storageID, ccdID, fcdID);
   cr.setMaxIterations( 10 );
   cr.setRelaxationModel( cr::HardContactSemiImplicitTimesteppingSolvers::ApproximateInelasticCoulombContactByDecoupling );
   cr.setRelaxationParameter( real_t(0.7) );
   cr.setGlobalLinearAcceleration( Vec3(0,0,0) );

We now have all our tools together and can start creating the simulation. First we have to specify the material parameters we want to use.

   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 );

And we have to initialize all body type ids.

   SetBodyTypeIDs<BodyTypeTuple>::execute();

Then we can set up the confining walls.

   createPlane(*globalBodyStorage, 0, Vec3(1,0,0), simulationDomain.minCorner(), material );
   createPlane(*globalBodyStorage, 0, Vec3(-1,0,0), simulationDomain.maxCorner(), material );
   createPlane(*globalBodyStorage, 0, Vec3(0,1,0), simulationDomain.minCorner(), material );
   createPlane(*globalBodyStorage, 0, Vec3(0,-1,0), simulationDomain.maxCorner(), material );
   createPlane(*globalBodyStorage, 0, Vec3(0,0,1), simulationDomain.minCorner(), material );
   createPlane(*globalBodyStorage, 0, Vec3(0,0,-1), simulationDomain.maxCorner(), material );

The gas particles are generated with the help of grid generators. These generators are iterators which facilitate the construction of particles on a regular grid. Currently grid_generator::SCIterator (for a simple cubic lattice) as well as grid_generator::HCPIterator (for a hexagonal close packing lattice) are supported. The spheres are create with createSphere which returns a SphereID of the created sphere. This SphereID acts like a pointer to Sphere and can be used like that. If for various reasons the sphere was not created the return value is NULL.

Attention

Before accessing the underlying sphere you should always check for NULL!

After you have initialized all particles you should synchronize the simulation to make sure that all information is distributed correctly. Two synchronization methods are available syncNextNeighbors() and syncShadowOwners().

   uint_t numParticles = uint_c(0);
   for (auto blkIt = forest->begin(); blkIt != forest->end(); ++blkIt)
   {
      IBlock & currentBlock = *blkIt;
      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 = createSphere( *globalBodyStorage, *forest, storageID, 0, *it, radius, material);
         Vec3 rndVel(math::realRandom<real_t>(-vMax, vMax), math::realRandom<real_t>(-vMax, vMax), math::realRandom<real_t>(-vMax, vMax));
         if (sp != nullptr) sp->setLinearVel(rndVel);
         if (sp != nullptr) ++numParticles;
      }
   }
   WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);
   syncNextNeighbors<BodyTypeTuple>(*forest, storageID);

Since the setup is finished now we can run the simulation loop. The simulation loop is as simple as:

   for (int i=0; i < simulationSteps; ++i)
   {
      if( i % 10 == 0 )
      {
         WALBERLA_LOG_PROGRESS_ON_ROOT( "Timestep " << i << " / " << simulationSteps );
      }
 
      cr.timestep( real_c(dt) );
      syncNextNeighbors<BodyTypeTuple>(*forest, storageID);
   }

cr::ICR::timestep() evolves your simulation in time. The subsequent synchronization keeps all particles that are known to more than one process in sync.

After the simulation is finished we can collect the results. In this case we only calculate the mean velocity of all particles. The particles can be easily accessed via the LocalBodyIterator. This iterator allows us to iterate through all particles and access their properties.

   Vec3 meanVelocity(0,0,0);
   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
   {
      for (auto bodyIt = LocalBodyIterator::begin(*blockIt, storageID); bodyIt != LocalBodyIterator::end(); ++bodyIt)
      {
        meanVelocity += bodyIt->getLinearVel();
      }
   }
   meanVelocity /= numParticles;
   WALBERLA_LOG_INFO( "mean velocity: " << meanVelocity );

Congratulation! You successfully created your first rigid body simulation. In the next tutorial we will look at possible extensions which can make your life easier.