Firedrake test

Helmholtz

We start by implementing the helmholtz equation on a 2d domain in order to learn how to express PDEs in the Firedrake methodology.

Let Ω\Omega be the unit square and use Γ\Gamma to denote the boundary.

The helmholtz equation with Neumann boundary conditions is ∇2u=u−f∇u⋅n=0. \begin{align*} \nabla^2 u &= u - f \\ \nabla u \cdot \mathbf{n} &= 0. \end{align*}

Let us use the galerkin formalism over a function space VV. Then let v∈Vv \in V be a test function.

We start by remembering how the variational problem is derived. Start with ∫Ω(u−∇⋅(∇u))v dx=∫Ωvf dx \int_\Omega (u - \nabla \cdot( \nabla u))v \intd{\mathbf{x}} = \int_\Omega vf \intd{\mathbf{x}} . Integration by parts in cartesian coordinates gives ∫Ω(∇⋅∇u)v dx=∫Ω∇v⋅∇u−∇⋅(v∇u) dx=div thm∫Ω∇v⋅∇u dx−∫Γv∇u⋅n dΓ \begin{align*} \int_\Omega (\nabla \cdot \nabla u) v \intd{\mathbf{x}} &= \int_\Omega \nabla v \cdot \nabla u - \nabla \cdot (v\nabla u) \intd{\mathbf{x}} \\ &\stackrel{\textrm{div thm}}{=}\int_\Omega \nabla v \cdot \nabla u \intd{\mathbf{x}} - \int_\Gamma v \nabla u \cdot \mathbf{n} \intd{\Gamma} \end{align*} and so we get ∫Ω∇u⋅∇v+uv dx=∫Ωvf dx+∫Γv∇u⋅n dΓ=∫Ωvf dx+0. \begin{align*} \int_\Omega \nabla u \cdot \nabla v + uv \intd{\mathbf{x}} &= \int_{\Omega} vf \intd{\mathrm{x}} + \int_\Gamma v \nabla u \cdot \mathbf{n} \intd{\Gamma} \\ &= \int_{\Omega} vf \intd{\mathbf{x}} + 0. \end{align*}

we can use the method of constructed solutions and set the ansatz u=cos⁡(2πx)cos⁡(2πy) u= \cos(2\pi x) \cos(2 \pi y) such that the boundary condition is satisified and f=(1.0+8.0π2)cos⁡(2πx)cos⁡(2πy)f = (1.0 + 8.0\pi^2)\cos(2\pi x)\cos(2\pi y).

Summary of how to use firedrake for linear variational problems

  • Construct mesh
  • Construct function space
  • Get trial/test functions
  • Construct governing equations
  • Projection is done via the interpolate method on a fd.Function object.
  • Pass to fd.solve which calls petsc

Note: fd.dx is a placeholder for integration over Rn\mathbb{R}^n. For variational problems Neumann boundary conditions are incorporated into the variational formulation. Determining the appropriate method of integrating over the boundary weakly must be determined.

Code:

import firedrake as fd
ne_x, ne_y = (10, 10)
mesh = fd.UnitSquareMesh(ne_x, ne_y)

cont_galerkin_deg = 1
cg_func_space = fd.FunctionSpace(mesh, "CG", cont_galerkin_deg)

trial_funcs = fd.TrialFunction(cg_func_space)
test_funcs = fd.TestFunction(cg_func_space) # there's jargon for when trial func space = test func space.
mesh_x, mesh_y = fd.SpatialCoordinate(mesh)
f = fd.Function(cg_func_space)
f.interpolate((1.0 + 8.0 * fd.pi**2) * fd.cos(mesh_x * fd.pi * 2) * fd.cos(mesh_y * fd.pi * 2))
a = (fd.inner(fd.grad(trial_funcs), fd.grad(test_funcs)) + fd.inner(trial_funcs,test_funcs)) * fd.dx
L = fd.inner(f, trial_funcs) * fd.dx

u = fd.Function(cg_func_space)

fd.solve(a==L, u, solver_parameters={'ksp_type': 'cg', 'pc_type': 'none'})

fd.File("/content/firedrake/helmholtz.pvd").write(u)

Parent post: