# 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: # # .. math:: # # \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 # # To facilitate our choice of time integrator, we make the substitution: # # .. math:: # # \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 # # We then form the weak form of the equation for :math:p. Find # :math:p \in V such that: # # .. math:: # # \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 :math:\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 :math:p # and :math:\phi by a half timestep. # # This time we created the mesh with Gmsh _: # # .. code-block:: bash # # 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 :class:.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 :class:.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 :math:\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 :math: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 += assemble(dt * inner(nabla_grad(v), nabla_grad(phi))*dx) / assemble(v*dx) # In the mass lumped case, we must now ensure that the resulting # solution for :math: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 :math:\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 :meth:~.File.write method, so that when # visualising the results Paraview will use it:: t += dt if step % 10 == 0: outfile.write(phi, time=t) # .. only:: html # # The following animation, produced in Paraview, illustrates the output of this simulation: # # .. container:: youtube # # .. youtube:: xhxvM1N8mDQ?modestbranding=1;controls=0;rel=0 # :width: 600px # # .. only:: latex # # An animation, produced in Paraview, illustrating the output of this simulation can be found on youtube _. # # # A python script version of this demo can be found here __. The gmsh input file is here __.