Howard Elman University of Maryland Vicki Howles John Shadid - - PowerPoint PPT Presentation

Howard Elman University of Maryland

2

• Vicki Howles
• David Kay
• Daniel Loghin
• Alison Ramage
• John Shadid
• David Silvester
• Ray Tuminaro
• Andy Wathen

3

Incompressible Navier-Stokes Equations

div grad ) grad (

2 t

= − = + ⋅ + ∇ − u f p u u u u ν α

α=0 ! steady state problem α=1 ! evolutionary problem

        =                 − f p u C B B F

T

Discretization and linearization Matrix equation x=b Goal: Robust general solution algorithms Easy to implement Derived from subsidiary building blocks Adaptible to a variety of scenarios (steady / evolutionary / Stokes / Boussinesq)

4

Outline

• 1. General approach

preconditioning for saddle point problems

• 2. Relation to traditional approaches

projection methods SIMPLE

• 3. Details for Navier-Stokes equations
• 4. Analytic / experimental results
• 5. Potential for more general problems

Preliminary: Steady Stokes Equations

div grad

2

= − = + ∇ − u f p u

, 0         =                 − f p u C B B A

T

Algebraic equation: Symmetric indefinite ! MINRES algorithm is applicable Preconditioning operator:

       

S

Q A

Generalized eigenvalue problem:

                =                 − p u Q A p u C B B A

S T

λ

A = discrete vector Laplacian Au+BTp = λ Au , λ ≠ 1 ) u= 1/( λ-1) A-1BTp Bu = λ QSp BA-1BTp=λ(λ-1) QSp Case C=0: Rusten &Winther, 1992 Silvester & Wathen, 1993

6

BA-1BTp=µ QSp, µ=λ(λ-1)

] ] [

• b

d

• a

c || || ) /( ) ( 1 ) /( ) ( 1 2 || ||

2 /

r ad bc ad bc r

k k

        + − ≤

[

Convergence bound for MINRES: Verfürth, 1984: For QS=pressure mass matrix, µ2 [as,bs] independent of discretization parameter h Under mapping µ λ=1§ (1+4µ)1/2, Computational requirements, for times a vector Poisson solve: can be approximated, e.g. with multigrid Mass matrix solve: cheap λ ½

1 −

       

S

Q A

7

Generalize to Navier-Stokes Equations

div grad ) grad (

) 1 ( ) ( ) 1 ( ) 1 ( ) ( ) 1 ( 2 ) 1 (

= − = + ⋅ + ∇ −

+ + + + + ∆ m m m m m m m t

u f p u u u u ν

α

Linearization via Picard iteration (slightly inaccurate notation):

      =             f p u B B F

T

Discretization Analogue of Stokes strategy: preconditioner

       

S

Q F

Same analysis ! BF-1BTp=µ QSp, µ=λ(λ-1) ) seek approximation QS to Schur complement N.B. Same question arises for other strategies for linearization

8

Under mapping µ λ=1§ (1+4µ)1/2, eigenvalues λ are clustered in two regions,

• ne on each side of imaginary axis

Suppose QS ¼ BF-1BT so that eigenvalues of BF-1BTp = µ QSp are tightly clustered. Can improve this: Observation: symmetry is important for Stokes solver MINRES: optimal with fixed cost per step ) need block diagonal preconditioner For Navier-Stokes: don’t have symmetry need Krylov subspace method for nonsymmetric matrices (e.g. GMRES)

9

        −

S T

Q B F

Alternative: block triangular preconditioner

                − =                 p u Q B F p u B B F

S T T

µ

! Generalized eigenvalue problem Fu+BTp = µ Fu , µ ≠ 1 ) u= -F-1BTp Bu = - µ QSp BF-1BTp=µ QSp µ ½ {1}[ Theorem (Fischer, Ramage, Silvester, Wathen): For preconditioned GMRES iteration, let p0 be arbitrary and u0 = F-1(f-BTu0) () r0=(0,w0)) (uk,pk)T be generated with block triangular preconditioner, (uk,pk)D be generated with block diagonal preconditioner. Then (u2k,p2k)D = (u2k+1,p2k+1)D = (uk,pk)T .

10

Computational requirements, to implement block triangular preconditioner: Solve Qs r = -q, then solve Fw = v-BTr The only difference from block diagonal solve: matrix-vector product BT r (negligible) For second step: convection-diffusion solve: can be approximated, e.g. with multigrid For first step: something new needed: QS

                − =        

