[PPT] - Optimization-based computational physics and high-order methods: PowerPoint Presentation

SLIDE 1

Optimization-based computational physics and high-order methods: from optimized analysis to design and data assimilation

Matthew J. Zahr† and Per-Olof Persson

CRD Postdoc Seminar Series Lawrence Berkeley National Laboratory, Berkeley, CA September 18, 2017

† Luis W. Alvarez Postdoctoral Fellow

Department of Mathematics Lawrence Berkeley National Laboratory 1 / 63

SLIDE 2

PDE optimization is ubiquitous in science and engineering

Design: Find system that optimizes performance metric, satisfies constraints Aerodynamic shape design of automobile Optimal flapping motion of micro aerial vehicle Control: Drive system to a desired state

2 / 63

SLIDE 3

PDE optimization is ubiquitous in science and engineering

Inverse problems: Infer the problem setup given solution observations Left: Material inversion – find inclusions from acoustic, structural measurements Right: Source inversion – find source of airborne contaminant from downstream measurements Full waveform inversion – estimate subsurface of Earth’s crust from acoustic measurements

3 / 63

SLIDE 4

Unsteady PDE-constrained optimization formulation

Goal: Find the solution of the unsteady PDE-constrained optimization problem minimize

U, µ

J (U, µ) subject to C(U, µ) ≤ 0 ∂U ∂t + ∇ · F (U, ∇U) = 0 in v(µ, t) where

U(x, t)

PDE solution

µ

design/control parameters

J (U, µ) =

Tf

T0

Γ

j(U, µ, t) dS dt

bjective function
C(U, µ) =

Tf

T0

Γ

c(U, µ, t) dS dt constraints

4 / 63

SLIDE 5

Nested approach to PDE-constrained optimization

Primal PDE Dual PDE Optimizer 5 / 63

SLIDE 6

Nested approach to PDE-constrained optimization

Primal PDE Dual PDE Optimizer µ 5 / 63

SLIDE 7

Nested approach to PDE-constrained optimization

Primal PDE Dual PDE Optimizer J (U, µ) 5 / 63

SLIDE 8

Nested approach to PDE-constrained optimization

Primal PDE Dual PDE Optimizer J (U, µ) µ U 5 / 63

SLIDE 9

Nested approach to PDE-constrained optimization

Primal PDE Dual PDE Optimizer J (U, µ)

dJ dµ (U, µ)

5 / 63

SLIDE 10

Highlights of globally high-order discretization

Arbitrary Lagrangian-Eulerian formulation:

Map, G(·, µ, t), from physical v(µ, t) to reference V ∂U X ∂t

X

+ ∇X · F X(U X, ∇XU X) = 0

Space discretization: discontinuous Galerkin

M ∂u ∂t = r(u, µ, t)

Time discretization: diagonally implicit RK

un = un−1 +

s

i=1

bikn,i Mkn,i = ∆tnr (un,i, µ, tn,i)

Quantity of interest: solver-consistency

F(u0, . . . , uNt, k1,1, . . . , kNt,s)

X1 X2 NdA V x1 x2 nda v G, g, vX

Mapping-Based ALE DG Discretization c1 a11 c2 a21 a22 . . . . . . . . . ... cs as1 as2 · · · ass b1 b2 · · · bs Butcher Tableau for DIRK

6 / 63

SLIDE 11

Adjoint method to efficiently compute gradients of QoI

Consider the fully discrete output functional F(un, kn,i, µ)
Represents either the objective function or a constraint
The total derivative with respect to the parameters µ, required in the context
f gradient-based optimization, takes the form

dF dµ = ∂F ∂µ +

Nt

n=0

∂F ∂un ∂un ∂µ +

Nt

n=1

s

i=1

∂F ∂kn,i ∂kn,i ∂µ

The sensitivities, ∂un

∂µ and ∂kn,i ∂µ , are expensive to compute, requiring the solution of nµ linear evolution equations

Adjoint method: alternative method for computing dF

dµ that require one linear evolution equation for each quantity of interest, F

7 / 63

SLIDE 12

Dissection of fully discrete adjoint equations

Linear evolution equations solved backward in time
Primal state/stage, un,i required at each state/stage of dual problem
Heavily dependent on chosen ouput

λNt = ∂F ∂uNt

T

λn−1 = λn + ∂F ∂un−1

T

+

s

i=1

∆tn ∂r ∂u (un,i, µ, tn−1 + ci∆tn)T κn,i M T κn,i = ∂F ∂uNt

T

+ biλn +

s

j=i

aji∆tn ∂r ∂u (un,j, µ, tn−1 + cj∆tn)T κn,j

Gradient reconstruction via dual variables

