`firedrake-adjoint`

¶*This example is modified from the equivalent dolfin-adjoint demo*

In this example, we will look at how to use `firedrake-adjoint`

to optimise for strong (Dirichlet) conditions in a steady problem. `firedrake-adjoint`

is a thin compatibility layer for the `dolfin-adjoint`

package, a python package to **automatically derive the discrete adjoint and tangent linear models** of forward problems written using Firedrake.

For the minimisation, we will need `scipy`

, which is installed with

```
pip install scipy
```

In the activated firedrake virtualenv.

As usual, we begin with some notebook magic (so that plots appear nicely) and importing Firedrake.

In [1]:

```
%matplotlib notebook
import matplotlib.pyplot as plt
from firedrake import *
```

`firedrake_adjoint`

package, which *overrides* much of the Firedrake interface so that an *annotated tape* for adjoint calculations can be built automatically.

In [2]:

```
from firedrake_adjoint import *
```

Now we will set up the problem. We consider minimising the compliance:

$$ \min_{g, u, p} \ \frac{1}{2}\int_{\Omega} \nabla u \cdot \nabla u\,\text{d}x + \frac{\alpha}{2} \int_{\Gamma_{\textrm{in}}} g^2\,\text{d}s $$subject to the Stokes equations $$ \begin{split}-\nu \Delta u + \nabla p &= 0 \quad \text{in $\Omega$}, \\ \nabla \cdot u &= 0 \quad \text{in $\Omega$}, \end{split} $$

and Dirichlet conditions

$$ \begin{split} u &= g \quad \text{on $\Gamma_\text{circ}$}, \\ u &= f \quad \text{on $\Gamma_\text{in}$}, \\ u &= 0 \quad \text{on $\Gamma_\text{top} \cup \Gamma_\text{bottom}$}, \\ p &= 0 \quad \text{on $\Gamma_\text{out}$}. \\ \end{split} $$Here, $u$ and $p$ are unknown velocity and pressure, $f$ is a prescribed inflow, $g$ is the control variable that we will optimise for and $\alpha$ is a regularisation parameter. This corresponds physically to minimising the loss of energy as heat by controlling the in/outflow on $\Gamma_\text{circ}$. The regularisation parameter penalises too many non-zero control values.

`Mesh`

constructor, passing the filename of the mesh in question.

In [3]:

```
mesh = Mesh("stokes-control.msh")
```

In [4]:

```
# NBVAL_IGNORE_OUTPUT
fig, axes = plt.subplots()
triplot(mesh, axes=axes)
axes.axis("off")
axes.set_aspect("equal")
axes.legend(loc="upper right");
```

The forward problem should be familiar by now, we create a mixed function space for $u$ and $p$ and set up trial and test functions. We specify a parabolic velocity at the inflow, and no-slip (zero velocity) conditions on the side walls. The zero-pressure outflow condition is enforced weakly.

Our variational formulation for Stokes reads as follows. Find $(u, p) \in V\times Q$ such that:

$$ \begin{align} \int_\Omega \nu \nabla u : \nabla v\,\text{d}x - \int_\Omega p \nabla \cdot v\,\text{d}x - \int_{\Gamma_{\text{circ}}} \nu (\nabla u \cdot n) \cdot v\,\text{d}s &= 0 \quad \forall v \in V,\\ - \int_\Omega q \nabla \cdot u\,\text{d}x &= 0 \quad \forall q \in Q.\\ u &= 0 \quad \text{on $\Gamma_{\text{top}} \cup \Gamma_{\text{bottom}}$},\\ u &= \left[\frac{y(10 - y)}{25}, 0\right]^T \quad \text{on $\Gamma_{\text{in}}$},\\ u &= g \quad \text{on $\Gamma_{\text{circ}}$},\\ \frac{\text{d}p}{\text{d}n} &= 0 \quad \text{on $\Gamma_{\text{out}}$}. \end{align} $$In [5]:

```
V = VectorFunctionSpace(mesh, "Lagrange", 2)
Q = FunctionSpace(mesh, "Lagrange", 1)
W = V*Q
v, q = TestFunctions(W)
u, p = TrialFunctions(W)
nu = Constant(1) # Viscosity coefficient
x, y = SpatialCoordinate(mesh)
u_inflow = as_vector([y*(10-y)/25.0, 0])
noslip = DirichletBC(W.sub(0), (0, 0), (3, 5))
inflow = DirichletBC(W.sub(0), interpolate(u_inflow, V), 1)
static_bcs = [inflow, noslip]
```

`Function`

which will hold the boundary values, we then build a `DirichletBC`

object as normal.

In [6]:

```
g = Function(V, name="Control")
controlled_bcs = [DirichletBC(W.sub(0), g, 4)]
bcs = static_bcs + controlled_bcs
```

Now we define the bilinear and linear forms.

In [7]:

```
a = nu*inner(grad(u), grad(v))*dx - inner(p, div(v))*dx - inner(q, div(u))*dx
L = Constant(0)*q*dx
```

`firedrake-adjoint`

annotates it. We'll also take a look at the solution.

In [8]:

```
w = Function(W)
solve(a == L, w, bcs=bcs, solver_parameters={"mat_type": "aij",
"ksp_type": "preonly",
"pc_type": "lu",
"pc_factor_shift_type": "inblocks"})
```

In [9]:

```
# NBVAL_IGNORE_OUTPUT
u_init, p_init = w.subfunctions
fig, axes = plt.subplots(nrows=2, sharex=True, sharey=True)
streamlines = streamplot(u_init, resolution=1/3, seed=0, axes=axes[0])
fig.colorbar(streamlines, ax=axes[0], fraction=0.046)
axes[0].set_aspect("equal")
axes[0].set_title("Velocity")
contours = tricontourf(p_init, 30, axes=axes[1])
fig.colorbar(contours, ax=axes[1], fraction=0.046)
axes[1].set_aspect("equal")
axes[1].set_title("Pressure");
```