[PPT] - Towards energy-stable DGSEM for Einsteins equations of general PowerPoint Presentation

SLIDE 1

Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form

ICERM Workshop, Brown University October 5, 2020 Kenneth Duru

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 1 / 37

SLIDE 2

The Einsteins equation

The harmonic coordinates xµ ≡

t, xi

= (t, x, y, z) gxµ := 1 √−g ∂a

−
−ggab∂bxµ

= 0, where g = det

gab

, gab is the inverse metric and g denotes the wave operator. In these coordinates Einstein’s equations reduce to 10 quasilinear wave equations for the metric tensor of the form ggab = Sab where Sab contains the non-linear terms that do not enter the principal part of the equation, and depends only on first derivatives of the metric. Now, the scalar wave equation given by gu := 1 √−g ∂a

−
−ggab∂bu
= 0,

with u a scalar field, has the same principal part as (2.2) and it therefore represent a fundamental model that allows us to test the numerical algorithms for the full non-linear gravitational problem.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 2 / 37

SLIDE 3

shifted wave equation

We consider the 1D shifted wave equation ∂ ∂t ∂u ∂t − a ∂u ∂x

− ∂

∂x

a

∂u ∂t − a ∂u ∂x

+ b ∂u

∂x

= 0,

−1 ≤ x ≤ 1, t ≥ 0. Here a, b, ∈ R, with b > 0.

B. Szil´

agyi, H.-O. Kreiss, and J. Winicour Phys. Rev. D 71, 104035

K. Mattsson and F

. Parisi / Commun. Comput. Phys., 7 (2010), pp. 103-137 a ≡ 0: scalar wave equation (acoustics)

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 3 / 37

SLIDE 4

shifted wave equation

We consider the 1D shifted wave equation ∂ ∂t ∂u ∂t − a ∂u ∂x

− ∂

∂x

a

∂u ∂t − a ∂u ∂x

+ b ∂u

∂x

= 0,

−1 ≤ x ≤ 1, t ≥ 0. Here a, b, ∈ R, with b > 0.

B. Szil´

agyi, H.-O. Kreiss, and J. Winicour Phys. Rev. D 71, 104035

K. Mattsson and F

. Parisi / Commun. Comput. Phys., 7 (2010), pp. 103-137 a ≡ 0: scalar wave equation (acoustics) c = b − a2 Solutions exhibits different characters when c > 0 and c < 0. ∂u ∂t − λ1 ∂ ∂x ∂u ∂t − λ2 ∂ ∂x

u = 0.

λ1 = a + √ b, λ1 = a − √ b, c > 0 : λ1 > 0, λ2 < 0, c < 0 : a > 0 : λ1 > 0, λ2 > 0, a < 0 : λ1 < 0, λ2 < 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 3 / 37

SLIDE 5

Energy estimate

Introduce the energy E(t) =

∂u

∂t − a ∂u ∂x

2

+ b

∂u

∂x

2

We have d dt E(t) ≤ α0E(t), α0 = max

| ∂a

∂x | + | a b ∂b ∂x | + | 1 b ∂b ∂t |

The constant coefficients Cauchy problem satisfies the energy equation

d dt E(t) = 0. We will consider boundary conditions such that d dt E(t) = BTs ≤ 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 4 / 37

SLIDE 6

Stable discrete approximation ?

Given: Wave equation + Initial Conditions + Boundary Conditions. Method of lines –discrete in space and continuous in time. Summation-By-Parts spectral difference operators (in space). Weak enforcement of boundary and inter-element conditions. Stability by energy methods (stable systems of ODEs). Fully discrete: Explicit Runge-Kutta.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 5 / 37

SLIDE 7

Summation-By-Parts (SBP) Vs Integration-By-Parts (IBP)

Continuous: ∂u/∂x ∂u ∂x , v

= −
u, ∂v

∂x

+ v(1)u(1) − v(−1)u(−1),

{IBP}. Discrete: Dxu ≈ ∂u/∂x Dxu, vH = −u, DxvH + vNuN − v1u1, {SBP}. Dx = H−1Q, Q + QT = B = diag ([−1, 0, 0, · · · , 0, 1]) , H = HT > 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 6 / 37

SLIDE 8

Summation-By-Parts (SBP) Vs Integration-By-Parts (IBP)

Continuous: ∂u/∂x ∂u ∂x , v

= −
u, ∂v

∂x

+ v(1)u(1) − v(−1)u(−1),

{IBP}. Discrete: Dxu ≈ ∂u/∂x Dxu, vH = −u, DxvH + vNuN − v1u1, {SBP}. Dx = H−1Q, Q + QT = B = diag ([−1, 0, 0, · · · , 0, 1]) , H = HT > 0. Continuous:

∂ ∂x

b(x) ∂u

∂x

∂

∂x

b(x) ∂u

∂x

, v
= −
b(x) ∂u

∂x , ∂v ∂x

+ v(1)b(1) ∂u

∂x (1) − v(−1)b(−1) ∂u ∂x (−1) {IBP}.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 6 / 37

SLIDE 9

Summation-By-Parts (SBP) Vs Integration-By-Parts (IBP)

Continuous: ∂u/∂x ∂u ∂x , v

= −
u, ∂v

∂x

+ v(1)u(1) − v(−1)u(−1),

{IBP}. Discrete: Dxu ≈ ∂u/∂x Dxu, vH = −u, DxvH + vNuN − v1u1, {SBP}. Dx = H−1Q, Q + QT = B = diag ([−1, 0, 0, · · · , 0, 1]) , H = HT > 0. Continuous:

∂ ∂x

b(x) ∂u

∂x

∂

∂x

b(x) ∂u

∂x

, v
= −
b(x) ∂u

∂x , ∂v ∂x

+ v(1)b(1) ∂u

∂x (1) − v(−1)b(−1) ∂u ∂x (−1) {IBP}. Discrete: D(b)

xx u ≈ ∂/∂x (b(x)∂u/∂x)

D(b)

xx u, vP = −vT M(b) x

u + vj=N(bSxu)j=N − vj=1(bSxu)j=1, uT M(b)

x

u ≥ 0, {SBP}.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 6 / 37

SLIDE 10

Summation-By-Parts (SBP) Vs Integration-By-Parts (IBP)

Continuous: ∂u/∂x ∂u ∂x , v

= −
u, ∂v

∂x

+ v(1)u(1) − v(−1)u(−1),

{IBP}. Discrete: Dxu ≈ ∂u/∂x Dxu, vH = −u, DxvH + vNuN − v1u1, {SBP}. Dx = H−1Q, Q + QT = B = diag ([−1, 0, 0, · · · , 0, 1]) , H = HT > 0. Continuous:

∂ ∂x

b(x) ∂u

∂x

∂

∂x

b(x) ∂u

∂x

, v
= −
b(x) ∂u

∂x , ∂v ∂x

+ v(1)b(1) ∂u

∂x (1) − v(−1)b(−1) ∂u ∂x (−1) {IBP}. Discrete: D(b)

xx u ≈ ∂/∂x (b(x)∂u/∂x)

D(b)

xx u, vP = −vT M(b) x

u + vj=N(bSxu)j=N − vj=1(bSxu)j=1, uT M(b)

x

u ≥ 0, {SBP}. vT M(b)

x

u ≈

b(x) ∂u

∂x , ∂v ∂x

,

(bSxv)j=N ≈ b(1) ∂u ∂x (1), (bSxv)j=1 ≈ b(−1) ∂u ∂x (−1). D(b)

xx = P−1

−M(b)

x

+ BSx

,

B = diag (−1, 0, 0, · · · , 0, 1) , P = PT > 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 6 / 37

SLIDE 11

Fully compatible SBP operators

Diagonal norms H (r = 1, 2, 3, 4, 5) Interior: 2r + boundary: r .

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 7 / 37

SLIDE 12

Fully compatible SBP operators

Diagonal norms H (r = 1, 2, 3, 4, 5) Interior: 2r + boundary: r .

Definition

Let Dxand D(b)

xx , defined above, denote SBP operators approximating ∂/∂x,

∂/∂x (b(x)∂/∂x), with 2r interior order of accuracy and r accurate boundary closures. Let the corresponding diagonal norm be denoted by H, P. The operators compatible SBP

perators if

P = Hx, M(b)

x

= DT

x HbDx + R(b) x ,

R(b)

x

=

R(b)

x

T , vT R(b)

x v ≥ 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 7 / 37

SLIDE 13

Fully compatible SBP operators

Diagonal norms H (r = 1, 2, 3, 4, 5) Interior: 2r + boundary: r .

Definition

Let Dxand D(b)

xx , defined above, denote SBP operators approximating ∂/∂x,

∂/∂x (b(x)∂/∂x), with 2r interior order of accuracy and r accurate boundary closures. Let the corresponding diagonal norm be denoted by H, P. The operators compatible SBP

perators if

P = Hx, M(b)

x

= DT

x HbDx + R(b) x ,

R(b)

x

=

R(b)

x

T , vT R(b)

x v ≥ 0

Definition

The operators fully compatible SBP operators if P = Hx, M(b)

x

= DT

x HbDx + R(b) x ,

R(b)

x

=

R(b)

x

T , vT R(b)

x v ≥ 0

Sx = Dx.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 7 / 37

SLIDE 14

Geometric flexibility: Duru and Virta JCP (2014)

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 8 / 37

SLIDE 15

Some remarks

Fully-compatible operators not sufficient deriving stable numerical methods for 1D shifted wave equation

K. Mattsson and F

. Parisi / Commun. Comput. Phys., 7 (2010), pp. 103-137 Even for the Cauchy problem

B. Szil´

agyi, H.-O. Kreiss, and J. Winicour Phys. Rev. D 71, 104035

Definition

The operators ultra-compatible SBP operators if P = Hx, M(b)

x

= DT

x HbDx + R(b) x ,

R(b)

x

≡ 0, Sx = Dx. A straightforward approach (use first derivative twice): ∂ ∂x

b(x) ∂u

∂x

≈ DxbDxu.

Yielded unsatisfactory results when c < 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University 9 / 37

SLIDE 16

DG operators

Consider the polynomial approximation of u(x) u(x) ≈ U(x) =

P+1

j=1

φj(x)Uj where Uj are degrees of freedom to be determined. 1

−1

φ(x) dU(x) dx dx =

P+1

j=1

Uj 1

−1

φ(x)φ′

j (x)dx.

φj(x) = Lj(x) – Lagrange polynomials = ⇒ Uj = u(xj) Galerkin approximation + Collocation + Quadrature 1

−1

Li(x) dU(x) dx dx =

P+1

j=1

QijUj = QU, The quadrature rule is exact for all polynomial integrand of degree ≤ 2P − 1! Qij =

N+1

m=1

wmLi(xm)L ′

j (xm) =

1

−1

Li(x)L ′

j (x)dx,

(1)

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University10 / 37

SLIDE 17

Q + QT = e(1)e(1)T − e(−1)e(−1)T , e(η) = [L1(η), L2(η), · · · , LP+1(η)]T . The matrices H, D ∈ R(P+1)×(P+1) are defined by D = H−1Q ≈ ∂ ∂x , H = ∆x 2 diag([w1, w2, · · · wP+1]). (2) D is a spectral difference approximation of the first derivative in one space dimension. ∂ ∂x

b(x) ∂u

∂x

≈ DbDU.

Note that DbD = H−1    −DT HbD +

e(1)e(1)T − e(−1)e(−1)T
B

D     .

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University11 / 37

SLIDE 18

Case 1: c > 0

Let us consider c > 0, this implies λ1 > 0 and λ2 < 0 for any a ∈ R. We impose the boundary conditions ∂u ∂t − λ1 ∂u ∂x

= 0,

x = −1 ∂u ∂t − λ2 ∂u ∂x

= 0,

x = 1 Define the energy E(t) =

∂u

∂t − a ∂u ∂x

2

+

√

b ∂u ∂x

2

. We have d dt E(t) = λ2 2 ∂u ∂t − λ1 ∂u ∂x 2

1 − λ1

2 ∂u ∂t − λ2 ∂u ∂x 2

−1 ≤ 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University12 / 37

SLIDE 19

Case 1: c > 0

