[PPT] - Discontinuous Galerkin Methods and Local Time-Stepping For Transient PowerPoint Presentation

SLIDE 1

Discontinuous Galerkin Methods and Local Time-Stepping For Transient Wave Propagation

Marcus Grote

joint work with:

J. Diaz, INRIA
T. Mitkova, A. Schneebeli, Univ. Basel
D. Sch¨
tzau, UBC

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 2

High-Order Explicit LTS Methods

Outline:

Introduction
The (damped) wave equation
CG and DG Finite Element discretizations
The damped case: Adams-Bashforth LTS methods
Numerical experiments
The undamped case: Leap-Frog LTS methods
Numerical experiments
Concluding remarks

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 3

Introduction

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 4

The (Damped) Wave Equation

Model problem (second-order form) utt + σut − ∇ · ( c ∇u) = f in Ω × (0, T) u = 0

n ∂Ω × (0, T)

u|t=0 = u0 , ut|t=0 = v0 in Ω

Ω ⊂ Rd bounded, σ(x) ≥ 0, c(x) > 0

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 5

The (Damped) Wave Equation

Model problem (second-order form) utt + σut − ∇ · ( c ∇u) = f in Ω × (0, T) u = 0

n ∂Ω × (0, T)

u|t=0 = u0 , ut|t=0 = v0 in Ω

Ω ⊂ Rd bounded, σ(x) ≥ 0, c(x) > 0

Weak formulation Find u ∈ C0(0, T; H1

0(Ω)) ∩ C1(0, T; L2(Ω)):

utt, v(H−1,H1

0) + (σut, v) + a(u, v) = (f, v),

∀v ∈ H1

0(Ω),

with a(u, v) = (c ∇u, ∇v)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 6

The (Damped) Wave Equation

Model problem (second-order form) utt + σut − ∇ · ( c ∇u) = f in Ω × (0, T) u = 0

n ∂Ω × (0, T)

u|t=0 = u0 , ut|t=0 = v0 in Ω

Ω ⊂ Rd bounded, σ(x) ≥ 0, c(x) > 0

Weak formulation Find u ∈ C0(0, T; H1

0(Ω)) ∩ C1(0, T; L2(Ω)):

utt, v(H−1,H1

0) + (σut, v) + a(u, v) = (f, v),

∀v ∈ H1

0(Ω),

with a(u, v) = (c ∇u, ∇v) Energy conservation For σ = 0, f = 0 the energy E[u](t), E[u](t) := 1 2

ut2 + a(u, u)
≡ const.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 7

Conforming Galerkin approximation

Finite element space V h ⊂ H1

0(Ω):

V h = {v ∈ H1

0(Ω) : v|K ∈ Sℓ(K),

K ∈ Th} ⇒ elements continuous across edges / faces

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 8

Conforming Galerkin approximation

Finite element space V h ⊂ H1

0(Ω):

V h = {v ∈ H1

0(Ω) : v|K ∈ Sℓ(K),

K ∈ Th} ⇒ elements continuous across edges / faces Semi-discrete Galerkin formulation: Find uh : [0, T] × V h → R such that (uh

tt, v) + (σuh t , v) + a(uh, v) = (f, v)

∀v ∈ V h, t ∈ (0, T) Conforming mass-lumped FEM: (Cohen-Joly-Roberts-Tordjman, SINUM, 2001) a(u, v) :=

K∈Th
K

c∇u · ∇v dx ⇒ The semi-discrete formulation automatically inherits the main properties from the continuous formulation.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 9

Properties of conforming FEM

For σ = 0 and f = 0, the energy E[uh](t),

E[uh](t) = 1 2

∂tuh2 + a(uh, uh)
,

is conserved.

For sufficiently smooth solutions we have the optimal error

estimates: (∂tu − ∂tuh)(t, .)L2 + ∇(u − uh)(t, .)L2 = O(hℓ), h → 0 (u − uh)(t, .)L2 = O(hℓ+1), see Dupont, Baker, SINUM, 1970’s

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 10

Discontinuous Galerkin (DG) methods

Finite element space with no inter-elemental continuity

constrains, i.e. V h ⊂ / H1

0(Ω):

V h = {v ∈ L2(Ω) : v|K ∈ Sℓ(K), K ∈ Th}

Local (elementwise) weak formulation
Numerical fluxes weakly enforce inter-elemental continuity

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 11

Discontinuous Galerkin (DG) methods

Finite element space with no inter-elemental continuity

constrains, i.e. V h ⊂ / H1

0(Ω):

V h = {v ∈ L2(Ω) : v|K ∈ Sℓ(K), K ∈ Th}

