```
# This magic makes plots appear in the browser
%matplotlib inline
import matplotlib.pyplot as plt
# Load Firedrake on Colab
try:
import firedrake
except ImportError:
!wget "https://github.com/g-adopt/tutorials/releases/latest/download/firedrake-install-real.sh" -O "/tmp/firedrake-install.sh" && bash "/tmp/firedrake-install.sh"
import firedrake
```

# The positive-definite Helmholtz equation¶

We start by considering the symmetric positive definite "Helmholtz" problem on a unit square domain $\Omega$ with boundary $\Gamma$. We seek to find the solution $u \in V$, where $V$ is some finite element space, satisfying the boundary value problem:

$$ -\nabla^2 u + u = f \text{ on } \Omega = [0, 1] \times [0, 1], \\ \nabla u \cdot \vec{n} = 0 \text{ on } \Gamma, $$

where $\vec{n}$ is the outward-pointing unit normal to $\Gamma$ and $f \in V$ is a known source function. Note that the Laplacian and mass terms have opposite signs, which makes this equation much more benign than the symmetric, indefinite operator $(\nabla^2 + I)$.

Now, we write the problem in its variational form by multiplying by a test function $v \in V$ and integrating by parts over the domain $\Omega$:

$$ \int_\Omega \nabla u\cdot\nabla v + uv\ \mathrm{d}x = \int_\Omega vf\ \mathrm{d}x + \underbrace{\int_\Gamma v \nabla u \cdot \vec{n}\, \mathrm{d} s}_{=0}. $$

Note that the boundary condition has been enforced weakly by removing the surface integral term resulting from the integration by parts.

We will now proceed to solve this problem using an $H^1$ conforming method. The full finite element problem reads as follows: find $u \in V$ such that:

$$ \int_\Omega \nabla u\cdot\nabla v + uv\ \mathrm{d}x = \int_\Omega vf\ \mathrm{d}x, \text{ for all } v \in V. $$

We can choose the source function $f$, so lets take it to be:

$$ f(x, y) = (1.0 + 8.0\pi^2)\cos(2\pi x)\cos(2\pi y). $$

This conveniently yields the analytical solution:

$$ u(x, y) = \cos(2\pi x)\cos(2\pi y). $$

However, we wish to employ this as an example for the finite element method, so lets go ahead and produce a numerical solution. First, we need our domain $\Omega$. Firedrake can read meshes from a number of popular mesh generators. In addition, it provides utility functions to create many "standard" meshes in one, two, and three dimensions. For a complete list of provided meshes, we can peruse the online Firedrake documentation. Unsurprisingly, amongst them is a `UnitSquareMesh`

.

To start, we must make Firedrake available in the notebook. We will import it through the `gadopt`

package. To save on typing, we will import all of the public API into the current namespace. For plotting within this notebook, we also explicitly import the `triplot`

function.

```
from gadopt import *
from firedrake.pyplot import triplot, tripcolor
```

Now we want to define our mesh. We already know that `UnitSquareMesh`

will give us an appropriate mesh of the domain $\Omega$, but what are the arguments to the constructor? The user-facing API is usually well-documented, and we can access this information via Python using the builtin `help`

command:

```
help(UnitSquareMesh)
```

Help on cython_function_or_method in module firedrake.utility_meshes: UnitSquareMesh(nx, ny, reorder=None, diagonal='left', quadrilateral=False, distribution_parameters=None, comm=<mpi4py.MPI.Intracomm object at 0x7266f1777e30>, name='firedrake_default', distribution_name=None, permutation_name=None) Generate a unit square mesh :arg nx: The number of cells in the x direction :arg ny: The number of cells in the y direction :kwarg quadrilateral: (optional), creates quadrilateral mesh. :kwarg reorder: (optional), should the mesh be reordered :kwarg distribution_parameters: options controlling mesh distribution, see :func:`.Mesh` for details. :kwarg comm: Optional communicator to build the mesh on. :kwarg name: Optional name of the mesh. :kwarg distribution_name: the name of parallel distribution used when checkpointing; if `None`, the name is automatically generated. :kwarg permutation_name: the name of entity permutation (reordering) used when checkpointing; if `None`, the name is automatically generated. The boundary edges in this mesh are numbered as follows: * 1: plane x == 0 * 2: plane x == 1 * 3: plane y == 0 * 4: plane y == 1

```
# In the notebook, we can also use a special "?" magic to pop up a help box
?UnitSquareMesh
```

We'll get to what some of the other arguments mean later, but right now we see that the first two allow us to specify the number of cells in the x-, and y-directions respectively.

```
mesh = UnitSquareMesh(10, 10)
```

What does this mesh look like? Firedrake has some built in plotting of meshes and functions, using matplotlib (that's why we had the magic matplotlib line at the top of this notebook). The built-in plotting functions correspond with those from matplotlib, like tripcolor, tricontour, tricontourf. For this plot, we'll also add a legend that will show us what the numeric ID is of each boundary segment. Visualizing the mesh with a legend is helpful when you have to apply boundary conditions to part of a mesh and forget how the segments are numbered.

```
# We test the output of the notebooks with out continuous integration system to
# make sure they run.
# Unfortunately we can't compare the plots from run to run, so the above line
# indicates that the output of this cell should be ignored for the purposes
# of testing.
fig, axes = plt.subplots()
triplot(mesh, axes=axes)
axes.legend();
```

Having selected a discretisation of our domain $\Omega$, we need to decide on the finite-dimensional function space $V_h \subset V$ in which we’d like to solve the problem. Since we are using an $H^1$ conforming method, the space of continuous piecewise defined polynomials of a fixed degree will be appropriate. As an example, let us use the space of piecewise linear polynomials that are continuous between elements:

```
V = FunctionSpace(mesh, "Lagrange", 1)
```

We now define the problem. We'll create the *symbolic* objects that correspond to test and trial spaces, and the linear and bilinear forms in our problem.

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

For the right hand side forcing, we'll use a UFL expression, incorporating information about the $x$ and $y$ coordinates. We make a symbolic representation of the coordinates in our mesh (these will be evaluated when we actually do the calculation).

```
x, y = SpatialCoordinate(mesh)
f = (1 + 8*pi*pi)*cos(2*pi*x)*cos(2*pi*y)
```

We can now define the bilinear and linear forms for the left and right hand sides of our equation respectively in UFL:

```
a = (dot(grad(v), grad(u)) + v * u) * dx
L = f * v * dx
```

Finally we are now ready to solve the equation. We define $u_h$ to be a function holding the solution:

```
uh = Function(V)
```

Since we know that the Helmholtz equation defines a symmetric problem, we instruct PETSc to employ the conjugate gradient method. We do not consider preconditioning, for now.

```
solve(a == L, uh, solver_parameters={'ksp_type': 'cg', 'pc_type': 'none'})
```

Let's have a look at the solution. You can use all the same arguments as for the matplotlib tripcolor function for things like changing the colormap, the min and max values, and so forth.