Automated shape differentiation in the Unified Form Language

Abstract

We discuss automating the calculation of weak shape derivatives in the Unified Form Language (ACM TOMS 40(2):9:1–9:37 2014) by introducing an appropriate additional step in the pullback from physical to reference space that computes Gâteaux derivatives with respect to the coordinate field. We illustrate the ease of use with several examples.

1 Introduction

Physical models often involve functionals that assign real values to the solutions of partial differential equations (PDEs). For instance, the compliance of a structure is a function of the solution to the elasticity equations, and the drag of a rigid obstacle immersed in a fluid is a function of the solution to the Navier-Stokes equations.

This type of functional depends on the PDE parameters, and it is often possible to compute the derivative of a functional with respect to a chosen set of parameters. This derivative can in turn be used to perform sensitivity analysis and to run optimization algorithms with respect to parameters in the PDE.

The shape of the physical domain that is part the PDE-model (like the shape of the rigid obstacle mentioned above) is a parameter that is not straightforward to handle. Although we can compute the shape derivative of a functional following shape calculus rules (Delfour and Zolésio 2011), this differentiation exercise is often tedious and error prone. In Schmidt (2018), Schmidt introduces the open-source library FEMorph: an automatic shape differentiation toolbox for the Unified Form Language (UFL, Alnæs et al. 2014). The library FEMorph is based on refactoring UFL expressions and applying shape calculus differentiation rules recursively. It can compute first- and second-order shape derivatives (both in so-called weak and strong form), and it has been successfully employed to solve shape optimization problems (Schmidt et al. 2018).

This article presents an alternative approach to automated shape differentiation. The key idea is to rely solely on pullbacks and standard Gâteaux derivatives. This approach is more generic and robust, because it does not require handling of special cases. In particular, it circumvents the implementation of shape calculus rules, and required only a minor modification of UFL, because UFL supports Gâteaux derivatives with respect to functions. As a result, UFL is now capable of shape differentiating any integral that can be expressed in it.

This article is organized as follows. In Section 2, we revisit shape calculus and describe how to shape differentiate using standard finite element software. In Section 3, we consider three test cases and show how to compute shape derivatives using Firedrake and UFL. In Section 4, we describe code validation of this new UFL feature. In Section 5, we solve a PDE-constrained shape optimization test case with a descent algorithm implemented in Firedrake and UFL. Finally, in Section 6, we summarize the contribution of this article.

2 Shape differentiation on the reference element

A shape functional is a map $\mathrm {J}:{\mathcal {U}}\to \mathbb {R}$ defined on the collection of domains ${\mathcal {U}}$ in $\mathbb {R}^{d}$, where $d\in \mathbb {N}$ denotes the space dimension.

We follow the perturbation of identity approach (Simon 1980) and for a vector field ${\mathbf {V}}\in W^{1,\infty }({\mathbb {R}^{d}},\mathbb {R}^{d})$, we consider the family of transformations {P_s}_s≥ 0 defined by P_s(x) = x + sV(x) for every ${\mathbf {x}}\in \mathbb {R}^{d}$. Note that P_s is a diffeomorphism for any sufficiently small $s\in \mathbb {R}^{+}$.

For a domain ${\Omega }\in {\mathcal {U}}$, let Ω_s : = P_s(Ω) and assume that ${\Omega }_{s}\in {\mathcal {U}}$ for s sufficiently small. The shape directional derivative of J at Ω in direction V is the derivative

$$ \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}]:= \lim_{s\searrow 0}\frac{\mathrm{J}({\Omega}_{s})-\mathrm{J}({\Omega})}{s} . $$

We say that J is shape differentiable in Ω if the directional derivative dJ(Ω)[V] exists for every direction ${\mathbf {V}}\in W^{1,\infty }({\mathbb {R}^{d}},\mathbb {R}^{d})$, and if the associated map $\mathrm {d\mathrm {J}}({\Omega }): W^{1,\infty }({\mathbb {R}^{d}},\mathbb {R}^{d})\mapsto \mathbb {R}$ is linear and continuous. In this case, the linear operator dJ(Ω) is called the shape derivative of J in Ω.

To illustrate the shape differentiation of a shape functional, we consider the prototypical example

