Tutorial - Useful Features

This tutorial will introduce some useful features of the waLBerla framework which can make your life easier.

Checkpointing

You can checkpoint the current state of your rigid body dynamics simulation at any point to restore it afterwards. First you have to store the current domain partitioning using blockforest::BlockForest::saveToFile().

   auto forest = blockforest::createBlockForest( math::AABB(0,0,0,60,60,60),
                                                 Vector3<uint_t>(2,2,2),                   // number of blocks
                                                 Vector3<bool>(false, false, false));      // periodicity
   forest->saveToFile("SerializeDeserialize.sbf");

Then you have to store the current simulation data using domain_decomposition::BlockStorage::saveBlockData().

   forest->saveBlockData("SerializeDeserialize.dump", storageID);

This will store all non global rigid bodies to the file system.

To load everything again you start by creating the blockforest::BlockForest. This time you will use a different constructor.

   auto forest = make_shared< BlockForest >( uint_c( MPIManager::instance()->rank() ), "SerializeDeserialize.sbf", true, false );

Instead of initializing the Storage BlockDatum like you normally would

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

you have to use domain_decomposition::BlockStorage::loadBlockData()

   auto storageID    = forest->loadBlockData("SerializeDeserialize.dump", createStorageDataHandling<BodyTuple>(), "Storage");

Unfortunately due to a misorder in the loading scheme you have to reload your coarse collision detection.

   for (auto blockIt = forest->begin(); blockIt != forest->end(); ++blockIt)
   {
      ccd::ICCD* ccd = blockIt->getData< ccd::ICCD >( ccdID );
      ccd->reloadBodies();
   }

Hopefully this gets fixed in the future. ;)

Attention

This method does not save global bodies nor solver settings. You have to take care to restore these settings on your own.

A fully working example can be found in the SerializeDeserialize.cpp test of the pe module.

VTK Output

For VTK Output you have to create vtk::VTKOutput objects. To output the domain partitioning use vtk::createVTKOutput_DomainDecomposition.

   auto vtkDomainOutput = vtk::createVTKOutput_DomainDecomposition( forest, "domain_decomposition", 1, "vtk_out", "simulation_step" );

To output all sphere particles use vtk::createVTKOutput_PointData in conjunction with SphereVtkOutput:

   auto vtkSphereHelper = make_shared<SphereVtkOutput>(storageID, *forest) ;
   auto vtkSphereOutput = vtk::createVTKOutput_PointData(vtkSphereHelper, "Bodies", 1, "vtk_out", "simulation_step", false, false);

Currently only spheres are supported for VTK output but you can easily write your own SphereVtkOutput and adapt it to the body you like.

To actually write something to disc call vtk::VTKOutput::write():

         vtkDomainOutput->write( );
         vtkSphereOutput->write( );

You can call this every time step if you want. The files will be automatically numbered so that ParaView can generate an animation.

Loading from Config

You can specify a config file as the first command line parameter. To access it you can use the Environment::config() function. You can access subblocks of the config with config::Config::getBlock().

   auto cfg = env.config();
   if (cfg == nullptr) WALBERLA_ABORT("No config specified!");
   const Config::BlockHandle configBlock  = cfg->getBlock( "LoadFromConfig" );

To get values from the config call config::Config::getParameter():

   real_t radius = configBlock.getParameter<real_t>("radius", real_c(0.4) );

Certain tasks already have predefined loading functions. You can for example directly create a BlockForest from the config file.

   shared_ptr<BlockForest> forest = blockforest::createBlockForestFromConfig( configBlock );

The corresponding block in the config file looks like:

simulationCorner < -15, -15, 0 >;
simulationDomain < 12, 23, 34 >;
blocks < 3, 4, 5 >;
isPeriodic < 0, 1, 0 >;

Also the HardContact solver can be configured directly from the config file:

   BlockDataID blockDataID;
   cr::HCSITS hcsits( globalBodyStorage, forest, blockDataID, blockDataID, blockDataID);
   configure(configBlock, hcsits);

The config file looks like:

HCSITSmaxIterations 123;
HCSITSRelaxationParameter 0.123;
HCSITSErrorReductionParameter 0.123;
HCSITSRelaxationModelStr ApproximateInelasticCoulombContactByDecoupling;
globalLinearAcceleration < 1, -2, 3 >;

Timing

To get additional information where your application spends its time you can use the WcTimingTree. It will give you a hierarchical view of the time used. Usage example:

   //WcTimingTree tt;
 
   tt.start("Initial Sync");
   syncCallWithoutTT();
   syncCallWithoutTT();
   tt.stop("Initial Sync");

Before you output the information you should collect all the information from all the processes if you are running in parallel.

   auto temp = tt.getReduced( );
   WALBERLA_ROOT_SECTION()
   {
      std::cout << temp;
   }

Many built-in functions like solver or synchronization methods come with an additional parameter where you can specify your timing tree. They will then include detailed information in your timing tree.

SQLite Output

waLBerla also supports SQLite databases for simulation data output. This can come in handy in parallel simulations as well as in data analysis. To store information in a SQLite database you have to fill three property maps depending on the type of information you want to store.

   std::map< std::string, walberla::int64_t >     integerProperties;
   std::map< std::string, double >      realProperties;
   std::map< std::string, std::string > stringProperties;

You can then dump the information to disc. timing::TimingPool and timing::TimingTree already have predefined save functions so you do not have to extract all the information yourself and save it in the property array.

   WALBERLA_ROOT_SECTION()
   {
      auto runId = sqlite::storeRunInSqliteDB( sqlFile, integerProperties, stringProperties, realProperties );
      sqlite::storeTimingPoolInSqliteDB( sqlFile, runId, *tpReduced, "Timeloop" );
      sqlite::storeTimingTreeInSqliteDB( sqlFile, runId, tt, "TimingTree" );
   }

Image Output

Using the pe::raytracing::Raytracer you can generate images of the simulation for each timestep and output them as PNG files. Setup the raytracer by reading a config object and optionally supply a shading and a visibility function, as it is done in the second pe tutorial. The shading function will be called during rendering for each body to apply user defined coloring to bodies, the visibility function to determine if a body should be visible or not.

   std::function<ShadingParameters (const BodyID body)> customShadingFunction = [](const BodyID body) {
      if (body->getTypeID() == Sphere::getStaticTypeID()) {
         return processRankDependentShadingParams(body).makeGlossy();
      }
      return defaultBodyTypeDependentShadingParams(body);
   };
   Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
                       cfg->getBlock("Raytracing"),
                       customShadingFunction);

Alternatively it may also be setup entirely in code:

Lighting lighting(Vec3(-12, 12, 12),
                  Color(1, 1, 1),
                  Color(1, 1, 1),
                  Color(real_t(0.4), real_t(0.4), real_t(0.4)));
Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
                    640, 480,
                    real_t(49.13), 2,
                    Vec3(-25, 10, 10), Vec3(-5, 10, 10), Vec3(0, 0, 1),
                    lighting,
                    Color(real_t(0.1), real_t(0.1), real_t(0.1)),
                    radius,
                    customShadingFunction);

After the configuration is done, images can be generated each timestep by calling Raytracer::generateImage<BodyTuple>() which will be output to the specified directory.

         raytracer.generateImage<BodyTuple>(size_t(i), &tt);

To hide certain bodies during rendering, the visibility function will be called with a BodyID as its sole argument and should return true if the object is supposed to be visible, false if not:

// [...]
std::function<bool (const BodyID body)> customVisibilityFunction = [](const BodyID body) {
  if (body->getTypeID() == Plane::getStaticTypeID()) {
     // hide all planes
     return false;
  }
  return true;
};
Raytracer raytracer(forest, storageID, globalBodyStorage, ccdID,
                   cfg->getBlock("Raytracing"),
                   customShadingFunction,
                   customVisibilityFunction);

For an overview over predefined shading functions, visit the file ShadingFunctions.h. For further information see the documentation for the classes pe::raytracing::Raytracer, pe::raytracing::Lighting and pe::raytracing::Color, the pe::raytracing::ShadingParameters struct and ShadingFunctions.h file may also be useful.