A Component-Based Reduced Basis Method for Many-Parameter Systems - - PowerPoint PPT Presentation

A Component-Based Reduced Basis Method for Many-Parameter Systems

David J. Knezevic Harvard University Institute for Applied Computational Science

Outline

• 1. Applications of the Reduced Basis Method

◮ Examples and Applications ◮ Parameter Estimation and Model Validation ◮ Toward More Parameters

• 2. Component-Based RB Method

◮ Motivation ◮ Formulation ◮ Akselos ◮ Numerical Results

• 1. Applications of the RB Method

Applications of the Reduced Basis Method

◮ Examples and Applications ◮ Parameter Estimation and Model Validation ◮ Toward More Parameters

Applications of the Reduced Basis Method

◮ Examples and Applications ◮ Parameter Estimation and Model Validation ◮ Toward More Parameters

Parametrized Partial Diff. Eqs.

Let F be a PDE operator, models some physical system (e.g. fluids, solids, acoustics, electromagnetics, ...) F depends on parameter vector µ ∈ RP: µ characterizes system (e.g. material properties, boundary conditions, geometry, ...) Given µ ∈ D ⊂ RP Find u(µ) ∈ X such that F(u(µ); µ) = 0, and output(s): ℓ(·; µ) : X → R, s(µ) = ℓ(u(µ); µ) ∈ R

Offline Implementation

RB for large-scale 3D PDEs requires high-performance implementation of “Offline stage” All results in this talk are based on:

◮ libMesh: Open source C++ parallel finite element library ◮ rbOOmit: Reduced Basis functionality within libMesh ◮ PETSc + SLEPc: Parallel linear algebra, eigensolver

Online Implementation

Once the Offline stage is complete, we can:

• 1. evaluate the RB model at any parameter value in the

pre-defined parameter range (µ ∈ D)

• 2. evaluate error bounds with respect to the underlying “truth”

finite element discretization

• 3. evaluate output quantities of interest
• 4. visualize solution field (linear combination of RB basis

functions) These Online computations are “lightweight”

◮ 1, 2, 3 are independent of FE discretization ◮ 4 depends on FE mesh, but still much cheaper than an FE

solve

Online Implementation: Smartphone App

Developed (with Phuong Huynh) open source RB smartphone app for Android to demonstrate efficiency of Online stage

Online Implementation: Smartphone App

Footbridge Example

Footbridge Example

Motivation:

◮ Develop an RB model for a Bridges to Prosperity bridge design ◮ RB model can be evaluated on a smartphone “in the field”

deck springer (I-beam) stiffener guardrails Los Montes Simple Span Bridge σd

d

(E Esp

sp,

, ρ ρsp

sp)

) wood td

d 3 3 3

E Ew

wo

• d

d =

= 11 11 (G (GP Pa) a) ρ ρw

wo

• od

d =

= 74 740 0 (k (kg/m g/m ) ) E Esteel

steel =

= 200 200 (G (GP Pa) a) ρ ρsteel

steel =

= 7 7, ,800 800 (kg/m (kg/m ) ) 2 2 ≤ ≤ t td

d ≤

≤ 20 20 (c (cm) m) 11 11 ≤ ≤ E Esp

sp ≤

≤ 200 200 (GP (GPa) a) 740 740 ≤ ≤ ρ ρsp

sp ≤

≤ 7 7,800 800 (kg/m (kg/m3

3)

) 20 20 ≤ ≤ σ σd

d ≤

≤ 1 1, ,000 000 (N/m (N/m2

2)

) ud ( ( E Esp

sp,

, ρ ρsp

sp,

, t td

d,

, σ σd

d)

) (m) (m)

Ou utpu tput t( (P Paramete arameter r) ) F Field ield ≡ ≡ Displ Displac acem emen ent( t(x x; ; P Paramete arameter r) ) P Paramete arameter r D Domain:

• main: P

P = = 4 4

“Truth” Finite element approximation

Starting point for RB: Introduce high-fidelity (“truth”) FE space X N ⊂ X, dim(X N ) = N For this footbridge mesh, N ≈ 1.5 × 106

Footbridge Example: Offline stage

