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

Discretization and Solution of and Solution of Discretization Convection- -Diffusion Problems Diffusion Problems Convection Howard Elman University of Maryland

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 1

The Convection-Diffusion Equation − ε ∇ + ⋅ ∇ = Ω ⊂ d = u w u f � d 2 in , 1 , 2 , 3 Boundary conditions: = ∂ Ω u g Ω ∂ Ω on D D N ∂ u ∂ Ω = ∂ Ω g on D ∂ n N N Inflow boundary: ∂Ω + = {x � ∂Ω | w � n > 0} Characteristic boundary: ∂Ω 0 = {x � ∂Ω | w � n = 0} Outflow boundary: ∂Ω - = {x � ∂Ω | w � n < 0} 2

The Convection-Diffusion Equation − ε ∇ + ⋅ ∇ = u w u f 2 Challenging / interesting case: ε 0 Reduced problem : w �� u = f , hyperbolic Streamlines : parameterized curves c(s) in Ω s.t. c(s) has tangent vector w(c(s)) on c d � ∇ ⋅ = = u w u c s f ( ) ( ( )) ds Solution u to reduced problem = solution to ODE If u ( s 0 ) � inflow boundary ∂Ω + , and u ( s 1 ) � ∂Ω , say outflow ∂Ω − , then boundary values are determined by ODE 3

Consequence − ε ∇ + ⋅ ∇ = u w u f 2 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: − ε + = u u ' ' ' 1 on ( 0 , 1 ) , = = u u ( 0 ) 0 , ( 1 ) 0 Solution = u x x ( ) � � − − ε x e / 1 − ε − − x e ( 1 / )( 1 ) � � ε − � e � 1 / 1 ≈ = x x except near 1 Solution to reduced equation 4

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

Conventions of Notation L = characteristic length scale in boundary e.g. length of inflow boundary x = x/L in normalized domain ˆ W = normalization for velocity (wind) w, e.g. w = W w *, where ||w * ||=1 In normalized variables: � � WL L � � 2 − ∇ + ⋅ ∇ = u � � w u f 2 � � ε ε � � * * * � � WL Peclet number , characterizes relative contributions P ≡ , ε of convection and diffusion 6

