Adjoint code development and optimization using automatic - - PowerPoint PPT Presentation

Adjoint code development and optimization using automatic differentiation (AD)

• Praveen. C

Computational and Theoretical Fluid Dynamics Division National Aerospace Laboratories Bangalore - 560 037

CTFD Seminar, March 31, 2006

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 1 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 2 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 3 / 56

Mathematical formulation

Constrained minimization problem min

α J(α, u)

subject to R(α, u) = 0 Find δα such that δJ < 0 δJ = ∂J ∂αδα + ∂J ∂u δu = ∂J ∂αδα + ∂J ∂u ∂u ∂αδα = Gδα Steepest descent δα = −ǫG⊤ δJ = −ǫGG⊤ = −ǫG2 < 0

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 4 / 56

Mathematical formulation

Constrained minimization problem min

α J(α, u)

subject to R(α, u) = 0 Find δα such that δJ < 0 δJ = ∂J ∂αδα + ∂J ∂u δu = ∂J ∂αδα + ∂J ∂u ∂u ∂αδα = Gδα Steepest descent δα = −ǫG⊤ δJ = −ǫGG⊤ = −ǫG2 < 0

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 5 / 56

Sensitivity approach

Linearized state equation ∂R ∂α δα + ∂R ∂u δu = 0

• r

∂R ∂u ∂u ∂α = −∂R ∂α Solve sensitivity equation iteratively ∂ ∂t ∂u ∂α + ∂R ∂u ∂u ∂α = −∂R ∂α Same eigenvalues as original system Costly if number of α are large

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 6 / 56

Adjoint approach

We have δJ = ∂J ∂αδα + ∂J ∂u δu and ∂R ∂α δα + ∂R ∂u δu = 0 Introduce a new unknown v δJ = ∂J ∂αδα + ∂J ∂u δu + v⊤ ∂R ∂α δα + ∂R ∂u δu

• =

∂J ∂α + v⊤ ∂R ∂α

• δα +

∂J ∂u + v⊤ ∂R ∂u

• δu

Adjoint equation ∂R ∂u ⊤ v = − ∂J ∂u ⊤ Iterative solution ∂v ∂t + ∂R ∂u ⊤ v = − ∂J ∂u ⊤

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 7 / 56

Continuous vs Discrete

Continuous approach:

Start with governing PDE Derive adjoint PDE and boundary conditions Discretize adjoint PDE and solve Must be re-derived whenever cost function changes Not all cost functions admissible

Discrete approach:

Start with discrete approximation like FVM Derive discrete adjoint equations Solve discrete adjoint equations True sensitivity of discrete solution

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 8 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 9 / 56

Techniques for computing gradients

Hand differentiation Finite difference method Complex variable method Automatic differentiation

Computer code to compute some function Chain rule of differentiation Generates a code to compute derivatives ADIFOR, ADOLC, ODYSEE, TAMC, TAF , TAPENADE

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 10 / 56

Differentiation: Example

A simple example f = (xy + sin x + 4)(3y2 + 6) Computer code, f = t10 t1 = x t2 = y t3 = t1t2 t4 = sin t1 t5 = t3 + t4 t6 = t5 + 4 t7 = t2

2

t8 = 3t7 t9 = t8 + 6 t10 = t6t9

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 11 / 56

Automatic Differentiation: F77 code

subroutine costfunc(x, y, f) t1 = x t2 = y t3 = t1*t2 t4 = sin(t1) t5 = t3 + t4 t6 = t5 + 4 t7 = t2**2 t8 = 3.0*t7 t9 = t8 + 6.0 t10 = t6*t9 f = t10 end

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 12 / 56

Differentiation: Direct mode

Some definitions Ik := {i : i < k and tk depends explicitly on ti} ˙ ti := ∂ti ∂α, ti,k := ∂ti ∂tk Apply chain rule of differentiation ˙ tk = ∂tk ∂α =

• i∈Ik

∂ti ∂α ∂tk ∂ti =

• i∈Ik

˙ titk,i k = 1, 2, . . . , 9, 10 α = x, or y

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 13 / 56

Differentiation: Direct mode

Apply chain rule of differentiation t1 = x ˙ t1 = ˙ x t2 = y ˙ t2 = ˙ y t3 = t1t2 ˙ t3 = ˙ t1t2 + t1˙ t2 t4 = sin(t1) ˙ t4 = cos(t1)˙ t1 t5 = t3 + t4 ˙ t5 = ˙ t3 + ˙ t4 t6 = t5 + 4 ˙ t6 = ˙ t5 t7 = t2

2

