AN AUTOMATED METHOD FOR SENSITIVITY ANALYSIS USING COMPLEX - - PDF document

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: nth 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) h + O(h)

• Central-difference:

f ′(x) ≈ f(x + h) − f(x − h) 2h + O(h2) Any finite-difference formula subject to subtractive cancellation.

– Joaquim Martins – January 12, 2000 – 3

Cauchy-Riemann Equations

• Extension of rational numbers to solve x2 = 2:

Define w = a + √ 2b where a, b are rational. Derivative of f(w): ∂f ∂w = ∂f ∂a = √ 2∂f ∂b.

• Extension of real numbers to solve x2 = −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 ∂y, ∂u ∂y = −∂v ∂x. Complex number really just one number.

– Joaquim Martins – January 12, 2000 – 4

Complex-Step Derivative Approximation

• From Cauchy-Riemann:

∂u ∂x = lim

h→0

v(x + i(y + h)) − v(x + iy) h . For a real functions of a real variable, y = 0, u(x) = f(x) and v(x) = 0, then, ∂f ∂x = lim

h→0

Im [f (x + ih)] h .

• From Taylor series expansion, step ih:

f(x+ih) = f(x)+ihf ′(x)−h2f ′′(x) 2! −ih3f ′′′(x) 3! +. . . ⇒ f ′(x) = Im [f(x + ih)] h + h2f ′′′(x) 3! + . . . No subtraction!

– Joaquim Martins – January 12, 2000 – 5

Simple Numerical Example

• Estimate derivative at x = 1.5 of the function,

f(x) = ex √ sin3x + cos3x

• ✥✧✦✩★✫✪✭✬✯✮✱✰✳✲✵✴✷✶✸✮✩✪

ε =

• f ′ − f ′

ref

• /
• f ′

ref

• – Joaquim Martins – January 12, 2000 –

6

Higher Derivative Approximations

• General form of Cauchy’s Integral,

f (n)(z) = n! 2πi

• Γ

f(ξ) (ξ − z)n+1dξ.

• Discretize using a trapezoidal-rule, integrate around

circle z + rei, f (n)(z) ≈ n! mr

m−1

• j=0

f

• z + rei2πj

m

• ei2πjn

m

.

• 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 ∂u ∂x = ∂v ∂y =

• −1

⇐ x < 0 +1 ⇐ x > 0 abs(x + iy) =

• −x − iy

⇐ x < 0 +x + iy ⇐ x ≥ 0 .

– Joaquim Martins – January 12, 2000 – 8

Overloaded Functions

abs ✘ abs(z) =

• −z

⇐ x < 0 +z ⇐ x ≥ 0 tan ✔ tan(z) = e−iz−eiz

e−iz+eiz

asin ✔ arcsin(z) = −i log[iz + (1 − z2)

1 2]

acos ✔ arccos(z) = −i log[z + (z2 − 1)

1 2]

atan ✔ arctan(z) = i

2 log

• 1−iz

1+iz

• sinh

✔ sinh(z) = ez−e−z

2

cosh ✔ cosh(z) = ez+e−z

2

tanh ✔ tanh(z) = ez−e−z

ez+e−z

dim ✘ dim(z1, z2) =

• z1 − z2

⇐ x1 > x2 ⇐ x1 ≤ x2 sign ✘ sign(z1, z2) =

• +|x1|

⇐ x2 ≥ 0 −|x1| ⇐ x2 < 0 max ✘ max(z1, z2) =

• z1

⇐ x1 ≥ x2 z2 ⇐ x1 < x2 min ✘ min(z1, z2) =

• z1

⇐ x1 ≤ x2 z2 ⇐ x1 > x2

– 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
• f 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

• verload

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.049760045804646 ADIFOR –39.049760045809059 Analytic –39.049760045805281 FD –39.049724352820375

• 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

RAE 2822 Airfoil α = 3.0o, M∞ = 0.70

1.20 0.80 0.40 0.00

• 0.40
• 0.80
• 1.20
• 1.60
• 2.00

Cp

+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + ++++++++++++++ + + +++++++++++++++++++ + + + + + + + + +

• 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

∂CD/∂M∞

10

−30

10

−20

10

−10

10

−8

10

−6

10

−4

10

−2

10

Step Size, h Normalized Error, ε

Complex−Step Finite−difference

– Joaquim Martins – January 12, 2000 – 15

Convergence of Sensitivity Estimates

50 100 150 10

−10

10

−5

10 10

5

10

10

Number of Iterations Residual

Complex−Step Finite−difference Average density residual

• Both estimates converge at same rate as solver.

– Joaquim Martins – January 12, 2000 – 16

Aerodynamic Sensitivity Results

Sensitivity Complex FD ∂CL/∂α 12.37054751691092 12.370726665267277 ∂CD/∂α 0.8602380269441042 0.86024234610682058 ∂CM/∂α –0.5026301652982372 –0.5026313670830616 ∂CL/∂M∞ 3.2499722985150532 3.2499447577549727 ∂CD/∂M∞ 0.43978671505102005 0.43978998888818932 ∂CM/∂M∞ –0.99037388394690016 –0.9903747405504145

• Imaginary part of result converged to the same number of digits

real part.

– Joaquim Martins – January 12, 2000 –

Drag Coefficient Sensitivity to Hicks-Henne “Bump” Functions

10 20 30 40 50 1 2 3 4 5 6 7 8 9 10 11

Design Points Drag Coefficient Sensitivity

Complex−Step Finite−difference

Method Time Memory Complex 1.00 1.00 FD 0.31 0.55

– Joaquim Martins – January 12, 2000 – 18

Conclusions

• The theory behind the application of the complex-

step method to real-world numerical algorithms was introduced.

• The implementation of the complex-step method

was successfully automated.

• The sensitivity results given by this method have

been further validated and shown to be very accurate.

• Unlike finite-differencing, the complex-step method

is relatively insensitive to the step size.

• The complex-step method is now an attractive

alternative to ADIFOR.

– Joaquim Martins – January 12, 2000 – 19

Further Work

• Continue to work on Complexify.py.

Need feedback from users.

• Apply the method to other solvers.
• More detailed stability and convergence analysis.
• Compare with CFD adjoint results.

– Joaquim Martins – January 12, 2000 – 20