[PPT] - On the usage of the Riemannian geometry framework for PDE PowerPoint Presentation

SLIDE 1

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Volker Schulz

On the usage of the Riemannian geometry framework for PDE constrained shape optimization

H. Zorn, K. Welker, M. Siebenborn University of Trier

Stephan Schmidt, Universität Würzburg Nicolas Gauger, TU Kaiserslautern Caslav Ilic, German Aerospace Center, Braunschweig TEAM:

SLIDE 2

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Outline

Motivating applications
shape Newton methods
shape SQP methods

SLIDE 3

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Aerodynamic Shape Optimization

VELA: Very Efficient Large

Aircraft

Design study for blended

wing-body configurations Joint work with DLR and Airbus

SLIDE 4

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Framework

Collaborations with German Aerospace Center

in Braunschweig and Airbus Germany within the projects: MEGADESIGN, MUNA, SPP1253, ComFliTe,

DGHPOPT

Flow solvers (Flower, Tau) are matrix-free and

based on multigrid methods featuring polynomial smoothers (aka RK) [for Euler-eq. cf. van Leer/Tai/Powell AIAA CFD 1989]

SLIDE 5

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Gradient-based optimization for

compute consistent y compute gradient w.r.t. p compute new p design

SLIDE 6

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Primal-dual SQP

equivalent linear-quadratic program

  Hyy Hyp c>

y

Hpy Hpp cp> cy cp     ∆y ∆p ∆λ   =  

ryL (y, p, λ) rpL (y, p, λ)

c(y, p, λ)  

min

∆y,∆p

1 2 ✓ ∆y ∆p ◆> H ✓ ∆y ∆p ◆

+ f >

y ∆y + f > p ∆p

s.t. cy∆y + cp∆p + c = 0

SLIDE 7

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

One-shot can be interpreted as an approximate reduced SQP method of the form:

where A Cy and B is a consistent approximation of the reduced Hessian. Convergence results for quadratic models in (Ito/Kunisch/ Schulz/Gherman, SIMAX 2010) Comparison with Griewank (ISNM 165, 2014)

Can be easily generalized to finitely many state constraints (lift)
Overall computational effort is less than 10 x effort for CFD solution.

SLIDE 8

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Modular one-shot approach to

Do few iterations in the forward solver
Do few iterations in the adjoint solver, i.e.

to the equation

Use the resulting gradient for some design

step, change grid accordingly

SLIDE 9

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

2D high lift configuration

3 part wing, M=0.17146
Shape and position of flap is

to be optimized

compr. Navier Stokes solver

TAU (DLR)

Reynolds number 14.7 * 10^6
Spalart-Almaras turbulence

model

Grid with 90 000 points
10 design variables
n=20
L-BFGS for B with 3 stages

SLIDE 10

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

design variables

SLIDE 11

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

2D High Lift results

SLIDE 12

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

… zoom

SLIDE 13

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Performance

33 Iterations until convergence (criterion:

norm of the reduced gradient)

From
To
Computing time: 4.5h compared to 65min

for one simulation only (=factor 4).

SLIDE 14

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Trouble with fine parameterizations

Sensitivities: solve #p linearized problems Adjoints: solve one adjoint problem by usage of an adjoint solver A supposedly trivial matrix vector product becomes a computational bottleneck

SLIDE 15

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Finite shape parametrization in many

industrial shape optimizations

– Pro: vector space setting, fits in CAD framework – Con: complexity inevitably increases with number of parameters, mesh sensitivities can become expensive, set of reachable shapes is restricted

Nonparametric approach built on shape calculus

– Pro: avoids cons of parametric approach, can be very efficient – Con: no longer vector space setting, theoretically more challenging

Parametric versus nonparametric

SLIDE 16

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Shape gradients for free node parametrizations

SLIDE 17

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Example: simple objective, no PDE constraint

Ωt = Tt,V(Ωo) Tt,V(x) = x + t · V(x)

directional derivative in direction :

V f (Ωt) =

Z

Ωt

g(x)dx d f (Ωt)[V] := d dt t=0 f (Ωt) = d dt t=0

Z

Ωt

g(x)dx

= d

dt t=0

Z

Ωo

g(Tt,V(x))|det(DTt,V(x))|dx

SLIDE 18

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

| {z }

divV(x)

(Gauss)

=

Z

Ωo

div(g(x)V(x))dx

=

Z

∂Ωo

V(x)

~

n(x) · g(x)dx

=

Z

Ωo

d dt t=0 g(Tt,V(x))|det(DTt,V(x))|dx

=

Z

Ωo

⇥g(x)V(x) + g(x) · tr(DV(x))dx

SLIDE 19

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Shape gradients for free node parametrizations

(Hadamard)

SLIDE 20

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Important shape derivatives

