数理解析研究所講究録989 巻 1997 年 78-91

78

Mesh Adaption and Automatic Differentiation for
Optimal Shape Design
Bijan Mohammadi
INRIA, BP 105, 78153 Le Chesnay, France
Bijan.Mohammadi@inria.fr
Abstract
A new approach for optimal shape design based on a CAD-free framework for shape andunstructured mesh deformations, automatic differentiation for the gradient computationand mesh adaption by metric control is presented. The CAD-free framework is shown
to be particularly convenient for optimization when the mesh connectivities and control
space size are variable during optimization. Constrained optimization for transonic regime
has been investigated.
Keywords: Shape Design, Unstructured Mesh Adaption, Automatic Differentiation.

1

Introduction

Local mesh adaption by metric control is a powerful tool for getting mesh independent
results at a reduced cost [1, 2, 3]. However, mesh adaption has never been used in optimization procedure because it implies a change of the desi space and mesh connectivities
after each adaption and this is difficult to take into account. Moreover, these procedures
are not differentiable at all. Therefore, optimization is usually performed on meshes with
fixed connectivities and on a given control space (i.e. no control parameters are createdduring optimization). However, as for direct problems, the need for adapted meshes isclear in design processes.
The aim of this paper is to show how to couple mesh adaption and shape optimization
tools. Here, not only the mesh in the computational domain, but also the number andrepartition of the control parameters over the shape will change during optimization.
Therefore, general shape and unstructured mesh deformation tools. In
particular, these tools have to be as much as possible CAD-free, in the sense that: theonly geometrical $entit?J$ available during
shonld be the mesh. Another reason
which ustifies this requirement is that CAD-based shape deformation tools
quite
complicated to use and also difficult to take into account when evaluating the Jacobian.
In [4] we showed that we can consider all the mesh points on the body (in both 2 and)asas far as we preserve the initial local regularity of the shape during
the design process. This makes possible and easy to consider all kind of shapes.
As we use amethod in our shape optimization, it is necessary to have a’
goodof the Jacobian of thecost function with respect of control$g\mathrm{n}$$\mathrm{a}\mathrm{l}\cdot \mathrm{e}\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{S}\mathrm{a}\mathrm{l}\backslash \mathrm{y}$$optimi_{\sim}^{\sim}ati_{\mathit{0}}n$$\mathrm{j}$$\mathrm{a}\mathrm{l}\cdot \mathrm{e}$$.3\mathrm{D})$$\mathrm{o}\mathrm{i}\mathrm{n}\mathrm{t}_{\mathrm{S}}$$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{o}\mathrm{l}_{1}$$\mathrm{g}_{1\partial \mathrm{d}\mathrm{i}\mathrm{e}\mathrm{n}}\mathrm{t}$’$\mathrm{e}\mathrm{v}\mathrm{a}1\iota\iota \mathrm{a}\{\mathrm{i}_{0}11$$\mathrm{c}\mathrm{l}\mathrm{i}\mathrm{S}\mathrm{C}\mathrm{l}\cdot \mathrm{e}\mathrm{t}\mathrm{e}$79parameters. It is clear that if the number of control parameters is large, an adjoint
method is necessary $[5, 6]$ . But classical adjoint methods are often built in continuous
level, after a linearization of the governing equations, and do not take into account for theinconsistencies of the numerical schemes (like numerical dissipation). In this, we use
automatic differentiation in reverse mode
to produce the Jacobian of the discrete
cost function. Hence, the non-consistencies of theschemes are naturally taken
into account in the Jacobian. Moreover, the cost of this evaluation is independent of thenumber of control parameters as in a classical adjoint method. Using this approach, wehave studied [4, 7, 8] the impact of using different gradients based on different numericalfluxes (Roe and Osher with and without a MUSCL type second order reconstruction[14, 17, 15] used for the solution of the Euler equations, on the optimization procedure.
We alsothat this approach does
suffer from a change in the nature of theequations (from parabolic to hyperbolic). Viscous turbulent configurations have also been
considered by tacking into account of a
two-equation turbulence model [9] and special
wall-laws (including pressure and convection effects [11, 10]) in the Jacobian. We noticed
that in an optimization procedure involving mesh deformation, wall-laws are more suitable
than low-Reynolds modeling (i.e. solving the flow up to the wall) as this former approachrequires much finer meshes for which conformal mesh deformation are hard to achieve.
In particular, we showed that once the difficulties coming from the nonlinearities of theoperators treated, it is easier to perform optimization for a turbulent flow configuration
than for the corresponding inviscid flow (i.e. at the sanle Mach number).
As this work includes several tools, only a short description of each of them is given. The
governing equations and the flow solver are described in section 2. Section 3 is devoted tothe shape and mesh deformation tools. A short description of our mesh adaption approachis given in section 4. Our adaptive optimization algorithm and related details including
the Jacobian
by AD in reverse mode are described in section 5. Section 6shows various results of adaptive shape optimization fronl transonic regime over airfoils.
$\mathrm{w}\mathrm{o}\mathrm{r}\mathrm{l}\backslash ’$$(\mathrm{A}\mathrm{D})$$\mathrm{n}\mathrm{u}\mathrm{n}\mathrm{l}\mathrm{e}\mathrm{r}\mathrm{i}_{\mathrm{C}}\mathrm{a}1$$)$$\mathrm{n}\mathrm{o}\grave{\mathrm{t}}$$\mathrm{s}\mathrm{h}\mathrm{o}\backslash 1^{\cdot}\mathrm{e}\mathrm{d}$$k-\epsilon$$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{u}\mathrm{p}\mathrm{U}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}$2Flow equations and solver