$$ \mathrm{J}({\Omega}_{s}) = {\int}_{{\Omega}_{s}} u_{{\mathbf{P}}_{s}} \mathrm{d}\boldsymbol{x} , $$

(1)

where $u_{{\mathbf {P}}_{s}}$ is a scalar function.¹ The subscript P_s highlights the possible dependence of $u_{{\mathbf {P}}_{s}}$ on the domain Ω_s.

The standard procedure to compute dJ is to employ transformation techniques and rewrite

$$ \mathrm{J}({\Omega}_{s})= {\int}_{\Omega} (u_{{\mathbf{P}}_{s}}\circ{\mathbf{P}}_{s}) \det({\mathbf{D}}{\mathbf{P}}_{s}) \mathrm{d}\boldsymbol{x} , $$

(2)

where DP_s denotes the Jacobian matrix of P_s and $u_{{\mathbf {P}}_{s}}\circ {\mathbf {P}}_{s}$ denotes the composition of $u_{{\mathbf {P}}_{s}}$ with P_s, that is, $(u_{{\mathbf {P}}_{s}}\circ {\mathbf {P}}_{s})({\mathbf {x}}) =u_{{\mathbf {P}}_{s}}({\mathbf {P}}_{s}({\mathbf {x}}))$ for every x ∈Ω. Note that $\det ({\mathbf {D}}{\mathbf {P}}_{s})>0$ for s sufficiently small. Then, by linearity of the integral, the shape derivative dJ is given by

$$ \begin{array}{@{}rcl@{}} \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] &= &{\int}_{\Omega} \text{d}_{\text{s}}((u_{{\mathbf{P}}_{s}}\circ{\mathbf{P}}_{s})\det{\mathbf{D}}{\mathbf{P}}_{s}) \mathrm{d}\boldsymbol{x} \\ &= &{\int}_{\Omega} \text{d}_{\text{s}}(u_{{\mathbf{P}}_{s}}\circ {\mathbf{P}}_{s}) + u_{{\mathbf{P}}_{0}} \text{div}({\mathbf{V}}) \mathrm{d}\boldsymbol{x} , \end{array} $$

(3)

where d˙s(⋅ ) denotes the derivative with respect to s at s = 0. The term $\text{d}_{\text{s}}(u_{{\mathbf {P}}_{s}}\circ {\mathbf {P}}_{s})$ is often called the material derivative (Berggren 2010). Its explicit formula depends on whether the function $u_{{\mathbf {P}}_{s}}$ does or does not dependent on Ω_s (see Section 3).

Next, we repeat the derivation of (3) in the context of finite elements and derive an alternative formula for dJ. Let $\{K_{i}\}_{i\in \mathcal {I}}$ be a partition of Ω such that $\dot \cup _{i} \overline K_{i} = \overline {\Omega }$ and such that the elements K_i are non-overlapping. Additionally, let $\{{\mathbf {F}}_i\}_{i\in \mathcal {I}}$ be a family of diffeomorphisms such that ${\mathbf {F}}_i(\hat {K}) = K_{i}$ for every $i\in \mathcal {I}$, where $\hat {K}$ denotes a reference element. This induces a partition $\{{\mathbf {P}}_{s}(K_{i})\}_{i\in \mathcal {I}}$ of Ω_s. To evaluate (1), standard finite element software rewrites it as

$$ \begin{array}{@{}rcl@{}} \mathrm{J}({\Omega}_{s}) &=& \sum\limits_{i \in \mathcal{I}}{\int}_{{\mathbf{P}}_{s}(K_{i})} u_{{\mathbf{P}}_{s}} \mathrm{d}\boldsymbol{x} \\ &= &\sum\limits_{i \in \mathcal{I}}{\int}_{\hat{K}} (u_{{\mathbf{P}}_{s}} \circ {\mathbf{P}}_{s}\circ {\mathbf{F}}_i) \vert\det({\mathbf{D}} ({\mathbf{P}}_{s}\circ {\mathbf{F}}_i))\vert \mathrm{d}{\hat{\boldsymbol{x}}} . \end{array} $$

(4)