Construct RB space XN ⊂ X N , N ≪ N, from solution snapshots at “greedily selected” parameters ξ1 ξ2 ξ3 . . . N = 33: Offline requires ≈ 2 hours (for 33 truth solves + extra preprocessing) on 64 processors

Footbridge Example: Online Stage

Then using RB model, footbridge problem can be solved (with error bounds) on smartphone in real-time Plot shows parameter-dependent output: Vertical displacement at bridge midpoint as function of bridge deck thickness (td)

Footbridge Example: Online Stage

We can also plot 3D solution fields on the phone1

1Uses Android’s implementation of OpenGL ES

Footbridge Example: Online Stage

Can also model dynamic response to harmonic forcing by solving parametrized Helmholtz equation, resonant frequency at ≈ 13Hz

CFD Example

CFD Example: Natural convection

We can also solve nonlinear PDEs with RB2, e.g. Boussinesq equations: Find u, p, T satisfying ∂u ∂t + 1 2 √ Gr Pr u · ∇u + √ Gr Pr∇p − ∇2u + √ Gr PrT ˆ g = 0 ∂T ∂t + 1 2 √ Gr Pr u · ∇T − 1 Pr∇2T = 0 ∇ · u = 0 Parameters:

◮ φ: direction of gravity, ˆ

g ≡ (− sin φ, 0, − cos φ)

◮ Gr: Grashof number (ratio of buoyancy to thermal diffusion) ◮ Pr = 0.71: Prandtl number (for air)

2Knezevic & Peterson, CMAME, 2011; Recent work on space-time

formulation for Boussinesq: Yano and Patera

CFD Example: Natural convection

3D domain: cross-section at y = 2.5 of the computational domain, three output regions Outputs: Average temperatures on D1, D2, D3

CFD Example: Natural convection

Finite element mesh: N = 241,520 degrees of freedom

Another example: Natural convection

FE solve requires 21.7 minutes on 128 processors, too long for many-query or real-time context x = 2.5 y = 2.5 (Gr, φ) = (6000, 0.2), t = 0.16 RB Offline stage: 42 hours on 128 processors, N = 90

CFD Example: Natural convection

RB Online on laptop:

◮ Output evaluation: 0.9 seconds ◮ Error bounds3: 18 seconds

0.05 0.1 0.15 0.2 0.25 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

(Gr, φ) = (6000, 0.2)

3Computational cost of error bounds scales as O(N4) for quadratic

nonlinearities

LabVIEW Example

Collaboration with National Instruments

Worked with engineers at National Instruments (NI) to build an RB plug-in for LabVIEW Goal: Evaluate RB models in LabVIEW on NI hardware for real-time control and data acquistion NI CompactRIO

Collaboration with National Instruments

RB model for 6-parameter thermal problem to simulate a 4-core CPU with thermal fin = ⇒ real-time PDE solves on CompactRIO

Applications of the Reduced Basis Method

◮ Offline/Online Implementations ◮ Parameter Estimation and Model Validation ◮ Toward More Parameters

Experimental Parameter Estimation

Experimental Parameter Estimation

Consider the parametrized physical system (transient heat transfer):

Macor top layer Macor bottom layer Artic Ag Thermal Grease Cu leads NiCr Wire X Y Y Z Z

Parameters: µ1: thickness of top Macor layer µ2: thermal conductivity of Arctic Ag Output: Temperature above NiCr on top surface

Experimental parameter estimation

• Prof. Ian Hunter’s group (MIT MechE) implemented this

experimental setup

Experimental parameter estimation

Goal: Determine µ1 (thickness of “top layer”) by fitting output of parametrized PDE to experimental measurements Nonlinear-least squares problem (µ → output is a nonlinear mapping), use Levenberg-Marquardt algorithm to find best fit Real-time response is desirable = ⇒ employ RB method

Experimental parameter estimation

Develop RB model for the system, transient solve requires approx. 0.01 seconds

50 100 150 200 5 10 15 20 25 30 35

RB temperature field Surface temp. output

Experimental parameter estimation

t = 0s t = 100s t = 200s

5 10 15 20 25 30 35 Run 1 Run 2 Best fit 50 100 150 200

time (s) T ฀ ฀ T T0

