21.4. Verlet time-stepping and Mass-lumping

For wave equations, explizit time-stepping methods are an interesting alternative. But still, one has to solve with the mass matrix \(M\). There are tricks to make the mass matrix cheaply invertible.

21.4.1. The Verlet method:

we start from the second-order wave equation

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

and discretize in space to obtain the ordinary differential equation

\[ M \frac{d^2 u}{dt^2} = -A u. \]

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

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

Initial conditions are discretized by setting

\[ u^0 = u_0 \qquad \text{and} \qquad u^{-1} = u_0 - \tau v_0 \]

21.4.2. Mass-lumping:

Choose a nodal basis, where the interpolation nodes correspond with 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: New higher-order mass-lumped tetrahedral elements for wave propagation modelling: https://arxiv.org/pdf/1803.10065.pdf (2018)

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.

Having a diagonal mass matrix allows a very cheap computation of the next timestep:

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

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(4): 
    mesh.Refine()
    maxh /= 2
mesh = Mesh(mesh)
# Draw (mesh)

order=2   # order=1 or order=2 is supported

tau = maxh / (5*order)
tend = 2
u0 = exp(-50**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())  

# or, in short:
# minv = fes.Mass(rho=1).Inverse()

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

scene = Draw(gfu, order=order, deformation=True, scale=3, euler_angles=[-60,0,-40])
sleep (2)
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()

For higher order version and local time-stepping see: Talk by Marcus Grote

Try examples from there.

21.4.3. Exercises

• Use the Gaussian wave packet as initial condition from before

• Try a problem with different refraction indices (e.g. \(n=1\) for \(x < 0\) and \(n=2\) for \(x \geq 0\)). It enters into the equation as

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

  For constant \(n\), the plane wave function

  \[ u(x,t) = {\mathcal Re} \left\{ e^{i k \cdot x - i \omega t} \right\}, \]

  with wave vector \(k\) and frequency \(\omega\) such that \(|k| = n \omega\) is a soluton to the wave equation on \(\Omega = {\mathbb R}^n\).

  Explain how wave packages are reflected at interfaces.

• Try problems in 3D