Tutorial - Confined Gas: Setting up a simple MESA-PD simulation

This tutorial will help you to set up a simple simulation using the MESA-PD module. The scenario for this tutorial is a particle gas (spheres) confined by solid walls.

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

#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>

The MESA-PD framework is comprised of different core components:

• collision_detection: collision_detection holds algorithms for determining collision information
• data: data holds all datastructures which can be used
• domain: domain contains the coupling to the BlockForest domain partitioning
• kernel: kernel holds all algorithms which act on particles
• mpi: mpi holds all functionality needed for synchronization
• vtk: vtk holds utility functions for writing the simulation to disk

#include <mesa_pd/collision_detection/AnalyticContactDetection.h>
#include <mesa_pd/data/LinkedCells.h>
#include <mesa_pd/data/ParticleAccessorWithShape.h>
#include <mesa_pd/data/ParticleStorage.h>
#include <mesa_pd/data/ShapeStorage.h>
#include <mesa_pd/domain/BlockForestDomain.h>
#include <mesa_pd/kernel/DoubleCast.h>
#include <mesa_pd/kernel/ExplicitEuler.h>
#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
#include <mesa_pd/kernel/ParticleSelector.h>
#include <mesa_pd/kernel/SpringDashpot.h>
#include <mesa_pd/mpi/ContactFilter.h>
#include <mesa_pd/mpi/ReduceProperty.h>
#include <mesa_pd/mpi/SyncNextNeighbors.h>
#include <mesa_pd/mpi/notifications/ForceTorqueNotification.h>
#include <mesa_pd/vtk/ParticleVtkOutput.h>

First, 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 * walberla::mpi::MPIManager::instance()->worldRank()) );
 
   real_t spacing          = real_c(1.0);
   real_t radius           = real_c(0.4);
   real_t density          = real_t(2707);
   real_t vMax             = real_c(1.0);
   int    simulationSteps  = 10;
   real_t dt               = real_c(0.01);

We will use the BlockForest datastructure for the domain partitioning of our simulation. 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 each direction of periodicity!

   auto forest = blockforest::createBlockForest( AABB(0,0,0,40,40,40), // simulation domain
                                                 Vector3<uint_t>(2,2,2), // blocks in each direction
                                                 Vector3<bool>(false, false, false) // periodicity
                                                 );
   if (!forest)
   {
      WALBERLA_LOG_INFO_ON_ROOT( "No BlockForest created ... exiting!");
      return EXIT_SUCCESS;
   }

Next, we initialize all data structures that we need for the simulation. The data::ParticleStorage is the data container for all particles. The data::ShapeStorage stores information about the shape of the particles. This information is separated from the particle due to that fact that most of the time there are many similar particles in one simulation. With this approach you do not need to store the information for every particle but only once.

Since all kernels are written against an abstract interface to access the particle properties we also need an intermediate accessor class. This class is in most cases data::ParticleAccessorWithShape.

Finally, to improve the complexity of the collision detection step, we need a data::LinkedCells acceleration data structure.

   auto ps = std::make_shared<data::ParticleStorage>(100);
   auto ss = std::make_shared<data::ShapeStorage>();
   data::ParticleAccessorWithShape accessor(ps, ss);
   data::LinkedCells     lc(localDomain.getExtended(spacing), spacing );

Now we are almost ready to initialize the particles for our simulation. For this purpose we will write two helper functions to create spheres and walls.

data::ParticleStorage::iterator createWall( data::ParticleStorage& ps,
                                            data::ShapeStorage& ss,
                                            const Vec3& pos,
                                            const Vec3& normal )
{
   using namespace walberla::mesa_pd::data::particle_flags;
 
   auto pl = ps.create(true);
   pl->setPosition          ( pos );
   pl->setInteractionRadius ( std::numeric_limits<real_t>::infinity() );
   pl->setShapeID           ( ss.create<data::HalfSpace>( normal ) );
   pl->setOwner             ( walberla::mpi::MPIManager::instance()->rank() );
   pl->setType              ( 1 );
   set(pl->getFlagsRef(), INFINITE);
   set(pl->getFlagsRef(), FIXED);
   set(pl->getFlagsRef(), NON_COMMUNICATING);
   return pl;
}
 
data::ParticleStorage::iterator createSphere( data::ParticleStorage& ps,
                                              const Vec3& pos,
                                              const real_t& radius,
                                              const size_t shapeID)
{
   auto sp = ps.create();
   sp->setPosition          ( pos );
   sp->setInteractionRadius ( radius );
   sp->setShapeID           ( shapeID );
   sp->setOwner             ( walberla::MPIManager::instance()->rank() );
   sp->setType              ( 0 );
   return sp;
}

