Discontinuous Galerkin for the Wave Equation

71. Discontinuous Galerkin for the Wave Equation#

We solve the first order wave equation by a matrix-free explicit DG method:

\[\begin{eqnarray*} \frac{\partial p}{\partial t} & = & \operatorname{div} u \\ \frac{\partial u}{\partial t} & = & \nabla p \end{eqnarray*}\]

We obtain the ODE

\[\begin{eqnarray*} M_p \dot{p} & = & -B^T u \\ M_u \dot{u} & = & B p, \end{eqnarray*}\]

with mass-matrices \(M_p\) and \(M_u\), and the discretization matrix \(B\) for the gradient, and \(-B^T\) for the divergence. The exact solution of the ODE converses energy:

\[\begin{eqnarray*} \frac{d}{dt} \left( \tfrac{1}{2} p^T M_p p + \tfrac{1}{2} u^T M_u u \right) & = & p^T M_p \frac{d}{dt} p + u^T M_u \frac{d}{dt} u \\ & = & -p^T B^T u + u^T B p = 0 \end{eqnarray*}\]

A matrix \(A\) is called symplectic if

\[ J A \]

is symmetric, where \(J = \left( \begin{array}{cc}0 & I \\ -I & 0 \end{array} \right)\). The right-hand side matrix \(\left( \begin{array}{cc}0 & -B^T \\ B & 0 \end{array} \right)\) is symplectic.

We want to want to use DG methods for reason of cheaply invertible mass matrices.

The DG bilinear-form for the gradient is:

\[ b(p,v) = \sum_{T} \Big\{ \int_T \nabla p \, v + \int_{\partial T} (\{ p \} - p) \, v_n \, ds \Big\}, \]

where \(\{ p \}\) is the average of \(p\) on the facet.

This form is consistent for the grad as well as the div operator in the following sense: If we define

\[ b(p,v) = \sum_{T} \Big\{ \int_T \nabla p \, v + \alpha \int_{\partial T} (\{ p \} - p) \, v_n \, ds \Big\} \]

for some \(\alpha \in {\mathbb R}\), the form is consistent for the gradient. If \(p\) is globally continuous, then \(p = \{ p \}\), and the boundary term cancels out, even for discontinuous test-functions \(v\).

From integration by parts we obtain

\[\begin{align*} b(p,v) &= \sum_{T} \Big\{ -\int_T p \operatorname{div} \, v + \alpha \int_{\partial T} \{ p \} \, v_n \, ds + (1-\alpha) \int_{\partial T} p v_n \, ds \Big\} \\ &= \Big\{ -\int_T p \operatorname{div} \, v + \alpha \int_{\partial T} p \, [ v_n ] \, ds + (1-\alpha) \int_{\partial T} p v_n \, ds \Big\}, \end{align*}\]

with the normal-jump \([v_n] := v_l n_l + v_r n_r\). If we choose \(\alpha = 1\), then the boundary terms cancel for normal-continuous \(v \), even for discontinuous test-functions \(p\).

Literature: Hesthaven+Warbuton: Nodal Discontinuous Galerkin Methods

71.1. Testing the differential operators#

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import unit_square

mesh = Mesh(unit_square.GenerateMesh(maxh=0.05))

order = 5
fes_pT = L2(mesh, order=order, all_dofs_together=True)
fes_pF = FacetFESpace(mesh, order=order)

fes_p = fes_pT * fes_pF
fes_u = VectorL2(mesh, order=order, piola=True)

# matrix computing average on facet from both sides: 
traceop = fes_pT.TraceOperator(fes_pF, average=True) 

(p,pF),(q,qF) = fes_p.TnT()
u,v = fes_u.TnT()
n = specialcf.normal(2)

B = BilinearForm(trialspace=fes_p, testspace=fes_u)
B += grad(p)*v * dx 
B += (pF-p)*(v*n) * dx(element_vb=BND)
B.Assemble();

compute \(u = \nabla p\) via \(u = M_u^{-1} B p\)

gfu = GridFunction(fes_u)
gfp = GridFunction(fes_p)

gfpT, gfpF = gfp.components
gfpT.Set( exp(-100*( (x-0.5)**2+(y-0.5)**2)))
gfpF.vec.data = traceop * gfpT.vec   # setting the mean value

Mu = BilinearForm(u*v*dx).Assemble()
Mp = BilinearForm(p*q*dx + pF*qF*dx(element_vb=BND)).Assemble()
invu = Mu.mat.Inverse()  # block-diagonal
invp = Mp.mat.Inverse()  

