5.2. DG Operators in NGSolve

We want to efficiently implement the mass and gradient bilinear forms from the previous section.

Spaces

We use L2 and VectorL2 spaces

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

geo = OCCGeometry(unit_square.shape, dim = 2)
mesh = Mesh(geo.GenerateMesh(maxh = 0.01))


order = 6
fes_p = L2(mesh, order = order)
fes_v = VectorL2(mesh, order = order-1)

p,p_ = fes_p.TnT()
v,v_ = fes_v.TnT()

print(fes_p.ndof, fes_v.ndof)

649656 974484

gfp = GridFunction(fes_p)
gfp.vec[100]=1.
Draw(gfp)

BaseWebGuiScene

Mass matrices

Mass Operators can be assembled

from time import time

now = time()
mass_p_inv = BilinearForm(p*p_*dx).Assemble().mat.Inverse()
mass_v_inv = BilinearForm(v*v_*dx).Assemble().mat.Inverse()
print("Assembling and factorization: ",time()-now)

tmp_p = mass_p_inv.CreateVector()
tmp_v = mass_v_inv.CreateVector()

tmp_p.SetRandom()
tmp_v.SetRandom()
now = time()
for i in range(20):
    tmp_p.data = mass_p_inv*tmp_p
    tmp_v.data = mass_v_inv*tmp_v
print("20 matrix applications: ",time()-now)

Assembling and factorization:  7.4811577796936035

20 matrix applications:  10.61510157585144

alternatively they can be applied matrix-free

from time import time
now = time()

mass_p_inv = fes_p.Mass(1).Inverse()
mass_v_inv = fes_v.Mass(1).Inverse()
print("Assembling and factorization: ",time()-now)
tmp_p = mass_p_inv.CreateVector()
tmp_v = mass_v_inv.CreateVector()
tmp_p.SetRandom()
tmp_v.SetRandom()
now = time()
for i in range(20):
    tmp_p.data = mass_p_inv*tmp_p
    tmp_v.data = mass_v_inv*tmp_v
    
print("20 matrix applications: ",time()-now)

Assembling and factorization:  0.05927157402038574

20 matrix applications:  0.5783483982086182

Gradients

Version with .Other()

fes = fes_p*fes_v

gf = GridFunction(fes)
gf.vec.SetRandom()


(P,V),(P_,V_) = fes.TnT()


n = specialcf.normal(mesh.dim)


now = time()
B = BilinearForm(fes, nonassemble=True)
B += grad(P)*V_ * dx + 0.5*(P.Other()-P)*(V_*n) * dx(element_boundary=True)
BT = BilinearForm(fes, nonassemble=True)
BT += grad(P_)*V * dx + 0.5*(P_.Other()-P_)*(V*n) * dx(element_boundary=True)
print("Assembling and factorization: ",time()-now)


now = time()
for i in range(20):
    gf.vec.data = B.mat*gf.vec
    
print("20 matrix applications: ",time()-now)

Assembling and factorization:  0.0006272792816162109

20 matrix applications:  3.8412270545959473

Version with TraceOperator

fes_tr = FacetFESpace(mesh,order=order)
pT = fes_tr.TrialFunction()


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



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

Btr = BilinearForm(trialspace=fes_tr, testspace=fes_v)
Btr += pT * (v_*n) * dx(element_boundary=True)
Btr.Assemble()

emb_p = fes.embeddings[0]
emb_v = fes.embeddings[1]

B = emb_v @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T
print("Assembling and factorization: ",time()-now)

now = time()
for i in range(20):
    gf.vec.data = B*gf.vec
    
print("20 matrix applications: ",time()-now)

Assembling and factorization:  1.4064881801605225

20 matrix applications:  1.1851904392242432

Geometry free TraceOperator

fes_v = VectorL2(mesh, order = order-1, piola = True)

v,v_ = fes_v.TnT()
now = time()
Bel = BilinearForm(trialspace=fes_p, testspace=fes_v, geom_free=True)
Bel += grad(p)*v_ * dx -p*(v_*n) * dx(element_boundary=True)
Bel.Assemble()

Btr = BilinearForm(trialspace=fes_tr, testspace=fes_v, geom_free=True)
Btr += 0.5 * pT * (v_*n) * dx(element_boundary=True)
Btr.Assemble()

emb_p = fes.embeddings[0]
emb_v = fes.embeddings[1]

B = emb_v @ (Bel.mat + Btr.mat @ traceop) @ emb_p.T
print("Assembling and factorization: ",time()-now)

now = time()
for i in range(20):
    gf.vec.data = B*gf.vec
    
print("20 matrix applications: ",time()-now)

Assembling and factorization:  0.043317556381225586

20 matrix applications:  0.2828328609466553