dF dµ = ∂F ∂µ + λ0

T ∂g

∂µ(µ) +

Nt

n=1

∆tn

s

i=1

κn,i

T ∂r

∂µ(un,i, µ, tn,i) [Zahr and Persson, 2016]

8 / 63

SLIDE 13

Optimal control, time-morphed geometry

Optimal Rigid Body Motion (RBM) and Time-Morphed Geometry (TMG), thrust = 2.5 Energy = 9.4096 Thrust = 0.1766 Energy = 4.9476 Thrust = 2.500 Energy = 4.6182 Thrust = 2.500 Initial Guess Optimal RBM Tx = 2.5 Optimal RBM/TMG Tx = 2.5

9 / 63

SLIDE 14

Energetically optimal flapping in three-dimensions

Energy = 1.4459e-01 Thrust = -1.1192e-01 Energy = 3.1378e-01 Thrust = 0.0000e+00

10 / 63

SLIDE 15

Extension: Parametrized time domain [Wang et al., 2017]

Parametrization of time domain, e.g., flapping frequency, leads to

parametrization of time discretization in fully discrete setting T(µ) = Nt∆t = ⇒ Nt = Nt(µ) or ∆t = ∆t(µ)

Choose ∆t = ∆t(µ) to avoid discrete changes
Does not change adjoint equations themselves, only reconstruction of gradient

from adjoint solution

11 / 63

SLIDE 16

Energetically optimal flapping vs. required thrust

Energy = 1.8445 Thrust = 0.06729 Energy = 0.21934 Thrust = 0.0000 Energy = 6.2869 Thrust = 2.5000 Initial Guess Optimal Tx = 0 Optimal Tx = 2.5

12 / 63

SLIDE 17

Super-resolution MR images through optimization

Space-time MRI data: noisy, low-resolution

In collaboration with research team at Lund University, working to use our

high-order optimization framework to generate “super-resolved” MR images

Idea: Match CFD parameters (material properties, boundary conditions) to

MRI data using optimization

Goal: visualize high-resolution flow and accurately compute quantities of

interest, i.e., wall shear stress

13 / 63

SLIDE 18

Phase I: synthetic data

3 14 2 6 Geometry and boundary conditions for synthetic MRI data assimilation setting. Boundary conditions: viscous wall ( ), parametrized inflow ( ), and outflow ( ). MRI data collected in the red shaded region.

14 / 63

SLIDE 19

Coarse MRI grid (24 × 36), 10 time samples, 3% noise

Synthetic MRI data d∗

i,n (top) and

computational representation of MRI data di,n (bottom) Reconstructed flow

15 / 63

SLIDE 20

Fine MRI grid (40 × 60), 20 time samples, 10% noise

Synthetic MRI data d∗

i,n (top) and

computational representation of MRI data di,n (bottom) Reconstructed flow

16 / 63

SLIDE 21

Stochastic PDE-constrained optimization formulation

minimize

µ∈Rnµ

E[J (u, µ, · )] subject to r(u; µ, ξ) = 0 ∀ξ ∈ Ξ

r : Rnu × Rnµ × Rnξ → Rnu

discretized stochastic PDE

J : Rnu × Rnµ × Rnξ → R

quantity of interest

u ∈ Rnu

PDE state vector

µ ∈ Rnµ

(deterministic) optimization parameters

ξ ∈ Rnξ

stochastic parameters

E[F] ≡
Ξ

F(ξ)ρ(ξ) dξ

17 / 63

SLIDE 22

Nested approach to stochastic PDE-constrained optimization

Ensemble of primal/dual PDE solves required at every optimization iteration

Primal PDE Dual PDE Optimizer 18 / 63

SLIDE 23

Nested approach to stochastic PDE-constrained optimization

Ensemble of primal/dual PDE solves required at every optimization iteration

Primal PDE Dual PDE Optimizer µ 18 / 63

SLIDE 24

Nested approach to stochastic PDE-constrained optimization

Ensemble of primal/dual PDE solves required at every optimization iteration

Primal PDE Dual PDE Optimizer E [J (u, µ, · )] 18 / 63

SLIDE 25

Nested approach to stochastic PDE-constrained optimization

Ensemble of primal/dual PDE solves required at every optimization iteration

Primal PDE Dual PDE Optimizer E [J (u, µ, · )] µ u 18 / 63

SLIDE 26

Nested approach to stochastic PDE-constrained optimization

Ensemble of primal/dual PDE solves required at every optimization iteration

Primal PDE Dual PDE Optimizer E [J (u, µ, · )] E

dJ

dµ (u, µ, · )

18 / 63

SLIDE 27

Proposed approach: managed inexactness

Replace expensive PDE with inexpensive approximation model