gfu.vec[:] = invu @ B.mat * gfp.vec
Draw (gfpT)
Draw (gfu[0], mesh);

vice versa: compute \(p = -\operatorname{div} u\) via \(p = M_u^{-1} B^T u\)

gfu.Set( (exp(-100*( (x-0.5)**2+(y-0.5)**2)), 0))

hv = (B.mat.T * gfu.vec).Evaluate()
hv[fes_p.Range(0)] += traceop.T * hv[fes_p.Range(1)]
gfp.vec[:] = invp * hv
Draw (gfpT);

71.2. Efficient implementation:#

  • due to orthogonal basis on the reference element, mass-matrices a diagonal, FESpace provides efficient operator for \(M^{-1}\). They are not diagonal anymore if elements are curved, or for non-constant coefficients. Efficient inverse anyway.

  • due to Piola-mapping of the vector-variable, element matrices are the same for all elements, and need to be stored just for the reference element (\(F\) the Jacobian, and \(J = \det F\)):

\[ \int_T u \cdot \nabla p \, dx = \int_{\hat T} J^{-1} F \hat u \cdot F^{-T} \nabla \hat p \, J d \hat x = \int_{\hat T} \hat u \cdot \nabla \hat p \, d\hat x \]

Thus, the huge memory requirement for storing matrices is completely avoided. On highly parallel processors the matrix-free version gets faster and faster in comparison to matrix-based implementations, since more and more computations can be done per memory-transfer.

71.3. Solving the wave equation:#

from ngsolve import *
from netgen.occ import *

from ngsolve.webgui import Draw
from netgen.webgui import Draw as DrawGeo

dim = 2

if dim==2:
    rect = MoveTo(-1,-1).Rectangle(2,2).Face()
    circ = Circle((0.5,0), 0.2).Face()
    shape = rect-circ
    DrawGeo(shape)
    geo = OCCGeometry(shape, dim=2)
    h = 0.05

else:
    box = Box((-1,-1,-1), (1,1,0))
    sp = Sphere((0.5,0,0), 0.2)
    shape = box-sp
    geo = OCCGeometry(shape)
    h = 0.1
    
    
mesh = Mesh( geo.GenerateMesh(maxh=h))
mesh.Curve(3)
Draw(mesh);

A new component is the TraceOperator:

Space provide geometry-free linear operators mapping form the element space to the face space. Face values can be averaged, or are summed up.

order = 4
fes_p = L2(mesh, order=order, all_dofs_together=True)
fes_u = VectorL2(mesh, order=order, piola=True)
fes_tr = FacetFESpace(mesh, order=order)

print ("ndof_p = ", fes_p.ndof, "+", fes_tr.ndof, ", ndof_u =", fes_u.ndof)

traceop = fes_p.TraceOperator(fes_tr, average=True) 

gfu = GridFunction(fes_u)
gfp = GridFunction(fes_p)
gftr = GridFunction(fes_tr)

gfp.Interpolate( exp(-400*(x**2+y**2+z**2)))
gftr.vec.data = traceop * gfp.vec

ndof_p =  50985 + 25955 , ndof_u = 101970

if dim == 2:
    Draw (gfp, order=3)
else:
    gftr.vec.data = traceop * gfp.vec
    Draw (gftr, draw_vol=False, order=3);

p = fes_p.TrialFunction()
v = fes_u.TestFunction()
phat = fes_tr.TrialFunction()

n = specialcf.normal(mesh.dim)

We define bilinear-forms for the element-wise \(p\), and for the facet-wise \(p_{trace}\), the test-function is \(v\). Thanks to the co-variant mapping of \(v\), both forms are independent of element-geometry, and only one element matrix is calculated for the reference element(s):

Bel = BilinearForm(trialspace=fes_p, testspace=fes_u, geom_free = True)
Bel += grad(p)*v * dx -p*(v*n) * dx(element_boundary=True)
Bel.Assemble()

Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u, geom_free = True)
Btr += phat * (v*n) *dx(element_boundary=True)
Btr.Assemble();

Combine linear operators:

B = Bel.mat + Btr.mat @ traceop

print (Bel.mat)
print (Bel.mat.GetOperatorInfo())

Sum of
Scale a = 1
Print base-matrix
Scale b = 1
Print base-matrix

