From 92227e717ce45eb1da9335516e6793cad9c3ec95 Mon Sep 17 00:00:00 2001
From: Sebastian Eibl <sebastian.eibl@fau.de>
Date: Fri, 4 Oct 2019 17:21:27 +0200
Subject: [PATCH] [ADD] first tutorial for mesa_pd
---
apps/tutorials/mesa_pd/01_LennardJones.cpp | 105 ---------
apps/tutorials/mesa_pd/01_MESAPD.cpp | 253 +++++++++++++++++++++
apps/tutorials/mesa_pd/01_MESAPD.dox | 71 ++++++
apps/tutorials/mesa_pd/CMakeLists.txt | 6 +-
doc/Mainpage.dox | 4 +
5 files changed, 331 insertions(+), 108 deletions(-)
delete mode 100644 apps/tutorials/mesa_pd/01_LennardJones.cpp
create mode 100644 apps/tutorials/mesa_pd/01_MESAPD.cpp
create mode 100644 apps/tutorials/mesa_pd/01_MESAPD.dox
diff --git a/apps/tutorials/mesa_pd/01_LennardJones.cpp b/apps/tutorials/mesa_pd/01_LennardJones.cpp
deleted file mode 100644
index 6ed34ccab..000000000
--- a/apps/tutorials/mesa_pd/01_LennardJones.cpp
+++ /dev/null
@@ -1,105 +0,0 @@
-//======================================================================================================================
-//
-// 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 01_LennardJones.cpp
-//! \author Sebastian Eibl <sebastian.eibl@fau.de>
-//
-//======================================================================================================================
-
-#include <mesa_pd/vtk/ParticleVtkOutput.h>
-
-#include <mesa_pd/data/LinkedCells.h>
-#include <mesa_pd/data/ParticleAccessor.h>
-#include <mesa_pd/data/ParticleStorage.h>
-
-#include <mesa_pd/kernel/DoubleCast.h>
-#include <mesa_pd/kernel/ExplicitEuler.h>
-#include <mesa_pd/kernel/ForceLJ.h>
-#include <mesa_pd/kernel/InsertParticleIntoLinkedCells.h>
-#include <mesa_pd/kernel/ParticleSelector.h>
-#include <mesa_pd/kernel/VelocityVerlet.h>
-
-#include <core/Environment.h>
-#include <core/grid_generator/SCIterator.h>
-#include <core/logging/Logging.h>
-#include <core/math/Random.h>
-#include <vtk/VTKOutput.h>
-
-#include <iostream>
-#include <memory>
-
-using namespace walberla;
-using namespace walberla::mesa_pd;
-
-int main( int argc, char ** argv )
-{
- Environment env(argc, argv);
- WALBERLA_UNUSED(env);
- mpi::MPIManager::instance()->useWorldComm();
-
- //domain setup
- const real_t spacing = real_t(1.0);
- math::AABB domain( Vec3(real_t(-0.5), real_t(-0.5), real_t(-0.5)),
- Vec3(real_t(+0.5), real_t(+0.5), real_t(+0.5)));
- domain.scale(real_t(20));
-
- //init data structures
- auto storage = std::make_shared<data::ParticleStorage> (100);
- data::LinkedCells linkedCells(domain.getScaled(2), real_t(1));
- data::ParticleAccessor ac(storage);
-
- //initialize particles
- for (auto it = grid_generator::SCIterator(domain, Vec3(spacing, spacing, spacing) * real_c(0.5), spacing);
- it != grid_generator::SCIterator();
- ++it)
- {
- data::Particle&& p = *storage->create();
- p.getPositionRef() = (*it);
-
- p.getLinearVelocityRef() = Vec3(math::realRandom(real_t(-1), real_t(+1)),
- math::realRandom(real_t(-1), real_t(+1)),
- math::realRandom(real_t(-1), real_t(+1)));
- }
-
- auto vtkOutput = make_shared<mesa_pd::vtk::ParticleVtkOutput>(storage) ;
- auto vtkWriter = walberla::vtk::createVTKOutput_PointData(vtkOutput, "Bodies", 1, "vtk", "simulation_step", false, false);
-
- // Init kernels
- kernel::ForceLJ lj(1);
- lj.setSigma(0,0,real_t(1));
- lj.setEpsilon(0,0,real_t(1));
- kernel::InsertParticleIntoLinkedCells ipilc;
- kernel::VelocityVerletPreForceUpdate vvPreForce( real_t(0.01) );
- kernel::VelocityVerletPostForceUpdate vvPostForce( real_t(0.01) );
-
- for (auto timestep = 0; timestep < 1000; ++timestep)
- {
- WALBERLA_LOG_DEVEL(timestep);
- linkedCells.clear();
- storage->forEachParticle(true, kernel::SelectAll(), ac, ipilc, ac, linkedCells);
- storage->forEachParticle(true, kernel::SelectLocal(), ac, vvPreForce, ac);
- linkedCells.forEachParticlePairHalf(true, kernel::SelectAll(), ac, lj, ac);
- const real_t coeff = real_t(0.2);
- storage->forEachParticle(true,
- kernel::SelectLocal(),
- ac,
- [coeff](const size_t idx, auto& access){ access.setForce(idx, -coeff*access.getPosition(idx) + access.getForce(idx)); },
- ac);
- storage->forEachParticle(true, kernel::SelectLocal(), ac, vvPostForce, ac);
- vtkWriter->write();
- }
-
- return EXIT_SUCCESS;
-}
diff --git a/apps/tutorials/mesa_pd/01_MESAPD.cpp b/apps/tutorials/mesa_pd/01_MESAPD.cpp
new file mode 100644
index 000000000..41f5c6867
--- /dev/null
+++ b/apps/tutorials/mesa_pd/01_MESAPD.cpp
@@ -0,0 +1,253 @@
+//======================================================================================================================
+//
+// 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 01_ConfinedGas.cpp
+//! \author Sebastian Eibl <sebastian.eibl@fau.de>
+//
+//======================================================================================================================
+
+//! [walberlaIncludes]
+#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>
+//! [walberlaIncludes]
+
+//! [mesapdIncludes]
+#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/ExplicitEulerWithShape.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>
+//! [mesapdIncludes]
+
+namespace walberla {
+namespace mesa_pd {
+
+//! [CreationHelper]
+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;
+}
+//! [CreationHelper]
+
+int main( int argc, char ** argv )
+{
+ //! [Parameters]
+ 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);
+ //! [Parameters]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** BLOCKFOREST ***");
+ // create forest
+ //! [BlockForest]
+ 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;
+ }
+ //! [BlockForest]
+
+ auto simulationDomain = forest->getDomain();
+ auto localDomain = forest->begin()->getAABB();
+ for (auto& blk : *forest)
+ {
+ localDomain.merge(blk.getAABB());
+ }
+
+ auto domain = domain::BlockForestDomain(forest);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** DATA STRUCTURES ***");
+ //! [DataStructures]
+ 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 );
+ //! [DataStructures]
+
+ const auto& generationDomain = simulationDomain; // simulationDomain.getExtended(-real_c(0.5) * spacing);
+ //! [Walls]
+ 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) );
+ //! [Walls]
+
+ //! [Spheres]
+ 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);
+ //! [Spheres]
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SETUP - END ***");
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** KERNELS ***");
+ //! [Kernels]
+ kernel::ExplicitEulerWithShape explicitEulerWithShape( 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;
+ //! [Kernels]
+
+ SNN(*ps, domain);
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - START ***");
+ //! [Loop]
+ 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, explicitEulerWithShape, accessor);
+
+ SNN(*ps, domain);
+ }
+ //! [Loop]
+ WALBERLA_LOG_INFO_ON_ROOT("*** SIMULATION - END ***");
+
+ WALBERLA_LOG_INFO_ON_ROOT("*** GETTING STATISTICAL INFORMATION ***");
+ //! [PostProcessing]
+ 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 /= numParticles;
+ WALBERLA_LOG_INFO_ON_ROOT( "mean velocity: " << meanVelocity );
+ //! [PostProcessing]
+
+ return EXIT_SUCCESS;
+}
+
+} // namespace mesa_pd
+} // namespace walberla
+
+int main( int argc, char* argv[] )
+{
+ return walberla::mesa_pd::main( argc, argv );
+}
diff --git a/apps/tutorials/mesa_pd/01_MESAPD.dox b/apps/tutorials/mesa_pd/01_MESAPD.dox
new file mode 100644
index 000000000..c3c1d3c81
--- /dev/null
+++ b/apps/tutorials/mesa_pd/01_MESAPD.dox
@@ -0,0 +1,71 @@
+namespace walberla {
+namespace mesa_pd {
+
+/**
+\page tutorial_mesapd_01 Tutorial - Confined Gas: Setting up a simple MESA-PD simulation
+
+This tutorial will help you to setup up a simple simulation using the MESA-PD module.
+The scenario for this tutorial a is particle gas (spheres) confined by solid walls.
+
+For this tutorial only a few waLBerla includes are needed, namely:
+\snippet 01_MESAPD.cpp walberlaIncludes
+
+The MESA-PD framework is comprised of different core components:
+- collision_detection: %collision_detection holds algorithms for determining collision information
+- \ref walberla::mesa_pd::data "data": data holds all datastructures which can be used
+- \ref walberla::mesa_pd::domain "domain": domain contains the coupling to the BlockForest domain partitioning
+- \ref walberla::mesa_pd::kernel "kernel": kernel holds all algorithms which act on particles
+- \ref walberla::mesa_pd::mpi "mpi": mpi holds all functionality needed for synchronization
+- \ref walberla::mesa_pd::vtk "vtk": vtk holds utility functions for writing the simulation to disk
+
+\snippet 01_MESAPD.cpp mesapdIncludes
+
+First, the waLBerla environment is initialized, the random number generator is seeded and some simulation parameters are set.
+\snippet 01_MESAPD.cpp Parameters
+
+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 \ref tutorial_basics_01 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!
+
+\snippet 01_MESAPD.cpp BlockForest
+
+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 on 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 a intermediate accessor class. This class is in most cases is the data::ParticleAccessorWithShape.
+
+Finally, to improve the complexity of the collision detection step, we need a data::LinkedCells acceleration data structure.
+
+\snippet 01_MESAPD.cpp DataStructure
+
+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.
+
+\snippet 01_MESAPD.cpp CreationHelper
+
+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.
+
+\snippet 01_MESAPD.cpp Walls
+
+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.
+
+\snippet 01_MESAPD.cpp Spheres
+
+Before we can run the simulation we have to set up all kernels that we will use.
+\snippet 01_MESAPD.cpp Kernels
+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.
+
+\snippet 01_MESAPD.cpp Loop
+
+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.
+
+\snippet 01_MESAPD.cpp PostProcessing
+
+Congratulation! You successfully created your first rigid particle simulation.
+
+*/
+
+}
+}
diff --git a/apps/tutorials/mesa_pd/CMakeLists.txt b/apps/tutorials/mesa_pd/CMakeLists.txt
index e74127bea..6c49f6996 100644
--- a/apps/tutorials/mesa_pd/CMakeLists.txt
+++ b/apps/tutorials/mesa_pd/CMakeLists.txt
@@ -1,5 +1,5 @@
waLBerla_link_files_to_builddir( *.cfg )
-waLBerla_add_executable ( NAME 01_LennardJones
- FILES 01_LennardJones.cpp
- DEPENDS core mesa_pd )
+waLBerla_add_executable ( NAME 01_MESA_PD
+ FILES 01_MESAPD.cpp
+ DEPENDS blockforest core mesa_pd )
diff --git a/doc/Mainpage.dox b/doc/Mainpage.dox
index 38d1b8777..d70a28f33 100644
--- a/doc/Mainpage.dox
+++ b/doc/Mainpage.dox
@@ -20,6 +20,10 @@ all the basic data strcutures and concepts of the framework.
- \ref tutorial_basics_03 \n
The third tutorial deals with communication, iterating field/grid/lattice data, and geometry setup.
+\subsection mesapd Modular and Extensible Software Architecture for Particle Dynamics
+
+- \ref tutorial_mesapd_01 \n
+
\subsection pe Rigid Body Dynamics (pe)
- \ref tutorial_pe_01 \n
--
GitLab