Space-Time Discontinuous Galerkin Discretizations for Linear - - PowerPoint PPT Presentation

www.kit.edu

Space-Time Discontinuous Galerkin Discretizations for Linear First-Order Hyperbolic Evolution Systems

Willy Dörfler | Stefan Findeisen | Christian Wieners

Institut für Angewandte und Numerische Mathematik

Outline

Linear hyperbolic first-order systems M∂tu + Au = f in [0, T], u(0) = u0 well-posedness in a space-time Hilbert space setting Discontinuous Galerkin method in space and Petrov-Galerkin method in time fully implicit adaptive space-time method for Lu = f with L = M∂t + A inf-sup stability of a discrete space-time Petrov-Galerkin setting Error Control duality based space-time error representation weighted residual error indicator A space-time multigrid strategy combined p−q−h subspace correction preconditioner smoothing in space and in time Applications linear transport equation Maxwell’s equations for polarized waves Acoustic waves (joint work with Daniel Ziegler)

Dörfler / Findeisen / Wieners: Space-time discontinuous Galerkin discretizations for linear first-order hyperbolic evolution systems. Comput. Methods Appl. Math. 2016. We gratefully acknowledge financial support by DFG through RTG 1294 and CRC 1173.

1

Linear Hyperbolic First-order Systems

We consider in the bounded domain Ω ⊂ RD and the time interval [0, T] M∂tu(t) + Au(t) = f(t) , t ∈ (0, T) , u(0) = u0 and a flux function F(v) = [B1v, . . . , BDv] with Bj ∈ L∞(Ω)J×J

sym such that

Av = div F(v) =

D

• d=1

∂d(Bdv) ∈ L2(Ω)J , v ∈ D(A) . Linear transport For ρ ∈ L∞(Ω) and q ∈ W1

∞(Ω)D with div q = 0 find u such that

ρ∂tu + div(uq) = f . Acoustic waves Find the pressure p and the velocity v (isotropic and homogeneous media) with ∂tp + div v = 0 , ∂tv + ∇p = 0 . Electro-magnetic waves Find an electric field E and magnetic field H such that ε∂tE − curl H = f , µ∂tH + curl E = 0 , div(εE) = ρ , div(µH) = 0 for given permeability µ and permittivity ε.

2

A Space-Time Setting

Let H ⊆ L2(Ω)J be a Hilbert space with (v, w)H = (Mv, w)0,Ω. Let D(A) ⊂ H be the domain of the operator A with (Av, v)0,Ω ≥ 0 in D(A). Consider L = M∂t + A on the space-time cylinder Q = Ω × (0, T). Let V = D(L) be the closure of

• u ∈ C1(0, T; D(A)): u(0) = 0
• with respect

to the weighted graph norm v2

V = (Mv, v)0,Q + (M−1Lv, Lv)0,Q.

Let W be the closure of L(V) in L2(0, T; H) with u2

W = (Mu, u)2 0,Q.

Lemma

For given f ∈ L2(Q)J a unique solution u ∈ V exists solving the variational problem (Lu, w)0,Q = (f, w)0,Q , w ∈ W For the proof we check inf-sup stability of b: V × W → R with b(v, w) = (Lv, w)0,Q. We observe vW ≤ 2T M−1LvW for all v ∈ V. This yields inf

v∈V sup w∈W

b(v, w) vVwW ≥ inf

v∈V

b(v, M−1Lv) vVM−1LvW = inf

v∈V

M−1LvW

• v2

W + M−1Lv2 W

≥ 1 √ 1 + 4T 2 .

3

A Petrov–Galerkin Space-Time Discretization

Let Q =

R∈R R be a decomposition of the space-time cylinder into space-time

cells R = K × I with K ⊂ Ω and I = (t−, t+) ⊂ (0, T). Choose locally Vh,R, Wh,R ⊂ L2(R)J with Wh,R ⊂ ∂tVh,R and set Vh =

• vh ∈ H1(0, T; H): vh(x, 0) = 0 for a.a. x ∈ Ω and vh,R = vh|R ∈ Vh,R
• ,

