[PPT] - On Energy Stable dG Approximation of the PML for the Wave Equation PowerPoint Presentation

SLIDE 1

On Energy Stable dG Approximation of the PML for the Wave Equation

Monash Workshop on Numerical Differential Equations and Applications February 2020 Kenneth Duru

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications1 /

SLIDE 2

Waves are everywhere

Simulations of seismic waves to quantify and assess earthquake risks and hazards.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications2 /

SLIDE 3

Wave propagation:

Forward Modeling

Efficient Time Domain Wave Propagation Tool Accurate and stable volume discretizations Efficient and reliable absorbing boundaries Accurate source generation

Efficient and scalable implementation on modern HPC platforms. My research has penetrated all components.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications3 /

SLIDE 4

Truncated Domain

Which boundary conditions ensure that numerical simulations converge to the infinite domain problem? (old but relevant: Engquist and Majda (1977)).

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications4 /

SLIDE 5

Reflections from boundaries

A solution of the acoustic pressure with a point source. 2 4 6 8 10 t[s]

0.4
0.2

0.2 0.4 p[MPa]

∆x = 5/9 ∆x = 5/27 Analytical

Classical absorbing boundary condition

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications5 /

SLIDE 6

Absorbing Layer

Equations must be perfectly matched: J. P . B´ erenger (1994).

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications6 /

SLIDE 7

Anechoic chamber

”non-reflective”, ”non-echoing”, ”echo-free”. A room designed to completely absorb reflections of either sound or electromagnetic waves.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications7 /

SLIDE 8

Complex coordinate stretching: Chew and Weedon (1994)

∂/∂x → 1/Sx∂/∂x, Sx := d˜ x/dx = 1 + dx(x)/s, dx ≥ 0. Simplifies PML construction for hyperbolic systems

−2 −1 1 2 3 −1 −0.5 0.5 1

UPML x

PML

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications8 /

SLIDE 9

Acoustic wave equation in first order form

1 κ ∂p ∂t + ∇ · v = 0, ρ ∂v ∂t + ∇p = 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications9 /

SLIDE 10

Acoustic wave equation in first order form

1 κ ∂p ∂t + ∇ · v = 0, ρ ∂v ∂t + ∇p = 0. (x, y, z) ∈ Ω = [−1, 1]3, BCs: 1 − rη 2 Zvη ∓ 1 + rη 2 p = 0, |rη| ≤ 1, at η = ±1.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications9 /

SLIDE 11

Acoustic wave equation in first order form

1 κ ∂p ∂t + ∇ · v = 0, ρ ∂v ∂t + ∇p = 0. (x, y, z) ∈ Ω = [−1, 1]3, BCs: 1 − rη 2 Zvη ∓ 1 + rη 2 p = 0, |rη| ≤ 1, at η = ±1.

dE dxdydz = 1 2   1 κ |p|2 + ρ

η=x,y,z

|vη|2   > 0, E(t) =

Ω

dE > 0.

d dt E(t) = −

∂Ω

p (n · v) dS = −

η=x,y,z

1

−1

1

−1

BT(η) dydzdx dη ≤ 0, BT(η) = 1 − |rη|2 Z |χ(−η)|2

η=−1 + 1 − |rη|2

Z |χ(+η)|2

η=1 ≥ 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications9 /

SLIDE 12

Derive PML in the Laplace domain

u(x, y, z, s) =

∞ e−stu (x, y, z, t) dt, s = a + ib, Re{s} = a > 0,

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications10

SLIDE 13

Derive PML in the Laplace domain

u(x, y, z, s) =

∞ e−stu (x, y, z, t) dt, s = a + ib, Re{s} = a > 0, 1 κ s p + ∇ · v = 0, ρs v + ∇ p = 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications10

SLIDE 14

Derive PML in the Laplace domain

u(x, y, z, s) =

∞ e−stu (x, y, z, t) dt, s = a + ib, Re{s} = a > 0, 1 κ s p + ∇ · v = 0, ρs v + ∇ p = 0. PML : ∂/∂η → 1/Sη∂/∂η, Sη = 1 + dη(η) s , dη(η) ≥ 0, η = x, y, z.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications10

SLIDE 15

Derive PML in the Laplace domain