We denote by $J(x, U(X))$ the cost function we want to ninimize under geometric andaerodynanlic constraints. Here, indicates the geometrical description of a configurationand the flow pattern $(U(x))$ around this shape is the solution of the steady Euler systemof fluid dynamics in conservation form:$x$$\nabla.(F(U))=0$ ,(1)where is the vector of conservative variables (i.e. $U=(\rho,$) andrepresents the advective operator. This system has 4 equations inand the system isclosed using the equation of state $p=p(\rho, T)$ . For seek of simplicity, we only consider theEuler equations but the extension is already available for viscous turbulent flows $[4, 7]$ .
The flow solver is based on aapproach on). A Roe [14] approximate Riemann solver has been usedinmeshes (togetherreconstruction and Van Albada limiters [15]. The time )$U$$\rho\tau\iota\rho(arrow,C_{v}T+\frac{1}{2}|\vec{u}|^{2}))t$$F$$2\mathrm{D}$$\mathrm{f}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{t}\mathrm{e}-\backslash " 01_{11}\mathrm{m}\mathrm{e}- \mathrm{G}\mathrm{a}\mathrm{l}\mathrm{e}1\backslash \mathrm{k}\mathrm{i}\mathrm{n}$$\mathrm{g}_{11}1\mathrm{a}\mathrm{r}$$\mathrm{f}\mathrm{e}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{h}\mathrm{e}\mathrm{d}\mathrm{l}\cdot \mathrm{a}$$\iota\iota \mathrm{n}\mathrm{S}\mathrm{t}\mathrm{l}\cdot n\mathrm{C}\mathrm{t}\iota \mathrm{l}\mathrm{r}\mathrm{e}\mathrm{C}\mathrm{l}\mathrm{t}\mathrm{l}\cdot \mathrm{i}\mathrm{a}\mathrm{l}\mathrm{l}-$$3\mathrm{D}$$\mathrm{w}\mathrm{i}\mathrm{t}\mathrm{h}-\backslash \mathrm{I}\mathrm{I}\mathrm{T}\mathrm{S}\mathrm{c}^{1}\mathrm{L}$$\mathrm{e}\mathrm{n}\mathrm{d}\mathrm{e}\mathrm{n}\mathrm{t}$$\mathfrak{c}1\mathrm{e}_{1}$80equ-. ation$\frac{\partial U}{\partial t}+\nabla.F(U)=0$,

is marched in time to a steady state. Our time discretization is based on a 4-stage Rungescheme. Inflow and outflow boundary conditions are of characteristic type [16]and for outflow boundaries care has been taken to correctly treat subsonic
.
$[12,
13]$
More technical details can be found in.
$\mathrm{I}\nwarrow^{r_{\mathfrak{U}}}\mathrm{t}\mathrm{t}\mathrm{a}$$\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{f}\mathrm{i}\mathrm{g}\mathrm{u}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}_{011}\mathrm{s}$3Shape and mesh deformation tools