−

q v Q B F r w

S T 1

Compute

11

        + −         =         −

− −

) (

1 1

C B BF B F I BF I C B B F

T T T

One more interpretation: )

        =         + −         −

− − −

I BF I C B BF B F C B B F

T T T 1 1 1

) (

1 −

        −

S T F

Q B Q

¼ Shows what is needed for stabilized discretization: QS ¼ BF-1BT + C

12

Relation to Projection Methods

“Classical” O(∆t) projection method (Chorin 1967, Temam 1969): Step 1:

( )

) ( 1 ) ( ) ( * 2 1 ) ( ) ( * 2 ) ( (*)

) grad ( ) grad (

m t m m t m m m

u u u f u I f u u u t u u

∆ ∆

+ ⋅ − = ∇ − = ⋅ + ∇ − ∆ − ν ν

In matrix form: (

)

) ( 1 ) ( * 1 m t m t

Mu Nu f u A M

∆ ∆

+ − = +ν

Step 2: In matrix form:

        =                 ∇ − ∇

∆ + + ∆ * 1 ) 1 ( ) 1 ( 1

u p u I

t m m t

        =                

∆ + + ∆ * 1 ) 1 ( ) 1 ( 1

Mu p u B B M

t m m T t

Performed via pressure-Poisson solve

13

Substitute u* from Step 1 into Step 2 :

( ) ( )

                − +

− ∆ − ∆ ∆

I B M I B M B B A M

T t T t t 1 1 1 1 1

ν

Perot, 1993 Contrast: update derived purely from linearization & discretization:

        + − =                 +

∆ + + ∆ ) ( 1 ) ( ) 1 ( ) 1 ( 1 m t m m m T t

Mu Nu f p u B B A M ν

( ) ( )

        +         + − +

− ∆ − ∆ ∆

I B A M I B A M B B A M

T t T t t 1 1 1 1 1

ν ν ν

( )( )

        + − =                 + +

∆ + + − ∆ ∆ ∆ ) ( 1 ) ( ) 1 ( ) 1 ( 1 1 1 1 m t m m m T t t t

Mu Nu f p u B B M A M A M ν ν

( )( )

) ( ) ( : Error

1 1 1 1

t O AB M t B M A M B

T T t t T

∆ = ∆ − = + −

− − ∆ ∆

ν ν

14

Dukowicz & Dvinsky 1992 Perot 1993 Quarteroni, Saleri &Veneziani 2000 Henriksen & Holman 2002 For higher order accuracy in time and related approaches:

15

Relation to SIMPLE

                − =        

− −

I B F I B BF B F B B F

T T T 1 1

                −

− −

I B F I B F B B Q

T T F

ˆ ˆ

1 1

¼ QF: approximate convection-diffusion solve F: diagonal part of F Many variants

^

Patankar & Spaulding, 1972

16

Perspective of New Approach

• Take on Schur complement directly
• Separate time discretization from algebraic algorithm
• Enable flexible treatment of time discretization, linearization

Allow choice of linearization Allow large time steps / CFL numbers if circumstances warrant

17

Approximation for the Schur Complement (I)

Suppose the gradient and convection-diffusion operators approximately commute (w=u(m)): Kay,Loghin, Wathen 2002

∇ ∇ ⋅ + ∇ − ≈ ∇ ⋅ + ∇ − ∇

u p

w w ) ( ) (

2 2

ν ν

Require pressure convection-diffusion operator Discrete analogue:

p p T u T T u u p p T u

M F B BM B BF B M F M F M B M

1 1 1 1 1 1 1 − − − − − − −

= ⇒ =

p

A

In practice: don’t have equality, instead

p p p S

M F A Q

1 −

≡

Requirements: Poisson solve

1 − p

A

Mass matrix solve

1 − p

M

+ Convection-diffusion solve

1 −

F

1

for

− S

Q

18

Evolutionary Equations

div grad ) grad (

2 t

= − = + ⋅ + ∇ − u f p u u u u ν div grad ) grad (

) 1 ( ) 1 ( ) 1 ( ) ( ) 1 ( 2 ) ( ) 1 (

= − = + ⋅ + ∇ − ∆ −

+ + + + + m m m m m m m

u f p u u u t u u ν

Backward Euler:

div ) ) grad ( ( grad ) ) grad ( (

) 1 ( ) ( ) ( ) ( 2 2 1 ) 1 ( ) 1 ( ) ( ) 1 ( 2 2 1 ) ( ) 1 (

= − ⋅ + ∇ − − = + ⋅ + ∇ − + ∆ −

+ + + + + m m m m m m m m m m

u u w u f p u w u t u u ν ν

Linearized 2nd order Crank-Nicolson (Simo & Armero):

) 1 ( ) ( ) (

5 . 5 . 1

−

− =

m m m

u u w

19

Matrix structure after discretization

N A M F B B F

v T

+ + =       ν α ,

Considerations are the same

,        

S T F

Q B Q

Preconditioner Convection-diffusion solve easier than for steady state

1 1 1

,

− − −

≈ = + + = ≈ ≈ S A F M Q N A M F S Q F Q

p p p S p p p p S F

ν α α

larger for CN than for BE For large α (small ∆t):

T h t d h t T u p p p S T h t T t d h p v

BB I h I B BM M F A Q BB B BF I N A M F

d d d

∆ ∆ − − ∆ − ∆

= = ≈ ⇒ + + = ) ( ) ( ~ ~

1 1 1

ν α

20

Approximation for the Schur Complement (II)

Elman 1999 Consider simple observation in linear algebra: let G, H be rectangular matrices Consider HT (GHT)-1G , maps R to range(HT) fixes range(HT)

G

n1 n2

1

n

HT (GHT)-1GT = I on range(HT) Take G=BF-1, H=B ) BT (BF-1BT)-1BF-1 = I on range(BT) BT (BF-1BT)-1B = F on range(F-1BT)

H

n1 n2

Suppose range(BT) ½ range(F-1BT) ) BT (BF-1BT)-1B = F on range(BT) (BBT) (BF-1BT)-1 (BBT)= BFBT (BF-1BT)-1¼ (BBT)-1(BFBT) (BBT)-1 ´ QS

• 1

21

Recapitulating: Two ideas under consideration

Preconditioners

        =                 − f p u C B B F

T

• 1. Fp preconditioner: QS = Ap Fp
• 1Mp

Requirements: Poisson solve, mass matrix solve, Fp on pressure space Decisions on boundary conditions

• 2. BFBt preconditioner: QS = (BBT) (BFBT)-1 (BBT)

Requirements: Two Poisson solves C=0

        −

S T F

Q B Q

for Requirement common to both: approximation of action of F-1 Overarching philosophy: subsidiary operations Poisson solve convection-diffusion solve are manageable

22

Benchmark problems

• 1. 2D Driven Cavity Problem
• 2. 2D Backward Facing Step
• 3. 3D Driven Cavity Problem

Various finite element / finite difference discretizations in space Backward Euler or Crank-Nicolson in time u1=u2=0, except u1=1 at top u1=u2=0, except u1=1-y2 at inflow ν ∂u1/∂x=p ∂u2/∂x=0 at outflow

23

Experiment 1: 2D driven cavity problem on [0,1]£[0,1]

Integrate from t=0 to t=1 Solve implicit systems with Fp-preconditioned GMRES Average iteration counts per linear solve

t ∆

1/8 1/16 1/32 1/64 6.9 5.6 4.0 2.9 8.4 6.9 5.1 3.6 9.3 8.1 6.2 4.3 9.9 8.6 6.9 5.0

1/40 1/80 1/160 1/320

ν

Discretization in space: marker-and-cell finite differences (Harlow & Welch), h=1/128 Discretization in time: backward Euler with various time steps 1/5000

ν

4-5 3-4

24

For same problem: GMRES behavior at t=1/4, 1/2, 3/4

h=1/64, ν = 1/160, 1/320

25

Experiment 2: 2D driven cavity flow on [-1,1]£[-1,1], Re=200

Discretization in space: Q2-Q1 finite elements Discretization in time: backward Euler with various time steps Iterations of GMRES at sample time / Picard step

26

Experiment 2, continued: Re=1000

27

Experiment 2, continued: Re=1000

28

Experiment 3: backward facing step, Re=200

Iterations of GMRES at sample time / Picard step Q2-Q1 fem spatial discretization Backward Euler time discretization

Residual reduction, ν=1/100, 32£96 grid Residual reduction, ν=1/100, 64£192 grid dt=1

