Discretization and Solution of and Solution of Discretization - - PowerPoint PPT Presentation

Discretization Discretization and Solution of and Solution of Convection Convection-

• Diffusion Problems

Diffusion Problems

Howard Elman University of Maryland

1

Overview

• 1. The convection-diffusion equation

Introduction and examples

• 2. Discretization strategies

Finite element methods Inadequacy of Galerkin methods Stabilization: streamline diffusion methods

• 3. Iterative solution algorithms

Krylov subspace methods Splitting methods Multigrid

2

The Convection-Diffusion Equation

3 , 2 , 1 , in

2

= ⊂ Ω = ∇ ⋅ + ∇ − d f u w u

d

• ε

Boundary conditions:

N N D D

g n u g u Ω ∂ = ∂ ∂ Ω ∂ =

• n
• n

D N

Ω ∂ Ω ∂ Ω

Inflow boundary: ∂Ω+ = {x∂Ω | w n > 0} Characteristic boundary: ∂Ω0 = {x∂Ω | w n = 0} Outflow boundary: ∂Ω- = {x∂Ω | w n < 0}

3

The Convection-Diffusion Equation

f u w u = ∇ ⋅ + ∇ −

2

ε

Challenging / interesting case: ε Reduced problem: w u = f , hyperbolic Streamlines: parameterized curves c(s) in Ω s.t. c(s) has tangent vector w(c(s)) on c

f s c u ds d w u = = ⋅ ∇

• ))

