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:

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.

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

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

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.

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.

PlaneID createPlane(BodyStorage &globalStorage, id_t uid, Vec3 normal, const Vec3 &gpos, MaterialID material)
Setup of a new Plane.
Definition: PlaneFactory.cpp:35
MaterialID createMaterial(const std::string &name, real_t density, real_t cor, real_t csf, real_t cdf, real_t poisson, real_t young, real_t stiffness, real_t dampingN, real_t dampingT)
Creating a new custom material.
Definition: Materials.cpp:161
static iterator< pe::BodyStorage::iterator > begin(IBlock &block, const BlockDataID &bodyStorageId)
Definition: BodyIterators.h:203
shared_ptr< BlockForest > createBlockForest(const AABB &domainAABB, const uint_t numberOfXBlocks, const uint_t numberOfYBlocks, const uint_t numberOfZBlocks, const uint_t numberOfXProcesses, const uint_t numberOfYProcesses, const uint_t numberOfZProcesses, const bool xPeriodic, const bool yPeriodic, const bool zPeriodic, const bool keepGlobalBlockInformation)
Function for creating a block forest that represents a uniform block grid.
Definition: Initialization.cpp:203
static iterator< pe::BodyStorage::iterator > end()
Definition: BodyIterators.h:219
SphereID createSphere(BodyStorage &globalStorage, BlockStorage &blocks, BlockDataID storageID, id_t uid, const Vec3 &gpos, real_t radius, MaterialID material, bool global, bool communicating, bool infiniteMass)
Setup of a new Sphere.
Definition: SphereFactory.cpp:35
Definition: Logging.h:669
Definition: Logging.h:697
void seedRandomGenerator(const std::mt19937::result_type &seed)
Definition: Random.cpp:43
Definition: DataTypes.h:232
math::Vector3< real_t > Vec3
Definition: Types.h:46
std::tuple< Sphere, Plane > BodyTypeTuple
Definition: 01_ConfinedGas.cpp:38
std::size_t uint_t
Definition: DataTypes.h:133
HardContactSemiImplicitTimesteppingSolvers HCSITS
Definition: HCSITS.h:273
float real_t
Definition: DataTypes.h:167
static void execute()
Initial setup of static type ids.
Definition: SetBodyTypeIDs.h:50
shared_ptr< HashGridsDataHandling > createHashGridsDataHandling(const shared_ptr< BodyStorage > &globalStorage, const BlockDataID &storageID)
Definition: HashGridsDataHandling.h:47
Materials::size_type MaterialID
Unique material ID.
Definition: Types.h:159
Header file which includes common pe headers!
GenericAABB< real_t > AABB
Definition: AABBFwd.h:33
Initialization of waLBerla.
Sphere * SphereID
Handle for a sphere primitive.
Definition: Types.h:119
@ ApproximateInelasticCoulombContactByDecoupling
Definition: HCSITS.h:104
#define WALBERLA_LOG_INFO(msg)
Definition: Logging.h:657