Let ${\mathbf {F}}_{i}^{-1}$ denote the inverse of F_i, that is, ${\mathbf {F}}_{i}^{-1}({\mathbf {F}}_{i}({\mathbf {x}}))= {\mathbf {x}}$ for every ${\mathbf {x}}\in \hat {K}$, and ${\mathbf {F}}_{i}({\mathbf {F}}_{i}^{-1}({\mathbf {x}}))= {\mathbf {x}}$ for every x ∈ K_i. Since ${\mathbf {P}}_{s} = ({\mathbf {P}}_{s}\circ {\mathbf {F}}_{i})\circ {\mathbf {F}}_{i}^{-1}$ and P_s ∘F_i = F_i + sV ∘F_i, we can rewrite (4) as follows:

$$ \begin{array}{@{}rcl@{}} \mathrm{J}({\Omega}_{s}) = \sum\limits_{i \in \mathcal{I}}{\int}_{\hat{K}} \left( u_{({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i)\circ{\mathbf{F}}_{i}^{-1}} \circ ({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i)\right) \\ \vert\det({\mathbf{D}} ({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i))\vert \mathrm{d}{\hat{\boldsymbol{x}}} . \end{array} $$

(5)

Let $\{g_{i}\}_{i\in \mathcal {I}}$ be the collection of maps defined by

$$ g_{i}({\mathbf{T}}) := (u_{{\mathbf{T}}\circ{\mathbf{F}}_{i}^{-1}}\circ{\mathbf{T}})\vert\det({\mathbf{D}}{\mathbf{T}})\vert . $$

Then, formula (5) can be rewritten as

$$ \mathrm{J}({\Omega}_{s}) = \sum\limits_{i \in \mathcal{I}}{\int}_{\hat{K}} g_{i}({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i) \mathrm{d}{\hat{\boldsymbol{x}}} , $$

and taking the derivative of (5) with respect to s implies

$$ \begin{array}{@{}rcl@{}} \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] = \sum\limits_{i \in \mathcal{I}}{\int}_{\hat{K}} \text{d}_{\text{s}}(g_{i}({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i)) \mathrm{d}{\hat{\boldsymbol{x}}} . \end{array} $$

(6)

Equation (6) gives an alternative and equivalent expression for the shape derivative (3). However, to derive formula (3), it is necessary to follow shape calculus rules by hand, which is often a tedious and error prone exercise. Equation (6), by contrast, can be derived automatically with finite element software. Indeed, to evaluate J(Ω), standard finite element software rewrites it as

$$ \mathrm{J}({\Omega}) = \sum\limits_{i \in \mathcal{I}}{\int}_{\hat{K}} g_{i}({\mathbf{F}}_i) . $$

In UFL, the maps $\{g_{i}\}_{i\in \mathcal {I}}$ are constructed symbolically and in an automated fashion. Therefore, it is possible to evaluate dJ(Ω)[V] by performing the steps necessary for the assembly of J(Ω) and, at the appropriate time, differentiating the maps $\{g_{i}\}_{i\in \mathcal {I}}$. To be precise, this differentiation corresponds to a standard Gâteaux directional derivative, because the integrand in (6) corresponds to the following limit

$$ \text{d}_{\text{s}}(g_{i}({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i))=\lim\limits_{s\searrow 0}\frac{g_{i}({\mathbf{F}}_i + s{\mathbf{V}} \circ {\mathbf{F}}_i)-g_{i}({\mathbf{F}}_i)}{s} , $$

which can be interpreted as the Gâteaux directional derivative of g_i at T = F_i in the direction V ∘F_i (Hinze et al. 2009, Def. 1.29). This viewpoint is important to correctly implement this differentiation step in the existing pipeline in UFL (see Fig. 1). We emphasize that this also enables computing higher order shape derivatives by simply taking higher order Gâteaux derivatives in (6).

Fig. 1

Symbolic workflow in UFL to transform integrals from physical to reference space

Remark 1

Lagrange finite element global basis functions are obtained by gluing local parametric basis functions, that is, basis functions $\{{b_{m}^{i}}\}_{m\in {}}$ defined only on K_i and of the form $b_{m} = \hat {b}_m\circ {\mathbf {F}}_i^{-1}$, where $\{\hat {b}_m\}_{m\in {}}$ is the set of reference local basis functions, which are defined only on the reference element $\hat {K}$. If V lives in a Lagrange finite element space built on the partitioning $\{K_{i}\}_{i\in \mathcal {I}}$, it is possible to evaluate dJ(T)[V] by computing the Gâteaux derivative in (6) in the direction of the reference local basis functions $\{\hat {b}_m\}_{m\in {}}$ (instead of in the direction V ∘F_i) and summing these values. This allows us to fully rely on the symbolic differentiation capabilities of UFL.

Remark 2

The approach does not rely on the element being affinely mapped, but extends to elements that are mapped using a Piola transform such as the Raviart-Thomas or Nedelec elements. However, the current implementation does not work for elements such as the Hermite element that require different pullbacks for point evaluation and derivative degrees of freedom.

3 Examples

In this section, we consider three examples based on (1) that cover most applications. For these examples, we give explicit expressions of dJ using (3) and (6) and show how to compute dJ with the finite element software Firedrake² (Rathgeber et al. 2016). To shorten the notation, we define V_i : = V ∘F_i.

Example 1

Let the integrand be independent of Ω, i.e., $u_{{\mathbf {P}}_{s}}= u$ for some function u. Then, the chain rule implies that $\text{d}_{\text{s}}(u_{{\mathbf {P}}_{s}}\circ {\mathbf {P}}_{s})=\text{d}_{\text{s}}(u\circ {\mathbf {P}}_{s})=\nabla u\cdot {\mathbf {V}}$. Recalling $\text{d}_{\text{s}}(\det {\mathbf {D}}{\mathbf {P}}_{s}) = \text {div}({\mathbf {V}})$, we conclude that (3) becomes

$$ \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] = {\int}_{\Omega} \nabla u \cdot {\mathbf{V}} + u\text{div}({\mathbf{V}}) \mathrm{d}\boldsymbol{x} . $$

(7)

On the other hand, inserting $u_{{\mathbf {P}}_{s}}=u$ into (6), we obtain the equivalent expression:

$$ \begin{array}{@{}rcl@{}} &&{\kern-.7pc}\mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] = \sum\limits_{i \in \mathcal{I}}{\int}_{\hat{K}}\left( \vphantom{{\mathbf{D}} {\mathbf{F}}_i^{-1}}{\mathbf{V}}_{i} \cdot (\nabla u \circ {\mathbf{F}}_i) \right.\\ &&{}\left.+ (u\circ {\mathbf{F}}_i) \text{tr}({\mathbf{D}} {\mathbf{V}}_{i}{\mathbf{D}} {\mathbf{F}}_i^{-1})\right)\vert\det({\mathbf{D}}{\mathbf{F}}_i)\vert \mathrm{d}{\hat{\boldsymbol{x}}} . \end{array} $$

(8)

Example code is shown in Listing 1. Functionals with domain independent integrands are used in applications including image segmentation (Hintermüller and Ring 2004) or, when u ≡ 1, to enforce volume constraints in shape optimization.

Listing 1

Firedrake code to compute dJ from example 1 when u(x, y) = x² + y² − 1

Example 2

Let $\{V_{h}({\Omega }_{s})\}_{s}$ be a family of scalar finite element spaces such that the global basis functions $\{{B_{s}^{i}}\}_{i}$ of V_h(Ω_s) are of the form ${B_{s}^{i}} = B^{i}\circ {\mathbf {P}}_{s}^{-1}$, where $\{B^{i}\}_{i}$ are basis functions of V_h(Ω) and ${\mathbf {P}}_{s}^{-1}$ is the inverse of ${\mathbf {P}}_{s}$, that is, ${\mathbf {P}}_{s}^{-1}({\mathbf {P}}_{s}({\mathbf {x}}))= {\mathbf {P}}_{s}({\mathbf {P}}_{s}^{-1}({\mathbf {x}}))={\mathbf {x}}$ for every ${\mathbf {x}}\in \mathbb {R}^{d}$. Let $v_{{\mathbf {P}}_{s}}\in V_{h}({\Omega }_{s})$ and $u_{{\mathbf {P}}_{s}}=v_{{\mathbf {P}}_{s}}+\| \nabla v_{{\mathbf {P}}_{s}} \|^{2}$. Since

$$ v_{{\mathbf{P}}_{s}} \circ {\mathbf{P}}_{s} + \|(\nabla v_{{\mathbf{P}}_{s}})\circ {\mathbf{P}}_{s}\|^{2} = v_{{\mathbf{P}}_{0}} + \|{\mathbf{D}}{\mathbf{P}}_{s}^{-T} \nabla v_{{\mathbf{P}}_{0}}\|^{2} , $$

(9)

equation (3) becomes

$$ \begin{array}{@{}rcl@{}} &&{\kern-.7pc}\mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] = {\int}_{\Omega} (v_{{\mathbf{P}}_{0}} + \| \nabla v_{{\mathbf{P}}_{0}}\|^{2}) \text{div}({\mathbf{V}})\\ &&{\kern-.8pc}-2 \nabla v_{{\mathbf{P}}_{0}} \cdot ({\mathbf{D}}{\mathbf{V}}^{T} \nabla v_{{\mathbf{P}}_{0}}) \mathrm{d}\boldsymbol{x} . \end{array} $$

(10)

On the other hand, note that for any $\hat {\mathbf {x}}\in \hat K$ and for $i\in \mathcal {I}$, it holds $v_{{\mathbf {P}}_{s}}({\mathbf {P}}_{s}({\mathbf {F}}_i({\hat {\boldsymbol {x}}}))) = v_{i}({\hat {\boldsymbol {x}}})$, where v_i is a linear combination of the local basis functions $\{\hat {b}_m\}_{m\in {}}$ defined on the reference element $\hat K$. Therefore,

$$ \begin{array}{@{}rcl@{}} &&{\kern-.7pc}v_{{\mathbf{P}}_{s}} \circ {\mathbf{P}}_{s} \circ {\mathbf{F}}_i + \|(\nabla v_{{\mathbf{P}}_{s}})\circ {\mathbf{P}}_{s} \circ {\mathbf{F}}_i\|^{2} \\ &&{\kern-.9pc}= \hat v_{i} + \|{\mathbf{D}}({\mathbf{P}}_{s} \circ {\mathbf{F}}_i)^{-T} \nabla \hat v_{i}\|^{2} \quad \text{on } \hat K , \end{array} $$

(11)

and (6) becomes

$$ \begin{array}{lll} \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] = \displaystyle\sum\limits_{i \in \mathcal{I}} \displaystyle{\int}_{\hat{K}} \left( (\hat v_{i} + \|{\mathbf{D}}{\mathbf{F}}_i^{-T} \nabla v_{i}\|^{2}) \text{tr}({\mathbf{D}}{\mathbf{V}}_{i} {\mathbf{D}}{\mathbf{F}}_i^{-1}) \right.\\ \left.- 2 \nabla \hat v_{i} \cdot {\mathbf{D}}{\mathbf{F}}_i^{-1}{\mathbf{D}}{\mathbf{V}}_{i}{\mathbf{D}}{\mathbf{F}}_i^{-1}{\mathbf{D}}{\mathbf{F}}_i^{-T}\nabla \hat v_{i}\right) \vert\det({\mathbf{D}}{\mathbf{F}}_i)\vert \mathrm{d}{\hat{\boldsymbol{x}}}. \end{array} $$

(12)

Listing 2 shows code for this case, using as $v_{{\mathbf {P}}_{0}}$, the piecewise affine Lagrange interpolant of $\sin \nolimits (x)\cos \nolimits (y)$.

Listing 2

Firedrake code to compute dJ from example 2. Note that in this case, v does not depend explicitly on x and y

Example 3

Let $u_{{\mathbf {P}}_{s}}$ be the finite element solution to the boundary value problem

$$ -{\Delta} u_{{\mathbf{P}}_{s}} + u_{{\mathbf{P}}_{s}} = f \quad \text{in } {\Omega}_{s} , \quad\nabla u_{{\mathbf{P}}_{s}} \cdot{\mathbf{n}} = 0\quad \text{on } \partial {\Omega}_{s} . $$

(13)

In this case, the functional (1) is said to be PDE-constrained, and computing its shape derivative is less straightforward. The standard procedure is to introduce an appropriate Lagrangian functional (Delfour and Zolésio 2011, Ch. 10, Sect. 5). For this example, the Lagrangian is

$$ \mathrm{L}_{s}(u_{{\mathbf{P}}_{s}},v_{s}) := \mathrm{J}({\Omega}_{s}) + \mathrm{e}_{s}(u_{{\mathbf{P}}_{s}},v_{s}), $$

(14)

where

$$ \mathrm{e}_{s}(u_{{\mathbf{P}}_{s}}, v_{s}) := {\int}_{{\Omega}_{s}}\nabla u_{{\mathbf{P}}_{s}}\cdot \nabla v_{s} + u_{{\mathbf{P}}_{s}}v_{s} - fv_{s} \mathrm{d}\boldsymbol{x} $$

(15)

stems from the weak formulation of the PDE constraint (13). The shape derivative dJ is equal to the shape derivative of $\mathrm {L}_{s}(u\circ {\mathbf {P}}_{s}^{-1},p\circ {\mathbf {P}}_{s}^{-1})$, where u is the solution to (15) for s = 0 and p ∈ V_h(Ω) is the solution to an adjoint boundary value problem. The shape derivative of $\mathrm {L}_{s}(u\circ {\mathbf {P}}_{s}^{-1},p\circ {\mathbf {P}}_{s}^{-1})$ can be computed as in example 2. The result is

$$ \begin{array}{@{}rcl@{}} \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] = {\int}_{\Omega} (u + \nabla u\cdot \nabla p + up - fp)\text{div}({\mathbf{V}})\\ - p\nabla f\cdot{\mathbf{V}} - \nabla u\cdot ({\mathbf{D}}{\mathbf{V}}+{\mathbf{D}}{\mathbf{V}}^{T})\nabla p \mathrm{d}\boldsymbol{x} . \end{array} $$

