Adjoint approach to optimization Praveen. C - - PowerPoint PPT Presentation

adjoint approach to optimization
SMART_READER_LITE
LIVE PREVIEW

Adjoint approach to optimization Praveen. C - - PowerPoint PPT Presentation

Adjoint approach to optimization Praveen. C praveen@math.tifrbng.res.in Tata Institute of Fundamental Research Center for Applicable Mathematics Bangalore 560065 http://math.tifrbng.res.in Health, Safety and Environment Group BARC 6-7


slide-1
SLIDE 1

Adjoint approach to optimization

  • Praveen. C

praveen@math.tifrbng.res.in

Tata Institute of Fundamental Research Center for Applicable Mathematics Bangalore 560065 http://math.tifrbng.res.in

Health, Safety and Environment Group BARC 6-7 October, 2010

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 1 / 58

slide-2
SLIDE 2

Outline

  • Minimum of a function
  • Constrained minimization
  • Finite difference approach
  • Adjoint approach
  • Automatic differentiation
  • Example
  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 2 / 58

slide-3
SLIDE 3

Minimum of a function

x f(x) x0 f'(x0)=0

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 3 / 58

slide-4
SLIDE 4

Steepest descent method

x1 x2 grad f(x1,x2)

xn+1 = xn − sn∇f(xn)

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 4 / 58

slide-5
SLIDE 5

Objectives and controls

  • Objective function

J(α) = J(α, u) mathematical representation of system performance

  • Control variables α

◮ Parametric controls α ∈ Rn ◮ Infinite dimensional controls α : X → Y ◮ Shape α ∈ set of admissible shapes

  • State variable u: solution of an ODE or PDE

R(α, u) = 0 ⇓ u = u(α)

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 5 / 58

slide-6
SLIDE 6

Gradient-based minimization: blackbox/FD approach

min

β∈RN I(β, Q(β))

  • Initialize β0, n = 0
  • For n = 0, . . . , Niter

◮ Solve R(βn, Qn) = 0 ◮ For j = 1, . . . , N ⋆ βn (j) = [βn 1 , . . . , βn j + ∆βj, . . . , βn N]⊤ ⋆ Solve R(βn (j), Qn (j)) = 0 ⋆ dI dβj ≈ I(βn

(j),Qn (j))−I(βn,Qn)

∆βj ◮ Steepest descent step

βn+1 = βn − sn dI dβ (βn)

Cost of FD-based steepest-descent

Cost = O(N + 1)Niter = O(N + 1)O(N) = O(N2)

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 6 / 58

slide-7
SLIDE 7

Accuracy of FD: Choice of step size

d dxf(x0) = f(x0 + δ) − f(x0) δ + O(δ) In principle, if we choose a small δ, we reduce the error. But computers have finite precision. Instead of f(x0) the computers gives f(x0) + O(ǫ) where ǫ = machine precision. [f(x0 + δ) + O(ǫ)] − [f(x0) + O(ǫ)] δ = f(x0 + δ) − f(x0) δ + C1 ǫ δ = d dxf(x0) + C2δ + C1 ǫ δ

  • Total error

, C1, C2 depend on f, x0, precision Total error is least when δ = δopt = C3 √ǫ, C3 =

  • C1/C2
  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 7 / 58

slide-8
SLIDE 8

Accuracy of FD: Choice of step size

Error δ δopt C1δ C2ǫ/δ

Precision ǫ δopt Single 10−8 10−4 Double 10−16 10−8 See: Brian J. McCartin, Seven Deadly Sins of Numerical Computation The American Mathematical Monthly, Vol. 105, No. 10 (Dec., 1998), pp. 929-941

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 8 / 58

slide-9
SLIDE 9

Drag gradient using FD (Samareh)

Errors in Finite Difference Approximations

Scaled Step Size % Error

10

  • 9

10

  • 8

10

  • 7

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10

  • 6

10

  • 5

10

  • 4

10

  • 3

10

  • 2

10

  • 1

10 10

1

Mid chord LE Sweep2 Tip chord LE sweep3 Twist (root) Twist (mid) Twist (tip)

Tip chord Mid chord LE Sweep2 LE Sweep3 Twist (root) Twist (mid) Twist (tip)

Roundoff Truncation

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 9 / 58

slide-10
SLIDE 10

Iterative problems

I(β, Q), where R(β, Q) = 0

  • Q is implicitly defined, require an iterative solution method.
  • Assume a Q0 and iterate Qn −

