[PPT] - Simulation du systme cardiovasculaire et modlisation rduite 24 PowerPoint Presentation

SLIDE 1

Jean-Frédéric Gerbeau & Damiano Lombardi

INRIA & UPMC Paris 6 France

Simulation du système cardiovasculaire et modélisation réduite

24 octobre 2013

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 2

Motivation: Blood flows

Output of interest

Ex: Pressure drop

2

Rapid prototyping

Boundary conditions (Windkessel), constitutive law coefficients,...

Moderate variability

Examples: pulmonary arteries

Astorino, JFG, Pantz, Traoré, CMAME 2009

Optimization

Ex: arterial by-pass optimization (Quarteroni-Rozza,

Sankaran-Marsden)

Aortic coarctation

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 3

Motivation: electrophysiology

3

Optimization

Ex: optimize multi-sites pacing Forward Inverse

Inverse problems

Electrocardiology

Potential in heart and the torso Electrocardiogram (ECG)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 4

Reduced Order Modeling

Many options:

Simplify the geometry and/or the physics

Examples:

1D model for blood flows (Formaggia, Peiró, ...)
Eikonal model in electrophysiology (Franzone, Sermesant, Frangi,...)
Keep the equations and the geometry,

but reduce the approximation space Examples:

Reduced Basis Method (Patera, Maday, Quarteroni, Rozza,...)
Proper Orthogonal Decomposition (POD), or Karhunen-Loève expansion

(Iollo, Farhat, Karniadakis, Kunisch, Gunzburger, Volkwein,...)

4

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 5

POD in a nutshell

5

Let (ϕi)i=1..n be a finite element basis
Singular Value Decomposition: S = Φ Σ ΨT , with Σ = diag(σ1, . . . , σp)
Reduced Order Model (ROM): search for Uh = PN

j=1 UjΦj such that:

d dt(Uh, Φi) + a(Uh, Φi) = (f, Φi), ∀i = 1..N

(Φ1, . . . , ΦN): N columns of Φ corresponding to the N largest σi, N n
Let S be the matrix whose columns are the Si, i = 1..p.
Full Order Model (FOM): search for uh = Pn

j=1 ujϕj such that:

d dt(uh, ϕi) + a(θ; uh, ϕi) = (f, ϕi), ∀i = 1..n

Compute p “snapshots” (i.e. solutions at p time instants, or parameters θ):

S1(u1

1, . . . , u1 n), . . . , Sp(up 1, . . . , up n)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 6

★ We are interested in problems with propagations ★ In the cardiovascular system:

cardiac electrophysiology
1D model for arterial networks

6

Motivation

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 7

7

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.2 0.2 0.4 0.6 0.8 1 1.2 x u t=0 exact POD(10)

5 10 15 20 25 30 35 40 45 50 10

−18

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10

n eigenvalues

POD with 10 modes

Simulation : D. Lombardi

Example 1: ∂u ∂t − ∂2u ∂x2 = 0 Eigenvalues

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 8

8

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1

x u

t=0 t=1

5 10 15 20 25 30 35 40 45 50 10

−8

10

−7

10

−6

10

−5

10

−4

10

−3

10

−2

10

−1

10

n eigenvalues

Simulation : D. Lombardi

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.2 0.2 0.4 0.6 0.8 1 1.2

x u

exact POD(10) 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.4 −0.2 0.2 0.4 0.6 0.8 1

x u

exact POD(50)

POD with 10 modes POD with 50 modes Example 2: ∂u ∂t + c∂u ∂x = 0 Eigenvalues

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 9

Basic idea

9

Instead of
use something like :
Questions :
how to compute the modes (reduced basis) ?
how to propagate them ?

u(x, t) =

N

X

j=1

βj(t)φj(x) u(x, t) =

N

X

j=1

βj(t)φj(x, t)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 10

Preliminary : Semi Classical Signal Analysis

Let u(x) be a non-negative signal and the “scattering” operator

10

Lχ(u)φ = −∆φ − χuφ

And approximate u(x) by the Deift-Trubowitz formula

˜ u(x) = χ−1

N−

X

m=1

κmφ2

m

Solve the eigenvalue problem

Lχ(u)φ = λφ where λm are the negative eigenvalues with κm = √−λm

Laleg, Crépeau, Sorine (2007 & 2012)

Choose χ > 0 to achieve a given accuracy
Remark: can be exact for “reflectionless” potentials

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 11

11

Reconstruction of aortic pressure

From: Laleg, Médigue, Papelier, Crépeau, et al. (2010)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 12

In some specific cases, solitons are better the Hilbert basis,

but there is no general result...

We will prefer the Hilbert basis !

12

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.5 1 1.5 2 2.5

x u

target eigenfunctions solitons

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −0.5 0.5 1 1.5 2 2.5 3

x u

target eigenfunctions solitons

5 modes 50 modes

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 13

Reduced basis definition

Schrödinger operator associated with the solution u :

13

Lχ(u)φ = −∆φ − χuφ

Approximate u(x) by
Solve the eigenvalue problem

Lχ(u)φm = λmφm ˜ u(x, t) =

NM

X

m=1

βm(t)φm(x, t), with βm = hu, φmi

(φm)m≥1 is a Hilbert basis in (L2(Ω), h·, ·i)

(χ > 0 is a given constant)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 14

Goal

Consider a general nonlinear evolution PDE :

14

⇢ ∂tu = F(u) u(0) = u0

For example: F(u) = ∆u + νu(1 − u)
Propagate the modes with an operator M (to be determined)
Compute the initial modes with L(u0)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 15

Lax Pairs in a nutshell

15

How do the eigenmodes evolve ?
Let Q(t) orthogonal such that φm(t) = Q(t)φm(0)

∂tφm(t) = M(t)φm(t)

Relation between L and M ?

L(t)φm(t) = λm(t)φm(t)

x1 x2 x1 x2 x3 x3

(∂tL + [L, M])φm = ∂tλmφm

[L, M] = LM − ML

Remark: Lax equation for integrable systems

∂tL + [L, M] = 0

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 16

Lax pair in a nutshell

Example: Korteweg de Vries equation:

16

∂tu + 6u∂xu + ∂3

xu = 0

Lax pair:

L(u)v = −∂2

xv − uv

M(u)v = 4∂3

xv + 6u∂xv + 3v∂xu

Example: 1-soliton: if

then: with: φ1(x) = sech(x) u0(x) = κ2

1

2 φ2

1(κ1x)

u(x, t) = κ2

1

2 φ2

1(κ1(x − κ2 1t))

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 17

17

Representation in the reduced basis

Matrix representation of the operators:

Mij = hMφj, φii Λij = hLφj, φii = diag{λi} u ≈ ˜ u =

NM

X

m=1

βmφm F(u) ≈ ˜ F(u) =

NM

X

m=1

γmφm

Approximation of the solution:

Θij = h ˜ F(u)φj, φii

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 18

18

∂tu = F(u) = ⇒ X ˙ βmφm + βm∂tφm = X γmφm Mφm = ⇒ ˙ β + Mβ = γ ˜ u =

NM

X

m=1

βmφm

Representation of the PDE

˜ F(u) =

NM

X

m=1

γmφm

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 19

19

Representation of the “Lax equation”

(∂tL + LM − ML)φm = ∂tλmφm Mmp(u) = χ λp − λm Θmp

if p 6= m and λp 6= λm( otherwise Mmp = 0)

˙ λm = −χΘmm dΛ dt + χΘ = ΛM − MΛ

( with Θij = h ˜ F(u)φj, φii)

−χ∂tu = −χF(u)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 20

20

