THE CONNECTION BETWEEN THE COMPLEX-STEP DERIVATIVE APPROXIMATION AND ALGORITHMIC DIFFERENTIATION Joaquim R. R. A. Martins Peter Sturdza Juan J. Alonso Department of Aeronautics and Astronautics Stanford University — Joaquim R. R. A. Martins — 1
Outline • Background on sensitivity analysis • Complex-step derivative approximation basics • Improvements and the connection to algorithmic differentiation • Implementations of the complex-step method (Fortran, C/C++) • Results – 2D Flow Solver (Fortran) – 3D High-Fidelity Aero-Structural Analysis (Fortran) – Supersonic Viscous/Inviscid Solver (Fortran, C++, Python) • Conclusions — Joaquim R. R. A. Martins — 2
Sensitivity Analysis • Motivation: – Gradient-based optimization requires accurate sensitivity information. • Typical approaches: – Finite-Difference: easy to implement, but lacks robustness and accuracy. – Algorithmic/Automatic/Computational Differentiation: accurate, ease of implementation varies. – Analytic Methods: efficient and accurate, but long development and implementation times. • Complex-Step Method: accurate and robust; easy to implement and maintain. — Joaquim R. R. A. Martins — 3
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 R. R. A. Martins — 4
Complex-Step Derivative Approximation If the function f ( x ) is analytic, it can be expanded as a Taylor series with a complex 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! f ′ ( x ) ≈ Im [ f ( x + ih )] ⇒ h No subtraction! — Joaquim R. R. A. Martins — 5
Why Improve the Complex-Step Method? • Improvements necessary because, � � � arcsin( z ) = − i log iz + 1 − z 2 , may yield a zero derivative... • How? If z = x + ih , where x = O (1) and h = O (10 − 20 ) then in the addition, iz + z = ( x − h ) + i ( x + h ) h vanishes when using finite precision arithmetic. • Would like to keep the real and imaginary parts separate. — Joaquim R. R. A. Martins — 6
Alternative Definition for “sin” • Complex definition of sine also problematic, sin( z ) = e iz − e − iz . 2 i • Complex trigonometric relation gives better alternative: sin( x + ih ) = sin( x ) cosh( h ) + i cos( x ) sinh( h ) . • Note that for small h this simplifies to, sin( x + ih ) ≈ sin( x ) + ih cos( x ) . — Joaquim R. R. A. Martins — 7
Alternative Definition for “arcsin” • From the standard complex definition, � � � arcsin( z ) = − i log iz + 1 − z 2 . • Need real and imaginary parts calculated separately. • Linearize in h about h = 0 , h arcsin( x + ih ) ≡ arcsin( x ) + i √ 1 − x 2 . — Joaquim R. R. A. Martins — 8
The Connection • The new function definitions in these examples can be generalized: f ( x + ih ) ≡ f ( x ) + ih∂f ( x ) ∂x . • The real part is the real function and the imaginary part is the derivative multiplied by h . • Defining functions this way, the complex-step method is the same as algorithmic differentiation. • For small enough step, using finite precision arithmetic, the complex-step method is the same as algorithmic differentiation. — Joaquim R. R. A. Martins — 9
Algorithmic Differentiation vs. Complex Step Look at a simple operation, e.g. f = x 1 x 2 , Algorithmic Complex-Step h 1 = 10 − 20 ∆ x 1 = 1 ∆ x 2 = 0 h 2 = 0 f = x 1 x 2 f = ( x 1 + ih 1 )( x 2 + ih 2 ) ∆ f = x 1 ∆ x 2 + x 2 ∆ x 1 f = x 1 x 2 − h 1 h 2 + i ( x 1 h 2 + x 2 h 1 ) d f/dx 1 = ∆ f d f/dx 1 = Im f/h Complex-step method computes one extra term. • Other functions are similar: – Superfluous calculations are made. – For h ≤ x × 10 − 20 they vanish but still affect speed. — Joaquim R. R. A. Martins — 10
Previously Unresolved Issues Large body of research on algorithmic differentiation can now be applied to the complex-step method: • Singularities: non-analytic points • if statements: piecewise function definitions • Convergence for iterative solvers • Other issues addressed by the algorithmic differentiation community — Joaquim R. R. A. Martins — 11
Automatic Differentiation Implementations • Algorithmic Differentiation – Source transformation (ADIFOR, ADIC): resulting code is unmaintainable. – Derived datatype and operator overloading (ADOL-F, ADOL-C): far fewer changes are necessary in source code, requires object-oriented language. • Complex Step: – Even fewer changes are required. – Resulting code is maintainable. – Can be easily implemented in any programming language that supports complex arithmetic. — Joaquim R. R. A. Martins — 12
Fortran Implementation • complexify.f90 : a module that defines additional functions and operators for complex arguments. • Complexify.py : Python script that makes necessary changes to source code, e.g., type declarations. • New improvements: – Script is more versatile: ∗ Compatible with many more platforms and compilers. ∗ Supports MPI based parallel implementations. ∗ Resolves some of the input and output issues. – Some of the function definitions were improved: tangent, inverse and hyperbolic trigonometric functions. — Joaquim R. R. A. Martins — 13
C/C++ Implementations • complexify.h : defines additional functions and operators for the complex-step method. • derivify.h : simple algorithmic differentiation. Defines a new type which contains the value and its derivative. • Compared run time with real-valued code: – Complexified version: ≈ × 3 – Algorithmic differentiation version: ≈ × 2 — Joaquim R. R. A. Martins — 14
2D Flow Solver: FLO82 -2.00 + + + + + + + + + + + + + + + + ++++++++++++++ + + + + + + + + + + + + + + + + -1.60 + + + + + -1.20 + + + -0.80 + +++++++++++++++++++ + Cp + -0.40 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.00 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 0.40 + + + + + + + + + 0.80 + + + + + + + 1.20 RAE 2822 Airfoil α = 3 . 0 o , M ∞ = 0 . 70 — Joaquim R. R. A. Martins — 15
Convergence of C D and ∂C D /∂M ∞ 0 10 C D ∂ C D / ∂ M ∞ −5 10 Relative Error, ε −10 10 −15 10 50 100 150 200 250 300 350 400 450 500 Iterations Imaginary part converges at similar rate but lagged. — Joaquim R. R. A. Martins — 16
Sensitivity Estimate vs. Step Size 0 10 Complex−Step Finite−difference −5 10 Relative Error, ε −10 10 −15 10 −10 −20 −30 10 10 10 Step Size, h Complex-step is accurate from h = 10 − 10 to h = 10 − 321 . — Joaquim R. R. A. Martins — 17
3D Aero-Structural Design Optimization Framework • Aerodynamics: SYN107-MB, a parallel, multiblock Euler flow solver. • Structures: detailed finite element model with plates and trusses. • Coupling: high-fidelity, consistent and conservative. • Geometry: centralized database for exchanges (jig shape, pressure distributions, displacements.) • Coupled-adjoint sensitivity analysis — Joaquim R. R. A. Martins — 18
Aero-Structural Model and Solution — Joaquim R. R. A. Martins — 19
Wing Structural Model Detail — Joaquim R. R. A. Martins — 20
Convergence of C D and ∂C D /∂b 1 C D 0 10 ∂ C D / ∂ b 1 −2 10 Reference Error, ε −4 10 −6 10 −8 10 100 200 300 400 500 600 700 800 Iterations — Joaquim R. R. A. Martins — 21
Sensitivity Estimate vs. Step Size 0 10 Complex−Step −1 Finite−difference 10 −2 10 Relative Error, ε −3 10 −4 10 −5 10 −6 10 −5 −10 −15 10 10 10 Step Size, h Finite-difference is practically useless. — Joaquim R. R. A. Martins — 22
Sensitivity of C D to Hicks-Henne Functions Complex−Step, h = 1 × 10 −20 0.15 Finite−Difference, h = 1 × 10 −2 0.1 ∂ C D / ∂ b i 0.05 0 −0.05 2 4 6 8 10 12 14 16 18 Shape variable, i Much effort expended for a finite-difference result this reasonable. — Joaquim R. R. A. Martins — 23
Supersonic Viscous/Inviscid Solver Tool for preliminary design of natural laminar flow supersonic aircraft • Transition prediction • Viscous and inviscid drag • Design optimization – Wing planform and airfoil design – Wing-Body intersection design — Joaquim R. R. A. Martins — 24
Supersonic Viscous/Inviscid Solver • Python wrapper defines geometry • CH GRID automatic grid generator – Wing only or wing-body – Complexified with our script • CFL3D calculates Euler solution – Version 6 includes complex-step – New improvements incorporated • C++ post-processor for the . . . • Quasi-3D boundary-layer solver – Laminar and turbulent – Transition prediction – C++ algorithmic differentiation • Python wrapper collects data and computes structural constraints — Joaquim R. R. A. Martins — 25
Recommend
More recommend