˙ t7 = 2t2˙ t2 t8 = 3t7 ˙ t8 = 3˙ t7 t9 = t8 + 6 ˙ t9 = ˙ t8 t10 = t6t9 ˙ t10 = ˙ t6t9 + t6˙ t9 ˙ x = 1, ˙ y = 0, ˙ t10 = ∂f

∂x and ˙

x = 0, ˙ y = 1, ˙ t10 = ∂f

∂y

tapenade -d -vars "x y" -outvars f costfunc.f

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 14 / 56

Automatic Differentiation: Direct mode

SUBROUTINE COSTFUNC_D(x, xd, y, yd, f, fd) t1d = xd t1 = x t2d = yd t2 = y t3d = t1d*t2 + t1*t2d t3 = t1*t2 t4d = t1d*COS(t1) t4 = SIN(t1) t5d = t3d + t4d t5 = t3 + t4 t6d = t5d t6 = t5 + 4 t7d = 2*t2*t2d t7 = t2**2 t8d = 3.0*t7d t8 = 3.0*t7 t9d = t8d t9 = t8 + 6.0 t10d = t6d*t9 + t6*t9d t10 = t6*t9 fd = t10d f = t10 END

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 15 / 56

Differentiation: Reverse mode

Some definitions Jk := {i : i > k and ti depends explicitly on tk} ¯ ti := ∂f ∂ti = ∂t10 ∂ti , ti,k := ∂ti ∂tk Apply chain rule of differentiation in reverse ¯ tk = ∂t10 ∂tk =

• i∈Jk

∂t10 ∂ti ∂ti ∂tk =

• i∈Jk

¯ titi,k k = 10, 9, . . . , 2, 1

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 16 / 56

Differentiation: Reverse mode

Apply chain rule of differentiation in reverse t1 = x ¯ t10 = 1 t2 = y ¯ t9 = ¯ t10t10,9 = t6 t3 = t1t2 ¯ t8 = ¯ t9t9,8 = t6 t4 = sin(t1) ¯ t7 = ¯ t8t8,7 = 3t6 t5 = t3 + t4 ¯ t6 = ¯ t10t10,6 = t9 t6 = t5 + 4 ¯ t5 = ¯ t6t6,5 = t9 t7 = t2

2

¯ t4 = ¯ t5t5,4 = t9 t8 = 3t7 ¯ t3 = ¯ t5t5,3 = t9 t9 = t8 + 6 ¯ t2 = ¯ t7t7,2 + ¯ t3t3,2 = 6t2t6 + t1t9 t10 = t6t9 ¯ t1 = ¯ t4t4,1 + ¯ t3t3,1 = t9 cos(t1) + t9t2 ¯ t1 = ∂f

∂x , ¯

t2 = ∂f

∂y

tapenade -b -vars "x y" -outvars f costfunc.f

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 17 / 56

Automatic Differentiation: Reverse mode

SUBROUTINE COSTFUNC_B(x, xb, y, yb, f, fb) t1 = x t2 = y t3 = t1*t2 t4 = SIN(t1) t5 = t3 + t4 t6 = t5 + 4 t7 = t2**2 t8 = 3.0*t7 t9 = t8 + 6.0 t10b = fb t6b = t9*t10b t9b = t6*t10b t8b = t9b t7b = 3.0*t8b t5b = t6b t3b = t5b t2b = t1*t3b + 2*t2*t7b t4b = t5b t1b = t2*t3b + COS(t1)*t4b yb = t2b xb = t1b fb = 0.0 END

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 18 / 56

Direct versus reverse AD

Direct mode: cost is proportional to number of independent variables Reverse mode: cost does not depend on the number of independent variables Reverse mode preferred when number of independent variables is large Scalar output f

Direct mode gives ∇f · s for given vector s Reverse mode gives ∇f

Vector output R ∈ Rn

Direct mode gives ∇R · s for given vector s Reverse mode gives (∇R)⊤ · s

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 19 / 56

Issues in reverse AD

Intermediate variables required in reverse order Some variables may be over-written Iterative solvers

Only fi nal solution required AD differentiates the iterative loop Intermediate solutions stored in stack Huge memory requirements Not practical for large problems

Piecemeal differentiation approach (Courty et al., Giles et al.):

Modular fl

• w solver

Adjoint solver written manually AD for differentiating the modules

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 20 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 21 / 56

Discrete 1-D adjoint equations

Finite volume residual for i’th cell Ri := Fi+1/2 − Fi−1/2 = 0 Numerical flux function F = F(X, Y), Fi+1/2 = F(Ui, Ui+1) Perturbation equation δRi =

∂ ∂X Fi+1/2δUi + ∂ ∂Y Fi+1/2δUi+1

− ∂

∂X Fi−1/2δUi−1 − ∂ ∂Y Fi−1/2δUi + ∂Ri ∂α δα