The position should be self explanatory. The interaction radius is used for a first check if two particles can collide. It is the radius of the bounding volume. The shapeID is the index into the data::ShapeStorage datastructure. The owner is the current owning process of this particle and the type identifies the class this particle belongs to. This type is later on used to determine which material parameters should be used for this particle.

   createWall(*ps, *ss, simulationDomain.minCorner(), Vec3(+1,0,0) );
   createWall(*ps, *ss, simulationDomain.maxCorner(), Vec3(-1,0,0) );
   createWall(*ps, *ss, simulationDomain.minCorner(), Vec3(0,+1,0) );
   createWall(*ps, *ss, simulationDomain.maxCorner(), Vec3(0,-1,0) );
   createWall(*ps, *ss, simulationDomain.minCorner(), Vec3(0,0,+1) );
   createWall(*ps, *ss, simulationDomain.maxCorner(), Vec3(0,0,-1) );

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.

   auto  sphereShape = ss->create<data::Sphere>( radius );
   ss->shapes[sphereShape]->updateMassAndInertia( density );
   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)
      {
         auto sp = createSphere( *ps, *it, radius, sphereShape);
         Vec3 rndVel(math::realRandom<real_t>(-vMax, vMax),
                     math::realRandom<real_t>(-vMax, vMax),
                     math::realRandom<real_t>(-vMax, vMax));
         sp->setLinearVelocity(rndVel);
         ++numParticles;
      }
   }
   walberla::mpi::reduceInplace(numParticles, walberla::mpi::SUM);
   WALBERLA_LOG_INFO_ON_ROOT("#particles created: " << numParticles);

Before we can run the simulation we have to set up all kernels that we will use.

   kernel::ExplicitEuler                 explicitEuler( dt );
   kernel::InsertParticleIntoLinkedCells ipilc;
   kernel::SpringDashpot                 dem( 2 );
   auto meffSphereSphere = real_t(0.5) * (real_c(4.0)/real_c(3.0) * math::pi) * radius * radius * radius * density;
   auto meffSphereWall   = real_t(1.0) * (real_c(4.0)/real_c(3.0) * math::pi) * radius * radius * radius * density;
   dem.setParametersFromCOR( 0, 0, real_t(0.9), real_t(20) * dt, meffSphereSphere );
   dem.setParametersFromCOR( 0, 1, real_t(0.9), real_t(20) * dt, meffSphereWall );
   collision_detection::AnalyticContactDetection              acd;
   kernel::DoubleCast                    double_cast;
   mpi::ContactFilter                    contact_filter;
   mpi::ReduceProperty                   RP;
   mpi::SyncNextNeighbors                SNN;

In this example we will use the explicit Euler integration method. We further need a kernel which updates our LinkedCell data structure and an interaction kernel. For the collision detection we select the analytic functions available within the framework. Together with all synchronization kernels we can now run the simulation.

Attention

Before you start the simulation loop you should synchronize once to make sure everything is in place.

   for (int i=0; i < simulationSteps; ++i)
   {
      if( i % 10 == 0 )
      {
         WALBERLA_LOG_PROGRESS_ON_ROOT( "Timestep " << i << " / " << simulationSteps );
      }
 
      lc.clear();
      ps->forEachParticle(true, kernel::SelectAll(), accessor, ipilc, accessor, lc);
 
      lc.forEachParticlePairHalf(true,
                                 kernel::SelectAll(),
                                 accessor,
                                 [&](const size_t idx1, const size_t idx2, auto& ac)
      {
         if (double_cast(idx1, idx2, ac, acd, ac ))
         {
            if (contact_filter(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), domain))
            {
               dem(acd.getIdx1(), acd.getIdx2(), ac, acd.getContactPoint(), acd.getContactNormal(), acd.getPenetrationDepth());
            }
         }
      },
      accessor );
 
      RP.operator()<ForceTorqueNotification>(*ps);
 
      ps->forEachParticle(true, kernel::SelectLocal(), accessor, explicitEuler, accessor);
 
      SNN(*ps, domain);
   }

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 data::ParticleStorage datastructure.

   auto meanVelocity = real_t(0);
 
   ps->forEachParticle(true,
                       kernel::SelectLocal(),
                       accessor,
                       [&meanVelocity](const size_t idx, auto& ac)
   {
      meanVelocity += length(ac.getLinearVelocity(idx));
   },
   accessor);
 
   walberla::mpi::reduceInplace(meanVelocity, walberla::mpi::SUM);
   meanVelocity /= real_c(numParticles);
   WALBERLA_LOG_INFO_ON_ROOT( "mean velocity: " << meanVelocity );

Congratulation! You successfully created your first rigid particle simulation.