dt=.1 ! CFL~ 83.2 dt=.1 ! CFL~ 41.6 For CFL=||u|| dt/h, ||u||¼26 )

29

Experiment 3, continued: backward facing step, Re=1000

dt=1 ! CFL~928 dt=.1 ! CFL~ 92.8 dt=1 ! CFL~464 dt=.1 ! CFL~ 46.4 For CFL = ||u||dt/h, ||u||¼29 )

30

Experiment 4: 3D driven cavity problem on [0,1]3

Marker-and-cell finite differences Pseudo-transient iteration: ten time steps at various CFL nos. and Re, h=1/64 Average iteration counts with Fp preconditioning to satisfy mild stopping criterion || rk|| · 10-2 ||f||, f=nonlinear residual CFL #

.1, .5, 1, 10, 50, 100 5000 10,000 50,000 2 5 6 9

Re

500 1000 5000 10,000 50,000 5 6 10

Iterations

31

Key aspect of computations:

Poisson solves: q = Ap

• 1p

Convection-diffusion solves: w = F-1v Each can be approximated using existing technology

• multigrid
• domain decomposition
• fast direct methods
• other iterative methods

required at each step

32

Experiment 4, continued:

Replace convection-diffusion solve and Poisson solve with multigrid approximations CFL #

50,000 9

Re

500 1000 50,000

Iterations Exact Inexact

12 3 Poisson 8 Conv-diff 10 13 3 Poisson 8 Conv-diff

33

Boundary conditions for preconditioners

To define operators Fp and Ap: need to “specify” boundary conditions on pressure space Derivation of preconditioners does not offer guidance Formulation of problem does: have convection-diffusion operator (-νr 2+(w¢r)) defined on pressure space, w=current velocity iterate No specific b.c. on pressures suggests “natural” condition ∂ p / ∂ n =0 But: flow problems require Dirichlet conditions on inflow boundary, where w¢n < 0

34

Therefore: formulate Fp using Dirichlet conditions p=0 on ∂ Ω- (inflow, w¢n<0) Neumann conditions ∂ p/∂ n =0 on ∂ Ω0 (characteristic, w¢n=0) ∂ Ω+ (outflow, w¢n>0) Comments:

• 1. Not really specifying values, just defining matrix Fp
• 2. Formulate Ap in compatible manner
• 3. This issue is important

but it only affects performance of solvers, not accuracy

35

For benchmark problems

• 1. 2D Driven Cavity Problem
• 2. 2D Backward Facing Step

u1=u2=0, except u1=1 at top u1=u2=0, except u1=1-y2 at inflow ν ∂u1/∂x=p ∂u2/∂x=0 at outflow u¢n´ 0 ) Fp defined using Neumann b.c u¢n < 0 at inflow ) Dirichlet b.c. there Otherwise Neumann

36

Analysis:

        − =         =

S T A T

Q B F Q B B F A ,

For solving AQA

• 1x = b using GMRES, assuming

AQA

• 1 =VΛ V-1

is diagonalizable:

|| || | ) ( | max min ) ( || ) ( || min || ||

) ( 1 ) ( 1 1 ) (

1

r p V r AQ p r

k AQ p A k p k

A k k

λ κ

σ λ

−

∈ = − =

≤ ≤

Here: For F_p preconditioning:

( ) ( )

|| || | ) ( |

1 1 2 1 2 || || 1 1

w c AQ c

t A w t

+ ≤ ≤

∆ − ∆ + ν ν ν

λ

) asymptotic convergence rate (||rk||/||r0||)1/k , k!1, independent of (small) h, ∆t pessimistic wrt ν (Loghin)

37

Generalizations:

Boussinesq equations: ! coefficient matrix “Ideal” preconditioner is For Picard iteration, H=0 and Schur complement is the same as for the Navier-Stokes equations

38

Add chemistry: molecular species with concentration Y

Add equation of form Coupled with ! coefficient matrix

39

Concluding remarks

Goal: develop strategies to handle linearized Navier-Stokes equations in a flexible manner

• Allow large time steps if stiffness is not critical
• Respect coupling of velocities and pressures
• Automatically adapt to handle different scenarios

(creeping flow, stiff systems, steady problems) Technical approach:

• Take advantage of saddle point structure of problem
• Develop preconditioners for Schur complement and

accompanying systems

This document was created with Win2PDF available at http://www.daneprairie.com. The unregistered version of Win2PDF is for evaluation or non-commercial use only.