u(x, y, z, s) =

∞ e−stu (x, y, z, t) dt, s = a + ib, Re{s} = a > 0, 1 κ s p + ∇ · v = 0, ρs v + ∇ p = 0. PML : ∂/∂η → 1/Sη∂/∂η, Sη = 1 + dη(η) s , dη(η) ≥ 0, η = x, y, z. 1 κ s p + ∇d · v = 0, ρs v + ∇d p = 0, ∇d = (1/Sx∂/∂x, 1/Sy∂/∂y, 1/Sz∂/∂z)T .

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications10

SLIDE 16

Time-domain PML

Auxiliary variables: s σ =

dx − dy

Sx Sy ∂ v ∂y , s ψ =

dx − dz

Sx Sz ∂ w ∂z . Modified PDE: 1 κ ∂p ∂t +dxp

+ ∇ · v−σ − ψ = 0,

ρ ∂v ∂t +dv

+ ∇p = 0,

d =    dx dy dz    .

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications11

SLIDE 17

Time-domain PML

Auxiliary variables: s σ =

dx − dy

Sx Sy ∂ v ∂y , s ψ =

dx − dz

Sx Sz ∂ w ∂z . Modified PDE: 1 κ ∂p ∂t +dxp

+ ∇ · v−σ − ψ = 0,

ρ ∂v ∂t +dv

+ ∇p = 0,

d =    dx dy dz    . Auxiliary differential equation: ODE ∂σ ∂t + dyσ

+ (dy − dx) ∂v

∂y = 0, ∂ψ ∂t + dzψ

+ (dz − dx) ∂w

∂z = 0, BCs: 1 − rη 2 Zvη ∓ 1 + rη 2 p = 0, at η = ±1.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications11

SLIDE 18

2D: SBP finite difference approximation

1000 2000 3000 4000 5000 10

−3

10

−2

10

−1

10 time Ezh 2nd−order 4th−order 6th−order

(b) Discrete PML

1000 2000 3000 4000 5000 10

−3

10

−2

10

−1

10 time Ezh 2nd−order 4th−order 6th−order

(c) No PML (dx = 0).

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications12

SLIDE 19

A nightmare for finite element and dG practitioners!

60
40
20

20 40 60 x[km] 10 20 30 40 50 y[km] t=150 s 0.02 0.04 0.06 0.08 0.1

GLL.

60
40
20

20 40 60 x[km] 10 20 30 40 50 y[km] t=35 s 0.02 0.04 0.06 0.08 0.1

GL.

60
40
20

20 40 60 x[km] 10 20 30 40 50 y[km] t=10 s 0.02 0.04 0.06 0.08 0.1

GLR.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications13

SLIDE 20

Some pioneering work: Cauchy problem

Abarbanel and Gottlieb (1997 & 98), Hesthaven et al. (1999), Collino and Tsogka (2001), B´ ecache et al. (2003): Geometric stability condition, Appel¨

, Hagstrom and Kreiss (2006),

Diaz and Joly (2006), Halpern et al. (2011), Duru and Kreiss (2012), ..., Duru (2016) Skelton, et al. (2007): Destabilizing effects of boundary conditions

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications14

SLIDE 21

Summary of literature: IBVP

Assume constant coefficients and consider:

1. IVP PML : (x, y) ∈ (−∞, ∞) × (−∞, ∞). Geometric stability condition B´

ecache et al.

2003. Classical Fourier analysis

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications15

SLIDE 22

Summary of literature: IBVP

Assume constant coefficients and consider:

1. IVP PML : (x, y) ∈ (−∞, ∞) × (−∞, ∞). Geometric stability condition B´

ecache et al.

2003. Classical Fourier analysis
2. IBVP PML : PML is asymmetric.

2a. Lower half–plane PML problem: (x, y) ∈ (−∞, ∞) × (−∞, 0), with LyU = 0, y = 0. 2b . Left half–plane PML problem: (x, y) ∈ (−∞, 0) × (−∞, ∞), with LxU = 0, x = 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications15

SLIDE 23

Summary of literature: IBVP

Assume constant coefficients and consider:

1. IVP PML : (x, y) ∈ (−∞, ∞) × (−∞, ∞). Geometric stability condition B´

