A Local Projection Stabilization method with shock capturing and - - PowerPoint PPT Presentation

A Local Projection Stabilization method with shock capturing and diagonal mass-matrix for solving non-stationary transport dominated problems

Friedhelm Schieweck1 and Piotr Skrzypacz 1

1Institut für Analysis und Numerik

Otto-von-Guericke-Universität Magdeburg

”Martin-60 workshop”, Dresden, 16-18th Dec. 2011

LPS-method with shock capturing and diagonal mass-matrix 1

Model Problem

1D model problem: Ω := (0, 1), 0 ≤ ε ≪ 1 ut − εuxx + b · ux + cu = f in Ω × (0, T), b.c. u(0, t) = g0(t) , u(1, t) = g1(t) for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω

LPS-method with shock capturing and diagonal mass-matrix 2

Model Problem

1D model problem: Ω := (0, 1), 0 ≤ ε ≪ 1 ut − εuxx + b · ux + cu = f in Ω × (0, T), b.c. u(0, t) = g0(t) , u(1, t) = g1(t) for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω asymptotic solution in the case b = 1, c = 0, f = 0 idea is the transformation: w(x, t) := u(x + bt, t) then PDE + initial condition are equivalent to wt − εwxx = in R × (0, T), w(x, 0) = ˜ u0(x) for x ∈ R (˜ u0 = extension of u0 by 0) ⇒ uε(x, t) = w(x − t, t) = 1 √ 4πεt 1 ˜ u0(y) exp

• −(x − t − y)2

4εt

• dy

LPS-method with shock capturing and diagonal mass-matrix 2

Quality of asymptotic solution uε

problem: Ω := (0, 1), 0 ≤ ε ≪ 1 ut − εuxx + ux = in Ω × (0, T), b.c. u(0, t) = 0 , u(1, t) = for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω

LPS-method with shock capturing and diagonal mass-matrix 3

Quality of asymptotic solution uε

problem: Ω := (0, 1), 0 ≤ ε ≪ 1 ut − εuxx + ux = in Ω × (0, T), b.c. u(0, t) = 0 , u(1, t) = for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω uε satisfies pde + initial condition, but boundary conditions are not satisfied; on the boundary ∂Ω it holds: |uε(x, t) − 0| ≤ Cεm ∀ m ∈ N, ∀ t ∈ (t1, t2)

LPS-method with shock capturing and diagonal mass-matrix 3

Quality of asymptotic solution uε

problem: Ω := (0, 1), 0 ≤ ε ≪ 1 ut − εuxx + ux = in Ω × (0, T), b.c. u(0, t) = 0 , u(1, t) = for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω uε satisfies pde + initial condition, but boundary conditions are not satisfied; on the boundary ∂Ω it holds: |uε(x, t) − 0| ≤ Cεm ∀ m ∈ N, ∀ t ∈ (t1, t2) from maximum principle we get: |uε(x, t) − u(x, t)| ≤ Cεm ∀ m ∈ N, ∀ (x, t) ∈ Ω × (t1, t2)

LPS-method with shock capturing and diagonal mass-matrix 3

Semi-discretization in space - part 1

1D model problem: Ω := (0, 1) ut − εuxx + b · ux + cu = f in Ω × (0, T), b.c. u(0, t) = g0(t) , u(1, t) = g1(t) for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω for 0 ≤ ε ≪ 1

LPS-method with shock capturing and diagonal mass-matrix 4

Semi-discretization in space - part 1

1D model problem: Ω := (0, 1) ut − εuxx + b · ux + cu = f in Ω × (0, T), b.c. u(0, t) = g0(t) , u(1, t) = g1(t) for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω for 0 ≤ ε ≪ 1 Galerkin part aG(u, v) := ε(ux, vx) + (b · ux, v) + (cu, v)

LPS-method with shock capturing and diagonal mass-matrix 4

Semi-discretization in space - part 1

1D model problem: Ω := (0, 1) ut − εuxx + b · ux + cu = f in Ω × (0, T), b.c. u(0, t) = g0(t) , u(1, t) = g1(t) for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω for 0 ≤ ε ≪ 1 Galerkin part aG(u, v) := ε(ux, vx) + (b · ux, v) + (cu, v) new semi-discretization: Find t → uh(t) ∈ Vh such that uh(0) = πhu0,

• dtuh(t), vh
• + ah
• ˜

uh(t), uh(t), vh

• = (f, vh)

∀vh ∈ Vh where ah(˜ uh, u, v) := aG(u, v) + aLPS(u, v) + aSC(˜ uh, u, v)

