hp -version discontinuous Galerkin methods for - - PowerPoint PPT Presentation

hp-version discontinuous Galerkin methods for advection-diffusion-reaction problems

• n polytopic meshes

Zhaonan Dong

Department of Mathematics, University of Leicester, UK Joint work with A. Cangiani, E. H. Georgoulis (Leicester) and P. Houston (Nottingham)

Polytopal Element Methods in Mathematics and Engineering Georgia Tech, October 26-28, 2015

1 / 29

Model Problem

Let Ω be a bounded open polyhedral domain in Rd, d = 2, 3. Consider the problem Lu ≡ −

d

• i,j=1

∂j(aij(x)∂iu) +

d

• i=1

bi(x)∂iu + c(x)u = f (x), (1) where c ∈ L∞(Ω), f ∈ L2(Ω), and b := (b1, b2, . . . , bd)⊤ ∈ [W 1

∞(Ω)]d.

a = {aij}d

i,j=1 such that aij ∈ L∞(¯

Ω), with ξ⊤a(x)ξ≥0 ∀ξ ∈ Rd, a.e x ∈ ¯ Ω. (2)

2 / 29

Model Problem

Let Ω be a bounded open polyhedral domain in Rd, d = 2, 3. Consider the problem Lu ≡ −

d

• i,j=1

∂j(aij(x)∂iu) +

d

• i=1

bi(x)∂iu + c(x)u = f (x), (1) where c ∈ L∞(Ω), f ∈ L2(Ω), and b := (b1, b2, . . . , bd)⊤ ∈ [W 1

∞(Ω)]d.

a = {aij}d

i,j=1 such that aij ∈ L∞(¯

Ω), with ξ⊤a(x)ξ≥0 ∀ξ ∈ Rd, a.e x ∈ ¯ Ω. (2) Under hypothesis (2), (1) is termed a partial differential equation with non-negative characteristic form. This class includes elliptic, parabolic, hyperbolic, mixed-type PDEs.

2 / 29

Model Problem

Let Γ signify the union of its (d − 1)-dimensional open faces of ∂Ω.

Ω

Γ+ Γ− ΓD ΓN ΓD b n⊤a(x)n = 0 n⊤a(x)n > 0

Problem (1) is well-posed under (standard) hypothesis: c0(x)2 := c(x) − 1 2∇ · b(x) ≥ γ0 > 0 a.e x ∈ Ω. (3)

3 / 29

Model Problem

Let Γ signify the union of its (d − 1)-dimensional open faces of ∂Ω.

Ω

Γ+ Γ− ΓD ΓN ΓD b n⊤a(x)n = 0 n⊤a(x)n > 0

Problem (1) is well-posed under (standard) hypothesis: c0(x)2 := c(x) − 1 2∇ · b(x) ≥ γ0 > 0 a.e x ∈ Ω. (3) Aim: Derive and analyse hp-version dG methods on extremely general polytopic meshes for this problem.

3 / 29

Challenges and their resolution

Classical hp-inverse estimates not sharp for polytopic elements with arbitrarily small (d − k)-dim faces, k = 1, . . . , d − 1. Inverse estimates dictate IPDG penalisation constant.

4 / 29

Challenges and their resolution

Classical hp-inverse estimates not sharp for polytopic elements with arbitrarily small (d − k)-dim faces, k = 1, . . . , d − 1. Inverse estimates dictate IPDG penalisation constant. → new sharp hp-inverse estimates for polytopic elements with arbitrarily small (d − k)-dim faces

4 / 29

Challenges and their resolution

Classical hp-inverse estimates not sharp for polytopic elements with arbitrarily small (d − k)-dim faces, k = 1, . . . , d − 1. Inverse estimates dictate IPDG penalisation constant. → new sharp hp-inverse estimates for polytopic elements with arbitrarily small (d − k)-dim faces No sharp hp-approximation results for L2-projector over polytopic