We describe the set of tools needed for shape and mesh deformation from a variation ofcontrol parameters:,$(x_{c}(n_{c}))$$\delta x_{c}arrow\delta x_{w}arrow\delta x_{m}$wheredenotes the discretization points on the geometry andthe internal
mesh nodes.
Forapplications, control points are usually fitted by splines. Splines have two
features:
1) They permit $n_{c}<<n_{w}$ .
2) They smooth the variations of control points when propagating to body discretization
points.
But in this work, following the solution (through a local metric),
control spacechanges during optimization as points may be created or removed on the body.
,the positions of the points have to be redefined after each adaption over the shape. Of
course, this is possible inbut much lnore difficult in. In all case, this increases thedifficulty by introducing another layer in the adaptive optimization. Moreover, generalsurface splines are quite complicated to handle, especially on unstructured meshes. To
avoid
difficulties and reduce the dependencies between two adaptions, wethe following general framework for shape deformation in 2 and:1. All the nodes on the shape are control points $(i.e. n_{c}=n_{w})$ .
2. To avoid oscillations, a smoothing opemtor is defined over the shape. This can be,
to solve the following system:
for instance, a few ’local ’ Jacobi$x_{w}(n_{w})$$x_{m}(l\mathrm{V})$$2\mathrm{D}$$\mathrm{t}\mathrm{l}\mathrm{l}\mathrm{e}$$\mathrm{T}\mathrm{h}\mathrm{e}\mathrm{r}\mathrm{e}\mathrm{f}_{0}1^{\cdot}\mathrm{e}$$2\mathrm{D}$$3\mathrm{D}$$3\mathrm{D}$$3\mathrm{D}$$\mathrm{i}\mathrm{n}\mathrm{t}_{1\mathrm{O}}.\mathrm{d}_{11}\mathrm{c}\mathrm{e}$$3\mathrm{D}$$itet^{\backslash }ati_{\mathit{0}}ns$$(I-\vee C\triangle)\delta_{\tilde{X}=}w\delta X_{w}$,(2)is the smoothed shape variation for the shape nodes andis the variation
given by tlle
tool. $By$ ’local’ we mean that if the predicted shape is locallysmooth, it rernains unchanged during this step.
In the past, we used this framework for optinlal shape design in 2 and $3\mathrm{D}[8,7,4,29]$ forcases with fixed connectivity and control space.
showed that configurations involving
several thousands of control points can be easily treated.
The importance of the smoothing step above can be understood by the following argument:
want } variation, if designs aof dimension $(n-1)$ ina domain.
Sobolev inclusions, we know that
. It is easy tounderstand that the
we use do not produce
variations$whe\uparrow\cdot e\delta.\tilde{\mathrm{t}}w$$\delta x_{w}$$optimi_{\sim}^{\sim}ation$$\mathrm{t}\prime \mathrm{V}\mathrm{e}$$\delta x_{w}\in C^{1}$ $(\Gamma)$$1\prime \mathrm{V}\mathrm{e}$$\mathrm{t}\mathrm{l}$$\mathrm{e}$$\Omega\in R^{\mathrm{n}}$$\mathrm{F}1^{\cdot}\mathrm{o}\mathrm{m}$$\mathrm{g}_{\mathrm{l}\mathrm{a}}\mathrm{d}\mathrm{i}\mathrm{e}\mathrm{n}\mathrm{t}_{1}\mathrm{n}\mathrm{e}\{\mathrm{h}\mathrm{o}\mathrm{d}$$\Gamma$$1\mathrm{n}\mathrm{a}\mathrm{n}\mathrm{i}\mathrm{f}_{0}1\mathrm{d}$$H^{n}(\Gamma)\subset C^{1}(\Gamma)$$\mathrm{n}\mathrm{e}\mathrm{c}\mathrm{e}\mathrm{s}\mathrm{s}\mathrm{a}1^{\backslash }\mathrm{i}1_{1,d}’\cdot c^{1}(\Gamma)$

81

and therefore we need to project them into
(an example of this is given in
picture (1) . This means that the projected variations
is solution of an elliptic$\delta x_{w}$$H^{n}(\Gamma)$$)$$(\delta\tilde{x}_{w})$Figure 1: Smoothed and unsmoothed shapes.
can see that the radient jumps
shocks and also prodnces a non-smooth shape in leading edge regions.
is the result. If we keep continue the shape becomes more andof the first iteration of themore irregular.
$l/Ve$$th\uparrow\cdot ough$$g$$Tll\dot{i}S$$optimi_{\sim}^{\sim}ation$system of degree . However, as we are using a
discretization, aorder elliptic
system is sufficient even inbecause the edges of the geometry (like a trailing edge)
are considered as being a geometrical constraint for the design. Therefore we project thevariations
only intoeven in.
known, we have to expand these variations
Once
the mesh. This is doneby solving a volunlic elasticity system which is also of the form of (2).
As these tools have to be differentiated, it is important to notice that both surfacic andvolumic elliptic systems are solved using iterative schemes. This is a real difficulty for AD
in reverse mode [8, 7, 4, 29] (see below).
$P^{1}$$n$$\mathrm{s}\mathrm{e}\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}\mathrm{d}$$3\mathrm{D}$$(\delta x_{w})$$H^{2}(\Gamma)$$3\mathrm{D}$$(\delta\tilde{x}_{w})$4$\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{l}\cdot \mathrm{a}\mathrm{l}\mathrm{l}$Mesh Adaption by Metric Control

We give a
description of our metric evaluation. The mesh generation technique is
not described here. It is based on a Delaunay algorithm $[22, 21]$ .
The key idea is to modify the scalar product used in the nlesh generator for distance,
surface (or volume) evaluations. Hence, the Delaunay approach will generate equilateralelements in a new nletric and not in the Euclidean one. This scalar product is based on
the evaluation of the Hessian of the variables of the problem. Indeed, for a
Lagrange
discretization of a variable , the interpolation error is bounded by:$\mathrm{b}_{1^{\backslash }}\mathrm{i}\mathrm{e}\mathrm{f}$$P^{1}$$?\ell$$\mathcal{E}=|u-\Pi_{h}u|_{0\leq cll^{2}}|D2\iota \mathrm{t}|_{0}$

where

$h$

is the element size,$\Pi_{h}$[the$P^{1}$

interpolation of(3),$\mathrm{t}\iota$

and$D^{2}u$its Hessian matrix.

82

This matrix being symmetric, we have:$D^{2}u=$$=\mathcal{R}\mathcal{R}^{-1}$,and its eigenvalues (always real). Using thiswhere is the eigenvectors matrix of:
information, we introduce the following metric tensor$R$$D^{2}\mathrm{c}\ell$$\lambda_{i}$$/\vee\iota$(4)$/\vee[=\mathcal{R}\mathcal{R}^{-1},$

where$\tilde{\lambda}_{i}=\min(\max(|\lambda_{i}|, \frac{1}{h_{\max}^{2}}),$$\frac{1}{h_{\min}^{2}})$

,

and
with
being the minimal and maximal edge lengths allowed in the mesh.
Now, if we generate, by a Delaunay procedure, an equilateral mesh with edges of length, the interpolation error is equi-distributed over the edgesof 1 in the metricprecisely, we haveof the mesh.
$h_{\max}$$h_{\min}$$\mathcal{E}$$/\mathrm{V}\{/(c\mathcal{E})$$a_{i}$$\backslash _{\perp}\mathrm{I}\mathrm{o}\mathrm{r}\mathrm{e}$$\frac{1}{c\mathcal{E}}a_{i}^{T}\mathrm{J}/Ia_{i}=1$.

(5)

VVe notice however that (3) leads to a global error while we would like to have a relative
one. We therefore use the following metric definition which takes into account not only
the dimension of the variables but also their magnitude:$\mathcal{M}=\frac{1}{c\mathcal{E}}|\frac{D^{2}(\Pi_{h}u)}{\max(|\Pi_{h}u_{h}|,\epsilon)}|_{0}$,(6)we have introduced the local value of the variable in the norm. is a cut-off which
defines the admitted difference of the order of magnitude between the larger and smallerscales we are looking for. Bellow , the phenomena is supposed impossible to be capturedby this grid and should therefore be modeled (as in turbulence modeling). This can also
be seen as looking for a nlore precise estimation in regions where the variable is small.
Another important consequence of this estimation is that it removes the dimensional
problems when intersecting metrics coming from different quantities.
Extension to:Forof equations, the previous approach leads to a metric for each variable andwe take theof all these metrics. Indeed, a metric can be seen through its
unit ball which
an ellipsoid in case of anisotropy (the unit ball of the Euclidean
metric is a sphere). lNIore precisely, for two metrics, we find an approximation of theirintersection by the following procedure:Letand$?,j=1,2$ the eigen-values and eigen-vectors of$j=1,2$ . Theintersection metricis defined by$\mathrm{w}1_{1}\mathrm{e}\mathrm{r}\mathrm{e}$$\epsilon$$\epsilon$$s\iota/ste\gamma nS$$\mathrm{s}\mathrm{y}_{\mathrm{S}}\mathrm{t}\mathrm{e}\mathrm{n}1\mathrm{s}$$\mathrm{i}\mathrm{n}\{\mathrm{e}\mathrm{l}\cdot \mathrm{s}\mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}$$\mathrm{b}_{\mathrm{e}\mathrm{C}\mathrm{o}\mathrm{n}}1\mathrm{e}\mathrm{s}$$\lambda_{i}^{i}$”$\tau_{i}^{j}$$\mathcal{M}_{i},$$(,\hat{\mathrm{W}})$$d^{\prime\hat{\vee}4}= \frac{/\hat{\vee}l_{1}+J\hat{\vee}\{2}{\underline{)}}.\cdot$

(7)

83

where
(resp.
) has the same eigen-vectors thaneigen-values defined by:$\hat{\mathcal{M}}_{2}$$\hat{\mathcal{M}}_{1}$, (resp.

$i=1,2$ .

,$\tilde{\lambda}_{i}^{1}=\max(\lambda^{1}, v_{i}\text{ノ}4il_{2}1\tau v^{1}i)$$\mathcal{M}_{1}$$\mathcal{M}_{2}$) but with
(8)

The previous algorithm is easy to extend to the case of several variables. As we said, one
advantage of (6) is to remove the problem of dimension when intersecting metrics comingfrom variables with different physical meaning and scale (for instance.pressure, densityand velocity).
Boundary layers and a priori discretization
The evaluation of the Hessian is done by the Green fornuula with Neumann boundary
condition. This does not lead to a suitable mesh for boundary layers and situations where
an a priori minimal discretization along the wall is required. Indeed, the position of thefirst layer of nodes can be quite irregular (this is important for viscous flows) and thediscretization along the wall can become too coarse for a good control space definition.
Another important ingredient therefore, is a mixed Dirichlet-Neumann boundary condition for the metric on wall nodes.
precisely, the eigenvectors for these nodes are
the normal and tangent unit vectors and the eigenvalues
are given by:$1\backslash \mathrm{I}_{0}\mathrm{r}\mathrm{e}$$(\lambda_{\mathcal{T}}, \lambda_{n})$$\lambda_{n}=?nax(\frac{1}{l\iota_{n}^{2}}, \lambda_{n}^{c})$,$\lambda_{\tau}=7na\backslash T(\frac{1}{h_{\tau}^{2}}, \lambda_{\eta}^{c})\sim$

’

where
and
prescribed values to monitor the minimal mesh refinement along
the wall in the normal and tangential directions andandare the computed values(coming
the Hessian). Along the wall the previous metric At $(x)$ is therefore replaced
by a new metric:,$h_{n}$$l\iota_{\tau}\mathrm{a}\mathrm{l}\cdot \mathrm{e}$$\lambda_{n}^{c}$$\lambda_{\tau}^{c}$$\mathrm{f}_{\Gamma \mathrm{O}}\mathrm{n}1$$\hat{\mathcal{M}}(x)$$\hat{\mathcal{M}}(x)=\tau\Lambda\tau^{-}1$

where
and

$\Lambda=dia_{\mathit{9}(}\lambda n’\lambda_{\tau})$$T=(??arrow(X),\vec{\mathcal{T}}(x))$.

From a practical point of view, for viscous applications
is small and function of the
Reynolds number while for inviscid cases
is taken large. This metric is slightly propagated through the flow by a smoothing operator.
In the past, we have applied these ingredients to inviscid and viscous laminar andturbulent flows involving shocks and boundary layers. Details on these techniques can befound in [1, 2, 3, 18].
$h_{n}$$h_{n}$

5

Optimization problem

We consider the following$\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{i}\mathrm{n}\mathrm{l}\mathrm{i}\mathrm{Z}\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}$problem:$\min_{x_{c}}J(X_{C}, U(x_{C}))$$E(x_{C}, [/(\vee.T_{C}))=0,$$\backslash \backslash \cdot \mathrm{h}\mathrm{e}\mathrm{l}\cdot \mathrm{e}’.\iota_{C}\in R^{7?\mathrm{c}}$$g_{1}(x_{c})\leq 0$

describe the control )direct (’geonletrical’, on$1$equation andconstraints.

$jC_{1,2}$$g_{2}(U(x_{C}))\leq 0$,

the flow, $E\in R^{N}$ the state) and indirect (on ( , such as a given lift)$\mathrm{a}\mathrm{r}\mathrm{a}\mathrm{m}\mathrm{e}\mathrm{t}\mathrm{e}\mathrm{l}\cdot \mathrm{S},$$\backslash ?_{C}$,,$U\in R^{N}$$U$$.T_{C}\mathrm{I}$84
Indirect constraints on have been taken into account in the cost function by penalty.
Geometrical constraints are of two types:),
1. shape deformation is allowed between two limiting surfaces (curves in
2. the volume of the shape can only increase.
The first constraint has been taken into account by projection and the second constraint
by penalty in the cost function.
$U$$2\mathrm{D}$5.1

The Gradient

As we said, in this approach all the nodes on the wall belong to the design space. This
makes out of reach an evaluation of the Jacobian by finite difference for instance. Introducing the Lagrangian $L=J+pE$ , by the optimality conditions we get:$\frac{\partial L}{\partial U}=\frac{\partial J}{\partial U}+p^{t}\frac{\partial E}{\partial U}=0$

whe.re to get

$p.$

’

(9)

,

we solve:$\partial E^{t}$$\partial J^{t}$$\overline{\partial U}\overline{\partial U}p=-$’which has the cost of one flow solution. And after substitution we have:$\frac{dJ}{dx_{\mathrm{c}}}=\frac{\partial L}{\partial x_{c}}=\frac{\partial J}{\partial\tau_{\mathrm{c}}}.+_{P}\frac{\partial E}{\partial x_{c}}t$.

. This can be
but we still need to compute
Here the cost is independent ofeasily done using AD in direct mode for instance [23]. In the present, to avoid developinga particular adjoint solver, we use AD in reverse mode.
$n_{c}$

5.2$\partial E/\partial x_{c}$AD in reverse mode

To get the Jacobian of the cost function with respect of control parameters, all the tools
described in the previous sections have to be differentiated. But the systems involved aresolved using iterative schemes. This is a real problem as the intermediate states have to
be stored for the adjoint computation $[7, 4]$ .
We use Odyss\’ee developed at INRIA [24, 25, 26]. In [27, 28, 29] the gradients obtained
by Odyss\’ee have been compared with those obtained using finite differences in similar
contexts. In this approach, the lines of the programs describing the relations between
the variation of the design variables and the cost function including the grid and the flow
equations are considered as constraints. We associate to each of them a Lagrange multipliers $(p)$ and we construct an augmented Lagrangian $(L)$ . The values of the Lagrange
multipliers are obtained from the condition that the first variations of with respect tointermediate variables vanish. The solution can be always obtained simply by back substitution (hence the notion of reverse mode). Once the Lagrange multipliers
evaluated,
the gradient of can be easily calculated. Technical details can be found in [4, 7, 29].
$L$$\mathrm{a}\mathrm{l}\cdot \mathrm{e}$

$L$

85
Limitations and keys points

5.2.1

As we said, in reverse mode a dual variable is introduced at each new substitution in the
program. This is a problem when the solver contains nested loops (for instance a time
integration loop with inside loops over nodes, segments and elements). The limitation of
the reverse mode comes therefore from the required nlemory for Lagrange multipliers andintermediate variables. More details about the limitations we encountered can be foundin [4, 7, 8, 29, 28]. For a efficient use of AD in reverse mode in our optimization process,the following key-remarks are essential:
1. For steady flows, the adjoint is independent of the intermediate solutions obtained
during time integration. Therefore, we only need to store one state if for a given shape
the steady solution is known.
2. Use interprocedural differentiation [26] which consists of replacing what is inside
nested loops by routines and differentiating these routines.
3. As we said, to compute the adjoint we always use the steady solution. Therefore, weneed to know when the adjoint convergence is sufficient for the optimization to converge.
Of course, we can integrate backward for a very long time but this will not be optimal. In
fact, it is sufficient for the norm of the gradient of the discrete cost function to be strictly
decreasing during optimization $[4, 7]$ ( being the optimization iteration):$n$$|\nabla J_{h}^{n+1}|<|\nabla J_{h}^{n}|$.

In other words, when computing the adjoint, we have two different numbers of time
step. One for the forward system (taken equal to one) and another for the adjoint(very large) and we leave the reverse time integration loop once the above criteria issatisfied.

$\mathrm{s}\mathrm{y}\mathrm{s}\mathrm{t}\mathrm{e}\ln$5.3

Adaptive optimization algorithm

The gradient is used in a gradient method to solve the optimization problem. Our minimization tool is quite simple. It is based on a gradient method with a fixed descent step.
But more sophisticated methods can be used as far as they do not require the solutionof large systems (proportional to ). Indeed, the number of nodes on the wall can bequite large, especially inapplications (several thousands in [4, 7]). Moreover, as inthis approach the number of control variables may vary during optimization due to mesh
adaption, the optinlization tool should not be sensible to this.
For seek of simplicity, the adaptive optimization algorithm is given with only one optimization after each adaption. But several iterations of optimization can be made on thesame nlesh. Indeed,use of the following criteria to detect the tinle when adaption canbe made:$n_{c}$$3\mathrm{D}$$\backslash \mathrm{v}\mathrm{e}$$J_{i}^{n} \leq J_{0}^{n}-\epsilon h^{2}\max$’where and denotethe adaption and the associated (on this level) optimization steps and is free. ?
is the maximum segment length allowed in the mesh andknown that all the segments in thehave the same length in the local metric M.
to notice that
Moreover, it isis usually chosen once for all. This criteria
has to be seen as a sufficient (not necessary) condition for the )
.
to
$n$

$i$

$\mathrm{r}\mathrm{e}\mathrm{s}1\supset \mathrm{e}\mathrm{c}\mathrm{t}\mathrm{i}_{\lambda}\cdot \mathrm{e}1\vee \mathrm{V}$$l$$\epsilon$$\mathrm{n}\mathrm{T}\mathrm{e}\mathrm{s}\mathrm{h}$$\backslash \backslash ’ \mathrm{e}$$\mathrm{i}\mathrm{n}\mathrm{t}\mathrm{e}\mathrm{l}\cdot \mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}g$$l_{l_{m\mathfrak{a}x}}$$\mathrm{t}\mathrm{i}_{1}\mathrm{n}\mathrm{i}_{\mathrm{Z}\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}}$$0_{1}$$\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{l}\backslash ;\mathrm{e}\vee \mathrm{r}\mathrm{g}\mathrm{c}$86Indeed, we use this criteria together with an a priori maximum number of optimization
iterations allowed between two adaptions. A theoretical background for the optimization with a variable discretization and justifications for the previous criteria under more
restrictive conditions are available in [30].
At step of adaption, we denote the mesh, the solution, the metric and the cost.
$J.(X_{ci}, U(X_{C}i))$ byand$i$$\mathcal{H}_{i},$$S_{i},$$J(x_{i})$$\mathcal{M}_{i}$$\mathrm{A}_{0\Gamma \mathrm{i}}\mathrm{t}\mathrm{h}\mathrm{n}1$$\mathcal{H}_{0},S_{0}$do, given,$i=0,$$..$, iadaptdefine the control space (wall nodes) :compute the gradient:$\mathcal{H}_{i}arrow x_{i}$$(x_{i}, \mathcal{H}_{i},S_{i})arrow\frac{dJ(x_{i})}{dx_{i}}$$. \tilde{x}_{i}=P_{i}(xi-\lambda\frac{dJ(x_{i})}{dx_{i}})$deforme the mesh :$(\tilde{x}_{i},\mathcal{H}_{i})arrow\tilde{\mathcal{H}}_{i}$,$(\tilde{\mathcal{H}}_{i}, s_{i})arrow\tilde{S}_{i}$$(\tilde{\mathcal{H}}_{i},\tilde{s}_{i})arrow \mathcal{M}_{i}$generate the new mesh using this metric:,,$(\tilde{\mathcal{H}}_{i}, \mathcal{M}_{i})arrow \mathcal{H}_{i+1}$

