Comparing DG and Nedelec finite element discretisations of the second-order time-domain Maxwell equation

Se ond-Order Time-Domain Maxwell Equation

D.Sármány a,

∗

,M.A.Bot hev a

,J.J.W.vanderVegt a

,J.G.Verwer b

a

DepartmentofAppliedMathemati s,UniversityofTwente,P.O.Box217,7500AEEns hede, theNetherlands

b

CenterforMathemati sandComputerS ien e,P.O.Box94079, 1090GBAmsterdam,theNetherlands

Abstra t

This arti le omparesthe dis ontinuous Galerkin nite element method (DG-FEM) with the

H(curl)

- onformingFEM inthedis retisationofthese ond-ordertime-domainMaxwellequationswithpossibly

nonzero ondu tivity term. While DG-FEM suers from an in reased number of degrees of freedom

omparedwith

H(curl)

- onformingFEM,ithasthe advantageof apurelyblo k-diagonalmassmatrix. Thismeansthat,aslongasanexpli ittime-integrations hemeisused,itisnolongerne essarytosolve

alinearsystem at ea h time step a lear advantageover

H

(curl)

- onforming FEM.It is knownthat DG-FEMgenerally favours high-ordermethods whereas

H

(curl)

- onforming FEM is moresuitable for low-order ones. The noveltywe provide in this work is adire t omparison of the performan e of the

twomethodswhenhierar hi

H

(curl)

- onformingbasisfun tionsareuseduptopolynomialorder

p

= 3

. Themotivation behindthis hoi eof basisfun tionsis itsgrowingimportan e inthedevelopmentof

p

-and

hp

-adaptiveFEMs.

Thefa t that weallowfor nonzero ondu tivityrequires spe ial attentionwithregardsto the

time-integrationmethodsappliedtothesemi-dis retesystems. High-orderpolynomialbasiswarrantstheuse

ofhigh-ordertime-integrations hemes,butexistinghigh-orders hemesmaysuerfromatoosevere

time-stepstabilityrestri tionasresultofthe ondu tivityterm. Weinvestigateseveralalternativesfromthe

pointof viewof a ura y,stability and omputationalwork. Finally, we arry outanumeri alFourier

analysis to study the dispersion and dissipation properties of the semi-dis rete DG-FEM s heme and

severalof thetime-integrationmethods. Itis instru tivein ourapproa hthat thedispersionand

dissi-pationpropertiesofthespatialdis retisationandthoseofthetime-integrationmethodsareinvestigated

separately,providingadditionalinsightintothetwodis retisationsteps.

Keywords:

H

(curl)

- onformingniteelementmethod, dis ontinuousGalerkinniteelementmethod, numeri altimeintegration,se ond-orderMaxwellwaveequation

1. Introdu tion

High-order nite element methods (FEM) are an in reasingly important te hnology in large-s ale

ele tromagneti simulationsthanks to theirability to ee tivelymodel omplexgeometri al stru tures

and long-time wave propagation. It has long been known that the standard

H

1

- onforming FEM for

ele tromagneti waves may result in non-physi al, spurious solutions. Instead, one may naturallyopt

forthe

H(curl)