ecache et al.

2003. Classical Fourier analysis
2. IBVP PML : PML is asymmetric.

2a. Lower half–plane PML problem: (x, y) ∈ (−∞, ∞) × (−∞, 0), with LyU = 0, y = 0. 2b . Left half–plane PML problem: (x, y) ∈ (−∞, 0) × (−∞, ∞), with LxU = 0, x = 0.

Assumption: In the absence of the PML, dx = 0, the IBVPs are stable in the sense of Kreiss. All eigenvalues are in the stable half of the complex plane. What will happen when the PML is present, dx > 0?

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications15

SLIDE 24

Summary of literature: IBVP

Assume constant coefficients and consider:

1. IVP PML : (x, y) ∈ (−∞, ∞) × (−∞, ∞). Geometric stability condition B´

ecache et al.

2003. Classical Fourier analysis
2. IBVP PML : PML is asymmetric.

2a. Lower half–plane PML problem: (x, y) ∈ (−∞, ∞) × (−∞, 0), with LyU = 0, y = 0. 2b . Left half–plane PML problem: (x, y) ∈ (−∞, 0) × (−∞, ∞), with LxU = 0, x = 0.

Assumption: In the absence of the PML, dx = 0, the IBVPs are stable in the sense of Kreiss. All eigenvalues are in the stable half of the complex plane. What will happen when the PML is present, dx > 0?

Theorem

If the geometric stability condition is satisfied then the PML IBVPs with dx ≥ 0 are asymptotically stable. The PML damping dx > 0 will move all eigenvalues further into the stable complex plane. Duru & Kreiss SINUM (2014), Duru SISC (2016), Duru et al. JCP (2016). PML is asymptotically stable ! Too technical to be extended to numerical approximations in multiple dimensions

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications15

SLIDE 25

Need to extend theory

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications16

SLIDE 26

Need to extend theory

1 Energy estimates + useful for designing stable numerical methods. 2 Energy estimate must account for BCs.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications16

SLIDE 27

Need to extend theory

1 Energy estimates + useful for designing stable numerical methods. 2 Energy estimate must account for BCs. 3 PML is asymmetric. Energy estimate in the time-domain technically difficult.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications16

SLIDE 28

Need to extend theory

1 Energy estimates + useful for designing stable numerical methods. 2 Energy estimate must account for BCs. 3 PML is asymmetric. Energy estimate in the time-domain technically difficult. 4 Energy estimate in the Laplace space. 5 Numerical boundary/inter-element procedures.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications16

SLIDE 29

Energy estimate in the Laplace space

Forcing:

fp, f, fσ, fψ

T with f = (fx, fy, fz).

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications17

SLIDE 30

Energy estimate in the Laplace space

Forcing:

fp, f, fσ, fψ

T with f = (fx, fy, fz). Laplace transform + Eliminate PML auxiliary variables 1 κ s p + ∇d · v = 1 κ

Fp,

ρs v + ∇d p = ρ f,

fη = 1

Sη

fη,
Fp =

1 Sx

fp −

κ sSySx

fσ −

κ sSzSx

fψ
,

BCs: 1 − rη 2 Z vη ∓ 1 + rη 2

p = 0,

at η = ±1.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications17

SLIDE 31

Energy estimate in the Laplace space

Forcing:

fp, f, fσ, fψ

T with f = (fx, fy, fz). Laplace transform + Eliminate PML auxiliary variables 1 κ s p + ∇d · v = 1 κ

Fp,

ρs v + ∇d p = ρ f,

fη = 1

Sη

fη,
Fp =

1 Sx

fp −

κ sSySx

fσ −

κ sSzSx

fψ
,

BCs: 1 − rη 2 Z vη ∓ 1 + rη 2

p = 0,

at η = ±1. 1 κ s2 p − ∇d · 1 ρ ∇d p

= s

κ

Fp − ∇d ·

f, BCs: 1 + rη 2Z s p ± 1 − rη 2 1 Sη ∂ p ∂η = 0, at η = ±1, s∗ Sηρ

p∗ ∂

p ∂η ≥ 0, η = −1; s∗ Sηρ

p∗ ∂

p ∂η ≤ 0, η = 1.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications17

SLIDE 32

