ToniLassila,AndreaManzoni,GianluigiRozza ReducedbasismethodforthereliablemodelreductionofNavier-Stokesequationsincardiovascularmodelling

(1)

Introduction Navier-Stokes equations Reduced basis approximation Application in cardiovascular modelling

Reduced basis method for the reliable model reduction of Navier-Stokes equations in cardiovascular modelling

Toni Lassila, Andrea Manzoni, Gianluigi Rozza CMCS - MATHICSE - ´ Ecole Polytechnique F´ ed´ erale de Lausanne

In collaboration with Alfio Quarteroni (EPFL & Politecnico di Milano)

Supported by the ERC-Mathcard Project (ERC-2008-AdG 227058), the Swiss National Science Foundation (Project 200021-122136), and the Emil Aaltonen Foundation

Model Reduction for Complex Dynamical Systems,

TU Berlin, December 2-4, 2010

(2)

1

Introduction

2

Navier-Stokes equations

3

Reduced basis approximation

4

Application in cardiovascular modelling

(3)

Challenges in modelling the human cardiovascular system

Human cardiovascular system is acomplex flow networkof different spatial and temporal scales.

When investigating fluid flow processes the flow geometries are changing over time. The geometric variation causes astrong nonlinearityin the equations.

Medical professionals are interested inaccurate simulationof spatial quantities, such as wall shear stresses at the location of a possible pathology.

Computational costscan become unacceptably high, especially if the objective is to model the entire network, and strategies to reduce numerical efforts and model order are being developed.

(4)

Parametric incompressible Navier-Stokes equations (steady case)

We consider the following model problem:

For a given parameter vectorµ∈D⊂R^P, findU(µ)∈X s.t.

a(U(µ),V;µ) =f(V;µ) ∀V∈X,∀µ∈D

whereU:= (u,p) andV:= (v,q) consist of the velocity field and the pressure, the product spaceX=V×Q⊂[H¹(Ω)]²×L²(Ω), and the problem consists of a linear parta0and a nonlinear (quadratic inU) parta1:

a(U,V;µ) :=a0(U,V;µ) +a1(U,U,V;µ) ∀U,V∈X,∀µ∈D

For example, if the parameter is simplyµ=ν(fluid viscosity), we have a0(U,V;µ) =

Z

Ω[µ∇u:∇v−pdiv(v)−qdiv(u)]dΩ a1(U,W,V) =

Z

Ωv·(u·∇)wdΩ + appropriate boundary conditions.

Typically we are interested in linear functionals of the field solutions (outputs) s(µ) :=`(U(µ)), i.e. need to find a reduced modeles(µ) that has is within certifiedtolerance of the actual outputs: |s(µ)−es(µ)|<TOL for allµ∈D.

(5)

Parametric incompressible Navier-Stokes equations (steady case)

a(U,V;µ) :=a0(U,V;µ) +a1(U,U,V;µ) ∀U,V∈X,∀µ∈D For example, if the parameter is simplyµ=ν(fluid viscosity), we have

a0(U,V;µ) = Z

Z

(6)

Parametric incompressible Navier-Stokes equations (steady case)

a(U,V;µ) :=a0(U,V;µ) +a1(U,U,V;µ) ∀U,V∈X,∀µ∈D For example, if the parameter is simplyµ=ν(fluid viscosity), we have

a0(U,V;µ) = Z

Z

(7)

Finite element approximation to the Navier-Stokes solution

Starting from initial guessU⁰, solve at each stepkof a FP iteration forU^ks.t.

a0(U^k,V;µ) +a1(U^k−1,U^k,V) =f(V) ∀V∈V×Q until convergence.

Stable discretization withP2/P1FE spaces for velocity and pressure Vh:={v∈C(Ω,R^d) :v|K∈[P2(K)]², ∀K∈Th} ⊂V Qh:={q∈C(Ω,R) :q|_K∈P1(K), ∀K∈Th} ⊂Q.

Galerkin projection in FE space: solve at each stepkforU_h^ks.t.

a0(U_h^k,Vh;µ) +a1(U_h^k−1,U_h^k,Vh) =f(Vh) ∀V_h∈Vh×Qh

until convergence.

Similar approach for the Newton’s method...

(8)

Reduced basis approximation of the finite element solution

1 Assumption: parametric manifold of FE solutionsMh⊂Xh

is 1) low dimensional and 2) depends smoothly onµ(valid for small Reynolds number)

2 Choose a representative set of parameter valuesµ¹, . . . ,µ^N

3 Snapshot solutionsuh(µ¹), . . . ,uh(µ^N) span a subspaceVh^N

for the velocity andph(µ¹), . . . ,ph(µ^N) span a subspaceQh^N

for the pressure

4 Galerkin reduced basis: givenµ∈D, findU_h^N(µ)∈X_h^N s.t.

a0(U_h^k,N,V_h^N;µ) +a1(U_h^k−1,N,U_h^k,V_h^N) =f(V_h^N) ∀V_h^N∈X_h^N

5 Adaptive sampling procedure (greedy algorithm) for the choice ofµ¹, . . . ,µ^N

M={U(µ)∈X;µ∈D} Mh={U_h(µ)∈Xh;µ∈D}

X_h^N= span{U_h(µⁱ), i= 1, . . . ,N}

Reliability / accuracy ?

1 is based on the quality of the sampling

2 relies on computable and rigorous a posteriori error estimator ∆N(µ): kU_h(µ)−U_h^N(µ)k_X≤∆_N(µ), |s(µ)−s^N(µ)| ≤∆^s_N(µ) =k`k_X0

h∆_N(µ)

(9)

Reduced basis approximation of the finite element solution

1 Assumption: parametric manifold of FE solutionsMh⊂Xh

is 1) low dimensional and 2) depends smoothly onµ(valid for small Reynolds number)

2 Choose a representative set of parameter valuesµ¹, . . . ,µ^N

3 Snapshot solutionsuh(µ¹), . . . ,uh(µ^N) span a subspaceVh^N

for the velocity andph(µ¹), . . . ,ph(µ^N) span a subspaceQh^N

for the pressure

4 Galerkin reduced basis: givenµ∈D, findU_h^N(µ)∈X_h^N s.t.

a0(U_h^k,N,V_h^N;µ) +a1(U_h^k−1,N,U_h^k,V_h^N) =f(V_h^N) ∀V_h^N∈X_h^N

5 Adaptive sampling procedure (greedy algorithm) for the choice ofµ¹, . . . ,µ^N

M={U(µ)∈X;µ∈D} Mh={U_h(µ)∈Xh;µ∈D} X_h^N= span{U_h(µⁱ), i= 1, . . . ,N}

Reliability / accuracy ?

1 is based on the quality of the sampling

2 relies on computable and rigorous a posteriori error estimator ∆N(µ):

kU_h(µ)−U_h^N(µ)k_X≤∆N(µ), |s(µ)−s^N(µ)| ≤∆^s_N(µ) =k`k_X0 h∆N(µ)

(10)

A posteriori error estimation of the reduced basis approximation

(Veroy-Patera 2005) If

τN

(µ)

<

1 and

βh

(µ)

>

0 there exists a unique solution

Uh

(µ) s.t.

||Uh

(µ)

−U_h^N

(µ)||

X≤

∆

N

(µ) =:

βh

(µ)

ρ(µ)

h

1−

p

1

−τN

(µ)

i Here:

βh(µ) is the Babuska inf-sup constant that needs to be estimated

Winf∈X_h sup

V∈X_h

da(Uh(µ);µ)(W,V)

||W||||V|| =βh(µ)>β0>0 for the Fr´echet derivative ofa(U,W,V) w.r.t first argument atUh

ρ(µ) is a Sobolev embedding constant that needs to be estimated τN(µ) :=^2ρ(µ)ε^N^(µ)

β_h(µ)² , whereεN(µ) :=||f(·;µ)−a(U_h^N,·;µ)||_X0

h is the RB residual

Note: for large viscosity we obtain the Stokes equations and the estimator simplifies to

||Uh

(µ)

−U_h^N

(µ)|| ≤ ∆

_N

(µ) =

εN

(µ)

βh

(µ)

key ingredients

(11)

A posteriori error estimation of the reduced basis approximation

(Veroy-Patera 2005) If

τN

(µ)

<

1 and

βh

(µ)

>

0 there exists a unique solution

Uh

(µ) s.t.

||Uh

(µ)

−U_h^N

(µ)||

X≤

∆

N

(µ) =:

βh

(µ)

ρ(µ)

h

1−

p

1

−τN

(µ)

i Here:

βh(µ) is the Babuska inf-sup constant that needs to be estimated

Winf∈X_h sup

V∈X_h

da(Uh(µ);µ)(W,V)

||W||||V|| =βh(µ)>β0>0 for the Fr´echet derivative ofa(U,W,V) w.r.t first argument atUh

ρ(µ) is a Sobolev embedding constant that needs to be estimated τN(µ) :=^2ρ(µ)ε^N^(µ)

β_h(µ)² , whereεN(µ) :=||f(·;µ)−a(U_h^N,·;µ)||_X0

h is the RB residual

Note: for large viscosity we obtain the Stokes equations and the estimator simplifies to

||Uh

(µ)

−U_h^N

(µ)|| ≤ ∆

_N

(µ) =

εN

(µ)

βh

(µ)

key ingredients

(12)

Generalization of reduced basis method to parametric geometries

Parametrized formulation on fixed reference domain (

Rozza 2009

)

Evaluate output of interests(µ) =`(Uh(µ))