interpolate the previous solution over the new mesh :compute the new solution over this mesh :,,update the solution over the deformed mesh:compute the metric:,,$(\tilde{\mathcal{H}}_{i},\tilde{S}_{i}, \mathcal{H}_{i+1})arrow\overline{S}_{i+1}$$(\mathcal{H}_{i+1},\overline{S}_{i}+1)arrow S_{i+1}$,,

end do.
In the previous algorithm, is the projection operator which changes after each adaption as the control space changes and$P_{i}$.

$\frac{dJ(\tau_{i})}{dx_{i}}.\cdot.=\frac{\partial J}{\partial x_{i}}+(\frac{\partial J}{\partial U})\tau(\frac{\partial U}{\partial x_{i}})$

6

Results

In this section we present some results for constrained optimization problems at varion a workstation making
ous Mach number. All these computations have beenwith less than 64 Megabytes of memory in simple precision. These conabout 10
figurations, including mesh adaptions, gradient computations and flow solutions requireless than 2 hours. Concerning memory requireinents, due to the optimization describedin the previous section, the reverse mode requires approximately 10 times more memory
the direct solver. Anof the nlemory required by the direct solver is given
by 50
words (1 word $=32$ or 64 bits depending on architectures), where $N$ is theof nodes. This is more thall what is necessary in a structured solver because that.
the data structmore complicated in an
involved$\mathrm{p}\mathrm{e}1^{\cdot}\mathrm{f}_{0}\mathrm{r}\mathrm{m}\mathrm{e}\mathrm{d}$$\mathrm{M}\dot{\mathrm{e}}\mathrm{g}\mathrm{a}\mathrm{f}\mathrm{l}\mathrm{o}\mathrm{p}$$\mathrm{t}1_{1}\mathrm{a}\mathrm{n}$$\mathrm{e}\mathrm{s}\mathrm{t}\mathrm{i}\mathrm{n}1\mathrm{a}\mathrm{t}\mathrm{i}_{\mathrm{o}\mathrm{n}}$$\cross N$$\mathrm{n}\iota\iota \mathrm{l}\mathrm{n}\mathrm{l}\supset \mathrm{e}\mathrm{r}$$\mathrm{a}\mathrm{l}\cdot \mathrm{e}\mathrm{n}\mathrm{l}\mathrm{t}\mathrm{c}\mathrm{h}$$\iota 1^{\cdot}\mathrm{e}\mathrm{s}$$\iota\iota \mathrm{n}\mathrm{S}\mathrm{t}\mathrm{r}\mathrm{t}\mathrm{l}\mathrm{C}\mathrm{t}_{\mathrm{t}\mathrm{r}}\mathrm{e}\mathrm{C}\mathrm{l}\mathrm{a}\mathrm{p}\mathrm{P}\mathrm{l}\cdot \mathrm{o}\mathrm{a}\mathrm{c}\mathrm{h}$87, the different penalty
In the following examples, when global constraints arecoefficients in the cost function are initially chosen for the different quantities involved tohave variations of the same order of magnitude. During optimization, they are reduced.
with the same ratio than the costwith constraints on the lift and the volume of the
The same drag reduction
airfoil has been considered for three different airfoils at different Mach regimes. The costfunction for all three cases is given by:$\mathrm{p}_{\mathrm{l}\mathrm{e}\mathrm{S}\mathrm{e}}\mathrm{n}\mathrm{t}$$\mathrm{f}_{\mathrm{t}\mathrm{t}}\mathrm{n}\mathrm{c}\mathrm{t}\mathrm{i}_{0}\mathrm{n}$$\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{b}\mathrm{l}\mathrm{e}\ln$$J(x)=C_{d}+\alpha|C_{l}-C_{l}^{0}|+\beta|Vol-V_{\mathit{0}}l0|$,

