10. 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)

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)
disc = Cylinder( Axes((0.5,0,0), X, Y), r=0.4, h=0.2)
shape = box-disc
geo = OCCGeometry(shape)

h = 0.1
        
mesh = Mesh( geo.GenerateMesh(maxh=h))
mesh.Curve(3)
Draw(mesh);

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 =  574630 + 521160 , ndof_u = 1723890

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.

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

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

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

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

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

cnt = 0
with TaskManager(): 
    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()

10.1. Time-stepping on the device

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

Sorry, no cuda device

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
with TaskManager(): 
    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()

print (op1.GetOperatorInfo())

ProductMatrix, h = 1723890, w = 574630
  ScaleMatrix, scale = 0.001, h = 1723890, w = 1723890
    SumMatrix, h = 1723890, w = 1723890
      SumMatrix, h = 1723890, w = 1723890
        BlockDiagonalMatrixSoA (bs = 3x3), h = 1723890, w = 1723890
        ProductMatrix, h = 1723890, w = 1723890
          ProductMatrix, h = 1723890, w = 33375
            Transpose, h = 1723890, w = 33375
              ConstantEBEMatrix (bs = 125x35), h = 33375, w = 1723890
            BlockDiagonalMatrixSoA (bs = 3x3), h = 33375, w = 33375
          ConstantEBEMatrix (bs = 125x35), h = 33375, w = 1723890
      ProductMatrix, h = 1723890, w = 1723890
        ProductMatrix, h = 1723890, w = 31125
          Transpose, h = 1723890, w = 31125
            ConstantEBEMatrix (bs = 125x35), h = 31125, w = 1723890
          BlockDiagonalMatrixSoA (bs = 3x3), h = 31125, w = 31125
        ConstantEBEMatrix (bs = 125x35), h = 31125, w = 1723890
  SumMatrix, h = 1723890, w = 574630
    SumMatrix, h = 1723890, w = 574630
      ConstantEBEMatrix (bs = 105x35), h = 1723890, w = 574630
      ConstantEBEMatrix (bs = 105x35), h = 1723890, w = 574630
    ProductMatrix, h = 1723890, w = 574630
      SumMatrix, h = 1723890, w = 521160
        ConstantEBEMatrix (bs = 105x60), h = 1723890, w = 521160
        ConstantEBEMatrix (bs = 105x60), h = 1723890, w = 521160
      ProductMatrix, h = 521160, w = 574630
        N4ngla14DiagonalMatrixIdEE, h = 521160, w = 521160
        SumMatrix, h = 521160, w = 574630
          ConstantEBEMatrix (bs = 60x35), h = 521160, w = 574630
          ConstantEBEMatrix (bs = 60x35), h = 521160, w = 574630

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 =  574630 + 1723890 , time per step = 0.03982951641082764

On the A-100:

ndof =  651035 + 1953105 , time per step = 0.0023579835891723634

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

Formulation of B. Kapidani, J. Schöberl, arxiv

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

from netgen.geom2d import SplineGeometry

geo = SplineGeometry()

xneg  =-0.43
xpos  = 0.43
yneg  =-0.48
ypos  = 0.48
wslab = 0.04
cringx = 0.0
cringy = 0.0
rring = 0.4
gap   = 0.005

pntx = [xneg,xpos]
pnty = [yneg,-rring-gap-wslab,-rring-gap,rring+gap,rring+gap+wslab,ypos]

pts = [geo.AddPoint(xi,yi) for yi in pnty for xi in pntx]

### geometry #######################################################
geo.Append (["line", pts[0], pts[1]], leftdomain=1, rightdomain=0)
geo.Append (["line", pts[1], pts[3]], leftdomain=1, rightdomain=0)
geo.Append (["line", pts[3], pts[2]], leftdomain=1, rightdomain=2)
geo.Append (["line", pts[2], pts[0]], leftdomain=1, rightdomain=0)

geo.Append (["line", pts[3], pts[5]], leftdomain=2, rightdomain=0,bc="normal_wg_rightbottom")
geo.Append (["line", pts[5], pts[4]], leftdomain=2, rightdomain=3)
geo.Append (["line", pts[4], pts[2]], leftdomain=2, rightdomain=0,bc="normal_wg_leftbottom")

geo.Append (["line", pts[5], pts[7]], leftdomain=3, rightdomain=0)
geo.Append (["line", pts[7], pts[6]], leftdomain=3, rightdomain=4)
geo.Append (["line", pts[6], pts[4]], leftdomain=3, rightdomain=0)