d dt du dt − aDu

− D
a

du dt − aDu

+ bDu
+ τNH−1e(1)eT (1)

du dt − λ2Du

+ τ0H−1e(−1)eT (−1)

du dt − λ1Du

= 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University13 / 37

SLIDE 20

Case 1: c > 0

d dt du dt − aDu

− D
a

du dt − aDu

+ bDu
+ τNH−1e(1)eT (1)

du dt − λ2Du

+ τ0H−1e(−1)eT (−1)

du dt − λ1Du

= 0

Introducing w such that du dt = aDu + w and d dt w = aDw − DT HT bDu + H−1 e(1)eT (1) − e(−1)eT (−1)

bDu

(3) − τNH−1e(1)eT (1) du dt − λ1Du

− τ0H−1e(−1)eT (−1)

du dt − λ2Du

= 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University13 / 37

SLIDE 21

Proof

Define A = DT bHD and the discrete energy E (t) =

w
2

H + uT Au > 0.

(4) We will prove τN = λ1 and τ0 = −λ2, gives d dt E (t) ≤ 0. (5)

d dt E (t) = wT Hwt + wT

t Hw + uT Aut + uT t Au

≤ − (τN − λ1) 2 du dt − λ2Du

−

τN λ1 − τN du dt − λ1Du 2 − τ2

N

2(λ1 − τN) du dt − λ1Du 2 − (τ0 + λ2) 2 du dt − λ1Du

−

τ0 λ1 − τ0 du dt − λ2Du 2 − τ2 2(λ2 + τ0) du dt − λ2Du 2 ≤ 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University14 / 37

SLIDE 22

Eigenvalues : a = 1, b = 2, c = 3 > 0

10
8
6
4
2

real

15
10
5

5 10 15 imaginary

P = 11

100
80
60
40
20

real

80
60
40
20

20 40 60 80 imaginary

P = 51

400
350
300
250
200
150
100
50

real

250
200
150
100
50

50 100 150 200 250 imaginary

P = 101

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University15 / 37

SLIDE 23

Case 2: c < 0

Let is us consider a > 0 and c < 0, this implies λ1 > λ2 > 0. We will need to impose 2 boundary conditions at x = 1, and there is no boundary conditions at x = −1. If a < 0 and c < 0, then λ1 < λ2 < 0, and the conditions will be reversed. We impose the boundary conditions u(x, t) = 0, ∂u ∂t − a ∂u ∂x = 0, x = 1

d dt    du dt − a   D − γ0 H−1e(1)eT (1)

u=0

   u    − D   a    du dt − a   D − γ0 H−1e(1)eT (1)

u=0

   u    + bDu    + τNH−1e(1)eT (1)a    du dt − a   D − γ0 H−1e(1)eT (1)

u=0

   u   

∂u

∂t −a ∂u ∂x =0

+ H−1 γNDT be(1)eT (1) + τ0e(1)eT (1)bH−1e(1)eT (1)

u
u=0

= 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University16 / 37

SLIDE 24

Case 2: c < 0

We set γ0 = τ0 = τN = 1 and γN = −1

d dt    du dt − a   D − H−1e(1)eT (1)

u=0

   u    −

D − H−1e(1)eT (1)


 a    du dt − a   D − H−1e(1)eT (1)

u=0

   u       − DbDu − H−1 DT be(1)eT (1) − e(1)eT (1)bH−1e(1)eT (1)

u
u=0

= 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University17 / 37

SLIDE 25

Case 2: c < 0

Introducing w such that du dt = a

D − H−1e(1)eT (1)
u + w

and d dt w = a

D − H−1e(1)eT (1)
w − H−1

D − H−1e(1)eT (1) T Hb

D − H−1e(1)eT (1)
u

(6) Define

D =
D − H−1e(1)eT (1)
,
A =

DT bH D Note that uT

H

D

+
H

D T u = −uT e(−1)eT (−1) + e(1)eT (1)

u = −u2(−1) − u2(1) ≤ 0.

We can define the discrete energy E (t) =

w
2

H + uT

Au > 0. (7)

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University18 / 37

SLIDE 26

Proof

