Automatic Differentiation-based perturbation methods for - PowerPoint PPT Presentation

Automatic Differentiation-based perturbation methods for uncertainties and errors Anca Belme, Massimiliano Martinelli, Laurent Hasco¨ et, Val´ erie Pascual, Alain Dervieux Tropics project http://www-sop.inria.fr/tropics/ INRIA Sophia-Antipolis Support of FP6-NODESIM-CFD 44` eme AAAF, Nantes, 23-25 mars 2009

Overview: 1. Progress in Automatic Differentiation, 2. Numerical error reduction, 3. Perturbation methods for uncertainties.

1. Progress in Automatic Differentiation The TAPENADE program differentiator: TAPENADE (115 000 lines of JAVA) transforms a FORTRAN program into a FORTRAN program which computes the derivative of the given FORTRAN program. TAPENADE is organised around an intermediate step, the analysis and differentation of an abstract imperative language. http://www-sop.inria.fr/tropics/ Differentiation Engine Imperative Language Analyzer (IL) (IL) Fortran77 parser Fortran77 printer Fortran95 parser Fortran95 printer C parser C printer Overall Architecture of tapenade

Differentiability/Differentiation modes Direct/Tangent mode: - Differentiated code computes ( c, δc ) �→ j ′ ( c ) .δc - Gˆ ateaux derivative of a functional - Computational cost factor: α T ≈ 4 - Chain rule applies in same order as initial code - Does not store intermediate variables Backward/Reverse mode: - Differentiated code computes ( c, δj ) �→ ( j ′ ( c )) ∗ .δj - gradient of a functional - Computational cost factor: α R ≈ 5 - chain rule applies in reverse order to initial code. - Needs storing intermediate variables in a first direct sweep .

Reverse differentiation strategies for iterative solvers Find W h such that | Ψ h ( W h ) | ≤ ε If the iteration for W h is stopped by a test with ε > 0 , then W h depends on initial conditions of solution algorithm, in a similar manner to an unsteady system. Also W h is only piecewise differentiable. It then is possible to apply R-AD to the whole programme, including iterative solution algorithms. Due to dependance to intermediate states they need be stored (or recomputed). For this case and for unsteady problems, efficients strategies of storage or recomputations can be introduced by user by means of “check point” con- trols. Assuming ε = 0 , state W h does not depend on initial conditions of solution algorithm and is differentiable with respect to c h . The iteration should be pointed as “independant” or parallel in order to avoid unnecessary storage of intermediate values.

1. Progress in Automatic Differentiation The TAPENADE program differentiator, recent advances: In its Version 3(alpha), TAPENADE takes into account user directives in- serted into sources. Directive II-LOOP allows for a more efficient reverse differentiation of par- allel loops by inhibiting storages useless in “parallel” case. Directive NOCHECKPOINT can be used for optimal checkpointing in reverse differentiation. TAPENADE lets the user specify finely which procedure calls must be checkpointed or not. TAPENADE V3 has been extended to deal with Fortran95 and with ANSI C. TAPENADE is optimised for successive reverse differentiations.

2. Numerical error reduction Abstract representation of the Partial Differential Equation: Ψ( W ) = 0 . Discretisation of the PDE: Ψ h ( W h ) = 0 ∈ R N W h ∈ R N , W h = [( W h ) i ] . Operators between continuous and discrete: R h : R N → V ⊂ L 2 (Ω) ∩ C 0 (¯ Ω) v h �→ R h v h → R N T h : V v �→ T h v W h ( x, y, z ) = ( R h W h )( x, y, z ) ; W − W h = small? corrector?

2. Numerical error reduction A priori estimate: Ψ h ( W ) − Ψ h ( W h ) = Ψ h ( W ) ⇒ W − W h ≈ [ ∂ Ψ h ] − 1 Ψ h ( W ) , ∂W h in practice: δW h = R h [ ∂ Ψ h ] − 1 (Ψ h − T h Ψ)( W ( h ) ) . ∂W h A posteriori estimate: Ψ( W ) − Ψ( W h ) = − Ψ( W h ) ⇒ W − W h ≈ − [ ∂ Ψ ∂W ] − 1 Ψ( W h ) , in practice: δW h = − R h [ ∂ Ψ h ] − 1 T h Ψ( R h W h ) , ∂W h quadratures formulas can be used for Ψ( R h W h ) .

2. Numerical error reduction Goal-oriented error: j ( W ) = ( g, W ) L 2 (Ω) s.t. Ψ( W ) = 0 ( ∂ Ψ ∂W ) ∗ p = g j h = ( g, R h W h ) L 2 (Ω) s.t. Ψ h ( W h ) = 0 g h = T h g ( ∂ Ψ h ) ∗ p h = g h ⇔ p h = [ ∂ Ψ h ] −∗ T h g h ∂ W h ∂W h p h = R h p h . A fundamental assumption of the present analysis is that this discrete adjoint is a good enough approximation of continuous adjoint.