Anisotropic sparse grids used for inexact integration of risk measures
Reduced-order models used for inexact PDE evaluations

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

m(µ)

19 / 63

SLIDE 28

Proposed approach: managed inexactness

Replace expensive PDE with inexpensive approximation model

Anisotropic sparse grids used for inexact integration of risk measures
Reduced-order models used for inexact PDE evaluations

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

m(µ)

HDM 19 / 63

SLIDE 29

Proposed approach: managed inexactness

Replace expensive PDE with inexpensive approximation model

Anisotropic sparse grids used for inexact integration of risk measures
Reduced-order models used for inexact PDE evaluations

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

m(µ)

HDM HDM 19 / 63

SLIDE 30

Proposed approach: managed inexactness

Replace expensive PDE with inexpensive approximation model

Anisotropic sparse grids used for inexact integration of risk measures
Reduced-order models used for inexact PDE evaluations

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

m(µ)

HDM HDM ROM 19 / 63

SLIDE 31

Proposed approach: managed inexactness

Replace expensive PDE with inexpensive approximation model

Anisotropic sparse grids used for inexact integration of risk measures
Reduced-order models used for inexact PDE evaluations

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

m(µ) Manage inexactness with trust region method

Embedded in globally convergent trust region method
Error indicators1 to account for all sources of inexactness
Refinement of approximation model using greedy algorithms

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

mk(µ) subject to ||µ − µk|| ≤ ∆k

1Must be computable and apply to general, nonlinear PDEs

20 / 63

SLIDE 32

First source of inexactness: anisotropic sparse grids

Stochastic collocation using anisotropic sparse grid nodes to approximate integral with summation minimize

u∈Rnu, µ∈Rnµ

E[J (u, µ, · )] subject to r(u, µ, ξ) = 0 ∀ξ ∈ Ξ

⇓

minimize

u∈Rnu, µ∈Rnµ

EI[J (u, µ, · )] subject to r(u, µ, ξ) = 0 ∀ξ ∈ ΞI [Kouri et al., 2013, Kouri et al., 2014]

21 / 63

SLIDE 33

Source of inexactness: anisotropic sparse grids

−1 −0.5 0 0.5 1 −1 −0.5 0.5 1 ξ1 ξ2 Quad rule ξ1 ⊗ ξ2 1 2 3 4 5 6 1 2 3 4 5 6 i1 i2 Index set Index set (I) – Neighbors (N(I)) –

22 / 63

SLIDE 34

Source of inexactness: anisotropic sparse grids

−1 −0.5 0 0.5 1 −1 −0.5 0.5 1 ξ1 ξ2 Quad rule ξ1 ⊗ ξ2 1 2 3 4 5 6 1 2 3 4 5 6 i1 i2 Index set Index set (I) – Neighbors (N(I)) –

23 / 63

SLIDE 35

Second source of inexactness: reduced-order models

Stochastic collocation of the reduced-order model over anisotropic sparse grid nodes used to approximate integral with cheap summation minimize

u∈Rnu, µ∈Rnµ

E[J (u, µ, · )] subject to r(u, µ, ξ) = 0 ∀ξ ∈ Ξ

⇓

minimize

u∈Rnu, µ∈Rnµ

EI[J (u, µ, · )] subject to r(u, µ, ξ) = 0 ∀ξ ∈ ΞI

⇓

minimize

ur∈Rku, µ∈Rnµ

EI[J (Φur, µ, · )] subject to ΦT r(Φur, µ, ξ) = 0 ∀ξ ∈ ΞI

24 / 63

SLIDE 36

Source of inexactness: projection-based model reduction

Model reduction ansatz: state vector lies in low-dimensional subspace

u ≈ Φur

Φ =
φ1

· · · φku

∈ Rnu×ku is the reduced (trial) basis (nu ≫ ku)
ur ∈ Rku are the reduced coordinates of u
Substitute into r(u, µ) = 0 and perform Galerkin projection

ΦT r(Φur, µ) = 0

25 / 63

SLIDE 37

Few global, data-driven basis functions v. many local ones

Instead of using traditional local

shape functions (e.g., FEM), use global shape functions

Instead of a-priori, analytical shape

functions, leverage data-rich computing environment by using data-driven modes

26 / 63

SLIDE 38

Proposed approach: managed inexactness

Replace expensive PDE with inexpensive approximation model

Anisotropic sparse grids used for inexact integration of risk measures
Reduced-order models used for inexact PDE evaluations

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

m(µ) Manage inexactness with trust region method

Embedded in globally convergent trust region method
Error indicators2 to account for all sources of inexactness
Refinement of approximation model using greedy algorithms

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

mk(µ) subject to ||µ − µk|| ≤ ∆k

2Must be computable and apply to general, nonlinear PDEs

27 / 63

SLIDE 39

Trust region ingredients for global convergence

Approximation models mk(µ) Error indicators ||∇F(µ) − ∇mk(µ)|| ≤ ξϕk(µ) ξ > 0 Adaptivity ϕk(µk) ≤ κϕ min{||∇mk(µk)|| , ∆k} Global convergence lim inf

k→∞

||∇F(µk)|| = 0

28 / 63

SLIDE 40

Trust region method: ROM/SG approximation model

Approximation models built on two sources of inexactness mk(µ) = EIk [J (Φkur(µ, · ), µ, · )] Error indicators that account for both sources of error ϕk(µ) = α1E1(µ; Ik, Φk) + α2E2(µ; Ik, Φk) + α3E4(µ; Ik, Φk) Reduced-order model errors E1(µ; I, Φ) = EI ∪ N (I) [||r(Φur(µ, ·), µ, · )||] E2(µ; I, Φ) = EI ∪ N (I)

rλ(Φur(µ, ·), Φλr(µ, · ), µ, · )
Sparse grid truncation errors

E4(µ; I, Φ) = EN (I) [||∇J (Φur(µ, · ), µ, · )||]

29 / 63

SLIDE 41

Adaptivity: Dimension-adaptive greedy method

while E4(Φ, I, µk) > κϕ 3α3 min{||∇mk(µk)|| , ∆k} do Refine index set: Dimension-adaptive sparse grids Ik ← Ik ∪ {j∗} where j∗ = arg max

j∈N (Ik)

Ej [||∇J (Φur(µ, · ), µ, · )||]

30 / 63

SLIDE 42

Adaptivity: Dimension-adaptive greedy method

while E4(Φ, I, µk) > κϕ 3α3 min{||∇mk(µk)|| , ∆k} do Refine index set: Dimension-adaptive sparse grids Ik ← Ik ∪ {j∗} where j∗ = arg max

j∈N (Ik)

Ej [||∇J (Φur(µ, · ), µ, · )||] Refine reduced-order basis: Greedy sampling while E1(Φ, I, µk) > κϕ 3α1 min{||∇mk(µk)|| , ∆k} do Φk ←

Φk

u(µk, ξ∗) λ(µk, ξ∗)

ξ∗ = arg max

ξ∈Ξj∗

ρ(ξ) ||r(Φkur(µk, ξ), µk, ξ)|| end while

30 / 63

SLIDE 43

Adaptivity: Dimension-adaptive greedy method

while E4(Φ, I, µk) > κϕ 3α3 min{||∇mk(µk)|| , ∆k} do Refine index set: Dimension-adaptive sparse grids Ik ← Ik ∪ {j∗} where j∗ = arg max

j∈N (Ik)

Ej [||∇J (Φur(µ, · ), µ, · )||] Refine reduced-order basis: Greedy sampling while E1(Φ, I, µk) > κϕ 3α1 min{||∇mk(µk)|| , ∆k} do Φk ←

Φk

u(µk, ξ∗) λ(µk, ξ∗)

ξ∗ = arg max

ξ∈Ξj∗

ρ(ξ) ||r(Φkur(µk, ξ), µk, ξ)|| end while while E2(Φ, I, µk) > κϕ 3α2 min{||∇mk(µk)|| , ∆k} do Φk ←

Φk

u(µk, ξ∗) λ(µk, ξ∗)

ξ∗ = arg max

ξ∈Ξj∗

ρ(ξ)

rλ(Φkur(µk, ξ), Φkλr(µk, ξ), µk, ξ)
end while

end while

30 / 63

SLIDE 44

Optimal control of steady Burgers’ equation

Optimization problem:

minimize

µ∈Rnµ

Ξ

ρ(ξ) 1 1 2(u(µ, ξ, x) − ¯ u(x))2 dx + α 2 1 z(µ, x)2 dx

dξ

where u(µ, ξ, x) solves −ν(ξ)∂xxu(µ, ξ, x) + u(µ, ξ, x)∂xu(µ, ξ, x) = z(µ, x) x ∈ (0, 1), ξ ∈ Ξ u(µ, ξ, 0) = d0(ξ) u(µ, ξ, 1) = d1(ξ)

Target state: ¯

u(x) ≡ 1

Stochastic Space: Ξ = [−1, 1]3, ρ(ξ)dξ = 2−3dξ

ν(ξ) = 10ξ1−2 d0(ξ) = 1 + ξ2 1000 d1(ξ) = ξ3 1000

Parametrization: z(µ, x) – cubic splines with 51 knots, nµ = 53

31 / 63

SLIDE 45