( ( ) (

Solution u to reduced problem = solution to ODE If u(s0) inflow boundary ∂Ω+, and u(s1)∂Ω, say outflow ∂Ω− , then boundary values are determined by ODE

4

Consequence

2

f u w u = ∇ ⋅ + ∇ −ε

For small ε, solution to convection-diffusion equation often has boundary layers, steep gradients near parts of ∂Ω. Αlso: discontinuities at inflow propagate into Ω along streamlines Simple (1D) example of first phenomenon:

1 near except 1 1 ) (

/ 1 / ) 1 )( / 1 (

= ≈

• −

− − =

− − −

x x e e e x x u

x x ε ε ε

Solution

) 1 ( , ) ( , ) 1 , (

• n

1 ' ' ' = = = + − u u u u ε

Solution to reduced equation

5

Additional Consequence These layers (steep gradients) are difficult to resolve with discretization

6

Conventions of Notation

f L u w WL u

• =

∇ ⋅

• +

∇ − ε ε

2 * * * 2

In normalized variables:

, ε WL P ≡

Peclet number, characterizes relative contributions

• f convection and diffusion

L = characteristic length scale in boundary e.g. length of inflow boundary = x/L in normalized domain W = normalization for velocity (wind) w, e.g. w = W w*, where ||w*||=1

x ˆ

7

Reference Problems

• −

− =

− − − ε ε / 2 / ) 1 (

1 1 ) , ( e e x y x u

y 2

f u w u = ∇ ⋅ + ∇ −ε

• 1. w=(0,1)

Dirichlet b.c. analytic solution

• 2. w=(0,(1+(x+1)2/4)

Neumann b.c. at outflow characteristic boundary layers

8

Reference Problems

2

f u w u = ∇ ⋅ + ∇ −ε

• 3. w: 30o left of vertical

interior layer from discontinuous b.c. downstream boundary layer

• 4. w=recirculating

(2y(1-x2),-2x(1-y2) characteristic boundary layers discontinuous b.c.

9

Weak Formulation

2

f u w u = ∇ ⋅ + ∇ −ε

}

• n

| { ) ( }

• n

| { ) ( ) ( ), ( all for s.t. ) ( Find

1 1 1 1 D E D D E N E E

v v H g v v H vg fv v u w v u H v H u

N

Ω ∂ = = Ω Ω ∂ = = Ω + = ∇ ⋅ + ∇ ⋅ ∇ Ω ∈ Ω ∈

• Ω

∂ Ω Ω

ε

Shorthand notation: a(u,v) =l(v) for all v Can show: a(u,u) ε ||u||2 = εΩ u u a(u,v) (ε+||w|| L) ||u|| ||v|| l(v) C ||v|| Lax-Milgram lemma existence and uniqueness

• f solution

10

Approximation by Finite Elements

a(uh,vh) =l(vh) for all vh Typically: finite element spaces are defined by low-order basis functions, e.g. — linear or quadratic functions on triangles — bilinear or biquadratic functions on quadrilaterals

• Ω

∂ Ω Ω Ω

+ = ∇ ⋅ + ∇ ⋅ ∇ ∈ ∈ ⊂ ⊂

N

N h h h h h h h h h E h E h E h E

g v v f v u w u u S v S u H S H S ) ( , all for such that find , , l dimensiona finite Given

1 1

ε

11

What happens in such cases?

Problem 1, accurate Problem 2, accurate Problem 1, inaccurate Problem 2, inaccurate

12

Explanations

small is if large , 1 1 , ) ( inf ) (

w

h

ε ε ε P WL v u u u

h S v w h

h E

+ = + = Γ − ∇ Γ ≤ − ∇

∈

• 1. Error analysis: discrete solution is quasi-optimal:
• 2. Mesh Peclet number:

2 2 ε Wh L Ph P

h

= =

If Ph>1, then — there are oscillations in the discrete solution — these become pronounced if mesh does not resolve layers — oscillations propagate into regions where solution is smooth — problem is most severe for exponential boundary layers

13

Revisit two examples

Problem 1, exponential layer, width ~ ε Problem 2, characteristic layer, width ~ ε1/2

14

Fix: The Streamline Diffusion Method

Petrov-Galerkin method: change the test functions Galerkin: a(uh,vh) =l(vh) for all vh Petrov-Galerkin: a(uh,vh+δ w vh) =l(vh+δ w vh) for all vh δ is a parameter Result: asd(uh,vh) =lsd(vh)

) ( ) ( ) ( ) ( ) ( ) )( ( ) ( ) , (

2

• Ω

∂ Ω Ω ∆ Ω Ω Ω

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

N k

N h h h h h h k h h h h h h h h h sd

g v w v v w f v f v l v w u v w u w v u w u u v u a δ δ δε δ ε

Streamline diffusion term 0 for linear/bilinear

15

The Streamline Diffusion Method Explained

Augment finite element space: Bh: bubble functions, with support local to element

h

B S S + =

h h

ˆ

We could pose the problem on the augmented space: find uh in s.t. a(uh,vh) =l(vh) for all vh in Then: decouple unknowns associated with bubble functions from system new problem on original grid Principle: augmented space places basis functions in layers not resolved by the grid

h

ˆ S

h

ˆ S

h

ˆ S

16

The Streamline Diffusion Method Explained

Under appropriate assumptions: this new problem is

) ( ) ( ) )( ( ) ( ) , (

• ∇

⋅ + = ∇ ⋅ ∇ ⋅ + ∇ ⋅ + ∇ ⋅ ∇ =

∆ Ω ∆ ∆ Ω Ω k h k h h h h k h h h h h h sd

v w f v f v l v w u w v u w u u v u a

k k k

δ δ ε

asd(uh,vh) =lsd(vh)

k determined from elimination of bubble functions

Streamline diffusion

17

Compare Galerkin and Streamline Diffusion

Middle: bilinear elements, Galerkin, 32 32 grid Top: accurate solution, ε=1/200 Bottom: bilinear elements, streamline diffusion, 32 32 grid

18

Error Bounds

) ( inf ) (

h

h S v w h

v u u u

h E

− ∇ Γ ≤ − ∇

∈

ε

For Galerkin: as noted earlier, quasi-optimality: More careful analysis: for linear/bilinear elements,

) (

2u

D Ch u u

h

≤ − ∇

Large in exponential boundary layers for small ε For streamline diffusion: use norm

( )

1/2 2 2

v w v v sd ∇ ⋅ + ∇ ≡ δ ε

2 2 / 3

u D Ch u u

sd h

≤ −

Then

19

These bounds do not tell the whole story

4.98e-7 4.98e-7 2.39 2.39 6464 (1) 1.11e-5 5.30e-2 3.23 3.81 3232 (2) 1.64e-5 1.48 4.01 4.91 1616 (4) 8.16e-7 3.25 4.34 5.62 88 (8) Str.Diff. Ω∗ Galerkin Ω* Str.Diff. Ω Galerkin Ω Grid (Ph ) For one example (Problem 1, ε=1/64), compare errors ||(u-uh)|| on Ω and Ω* = (-1,1)(-1,3/4) (to exclude boundary layer)

20

Choice of parameter δ δ δ δ

• ∆

∆ Ω Ω

∇ ⋅ ∇ ⋅ + ∇ ⋅ + ∇ ⋅ ∇ =

k k

h h k h h h h h h sd

v w u w v u w u u v u a ) )( ( ) ( ) , ( δ ε

Made element-wise:

( )

• ≤

> − = 1 if 1 if / 1 1 | | 2

k h k h k h k k k

P P P w h δ

21

Matrix Properties

n i l a u S S

j j i j j j j h h E n n n j h n j j

,..., 2 , 1 ), ( u ) , ( such that } {u find : becomes problem , u is function element Finite for } by{ extended , for } { basis a Given

j j 1 1

= = =

• ∂

+ + =

ϕ ϕ ϕ ϕ ϕ ϕ

• r asd

Leads to matrix equation Fu=f, F=ε A + N (+ S)

22

Matrix Properties

Matrix equation Fu=f, F=ε A + N (+ S) A=[aij], aij=Ω φjφi, , discrete Laplacian, symmetric positive-definite N=[nij], nij=Ω (wφj)φi , discrete convection operator, skew-symmetric (N=-NT) S=[sij], sij=Ω (wφj)(wφi) , discrete streamline upwinding operator, positive semi-definite

23

End of Part I

Next: how to solve Fu=f ?

24

Iterative Solution Algorithms: Krylov Subspace Methods

System Fu=f

• F is a nonsymmetric matrix, so an appropriate Krylov subspace

method is needed

• Examples:
• GMRES
• GMRES(k) restarted
• BiCGSTAB
• BiCGSTAB(l)
• Our choices:
• Full GMRES for optimal algorithm, or
• BiCGSTAB(2) for suboptimal

25

Properties of Krylov Subspace Methods

Drawback of GMRES: work & storage requirements at step k are proportional to kN BiCGSTAB: Fixed cost per step, independent of k Drawback: No convergence analysis Variant: BiCGSTAB(l), more robust for complex eigenvalues, somewhat higher cost per step (l=2), but still fixed

26

Convergence of GMRES

GMRES: Starting with u0, with residual r0=f-Fu0, computes uk span{r0, Fr0, …,Fk-1r0} for which rk = f-Fuk satisfies ||rk|| = minpk(0)=1 || pk(F)r0||. Consequence: Theorem: For diagonalizable F=VΛV-1, ||rk|| ||V|| ||V-1|| minpk(0)=1 maxλσ(F)|pk(λ)| ||r0||. (||rk||/||r0||)1/k (||V|| ||V-1||)1/k (minpk(0)=1}maxλσ(F)|pk(λ)|)1/k

ρ ˆ

Loosely speaking: residual is reduced by factor of at each step Want eigenvalues to lie in compact set

ρ ˆ

27

Convergence of GMRES

) ( 1 1 ˆ

1 2 2 2 F F R

Q a d d a a

−

= ≤ − + − + ≈ ρ ρ

1

. .

1+d 1+a

Size of convergence factor ρ

ˆ

28

Key for Fast Convergence: Preconditioning Splitting operators

Seek QF F such that

• the approximation is good, and
• it is inexpensive to apply the action of Q-1 to a vector

Splitting: F=QF-RF stationary iteration Error ek = u-uk satisfies ) (

1 1

f u R Q u

k F F k

+ =

− +

( )

) ( ) ( / ) ( ) ( ) )( ( ) (

1 / 1 1 / 1 1 1 1 1 1

F Q I F Q I e e e F Q I e e F Q I e u u F Q I u u R Q u u

F k k F k k k F k k F k k F k F F k − − − − − − +

− ≈ − ≤ − ≤ − =

• −

− = − = − ρ

29

Preconditioning / Splitting operators

Thus: want to be as small as possible Equivalently: eigenvalues of = eigenvalues of ) (

1F

Q I

F −

− ρ F QF

1 − 1 − F

FQ as close to 1 as possible

1

) (

1 F F R

Q− ρ This is similar to the requirement for rapid convergence of GMRES Solve u Q u f u FQ f Q Fu Q

F F F F

ˆ , ˆ

• r

1 1 1 1 − − − −

= = =

30

Examples of splitting operators

=

F

Q lower triangle of

A

Gauss-Seidel =

F

Q block lower triangle of

A

Line Gauss-Seidel

F F F

U L Q = Symmetric versions Incomplete LU factorization

F

Q LU F = ≈ Comments:

• All depend on ordering of underlying grid
• Symmetric versions (symmetric GS, ILU)

take some account of underlying flow

• Line/block versions can handle irregular grids

31

Convergence Analysis (Parter & Steuerwalt)

Seek maximal eigenvalue of u u R Q

F F

λ =

−1

• r

u R u Q

F F

= λ Subtract

u RF λ from both sides

u R h u R u R Q

F F F F

h

) ( ) (

2 2

1 1

• =
• =

−

− − λ λ λ λ

u R h Fu

F h

) (

2

µ = Suggests relation to

u u µ =

• u

w u u ∇ ⋅ + ∇ − =

2

ε to be determined

• For many examples of splittings:

) ( ), , ( ) , (

2

x r r v u r v u R h

h h F

= ≈

• perator"

tion multiplica weak " a is

2 F

R h and (2)

• f

eigenvalue minimal

) (

= → µ µh (1) (2) (defines )

32

are known:

• +
• +

=

• +
• +

+ = =

2 2 2 ) ( 2 2 2 2 2 2 / 2 / ) , (

2 2 2 2 2 ) ( ) sin( ) sin(

2 1

ε ε π ε µ ε ε π ε µ π π

y x y x jk y w x w k j

w w r w w k j r y k e x j e u

Consequence:

2 ) ( 1

1 ) ( h R Q

F F

µ ρ − =

−

For model problems: (i) r is constant (will demonstrate in a moment) (ii) on square domains, eigenvalues, eigenvectors of u u µ =

To find r: consider centered finite differences

• +

− 2 h wx ε

• −

− 2 h wx ε

• −

− 2 h wy ε

• +

− 2 h wy ε ε 4

For horizontal line Jacobi splitting, = block tridiagonal:

F

Q ) ( 2 2 2 ] [

2 1 , 1 ,

h O u u h w u h w u R

ij j i y j i y ij F

+ =

• +

+

• −

=

− +

ε ε ε

, 2ε =

• r

Key point: convection terms lead to smaller convergence factors

2 2 2 2

2 2 2 2 1 1 h w w

y x

• +
• +

− = ε ε π ρ

34

Comments / Extensions

• Similar results obtained from matrix/Fourier analysis
• Young:
• “Multi-line” (k-line) splittings

r =

• Can extend to other splittings via matrix comparison theorems

(Varga-Wonicki):

k / 2ε

2

• perator)

Jacobi (line

• perator)

Seidel

• Gauss

(line ρ ρ =

) ( ) (

1

• 1

1 2

• 1

2 1 1 1 2

R Q R Q Q Q ρ ρ ≤

• ≥

− −

35

x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x

6 5 4 3 2 1

Limitations of analysis above:

It does not discriminate among different orderings

Natural ordering of grid Left-to-right, bottom-to-top Plus resulting matrix structure: Horizontal line red-black Ordering and matrix structure:

x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x x

6 3 5 2 4 1

Young theory: spectral radii (Jacobi or Gauss-Seidel) independent of ordering Performance of GS: depends on ordering

36

Example: Problem 1

60 P elements, linear piecewise , ) 1 , (

2 2

= = + ∇ − f u u

y

ε

Four solution strategies: line Gauss-Seidel iteration with natural line ordering, following the flow (bottom-to-top) natural line ordering, against the flow (top-to-bottom) red-black line ordering, with the flow red-black line ordering, against the flow

10-5 10-4 10-3 10-2 10-1 100 101 102 10 20 30 40 50 60

*

• *
• *
• k

||e^(k)|| Natural with flow Red- black Natural against flow With Against

37

Ordering effects

• =

− =

− 1

) ( satisfies Error e R Q e u u e

k F F k k k 1 1

) ( ) ( e R Q e R Q e

k F F k F F k − −

≤ =

• “Classical” analysis only provides

insight in asymptotic sense: 1 ) ( lim

/ 1 1

=

− ∞ →

ρ

k k F F k

R Q

• E. & Chernesky:

bounds for for 1D problems

k F F R

Q ) (

1 −

38

Practical consequences

For nonconstant flows: inherent latencies if sweeps don't follow flow Possible fixes:

• flow-directed orderings (Bey & Wittum, Kellogg, Hackbusch, Xu)
• iterations based on multi-directional sweeps

2D version: Speeds convergence when recirculations are present ) ( ) ( ) ( ) (

4 / 3 1 4 4 / 3 1 k 2 / 1 1 3 2 / 1 3/4 k 4 / 1 1 2 4 / 1 1/2 k 1 1 1/4 k + − + + + − + + + − + + − +

− + = − + = − + = − + =

k k k k k k k k

Fu f Q u u Fu f Q u u Fu f Q u u Fu f Q u u

Contours of stream function

39

Summarizing with an experiment:

Eigenvalues of line-GS preconditioned operator, vertical flow, P=40, h=1/32 Asymptotic convergence rate is faster with Krylov acceleration However: does not overcome latencies

40

Multigrid

Flow-following methods are effective for convection-dominated problems: But: ultimately, solvers discussed above are mesh dependent GMRES performance for Problem 4

41

Multigrid

V-cycle multigrid: end iteration) next for (update h) (postsmoot ) ( steps, for update) and correction (prolong ˆ ˆ ˆ problem coarse to system multigrid apply residual) (restrict ) ( r ˆ ) (presmooth ) ( steps, for e convergenc until for Choose

1 1 1 2 T 1 1 i i F i F i i i h i F i F i

u u f Q u F Q I u m e P u u r e F Fu f P f Q u F Q I u k i u ← + − ← + ← = − = + − ← =

+ − − − −

42

Bottom Line: Performance

4 3 2 1 4 3 2 1 Grid 2 3 3 3 2 5 3 3 128 128 3 3 3 3 2 4 3 3 6464 5 4 3 3 2 4 3 3 32 32 6 6 3 3 3 3 3 4 16 16 Example Example ε=1/25 ε=1/200 Multigrid iterations for ||rk||/||r0||<10-6

43

For this to happen:

) ( : Smoothing 1.

1 1

f Q u F Q I u

F i F i − −

+ − ← Two things have to be done correctly r e F

h

ˆ ˆ : solve grid Coarse 2.

2

= Smoother must take underlying flow into account For results above for Problem 4 (recirculating wind): smoother is 4-directional Coarse grid operator must be stable Even if fine grid is “fine enough,” coarse grid operators should include streamline diffusion

44

Example / effect

• f smoother

After one four-directional Gauss-Seidel step After four one-directional Gauss-Seidel steps

45

Concluding Remarks

• Discretization requires stabilization for convection-dominated

problems

• The best solution algorithms combine
• general techniques of iterative methods
• splitting strategies coupled to the underlying physics
• stabilization when needed

46

References

• J. J. H. Miller, E. O’Riordan and G. I. Shishkin,

Fitted Numerical Methods for Singularly Perturbed Problems, World Scientific, 1995.

• K. W. Morton, Numerical Solution of Convection-Diffusion

Problems, Chapman & Hall, 1996.

• H.-G. Roos, M. Stynes and L. Tobiska, Numerical Methods

for Singularly Perturbed Differential Equations, Springer, 1996.

• H. C. Elman, D. J. Silvester and A. J. Wathen, Finite Elements

and Fast Iterative Solvers, Oxford University Press, 2005.