We will prove d dt E (t) ≤ 0. (8) d dt E (t) = wT Hwt + wT

t Hw + uT

AUt + uT

t

Au = awT

H

D T +

H

D

w − wT

Au − uT Aw + auT A Du + uT Aw + a

Du

T Au + wT Au = −a(w2(−1) + w2(1)) + ab

Du

T H D

+
H

D T

Du
= −a(w2(−1) + w2(1)) − ab
Du

T e(−1)eT (−1) + e(1)eT (1)

Du
≤ 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University19 / 37

SLIDE 27

Eigenvalues: a = −2, b = 1, c = −3 < 0

35
30
25
20
15
10
5

real

60
40
20

20 40 60 imaginary

P = 11

700
600
500
400
300
200
100

real

1500
1000
500

500 1000 1500 imaginary

P = 51

3000
2500
2000
1500
1000
500

real

5000
4000
3000
2000
1000

1000 2000 3000 4000 5000 imaginary

P = 101

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University20 / 37

SLIDE 28

Interface treatments

The goal is to derive a conservative discontinuous formulation of the scalar shifted wave equation to patch elements together. Let u+ denote the solution in the positive real line x ≥ 0 and u− denote the solution in the negative real line x ≤ 0. ∂ ∂t ∂u− ∂t − a ∂u− ∂x

− ∂

∂x

a

∂u− ∂t − a ∂u− ∂x

+ b ∂u−

∂x

= 0,

x < 0, ∂ ∂t ∂u+ ∂t − a ∂u+ ∂x

− ∂

∂x

a

∂u+ ∂t − a ∂u− ∂x

+ b ∂u+

∂x

= 0,

x > 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University21 / 37

SLIDE 29

Interface treatments

The goal is to derive a conservative discontinuous formulation of the scalar shifted wave equation to patch elements together. Let u+ denote the solution in the positive real line x ≥ 0 and u− denote the solution in the negative real line x ≤ 0. ∂ ∂t ∂u− ∂t − a ∂u− ∂x

− ∂

∂x

a

∂u− ∂t − a ∂u− ∂x

+ b ∂u−

∂x

= 0,

x < 0, ∂ ∂t ∂u+ ∂t − a ∂u+ ∂x

− ∂

∂x

a

∂u+ ∂t − a ∂u− ∂x

+ b ∂u+

∂x

= 0,

x > 0. Multiply the equation by a sufficiently smooth function φ(x) and integrate over the domain, we have Collecting contributions from both sides of the interface gives

φ(x) ∂

∂t

∂u−

∂t − a ∂u− ∂x

dx +
∂φ(x)

∂x

a
∂u−

∂t − a ∂u− ∂x

+ b ∂u−

∂x

dx+
φ(x) ∂

∂t

∂u+

∂t − a ∂u+ ∂x

dx +
∂φ(x)

∂x

a
∂u+

∂t − a ∂u+ ∂x

+ b ∂u+

∂x

dx − φ(x)[

[F(u)] ] = 0,

where F(u) is the flux given by F(u) =

a
∂u

∂t − a ∂u ∂x

+ b ∂u

∂x

.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University21 / 37

SLIDE 30

Interface treatments

The flux has two essential components:

1. a
∂u

∂t − a ∂u ∂x

the contribution form the advective transport when a = 0,
2. b ∂u

∂x the contribution from the ”expanding pressure wave”.

Any conserved quantity propagated by the PDE must ensure that the jump in the flux vanishes [ [F(u)] ] = 0. This will also be useful when designing stable and conservative schemes. [ [u] ] := u+ − u− = 0 = ⇒ [ [∂u ∂t ] ] = ∂u+ ∂t − ∂u− ∂t = 0, and [ [∂u ∂x ] ] = ∂u+ ∂x − ∂u− ∂x = 0 (9) d dt

E−(t) + E+(t)
= 0.

(10)

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University22 / 37

SLIDE 31

Interface treatments

Introduce the interface matrix

B =
e(1)eT (1)

−e(1)eT (−1) e(−1)eT (1) −e(−1)eT (−1)

with
v−

v+ T

B
u−

u+

=
v−(1) − v+(−1)

u−(1) − u+(−1)

