68. Mass-lumping and Local time-stepping

Mass-lumping: Choose a nodal basis, where the interpolation nodes correspond to integration points. Thus, the mass-matrix becomes diagonal.

Well know for first order elements.

More recent for higher order methods:

• G. Cohen, P. Joly, J. E. Roberts, and N. Tordman: Higher order triangular finite elements with mass lumping for the wave equation , SINUM 38(6), pp 2047-2078 (2000)

• S. Geevers, W.A. Mulder, and J.J.W. van der Vegt: Ner higher-order mass-lumped tetrahedral elements for wave propagation modelling: https://arxiv.org/pdf/1803.10065.pdf

Basic idea for higher order triangles:

\(P^2\) triangle with the 6 basis functions does not work with mass-lumping, however adding the interior bubble works ! Then we have 7 integration points on the triangle.

from ngsolve import *
from ngsolve.webgui import Draw
from netgen.occ import unit_square
from time import sleep, time

maxh = 0.15
mesh = unit_square.GenerateMesh(maxh=maxh)
for l in range(5): 
    mesh.Refine()
    maxh /= 2
mesh = Mesh(mesh)
# Draw (mesh)

order=2

tau = maxh / (5*order)
tend = 2
u0 = exp(-100**2*( (x-0.5)**2 + (y-0.5)**2))
v0 = 0

fes = H1LumpingFESpace(mesh, order=order)  
u,v = fes.TnT()

mform = u*v*dx(intrules=fes.GetIntegrationRules())
aform = grad(u)*grad(v)*dx

m = BilinearForm(mform, diagonal=True).Assemble()
a = BilinearForm(aform).Assemble()
minv = m.mat.Inverse(fes.FreeDofs())  

68.1. The Verlet method:

the second-order wave equation

\[ \frac{\partial^2 u}{\partial t^2} = \Delta u \]

is discretized with mass matrix \(M\) and stiffness matrix \(A\) as

\[ M \frac{\partial^2 u}{\partial t^2} = -A u. \]

The time-derivative is approximated by a second order finite difference stencil:

\[ M \frac{u^{n+1} - 2 u^n + u^{n-1}}{\tau^2} = -A u^n \]

gfu = GridFunction(fes)
gfu.Set(u0)

scene = Draw(gfu, order=1, deformation=True)
# sleep (3)
unew = gfu.vec.CreateVector()
uold = gfu.vec.CreateVector()
uold.data = gfu.vec

with TaskManager(): 
    for n in range(int(tend/tau)):
        unew.data = 2*gfu.vec - uold 
        unew.data -= tau**2 * minv@a.mat * gfu.vec
        uold.data = gfu.vec
        gfu.vec.data = unew.data
        if n % 100 == 0:
            scene.Redraw()

scene.Redraw()

68.2. Geometry with local details:

Use large time-steps for large elements, and perform smaller time-steps where the mesh is finer.

J. Diaz and M.J. Grote: Energy conserving explicit local time stepping for second-order wave equations, SISC 31(3), pp 1985-2014 (2009)

from netgen.occ import *
rect = MoveTo(-1,-1).Rectangle(2.5,2).Face()
hole = MoveTo(0.5,0.01).Rectangle(0.001,0.8).Face() + \
    MoveTo(0.5,-0.81).Rectangle(0.001,0.8).Face()
local = Circle((0.5,0.8),0.1).Face() + \
    Circle((0.5,0),0.1).Face() + \
    Circle((0.5,-0.8),0.1).Face()
large = rect-local-hole
small = local-hole
large.faces.name="large"
small.faces.name="small"

shape = Glue ([large,small])
geo = OCCGeometry(shape, dim=2)
mesh = Mesh(geo.GenerateMesh(maxh=0.02, grading=0.5))
Draw (mesh);

tau = 0.01
tend = 1
u0 = exp(-10**2*( x**2 + y**2))
v0 = 0
substeps = 20

fes = H1(mesh, order=1)
u,v = fes.TnT()

localdofs = fes.GetDofs(mesh.Materials("small"))
print ("local dofs: ", localdofs.NumSet(),"/",len(localdofs))
Ps = Projector(localdofs, True)   # projection to small
Pl = Projector(localdofs, False)  # projection to large

lumping = IntegrationRule( [(0,0),(1,0),(0,1)], [1/6, 1/6, 1/6])
mform = u*v*dx(intrules = { TRIG: lumping })
aform = grad(u)*grad(v)*dx

mmat = BilinearForm(mform).Assemble().mat
amat = BilinearForm(aform).Assemble().mat
minva = mmat.CreateSmoother()@amat
# APs = minva@Ps    # operator composition
APs = mmat.CreateSmoother()@(amat@Ps.CreateSparseMatrix()).DeleteZeroElements(1e-12)
APl = minva@Pl

local dofs:  563 / 14351

gfu = GridFunction(fes)
gfu.Set(u0)

scene = Draw(gfu, order=1, deformation=True)
sleep (3)
unew = gfu.vec.CreateVector()
uold = gfu.vec.CreateVector()
z = gfu.vec.CreateVector()
znew = gfu.vec.CreateVector()
zold = gfu.vec.CreateVector()
w = gfu.vec.CreateVector()
uold.data = gfu.vec

with TaskManager(): # pajetrace=10**8):
    for n in range(int(tend/tau)):
        w.data = APl * gfu.vec
        zold.data = gfu.vec
        z.data = zold - (tau/substeps)**2/2*(w+APs*zold)
        for m in range(1, substeps):
            znew.data = 2*z-zold
            znew.data -= (tau/substeps)**2*(w+APs*z)
            zold,z,znew = z,znew,zold
        unew.data = 2*z-uold
        uold.data = gfu.vec
        gfu.vec.data = unew.data
        if n % 10 == 0:
            scene.Redraw()

print ("nze", amat.nze)
asmall =  (amat@Ps.CreateSparseMatrix()).DeleteZeroElements(1e-12)
print ("nzes", asmall.nze)

nze 99169
nzes 3777