Wh =

• wh ∈ L2(0, T; H): wh,R = wh|R ∈ Wh,R
• .

Tensor-product space-time meshes Consider a fixed mesh K in space and a time series 0 = t0 < t1 < · · · < tN = T, i.e., R =

K∈K

N

n=1 K × (tn−1, tn). Then, set Wh,R =

• Pp(K) ⊗ Pq−1

J and Vh =

• vh ∈ H1(0, T; H) : vh(x, 0) = 0 and for (x, t) ∈ R = K × (tn−1, tn)

vh(x, t) = tn − t tn − tn−1 vh(x, tn−1) + t − tn−1 tn − tn−1 wR,h(x, t) , wh ∈ Wh,R

• .

General space-time meshes Set Wh,R =

• PpR(K) ⊗ PqR−1

J, and define Vh,R =

• vh,R ∈ L2(R)J : vh,R(x, t) = t+ − t

t+ − t− vh(x, t−) + t − t− t+ − t− wh,R(x, t) , vh ∈ Vh|[0,t−] , wh,R ∈ Wh,R , (x, t) ∈ R = K × (t−, t+)

• .

4

A Discontinuous Galerkin Approximation

We define the discontinuous Galerkin operator Ah ∈ L(Vh, Wh) by

• Ahvh, wh)0,Q =
• R=K×I
• div F(vh,R), wh,R
• 0,R

+

• f∈FK
• nK · (Fnum

K

(vh) − F(vh,R)), wh,R

• 0,f×I
• using the upwind flux Fnum

K

(vh) on every face f ⊂ ∂K. Let Πh : W → Wh be the L2-projection, and define dT(t) = T − t.

Lemma

Let Lh = M∂t + Ah and f ∈ L2(Q)J, and assume that

• M∂tvh, dTvh
• 0,Q ≤
• Lhvh, dTΠhvh
• 0,Q ,

vh ∈ Vh . (1) Then, a unique discrete solution uh ∈ Vh exists solving (Lhuh, vh)0,Q = (f, vh)0,Q , vh ∈ Wh . The proof relies on discrete inf-sup stability and the identity for vh ∈ Vh vh2

W = −

T

• Mvh(t), vh(t)
• 0,Ω∂tdT(t) dt = 2
• M∂tvh, dTvh
• 0,Q .

5

A Discontinuous Galerkin Approximation

Theorem

Assume that (1) holds and set vh2

Vh = vh2 W + M−1 h Lhvh2

• W. Then, we have

u − uhVh ≤

• 1 +
• 4T 2 + 1
• inf

vh∈Vh

u − vhVh . If in addition the solution is sufficiently smooth, we obtain the a priori error estimate u − uhVh ≤ C

• △tq + △xp

∂q+1

t

u0,Q + Dp+1u0,Q

• for △t, △x and p, q ≥ 1 with △t ≥ t+ − t−, △x ≥ diam K, p ≤ pR and q ≤ qR.

Lemma

In the case of tensor product space-time discretizations, the condition (1) is satisfied, and we have for vh ∈ Vh Πh∂tvh = ∂tvh ,

• Mh∂tvh, dTvh
• 0,Q ≤
• Mh∂tvh, dTΠhvh
• 0,Q , 0 ≤
• Ahvh, dTΠhvh
• 0,Q .

The proof relies on

• t∂tλk, λk
• 0 = k for orthonormal Legendre polynomials λk.

6

A Dual Space-Time Error Representation

Let L∗ = −L be the adjoint operator defined on the adjoint Hilbert space V ∗ =

• w ∈ W : there exists g ∈ W such that (Lv, w)0,Q = (v, g)0,Q for all v ∈ V
• .

Note that we have w(T) = 0 for w ∈ V ∗. For a given error functional E ∈ W ′ we introduce the dual solution u∗ ∈ V ∗ with (w, L∗u∗)0,Q = E, w , w ∈ W .

Lemma

