Requirements for iterative linear equation solvers in problems of - PowerPoint PPT Presentation

Requirements for iterative linear equation solvers in problems of solid- and structural mechanics Reijo Kouhia Tampere University http://www.tut.fi/rakmek/personnel/kouhia/ MTU60: Let’s make it sparse, Prague, Czech Republic, Jan 31 - Feb 2, 2019

Happy birthday Mirek! Ver¨ aj¨ am¨ aki May 27, 2000.

Some history ◮ I became interested in iterative linear solvers around 1995. ◮ After hitting my head to the wall with X preconditioners I found the paper: M. Benzi, C.D. Meyer, M. T˚ uma, A sparse approximate inverse preconditioner for the CG method, SIAM J. Sci. Comp. , 17 ,1135-1149 (1996) and I felt fascinated about the idea. ◮ I sent a mail to Michele (autumn 1996?) and then everything started to evolve fast. ◮ Mirek, Michele and I wrote our first journal paper and it was published 1998 in Commun. Numer. Methods Eng. ◮ I met Mirek first time face to face in Prague May 1999 and Michele in Finland July 2004. ◮ Mirek visited Helsinki in May 2000.

It started with shells Bending dominated loading, 30 × 30 uniform mesh, stabilized quadrilateral MITC4 shell elements 5489 unknowns. Convergence of the PCG depends highly on the “thinness” of the shell, here t/R = 0 . 1 , 0 . 01 , 0 . 001 and convergence with the IC(0)-preconditioner shown below. M. Benzi, RK, M. T˚ uma, An assessment of some preconditioning techniques in shell problems, Communications in Numerical Methods in Engineering , 14 , 1998, 897-906

In 1999 it continued with solids etc. Block versions of the stabilized AINV Discussion about artificial rhs-vectors Ver¨ aj¨ am¨ aki May 27, 2000. M. Benzi, RK, M. T˚ uma, Stabilized and block approximate inverse preconditioners for problems in solid and structural mechanics, Computer Methods in Applied Mechanics and Engineering , 190 , 2001, 6533-6554

Why I am interested in shells? - Shells are prone to stability problems R/L = 4 (a) R/t = 100 (b) R/t = 400 P L R θ θ

Path-following methods The non-linear mapping f defines the equilibrium path in the displacement q and load parameter λ space: f ( q , λ ) ≡ r ( q ) − λ p r ( q ) = 0 (1) and constitutes the balance between internal- and external forces. To overcome limit points an additional equation is usually augmented to (1), thus we have to solve � f ( q , λ ) = r ( q ) − λ p r ( q ) = 0 (2) c ( q , λ ) = 0 . Newton’s linearisation step results in A δ x = − g , (3) where � K � � δ q � � � , K = f ′ = ∂ f c = ∂c ∂ q , e = ∂c p r f A = , δ x = , g = ∂ q , ∂λ. c T e δλ c

Requirements for linear solvers in non-linear analysis ◮ Be able to handle indefinite matrices. ◮ Symmetry (if existing) of the Jacobian matrix of f should be exploitable. ◮ Information about definiteness (positive definite/indefinite) of the matrix should be extractable.

Stability of equilibrium state ? Singularity test functions (STF) ◮ determinant ◮ smallest eigenvalue ◮ smallest pivot ◮ current stiffness parameter Requirements for the STF ◮ mesh size independent ◮ “good predictor” ◮ indicates all kind of singular behaviour ◮ easy and cheap to evaluate

Check for the positivity of a matrix A when applying preconditioned Krylov subspace methods to solve Ax = b with accuracy � r i � = � b − Ax i � < ε r � b − Ax 0 � + ε a . ◮ λ i is an eigenvalue of A 1. form M (or directly M − 1 ), ◮ λ 1 = min λ i 2. initialize r 0 = b − Ax 0 , apply ◮ | λ 1 | < δ preconditioner d 0 = M − 1 r 0 , compute τ 0 = r T 0 d 0 , Conjecture 3. iterate i = 0 , 1 , 2 , ... until convergence: If A is indefinite a i = τ i / d T 3.1 compute: s = Ad i , i s , then there exists with probability 1 − η ( δ, ε r , ε a ) 3.2 update: x i +1 = x i + a i d i , at least one iterate k such that d T k Ad k < 0 . r i +1 = r i − a i s , 3.3 apply preconditioner: z = M − 1 r i +1 , 3.4 compute: τ i +1 = r T i +1 z , b i = τ i +1 /τ i , To be usefull: η ≪ 1 and Ad k should be 3.5 update: d i +1 = z + b i d i . available

It works with PCG but not with preconditioned Bi-CGSTAB - Why? 1. form M (or directly M − 1 ), r T r 0 , set d 0 = r 0 , 2. compute r 0 = b − Ax 0 , choose ˜ r , compute ρ 0 = ˜ 3. iterate i = 0 , 1 , 2 , ... , until converged: 3.1 apply preconditioner: z = M − 1 d i , r T v i , 3.2 compute: v i = Az , a i = ρ i / ˜ s = r i − a i v i , s = M − 1 s , 3.3 apply preconditioner: ˜ ω i = w T s / w T w , 3.4 compute: w = A ˜ s , 3.5 update: x i +1 = x i + a i z + ω i ˜ s , r i +1 = s − ω i w , 3.6 compute: � r i +1 � 2 , r T r i +1 and 3.7 compute: ρ i +1 = ˜ b i +1 = ( ρ i +1 /ρ i )( a i /ω i ) , 3.8 update: d i +1 = r i +1 + b i +1 ( d i − ω i v i ) .