d(

Z

∂Ω

g(s)ds)[W] =

Z

∂Ω

( ∂g

∂~ n(s) + κ(s)g(s))(W,~ n)ds

d(

Z

∂Ω

h(s)>n(s)ds)[W] =

Z

∂Ω

div(h(s))(W, n)ds

(cf. Delfour/Zolesio, 2001) (cf. Stephan Schmidt, 2010)

d(

Z

Ω

g(x)dx)[W] =

Z

∂Ω

g(s)(W, n)ds

SLIDE 21

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Further names to be mentioned in the field of shape calculus: Zolesio, Haslinger, Sokolowski, Pironneau, Mohammadi, Delfour, Neittanmäki, Berggren, Hintermüller, Ring, Eppler, Harbrecht, Zuazua, Sturm…

SLIDE 22

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Central observation for PDE constraints (Correa-Seeger)

Thus, we just have to build the Lagrangian and perform the (partial) shape differentiation

d f (y(Ω), Ω)[V] = dL(y, λ, Ω)[V]

SLIDE 23

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Example: drag in volume formulation (dissipation of kinetic energy into heat)

SLIDE 24

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Shape derivatives

(Schmidt/S.: Control and Cybernetics, 2010)

SLIDE 25

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Usage of shape derivatives alone may lead to unphysical

geometries

Shape Hessian approximations help to

– „smooth“ gradients – Speed up convergence in the fashion of a Newton-like method – Give potentially mesh independence

Potential problem: too much freedom

SLIDE 26

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Shape Hessians

Can become rather cumbersome and are difficult

to interprete the operators (Eppler/S./Schmidt JOTA

2008)

Fourier mode analysis gives an idea of the

general structure (S./Schmidt, SICON 2009)

d2J(u, Ω)[V1, V2]

SLIDE 27

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Performance for Navier-Stokes

12 vs. 200: 96% less iterations

speed

SLIDE 28

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Optimal non-parametric design for Euler flow in TAU (DLR)

From NACA0012 to Haack Ogive Mach 2.0 strong detached bow shock transformed to weak Drag reduction 45% pressure

SLIDE 29

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Overall performance

Wall clock time reduced by 99% (2.77h versus 100s)

SLIDE 30

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Vela aircraft

VELA: Very efficient large

aircraft

Design study for blended

wing-body configurations

115,673 surface nodes to

be optimized

Planform constant

[Schmidt/Ilic/Gauger/S. AIAA Journal 2013]

SLIDE 31

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

SLIDE 32

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

SLIDE 33

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

A closer look at the shape Hessian

Symmetry
Taylor series

⇒ Sufficient conditions ⇒ Quadratic convergence of Newton method

No

?

d(d f (Ω)[W])[V] =

Z

∂Ω( ∂g

∂~ n + κcg) hW,~ ni hV,~ ni + g hDW V,~ ni ds

SLIDE 34

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Hausdorff distance
Distance concepts

„Morphing“: Riemannian length of shortest connecting path

Ω1 Ω2

SLIDE 35

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

The Riemannian metric of Michor and Mumford (2006)

Be(S1,

2) = Emb(S1, 2)/Diff(S1)

Shape set is a manifold with tangent space

GA(h, k) =

Z

c(1 + Aκ2 c)αβds ,

A> 0

Scalar product defines a Riemannian manifold

TcBe ∼

= {h | h = α~

n, α ∈ C∞(S1,

)}

SA(h, k) =

Z

c αβ + Aα0β0ds

SLIDE 36

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

A Riemannian view on shape

ptimization

(arXiv:1203.1493)

Defining the action of a vector field as the shape

derivative, we can unleash the Riemannian structure on shape optimization

Optimization on manifolds can be performed as in [Absil

2008] for matrix manifolds.

h( f )(c) =d f (Ω)[V] , V = h~ n , c = ∂Ω

SLIDE 37

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Consequences

Symmetric Riemannian shape Hessian
Taylor series expansion

GA(Hessf (Ω)[V], W) :=GA(rVgradf (Ω), W)

=d(d f (Ω)[W])[V] d f (Ω)[rVW] f (expΩ(h)) = f (Ω) + GA(gradf (Ω), h) + 1 2GA(Hessf (Ω)h, h) + O(khk3)

SLIDE 38

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Newton convergence

Affine covariant contraction result (~Bock/Deuflhard)

Geodesic
With parallel transport
Lipschitz condition 1
Lipschitz condition 2
Contraction property as consequence

γ : [0, 1] ! M, t 7! γ(t) := expx(t∆x)

Pα,β : Tγ(α)M → Tγ(β)M

kM(y)(Pt,1J(γ(t))P0,t P0,1J(x))∆x)k  ωtk∆xk2

kM(y)P0,1(F(x) + J(x)∆x)k =: ˜