If the dual solution is sufficiently smooth, we have for all wh ∈ Wh E, u − uh =

• R=K×I∈R
• f − M∂tuh − div F(uh), u∗ − wh
• 0,R

+

• nK · (F(uh) − Fnum

K

(uh)), u∗ − wh

• 0,∂K×I
• .

Inserting the discrete dual solution u∗

h ∈ Wh with (vh, L∗ hu∗ h)0,Q = E, vh for vh ∈ Vh

and a higher-order recovery Ihu∗

h defines the weighted residual error indicator

ηR = f − M∂tuh − div F(uh)0,Ru∗

h − Ihu∗ h0,R

+ nK · (F(uh) − Fnum(uh))0,∂K×Iu∗

h − Ihu∗ h0,∂K×I .

7

A Multigrid Preconditioner

Let R0,0 be the coarse space-time mesh, and let Rl,k be the mesh obtained by l = 1, . . . , lmax refinements in space and k = 1, . . . , kmax refinements in time.

BGS

l,k

Rl,k

l−1,k

BGS

l−1,k

BJ

0,k

R0,k

0,k−1

BJ

0,k−1

BML

0,0

BJ

0,k−1

P0,k

0,k−1

BJ

0,k

BGS

l−1,k

Pl,k

l−1,k

BGS

l,k

(l, k) (l − 1, k) (0, k) (0, k − 1) (0, 0)

Ll,k ∈ L(Vl,k, Wl,k) approximates L on Rl,k the prolongation Pl,k

l−1,k ∈ L(Vl−1,k, Vl,k) represents the natural injection

the restriction Rl,k

l−1,k ∈ L(Wl,k, Wl−1,k) represents the L2 projection

BJ

l,k = θl,k block_diag(Ll,k)−1 is the block-Jacobi smoother with damping θl,k

BGS

l,k = θl,k

• block_lower(Ll,k) + block_diag(Ll,k)

−1 Gauss-Seidel smoother

8

Numerical Test: Rotating Cone

∂tu + div(uq) = 0 in Ω = (−0.5, 0.5)2 for t ∈ (0, 1) with q(x1, x2) = 2π(−x2, x1)

9

Numerical Test: Rotating Cone

k = 1 k = 2 k = 3 k = 4 k = 5 k = 6 degrees of freedom on the space-time mesh Rl,k with polynomial degrees (p, q) = (2, 2) l = 1 768 × 16 768 × 32 768 × 64 768 × 128 768 × 256 768 × 512 l = 2 3 072 × 16 3 072 × 32 3 072 × 64 3 072 × 128 3 072 × 256 3 072 × 512 l = 3 12 288 × 16 12 288 × 32 12 288 × 64 12 288 × 128 12 288 × 256 12 288 × 512 l = 4 49 152 × 16 49 152 × 32 49 152 × 64 49 152 × 128 49 152 × 256 49 152 × 512 two-level iteration in time with Jacobi smoothing (νl,k = 2, θl,k = 0.5) l = 1 26 (4.97e-1) 24 (4.65e-1) 9 (1.31e-1) 7 (5.87e-2) 7 (5.81e-2) 7 (5.77e-2) l = 2 38 (6.09e-1) 32 (5.65e-1) 7 (6.88e-2) 7 (5.22e-2) 7 (5.30e-2) l = 3 57 (7.19e-1) 48 (6.83e-1) 6 (4.52e-2) 6 (4.30e-2) l = 4 106 (8.40e-1) 94 (8.20e-1) 6 (3.51e-2) two-level iteration in space with Gauss-Seidel smoothing (νl,k = 5) l = 1 4 (1.31e-4) 4 (1.36e-4) 4 (1.85e-4) 4 (1.79e-4) 4 (1.68e-4) 4 (1.68e-4) l = 2 5 (1.35e-2) 5 (5.50e-3) 5 (4.72e-3) 5 (5.50e-3) 5 (5.28e-3) 5 (4.94e-3) l = 3 8 (7.57e-2) 7 (3.63e-2) 7 (2.26e-2) 6 (4.24e-2) 6 (4.07e-2) 6 (3.89e-2) l = 4 15 (2.70e-1) 11 (1.61e-1) 10 (1.34e-1) 9 (1.27e-1) 9 (1.24e-1) 9 (1.20e-1) full space-time multilevel method l = 1 4 (1.25e-4) 4 (1.26e-4) 4 (1.80e-4) 4 (1.92e-4) 4 (1.86e-4) 4 (1.75e-4) l = 2 5 (1.35e-2) 5 (5.50e-3) 5 (4.71e-3) 5 (5.50e-3) 5 (5.28e-3) 5 (4.94e-3) l = 3 8 (7.57e-2) 7 (3.63e-2) 7 (2.25e-2) 6 (4.24e-2) 6 (4.07e-2) 6 (3.89e-2) l = 4 15 (2.73e-1) 11 (1.61e-1) 10 (1.34e-1) 9 (1.27e-1) 9 (1.24e-1) 9 (1.20e-1)