Local (elementwise) weak formulation
Numerical fluxes weakly enforce inter-elemental continuity

+ Flexibility in mesh-design (non-matching grids) + easily handles varying polynomial degree (hp-adaptivity) + (Block-)Diagonal mass matrix ⇒ fully explicit time stepping!

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 12

Discontinuous Galerkin (DG) methods

Finite element space with no inter-elemental continuity

constrains, i.e. V h ⊂ / H1

0(Ω):

V h = {v ∈ L2(Ω) : v|K ∈ Sℓ(K), K ∈ Th}

Local (elementwise) weak formulation
Numerical fluxes weakly enforce inter-elemental continuity

+ Flexibility in mesh-design (non-matching grids) + easily handles varying polynomial degree (hp-adaptivity) + (Block-)Diagonal mass matrix ⇒ fully explicit time stepping!

increase the number of degrees of freedom
increase the condition number of the stiffness matrix

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 13

DG discretization: notation

Th: shape-regular mesh K: tetrahedra or hexahedra Fh: faces of K ∈ Th f ∈ Fh, f = K .

+ ∩ K

.

−

n±: outward normal on ∂K± v±: traces of v from K±

K+ K! n! n+

Jumps: [ [v] ] := n+v+ + n−v− Averages: { {v} } := (v+ + v−)/2

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 14

Interior Penalty (IP) formulation

u, v ∈ Vh:

ah(u, v) := (c∇hu, ∇hv) −

Fh

[ [v] ] · { {c ∇hu} } ds u, v ∈ Vh := {v ∈ L2(Ω) : v|K ∈ Sℓ(K) ∀K ∈ Th}

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 15

Interior Penalty (IP) formulation

u, v ∈ Vh:

ah(u, v) := (c∇hu, ∇hv) −

Fh

[ [v] ] · { {c ∇hu} } ds −

Fh

[ [u] ] · { {c ∇hv} } ds u, v ∈ Vh := {v ∈ L2(Ω) : v|K ∈ Sℓ(K) ∀K ∈ Th}

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 16

Interior Penalty (IP) formulation

u, v ∈ Vh:

ah(u, v) := (c∇hu, ∇hv) −

Fh

[ [v] ] · { {c ∇hu} } ds −

Fh

[ [u] ] · { {c ∇hv} } ds +

Fh

a [ [u] ] · [ [v] ] ds. u, v ∈ Vh := {v ∈ L2(Ω) : v|K ∈ Sℓ(K) ∀K ∈ Th}

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 17

Interior Penalty (IP) formulation

u, v ∈ Vh:

ah(u, v) := (c∇hu, ∇hv) −

Fh

[ [v] ] · { {c ∇hu} } ds −

Fh

[ [u] ] · { {c ∇hv} } ds +

Fh

a [ [u] ] · [ [v] ] ds. u, v ∈ Vh := {v ∈ L2(Ω) : v|K ∈ Sℓ(K) ∀K ∈ Th}

Interior penalty stabilization function

a := α c h−1 ∈ L∞(Fh) α > 0: IP stabilization parameter h|F := min{hK+, hK−}, hK diameter of element K c|F := max{c|K+(x), c|K−(x)}

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 18

Semi-discrete IP-DG approximation

Find uh : [0, T] × V h → R such that (uh

tt, v) + (σuh t , v) + ah(uh, v) = (f, v)

∀v ∈ V h, t ∈ (0, T), uh|t=0 = Πhu0, uh

t |t=0 = Πhv0.

Πh : L2-projection onto V h The DG-bilinear form ah(u, v) is

symmetric, continuous, and coercive, for α ≥ αmin,

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 19

Semi-discrete IP-DG approximation

Find uh : [0, T] × V h → R such that (uh

tt, v) + (σuh t , v) + ah(uh, v) = (f, v)

∀v ∈ V h, t ∈ (0, T), uh|t=0 = Πhu0, uh

t |t=0 = Πhv0.

Πh : L2-projection onto V h The DG-bilinear form ah(u, v) is

symmetric, continuous, and coercive, for α ≥ αmin,
and the solution conserves the (discrete) energy

Eh(t) := 1 2uh

t (t)2 0 + 1

2ah(uh(t), uh(t)) if σ = 0 and f = 0.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 20

IP-DG for second-order wave equations

Symmetric bilinear form ah ⇒ conservation of discrete

energy (G., Schneebeli, Sch¨

tzau:

SINUM, 2006; JCAM 2007; IMA J. Num. Anal. 2008) (G. and Sch¨

tzau: J. Sc. Comp., 2009)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 21

IP-DG for second-order wave equations

Symmetric bilinear form ah ⇒ conservation of discrete

energy