• meshes. L2-projector key for first order terms Houston, Schwab & S¨

uli (’02).

4 / 29

Challenges and their resolution

Classical hp-inverse estimates not sharp for polytopic elements with arbitrarily small (d − k)-dim faces, k = 1, . . . , d − 1. Inverse estimates dictate IPDG penalisation constant. → new sharp hp-inverse estimates for polytopic elements with arbitrarily small (d − k)-dim faces No sharp hp-approximation results for L2-projector over polytopic

• meshes. L2-projector key for first order terms Houston, Schwab & S¨

uli (’02).

→ error analysis via hp-inf-sup stability on stronger norms.

` a la Johnson & Pitk¨ aranta (’86), Buffa, Hughes & Sangalli (’06), Ayuso & Marini (’09),...

4 / 29

The Finite Element Space for DGFEMs

DGFEMs are not restricted to employing standard polynomial spaces mapped from a reference frame. In this work, the FE space are constructed with following properties: Polynomial basis are defined over physical meshes On each element κ, we employ the total degree space Pp(κ).

5 / 29

The Finite Element Space for DGFEMs

DGFEMs are not restricted to employing standard polynomial spaces mapped from a reference frame. In this work, the FE space are constructed with following properties: Polynomial basis are defined over physical meshes On each element κ, we employ the total degree space Pp(κ). Definition of Finite Element Space: mesh T : disjoint subdivision of Ω into polygons/polyhedra; discretisation parameters:

◮ mesh size hκ := diam(κ), κ ∈ T ; ◮ polynomial degree vector p : (pκ ∈ N, κ ∈ T )

The DG finite element space Sp

T with respect to T and p defined by:

Sp

T := {u ∈ L2(Ω) : u|κ ∈ Ppκ(κ), κ ∈ T },

5 / 29

Notation for DGFEMs

Average { {·} } of v and q at x ∈ F: { {v} } := 1

2(vi + vj),

F ⊂ Γint vi, F ⊂ Γ , { {q} } := 1

2(qi + qj),

F ⊂ Γint qi, F ⊂ Γ . Jump of v and q at x ∈ F: [ [v] ] := vi ni + vj nj, F ⊂ Γint vi ni, F ⊂ Γ , [ [q] ] := qi · ni + qj · nj, F ⊂ Γint qi · ni, F ⊂ Γ . Upwind jump of (scalar-valued) function v across a face F: ⌊v⌋ := v+

κ − v− κ ,

F ⊂ ∂−κ\Γ v+

κ ,

F ⊂ Γ v+

κ and v− κ upwind and downwind w.r.t b values.

6 / 29

Weak formulation of DGFEM

For simplicity of presentation, we assume a ∈ [S0

T ]d×d sym .

and that b · ∇ξ ∈ Sp

T , for ξ ∈ Sp T . (important for p-version bounds.)

7 / 29

Weak formulation of DGFEM

For simplicity of presentation, we assume a ∈ [S0

T ]d×d sym .

and that b · ∇ξ ∈ Sp

T , for ξ ∈ Sp T . (important for p-version bounds.)

Then the IP DGFEM is given by: find uh ∈ Sp

T such that

B(uh, vh) = ℓ(vh), (4) for all vh ∈ Sp

T . Here, the bilinear form B(·, ·) : Sp T × Sp T → R is defined as

the sum of two parts B(u, v) := Bar(u, v) + Bd(u, v),

7 / 29

Weak formulation of DGFEM

The first part Bar accounts for the advection and reaction terms: Bar(u, v) :=

• κ∈T
• κ
• b · ∇u + cu
• v dx

−

• κ∈T
• ∂−κ\Γ

(b · n)⌊u⌋v+ ds −

• κ∈T
• ∂−κ∩(ΓD∪Γ−)

(b · n)u+v+ ds, and the second part Bd is responsible for the diffusion term: Bd(u, v) :=

• κ∈T
• κ

a∇u · ∇v dx +

• Γint∪ΓD

σ[ [u] ] · [ [v] ] ds −

• Γint∪ΓD
• {

{a∇hu} } · [ [v] ] + { {a∇hv} } · [ [u] ]

• ds.

8 / 29

Weak formulation of DGFEM

The linear functional ℓ : Sp

T → R is defined by

ℓ(v) :=

• κ∈T
• κ

fv dx −

• κ∈T
• ∂−κ∩(ΓD∪Γ−)

(b · n)gDv+ ds −

• ΓD

gD

• (a∇hv) · n − σv
• ds +
• ΓN

gNv ds. Remark: The discontinuity-penalization parameter σ ∈ L∞(Γint ∪ ΓD) will be defined later.

9 / 29

hp-inverse estimate from trace

F κF

♭

For v ∈ Pp(κ), a classical inverse estimate on the simplex κF

♭ gives:

v2

L2(F) p2|F|

|κF

♭ | v2 L2(κF

♭ ) p2|F|

|κF

♭ | v2 L2(κ).

Note: This is not sharp when |F| → 0 !

10 / 29

hp-inverse estimate from trace

F κF

♭

For v ∈ Pp(κ), a classical inverse estimate on the simplex κF

♭ gives:

v2

L2(F) p2|F|

|κF

♭ | v2 L2(κF

♭ ) p2|F|

|κF

♭ | v2 L2(κ).

Note: This is not sharp when |F| → 0 ! Using some developments from Georgoulis (’2008), we can show the following inverse estimate which remains sharp when |F| → 0:

10 / 29

Inverse estimate lemma

Lemma (Trace-inverse estimate)

If κ ∈ ˜ T then for each v ∈ Pp(κ), we have the inverse estimate v2

L2(F) ≤ CINV(p, κ, F)p2|F|

|κ| v2

L2(κ),

(5) with CINV(p, κ, F) := Cinv          min

• |κ|

supκF

♭ ⊂κ |κF

♭ |, p2d

• ,

if κ ∈ ˜ T , |κ| supκF

♭ ⊂κ |κF

♭ |,

if κ ∈ T \ ˜ T ,

11 / 29

Inverse estimate lemma

Lemma (Trace-inverse estimate)

If κ ∈ ˜ T then for each v ∈ Pp(κ), we have the inverse estimate v2

L2(F) ≤ CINV(p, κ, F)p2|F|

|κ| v2

L2(κ),

(5) with CINV(p, κ, F) := Cinv          min

• |κ|

supκF

♭ ⊂κ |κF

♭ |, p2d

• ,

if κ ∈ ˜ T , |κ| supκF

♭ ⊂κ |κF

♭ |,

if κ ∈ T \ ˜ T ,

Lemma (H1-inverse estimate)

If for each v ∈ Pp(κ), we have the inverse estimate ∇v2

L2(κ) ≤ ˜

Cinv p4 h2

κ

v2

L2(κ),

(6) where the constant ˜ Cinv is independent of discretization parameter hκ, pκ.

11 / 29

Main Result

Let penalisation σ : Γ\ΓN → R+ be defined facewise by σ(x) :=        Cσ max

κ∈{κ1,κ2}

• CINV(pκ, κ, F) ¯

aκp2

κ|F|

|κ|

• ,

x ∈ F ⊂ Γint, F = ∂κ1 ∩ ∂κ2, CσCINV(pκ, κ, F) ¯ aκp2|F| |κ| , x ∈ F ⊂ ΓD, F = ∂κ ∩ ΓD, (7) with Cσ > 0 large enough, and independent of p, |F|, and |κ|.

12 / 29

Main Result

Let penalisation σ : Γ\ΓN → R+ be defined facewise by σ(x) :=        Cσ max

κ∈{κ1,κ2}

• CINV(pκ, κ, F) ¯

aκp2

κ|F|

|κ|

• ,

x ∈ F ⊂ Γint, F = ∂κ1 ∩ ∂κ2, CσCINV(pκ, κ, F) ¯ aκp2|F| |κ| , x ∈ F ⊂ ΓD, F = ∂κ ∩ ΓD, (7) with Cσ > 0 large enough, and independent of p, |F|, and |κ|. Introducing the streamline-DG norm: |v|2

s := |v|2 DG +

• κ∈T

τκb · ∇v2

L2(κ),

(8) where |v|2

DG

:=

• κ∈T
• c0v2

L2(κ) + 1

2v +2

∂−κ∩(ΓD∪Γ−) + 1

2v + − v −2

∂−κ\Γ

+ 1 2v +2

∂+κ∩Γ

• +
• κ∈T

√a∇v2

L2(κ) +

• Γint∪ΓD

σ|[ [v] ]|2 ds, with τκ depending on discretisation parameters.

12 / 29

Main Result

Theorem

There exist a positive constant Λs which is independent of the mesh size h and polynomial degree p, such that: inf

ν∈Sp

T \{0}

sup

µ∈Sp

T \{0}

˜ B(ν, µ) |ν|s|µ|s ≥ Λs. (9) Test function µ = Tν = ν + α

κ∈T τκ(b · ∇ν) ≡ ν + ανs.

Theorem

For each κ ∈ T , such that Eu|K ∈ Hlκ(K), where K ∈ T♯ with κ ⊂ K and sκ = min{pκ + 1, lκ}. Then, the following bound holds: |u − uh|2

s ≤ C

• κ∈T

h2sκ

κ

p2lκ

κ

• Gκ(F, Cm, pκ, τκ) + Dκ(F, CINV, Cm, pκ)
• Eu2

Hlκ(K),

13 / 29

Main Result

If hκ = h, pκ = p, and diam(F) ∼ hκ, h⊥

κ

∼ hκ, F ⊂ ∂κ, κ ∈ T . For pure hyperbolic problem (a = 0): |u − uh|ar ≤ |u − uh|s ≤ C hs− 1

2

pl−1 uHl(Ω). Which is h optimal and p sub-optimal by 1/2 order. For pure elliptic problem (b = 0): |u − uh|d ≤ C hs−1 pl− 3

2

uHl(Ω). Which is h optimal and p sub-optimal by 1/2 order.

14 / 29

Numerical example (Hyperbolic)

For Ω = [−1, 1]2, b = (2 − y2, 2 − x), and c = 1 + (1 + x)(1 + y)2 the forcing function f is chosen to satisfy the analytic solution u(x, y) = 1 + sin(π(1 + x)(1 + y)2/8).

−1 −0.5 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

256 standard polygonal mesh

Figure: 256 polygonal mesh.

15 / 29

Numerical example (Hyperbolic)

Dof1/2

101 102 103

||u-uh ||L2(Ω)

10-15 10-10 10-5 100 p=1 p=2 p=3 p=4 p=5 p=6

DGFEM DGFEM(P) DGFEM(Q)

Dof1/2

101 102 103

||| u-uh|||DG

10-10 10-5 100 p=1 p=2 p=3 p=4 p=5 p=6

DGFEM DGFEM(P) DGFEM(Q)

(a) (b)

Figure: Convergence of DGFEM under h–refinement for p = 1, 2, 3, 4, 5, 6. (a) u − uhL2(Ω) ∼ hp+1; (b) |u − uh|DG ∼ hp+1/2

16 / 29

Numerical example (Hyperbolic)

Dof1/2

20 40 60 80 100 120 140

||u-uh ||L2(Ω)

10-15 10-10 10-5 100 64 Elements 256 Elements 64 Elements 256 Elements

DGFEM DGFEM(P) DGFEM(Q)

Dof1/2

20 40 60 80 100 120 140

||| u-uh|||DG

10-14 10-12 10-10 10-8 10-6 10-4 10-2 100 64 Elements 256 Elements

DGFEM DGFEM(P) DGFEM(Q)

Figure: Convergence of DGFEM under p–refinement.

17 / 29

Numerical example (Mixed type)

We consider the PDE problem: −x2uyy + ux + u = 0 for − 1 ≤ x ≤ 1, y > 0; ux + u = 0 for − 1 ≤ x ≤ 1, y ≤ 0, which has analytic solution: u(x, y) =

• sin( 1

2π(1 + y)) exp(−(x + π2x3 12 ))

for − 1 ≤ x ≤ 1, y > 0; sin( 1

2π(1 + y)) exp(−x)

for − 1 ≤ x ≤ 1, y ≤ 0, This problem is hyperbolic in the area y ≤ 0 and is parabolic in the area y > 0. The accuracy of the method is influenced by whether the mesh contains the line y = 0 in its skeleton or not.

18 / 29

Numerical example (Mixed type)

−1 −0.5 0.5 1 −1 −0.8 −0.6 −0.4 −0.2 0.2 0.4 0.6 0.8 1

256 modified polygonal mesh

Figure: 256 modified polygonal mesh.

19 / 29

Numerical example (Mixed type)

Dof1/2

20 40 60 80 100 120 140

||u-uh ||L2(Ω)

10-12 10-10 10-8 10-6 10-4 10-2 100 64 Elements 256 Elements

DGFEM DGFEM(P) DGFEM(Q)

Dof1/2

20 40 60 80 100 120 140

||| u-uh|||DG

10-12 10-10 10-8 10-6 10-4 10-2 100 64 Elements 256 Elements

DGFEM DGFEM(P) DGFEM(Q)

Figure: Convergence of DGFEM under p–refinement.

20 / 29

Numerical example (3D)

Let Ω be (0, 1)3 and set a ≡ 0, b = (−y, z, x), c = xy2z, f is then selected so that the analytical solution to (1) is u(x, y, z) = 1 + sin(πxy2z/8). (a) (b)

Figure: 3D example: Agglomerated meshes. (a) 64 elements; (b) 32768 elements.

21 / 29

3D Numerical example

dof1/3 101 102 ||u-uh||L2(Ω) 10-12 10-10 10-8 10-6 10-4 10-2 p=1 p=2 p=3 p=4

DGFEM DGFEM(P)

dof1/3

101 102

|||u-uh|||DG

10-10 10-8 10-6 10-4 10-2

p=1 p=2 p=3 p=4

DGFEM DGFEM(P)

Figure: Convergence of DGFEM under h–refinement for p = 1, 2, 3, 4.

22 / 29

3D Numerical example

dof1/3 10 20 30 40 50 60 70 ||u-uh||L2(Ω) 10-10 10-8 10-6 10-4 10-2 64 Elements 512 Elements 4096 Elements

DGFEM DGFEM(P) DGFEM(Q)

dof1/3

10 20 30 40 50 60 70

|||u-uh|||DG

10-10 10-8 10-6 10-4 10-2

64 Elements 512 Elements 4096 Elements

DGFEM DGFEM(P) DGFEM(Q)

Figure: Convergence of DGFEM under p–refinement.

23 / 29

Numerical example – parabolic problem

Let Ω (0, 1)2, t ∈ J = (0, 1) and we consider following heat equation: ut − ∆u = f . f is then selected so that the analytical solution to (1) is u(x, y, t) = sin(20πt)e−5((x−0.5)2+(y−0.5)2), Numerical Experiments:

DG(P)-Pp basis over 3D space-time elements (poly/rect spatial meshes) DG(PQ)-Pp basis over 2D rect spatial meshes tensor product temporal basis DG(Q)-Qp basis over 3D space-time elements with rect spatial meshes FEM(Q)-Qp basis over 3D space-time elements with rect spatial meshes

Red denotes space-time DG schemes, while blue denotes DG time-stepping schemes with FEM in space.

24 / 29

Numerical example – parabolic problem

10

1

10

2

10

3

10

−10

10

−9

10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10

Dof1/3 || u−uh||L2((0,T),H1(Ω)) L2((0,T),H1(Ω)) norm error under h refinement

DG rect P2 slope 2.0118 DG rect P5 slope 5.143 DG rect P6 slope 6.2208 DG poly P2 slope 1.9935 DG poly P5 slope 5.0011 DG poly P6 slope 6.02102 DG rect PQ1 slope 1.0077 DG rect PQ4 slope 4.0146 DG rect Q1 slope 1.0059 DG rect Q4 slope 4.0006 FEM rect Q1 slope 1.0363 FEM rect Q4 slope 4.024

10

1

10

2

10

3

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10

Dof1/3 || u−uh||L∞((0,T),L2(Ω)) L∞((0,T),L2(Ω)) norm error under h refinement

DG rect P2 slope 2.8031 DG rect P5 slope 5.6289 DG rect P6 slope 6.6693 DG poly P2 slope 2.8084 DG poly P5 slope 5.5871 DG rect P6 slope 6.6630 DG rect PQ1 slope 1.9968 DG rect PQ4 slope 4.9904 DG rect Q1 slope 1.9864 DG rect Q4 slope 4.9903 FEM rect Q1 slope 2.021 FEM rect Q4 slope 5.0084

Figure: space-time h–refinement.

25 / 29

Numerical example – parabolic problem

20 40 60 80 100 120 140 10

−10

10

−8

10

−6

10

−4

10

−2

10 Dof1/3 || u−uh||L2((0,T),H1(Ω)) p refinement 64 rect DG(P) 80 time steps 64 poly DG(P) 80 time steps 64 rect DG(PQ) 80 time steps 64 rect DG(Q) 80 time steps 64 rect FEM(Q) 80 time steps 20 40 60 80 100 120 140 10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 Dof1/3 || u−uh||L∞((0,T),L2(Ω)) p refinement 64 rect DG(P) 80 time steps 64 poly DG(P) 80 time steps 64 rect DG(PQ) 80 time steps 64 rect DG(Q) 80 time steps 64 rect FEM(Q) 80 time steps

Figure: p–refinement, T = 1.

26 / 29

Numerical example – parabolic problem

50 100 150 200 250 300 350 400 10

−10

10

−8

10

−6

10

−4

10

−2

10 10

2

Dof1/3 || u−uh||L2((0,T),H1(Ω)) p refinement

64 rect DG(P) 3200 time steps 64 rect FEM(Q) 3200 time steps 64 poly DG(P) 3200 time steps 50 100 150 200 250 300 350 400 10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10

Dof1/3 || u−uh||L∞((0,T),L2(Ω)) p refinement

64 rect DG(P) 3200 time steps 64 rect FEM(Q) 3200 time steps 64 poly DG(P) 3200 time steps

Figure: p–refinement, T = 40.

27 / 29

Current Work

A-posteriori error estimates for hp-version DG on general polytopic meshes with degenerating faces. hp-adaptive DG method on agglomeration meshes for evolution problems. Explore relationship between Dof against error for hp-version DG method with different basis.

28 / 29

• A. Cangiani, E. H. Georgoulis, and P. Houston.

hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. Mathematical Models and Methods in Applied Sciences, 24(10):2009–2041, 2014.

• A. Cangiani, Z.Dong, E. H. Georgoulis, and P. Houston.

hp-version discontinuous Galerkin methods for advection-diffusion-reaction problems on polytopic meshes. M2AN, Published online 2015.

• A. Cangiani, Z.Dong and E. H. Georgoulis.

hp-version space time discontinuous Galerkin methods for parabolic problems on polytopic meshes. In preparation.

• A. Cangiani, Z.Dong, E. H. Georgoulis, and P. Houston.

hp-version discontinuous Galerkin methods on polygonal and polyhedral meshes. Springer, in preparation.

29 / 29