iteration steps and averaged rates for a residual reduction by the factor 10−8

10

Numerical Test: Rotating Cone

The adaptive solution with 3 303 810 degrees of freedom computes the same solution as the uniform computation on 524 288 = 4 096 × 128 space-time cells and 31 703 040 degrees of freedom. solution sliced at times t = 0, 0.3, 0.6, 1 location of the highest polynomial degrees

11

Numerical Test: A Double-slit Configuration

We consider an interfering wave in a double-slit experiment with reflecting boundary conditions constant material parameters µ = 1, ε = 1 wavelength λ = 0.5 slit gap a = 1.25 and slit width b = 0.25 the error functional E(e) := |S|−1

• S

ey d(x, t) with region of interest S = (5.5, 6) × (0, 2) × (0, 8) analytic diffraction pattern Idiff(α) = I0sinc2π λb sin α

• cos2 π

λa sin α

• in dependence of the observation angle α

triangular mesh in domain Ω with region of interest S

12

Numerical Test: A Double-slit Configuration

time evolution of an interfering wave with variable polynomial degrees pR

13

Numerical Test: Adaptivity

uniform mesh adaptive mesh polynomial degree pR

14

Numerical Test: Parallel Performance on ForHLR

uniform refinement adaptive refinement level (p, q) #DoFs steps (rate) |El − Eex| #DoFs (effort) steps (rate) |El − Eex| l = 1 (1, 1) 5 477 184 10 (0.11) 3.1303e-1 5 477 184 10 (0.11) 3.1303e-1 l = 2 (2, 2) 21 908 736 17 (0.28) 2.6220e-2 9 645 930 (44%) 13 (0.23) 2.6230e-2 l = 3 (3, 3) 54 771 840 24 (0.45) 1.4502e-3 17 309 043 (32%) 18 (0.34) 1.4502e-3 l = 4 (4, 4) 109 543 680 34 (0.57) 8.0209e-5 29 064 348 (27%) 23 (0.44) 8.0209e-5

Uniform vs. adaptive refinement

• n 606 208 = 2 368 × 256 space-time cells

distributed to 256 processes and

• n 4 849 664 = 9 472 × 512 space-time cells

distributed to 1024 processes, respectively.

uniform refinement adaptive refinement level (p, q) #DoFs steps (rate) |El − Eex| #DoFs (effort) steps (rate) |El − Eex| l = 1 (1, 1) 43 732 224 17 (0.31) 9.7373e-2 43 732 224 17 (0.31) 9.7373e-2 l = 2 (2, 2) 174 928 896 31 (0.54) 1.7630e-3 68 437 899 (39%) 21 (0.40) 1.7530e-3 l = 3 (3, 3) 437 322 240

• ut of memory

115 207 920 (26%) 28 (0.51) 7.3043e-5 l = 4 (4, 4) 874 644 480

• ut of memory

184 208 094 (21%) 37 (0.58) 3.0435e-6

Computing time for the largest computation: approx. 3 seconds per time step

15