→ Qn+1 until ||R(β, Qn)|| ≤ TOL

  • If TOL is too small, need too many iterations
  • Many problems, we cannot reduce ||R(β, Qn)|| to small values
  • This means that numerical value of I will be noisy

Finite difference will contain too much error, and is useless RAE5243 airfoil, Mach=0.68, Re=19 million, AOA=2.5 deg. iter Lift Drag 41496 0.824485788042416 1.627593747613790E-002 41497 0.824485782714867 1.627593516695762E-002 41498 0.824485777387834 1.627593285794193E-002 41499 0.824485772061306 1.627593054909022E-002 41500 0.824485766735297 1.627592824040324E-002

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 10 / 58

slide-11
SLIDE 11

Complex variable method

f(x0 + iδ) = f(x0) + iδf′(x0) + O(δ2) + iO(δ3) f′(x0) = 1 δ imag[f(x0 + iδ)] + O(δ2)

  • No roundoff error
  • We can take δ to be very small, δ = 10−20
  • Can be easily implemented

◮ fortran: redeclare real variables as complex ◮ matlab: no change

  • Iterative problems: β −

→ β + i∆β

◮ Obtain ˜

Q = Q(β + i∆β) by solving R(β + i∆β, ˜ Q) = 0

◮ Then gradient