Best fit for µ1: 1.66mm (correct to three significant digits!), requires 1 second

In Situ Model Validation

In Situ Model Validation

Suppose we have a Physical System (PS) which provides (noisy) measurement data at times tm, 1 ≤ m ≤ M Suppose we have a proposed parametrized PDE model that represents PS Model Validation Question: For which set of parameters A ⊂ D is the proposed PDE consistent with PS? PDE (RB) Accept PS ˆ

µ ∈ D

PDE(ˆ

µ)?

sn(tm;ˆ

µ)

∆s

n(tm;ˆ µ)

ˆ

µ ∈ A

ˆ

µ ∈ A

Yes No sPS(tm) + noise

In Situ Model Validation

We employ a frequentistic approach:4 Perform an independent hypothesis test at each candidate ˆ µ to determine if ˆ µ ∈ A

◮ Hypothesis test (introduce null hypothesis, use confidence

regions, etc): Reject ˆ µ if RB(ˆ µ) vs. PS(ˆ µ) misfit ≥ ∆s

N(ˆ

µ)

◮ If hypothesis test with confidence level γ rejects ˆ

µ, then probability that PDE(ˆ µ) is consistent with PS is ≤ (1 − γ)

◮ If we reject all parameters, then we conclude that the PDE

model is invalid and should be rejected

4Huynh, Knezevic, Patera, CMAME 2012; Related prior work: Balci:1982,

McFarland:2008

In Situ Model Validation: FRP

Example application:

◮ Fiber Reinforced Polymer (FRP) is attached to concrete to

provide structural reinforcement

◮ Important to detect FRP debonding cracks, which can

compromise structural integrity

In Situ Model Validation: FRP

Debonding cracks can be detected via thermal tests

◮ Apply heat to FRP ◮ Deduce based on surface temperature measurement if there is

a crack between concrete and FRP

In Situ Model Validation: FRP

X Y Z

(κRFP

RFP)

FRP Laminate

(hc

c)

WIND

Proposed PDE model (heat equation): Find T(t; µ) ∈ X such that

• Ω

∂T ∂t v+

2

• ℓ=1
• Ωℓ

κℓ∇T·∇v+hc

• ∂ΩFRP

T v =

• ∂Ωin

qin v, ∀v ∈ X s(t; µ) = 1 |Ωout|

• Ωout

T(t; µ)

In Situ Model Validation: FRP

Develop 3D FEM/RB model, N = 98,009 vs. N = 51

In Situ Model Validation: FRP

Generate “experimental data”: add 10 different Gaussian noise realizations to FE output for µ∗ = (5, 15) Compute A on 100 × 100 grid (confidence level γ = 0.95, noise ≈ 1%)

1 2 3 4 5 6 7 8 9 10 5 10 15 20

κ κFRP

FRP

h hc

c

1 2 3 4 5 6 7 8 9 10 5 10 15 20

κFRP

FRP

h hc

c

N = 20 N = 51 N = 20: for many ˆ µ we can’t distinguish between RB error and misfit

In Situ Model Validation: FRP

Introduce a “debonding crack” between FRP and concrete, still µ∗ = (5, 15)

X Y Z

(κRFP

RFP)

FRP Laminate

(hc)

WIND DELAMINATION CRACK

1 2 3 4 5 6 7 8 9 10 5 10 15 20

κ κFRP

FRP

h hc

c

Conclusions:

◮ A is empty: Reject proposed PDE model ◮ We have detected unmodeled physics wrt FE model with 95%

confidence

Applications of the Reduced Basis Method

◮ Offline/Online Implementations ◮ Parameter Estimation and Model Validation ◮ Toward More Parameters

Toward More Parameters

Examples presented so far only had a few parameters, but Greedy algorithm can be effective for problems with more parameters For example, we consider a 3D heat transfer problem with 27 parameters (conductivities in 3 × 3 × 3 grid) and N = 241, 520

Toward More Parameters

Thermal problem with 27 parameters:5

◮ Greedy: With ntrain = 106, require N = 100 to satisfy error

tolerance of ǫ = 0.005

◮ Perform Offline computations on 512 processors (≈ 31 hours) ◮ Training set still very sparse in 26-dimensional space

