Tensor Train Accelerated Solvers for Nonsmooth Rigid Body Dynamics

Abstract

In the last two decades, increased need for high-fidelity simulations of the time evolution and propagation of forces in granular media has spurred a renewed interest in the discrete element method (DEM) modeling of frictional contact. Force penalty methods, while economic and widely accessible, introduce artificial stiffness, requiring small time steps to retain numerical stability. Optimization-based methods, which enforce contacts geometrically through complementarity constraints leading to a differential variational inequality problem (DVI), allow for the use of larger time steps at the expense of solving a nonlinear complementarity problem (NCP) each time step. We review the latest efforts to produce solvers for this NCP, focusing on its relaxation to a cone complementarity problem (CCP) and solution via an equivalent quadratic optimization problem with conic constraints. We distinguish between first order methods, which use only gradient information and are thus linearly convergent and second order methods, which rely on a Newton type step to gain quadratic convergence and are typically more robust and problem-independent. However, they require the approximate solution of large sparse linear systems, thus losing their competitive advantages in large scale problems due to computational cost.In this work, we propose a novel acceleration for the solution of Newton step linear systems in second order methods using low-rank compression based fast direct solvers, leveraging on recent direct solver techniques for structured linear systems arising from differential and integral equations. We employ the Quantized Tensor Train (QTT) decomposition to produce efficient approximate representations of the system matrix and its inverse. This provides a versatile and robust framework to accelerate its solution using this inverse in a direct or a preconditioned iterative method. We demonstrate compressibility of the Newton step matrices in Primal Dual Interior Point (PDIP) methods as applied to the multibody dynamics problem. Using a number of numerical tests, we demonstrate that this approach displays sublinear scaling of precomputation costs, may be efficiently updated across Newton iterations as well as across simulation time steps, and leads to a fast, optimal complexity solution of the Newton step. This allows our method to gain an order of magnitude speedups over state-of-the-art preconditioning techniques for moderate to large-scale systems, hence mitigating the computational bottleneck of second order methods.

On convex numerical schemes for inelastic contacts with friction

“…As in the previous section, one can follow [6,18,19,39], and discretize the normal constraint using (8). Here, following again [27] and [5], we consider the Taylor expansion of the constraint (10) in order to stabilize the scheme. Then, to solve (21,25), a natural Euler-based time-stepping discretization is…”

“…Improving the available algorithms to solve the friction problem is still an active domain of research. One can cite for example the recent works [10,13] where the authors propose a method to accelerate the Newton step in second-order methods.…”

“…This allows the use of several existing optimization solvers. For example, it can be solved by accelerated projected gradient descent methods [13,28], various spectral and Krylov subspace methods [16] or interior point methods [3,10,20,37]. To avoid the convex relaxation of the problem, the authors in [2] propose to iterate a fixed point functional to solve the original LCP problem.…”

Tensor-Train Compression of Discrete Element Method Simulation Data

“…friction), leading to highfidelity and robust predictions. In recent years, substantial progress has been made to address the algorithmic challenges in DEM simulations: in fast optimization solvers [1,2,3,4], parallel computing frameworks [5,6,7,8,9,10], and multiscale modeling [11,12,13,14]. Unfortunately, for large-scale problems such as those in vehicle-terrain simulation (see Figure 1), each run can still take days or even weeks to complete.…”

“…For instance, the entry F (11) is represented as F {3} (3, 1, 2). Assuming a column-major order for multidimensional arrays, this is essentially a reshaping: F {3} = reshape(F,[4,2,2]).…”