LPS-method with shock capturing and diagonal mass-matrix 4

Semi-discretization in space - part 2

(time exact) ansatz: uh(t) =

Nh

• j=1

Uj(t) ∈ R ϕj , uh(t) ∈ Vh , Nh = dim (Vh)

LPS-method with shock capturing and diagonal mass-matrix 5

Semi-discretization in space - part 2

(time exact) ansatz: uh(t) =

Nh

• j=1

Uj(t) ∈ R ϕj , uh(t) ∈ Vh , Nh = dim (Vh) Find t → uh(t) ∈ Vh such that uh(0) = πhu0, (dtuh(t), vh) + ah(˜ uh(t), uh(t), vh) = (f(t), vh) ∀vh ∈ Vh

LPS-method with shock capturing and diagonal mass-matrix 5

Semi-discretization in space - part 2

(time exact) ansatz: uh(t) =

Nh

• j=1

Uj(t) ∈ R ϕj , uh(t) ∈ Vh , Nh = dim (Vh) Find t → uh(t) ∈ Vh such that uh(0) = πhu0, (dtuh(t), vh) + ah(˜ uh(t), uh(t), vh) = (f(t), vh) ∀vh ∈ Vh equivalent to nonlinear ODE-system (time exact): Find t → U(t) :=

• Uj(t)
• ∈ V := RNh such that

MdtU(t) = F

• t, U(t)
• ,

U(0) = U0 where

• F
• t, U(t)
• i =
• f(t), ϕi
• − a
• uh( ˜

U(t)), uh(U(t)), ϕi

• LPS-method with shock capturing and diagonal mass-matrix

5

Time discretization by dG(1)

0 = t0 < t1 < . . . < tN = T, I = ∪N

n=1In,

In := (tn−1, tn] U(t) := (Uj(t)) ∈ RNh =: V tn,1, tn,2 ∈ (tn−1, tn] from 2-point Gauss-Radau formula Uj

n = uτ(tn,j),

j = 1, 2 dG(1)-method: uτ

• In(t) = U1

nφn,1(t) + U2 nφn,2(t)

Find U1

n, U2 n ∈ V s.t.

3

4M + τn 2 A1 1 4M

−9

4M 5 4M + τn 2 A2

U1

n

U2

n

• =
• MU0

n + τn 2 f(tn,1)

−MU0

n + τn 2 f(tn,2)

• matrix Aj ∼ ah(uh( ˜

U(tn,j)), ·, ·) = ah(uh(U(tn−1)), ·, ·)

LPS-method with shock capturing and diagonal mass-matrix 6

desirable FE-spaces

wish for FE-space:

∃ L2-orthogonal local basis ⇒ mass-matrix M is diagonal !!! example: non-conforming Crouzeix-Raviart element (P1/P0)

LPS-method with shock capturing and diagonal mass-matrix 7

desirable FE-spaces

wish for FE-space:

∃ L2-orthogonal local basis ⇒ mass-matrix M is diagonal !!! example: non-conforming Crouzeix-Raviart element (P1/P0)

idea: construct in 1D an L2-orthogonal basis on ˆ K = [−1, 1]

hierarchical basis of P2 : ˜ ϕ1(x) = 1

2(1 − x) ,

˜ ϕ2(x) = 1

2(1 + x) ,

ϕ3(x) = 1 − x2 new enriched bubble-function of P5 : (theory see [Ma/Sk/To ’07]) ϕ4(x) = (1 − x2)x(1 + cx2), c ∈ R suitably chosen !!!

LPS-method with shock capturing and diagonal mass-matrix 7

desirable FE-spaces

wish for FE-space:

∃ L2-orthogonal local basis ⇒ mass-matrix M is diagonal !!! example: non-conforming Crouzeix-Raviart element (P1/P0)

idea: construct in 1D an L2-orthogonal basis on ˆ K = [−1, 1]

hierarchical basis of P2 : ˜ ϕ1(x) = 1

2(1 − x) ,

˜ ϕ2(x) = 1

2(1 + x) ,

ϕ3(x) = 1 − x2 new enriched bubble-function of P5 : (theory see [Ma/Sk/To ’07]) ϕ4(x) = (1 − x2)x(1 + cx2), c ∈ R suitably chosen !!! ansatz: ϕi(x) = ˜ ϕi(x) + aiϕ3(x) + biϕ4(x) i = 1, 2 determine c, ai, bi such that: (ϕj, ϕi) = 0 ∀ i = j

LPS-method with shock capturing and diagonal mass-matrix 7