geo.Append (["line", pts[7], pts[9]], leftdomain=4, rightdomain=0,bc="normal_wg_righttop")
geo.Append (["line", pts[9], pts[8]], leftdomain=4, rightdomain=5)
geo.Append (["line", pts[8], pts[6]], leftdomain=4, rightdomain=0,bc="normal_wg_lefttop")

geo.Append (["line", pts[9], pts[11]], leftdomain=5, rightdomain=0)
geo.Append (["line", pts[11], pts[10]], leftdomain=5, rightdomain=0)
geo.Append (["line", pts[10], pts[8]], leftdomain=5, rightdomain=0)

geo.AddCircle(c=(cringx,cringy), r=rring, leftdomain=6, rightdomain=3)
geo.AddCircle(c=(cringx,cringy), r=rring-wslab, leftdomain=7, rightdomain=6)

for i in (1,3,5,7):
    geo.SetMaterial(i, "air")
for i in (2,4,6):
    geo.SetMaterial(i, "eps_nine")
data = geo.CreatePML(0.05)
normals = data["normals"]


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

eps_r = {"air" : 1, "eps_nine" : 3**3}

for mat in mesh.GetMaterials():
    if mat.startswith("pml_normal_wg"):
        eps_r[mat] = eps_r["eps_nine"]

for mat in mesh.GetMaterials():
    if mat not in eps_r:
        eps_r[mat] = eps_r["air"]

### Parameters for Source field ######
wavelength = 1.542
fcen       = 5/wavelength
df         = 0.1
tpeak      = 1
order = 3

mesh.Curve(order)
fes_facet = FacetFESpace(mesh, order=order+1)
gfsource = GridFunction(fes_facet)

source_cf =  sin( (pi/wslab)*(y-pnty[3]) ) 
gfsource.Set(source_cf,definedon=mesh.Boundaries("normal_wg_lefttop"))

fes_u = VectorL2(mesh, order=order, piola=True, order_policy=ORDER_POLICY.CONSTANT)
fes_p = L2(mesh, order=order+1, all_dofs_together=True, order_policy=ORDER_POLICY.CONSTANT)
fes_tr = FacetFESpace(mesh, order=order+1)
fes_hdiv = HDiv(mesh, order=order+1, orderinner=1)

n = specialcf.normal(2) 

p,q = fes_p.TnT()
u,v = fes_u.TnT()
ptr,qtr = fes_tr.TnT()
uhdiv,vhdiv = fes_hdiv.TnT()

Bel = BilinearForm(grad(p)*v*dx - p*(v*n)*dx(element_boundary=True), geom_free=True).Assemble()
Btr = BilinearForm(0.5*ptr*(n*v)*dx(element_boundary=True), geom_free=True).Assemble()
Bstab = BilinearForm(p*(vhdiv*n)*dx(element_boundary=True),  geom_free=True).Assemble()


nvec = { mat : ((normals[mat][0], normals[mat][1]) if mat in normals else (0,0)) for mat in mesh.GetMaterials() }

cfn = CF( [CF(nvec[mat]) for mat in mesh.GetMaterials()])
cft = CF( ( cfn[1], -cfn[0] ) )

pml1d = mesh.Materials("pml_default.*|pml_normal.*")
eps = CF([eps_r[mat] for mat in mesh.GetMaterials()])


fes = fes_p*fes_p*fes_u*fes_u*fes_hdiv
emb_p, emb_phat, emb_u, emb_uhat, emb_hdiv = fes.embeddings

# gradient operator
traceop = fes_p.TraceOperator(fes_tr, False)
fullB = emb_u @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T + emb_hdiv@Bstab.mat@emb_p.T

# mass matrices
invmassp = fes_p.Mass(eps).Inverse()
invmassu = fes_u.Mass(Id(mesh.dim)).Inverse()
Mstab = BilinearForm(uhdiv*vhdiv*dx, diagonal=True).Assemble()
Mstabinv = Mstab.mat.Inverse()


invp = emb_p @ invmassp @ emb_p.T + emb_phat @ invmassp @ emb_phat.T
invu = emb_u @ invmassu @ emb_u.T + emb_uhat @ invmassu @ emb_uhat.T + emb_hdiv@Mstabinv@emb_hdiv.T


# damping matrices
dampp1 = fes_p.Mass (eps, definedon=pml1d)
dampp2 = fes_p.Mass (eps, definedon=mesh.Materials("pml_corner"))
dampu1 = fes_u.Mass (OuterProduct(cfn,cfn), definedon=pml1d)
dampu2 = fes_u.Mass (OuterProduct(cft,cft), definedon=pml1d)