I′(β) ≈ 1 ∆β imag[I(β + i∆β, Q(β + i∆β)]

  • Computational cost is O(N2) or higher (due to complex

arithmetic)

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 11 / 58

slide-12
SLIDE 12

Objectives and controls

  • Objective function

J(α) = J(α, u) mathematical representation of system performance

  • Control variables α

◮ Parametric controls α ∈ Rn ◮ Infinite dimensional controls α : X → Y ◮ Shape α ∈ set of admissible shapes

  • State variable u: solution of an ODE or PDE

R(α, u) = 0

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 12 / 58

slide-13
SLIDE 13

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 ∂αδα = ∂J ∂α + ∂J ∂u ∂u ∂α

  • δα =: Gδα
  • Steepest descent

δα = −ǫG⊤ δJ = −ǫGG⊤ = −ǫG2 < 0

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 13 / 58

slide-14
SLIDE 14

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 ∂αδα = ∂J ∂α + ∂J ∂u ∂u ∂α

  • δα =: Gδα
  • Steepest descent

δα = −ǫG⊤ δJ = −ǫGG⊤ = −ǫG2 < 0

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 14 / 58

slide-15
SLIDE 15

Sensitivity approach

  • Linearized state equation R(α, u) = 0

∂R ∂α δα + ∂R ∂u δu = 0

  • r

∂R ∂u ∂u ∂α = −∂R ∂α

  • Solve sensitivity equation iteratively, e.g.,

∂ ∂t ∂u ∂α + ∂R ∂u ∂u ∂α = −∂R ∂α

  • Gradient

dJ dα = ∂J ∂α + ∂J ∂u ∂u ∂α

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 15 / 58

slide-16
SLIDE 16

Sensitivity approach: Computational cost

  • n design variables: α = (α1, . . . , αn)
  • Solve primal problem R(α, u) = 0 to get u(α)
  • For i = 1, . . . , n

◮ Solve sensitivity equation wrt αi

∂R ∂u ∂u ∂αi = − ∂R ∂αi

◮ Compute derivative wrt αi

dJ dαi = ∂J ∂αi + ∂J ∂u ∂u ∂αi

  • One primal equation, n sensitivity equations

Computational cost = n + 1

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 16 / 58

slide-17
SLIDE 17

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 17 / 58

slide-18
SLIDE 18

Adjoint approach: Computational cost

  • n design variables: α = (α1, . . . , αn)
  • Solve primal problem R(α, u) = 0 to get u(α)
  • Solve adjoint problem

∂R ∂u ⊤ v = − ∂J ∂u ⊤

  • For i = 1, . . . , n

◮ Compute derivative wrt αi

dJ dαi = ∂J ∂αi + v⊤ ∂R ∂αi

  • One primal equation, one adjoint equation

Computational cost = 2, independent of n

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 18 / 58

slide-19
SLIDE 19

Adjoint: Two approaches

Continuous or Discrete or differentiate-discretize discretize-differentiate

PDE Adjoint PDE Discrete adjoint PDE Discrete PDE Discrete adjoint

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 19 / 58

slide-20
SLIDE 20

Techniques for computing gradients

  • Hand differentiation
  • Finite difference method
  • Complex variable method
  • Automatic Differentiation (AD)

◮ Computer code to compute J(α, u) and R(α, u) ◮ Chain rule of differentiation ◮ Generates a code to compute derivatives ◮ ADIFOR, ADOLC, ODYSEE, TAMC, TAF, TAPENADE

see http://www.autodiff.org Code for J(α, u) AD tool Code for

∂J ∂α

α

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 20 / 58

slide-21
SLIDE 21

Derivatives

  • Given a program P computing a function F

F : Rm → Rn X → Y

  • build a program that computes derivatives of F
  • X : independent variables
  • Y : dependent variables
  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 21 / 58

slide-22
SLIDE 22

Derivatives

  • Jacobian matrix: J =

∂Yj

∂Xi

  • Directional or tangent derivative

˙ Y = J ˙ X

  • Adjoint mode

¯ X = J⊤ ¯ Y

  • Gradients (n = 1 output)

J = ∂Y ∂Xi

  • = ∇Y
  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 22 / 58

slide-23
SLIDE 23

Forward differentiation

  • Program P is a sequence of computations Fk
  • T0 = X, given
  • k’th line

tk = Fk(Tk−1), Tk = Tk−1 tk

  • Function is a composition

F = Fp ◦ Fp−1 ◦ . . . ◦ F2 ◦ F1

  • Chain rule

˙ Y = F ′(X) ˙ X = F ′

p(Tp−1)F ′ p−1(Tp−2) . . . F ′ 2(T1)F ′ 1(T0) ˙

X X, ˙ X → Y, ˙ Y

  • cost( ˙

Y ) = 4 ∗ cost(Y )

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 23 / 58

slide-24
SLIDE 24

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 24 / 58

slide-25
SLIDE 25

F77 code: costfunc.f

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 25 / 58

slide-26
SLIDE 26

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 26 / 58

slide-27
SLIDE 27

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 27 / 58

slide-28
SLIDE 28

Backward differentiation

  • Program P is a sequence of computations Fk
  • T0 = X, given
  • k’th line

tk = Fk(Tk−1), Tk = Tk−1 tk

  • Function is a composition

F = Fp ◦ Fp−1 ◦ . . . ◦ F2 ◦ F1

  • Chain rule

¯ X = [F ′(X)]⊤ ¯ Y = [F ′

1(T0)]⊤[F ′ 2(T1)]⊤ . . . [F ′ p−1(Tp−2)]⊤[F ′ p(Tp−1)]⊤ ¯

Y X, ¯ Y → ¯ X

  • cost( ¯

X) = 4 ∗ cost(Y )

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 28 / 58

slide-29
SLIDE 29

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 29 / 58

slide-30
SLIDE 30

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 30 / 58

slide-31
SLIDE 31

Direct versus reverse AD

F : Rm → Rn

  • Direct mode

cost(J) = m ∗ 4 ∗ cost(P)

  • Reverse mode

cost(J) = n ∗ 4 ∗ cost(P)

  • Scalar output F ∈ R, n = 1

◮ Direct mode gives ∇F · ˙

X for given vector ˙ X

◮ Reverse mode gives ∇F, hence preferred

  • Vector output F ∈ Rn

◮ Direct mode gives ∇F · ˙

X for given vector ˙ X use for sensitivity equation approach

◮ Reverse mode gives (∇F)⊤ · ¯

Y use for adjoint approach

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 31 / 58

slide-32
SLIDE 32

Black-box AD

  • cost depends on alpha

call ComputeCost(alpha, cost)

  • Subroutine for cost function

subroutine ComputeCost(alpha, cost) call SolveState(alpha, u) call CostFun(alpha, u, cost) end

  • Reverse differentiation using AD

tapenade -backward \

  • head ComputeCost \
  • vars alpha \
  • outvars cost \

ComputeCost.f SolveState.f CostFun.f

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 32 / 58

slide-33
SLIDE 33

Black-box AD

  • Diffentiated subroutines:

ComputeCost b.f, SolveState b.f, CostFun b.f subroutine ComputeCost_b(alpha,alphab,cost,costb) call SolveState(alpha, u) call CostFun_b(alpha, alphab, u, ub, cost, costb) call SolveState_b(alpha, alphab, u, ub) end

  • Compute gradient using

costb = 1.0 call ComputeCost_b(alpha, alphab, cost, costb)

  • Gradient given by alphab

alphab = ∂(cost) ∂(alpha)

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 33 / 58

slide-34
SLIDE 34

AD for iterative problems

  • Solve state equation R(α, u) = 0 as steady state of

du dt + R(α, u) = 0, u(0) = uo

  • State solver

subroutine SolveState(alpha, u) u = 0.0 do call Residue(alpha, u, res) u = u - dt * res while( abs(res) > TOL) end subroutine

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 34 / 58

slide-35
SLIDE 35

AD for iterative problems

Adjoint solver, hand written subroutine SolveState_b(alpha,alphab,u,ub) resb = 0.0 do ub1 = 0.0 call Residue_bu(alpha,u,ub1,res,resb) resb = resb - (ub + ub1) while( abs(ub+ub1) .gt. 1.0e-5) call Residue_ba(alpha,alphab,u,res,resb) end subroutine resb = Adjoint variable ub = ∂J ∂u, ub1 = ∂R ∂u ⊤ v

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 35 / 58

slide-36
SLIDE 36

Adjoint iterative scheme

  • Forward iterations linearly stable

un+1 = un − ∆t R(α, un), ∆t < S(σ(R′))

  • Adjoint iteration

vn+1 = vn − ∆t

  • [R′(α, u∞)]⊤vn + ∂J

∂u

  • [R′]⊤ has same eigenvalues as R′ =

⇒ adjoint iterations stable under same condition on ∆t

  • Preconditioner for adjoint = (preconditioner for primal problem)⊤
  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 36 / 58

slide-37
SLIDE 37

AD tools

  • Source transformation (ST) and operator overloading (OO)
  • See http://www.autodiff.org

Tool Modes Method Lang ADIFOR F ST Fortran ADOLC F/B OO C TAPENADE F/B ST Fortran/C TAF/TAMC F/B ST Fortran/C MAD F OO Matlab

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 37 / 58

slide-38
SLIDE 38

Quasi 1-D flow

Inflow Outflow h(x) x

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 38 / 58

slide-39
SLIDE 39

Quasi 1-D flow

  • Quasi 1-D flow in a duct

∂ ∂t(hU) + ∂ ∂x(hf) = dh dxP, 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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 39 / 58

slide-40
SLIDE 40

Quasi 1-D flow

1 i i − 1 / 2 i + 1 / 2 N Inflow face Outflow face

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 40 / 58

slide-41
SLIDE 41

Quasi 1-D flow

  • Finite volume scheme for cell Ci = [xi−1/2, xi+1/2]

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 41 / 58

slide-42
SLIDE 42

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 42 / 58

slide-43
SLIDE 43

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 43 / 58

slide-44
SLIDE 44

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 44 / 58

slide-45
SLIDE 45

Adjoint density

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

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 45 / 58

slide-46
SLIDE 46

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 46 / 58

slide-47
SLIDE 47

Shape gradient

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

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 47 / 58

slide-48
SLIDE 48

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 48 / 58

slide-49
SLIDE 49

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 49 / 58

slide-50
SLIDE 50

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 50 / 58

slide-51
SLIDE 51

Gradient smoothing

  • Non-smooth gradients G especially in the presence of shocks
  • Smooth using an elliptic equation
  • 1 − ǫ(x) d2

dx2

  • ¯

G(x) = G(x) ǫ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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 51 / 58

slide-52
SLIDE 52

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 52 / 58

slide-53
SLIDE 53

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 53 / 58

slide-54
SLIDE 54

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 54 / 58

slide-55
SLIDE 55

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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 55 / 58

slide-56
SLIDE 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

(TIFR-CAM) Optimization BARC, 6 Oct 2010 56 / 58

slide-57
SLIDE 57

Quasi 1-D optimization: Adjoint density

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

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 57 / 58

slide-58
SLIDE 58

Example codes

  • 1-D example: nozzle flow (TAPENADE)

http://cfdlab.googlecode.com

  • 1-D example: nozzle flow (ADOLC)

http://cfdlab.googlecode.com

  • 2-D example: unstructured grid Euler solver (TAPENADE)

http://euler2d.sourceforge.net

  • 2/3-D example: structured grid Euler solver (TAPENADE)

http://nuwtun.berlios.de

  • Praveen. C

(TIFR-CAM) Optimization BARC, 6 Oct 2010 58 / 58