mass matrix (block-) diagonal ⇒ fully explicit time

integration (G., Schneebeli, Sch¨

tzau:

SINUM, 2006; JCAM 2007; IMA J. Num. Anal. 2008) (G. and Sch¨

tzau: J. Sc. Comp., 2009)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 22

IP-DG for second-order wave equations

Symmetric bilinear form ah ⇒ conservation of discrete

energy

mass matrix (block-) diagonal ⇒ fully explicit time

integration

semi-discrete convergence analysis:
optimal a-priori error bounds in energy norm
optimal a-priori error bounds in L2 norm

(G., Schneebeli, Sch¨

tzau:

SINUM, 2006; JCAM 2007; IMA J. Num. Anal. 2008) (G. and Sch¨

tzau: J. Sc. Comp., 2009)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 23

IP-DG for second-order wave equations

Symmetric bilinear form ah ⇒ conservation of discrete

energy

mass matrix (block-) diagonal ⇒ fully explicit time

integration

semi-discrete convergence analysis:
optimal a-priori error bounds in energy norm
optimal a-priori error bounds in L2 norm
fully discrete convergence analysis (L-F in time):
optimal a-priori error bounds in L2 norm
CFL condition: ∆t ≤ CSh / √α cmax

(G., Schneebeli, Sch¨

tzau:

SINUM, 2006; JCAM 2007; IMA J. Num. Anal. 2008) (G. and Sch¨

tzau: J. Sc. Comp., 2009)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 24

IP-DG for second-order wave equations

Symmetric bilinear form ah ⇒ conservation of discrete

energy

mass matrix (block-) diagonal ⇒ fully explicit time

integration

semi-discrete convergence analysis:
optimal a-priori error bounds in energy norm
optimal a-priori error bounds in L2 norm
fully discrete convergence analysis (L-F in time):
optimal a-priori error bounds in L2 norm
CFL condition: ∆t ≤ CSh / √α cmax
Optimal convergence also proved for Maxwell’s Equations.

(G., Schneebeli, Sch¨

tzau:

SINUM, 2006; JCAM 2007; IMA J. Num. Anal. 2008) (G. and Sch¨

tzau: J. Sc. Comp., 2009)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 25

The (Damped) Wave Equation

Model problem (first-order form, v := ut and w := −∇u) vt + σv + ∇ · (c w) = f in Ω × (0, T) wt + ∇v = 0 in Ω × (0, T) v = 0 , w = 0

n ∂Ω × (0, T)

v|t=0 = v0 , w|t=0 = −∇u0 in Ω qt + Σ q + ∇ · F(q) = S with q = (v, w)t

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 26

Nodal DG FE Formulation

Find qh : [0, T] × Vh → R such that (qh

t , ψ)+(Σ qh, ψ)+aDG(qh, ψ) = (S, ψ)

∀ ψ ∈ Vh , t ∈ (0, T) .

Nodal DG FEM: (Hesthaven-Warburton, Springer, 2008)

aDG(q, ψ) :=

K∈Th
K

(∇ · F(q)) · ψ dx −

e∈Eh
e

(n · F(q) − (n · F(q))∗) · ψ dA Here, (n · F(q))∗ denotes a suitably chosen numerical flux in the unit normal direction n.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 27

Fully Discrete FE Formulation

The discretization in space leads to either system of ODE’s: M d2U dt2 (t) + Mσ dU dt (t) + K U(t) = F(t) , t ∈ (0, T)

r

MdQ dt (t) + Mσ Q(t) + K Q(t) = F(t) , t ∈ (0, T) . The stiffness matrix K and the mass matrices M, Mσ are SPD. Moreover, the mass matrix M is (block-)diagonal ⇒ computing M−1 is cheap ⇒ fully explicit time-stepping!

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 28

Fully Discrete FE Formulation

The discretization in space leads to either system of ODE’s: M d2U dt2 (t) + Mσ dU dt (t) + K U(t) = F(t) , t ∈ (0, T)

r

MdQ dt (t) + Mσ Q(t) + K Q(t) = F(t) , t ∈ (0, T) . The stiffness matrix K and the mass matrices M, Mσ are SPD. Moreover, the mass matrix M is (block-)diagonal ⇒ computing M−1 is cheap ⇒ fully explicit time-stepping! adaptivity, small geometric features ⇓ locally refined meshes ⇓ CFL condition for explicit time-stepping ∆t ≤ C h, h = min

T ∈Th hT Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 29

Local Time-Stepping: Previous Work

Collino-Fouquet-Joly, Numer. Math. 2003, JCP 2006

Piperno, M2AN 2006; Cohen-Montseny-Pernet-Ferri´ eres, JCP 2008