Reference Problems − ε ∇ + ⋅ ∇ = u w u f 2 1. w =(0,1) Dirichlet b.c. analytic solution � � − − ε − y e ( 1 ) / 1 = u x y x � � ( , ) − ε − e � � 2 / 1 2. w =(0,(1+(x+1) 2 /4) Neumann b.c. at outflow characteristic boundary layers 7

Reference Problems − ε ∇ + ⋅ ∇ = u w u f 2 3. w : 30 o left of vertical interior layer from discontinuous b.c. downstream boundary layer 4. w= recirculating (2y(1-x 2 ),-2x(1-y 2 ) characteristic boundary layers discontinuous b.c. 8

Weak Formulation − ε ∇ + ⋅ ∇ = u w u f 2 ∈ Ω ∈ Ω u H 1 v H 1 Find ( ) s.t. for all ( ), E E 0 � � � ε ∇ ⋅ ∇ + ⋅ ∇ = + u v w u v fv vg ( ) N Ω Ω ∂ Ω N Ω = = ∂ Ω 1 H v v g ( ) { | on } E D D Ω = = ∂ Ω H 1 v v ( ) { | 0 on } E D 0 Shorthand notation: a ( u , v ) = l ( v ) for all v Can show: a ( u,u ) � ε || � u || 2 = ε � Ω � u �� u Lax-Milgram lemma existence and uniqueness a ( u,v ) � ( ε +|| w || � L) || � u || || � v|| of solution l ( v ) � C || � v|| 9

Approximation by Finite Elements h ⊂ h ⊂ S H S H 1 1 Given finite dimensiona l , , E E E 0 0 ∈ h ∈ h u S v S find such that for all , h E h 0 � � � � ε ∇ ⋅ ∇ + ⋅ ∇ = + u u w u v f v v g ( ) h h h h h h N Ω Ω Ω ∂ Ω N a ( u h , v h ) = l ( v h ) for all v h 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 10

What happens in such cases? Problem 1, accurate Problem 1, inaccurate Problem 2, accurate Problem 2, inaccurate 11

Explanations 1. Error analysis: discrete solution is quasi-optimal : Γ ∇ − ≤ ∇ − u u w u v ( ) inf ( ) , ε h h h ∈ v S E h WL Γ = + = + ε P 1 1 , large if is small ε w Ph Wh = = P 2. Mesh Peclet number : ε h L 2 2 If P h >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 12

Revisit two examples Problem 1, exponential layer, width ~ ε Problem 2, characteristic layer, width ~ ε 1/2 13

Fix: The Streamline Diffusion Method Petrov-Galerkin method : change the test functions Galerkin: a ( u h , v h ) = l ( v h ) for all v h Petrov-Galerkin: a ( u h ,v h + δ w �� v h ) = l ( v h + δ w �� v h ) for all v h δ is a parameter Result: a sd ( u h , v h ) = l sd ( v h ) Streamline diffusion term � � � = ε ∇ ⋅ ∇ + ⋅ ∇ + δ ⋅ ∇ ⋅ ∇ a u v u u w u v w u w v ( , ) ( ) ( )( ) sd h h h h h h h h Ω Ω Ω � � − δε ∇ ⋅ ∇ 2 u w v ( ) ( ) 0 for linear/bilinear h h k ∆ k � � � = + δ ⋅ ∇ + + δ ⋅ ∇ l v f v f w v v w v g ( ) ( ) ( ) h h h h h N Ω Ω ∂ Ω N 14

The Streamline Diffusion Method Explained Augment finite element space: ˆ = + h S S B h h B h : bubble functions, with support local to element ˆ S h Principle: augmented space places basis functions in layers not resolved by the grid We could pose the problem on the augmented space: ˆ ˆ find u h in s.t. a(u h ,v h ) =l(v h ) for all v h in S S h h Then: decouple unknowns associated with bubble functions from system new problem on original grid 15

The Streamline Diffusion Method Explained Under appropriate assumptions: this new problem is a sd ( u h , v h ) = l sd ( v h ) � � = ε ∇ ⋅ ∇ + ⋅ ∇ a u v u u w u v ( , ) ( ) sd h h h h h h Ω Ω � � + δ ⋅ ∇ ⋅ ∇ w u w v ( )( ) k h h ∆ ∆ k k � � � = + δ ⋅ ∇ l v f v f w v ( ) ( ) h h k h k Ω ∆ k k determined from elimination of bubble functions Streamline diffusion 16

Compare Galerkin and Streamline Diffusion Top: accurate solution, ε =1/200 Middle: bilinear elements, Galerkin, 32 � 32 grid Bottom: bilinear elements, streamline diffusion, 32 � 32 grid 17

Error Bounds For Galerkin : as noted earlier, quasi-optimality: Γ ∇ − ≤ ∇ − u u w u v ( ) inf ( ) ε h h h ∈ v S E h More careful analysis : for linear/bilinear elements, ∇ − ≤ u u Ch D 2 u Large in exponential ( ) h boundary layers for small ε For streamline diffusion: use norm ( ) 1/2 2 2 ≡ ε ∇ + δ ⋅ ∇ v sd v w v Then − ≤ u u Ch D u 3 / 2 2 h sd 18

These bounds do not tell the whole story For one example (Problem 1, ε =1/64), compare errors || � (u-u h )|| on Ω and Ω * = (-1,1) � (-1,3/4) (to exclude boundary layer) Grid Galerkin Str.Diff. Galerkin Str.Diff. Ω Ω * Ω Ω ∗ (P h ) 8 � 8 5.62 4.34 3.25 8.16e-7 (8) 16 � 16 4.91 4.01 1.48 1.64e-5 (4) 32 � 32 3.81 3.23 5.30e-2 1.11e-5 (2) 64 � 64 2.39 2.39 4.98e-7 4.98e-7 (1) 19

Choice of parameter δ δ δ δ � � = ε ∇ ⋅ ∇ + ⋅ ∇ a u v u u w u v ( , ) ( ) Made element-wise: sd h h h h h h Ω Ω � � + δ ⋅ ∇ ⋅ ∇ w u w v ( )( ) k h h ∆ ∆ k k ( ) � h − k k > k P P � 1 1 / if 1 δ = w h h � 2 | | k k � k ≤ P � 0 if 1 h 20

Matrix Properties + n h n n h ϕ ϕ S S ∂ Given a basis { } for , extended by{ } for = j j j + E n 1 0 1 � = ϕ u Finite element function is u , problem h j j j becomes : find {u } such that j � ϕ ϕ = ϕ = a l i n ( , ) u ( ), 1 , 2 ,..., j i j j j or a sd Leads to matrix equation F u = f, F= ε A + N ( + S ) 21

Matrix Properties Matrix equation F u = f, F= ε A + N ( + S ) A= [ a ij ], a ij = � Ω � φ j �� φ i, , discrete Laplacian, symmetric positive-definite N= [ n ij ], n ij = � Ω ( w �� φ j ) φ i , discrete convection operator, skew-symmetric ( N=-N T ) S= [ s ij ], s ij = � Ω ( w �� φ j )( w �� φ i ) , discrete streamline upwinding operator, positive semi-definite 22

End of Part I Next: how to solve F u = f ? 23

Iterative Solution Algorithms: Krylov Subspace Methods System F u = 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 24