basis functions of the new Q+

2 element in 1D

−1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1 −0.4 −0.2 0.2 0.4 0.6 0.8 1 x L2−orthogonal basis for LPS−Q2

+

φ1 φ2 φ3 φ4 φ1 φ2 φ3 φ4

LPS-method with shock capturing and diagonal mass-matrix 8

tensor product basis of the new Q+

2 element in 2D

−1 −0.5 0.5 1 −1 −0.5 0.5 1 −0.2 0.2 0.4 0.6 0.8 1 1.2 x φ 3 y −1 −0.5 0.5 1 −1 −0.5 0.5 1 −0.2 0.2 0.4 0.6 0.8 1 1.2 x φ 7 y −1 −0.5 0.5 1 −1 −0.5 0.5 1 −0.4 −0.3 −0.2 −0.1 0.1 0.2 0.3 x φ 8 y −1 −0.5 0.5 1 −1 −0.5 0.5 1 −0.1 −0.05 0.05 0.1 x φ16 y

LPS-method with shock capturing and diagonal mass-matrix 9

LPS = Local Projection Stabilization

Galerkin bilinear form: aG(u, v) := ε (∇u, ∇v) + ((b · ∇)u, v) + (cu, v) LPS stabilization: aLPS(u, v) := c0

• K

hK

• ∇u − πh(∇u), ∇v − πh(∇v)
• K

where πh : L2(Ω) → Dh := Pdc

1

(local projection space)

LPS-method with shock capturing and diagonal mass-matrix 10

LPS = Local Projection Stabilization

Galerkin bilinear form: aG(u, v) := ε (∇u, ∇v) + ((b · ∇)u, v) + (cu, v) LPS stabilization: aLPS(u, v) := c0

• K

hK

• ∇u − πh(∇u), ∇v − πh(∇v)
• K

where πh : L2(Ω) → Dh := Pdc

1

(local projection space) sufficient for theory: local inf-sup-condition 0 < c2 ≤ inf

wL∈Dh

sup

vh∈Bh(Vh)

(wL, vh)K wL0,Kvh0,K ∀ K. . . . is fulfilled for our new Q+

2 -element

LPS-method with shock capturing and diagonal mass-matrix 10

LPS = Local Projection Stabilization

Galerkin bilinear form: aG(u, v) := ε (∇u, ∇v) + ((b · ∇)u, v) + (cu, v) LPS stabilization: aLPS(u, v) := c0

• K

hK

• ∇u − πh(∇u), ∇v − πh(∇v)
• K

where πh : L2(Ω) → Dh := Pdc

1

(local projection space) sufficient for theory: local inf-sup-condition 0 < c2 ≤ inf

wL∈Dh

sup

vh∈Bh(Vh)

(wL, vh)K wL0,Kvh0,K ∀ K. . . . is fulfilled for our new Q+

2 -element

main advantage of LPS compared with SDFEM: no coupling with dtu term

LPS-method with shock capturing and diagonal mass-matrix 10

exploiting diagonal mass-matrix M

semi-discrete conv.-diff. problem: Mdtu(t) = f(t) − Au(t) dG(1) in Q+

2 basis:

find U1

n, U2 n ∈ V s.t.

• αM + τn

2 A1

βM γM δM + τn

2 A2

U1

n

U2

n

• =
• R1

n

R2

n

• eliminate U1

n from eq. (2):

U1

n = 1 γ M−1R2 n − 1 γ (δI + τn 2 M−1A2)U2 n

LPS-method with shock capturing and diagonal mass-matrix 11

exploiting diagonal mass-matrix M

semi-discrete conv.-diff. problem: Mdtu(t) = f(t) − Au(t) dG(1) in Q+

2 basis:

find U1

n, U2 n ∈ V s.t.

• αM + τn

2 A1

βM γM δM + τn

2 A2

U1

n

U2

n

• =
• R1

n

R2

n

• eliminate U1

n from eq. (2):

U1

n = 1 γ M−1R2 n − 1 γ (δI + τn 2 M−1A2)U2 n

insert in eq. (1) ⇒ to solve: ˜ A U2

n = ˜

R ⇒ U1

n = . . .

˜ A = βI − 1

γ

• αI + τn

2 M−1A1

• δI + τn

2 M−1A2

• sparse !!!

˜ R = M−1R1

n − 1 γ

• αI + τn

2 M−1A1

• M−1R2

n

LPS-method with shock capturing and diagonal mass-matrix 11

Shock capturing based on postprocessing

Observation: LPS-solution is still deteriorated by spurious