First-order wave and Maxwell system
Second-order accuracy and conservation of energy
Dumbser-K¨

aser-Toro, Geophys. J. Int. 2007 Dumbser-Munz-Schneider-Taube, Int. J. Numer. Model. 2009

LTS ADER-DG schemes
Constantinescu-Sandu, SISC 2009
Multirate explicit Adams methods for conservation laws
Limited to second order accuracy
Diaz-G., SISC 2009
Second-order wave equation (σ = 0)
LTS-LF of arbitrarily high accuracy; conservation of energy
G.-Mitkova, JCAM 2010
Maxwell’s equations in second-order form
σ = 0: LTS-LF, arbitrarily high accuracy
σ > 0: LTS-LF/CN, second order accuracy

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 30

High-Order Explicit LTS

We consider the semi-discrete problems

M d2U dt2 (t)+Mσ dU dt (t)+K U(t) = F(t) resp. MdQ dt (t)+Mσ Q(t)+K Q(t) = F(t)

and rewrite them as

d2Z dt2 (t) + DdZ dt (t) + AZ(t) = R(t) resp. de Z dt (t) + D e Z(t) + A e Z(t) = R(t) Z=M

1 2 U, e

Z=M

1 2 Q, D=M− 1 2 MσM− 1 2 , A=M− 1 2 KM− 1 2 , R=M− 1 2 F.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 31

High-Order Explicit LTS

We consider the semi-discrete problems

M d2U dt2 (t)+Mσ dU dt (t)+K U(t) = F(t) resp. MdQ dt (t)+Mσ Q(t)+K Q(t) = F(t)

and rewrite them as

d2Z dt2 (t) + DdZ dt (t) + AZ(t) = R(t) resp. de Z dt (t) + D e Z(t) + A e Z(t) = R(t) Z=M

1 2 U, e

Z=M

1 2 Q, D=M− 1 2 MσM− 1 2 , A=M− 1 2 KM− 1 2 , R=M− 1 2 F.

Goal: Derive explicit LTS schemes of arbitrarily high accuracy starting from explicit multi-step Adams-Bashforth methods for dY dt (t) = BY(t) , t ∈ (0, T) (R(t) ≡ 0)

Y(t) = „ Z(t), dZ dt (t) «T , B = „ I −A −D «

r Y(t) = e

Z(t) , B = −A−D .

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 32

High-Order Explicit LTS

Classical k-step Adams-Bashforth method We consider dY dt (t) = BY(t) in integrated form Y(tn + ξ∆t) = Y(tn) + tn+ξ∆t

tn

BY(t) dt , 0 < ξ ≤ 1 , and replace Y(t) by the polynomial interpolant through (ti, yi), i = n − k + 1, . . . , n. Then, we have Y(tn + ξ∆t) ≈ Yn+ξ = Yn + ∆tB

k−1

j=0

γj(ξ)∇jYn , γj(ξ) ∈ Pj+1 , ∇0Yn = Yn , ∇j+1Yn = ∇jYn − ∇jYn−1 .

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 33

High-Order Explicit LTS

LTS-ABk(p) method Let us now split Y in two parts Y =

Ycoarse
+
Yfine
= (I − P)Y + PY, with P2 = P.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 34

High-Order Explicit LTS

LTS-ABk(p) method Let us now split Y in two parts Y =

Ycoarse
+
Yfine
= (I − P)Y + PY, with P2 = P.

Then, we have d dtY = B(I − P)Y + BPY

r

Y(tn + ξ∆t) = Y(tn) + Z tn+ξ∆t

tn

B(I − P)Y(t) dt + Z tn+ξ∆t

tn

BPY(t) dt .

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 35

High-Order Explicit LTS

Y(tn + ξ∆t) = Y(tn) + Z tn+ξ∆t

tn

B(I − P)Y(t) dt + Z tn+ξ∆t

tn

BPY(t) dt ≈ Yn + ∆tB(I − P)

k−1

X

j=0

γj(ξ)∇jYn + Z ξ∆t BP e Y(τ) dτ Where e Y is the solution of ˛ ˛ ˛ ˛ ˛ ˛ ˛ e Y(0) = Yn d dτ e Y(τ) = B(I − P)

k−1

X

j=0

d dτ γj “ τ ∆t ” ∇jYn + BP e Y(τ)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 36

High-Order Explicit LTS

Y(tn + ξ∆t) = Y(tn) + Z tn+ξ∆t

tn

B(I − P)Y(t) dt + Z tn+ξ∆t

tn

BPY(t) dt ≈ Yn + ∆tB(I − P)

k−1

X

j=0