(16)

For this example, we omit the equivalent formula on the reference element because of its length. However, as Listing 3 shows, UFL removes the tedium of deriving the shape derivative, and we can easily compute dJ.

Listing 3

Firedrake code to compute dJ from example 3 when f(x, y) = xy in (15)

Remark 3

With appropriate modifications, the same code can be use for functionals constrained to boundary value problems with Neumann or Dirichlet boundary conditions. For the Neumann case, it is sufficient to add the Neumann forcing term in line 6 of Listing 3. For the Dirichlet case, one needs to replace u with u+g in lines 6 and 7 (where g is the function defined in terms of X that describes the Dirichlet boundary condition) and impose homogeneous Dirichlet boundary conditions in lines 8 and 9.

Remark 4

To evaluate the action of the shape Hessian of a PDE-constrained functional, one can follow the instructions given in Hinze et al. (2009, p. 65). Note that by computing shape derivatives as in (6), it is straightforward to combine shape derivatives of $\mathrm {L}_{s}(u\circ {\mathbf {P}}_{s}^{-1},p\circ {\mathbf {P}}_{s}^{-1})$ with standard Gâteaux derivatives with respect to $u\circ {\mathbf {P}}_{s}^{-1}$.

4 Code validation

We validate our implementation by testing that the Taylor expansions truncated to first and second order satisfy the asymptotic conditions