,
v−

v+ T

B
u−

u+

=
v−(1) + v+(−1)

u−(1) − u+(−1)

.

Note that

BT =
e(1)eT (1)

e(1)eT (−1) −e(−1)eT (1) −e(−1)eT (−1)

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University23 / 37

SLIDE 32

Interface treatments

A consistent discontinuous Galerkin spectral difference approximation of the shifted wave equation is

d dt          

du− dt du+ dt

  − a        

D

D

− γ0
H

H −1

B
continuity of solution

       

u−

u+



       −

D

D



       a          

du− dt du+ dt

  − a        

D

D

− γ0
H

H −1

B
continuity of solution

       

u−

u+



       + b

D

D u− u+



       + τN

H

H −1

B

        a          

du− dt du+ dt

  − a        

D

D

− γ0
H

H −1

B
continuity of solution

       

u−

u+



       + b

D

D u− u+



      

continuity of flux

+

H

H −1  γN

DT b

DT b

B + τ0

BT b

H

H −1

B

 

u−

u+

continuity of solution

= 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University24 / 37

SLIDE 33

Interface treatment

The blue penalty terms are flux corrections due to advective transport, when a = 0, and the remaining penalty terms are flux corrections due to expanding pressure waves. Next, we will determine the penalty parameters such that the scheme is energy stable. We set τN = 1/2, γN = −1/2, γ0 = 1/2 and τ0 = 1/4 obtaining

d dt    

du− dt du+ dt

  − a  

D

D

−

1 2

H

H −1

B

 

u−

u+   −

D

D  a    

du− dt du+ dt

  − a  

D

D

−

1 2

H

H −1

B

 

u−

u+   + b

D

D u− u+   + 1 2

H

H −1

B

 a    

du− dt du+ dt

  − a  

D

D

−

1 2

H

H −1

B

 

u−

u+   + b

D

D u− u+   − 1 2

H

H −1  

DT b

DT b

B −

1 2

BT b
H

H −1

B

 

u−

u+

= 0

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University25 / 37

SLIDE 34

Interface treatment

Introduce

D =

 

D

D

− 1

2

H

H −1

B

  U =

u−

u+

,

W =

w−

w+

,

H =

H

H

d

dt U = a DU + W d dt W = a DW − H−1AU A = DT bH D

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University26 / 37

SLIDE 35

Interface treatment

We will use WT

H

D

+
H

D T W = 0, ∀W = 0, We can define the discrete energy E (t) =

W
2

H + UT AU > 0.

We will prove d dt E (t) = 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University27 / 37

SLIDE 36

Proof

d dt E (t) = WT HWt + WT

t HW + UT AUt + UT t AU

= aWT

H

D T +

H

D

W − WT AU − UT AW

+ aUT A DU + UT AW + a

DU

T AU + WT AU = ab

DU

T H D

+
H

D T

DU
= 0.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University28 / 37

SLIDE 37

Eigenvalues : a = 1, b = 2, c = 1 > 0

Number of elements = 11

30
25
20
15
10
5

real

800
600
400
200

200 400 600 800 imaginary

P = 11

800
700
600
500
400
300
200
100

real

1.5
1
0.5

0.5 1 1.5 imaginary 104

P = 51

3000
2500
2000
1500
1000
500

real

6
4
2

2 4 6 imaginary 104

P = 101

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University29 / 37

SLIDE 38

Eigenvalues: a = −2, b = 1, c = −3 < 0

Number of elements = 11

400
350
300
250
200
150
100
50

real

1500
1000
500

500 1000 1500 imaginary

P = 11

8000
7000
6000
5000
4000
3000
2000
1000

real

3
2
1

1 2 3 imaginary 104

P = 51

3
2.5
2
1.5
1
0.5

real 104

1.5
1
0.5

0.5 1 1.5 imaginary 105

P = 101

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University30 / 37

SLIDE 39

Accuracy

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t 10-12 10-10 10-8 10-6 10-4 10-2 100 error 21-elements 41-elements 81-elements 161-elements 321-elements

GLL.

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t 10-12 10-10 10-8 10-6 10-4 10-2 100 error 21-elements 41-elements 81-elements 161-elements 321-elements

