AN AUTOMATED METHOD FOR SENSITIVITY ANALYSIS USING COMPLEX VARIABLES Joaquim R. R. A. Martins Ilan M. Kroo Juan J. Alonso Department of Aeronautics and Astronautics Stanford University January 12, 2000 – Joaquim Martins – January 12, 2000 –
Outline • Background on complex-step and other sensitivity analysis methods • Theory – First derivative approximations – Higher derivative approximations • Implementation of the complex-step in algorithms – Fortran functions and operators – Automatic implementation • Results – Structural sensitivities – Aerodynamic sensitivities • Conclusions – Joaquim Martins – January 12, 2000 – 1
Sensitivity Analysis Methods • Finite-Difference Methods • ADIFOR • Analytic Methods • Complex-Step: – Lyness & Moler, 1967: n th derivative approximation by integration in the complex plane – Squire & Trapp, 1998: Simple formula for first derivative, f ′ ≈ Im [ f ( x + ih )] /h – Newman, Anderson & Whitfield, 1998: Aero-structural code – Anderson, Whitfield & Nielsen, 1999: 3D Navier-Stokes with turbulence model – Joaquim Martins – January 12, 2000 – 2
Finite-Difference Derivative Approximations From Taylor series expansion. • Forward-difference: f ′ ( x ) ≈ f ( x + h ) − f ( x ) + O ( h ) h • Central-difference: f ′ ( x ) ≈ f ( x + h ) − f ( x − h ) + O ( h 2 ) 2 h Any finite-difference formula subject to subtractive cancellation . – Joaquim Martins – January 12, 2000 – 3
Cauchy-Riemann Equations • Extension of rational numbers to solve x 2 = 2 : √ Define w = a + 2 b where a , b are rational. Derivative of f ( w ) : √ ∂w = ∂f ∂f 2 ∂f ∂a = ∂b. • Extension of real numbers to solve x 2 = − 1 : Define z = x + iy , where x , y are real. Derivative of f ( z ) : ∂f ∂z = ∂f ∂y = i∂f ∂x . Define f ( z ) = u ( z ) + iv ( z ) , ∂u ∂x = ∂v ∂u ∂y = − ∂v ∂y, ∂x. Complex number really just one number . – Joaquim Martins – January 12, 2000 – 4
Complex-Step Derivative Approximation • From Cauchy-Riemann: ∂u v ( x + i ( y + h )) − v ( x + iy ) ∂x = lim . h h → 0 For a real functions of a real variable, y = 0 , u ( x ) = f ( x ) and v ( x ) = 0 , then, ∂f Im [ f ( x + ih )] ∂x = lim . h h → 0 • From Taylor series expansion, step ih : f ( x + ih ) = f ( x )+ ihf ′ ( x ) − h 2 f ′′ ( x ) − ih 3 f ′′′ ( x ) + . . . 2! 3! ⇒ f ′ ( x ) = Im [ f ( x + ih )] + h 2 f ′′′ ( x ) + . . . h 3! No subtraction! – Joaquim Martins – January 12, 2000 – 5
✖ ✔✕ ✒ ✔ ✔ ✢ ✜ ✛ Simple Numerical Example • Estimate derivative at x = 1 . 5 of the function, e x f ( x ) = √ sin 3 x + cos 3 x ✥✧✦✩★✫✪✭✬✯✮✱✰✳✲✵✴✷✶✸✮✩✪ ✹✻✺✩✼✾✽❀✿✩✼❂❁✭❃❅❄❇❆❉❈❊❈❂❋✩✼●❋☞❍✱■❏❋ ❑❇▲✩▼❖◆◗P❂❘✩❙❯❚❅❱❇❲❉❳❊❳❂▲✩P❂▲✩▼✱❨❏▲ � ✔✤✣ ✗✘✚✙ ✑✓✒ �✂✁☎✄✝✆✞�✠✟☛✡☞✄✍✌✏✎ � f ′ − f ′ � / � � � � � f ′ ε = � ref ref – Joaquim Martins – January 12, 2000 – 6
Higher Derivative Approximations • General form of Cauchy’s Integral, f ( n ) ( z ) = n ! � f ( ξ ) ( ξ − z ) n +1 dξ. 2 πi Γ • Discretize using a trapezoidal-rule, integrate around circle z + re i , � z + re i 2 πj � m − 1 f m f ( n ) ( z ) ≈ n ! � . e i 2 πjn mr m j =0 • Fixed r , use more points for better approximation. Bound on error. • With finite-difference, have to decrease h for better accuracy. • First derivative formula can be derived from this... ...but it is the only one that does not involve subtraction. – Joaquim Martins – January 12, 2000 – 7
Complex Functions and Operators in Fortran • Relational operators – Used with if statements to direct the execution thread. – Complex algorithm must follow same thread. – Therefore, compare only the real parts. – Also, max , min , etc. • Arithmetic functions and operators: – Most of these have a mathematical standard definition that is analytic. – Some of them are implemented in Fortran. – Exception: abs � ∂x = ∂v ∂u − 1 ⇐ x < 0 ∂y = +1 ⇐ x > 0 � − x − iy ⇐ x < 0 abs ( x + iy ) = ⇐ x ≥ 0 . + x + iy – Joaquim Martins – January 12, 2000 – 8
Overloaded Functions � − z ⇐ x < 0 ✘ abs( z ) = abs + z ⇐ x ≥ 0 tan( z ) = e − iz − e iz ✔ tan e − iz + e iz 1 arcsin( z ) = − i log[ iz + (1 − z 2 ) ✔ 2 ] asin arccos( z ) = − i log[ z + ( z 2 − 1) 1 ✔ 2 ] acos � � 1 − iz arctan( z ) = i ✔ 2 log atan 1+ iz sinh( z ) = e z − e − z ✔ sinh 2 cosh( z ) = e z + e − z ✔ cosh 2 tanh( z ) = e z − e − z ✔ tanh e z + e − z � z 1 − z 2 ⇐ x 1 > x 2 ✘ dim( z 1 , z 2 ) = dim 0 ⇐ x 1 ≤ x 2 � + | x 1 | ⇐ x 2 ≥ 0 ✘ sign( z 1 , z 2 ) = sign −| x 1 | ⇐ x 2 < 0 � z 1 ⇐ x 1 ≥ x 2 ✘ max( z 1 , z 2 ) = max z 2 ⇐ x 1 < x 2 � z 1 ⇐ x 1 ≤ x 2 ✘ min( z 1 , z 2 ) = min z 2 ⇐ x 1 > x 2 – Joaquim Martins – January 12, 2000 – 9
Automatic Implementation • Cookbook procedure: – Substitute all real type variable declarations with complex declarations. – Define all functions and operators that are not defined for complex arguments. – A complex-step can then be added to the desired variable and the derivative can be estimated by f ′ ≈ Im[ f ( x + ih )] /h . • Fortran 77: write new subroutines, substitute some of the intrinsic function calls by the subroutine names, e.g. abs by c abs . But ... need to know variable types in original code. • Fortran 90: can overload intrinsic functions and operators, including comparison operators. Compiler knows variable types and chooses correct version of the function or operator. • Complexify.py : Python script processes code. complexify.f90 : overloaded definitions. – Joaquim Martins – January 12, 2000 – 10
Structural Sensitivities: Finite Element Solver • FESMEH, finite element solver used in MDO. • Triangular plates and truss elements. • Wing modeled with spars, ribs and skin. • Spars and ribs modeled with plates and trusses. – Joaquim Martins – January 12, 2000 – 11
✔ ✜ ✒ ✔ ✔✕ ✖ ✢ ✛ Sensitivity Estimates vs. Step Size • Sensitivity of truss stress w.r.t. truss area. ✥✧✦✩★✫✪✭✬✯✮✱✰✳✲✵✴✩✶✷✮✩✪ ✸✩✹✻✺✭✹✽✼✿✾✭❀❂❁✧✹❄❃❅❃✿✾✩❆✿✾✩✺✱❇❈✾ � ✔✤✣ ✗✘✚✙ ✑✓✒ �✂✁☎✄✝✆✞�✠✟☛✡☞✄✍✌✏✎ • Finite-difference: reasonable estimate, if you’re lucky... • Complex-step: practically insensitive to step size. – Joaquim Martins – January 12, 2000 – 12
Computational Accuracy and Cost Method Sample Sensitivity Complex –39.04976004580 4646 ADIFOR –39.04976004580 9059 Analytic –39.04976004580 5281 FD –39.0497 24352820375 • Finite-difference is worst. • Other methods achieve the solver’s precision. Method Time Memory Complex 1.00 1.00 ADIFOR 2.33 8.09 Analytic 0.58 2.42 FD 0.88 0.72 • Analytic: best, but not easy to implement. • ADIFOR: costly. • Complex-step: good compromise. • Caveat: ratios depend on problem. – Joaquim Martins – January 12, 2000 – 13
Aerodynamic Sensitivities: FLO82 -2.00 + + + + + + + + + + + + + + + + ++++++++++++++ + + + + + + + + + + + + + + + + -1.60 + + + + + -1.20 + + + -0.80 + +++++++++++++++++++ + Cp + -0.40 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.00 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.40 + + + + + + + + + 0.80 + + + + + + + RAE 2822 Airfoil 1.20 α = 3 . 0 o , M ∞ = 0 . 70 • 2D, cell-centered, finite volume Euler solver. • Multigrid. • Implicit residual smoothing. • Artificial dissipation options: JST, CUSP, ECUSP and HCUSP. • Same numerical schemes as FLO87 and FLO107. – Joaquim Martins – January 12, 2000 – 14
Sensitivity Estimates vs. Step Size ∂C D /∂M ∞ 0 10 Complex−Step −2 Normalized Error, ε 10 Finite−difference −4 10 −6 10 −8 10 −10 −20 −30 10 10 10 Step Size, h – Joaquim Martins – January 12, 2000 – 15
Convergence of Sensitivity Estimates 10 10 Complex−Step Finite−difference 5 10 Average density residual Residual 0 10 −5 10 −10 10 0 50 100 150 Number of Iterations • Both estimates converge at same rate as solver. – Joaquim Martins – January 12, 2000 – 16
Aerodynamic Sensitivity Results Sensitivity Complex FD ∂C L /∂α 12.370547 51691092 12.370 726665267277 ∂C D /∂α 0.8602380 269441042 0.8602 4234610682058 ∂C M /∂α –0.502630 1652982372 –0.50263 13670830616 ∂C L /∂M ∞ 3.249972 2985150532 3.2499 447577549727 ∂C D /∂M ∞ 0.439786 71505102005 0.43978 998888818932 ∂C M /∂M ∞ –0.990373 88394690016 –0.99037 47405504145 • Imaginary part of result converged to the same number of digits real part. – Joaquim Martins – January 12, 2000 –
Recommend
More recommend
