[PPT] - A posteriori analysis of discontinuous Galerkin schemes for systems PowerPoint Presentation

SLIDE 1

A posteriori analysis of discontinuous Galerkin schemes for systems of hyperbolic conservation laws

Jan Giesselmann

joint work with Ch. Makridakis (Univ. of Sussex), T. Pryer (Univ. of Reading) Sino-German Symposium on Modern Numerical Methods for Compressible Fluid Flows and Related Problems

May 21 - 27 2014

1 / 22

SLIDE 2

Outline

1. Introduction
2. Relative entropy
3. Semi-discrete DG schemes
4. Reconstruction approach
5. Error estimate
6. Numerical examples

2 / 22

SLIDE 3

Introduction

We study the approximation of (systems of) hyperbolic conservation laws in one space dimension ut + f(u)x = 0 in (0, ∞) × I, (CLAW) where u : (0, ∞) × I → U ⊂ Rn and f ∈ C2(U, Rn), by semi-discrete discontinuous Galerkin schemes. Initial conditions: u(0, ·) = u0. Periodic boundary conditions.

3 / 22

SLIDE 4

Introduction

We study the approximation of (systems of) hyperbolic conservation laws in one space dimension ut + f(u)x = 0 in (0, ∞) × I, (CLAW) where u : (0, ∞) × I → U ⊂ Rn and f ∈ C2(U, Rn), by semi-discrete discontinuous Galerkin schemes. Initial conditions: u(0, ·) = u0. Periodic boundary conditions. Entropy/ entropy flux pair: η : U → R, q : U → R such that ∇ηDf = ∇q. (∗) Definition: Entropy solution A weak solution of (CLAW) is called an entropy solution with respect to the entropy/ entropy flux pair (η, q) provided it satisfies η(u)t + q(u)x ≤ 0 in the sense of distributions.

3 / 22

SLIDE 5

Known results

A posteriori error estimates: Ohlberger, Kröner ’00: Upwind finite volume schemes for scalar conservation laws. Hartmann, Houston ’02: Space-time DG schemes for systems of conservation laws. Jovanovic, Rohde ’05: Finite volume schemes for Friedrichs systems. Dedner, Ohlberger, Makridakis ’07: Runge-Kutta DG schemes for scalar conservation laws. A priori error estimates for DG schemes: Zhang, Shu ’04, ’10: For RKDG2, RKDG3 schemes the spatial error is O(hp+ 1

2 ). It is O(hp+1) for upwind flux.

4 / 22

SLIDE 6

Aim

Construct numerical solution uh

t + f(uh)x = Rh

with residual Rh. Construct explicit estimator E = E (uh) such that u − uhL∞(0,T;L2(I)) ≤ E (uh). E (uh) should be of the same order as the error.

5 / 22

SLIDE 7

Remarks on entropies

In the scalar case every function η : U → R is an entropy. In the scalar case entropy solutions (satisfying the entropy inequality for all convex entropies) are unique. For systems of hyperbolic conservation laws there is usually only one (physically motivated) entropy/entropy flux pair. (For systems) entropy solutions need not be unique. In most cases the entropy is convex. But there are important cases where it is not. In the sequel we consider a system endowed with a strictly convex entropy.

6 / 22

SLIDE 8

Relative entropy

In which way does the entropy inequality ensure stability? Relative entropy & relative entropy flux; (Dafermos ’79, Di Perna ’79): Let (η, q) be an entropy/entropy flux pair. η(u|v) := η(u) − η(v) − ∇η(v)(u − v) q(u|v) := q(u) − q(v) − ∇η(v)(f(u) − f(v)). For K ⊂ U convex and compact, η strictly convex, there exists c > 0 such that η(u|v) ≥ c|u − v|2 ∀ u, v ∈ K.

7 / 22

SLIDE 9

Relative entropy