(e.g. 226 ≈ 6.7 × 107) = ⇒ a posteriori error bounds crucial

5Knezevic & Peterson, CMAME 2011

hp-Greedy Algorithm

With more parameters or more complicated parameter-dependence, we need more RB basis functions to “cover” D This leads to more expensive Offline and more expensive Online To control the Online cost, natural idea is to:

◮ Adaptively subdivide the parameter domain ◮ Perform separate Greedy algorithm on each subdomain =

⇒ reduces N This is referred to as hp-Greedy Algorithm, in analogy to hp finite elements6

6See Ph.D. thesis and papers by Eftang; Stamm et al.

hp-Greedy Algorithm

For the Boussinesq example introduced earlier:7

◮ Adaptive subdivision into 45 subdomains ◮ 8× speedup of Online phase

20 40 60 80 10

−3

10

−2

10

−1

10 10

1

10

2

N ǫmax

N,M

4000 4500 5000 5500 6000 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2 Gr φ

7Eftang, Knezevic, Patera, MCMDS 2011

• 2. Component-Based RB Method

Component-Based RB Method

◮ Motivation ◮ Formulation ◮ Akselos ◮ Numerical Results

Component-Based RB Method

◮ Motivation ◮ Formulation ◮ Akselos ◮ Numerical Results

Component-Based RB Method

“Standard” RB Method has some significant limitations Many-parameter problems (e.g. 10 parameters) are difficult

◮ Very large Greedy training sets, many RB basis functions ◮ Quickly becomes prohibitively expensive in both “Offline” and

“Online” stages Limited modeling flexibility

◮ Only “parametric variations” are permitted ◮ Models must be mostly pre-determined, major

reconfigurations are not possible We can address these issues by combining RB with Domain Decomposition (DD) to obtain a component-based approach

Component-Based RB Method

Some approaches for combining model reduction with DD:

◮ Maday, Rønquist: RB Element Method, Lagrange multipliers

to glue non-conforming RB approximations

◮ Chen, Hesthaven, Maday: Seamless Reduced Basis Element

Method, DG formulation to avoid Lagrange multipliers

◮ Nguyen: Multiscale Reduced Basis method, similar to MsFEM

but uses RB for cell problems

◮ Iapichino, Quarteroni, Rozza: Reduced Basis Hybrid method,

couples components via coarse grid and Lagrange multipliers

◮ Heinkenschloss et al.: Balanced truncation model reduction

• n subdomains coupled to FE

◮ Craig-Bampton, Component Mode Synthesis: Widely used in

industry for non-parametrized eigenvalue or dynamic problems

Component-Based RB Method

We propose an alternative approach — more like “domain assembly” than “domain decomposition”8 Key features:

◮ Enables assembly of arbitrary systems from a set of

pre-computed parametrized components

◮ Efficiently handles systems with many components (and hence

many parameters)

◮ A posteriori error bounds can be computed

8Huynh, Knezevic, Patera, “A Static Condensation Reduced Basis Element

Method: Approximation and A Posteriori Error Estimation”, M2AN

Model Assembly

Acoustic muffler component library:9

µ 1

1 ≡

≡ k k

µ 1

1 ≡

≡ k k

9Huynh, Knezevic, Patera, “A Static Condensation Reduced Basis Element

Method: Complex Problems”, accepted to CMAME

Model Assembly

We can assemble many different models from these components... Question: How do we efficiently compute the PDE solution in the assembled domains?

Component-Based RB Method

◮ Motivation ◮ Formulation ◮ Akselos ◮ Numerical Results

Truth Formulation: Linear Elliptic PDE

Let X N be a “truth” FE space on global domain Ω, where dim(X N ) = N is large Consider a parametrized elliptic PDE: For µ ∈ D, find uN (µ) ∈ X N such that: a(uN (µ), v; µ) = f (v; µ), ∀v ∈ X N We then decompose this system level formulation into components using static condensation

Truth Static Condensation

Subdivide Ω into a set of components and define interface functions ψk,COM local to each component, COM:

◮ Compute a set of modes on each port10 ◮ Solve a Laplace problem for each mode to “elliptically lift”

