## 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.

## References

References
1.
Cundall
,
P. A.
, and
Strack
,
O. D.
,
1979
, “
A Discrete Numerical Model for Granular Assemblies
,”
Geotechnique
,
29
(
1
), pp.
47
65
.10.1680/geot.1979.29.1.47
2.
Pazouki
,
A.
,
Kwarta
,
M.
,
Williams
,
K.
,
Likos
,
W.
,
Serban
,
R.
,
Jayakumar
,
P.
, and
Negrut
,
D.
,
2017
, “
Compliant Contact Versus Rigid Contact: A Comparison in the Context of Granular Dynamics
,”
Phys. Rev. E
,
96
(
4
), p.
042905
.10.1103/PhysRevE.96.042905
3.
Anitescu
,
M.
,
2006
, “
Optimization-Based Simulation of Nonsmooth Rigid Multibody Dynamics
,”
Math. Program.
,
105
(
1
), pp.
113
143
.10.1007/s10107-005-0590-7
4.
Tasora
,
A.
, and
Anitescu
,
M.
,
2010
, “
A Convex Complementarity Approach for Simulating Large Granular Flows
,”
ASME J. Comput. Nonlinear Dyn.
,
5
(
3
), p.
031004
.10.1115/1.4001371
5.
Anitescu
,
M.
, and
Potra
,
F. A.
,
1997
, “
Formulating Dynamic Multi-Rigid-Body Contact Problems With Friction as Solvable Linear Complementarity Problems
,”
Nonlinear Dyn.
,
14
(
3
), pp.
231
247
.10.1023/A:1008292328909
6.
Stewart
,
D. E.
, and
Trinkle
,
J. C.
,
1996
, “
An Implicit Time-Stepping Scheme for Rigid Body Dynamics With Inelastic Collisions and Coulomb Friction
,”
Int. J. Numer. Methods Eng.
,
39
(
15
), pp.
2673
2691
.10.1002/(SICI)1097-0207(19960815)39:15<2673::AID-NME972>3.0.CO;2-I
7.
Fang
,
L.
,
2015
, “
A Primal-Dual Interior Point Method for Solving Multibody Dynamics Problems With Frictional Contact
8.
Heyn
,
T.
,
Anitescu
,
M.
,
Tasora
,
A.
, and
Negrut
,
D.
,
2013
, “
Using Krylov Subspace and Spectral Methods for Solving Complementarity Problems in Many-Body Contact Dynamics Simulation
,”
Int. J. Numer. Methods Eng.
,
95
(
7
), pp.
541
561
.10.1002/nme.4513
9.
Jan
,
K.
,
2015
,
Simulating Granular Material Using Nonsmooth Time-Stepping and a Matrix-Free Interior Point Method
,
Fraunhofer Verlag
, Stuttgart, Germany.
10.
Mazhar
,
H.
,
Heyn
,
T.
,
Negrut
,
D.
, and
Tasora
,
A.
,
2015
, “
Using Nesterov's Method to Accelerate Multibody Dynamics With Friction and Contact
,”
ACM Trans. Graph. (TOG)
,
34
(
3
), p.
32
.10.1145/2735627
11.
Petra
,
C.
,
Gavrea
,
B.
,
Anitescu
,
M.
, and
Potra
,
F.
,
2009
, “
A Computational Study of the Use of an Optimization-Based Method for Simulating Large Multibody Systems?
,”
Optim. Methods Software
,
24
(
6
), pp.
871
894
.10.1080/10556780902806094
12.
Daniel
,
M.
,
Luning
,
F.
,
Paramsothy
,
J.
, and
Dan
,
N.
,
2017
, “
A Comparison of Numerical Methods for Solving Multibody Dynamics Problems With Frictional Contact Modeled Via Differential Variational Inequalities
,”
Comput. Methods Appl. Mech. Eng.
,
320
, pp.
668
693
.10.1016/j.cma.2017.03.010
13.
Benzi
,
M.
,
2002
, “
Preconditioning Techniques for Large Linear Systems: A Survey
,”
J. Comput. Phys.
,
182
(
2
), pp.
418
477
.10.1006/jcph.2002.7176
14.
,
Y.
,
2003
,
Iterative Methods for Sparse Linear Systems
, Vol.
82
,
.
15.
Li
,
A.
,
Serban
,
R.
, and
Negrut
,
D.
,
2017
, “
Analysis of a Splitting Approach for the Parallel Solution of Linear Systems on GPU Cards
,”
SIAM J. Sci. Comput.
,
39
(
3
), pp.
C215
C237
.10.1137/15M1039523
16.
Corona
,
E.
,
Rahimian
,
A.
, and
Zorin
,
D.
,
2015
, “
A Tensor-Train Accelerated Solver for Integral Equations in Complex Geometries
,” preprint
arXiv: 1511.06029
.https://arxiv.org/abs/1511.06029
17.
Ambikasaran
,
S.
, and
Darve
,
E.
,
2013
, “
An $O$(N Log N) Fast Direct Solver for Partial Hierarchically Semi-Separable Matrices
,”
J. Sci. Comput.
,
57
(
3
), pp.
477
501
.10.1007/s10915-013-9714-z
18.
Chandrasekaran
,
S.
,
Dewilde
,
P.
,
Gu
,
M.
,
Lyons
,
W.
, and
Pals
,
T.
,
2007
, “
A Fast Solver for HSS Representations Via Sparse Matrices
,”
SIAM J. Matrix Anal. Appl.
,
29
(
1
), pp.
67
81
.10.1137/050639028
19.
Chandrasekaran
,
S.
,
Gu
,
M.
, and
Pals
,
T.
,
2006
, “
A Fast ULV Decomposition Solver for Hierarchically Semiseparable Representations
,”
SIAM J. Matrix Anal. Appl.
,
28
(
3
), pp.
603
622
.10.1137/S0895479803436652
20.
Gillman
,
A.
,
2011
, “
Fast Direct Solvers for Elliptic Partial Differential Equations
,” Ph.D. thesis, University of Colorado, Boulder, CO.
21.
Kenneth
,
L. H.
, and
Lexing
,
Y.
,
2015
, “
Hierarchical Interpolative Factorization for Elliptic Operators: Differential Equations
,”
Commun. Pure Appl. Math.
,
69
(
8
), pp.
1415
1451
.10.1002/cpa.21582
22.
Xia
,
J.
,
Chandrasekaran
,
S.
,
Gu
,
M.
, and
Li
,
X. S.
,
2010
, “
Superfast Multifrontal Method for Large Structured Linear Systems of Equations
,”
SIAM J. Matrix Anal. Appl.
,
31
(
3
), pp.
1382
1411
.10.1137/09074543X
23.
Xia
,
J.
,
Chandrasekaran
,
S.
,
Gu
,
M.
, and
Li
,
X. S.
,
2010
, “
Fast Algorithms for Hierarchically Semiseparable Matrices
,”
Numer. Linear Algebra Appl.
,
17
(
6
), pp.
953
976
.10.1002/nla.691
24.
Bebendorf
,
M.
,
2008
, “
A Means to Efficiently Solve Elliptic Boundary Value Problems
,”
Hierarchical Matrices
(Lecture Notes in Computational Science and Engineering, Vol.
63
),
Springer-Verlag
,
Berlin
.
25.
Börm
,
S.
,
Grasedyck
,
L.
, and
Hackbusch
,
W.
,
2003
,
Hierarchical Matrices
(Lecture Notes, Vol.
21
), Max-Planck-Institut für Mathematik, Leipzig, Germany.
26.
Börm
,
S.
,
2010
, Efficient Numerical Methods for Non-Local Operators: $H$2-Matrix Compression, Algorithms and Analysis (EMS Tracts in Mathematics, Vol.
14
),
European Mathematical Society (EMS)
,
Zürich, Switzerland
.
27.
Coulier
,
P.
,
Pouransari
,
H.
, and
Darve
,
E.
,
2017
, “
The Inverse Fast Multipole Method: Using a Fast Approximate Direct Solver as a Preconditioner for Dense Linear Systems
,”
SIAM J. Sci. Comput.
,
39
(
3
), pp.
A761
A796
.10.1137/15M1034477
28.
Pouransari
,
H.
,
Coulier
,
P.
, and
Darve
,
E.
,
2017
, “
Fast Hierarchical Solvers for Sparse Matrices Using Extended Sparsification and Low-Rank Approximation
,”
SIAM J. Sci. Comput.
,
39
(
3
), pp.
A797
A830
.10.1137/15M1046939
29.
Stewart
,
D. E.
,
2000
, “
Rigid-Body Dynamics With Friction and Impact
,”
SIAM Rev.
,
42
(
1
), pp.
3
39
.10.1137/S0036144599360110
30.
Barzilai
,
J.
, and
Borwein
,
J. M.
,
1988
, “
,”
IMA J. Numer. Anal.
,
8
(
1
), pp.
141
148
.10.1093/imanum/8.1.141
31.
Andersen
,
M.
,
Dahl
,
J.
,
Liu
,
Z.
, and
Vandenberghe
,
L.
,
2011
, “
Interior-Point Methods for Large-Scale Cone Programming
,”
Optim. Mach. Learn.
,
5583
, pp.
1
26
.http://people.compute.dtu.dk/~mskan/publications/mlbook.pdf
32.
Kučera
,
R.
,
Machalová
,
J.
,
Netuka
,
H.
, and
Ženčák
,
P.
,
2013
, “
An Interior-Point Algorithm for the Minimization Arising From 3D Contact Problems With Friction
,”
Optim. Methods Software
,
28
(
6
), pp.
1195
1217
.10.1080/10556788.2012.684352
33.
Kleinert
,
J.
,
Simeon
,
B.
, and
Obermayr
,
M.
,
2014
, “
An Inexact Interior Point Method for the Large-Scale Simulation of Granular Material
,”
Comput. Methods Appl. Mech. Eng.
,
278
, pp.
567
598
.10.1016/j.cma.2014.06.009
34.
Mangoni
,
D.
,
Tasora
,
A.
, and
Garziera
,
R.
,
2018
, “
A Primal–Dual Predictor–Corrector Interior Point Method for Non-Smooth Contact Dynamics
,”
Comput. Methods Appl. Mech. Eng.
,
330
, pp.
351
367
.10.1016/j.cma.2017.10.030
35.
Nocedal
,
J.
, and
Wright
,
S.
,
2006
,
Nonlinear Equations
,
Springer
, New York.
36.
Oseledets
,
I. V.
, and
Tyrtyshnikov
,
E.
,
2010
, “
TT-Cross Approximation for Multidimensional Arrays
,”
Linear Algebra Its Appl.
,
432
(
1
), pp.
70
88
.10.1016/j.laa.2009.07.024
37.
Kazeev
,
V. A.
, and
Khoromskij
,
B.
,
2012
, “
Low-Rank Explicit QTT Representation of the Laplace Operator and Its Inverse
,”
SIAM J. Matrix Anal. Appl.
,
33
(
3
), pp.
742
758
.10.1137/100820479
38.
Oseledets
,
I. V.
,
2010
, “
Approximation of 2d × 2d Matrices Using Tensor Decomposition
,”
SIAM J. Matrix Anal. Appl.
,
31
(
4
), pp.
2130
2145
.10.1137/090757861
39.
Oseledets
,
I. V.
,
Tyrtyshnikov
,
E.
, and
Zamarashkin
,
N.
,
2011
, “
Tensor-Train Ranks for Matrices and Their Inverses
,”
Comput. Methods Appl. Math.
,
11
(
3
), pp.
394
403
.10.2478/cmam-2011-0022
40.
Bern
,
M.
,
Eppstein
,
D.
, and
Teng
,
S.-H.
,
1999
, “
Parallel Construction of Quadtrees and Quality Triangulations
,”
Int. J. Comput. Geom. Appl.
,
9
(
06
), pp.
517
532
.10.1142/S0218195999000303
41.
Dolgov
,
S.
, and
Savostyanov
,
D. V.
,
2013
, “
Alternating Minimal Energy Methods for Linear Systems in Higher Dimensions—Part I: SPD Systems
,” Preprint
arXiv:1301.6068
.https://arxiv.org/abs/1301.6068
42.
Dolgov
,
S.
, and
Savostyanov
,
D. V.
,
2013
, “
Alternating Minimal Energy Methods for Linear Systems in Higher Dimensions—Part II: Faster Algorithm and Application to Nonsymmetric Systems
,” Preprint
arXiv: 1304.1222
.https://arxiv.org/abs/1304.1222
43.
Oseledets
,
I. V.
, and
Dolgov
,
S. V.
,
2012
, “
Solution of Linear Systems and Matrix Inversion in the TT-Format
,”
SIAM J. Sci. Comput.
,
34
(
5
), pp.
A2718
A2739
.10.1137/110833142
44.
Melanz
,
D.
,
Jayakumar
,
P.
, and
Negrut
,
D.
,
2016
, “
Experimental Validation of a Differential Variational Inequality-Based Approach for Handling Friction and Contact in Vehicle/Granular-Terrain Interaction
,”
J. Terramechanics
,
65
, pp.
1
13
.10.1016/j.jterra.2016.01.004
45.
Tasora
,
A.
,
Serban
,
R.
,
Mazhar
,
H.
,
Pazouki
,
A.
,
Melanz
,
D.
,
Fleischmann
,
J.
,
Taylor
,
M.
,
Sugiyama
,
H.
, and
Negrut
,
D.
,
2015
, “
Chrono: An Open Source Multi-Physics Dynamics Engine
,”
International Conference on High Performance Computing in Science and Engineering
(
HPCSE
),
Solan, Czech Republic
,
May 25–28
, pp.
19
49
.https://www.researchgate.net/publication/303761434_Chrono_An_Open_Source_Multi-physics_Dynamics_Engine
46.
Oseledets
,
I. V.
,
2012
, “TT-Toolbox 2.2,” Skoltech, Moscow, Russia.