Theorem (Dafermos): Let a system of conservation laws be endowed with a strictly convex entropy/entropy flux pair (η, q). Let u be an entropy solution and v a Lipschitz solution. Assume that u and v take values in some convex and compact K ⊂ U. Then there exist a, b > 0 such that for t > 0 u(t, ·) − v(t, ·)L2(I) ≤ aebtu0 − v0L2(I). b depends on the Lipschitz constant of v. Remark 1: Uniqueness of Lipschitz solutions in the class of entropy solutions. Remark 2: In the scalar case the abundance of entropies leads to L1-contraction: u(t, ·) − v(t, ·)L1(I) ≤ u(s, ·) − v(s, ·)L1(I) for t > s for entropy solutions u, v.

8 / 22

SLIDE 10

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

9 / 22

SLIDE 11

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

9 / 22

SLIDE 12

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

η(u|v)t + q(u|v)x ≤ −vT

t D2η(v)(u − v) − vT x D2η(v)(f(u) − f(v))

9 / 22

SLIDE 13

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

η(u|v)t + q(u|v)x ≤ −vT

t D2η(v)(u − v) − vT x D2η(v)(f(u) − f(v))

η(u|v)t + q(u|v)x ≤ (Df(v)vx)TD2η(v)(u − v) − vT

x D2η(v)(f(u) − f(v))

9 / 22

SLIDE 14

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

η(u|v)t + q(u|v)x ≤ −vT

t D2η(v)(u − v) − vT x D2η(v)(f(u) − f(v))

η(u|v)t + q(u|v)x ≤ (Df(v)vx)TD2η(v)(u − v) − vT

x D2η(v)(f(u) − f(v))

Existence of an entropy flux implies (Df)TD2η = D2ηDf

9 / 22

SLIDE 15

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

η(u|v)t + q(u|v)x ≤ −vT

t D2η(v)(u − v) − vT x D2η(v)(f(u) − f(v))

η(u|v)t + q(u|v)x ≤ (Df(v)vx)TD2η(v)(u − v) − vT

x D2η(v)(f(u) − f(v))

Existence of an entropy flux implies (Df)TD2η = D2ηDf η(u|v)t + q(u|v)x = −vT

x D2η(v)

f(u) − f(v) − Df(v)(u − v)
9 / 22

SLIDE 16

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

η(u|v)t + q(u|v)x ≤ −vT

t D2η(v)(u − v) − vT x D2η(v)(f(u) − f(v))

η(u|v)t + q(u|v)x ≤ (Df(v)vx)TD2η(v)(u − v) − vT

x D2η(v)(f(u) − f(v))

Existence of an entropy flux implies (Df)TD2η = D2ηDf η(u|v)t + q(u|v)x = −vT

x D2η(v)

f(u) − f(v) − Df(v)(u − v)
d

d t

I η(u|v) d x ≤ cu − v2

L2 ≤ b

I η(u|v) d x

9 / 22

SLIDE 17

Idea of the proof (formal)

For the rigorous proof, see Dafermos’ book. η(u|v)t = η(u)t − vT

t D2η(v)(u − v) − ∇η(v)ut

q(u|v)x = q(u)x − vT

x D2η(v)(f(u) − f(v)) − ∇η(v)f(u)x

η(u|v)t + q(u|v)x ≤ −vT

t D2η(v)(u − v) − vT x D2η(v)(f(u) − f(v))

η(u|v)t + q(u|v)x ≤ (Df(v)vx)TD2η(v)(u − v) − vT

x D2η(v)(f(u) − f(v))

Existence of an entropy flux implies (Df)TD2η = D2ηDf η(u|v)t + q(u|v)x = −vT

x D2η(v)

f(u) − f(v) − Df(v)(u − v)
d

d t

I η(u|v) d x ≤ cu − v2

L2 ≤ b

I η(u|v) d x

cu(t, ·) − v(t, ·)2

L2 ≤

I η(u(t, ·)|v(t, ·)) d x ≤
I η(u0|v0) d x ebt

9 / 22

SLIDE 18