Energy estimate in the Laplace space

Re (sSη)∗ Sη

= a + ǫη(s, dη),

ǫη(s, dη) := 2dηb2 |sSη|2 ≥ 0. Define the energy

E2

p (s, dη) =

s

p

2

1/κ +

η=x,y,z
1

Sη ∂ p ∂η

2

1/ρ > 0,

E2

f (s, dη) =

s

Fp

2

1/κ +

η=x,y,z
1

Sη ∂ fη ∂η

2

κ > 0.

Theorem

Consider the PML IBVP with s = 0, Re{s} = a > 0 and piecewise constant dη ≥ 0.

a E2

p (s, dη) +

η=x,y,z
1

Sη ∂ p ∂η

2

ǫη/ρ + BT (s, dη) ≤ 2

Ep (s, dη) Ef (s, dη), BT (s, dη) = 1

−1

1

−1

η=x,y,z

s∗ Sηρ

p∗ ∂

p ∂η

η=−1 −

s∗ Sηρ

p∗ ∂

p ∂η

η=1

dxdydz dη ≥ 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications18

SLIDE 33

Well-posedness & Stability

Theorem

The problem, with dη ≥ 0, is asymptotically stable in the sense that no exponentially growing solutions are supported. Proof: Make the ansatz Q(x, y, z, t) = est Q(x, y, z), with s = a + ib, a E2

p (s, dη) +

η=x,y,z
1

Sη ∂ p ∂η

2

ǫη/ρ + BT (s, dη) = 0.

With Re{s} = a > 0 the expression in the left hand side is positive. The conclusion is that

Q(x, y, z) ≡ 0.
Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications19

SLIDE 34

Well-posedness & Stability

Theorem

The problem, with dη ≥ 0, is asymptotically stable in the sense that no exponentially growing solutions are supported. Proof: Make the ansatz Q(x, y, z, t) = est Q(x, y, z), with s = a + ib, a E2

p (s, dη) +

η=x,y,z
1

Sη ∂ p ∂η

2

ǫη/ρ + BT (s, dη) = 0.

With Re{s} = a > 0 the expression in the left hand side is positive. The conclusion is that

Q(x, y, z) ≡ 0.
Theorem

Let the energy norms E2

p (t, dη) > 0, E2 f (t, dη) > 0. For any a > 0 and T > 0

T e−2atE2

p (t, dη) dt ≤ 4

a2 T e−2atE2

f (t, dη) dt,

a > 0. Key: Construct a discrete approximation that, as far as possible, mimics the energy estimate.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications19

SLIDE 35

Dicretize domain

Discretize the domain (x, y, z) ∈ Ω = ∪Ωlmn with Ωlmn = [xl, xl+1] × [ym, ym+1] × [zn, zn+1], (x, y, z) ← → (q, r, s) ∈ Ω = [−1, 1]3:

x = xl + ∆xl 2 (1 + q) , y = ym + ∆ym 2 (1 + r) , z = zn + ∆zn 2 (1 + s) ,

with

J = ∆xl 2 ∆ym 2 ∆zn 2 > 0, ∆xl = xl+1 − xl, ∆ym = ym+1 − ym, ∆zn = yn+1 − yn.

Ω

f(x, y, z)dxdydz =

K

k=1

L

l=1

M

m=1
Ωklm

f(x, y, z)dxdydz =

K

k=1

L

l=1

M

m=1
Ω

f(q, r, s)Jdqdrds.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications20

SLIDE 36

Dicretize domain

Discretize the domain (x, y, z) ∈ Ω = ∪Ωlmn with Ωlmn = [xl, xl+1] × [ym, ym+1] × [zn, zn+1], (x, y, z) ← → (q, r, s) ∈ Ω = [−1, 1]3:

x = xl + ∆xl 2 (1 + q) , y = ym + ∆ym 2 (1 + r) , z = zn + ∆zn 2 (1 + s) ,

with

J = ∆xl 2 ∆ym 2 ∆zn 2 > 0, ∆xl = xl+1 − xl, ∆ym = ym+1 − ym, ∆zn = yn+1 − yn.

Ω

f(x, y, z)dxdydz =

K

k=1

L

l=1

M

m=1
Ωklm