Optimal control and statistics

0.2 0.4 0.6 0.8 1 2 x z(µ, x) 0.2 0.4 0.6 0.8 1 1 x E[u(µ, · , x)] Optimal control and corresponding mean state ( ) ± one ( ) and two ( ) standard deviations

32 / 63

SLIDE 46

Global convergence without pointwise agreement

1 2 3 4 10−8 10−5 10−2 Major iteration ( ) |F(µk) − F(µ∗)| ( ) |F(ˆ µk) − F(µ∗)| ( ) |mk(µk) − F(µ∗)| ( ) |mk(ˆ µk) − F(µ∗)|

F(µk) mk(µk) F(ˆ µk) mk(ˆ µk) ||∇F(µk)|| ρk Success? 6.6506e-02 7.2694e-02 5.3655e-02 5.9922e-02 2.2959e-02 1.0257e+00 1.0000e+00 5.3655e-02 5.9593e-02 5.0783e-02 5.7152e-02 2.3424e-03 9.7512e-01 1.0000e+00 5.0783e-02 5.0670e-02 5.0412e-02 5.0292e-02 1.9724e-03 9.8351e-01 1.0000e+00 5.0412e-02 5.0292e-02 5.0405e-02 5.0284e-02 9.2654e-05 8.7479e-01 1.0000e+00 5.0405e-02 5.0404e-02 5.0403e-02 5.0401e-02 8.3139e-05 9.9946e-01 1.0000e+00 5.0403e-02 5.0401e-02

2.2846e-06
Convergence history of trust region method built on two-level approximation

33 / 63

SLIDE 47

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−5 10−4 10−3 10−2 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

34 / 63

SLIDE 48

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−5 10−4 10−3 10−2 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

34 / 63

SLIDE 49

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−5 10−4 10−3 10−2 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

34 / 63

SLIDE 50

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−5 10−4 10−3 10−2 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

34 / 63

SLIDE 51

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−5 10−4 10−3 10−2 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

34 / 63

SLIDE 52

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−5 10−4 10−3 10−2 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

34 / 63

SLIDE 53

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis

35 / 63

SLIDE 54

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis

35 / 63

SLIDE 55

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis

35 / 63

SLIDE 56

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis

35 / 63

SLIDE 57

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations

35 / 63

SLIDE 58

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations

35 / 63

SLIDE 59

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations

35 / 63

SLIDE 60

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations

35 / 63

SLIDE 61

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations Proposed solution align features of solution basis with features in the solution using

ptimization formulation and solver

35 / 63

SLIDE 62

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations Proposed solution align features of solution basis with features in the solution using

ptimization formulation and solver

35 / 63

SLIDE 63

Optimization beyond design/control: high-order shock resolution

Fundamental issue: interpolate discontinuity with polynomial basis Exising solutions: limiting, adaptive refinement, artificial viscosity usually result in order reduction or very fine discretizations Proposed solution align features of solution basis with features in the solution using

ptimization formulation and solver

35 / 63

SLIDE 64

Shock tracking optimization formulation

Consider the spatial discretization of the conservation law

∇X · F (U; X) = 0 → r(u; x) = 0

U, X are the continuous state vector and coordinate
x contains the coordinates of the mesh nodes
u contains the discrete state vector corresponding to U at the mesh nodes
Shock tracking formulation

minimize

u, x

f(u, x) subject to r(u; x) = 0 Key assumption: Solution basis supports discontinuties along element edges, i.e., discontinuous Galerkin, finite volume

36 / 63

SLIDE 65

Shock tracking objective function

Requirements on objective function

btains minimum when mesh

edge aligned with shock and monotonically decreases to minimum in (large) neighborhood f(u; x) = fshk(u; x) + αfmsh(x) fshk(u, x) =

ne

e=1
Ωe(x)

|u − ¯ u|2 dV fmsh(x) =

ne

e=1

ne

q

k=1
tr GT G

det G

0.64

0.65 0.66 0.67 2 4 6 8 ·10−4 position of mesh edge closest to shock

Objective function as an element edge is smoothly swept across shock location for: fshk(u, x) ( ), residual-based objective ( ), and Rankine-Hugniot-based objective ( ).

37 / 63

SLIDE 66

Full space optimization solver for shock tracking

Cannot use nested approach to PDE optimization because it requires solving r(u; x) = 0 for x = x∗ = ⇒ crash

Full space approach: u → u∗ and x → x∗ simultaneously
Define Lagrangian

L(u, x, λ) = f(u; x) − λT r(u; x)

First-order optimality (KKT) conditions for full space optimization problem