The case with residuals

Let v solve a perturbed problem: ut + f(u)x = 0, vt + f(v)x = R. Then,

I

η(u|v)t d x =

I

−

f(u) − f(v) − Df(v)(u − v)

TD2η(v)vx − RTD2η(v)(u − v) d x, such that

I

η(u|v)t d x ≤ bu − v2

L2(I) + R2 L2(I).

10 / 22

SLIDE 19

A first estimate

Proposition Let u be an entropy solution of (CLAW) and v a Lipschitz solution of a perturbed problem with residual R. Let u, v take values in a compact, convex K ⊂ U. Then there exist a, b > 0 such that for t > 0 u(t, ·) − v(t, ·)2

L2(I) ≤ a

u0 − v02

L2(I) + R2 L2([0,t)×I)

ebt.

b is related to the Lipschitz constant of v.

11 / 22

SLIDE 20

Semi-discrete DG scheme

For I = (0, 1) choose partition 0 = x0 < x1 < · · · < xN = 1. Vp space of piecewise polynomials of degree ≤ p. Find uh ∈ C 1((0, ∞), Vd

p) such that

I

uh

t · φ − f(u) · φx d x + N−1

n=0

F(uh(x−

n ), uh(x+ n )) · [φ]n = 0

for all φ ∈ Vd

p with [φ]n := φ(x− n ) − φ(x+ n ).

F : U × U → Rn is the numerical flux function. We impose that there exist L > 0 and w : U × U → U such that F(a, b) = f(w(a, b)) |w(a, b) − a| + |w(a, b) − b| ≤ L|a − b| ∀a, b ∈ U.

12 / 22

SLIDE 21

Condition on the numerical flux

For Burger’s equation our condition is not satisfied for the (local and global) Lax-Friedrichs flux. For Burger’s equation it is satisfied for the Roe flux with w(a, b) =

1

2a2(1 + sgn(a − b)) + 1 2b2(1 − sgn(a − b)) It is satisfied for exact Riemann solvers. -> Godunov scheme. It is satisfied for approximate Riemann solvers computing the flux by evaluating the flux function on a state of the approximate Riemann solution.

13 / 22

SLIDE 22

The need for reconstruction

In employing the relative entropy framework for estimating the difference between exact and approximate solution there are two key challenges: uh(t, ·), u(t, ·) are not Lipschitz continuous. uh

t + f(uh)x =: Rh is just measure valued, Rh(t, ·) ∈ L2(I).

14 / 22

SLIDE 23

The need for reconstruction

In employing the relative entropy framework for estimating the difference between exact and approximate solution there are two key challenges: uh(t, ·), u(t, ·) are not Lipschitz continuous. uh

t + f(uh)x =: Rh is just measure valued, Rh(t, ·) ∈ L2(I).

Basic idea of reconstruction: Introduce intermediate quantity ^ u such that ^ u − uh can be explicitly bounded and is of optimal order. ^ u is Lipschitz & satisfies a perturbed PDE with an L2 residual. Use u − uh ≤ u − ^ u + ^ u − uh.

14 / 22

SLIDE 24

The reconstructions

We employ a reconstruction technique used for DG-time-discretisations by Makridakis & Nochetto ’06. Let In denote the n-th sub-interval In := [xn, xn+1].

15 / 22

SLIDE 25

The reconstructions

We employ a reconstruction technique used for DG-time-discretisations by Makridakis & Nochetto ’06. Let In denote the n-th sub-interval In := [xn, xn+1]. We define ^ u ∈ C 1((0, ∞), Vd

p+1) by

I

(^ u − uh) · φ d x = 0 for all φ ∈ Vd

p−1

^ u(x±

n ) = w(uh(x− n ), uh(x+ n )).

15 / 22

SLIDE 26

The reconstructions

We employ a reconstruction technique used for DG-time-discretisations by Makridakis & Nochetto ’06. Let In denote the n-th sub-interval In := [xn, xn+1]. We define ^ u ∈ C 1((0, ∞), Vd

p+1) by

