Preconditioning of Elliptic Saddle Point Systems by Substructuring - - PowerPoint PPT Presentation

Preconditioning of Elliptic Saddle Point Systems by Substructuring and a Penalty Approach

16 16th

th International Conference on Domain

International Conference on Domain Decomposition Methods Decomposition Methods January 12-15, 2005

Clark R. Dohrmann Structural Dynamics Research Department Sandia National Laboratories Albuquerque, New Mexico

Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Department of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

• DD16 scientific & organizing committees
• Sandia

– Rich Lehoucq – Pavel Bochev – Kendall Pierson – Garth Reese

• University

– Jan Mandel

Thanks

Overview

• Saddle Point Preconditioner

– related penalized problem – positive real eigenvalues – conjugate gradients

• BDDC Preconditioner

– overview – modified constraints

• Examples

– H(grad), H(div), H(curl)

Motivation

• riginal system:

penalized (regularized) system: Axelsson (1979) ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − g f p u C B B A

T

~ C C C

T

~ ~ , ~ = > ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − g f p u C B B A

T

A > 0 on kernel of B, C ≥ 0 AT = A, CT = C, B full rank

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −

p u p u T

r r z z C B B A ~

) ~ ( ) ( ~

1 1 1 p T u A u p u p

r C B r S z r Bz C z

− − −

+ = − =

B C B A S

T A 1

~− + = exact solution of penalized system primal rather than dual Schur complement considered is ) ~ (

1 ~ T C

B BA C S

−

+ − =

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + ⇒ f p u B B B B A

T T

ρ f u B B A

T

= + ⇒ ) ( ρ

Some Connections

2 / ) ( ) ( 2 / Bu Bu f u Au u G

T T T

ρ + − = penalty method for C = 0 and g = 0: ρ / ~ I C = matrix same as SA for Bu p G L

T

+ = augmented Lagrangian:

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ −

p u p u T

r r z z C B B A ~ ) ( ~ 1

p u p

r Bz C z − =

−

regularized constraint preconditioning: Benzi survey (2005) regularized constraint equation satisfied exactly

Penalty Preconditioner

) ( ~ ) ~ (

1 1 1 p u p p T u A u

r Bz C z r C B r S z − = + =

− − −

B C B A S

T A 1

~− + = recall exact solution of penalized system ⇒ preconditioner: ) ( ~ ) ~ (

1 1 1 p u p p T u u

r Bz C z r C B r M z − = + =

− − −

M preconditioner for SA

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛

− − − − − p u T p u

r r I C B I C M I B C I z z

• 1

1 1 1 1

~ ~ ~ M ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − = C B B A

T

A matrix representation: A M 1

−

question: what about spectrum of

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = C C M S A ~ H z z M A λ =

z z H A HM λ =

−

• symmetric

1

eigenproblem: introduce: Bramble & Pasciak (1988), Klawonn (1998) eigenvalues same as those for can show H > 0 ⇒ HM-1A > 0

equivalent linear system: HM-1Aw = HM-1b

CG Connection (B&P too)

• riginal linear system: Aw = b

⇒ solve using pcg with H as preconditioner eigenvalues of preconditioned system same as those of M-1A H not available, but not a problem ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = C C M S A ~ H

k k k k k k k k k k k k k k k k

p r r p z z p r r p w w A HM A M A A

1 1 1 1 1 1

~ ~

− − − − − −

− = − = − = + = α α α α ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =

− − − −

) ( ~ ) ~ (

1 1 1 1 p u p T u p u

a Bd C a C B a M d d a M ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − + − =

− − p p u p T u u A

Cd a Bd a C B a d S a ) ~ (

1 1

HM required calculations: recurrences:

Theory

2 1 1

) ( δ σ δ ≤ ≤

− A

M

m T T T T T n T A T T

p p B BM p p C p p B BM p u Mu u u S u Mu u ℜ ∈ ∀ ≤ ≤ ℜ ∈ ∀ ≤ ≤

− − 1 2 1 1 2 1

~ γ γ α α

Theorem (C = 0): Given α1 > 1, γ1 > 0, and and σ1 > 0, σ2 > 0, σ1 + σ2 = 1 eigenvalues satisfy where

{ } { }

1 2 1 2 2 2 2 2 1 2 1 2 1

/ ) / 2 ( , 2 max ) /( ), / ( min γ α σ σ α δ γ α σ α α σ δ − − = =

as ε → 0

( )

2 / 3 ) ( 2 / /

2 1 2 1

α σ α α ≤ ≤

− A

M

Simplification for

C B BS

T A

~

1

→

− m T T T T T n T A T T

p p B BM p p C p p B BM p u Mu u u S u Mu u ℜ ∈ ∀ ≤ ≤ ℜ ∈ ∀ ≤ ≤

− − 1 2 1 1 2 1

~ γ γ α α recall with α1 > 1 ⇒ γ1 bounded below by 1/α2 and γ2 bounded above by 1/α1 ⇒

~ C C ε =

two goals:

• 1. Ensure α1 > 1 (i.e. SA – M > 0)
• 2. Minimize α2/α1 (M good preconditioner for SA)

potential issues:

• 1. Scaling preconditioner M to satisfy Goal 1
• 2. SA becomes very poorly conditioned as ε → 0

conclusion: effective preconditioner for SA is essential

Recap of Penalty Preconditioner

• Based on approximate solution of related

penalized problem

• Conjugate gradients can be used to solve saddle

point system

• Theory provides conditions for scalability
• Effectiveness hinges on preconditioner for primal

Schur complement SA

Ref: C. R. Dohrmann and R. B. Lehoucq, “A primal based penalty preconditioner for elliptic saddle point systems,” Sandia National Laboratories, Technical report SAND 2004-5964J.

BDDC Overview

• Primal Substructuring Preconditioner

– no coarse triangulation needed – recent theory by Mandel et al – multilevel extensions straightforward

• Numerical Properties

– C(1+log(H/h))2 condition number bounds – preconditioned eigenvalues all ≥ 1 (woo hoo!)

• eigenvalues identical to FETI-DP

– local and global problems only require sparse solver for definite systems

M7

FE discretization

• f pde

elements partitioned into substructures solve for unknowns

• n sub boundaries

(Schur complement)

Building Blocks

Additive Schwarz:

– Coarse grid correction v1 – Substructure correction v2 – Static condensation correction v3

FETI-DP counterparts:

k T I cc I

rc rc

F K F v λ

1

* 1

−

↔

∑

=

−

↔

s T

N s k s r s rr s r

B K B v

1 2

1

λ ↔

3

v

Dirichlet preconditioner

3 2 1 1

v v v r M + + =

−

Coarse Grid and Substructure Problems

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ Λ Φ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ I Q Q K

i i i T i i

i i T i ci

K K Φ Φ =

⇒

• Kci coarse element matrix
• assemble Kci ⇒ Kc
• Kc positive definite

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛

2 2 i i i i T i i

r v Q Q K λ

substructure problem: coarse problem:

Mixed Formulation of Linear Elasticity

f u B C B A f p u C B B A

T T

= + ⇒ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ −

−

) (

1

• A + BTC-1B has same sparsity as A (discontinuous pressure p)
• BDD preconditioner with enriched coarse space investigated

by Goldfeld (2003)

• related work for incompressible problems in Pavarino and

Widlund (2002) and Li (2001), see also M7

B M B A B C B A

p T T 1 1 − −

+ = + ⇒ λ µ Let A = µA0 and C = (1/λ)Mp where ) 2 1 )( 1 ( ) 1 ( 2 ν ν ν λ ν µ − + = + = E E recall f u B C B A

T

= +

−

) (

1

notice as ν → ½ that λ→ ∞ ⇒ condition number → ∞

2D Plane Strain Example

6.2e3 0.4999999 6.2e2 0.499999 63 0.49999 7.2 0.4999 2.2 0.499 2.1 0.49 1.8 0.4 2.0 0.3

κ ν condition number estimates for 4 substructure problem with H/h = 8

κ ≈ 1/(1 − 2ν) for ν near ½ Q: why does κ → ∞ as ν → ½ for BDDC preconditioner?

• f preconditioned system

• Explanation

– Movement of substucture boundary nodes is weighted average of coarse grid and substructure corrections (v1 and v2) from neighboring subs – If substructure volume changes, then strain energy of substructure → ∞ as ν → ½ – cg step length α → 0 since cg minimizes energy

• Solution

– Modify “standard” BDDC constraint equations to enable enforcement of zero volume change – Same number of constraints in 2D slightly more in 3D

condition number estimates for 4 substructure problem with H/h = 8

modified

• riginal

2.3 2.3 2.3 2.3 2.3 2.2 1.8 2.0 6.2e3 0.4999999 6.2e2 0.499999 63 0.49999 7.2 0.4999 2.2 0.499 2.1 0.49 1.8 0.4 2.0 0.3

κ ν

Recap of BDDC Preconditioner

Ref: C. R. Dohrmann, “A substructuring preconditioner for nearly incompressible elasticity problems,” Sandia National Laboratories, Technical report SAND 2004-5393J.

• Performance sensitive to values of ν near ½ if

“standard” constraints used

• Sensitivity caused by substructure volume

changes for nearly incompressible materials

• Simple modification of constraints can effectively

accommodate problems with ν near ½

Examples

• Problem Types

– H(grad): incompressible elasticity – H(div): Darcy’s problem – H(curl): magnetostatics

• Preconditioning Approaches

– Penalty/BDDC for H(grad) & H(div) – BDDC for stabilized H(curl) problem

2D Structured Meshes

H(grad): Q2 - P1 H(grad): P2+ - P1 H(div) & H(curl): RT0 & Ned

2D Unstructured Mesh

2937 elements

3D Unstructured Mesh

10194 elements

Incompressible Elasticity

mixed variational formulation:

q udx q w wdx f wdx p dx v u ∀ ∫ = ⋅ ∇ ∀ ∫ ∫ ∫ ⋅ = ⋅ ∇ +

Ω Ω Ω Ω

) ( : ) ( 2 ε ε µ

where εi,j(u) = (ui,j + uj,i)/2, u,w ∈ H1(Ω), p,q ∈ L2(Ω), µ = 1 recall:

⎟ ⎠ ⎞ ⎜ ⎝ ⎛ = B B A

T

A

penalty preconditioner: = (2,2) block of A for ν < ½ C ~ −

2D Incompressible Plane Strain Example

9 (2.6) 9 0.49999 9 (2.7) 9 0.4999 10 (2.7) 10 0.499 11 (3.0) 11 0.49 17 (7.2) 15 0.4 23 (16) 19 0.3

PCG GMRES ν iterations for rtol = 10-6 for 16 substructure problem with H/h = 8

Note: ν is Poisson ratio used to define in penalty preconditioner C ~

Other Saddle Point Preconditioners

block diagonal: Fortin, Silvester, Wathen, Klawonn, …

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =

p A D

M M M ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − =

p T A T

M B M M

MA: BDDC preconditioner for A, Mp: pressure mass matrix block triangular: Elman, Silvester, Klawonn, … ⇒ MINRES ⇒ GMRES

2D Structured Mesh Comparison

30 47 11 (3.1) 10 256 30 45 11 (3.1) 10 196 29 42 10 (3.1) 10 144 28 40 10 (3.0) 10 100 26 38 10 (2.9) 9 64 23 35 9 (2.6) 9 36 20 30 8 (2.1) 8 16 16 26 6 (1.8) 6 4 GMRES GMRES PCG GMRES

MT MD M N iterations for rtol = 10-6 for N substructure problem with H/h = 4 Note: ν = 0.49999 for penalty preconditioner

Unstructured Mesh Comparisons

95 140 27 (37) 23 3D 34 53 11 (3.2) 11 2D GMRES GMRES PCG GMRES

MT MD M dim Note: ν = 0.49999 in 2D and ν = 0.4999 in 3D for penalty preconditioner

Darcy’s Problem

governing equations: compatibility condition: penalty preconditioner: chosen as εDs and BDDC for SA C ~ u = −K∇ p in Ω ∇⋅ u = f in Ω u ⋅ n = 0 on ∂Ω ∫Ω fdx = 0 discretization: lowest-order R-T simplicial elements

2D Darcy’s Problem Example

5 (1.5) 5 0.00001 7 (1.9) 6 0.0001 10 (5.8) 10 0.001 20 (34) 18 0.01 47 (2.1e2) 38 0.1 100 (1.9e3) 91 1

PCG GMRES ε iterations for rtol = 10-6 for 16 substructure problem with H/h = 8, K = I

Note: in penalty preconditioner

s

D C ε = ~

⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =

p div D

M M

1

M ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ =

Sp A D

M M

2

M

Other Saddle Point Preconditioners

block diagonal 1: norm equivalence (Klawonn, 1995) block diagonal 2: some others: balancing Neumann-Neumann (BDD),

• verlapping methods (SPD reduction for div free space)

∫ ∫ ⋅ ∇ ⋅ ∇ + ⋅ =

Ω Ω

dx w u wdx u w u

div

) )( ( ) , (

T

B BA Sp

1 −

=

2D Structured Mesh Comparison

10 7 (2.0) 6 256 10 7 (1.9) 6 196 10 7 (1.9) 6 144 10 6 (1.8) 6 100 10 6 (1.7) 6 64 10 6 (1.6) 5 36 8 4 (1.2) 4 16 5 2 (1.01) 2 4 GMRES PCG GMRES

MD1 M N iterations for rtol = 10-6 for N substructure problem with H/h = 4, K = I

5

10

−

= ε

similar results for 3D

2D Structured Mesh (H/h dependence)

6 (2.09) 20 5 (1.92) 16 5 (1.74) 12 5 (1.53) 8 4 (1.24) 4

PCG H/h iterations for rtol = 10-6 for 9 substructure problem with K = I similar results for 3D and K ≠ constant

Unstructured Mesh Comparisons

17 8 (2.5) 8 3D 12 5 (1.3) 5 2D GMRES PCG GMRES

MD1 M dimen Note: ε = 10-5 for penalty preconditioner

Magnetostatics (B = ∇ × u)

mixed variational form with Coulomb gauge (∇ ⋅ u = 0): ∇ ⋅ J = 0 ⇒ ∫J ⋅ ∇p = 0 and w = ∇p ⇒ p = 0 ⇒

) ( ) ( ) )( / 1 (

0 curl

H w wdx J pdx w dx w u ∈ ∀ ∫ ∫ ⋅ = ∇ ⋅ + ∫ × ∇ ⋅ × ∇

Ω Ω Ω

µ ∫ ∈ ∀ = ∇ ⋅

Ω 1

H q qdx u

Note: ∇ ⋅ u was integrated by parts. Bummer, but … Acurlx = b with Acurl ≥ 0, but system consistent and B unique

Q: How to solve Acurlx = b with Acurl ≥ 0? One option: Solve in space restricted to ∇ ⋅ u = 0 ⇒ need basis for this space or back to saddle point system Another option: Solve Acurl x = b using CG with preconditioner for Ap Ap = Acurl + Adiv Advantage: can apply preconditioners for H(grad) problems!

Ref: C. R. Dohrmann, “Preconditioning of curl-curl equations by a penalty approach,” Sandia National Laboratories, Technical report, in preparation.

Example

128 elements, 208 edges, 81 nodes Acurl Ap

80 zevals before 49 zevals after 2 zevals before 0 zevals after

where σ > 0 ⇒ (Acurl + σAmass)x = b, σ → 0 equivalance, and precondition:

Some Other Options

• 1. precondition saddle point system: plenary L13 (Zou)
• 2. introduce regularized problem: Rietzinger & Schöberl (2002)

∫ ∫ ⋅ = ⋅ + ∫ × ∇ ⋅ × ∇

Ω Ω Ω

wdx J wdx u dx w u σ µ ) ( ) )( / 1 (

• multigrid (Hiptmair, Arnold, Falk, Winther, R&S, …)
• overlapping (Toselli, Hiptmair)
• substructuring (Toselli, Widlund, Wohlmuth, Hu, Zou)
• FETI-DP: plenary L11 (Toselli)

Recall Ap = Acurl + Adiv. Apply BCs and consider

) ( 1

curl curl 1

A R x x A x x A x

p T T

∈ ∀ ≤ ≤ α

Note: max α1 is smallest nonzero eig of Acurlw = λ Apw

0.7964 1/24 0.7966 1/20 0.7970 1/16 0.7977 1/12 0.7994 1/8 0.8064 1/4 α1 h

Note: PCG convergence depends

• n 1/α1 for exact solves w/ Ap

2D Structured Mesh

BDDC preconditioner for Ap with H/h = 4, rtol = 10-10

2.3 2.2 2.1 2.1 1.9 cond 15 15 14 13 12 iter 0.7964 1/24 0.7966 1/20 0.7970 1/16 0.7977 1/12 0.7994 1/8 0.8064 1/4 α1 h

Similar results for 3D problems, OS too

Recap of Talk

• Penalty Preconditioner

– approximate solution of related penalized problem – symmetric indefinite, but σ(M-1A) real and positive – CG for saddle point systems

• BDDC Preconditioner

– well suited for use with penalty preconditioner – constraint modification for near incompressibility – applicable to problems in H(grad), H(div), and H(curl)

• Divergence Stabilization

– permits use of H(grad) preconditioners for curl-curl