κ(x)k∆xk , ˜ κ(x)  κ < 1

k∆xk+1k  δkk∆xkk = (κ + ω

2 k∆xkk)k∆xkk

SLIDE 39

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Newton convergence

Riemannian variants of quadratic convergence

results are possible.

In particular, if we eliminate expressions from the

Hessian, which are zero at the solution, still quadratic convergence is achieved.

SLIDE 40

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Optimization algorithms

Ωk+1 = expΩk(αdk) dk = −gradf (Ωk) or dk = −Hess(Ωk)−1gradf (Ωk)

∂Ωk+1 = ∂Ωk + αdk

Retraction Steepest descent linear conv. Newton method quadratic conv.

SLIDE 41

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Back to the simple example

min

Ω

f (Ω) :=

Z

Ω

g(x)dx with g(x) =x>x 1

Riemannian shape Hessian

2

g

GA(Hessf (Ω)[V], W) = d(d f (Ω)[W])[V] d f (Ω)[rVW]

=

Z

∂Ω( ∂g

∂~ n + κc 2 g gAκ3

c

1 + Aκ2

c

)) hV,~

ni hW,~ ni ds

Z

∂Ω gAκc (hV,~

ni hW,~ ni)ττ ds

SLIDE 42

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Iteration objective line search

2.66624

3.4135 0.69 1

4.18749

7.2458 E-02 0.50 2

4.18879

6.4371 E-04 0.50 3

4.18879

7.0061 E-08

Performance of optimization algorithms

Steepest descent with exact linesearch

dA(Ωk, ˆ Ω)

ˆ Ω

Equivalent to Newton method

> example too simple

SLIDE 43

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

g(x) = x2

1 + 4x2 2

Iteration SD-objective SD-line search NM-objective NM-line search

0.5571

.9222E+00 .5100E+00

0.5571

.9222E+00 0.6400 1

0.7630

.2142E+00 .3100E+00

0.7775

.1359E+00 1.0000 2

0.7830

.5882E-01 .3700E+00

0.7854

.4265E-02 1.0000 3

0.7852

.1662E-01 .3200E+00

0.7854

.1037E-04 1.0000 4

0.7854

.4459E-02 .3500E+00

0.7854

.2590E-09 5

0.7854

.1356E-02 .3300E+00 6

0.7854

.3951E-03 .3400E+00 7

0.7854

.1218E-03 .3300E+00 8

0.7854

.3643E-04 .3500E+00 9

0.7854

.1194E-04 .3200E+00 10

0.7854

.3421E-05 .3500E+00 11

0.7854

.1137E-05 .3200E+00 12

0.7854

.3305E-06 .3500E+00 13

0.7854

.1109E-06 .3000E+00 14

0.7854

.2894E-07

dSC(Ωk, ˆ Ω)

dNM(Ωk, ˆ Ω)

SLIDE 44

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Advantages of the Riemannian view

Symmetric shape Hessian variant
Sufficient optimality conditions

– Example: is coercive

Analysis of convergence order of shape
ptimization methods

– Quadratic convergence for Newton method can be observed

Hessf ( ˆ Ω) = 2 · id

SLIDE 45

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Related approaches

Hintermüller/Ring 2003/2004: special choice of

perturbations

Hintermüller 2005: convergence from descent

property

Ring/Wirth 2012: convergence theory on

Riemannian manifolds

Sundamoorthi/Mennucci/Soatto/Yezzi 2011:

Sobolev-type metric on shape space

Harbrecht/Eppler: well-posedness studies for

star shaped domains

SLIDE 46

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Towards a Lagrange-Newton approach in shape calculus

model problem:

min

u

J(y, u) ⌘ 1 2

Z

Ω(u)(y(x) ¯

y(x))2dx + µ

Z

u 1ds

s.t. 4y(x) = f (x) , 8x 2 Ω(u) y(x) = 0 , 8x 2 ∂Ω(u)

