\documentclass[twoside]{article} \usepackage{epsfig} \pagestyle{myheadings} \markboth{ Use of discrete sensitivity analysis }{ Clarence O. E. Burg } \begin{document} \setcounter{page}{13} \title{\vspace{-1in}\parbox{\linewidth}{\footnotesize\noindent Fourth Mississippi State Conference on Differential Equations and \newline Computational Simulations, Electronic Journal of Differential Equations, \newline Conference 03, 1999, pp 13--27. \newline http://ejde.math.swt.edu or http://ejde.math.unt.edu \newline ftp ejde.math.swt.edu or ejde.math.unt.edu (login: ftp)} \vspace{\bigskipamount} \\ % Use of discrete sensitivity analysis to transform explicit simulation codes into design optimization codes \thanks{ {\em Mathematics Subject Classifications:} 76N25, 49Q12. \hfil\break\indent {\em Key words:} design optimization, sensitivity analysis, adjoint methods, \hfil\break\indent computational fluid dynamics. \hfil\break\indent \copyright 2000 Southwest Texas State University and University of North Texas. \hfil\break\indent Published July 10, 2000. } } \date{} \author{ Clarence O. E. Burg } \maketitle \begin{abstract} Sensitivity analysis is often used in high fidelity numerical optimization to estimate design space derivatives efficiently. Typically, explicit codes are combined with the adjoint formulation of continuous sensitivity analysis, which requires the derivation and solution of the adjoint equations along with appropriate boundary conditions. However, for implicit codes, which already calculate the Jacobian matrix of the discretized governing equations, the discrete approach of sensitivity analysis is relatively easy to implement. Using the complex Taylor's series expansion method to generate derivatives, a highly accurate approximation to the Jacobian matrix can be generated for implicit or explicit codes, allowing uniform application of discrete sensitivity analysis to both implicit and explicit codes. \end{abstract} \section{Introduction} High fidelity numerical simulation codes have been developed to solve a wide variety of partial differential equations, including the Euler and Navier-Stokes equations for compressible or incompressible flow, the shallow water equations for flow in the ocean, rivers and channels and the porous media equation for groundwater flow in underground aquifers. The primary goal of these codes has been to simulate the flow for a particular geometry and/or parameter distribution. The codes are used in conjunction with scale models and field measurements for evaluation and analysis of a particular design. Then, the design or the parameters are changed based on the knowledge and experience of the researcher, and a new simulation is run. This trial-and-error approach to design is inefficient and does not take full advantage of the computational tools. By changing the simulation code into a gradient-based optimization code, the results from the flow simulation are used to direct the search and can efficiently locate much improved designs. For instance, in the inverse design of airfoils, the designer determines the desired pressure distribution based on certain design criteria and uses an Euler or Navier-Stokes code to find an airfoil whose pressure distribution closely matches the target distribution. Thus, the objective function to be minimized could be a least squares function between the target pressures and the actual pressures for a particular design, and the design parameters would be the shape of the airfoil and the flight conditions, such as Mach number and angle of attack. By varying these parameters, the designer hopes to minimize the objective function. For notation, $\vec{\beta}$ is the vector of design parameters, $F(\vec{\beta})$ is the objective function, $\frac{\partial F}{\partial \beta_i}$ is the design space derivative of $F$ with respect to $\beta_i$ and $\nabla_{\vec{\beta}} F$ is the design space gradient of $F$ with respect to the vector of design variables $\vec{\beta}$. These simulation codes can be used to estimate design space derivatives by using one-sided finite differences \begin{equation} \frac{\partial F(\vec{\beta})}{\partial \beta_i} \approx \frac{F(\vec{\beta}+e_i \Delta \beta) - F(\vec{\beta})}{\Delta \beta} \end{equation} or central differences \begin{equation} \frac{\partial F(\vec{\beta})}{\partial \beta_i} \approx \frac{F(\vec{\beta}+e_i \Delta \beta) - F(\vec{\beta}-e_i \Delta \beta)}{2\Delta \beta} \end{equation} where $e_i$ is a unit vector in the $i$th direction. These finite difference formulae for estimating the design space gradient are easy to implement. However, for one-sided finite differences, an $N$ additional steady-state simulations are required where $N$ is the number of design parameters, and when central differences are used, an additional $2 N$ steady-state simulations are needed. Since a high-fidelity flow solver is being used to calculate the steady-state flow variables, the computational cost of evaluating the objective function is quite large. Thus, the finite difference method is computationally expensive. Furthermore, the accuracy of the method is highly dependent on the size of the perturbation $\Delta \beta$. If the perturbation is too large, the non-linearity of the objective function will dominate; if it is too small, round-off error reduces the accuracy of the derivative. Finally, if the solutions are not strongly converged, then the accuracy of the function values will be limited. In this case, the function values being used in these formulae must be taken from identical locations in the solution process to have any chance of being accurate. Other methods being investigated to generate numerically exact derivatives include automatic differentiation such as ADIFOR and the complex Taylor's series expansion method. These two methods are relatively easy to implement but are approximately as computationally expensive as finite differences, although current research efforts are being investigated to reduce their computational cost for steady-state problems. Table 1 shows a comparison of these methods in regards to their computational cost, accuracy of derivative and ease of implementation. \begin{center} \begin{tabular}{|c|c|c|c|} \hline Derivative & Cost of & Accuracy of & Ease of \\ Method & Comput. & Derivative & Implement. \\ \hline \hline Finite Difference & Expensive & Moderate & Very Easy \\ \hline Automatic Differentiation & Expensive & Num. Exact & Easy \\ \hline Complex Taylor's Series & & & \\ Expansion Method & Expensive & Num. Exact & Easy \\ \hline Continuous Sensitivity Analysis & & & \\ for Explicit Codes & Cheap & Moderate & Very Difficult \\ \hline Discrete Sensitivity Analysis & & & \\ for Implicit Codes & Cheap & Highly & Moderate \\ \hline Discrete Sensitivity Analysis & & & \\ for Explicit Codes & Cheap & Moderate & Difficult \\ \hline \end{tabular} Table 1. Costs and Benefits of Various Methods \\ to Generate Design Space Derivatives. \end{center} Sensitivity analysis can be used to reduce the computational cost of estimating the design space gradient, by taking advantage of derivative information obtained via differentiation of the governing system of equations. Sensitivity analysis can be divided into two different approaches - the continuous approach and the discrete approach. In the continuous approach, the system of governing differential equations is differentiated to form a separate set of continuous adjoint equations. In the discrete approach, the system of governing equations is discretized first, creating a system of discretized equations $W(Q)=0$. These discretized equations are differentiated to form a system of discrete adjoint equations. For an overview of sensitivity analysis as applied to complex aerodynamic applications, see Newman, etal \cite{n2}. Figure 1 gives a breakdown of the various formulations within sensitivity analysis. The adjoint formulations of both continuous and discrete sensitivity analysis are typically used due to the need to calculate the adjoint variable once per objective function regardless of the number of design variables. Each formulation can be implemented within an explicit or an implicit code; however, traditionally, discrete sensitivity analysis has been used for implicit codes while continuous sensitivity analysis has been used for explicit codes. \begin{center} \epsfig{file=sensitivity.eps, height=1.4in} Figure 1. Breakdown of Formulations within Sensitivity Analysis. \end{center} For explicit codes, the continuous adjoint equations can be discretized and solved in any consistent fashion; however, Shubin and Frank \cite{s1} showed that to achieve the best agreement between the derivatives produced via the continuous approach with the finite difference derivatives, the same discretization method should be employed to solve the continuous adjoint equations as was used to solve the continuous governing equations. Reuther \cite{r1}, in his dissertation research, demonstrated these techniques for the two-dimensional Euler equations as applied to airfoil design. Once the adjoint equations are derived and successfully implemented, they are driven to steady-state, which takes approximately twice as long as a flow simulation, according to Reuther. Since the adjoint equations must be solved only once for each objective function regardless of the number of design variables and since the number of objective functions is often quite small in comparison to the number of design variables for single discipline designs, the computational savings of using the adjoint formulation as opposed to finite differences can be quite large. More recently, Nadarajah and Jameson \cite{n1} showed that the discrete approach of sensitivity analysis could be applied to explicit codes by differentiating the discretized equations solved via the explicit code. However, in his derivations, Nadarajah froze certain terms to reduce the complexity of the differentiated equations, which resulted in a loss of accuracy between the finite difference derivative and the derivative calculated via his approach. Furthermore, the computational cost of solving the discrete adjoint equations via an explicit method is quite similar to solving the discrete governing equations. Thus, for explicit codes, both the continuous and discrete approaches of sensitivity analysis can be used, although the complexity of the derivations, the computational cost and the inaccuracy between the derivative results are limitations for these codes. For implicit codes, both the continuous and discrete approaches can also be employed. However, because implicit codes already calculate the Jacobian matrix of the discretized governing equations with respect to the flow variables, the discrete approach is quite easy to implement; whereas, for the continuous approach, the Jacobian of the discrete adjoint equations would be needed. The focus of the research presented herein is to develop a method for easily generating the Jacobian of the discretized governing equations as implemented within either an implicit or an explicit code, so that discrete sensitivity analysis can be applied uniformly to both solution methodologies. One final distinction between the continuous and discrete approaches of sensitivity analysis should be emphasized. The derivatives generated by the continuous approach are not ``discretely adjoint'' to the discretized governing equations. As a result, the derivatives produced by the continuous approach are generally not as accurate as those produced by the discrete approach. The error is a result of the order of discretization and can be summarized by Figure 2. In discrete sensitivity analysis, the discretization error associated with the adjoint variable is consistent with the discretization error generated by solving the discretized governing equations. However, for continuous sensitivity analysis, the discretization error results from solving the discretized adjoint equations, which does not have a direct relationship to the discretization error coming from the governing equations. For discrete simulations, the adjoint variable $\lambda$ generated by discrete sensitivity analysis will be discretely adjoint to the discretized governing equations. Since the derivatives are consistent with the discrete solution, these derivatives will be in better agreement with derivatives produced via finite differences or the complex Taylor's series expansion method than those derivatives produced by the continuous approach. However, as the mesh size is reduced, the discretization error will tend towards zero, and the derivatives produced by the two approaches will agree. \begin{center} \epsfig{file=discret.eps, height=2.7in} Figure 2: Effects of Discretization Error on Accuracy of Adjoint Variable. \end{center} For discrete sensitivity analysis, one of the primary factors affecting the accuracy of the design space derivatives is the accuracy of the Jacobian matrix. To obtain each term in this matrix by hand-differentiation is also quite challenging and error-prone. However, the complex Taylor's series expansion method can be used to overcome the difficulties associated with determining the Jacobian matrix, by using complex arithmetic to generate highly accurate derivatives without the need for hand-differentiation. Since the steady-state discretized equations are the same for explicit and implicit codes, the complex Taylor's series expansion method can be used to generate highly accurate Jacobian matrix for both explicit and implicit codes. Using this method, highly accurate Jacobian matrices can be generated for explicit and implicit codes, and discrete sensitivity analysis can be uniformly applied to both solution methodologies. Sensitivity analysis is presented in detail in the next section. The complex Taylor's series expansion method is discussed in section 3. In section 4, transformation of the Euler simulation code is explained, and in section 5, the resulting optimization code is discussed and is applied to an inverse design problem. \section{Sensitivity Analysis} The objective function $F$ is explicitly a function of the flow variables $Q$ and possibly of the grid $\chi$ or the design variables $\vec{\beta}$, so $F(Q(\vec{\beta}),\chi(\vec{\beta}), \vec{\beta})$. Thus, the design space derivative of $F$ with respect to the design variable $\beta_i$ can be expressed as \begin{equation} \frac{dF}{d\beta_i} = \frac{\partial F}{\partial Q} \frac{\partial Q}{\partial \beta_i} + \frac{\partial F}{\partial \chi} \frac{\partial \chi}{\partial \beta_i} + \frac{\partial F}{\partial \beta_i} \end{equation} Since $F$ is explicitly dependent on $Q$, $\chi$ and $\beta$, the terms $\frac{\partial F}{\partial Q}$, $\frac{\partial F}{\partial \chi}$ and $\frac{\partial F}{\partial \beta_i}$ can be calculated easily by hand. The term $\frac{\partial \chi}{\partial \beta_i}$ can be estimated via finite differencing the results of the grid generation code. Unfortunately, the term $\frac{\partial Q}{\partial \beta_i}$ can not be estimated via finite differences without an additional steady-state simulation. ($\frac{dF}{d\beta_i}$ is actually a partial derivative because $F$ depends on several design variables; however, due to the limitations of notation, it is written as the total derivative, being the sum of the explicit and implicit dependencies of $F$ on the design variable $\beta_i$.) The system of governing equations can be expressed as $W(Q(\vec{\beta}),\chi({\vec{\beta}}), \vec{\beta})=0$. Thus, \begin{equation} \frac{dW}{d\beta_i} = \frac{\partial W}{\partial Q} \frac{\partial Q}{\partial \beta_i} + \frac{\partial W}{\partial \chi} \frac{\partial \chi}{\partial \beta_i} + \frac{\partial W}{\partial \beta_i}=0 \end{equation} Again, since $W$ is explicitly dependent on $Q$, $\chi$ and $\beta$, the associated terms $\frac{\partial W}{\partial Q}$, $\frac{\partial W}{\partial \chi}$ and $\frac{\partial W}{\partial \beta_i}$ can be calculated efficiently. The term $\frac{\partial \chi}{\partial \beta_i}$ is already known via finite differences. Hence, the only unknown term in this equation is $\frac{\partial Q}{\partial \beta_i}$. By multiplying equation (4) by the adjoint variable $\lambda$ and adding to equation (3), we get \begin{equation} \frac{\partial F}{\partial \beta_i} = \left(\frac{\partial F}{\partial Q} +\lambda^T \frac{\partial W}{\partial Q} \right) \frac{\partial Q}{\partial \beta_i} + \left(\frac{\partial F}{\partial \chi} + \lambda^T \frac{\partial W}{\partial \chi} \right) \frac{\partial \chi}{\partial \beta_i} + \left(\frac{\partial F}{\partial \beta_i} + \lambda^T \frac{\partial W}{\partial \beta_i}\right) \end{equation} By choosing $\lambda$ such that \begin{equation} \frac{\partial F}{\partial Q} +\lambda^T \frac{\partial W}{\partial Q} = 0 \end{equation} or \begin{equation} \frac{\partial W}{\partial Q}^T\lambda = -\frac{\partial F}{\partial Q}^T \end{equation} the need to calculate the term $ \frac{\partial Q}{\partial \beta_i}$ is removed, and the design space derivative can be calculated without the need of any additional steady-state simulations. Equation (7) does not depend on the design variables; hence $\lambda$ is the same for all the design variables, and equation (7) must only be solved once, regardless of the number of design variables. However, equation (7) is dependent on the objective function $F$; thus, this equation must be solved for each function for which a derivative is needed. In the discrete approach, the terms $F$ and $W$ represent the discretized objective function and governing equations, whereas in the continuous approach, these terms are the continuous objective function and governing differential equations. For the objective function used herein, $F$ is only explicitly dependent on the flow variables $Q$; hence, after solving for $\lambda$, the design space derivative is \begin{equation} \frac{\partial F}{\partial \beta_i} = \lambda^T \left(\frac{\partial W}{\partial \chi} \frac{\partial \chi}{\partial \beta_i} + \frac{\partial W}{\partial \beta_i} \right) = \lambda^T \frac{dW}{d\beta_i}\bigg{|}_{\mbox{Q fixed}} \end{equation} To study the differences in implementation between the two approaches, a review of the differences between implicit and explicit simulation codes is appropriate. For explicit codes, the flow variables are updated via an equation similar to \begin{equation} Q^{n+1} = Q^n + f_{explicit}(Q^n,dt^n) \end{equation} whereas for implicit codes, the flow variables are updated via \begin{equation} Q^{n+1} = Q^n + f_{implicit}(Q^{n+1},Q^n,dt^n) \end{equation} where $Q^{n+1}$ is the vector of flow variables at time level $n+1$, $Q^n$ is the vector of flow variables at time level $n$, which are known, and $f_{explicit}$ and $f_{implicit}$ are functions that define the discretized spatial derivative. If higher order discretizations of the temporal derivatives are used or if fractional step algorithms are used, equations (9) and (10) are modified to account for these changes. For explicit codes, the function $f_{explicit}$ depends only on known values, whereas for implicit codes, the unknowns $Q^{n+1}$ are included within the update function $f_{implicit}$. For implicit codes, equation (10) can be recast as \begin{equation} W(Q^{n+1},Q^n,dt^n) = Q^{n+1} - Q^n - f_{implicit}(Q^{n+1},Q^n,dt^n) = 0 \end{equation} and can be solved iteratively via \begin{equation} \frac{\partial W}{\partial Q^{n+1}} \Delta Q^{m} = - W(Q^{n,m},Q^n,dt^n) \end{equation} where $Q^{n,m+1} = Q^{n,m} + \Delta Q^m$. This equation is solved iteratively until the residual vector $W(Q^{n,m},Q^n,dt^n)$ is sufficiently small, at which point $Q^{n+1} = Q^{n,m}$. The matrix $\frac{\partial W}{\partial Q^{n+1}}$ is called the Jacobian matrix of the residual vector with respect to the flow variables. Since equation (12) is solved iteratively in delta form, the Jacobian matrix does not need to be exact for the flow solver as long as $W$ is driven to zero. At steady-state, $Q^{n+1} = Q^n$, so $f_{explicit}=0$ and $f_{implicit}=0$. Furthermore, after using the identity that $Q^{n+1} = Q^n$, the Jacobian matrices can be expressed as \begin{equation} \frac{\partial W}{\partial Q^{n+1}} = - \frac{\partial f_{implicit}}{\partial Q^{n+1}} \end{equation} and \begin{equation} \frac{\partial W}{\partial Q^n} = -\frac{\partial f_{explicit}}{\partial Q^n} \end{equation} Typically, $\frac{\partial f_{explicit}}{\partial Q^n}$ is not calculated in explicit codes because it is not needed for the flow solver. However, if this matrix were available, discrete sensitivity analysis could be applied to explicit codes in the same way as it is applied to implicit codes. Since the Jacobian matrix is used in implicit codes, it is relatively straight-forward to convert implicit flow solvers into optimization codes using discrete sensitivity analysis, although the accuracy of the resulting design space derivatives will be limited by the accuracy of the Jacobian matrix. By using the complex Taylor's series expansion method, the exact Jacobian can be generated regardless of the complexity of the turbulence models or spatial discretization methods, and highly accurate design space derivatives can be achieved. To avoid the difficulties associated with the continuous approach of sensitivity analysis, the complex Taylor's series expansion method can be applied to explicit codes at steady-state to calculate the Jacobian matrix $\frac{\partial f_{explicit}}{\partial Q^n}$. Hence, discrete sensitivity analysis can be used in conjunction with explicit codes to generate the design space derivatives efficiently without the need for deriving and implementing the adjoint equations as well as driving these equations to steady-state. \section{Complex Taylor's Series Expansion Method} The complex Taylor's series expansion method has been recently used by Squire and Trapp \cite{s2} to calculate the derivatives of real-valued functions and by Newman etal \cite{n3} to estimate derivatives for an aero-structural design problem. This simple method can be applied to any numerical code that yields real-valued functional information and is based on the Taylor's series expansion of $F(x+i\Delta x)$ or \begin{equation} F(x+i\Delta x) = F(x) + i F'(x) \Delta x - \frac{F''(x) \Delta x^2}{2!} - i \frac{F^{(3)}(x) \Delta x^3}{3!} + O(\Delta x^4) \end{equation} Since $F(x)$ is real valued, the Taylor's series alternates between real and imaginary terms and can be decomposed as \begin{eqnarray} \lefteqn{F(x+i\Delta x) }\\ &=& \left(F(x) - \frac{F''(x) \Delta x^2}{2!} + O(\Delta x^4), F'(x) \Delta x - \frac{F^{(3)}(x) \Delta x^3}{3!} + O(\Delta x^5)\right)\nonumber \end{eqnarray} Hence, \begin{equation} Im(F(x+i\Delta x)) = F'(x) \Delta x - \frac{F^{(3)}(x) \Delta x^3}{3!} + O(\Delta x^5) \end{equation} or \begin{equation} F'(x) = \frac{Im(F(x+i\Delta x))}{\Delta x} + \frac{F^{(3)}(x) \Delta x^2}{3!} + O(\Delta x^4) \end{equation} This method is second order accurate as is the central finite difference equations, but as there is no subtraction error, there are no round-off errors, and $\Delta x$ can be as small as the computer will allow. Hence, choosing $\Delta x$ such that the higher order terms are smaller than machine precision, we have \begin{equation} F'(x) = \frac{Im(F(x+i\Delta x))}{\Delta x} \end{equation} Thus, numerically exact Jacobians can be estimated via the application of the complex Taylor's series expansion method to the residual vector $W(Q^n,\Delta t^n)$. For a two-dimensional, structured grid, using a first-order spatial discretization, the equations at an interior node $(i,j)$ will be dependent on the values at $(i,j)$, $(i-1,j)$, $(i+1,j)$, $(i,j-1)$ and $(i,j+1)$, as shown in Figure 3. For the two-dimensional Euler equations, there are four differential equations to be solved and four flow variables to be determined at each node. Hence, there are four discretized equations associated with each node, and each equation can be dependent on as many as twenty different flow variables. \begin{center} \epsfig{file=gdepend.eps, height=1.50in} Figure 3. Grid Dependencies for Two-Dimensional, First-Order Scheme. \end{center} Repeated application of the complex Taylor's series expansion method to each discretized equation represented in the residual vector $W$ will generate the values of the twenty derivative terms associated with a row in the Jacobian matrix, without any need for hand-differentiation. Hence, any turbulence model or spatial discretization scheme can be used, regardless of its complexity, as long as the derivatives exist. Furthermore, for unstructured grids, or for higher-order spatial discretizations, the only added complexity is the connectivity stencil. In other words, to generate the exact Jacobian matrix, the dependencies of the vector of discretized equations with the flow variables must be known, but the knowledge of the particulars of these dependencies is unnecessary as the complex Taylor's series expansion method will handle these complexities. \section{Application to Euler Code} The explicit flow simulation code that has been transformed into a design optimization code solves the two-dimensional Euler equations for flow around the NACA64A006 airfoil, using a structured C-grid. A wide variety of solution methods are available within the code, including implicit methods and fractional step methods; however, these subroutines have been deleted, leaving only the explicit implementation of the first order Flux-Vector-Split method of Steger and Warming \cite{s3}. Since knowledge of the implementation of this method is not used in the transformation of the simulation code, a discussion of this method will not be presented here - the key element is that this method is explicit, and the stencil is the same as shown in Figure 3. \begin{center} \epsfig{file=derivatives.eps, height=4.0in} Figure 4. Flow-Chart for Algorithm to Estimate Derivatives \end{center} The algorithm for using discrete sensitivity analysis to estimate the design space derivatives is presented in Figure 4. The grid is determined by perturbing the shape of the NACA64A006 airfoil and propagating the perturbation into the two-dimensional grid. The shape of the airfoil was controlled by a B-spline curve with 14 control points. Ten of these control points were determined by the ten design variables, with the additional control points controlling the smoothness at the ends of the curve. Given a set of design variables, the B-spline curve was determined for each node on the surface of the airfoil. The value of the curve at each node determined the displacement in the y-direction that the new airfoil was from the NACA64A006 airfoil. Hence, if the B-spline curve's value at node 98 was +0.0012, then the y-location of the surface at node 98 for the new airfoil was the value for the NACA airfoil plus 0.0012. Once the displacement values were determined for the surface of the airfoil, these values were propagated into the grid by solving Laplace's equation for the NACA grid where the boundary conditions were the displacement values. Once the values were propagated into the interior of the grid, the y-values for the new grid were determined by adding the value at a node to its original value. The complex Taylor's series expansion method is used to generate the Jacobian matrix $\frac{\partial W}{\partial Q}$. The vector $\frac{\partial F}{\partial Q}$ is calculated via hand differentiation of the function $F$ which is explicitly a function of the flow variables $Q$. The vector $\frac{dW}{d\beta_i}\big{|}_{\mbox{Q fixed}}$ is calculated via central differences or \begin{equation} \frac{dW}{d\beta_i} = \frac{W(Q(\vec{\beta}),\chi(\vec{\beta} +e_i \Delta \beta_i),\vec{\beta}+e_i \Delta \beta_i) - W(Q(\vec{\beta}),\chi(\vec{\beta}-e_i \Delta \beta_i),\vec{\beta}-e_i \Delta \beta_i)}{2\Delta \beta_i} \end{equation} This calculation requires the formation of the grids associated with $\chi(\vec{\beta}+e_i \Delta \beta_i)$ and $\chi(\vec{\beta}-e_i \Delta \beta_i)$. \begin{center} \epsfig{file=cgrid.eps, height=2.3in} Figure 5. Typical C-grid around an airfoil. \end{center} To generate the Jacobian matrix, the nodal dependencies of the residual vector must be analyzed. The NACA64A006 grid is a C-grid, similar to the grid shown in Figure 5. For the interior nodes, the nodal connectivity stencil is the same as shown in Figure 3, resulting in a banded, block structure for the Jacobian matrix. For the far-field boundary and impermeable boundary, the connectivity stencil is a subset of the stencil for the interior nodes, so the banded, block structure of the Jacobian matrix is unchanged at these boundaries. However, for the re-entrant boundary, the nodal connectivity reaches to the other side of the grid, destroying the banded, block structure of the Jacobian. Thus, to solve the matrix equation (7), an iterative solver is employed. Unfortunately, this solver does not converge rapidly, and the solution is not highly converged. Finally, to update the design variables, the BFGS update method has been employed. This method gradually builds an approximation to the Hessian matrix by updating the approximate Hessian matrix with gradient and design variable information. One additional function evaluation is performed in the search direction to determine a more appropriate step size. See Gill, Murray and Wright \cite{g1} for more details on this method. \section{Results} The design problem is an inverse design where the pressure distribution for the target design is specified. The design variables for the target design are (-0.001, -0.003, 0.001, 0.005, 0.01, 0.01, 0.02, 0.03, 0.02, 0.01). The steady-state solution for this design is approximated, and the pressure distribution along the surface of the airfoil are stored in $Cp^{target}$. The mach number is 0.8700, which results in transonic flow across the airfoil, and the angle of attack for the simulation is +3.2 degrees. The objective function is \begin{equation} F(\vec{\beta}^n) = \sum_{i=itl+1}^{itu} \left( Cp_i(\vec{\beta}^n) - Cp_i^{target}\right)^2 \end{equation} where $itl$ and $itu$ specify the initial and final nodes for the surface of the airfoil in the C grid. The initial set of design variables is a vector of 0.0's, which represents the unmodified NACA airfoil. To analyze the accuracy of the design space derivatives, the derivatives for this set of design variables are calculated via the complex Taylor's series expansion method applied to the objective function. Hence, these derivatives are numerically exact, although the computational cost is quite large. These exact derivatives are compared with the derivatives generated by discrete sensitivity analysis and by one-sided finite differences and are presented in Table 2. \begin{center} \begin{tabular}{|c|c|c|c|} \hline Design & Complex Taylor's & Finite & Adjoint Form. of Discrete \\ Variables & Series Expansion & Differences & Sensitivity Analysis \\ \hline \hline 1 & -58.918524793680 & -58.918108 & -58.925014473517 \\ \hline 2 & 2.1553606507097 & 2.155659 & 2.1527299069564 \\ \hline 3 & 9.3080663602748 & 9.308617 & 9.3059170047628 \\ \hline 4 & -4.4977859312973 & -4.496398 & -4.5001198714489 \\ \hline 5 & 3.1437212541549 & 3.146781 & 3.1419623720023 \\ \hline 6 & 587.33976898291 & 587.343377 & 587.35094264658 \\ \hline 7 & -457.04698505623 & -457.044045 & -457.04649298139 \\ \hline 8 & -431.92433510382 & -431.923165 & -431.98154907862 \\ \hline 9 & -60.885888964239 & -60.885054 & -60.878678137755 \\ \hline 10 & -128.31651602328 & -128.315028 & -128.20253621961 \\ \hline Comp. Cost & 25870.7 sec. & 4136.0 sec. & 159.9 sec. \\ \hline \end{tabular} Table 2. Comparison of Design Space Derivatives. \end{center} Due to the errors introduced into the design space derivatives by the inexact solution of equation (7), the derivatives produced by discrete sensitivity analysis are not quite as accurate as those produced by finite differences, although they agree with the exact derivatives to at least three significant digits. However, the computational cost of such calculations, given in the bottom row, show the tremendous savings of discrete sensitivity analysis over the other two methods. In Table 3, the computational cost for each component in a design iteration using discrete sensitivity analysis is given. Since there are 10 design variables, the cost of estimating the design space gradient via finite differences is approximately 4100 seconds, which is substantially more than the cost associated with discrete sensitivity analysis. Furthermore, as the number of design variables increases, the computational cost of using finite differences scales with the number of design variables, whereas the cost of discrete sensitivity analysis grows marginally, because the majority of the cost lies in determining the adjoint vector, which must only be calculated once, regardless of the number of design variables. \begin{center} \begin{tabular}{|c|c|} \hline Component & Cost \\ \hline \hline Generating 1 Steady-State Solution & 413.8 sec. \\ \hline Generating $\nabla_{\vec{\beta}} F$ Via Discrete Sensitivity Analysis & 159.9 sec. \\ \hline Updating Design Variables (includes an additional solution) & 413.6 sec. \\\hline Total Cost of Design Iteration Via Discrete Sensitivity Analysis & 987.3 sec.\\ \hline \end{tabular} Table 3. Computational Costs for Components of Design Process. \end{center} Figure 6 shows the convergence history of the objective function values for the BFGS update method, using one additional function evaluation. The method is slow at improving the design variables for the first fifteen design iterations; however, after the twentieth design iteration, the objective function decreases greatly, stabilizing near $2.3E-10$. By using the BFGS update method, superlinear convergence results have been obtained, which reduces the number of design iterations and hence the total computational cost of the design process. The design variables corresponding to iteration \#23 are in agreement with the target set of design variables to at least seven decimal places. Hence, one can conclude that the design space derivatives are accurate enough to drive the design variables quite close to the target. \begin{center} \epsfig{file=euler.eps, height=2.4in} Figure 6. Optimization History for BFGS Update Method. \end{center} \section{Conclusions} Discrete sensitivity analysis in conjunction with the complex Taylor's series expansion method has been used to transform an explicit, two-dimensional Euler solver into a design optimization code. As has been derived in many other papers, the formulation of discrete sensitivity analysis presented herein requires the Jacobian of the discretized system of equations, which is not generated or used in the explicit code. The complex Taylor's series expansion method is used to generate the exact Jacobian matrix, once the steady-state solution has been obtained. Thus, there is no need to derive or solve the continuous adjoint equations as is typically done when applying the continuous approach of sensitivity analysis to explicit codes. The only information needed about the discretized equations is the connectivity stencil, showing the dependence of the discretized equations on the nodes within the grid. As a result, the complex Taylor's series expansion method should be easily applicable to complicated turbulence models and higher-order schemes. This method has been applied to an explicit code that solves the two-dimensional Euler equations for flow around an airfoil. The resulting design space derivatives have been compared with the numerically exact derivatives, agreeing to at least three significant digits. These design space derivatives are also accurate enough to drive the design variables towards the target set of design variables, with the converged set of design variables agreeing with the target set to at least seven decimal places. Furthermore, the computational cost of estimating these derivatives for the discrete sensitivity analysis is substantially less than for finite differences. The explicit code is only first-order in space and does not use fractional steps. More research is necessary to address the difficulties of using higher-order spatial discretizations, which increase the connectivity stencil and hence the computational cost of solving equation (7). Furthermore, more study is necessary to apply discrete sensitivity analysis to fractional step methods, such as a predictor-corrector method. \paragraph{Acknowledgments.} The author would like to thank Dr. Mark Janus for providing the explicit Euler code used in this study. Furthermore, the author would like to thank Dr. David Huddleston for introducing the author to the concept of sensitivity analysis and Dr. David Whitfield for introducing the author to the complex Taylor's series expansion method. \begin{thebibliography}{0} \bibitem{g1} Gill, P. E., Murray, W., and Wright, M. H., {\it{Practical Optimization}}, Academic Press Ltd., San Diego, 1981. \bibitem{n1} Nadarajah, S., and Jameson, A., ``A Comparison of the Continuous and Discrete Adjoint Approach to Automatic Aerodynamic Optimization'', AIAA Paper 2000-0667. \bibitem{n2} Newman, J. C. III, Taylor, A. C. III, Barnwell, R. W., Newman, P. A., and Hou, G. J. W., ``Overview of Sensitivity Analysis and Shape Optimization for Complex Aerodynamic Configurations'', {\it{AIAA J. of Aircraft}}, Vol. 36, No.1, Jan.-Feb., 1999., pp. 87-96. \bibitem{n3} Newman, J. C. III, Anderson, K. W., and Whitfield, D. L., ``Multidisciplinary Sensitivity Derivatives Using Complex Variables'', Mississippi State University Publication, MSSU-EIRS-ERC-98-08, July, 1998. \bibitem{r1} Reuther, J. J., ``Aerodynamic Shape Optimization Using Control Theory'', Dissertation, University of California, Davis, 1996. \bibitem{s1} Shubin, G. R., and Frank, P. D., ``A Comparison of Two Closely-Related Approaches to Aerodynamic Design Optimization'', Third International Conf. on Inverse Design Concepts and Opt. in Eng. Sciences (ICIDES-III), Washington, D. C., Oct. 23-25, 1991. \bibitem{s2} Squire, W., and Trapp, G., ``Using Complex Variables to Estimate Derivatives of Real Functions'', {\it{Soc. Ind. Appl. Math.}}, Vol. 40, No.1, pp.110-112, March, 1998. \bibitem{s3} Steger, J. L., and Warming, R. F., ``Flux Vector Splitting of the Inviscid Gas Dynamic Equations with Applications to Finite-Difference Methods'', {\it{J. Comp. Phys.}}, Vol. 40, 1981, pp. 263-293. \end{thebibliography} \noindent{\sc Clarence O. E. Burg}\\ Research Engineer, Computational Simulation and Design Center \\ Engineering Research Center \\ Mississippi State University, Mississippi State, MS, USA \\ email: burg@erc.msstate.edu \end{document}