γj(ξ)∇jYn + Z ξ∆t BP e Y(τ) dτ Where e Y is the solution of ˛ ˛ ˛ ˛ ˛ ˛ ˛ e Y(0) = Yn d dτ e Y(τ) = B(I − P)

k−1

X

j=0

d dτ γj “ τ ∆t ” ∇jYn + BP e Y(τ) e Y(ξ∆t) = Yn + ∆tB(I − P)

k−1

X

j=0

γj(ξ)∇jYn + Z ξ∆t BP e Y(τ) dτ Y(tn + ξ∆t) ≈ e Y(ξ∆t) = ⇒ Yn+1 = e Y(∆t)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 37

High-Order Explicit LTS

Y(tn + ξ∆t) = Y(tn) + Z tn+ξ∆t

tn

B(I − P)Y(t) dt + Z tn+ξ∆t

tn

BPY(t) dt ≈ Yn + ∆tB(I − P)

k−1

X

j=0

γj(ξ)∇jYn + Z ξ∆t BP e Y(τ) dτ Where e Y is the solution of ˛ ˛ ˛ ˛ ˛ ˛ ˛ e Y(0) = Yn d dτ e Y(τ) = B(I − P)

k−1

X

j=0

d dτ γj “ τ ∆t ” ∇jYn + BP e Y(τ) e Y(ξ∆t) = Yn + ∆tB(I − P)

k−1

X

j=0

γj(ξ)∇jYn + Z ξ∆t BP e Y(τ) dτ Y(tn + ξ∆t) ≈ e Y(ξ∆t) = ⇒ Yn+1 = e Y(∆t) We solve the problem ( e Y) from τ = 0 to τ = ∆t, using a k-step Adams-Bashforth method with ∆τ = ∆t/p. THEOREM: The LTS-ABk(p) scheme is kth-order accurate in time.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 38

High-Order Explicit LTS: Algorithm AB3(2)

AB3: Yn+1 = Yn + ∆tB »23 12Yn−16 12Yn−1 + 5 12Yn−2 – LTS-AB3(2): k = 3, p = 2 Wj := B(I−P)Yj, j = n, n−1, n−2 , e Y0 := Yn, e Y− 1

2 := Yn− 1 2 , e

Y−1 := Yn−1

tn−1 tn−1/2 tn tn+1/2 tn+1

Yn−1 Yn

Y−1/2
Y0

tn−2

Yn−2

Y−1

Solve from 0 to ∆t with ∆τ = ∆t

2

d dτ e Y(τ) =B(I − P)

k−1

X

j=0

d dτ γj “ τ ∆t ” ∇jYn + BP e Y(τ)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 39

High-Order Explicit LTS: Algorithm AB3(2)

AB3: Yn+1 = Yn + ∆tB »23 12Yn−16 12Yn−1 + 5 12Yn−2 – LTS-AB3(2): k = 3, p = 2 Wj := B(I−P)Yj, j = n, n−1, n−2 , e Y0 := Yn, e Y− 1

2 := Yn− 1 2 , e

Y−1 := Yn−1 e Y

1 2 = e

Y0+∆t 2 »17 12Wn− 7 12Wn−1+ 2 12Wn−2 – +∆t 2 BP »23 12 e Y0−16 12 e Y− 1

2 + 5

12 e Y−1 –

tn−1 tn−1/2 tn tn+1/2 tn+1

Yn−1 Yn

Y−1/2
Y0

tn−2

Yn−2

Y−1
Y1/2

Integration from 0 to ∆t/2 of d dτ e Y(τ) =B(I − P)

k−1

X

j=0

d dτ γj “ τ ∆t ” ∇jYn + BP e Y(τ)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 40

High-Order Explicit LTS: Algorithm AB3(2)

AB3: Yn+1 = Yn + ∆tB »23 12Yn−16 12Yn−1 + 5 12Yn−2 – LTS-AB3(2): k = 3, p = 2 Wj := B(I−P)Yj, j = n, n−1, n−2 , e Y0 := Yn, e Y− 1

2 := Yn− 1 2 , e

Y−1 := Yn−1 e Y

1 2 = e

Y0+∆t 2 »17 12Wn− 7 12Wn−1+ 2 12Wn−2 – +∆t 2 BP »23 12 e Y0−16 12 e Y− 1

2 + 5

12 e Y−1 – e Y1 = e Y

1 2 +∆t

2 »29 12Wn−25 12Wn−1+ 8 12Wn−2 – +∆t 2 BP »23 12 e Y

1 2 −16

12 e Y0+ 5 12 e Y− 1

2

–

tn−1 tn−1/2 tn tn+1/2 tn+1

Yn−1 Yn

Y−1/2
Y0

tn−2

Yn−2