GL.

0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 t 10-12 10-10 10-8 10-6 10-4 10-2 100 error 21-elements 41-elements 81-elements 161-elements 321-elements

GLR.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University31 / 37

SLIDE 40

spectral convergence

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

Figure: Spectral convergence polynomial degrees P = {2, 4, 6, 8, 10}.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University32 / 37

SLIDE 41

Convergence rate

Table: GL.

∆x error rate 0.1 5.3624×10−3

0.05

1.3727×10−4 5.2877 0.025 6.1924×10−6 4.4704 0.0125 2.4207×10−7 4.6770 0.00625 8.2462×10−9 4.8756

Table: GLL.

∆x error rate 0.1 5.2286×10−3

0.05

1.1270×10−4 5.5359 0.025 4.42875×10−7 7.9913 0.0125 4.9542×10−9 6.4821 0.00625 7.0721×10−11 6.1303

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University33 / 37

SLIDE 42

Example 1: a = −1, b = 2, c = 1 > 0

λ1 = −1 + √ 2 > 0, λ2 = −1 − √ 2 < 0 3 elements, P = 101 polynomial degree 10 PPW

5

5 10 15 x

1

1 2 3 4 5 6 7 8 9 amplitude t=0.72553

5

5 10 15 x 1 2 3 4 5 6 7 8 9 10 amplitude t=2.1775

5

5 10 15 x 1 2 3 4 5 6 7 8 9 10 amplitude t=3.6295

5

5 10 15 x 1 2 3 4 5 6 7 8 9 amplitude t=5.0814

5

5 10 15 x

1

1 2 3 4 5 6 7 8 9 amplitude t=10.8893

4
2

2 4 6 8 10 12 14 x

1.5
1
0.5

0.5 1 amplitude 10-4 t=16.6972

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University34 / 37

SLIDE 43

Example 2: a = −2, b = 1, c = −3 < 0

λ1 = −1 < 0, λ2 = −3 < 0 3 elements, P = 101 polynomial degree 10 PPW

5

5 10 15 x

5

5 10 15 amplitude t=0.72553

5

5 10 15 x

5

5 10 15 amplitude t=2.1775

5

5 10 15 x

5

5 10 15 amplitude t=3.6295

5

5 10 15 x

2

2 4 6 8 10 12 14 16 amplitude t=5.0814

5

5 10 15 x

2

2 4 6 8 10 12 14 16 amplitude t=10.8893

5

5 10 15 x

3
2
1

1 2 3 amplitude 10-4 t=16.6972

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University35 / 37

SLIDE 44

61 elements, P = 4 polynomial degree

Keeping the degrees of freedom constant: 10 PPW

5

5 10 15 x

5

5 10 15 amplitude t=0.72553

5

5 10 15 x

5

5 10 15 amplitude t=2.1775

5

5 10 15 x

5

5 10 15 amplitude t=3.6295

5

5 10 15 x 2 4 6 8 10 12 14 amplitude t=5.0814

5

5 10 15 x

2

2 4 6 8 10 12 14 16 amplitude t=10.8893

5

5 10 15 x

0.4
0.3
0.2
0.1

0.1 0.2 0.3 0.4 0.5 amplitude t=16.6972

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University36 / 37

SLIDE 45

Summary and outlook

Initial developments of an efficient and robust numerical method for the Einstein’s equations modelling gravitational waves in second order form. The new method combines the advantages and central idea of two very successful numerical techniques, the summation-by-parts finite difference methods and the discontinuous Galerkin methods. We consider a 1D model problem, the shifted wave in second order form in one space dimension, and derive compatible derivative operators and numerical fluxes for the different components of the PDE flux terms. The numerical flux allows for the derivation of a provably energy-conserving and arbitrarily high-order accurate discontinuous Galerkin spectral element method for the 1D shifted wave in second order form. The method can be extended to higher space dimensions using tensor product elements. Presented the proof of numerical stability and numerical experiments verifying the accuracy and stability properties of the method.

Kenneth Duru: Towards energy-stable DGSEM for Einstein’s equations of general relativity in second order form— ICERM Workshop, Brown University37 / 37