• An Algebraic Multigrid Algorithm based on Smoothed ...


  •   
  • FileName: slides_beamer3.pdf [read-online]
    • Abstract: Newton’s method for solving the LCP. Non-Smooth Newton’s Method. Algebraic Multigrid Preconditioner based on Smooth Aggregation ... Newton’s method for solving the LCP. Non-Smooth Newton’s Method. Algebraic Multigrid ...

Download the ebook

Geometric Multigrid, Algebraic Multigrid
Rigid Multibody Dynamics with Frictional Contacts
Newton’s method for solving the LCP
An Algebraic Multigrid Algorithm based on
Smoothed Aggregation for Preconditioning a Large
Sparse System of Linear Equations
R. Ortiz
Department of Mathematics
The University of Iowa
Numerical Analysis Seminar - Oct. 17 2006
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Rigid Multibody Dynamics with Frictional Contacts
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
Geometric Multigrid (GMG or just MG)
Originally used for solving the matrix equations that arise in a wide
variety of applications, including discretized PDEs.
It automatically construct a sequence of increasingly smaller matrix
problems that hopefully enables efficient resolution of all scales present
in the solution
The methodology is based on measuring how the algebraically smooth
error value at one point depends in its value at another.
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
Boundary Value Problem Example
Consider the following boundary value problem
−u (x) = f (x), 0 < x < 1,
u(0) = u(1) = 0
Create Grid: The interval [0, 1] is discretized using centered difference
approximations, using equally spaced n + 2 points.
xi = i × h, i = 0, . . . , n + 1,
where h = 1/(n + 1).
Replacing u (xi ) by [u(xi−1 ) − 2u(xi ) + u(xi+1 )]/h2 , results in the
system Ax = b, where
   
2 −1 f (x0 )
 −1 2 −1   . 
   
 . . , b = h2 
  . 
A= 

 . . . 


 . 

 . . −1   . 
−1 2 f (xn−1 )
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
MG Motivation
MG can be motivated by taking an in-depth look at a simple iterative
scheme such as Richardson iteration.
uk+1 = (I − ωA)uk + ωb
Convergence takes place if 0 < ω < 2/ρ(A). Let Mω = I − ωA.
If u∗ is the exact solution, then ek = (Mω )k e0 .
If we expand the vector e0 = ∑n ξk wk
k=1
n
ek = ∑ (1 − λk /γ)j ξk wk
k=1
Ricardo Ortiz AMG SA
Note that each component is reduced by (1 − λk /γ)j , where ρ(A) ≤ γ.
The smallest converging component correspond to the smallest
eigenvalue λ1 .
θk
For the model problem above we have γ = 4, λk = 4 sin2 ( ) , and
2

wk = [sin θk , sin(2θk ), . . . , sin(nθk )]T , for θk =
n+1
So
λ1 π πh 2
1− = 1 − sin2 ≈ 1−( ) = 1 − O(h2 )
γ 2(n + 1) 2
The basic observation on which MG methods are founded is that
convergence is not the same for all components.
Half of the error components see a very good decrease. This is the case
of the high frequency components k > n/2.
kπ kπ 1
ηk = 1 − sin2 = cos2 ≤
2(n + 1) 2(n + 1) 2
Two important observations:
The oscillatory components undergo better than 1/2 reduction at each
step of the iteration.
This factor is independent of h
So we can only reduce efficiently half of the components, so how do we
reduce the other components?
Here is where we introduce a coarse-grid problem.
11
0 0 Fine Grid 1
0
1 11
0 00
Coarse Grid
11
00
1
0 11
00 5 6 7 0 0
0 1 2 3 4
10
8 1 1 11
00
2 3 11
00
4
111
000
1
0 11
00
Figure: k = 2 mode on a fine grid (n = 7) and a coarse grid (n = 3)
Some of the modes that where smooth on the finer grid will become
oscillatory in the coarse grid.
And at the same time the oscillatory modes present in the fine mesh are
not longer represented in the coarse grid.
The iteration fails to make progress on the fine grid when the only
components left are those associated with smooth modes.
MG strategies does not attempt to eliminate these components on the
fine grid, instead they first move down a to a coarser grid where smooth
modes are translated in to oscillatory ones.
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
AMG
GMG depend on the availability of an underlying mesh, in addition its
performance deteriorates for problems with anisotropic or discontinuous
coefficients
It is difficult also to define GMG on domains that are not rectangular
Initially introduced in 1984 by Brant A, McCormick SF, Ruge JW
AMG offer the efficiency of GMG methods but without reliance on
knowledge of the grid geometry or of variation in the coefficient of the
differential operator
For this reasons, AMG algorithms are often the method of choice for
solving linear systems arising from discretizations over unstructured
meshes
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Multigrid Methods
Rigid Multibody Dynamics with Frictional Contacts
Algebraic Multigrid (or AMG)
Newton’s method for solving the LCP
Algorithm
vk+1 ← MGk+1 (bf , vk , µ1 , µ2 , γ)
1: Let Rl : Rnl → Rnl , l = 1, · · · , L − 1 be the given smothers and µi , γ > 0 be
a given smoothing and cycle parameters, respectivelly
2: Set MG = MG1 where MGl (·, ·) for l = 1, · · · , L − 1 is defined by
3: Pre-smoothing: Perform µ1 iterations of xl ← (I − Rl AL )xl + Rl bl
4: Coarse grid correction:
5: Set bl+1 = (Pl )(bl − Al xl )
l+1
6: if l + 1 = L then
7: solve AL xL = bL by a direct method
8: else
9: Set xl+1 = 0 and perform γ iterations of xl+1 ← MGl+1 (xl+1 , bl+1 )
10: end if
11: Correct the solution on level l by xl ← xl + Pl xl+1l+1
12: Post-smoothing: Perform µ2 iterations of xl ← (I − Rl AL )xl + Rl bl
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Frictionless Case
Define the feasible set as C = {q | ϕ(q) ≥ 0}, with ϕ : ℜn → ℜ.
For a given particle at a point q ∈ C
1 Kinetic energy: 1 vT M(q)v, where M(q) is the generalized mass matrix
2
and v = q, the generalized velocity vector.
˙
2 Potential energy: V(q)
Lagrangian
1
L (q, v, λ) = 2 vT M(q)v − V(q) + N ϕ(q)
Where N is the a Lagrange multiplier.
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Frictionless Case
Define the feasible set as C = {q | ϕ(q) ≥ 0}, with ϕ : ℜn → ℜ.
For a given particle at a point q ∈ C
1 Kinetic energy: 1 vT M(q)v, where M(q) is the generalized mass matrix
2
and v = q, the generalized velocity vector.
˙
2 Potential energy: V(q)
Lagrangian
1
L (q, v, λ) = 2 vT M(q)v − V(q) + N ϕ(q)
Where N is the a Lagrange multiplier.
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Frictionless Case
Euler-Lagrange Equation of Motion with q ∈ ∂C
equation dv
M(q) = k(q, v) − ∇V(q) + N∇ϕ(q)
d ∂L ∂L dt
− =0
dt ∂v ∂q
0 ≤ N ⊥ ϕ(q) ≥ 0
1 ∂Mir ∂Mis ∂Mrs
Here ki (q, v) = − + − vr vs
2 ∑ ∂qs
r,s ∂qr ∂qi
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Frictional Case (Coulomb Friction)
N∇ϕ(q)
The friction force inside this polyhedral
friction cone can be represented by
D(q)β, β ≥ 0, eT β ≤ µN
Now we use the maximum dissipation
principle(Erdmann, 1994)
β minimizes vT D(q)β
subject to β≥0
eT β ≤ µN
Figure: Polyhedral friction cone
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Fictional Case
Using the Kuhn-Tucker conditions we reformulate this LP as
Complementarity Equation of Motion, here again assume
conditions q ∈ ∂C
T dv
0 ≤ β ⊥ D(q) v + λe ≥ 0
M(q) = k(q, v) − ∇V(q) + N∇ϕ(q) + D(q)β
dt
0 ≤ λ ⊥ µN + eT β ≥ 0
0 ≤ N ⊥ ϕ(q) ≥ 0
0 ≤ β ⊥ D(q)T v + λe ≥ 0
0 ≤ λ ⊥ µN + eT β ≥ 0
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Implicit Euler Method
Now we set up a time-stepping method which allows for impulsive forces
Let h = t +1 − t and q ≈ q(t ), v ≈ v(t )
Our final frictional problem becomes: Given q , v ; find
q +1 , v +1 , β ,N ,λ
+1
M(q)(v − v ) = h k(q , v ) − ∇V(q ) + h∇ϕ(q )T N
+hD(q )β
+1
0≤N ⊥ ∇ϕ(q )T γv + (1 − γ)v ≥0
+1
0≤β ⊥ D(q )T γv + (1 − γ)v + Eλ ≥ 0
0≤λ ⊥ µN − ET β ≥ 0
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Equations of Motion
Rigid Multibody Dynamics with Frictional Contacts
Time-Stepping
Newton’s method for solving the LCP
Implicit Euler Method
We can rearrange and substitute so it looks like a big LCP of the in the
form
0 ≤ z ⊥ Az + b ≥ 0
The LCP matrix
∇ϕ(q )T M(q )−1 ∇ϕ(q ) ∇ϕ(q )T M(q )−1 D(q ) 0
 