{

Tijk

Representation of the “multiplication by F(u)”

Θij = h ˜ F(u)φj, φii =

NM

X

k=1

γkhφkφj, φii ˙ Tijk = h∂tφkφj, φii + hφk∂tφj, φii + hφkφj, ∂tφii

{M, T}(3)

ijk =

X (MliTljk + MljTilk + MlkTijl)

Θij =

NM

X

k=1

γkTijk ˙ Tijk = {M, T}(3)

ijk

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 21

Summary

21

Full Order Space Reduced Order Space u βm F(u) γm Lχ(u) Λ = diag(λm) M(u) Mmp F(u)· Θmp (∂tLχ + [L, M])ϕm = ∂tλmϕm

dΛ dt + χΘ = ΛM − MΛ

∂tu = F(u) ˙ β + Mβ = γ

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 22

Summary

22

A relation between β and γ has to be derived for the specific PDE

(i) ˙ β + Mβ = γ (ii) Θij =

NM

X

k=1

γkTijk (iii) ˙ λm = −χΘmm (iv) Mmp(u) = χ λp − λm Θmp (v) ˙ Tijk = {M, T}(3)

ijk

Generic relations between β, γ, T, M, λ, θ

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 23

Example : Fisher-Kolmogorov

The Fisher-Kolmogorov equation:

23

∂u ∂t − ∆u = νu(1 − u)

Using the approximation of u :

NM

X

j=1

( ˙ βjφj + βj∂tφj − βj∆φj) = ν

NM

X

j=1

βjφj − ν

NM

X

j,k=1

βjβkφiφk

Hence, the relation between β and γ:

(vi) γi = (ν − λi)βi − (ν + χ)

NM

X

j,k=1

Tijkβjβk

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 24

24

           ∂u ∂t = ∆u + νu(1 − u) in Ω, ∂u ∂n = 0 on ∂Ω, u0(x, y) = e−50((x−0.5)2+(y−0.25)2) t=0

Example 1 : FKPP

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 25

25

ALP (40 modes) FEM (5700 DOF)

Example 1 : FKPP

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 26

26

ε2

L2 = 1

T Z T ku uALP )k2

L2

kuk2

L2

dt εT = ku(T) uALP (T)kL2 ku(T)kL2

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 27

Example 1 : FKPP

27

11350 vertices

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 28

28

Example 1 : FKPP

FEM (11350 dofs) ALP (30 modes)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 29

29

Example 2: Electrophysiology

8 > < > : Am ✓ Cm ∂u ∂t + Iion(u, w) ◆ div(σIru) = AmIapp ∂w ∂t g(vm, w) =

“Monodomain equation”:
FitzHugh Nagumo ionic current:

⇢ Iion(u, w) = su(u − a)(u − 1) + w g(u, w) = ✏(u − w)

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 30

30

Example 2: Electrophysiology

FEM ALP

25 modes 6000 dof

t=0 ms t=125 ms t=250 ms t=375 ms

Monodomain equation with Fitzhugh-Nagumo model

“Infarcted” area

Simulation : E. Schenone

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 31

Example 3 : arterial network

55 arteries
FEM, Taylor-Galerkin
6000 DOF
1 cycle: about 1 minute

31

∂tA + ∂x(Au) = 0

∂tu + ∂x ✓u2 2 + p ◆ = 0

p = pe + β ρ ⇣ A1/2 − A1/2 ⌘

1 2 3

Inlet: Q
Outlet: RCR
Bifurcations :

(Au)1 = (Au)2 + (Au)3 1 2ρu2

1 + p1 = 1

2ρu2

2,3 + p2,3

Simulation : D. Lombardi

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 32

32

a) b) c) d)

3 or 4 points per branch
10 modes
1 cycle : about 1 second
ALP: Temporal signal with Schrödinger eigenfunctions

∂τq + ∂ξA = 0 ∂τπ + ∂ξu = 0

τ = x, ξ = t

∂tA + ∂x(Au) = 0 ∂tu + ∂x ✓u2 2 + p ◆ = 0

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 33

33

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −0.01 0.01 0.02 0.03 0.04 0.05 0.06

t u

FEM ROM 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 −0.01 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08

t u

reference ROM 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16

t u

reference ROM 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.02 0.04 0.06 0.08 0.1 0.12

t u

FEM ROM

Coronary Femoral artery Aorta inlet Aorta outlet

COLLOQUE EDP-NORMANDIE CAEN 2013

SLIDE 34

Perspectives

Still many things to test and to do

★ Try other operators than Schrödinger ★ Adaptation of the number of modes ★ Error estimator ★ Stability analysis ★ Non-polynomial nonlinearities

In progress:

★ Inverse problems for arterial networks ★ Cardiac electrophysiology (with Elisa Schenone)

34