$$ \mathrm\delta_{1}(\mathrm{J}, s) = O(s^{2}) \quad \text{and} \quad \mathrm\delta_{2}(\mathrm{J}, s) = O(s^{3}) , $$

(17)

where

$$ \mathrm\delta_{1}(\mathrm{J}, s) := \|\mathrm{J}({{\Omega}_{s}}) - \mathrm{J}({\Omega}) - s \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}]\| $$

and

$$ \begin{array}{@{}rcl@{}} \mathrm\delta_{2}(\mathrm{J}, s) := \|\mathrm{J}({{\Omega}_{s}}) &-& \mathrm{J}({\Omega})\\ &-& s \mathrm{d\mathrm{J}}({\Omega})[{\mathbf{V}}] - \frac{1}{2} s^{2} \mathrm{d^2\mathrm{J}}({\Omega})[{\mathbf{V}}, {\mathbf{V}}]\| . \end{array} $$

In Fig. 2, we plot the values of δ₁ and δ₂ for s = 2^− 1,2^− 2, $\dots , 2^{-10}$, and J as in examples 1 and 3 from the previous section (we denote these functionals J₁ and J₂ respectively). The vector field V is chosen randomly. This experiment clearly displays the asymptotic rates predicted by (17).

Fig. 2

Taylor test for examples 1 and 3. The convergence rates match the expected convergence