Y−1
Y1/2
Y1

Yn+1

Integration from ∆t/2 to ∆t of d dτ e Y(τ) =B(I − P)

k−1

X

j=0

d dτ γj “ τ ∆t ” ∇jYn + BP e Y(τ)

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 41

Numerical Experiments 1D: LTS-AB3(p)

Computational domain:

Ω = [0, 6] , Ωcoarse = [0, 2] ∪ [4, 6] , Ωfine = [2, 4] , p = 2, 4, 13

Exact solution: (c ≡ 1)

u(x, t) = 2e− σt

2

√ 4π2 − σ2 sin(πx) sin t 2

4π2 − σ2
Homogeneous source data
Homogeneous Dirichlet boundary condition
Damping parameter: σ ≡ 0.1
Space discretization with P2 mass-lumped, IP-DG and

nodal DG FEM

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 42

Numerical Experiments 1D: LTS-AB3(p)

local refinement p = 2, ∆t = ∆tAB3 (optimal CFL condition)

∆tAB3 denotes the largest time-step allowed by the classical AB3 scheme solution at T = 10

10

−2

10

−1

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

meshsize hcoarse L2−error

cont. FEM

IP−DG nodal DG hcoarse

3

L2-error at T = 10

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 43

Numerical Experiments 1D: LTS-AB3(p)

local refinement p = 4, ∆t = ∆tAB3 (optimal CFL condition)

∆tAB3 denotes the largest time-step allowed by the classical AB3 scheme solution at T = 10

10

−2

10

−1

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

meshsize hcoarse L2−error

cont. FEM

IP−DG nodal DG hcoarse

3

L2-error at T = 10

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 44

Numerical Experiments 1D: LTS-AB3(p)

local refinement p = 13, ∆t = ∆tAB3 (optimal CFL condition)

∆tAB3 denotes the largest time-step allowed by the classical AB3 scheme solution at T = 10

10

−2

10

−1

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

meshsize hcoarse L2−error

cont. FEM

IP−DG nodal DG hcoarse

3

L2-error at T = 10

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 45

Numerical experiments 2D: Roof mounted antenna

Computational domain Ω × (0, 1.3): Ω is rectangular of size

[0, 3] × [0, 1] adjacent to a roof mounted antenna

Mesh with ratio of local refinement p = 7 resolves the small

geometric features of the antenna

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 46

Numerical experiments 2D: Roof mounted antenna

Computational domain Ω × (0, 1.3): Ω is rectangular of size

[0, 3] × [0, 1] adjacent to a roof mounted antenna

Mesh with ratio of local refinement p = 7 resolves the small

geometric features of the antenna

Homogeneous source data
Homogeneous Neumann BC
Model parameters: c ≡ 1, σ = 0.1
A plane wave is excited through the initial conditions
LTS-AB3(7) combined with P2 mass-lumped FEM
Initial time steps: Runge-Kutta method of order 4

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 47

Numerical Experiments 2D: Roof mounted antenna

t = 0 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 48

Numerical Experiments 2D: Roof mounted antenna

t = 0.1 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 49

Numerical Experiments 2D: Roof mounted antenna

t = 0.2 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 50

Numerical Experiments 2D: Roof mounted antenna

t = 0.4 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 51

Numerical Experiments 2D: Roof mounted antenna

t = 0.55 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 52

Numerical Experiments 2D: Roof mounted antenna

t = 0.71 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 53

Numerical Experiments 2D: Roof mounted antenna

t = 0.82 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 54

Numerical Experiments 2D: Roof mounted antenna

t = 1 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 55

Numerical Experiments 2D: Roof mounted antenna

t = 1.2 Roof mounted antenna: numerical solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 56

The undamped case (σ = 0): LTS L-F methods

The spatial discretization of the wave equation with no damping, either with conforming FE with mass lumping or with the IP-DG method, yields M d2 dt2 U + KU = F

M and K are symmetric positive definite matrices
M is (block-) diagonal

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 57

Local Time Stepping: L-F methods

Let F = 0 and multiply by M−1/2. This yields d2 dt2 Y + AY = 0, Y = M1/2U, A = M−1/2KM−1/2 with A spd and still sparse. Then [ ˙ Y 2 + Y ⊤AY ]/2 = const.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 58

Local Time Stepping: L-F methods

Let F = 0 and multiply by M−1/2. This yields d2 dt2 Y + AY = 0, Y = M1/2U, A = M−1/2KM−1/2 with A spd and still sparse. Then [ ˙ Y 2 + Y ⊤AY ]/2 = const. The exact solution satisfies: Y (t+∆t)−2Y (t)+Y (t−∆t) = −∆t2 1

−1

(1−|τ|)AY (t+τ∆t) dτ.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 59

Local Time Stepping: L-F methods

Let F = 0 and multiply by M−1/2. This yields d2 dt2 Y + AY = 0, Y = M1/2U, A = M−1/2KM−1/2 with A spd and still sparse. Then [ ˙ Y 2 + Y ⊤AY ]/2 = const. The exact solution satisfies: Y (t+∆t)−2Y (t)+Y (t−∆t) = −∆t2 1

−1

(1−|τ|)AY (t+τ∆t) dτ. Replacing Y (t + τ∆t) ≈ Y (t) yields the leap-frog scheme Y (t + ∆t) − 2Y (t) + Y (t − ∆t) = −∆t2AY (t) + O(∆t4).

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 60

Local Time Stepping: L-F methods

Again, split Y as: Y = Y coarse Y fine

−1

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

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 61

Local Time Stepping: L-F methods

Again, split Y as: Y =

Y c
+

Y f

= (I − P)Y + PY

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

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 62

LTS L-F methods: Algorithm

Can show Y (t+∆t)−2Y (t)+Y (t−∆t) = 2(Q(∆t)−Q(0))+O(∆t4) where

Q(0) = Y n

Q′(0) = 0 d2 dτ 2 Q(τ) + A(I − P)Y n + APQ(τ) = 0 from τ = 0 to τ = ∆t, solved with L-F, with ∆τ = ∆t p .

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 63

LTS L-F methods: Algorithm

Can show Y (t+∆t)−2Y (t)+Y (t−∆t) = 2(Q(∆t)−Q(0))+O(∆t4) where

Q(0) = Y n

Q′(0) = 0 d2 dτ 2 Q(τ) + A(I − P)Y n + APQ(τ) = 0 from τ = 0 to τ = ∆t, solved with L-F, with ∆τ = ∆t p . Then Y n+1 = 2Y n − Y n−1 + 2(Q(∆t) − Q(0))

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 64

LTS L-F methods: Algorithm

Can show Y (t+∆t)−2Y (t)+Y (t−∆t) = 2(Q(∆t)−Q(0))+O(∆t4) where

Q(0) = Y n

Q′(0) = 0 d2 dτ 2 Q(τ) + A(I − P)Y n + APQ(τ) = 0 from τ = 0 to τ = ∆t, solved with L-F, with ∆τ = ∆t p . Then Y n+1 = 2Y n − Y n−1 + 2(Q(∆t) − Q(0)) This algorithm requires only one multiplication by A(I − P) and p multiplications by AP.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 65

LTS L-F methods: energy conservation

This algorithm is equivalent to Y n+1 = 2Y n − Y n−1 − ∆t2ApY n with Ap = A − 2 p2

p−1

j=1

∆t p 2j αp

j(AP)jA symmetric.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 66

LTS L-F methods: energy conservation

This algorithm is equivalent to Y n+1 = 2Y n − Y n−1 − ∆t2ApY n with Ap = A − 2 p2

p−1

j=1

∆t p 2j αp

j(AP)jA symmetric.

Hence it conserves the discrete energy: En+ 1

2

= 1 2

I − ∆t2

4 Ap Y n+1 − Y n ∆t , Y n+1 − Y n ∆t

+

1 2

Ap

Y n+1 + Y n 2 , Y n+1 + Y n 2

Marcus Grote

RICAM workshop, Linz, 21–25.11.2011

SLIDE 67

Experiments in 2D: P3-elements & 4th order in time

hcoarse = 0.0125, hfine = 7.62.10−4 ≈ hcoarse/17

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 68

Experiments in 2D: P3-elements & 4th order in time

hcoarse = 0.0125, hfine = 7.62.10−4 ≈ hcoarse/17

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 69

Experiments in 2D: P3-elements & 4th order in time

hcoarse = 0.0125, hfine = 7.62.10−4 ≈ hcoarse/17

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 70

Local Time-Stepping: L-F methods

Example: (Diaz-G., SISC 2009)

Ω = [−0.5 ; 0.5]2 with Neumann boundary condition,

width of the slot: 0.004

IP-DG, P3-elements with α = 11
LTS-LF of order 4, using modified equation approach

hfine ≈ hcoarse/17 solution at t = 0.09

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 71

Local Time-Stepping: L-F methods

Example: (Diaz-G., SISC 2009)

Ω = [−0.5 ; 0.5]2 with Neumann boundary condition,

width of the slot: 0.004

IP-DG, P3-elements with α = 11
LTS-LF of order 4, using modified equation approach

hfine ≈ hcoarse/17 solution at t = 0.17

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 72

Local Time-Stepping: L-F methods