= 0 Introduce adjoint variable Vi for i’th cell δJ = ∂J ∂αδα +

• i

∂J ∂Ui δUi +

• i

V ⊤

i δRi

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 22 / 56

Discrete adjoint equations

Collecting terms containing δUi δJ =

• i

∂J ∂Ui + V ⊤

i−1

∂ ∂Y Fi−1/2 + V ⊤

i

∂ ∂X Fi+1/2 − ∂ ∂Y Fi−1/2

• −

V ⊤

i+1

∂ ∂X Fi+1/2

• δUi + . . .

Adjoint equation for i’th cell ∂J ∂Ui ⊤ + ∂ ∂Y Fi−1/2 ⊤ Vi−1 + ∂ ∂X Fi+1/2 − ∂ ∂Y Fi−1/2 ⊤ Vi − ∂ ∂X Fi+1/2 ⊤ Vi+1 = 0

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 23 / 56

Example fl

• w solver

While u is not converged res = 0.0 fluxinflow(u(1), res(1)) do i=1,N-1 fluxinterior(u(i), u(i+1), res(i), res(i+1)) enddo fluxoutflow(u(N), res(N)) do i=1,N u(i) = u(i) – (dt/dx)*res(i) enddo endwhile cost=0.0 do i=1,N costfunc(u(i), cost) enddo

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 24 / 56

Example adjoint solver

costb=1.0 do i=1,N costfunc_b(u(i), ub1(i), cost, costb) enddo v=0.0 While v is not converged ub2 = 0.0 fluxinflow_b(u(1), ub2(1), res(1), v(1)) do i=1,N-1 fluxinterior_b(u(i), ub2(i), u(i+1), ub2(i+1), res(i), v(i), res(i+1), v(i+1)) enddo fluxoutflow_b(u(N), ub2(N), res(N), v(N)) do i=1,N v(i) = v(i) – (dt/dx)*(ub1(i) + ub2(i)) enddo endwhile

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 25 / 56

Quasi 1-D fl

• w

Quasi 1-D flow in a duct ∂ ∂t (hU) + ∂ ∂x (hf) = dh dx P, x ∈ (a, b) t > 0 U =   ρ ρu E   , f =   ρu p + ρu2 (E + p)u   , P =   p   h(x) = cross-section height of duct Inverse design: find shape h to get pressure distribution p∗ Optimization problem: find the shape h which minimizes J = b

a

(p − p∗)2dx

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 26 / 56

Quasi 1-D fl

• w

Finite volume scheme hi dUi dt + hi+1/2Fi+1/2 − hi−1/2Fi−1/2 ∆x = (hi+1/2 − hi−1/2) ∆x Pi Discrete cost function J =

N

• i=1

(pi − p∗

i )2

Control variables h1/2, h1+1/2, . . . , hi+1/2, . . . , hN+1/2 N = 100

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 27 / 56

Duct shape

2 4 6 8 10 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 x h Duct shape Target shape Current shape

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 28 / 56

Target pressure distribution p∗

2 4 6 8 10 0.5 1 1.5 2 2.5 x Pressure Target pressure AUSMDV KFVS LF

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 29 / 56

Current pressure distribution

2 4 6 8 10 0.5 1 1.5 2 2.5 x Pressure Starting pressure AUSMDV KFVS LF

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 30 / 56

Adjoint density

2 4 6 8 10 −35 −30 −25 −20 −15 −10 −5 5 x Adjoint density AUSMDV KFVS LF

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 31 / 56

Convergence history: Explicit Euler

1000 2000 3000 4000 5000 6000 7000 −10 −9 −8 −7 −6 −5 −4 −3 −2 −1 No of iterations Residue Convergence history with AUSMDV flux Flow Adjoint

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 32 / 56

Shape gradient

2 4 6 8 10 5 10 15 20 25 30 35 x Gradient Gradient AUSMDV KFVS LF

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 33 / 56

Shape gradient

2 2.5 3 3.5 4 4.5 5 5.5 6 −5 5 10 15 20 25 30 35 40 x Gradient Gradient AUSMDV KFVS LF

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 34 / 56

Validation of Shape gradient

2 4 6 8 10 5 10 15 20 25 30 35 x Gradient Gradient with AUSMDV flux A B C

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 35 / 56

Validation of shape gradient

∂J ∂h ≈ J(h + ∆h) − J(h − ∆h) 2∆h ∆h A B C 0.01 0.4191069499 35.18452823 2.545316345 0.001 0.4231223499 36.10982621 2.556461900 0.0001 0.4231624999 36.11933154 2.556573499 0.00001 0.4231599998 36.11942125 2.556575000 0.000001 0.4231999994 36.11942305 2.556550001 0.0000001 0.4229999817 36.11942329 2.556499971 AD 0.4231628330 36.11941951 2.556574450

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 36 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 37 / 56

