5. Local Implicit Time Stepping For The Wave Equation

Rafael Dorigo

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*}\]

form a simple DG version with central fluxes. The discrete gradient \(B\) is defined by the bilinear-form

(5.1)\[\begin{align} b(p,v) = \sum_{T} \Big\{ \int_T \nabla p \, v + \int_{\partial T} (\{ p \} - p) \, v_n \, ds \Big\} \end{align}\]

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import *
import numpy as np
from tqdm import tqdm

Simple trumpet with a two small holes to create thin elements @shirnschall

airBoxW = 12
airBoxH = 12
outletAngleDeg = 60
tubeW = 4
tubeH = 1
holeW = 0.001
holeH = 0.05
def create_trumpet():    
    rect = MoveTo(0,-airBoxH/2).Rectangle(airBoxW,airBoxH).Face()
    tube = MoveTo(-tubeW,-tubeH/2).Rectangle(tubeW,tubeH).Face()
    outlet = Face(MoveTo(0,-tubeH/2).Rotate(-outletAngleDeg/2).Line(3).Rotate(-90).Line(.25).Rotate(-90).Line(3).Close().Wire()).Reversed()
    hole = MoveTo(-tubeW/2, tubeH/2 - holeH).Rectangle(holeW, holeH).Face()
    hole2 = MoveTo(-tubeW/2 - 2 * holeW, tubeH/2 - holeH).Rectangle(holeW, holeH).Face()

    geo = rect + tube - outlet - outlet.Mirror(Axis((0,0,0),X)) - hole - hole2
    return geo

dim = 2
h = 0.2

shape = create_trumpet()
shape = OCCGeometry(shape, dim = dim)

mesh = Mesh(shape.GenerateMesh(maxh = h))
Draw(mesh)
print(f"Number of elements: {mesh.ne}")

Number of elements: 8691

order = 2
fes_p = L2(mesh, order=order+1, dgjumps=True) 
fes_u = VectorL2(mesh, order=order, piola=True, dgjumps=True)
fes_tr = FacetFESpace(mesh, order=order+1)
fes = fes_p * fes_u

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

gfu = GridFunction(fes_u)
gfp = GridFunction(fes_p)
gfp.Set( exp(-10*((x+tubeW)**2+y**2+z**2)))

p, u = fes.TrialFunction()
q, v = fes.TestFunction()
phat = fes_tr.TrialFunction()

n = specialcf.normal(mesh.dim)
dS = dx(skeleton=True)   
def jump(p): return p.Other()-p
def avgn(v): return 0.5*(v*n-v.Other()*n.Other())

embp, embu = fes.embeddings

dt = 0.2 * h / (order+1)**2
print ("dt = ", dt)

dt =  0.004444444444444445

print(f"fes ndofs: {fes.ndof}")
print(f"fes_u ndofs: {fes_u.ndof}")
print(f"fes_p ndofs: {fes_p.ndof}")

fes ndofs: 191202
fes_u ndofs: 104292
fes_p ndofs: 86910

Draw (gfp, order=3)

BaseWebGuiScene

5.1. Mesh splitting

F = specialcf.JacobianMatrix(mesh.dim)
Finv = Inv(F)
detF = Det(F)
Norm_Finv = Norm(Finv)
el_norms = Integrate(Norm_Finv*1/detF, mesh, element_wise=True)
el_norms_numpy = np.array(el_norms)

sorted_el_norms = -np.sort(-el_norms_numpy) #sort descending
ref_norm = el_norms_numpy.mean()# using elements > mean 
impl_els = np.where(el_norms_numpy > ref_norm, 1, 0)
    
print(f'min = {min(el_norms_numpy)}')
print(f'max = {max(el_norms_numpy)}')
print(f'mean = {el_norms_numpy.mean()}')
print(f'median = {np.median(el_norms_numpy)}')
print(f'ref_norm = {ref_norm}')
print(f'implicit / total elements: {np.count_nonzero(impl_els)}/{len(impl_els)}')

min = 2.8500061904762743
max = 832.1624920872001
mean = 14.553329401327039
median = 4.082563397670109
ref_norm = 14.553329401327039
implicit / total elements: 389/8691

ba_implicit_els = BitArray(mesh.ne)
ba_explicit_els = BitArray(mesh.ne)
ba_interface_edges = BitArray(mesh.nedge)
ba_explicit_edges = BitArray(mesh.nedge)
ba_implicit_edges = BitArray(mesh.nedge)

ba_implicit_els[:] = 0
ba_explicit_els[:] = 0
ba_interface_edges[:] = 0
ba_explicit_edges[:] = 0
ba_implicit_edges[:] = 0

for el in mesh.Elements():
    if impl_els[el.nr] == 1:
        ba_implicit_els[el.nr] = 1
        for e in el.edges:
            ba_implicit_edges[e.nr] = 1
    else:
        ba_explicit_els[el.nr] = 1
        for e in el.edges:
            ba_explicit_edges[e.nr] = 1
        
ba_interface_edges = ba_explicit_edges & ba_implicit_edges     

ba_local_implicit_dofs = BitArray(fes.ndof)
ba_local_implicit_dofs[:] = 0