Example: (Diaz-G., SISC 2009)

Ω = [−0.5 ; 0.5]2 with Neumann boundary condition,

width of the slot: 0.004

IP-DG, P3-elements with α = 11
LTS-LF of order 4, using modified equation approach

hfine ≈ hcoarse/17 solution at t = 0.35

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 73

Local Time-Stepping: L-F methods

Example: (Diaz-G., SISC 2009)

Ω = [−0.5 ; 0.5]2 with Neumann boundary condition,

width of the slot: 0.004

IP-DG, P3-elements with α = 11
LTS-LF of order 4, using modified equation approach

hfine ≈ hcoarse/17 solution at t = 0.44

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 74

Local Time-Stepping: L-F methods

Example: (Diaz-G., SISC 2009)

Ω = [−0.5 ; 0.5]2 with Neumann boundary condition,

width of the slot: 0.004

IP-DG, P3-elements with α = 11
LTS-LF of order 4, using modified equation approach

hfine ≈ hcoarse/17 solution at t = 0.61

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 75

Local Time-Stepping: L-F methods

Example: (Diaz-G., SISC 2009)

Ω = [−0.5 ; 0.5]2 with Neumann boundary condition,

width of the slot: 0.004

IP-DG, P3-elements with α = 11
LTS-LF of order 4, using modified equation approach

hfine ≈ hcoarse/17 solution at t = 0.70

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 76

Multi-level local time stepping (L-F & DG-P 3)

Computational Domain

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 77

Multi-level local time stepping (L-F & DG-P 3)

Mesh

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 78

Multi-level local time stepping (L-F & DG-P 3)

Top left corner p = 3

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 79

Multi-level local time stepping (L-F & DG-P 3)

Bottom left corner p1 = 3, p2 = 3

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 80

Multi-level local time stepping (L-F & DG-P 3)

Bottom right corner p1 = 2, p2 = 5, p3 = 4

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 81

Multi-level local time stepping (L-F & DG-P 3)

Top right corner p1 = 3, p2 = 3, p3 = 3, p4 = 4

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 82

Multi-level local time stepping (L-F & DG-P 3)

∆t = ∆topt

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 83

Multi-level local time stepping (L-F & DG-P 3)

Solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 84

Multi-level local time stepping (L-F & DG-P 3)

Solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 85

Multi-level local time stepping (L-F & DG-P 3)

Solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 86

Multi-level local time stepping (L-F & DG-P 3)

Solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 87

Multi-level local time stepping (L-F & DG-P 3)

Solution

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 88

Concluding Remarks

Wave equation with (or without) damping
Mass-lumped/IP-DG/nodal DG FE discretization ⇒

block-diagonal mass matrix ⇒ explicit time integration

Proposed explicit LTS schemes of arbitrarily high accuracy:
With damping (σ ≤ 0): LTS AB(k) methods

⇒ optimal CFL for k ≥ 3 (arXiv:1109.4480v1, 2011).

With no damping (σ = 0): LTS L-F type methods

⇒ discrete energy conserved (SISC 2009, JCAM 2010).

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 89

Concluding Remarks

Wave equation with (or without) damping
Mass-lumped/IP-DG/nodal DG FE discretization ⇒

block-diagonal mass matrix ⇒ explicit time integration

Proposed explicit LTS schemes of arbitrarily high accuracy:
With damping (σ ≤ 0): LTS AB(k) methods

⇒ optimal CFL for k ≥ 3 (arXiv:1109.4480v1, 2011).

With no damping (σ = 0): LTS L-F type methods

⇒ discrete energy conserved (SISC 2009, JCAM 2010).

Current work:

Exponential LTS-ABk schemes (Hochbruck & Ostermann)

⇒ more expensive, but additional work independent of p.

Explicit LTS-RK schemes.

Marcus Grote RICAM workshop, Linz, 21–25.11.2011

SLIDE 90

Concluding Remarks

Wave equation with (or without) damping
Mass-lumped/IP-DG/nodal DG FE discretization ⇒

block-diagonal mass matrix ⇒ explicit time integration

Proposed explicit LTS schemes of arbitrarily high accuracy:
With damping (σ ≤ 0): LTS AB(k) methods

⇒ optimal CFL for k ≥ 3 (arXiv:1109.4480v1, 2011).

With no damping (σ = 0): LTS L-F type methods

⇒ discrete energy conserved (SISC 2009, JCAM 2010).

Current work:

Exponential LTS-ABk schemes (Hochbruck & Ostermann)

⇒ more expensive, but additional work independent of p.

Explicit LTS-RK schemes.

Thank you for your attention!

Marcus Grote RICAM workshop, Linz, 21–25.11.2011