s.t. U_h= (u_h(µ),p(µ))∈Vh(bΩ)×Qh(bΩ) solves a(Uh(µ),Vh;µ) =f(Vh;µ) ∀Vh∈X(Ω)b

a((v,p),(w,q);µ) = Z

Ωb

∂v

∂xi

νij(x,µ)∂w

∂xj

dΩ− Z

Ωbpχij(x,µ)∂wj

∂xi

dΩ− Z

Ωbqχij(x,µ)∂vj

∂xi

dΩ,

The parametrized (original) domain Ω(µ) is the image of a reference domainΩb through aparametric mappingT(·;µ) :Ωb→Ω(µ)

One possible parametrization using free-form deformations (L.-Rozza 2009) Transformation tensors (J_T=J_T(x,µ) = Jacobian ofT(x,µ))

ν(x,µ) =J⁻¹_T ν^oJ^−T_T |J_T| and χ(x,µ) =J⁻¹_T χ^o|J_T| Problem reduced to aparametric PDEssystem onΩ (reference domain)b

(13)

Reduced basis offline/online computational framework

Offline stageinvolves precomputation of structures required for the certified error estimate and choice of the reduced basis functions.

Online stagehas complexity only depending onNand allows evaluation of outputs(µ) for anyµ∈D with a certified error bounds.

(14)

Shape Optimization of Aorto-Coronaric Bypass Grafts

Shape optimization of cardiovascular geometries helps to avoid post-surgical complications

Local fluid patterns (vorticity) and wall shear stress are strictly related to the development of cardiovascular diseases

Shape optimization problem[Agoshkov-Quarteroni-Rozza 2006] minJ(Ω;v) s.t.











−ν∆v+∇p=f inΩ

∇·v= 0 inΩ

v=vg onΓD:=∂Ω\Γ_out,

−p·n+ν∂v

∂n=0 onΓout

Jo(Ω;v) = Z

Ω^df

|∇×v|²dΩ, Jo(Ω;v) =− Z

∂Ων∂v

∂n·tdΓo

Pictures taken from: Lei et al., J Vasc. Surg. 25(1997),637-646; Loth et al., Annu. Rev. Fluid Mech. 40(2008),367-393.

(15)

Shape Optimization of Aorto-Coronaric Bypass Grafts

Shape optimization of cardiovascular geometries helps to avoid post-surgical complications

Local fluid patterns (vorticity) and wall shear stress are strictly related to the development of cardiovascular diseases

Shape optimization problem[Agoshkov-Quarteroni-Rozza 2006]

minJ(Ω;v) s.t.











−ν∆v+∇p=f inΩ

∇·v= 0 inΩ

v=vg onΓD:=∂Ω\Γout,

−p·n+ν∂v

∂n=0 onΓout

Jo(Ω;v) = Z

Ω^df

|∇×v|²dΩ, Jo(Ω;v) =− Z

∂Ων∂v

∂n·tdΓo

Pictures taken from: Lei et al., J Vasc. Surg. 25(1997),637-646; Loth et al., Annu. Rev. Fluid Mech. 40(2008),367-393.

(16)

Shape optimization of aorto-coronaric bypass grafts

A possible free-form deformation approach

[Manzoni-Quarteroni-Rozza 2010]

Several analyses show a deep impact of the graft-artery diameter ratio Φ and anastomotic angleαon shear stress and vorticity distributions

Oscillatory shear stress with different graft-artery diameter ratios Φ and anastomotic anglesα.

Picture taken from F.L. Xiong, C.K. Chong, Med. Eng. & Phys. 30(2008),311-320.

In order to get a low-dimensional FFD parametrization we need to maximize the influence of the control points by placing them close to the sensitive regions

8 parameters (7 vertical•and 1 horizontal•displacements) to control the anastomotic angle, the graft-artery diameter ratio, the upper side, the lower wall

(17)

Shape optimization of aorto-coronaric bypass grafts

A possible free-form deformation approach

Several analyses show a deep impact of the graft-artery diameter ratio Φ and anastomotic angleαon shear stress and vorticity distributions

Oscillatory shear stress with different graft-artery diameter ratios Φ and anastomotic anglesα.

Picture taken from F.L. Xiong, C.K. Chong, Med. Eng. & Phys. 30(2008),311-320.

In order to get a low-dimensional FFD parametrization we need to maximize the influence of the control points by placing them close to the sensitive regions

8 parameters (7 vertical•and 1 horizontal•displacements) to control the anastomotic angle, the graft-artery diameter ratio, the upper side, the lower wall

(18)

Shape optimization of aorto-coronaric bypass grafts

RB approximation space construction

0 5 10 15 20 25

10⁻³ 10⁻² 10⁻¹ 10⁰

RB space dimension N

!N(µ)