A =  D(q )T M(q )−1 ∇ϕ(q ) D(q )T M(q )−1 D(q ) E 
µ −ET 0
Note here A is not symmetric
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
LCP Solvers
1 There are three types of LCP solvers
2 Pivoting methods:
Lemke’s pivot method and its variants(Lemke 1964)
Cottle-Dantzig’s principal pivot method (Cottle 1968)
Keller’s principal pivot method (Keller 1973)
among others...
3 Newton Methods: Line search methods
4 Iterative Methods: Splitting methods
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Newton Method
We now implement a non-smooth Newton scheme to solve the LCP
0 ≤ z ⊥ Az − b ≥ 0
Define H(z) = min(Az − b, z) for the frictionless case and understood
component-wise
A vector z is a solution of the LCP if and only if H(z) = 0
So we translate the LCP to a zero finding problem that we can solve with
a Newton method
1 A function H : ℜn → ℜn is said to be B-differentiable at a point z ∈ ℜn if
H is Lipschitz continuous on a neighborhood of z and
Directionally differentiable at z along a direction d ∈ ℜn , denoted as
∇H(z; d).
2 We extend the classical Newton method by replacing the linear Newton
equation by one obtained from the use of the B-derivative.
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Algorithm(J.S. Pang 1990)
B-differentiable Newton method for z ← NSN(A, z0 , b, α, β, ε)
solving H(z) = 0
1: Initialize β, α, ε > 0
Θ(z) = 1 H(z)T H(z) is the merit
2 2: Solve for dz in
function
H(zl ) + H (zl )dzl = 0
∇H(z) is the Jacobian matrix defined 3: Determine step-size
as 4: if Θ(zl + βm dzl ) ≤
If (Az − b)i < zi , then (1 − 2αβm )Θ(zl ) then
5: Set zl+1 ≡ zl + βm dzl
∇H(i, :) = A(i, :)
6: end if
If (Az − b)i ≥ zi , then 7: if Θ(zl+1 ) ≤ ε then
8: Terminate
∇H(i, :) = Id(i, :)
9: else
10: Return to Step 2
11: end if
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Outline
1 Multilevel Methods
Multigrid Methods
Algebraic Multigrid (or AMG)
2 Rigid Multibody Dynamics Model
Equations of Motion
Time-Stepping Scheme
3 Newton’s method for solving the LCP
Non-Smooth Newton’s Method
Algebraic Multigrid Preconditioner based on Smooth Aggregation
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Algebraic Multigrid
Algebraic multigrid do not require geometric grid information so they are
appropriate for Multibody Dynamics
We use an AMG to precondition the system
∇H(z )dz = −H(z )
The system matrix takes a special form
A A
∇H(z ) =
0 I
Where A = PT AP and symmetric positive semi-definite
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Coarsening Strategy for the Contacts
The strategy consist in constructing the multigrid operators based on a
clustering of nodes in the contact graph
For any multigrid algorithm, the same fundamental components are
required
There must be a sequence of grids
Intergrid transfer operators (Restriction and Interpolation)
A relaxation or smoothing operator (Jacobi, GS, SOR, etc.)
Coarse-grid version of the fine-grid operator
A solver for solving a small system in the coarsest grid
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Clustering Strategy
Here is a 8 bodies configuration making 8 contacts with each other
Black dots represent constraints between two bodies
7
7
C4
8 8
3 3
2 2
C1
6
6
C3
1
1
4 C2 5
4 5
(b) Clustering of nodes, coarse
(a) Node configuration in finer level
level
Figure: Clustering Strategy
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Smoothed Aggregation (SA)
During the last decade, SA has emerge as an efficient AMG solver, in
particular from problems that arise from 3D thin body elasticity
As with regular AMG, SA construct intergrid operators by doing certain
assumptions about the nature of the algebraically smooth error
Usually for SA, knowledge of the near-kernel is needed. An improved
SA, Adaptive SA deal with the case where we don’t have that
information available (Brezina et al. 2005).
The hierarchy of coarse level matrices is defined by
Al+1 = (Sl Pl )T Al Sl Pl , A1 = A
l+1 l+1
P1 = P1 P2 . . . Pl−1 is known as the composite tentative prolongator
l+1 2 3 l
and Sl is the prolongator smoother defined by
4 −1
Sl = I − M Al
3λl l
where Ml = (P1 )T P1 and λl ≥ ρ(Ml−1 Al )
l l
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Clustering Strategy
Suppose we are dealing with the frictionless case for the moment, then
we can write the system matrix as
A = JM −1 J T
Each row of J correspond to a linear constraint between any two bodies
in the configuration
Use aggregation to construct the coarse-grid operators
We define the interpolation operator as
1 if i ∈ Cjl
(Pl
l+1 )ij =
0 otherwise
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Multigrid Algorithm - V-Cycle Scheme
vk+1 ← AMG(bf , vk , µ1 , µ2 )
1: Define l = 0 and set A0 = A, b0 = b and x0 = vk
2: while l < levels − 1 do
3: Relax µ1 times on Al xl = bl
4: Calculate residual rl = Al xl − bl
5: Set bl+1 = (Sl Pl )T rl
l+1
6: Set l = l + 1
7: end while
8: Solve the system Alevels xlevels = blevels /// Solve coarsest level
9: while l ≥ 0 do
10: Set xl = xl + Sl Pl xl+1
l+1
11: Relax µ2 times on Al xl = bl
12: end while
13: Set vk = x0
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Implementation
We have implemented this algorithms in Matlab and C++.
This algorithm is for the frictionless case, although friction will be
considered in near future.
Ricardo Ortiz AMG SA
Geometric Multigrid, Algebraic Multigrid
Non-Smooth Newton Method
Rigid Multibody Dynamics with Frictional Contacts
AMG iterative solver
Newton’s method for solving the LCP
Further Reading I
Y. Saad.
Iterative Methods for Linear Systems.
SIAM, 2003.
WL Brigs, VE Henson, SF McCormick.
A multigrid tutorial.
SIAM, 2000.
S. Oliveira, D. Stewart.
Writting scientific software.
Cambridge, 2006.
Petr Vanek , Marian Brezina , Jan Mandel.
Convergence of algebraic multigrid based on smoothed aggregation.
Numer. Math. (2001) 88: 559-579.
Ricardo Ortiz AMG SA


Use: 0.2986