∇uL(u∗, x∗, λ∗) = 0, ∇xL(u∗, x∗, λ∗) = 0, ∇λL(u∗, x∗, λ∗) = 0

Apply (quasi-)Newton method3 to solve nonlinear KKT system for u∗, x∗, λ∗

3usually requires globalization such as linesearch or trust-region

38 / 63

SLIDE 67

Nozzle flow: quasi-1d Euler equations

0.5 0.5 1 0.6 A(x) Geometry and boundary conditions for nozzle flow. Boundary conditions: inviscid wall ( ), inflow ( ), outflow ( ).

39 / 63

SLIDE 68

Resolution of 1d transonic flow with only 4 quartic elements

0.2 0.4 0.6 0.8 1 0.5 1 x ρ(x) 0.2 0.4 0.6 0.8 1 x The solution of the quasi-1d Euler equations using: 300 linear elements ( ) and 4 quartic elements ( ) on a mesh not aligned (left) and aligned (right) with the shock.

40 / 63

SLIDE 69

Supersonic flow around cylinder: 2D Euler equations

5 4 1 Geometry and boundary conditions for supersonic flow around cylinder. Boundary conditions: inviscid wall/symmetry condition ( ) and farfield ( ).

41 / 63

SLIDE 70

Resolution of 2D supersonic flow with only 67 quadratic elements

The solution of the 2d Euler equations using: 67 quadratic elements on a mesh not aligned with the shock (left), 67 linear elements on a mesh aligned with the shock (middle), 67 quadratic elements on a mesh aligned with the shock (right).

42 / 63

SLIDE 71

Convergence to optimal solution and mesh

43 / 63

SLIDE 72

PDE-constrained optimization for design/control and beyond

Globally high-order numerical method and

adjoint-based gradient computations for efficient design and data assimilation

energetically optimal flapping, energy harvesting

mechanisms, super-resolution MRI

Globally convergent multifidelity framework for

PDE-constrained optimization under uncertainty

risk-averse flow control
Optimization-based shock tracking

framework for highly resolved supersonic flows

n extremely coarse meshes

44 / 63

SLIDE 73

References I

Kouri, D. P., Heinkenschloss, M., Ridzal, D., and van Bloemen Waanders,

B. G. (2013).

A trust-region algorithm with adaptive stochastic collocation for PDE optimization under uncertainty. SIAM Journal on Scientific Computing, 35(4):A1847–A1879. Kouri, D. P., Heinkenschloss, M., Ridzal, D., and van Bloemen Waanders,

B. G. (2014).

Inexact objective function evaluations in a trust-region algorithm for PDE-constrained optimization under uncertainty. SIAM Journal on Scientific Computing, 36(6):A3011–A3029. Wang, J., Zahr, M. J., and Persson, P.-O. (6/5/2017 – 6/9/2017). Energetically optimal flapping flight based on a fully discrete adjoint method with explicit treatment of flapping frequency. In Proc. of the 23rd AIAA Computational Fluid Dynamics Conference, Denver, Colorado. American Institute of Aeronautics and Astronautics.

45 / 63

SLIDE 74

References II

Zahr, M. J., Huang, Z., and Persson, P.-O. (1/8/2018 – 1/12/2018). Adjoint-based optimization of time-dependent fluid-structure systems using a high-order discontinuous Galerkin discretization. In AIAA Science and Technology Forum and Exposition (SciTech2018), Kissimmee, Florida. American Institute of Aeronautics and Astronautics. Zahr, M. J. and Persson, P.-O. (2016). An adjoint method for a high-order discretization of deforming domain conservation laws for optimization of flow problems. Journal of Computational Physics. Zahr, M. J., Persson, P.-O., and Wilkening, J. (2016). A fully discrete adjoint method for optimization of flow problems on deforming domains with time-periodicity constraints. Computers & Fluids.

46 / 63

SLIDE 75

PDE optimization is ubiquitous in science and engineering

Control: Drive system to a desired state Boundary flow control Metamaterial cloaking – electromagnetic invisibility

47 / 63

SLIDE 76

High-order discretization of PDE-constrained optimization

Continuous PDE-constrained optimization problem

minimize

U, µ

J (U, µ) subject to C(U, µ) ≤ 0 ∂U ∂t + ∇ · F (U, ∇U) = 0 in v(µ, t)

Fully discrete PDE-constrained optimization problem

minimize

u0, ..., uNt∈RNu, k1,1, ..., kNt,s∈RNu, µ∈Rnµ

J(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) subject to C(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) ≤ 0 u0 − g(µ) = 0 un − un−1 −

s

i=1

bikn,i = 0 Mkn,i − ∆tnr (un,i, µ, tn,i) = 0

48 / 63

SLIDE 77

Adjoint equation derivation: outline