andare the actual and initial lift coefficientsis the drag coefficient,where
the actual and initial airfoil volumes. Our aim is therefore to reduce
and $Vol$ andthe drag while the lift and the airfoil volume are kept unchanged. To be sure that thesolutions are mesh-independent, a final computation has been done on the final shapeuntil convergence.
are in free access under
and the mesh generator
The flow solver$C_{d}$$C_{l}\text{ノ}$$C_{l}^{0}$$Vol_{0}$$\mathrm{B}\mathrm{L}2\mathrm{D}$$\mathrm{N}\mathrm{S}\mathrm{C}2\mathrm{I}\backslash \mathrm{E}\vee$$\mathrm{f}\mathrm{t}\mathrm{p}://\mathrm{f}\mathrm{t}\mathrm{p}$6.1. inria.

$\mathrm{f}\mathrm{r}/\mathrm{I}\mathrm{N}\mathrm{R}\mathrm{I}\mathrm{A}/\mathrm{p}\mathrm{r}\mathrm{o}\mathrm{i}\mathrm{e}\mathrm{C}\mathrm{t}\mathrm{S}/\mathrm{G}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}/$.

Drag reduction for a transonic flow

The initial airfoil is the NACA 0012. The design takes place at Mach number 0.754 and2 degrees of incidence. The drag coefficient has been reduced by more than 10 percent
while the lift and volume slightly increased.

