What's behind MWDiEM?
- The time integration of the Newton equations was achieved by a quaternion based Verlet algorithm [1]. As this part of the algorithm always acts on nearly if not all the particles separately, this algorithm was ported to the Graphical Processing Unit (GPU). Other tested algorithm included another version of the Verlet algorithm, and a Rodrigues rotation vector method instead of quaternions, however these have proven to be less efficient [2].
- The final version of the implemented contact detection algorithm has three subparts. First a coarse contact detection is made, identifying possible interaction pairs and is based on a special type of hashing [3]. Second after considering and testing algorithms based on the discretization of parts of the polyhedra onto a finite mesh, and different types of algorithms based on the separating axis theorem, the final choice fell on a Gilbert–Johnson–Keerthi based algorithm [4], that was slightly modified to increase its efficiency even more. This choice of this algorithm makes it even possible to later easily rewrite the code to be able to handle any convex particles, not just polyhedra. This algorithm was also ported to the GPU. The third phase is the determination of the penetration depth and penetration normal. After careful consideration an Expanding Polytope Algorithm (EPA) was chosen. This part was also ported to the GPU, and currently only determines the penetration depth. In the future it will be redesigned to provide the full penetration volume, resulting in more exact force estimates that will lead to even better simulation results.
- The particle shapes possible range from tetrahedrons to arbitrary polyhedra and the option of adding spheres into the system is also possible. As mentioned above, with minor modification the code could be easily modified to model arbitrary convex shapes altogether.
- Possible fracture of the particles was implemented. Breaking threshold can be set manually or depending on the local shear rate. From simply halving the particles, to it losing a small chunk or even pulverizing in one single instance are all possible, however one should note that this gravely influences the performance of the code as it drasticly increases the number of the particles.
- The code can also efficiently run on a good PC with a higher end GPU card, which makes testing simulations especially convenient. For smaller systems it can be used without a GPU card, solely using the CPU of a PC, this is thanks to the fact that all algorithms were first developed to work on the CPU and only then were they rewritten to work on the GPU.
- The code is prepared for Geographic Information Systems (GIS) implementation, meaning that it can import raster files for the base geometry, and also the full system topography can be extracted at any time from the simulations to be able to be loaded into a GIS program.
- The use of different basal topographies is implemented. Bases made of particles “glued together” (meaning constraints between sets of particles), or bases consisting of a large number of planes, or triangulated planes is possible. With this even very complex basal or system topographies are possible.
- Basal entrainment is also possible by setting the force criteria at which particles dislocate from the base.
- Obstacles may be included in the path of the movement (e.g. a dam or columns) with the same method as in point (2).
- Different slide initiation/slope failure methods are implemented, like dams, vibrations or forced breakup of a block of mass.
- Possible to use data acquired from r.avaflow as an initial condition.
- The inter-particle force law is currently a linear spring-dashpot based contact model [3], however the code is written in a modulated way, so other contact theories (e.g. with cohesion or for very high speed collisions found in falls or very rapid, almost gaseous flows) are easily applicable in the future.
- Most of the visualization is based on the Paraview and Visit packages.
The code is uploaded to Github, is open source and freely downloadable for anyone. A simple manual has been prepared there. Even though the code is ready for use, there are number of points that could make the code even more efficient and accurate and will be developed in the future. One is rewriting the EPA algorithm part to not just calculate the penetration normal, but the full overlap volume for more precise force calculations, while the other considers more efficient GPU-CPU memory management tasks, as currently these take the up most of the simulation time, and a third one will be implementing more efficient wall-particle interactions.