# 5.5.2 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*}

Using DG discretization in space we obtain the ODE system \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 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.

The computational advantages of a proper version of DG is:

• universal element-matrices for the differntial operator $$B$$

• cheaply invertible mass matrices (orthogonal polynomials, sum-factorization)

[1]:

from ngsolve import *
from netgen.occ import *
from time import time

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

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

[2]:

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 =  650405 + 585960 , ndof_u = 1951215

[3]:

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

n = specialcf.normal(mesh.dim)

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

[4]:

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();

B = Bel.mat + Btr.mat @ traceop

[5]:

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

[6]:

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

scene = Draw (gftr, draw_vol=False, order=3);

t = 0
dt = 0.3 * h / (order+1)**2
tend = 1

op1 = dt * invmassu @ B
op2 = dt * invmassp @ B.T

cnt = 0
while t < tend:
t = t+dt
gfu.vec.data += op1 * gfp.vec
gfp.vec.data -= op2 * gfu.vec
cnt = cnt+1
if cnt%10 == 0:
gftr.vec.data = traceop * gfp.vec
scene.Redraw()


## Time-stepping on the device¶

[7]:

try:
import ngsolve.ngscuda
except:
print ("Sorry, no cuda device")

Sorry, no cuda device

[8]:

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

scene = Draw (gftr, draw_vol=False, order=3);

t = 0
dt = 0.5 * h / (order+1)**2 / 2
tend = 0.1

op1 = (dt * invmassu @ B).CreateDeviceMatrix()
op2 = (dt * invmassp @ B.T).CreateDeviceMatrix()

devu = gfu.vec.CreateDeviceVector(copy=True)
devp = gfp.vec.CreateDeviceVector(copy=True)

cnt = 0
while t < tend:
t = t+dt
devu += op1 * devp
devp -= op2 * devu
cnt = cnt+1
if cnt%10 == 0:
gfp.vec.data = devp
gftr.vec.data = traceop * gfp.vec
scene.Redraw()

[9]:

print (op1.GetOperatorInfo())

ProductMatrix, h = 1951215, w = 650405
ScaleMatrix, scale = 0.001, h = 1951215, w = 1951215
SumMatrix, h = 1951215, w = 1951215
SumMatrix, h = 1951215, w = 1951215
BlockDiagonalMatrixSoA (bs = 3x3), h = 1951215, w = 1951215
ProductMatrix, h = 1951215, w = 1951215
ProductMatrix, h = 1951215, w = 18000
Transpose, h = 1951215, w = 18000
ConstantEBEMatrix (bs = 125x35), h = 18000, w = 1951215
BlockDiagonalMatrixSoA (bs = 3x3), h = 18000, w = 18000
ConstantEBEMatrix (bs = 125x35), h = 18000, w = 1951215
ProductMatrix, h = 1951215, w = 1951215
ProductMatrix, h = 1951215, w = 24000
Transpose, h = 1951215, w = 24000
ConstantEBEMatrix (bs = 125x35), h = 24000, w = 1951215
BlockDiagonalMatrixSoA (bs = 3x3), h = 24000, w = 24000
ConstantEBEMatrix (bs = 125x35), h = 24000, w = 1951215
SumMatrix, h = 1951215, w = 650405
SumMatrix, h = 1951215, w = 650405
ConstantEBEMatrix (bs = 105x35), h = 1951215, w = 650405
ConstantEBEMatrix (bs = 105x35), h = 1951215, w = 650405
ProductMatrix, h = 1951215, w = 650405
SumMatrix, h = 1951215, w = 585960
ConstantEBEMatrix (bs = 105x60), h = 1951215, w = 585960
ConstantEBEMatrix (bs = 105x60), h = 1951215, w = 585960
ProductMatrix, h = 585960, w = 650405
N4ngla14DiagonalMatrixIdEE, h = 585960, w = 585960
SumMatrix, h = 585960, w = 650405
ConstantEBEMatrix (bs = 60x35), h = 585960, w = 650405
ConstantEBEMatrix (bs = 60x35), h = 585960, w = 650405


[10]:

ts = time()
steps = 10
for i in range(steps):
devu += op1 * devp
devp -= op2 * devu
te = time()
print ("ndof = ", gfp.space.ndof, "+", gfu.space.ndof, ", time per step =", (te-ts)/steps)

ndof =  650405 + 1951215 , time per step = 0.15075998306274413


On the A-100:

ndof = 651035 + 1953105 , time per step = 0.0023579835891723634

## Time-domain PML¶

PML (perfectly matched layers) is a method for approximating outgoing waves on a truncated domain. Its time domain version leads to additional field variables in the PML region, which are coupled via time-derivatives only.

[11]:

from ring_resonator_import import *
from ngsolve import *
from ngsolve.webgui import Draw

define envelope

WARNING: kwarg 'orderinner' is an undocumented flags option for class <class 'ngsolve.comp.HDiv'>, maybe there is a typo?

[12]:

Draw (mesh);

[13]:

scene = Draw (gfu.components[0], order=3, min=-0.05, max=0.05, autoscale=False)

tend = 15
i = 0
while t < tend:
w.data = -fullB.T * gfu.vec
w.data += Envelope(t)*Lsrc.vec
w.data -= sigma * dampingp * gfu.vec
w.data -= emb_p @ Bstab.mat.T * gfstab.vec

gfu.vec.data += tau * invp * w

w.data = fullB * gfu.vec
w.data -= sigma * dampingu * gfu.vec
gfu.vec.data += tau * invu * w

gfstab.vec.data += tau * Mstabinv@Bstab.mat@emb_p.T * gfu.vec

t += tau
i += 1
if i%20 == 0:
scene.Redraw()


The differential operators and coupling terms to the pml - variables are represented via linear operators;

[14]:

print (fullB.GetOperatorInfo())
print (dampingp.GetOperatorInfo())

SumMatrix, h = 64960, w = 64960
SumMatrix, h = 64960, w = 64960
ConstantEBEMatrix (bs = 20x15), h = 64960, w = 64960
ConstantEBEMatrix (bs = 20x15), h = 64960, w = 64960
EmbeddedTransposeMatrix, h = 64960, w = 64960
EmbeddedMatrix, h = 64960, w = 13920
ProductMatrix, h = 18560, w = 13920
SumMatrix, h = 18560, w = 7165
ConstantEBEMatrix (bs = 20x15), h = 18560, w = 7165
ConstantEBEMatrix (bs = 20x15), h = 18560, w = 7165
SumMatrix, h = 7165, w = 13920
ConstantEBEMatrix (bs = 15x15), h = 7165, w = 13920
ConstantEBEMatrix (bs = 15x15), h = 7165, w = 13920

SumMatrix, h = 64960, w = 64960
SumMatrix, h = 64960, w = 64960
EmbeddedTransposeMatrix, h = 64960, w = 64960
EmbeddedMatrix, h = 64960, w = 13920
N6ngcomp9ApplyMassE, h = 13920, w = 13920
ProductMatrix, h = 64960, w = 64960
EmbeddedMatrix, h = 64960, w = 13920
N6ngcomp9ApplyMassE, h = 13920, w = 13920
SumMatrix, h = 13920, w = 64960
ScaleMatrix, scale = 2, h = 13920, w = 64960
N4ngla18EmbeddingTransposeE, h = 13920, w = 64960
N4ngla18EmbeddingTransposeE, h = 13920, w = 64960
EmbeddedTransposeMatrix, h = 64960, w = 64960
EmbeddedMatrix, h = 64960, w = 13920
N6ngcomp9ApplyMassE, h = 13920, w = 13920


[ ]: