The energy based discontinuous Galerkin method Daniel Appel Applied - - PowerPoint PPT Presentation

The energy based discontinuous Galerkin method

Daniel Appelö Applied Mathematics University of Colorado, Boulder

Outline

ICERM 2018

• Applications / interests
• Why second order?
• Why energy based?
• Methods

Please Ask Questions Throughout

Contributors 
 Thomas Hagstrom, Southern Methodist U. Fatemeh Pourahmadian, CU, Boulder. Olof Runborg, Royal Inst. of Technology. Siyang Wang, Chalmers Inst. of Technology.

ICERM 2018

Micropolar Elasticity

• Classic Hooke’s theory: Continuum description with

balance of force and stress. Symmetric tensors.

• Voigt (1877), Cosserat brothers (1909): Asymmetric theory

with force-stress and couple-stress balance.

• Micropolar theory of Eringen-Nowacki: displacements and

rotations.

• Allows for construction of phononic materials.
• Acoustic properties depends on grain-to-grain 

• forces. These can be controlled using EM fields.

ρvtt = (µ + κ)vxx − κθx, ρJθtt = γθxx + κ(vx − 2θ).

1 2 3 4 5 6 1 2 3 4 5 6

Systems of wave equations (dispersive)

ICERM 2018

Interests in numerical homogenization and beyond

• Developing wave solvers that are accurate and

geometrically flexible (high order dG with upwind fluxes).

• Methods that couple well with micro-models. 


In particular micro models described in terms of (discrete) energy balance (Heterogenous Multiscale Models).

• Methods that allow the macro-models to have anisotropy,

dispersive and non-linear effects.

• Fast solvers that can be used for direct simulation of

Maxwell / acoustic cell problems.

• Variable order methods that can be used together with

multi-level algorithms in UQ, inverse problems and frequency domain solvers.

• Taming of the poor CFL constraints of dG.

ICERM 2018

Motivation: Why solve in second order form?

• “There are three kinds of mathematicians; those who can

count and those who cannot… ” .

• Maxwell in 3D has 12 derivative matrices.
• Formulated as a second order problem there can be as few

as 3 derivative matrices (although for complex materials the saving may be less).

• Many problems in solid mechanics are formulated in terms
• f energy densities - this fits well with our formulation.
• The energy based formulation naturally incorporates

upwind fluxes for second order formulations.

• The method is parameter free.

ICERM 2018

Energy based discontinuous Galerkin method

Plan:

• What energy are we talking about?
• The general formulation…
• What the general formulation really is…
• SIPDG.
• EDG with central flux / Baumann-Oden form.
• Full EDG (upwind / alt. / central flux).
• Extension to SBP-SAT finite differences.

ICERM 2018

What energy are we talking about?

E(t) =

• Ω

1 2 ⏐ ⏐ ⏐ ⏐ ∂u ∂t ⏐ ⏐ ⏐ ⏐

2

+ G(u, ∇u, x).

The energy we consider is a non-negative functional of the kinetic and potential energy density. For example:

ICERM 2018

What energy are we talking about?

E(t) =

• Ω

1 2 ⏐ ⏐ ⏐ ⏐ ∂u ∂t ⏐ ⏐ ⏐ ⏐

2

+ G(u, ∇u, x).

The energy we consider is a non-negative functional of the kinetic and potential energy density. For example: The equations we treat result from the Euler -Lagrange equations derived from the action principle associated with the Lagrangian

n be identified as the n 1

2

⏐ ⏐ ∂u

∂t

⏐ ⏐2 − G − u · f

ICERM 2018

The Euler-Lagrange equations are

∂2ui ∂t2 =

• k

∂ ∂xk ∂G ∂ui,k

• − ∂G

∂ui + fi,

∂t − ∂vi ∂t −

• k

∂ ∂xk ∂G ∂ui,k

• + ∂G

∂ui = fi, ∂ui ∂t − vi = 0,

We prefer to evolve “displacement” and “velocity”

ICERM 2018

The Euler-Lagrange equations are

∂2ui ∂t2 =

• k

∂ ∂xk ∂G ∂ui,k

• − ∂G

∂ui + fi,

d dt

• Ωj

1 2|v|2 + G =

• Ωj

v · f +

• ∂Ωj
• i,k

vi ∂G ∂ui,k nk, .

∂t − ∂vi ∂t −

• k

∂ ∂xk ∂G ∂ui,k

• + ∂G

∂ui = fi, ∂ui ∂t − vi = 0,

We prefer to evolve “displacement” and “velocity” Seek methods that automagically satisfies

ICERM 2018

B0ring testing for velocity

∂t − ∂vi ∂t −

• k

∂ ∂xk ∂G ∂ui,k

• + ∂G

∂ui = fi,

ICERM 2018

Boring testing for velocity

∂t − ∂vi ∂t −

• k

∂ ∂xk ∂G ∂ui,k

• + ∂G

∂ui = fi,

What about the displacement equation?

∂ui ∂t − vi = 0,

ICERM 2018

∂t − ∂vi ∂t −

• k

∂ ∂xk ∂G ∂ui,k

• + ∂G

∂ui = fi,

d dt

• Ωj

1 2|v|2 + G =

• Ωj

v · f +

• ∂Ωj
• i,k

vi ∂G ∂ui,k nk, .

Note that we are halfway there!

Boring testing for velocity

ICERM 2018

Santa’s test

∂ui ∂t − vi = 0,

Test the displacement with the variation of the potential energy

ICERM 2018

Santa’s test

∂ui ∂t − vi = 0,

ICERM 2018

Use solution as test function, integrate by parts, 
 add up to get a Christmas present!

ICERM 2018

Use solution as test function, integrate by parts, 
 add up to get a Christmas present!

Theorem 1. Suppose U h(t) and the fluxes v∗, w∗ are given. Suppose further that (2.2) is linear. Then

dUh dt

satisfying Problem 1 is uniquely determined and the energy identity for Eh(t) =

j Eh j (t)

(2.16) dEh dt =

• i,j
• Ωj

vh

i fi(x, t) +

• i,j
• ∂Ωj
• k

nk

• v∗

i − vh i

∂G ∂ui,k (uh, ∇uh, x) + vh

i w∗ i,k

• ,

is satisfied.

h

Element-wise energy identity

Eh

j (t) =

• Ωj

1 2 ⏐ ⏐vh⏐ ⏐2 + G(uh, ∇uh, x)

Where

ICERM 2018

Accepts upwind and averaged fluxes

For the general parametrization We get element-wise contributions to the energy

ICERM 2018

Accepts upwind and averaged fluxes

For the general parametrization We get element-wise contributions to the energy

v∗

i = vh i,1,

w∗

i,k =

∂G ∂ui,k (uh

2, ∇uh 2, x).

v∗

i = {{vh i }} − ζi

2 [[D∇uiGh]], {{ }} − 2 w∗

i,k = − 1

2ζi [[vh

i ]]k + {{ ∂Gh

∂ui,k }}.

Alternating flux Sommerfeld flux Two attractive choices are:

ICERM 2018

Examples of equations that fit into the framework

σij = λuk,k δij + µ(ui,j + uj,i) + κ

• (ui,j − uj,i)/2 − εkijϕk
• mij = αϕk,kδij + β + γ

2 (ϕi,j + ϕj,i) + β − γ 2 (ϕi,j − ϕj,i)

mji,j + εijkσjk − ρJ ¨ ϕi = 0 σji,j − ρ¨ ui = 0

Et = U, µ0ε0Ut = ∆E − µ0Z, Jt = Z, Zt = −ΓeZ + ε0ω2

peU,

Ht = V, µ0ε0Vt = ∆H − ε0Y, Kt = Y, Yt = −ΓmY + µ0ω2

pmV,

E = 1 2

• ε0µ0
• U2 + V2

+ (∇ × E)2 + (∇ × H)2 + µ0 ε0ω2

pe

Z2 + ε0 µ0ω2

pm

Y2

• .

− G

Micropolar elasticity Maxwell Drude material

ICERM 2018

“Ok this is very general, but how do I implement it?!”

ICERM 2018

Back to basics SIPDG for wave.

utt − ∇ · (c∇u) = f, in J × Ω, u = 0,

• n J × ∂Ω,

(uh

tt, U) + ah(uh, U)

= (f, U) ∀ U ∈ Vh, t ∈ J,

Test as usual, U is the test function.

ICERM 2018

Back to basics SIPDG for wave.

utt − ∇ · (c∇u) = f, in J × Ω, u = 0,

• n J × ∂Ω,

(uh

tt, U) + ah(uh, U)

= (f, U) ∀ U ∈ Vh, t ∈ J,

ah(u, w) =

• K∈Th
• K

c∇u · ∇wdx −

• F ∈Fh
• F

[[u]] · {{c∇w}}dA −

• F ∈Fh
• F

[[w]] · {{c∇u}}dA + γ h

• F ∈Fh
• F

[[u]] · [[w]]dA.

Bilinear form: Test as usual, U is the test function.

ICERM 2018

Back to basics SIPDG for wave.

utt − ∇ · (c∇u) = f, in J × Ω, u = 0,

• n J × ∂Ω,

(uh

tt, U) + ah(uh, U)

= (f, U) ∀ U ∈ Vh, t ∈ J,

ah(u, w) =

• K∈Th
• K

c∇u · ∇wdx −

• F ∈Fh
• F

[[u]] · {{c∇w}}dA −

• F ∈Fh
• F

[[w]] · {{c∇u}}dA + γ h

• F ∈Fh
• F

[[u]] · [[w]]dA.

Test as usual, U is the test function. Bilinear form: I.b.p, symmetry, penalty (stabilization). Energy conserving - optimal convergence - known bf.s.

ICERM 2018

ut − v = 0, in J × Ω, vt − ∇ · (c∇u) = f, in J × Ω, First order in time. (vh

t , V ) + bh(uh, V ) = (f, V ),

U = 0, ∀ V ∈ Vh, t ∈ J, Velocity equation looks as in SIPDG (b is not the same).

Conservative Energy DG

ICERM 2018

ut − v = 0, in J × Ω, vt − ∇ · (c∇u) = f, in J × Ω, First order in time. (vh

t , V ) + bh(uh, V ) = (f, V ),

U = 0, ∀ V ∈ Vh, t ∈ J,

) = 1 2 d dt

• K∈Th
• K

(vh)2 + c∇uh · ∇uhdx = (f, vh).

Velocity equation looks as in SIPDG (b is not the same). We would like this type of estimate! What is missing above?

Conservative Energy DG

ICERM 2018

(vh

t , V ) + bh(uh, V ) = (f, V ),

U = 0, ∀ V ∈ Vh, t ∈ J,

Conservative Energy DG

dh(uh

t , U) − bh(U, vh) = 0,

V = 0, ∀ U ∈ Vh, t ∈ J.

Time for Christmas gifus!

ICERM 2018

(vh

t , V ) + bh(uh, V ) = (f, V ),

U = 0, ∀ V ∈ Vh, t ∈ J,

Conservative Energy DG

dh(uh

t , U) − bh(U, vh) = 0,

V = 0, ∀ U ∈ Vh, t ∈ J.

Time for Christmas gifus!

bilinear forms bh, and dh are dh(u, w) =

• K∈Th
• K

c∇u · ∇wdx, bh(u, w) =

• K∈Th
• K

c∇u · ∇wdx −

• F ∈Fh
• F

{{c∇u}} · [[w]]dA.

With the bilinear forms being.

ICERM 2018

(vh

t , V ) + bh(uh, V ) = (f, V ),

U = 0, ∀ V ∈ Vh, t ∈ J,

Conservative Energy DG

dh(uh

t , U) − bh(U, vh) = 0,

V = 0, ∀ U ∈ Vh, t ∈ J.

Time for Christmas gifus!

bilinear forms bh, and dh are dh(u, w) =

• K∈Th
• K

c∇u · ∇wdx, bh(u, w) =

• K∈Th
• K

c∇u · ∇wdx −

• F ∈Fh
• F

{{c∇u}} · [[w]]dA.

dx = (f, vh). (vh

t , vh) + dh(uh t , uh) =

With the bilinear forms being. And we indeed get the energy identity.

ICERM 2018

The equation for u is

Evolving the displacement

dh(uh

t − vh, U) = −

• F ∈Fh
• F

{{c∇U}} · [[vh]]dA,

dh(u, w) =

• K∈Th
• K

c∇u · ∇wdx,

S(uh

t − vh) = L1vh. The form d is the “stiffness matrix” , which is singular So how should we evolve this?

ICERM 2018

Evolving the displacement

S(uh

t − vh) = L1vh.

uh

t = vh + S†L1vh.

• K∈Th
• K

N(uh

t − vh)dx = 0,

) = 1 2 d dt

• K∈Th
• K

(vh)2 + c∇uh · ∇uhdx

The energy does not see constants. So we add moments of its null space. Equivalently - compute Moore - Penrose inverse

ICERM 2018

Inverting the stiffness

S+ = (S + N)−1 − N.

hat g = (1, 0, . . . , 0)T .

is g = 1

n(1, 1, . . . , 1)T

sy to check that =

0) . If we use t N = ggT Given a choice of basis we know the null vector. Modal. Nodal. Not hard to show that

S†LR = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ . . . 1/2 1/2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ , S†LL = ⎡ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣ . . . −1/2 1/2 ⎤ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦ .

uh

t = vh + S†L1vh.

Sometimes we can do even better! Legendre basis in 1D

ICERM 2018

Upwind fluxes

Add appropriate dissipative but consistent terms

dh(uh

t , U) − dh(vh, U)

= −

• F ∈Fh
• F

{{c∇U}} · [[vh]] + [[∂U ∂n ]][[∂uh ∂n ]]dA, (vh

t , V ) + dh(uh, V )

= +

• F ∈Fh
• F

{{c∇uh}} · [[V ]] − [[V ]][[vh]] dA,

The upwind fluxes yields optimal convergence on general meshes (observed) and in 1D (theorem). Alternating fluxes are easily incorporated by replacing arithmetic averages with weighted averages

ICERM 2018

Conclusion - Energy based DG

• Parameter free stability.
• Naturally introduces upwinding.
• Can be expressed in commonly used bilinear forms for easy

implementation in open source FEM libraries (MFEM, deal.ii, Fenics).

• Implementation in MFEM took one day, 


including learning MFEM…

• The “add and subtract” to get an energy 


estimate is easily extended to SBP-SAT, 
 Galerkin differences, etc.

• Nonlinear problems can be handled as well.
• The u,v formulation lends itself to staggering.