Discontinuous Galerkin for the Wave Equation

69. 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

69.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);

69.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.

69.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 =  50835 + 25880 , ndof_u = 101670

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 = 101670, w = 50835
  ConstantEBEMatrix (bs = 30x15), h = 101670, w = 50835
  ConstantEBEMatrix (bs = 30x15), h = 101670, w = 50835

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 = 20
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

---------------------------------------------------------------------------
KeyboardInterrupt                         Traceback (most recent call last)
Cell In[12], line 16
     14 cnt = 0
     15 with TaskManager(): 
---> 16     while t < tend:
     17         t = t+dt
     18         gfu.vec.data += dt * invmassu @ B * gfp.vec

KeyboardInterrupt:

69.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 > 50000:
    raise Exception("come back with coarser mesh")

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)

lam, vecs = eigensys

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

help (PINVIT)

69.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