2. Numerical error reduction A posteriori goal-oriented errors and correctors Field corrector δW h = − R h [ ∂ Ψ h ] − 1 T h Ψ( W h ) ∂W h Direct-linearized goal-oriented corrector  g , R h [ ∂ Ψ h   ] − 1 T h Ψ( W h ) δ 1 j = −    ∂W h L 2 (Ω) Adjoint-based goal-oriented correctors (Giles-Pierce) δ 1 j = − ( p h , T h Ψ( W h )) L 2 (Ω)

2. Numerical error reduction An equivalence: Assuming that we have chosen: T h = R ∗ h , then direct-linearised and adjoint-based correctors are identical  g , R h [ ∂ Ψ h  [ ∂ Ψ h     ] − 1 T h Ψ( W h ) ] −∗ R ∗ − = − h g , T h Ψ( W h )       ∂W h ∂W h L 2 (Ω) R N h [ ∂ Ψ h   ] −∗ R ∗  T ∗ = − h g , Ψ( W h )    ∂W h L 2 (Ω)

2. Numerical error reduction A preliminary numerical example 3D Euler flow around a wing with NACA0012 section. 1. Direct calculation 2. Direct calculation + Correction CPU is 30% larger, error on entropy deviation is 10 times smaller. Entropy spurious generation for a direct computation of a steady flow and for a corrected one.

2. Numerical error reduction by correction, conclud’d Which A posteriori goal-oriented errors and correctors Direct-linearized goal-oriented corrector  g , R h [ ∂ Ψ h   ] − 1 T h Ψ( W h ) δ 1 j = −    ∂W h L 2 (Ω) Adjoint-based goal-oriented correctors δ 1 j = − ( p h , T h Ψ( W h )) L 2 (Ω) Adjoints, as compared with linearised perturbations, are slightly more com- plex to compute for steady PDE, and much more complex for unsteady PDE’s (despite TAPENADE is well adapted for this exercice). For com- puting correctors and for several mesh adaptation applications, the “Direct- linearized goal-oriented corrector” can be used for an identical result.

2. A priori numerical error reduction, mesh adaption Minimizing P1-Interpolation error (Loseille, PhD, Paris VI, 2008) 2  | ∂ 2 u ξ + | ∂ 2 u   || u − π M u || 2 ∂ξ 2 | .m 2 ∂η 2 | .m 2 � = dxdy    η   where ξ and η are directions of diagonalization of the Hessian of u . min E M M under the constraint N M = N. This can be solved analytically.

2. A priori numerical error reduction, mesh adaption Optimal Metric for interpolation − 5 / 6 | ∂ 2 u 1 / 6   | ∂ 2 u ∂η 2 | ∂ξ 2 | 0 M opt = C   N R − 1    R . (1)     − 5 / 6 | ∂ 2 u 1 / 6 | ∂ 2 u   0 ∂ξ 2 | ∂η 2 |    � 2 6 dxdy . | ∂ 2 u ∂ξ 2 | . | ∂ 2 u � with: C = ∂η 2 | � But this only minimises an interpolation error...

2. A priori numerical error reduction, mesh adaption Approximation a priori error (with Alauzet & Loseille) Steady Euler equations: � 5 , ∀ φ ∈ V, H 1 (Ω) � W ∈ V = Γ φ ˆ � � (Ψ( W ) , φ ) = Ω φ ∇ . F ( W ) d Ω − F ( W ) . n d Γ = 0 where Γ = ∂ Ω and ˆ F is B.C. fluxes. Mixed-Element-Volume appoximation: Γ h φ h Π h ˆ � � ∀ φ h ∈ V h , Ω h φ h ∇ . Π h F ( W h ) d Ω h − F ( W h ) . n d Γ h = � − Ω h φ h D h ( W h ) d Ω h , where Π h is the usual elementwise linear interpolation and where D h holds for a numerical dissipation term.

2. A priori numerical error reduction, mesh adaption A priori adjoint-based error estimate: ∗  ∂ Ψ   ( g, W − W h ) ≈ ((Ψ h − Ψ)( W ) , P ) , with P = g,   ∂W  the optimal mesh is obtain after some transformations by solving: � Find M opt = Argmin M Ω |∇ P ||F ( W ) − π M F ( W ) | d Ω Γ | P || ( ¯ F ( W ) − π M ¯ � + F ( W )) . n | d Γ (2) under the constraint C ( M ) = N . Solved analytically as interpolation case. Remark: The adjoint-based formulation is compulsory.

Application to sonic boom : Hessian-based (Mach L2)

Goal-oriented + Hessian-based (“foot print” funct.)

Pressure track at ground

3. Hessian-based perturbation methods for uncertainties A differentiation problem: Given a computer program computing a functional j ( c ) for c ∈ R n j ( c ) = J ( c, W ( c )) , such that Ψ( c, W ) = 0 we want, applying a source-to-source programme differentiation soft- ware , viz. TAPENADE , to get a computer program computing the sec- ond derivatives: d 2 j dc 2 = ( H ii ) . The solution W ( c ) ∈ R N of state equation Ψ( c, W ) = 0 is supposed here to be steady-state. It can be obtained by explicit or implicit pseudo-time advancing techniques. Application to an unsteady problem can be treated in a similar manner, but with important extra difficulties.