Gradient smoothing

Non-smooth gradients especially in the presence of shocks Smooth using an elliptic equation

• 1 − ǫi

d2 dx2

• ¯

gi = gi ǫi = {|gi+1 − gi| + |gi − gi−1|}Li Li = |gi+1 − 2gi + gi−1| max(|gi+1 − gi| + |gi − gi−1|, TOL) Finite difference with Jacobi iterations

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 38 / 56

Gradient smoothing

2 4 6 8 10 5 10 15 20 25 30 35 x Gradient Gradient using AUSMDV flux Original gradient Smoothed gradient

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 39 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 40 / 56

Quasi 1-D optimization: Shape

1 2 3 4 5 6 7 8 9 10 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 x h Target Initial

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 41 / 56

Quasi 1-D optimization: Final shape

1 2 3 4 5 6 7 8 9 10 1 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9 x h Target Initial Final

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 42 / 56

Quasi 1-D optimization: Pressure

1 2 3 4 5 6 7 8 9 10 0.5 1 1.5 2 2.5 x Pressure Target Initial Final

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 43 / 56

Quasi 1-D optimization: Convergence

5 10 15 20 25 30 35 40 45 50 10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10

• No. of iterations

Cost/Gradient Cost Gradient

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 44 / 56

Quasi 1-D optimization: Adjoint density

2 4 6 8 10 −25 −20 −15 −10 −5 x Adjoint density Initial Final

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 45 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 46 / 56

2-D fi nite volume fl

• w solver

Unstructured triangular grid Vertex-centroid scheme (Frink, Jameson) Vertex values obtained using inverse-distance averaging with correction for linear consistency

i p

R L

up =

• i

wi ri ui

• i

wi ri

Convergence acceleration using matrix-free LUSGS with spectral radius approximation

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 47 / 56

2-D fi nite volume adjoint solver

Vertex averaging including weights are differentiated For transonic flow limiter is also differentiated Adjoint solver requires 38% more memory compared to flow solver Time per adjoint iteration = Twice the time per flow solver iteration Convergence acceleration using LUSGS with spectral-radius approximation Automated using Makefile Adjoint validation for J(α) = Cl

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 48 / 56

Grid for NACA-0012: 160 × 40

Number of vertices = 6560 Number of triangles = 12800

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 49 / 56

Validation of adjoints: Lift-curve slope, subsonic fl

• w

0.5 1 1.5 2 2.5 3 0.05 0.1 0.15 0.2 0.25 0.3 0.35 Angle of attack in degrees Cl

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 50 / 56

Validation of adjoints: Lift-curve slope, transonic fl

• w

0.5 1 1.5 2 2.5 3 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 Angle of attack in degrees Cl

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 51 / 56

Lift-curve slope: Clα

AOA M∞ = 0.4 M∞ = 0.8 Clα Clα/Cth

lα

Clα Clα/Cth

lα

0.0 7.2095 1.0516 13.7506 1.3130 0.5 7.2069 1.0512 14.0556 1.3422 1.0 7.2000 1.0502 13.9544 1.3325 1.5 7.1901 1.0487 13.2116 1.2616 2.0 7.1789 1.0471 12.8314 1.2253 2.5 7.1675 1.0455 11.4077 1.0893 3.0 7.1562 1.0438 10.7135 1.0230

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 52 / 56

Convergence history: α = 3o

M∞ = 0.4 M∞ = 0.8

1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 1 10 100 2000 4000 6000 8000 10000 12000 Residue

• No. of iterations

Flow Adjoint 1e-07 1e-06 1e-05 1e-04 0.001 0.01 0.1 1 10 1000 2000 3000 4000 5000 6000 7000 8000 9000 10000 Residue

• No. of iterations

Flow Adjoint

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 53 / 56

Adjoint density for lift coeffi cient: α = 3

• M∞ = 0.4

M∞ = 0.8

Adjoint_Density, min = -57.979, max = 58.0587 Adjoint_Density, min = -140.596, max = 134.291

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 54 / 56

Outline

1

Mathematical formulation

2

Computing gradients

3

Quasi 1-D flow

4

Gradient smoothing

5

Quasi 1-D optimization: Pressure matching

6

2-D adjoint Euler solver

7

Future work

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 55 / 56

Future work

How to deform the grid ? How to account for grid perturbation ? Smoothing of gradients in 2-D Application to 2-D optimization problems Extension to viscous flows Grid adaptation

• Praveen. C (CTFD, NAL)

Adjoints and optimization CTFD Seminar 56 / 56