I

(^ u − uh) · φ d x = 0 for all φ ∈ Vd

p−1

^ u(x±

n ) = w(uh(x− n ), uh(x+ n )).

We reconstruct the flux by ^ f ∈ C 0((0, ∞), Vd

p+1) such that

I

^ fx · φ d x = −

I

f(uh) · φx d x +

N−1

n=0

F(uh(x−

n ), uh(x+ n )) · [φ]n

for all φ ∈ Vd

p and

^ f(x−

n ) = F(uh(x− n ), uh(x+ n )).

15 / 22

SLIDE 27

Properties of the reconstructions

^ u,^ f are uniquely determined and continuous. Indeed, they are Lipschitz continuous for every t. ^ u,^ f can be computed explicitly and locally. By construction:

I

uh

t · φ +^

fx · φ d x = 0 for all φ ∈ Vd

p.

This implies uh

t +^

fx = 0 a.e. Thus, ^ ut + f(^ u)x = (^ u − uh)t + (f(^ u) −^ f)x =: Rh such that Rh(t, ·) ∈ L2(I) for every t. = ⇒ We can use the relative entropy framework to bound u − ^ uL2(I).

16 / 22

SLIDE 28

Representation of the reconstruction ^ u

Let ln

p and ln p+1 be Legendre polynomials of degree p, p + 1 rescaled to

[xn, xn+1] then for x ∈ [xn, xn+1] 2(^ u − uh)(x) =

(−1)p(w(uh(x−

n ), uh(x+ n )) − uh(x+ n ))

+ w(uh(x−

n+1), uh(x+ n+1)) − uh(x− n+1)

ln

p (x)

+

(−1)p+1(w(uh(x−

n ), uh(x+ n )) − uh(x+ n ))

+ w(uh(x−

n+1), uh(x+ n+1)) − uh(x− n+1)

ln

p+1(x)

such that ^ u − uh2

L2(I) ≤ L2 n

(xn+1 − xn)

uh

n

2 +
uh

n+1

2

. Arguments of Makridakis & Nochetto ’06 indicate that the right hand side is of optimal order. Analogously we get a bound for ^ ut − uh

t 2 L2(I).

17 / 22

SLIDE 29

Estimate for (^ f − f(^ u))x

More involved, as the x derivative might cause loss of one order. Here it is crucial that f(^ u(xn)) = ^ f(xn) for all n. This is the reason for our assumption on the numerical flux. (^ f − f(^ u))xL2 ≤ C

N−1

n=0

hn(an + bn) with an :=

uh

n

2 +
uh

n+1

2
uh

n

+
uh

n+1

hn

+

∂xuh
L∞(In)
bn :=
p
k=0

p + 1 k

hp+1

n

∂k+1

x

uh

L∞(In) + hp−k

n

uh

n

+
uh

n+1

×
∂p+1−k

x

Df(uh)

2

.

18 / 22

SLIDE 30

Error estimate

Theorem (JG, Makridakis, Pryer; ’14) Let u be an entropy solution of (CLAW) and uh a solution of a semi-discrete dG scheme with reconstruction ^ u. Let u, ^ u take values in a compact, convex K ⊂ U. Then there exist a, b > 0 such that for t > 0 u(t, ·) − uh(t, ·)2

L2(I) ≤ ^

u(t, ·) − uh(t, ·)2

L2(I)

+ a

u0 − ^

u02

L2(I) + Rh2 L2([0,t)×I)

ebt.

(∗) ^ u(t, ·) − uh(t, ·)L2(I), Rh, a, b are explicitly computable. b depends on the Lipschitz constant of ^

u. It is expected to be

bounded for h → 0, if u is Lipschitz. We expect the right hand side of (∗) to be of the order of the error, if u is Lipschitz. We expect the right hand side of (∗) to blow up for h → 0, if u is discontinuous.

19 / 22

SLIDE 31

