Automated shape differentiation in the Unified Form Language

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.


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 [5], this differentiation exercise is often tedious and error prone. In [15], Schmidt introduces the open-source library FEMorph: an automatic shape differentiation toolbox for the Unified Form Language (UFL, [2]). 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 [16].
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 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. arXiv:1808.08083v2 [math.NA] 4 Apr 2019 2 Shape differentiation on the reference element A shape functional is a map J : U → R defined on the collection of domains U in R d , where d ∈ N denotes the space dimension. We follow the perturbation of identity approach [17] and for a vector field V ∈ W 1,∞ (R d , R d ), we consider the family of transformations {P s } s≥0 defined by P s (x) = x + sV(x) for every x ∈ R d . Note that P s is a diffeomorphism for any sufficiently small s ∈ R + .
For a domain Ω ∈ U, let Ω s := P s (Ω) and assume that Ω s ∈ U for s sufficiently small. The shape directional derivative of J at Ω in direction V is the derivative We say that J is shape differentiable in Ω if the directional derivative dJ(Ω)[V] exists for every direction V ∈ W 1,∞ (R d , R d ), and if the associated map dJ(Ω) : W 1,∞ (R d , R d ) → 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 where u Ps is a scalar function 1 . The subscript P s highlights the possible dependence of u Ps on the domain Ω s . The standard procedure to compute dJ is to employ transformation techniques and rewrite where DP s denotes the Jacobian matrix of P s , and u Ps • P s denotes the composition of u Ps with P s , that is, (u Ps • P s )(x) = u Ps (P s (x)) for every x ∈ Ω. Note that det(DP s ) > 0 for s sufficiently small. Then, by linearity of the integral, the shape derivative dJ is given by where d s ( · ) denotes the derivative with respect to s at s = 0. The term d s (u Ps •P s ) is often called the material derivative [4]. Its explicit formula depends on whether the function u Ps 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∈I be a partition of Ω such thaṫ ∪ i K i = Ω and such that the elements K i are nonoverlapping. Additionally, let {F i } i∈I be a family of diffeomorphisms such that F i (K) = K i for every i ∈ I, whereK denotes a reference element. This induces a partition {P s (K i )} i∈I of Ω s . To evaluate (1), standard finite element software rewrites it as Let {g i } i∈I be the collection of maps defined by Then, formula (5) can be rewritten as and taking the derivative of (5) with respect to s implies 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 In UFL the maps {g i } i∈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∈I . To be precise, this differentiation corresponds to a standard Gâteaux directional derivative, because the integrand in (6) corresponds to the following limit which can be interpreted as the Gâteaux directional derivative of This viewpoint is important to correctly implement this differentiation step in the existing pipeline in UFL (see Figure 1). We emphasize that this also enables computing higher-order shape derivatives by simply taking higher-order Gâteaux derivatives in (6).

Input: integrand in physical space
Estimate quadrature degree Pullback to reference element New: calculate derivatives with respect to coordinates Output: integrand in reference space

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, it does not work for elements such as the Hermite element that require different pullbacks for point evaluation and derivative degrees of freedom.

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 2 [13]. To shorten the notation, we define Example 1 Let the integrand be independent of Ω, i.e. u Ps = u for some function u. Then, the chain rule im- On the other hand, inserting u Ps = u into (6), we obtain the equivalent expression: Example code is shown in Listing 1. Functionals with domain independent integrands are used in applications including image segmentation [7] or, when u ≡ 1, to enforce volume constraints in shape optimization.
Example 2 Let {V h (Ω s )} s be a family of scalar finite element spaces such that the global basis functions 2 Examples using FEniCS [11,12] will be almost identical, modulo small differences in setting up initial conditions On the other hand, note that for anyx ∈K and for i ∈ I, it holds v Ps (P s (F i (x))) = v i (x), where v i is a linear combination of the local basis functions {b m } m∈M defined on the reference elementK. Therefore, and (6) becomes Listing 2 shows code for this case, using as v P0 the piecewise affine Lagrange interpolant of sin(x) cos(y).
Example 3 Let u Ps be the finite element solution to the boundary value problem −∆u Ps + u Ps = f in Ω s , ∇u Ps · n = 0 on ∂Ω s .
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 [5, Ch. 10, Sect. 5]. For this example, the Lagrangian is where e s (u Ps , v s ) := Ωs ∇u Ps · ∇v s + u Ps v s − f v s dx (15) stems from the weak formulation of the PDE constraint (13). The shape derivative dJ is equal to the shape derivative of L s (u•P −1 s , p•P −1 s ), 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 L s (u • P −1 s , p • P −1 s ) can be computed as in Example 2. The result is 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.

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 [8, p. 65]. Note that by computing shape derivatives as in (6), it is straightforward to combine shape derivatives of L s (u•P −1 s , p•P −1 s ) with standard Gâteaux derivatives with respect to u • P −1 s .

Code validation
We validate our implementation by testing that the Taylor expansions truncated to first and second order satisfy the asymptotic conditions In Figure 4, we plot the values of δ 1 and δ 2 for s = 2 −1 , 2 −2 , . . . , 2 −10 , and J as in Examples 1 and 3 from the previous section (we denote these functionals J 1 and J 2 respectively). The vector field V is chosen randomly. This experiment clearly displays the asymptotic rates predicted by (17). We have repeated this numerical ex- Taylor test periment 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 [18].

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 [14,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 Figure 3 (top). The pipe contains viscous fluid (with viscosity µ), which flows in from the left and is modelled 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 −ν∆u + u∇u + ∇p = 0 in Ω , div u = 0 in Ω , Here g is given by a Poiseuille flow at the inlet and is zero on the walls of the pipe The goal is to modify the central region of the pipe so that the shape functional J(Ω) = Ω ν∇u : ∇u dx 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 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 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 Figure 3 (bottom), the convergence history is in Figure  4. These results are compatible with those in [14, Sect. 6.2.3] and clearly indicate the success of the shape optimization algorithm.
Description of Listing 4 In line 2-4, we load the finite element mesh pipe.msh and extract the vertex coordinates. This mesh is generated with Gmsh [6] and is available as part of [18]. 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 [3] 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 [1]

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 [10], as well as compute second order derivatives of J and implement (quasi-)Newton methods [15]. 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.

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 non-linear PDEs. During shape differentiation, our code treats finite element functions and global functions differently. This behavior is correct and necessary to handle PDEconstraints 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 [4,9,19] that the weak form is superior when the state and the adjoint equations are discretized by finite elements.