f(x, y, z)dxdydz =

K

k=1

L

l=1

M

m=1
Ω

f(q, r, s)Jdqdrds. Physics based flux fluctuations:

F x := Z 2

vx −

vx

+ 1

2

p −

p

= 0,

x = xl, Gx := Z 2

vx −

vx

− 1

2

p −

p

= 0,

x = xl+1. (1)

vx,

p encode the solution of the IBVP at element boundaries.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications20

SLIDE 37

Integral form

Ω

φp 1 κ ∂p ∂t + dxp

+ ∇ · v − σ − ψ
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φp Z F η|η=−1 − φp Z Gη|η=1 dqdrds dη ,

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications21

SLIDE 38

Integral form

Ω

φp 1 κ ∂p ∂t + dxp

+ ∇ · v − σ − ψ
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φp Z F η|η=−1 − φp Z Gη|η=1 dqdrds dη ,

Ω

θT

ρ

∂v ∂t + dv

+ ∇p
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

θT nF η|η=−1 + θT nGη|η=1

dqdrds dη ,

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications21

SLIDE 39

Integral form

Ω

φp 1 κ ∂p ∂t + dxp

+ ∇ · v − σ − ψ
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φp Z F η|η=−1 − φp Z Gη|η=1 dqdrds dη ,

Ω

θT

ρ

∂v ∂t + dv

+ ∇p
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

θT nF η|η=−1 + θT nGη|η=1

dqdrds dη ,

Ω

φσ ∂σ ∂t + dyσ + (dy − dx) 2 ∆y ∂vy ∂r

Jdqdrds

= − ωy (dy − dx)

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φσ Z nyF η|η=−1 − φσ Z nyGη|η=1 dqdrds dη

PML stabilizing flux fluctuation

,

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications21

SLIDE 40

Integral form

Ω

φp 1 κ ∂p ∂t + dxp

+ ∇ · v − σ − ψ
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φp Z F η|η=−1 − φp Z Gη|η=1 dqdrds dη ,

Ω

θT

ρ

∂v ∂t + dv

+ ∇p
Jdqdrds

= −

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

θT nF η|η=−1 + θT nGη|η=1

dqdrds dη ,

Ω

φσ ∂σ ∂t + dyσ + (dy − dx) 2 ∆y ∂vy ∂r

Jdqdrds

= − ωy (dy − dx)

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φσ Z nyF η|η=−1 − φσ Z nyGη|η=1 dqdrds dη

PML stabilizing flux fluctuation

,

Ω

φψ ∂ψ ∂t + dzψ + (dz − dx) 2 ∆z ∂vz ∂s

Jdqdrds

= − ωz (dz − dx)

η=q,r,s

1

−1

1

−1

J

η2

x + η2 y + η2 z

φψ Z nzF η|η=−1 − φψ Z nzGη|η=1 dqdrds dη

PML stabilizing flux fluctuation

.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications21

SLIDE 41

The dG Approximation

Galerkin approximation + Polynomial interpolants + GL, GLL, GLR, Collocation nodes:

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications22

SLIDE 42

The dG Approximation

Galerkin approximation + Polynomial interpolants + GL, GLL, GLR, Collocation nodes:

κ−1 dp(t) dt + dxp(t)

+ ∇D · v(t) + σ(t) + ψ(t) = −
η=x,y,z

H−1

η

eη(−1) Z F η − eη(1) Z Gη

,

ρ dv(t) dt + dv(t)

+ ∇Dp(t) = −
η=x,y,z

H−1

η

eη(−1) Z nF η − eη(1) Z nGη

,

dσ(t) dt + dyσ(t)

+ (dy − dx)Dyvy(t) = − ωy (dy − dx)
η=x,y,z

H−1

η

eη(−1) Z nyF η − eη(1) Z nyGη

PML stabilizing flux fluctuation

, dψ(t) dt + dyψ(t)

+ (dz − dx)Dzvz(t) = − ωz (dz − dx)
η=x,y,z

H−1

η

eη(−1) Z nzF η − eη(1) Z nzGη

PML stabilizing flux fluctuation

.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications22

SLIDE 43

The dG Approximation

Galerkin approximation + Polynomial interpolants + GL, GLL, GLR, Collocation nodes:

