12 Implicit Spatial Discretization for Advection-Diffusion-Reaction - - PowerPoint PPT Presentation

12

◭ ◭ ◭ ◮ ◮ ◮

1/35

◭ ◭ ◭ ◮ ◮ ◮

1/35

Implicit Spatial Discretization for Advection-Diffusion-Reaction Equation

Kundan Kumar 10-Dec-2008

12

◭ ◭ ◭ ◮ ◮ ◮

2/35

◭ ◭ ◭ ◮ ◮ ◮

2/35

Introduction Applications of Advection-Diffusion Reaction Equations

Chemical Vapor Deposition

12

◭ ◭ ◭ ◮ ◮ ◮

3/35

◭ ◭ ◭ ◮ ◮ ◮

3/35

Introduction

Setting:

• Advection-Diffusion-Reaction Equation
• φt + uφx = ǫφxx + s(x, t),
• Advection Velocity : u
• Diffusion Coefficient : ǫ
• Source term : s(x, t)

s(x, t) = b2ǫ cos(b(x − ut)).

12

◭ ◭ ◭ ◮ ◮ ◮

4/35

◭ ◭ ◭ ◮ ◮ ◮

4/35

Introduction

Setting:

• Exact Solution:

φ = cos(b(x − ut)) + exp(−a2ǫt) cos(a(x − ut)).

• Dirichlet Boundary Conditions.
• Initial Condition, φ(x, t) at t = 0.

12

◭ ◭ ◭ ◮ ◮ ◮

5/35

◭ ◭ ◭ ◮ ◮ ◮

5/35

Contents

1 Discretization 6 1.1 Order Condition . . . . . . . . . . . . . . . . . . . . . . . . 8 2 Examples 10 3 Stability 16 4 Time Integration Aspect 18 5 Numerical Computations 19 6 Conclusion 35

12

◭ ◭ ◭ ◮ ◮ ◮

6/35

◭ ◭ ◭ ◮ ◮ ◮

6/35

1. Discretization

φt + uφx = ǫφxx + s(x, t)

Discretization:

1

• k=−1

βkw′

j+k(t) = h−2 1

• k=−1

αkwj+k(t) +

1

• k=−1

βkgj+k(t) wj(t) ≈ φ(xj, t); gj(t) = s(xj, t);

1

• k=−1

βk = 1.

12

◭ ◭ ◭ ◮ ◮ ◮

7/35

◭ ◭ ◭ ◮ ◮ ◮

7/35

Discretization

Vector Notation:

Bw′(t) = Aw(t) + Bg(t), A = (aij) = (h−2αj−i) B = (bij) = (βj−i).

Define:

ξk = (−1)kα−1 + α1, ηk = (−1)kβ−1 + β1.

12

◭ ◭ ◭ ◮ ◮ ◮

8/35

◭ ◭ ◭ ◮ ◮ ◮

8/35

1.1. Order Condition

Let φh be the restriction of the exact solution φ to the grid. Spatial truncation error:

σh(t) = Bφ′

h(t) − Aφh(t) − Bg(t).

Truncation error in a point (xj, t) equals:

σh,j(t) = h−2(C0φ + hC1φx + h2C2φxx + h3C3φxxx + ...)|(xj,t)

Order Condition: The discretization has order q if:

σh = O(hq),

translates to:

Ck = O(hq+2−k), k = 0, 1, · · · , q + 2.

12

◭ ◭ ◭ ◮ ◮ ◮

9/35

◭ ◭ ◭ ◮ ◮ ◮

9/35

Order Condition

Error coefficients:

C0 = −ξ0, C1 = −ξ1 − uhη0, Ck = −1 k! (ξk + kuhηk−1 − k(k − 1)ǫηk−2); k ≥ 2.

where,

ξk = (−1)kα−1 + α1, ηk = (−1)kβ−1 + β1.

Use the order condition to determine αj and βj.

12

◭ ◭ ◭ ◮ ◮ ◮

10/35

◭ ◭ ◭ ◮ ◮ ◮

10/35

2. Examples

Explicit Central Difference

w′

j = u

2h(wj−1 − wj+1) + ǫ h2(wj−1 − 2wj + wj+1) + gj,

Implicit Central Difference

1 6(w′

j−1 + 4w′ j + w′ j+1) =

u 2h(wj−1 − wj+1) + ǫ h2(wj−1 − 2wj + wj+1) + 1 6(gj−1 + 4gj + gj+1)

12

◭ ◭ ◭ ◮ ◮ ◮

11/35

◭ ◭ ◭ ◮ ◮ ◮

11/35

Examples

Define:

µ = uh/ǫ (Peclet Number).

Explicit Adaptive Upwinding

w′

j = u

2h(wj−1 − wj+1) + ǫ + 0.5uhκ h2 (wj−1 − 2wj + wj+1) + gj,

Where κ is defined as:

κ = max(0, 1 − 2/µ).

12

◭ ◭ ◭ ◮ ◮ ◮

12/35

◭ ◭ ◭ ◮ ◮ ◮

12/35

Examples

Implicit Adaptive Upwinding

1 2κw′

j−1 + (1 − 1

2κ)w′

j =

u 2h(wj−1 − wj+1) + ǫ + 0.5uhκ h2 (wj−1 − 2wj + wj+1) + 1 2κgj−1 + (1 − 1 2κ)gj.

12

◭ ◭ ◭ ◮ ◮ ◮

13/35

◭ ◭ ◭ ◮ ◮ ◮

13/35

Examples

Peclet Number µ:

µ = uh/ǫ.

Explicit Exponential Fitting

w′

j = 1

• k=−1

αkwj+k + gj, α−1 = uh exp(µ) exp(µ) − 1, α1 = uh 1 exp(µ) − 1, α0 = −(α1 + α−1).

12

◭ ◭ ◭ ◮ ◮ ◮

14/35

◭ ◭ ◭ ◮ ◮ ◮

14/35

Implicit Exponential Fitting

β−1w′

j−1 + β0w′ j + β1w′ j+1 = 1

• k=−1

αkwj+k + β−1gj−1 + β0gj + β1gj+1.

where

β−1 = 1 2

• exp(µ)

exp(µ) − 1 − 1 µ

• ,

β0 = 1 2, β1 = 1 2 1 µ − 1 exp(µ) − 1

• .

12

◭ ◭ ◭ ◮ ◮ ◮

15/35

◭ ◭ ◭ ◮ ◮ ◮

15/35

Examples

Compact Schemes:

α−1 = ǫ + 1 2uh − uh(β1 − β−1), α1 = ǫ − 1 2uh − uh(β1 − β−1), α0 = −(α−1 + α1), β−1 = 1 γ(6 + 3µ − µ2), β0 = 1 γ(60 − 4µ2), β1 = 1 γ(6 − 3µ − µ2)

and γ is a scaling factor given by:

γ = 72 − 6µ2.

12

◭ ◭ ◭ ◮ ◮ ◮

16/35

◭ ◭ ◭ ◮ ◮ ◮

16/35

3. Stability

Requirement:

|| exp(tB−1A)|| ≤ C,

for all t > 0. We can write:

A = V diag(ak)V −1, B = V diag(bk)V −1,

with ak, bk eigenvalues of A, B respectively. Define global error e(t):

e(t) = φh(t) − w(t), ˆ e(t) = V −1e(t)

Discretization error σh(t):

σh(t) = Bφ′

h(t) − Aφh(t) − Bg(t),

ˆ σh(t) = V −1σh(t).

12

◭ ◭ ◭ ◮ ◮ ◮

17/35

◭ ◭ ◭ ◮ ◮ ◮

17/35

Stability

The error equation then reads:

bk d dtˆ e(t) = akˆ e(t) + ˆ σh(t).

Stability if:

Re(ak/bk) ≤ 0

and

|ak| + |bk| > 0.

