71. Discontinuous Galerkin for the Wave Equation#
We solve the first order wave equation by a matrix-free explicit DG method:
We obtain the ODE
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:
A matrix \(A\) is called symplectic if
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:
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
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
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\)):
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 = 51015 + 25970 , ndof_u = 102030
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 = 102030, w = 51015
ConstantEBEMatrix (blocks = 1765, bs = 30x15), h = 102030, w = 51015
ConstantEBEMatrix (blocks = 1636, bs = 30x15), h = 102030, w = 51015
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 24
20 cnt = cnt+1
21 if cnt%10 == 0:
22 if dim == 3:
23 gftr.vec.data = traceop * gfp.vec
---> 24 scene.Redraw()
File /Applications/Netgen.app/Contents/Resources/lib/python3.14/site-packages/netgen/webgui.py:259, in WebGLScene.Redraw(self, *args, **kwargs)
257 self.args = new_scene.args
258 self.kwargs = new_scene.kwargs
--> 259 super().Redraw()
File /Library/Frameworks/Python.framework/Versions/3.14/lib/python3.14/site-packages/webgui_jupyter_widgets/widget.py:81, in BaseWebGuiScene.Redraw(self)
79 def Redraw(self):
80 self.encoding='binary'
---> 81 self.widget.value = self.GetData(set_minmax=False)
File /Applications/Netgen.app/Contents/Resources/lib/python3.14/site-packages/netgen/webgui.py:266, in WebGLScene.GetData(self, set_minmax)
264 d = None
265 if type(self.obj) in _registered_draw_types:
--> 266 d = _registered_draw_types[typ](self.obj, self.args, self.kwargs)
267 else:
268 import inspect
File /Applications/Netgen.app/Contents/Resources/lib/python3.14/site-packages/ngsolve/webgui.py:163, in GetData(obj, args, kwargs)
161 if 'intpoints' not in kwargs:
162 kwargs['intpoints'] = None
--> 163 d = BuildRenderData(mesh, cf, order=kwargs['order'], draw_surf=kwargs['draw_surf'], draw_vol=kwargs['draw_vol'], intpoints=kwargs['intpoints'], deformation=kwargs['deformation'], regions=regions, objects=kwargs['objects'], nodal_p1=kwargs['nodal_p1'], settings=kwargs['settings'])
165 if isinstance(cf, ngs.GridFunction) and len(cf.vecs)>1:
166 # multidim gridfunction - generate data for each component
167 gf = ngs.GridFunction(cf.space)
File /Applications/Netgen.app/Contents/Resources/lib/python3.14/site-packages/ngsolve/webgui.py:390, in BuildRenderData(mesh, func, order, draw_surf, draw_vol, intpoints, deformation, regions, objects, nodal_p1, encoding, settings)
387 timer2map.Stop()
388 pmat = cf_surf(pts)
--> 390 pmima = updatePMinMax(pmat)
392 timermult.Start()
393 pmat = pmat.reshape(-1, og+1, cf.dim)
File /Applications/Netgen.app/Contents/Resources/lib/python3.14/site-packages/ngsolve/webgui.py:10, in updatePMinMax(pmat, pmima)
9 def updatePMinMax( pmat, pmima=None ):
---> 10 pmima_new = [ngs.Vector(pmat[:,i], copy=False).MinMax(ignore_inf=True) for i in range(3)]
12 if pmima is not None:
13 for i in range(3):
KeyboardInterrupt:
71.4. Eigenvalues of the discretized Laplace-operator#
Combining operators we obtain a discretization for the Laplace-operator:
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)
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.