κ−1 dp(t) dt + dxp(t)

+ ∇D · v(t) + σ(t) + ψ(t) = −
η=x,y,z

H−1

η

eη(−1) Z F η − eη(1) Z Gη

,

ρ dv(t) dt + dv(t)

+ ∇Dp(t) = −
η=x,y,z

H−1

η

eη(−1) Z nF η − eη(1) Z nGη

,

dσ(t) dt + dyσ(t)

+ (dy − dx)Dyvy(t) = − ωy (dy − dx)
η=x,y,z

H−1

η

eη(−1) Z nyF η − eη(1) Z nyGη

PML stabilizing flux fluctuation

, dψ(t) dt + dyψ(t)

+ (dz − dx)Dzvz(t) = − ωz (dz − dx)
η=x,y,z

H−1

η

eη(−1) Z nzF η − eη(1) Z nzGη

PML stabilizing flux fluctuation

. ∇D = (Dx, Dy, Dz)T , Dx = 2 ∆x (D ⊗ I ⊗ I) , Hx = ∆x 2 (H ⊗ I ⊗ I) , ex(η) = (e(η) ⊗ I ⊗ I) , D = H−1A ≈ ∂ ∂q , H = diag[h1, h2, · · · , hP+1], Aij =

P+1

m=1

hmLi(qm)L ′

j (qm) =

1

−1

Li(q)L ′

j (q)dq,

e(η) = [Li(η), Li(η), · · · , LP+1(η)]T .

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications22

SLIDE 44

Stability of the discrete undamped problem

Introduce the elemental energy density

dE (q, r, s, t) = 1 2   1 κ(q, r, s) |p(q, r, s, t)|2 + ρ(q, r, s)

η=x,y,z
|vη(q, r, s, t)|2

 

and the corresponding semi-discrete energy

E (t) =

l=1
m=1
n=1

P+1

i=1

P+1

j=1

P+1

k=1

dE (qi, rj, sk, t)Jhihjhk, J = ∆xl 2 ∆ym 2 ∆zn 2 .

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications23

SLIDE 45

Stability of the discrete undamped problem

Introduce the elemental energy density

dE (q, r, s, t) = 1 2   1 κ(q, r, s) |p(q, r, s, t)|2 + ρ(q, r, s)

η=x,y,z
|vη(q, r, s, t)|2

 

and the corresponding semi-discrete energy

E (t) =

l=1
m=1
n=1

P+1

i=1

P+1

j=1

P+1

k=1

dE (qi, rj, sk, t)Jhihjhk, J = ∆xl 2 ∆ym 2 ∆zn 2 .

Theorem

Consider the semi-discrete approximation. When all the damping vanish, dη = 0, the solution

f the semi-discrete approximation satisfies the energy identity

d dt E (t) = −

L

l=1

M

m=1

N

n=1

 

η=q,r,s

N+1

i=1

N+1

j=1
J
η2

x + η2 y + η2 z

1 Z |Fη|2 + 1 Z |Gη|2 + BT(η)

ij

hi hj  

lmn

≤ 0.

The numerical approximation is asymptotically stable! Result does not translate to the PML when dη > 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications23

SLIDE 46

Stability of the discrete PML problem

Laplace transform + Eliminating the PML auxiliary variables:

1 κ s p + ∇D · v = 1 κ

fp −
η=x,y,z

1 Sη H−1

η

eη(−1) Z F η − eη(1) Z Gη

+
η=y,z

(1 − ωη) (dη − dx) sSηSx H−1

η

eη(−1) Z F η − eη(1) Z Gη

,

ρs v + ∇D p = ρ

f −
η=x,y,z

1 Sη H−1

η

eη(−1) Z nF η − eη(1) Z nGη

,
∇D = (1/SxDx, 1/SyDy, 1/SzDz)T .

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications24

SLIDE 47

Stability of the discrete PML problem

Laplace transform + Eliminating the PML auxiliary variables:

1 κ s p + ∇D · v = 1 κ

fp −
η=x,y,z

1 Sη H−1

η

eη(−1) Z F η − eη(1) Z Gη

+
η=y,z

(1 − ωη) (dη − dx) sSηSx H−1

η

eη(−1) Z F η − eη(1) Z Gη

,

ρs v + ∇D p = ρ

f −
η=x,y,z

1 Sη H−1

η

eη(−1) Z nF η − eη(1) Z nGη

,
∇D = (1/SxDx, 1/SyDy, 1/SzDz)T .

Consider a single element and introduce the modified discrete operators

Dη =

1 Sη

Dη + 1 + rη

2 H−1

η

(Bη (−1, −1) − Bη (1, 1))

,
Hη = H
I + (1 − rη)c

2sSη H−1

η

(Bη(−1, −1) + Bη(1, 1)) −1 .

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications24

SLIDE 48

Stability of the discrete PML problem

Laplace transform + Eliminating the PML auxiliary variables:

1 κ s p + ∇D · v = 1 κ

fp −
η=x,y,z

1 Sη H−1

η

eη(−1) Z F η − eη(1) Z Gη

+
η=y,z

(1 − ωη) (dη − dx) sSηSx H−1

η

eη(−1) Z F η − eη(1) Z Gη

,

ρs v + ∇D p = ρ

f −
η=x,y,z

1 Sη H−1

η

eη(−1) Z nF η − eη(1) Z nGη

,
∇D = (1/SxDx, 1/SyDy, 1/SzDz)T .

Consider a single element and introduce the modified discrete operators

Dη =

1 Sη

Dη + 1 + rη

2 H−1

η

(Bη (−1, −1) − Bη (1, 1))

,
Hη = H
I + (1 − rη)c

2sSη H−1

η

(Bη(−1, −1) + Bη(1, 1)) −1 .

Eliminate the velocity fields:

s∗Hsκ−1s p +

η

1 Sη

Dη

† (s∗S∗

η)

ρSη

Hη

1 Sη

Dη
p

+ |s|2

η

1 + rη 2ZSη HHη

−1 (Bη(−1, −1) + Bη(1, 1))

p = s∗Hκ−1s Fp −

η

s∗H 1 Sη

D0η
fη.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications24

SLIDE 49

Discrete energy estimate in the Laplace space

Introduce the discrete scalar product

u,

v

H =

v†H u, (2) and the energy norms

E 2

p (s) =

s

p, s p

H/κ +
η=x,y,z

1 Sη

Dη

p, 1 Sη

Dη

p

Hη/ρ > 0,

(3)

E 2

f (s) =

s

Fp, s Fp

H/κ +
η=x,y,z

1 Sη

D0η
fη, 1

Sη

D0η
fη
κH > 0.

(4)

Theorem

Consider the one element dG approximation of the PML in the Laplace space with Re{s} = a > 0 and constant damping dη ≥ 0. If ωη = 1, then we have

a E 2

p (s) +

η=x,y,z

1 Sη

Dη

p, 1 Sη

Dη

p

ǫη

Hη/ρ +

η

1 − rη ρ BT(η)

num + BT(s) ≤ 2

Ep(s) Ef (s), BT(s) = |s|2

η

Re 1 Sη 1 + rη 2Z

p†

HH−1

η

(Bη(−1, −1) + Bη(1, 1))

p ≥ 0.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications25

SLIDE 50

2D: Stable approximations

60
40
20

20 40 60 x[km] 10 20 30 40 50 y[km] t=500 s 1 2 3

×10 -5

GLL.

60
40
20

20 40 60 x[km] 10 20 30 40 50 y[km] t=500 s 1 2 3

×10 -5

GL.

60
40
20

20 40 60 x[km] 10 20 30 40 50 y[km] t=500 s 1 2 3

×10 -5

GLR.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications26

SLIDE 51

2D: Stable approximations

100 200 300 400 500 t[s] 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 L∞-norm

ω = 0 ω = 1

GLL.

100 200 300 400 500 t[s] 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 L∞-norm

ω = 0 ω = 1

GL.

100 200 300 400 500 t[s] 10 -4 10 -3 10 -2 10 -1 10 0 10 1 10 2 10 3 10 4 10 5 L∞-norm

ω = 0 ω = 1

GLR.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications27

SLIDE 52

h- and p-convergence

2 4 6 8 10 ∆x 10 -10 10 -8 10 -6 10 -4 10 -2 error GLL GL GLR

C∆x5

h-convergence.