Numerical example

For Burgers equation on the interval (−π, π) with initial data u0(x) = sin (x) the exact solution (pre-shock) is given by u(x, t) = −2

∞

k=1

Jk(kt) kt sin (kx) , where Jk is the k − th Bessel function. Instead of computing the full estimator we compute

n

(xn+1 − xn)

uh

n

2 +
uh

n+1

2

, which is of the same order. We used the Engquist-Osher flux, which satisfies our condition on the numerical flux. Time discretisation is RK4 with τ = h/4.

20 / 22

SLIDE 32

Numerical results

Numerical results for p = 1; EOC: Experimental order of convergence, EI: efficiency index N

u − uh
L∞(L2)

EOC E L∞(L2) EOC EI 8 2.3336e-01 0.000 3.5500e-01 0.000 1.521 16 8.6657e-02 1.429 1.3541e-01 1.390 1.563 32 3.1863e-02 1.443 5.1422e-02 1.397 1.614 64 1.1753e-02 1.439 1.9416e-02 1.405 1.652 128 4.2916e-03 1.453 7.1950e-03 1.432 1.677 256 1.5501e-03 1.469 2.6220e-03 1.456 1.692 512 5.5526e-04 1.481 9.4403e-04 1.474 1.700 1024 1.9779e-04 1.489 3.3723e-04 1.485 1.705 2048 7.0216e-05 1.494 1.1990e-04 1.492 1.708 4096 2.4879e-05 1.497 4.2518e-05 1.496 1.709

21 / 22

SLIDE 33

Numerical results

Numerical results for p = 2 EOC: Experimental order of convergence, EI: efficiency index N

u − uh
L∞(L2)

EOC E L∞(L2) EOC EI 8 2.2135e-02 0.000 9.5664e-02 0.000 0.000 16 2.7472e-03 3.010 1.4455e-02 2.726 5.262 32 3.4615e-04 2.988 2.1409e-03 2.755 6.185 64 4.4617e-05 2.956 3.3881e-04 2.660 7.594 128 5.6079e-06 2.992 5.1629e-05 2.714 9.207 256 7.0465e-07 2.992 7.4392e-06 2.795 10.557 512 8.8207e-08 2.998 1.0230e-06 2.862 11.598 1024 1.1040e-08 2.998 1.3598e-07 2.911 12.317

21 / 22

SLIDE 34

Numerical results

Numerical results for p = 3 EOC: Experimental order of convergence, EI: efficiency index N

u − uh
L∞(L2)

EOC E L∞(L2) EOC EI 8 2.4639e-03 0.000 1.2564e-02 0.000 0.000 16 2.4557e-04 3.327 2.6343e-03 2.254 10.727 32 2.2211e-05 3.467 3.3904e-04 2.958 15.264 64 1.9772e-06 3.490 3.4219e-05 3.309 17.307 128 1.7515e-07 3.497 3.2659e-06 3.389 18.647 256 1.5491e-08 3.499 3.0403e-07 3.425 19.626

21 / 22

SLIDE 35

Summary and outlook

Summary: We have derived an a posteriori error estimate for DG schemes approximating hyperbolic conservation laws with convex entropies. In case the solution is at least Lipschitz the estimator is expected to be of optimal order. The estimate also holds for discontinuous solutions, but in that case the estimator blows up for h → 0.

22 / 22

SLIDE 36

Summary and outlook

Summary: We have derived an a posteriori error estimate for DG schemes approximating hyperbolic conservation laws with convex entropies. In case the solution is at least Lipschitz the estimator is expected to be of optimal order. The estimate also holds for discontinuous solutions, but in that case the estimator blows up for h → 0. Outlook: Use known results on Runge-Kutta methods to obtain estimates for fully discrete RKDG schemes. Extension to several space dimensions (-> only reconstruction). Extension to higher order problems (e.g. NSK). Extension to non-convex entropies -> involution. Include e.g. flux limiters.

22 / 22