- onformingFEM pioneered by Nédéle [1, 2℄ and Bossavit[3, 4℄. Ithas theadvantage of mimi king the geometri al properties of the Maxwell equations at the dis rete level. However, in

time-domain omputationsitrequiressolvinglinearsystemswithmassmatri esevenifanexpli it

time-integration method is used. One attra tive alternative also free of spurious solutions under ertain

onditionsis the dis ontinuous Galerkin FEM (DG-FEM) [5, 6, 7℄, where the resulting massmatrix

isblo k-diagonalandthereforethe omputational ostofitsinversionisnegligible. Butthisadditional

exibility omesata ost. ThenumberofdegreesoffreedominDGdis retisationsishigherthanthatin

the

H

(curl)

- onformingdis retisation, although thedieren e de reasesasthepolynomialorder in the

∗

Correspondingauthor

Emailaddresses: d.sarmanymath.utwente.n l(D.Sármány),m.a.bot hevmath.utwente .nl (M.A.Bot hev),

forbothmethodswhenamesh of320tetrahedraandthird-orderpolynomialsareused.

0

1000 2000 3000 4000 5000 6000

0

1000

2000

3000

4000

5000

6000

nz = 469124

0

5000

10000

15000

0

2000

4000

6000

8000

10000

12000

14000

16000

18000

nz = 795328

Figure1: Sparsity patternofthemassmatrixfor

H(curl)

- onformingFEM (left)andDG-FEM(right) forameshwith320elements. Third-orderpolynomialsareused,whi hmeansthatthesizeoftheblo ks

intherightplotis60

×

60. Notethedieren einsizebetweenthetwomatri es.

Sothereappearstobeatrade-obetweenthetwomethodsintime-domain omputations. Ingeneral,

the

H(curl)

- onformingapproa his moree ientwithlow-orderpolynomialsandDG-FEMwith high-order ones. The expe ted break-evenpoint depends on anumber of fa tors, su h asthe onditioning

andsparsityofthemassand stinessmatri es intheresultingsemi-dis retesystems. Asanovelty,the

fo usofthisworkistoprovidea omparisonofthe omputationalperforman eofthetwomethodswhen

hierar hi

H(curl)

- onforming basis fun tions [8, 9℄ are used on tetrahedral meshes. The motivation behindthis hoi e isthat these basisfun tionsplayanevermoreimportantrole in thedevelopmentof

p

-and

hp

-adaptivemethods[10℄fortheMaxwellequation.

Throughout the arti le, the dierentdis retisation te hniques are applied to the three-dimensional

Maxwellequationsinthese ond-ordertime-dependentform,

ε

r

∂

2

_E

∂t

2

+ σ

∂E

∂t

+ ∇ × µ

−1

r

∇ × E
= −

∂J

∂t

,

(1)

withhomogeneousboundary onditions

n

× E = 0

. Allquantitiesaredimensionless 1

in(1), where

E

is theele tri eldand

J

is theele tri urrentdensity. Thevalues

σ

,

ε

r

and

µ

r

areassumedto be time-independent onstant s alars,andthey respe tivelydenote ondu tivity, relative diele tri permittivity

andrelativemagneti permeability. Ifthedomainislledwithnon ondu tivematerial,thedampingterm

σ

∂

_∂t

E

is absent. If, inaddition, thesour eterm

−∂

t

J

isalsotakento bezero,wehavethe onservative Maxwellwaveequation.

Followingthemethodoflines, werstdis retisethespatial operators,usingthe

H(curl)

- onforming FEM or theDG-FEM. In either ase, we arriveat asemi-dis rete system in the form of se ond-order

ordinarydierentialequations(ODEs)in

R

n

,

M

ε

u

′′

+ M

σ

u

′

+ S

µ

u = j,

(2)

where

u

istheunknownve torof

N

s alar oe ientsasso iatedwiththeapproximationoftheele tri eld

E

. Thesour e term

j

isthe proje tionof

−∂

t

J

onto thenite-elementspa e andin generalmay also ontainboundarydata. Forsimpli ity,however,werestri tourselvesto thehomogeneousDiri hlet

1

We an derivethe dimensionlessformbyusingthe s alings

x

= ˜

x

/ ˜

L

,

t = ˜

t/( ˜

L/˜

c0)

,

E

= ˜

E

/( ˜

Z0

H0

˜

)

,

H

= ˜

H

/ ˜

H0

,

J

_{= ˜}

J

_{/( ˜}

_{H0/ ˜}

_L)

and

σ = ˜

J ˜

L ˜

Z0/ ˜

E

,withtildedenotingthe dimensionalquantities. Here

L

˜

isthe referen elength,

c0

˜

=

(˜

µ

0

ε

˜

0

)

−

1

/2

isthespeedoflightinva uum,

H

isthemagneti eld(eliminatedin(1)),

H

˜

0

isthe referen emagneti eld strengthand

Z

˜

0

= (i˜

ω ˜

µ

0

/(˜

σ + i˜

ω˜

ε

0

))

1

/2

istheintrinsi impedan e,with

ω

˜

beingtheangularfrequen yand

i

theimaginary unit.

boundary ondition,

n

× E = 0

,inthisarti le. Ea htermintheleft-handsideof(2) orrespondstothe respe tivetermin theleft-handside of(1). Themassmatrix

M

ε

issymmetri positivedeniteandthe ondu tivitymatrix

M

σ

is symmetri positive semi-denite. In addition, for onstants alars

σ

and

ε

r

thematri es

M

ε

and

M

σ

are identi alupto a onstant. The stinessmatrix

S

µ

isthedis retisationof thewavetermandissymmetri positivesemi-denite.

Convergen e resultsfor the

H(curl)

- onformingsemi-dis rete approximation (2) are relatively well-established[11, 12℄. Resultson thesemi-dis rete DGdis retisation aremorere ent: energy-norm

esti-mates [13℄ and

L

2

-estimates [14℄ havebeen derived for theMaxwellequations; optimal errorestimates

for the fully-dis rete se ond-order s alar wave equation have been provided in [15℄; and a promising

energy- onservinglo al-timesteppings hemehasbeendevelopedin [16℄.

A vital feature of (1) and (2) from the point of view of time integration is that it in ludes the

ondu tivity

σ

. Evenmoderatevaluesof

σ

mayresultin aprohibitivelysmalltimestepformanyofthe populartime-integrations hemes. Therefore,wepayspe ialattentiontotime-integrationmethodsthat

treatthe ondu tivitymassmatrix

M

σ

inanimpli itway. Manyofsu hmethodsandothersdis ussedin thisarti lehavebeenpreviouslystudiedin[17℄forthesystemofrst-orderMaxwellequationsdis retised

bythelowest-order

H(curl)

- onformingelements. Seealso[18℄formoredetails on ompositionmethods forthe ondu tion-freeMaxwellequations.

Thesemi-dis retesystem(2) onserves(dis rete)energyforthespatialdis retisationsdis ussedhere,

sin ethese arebothsymmetri . Hen e,usinganenergy- onservativetime-integrationmethod resultsin

a onservativefully-dis retes heme. Weinvestigate thedispersionand dissipationerrorof thes hemes

in two steps. First,wedetermine thedispersionerrorof the semi-dis retes heme bysolvingthe

time-harmoni eigenvalueproblem orrespondingtothesemi-dis retesystem. Se ond,we anthenapplyany

giventime-integrations hemeto asimple,but equivalent,model problemthat in ludestheinformation

ofthesemi-dis retenumeri alfrequen y,andthusdenethedispersion(and,ifthereisany,dissipation)

errorof thetime-integration method. This approa hshowsifthe dispersionerroris dominated bythe

spatialortemporaldis retisationapie e ofinformationthat mayproveusefulin de idingwhether or

nottogoforhigh-ordertime-integrations hemes.

The omputational performan e of the

H(curl)

- onformingmethod hinges to a great degree on ef- iently solvingthe linearsystem with themass matrix. A numberof advan ed te hniques have been

proposedre ently,in ludingmasslumping[19,20,21℄,theexpli it omputationofanapproximatesparse

inversemass matrix[22℄,orthe onstru tionofspe ialpre onditioners. These approa hes,however,do

notin their urrentstatesprovideageneralframework andtherefore annotbeextended to high-order

dis retisationsinastraightforwardmanner. Thatisthereasonwhyinthisarti leweresorttostandard

pre onditioners. It is of ourse also possibleto use sparse dire t solversbut in test problemswefound

thattheyaretoomemorydemanding forlarge-s alethree-dimensional omputations.

The remaining part of the arti le is organised as follows. The weak formulations of the

H(curl)

- onformingFEM andtheDG-FEMaregivenin Se tion2. Thesemi-dis retesystemarisingfromeither

of the spatial dis retisationsis analysed in Se tion 3, while we briey des ribe a number of the most

widely-usedtime-integrationmethodsinSe tion4. Numeri alexamplesthat omparethe omputational

performan eofthetwoniteelementapproa hesarepresentedinSe tion5,wherewetestbothlow-order

andhigh-orderapproximations. Se tion 6 on ludesthearti lewithnalremarks.

2. The weak formulation

Before we present theweak formulationsthat result from the

H(curl)

- onforming and theDG dis- retisations, we introdu ethe tessellation

T

h

that partitions thepolyhedral domain

Ω ⊂ R

3

into aset

of tetrahedra

{K}

. Throughout the arti le we assume that the mesh is shape-regularand that ea h tetrahedronis straight-sided. The notations

F

h

,

F

i

h

and

F

b

h

stand respe tively for the set of all fa es

{F }

,thesetofallinternalfa es, andtheset ofallboundaryfa es. Onthe omputationaldomain

Ω

,wedene thespa es

H(curl; Ω) :=

n

u

_∈

L

2

(Ω)



3

: ∇ × u ∈

L

2

_(Ω)



3

o

,

H

0

(curl; Ω) :=

n

u

_{∈ H(curl; Ω)}







n

× u = 0

on

∂Ω

o

,

andthe

L

2

innerprodu t

(·, ·)

(u, v) =

Z

Ω

u

_{· v dV.}

The ontinuous weak formulation of (1) now reads asfollows: Find

E

∈ H

0

(curl, Ω)

su h that

∀w ∈

H

0

(curl, Ω)

therelation

∂

2

∂t

2

(ε

r

E, w) +

∂

∂t

(σE, w) + µ

−1

r

∇ × E, ∇ × w
= −

 ∂J

∂t

, w



(3)

issatised. Seee.g. [23,12℄.

2.1. Weakformulationof the globally

H(curl)

- onforming dis retisation

Inorder todis retise (3), werstintrodu ethe niteelementspa easso iatedwiththetessellation

T

h

. Let

P

p

(K)

bethe spa eof polynomialsof degreeat most

p ≥ 1

on

K ∈ T

h

. Overea helement

K

the

H(curl)

- onformingpolynomialspa eisdened as

Q

p

=

n

u

_{∈ [P}

p

(K)]

3

; u

T

|

F

K

i

∈

P

p

(F

K

i

)



2

; u · τ

j

|

e

K

j

∈ P

p

(e

K

j

)

o

,

(4) where

F

K

i

,

i = 1, 2, 3, 4

arethefa es of theelement;

e

K

j

,

j = 1, 2, 3, 4, 5, 6

aretheedges oftheelement;

u

T

is the tangential omponent of

u

; and

τ

j

is the dire ted tangential ve tor on edge

e

K

j

. For the onstru tionof

Q

p

,weuseaset of

H(curl)

- onforminghierar hi basisfun tions[8,9℄. Next,weintrodu ethedis retespa eofglobally

H(curl)

- onformingfun tions

Υ

p

h

:=

n

υ

_{∈ [H}

0

(curl, Ω)]

3







υ

|

K

∈ Q

p

, ∀K ∈ T

h

o

,

and let the set of basis fun tions

{ψ

i

}

span the spa e

Υ

p

h

. See [12℄ for a detailed dis ussion on both ontinuousanddis rete

H(curl)

- onformingspa es. We anthenapproximatetheele tri eld

E

as

E

_{≈ E}

h

=

X

i

u

i

(t)ψ

i

(x),

(5)

fromwhi hthedis reteweakformulationreadsasfollows: Find

E

h

∈ Υ

p

h

su hthat

∀φ ∈ Υ

p

h

therelation

∂

2

∂t

2

(ε

r

E

h

, φ) +

∂

∂t

(σE

h

, φ) + µ

−1

r

∇ × E

h

, ∇ × φ
= −

 ∂J

∂t

, φ



(6)

issatised. Notethat(6)issatisedifandonlyifitissatisedforeverybasisfun tion

ψ

i

, i = 1, . . . , N

, with

N

beingtheglobalnumberofdegreesoffreedom. Asaresult,substitutionof(5)into(6)yieldsthe semi-dis retesystem(2)with

[M

ε

]

ij

= ε

r

ψ

i

, ψ

j

,

[S

µ

]

_ij

= µ

−1

r

∇ × ψ

i

, ∇ × ψ

j

,

[M

σ

]

ij

= σψ

i

, ψ

j

,

[j]

i

= −

 ∂J

∂t

, ψ

i



.

Ea hoftheabovematri es

M

ε

,

M

σ

and

S

µ

hasalargenumberentriesfarothediagonal,in reasing the omputational ostforbothexpli itandimpli ittime-integrationmethods.

2.2. Weakformulationof DG-FEM

In ontrast to the

H(curl)

- onforming dis retisation, in DG-FEM we are looking for the dis rete solutionin thespa e

Σ

p

_h

:=

n

σ

_∈

L

2

(Ω)



3







σ

|

K

∈ Q

p

, ∀K ∈ T

h

o

.

Thatis,weallowthepolynomialfun tionstobefullydis ontinuousa rosselementinterfa esandassume

that the set of basis fun tions

{ψ

i

}

now span the spa e

Σ

p

h

. Instead of enfor ing ontinuity of the tangential omponents, the information between elements is now oupled through the numeri al ux

needto introdu emorenotation.

Consideran interfa e

F ∈ F

h

betweenelement

K

L

andelement

K

R

, and let

n

L

and

n

R

represent

theirrespe tiveoutwardpointingnormalve tors. Wedene thetangentialjumpandtheaverageofthe

quantity

u

a rossinterfa e

F

as

[[u]]

_T

= n

L

× u

L

+ n

R

× u

R

and

{{u}} = u

L

_{+ u}

R

_/2,

respe tively. Here

u

L

and

u

R

are thevalues of thetra eof

u

at

∂K

L

and

∂K

R

, respe tively. At the

boundary

Γ

,weset

{{u}} = u

and

[[u]]

T

= n × u

. Wefurthermoreintrodu ethegloballifting operator

R(u) :

L

2

_(F

h

)



3

→ Σ

p

h

as

(R(u), v)

Ω

=

Z

F

h

u

_{· {{v}} dA,}

_{∀v ∈ Σ}

p

_h

,

(7)

and,foragivenfa e

F ∈ F

h

,thelo alliftingoperator

R

F

(u) :

L

2

_{(F )}



3

→ Σ

p

h

as

(R

F

(u), v)

Ω

=

Z

F

u

_{· {{v}} dA,}

_{∀v ∈ Σ}

p

_h

.

(8)

Note that

R

F

(u)

vanishes outside the elements onne ted to the fa e

F

so that for a given element

K ∈ T

h

wehavetherelation

R(u) =

X

F ∈F

h

R

F

(u),

∀u ∈

L

2

(F

h

)



3

.

(9)

The dis rete weak formulation for DG-FEM now reads as follows [25℄: Find

E

h

∈ Σ

p

h

su h that

∀φ ∈ Σ

p

h

therelation

∂

2

∂t

2

(ε

r

E

h

, φ) +

∂

∂t

(σE

h

, φ) + µ

−1

r

∇

h

× E

h

, ∇

h

× φ

−

Z

F

h

[[E

h

]]

_T

· {{∇

h

× φ}} dA −

Z

F

h

{{∇

h

× E

h

}} · [[φ]]

_T

dA

+

X

F ∈F

h

C

F

(R

F

([[E]]

T

), R

F

([[φ]]

T

))

Ω

= −

 ∂J

∂t

, φ



(10)

issatised,wheretheoperator

∇

h

denotestheelementwiseappli ationof

∇

. Forstability,the onstant

C

F

hastosatisfythe ondition[25℄

C

F

≥ n

f

C

1

+ min

 1

2

,

1

C

2



,

where

C

1

and

C

2

are positive onstantsand

n

f

denotesthenumberofsides ofea h element,that isfor tetrahedra

n

f

= 4

. Asa onsequen e,the onstant

C

F

isindependentofboththepolynomialorderand themesh size.

Again,(10) issatisedifand only ifitis satisedforeverybasis fun tion

ψ

i

, i = 1, . . . , N

, with

N

beingtheglobalnumberofdegreesoffreedom. Substitutionof

E

≈ E

h

=

P

i

u

i

(t)ψ

i

(x)

into(10)yields thesemi-dis retesystem(2)with

[M

ε

]

ij

= ε

r

ψ

i

, ψ

j

,

[M

σ

]

ij

= σψ

i

, ψ

j

,

[j]

i

= −

 ∂J

∂t

, ψ

i



,

[S

µ

]

ij

= µ

−1

r

∇

h

× ψ

i

, ∇

h

× ψ

j

−

Z

F

h

[[ψ

i

]]

T

·



∇

h

× ψ

j

dA

−

Z

F

h

{{∇

h

× ψ

i

}} ·

ψ

j



T

dA +

X

F ∈F

h

C

F



R

F

([[ψ

i

]]

T

), R

F

(ψ

j



T

)



Ω

.

Thematri es

M

ε

and

M

σ

arenowblo k-diagonalwiththeelementwisematri esbeingtheblo ks. How-ever,thestinessmatrix

S

µ

hasstillmanyentriesfarothediagonalbe auseofthefa eintegralsinits onstru tion. That iswhy, DGin generalwarrantstheuseof expli ittime-integrations hemes butnot

impli itones.

We emphasise that (10) is only one of many possible formulations of DG-FEM, depending on the

numeri aluxone hoosestouse. Theonewehaveintrodu edhereisbasedonthenumeri alux from

[26℄(seealso[27℄),andwasanalysedin detailforthetime-harmoni Maxwellequationsin[25℄. Seealso

[24℄foranoverviewofDG-FEMmethodsforellipti problemsandforalargenumberofpossible hoi es

forthenumeri alux.

2.3. Theenergynorm

Convergen e results for FEMs are generally derived not only in the

L

2

-norm but also in a norm

asso iated with the dis rete energy of the approximation [13, 15℄. These are dened for the

H(curl)

- onformingandDGdis retisationsas

kvk

2

H(curl)

= kvk

2

+ k∇ × vk

2

and

kvk

2

DG

= kvk

2

+ k∇

h

× vk

2

+ k

h

−

1

2

[[v]]

T

k

2

F

h

,

respe tively. Intheabovedenition,

k · k

F

h

denotesthe

L

2

_(F)

norm andh

(x) = h

F

isthediameter of fa e

F

ontaining

x

. We note that thetwodenitions ofthe energynorm are a tuallyidenti al as

∇

h

be omes

∇

and

[[·]]

T

vanishesif

H(curl)

- onformingdis retisationisused.

3. Stability ofthe semi-dis rete system

To arryoutabasi stabilityanalysis,wersttransform(2)intoarst-ordersystemofODEs,

u

′

_{= v,}

(11)

M

ε

v

′

+ M

σ

v + S

µ

u = j.

Re allthat

S

µ

is symmetri andtherefore using theinner-produ tnotation fordis rete ve torswe havetheproperty

d

dt

v

T

_M

ε

v + u

T

S

µ

u
=

dv

T

dt

M

ε

v + v

T

_M

ε

dv

dt

+

du

T

dt

S

µ

u + u

T

_S

µ

du

dt

=

2v

T

(−M

σ

v − S

µ

u + j) + 2v

T

S

µ

u = 2v

T

j − 2v

T

M

σ

v.

(12)

If

j = 0

,thisentailsstability,that is

d

dt

v

T

_M

ε

v + u

T

S

µ

u

= −2v

T

M

σ

v ≤ 0,

sin e,for onstant

σ

,thematrix

M

σ

ispositivedeniteif

σ > 0

and

M

σ

= 0

if

σ = 0

. Therefore,if

σ = 0

inadditionto

j = 0

,(12)shows onservation.

Inordertouseastabilitytestmodelintrodu edlaterinthisse tion,wetransform(11)toanequivalent

expli it form. To do so, we multiply the rst equation in (11) with

M

ε

and introdu e the Cholesky fa torisation

LL

T

_{= M}

ε

. Thenewvariables

v = L

˜

T

_v

and

u = L

˜

T

_u

thensatisfythesystem

 ˜

u

′

˜

v

′



=



0

I

− ˜

S

µ

− ˜

M

σ

  ˜

u

˜

v



+

0

_˜j



,

(13) where

˜j = L

−1

_j,

_S

˜

_µ

_{= L}

−1

_S

µ

L

−T

,

M

˜

σ

= L

−1

M

σ

L

−T

.

Sin eboththe ondu tivity oe ient

σ

andthe permittivity oe ient

ε

r

are onstants alarsin (1), thematrix

M

˜

σ

in(13)isthe onstantdiagonalmatrix

˜

M

σ

= γI,

γ =

σ

ε

r

.

From this we an derivea two-by-twosystem through whi h stability of time-integration methods for

(11) anbeexamined.

Thematrix

S

˜

µ

issymmetri positivesemi-denite soit anbede omposedas

S

˜

µ

= UΛU

T

,where

Λ

isadiagonalmatrixwiththeeigenvaluesof

S

˜

µ

onitsdiagonal

λ

1

≥ λ

2

≥ . . . ≥ λ

r

≥ λ

r+1

= λ

r+1

= . . . = λ

n

= 0,

where

r

is therankof thematrix. Thematrix

U

is orthogonaland its olumns arethe eigenve torsof

˜

S

µ

. Usingapermutation matrix

P

,wehave

A =



0

I

− ˜

S

µ

− ˜

M

σ



=



0

U U

T

−UΛU

T

_−γI



=

U

0

0

U

 

0

I

−Λ −γI

 U

T

₀

0

U

T



=

U

0

0

U



PΛ

P

P

T

U

T

₀

0

U

T



,

(14)

where

Λ

P

isablo k-diagonalmatrixwithtwo-by-twoblo ks



0

1

−λ

k

−γ



,

k = 1, . . . , N.

(15)

Thisallowsusto statethefollowingproposition.

Proposition1. Assumethat

σ

and

ε

r

ares alar and

γ = σ/ε

r

. Then the matrix

A

has (i)

n − r

zeroeigenvalues,

(ii)

n − r

eigenvalueswhi h equal

−γ

,

(iii)

2r

eigenvalues whi hare

−γ ±

pγ

2

_{− 4λ}

k

2

,

k = 1, . . . , r.

Thus, theorthogonaltransformation

V ≡ (

U 0

0 U

) P

de ouples (13)into

r

two-by-twosystems

 ˆ

u

′

ˆ

v

′



=

 0

1

−λ −γ

  ˆ

u

ˆ

v



+

0

_ˆj



,

with

λ = λ

k

> 0

,

k = 1, . . . , r

,and

n − r

two-by-twosystems

 ˆ

u

′

ˆ

v

′



=

0

1

0

−γ

  ˆ

u

ˆ

v



+

0

_ˆj



.

Forthestabilityanalysis,wemaynegle tthesour etermandthusarriveatthetwo-by-twostabilitytest

model

 ˆ

u

′

ˆ

v

′



=

 0

1

−λ −γ

  ˆ

u

ˆ

v



,

λ ≥ 0,

γ ≥ 0.

(16)

Theattra tivefeature of this formulationis that stability for thetest model (16) indu es stability for

(11)in thenormgeneratedbytheinnerprodu tin (12).

Useful though equation (16) is, it is important to emphasise that the derivation of (16) requires

Probably the most popular time-integration methods to use in ombination with high-order DG

methods are high-order Runge-Kutta methods, giving rise to what are olle tively alled the

Runge-KuttaDG (RKDG) methods [5℄. For ontinuous and

H(curl)

- onforming FEMs geometri integrators arealsowidely usedthanks totheir ability to onservesymple ti ity

2

at thedis retelevel[18℄. In this

se tion,webrieyre allthe onstru tionofthese twofamiliesofmethodsandwealsodis usslo al and

globalRi hardsonextrapolation.

Thehighest-order polynomial weuse within thenite elementmethods is

p = 3

. Forboththe DG andthe

H(curl)

- onformingmethods,this orrespondstofourth-order onvergen eforthesemi-dis rete system(2) provided that thesolution issmooth[13, 14, 11, 12℄. Therefore, we now onlydis uss

time-integrationmethodsthat arealsoat mostfourth-ordera urate. Extensionto higherorder,however,is

usuallystraightforward.

Forinvestigatingthepropertiesofanygiventime-integrationmethod,let

τ

denotethetime-stepsize andintrodu e

z

λ

= τ

√

λ

and

z

γ

= τ γ

. 3

Thestabilityofthetime-integrationmethod anthen,ingeneral,

bebestinspe tedthroughthe(numeri allydetermined)stabilityregion

S = {(z

λ

, z

γ

) : z

λ

, z

γ

≥ 0

with

|µ| < 1, µ

eigenvaluesoftheampli ationoperator

}

asso iatedwiththetest model(16).

4.1. Runge-Kuttamethods

Out of the many dierent typesof Runge-Kuttamethods, strong-stability-preserving Runge-Kutta

methods (SSPRK) [5℄ are parti ularly well suited for the time integration of semi-dis rete hyperboli

problems.

With thedenition oftheinitial values

U

0

= u

n

and

V

0

= v

n

forthetime stepfrom

t

n

to

t

n+1

,the general

s

-stageSSPRK s hemefor(11)reads

U

k

=

k−1

X

l=0

(α

kl

U

l

+ τ β

kl

V

l

) ,

M

ε

V

k

=

k−1

X

l=0

(α

kl

V

l

+ τ β

kl

(−S

µ

U

l

− M

σ

V

l

+ j(t

l

))) ,

u

n+1

= U

s

,

v

n+1

= V

s

,

(17)

where

k = 1, . . . , s

while

α

kl

and

β

kl

are the oe ientsin the SSPRK method. Applying (17)to the testequation(16),theampli ationoperator

M

s

ssp

ofan

s

-stageSSPRKmethodis

M

k

ssp

=

k−1

X

l=1

B

kl

M

l−1

ssp

with

B

kl

= α

kl

1 0

0

1



+ β

kl

 0

1

z

2

λ

−z

γ



,

(18)

where, again,

k = 1, . . . , s

and

M

0

ssp

is the identity. We show the stability regions of several SSPRK s hemes in Figure 2, where we refer to an

s

-stage

p

th-order SSPRK method as SSPRK(

s

,

p

). It is importantto emphasisethat any(standard aswell as`nonstandard'su hasSSP) expli it

s

-stage

p

th-orderRKmethodswith

s = p

hasthesameampli ationmatrix. Inthose asesthe hoi eof oe ients onlydeterminestheSSPpropertyandnottheshapeofthestabilityregion. Forthe aseswhen

s = p + 1

, we display the stability regions of the methods that were derived and analysed in [28, 29, 30℄. The

plotssuggestthat in reasingthenumberof stages,whilekeepingthepolynomialorderxed, resultsin

2

Thepreservationofsymple ti ityisimportantbe auseitisrelatedtoenergy.Morepre isely,forsymple ti integrators

theerror intotalenergywillremainwithina ertainmarginthroughouttheentiretimeintegration.

3

Thesevaluesappearinanaturalwayintheampli ationmatri esofmosttime-integrationmethodsdes ribedlaterin

of introdu ing an additional stage at ea h time step. This is in line with known results forthe linear

adve tionequation[28,29℄. Nevertheless,expli itSSPRKmethodstreatthe ondu tionterm,aswellas

thewaveterm,expli itly,whi hentailsatime-step onditionthat istoorestri tiveevenformoderately

ondu tivematerials(see nextse tion). Notethatthese ond-ordermethodsareonlystable for

z

γ

> 0

, i.e. for

σ > 0

.

4.2. Compositionmethods

Composition methods [31, 32, 33℄ are espe ially suitablefor geometri integration [18℄ andthus for

thetimeintegrationofrst-orderHamiltoniansystems. Ourdes riptionofthe ompositionmethodshere

stri tlyfollowsthatin [17℄andwereferto thatworkformoredetails.

These ond-order omposition methodfor(11)isdened as

u

n+1/2

− u

n

τ

=

1

2

v

n

,

M

ε

v

n+1

− v

n

τ

= −S

µ

u

n+1/2

−

1

2

M

σ

(v

n

+ v

n+1

) +

1

2

(j(t

n

) + j(t

n+1

)) ,

u

n+1

− u

n+1/2

τ

=

1

2

v

n+1

,

(19)

whi h isakin to theubiquitousleapfrog s heme, withthe onlydieren ebeingin the treatmentof the

sour eterm( f. [34℄). Ifappliedtothetestmodel (16),ithastheampli ationmatrix

M

co2

=

1 −

1

2

z

2

λ

+

1

2

z

γ

1 −

1

4

z

2

λ

−z

2

λ

1 −

1

2

z

2

λ

−

1

2

z

γ



,

(20)

whi h entailsthestabilityproperties:

z

λ

≤ 2

if

z

γ

= 0

and

z

λ

< 2

if

z

γ

> 0

. An attra tivefeature of thismethodoverexpli itRKmethodsisthat itisun onditionallystablewithrespe ttothe ondu tion

term.

Inprin iple,itispossibleto onstru tanarbitraryhigh-order ompositionmethod[31℄. Inthisarti le,

however,weareonlyinterestedinatmostfourth-ordera uratemethodssowewillnowonlydis ussthe

fourth-order omposition method. Wedene the initial valuesfor theinner time step as

U

0

= u

n

and

V

0

= v

n

,timelevels

t

u

,

t

v

for

u

,

v

and oe ients

β

0

= α

0

= 0,

β

1

= α

5

=

14−

√

19

108

,

β

2

= α

4

=

−23−20

√

19

270

,

β

3

= α

3

=

1

₅

, β

4

= α

2

=

−2+10

√

19

135

,

β

5

= α

1

=

146+5

√

19

540

.

Thefourth-order ompositionmethod [31,32,17℄for(11)nowreads

U

k

− U

k−1

τ

= (β

k

+ α

k−1

) V

k−1

,

M

ε

V

k

− V

k−1

τ

= β

k

−S

µ

U

k

− M

σ

V

k−1

+ j(t

v

k−1

)
+ α

k

(−S

µ

U

k

− M

σ

V

k

+ j(t

v

k

)) ,

v

n+1

= V

s

,

u

n+1

= U

s

+ α

s

τ V

s

,

(21)

where

k = 1, . . . , s

,

s = 5

is the number of internal time levels, and

t

v

k

= t

n

+ (˜

α

k

+ ˜

β

k

)τ

and

t

u

k

=

t

n

+ (˜

α

k−1

+ ˜

β

k

)τ

with the oe ients

α

˜

k

= α

1

+ · · · + α

k

and

β

˜

k

= β

1

+ · · · + β

k

. Theampli ationoperatorof(21)whenapplied to(16)isthen

1

Y

k=5

1

1 + α

k

z

γ



1 + α

k

z

γ

(1 + α

k

z

γ

) (α

k−1

+ β

k

)

− (α

k

+ β

k

) z

2

λ

1 − β

k

z

γ

− (β

k

+ α

k−1

) (β

k

+ α

k

) z

λ

2



.

(22)

Animportantpropertyofanyfourth-order omposition methodisthatitinevitably ontainsanegative

oe ient,whi hin our ase is

α

4

= β

2

. This entailsastabilityrestri tionthat is onditionalevenfor an impli itly treated ondu tion term. This is illustrated in Figure 3, where parts of the upper right

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

(a)SSPRK(

2

,

2

)

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

(b)SSPRK(

3

,

2

)

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

( )SSPRK(

3

,

3

)

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

(d)SSPRK(

4

,

3

)

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

(e)SSPRK(

4

,

4

)

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

(f)SSPRK(

5

,

4

)

Figure2:Stabilityregions(shadedareas)forseveralexpli itSSPRK(

s

,

p

)methods,where

s

isthenumber ofstagesand

p

istheorderofthemethod. Notethatallexpli itRKmethodswith

s = p

havethesame stabilityregionsasSSPRK(

s

,

p

)with

s = p

(left olumn).

halfof thestability regionfor(21)isshown. Stability isguaranteedaslongas

z

γ

< 2.4

and

z

λ

< 3

, or equivalently,if

τ < 2.4/γ

and

τ < 3/

√

λ

.

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

1

2

3

4

0

0.5

1

1.5

2

2.5

3

3.5

4

Figure3: Stabilityregion(shaded area)of thefourth-order ompositionmethod. Therightplot zooms

inontheregionwhere stabilityisguaranteed. (Cf. Figure5.1 in[17℄.)

4.3. Fourth-orderglobal Ri hardson extrapolation

As already mentioned in the previous se tion, when

σ > 0

the stability ondition may be very restri tive even for moderately ondu tive materials. In these ases, high-order omposition methods

andSSPRKmethodsarenot ompetitive. Instead,onewouldprefertouseexpli itmethodswhi htreat

the ondu tionterminanun onditionallystablemanner. Sin ethese ond-order ompositionmethodis

su hamethod,extendingittohigherorderthroughRi hardsonextrapolationisanobviousalternative.

Wereferto[17℄ foradetailed dis ussiononthestabilitypropertiesof thefourth-orderlo aland global

versionsoftheRi hardsonextrapolation. Herewerstre allthe onstru tionofthefourth-orderglobal

Ri hardsonextrapolation(GEX4)

u

gex4

τ

=

4

3

u

co2

τ

2

−

1

3

u

co2

τ

,

(23) where

u

co2

τ

2

and

u

co2

τ

denotetheresultsat naltime omputed bythese ond-order omposition method withtimesteps

τ

2

and

τ

, respe tively. Sin e extrapolationonlytakespla e on eatthenal timeofthe integration,thismethodhasthesamestabilitypropertiesasthese ond-order ompositionmethod. Note

thatitonlyneedsthree timesasmu h omputationalworkpertimestep.

Forlongtimeintegrationandintheabsen eofdamping,globalextrapolationmaynotbesu iently

ee tiveinannihilatingleadingerrorterms. Inthese ases,thelo alversionofRi hardsonextrapolation

when theextrapolationis performedatea htime stepisusually morebene ial. Thelo alversion

of(23)is, however,notun onditionallystable with respe tto

z

γ

. Instead,we anusethefourth-order lo alextrapolation(LEX4)denedas

u

lex4

τ

=

9

8

u

co2

τ /3

−

1

8

u

co2

τ

,

(24)

where thework per time stepis approximatelyfour times asmu h asthat ofCO2. Theampli ation

operatorofLEX4forthetestmodel(16)reads

9

8

M

3

co2

(z

λ

/3, z

γ

/3) −

1

8

M

co2

(z

λ

, z

γ

),

(25)

where

M

co2

(z

λ

, z

γ

)

denotes the ampli ation operator (20) of CO2. Figure 4 shows the asso iated stability region

S

, whi h indi ates an approximate stability interval

0 ≤ z

λ

≤ 2.85

and un onditional stabilityfor

z

γ

.

5. Numeri alexperiments

In this se tion, we perform numeri al tests to establish the onvergen e rates of the fully dis rete

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

2

4

6

8

10

0

2

4

6

8

10

z

_γ

= τ γ

z

λ

=

τ

λ

1/2

0

1

2

3

4

0

0.5

1

1.5

2

2.5

3

3.5

4

Figure4:Stabilityregion(shadedarea)ofthefourth-orderlo alRi hardsonextrapolation(24). Theright

plotzoomsinontheregionwherestabilityforwaveterm

z

λ

isguaranteed. Stabilityforthe ondu tion term

z

γ

isun onditional.

also arry out a numeri al dispersion analysis of the semi- and fully dis rete system with DG spatial

dis retisation. This isdonein thefollowingway: i)solvethetime-harmoni eigenvalueproblem, whi h

orresponds to the semi-dis rete systemwith Fouriermode initial onditions; ii) apply a hosen

time-integrationmethod to thetest model (16)with the omputedsemi-dis rete numeri al frequen y. This

approa hhastwomainadvantagesoversimplysolvingtheeigenvalueproblemthatresultsfromapplying

the ampli ation matrix dire tly to (11). First, it is more e ient be ause we solve an eigenvalue

problemthat issmallerandalwayssymmetri . Se ond,itmakesitpossibletostudythedispersion(and

dissipation)propertiesofthetime-integrations hemeseparatelyfromthoseofthesemi-dis retes heme.

5.1. Convergen eand omparison of performan e

Weuseasimpletestexampletoillustratethenumeri alperforman eofthetwospatialdis retisation

te hniquesdes ribedinSe tion2. Forbothmethods,thepredi ted onvergen erateofthesemi-dis rete

systemis

O(h

p+1

₎

inthe

L

2

_(Ω)

normand

O(h

p

₎

intheenergynormforsmoothsolutions;seeforexample

[12,35, 25, 13, 14℄. It isthus naturalto hoosethetime-integration methodsu h thatit guaranteesat

leastthe sameorderof onvergen e. Therefore, ifthepolynomialorder in the FEM isat most one we

usethe se ond-order omposition method; if the polynomial order is two orthree weapply oneof the

possiblefourth-ordermethodsdes ribedinSe tion4.

Thenumeri al testsare implementedin

hp

GEM 4

[36℄, ageneralnite elementpa kagesuitable for

solvingavarietyof physi al problems in uid dynami sand ele tromagnetism. To integratethe

semi-dis rete system in time we use PETS [37℄, whi h is parti ularly e ient in omputing matrix-ve tor

multipli ationsandsolvinglinearsystemsforlargesparsematri es. Asastopping riterionin thelinear

solverforthe

H(curl)

- onformingmethod, wesetthetoleran eat

tol

= 10

−8

.

Intheexample,we onsider(1)inthe ubi domain

Ω = (0, 1)

3

. Wedenethetime-independenteld

¯

E(x, y, z) =





sin(πy) sin(πz)

sin(πz) sin(πx)

sin(πx) sin(πy)





and hoosethesour etermto be

−

∂J

_∂t

= ε

r

η

′′

(t) + ση

′

(t) + 2π

2

η(t)

E(x, y, z).

¯

4

E

(t, x, y, z) = η(t) ¯

E

(x, y, z),

η(t) =

3

X

k=1

cos ω

k

t,

(26)

with

ω

1

= 1

,

ω

2

= 1/2

,

ω

3

= 1/3

,

ε

r

= 1

. Weeitherset

σ = 0

or

σ = 60π

, andweintegrateuntil nal time

T

end

= 12π

(exa tly onetimeperiod).

Whentheglobally

H(curl)

- onformingdis retisation[12,11℄isused,weneedtosolvealinearsystem at ea h time step. Sin e the matrix

M

ε

is positive denite, a natural hoi e of linear solver is the pre onditioned onjugate gradient (PCG) method. For simpli ity, we apply the in omplete Cholesky

pre onditionerforallmeshesandpolynomialorders. Weemphasisethat pre onditioningis notanissue

fortheDGmethod sin ewe ansimplyinverttheblo k-diagonalmassmatrixat negligible ost.

As arstexample, we run(26)with

σ = 0

and nal time

T

end

= 12π

, onasequen e of stru tured mesheswith

N

el

= 5, 40, 320, 2560, 20480, 163840

elements. Inea hmesh the largestfa e diameter

h

is exa tlyhalfthatofthepreviousmesh. Weplotthe onvergen eratesinFigure5inboththe

L

2

_(Ω)

-norm

andtheenergynormforpolynomialorders

p = 1, 2, 3

. Notethatthe onvergen eisshownasafun tion of degrees of freedom, whi h is equivalent to showing the onvergen e as fun tion of

1/h

in the DG ase. However,in the

H(curl)

- onforming asethereis somedieren ebetweenthetwo,asthe number of degrees of freedom generally in rease slightly more than eightfold when

h

is halved. Nevertheless, we ansee that the expe ted onvergen e rates are a hieved asymptoti ally for both the DG and the

H(curl)

- onformingmethods. We analsoobservethatittakesfewerdegreesoffreedomforthe

H(curl)

- onformingdis retisation to rea h a given a ura y. Furthermore, we an onrm the well-established

observationthattheuseofhigh-orderapproximationspayso(atleastforsmoothsolutions)intermsof

a ura yperdegreesoffreedom.

Togainfurtherinsightintothe omputational ostsofthetimeintegration,weshowtheperforman e

of the DG method in Tables I and III; and that of the

H(curl)

- onforming method in Tables II and IV. In this parti ular example, we use astru tured mesh with 320elementsand an unstru tured one

with432elements. 5

Althoughthe a ura yofthetwomethodsis omparable,the omputational osts

arenotandthepattern hangesdramati allyastheorderin reases. Thetotalnumberofmatrix-ve tor

multipli ations(matve s)neededtointegrateuntil

T

end

isalwayshigherforthe

H(curl)

- onforming ase thanforthe DGmethod. Thisis notsurprising giventhat at ea h timestep alinear systemhasto be

solved. However,this seeminglyunfavourablepropertydoesnotmanifestsitselfinlonger omputational

time for

p = 1

and

p = 2

on stru tured meshes, thanks in part to the smaller size of the system and in partto aweakertime-steprestri tionin the

H(curl)

- onformingFEM.The situation isdierentfor

p = 3

. Here,thein reasednumberofmatve stranslatesreadilyintomoreCPUtimeonbothstru tured and unstru tured meshes. This is partly be ause of a trade-o between the onditioning of the mass

matrixandtheuseofthehierar hi basis. Massmatri esbasedonhierar hi basestendtoberelatively

badly onditioned. This does not inuen e the performan e of the DG method. But it renders the

H(curl)

- onformingmethod less ee tivebe ause thenumberof iterationsin solving thelinearsystem at ea h time step growssigni antlywith the polynomial order. This ee t is even morepronoun ed

onunstru tured meshes, where DG performs slightlybetterfor

p = 2

already and where the

H(curl)

- onforming omputationfor

p = 3

is ex essivelylongwhi h isone reasonwhyweonly ompleted one ofthem.

The hoi eofthetime-integrationmethoddoesnotinuen e the omputationalresultsmu hin this

example. Nevertheless,LEX4appearstobethemoste ientthankstothebalan ebetweentheallowable

time-stepsize and the omputationalworkneeded per timestep. Wealso note that forthis parti ular

meshthe useof fourth-ordertime-integration methodsmaynotbene essaryevenfor

p = 2, 3

. This is solelybe ausethespatial errorisnotyet intheasymptoti regime andthereforedominates. Wetakea

loserlook atthisshortlyin termsofnumeri aldispersionanddissipation.

Onstru turedmeshes, werepeatexample(26)with ondu tivity

σ = 60π

,whi h orrespondstothe dimensionalvalue

σ = 0.5 S m

˜

−1

,typi alofthehumanabdomen. The onvergen eresultsare shownin

5

Amesh of 320 or432 tetrahedra issu ient to ompare the dierentmethods fromthe pointof view of a ura y

and omputationalwork. Anermeshwouldnaturallygiveamorea uratesolutionbuttherelativeperforman eofthe

TableI:Computational ostsoftheDGmethod forexample(26)with

σ = 0

. Astru turedmeshof 320 elementsisusedwithspatialpolynomialorders

p = 1, 2, 3

.

method #DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 3840 1.2174e-01 4526 0.0167 8s

p = 2

CO2 9600 1.1696e-02 7542 0.0100 114s

p = 2

GEX4 9600 1.2303e-02 22624 0.0100 342s

p = 3

CO4 19200 7.0432e-04 35192 0.0107 2013s

p = 3

GEX4 19200 9.0148e-04 31672 0.0071 1863s

p = 3

LEX4 19200 6.1762e-04 28154 0.0107 1623s

Table II: Computational osts of the

H(curl)

- onforming method for example (26) with

σ = 0

. A stru turedmeshof320elementsisusedwithspatialpolynomialorders

p = 1, 2, 3

.

method #DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 504 2.7283e-01 7783 0.0417 2s

p = 2

CO2 2388 1.3642e-02 59201 0.0250 87s

p = 2

GEX4 2388 1.2942e-02 180067 0.0250 264s

p = 3

CO4 6640 7.4117e-04 817880 0.0268 19683s

p = 3

GEX4 6640 7.7523e-04 736276 0.0179 17492s

p = 3

LEX4 6640 8.5234e-04 653611 0.0268 15510s

TableIII:Computational ostsoftheDG methodfor example(26)with

σ = 0

. An unstru tured mesh of432elementsisusedwithspatialpolynomialorders

p = 1, 2, 3

.

method # DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 5184 1.9583e-01 11878 0.00635 41s

p = 2

CO2 12960 1.3396e-02 19796 0.00381 429s

p = 2

GEX4 12960 1.4316e-02 59384 0.00381 1263s

p = 3

CO4 25920 1.4311e-03 92372 0.00408 7585s

p = 3

GEX4 25920 1.5558e-03 83134 0.00272 6749s

p = 3

LEX4 25920 1.4038e-03 73898 0.00408 5909s

Table IV: Computational osts of the

H(curl)

- onforming method for example (26) with

σ = 0

. An unstru turedmeshof432elementsisusedwith spatialpolynomialorders

p = 1, 2, 3

.

method # DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 744 1.9113e-01 33691 0.02380 18s

p = 2

CO2 3420 1.7294e-02 169129 0.01429 963s

p = 2

GEX4 3420 1.8860e-02 406784 0.01429 2778s

p = 3

CO4 9360  >1e+07 0.01530 >5e+05s

p = 3

GEX4 9360 1.9676e-03 21337490 0.01020 782647s

p = 3

LEX4 9360  >1e+07 0.01530 >5e+05s

Figure6,fromwhi hitappearsthattheyaresimilartothenon ondu tive aseex eptthatoptimalrates

of onvergen earerea hedsooner. On unstru turedmeshes,theexampleisrepeated with ondu tivity

σ = 450π

, a value more typi al of seawater. See Table V for the ondu tivity of asmall sele tion of materials(sour e: en.wikipedia.org/wiki/Ele tri al_ ondu tivity).

The omputationalwork,depi ted inTablesVIIX,also showsasimilarpattern tothe

ondu tion-free ase,ex eptwhen thefourth-order omposition methodis used. Inthat ase,the ondu tion term

poses a stri ter time-step size than the wave term and in reases the number of time steps and thus

10

2

10

4

10

6

10

−2

10

−1

10

0

p = 1, convergence in the L

2

norm

Degrees of Freedom

L

2

error

DG

H(curl)

order 2

10

2

10

4

10

6

10

0

10

1

p = 1, convergence in the energy norm

Degrees of Freedom

H(curl,

Ω

) / DG error

DG

H(curl)

order 1

10

2

10

4

10

−3

10

−2

10

−1

10

0

p = 2, convergence in the L

2

norm

Degrees of Freedom

L

2

error

DG

H(curl)

order 3

10

2

10

4

10

−1

10

0

10

1

p = 2, convergence in the energy norm

Degrees of Freedom

H(curl,

Ω

) / DG error

DG

H(curl)

order 2

10

2

10

4

10

−2

10

−1

10

0

10

1

p = 3, convergence in the energy norm

Degrees of Freedom

H(curl,

Ω

) / DG error

DG

H(curl)

order 4

10

2

10

4

10

−4

10

−3

10

−2

10

−1

10

0

p = 3, convergence in the L

2

norm

Degrees of Freedom

L

2

error

DG

H(curl)

order 3

Figure 5: Convergen e plotsin the

L

2

-norm (left olumn) and in the energy norm (right olumn) for

testexample (26)with

σ = 0

. Inea h plotthe onvergen e ratesof theDG method and the

H(curl)

- onformingmethod areshownalongwiththeexpe tedorderof onvergen e.

H(curl)

- onformingdis retisationbe ausethestinessmatrixintheDGmethodhasasigni antlylarger spe tralradius(andthereforeitstilldeterminesthestability ondition). Ontheunstru turedmeshwith

(

S m

−1

). For the dimensionless value a multipli ation by

120π

is needed. Sour e: en.wikipedia.org/wiki/Ele tri al_ ondu tivity.

Material Condu tivity(

S m

−1

) Note

Silver 63.0e+06 Bestele tri al ondu tor

Copper 59.6e+06

Gold 45.2e+06 Commonlyusedinele tri al onta ts

Aluminium 37.8e+06

Seawater 4.8 Foraveragesalinityof35g/kg

HumanBody 0.0061.5 Variesfromboneto erebrospinaluids

Drinkingwater 0.00050.05

Deionisedwater 5.5e-06 Lowestvalue,withmonoatomi gasespresent

Air 5e-15 Variesslightlydependingonhumidity

432elements and

σ = 450π

, however, italready ae ts theDG dis retisation too. This indi ates that largevaluesof

σ

prohibittheuseoffourth-order(or,indeed,anyhigh-order) ompositionmethods,aswell asexpli itRKmethods, su h as SSPRK. Instead, Ri hardson extrapolationbased onthe se ond-order

ompositionmethod maybeused sin etheyareun onditionallystablewithrespe ttothe ondu tivity

term. Similarlyto the ondu tion-free asewekilledthe

H(curl)

- onforming omputationsafter almost sixdaysalreadysigni antlymorethanwhattheDG omputationstake.

TableVI:Computational ostsoftheDGmethodforexample(26)with

σ = 60π

. Ameshof320elements isusedwithspatialpolynomialorders

p = 1, 2, 3

.

method #DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 3840 6.6817e-02 4526 0.0167 8s

p = 2

CO2 9600 8.4244e-03 7542 0.0100 113s

p = 2

GEX4 9600 8.4243e-03 22624 0.0100 341s

p = 3

CO4 19200 5.5619e-04 35192 0.0107 2012s

p = 3

GEX4 19200 5.5612e-04 31672 0.0071 1864s

p = 3

LEX4 19200 5.5612e-04 28154 0.0107 1623s

Table VII: Computational ostsof the

H(curl)

- onforming method forexample (26)with

σ = 60π

. A meshof320elementsisusedwithspatial polynomialorders

p = 1, 2, 3

.

method #DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 504 1.1789e-01 6253 0.0417 1s

p = 2

CO2 2388 1.2315e-02 56301 0.0250 82s

p = 2

GEX4 2388 1.2314e-02 166303 0.0250 247s

p = 3

CO4 6640 7.3357e-04 1717157 0.0127 40862s

p = 3

GEX4 6640 7.3358e-04 734732 0.0179 17472s

p = 3

LEX4 6640 7.3358e-04 653024 0.0268 15498s

5.2. Numeri al dispersionanalysis

Toinvestigatethedispersionanddissipationpropertiesofthefullydis retes hemes,we onsiderthe

semi-dis retesystem(11)with

σ = 0

and

j = 0

,

u

′

v

′



= A

u

_v



with

A =



0

I

−M

−1

ε

S

µ

0



,

(27)

10

2

10

4

10

6

10

−2

10

−1

p = 1, convergence in the L

2

norm

Degrees of Freedom

L

2

error

DG

H(curl)

order 2

10

2

10

4

10

6

10

0

p = 1, convergence in the energy norm

Degrees of Freedom

H(curl,

Ω

) / DG error

DG

H(curl)

order 1

10

2

10

4

10

−3

10

−2

10

−1

10

0

p = 2, convergence in the L

2

norm

Degrees of Freedom

L

2

error

DG

H(curl)

order 3

10

2

10

4

10

−1

10

0

10

1

p = 2, convergence in the energy norm

Degrees of Freedom

H(curl,

Ω

) / DG error

DG

H(curl)

order 2

10

2

10

3

10

4

10

−1

10

0

p = 3, convergence in the energy norm

Degrees of Freedom

H(curl,

Ω

) / DG error

DG

H(curl)

order 4

10

2

10

3

10

4

10

−3

10

−2

10

−1

p = 3, convergence in the L

2

norm

Degrees of Freedom

L

2

error

DG

H(curl)

order 3

Figure 6: Convergen e plotsin the

L

2

-norm (left olumn) and in the energy norm (right olumn) for

testexample(26)with

σ = 60π

. Inea hplotthe onvergen eratesoftheDGmethodandthe

H(curl)

- onformingmethod areshownalongwiththeexpe tedorderof onvergen e.

andassumeaplanewaveexa tsolution

TableVIII:Computational osts oftheDG method for example(26)with

σ = 450π

. An unstru tured meshof432elementsisusedwithspatial polynomialorders

p = 1, 2, 3

.

method # DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 5184 4.6168e-02 11878 0.00635 41s

p = 2

CO2 12960 8.1650e-03 19796 0.00381 422s

p = 2

GEX4 12960 8.1650e-03 59384 0.00381 1240s

p = 3

CO4 25920 8.5690e-04 222072 0.00170 17606s

p = 3

GEX4 25920 8.5671e-04 83134 0.00272 6618s

p = 3

LEX4 25920 8.5690e-04 73898 0.00136 5906s

TableIX:Computational ostsof the

H(curl)

- onformingmethod forexample(26)with

σ = 450π

. An unstru turedmeshof432elementsisusedwith spatialpolynomialorders

p = 1, 2, 3

.

method # DoF

L

2

_(Ω)

error #matve s

τ

CPUtime

p = 1

CO2 744 1.2498e-01 32049 0.02380 17s

p = 2

CO2 3420 1.3102e-02 163451 0.01429 938s

p = 2

GEX4 3420 1.3102e-02 490287 0.01429 2677s

p = 3

CO4 9360  >1e+07 0.01530 >5e+05s

p = 3

GEX4 9360  >1e+07 0.01020 >5e+05s

p = 3

LEX4 9360  >1e+07 0.01530 >5e+05s

with periodi boundary onditions and

E

ˆ

= 1

. In (28),

i

2

_{= −1}

,

ω

denotes the angular frequen y,

k = (k

x

, k

y

, k

z

)

T

isthewavenumber. Betweenthesequantitiesthe(exa t)dispersionrelation

ω

2

_{= k}

2

_/c

2

holdswith

k

2

_{= k}

2

x

+ k

2

y

+ k

z

2

andwith

c = 1/ (ε

r

µ

r

)

1/2

, whi histhespeedoflight.

Asarst step,weproje ttheexa tinitial onditions

E

(x, 0)

and

∂

t

E

(x, 0)

ontothenite-element spa e

E

_h

j

(0) = E(x, 0), ψ

j

Ω

,

j = 1 . . . N,

d

dt

E

j

h

(0) = ∂

t

E(x, 0), ψ

j

Ω

,

j = 1 . . . N.

(29)

We an now obtain the initial onditions for (27) through the relations

u

0

= u(0) = M

−1

ε

r

E

h

(0)

and

v

0

= v(0) = u

′

(0) = M

ε

−1

r

d

dt

E

h

(0)

. Thetime-exa tdis reteFouriermodeattimelevel

nτ

isthendened as

u

n

v

n



= ν

n

u

0

v

0



with

ν

n

_{= e}

−iω

h

nτ

_,

(30) where

ν

n

istheexa tampli ationfa torand

ω

h

isthesemi-dis retenumeri alfrequen y. Torstseetheimpa tofthespa edis retisationonly,we onsider thesemi-dis reteequation

M

ε

u

′′

+ S

µ

u = 0

(31)

withperiodi boundary onditionsandaplanewaveinitial ondition(28). Inthis ase,(31)isequivalent

tothedis retetime-harmoni Maxwelleigenvalueproblem

S

µ

u − ω

h

2

M

ε

u = 0

(32)

withperiodi boundary onditions. All semi-dis rete eigenvalues

ω

2

h

of (32)are real and non-negative, whi h entailsthat the spa e dis retisationimposes no dissipation. In Table X, we show thenumeri al

frequen iesofthespatialDGdis retisationfortheFouriermodewith

k

x

= 2π, k

y

= −2π, k

z

= 0

,i.e.with exa tangular frequen y

ω

ex

=

√

8π

. Thenumberof elementsforea hmesh is

N

el

= 5(

1

h

)

3

and inea h elementthelo alnumberofdegreesoffreedomis

1