2 4 6 8 polynomial degree 10 -10 10 -8 10 -6 10 -4 10 -2 error GLL GL GLR

p-convergence. dx (x) =    if |x| ≤ 50 km, d0 |x|−50

δ

3 if |x| ≥ 50 km, (5) d0 = 4c 2δ ln 1 tol , tol = C0 1 δ ∆x P + 1 P+1 . (6)

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications28

SLIDE 53

Stable PML boundaries

1 2 3 4 X Axis 4 Y Axis 4 Y Axis 1 2 3 4 Z Axis 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Z Axis 1 2 3 4 X Axis

t = 1.7 s

1 2 3 4 X Axis 4 Y Axis 4 Y Axis 1 2 3 4 Z Axis 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Z Axis 1 2 3 4 X Axis

t = 2.5 s

1 2 3 4 X Axis 4 Y Axis 4 Y Axis 1 2 3 4 Z Axis 0.5 1 1.5 2 2.5 3 3.5 4 4.5 Z Axis 1 2 3 4 X Axis

t = 3.0 s

2 4 6 8 10 t[s]

0.4
0.2

0.2 0.4 p[MPa]

∆x = 5/9 ∆x = 5/27 Analytical

Absorbing boundary condition

2 4 6 8 10 t[s]

0.4
0.2

0.2 0.4 p[MPa]

∆x = 5/9 ∆x = 5/27 Analytical

PML

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications29

SLIDE 54

Extensions to linear elasticity

1 2 3

time [s]

0.6
0.4
0.2

0.2 0.4 0.6

velocity [m/s]

PML ABC analytical Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications30

SLIDE 55

LOH1: A Seismology benchmark problem

time [s] 15 10 5 5 10 15 vx [m/s] analytical ABC PML 15 10 5 5 10 15 vy [m/s] 1 2 3 4 5 6 7 8 9 time [s] 7.5 5.0 2.5 0.0 2.5 5.0 7.5 vz [m/s]

Receiver 4

time [s] 4 2 2 4 vx [m/s] analytical ABC PML 4 2 2 4 vy [m/s] 1 2 3 4 5 6 7 8 9 time [s] 2 1 1 2 vz [m/s]

Receiver 5

time [s] 0.5 0.0 0.5 vx [m/s] analytical ABC PML 0.5 0.0 0.5 vy [m/s] 1 2 3 4 5 6 7 8 9 time [s] 1.0 0.5 0.0 0.5 1.0 vz [m/s]

Receiver 6

time [s] 1.5 1.0 0.5 0.0 0.5 1.0 1.5 vx [m/s] analytical ABC PML 1.5 1.0 0.5 0.0 0.5 1.0 1.5 vy [m/s] 1 2 3 4 5 6 7 8 9 time [s] 1.0 0.5 0.0 0.5 1.0 vz [m/s]

Receiver 9

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications31

SLIDE 56

ExaHyPE: Exa-scale Hyperbolic PDE simulation Engine

Exascale Spacetree ADER-DG

TUM TRE DUR Seismic Astro

Software for next generation super computers (1018 flops/sec), Exascale HPC scalability, Energy efficiency, etc.

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications32

SLIDE 57

Summary

PML for acoustics equation is well-posed and asymptotically stable. Numerical flux procedures can introduce instability. Stable numerical flux procedures can be designed by mimicking continuous energy estimate. Ideas has been extended to linear elasticity. Initial ideas (2D SBP FDM): K. Duru SIAM J. Sci. Comput., 38(2016), A1171-A1194.

K. Duru, A.-A. Gabriel and G. Kreiss, Computer Methods in Applied Mechanics and

Engineering, 350(2019), 898–937.

K. Duru, et al., Extensions to linear elastodynamics, Submitted to Numerische Mathematik

(2019).

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications33

SLIDE 58

Support & Acknowledgements

Prof. Gunilla Kreiss, Uppsala University Sweden.
Prof. Eric M. Dunham, Stanford University, CA.
Dr. Alice-Agnes Gabriel, LMU Munich.

ExaHyPE team: Heinz-Otto Kreiss

Kenneth Duru: On Energy Stable dG Approximation of the PML for the Wave Equation— Monash Workshop on Numerical Differential Equations and Applications34