# 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 = File("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:


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.