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{align}\begin{aligned}\frac{\partial^2\phi}{\partial t^2} - \nabla^2 \phi = 0\\\nabla \phi \cdot n = 0 \ \textrm{on}\ \Gamma_N\\\phi = \frac{1}{10\pi}\cos(10\pi t) \ \textrm{on}\ \Gamma_D\end{aligned}\end{align} \]

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

\[ \begin{align}\begin{aligned}\frac{\partial\phi}{\partial t} = - p\\\frac{\partial p}{\partial t} + \nabla^2 \phi = 0\\\nabla \phi \cdot n = 0 \ \textrm{on}\ \Gamma_N\\p = \sin(10\pi t) \ \textrm{on}\ \Gamma_D\end{aligned}\end{align} \]

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

\[\int_\Omega \frac{\partial p}{\partial t} v\,\mathrm d x = \int_\Omega \nabla\phi\cdot\nabla v\,\mathrm d x \quad \forall v \in V\]

For a suitable function space V. Note that the absence of spatial derivatives in the equation for \(\phi\) 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 \(\phi\) 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 \(\phi\) 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 \(\phi\) 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.