We have repeated this numerical experiment for many other test cases, including functionals that are not linear in u, functionals given by integrals over ∂Ω involving the normal n, and functionals that are constrained to linear and nonlinear boundary value problems with nonconstant right-hand sides and nonconstant Neumann and Dirichlet boundary conditions. In every instance, we have observed the asymptotic rates predicted by (17). The code for these numerical experiments is available at “Software used in ‘Automated shape differentiation in the Unified Form Language’ ()”.

5 Shape optimization of a pipe

In this section, we show how to use Firedrake and the new UFL capability to code a PDE-constrained shape optimization algorithm. As test case, we consider the optimization of a pipe to minimize the dissipation of kinetic energy of the fluid into heat. This example is taken from Schmidt (2010, Sect. 6.2.3). To simplify the exposition, we use a very simple optimization strategy. At the end of the section, we will comment on possible improvements.

The initial design of the pipe is shown in Fig. 3 (top). The pipe contains viscous fluid (with viscosity ν), which flows in from the left and is modeled using the incompressible Navier-Stokes equations. To be precise, let Ω be the shape of the pipe, Γ ⊂ ∂Ω be the outflow boundary of the pipe (that is, the end of the pipe on the right), and u and p be the velocity and the pressure of the fluid, respectively. Then, u and p satisfy