Define auxiliary PDE-constrained optimization problem

minimize

u0, ..., uNt∈RNu, k1,1, ..., kNt,s∈RNu

F(u0, . . . , uNt, k1,1, . . . , kNt,s, µ) subject to R0 = u0 − g(µ) = 0 Rn = un − un−1 −

s

i=1

bikn,i = 0 Rn,i = Mkn,i − ∆tnr (un,i, µ, tn,i) = 0

Define Lagrangian

L(un, kn,i, λn, κn,i) = F − λ0

T R0 − Nt

n=1

λn

T Rn − Nt

n=1

s

i=1

κn,i

T Rn,i

The solution of the optimization problem is given by the

Karush-Kuhn-Tucker (KKT) sytem ∂L ∂un = 0, ∂L ∂kn,i = 0, ∂L ∂λn = 0, ∂L ∂κn,i = 0

49 / 63

SLIDE 78

Extension: constraint requiring time-periodicity [Zahr et al., 2016]

Optimization of cyclic problems requires finding time-periodic solution of PDE
Necessary for physical relevance and avoid transients that may lead to crash

minimize

U, µ

F(U, µ) subject to U(x, 0) = U(x, T) ∂U ∂t + ∇ · F (U, ∇U) = 0 λNt = λ0 + ∂F ∂uNt

T

λn−1 = λn + ∂F ∂un−1

T

+

s

i=1

∆tn ∂rn,i ∂u

T

κn,i M T κn,i = ∂F ∂uNt

T

+ biλn +

s

j=i

aji∆tn ∂rn,i ∂u

T

κn,j

2 4 −60 −40 −20 time power 2 4 −4 −2 time power Time history of power on airfoil of flow initialized from steady-state ( ) and from a time-periodic solution ( )

50 / 63

SLIDE 79

Energetically optimal flapping vs. required thrust: QoI

0.5 1 1.5 2 2.5 2 4 6 W∗ 0.5 1 1.5 2 2.5 0.1 0.15 0.2 0.25 f ∗ 0.5 1 1.5 2 2.5 1.2 1.4 1.6 1.8 2 ¯ Tx y∗

max

0.5 1 1.5 2 2.5 20 40 60 ¯ Tx θ∗

max

The optimal flapping energy (W ∗), frequency (f ∗), maximum heaving amplitude (y∗

max),

and maximum pitching amplitude (θ∗

max) as a function of the thrust constraint ¯

Tx.

51 / 63

SLIDE 80

Extension: Multiphysics problems [Zahr et al., 2018]

For problems that involve the interaction of multiple types of physical

phenomena, no changes required if monolithic system considered M 0 ˙ u0 = r0(u0, c0(u0, u1)) M 1 ˙ u1 = r1(u1, c1(u0, u1))

However, to solve in partitioned manner and achieve high-order, split as

follows and apply implicit-explicit Runge-Kutta M 0 ˙ u0 = r0(u0, c0(u0, u1)) M 1 ˙ u1 = r1(u1, ˜ c1) + (r1(u1, c1(u0, u1)) − r1(u1, ˜ c1))

Adjoint equations inherit explicit-implicit structure

52 / 63

SLIDE 81

Optimal energy harvesting from foil-damper system

Goal: Maximize energy harvested from foil-damper system maximize

µ

1 T T (c˙ h2(us) − Mz(uf) ˙ θ(µ, t)) dt

Fluid: Isentropic Navier-Stokes on deforming domain (ALE)
Structure: Force balance in y-direction between foil and damper
Motion driven by imposed θ(µ, t) = µ1 cos(2πft)

c θ(µ, t) h(us) µ∗

1 ≈ 45◦ 53 / 63

SLIDE 82

MRI data assimilation formulation

d∗

i,n : MRI measurement taken in voxel i at the nth time sample

di,n(U, µ): computational representation of d∗

i,n

di,n(U, µ) = T

V

wi,n(x, t) · U(x, t) dV dt wi,n(x, t) = χs(x; xi, ∆x)χt(t; tn, ∆t) χt(s; c, w) = 1 1 + e−(s−(c−0.5w))/σ − 1 1 + e−(s−(c+0.5w))/σ χs(x; c, w) = χt(x1; c1, w1)χt(x2; c2, w2)χt(x3; c3, w3)

xi - center of ith MRI voxel
tn time instance of n MRI sample
∆x - size of MRI voxel in each dimension
∆t sampling interval in time

minimize

U, µ nxyz

i=1

nt

n=1

αi,n 2

di,n(U, µ) − d∗

i,n

2

2 54 / 63

SLIDE 83

Coarse MRI grid (24 × 36), 20 time samples, 3% noise

Synthetic MRI data d∗