Acoustic Waves: The Tunnel Experiment

t = 0.0 t = 0.6 t = 1.2 t = 1.8 t = 2.4 t = 3.0

16

Acoustic Waves: Estimator without dual Regularity

Let u ∈ V be the solution of Lu = f, let uh ∈ Vh be its approximation, and let u∗ ∈ V ∗ by the dual solution with (w, L∗u∗)0,Q = E, w for w ∈ W. Then we obtain for all wh ∈ Wh ∩ V ∗ the error representation E, u − uh =

• R=I×K∈R
• f − M∂tuh,R −div F(uh,R), u∗− wh
• 0,R+
• nK · F(uh,R), u∗− wh
• 0,I×∂K

and for acoustic waves

• E, (p − ph, v − vh)
• =
• R=I×K∈R
• − ∂tph,R + div vh,R, p∗ − ϕh
• 0,R +
• b − ∂tvh,R + ∇ph,R, v∗ − ψh
• 0,R

+ 1 2

• f∈FK
• [vh]K,f · nK, p∗ − ϕh
• 0,I×f +
• [ph]K,f, nK · (v∗ − ψh)
• 0,I×f
• .

Inserting the L2 projection Qh in K we define the local error estimate ηR =

• − κ−1∂tph,R + div vh,R
• 0,Rh1/2

K

• [Qhp∗

h]K

• 0,I×∂K

+

• b − ρ∂tvh,R + ∇ph,R
• 0,Rh1/2

K

• [Qhv∗

h]K · nK

• 0,I×∂K

+ 1 2

• f∈FK
• [vh]K,f0,I×f
• [Qhp∗

h]K,f

• 0,I×f +
• [σh]K,fnK0,I×f
• [Qhv∗

h]K,f · nK

• 0,I×f
• .

17

Acoustic Waves: The Tunnel Experiment

t = 0.0 t = 0.6 t = 1.2 t = 1.8 t = 2.4 t = 3.0

18

Acoustic Waves: The Tunnel Experiment

level (p, q) #DoFs (effort) mg steps (rate) El(uh) |El(uh) − Eex| uniform refinement l = 0 (1, 1) 801 792 8 (1.49e-1) 4.4925e-3 1.4647e-3 l = 1 (2, 2) 3 207 168 13 (2.95e-1) 5.8966e-3 6.0637e-5 l = 2 (3, 3) 8 017 920 20 (4.25e-1) 5.9523e-3 4.8958e-6 l = 3 (4, 4) 16 035 840 28 (5.79e-1) 5.9568e-3 3.9528e-7 estimator including numerical flux higher order interpolation l = 0 801 792 8 (1.49e-1) 4.4925e-3 1.4647e-3 l = 1 1 395 819 (44%) 13 (2.44e-1) 5.8907e-3 6.6512e-5 l = 2 2 469 042 (31%) 18 (3.85e-1) 5.9524e-3 4.7504e-6 l = 3 4 271 634 (27%) 28 (5.64e-1) 5.9566e-3 5.6460e-7 estimator without dual regularity assumptions l = 0 801 792 8 (1.49e-1) 4.4925e-3 1.4647e-3 l = 1 1 514 430 (47%) 13 (2.80e-1) 5.8966e-3 6.0647e-5 l = 2 2 807 919 (35%) 20 (4.27e-1) 5.9522e-3 4.9620e-6 l = 3 4 931 487 (31%) 27 (5.66e-1) 5.9567e-3 4.8728e-7 Uniform vs. adaptive refinement on 88 608 = 1 846 × 48 space-time cells on 64 procs (ϑ = 6e-4, extrapolated value of the goal functional Eex = 5.9572e-3).

19

Acoustic Waves: The Tunnel Experiment

106 107 10−7 10−6 10−5 10−4 10−3 10−2 (1, 1) (2, 2) (3, 3) (4, 4) #DoFs |El(uh) − Eex| uniform refinement estimator with numerical flux alternative adaptive refinement

20

Acoustic Waves: The Tunnel Experiment

21