$$ \begin{array}{@{}rcl@{}} -\nu {\Delta} {\mathbf{u}} + {\mathbf{u}} \nabla {\mathbf{u}} + \nabla p &= 0 & ~~~~~~~~~~~~~~\text{in } {\Omega} ,\\ \text{div}~ {\mathbf{u}} &= 0& ~~~~~~~~~~~~~~\text{in } {\Omega} ,\\ {\mathbf{u}} &= {\mathbf{g}} &~~~~~\text{on } \partial{\Omega}\setminus {\Gamma} ,\\ p{\mathbf{n}} - \nu \nabla u\cdot {\mathbf{n}} & = 0 &~~~~~~~~~~~~~~\text{on } {\Gamma} . \end{array} $$

Here, g is given by a Poiseuille flow at the inlet and is zero on the walls of the pipe

Fig. 3

Initial (top) and optimized (bottom) shape of a pipe connecting a given inflow and outflow

The goal is to modify the central region of the pipe so that the shape functional

$$ \mathrm{J}({\Omega}) = {\int}_{\Omega} \nu \nabla{\mathbf{u}} : \nabla{\mathbf{u}} \mathrm{d}\boldsymbol{x} $$

is minimized. To solve this shape optimization problem, we parametrize the initial design with a polygonal mesh and update the node coordinates using a descent direction optimization algorithm with a fixed step size. As descent directions, we use Riesz representatives of the shape gradient with respect to the inner product induced by the Laplacian, i.e., at each step the deformation is given by the solution to

$$ \begin{array}{@{}rcl@{}} -{\Delta} {\mathbf{V}} &=& -\mathrm{d\mathrm{J}}({\Omega}) \quad\text{in } {\Omega}\\ {\mathbf{V}} &=& 0 \qquad\quad{\kern7pt}\text{on fixed boundaries.} \end{array} $$

(18)

This approach is also known as Laplace smoothing. To avoid degenerate results, we penalize changes of the pipe volume. The whole algorithm, comprising of state and adjoint equations and shape derivatives, is contained in Listing 4 and described in detail in the following paragraph. The optimized shape is displayed in Fig. 3 (bottom), the convergence history is in Fig. 4. These results are compatible with those in Schmidt (2010, Sect. 6.2.3) and clearly indicate the success of the shape optimization algorithm.

Listing 4

Firedrake code to optimize the shape of a pipe and minimize the dissipation of kinetic energy into heat. Lines 26 and 28 use the newly developed automatic shape differentiation

Description of Listing 4

