An adaptive hierarchical domain decomposition method for parallel contact dynamics simulations of granular materials

Abstract

A fully parallel version of the Contact Dynamics (CD) method is presented in this paper. For large enough systems, 100% efficiency has been demonstrated for up to 256 processors using a hierarchical domain decomposition with dynamic load balancing. The iterative scheme to calculate the contact forces is left domain-wise sequential, with data exchange after each iteration step, which ensures its stability. The number of additional iterations required for convergence by the partially parallel updates at the domain boundaries becomes negligible with increasing number of particles, which allows for an effective parallelization. Compared to the sequential implementation, we found no influence of the parallelization on simulation results. decomposition leads to causality problems. The algorithm presented in [15] conserves causality by reverting to an older state when violated. The best efficiency reached so far is a speedup proportional to the square root of the number of processors [15].In contrast to ED, lasting contacts between rigid bodies are considered in the realm of (multi)-rigid-body dynamics. Common to all its realizations is the treatment of contact forces as constraint forces, preventing interpenetration and, to a certain extent in the case of frictional contacts, sliding. When applying the rigid body modelling to problems like e.g. robotics [16,17] or granular media [18,19,20,21], different algorithms can in principle be used. Approximations with respect to the constraint of dry Coulomb friction enable the usage of powerful standard techniques for linear complementary problems (LCP) [22]. Other algorithms keep the isotropic friction cone, using a solver based on a modified time stepping scheme leading to a cone complementary problem (CCP) for the simulation of frictional contact dynamics [23]. Other approximations, leading to fast frictional dynamics (FFD) [24], yield a computational cost being only linear in the number of contacts and thus allow for impressive system sizes in terms of the number of particles. For investigations of e.g. the stress field in granular media, these approximations are prohibitive, though, and thus the non-smooth contact dynamics method [7], or commonly just contact dynamics, is widely employed. We will sketch the principle of this iterative procedure in section 2.1. Parallelization of the FFD method is straightforward and efficient [25,26], on the other hand, the parallel version suffers also from the undesired approximations. The parallel implementation of the CCP algorithm by the use of the Graphics Processing Unit (GPU) for large-scale multibody dynamics simulations is presented in [27]. In the present work we investigate the impact of the parallelization on the numerical solution of the CD method going beyond [25,26,27].Providing a parallel CD code is motivated by the need for large-scale simulations of dense granular systems of hard particles. The computation time even scales as O(N 1+2/d ) with the number of particles in CD [8] (d is the dimension of the system), while it...

A parallel version of the contact dynamics method

“…The efficiency break down at N p =24 is attributed to the special architecture of the distributed memory cluster used for simulations (Cray-XT6m with 24 cores per board), meaning that the speed of the inter-board communications is much slower than the intra-board one. We also note that the impact of paralellization on the physical properties of the system have been carefully investigated [13], and no significant signature of parallel update is observed.…”

“…Finally, we test the efficiency of the parallel algorithm in homogeneous dense packings of 500, 8000, and 10 6 particles [13]. We mainly focus on the performance of the contact force iteration loop as the most time consuming part of the contact dynamics simulation.…”

“…A parallel update scheme (where the contact forces are updated based on the old information, and the new forces are replaced only at the end of the iteration step) is unstable. A mixture of the two schemes requires more iteration steps to achieve a given accuracy level compared with the sequential scheme, and becomes unstable when the fraction of sequentially updated contacts is less than 30−40% of the contacts [13]. This is an important point for parallelization of the CD method because the domain decomposition method allows for a sequential update only in the bulk of each domain, while the boundary contacts are updated in a parallel way.…”

High performance computing of discrete nonsmooth contact dynamics with domain decomposition

“…To address these issues, several answers had been proposed, for instance, several time stepping schemes [1,13], and parallelization of solvers [16,27,28,29].…”

“…Up to our knowledge, the only solution methods coupling a domain decomposition method with non-regular implicit contact dynamics for granular models are available in [5,23,29]. The approach described in [23] mainly uses an asynchronous distributed-memory Non-Linear Gauss-Seidel solver (close to the one used in [17] for a synchronous solver and in [18] for an asynchronous version, both for shared-memory architectures.…”