SumMatrix, h = 101970, w = 50985
  ConstantEBEMatrix (blocks = 1764, bs = 30x15), h = 101970, w = 50985
  ConstantEBEMatrix (blocks = 1635, bs = 30x15), h = 101970, w = 50985

Inverse mass matrices: either (block)diagonal, or operator application via sum factorization:

invmassp = fes_p.Mass(1).Inverse()
invmassu = fes_u.Mass(1).Inverse()

gfp.Interpolate( exp(-100*(x**2+y**2+z**2)))
gfu.vec[:] = 0

if dim == 2:
    scene = Draw (gfp, order=3, deformation=True);
else:
    scene = Draw (gftr, draw_vol=False, order=3);

t = 0
tend = 1
dt = 0.5 * h / (order+1)**2 / 2
print ("dt = ", dt)

cnt = 0
with TaskManager(): 
    while t < tend:
        t = t+dt
        gfu.vec.data += dt * invmassu @ B * gfp.vec
        gfp.vec.data -= dt * invmassp @ B.T * gfu.vec
        cnt = cnt+1
        if cnt%10 == 0:
            if dim == 3:
                gftr.vec.data = traceop * gfp.vec
            scene.Redraw()

dt =  0.0005

71.4. Eigenvalues of the discretized Laplace-operator#

Combining operators we obtain a discretization for the Laplace-operator:

\[ -\Delta_h = B^T M_u^{-1} B, \]

since \(B \approx \nabla\) and \(B^T \approx -\operatorname{div}\).

If the fe-space for \(u\) is too small, then \(\Delta_h\) gets additional zero - eigenvalues. We check that numerically by computing eigenvalues of \(-\Delta_h p = \lambda M_p p\).

Laplace = B.T @ invmassu @ B
mass = fes_p.Mass(1)

fes_p = L2(mesh, order=order, all_dofs_together=True, dgjumps=True)

if fes_p.ndof > 100000:
    raise Exception("come back with coarser mesh, ndof=", fes_p.ndof)

bfpre = BilinearForm(fes_p)
p,q = fes_p.TnT()
bfpre += grad(p)*grad(q)*dx + p*q*dx
bfpre += 1/h*(p-p.Other())*(q-q.Other()) * dx(skeleton=True)

with TaskManager():
    bfpre.Assemble()
    inv = bfpre.mat.Inverse(inverse="sparsecholesky")

from ngsolve.solvers import PINVIT
with TaskManager():
    eigensys = PINVIT(Laplace, mass, inv, num=8, maxit=50)

0 : [0.0019610691782318895, 12.530097553318297, 30.071959767920585, 44.410124391680235, 111.27047830136578, 125.87981867431681, 373.26816082154824, 519.634322447681]
1 : [0.0014151881769126334, 5.225415932359963, 7.3575826246171525, 13.211641321797787, 27.35157384116149, 34.456917882418836, 60.60766731363219, 65.05584265387385]
2 : [0.000975069231975477, 3.5262200551942495, 3.544005806596869, 6.649244667185429, 12.718358640298115, 16.24648897104735, 22.07147983053837, 26.227146148469842]
3 : [0.0006968604580658802, 2.9254019154065816, 3.0620927099610853, 5.961604255858556, 9.910196187032746, 13.50285365901096, 16.230984021795205, 18.12233489466299]
4 : [0.0005608189260820888, 2.7443252539524443, 2.8799540121040303, 5.661322991398223, 9.591059349873706, 12.634629366864042, 15.292915463017845, 16.520527971728463]
5 : [0.00046579731589452083, 2.6413663286264972, 2.7710907429403915, 5.461482941176028, 9.485224986561821, 12.128775079123026, 14.835902724967195, 15.732350771987893]
6 : [0.00039062189751259827, 2.560229600324007, 2.6950553221812146, 5.31824498415214, 9.409109730306385, 11.762861788003466, 14.495030612220166, 15.155592893668446]
7 : [0.0003300282622792194, 2.495367935557823, 2.63842291433198, 5.20780681484817, 9.345895678381895, 11.493777614461361, 14.239807143091744, 14.71776355719705]
8 : [0.0002805184506089352, 2.4432129136664176, 2.595150243934314, 5.1216388962948525, 9.291695345491545, 11.288095429922729, 14.035485037980463, 14.368406226609025]
9 : [0.00023972076349291322, 2.401035962821033, 2.5610751331099992, 5.052814774189954, 9.244676238525523, 11.127545194369905, 13.871280492089983, 14.088721043899092]

10 : [0.00020588732467636457, 2.366619149955629, 2.5339243506803264, 4.997389855877352, 9.20365998231166, 11.000218971727886, 13.734454124722511, 13.8584017339499]
11 : [0.00017767518756896892, 2.338361829670346, 2.5119087716222763, 4.952190394289993, 9.167616792450728, 10.89712959535109, 13.616459580820948, 13.671813205976868]
12 : [0.0001540423652007051, 2.314990627788917, 2.4939354627155734, 4.9151063538243704, 9.135867883250839, 10.81277048322605, 13.489223186970888, 13.541000262646078]
13 : [0.0001341620838840585, 2.2955553680857257, 2.4790891037108436, 4.884414799149257, 9.10776084785462, 10.742685293569647, 13.356336185945942, 13.454646768918472]
14 : [0.00011737517461220467, 2.279297361426251, 2.4667589967148507, 4.858878072701796, 9.082829128122407, 10.683958900723692, 13.23831755283918, 13.383674940455194]
15 : [0.00010314731601040092, 2.265631259380973, 2.4564296169864144, 4.83748608482573, 9.060646165949413, 10.63425284613956, 13.135544801782142, 13.323112581124242]
16 : [9.104558928875343e-05, 2.25408650260681, 2.4477326693537402, 4.81947364148078, 9.040869847901186, 10.591880670150983, 13.045516460669337, 13.270523895611174]
17 : [8.071484085630474e-05, 2.244290093649316, 2.4403593426127297, 4.8042212186230335, 9.023201644855128, 10.555506325896008, 12.966537965132005, 13.22467522592536]
18 : [7.186402839649455e-05, 2.235940267447839, 2.4340769810501515, 4.791238438812808, 9.007383664486253, 10.524090725741827, 12.896798633220365, 13.184333532834613]
19 : [6.425308585830436e-05, 2.22879289299042, 2.428692738661899, 4.780133588365316, 8.993197076186522, 10.496813626792058, 12.835085726278352, 13.14874821596208]
20 : [5.768348257520841e-05, 2.222649362107733, 2.424054906734523, 4.770584039884224, 8.980445622046407, 10.473004799002306, 12.780195570063217, 13.117144164013688]
21 : [5.1991081281590905e-05, 2.217346510800938, 2.4200394081047225, 4.762336368277239, 8.968965438349066, 10.45213431639568, 12.731272262274818, 13.089011985370156]
22 : [4.7039084180760346e-05, 2.212750942961908, 2.4165449660720277, 4.755173965781026, 8.9586067011406, 10.433754477326941, 12.687490316751461, 13.063836164171157]

23 : [4.271424589078034e-05, 2.208751765162773, 2.4134898843027406, 4.748929737271881, 8.949245458839577, 10.4175102423239, 12.648234029103469, 13.041256864928334]
24 : [3.89215341158034e-05, 2.205257913857816, 2.410805305541306, 4.743455863587409, 8.94076706856538, 10.403094016632934, 12.612919030581711, 13.020919263218294]
25 : [3.558218248694774e-05, 2.2021929581480886, 2.4084363909829936, 4.738640307100792, 8.933077112268858, 10.39026124479051, 12.581093904109482, 13.002562734784519]
26 : [3.262967697767145e-05, 2.199493889099703, 2.406335619782723, 4.734380793340841, 8.926087570170047, 10.378795547355319, 12.552334043399304, 12.985936286224]
27 : [3.0008807435174038e-05, 2.1971073804346646, 2.404465542771998, 4.73060101988937, 8.919726032121678, 10.368524608977475, 12.5263029503987, 12.970847530430802]
28 : [2.7672672373998064e-05, 2.194989319982176, 2.4027929355688475, 4.72722934888614, 8.91392459510606, 10.35929308658514, 12.50268549394198, 12.957114675391006]
29 : [2.5582270273592196e-05, 2.193102085513224, 2.4012918070189944, 4.7242130272305465, 8.908627319698214, 10.350977172422636, 12.48122720423944, 12.944593439617545]
30 : [2.370427095165695e-05, 2.191414459878113, 2.399938646795823, 4.721501289108385, 8.903781506645101, 10.343463326093982, 12.461690048094432, 12.933149361339717]
31 : [2.2010896203963964e-05, 2.1898996194818907, 2.3987151275149, 4.719057113665783, 8.899343599496435, 10.336661033846253, 12.443879138747146, 12.922672717232004]
32 : [2.0478260136062454e-05, 2.188535241919751, 2.3976044142282626, 4.716844106608183, 8.895272415485913, 10.330485968322035, 12.427611968155642, 12.913062163107492]
33 : [1.908640095139863e-05, 2.187302004620112, 2.3965933755399456, 4.71483589108621, 8.891533759521552, 10.324870904893558, 12.412737588184285, 12.904233187221493]
34 : [1.7818041102450526e-05, 2.1861837778986075, 2.395669808147888, 4.7130060616738065, 8.888095207735471, 10.319752298441507, 12.399114237596686, 12.896108100713224]
35 : [1.66586928937571e-05, 2.185166497058788, 2.3948241495056495, 4.71133553100744, 8.884929688259712, 10.315079492339123, 12.386623799768886, 12.888621030718166]
36 : [1.5595728533079688e-05, 2.18423838039716, 2.3940474383307215, 4.709804898287155, 8.882011475494716, 10.310803965065515, 12.375154905883043, 12.88171150967329]

37 : [1.4618515865450911e-05, 2.183389079888074, 2.3933325953744164, 4.708400128656904, 8.879318957385712, 10.306887028688342, 12.364614257300868, 12.875327609833588]
38 : [1.3717718213876082e-05, 2.1826098912180214, 2.3926729483373967, 4.707106777194164, 8.876831564457603, 10.303291182707614, 12.35491348529469, 12.869421626528611]
39 : [1.2885435632537004e-05, 2.1818931148861314, 2.3920631693492993, 4.705914354197892, 8.874531902337274, 10.299986508260899, 12.345978246334381, 12.863952181901217]
40 : [1.2114676295632756e-05, 2.181232242249903, 2.391498216683536, 4.704811954121162, 8.872403399581287, 10.296943700067796, 12.337737807181698, 12.85888119036187]
41 : [1.139948977370347e-05, 2.180621476450127, 2.3909740152140686, 4.703791606539386, 8.870431946554765, 10.29413934900636, 12.330132438195784, 12.854175340848009]
42 : [1.0734567447098367e-05, 2.1800558885312946, 2.390486696201572, 4.702844959397771, 8.868604095873494, 10.291550321162465, 12.323105089631335, 12.849803912039741]
43 : [1.0115360984259937e-05, 2.179531059836302, 2.3900330938111454, 4.701965856437029, 8.866908317340382, 10.289158099964776, 12.316607398202386, 12.845739826148193]
44 : [9.537780424779887e-06, 2.1790432081308464, 2.389610192363181, 4.701147808411721, 8.865333623942625, 10.28694425363793, 12.310593044331622, 12.841958081125302]
45 : [8.998295148530828e-06, 2.1785889218796943, 2.389215495994879, 4.700385983783346, 8.863870527880664, 10.28489398593868, 12.305022580235324, 12.838436466564913]
46 : [8.493706188243065e-06, 2.178165259325979, 2.3888466170569984, 4.699675268611956, 8.86250999494075, 10.282992486828169, 12.299858178682207, 12.835154478544357]
47 : [8.021229520772934e-06, 2.17776955183001, 2.3885015607039146, 4.6990118129178855, 8.861244167695043, 10.281227820654832, 12.295067443190675, 12.832093747763125]
48 : [7.5783248287882935e-06, 2.177399480200587, 2.3881784080586255, 4.698391531262507, 8.860065573511271, 10.279587994787759, 12.290619320725037, 12.8292373461617]
49 : [7.1627624544712525e-06, 2.177052929397117, 2.3878755443490745, 4.6978113120119, 8.858967666764388, 10.278063296596175, 12.286487037562203, 12.82656996756256]

lam, vecs = eigensys

gfp = GridFunction(fes_p)
gfp.vec.data = vecs[2]
Draw (gfp)

BaseWebGuiScene

help (PINVIT)

Help on function PINVIT in module ngsolve.eigenvalues:

PINVIT(mata, matm, pre, num=1, maxit=20, printrates=True, GramSchmidt=True)
    preconditioned inverse iteration

71.5. Exercise:#

  • check conservation of energy over time

  • design some wave-guide:
    A thin strip with (e.g. 2-times) higher refraction index, i.e. the coefficient in the coefficient of \(M_p\). Start with a peak in the wave-guide, it should travel within the strip.

Contents