In lines 2–4, we load the finite element mesh pipe.msh and extract the vertex coordinates. This mesh is generated with Gmsh (Geuzaine and Remacle 2009) and is available as part of “Software used in ‘Automated shape differentiation in the Unified Form Language’ (2019)”. Lines 5–8 define the Gramian matrix of the inner product employed to compute descent directions. In lines 9–14, we define the space of P2-P1 Taylor-Hood finite elements, which we use to discretize the weak formulation of the Navier-Stokes equations, and set up the functions containing the solutions to the state and adjoint equation as well as the test functions for the weak form. In lines 15–22, we define the weak formulation of the Navier-Stokes equations and certain parameters to prescribe the use of the MUMPS direct solver (Amestoy et al. 2000) to solve the linearized equations. In lines 23–29, we define the shape functional J, the functional describing the volume of the shape, as well as the Lagrangian and its derivative. In particular, note that the shape derivative of the Lagrangian can be computed with the simple command dL = derivative(L, X) in line 28. Without the new automatic shape differentiation capability in UFL, line 28 would have to be replaced with the following formula

dL = -inner(nu*grad(u)*grad(W), grad(v))*dx -inner(nu*grad(u), grad(v)*grad(W))*dx -inner(v, grad(u)*grad(W)*u)*dx +tr(grad(v)*grad(W))*p*dx -tr(grad(u)*grad(W))*q*dx +div(W)*inner(nu*grad(u), grad(v))*dx -div(W)*inner(div(v), p)*dx +div(W)*inner(div(u), q)*dx +div(W)*inner(v, grad(u)*u)*dx +nu*inner(grad(u), grad(u))*div(W)*dx -2*nu*inner(grad(u)*grad(W), grad(u)) *dx

In lines 30–35, we set up a function that updates the solution to the state and the adjoint equations. We emphasize that this shape optimization problem is not self-adjoint and that UFL derives the adjoint equation automatically. Note that, whenever the function solve_state_and_adjoint is called, the new values of the velocity u are stored in the file u.pvd (which can visualized using Paraview (Ahrens et al. 2005)). Finally, lines 36–46 contain the optimization algorithm: for 100 iterations, we compute the shape derivative and penalize volume changes (lines 38–40), compute the descent direction (lines 41–42), update the domain (line 45), and update the state and adjoint solutions (line 46).

Fig. 4

The value of the objective (plotted in linear scale) is reduced from approximately 0.61298 to 0.40506. The H¹-norm of the gradient (plotted in logarithmic scale) is reduced from 0.487274 to 0.000870

Remark 5

The optimization algorithm of Listing 4 is based on a simple optimization strategy and can be improved in several ways, at the mere cost of adding lines of code. For instance, instead of using a fixed step size and a fixed number of iterations, one could implement an adaptive step-size selection and stopping criteria. Additionally, one could experiment with different inner products to define descent directions (Iglesias et al. 2017), as well as compute second-order derivatives of J and implement (Quasi-)Newton methods (Schmidt 2018). Despite the room for improvement, we would like to stress that Listing 4 can be readily used for a 3D problem by simply passing a 3D mesh and changing the inflow boundary condition in line 18.

6 Discussion

We have presented a new and equivalent formulation of shape derivatives in the context of finite elements as Gâteaux derivatives on the reference element. While the formulation applies to finite elements in general, we have implemented this new approach in UFL due to its extensive support for symbolic calculations. This new UFL capability allows computing shape derivatives of functionals that are defined as volume or boundary integrals, and that are constrained to linear and nonlinear PDEs. During shape differentiation, our code treats finite element functions and global functions differently. This behavior is correct and necessary to handle PDE-constraints properly. In combination with a finite element software package, such as FEniCS or Firedrake, that takes as input UFL, this enables the entirely automated shape differentiation of functionals subject to boundary value problems. This notably simplifies tackling PDE-constrained shape optimization problems.

Compared to the existing shape differentiation toolbox FEMorph, our code does not compute shape derivatives in strong form because it neither relies on shape calculus differentiation rules nor performs integration by parts. However, in practice, we do not consider this a limitation as it has been shown in Hiptmair et al. (2015), Berggren (2010), and Zhu (2018) that the weak form is superior when the state and the adjoint equations are discretized by finite elements.

7 Replication of results

The code for the numerical experiments is available at Software used in ‘Automated shape differentiation in the Unified Form Language’ (2019).