f (x) ≡ ( f1(x) = const. , ∀x ∈ Ω1(u) f2(x) = const. , ∀x ∈ Ω2(u)

[S. /Siebenborn/Welker, 2014 arXiv:1405.3266]

SLIDE 47

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Weak formulation

Gives rise to the Lagrangian

min

u

1 2

Z

Ω(u)(y(x) ¯

y(x))2dx + µ

Z

u 1ds

s.t. au(y, p) = bu(p) , 8p 2 H1

0(Ω(u))

where au(y, p) :=

Z

Ω(u) ry(x)>rp(x)dx

bu(p) :=

Z

Ω(u) f (x)p(x)dx

L(y, u, p) := J(y, u) + au(y, p) − bu(p)

SLIDE 48

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

The adjoint problem

allows for a comfortable computation of the shape derivative (cf. Cea 1980‘s, Correa/Seeger 1985, Delfour/Zolesio 2012) Z

Ω

rz>rp dx =

Z

Ω

(y ¯

y)z dx ,

8 ¯

y 2 H1

0(Ω(u))

dJ(y(u), u)[V] = dL(y, u, p)[V]

SLIDE 49

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Resulting KKT-conditions

(adjoint) (design) (state) What about applying Newton‘s method to these equations?

or even multigrid?

But the function spaces for y and p vary with the shape!

∂ ∂yL(y, u, p) = 0 , ∀y dL(y, u, p)[V] = 0 , ∀V ∂ ∂pL(y, u, ¯ p) = 0 , ∀p

SLIDE 50

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Remedy: yet another structure from differential geometry: vector bundles

u := Γ ∈ N := B0

e ([0, 1], 2) := Emb0([0, 1], 2)/Diff([0, 1])

SLIDE 51

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

E := {(H(u), u) | u ∈ N } τi : π−1(Ui) → H0 × Ui τi(u) : π−1(u) = H(u) → H0

Vector bundle local diffeomorphisms: domain deformation

SLIDE 52

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Variational Newton step h

G(HessL(y, u, p)h, k) = G(rL(y, u, p), k) , 8k 2 TE

After omitting terms in the Hessian, which are zero at the solution, this can be identified with the linear-quadratic

ptimal control problem:

min

(z,w)

Z

Ω(u)

z2 2 + (y ¯ y)zdx + µ

Z

u κ1wds + 1

2

Z

u ( f2 f1) ∂p

∂n1 w2 + µ ✓∂w ∂τ ◆2 ds s.t.

Z

Ω(u) rz>r ¯

qdx +

Z

u ( f2 f1) ¯

qwds =

Z

Ω(u) ry>r ¯

qdx +

Z

Ω(u) f ¯

qdx ,

8 ¯

q 2 H1

0(Ω(u))

> can be efficiently solved by multigrid methods (e.g. Borzi/
S. 2009/2012, or Ascher, or Zulehner, or Nash)

SLIDE 53

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Strong form of QP

  Hyy Hyu c>

y

Huy Huu c>

u

cy cu     hy hu hp   =  

ryL (y, u, p) ruL (y, u, p)

c(y, u, λ)  

min

hy,hu

1 2 ✓ hy hu ◆> H ✓ hy hu ◆

+ f >

y hy + f > p hu

s.t. cyhy + cuhu + c = 0

min

(z,w)

Z

Ω(u)

z2 2 + (y ¯ y)zdx+µ

Z

u κ1wds + 1

2

Z

u ( f2 f1) ∂p

∂n1 w2 + µ ✓∂w ∂τ ◆2 ds s.t. 4z = 4y + f1 in Ω1(u)

4z = 4y + f2

in Ω2(u) ∂z ∂n1

= f1w

n u

∂z ∂n2

= f2w

n u

z = 0

n ∂Ω(u)

SLIDE 54

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Numerical results

[0] [1] [2]

SLIDE 55

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Quadratic convergence?

Distance to solution on several grids
Major advantages over steepest descent:

– natural step scaling of 1 (vs. 10000) – Optimal control technology (multigrid!) can be exploited

It.-No. Ω1

h

Ω2

h

Ω3

h

0.0705945 0.070637 0.0706476 1 0.0043115 0.004104 0.0040465 2 0.0003941 0.000104 0.0000645

SLIDE 56

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Parabolic interface problem

SLIDE 57

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

s.t.

SLIDE 58

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Pantz (2005) notes that Cea has to be

modified here

Very similar to results by Ito/Kunisch 2008

and Paganini 2014 for the elliptic case

Shape derivative

SLIDE 59

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Adaptive meshes

SLIDE 60

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Riemannian Quasi-Newton

SLIDE 61

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Riemannian limited BFGS

SLIDE 62

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

L-BFGS

2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40 0.0001 0.001 0.01 0.1 1 full BFGS l5-BFGS l2-BFGS gradient 3 6 9 12 15 18 21 24 27 30 33 36 39 42 45 48 0.0001 0.001 0.01 0.1 1 BFGS 100k BFGS 3k gradient 100k gradient 3k

S./Siebenborn/Welker, 2014, arXiv:1409.3464 3D Scalability study Up to 320M el. on 64k processors (subm. CVS 2015)

SLIDE 63

Volker Schulz

ESI w ESI wor

rkshop

shop, February 16, 2015

Conclusions

Many applications profit from a non-parametric

approach.

The Riemannian perspective enables SQP-like shape
ptimization methods.
Several open questions remain

– what is the appropriate Riemannian metric? – C-infinity is far too smooth for practical applications – What about completions in an appropriate norm? – What features of the exponential map are essential for approximations (retractions)? – Large shape deformations tend to develop kinks