into interior ψ1,COM ψ2,COM ψ3,COM ψ4,COM . . .

10Default: number of modes = number of FE dofs on port

Truth Static Condensation

Let X N

COM be the FE “bubble space” on COM: ◮ restriction of X N to COM ◮ vanishes on all local ports of COM

Express solution as a sum of “bubble functions” bN

COM(µ) ∈ X N COM

and “interface functions” with interface coefficients Uk(µ): uN (µ) =

• COM

bN

COM(µ) +

• “global ports”
• k

Uk(µ)ψk Apply static condensation11:

• 1. Solve for the bubble bN

COM(µ) on each component

• 2. Assemble and solve the remaining system for U(µ) ∈ RnSC

11AKA substructuring AKA block Gauss elimination AKA Schur complement

Solving for Bubble Functions

Set test functions to v ∈ X N

COM in a(uN (µ), v; µ) = f (v; µ) to get

COM-local bubble problems COM-local problem for bN

COM(µ) ∈ X N COM:

a(bN

COM(µ), v; µ) = f (v; µ) −

• “COM ports”
• k

Uk(µ) a(ψk, v; µ), for all v ∈ X N

COM

Solving for Bubble Functions

From linearity of a we have: bN

COM(µ) = bN f ,COM(µ) +

• “COM ports”
• k

Uk(µ) bN

k,COM(µ)

where bN

k,COM(µ), bN f ,COM(µ) satisfy

a(bN

k,COM(µ), v; µ)

= −a(ψk,COM, v; µ), ∀v ∈ X N

COM

a(bN

f ,COM(µ), v; µ)

= f (v; µ), ∀v ∈ X N

COM

Solving for U(µ)

Once the bubble solves are complete, we:

◮ Substitute our representation of uN (µ) into

a(uN (µ), v; µ) = f (v; µ)

◮ Test on the interface functions

This yields an nSC × nSC static condensation system for U(µ) ∈ RnSC: A(µ)U(µ) = F(µ) Typically nSC ≪ N hence this nSC × nSC solve is usually very fast, GOOD!

Truth Static Condensation

However, to assemble A(µ), F(µ) we need a component-local FE solve for each interface function on each component For large systems with many components this requires many FE solves, SLOW!

SCRBE Method

It’s natural to accelerate this procedure by replacing FE bubble solves by RB bubble solves This yields the Static Condensation Reduced Basis Element (SCRBE) Method Using RB bubble solves, we obtain an approximate static condensation system:

• A(µ)

U(µ) = F(µ) Assembly of SCRBE system only requires RB calculations, hence

• rders of magnitude faster than assembly of “truth” system

SCRBE Error Bounds

A posteriori error bounds can be developed based on matrix perturbation:

• AδU = δF − δA

U − δAδU, where δA = A(µ) − A(µ), δF = F(µ) − F(µ), δU = U(µ) − U(µ) Summing standard RB error bounds on each component gives: δF2 ≤ σ1(µ), δA2 ≤ σ2(µ) Then, for any µ we have: U(µ) − U(µ)2 ≤ σ1(µ) + σ2(µ) U(µ)2 λmin( A(µ)) − σ2(µ) ≡ ∆U(µ)

SCRBE: Port Modes

If we use all modes on a port, we typically get very high accuracy (only error is due to RB bubbles) However, using all modes can be expensive! O(103) degrees-of-freedom per port is typical for 3D problems = ⇒ significant computational cost, large storage requirements But there is a very natural solution to this issue: Use a subset of the port modes!

SCRBE: Port Reduction

Port Reduction Option 1: Truncated eigenmodes

◮ Use a subset of discrete eigenmodes of the Laplacian on the

ports Port Reduction Option 2: Empirical modes

◮ Compute a (large) set of “representative” port data by

extracting slices from FE results

◮ Subdivide into training set and test set ◮ Perform POD on the training set =

⇒ empirical modes

◮ Ensure maximum approximation error for empirical modes on

test set is below a specified tolerance

SCRBE: Port Reduction

Port reduction is essential in order to apply SCRBE to large-scale 3D problems Currently an active area of research:

◮ JL Eftang, DBP Huynh, DJ Knezevic, EM Rønquist, AT Patera. Port

reduction in static condensation, MATHMOD 2012.

◮ JL Eftang, AT Patera. Port Reduction in Component-Based Static

Condensation for Parametrized Problems, submitted to IJNME, 2012.

Key issues:

◮ A posteriori error analysis — need to account for “port

residual”

◮ Strategies for generating empirical modes

Component-Based RB Method

◮ Motivation ◮ Formulation ◮ Akselos ◮ Numerical Results

Akselos

Mission: Empower engineers (non-experts in computational methods) to perform large-scale, high-fidelity simulations Key enabler: Component-Based Reduced Basis Method

◮ User friendly: Components provide a simple “wrapper” for FE,

hides details of CAD, meshing, boundary conditions, etc

◮ Design flexibility: “Parametrized lego blocks”, enables

assembly of a huge number of parametrized models

◮ Model Re-use: Web-based library of RB components that can

be expanded over time to address a wide range of applications

Akselos

Akselos is building:

◮ Library of Components ◮ Library of Models ◮ GUI for assembling models, calling solver, visualizing results ◮ Cloud-based solver

Akselos

Akselos is building:

◮ Library of Components ◮ Library of Models ◮ GUI for assembling models, calling solver, visualizing results ◮ Cloud-based solver

Akselos

Akselos is building:

◮ Library of Components ◮ Library of Models ◮ GUI for assembling models, calling solver, visualizing results12 ◮ Cloud-based solver

12Python, Qt, OpenGL

Akselos

Akselos is building:

◮ Library of Components ◮ Library of Models ◮ GUI for assembling models, calling solver, visualizing results ◮ Cloud-based solver13

PETSc

13MPI/multithreaded

Akselos

Commercial applications:

◮ Commenced first project with a major engineering company,

in talks with several others Educational applications:

◮ Developing a simulation platform for education ◮ Integrating with edX (college courses) and ck12 (K-12 STEM)

Component-Based RB Method

◮ Motivation ◮ Formulation ◮ Akselos ◮ Numerical Results

Helmholtz Acoustics

Helmholtz Acoustics

Helmholtz equation (complex valued) on Ω:14 (1 + ikε)∆u(µ) + k2u(µ) = 0 Incoming wave on Γin: iku + ∂u ∂n = 2ik Non-reflection (radiation) on Γout:

• ik +

1 Rrad

• u + ∂u

∂n = ik Output of interest: RC(µ) =

• 1

|Γin|

• Γin u − 1
• 14ε models energy dissipation, here we use ε = 10−3

Helmholtz Acoustics: Without Port Reduction

k = 0.9 FE N ≈ 245, 000 180s per solve SCRBE nSC = 3, 879 6.8s per solve

Helmholtz Acoustics: Without Port Reduction

System a System b

0.6 0.7 0.8 0.9 0.2 0.4 0.6 0.8 1 k

• RC(µ)

Systema

5

Systemb

5

Helmholtz Acoustics: With Port Reduction

Helmholtz Acoustics: With Port Reduction

k = 1.2 k = 1.4 FE N ≈ 304, 261 219s per solve SCRBE (truncated eigen) nSC = 480 2.0s per solve SCRBE (empirical modes) nSC = 288 1.2s per solve

Helmholtz Acoustics: With Port Reduction

1.1 1.2 1.3 1.4 0.2 0.4 0.6 0.8 1 k RC(µ) scRBE−POD scRBE FE 1.35 1.36 1.37 1.38 1.39 1.4 0.7 0.75 0.8 0.85 k RC(µ) scRBE−POD scRBE FE

Linear Elasticity

Linear Elasticity

Equilibrium linear elasticity on Ω: ∂ ∂xj

• Cijkl(µ)∂uk(µ)

xl

• = fi(µ)

Displacement boundary conditions on Γd: u = gd Traction boundary conditions on Γt: σ · n = gt

Linear Elasticity

Linear Elasticity

Linear Elasticity

Currently developing a demo requested by a mining company: Structural assessment of a shiploader

Linear Elasticity

SCRBE simulation of “shuttle” part of the shiploader (the rest of the structure is coming soon...)