{'fraction_non_alphanumeric': 0.02957050799777567, 'fraction_numerical': 0.010745477740342154, 'mean_word_length': 4.696198993851314, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Combining Sobolev Smoothing with Parameterized Shape Optimization

On the one hand, Sobolev gradient smoothing can considerably improve the performance of aerodynamic shape optimization and prevent issues with regularity. On the other hand, Sobolev smoothing can also be interpreted as an approximation for the shape Hessian. This paper demonstrates, how Sobolev smoothing, interpreted as a shape Hessian approximation, offers considerable benefits, although the parameterization is smooth in itself already. Such an approach is especially beneficial in the context of simultaneous analysis and design, where we deal with inexact flow and adjoint solutions, also called One Shot optimization. Furthermore, the incorporation of the parameterization allows for direct application to engineering test cases, where shapes are always described by a CAD model. The new methodology presented in this paper is used for reference test cases from aerodynamic shape optimization and performance improvements in comparison to a classical Quasi-Newton scheme are shown.


Introduction
It is well-known from the literature that Sobolev gradient smoothing improves the performance of aerodynamic shape optimization considerably [1,2], because it counteracts the so-called loss of regularity [3,4]. This is a phenomenon, where the smoothness of the unknown shape decreases progressively during the optimization. Naturally, using a very large design space, such as all surface mesh vertex positions, which is also called a free node parameterization, is more prone to this problem. There are several interpretations on why this occurs, one is not explicitly demanding the regularity of the search space as a constraint. Hence, using a coarser parameterization, which only produces shapes of sufficient regularity, circumvents this problem [5]. The most popular solution to solve this problem is to apply an elliptic partial differential equation to smooth the gradient [6,7]. This is equivalent to a reinterpretation of the gradient in a search space of higher regularity. Algorithmically, the resulting Sobolev optimization scheme corresponds to an approximative Newton scheme, where the reduced Hessian is approximated Ω Γ n v Figure 1: Flow domain and design shape Consider a vector of design parameters p ∈ S p ⊂ R np , describing the shape of Γ, where S p is a set of parameters resulting in feasible designs. These would normally be the coordinates of a shape parameterization. For example, Hicks-Henne functions [17], to deform the surface, or Free-Form Deformation (FFD) points [18], where control points of a regular box around the shape are deformed and the surface coordinates are moved using a smooth spline interpolation. In a broader context, any coordinates of a CAD-based representation of our geometry may be used here.
Definition 1 (discretized shape optimization problem). Let m ∈ R d·nm be the coordinates of a d-dimensional mesh triangulation of the domain Ω, F : R nu × R d·nm → R be a discretized objective function, M : R np → R d·nm be a parameterization of the computational mesh and H : R nu × R d·nm → R nu a discretization of the flow equations on Ω. Then is the discretized shape optimization problem.
Remark: For the rest of this paper we assume a simplification. The flow equation can be written in a fixed point form H(u, m) = 0 ⇔ u = G(u, m). For steady-state solutions, this is just a stationary point in the time marching scheme of the used finite volume solver. Now, we can define the Lagrangian function L : R nu × R n λ × R d·nm → R, We can now formulate Karush Note, that we use the expression D u for the Jacobian matrix of a function with respect to variable u. These are well-known first-order optimality conditions and a variety of gradientbased optimization techniques is available to solve those equations. In this paper, a discrete adjoint equation is set up using algorithmic differentiation, to ensure that the adjoint equation is always consistent to the flow and turbulence model on a discrete level [20]. Since we assume the existence and uniqueness of the flow solution, the value of the objective function is completely determined by the design parameters. This means, that it is possible to formulate Definition 1 as an unconstrained optimization problem min p∈SpF (p), whereF (p) := F (u(M (p)), M (p)). We use this formulation to calculate expressions for the gradient and the Hessian of the objective w.r.t. the design parameters p, in Section 3. It is important to remark thatF (p) always implicitly assumes the flow equations to be fulfilled. Additionally, if one assumes the adjoint equation to be fulfilled it is possible to show the following lemma. Therefore, the derivative ofF w.r.t. p is a reduced derivative via the implicit function theorem.

Hessians and Parameterization
The shape Hessian can be defined either in a continuous or in a discrete setting. In the introduction, we discussed how in engineering applications the geometry is usually not described by a shape Γ. Instead, a parameterization via a Computer-aided design (CAD) software is used. This allows the designer to apply changes in a more intuitive way, i.e., change the thickness or length of a component instead of moving individual mesh nodes. Also, for industrial applications such CAD-based representations are the common standard, as they allow multiple simulations with varying software tools, e.g., flow, elasticity, or heat transfer simulations, independent of certain mesh file formats. It is important to obtain optimized designs in these terms to ensure that they can be manufactured. For our application, this means we are dealing with a set of design parameters p ∈ R np , as in Definition 1. If we want to discretize the PDE from Theorem 3, we have two paradigms.
1. The mesh for the CFD computation is our discretization of the shape, implying that we need to formulate a discrete version of the shape Hessian on this level.
2. We are interested in sensitivities w.r.t. the parameters p since we rely on optimizing those values for industrial application.
Traditionally, one has worked with the second-order derivatives of the Lagrangian w.r.t. the mesh coordinates as the shape Hessian, i.e., D mm L(u, λ, m). We now discuss how to translate this into the situation of a parameterized mesh deformation, where we need to calculate D ppF (p), in such a way that it is compatible with Sobolev smoothing.
Theorem 1 (discrete parameterization of the shape Hessian). Assume thatF : R np → R, p → y is twice continuously differentiable, for p ∈ S p ⊂ R np in the space of possible parameters, then for the second order derivatives ofF the following equation holds Here, m k denotes the k-th component of the respective vector m = M (p).
Proof. For the gradient ofF , Lemma 1 gives a connection to the derivatives of L. Application of the chain rule to the i-th component yields Once again, the adjoint calculus with the Lagrangian is crucial to remove dependencies on the flow state u and to transform total into partial derivatives. Now take a look at the i, j-th component of the second-order derivative, i.e., the Hessian matrix. d 2 dp i dp jF (p) = d dp i d dp jF (p) Applying the product rule to each component of this sum results in d 2 dp i dp jF For the mixed second order derivative d∂ dp i ∂m k L(u, λ, m) the chain rule can be applied again to replace p i with terms of m.
Reordering the scalar multiplications yields d 2 dp i dp jF The appearing sums are just matrix-vector multiplications. This means that the overall expression is Overall, this concludes the proof.
The above can be interpreted as a special case of the generalized Faà di Bruno formula [13,14]. Theorem 1 connects the second-order derivatives with respect to the mesh m and the parameters p. The second term of the right-hand side of Equation (6) can be neglected in most cases, since many parameterizations used in industrial and engineering applications are linear in the design parameters. We will discuss the exact form of these derivatives and their role in more details in Section 5. This is a useful feature in a design process, where gradual changes are applied to study the behavior of the system. Naturally, one wants this to be equivalent to one large overall change, i.e., linear in the design parameter. Overall, it is reasonable to approximate the Hessian by the first term of the sum in Equation (6) only. To do a Newton step, all we are missing is the first-order derivative ofF (p). It can be expressed simply by multiplication of Jacobian matrices, as seen in Lemma 1, together with the chain rule. D pF (p) = D m L(u, λ, M (p))D p M (p) By using Equations (6) and (7) in a Newton step for the parameterized problem, i.e., is now transformed into Here D mm L(u, λ, M (p)) is the reduced, discrete Hessian with respect to the mesh. Together this gives an approximation for the Hessian with respect to the design parameters, thereby eliminating the conflicting objectives described above. Any suitable approximation of the shape Hessian on the mesh B ≈ D mm L(u, λ, M (p)) can now be applied for a Quasi-Newton step w in the parameter space via the formula

Sobolev smoothing on function spaces
The idea of Sobolev gradient smoothing is based on the connection between gradients in different function spaces. It is well known from functional analysis, that the gradient as a mathematical object is dependent on the scalar product of the underlying Hilbert space. To illustrate this point further, we will introduce the relevant spaces for our application [21]. Note, that we present the method for scalar-valued functions. When dealing with vector-valued functions, Sobolev smoothing can be generalized naturally by a componentwise application.
Definition 2 (function spaces). We define the scalar product (11) and call the vector space the space of square-integrable functions. In addition, we define the scalar product where ∇f (x) is the weak derivative of f (x) and the vector space the Sobolev space of weakly differentiable functions.
Furthermore, we know that H 1 (Ω) ⊂ L 2 (Ω) and that the directional derivative in an arbitrary direction v ∈ H 1 (Ω), defined as is a unique linear mapping on H 1 (Ω). Indeed by Riesz representation theorem, we know that the directional derivative and the gradient are directly linked together.
Theorem 2 (Riesz representation of the gradient). Let f : H → R be a differentiable function on a Hilbert space H, with inner product ·, · , then there exists a unique representative ∇f such that We call this representative the gradient of f.
One should note that in the case of a finite-dimensional, real-valued vector space, i.e., R n , with the standard scalar product, this leads to the well-known connection between the Jacobian matrix and the gradient.
Furthermore, since the gradient depends on the underlying Hilbert space we can mark this in the notation. ∇ L f : gradient in L 2 (Ω) The following Theorem 3 demonstrates the connection between the two gradients and introduces the basic equation for Sobolev smoothing.

Theorem 3 (Sobolev Smoothing).
For a function f ∈ H 1 (Ω), which fulfills the following equation holds Proof. Note that in Equation (20) the differential operator acts componentwise on the gradient. By Theorem (2), we know that the directional derivative can be expressed in terms of a scalar product with the gradient.
For the second integral in the equation, we use integration by parts with a zero Neumann boundary as stated in Equation (19).
This can be rewritten using scalar products, resulting in Since the directional derivative is a unique linear mapping on H 1 (Ω) and L 2 (Ω), and H 1 (Ω) ⊂ L 2 (Ω), by Riesz representation theorem we get the following equation Plugin this into the previous equation we get the following statement which implies the terms are equal outside of sets with measure 0 The connection to Newton methods becomes apparent if we take a look at the update step for an iterative optimizer. Assume we want to use Newton's method for minimizing f (x), then a single step looks like On the other hand, if we do a steepest descent step with the Sobolev gradient the resulting formula is Combined with Equation (20) we get that This motivates an approximation of the Hessian by Sobolev smoothing. How to proceed depends on whether one considers the motion of the volume vertices as independent, or as dependent on the motion of the boundary vertices only. Normally, the analytic Hessian only depends on the surface deformation, as movement inside the flow volume does not change the objective function [8]. Therefore, a discretized volume approximation will have a rank deficit and be ill-posed. However, when replacing the Hessian with a smoothing operator it seems possible to either smooth the sensitivities utilizing a volumetric operator or to smooth along the surface only. See the next Section 5 for more details on the surface and volume formulations.
Here, we should remark that this also influences the choice of boundary conditions. If we perform the Sobolev smoothing on the design surface Γ, we can keep the Neumann boundary condition from Theorem 3. If the Sobolev smoothing is performed on the whole volume mesh, one usually asks for ∇ H f = 0 on the outer freestream boundary, since only the design boundary should be deformed. Notice that this is a stronger assumption than having homogeneous Neumann boundaries everywhere. We demonstrate the formulation of the Sobolev smoothing algorithm on the surface by using the Laplace-Beltrami operator. Assume we have a discretization of this operator on the surface nodes (I − Γ ), then we can use it to approximate the reduced Hessian of the Lagrangian function If we combine this with Equation (10), we get the formulation for Sobolev smoothing on the parameter level This formula offers a way to formulate an approximated Newton scheme w.r.t the parameters, even though the design parameters do not necessarily come from a Hilbert space structure.
In addition, we can introduce weights 1 , 2 into the equation, to better adapt the operator to a given test case.
For ease of notation, we will denote the entire left-hand side matrix from this equation by At this point, we would also like to mention similar smoothing approaches found throughout the literature. Most of them consider applying a partial differential equation to the search direction to change the space in which we perform a descent step and thereby achieve a higher regularity [22]. One example is the paper of Kusch, et.al. [9] on the drag minimization problem for Euler equations, where the authors derive the operator symbol of the Hessian matrix. Another example is the work of Ta'asan [10]. While this only deals with potential flows and inviscid flows, e.g., situations governed by Euler equations, it offers detailed and deep insight into the theoretical side of the problem. A more recent paper by Müller, Kühl, et.al. [7] suggests solving a p-Laplace problem as a relaxation for the steepest descent step. This is done in a free node optimization setting and can deal effectively with situations where the optimum might contain points or edges in the shape. Finally, we also have to consider the choice of the weights 1 , 2 in Equation (26). One idea is to optimize them as additional parameters in an optimization algorithm, an idea that was proposed by Jameson [23]. A more in-depth mathematical analysis of the -values and a strategy to compute them is discussed by Kusch, et.al. in [24]. All of these approaches are, however, hard to compute in a practical implementation. Therefore, the values are usually determined by a parameter study.

Formulation for Surface and Volume Mesh
In the previous Section 4, we already stated that the Sobolev smoothing method might be formulated on the surface or the volume mesh. Here, we will extend on this point and show some implications of this choice. Recall the smoothing Equation (10), Depending on the formulation of the Hessian approximation B, we might get different formulations. In the case of a volume formulation, assume that the parameterization M consists of two parts. First, the surface parameterization and second, the mesh deformation from surface to volume Together, they form the complete design parameterization In Equation 6, the Hessian matrix of an individual mesh component w.r.t the parameters D pp M k appears. These component functions M k are combined mappings Therefore, we can apply the Faà di Bruno formula again to get Inserting Equation (33) into the original result from Equation (6) yields a formulation for the complete second derivatives, In many applications, this formulation can be simplified significantly. Many common surface parameterizations, e.g., the FFD boxes [18] or Hicks-Henne functions [17] used in this paper in Section 7, are linear. Therefore, the second order derivative terms of the surface parameterization D pp S l (p) will vanish. This leaves the formulation with some additional cross derivatives for the volume formulation, In addition, the mesh deformation between the surface and volume nodes is oftentimes linear as well. For example, the SU2 framework [25] used in this paper computes the deformation with a finite element stiffness approach. See Appendix A for more information on the implementation. Numerically, this means assembling a stiffness matrix and solving the corresponding linear system. From this, we can conclude that the second-order derivatives of the mesh deformation D ss V k (s) will vanish as well if an appropriate deformation method is used.
In conclusion, we can neglect the second-order derivatives of the mesh parameterization M (p) and formulate a parameterized Sobolev smoothing algorithm based on the first term in Equation (34). For the surface formulation, we end up with leading to the smoothing equation Accordingly, we get a similar equation for the volume case, where denotes the d-dimensional Laplace operator. For more complicated CAD models these linearity assumptions might not hold. In such a case, it would be necessary to compute the additional terms from Equation (34) and use them in our approximated Newton step.

Reduced SQP Optimizer
In Section 3, we saw how to smooth the gradient for the discrete shape optimization problem from Definition 1. However, we are still missing one important additional complication. In literally any industrially relevant aerodynamic application we have to deal with additional constraints in the optimization process, other than the mesh and flow equations. In this subsection, we will introduce the handling of additional constraints in the optimization process and demonstrate how our introduced Hessian approximation, by the Sobolev smoothing of gradients, fits naturally into a reduced sequential quadratic programming (SQP) framework [19,Chapter 18]. Examples of such constraints could be keeping certain lift or pitching moments for wings, or similar conditions. Definition 3 (constrained shape optimization problem). Use the terms from Definition 1 and let E : R nu × R d·nm → R n E be a set of additional constraints, then we call the constrained shape optimization problem.
With this we get an extended Lagrangian function L Ext : Here, the multiplier ν weights how much we value improvement in the objective against holding the constraints, and we can derive corresponding KKT conditions.
We want to solve these KKT conditions and a natural way to solve this system of nonlinear equations is to apply a Newton step [26].
In this equation, we have the Hessian of the Lagrangian function as the upper left part of the matrix. Also, we still have a dependence on u and p. Solving the adjoint problem and applying Lemma 1 we can remove the dependency on u and end up with reduced derivatives w.r.t. pD and simplify Equation (42) to Here we make use of the Hessian approximation with B defined from Equation (26). The relation to classical SQP methods becomes apparent if we take a look at the following quadratic optimization problem.
This quadratic approximation of the Lagrangian combined with linear constraints is known as the basic idea of SQP methods. The quadratic optimization problem has a unique solution (v, µ), fulfilling the relationD If we subtractD p Eµ from both sides we have Since we use the approximationD pp L Ext ≈ B and since the involved matrices are invertible we can identify the solutions of equations 44 and 48. This gives us At last, all we are missing is a way to calculateD p E(u, m). This gradient can be calculated by the adjoint solver implementation, similar to the gradient of the objective function. Therefore, we need additional calls to the adjoint solver for each constraint and adjoint constraint states To finally summarize our method, we have introduced a parameterized, discrete shape Hessian approximation, associated it with the Sobolev smoothing of free node gradients, shown how the Laplace-Beltrami operator can be used as a preconditioner, and demonstrated that this is the same as using a reduced SQP method with an approximated Hessian. With all those Algorithm 6.1 SQP method for equality constraint optimization input Initial design variables p 0 , set iteration counter i = 0

the design equations and compute reduced gradients
parts in place, Algorithm 6.1 shows the overall optimization algorithm. As the final step in this Section, we will state what an expansion of the algorithm for inequality constraints looks like. It is entirely possible to incorporate additional inequality constraints as well. We will not give the full analysis here, since it is beyond the scope of this paper. However, since we have an SQP method, we can just use the available adaptation from the literature [19,Chapter 18]. Assume that we have an additional set of constraints C : then the resulting optimization method is shown in Algorithm 6.2.

Hybrid Laplace-Beltrami operator
Having formulated the optimization algorithm, we have to reconsider the choice of the approximation forD pp L Ext (u, λ, ν, M (p)). Here, we will discuss two different aspects, not covered until now, to formulate our final expression for B. First, when dealing with additional constraints we saw how to extend the Lagrangian and the KKT conditions accordingly We propose using the parameterized formulation of Sobolev gradient reinterpretation, as described by Equation (26). Mathematically this means that Equation (27) is an approximation Algorithm 6.2 SQP method for mixed constraint optimization input Initial design variables p 0 , set iteration counter i = 0

the design equations and compute reduced gradients
of the Hessian of the original Lagrangian without additional constraints.
Neglecting the Hessian of the constraints will negatively affect the accuracy of our approximation. Furthermore, this is the only dependence of the optimizer on the nonlinearity of the constraint. Numerical experiments have shown that large design updates with this formulation have trouble abiding by the constraints, since the linearization in Equation (44) is no longer valid. Due to this, we propose the introduction of a regularization term The constant factor is ideally chosen as c = ν TD pp E(u, M (p)) . However, since computing this terms is as computationally expensive as computing the entire Hessian matrix, we need a simpler heuristic to determine c in practice.
Second, we have to ask ourselves about the structure and properties of the Laplace-Beltrami operator It can be split into two parts 1. The identity I is the canonical term to regularize the expression and stands for the standard L 2 scalar product when thinking about gradient representation.
2. The Laplace operator is how higher-order derivative information is introduced into Sobolev smoothing. It has strong smoothing properties onto the sensitivities. In the derivation of the gradient reinterpretation, we showed its connection to the H 1 scalar product.
To formulate the Laplace operator we need a spatial structure, i.e., an embedding in R 2 or R 3 . The design parameter space does not have this property, so it must be assembled on the computational mesh and then projected. In contrast to this, the identity can be defined as the componentwise identity for each of the design parameters. Combining both viewpoints we construct the system matrix as follows This expression is referred to as the hybrid Laplace-Beltrami operator. It is motivated by the two considerations above and is a fast and computationally cheap way to regularize our Hessian approximation. We already stated in Section 4 that optimal weights 1 , 2 could be computed analytically. However, the necessary terms are oftentimes too expensive to evaluate in practice. For practical applications of the algorithm, we found that a simple heuristic that ensures that the hybrid Laplace-Beltrami operator (55) remains a symmetric positive definite matrix works well enough.

One Shot
The One Shot method is a framework for the simultaneous solution of flow, adjoint, and design equations. Mathematically, it is a Newton-type method to solve the complete nonlinear system of equations in the KKT condition at once [27]. We can understand this better when we view it as a natural extension of the adjoint framework. First, start with a coupled fixed point iteration for the flow and adjoint equations, known as the piggyback iteration The following notation can be significantly simplified by introducing the so-called shifted Lagrangian function For the classical One Shot algorithm, the intermediate adjoint solution is directly used to update the mesh coordinates.
Such an iteration solves all three KKT conditions at once, hence the name One Shot method. Naturally, the convergence of this coupled iteration is strongly dependent on an appropriate choice of the preconditioners P i . Extensive convergence analysis for this formulation has been done in the past. We refer the reader to the works of Hamdi and Griewank [27], where the method is treated as a Quasi-Newton method for the KKT system.

Algorithm 6.3 Multistep One Shot with mesh update
input Initial values p 0 for i = 0, 1, . . . , For practical applications, Equation (58) makes for a poor algorithm, due to a couple of reasons.
• Updating the design only after several iterations would intuitively increase the accuracy of the involved functions and gradients. Such algorithms are known as multistep One Shot algorithms.
• For our purposes the mesh parameterization has to be incorporated in the update procedure as well. In industrial applications updating the geometry can be quite expensive in terms of runtime. The implementation has to project mesh sensitivities back onto the design parameters and additionally call CAD tools, mesh generators, or similar necessary software to compute an updated geometry.
Taking the arguments above into consideration the parameterized multistep One Shot algorithm 6.3 can be formulated. However, this method is still lacking a crucial component. For aerodynamic shape optimization, we need a way to incorporate additional constraints into the optimization, as introduced in Subsection 6.1. There we took the classical approach of computing converged flow and adjoint solutions and then performing a Quasi-Newton step for the design equation, and replaced it with a reduced SQP optimizer. We can combine this with the One Shot algorithm 6.3 by handing intermediate flow and adjoint states to the quadratic subproblem in the SQP algorithm. Combining Algorithms 6.2 and 6.3 into one method leads to the One Shot optimization algorithm 6.4 used in this paper.

NACA 0012
In this section, we will show results for our method applied to standard test cases from aerodynamic shape optimization. The first test case is the well-known NACA 0012 airfoil. This geometry dates back to the first half of the 20th century [28] and is one of the most widespread test cases for CFD in general. For the flow simulation, we use Euler equations, where the conditions are a Mach number of M = 0.8, an angle of attack α = 1.25 • , and standard air, resulting in a typical shock appearing on the upper side of the airfoil. The computational mesh contains 5233 points, 10216 elements, and two boundary markers, one for the airfoil and farfield boundary respectively. The resulting flow field can be seen in Figure 2. In fact, the drag value in this Euler simulation results from the shock. Therefore, assuming one wants to minimize the drag coefficient it is intuitively clear that the shape optimization must remove the shock. As a first test, we validate the finite elements solver for the Laplace-Beltrami operator by applying it to the surface sensitivities for the drag coefficient c D , according to Equation (20). In Figure 3, the upper part shows the L 2 gradient ∇ L c D , while the lower part shows the reinterpreted H 1 Sobolev gradient ∇ H c D on the surface. Both represent the directional derivative for their respective scalar products. As expected, Sobolev reinterpretation shows the classic properties of a Laplace-Beltrami operator, e.g., high sensitivities around the trailing edge are smoothed out. Of course, the main idea of this paper is to apply the Sobolev smoothing method to the design parameters. For this purpose, the airfoil is parameterized by 38 Hicks-Henne bump functions. 19 to deform the upper and 19 to deform the lower curve. The peaks of these functions are distributed equally along the chord at relative lengths of 0.05, 0.10, . . . , 0.95. For the optimization, we minimize the drag coefficient c D , in drag counts (1 count = 10 −4 ), such that the lift coefficient of 32.69 lift counts (1 count = 10 −2 ) is kept constant. The results of the optimization can be seen in Figure 4. We compare the optimization algorithm 6.1 with our approximated Hessian matrix to a simple constrained gradient descent method, with projection to the constraints. Also, to better understand how the choice of 1 and 2 effects the optimization, we test two different settings. In [24], it was suggested to choose a ratio 1 = 10 · 2 between the two parameters. This is test case dependent and a parameter study for the NACA 0012 case resulted in the choice of 1 = 1.0, 2 = 0.0625. To demonstrate this is a good choice, we show a second setting using 1 = 1.0, 2 = 0.625. With these settings, the Laplace-Beltrami operator is assembled on the design surface according to Equation (27). Figure 4: Comparison of the shape optimization for the NACA 0012 test case All methods achieve a similar level of improvement in the objective function, with the Sobolev reinterpretation of the gradients yielding a considerably faster convergence than gradient descent. For all methods, the final c D values for the optimized shapes are about 13 drag counts, which is within the range of spurious drag. The deviation in the lift constraint is more pronounced for the Sobolev methods. However, it stays within 0.6 lift counts of the target lift, which is within the numerical accuracy. Comparing the two Sobolev variants to each other, the first set of parameters is noticeably faster in the decrease of the objective than the second, at the cost of less than one lift count. Due to the inviscid nature of this test case, optimal solutions have to be shock free. The computed drag values in those shock free situation stem from numerical stabilization. To this end, Figure 5 shows the shape and pressure field for the starting geometry (left), the result of our Sobolev smoothed SQP optimization with 1 = 1.0, 2 = 0.0625 (center), and gradient descent (right). Each optimization algorithm efficiently removes the shock on the upper side of the airfoil and thereby the shock-induced drag. Upon closer inspection, one could argue that a negligible shock remains for the gradient descent algorithm. This is in correspondence to the spurious nature of the drag as discussed above.

ONERA M6
The NACA 0012 test in Section 7.1 gave a proof of concept for our algorithm on a wellunderstood reference test case. However, we are interested in applying preconditioning to more involved optimization problems. To demonstrate the potential for more complex geometries, and especially other parameterizations with higher numbers of design parameters, we take a look at the ONERA M6 test case. This is one  Figure 7. Wind tunnel data for the pressure coefficients on seven spanwise cross-sections are stated in the reference database [29]. Their positions in percentages of the total length are listed in Table 1 below. For each position 34 experimental measurements at different locations, along the cross-section, are given.  For this test, we use a structured C-type mesh for the ONERA M6 wing, consisting of 306, 577 points and 294, 912 rectangular elements. Around the wing the elements follow a logarithmic  For the flow simulation, the Reynolds-averaged Navier-Stokes (RANS) equations with a shear stress transport (SST) turbulence model [30] are solved with a Jameson-Schmidt-Turkel (JST) scheme [31]. Artificial dissipation coefficients for this test were set to d (2) = 1 2 and d (4) = 1 128 . For time integration we use a simple implicit Euler scheme. The parameterization is a three-dimensional FFD box around the wing, which has 10 cells in x-direction, 8 cells in y-direction, and 1 cell in z-direction, with the control points being only allowed to move in z-direction. Other deformations are suppressed to keep the span and chord length of the wing constant. Overall, this gives us 198 design parameters for the optimization. We start with a validation of the mesh and solver settings. A comparison between the CFD solution and experimental wind tunnel data for the original ONERA M6 test case is plotted in Figure 8. The CFD simulation can predict the shock position on the upper wing surface well and shows a good correlation to the experimental data. The accuracy is comparable to other relevant numerical studies featuring the same test case [32,33]. For cross-sections at 20%, 44%, 65%, and 90% simulation results correspond well with the experimental data.
Furthermore, both the SST and SA turbulence models yield very similar results. Nonetheless, the flow solver has some issues predicting the exact behavior around the merging double shock at 80% of the length and on the wingtip at 99%, due to limitations in the turbulence model. More accurate solutions would require the use of a more involved Reynolds stress tensor model, and for reference, we refer the reader to the work of Jakirlić, Eisfeld, et. al. [34]. In total, the test case setting gives us a decent resolution of the physical flow and is a stable basis for shape optimization. For optimization, we minimize the drag coefficient c D , with respect to a mix of equality and inequality constraints. The equality constraint keeps the lift coefficient of the original design c L = 24.99 lift counts and as inequality constraints, the minimal airfoil thickness is kept at the 5 different cross-sections listed in Table 2  Algorithm 6.2, with our approximated Hessian matrix, together with two other optimizers. First, a simple constrained gradient descent method with projection to the constraints, which can be obtained by simply using B = I in Algorithm 6.2. Second, we use the SLSQP optimization algorithm as a reference point [35]. This is the standard Quasi-Newton algorithm within the python SciPy optimize package [36]. It can be seen as a typical example of an SQP algorithm using gradient-based, iterative updates, e.g., a BFGS formula, to determine an approximated Hessian. For the reduced SQP algorithm, we use two different ways to assemble the hybrid Laplace-Beltrami operator. One on the surface mesh with weights 1 = 56.9, 2 = 0.9, and 3 = 0.1. For the other we assemble everything in the volume mesh with weights 1 = 0.0, 2 = 7.1, and  Results of the optimization can be seen in Figure 9. The objective and the constraints are plotted over the course of the optimization. We can observe that all three methods achieve an improvement in the objective function, from a starting value of 143.9 drag counts, albeit at different levels. The exact resulting values are shown in Table 3. Using reduced SQP with the smoothing matrix assembled on the surface results in an optimized drag value of 105.3 drag counts, this represents a 26.8% decrease in the objective function, while smoothing on the volume mesh results in a comparatively smaller improvement to 108.1 drag counts, which is a 24.89% decrease. In comparison, gradient descent results in a drag coefficient of 105.6 drag counts, i.e., a 26.6% drop. All of these optimizations result in very similar values, while the SLSQP algorithm can achieve an even smaller objective function value of c D = 102.4 drag counts, which represents a 28.85% decrease. At the same time, all of the optimization algorithms adhere to the optimization constraints very well. The lift constraint is kept up to a small error margin of less than 2.1%. Geometric inequality constraints are also kept and can be observed to be active constraints in the achieved minimum as well. Furthermore, we can observe how the optimizer oscillates between steps improving the objective function and steps increasing the feasibility, e.g., adherence to the constraints. This is a typical behavior for reduced SQP methods with approximated Hessian matrices, see for example [37]. They aim at reducing the Lagrangian function which may lead to some steps increasing the objective function to obtain better feasibility in the constraints. We can take a look at the flow fields around the final designs to achieve a better understanding of the minima found by the different optimization algorithms. A comparison is shown in Figure 10. In the upper left, we can see the original wing with the characteristic lambda shock on the surface. All of the used optimizers can remove most of the double shock system. This is important since shock-induced drag is an important source of the overall drag in this test case. A second visualization of the computed designs is shown in Figure 11, namely wing crosssections at 65% of the length, together with the surrounding flow fields. This gives an impression of the thickness and surface curvature for the different optimized designs. We can observe that the SLSQP algorithm results in a thinner leading edge of the wing, while the other optimizations result in rounder leading edges. This is a sign of the different methods finding different local minima, a common phenomenon in nonlinear optimization. Overall, we can observe that the application of our optimization method, with smoothing as a Hessian approximation, yields a more regular search direction and a faster rate of convergence for the optimizer if compared to gradient descent. Here, we can observe the effect of changing from a L 2 to a H 1 gradient, i.e., a reinterpretation of the derivatives in a different scalar product.   One Shot preconditioned One Shot grad. desc. time skylake (sec) 41280.76 39346.5 Another important factor to take into account is the runtime necessary to achieve the presented optimization results. The computational costs of the Hessian computation are negligible compared to the time needed to solve the flow and (multiple) adjoint simulations, as we will see now. For this, we run all algorithms on the same setting on a high-performance computing (HPC) cluster. The test architecture are two Intel skylake XEON SP 6126 nodes with 12 cores and 96 GB of RAM each, where the simulations are executed on 48 parallel MPI processes. We perform a number of 25 optimizer steps, or less if convergence is reached before. The measured runtimes can be seen in Table 4. The SciPy optimizer excels in this setting since it requires only 18 steps to converge and is thus much faster in overall time. Also applying the reduced SQP method with surface smoothing shows are very competitive performance in terms of runtime. It should be also taken into account that the flow and adjoint solutions are used as restart values in the next optimization. Therefore, the extra cost for the Hessian computation for the reduced SQP algorithm, compared to gradient descent, can be completely offset by smoother deformations leading to faster convergence in the next step.
Next, we assess the performance of our preconditioning method when applied in the context of One Shot optimization. Here, the hybrid Laplace-Beltrami operator from Equation (55) takes the place of the Hessian preconditioner in the One Shot algorithm 6.4. Before we present the results there are some settings for the algorithm to be considered. First, we must choose an adequate number of inner piggyback steps between design updates. Obviously, this will depend on the exact test case, to ensure that the flow and adjoints are recovered after each design update. For the presented ONERA M6 test case using 10 piggyback steps shows the desired behavior. Linked to this is the maximum allowed update step size. While the preconditioning itself enforces a scaling on the gradient, numerical experiments have shown the importance of further restrictions. Together with the number of piggyback steps, this choice helps to achieve a good combination of adjoint stability and design progress. We will assemble the Laplace-Beltrami operator on the surface mesh, since it is computationally cheap and showed better performance with converged gradients. As weights we use the same setting as previously, i.e., 1 = 56.9, 2 = 0.9, and 3 = 0.1. Once again we compare our approach to gradient descent and SLSQP. The gradient descent algorithm in the context of One Shot needs preconditioning as well to remain stable. Here, we restrict the max step size, which is equivalent to using a multiple of the identity as a preconditioner. The SLSQP algorithm is easy to set up as well, since it only requires handing the function values and gradients computed by piggyback to the SciPy interface.
We can see the results for One Shot optimization in Figure 12. From the plot, we can draw several observations. First, our Hessian approximation, shown in dark green, performs a successful optimization and achieves a minimal drag coefficient as low as 100.5 drag counts, with our proposed smoothing preconditioner. This represents a 30.14% improvement, while having c L = 24.985 lift counts, which is ca. 0.02% deviation from the target lift. Here, two settings for the steep sizes are shown. A heuristic is used to limit the maximal design update to the values given in Figure 12. As can be seen by the light green line, limiting the value from 5e −3 to 2.5e −3 slows the algorithm, but gives a similar result. Second, the gradient descent method, shown in blue, also performs a successful optimization, in accordance with the theory stating that for small enough update steps the One Shot algorithm must be convergent. The optimization results in an average drag coefficient of 101.6 drag counts, or 29.38% improvement, while showing small oscillating behavior in the end. Comparing the dark green and blue lines in Figure 12, it can be observed that constant preconditioners lead to higher oscillations for equal step sizes. Also, they tend to diverge faster if the maximal design updates are not limited carefully, as shown by the magenta line in Figure 12. This means that the use of an approximated, reduced shape Hessian as a preconditioner will allow the One Shot algorithm to perform larger design updates and result in a faster convergence to an optimal design. While at the same time keeping the additional time to assemble the Sobolev smoothing operator low, as can be seen from the runtimes measured in Table 5. The third observation is that a direct application of the SLSQP algorithm diverges after a few optimization steps. This can be explained as well since it is well known that BFGS updates will suffer from a high number of matrix updates. A problem, which is further increased by using approximated One Shot gradients of questionable precision. Two more general points to notice are how convergent One Shot optimizations adhere to all design constraints and how in general One Shot needs more optimization iterations than working with converged values. This is natural since we intend to do more, and in tendency smaller, design updates, albeit each at a fraction of the price for computing exact function and gradient values. After considering the optimization progress, we will investigate the resulting wing profiles in more detail. In Figure 13, the pressure coefficient on the surface is shown in relation to the original distribution. As we can observe, the One Shot method efficiently removes the characteristic lambda shock. Both optimizers remove the shock, thereby reducing the induced drag from the system. In Figure 14, we can see the cross-section at 65% of the wing length. Here, we can observe how the double shock gets reduced, and it is also worth noting how the final cross-sections are very close to the one found by the SLSQP algorithm, with converged values, from Figure 11. Next, we take a look at the runtime of the One Shot algorithm. All computations are done on an HPC cluster and we use the same setting as in the previous runtime test. The measured times can be found in Table 5. For this test, we perform 250 One Shot design updates with each algorithm. As we can see, the use of a One Shot approach leads to significant improvements in runtime. This proves that performing a higher number of cheap design updates can be advantageous compared to doing a few costly steps. To better summarize this point, we can calculate the so-called retardation factor. That is the quotient between the time for the complete optimization and one full convergence of the flow solver. Since converging the flow simulation requires 2224.1 sec. on the two skylake nodes, we end up with the factors shown in Table 6. The numbers presented here clearly show that applying One Shot can give crucial advantages in runtime. The presented Hessian approximation techniques can be used in this context as efficient preconditioning to ensure fast and stable convergence of the One Shot algorithm. All while keeping the computational cost low. If we compare the time retardation factors, we can conclude that the One Shot optimization algorithm achieves a full optimization in between 18 − 20 times the wall clock time needed for a simple flow simulation. At the same time, when working with exact gradients, even the competitive SLSQP implementation requires ca. 77 times the runtime at best. The retardation factors reveal that One Shot optimization with our proposed Sobolev smoothing is between 3.98 − 4.19 times more efficient, in terms of runtime than established Quasi-Newton methods with exact derivatives. For One Shot optimizers it is common practice to calculate retardation factors w.r.t. the iteration count, too. Here, the base value are the 15371 iterations necessary to fully converge the flow solver. For the optimizer, we need to converge the flow and adjoint solvers to a certain level, in order to gain initial restart solutions. Then for the given setting the One Shot optimization must perform 1 flow and 2 adjoint steps during each piggyback iteration and at the end, the flow is fully converged again. This results in the iteration retardation factors stated in Table 6. The discrepancy between the time and iteration retardation is largely due to the additional restarts and file I/O our implementation has to perform. Overall, we can conclude that the presented smoothing techniques are an effective preconditioning method for One Shot algorithms and result in a highly effective optimization strategy.

Conclusion
In this paper, we have shown how Sobolev smoothing further improves performance in aerodynamic shape optimization, even in situations where a smooth parameterization is used anyway. This is achieved by interpreting Sobolev smoothing as an approximation of the reduced shape Hessian, which is then incorporated into the parameterized setting via a tailor-made adaptation of the Faà di Bruno formula.
Although valid for all kinds of parameterizations, the resulting expressions based on our modified Faà di Bruno formula are particularly easy to implement for linear parameterizations. This makes the novel methodology presented herein convenient for a wide variety of industrial and engineering problems, since most of the common parameterizations, such as Hicks-Henne functions or FFD boxes, are typically linear. Furthermore, our modified Faà di Bruno formula can immediately be used as the KKT system of an approximated Newton-type optimizer, such as One Shot. The Laplace-Beltrami operator for Sobolev smoothing is extended by coefficients for zeroth and second-order smoothing. We compare this One Shot optimizer with a state-of-the-art SLSQP optimizer for industry-sized problems. This is demonstrated by a shape optimization of an ONERA M6 wing for drag reduction, subject to constant lift and thickness. For comparable gain in drag reduction, we achieve a time retardation factor of 18.56 with our One Shot approach, compared to 77.81 with the SLSQP optimizer. This represents a decrease factor of 4.19 in terms of runtime.

Acknowledgments
The authors acknowledge the support by Lisa Kusch concerning discussions on One Shot optimization and by Payam Dehpanah (both Scientific Computing Group, TU Kaiserslautern) for creating computational meshes for the numerical test cases. In addition, the authors acknowledge the productive discussions with Bernhard Eisfeld (DLR Braunschweig) and his valuable feedback on the ONERA M6 test case used in this paper. The presented computations were performed with resources provided by the Computing Center (RHRK) at TU Kaiserslautern.

Appendix A: Implementation
The methods discussed in this paper were implemented within the SU2 solver [25]. SU2 is an open-source multiphysics simulation and design software. For the scope of our work, the most important components are the available design optimization capabilities and the included discrete adjoint solver [20]. The adjoint solver is based on algorithmic differentiation, via operator overloading, and ensures that the adjoint solver is consistent for the governing flow and turbulence equations. We implement an efficient Sobolev smoothing approach, as introduced in previous sections, by providing additional functionality in two main areas.
1. An extension of the discrete adjoint flow solver to enable the assembly, solution, and output of the approximated Hessian matrix. For this, the Laplace-Beltrami operator is constructed on the mesh level and multiplications with the parameterization Jacobian are calculated.
2. The addition of an SQP optimizer capable of working with approximated Hessian matrices. Together with the existing design capabilities and the new Hessian approximation, this allows for efficient optimization. These capabilities are added to the existing python optimization toolbox FADO.

Gradient Smoothing Solver
The core of the presented implementation is a gradient smoothing solver capable of constructing the Hessian approximation from Equation (27) and handing it to an optimizer. To explain the key ideas of the implementation first, we need to introduce the basic structure of the SU2 C++ code. The current layout is based on a hierarchy of classes, which fulfill different roles for the multiphysics simulation. The general class hierarchy of SU2 is shown in Figure 15. For our purposes the two most important levels are: 1. driver classes: These are high-level classes derived from a common interface named 'CDriver'. They determine the kind of physical problem and options can include whether this is a flow, structural, or heat transport problem, or an adjoint run. They also handle the spatial layout of the problem being all on a single zone, e.g., only flow solutions, or having multiple zones, e.g., fluid-structure interactions. Finally, they control the temporal nature of the involved physics, which could be a steady-state problem with pseudo-time-stepping or a time-dependent problem. Here, we are interested in the 'CDiscAdjSingelzoneDriver' class, which runs the discrete adjoint simulation of a flow problem, thus computing the sensitivities of an objective function with respect to the mesh nodes.
2. solver classes: These classes implement a specific numerical solver and represent a certain set of governing equations. This could be flow, turbulence, structure, or any other kind of partial differential equation. All solver classes are based on a common interface named 'CSolver'. For our purposes, we are interested in performing operations on the solution of the adjoint flow solver class 'CDiscAdjSolver'.
Considering the structure, it is clear that to implement a Sobolev gradient treatment and Hessian approximation, we need to adapt the existing discrete adjoint framework on those two levels. First, we implement a new solver class 'CGradientSmoothingSolver', which has to assemble and solve the discretized PDE from Equation (26). From the structure of the Laplace-Beltrami operator, it is obvious that finite elements can be used to get an efficient discrete representation. This is done by extending the existing finite element implementation in SU2, provided by the structural mechanic features of the code. Several new functions are necessary to overcome two main hurdles. First, we need to extend the numerical integration to use second-order Gauss quadrature. This is necessary to accurately represent the identity with finite elements. Second, the existing elements and transformations are extended to deal with surface cells. For example, a two-dimensional surface cell is embedded in a three-dimensional mesh, meaning the transformation to reference elements maps R 3 → R 2 and this change in dimension requires some extensions to the existing implementation. Overall, this solver assembles the discrete equivalent of the Laplace-Beltrami operator on the computational mesh.  Next, the parameterization Jacobian D p M (p) and its transposed must be available in some form. Normally, all driver and solver classes exist within an executable SU2_CFD_AD and the 'CDiscAdjSinglezoneDriver' only computes sensitivities with respect to the mesh coordinates. The chain rule evaluation to get sensitivities w.r.t. design parameters is then handled by a different executable named SU2_DOT_AD, as is depicted in Figure 16. This modular setup allows the adjoint solvers to be independent of the choice of parameterization, yet for our implementation, it is necessary to transfer functions for the projection of sensitivities from SU2_DOT_AD. Our implementation enables the computations of all the steps inside of a single instance of the class 'CGradientSmoothingSolver', including the computation of the Jacobian matrix of the mesh parameterization. Reviewing Equation (26) carefully shows that the full matrices are not needed at all.Consider the involved dimensions, D p M (p) ∈ R (d·nm)×np , ( 1 I − 2 ) ∈ R (d·nm)×(d·nm) , and B ∈ R np×np . We expect n p << n m and therefore it is much more efficient to only construct the complete matrix B. If we have a way to compute arbitrary scalar products with D p M (p) and D p M (p) T this can be achieved. Algorithmic Differentiation (AD) is an efficient way to calculate derivatives [38], and is especially efficient if only matrix-vector products with Jacobians are required. AD offers two basic options called the forward and reverse mode. Given a function M : R np → R d·nm and arbitrary vectorṡ x ∈ R np ,ȳ ∈ R d·nm we can computė y = D p M (p)ẋ (forward mode) x = D p M (p) Tȳ (reverse mode). (A.1) In our implementation, we use the AD tool CoDiPack [39] to efficiently calculate the derivatives. CoDiPack is already the standard AD tool in SU2 and is heavily utilized to solve the adjoint equation. Here, we use the available CoDiPack features to efficiently evaluate the two expressions in Equation (A.1). With all these functions in place, we can assemble and solve Equation (26) in SU2, as well as extract the left-hand side matrix for use as a Hessian approximation. Finally, the existing 'CDiscAdjSinglezoneDriver' class is extended in several ways. We add an instance of the new 'CGradientSmoothingSolver', and couple it with the discrete adjoint solver as a postprocessing step. Together, they enable us to calculate D p L and B with one adjoint simulation run.

Piggyback Driver
The One Shot algorithm 6.4 is built around the idea of a piggyback iteration. This means we compute a coupled iteration of flow and adjoints, where the current flow state u j is used to compute the current adjoint state λ j . While we will not go into specific details about the implementation, it is worth pointing out the key idea and how this can be set up. Within SU2 the general adjoint logic is included within the driver classes, derived from a common interface 'CDriver', see Figure 15 for the class hierarchy. Driver classes represent a certain type of physical problem, e.g., flow, structure, adjoint, etc., and the driver controls the time-stepping process and calls the involved solvers. For many adjoint solvers, including the SU2 implementation, adjoints are computed using the idea of reverse accumulation. The algorithm is based on a theorem by Christianson regarding AD of fixed-point iterations [40]. Usually, when applying the reverse mode of AD to an iterative procedure, e.g., u j = G(u j , m), one needs to store all intermediate states u 1 , . . . , u J for the derivative computation. Of course, this results in unfeasible memory overheads. The theorem states that for a converged fixed point solution u * it is sufficient to store only this fixed point to compute the exact derivatives.
flow solver u j = G(u j , m i ) In terms of data flow, the situation is depicted in Figure 17. First, an instance of the flow driver executes a pseudo-time-stepping to solve the flow equation, then the resulting solution u * is handed to the adjoint driver. There an iteration is set up using reverse-mode AD and then executed until a converged adjoint state λ * is reached. For our implementation, we introduced a new piggyback driver class named 'COneShotSingle-zoneDriver'. Following the idea presented in Equation (56) the data flow is changed to the layout depicted in Figure 18. The new driver performs a piggyback iteration, by doing one step of the flow solver and then handing the intermediate result u j to the adjoint solver to do one adjoint step using reverse-mode AD, and returns λ j . Those two steps are then repeated to perform a combined fixed point iteration. Using these capabilities, we are able to compute all necessary reduced gradientsD p F,D p E k ,D p C l from Algorithm 6.4 with the SU2 adjoint framework. Please keep in mind that these principles apply to all AD-based adjoint solvers when used for One Shot optimization. By using a consistent piggyback implementation we ensure that function and gradient values, while inaccurate, are still computed consistently with one another.

Optimizer implementation
In Section 6, we introduced an SQP optimization strategy using Sobolev gradient smoothing for design parameters to compute an approximated Hessian. Here, we will discuss the implementation of Algorithms 6.1, 6.2, and 6.4. The flow and adjoint simulations in this paper are done using the SU2 framework. SU2 already provides a number of python scripts to simplify various tasks, by a python interface. These scripts allow the computation of flow and adjoint simulations, using the different SU2 C++ executables, and the evaluation of functionals and gradients, respectively. The results can be handed to the SciPy optimization library [36] by the interface. However, the provided optimization capabilities are insufficient to run Algorithms 6.2 and 6.4 for several reasons. The major concern is that we want to provide an approximated Hessian to the SQP optimizer, while most software packages either expect an exact Hessian or numerically approximate the Hessian themselves. To be able to provide our Hessian approximation we decided to implement the core iterative loop of the SQP method ourselves and use existing packages only to solve the individual steps of the optimization algorithms. To achieve this we use the FADO toolbox 1 . This toolbox provides a python interface to call SU2 simulations, to store and control different designs and meshes, or to adapt config options. It is deliberately open for the addition of different optimizers and we added an implementation of the basic SQP algorithm to run our process. In addition, in Algorithms 6.2 and 6.4 we have to deal with a more complex constrained quadratic program. For this, we use the CVXOPT python package 2 for convex programming [41,42]. The Hessian approximation in Equation (27) is symmetric, positive defi-nite, which implies that the quadratic subproblems for SQP are convex, so this is a valid choice.