7

Concluding Remarks

A new approach for optimal shape design on variable control spaces and on meshes withconnectivities has been presented. The reverse mode of automatic differentiation
has been used to produce the Jacobian of the discrete operators. In this approach, themesh is part of the optimization procedure and function of the solution. This is an easy
way to avoid mesh dependencies in the optimization as well as in the direct problem.
Acknowledgments. The author
The author would like to thank professor $0$ . Pironneau for hiswould like also to thank N. Rostaing-Schmidt and C. Faure of INRIA Sophia Antipolis
for their assistance in using Odyss\’ee.
$\mathrm{v}\mathrm{a}\mathrm{l}\cdot \mathrm{i}\mathrm{a}\mathrm{b}\mathrm{l}\mathrm{e}$$\mathrm{s}\mathrm{u}\mathrm{p}\mathrm{p}\mathrm{o}\mathrm{r}\dot{\mathrm{t}}$References[1] H. Borouchaki. M.J. Castro-Diaz, P.L. George, F. Hecht, B. Mohammadi (1996),
Anisotropic adaptive mesh generation in two dimensions for CFD, 5th Inter. Conf.
on Numerical Grid Generation in Computational Field Simulations, Mississipi State

Univ.
$Anisot\Gamma opic$ Grid Adaption[2] M. Castro-Diaz, F. Hecht and B. Mohammadi.
Flows Simulations”, to appear infor Inviscid and$(1995)\text{ノ}.$$l^{r}\prime i_{S}cous$”$\mathrm{J}\mathrm{C}’|\mathrm{P}$88[3] H. Borouchaki, P.L.George, B. Mohammadi, (1996) Delaunay lVIesh Generation
Specifications. Part II.. Applications, Finite Element in Analysis anderned byissue
on mesh adaption.
specialDesign,$G_{ov}-$$P\prime Iet?’

ic$

Optimal[4] B. Mohammadi
Design Procedure for Inviscid and Viscous Turbulent Flows, submitted to IJNMF.
$(1996),A\wedge^{\gamma}ew$$Sl\iota ape$[5] 0. Pironneau (1984), Optimal Shape Design for Elliptic Systerns, Springer-Verlag,
New York.

[6]. A. JanTeSon (1994), Optimum
Report 803, Von$\mathrm{I}\backslash ^{r}\mathrm{a}\mathrm{r}\mathrm{m}\mathrm{a}\mathrm{n}$$Ae?^{\backslash }odyn$

amic Design via Boundary Control, AGARD

Institute Courses.

[7] B. Mohammadi (1996), Optimal Shape $Design_{f}$ Reverse Iode of Autornatic Differentiation and Turbulence, 35 th AIAA conference, Reno.
$\Lambda$[8] M. Hafez, B. Mohammadi, 0. Pironneau (1996), Optimum Shape Design using Automatic Differentiation in Reverse Mode , ICNMFD conference, Monterey.
[9] B.E. Launder and D.B. Spalding $(197‘ 2..),\Lambda/Iathemati_{C}al$ AIodels of Turbulence, Academic Press.

[10] B. Mohammadi and 0. Pironneau $(1994),Analy_{3}\backslash is$ of theModel, WILEY.

-Epsilon Turbulence$I\mathrm{i}’$[11] B. Mohammadi, 0. Pironneau (1996), Unsteady Separated Turbulent Flows CompuModel, Submitted to Computer Method in Applied
tations with MVall-Laws and.
Mechanics and Engineering.
$k-\epsilon$[12] A. Dervieux $(19\mathrm{s}5),steadyEuler$ Simulations using Unstructured Meshes, VKI lecture
series, 1884-04.
[13] B. Mohammadi (1994), CFD with
No.164.

$NSC\mathit{2}I\mathrm{i}E’$: an User Guide, Technical report INRIA[14] P.L.Roe (1981), Approximate Riemann $SolverS_{f}$ Parameters Vectors and Difference$SchemeS_{;}$