Number of FE dofNv+Np 35997 Lattice FFD control pointsPi,j 5×6 Number of design variablesP 8 Number of RB functionsN 22 Error tolerance RB greedyε_tol^RB 5×10⁻³ Affine operator componentsQ 222

Error estimation (energy norm) for RB space construction (greedy procedure) and selected snapshots

0 5 10 15 20

0 0.2 0.4 0.6 0.8

1 µ₅

0 5 10 15 20

−0.2

−0.15

−0.1

−0.05 0 0.05 0.1 0.15 0.2

µ₁ µ₂ µ₃ µ₄ µ₆ µ₇ µ₈

Reduction in linear system dimension 500:1

Computational speedup (single flow simulation) 107 Reduction in parametric complexity w.r.t. explicit nodal deformation 102:1

(19)

Shape optimization of aorto-coronaric bypass grafts

Vorticity Minimization (downfield region)

Automatic iterative minimization procedure (sequential quadratic programming) Vorticity evaluation by using the reduced basis velocity at each step

Optimized bypass anastomosis and Stokes flow (velocity magnitude and pressure) Optimal (black) and unperturbed (grey) configurations of FFD parametrization Vorticity magnitude for the unperturbed (left) and optimal (right) configuration

(20)

Reduction of time-dependent Navier-Stokes (ongoing work)

One-dimensional prototype problem (viscous Burgers’ equation): find

u(µ)∈L²

(0,T ;V )

∩C⁰

([0,T ];L

²

(Ω)) s.t.

d

dt

(u(µ),v ) +

µ Z

Ωux

(µ)v

xdx−¹₂ Z

Ωu²

(µ)v

xdx

=

f

(v )

∀v∈V

where Ω = (0,1) and

V ⊂H¹

(Ω).

Time discretization with implicit Euler =

⇒

time-discrete equations Spatial discretization with FEM + reduced basis reduction as before Stability constant

ρN

:= inf

v∈V

da(u_h^N,k

)(v

,v

)

||v||X

not necessarily positive!

A posteriori estimator (Nguyen-Rozza-Patera 2009) for

k

= 1, . . . ,

_∆t^T

||u^k_h

(µ)−

u_h^N,k

(µ)|| ≤ 4

^k_N

(µ) :=

v u u t

∆t µ ∑^k_m=1

ε_N²

(t

^m

;µ)

∏^m−1_j=1

(1 + ∆tρ

N

(t

^j

;µ))

∏^k_m=1

(1 + ∆tρ

N

(t

^m

;

µ))

BUT the error bound grows exponentially in time

(21)

Conclusions

Reduced basis methods a reliable MOR method for parametric PDEs Parameters can also describe the (variable) flow geometry

Certified error bounds for spatial outputs of reduced field variables Extensions to noncoercive (Stokes) and nonlinear (Navier-Stokes) cases Future work

Time-dependent Navier-Stokes, improved error estimates Reduction of coupled multiphysics problems

Parameter identification and inverse problems

(22)

References

V.I. Agoshkov, A. Quarteroni and G. Rozza. Shape design in aorto-coronaric bypass using perturbation theory. SIAM J. Num. Anal. 44(1), 367–384, 2006.

D.B.P. Huynh, D.J. Knezevic, Y. Chen, J.S. Hesthaven, and A.T. Patera. A natural-norm successive constraint method for inf-sup lower bounds. Comput. Methods Appl. Mech. Engrg. 199:1963-1975, 2010.

T. L. and G. Rozza. Parametric free-form shape design with PDE models and reduced basis method.

Comput. Methods Appl. Mech. Engrg. 199(23-24):1583–1592, 2010.

A. Manzoni, A. Quarteroni, G. Rozza. Shape optimization for viscous flows by reduced basis methods and free-form deformation, J. Comp. Phys., submitted, 2010.

A. Quarteroni and G. Rozza. Numerical solution of parametrized Navier-Stokes equations by reduced basis methods. Numer. Methods Partial Differential Equations, 23(4):923–948, 2007.

N.-C. Nguyen, G. Rozza and A.T. Patera. Reduced basis approximation and a posteriori error estimation for the time-dependent viscous Burgers’ equation. Calcolo 46(3):157–185, 2009.

G. Rozza, D.B.P. Huynh, A.T. Patera. Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive PDEs. Arch. Comput. Methods Engrg.,15: 229–275, 2008.

G. Rozza and K. Veroy. On the stability of the reduced basis method for Stokes equations in parametrized domains. Comput. Methods Appl. Mech. Engrg., 196(7):1244–1260, 2007.

K. Veroy and A.T. Patera. Certified real-time solution of the parametrized steady incompressible Navier-Stokes equations; rigorous reduced-basis a posteriori error bounds. Int. J. Numer. Methods Fluids 47(89), 773-788, 2005.