Result: For the three point scheme considered with C0 = C1 = 0, C2 =

O(h), and assume that: h−2|α0| + |β0 − 1 2| > 0,

then the stability condition holds iff:

2ah(β1 − β−1) ≥ α0,

and

α0(1 − 2β0) ≥ 0.

12

◭ ◭ ◭ ◮ ◮ ◮

18/35

◭ ◭ ◭ ◮ ◮ ◮

18/35

4. Time Integration Aspect

Ode system:

Bw′(t) = Aw + Bg(t).

Define:

F(t, w) = Aw(t) + Bg(t).

• We use the θ method (with θ = 0.5:

Bwn+1 = Bwn + 0.5τF(tn, wn) + 0.5τF(tn+1, wn+1).

• With Explicit method, there is some amount of ’implicitness’!.
• Stability conditions in general become more stringent in case of implicit

discretization method.

• For an implicit A-stable ODE method for time stepping, little difference

between the two methods.

12

◭ ◭ ◭ ◮ ◮ ◮

19/35

◭ ◭ ◭ ◮ ◮ ◮

19/35

5. Numerical Computations

Error for Implicit vs Explicit Central Difference

12

◭ ◭ ◭ ◮ ◮ ◮

20/35

◭ ◭ ◭ ◮ ◮ ◮

20/35

Implicit vs Explicit Adaptive Upwinding

12

◭ ◭ ◭ ◮ ◮ ◮

21/35

◭ ◭ ◭ ◮ ◮ ◮

21/35

Implicit vs Explicit Exponential Fitting

12

◭ ◭ ◭ ◮ ◮ ◮

22/35

◭ ◭ ◭ ◮ ◮ ◮

22/35

Implicit vs Explicit

12

◭ ◭ ◭ ◮ ◮ ◮

23/35

◭ ◭ ◭ ◮ ◮ ◮

23/35

Implicit vs Explicit Central Difference

12

◭ ◭ ◭ ◮ ◮ ◮

24/35

◭ ◭ ◭ ◮ ◮ ◮

24/35

Implicit vs Explicit Adaptive Upwinding

12

◭ ◭ ◭ ◮ ◮ ◮

25/35

◭ ◭ ◭ ◮ ◮ ◮

25/35

Implicit vs Explicit Exponential Fitting

12

◭ ◭ ◭ ◮ ◮ ◮

26/35

◭ ◭ ◭ ◮ ◮ ◮

26/35

12

◭ ◭ ◭ ◮ ◮ ◮

27/35

◭ ◭ ◭ ◮ ◮ ◮

27/35

12

◭ ◭ ◭ ◮ ◮ ◮

28/35

◭ ◭ ◭ ◮ ◮ ◮

28/35

12

◭ ◭ ◭ ◮ ◮ ◮

29/35

◭ ◭ ◭ ◮ ◮ ◮

29/35

12

◭ ◭ ◭ ◮ ◮ ◮

30/35

◭ ◭ ◭ ◮ ◮ ◮

30/35

12

◭ ◭ ◭ ◮ ◮ ◮

31/35

◭ ◭ ◭ ◮ ◮ ◮

31/35

12

◭ ◭ ◭ ◮ ◮ ◮

32/35

◭ ◭ ◭ ◮ ◮ ◮

32/35

12

◭ ◭ ◭ ◮ ◮ ◮

33/35

◭ ◭ ◭ ◮ ◮ ◮

33/35

12

◭ ◭ ◭ ◮ ◮ ◮

34/35

◭ ◭ ◭ ◮ ◮ ◮

34/35

12

◭ ◭ ◭ ◮ ◮ ◮

35/35

◭ ◭ ◭ ◮ ◮ ◮

35/35

6. Conclusion

When do we use Implicit Spatial Discretization?

• To achieve higher order without using wider stencils.
• To reduce the artificial oscillations in the numerical solution.
• Provides extra degrees of freedom for the numerical scheme.

Disadvantages

• Positivity may be lost.
• Stringent conditions for explicit time integration methods.