dampingu = emb_u @ dampu1 @ emb_u.T + (-emb_u + emb_uhat) @ dampu2 @ (emb_u.T + emb_uhat.T)
dampingp = emb_p @ dampp1 @ emb_p.T + emb_p @ dampp2 @ (2*emb_p.T-emb_phat.T) + emb_phat @ dampp2 @ emb_p.T


# source term
Lsrc  = LinearForm(gfsource*q*dx(element_boundary=True)).Assemble()
srcvec = emb_p * (invmassp*Lsrc.vec).Evaluate()

# time-envelope:
def Envelope(t):
    if abs((t-tpeak)/tpeak) < 1:
        return (2*exp(1)/sqrt(pi))*sin(2*pi*fcen*t)*exp (-1/(1-((t-tpeak)/tpeak)**2))
    else:
        return 0

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

t = 0
tend = 12
tau = 2e-4
i = 0
sigma = 10   # pml damping parameter

op1 = invp@(-fullB.T-sigma*dampingp) 
op2 = invu@(fullB-sigma*dampingu)
with TaskManager(): 
    while t < tend:

        gfu.vec.data += tau*Envelope(t)*srcvec
        gfu.vec.data += tau*op1*gfu.vec
        gfu.vec.data += tau*op2*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;

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

SumMatrix, h = 72125, w = 72125
  EmbeddedTransposeMatrix, h = 72125, w = 72125
    EmbeddedMatrix, h = 72125, w = 13920
      SumMatrix, h = 18560, w = 13920
        SumMatrix, h = 18560, w = 13920
          ConstantEBEMatrix (bs = 20x15), h = 18560, w = 13920
          ConstantEBEMatrix (bs = 20x15), h = 18560, 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
  EmbeddedTransposeMatrix, h = 72125, w = 72125
    EmbeddedMatrix, h = 72125, w = 13920
      SumMatrix, h = 7165, w = 13920
        ConstantEBEMatrix (bs = 15x15), h = 7165, w = 13920
        ConstantEBEMatrix (bs = 15x15), h = 7165, w = 13920

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

print (op1.GetOperatorInfo())

ProductMatrix, h = 72125, w = 72125
  SumMatrix, h = 72125, w = 72125
    EmbeddedTransposeMatrix, h = 72125, w = 72125
      EmbeddedMatrix, h = 72125, w = 13920
        N6ngcomp9ApplyMassE, h = 13920, w = 13920
    EmbeddedTransposeMatrix, h = 72125, w = 72125
      EmbeddedMatrix, h = 72125, w = 13920
        N6ngcomp9ApplyMassE, h = 13920, w = 13920
  SumMatrix, h = 72125, w = 72125
    ScaleMatrix, scale = -1, h = 72125, w = 72125
      Transpose, h = 72125, w = 72125
        SumMatrix, h = 72125, w = 72125
          EmbeddedTransposeMatrix, h = 72125, w = 72125
            EmbeddedMatrix, h = 72125, w = 13920
              SumMatrix, h = 18560, w = 13920
                SumMatrix, h = 18560, w = 13920
                  ConstantEBEMatrix (bs = 20x15), h = 18560, w = 13920
                  ConstantEBEMatrix (bs = 20x15), h = 18560, 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
          EmbeddedTransposeMatrix, h = 72125, w = 72125
            EmbeddedMatrix, h = 72125, w = 13920
              SumMatrix, h = 7165, w = 13920
                ConstantEBEMatrix (bs = 15x15), h = 7165, w = 13920
                ConstantEBEMatrix (bs = 15x15), h = 7165, w = 13920
    ScaleMatrix, scale = 10, h = 72125, w = 72125
      SumMatrix, h = 72125, w = 72125
        SumMatrix, h = 72125, w = 72125
          EmbeddedTransposeMatrix, h = 72125, w = 72125
            EmbeddedMatrix, h = 72125, w = 13920
              N6ngcomp9ApplyMassE, h = 13920, w = 13920
          ProductMatrix, h = 72125, w = 72125
            EmbeddedMatrix, h = 72125, w = 13920
              N6ngcomp9ApplyMassE, h = 13920, w = 13920
            SumMatrix, h = 13920, w = 72125
              ScaleMatrix, scale = 2, h = 13920, w = 72125
                N4ngla18EmbeddingTransposeE, h = 13920, w = 72125
              N4ngla18EmbeddingTransposeE, h = 13920, w = 72125
        EmbeddedTransposeMatrix, h = 72125, w = 72125
          EmbeddedMatrix, h = 72125, w = 13920
            N6ngcomp9ApplyMassE, h = 13920, w = 13920