i,n (top) and

computational representation of MRI data di,n (bottom) Reconstructed flow

55 / 63

SLIDE 84

Fine MRI grid (40 × 60), 20 time samples, 3% noise

Synthetic MRI data d∗

i,n (top) and

computational representation of MRI data di,n (bottom) Reconstructed flow

56 / 63

SLIDE 85

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 86

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 87

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 88

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 89

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 90

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 91

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 92

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 93

Trust region framework for optimization with ROMs

µ-space

57 / 63

SLIDE 94

Trust region ingredients for global convergence

minimize

µ∈Rnµ

F(µ) − → minimize

µ∈Rnµ

mk(µ) subject to ||µ − µk|| ≤ ∆k Approximation models mk(µ), ψk(µ) Error indicators ||∇F(µ) − ∇mk(µ)|| ≤ ξϕk(µ) ξ > 0 |F(µk) − F(µ) + ψk(µ) − ψk(µk)| ≤ σθk(µ) σ > 0 Adaptivity ϕk(µk) ≤ κϕ min{||∇mk(µk)|| , ∆k} θk(ˆ µk)ω ≤ η min{mk(µk) − mk(ˆ µk), rk}

58 / 63

SLIDE 95

Trust region method with inexact gradients and objective

1: Model update: Choose model mk and error indicator ϕk

ϕk(µk) ≤ κϕ min{||∇mk(µk)|| , ∆k}

2: Step computation: Approximately solve the trust region subproblem

ˆ µk = arg min

µ∈Rnµ mk(µ)

subject to ||µ − µk|| ≤ ∆k

3: Step acceptance: Compute approximation of actual-to-predicted reduction

ρk = ψk(µk) − ψk(ˆ µk) mk(µk) − mk(ˆ µk) if ρk ≥ η1 then µk+1 = ˆ µk else µk+1 = µk end if

4: Trust region update:

if ρk ≤ η1 then ∆k+1 ∈ (0, γ ||ˆ µk − µk||)] end if if ρk ∈ (η1, η2) then ∆k+1 ∈ [γ ||ˆ µk − µk|| , ∆k] end if if ρk ≥ η2 then ∆k+1 ∈ [∆k, ∆max] end if

59 / 63

SLIDE 96

Final requirement for convergence: Adaptivity

With the approximation model, mk(µ), and gradient error indicator, ϕk(µ) mk(µ) = EIk [J (Φkur(µ, · ), µ, · )] ϕk(µ) = α1E1(µ; Ik, Φk) + α2E2(µ; Ik, Φk) + α3E4(µ; Ik, Φk) the sparse grid Ik and reduced-order basis Φk must be constructed such that the gradient condition holds ϕk(µk) ≤ κϕ min{||∇mk(µk)|| , ∆k} Define dimension-adaptive greedy method to target each source of error such that the stronger conditions hold E1(µk; I, Φ) ≤ κϕ 3α1 min{||∇mk(µk)|| , ∆k} E2(µk; I, Φ) ≤ κϕ 3α2 min{||∇mk(µk)|| , ∆k} E4(µk; I, Φ) ≤ κϕ 3α3 min{||∇mk(µk)|| , ∆k}

60 / 63

SLIDE 97

Backward facing step: minimize recirculation

0.2 1.0 0.5 0.5 Geometry and boundary conditions for backward facing step. Boundary conditions: viscous wall ( ), parametrized inflow(µ) ( ), stochastic inflow(ξ) ( ), outflow ( ). Vorticity magnitude minimized in red shaded region.

61 / 63

SLIDE 98

Mean vorticity corresponding to no inflow (left) and optimal inflow (right) along parametrized boundary.

62 / 63

SLIDE 99

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−7 10−5 10−3 10−1 101 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

63 / 63

SLIDE 100

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−7 10−5 10−3 10−1 101 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

63 / 63

SLIDE 101

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−7 10−5 10−3 10−1 101 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

63 / 63

SLIDE 102

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−7 10−5 10−3 10−1 101 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

63 / 63

SLIDE 103

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−7 10−5 10−3 10−1 101 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

63 / 63

SLIDE 104

Significant reduction in cost, if ROM only 10× faster than HDM

Cost = nHdmPrim + 0.5 × nHdmAdj + τ −1 × (nRomPrim + 0.5 × nRomAdj) 101 102 103 104 105 10−7 10−5 10−3 10−1 101 Cost |F(µ) − F(µ∗)| 5-level isotropic SG ( ), dimension-adaptive SG [Kouri et al., 2014] ( ), and proposed ROM/SG for τ = 1 ( ), τ = 10 ( ), τ = 100 ( ), τ = ∞ ( )

63 / 63