• scillations at the boundary layer.

LPS-method with shock capturing and diagonal mass-matrix 12

Shock capturing based on postprocessing

Observation: LPS-solution is still deteriorated by spurious

• scillations at the boundary layer.

Remedy: use additional nonlinear shock capturing term aSC(˜ uh, u, v) =

• K∈T ∗

h (˜

uh)

γK(˜ uh)(ux, vx)K where γK(˜ uh) := γ0hK ˜ uh − R1

hπ0 h˜

uhL∞(K) and R1

hπ0 h : Vh → P1

LPS-method with shock capturing and diagonal mass-matrix 12

Shock capturing based on postprocessing

Observation: LPS-solution is still deteriorated by spurious

• scillations at the boundary layer.

Remedy: use additional nonlinear shock capturing term aSC(˜ uh, u, v) =

• K∈T ∗

h (˜

uh)

γK(˜ uh)(ux, vx)K where γK(˜ uh) := γ0hK ˜ uh − R1

hπ0 h˜

uhL∞(K) and R1

hπ0 h : Vh → P1

R1

hπ0 h is a low order postprocessing based on the transfer operator

Rr

h of [Sch. 2000, preprint],

T ∗

h (˜

uh) =set of “bad cells“ (with high postprocessing error ˜ uh − R1

hπ0 h˜

uhL∞(K) )

LPS-method with shock capturing and diagonal mass-matrix 12

Postprocessing operator

Example of P1-postprocessing:

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 x u ihu π0

hu

Rh

1π0 hu

let π0

h : L2(Ω) → Pdisc

, R1

h : Pdisc

→ P1 ∩ H1

0(Ω)

LPS-method with shock capturing and diagonal mass-matrix 13

Postprocessing operator

Example of P1-postprocessing:

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.2 0.4 0.6 0.8 1 x u ihu π0

hu

Rh

1π0 hu

let π0

h : L2(Ω) → Pdisc

, R1

h : Pdisc

→ P1 ∩ H1

0(Ω)

let uh := ihu ∈ P1, ¯ uh = π0

huh ∈ Pdisc

and up

h = R1 h¯

uh ∈ P1 ∩ H1

0(Ω) ,

u − up

h0 ≤ Ch2

⇒ for smooth uh-parts the post-proc. error is small

LPS-method with shock capturing and diagonal mass-matrix 13

A 1D example for time-dependent transport

problem: Ω := (0, 1), 0 ≤ ε ≪ 1 ut − εuxx + ux = in Ω × (0, T), b.c. u(0, t) = 0 , u(1, t) = for t ∈ (0, T) i.c. u(x, 0) = u0(x) for x ∈ Ω h = (1/10)2−ℓ, τ = ch, ε = 0, Ω0(t) = (0.1 + t, 0.4 + t) ℓ maxn (u − uh,τ)(tn)L2(Ω0(tn))

• rder

1 6.513e − 04 2 1.539e − 05 5.40e + 00 3 1.496e − 06 3.36e + 00 4 1.884e − 07 2.99e + 00 5 2.355e − 08 3.00e + 00

LPS-method with shock capturing and diagonal mass-matrix 14

A 2D example for time-dependent transport

development of a shock: Ω = (0, 1)2, ν = 0, h = 1

32, τ = 1 128

convection-field: (velocity, wind) find c : Ω × [0, T] → R such that dtc(t) + v · ∇c(t) = 0, Ω × (0, T) c(t) = g(t), ∂Ω × [0, T] c(0) = 0, Ω × {0} animation: exact solution vs. LPS solution animation: exact error vs. error indicator (by post-processing) animation: exact solution vs. LPS+(shock-capturing) solution

LPS-method with shock capturing and diagonal mass-matrix 15

• nly for Martin

LPS-method with shock capturing and diagonal mass-matrix 16

• nly for Martin

Theorem [properties of Martin Stynes]

Martin Stynes is: honorable, a good friend, a good mathematician.

LPS-method with shock capturing and diagonal mass-matrix 16

• nly for Martin

Theorem [properties of Martin Stynes]

Martin Stynes is: honorable, a good friend, a good mathematician. Proof : It follows obviously from our experience.

• LPS-method with shock capturing and diagonal mass-matrix

16

• nly for Martin

Theorem [properties of Martin Stynes]

Martin Stynes is: honorable, a good friend, a good mathematician. Proof : It follows obviously from our experience.

• Martin, all the best in the period until the

”Martin-70 workshop” !!!

LPS-method with shock capturing and diagonal mass-matrix 16