J.C.P. Vol.43.

[15] G.D.Van Albada, B. Van Leer $(1984),Flux$ Vector Splitting andfor the Euler Equations, ICASE 84-27.

$Runge-I\mathrm{i}^{r}utta$

Methods[16] J. Steger, R.F. Warming (1983), Flux Vector Splitting for the Inviscid gas dynamic
with Applications to Finite-Difference Methods, J. Comp. Phys. 40, pp: 263-293.

[17] S. Osher, S. Chakravarthy (1982), Upwind Difference Schemes for the Hyperbolic
systems of conservation laws, mathematics of Computation.
[18] F. Hecht, B. Mohammadi $(1997),\Lambda IeSh$ Adaption by MetricPhenomena and Turbulence, submitted to the 35th AIAA, Reno.

$c_{ont},\backslash ol$

for Multi-scale$\Lambda$

89[19] D. Vandromme (1983), Contribution \‘a la
et la pr\’ediction d’\’ecoulements
turbulents \‘a masse volumique $variable_{f}$ Ph.D. thesis, University of Lille.
$mode^{\text{ノ}}li_{S}ation$[20] R. Struijs, H. Deconinck, P. de Palma, P. Roe, G.G.Powel (1991), Progress on AIultidimensional Upwind Euler Solvers for Unstructured Grids, AIAA paper 91-1550.
[21] P.L. George, F. Hecht, E. Saltel (1990), Fully automatic mesh generator for, Impact of Comp. in Sci. and Eng., 2, pp.187-218.
mains of any$\mathit{3}dd_{\mathit{0}}-$$sl\iota ape$[22] P.L. George (1991), Automatic mesh generation. Applications to finite element
method, Wiley.
[23] J.M. Mal\’e, B. Mohammadi, N. Rostaing-Schmidt $(1996),Direct$ and Reverse AIodes
of Automatic $Differe\uparrow?tiation$ of Programs for Inverse Problems: Application to Optimum Shapes Design Proc. 2nd Int. SIAM Workshop on Compu, tational Differentiation, Santafe.
[24] J.C. Gilbert, G. Le Vey, J. Masse $(1991),La$ diff\’erentiation automatique de fonctions
par des programmes, Rapport de Recherche INRIA 1557.
$repre\text{ノ}sent\acute{e}eS$[25] N. Rostaing-Schmidt (1993),
d’optimisation en

automatique: Application \‘a un probl\‘eme, Ph.D. Thesis, University of Nice.

$Diffe^{\text{ノ}}\Gamma entiation$$m\acute{e}t_{6}’orologie$[26] C. Faure (1996), Splitting of Algebraic Expressions for Automatic $Differentiation_{f}$ In
proc. of the Second SIAM Inter. Workshop on
Differentiation, SantaFe.
$\mathrm{C}_{\mathrm{o}\mathrm{n}1}\mathrm{p}_{\mathrm{U}}\mathrm{t}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{a}1$[27] B. Mohammadi, 0. Pironneau $(1996),New$ Progress in $O_{l},timum$ Shapes $Desig?\iota$ , CFD
Review 95, Hafez-Oshima Eds, J. Wiley Pub.
[28] B. Mohamnladi (1996), Automatic Differentiation and Nonlinear PDE, ECCOMAS,
Paris.

[29] B. Mohammadi $(1996),Diff\acute{e}rentiati_{on}$ Automatique par Programme et Optimisationde Formes A\’erodynamiques, MATAPLI 07/96.
[30] E. Polack (1993), On tlre Use of Consistent Approximations in the Solutions of SemiInfinite $Optimi^{\sim}\sim ation$ and Optimal Control Problems, Math. Progranuming, Series B,
Vol. 62, No. 2, pp: 385-414.

90/$\nearrow$$\gamma\iota$’$(^{\mathrm{A}_{/}}\sim\grave{\overline{\mathrm{C}^{\overline{(}}}}//,\varphi\backslash /\nearrow_{\sim}1\nearrow\sim\backslash -$$.d\backslash$$\rangle$-$(\vee\backslash$

’

’.//

$/\sim’\backslash \prime’/\backslash \sim\backslash$$\backslash \backslash \vee^{arrow/}\prime’/-$$-_{-\prime}\backslash //$$\backslash .\sim.<\backslash -’\backslash ’\backslash \backslash \simeq\backslash \backslash \sim’\backslash \nearrow_{\backslash }/\backslash ’\backslash \backslash -^{--}\backslash \backslash /$$/–\sim/\sim/\angle\backslash \backslash b’\prime^{4}.\tau$,$.\backslash \backslash \acute{\acute{k}}- h\prime’\backslash \backslash \grave{\mathfrak{x}}’\backslash -\backslash /\backslash \backslash \sim\backslash \backslash$$\backslash$$\backslash -\vee/\backslash \cdot.’.\mathrm{t}\mathrm{k}_{\prec}^{\backslash \prime}\mathrm{i}‘--/,/,’-\wedge’\sim.\backslash \tau\backslash .\backslash ’\.\backslash ’\backslash ’\backslash ,\grave{j}’\backslash \underline{\mathrm{K}}\wedge^{\backslash }\vee^{-}\sim’\backslash ’/-\backslash /\vee^{\backslash \backslash \backslash }-\wedge\backslash \backslash \sim//$$/-$$\backslash //\sim\sim-^{-\prime’}\backslash \sim\backslash \backslash \backslash /\backslash /^{\backslash }$91$/fl_{k}$$\mathrm{c}$,$\mathfrak{S}$