Double slit experiment

Here we solve a linear wave equation using an explicit timestepping scheme. This example demonstrates the use of an externally generated mesh, pointwise operations on Functions, and a time varying boundary condition. The strong form of the equation we set out to solve is:

$$\begin{array}{r}\begin{array}{c}\frac{{\partial }^2\varphi }{\partial t^2}-{\mathrm{\nabla }}^2\varphi =0\\ \mathrm{\nabla }\varphi \cdot n=0\text{on}{\mathrm{\Gamma }}_N\\ \varphi =\frac{1}{10\pi }\cos (10\pi t)\text{on}{\mathrm{\Gamma }}_D\end{array}\end{array}$$

To facilitate our choice of time integrator, we make the substitution:

$$\begin{array}{r}\begin{array}{c}\frac{\partial \varphi }{\partial t}=-p\\ \frac{\partial p}{\partial t}+{\mathrm{\nabla }}^2\varphi =0\\ \mathrm{\nabla }\varphi \cdot n=0\text{on}{\mathrm{\Gamma }}_N\\ p=\sin (10\pi t)\text{on}{\mathrm{\Gamma }}_D\end{array}\end{array}$$

We then form the weak form of the equation for $p$. Find $p\in V$ such that:

$${\int }_{\mathrm{\Omega }}\frac{\partial p}{\partial t}v\mathrm{d}x={\int }_{\mathrm{\Omega }}\mathrm{\nabla }\varphi \cdot \mathrm{\nabla }v\mathrm{d}x\mathrm{\forall }v\in V$$

For a suitable function space V. Note that the absence of spatial derivatives in the equation for $\varphi$ makes the weak form of this equation equivalent to the strong form so we will solve it pointwise.

In time we use a simple symplectic method in which we offset $p$ and $\varphi$ by a half timestep.

This time we created the mesh with Gmsh:

gmsh -2 wave_tank.geo

We can then start our Python script and load this mesh:

from firedrake import *
mesh = Mesh("wave_tank.msh")

We choose a degree 1 continuous function space, and set up the function space and functions. Setting the name parameter when constructing Function objects will set the name used in the output file:

V = FunctionSpace(mesh, 'Lagrange', 1)
p = Function(V, name="p")
phi = Function(V, name="phi")

u = TrialFunction(V)
v = TestFunction(V)

Output the initial conditions:

outfile = VTKFile("out.pvd")
outfile.write(phi)

We next establish a boundary condition object. Since we have time-dependent boundary conditions, we first create a Constant to hold the value and use that:

bcval = Constant(0.0)
bc = DirichletBC(V, bcval, 1)

Now we set the timestepping variables:

T = 10.
dt = 0.001
t = 0
step = 0

Finally we set a flag indicating whether we wish to perform mass-lumping in the timestepping scheme:

lump_mass = True

Now we are ready to start the timestepping loop:

while t <= T:
    step += 1

Update the boundary condition value for this timestep:

bcval.assign(sin(2*pi*5*t))

Step forward $\varphi$ by half a timestep. Since this does not involve a matrix inversion, this is implemented as a pointwise operation:

phi -= dt / 2 * p

Now step forward $p$. This is an explicit timestepping scheme which only requires the inversion of a mass matrix. We have two options at this point, we may either lump the mass, which reduces the inversion to a pointwise division:

if lump_mass:
    p.dat.data[:] += assemble(dt * inner(nabla_grad(v), nabla_grad(phi))*dx).dat.data_ro / assemble(v*dx).dat.data_ro

In the mass lumped case, we must now ensure that the resulting solution for $p$ satisfies the boundary conditions:

bc.apply(p)

Alternatively, we can invert the mass matrix using a linear solver:

else:
    solve(u * v * dx == v * p * dx + dt * inner(grad(v), grad(phi)) * dx,
          p, bcs=bc, solver_parameters={'ksp_type': 'cg',
                                        'pc_type': 'sor',
                                        'pc_sor_symmetric': True})

Step forward $\varphi$ by the second half timestep:

phi -= dt / 2 * p

Advance time and output as appropriate, note how we pass the current timestep value into the write() method, so that when visualising the results Paraview will use it:

t += dt
if step % 10 == 0:
    outfile.write(phi, time=t)

The following animation, produced in Paraview, illustrates the output of this simulation:

A python script version of this demo can be found here. The gmsh input file is here.