for el in mesh.Elements():
    if ba_implicit_els[el.nr] == 1:
        for nr in fes.GetDofNrs(el):
            ba_local_implicit_dofs[nr] = 1
    if ba_explicit_els[el.nr] == 1:
        for e in el.edges:
            if ba_interface_edges[e.nr] == 1:
                for nr in fes.GetDofNrs(el):
                    ba_local_implicit_dofs[nr] = 1

ba_local_implicit_dofs_u = BitArray(fes_u.ndof)
ba_local_implicit_dofs_u[:] = 0

for el in mesh.Elements():
    if ba_implicit_els[el.nr] == 1:
        for nr in fes_u.GetDofNrs(el):
            ba_local_implicit_dofs_u[nr] = 1
    if ba_explicit_els[el.nr] == 1:
        for e in el.edges:
            if ba_interface_edges[e.nr] == 1:
                for nr in fes_u.GetDofNrs(el):
                    ba_local_implicit_dofs_u[nr] = 1

5.2. Time stepping combination of Verlet and Crank-Nicolson

(5.2)\[\begin{align} p_{h}^{n+\frac{1}{2}} - p_{h}^{n} &= - \frac{\tau}{2} M_p^{-1} B^{\top} u_{h}^{n} \\ \left(I + \frac{\tau^2}{4} M_u^{-1} B_i M_p^{-1} B_i^{\top} \right) u_{h}^{n+1} &= u_{h}^{n} + \tau M_u^{-1} B_e p_{h}^{n+\frac{1}{2}} + \frac{\tau}{2} M_u^{-1} B_i (p_{h}^{n+\frac{1}{2}} + p_{h}^{n}) \\ p_{h}^{n+1}- p_{h}^{n+\frac{1}{2}} &= - \frac{\tau}{2} M_p^{-1} B^{\top} u_{h}^{n+1} \end{align}\]

A = BilinearForm(fes)
A += -p*q*dx +u*v*dx 
A += dt/2*(grad(p)*v + grad(q)*u)*dx
A += dt/2*(jump(p)*avgn(v)+jump(q)*avgn(u)) * dS
A.Assemble()

p1, q1 = fes_p.TnT()
u1, v1 = fes_u.TnT()
phat1, qhat1 = fes_tr.TnT()
Bel = BilinearForm(trialspace=fes_p, testspace=fes_u, geom_free = True)
Bel += grad(p1)*v1 * dx -p1*(v1*n) * dx(element_boundary=True)
Bel.Assemble()
Btr = BilinearForm(trialspace=fes_tr, testspace=fes_u, geom_free = True)
Btr += phat1 * (v1*n) *dx(element_boundary=True)
Btr.Assemble();

B = Bel.mat + Btr.mat @ traceop
B_T = B.T

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

Ps = Projector(ba_local_implicit_dofs_u, True)   # projection to small
Pl = Projector(ba_local_implicit_dofs_u, False)  # projection to large
B_e = Pl @ B
B_i = Ps @ B

print (f"local implicit dofs: {ba_local_implicit_dofs.NumSet()} / {len(ba_local_implicit_dofs)}")
print ("local implicit dofs of fes u: ", ba_local_implicit_dofs_u.NumSet(),"/",len(ba_local_implicit_dofs_u))

local implicit dofs: 8756 / 191202
local implicit dofs of fes u:  4776 / 104292

invA = A.mat.Inverse(freedofs=ba_local_implicit_dofs, inverse="sparsecholesky")

# delete non zeros elements of matrix to speed up matrix multiplication
Anze = A.mat.DeleteZeroElements(10e-12)
invAnze = Anze.Inverse(freedofs=ba_local_implicit_dofs, inverse="sparsecholesky")
invmstar = embu.T @ invAnze @ embu
mstarloc = massu + dt*dt/4*B_i @ invmassp @ B_T

invmassuB = invmassu @ B
invmasspB_T = invmassp @ B_T
invmstar_mstar = invmstar @ mstarloc

# Visulaize implicit dofs
ba_gfu = BitArray(fes.ndof)
ba_gfu[:] = 0
gfuim = GridFunction(fes)
gfuim.vec[:] = 0
for i in range(len(ba_local_implicit_dofs)):
    if ba_local_implicit_dofs[i] == 1:
        gfuim.vec.data[i] = 1
        
Draw(gfuim.components[0], mesh)

BaseWebGuiScene

gfp.Set( exp(-10*((x+tubeW)**2+y**2+z**2)))
gfu.vec[:] = 0
gfp_halfstep = gfp.vec.CreateVector()
gfuold = gfu.vec.CreateVector()
res = gfu.vec.CreateVector()

scene = Draw (gfp, order=3, deformation=True)
    
t = 0
tend = 10
cnt = 0
loop = [i*dt for i in range(int(tend/dt))]
with TaskManager():
    for i in loop:
        gfp_halfstep.data = gfp.vec - dt/2 * invmassp @ B.T * gfu.vec
        
        res.data = dt * B_e * gfp_halfstep + dt/2 * B_i * (gfp_halfstep + gfp.vec) + massu * gfu.vec
        gfuold.data = invmassu * res 
        gfu.vec.data = gfuold
        gfu.vec.data += invmstar * (res - mstarloc * gfuold)
        gfp.vec.data = gfp_halfstep - dt/2 * invmassp @ B.T * gfu.vec
        cnt = cnt+1
        if cnt%20 == 0:
            scene.Redraw()