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

COLLOQUE EDP-NORMANDIE CAEN 2013 Simulation du système cardiovasculaire et modélisation réduite 24 octobre 2013 Jean-Frédéric Gerbeau & Damiano Lombardi INRIA & UPMC Paris 6 France

COLLOQUE EDP-NORMANDIE CAEN 2013 Motivation: Blood flows • Output of interest Ex: Pressure drop • Optimization Ex: arterial by-pass optimization (Quarteroni-Rozza, Aortic coarctation Sankaran-Marsden) • Rapid prototyping Boundary conditions (Windkessel), constitutive law coefficients,... • Moderate variability Examples: pulmonary arteries Astorino, JFG, Pantz, Traoré, CMAME 2009 2

COLLOQUE EDP-NORMANDIE CAEN 2013 Motivation: electrophysiology • Optimization Ex: optimize multi-sites pacing • Inverse problems Electrocardiology Forward Inverse Potential in heart and the torso Electrocardiogram (ECG) 3

COLLOQUE EDP-NORMANDIE CAEN 2013 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 POD in a nutshell • Let ( ϕ i ) i =1 ..n be a finite element basis • Full Order Model (FOM): search for u h = P n j =1 u j ϕ j such that: d dt ( u h , ϕ i ) + a ( θ ; u h , ϕ i ) = ( f, ϕ i ) , ∀ i = 1 ..n • Compute p “snapshots” (i.e. solutions at p time instants, or parameters θ ): S 1 ( u 1 1 , . . . , u 1 n ) , . . . , S p ( u p 1 , . . . , u p n ) • Let S be the matrix whose columns are the S i , i = 1 ..p . • Singular Value Decomposition: S = Φ Σ Ψ T , with Σ = diag ( σ 1 , . . . , σ p ) • ( Φ 1 , . . . , Φ N ) : N columns of Φ corresponding to the N largest σ i , N � n • Reduced Order Model (ROM): search for U h = P N j =1 U j Φ j such that: d dt ( U h , Φ i ) + a ( U h , Φ i ) = ( f, Φ i ) , ∀ i = 1 ..N 5

COLLOQUE EDP-NORMANDIE CAEN 2013 Motivation ★ We are interested in problems with propagations ★ In the cardiovascular system: - cardiac electrophysiology - 1D model for arterial networks 6

COLLOQUE EDP-NORMANDIE CAEN 2013 ∂ t − ∂ 2 u Example 1: ∂ u ∂ x 2 = 0 0 10 1.2 t=0 exact − 2 10 POD(10) 1 − 4 10 0.8 − 6 10 eigenvalues 0.6 − 8 10 u − 10 10 0.4 − 12 10 0.2 − 14 10 0 − 16 10 − 18 − 0.2 10 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 5 10 15 20 25 30 35 40 45 50 x n POD with 10 modes Eigenvalues Simulation : D. Lombardi 7

COLLOQUE EDP-NORMANDIE CAEN 2013 Example 2: ∂ u ∂ t + c ∂ u ∂ x = 0 1 0 10 0.9 t=0 t=1 − 1 10 0.8 − 2 0.7 10 0.6 − 3 10 eigenvalues 0.5 u − 4 10 0.4 − 5 10 0.3 0.2 − 6 10 0.1 − 7 10 0 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 − 8 x 10 0 5 10 15 20 25 30 35 40 45 50 n Eigenvalues 1.2 exact POD(10) 1 exact POD(50) 1 0.8 0.8 0.6 0.6 0.4 u u 0.4 0.2 0.2 0 0 − 0.2 − 0.2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x − 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x POD with 10 modes POD with 50 modes Simulation : D. Lombardi 8

COLLOQUE EDP-NORMANDIE CAEN 2013 Basic idea • Instead of N X u ( x, t ) = β j ( t ) φ j ( x ) j =1 • use something like : N X u ( x, t ) = β j ( t ) φ j ( x, t ) j =1 • Questions : • how to compute the modes (reduced basis) ? • how to propagate them ? 9

COLLOQUE EDP-NORMANDIE CAEN 2013 Preliminary : Semi Classical Signal Analysis Laleg, Crépeau, Sorine (2007 & 2012) • Let u(x) be a non-negative signal and the “scattering” operator L χ ( u ) φ = − ∆ φ − χ u φ • Solve the eigenvalue problem L χ ( u ) φ = λφ • And approximate u(x) by the Deift-Trubowitz formula N − X u ( x ) = χ − 1 κ m φ 2 ˜ m m =1 with κ m = √− λ m where λ m are the negative eigenvalues • Choose χ > 0 to achieve a given accuracy • Remark: can be exact for “reflectionless” potentials 10

COLLOQUE EDP-NORMANDIE CAEN 2013 Reconstruction of aortic pressure From: Laleg, Médigue, Papelier, Crépeau, et al. (2010) 11

COLLOQUE EDP-NORMANDIE CAEN 2013 • In some specific cases, solitons are better the Hilbert basis, but there is no general result... 3 2.5 target target eigenfunctions eigenfunctions solitons 2.5 solitons 2 2 1.5 1.5 u u 1 1 0.5 0.5 0 0 − 0.5 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 x x 5 modes 50 modes • We will prefer the Hilbert basis ! 12

COLLOQUE EDP-NORMANDIE CAEN 2013 Reduced basis definition • Schrödinger operator associated with the solution u : L χ ( u ) φ = − ∆ φ − χ u φ ( χ > 0 is a given constant) • Solve the eigenvalue problem L χ ( u ) φ m = λ m φ m • ( φ m ) m ≥ 1 is a Hilbert basis in ( L 2 ( Ω ) , h · , · i ) • Approximate u(x) by N M X u ( x, t ) = ˜ β m ( t ) φ m ( x, t ) , with β m = h u, φ m i m =1 13

COLLOQUE EDP-NORMANDIE CAEN 2013 Goal • Consider a general nonlinear evolution PDE : ⇢ = F ( u ) ∂ t u u (0) = u 0 • For example: F ( u ) = ∆ u + ν u (1 − u ) • Compute the initial modes with L ( u 0 ) • Propagate the modes with an operator M (to be determined) 14

COLLOQUE EDP-NORMANDIE CAEN 2013 Lax Pairs in a nutshell L ( t ) φ m ( t ) = λ m ( t ) φ m ( t ) x 1 x 1 • How do the eigenmodes evolve ? x 2 x 2 x 3 x 3 • Let Q ( t ) orthogonal such that φ m ( t ) = Q ( t ) φ m (0) ∂ t φ m ( t ) = M ( t ) φ m ( t ) • Relation between L and M ? ( ∂ t L + [ L , M ]) φ m = ∂ t λ m φ m [ L , M ] = LM − ML • Remark: Lax equation for integrable systems ∂ t L + [ L , M ] = 0 15

COLLOQUE EDP-NORMANDIE CAEN 2013 Lax pair in a nutshell • Example: Korteweg de Vries equation: ∂ t u + 6 u ∂ x u + ∂ 3 x u = 0 • Lax pair: L ( u ) v = − ∂ 2 x v − uv M ( u ) v = 4 ∂ 3 x v + 6 u ∂ x v + 3 v ∂ x u u 0 ( x ) = κ 2 • Example: 1-soliton: if 2 φ 2 1 1 ( κ 1 x ) u ( x, t ) = κ 2 2 φ 2 1 1 ( κ 1 ( x − κ 2 1 t )) then: with: φ 1 ( x ) = sech(x) 16

COLLOQUE EDP-NORMANDIE CAEN 2013 Representation in the reduced basis • Approximation of the solution: N M X β m φ m u ≈ ˜ u = m =1 N M F ( u ) ≈ ˜ X F ( u ) = γ m φ m m =1 • Matrix representation of the operators: M ij = h M φ j , φ i i Λ ij = h L φ j , φ i i = diag { λ i } Θ ij = h ˜ F ( u ) φ j , φ i i 17

COLLOQUE EDP-NORMANDIE CAEN 2013 Representation of the PDE N M N M ˜ X X F ( u ) = β m φ m γ m φ m u = ˜ m =1 m =1 M φ m X ˙ X ∂ t u = F ( u ) = β m φ m + β m ∂ t φ m = γ m φ m ⇒ ˙ β + M β = γ = ⇒ 18

COLLOQUE EDP-NORMANDIE CAEN 2013 Representation of the “Lax equation” − χ∂ t u = − χ F ( u ) ( ∂ t L + LM − ML ) φ m = ∂ t λ m φ m ˙ λ m = − χ Θ mm d Λ dt + χ Θ = Λ M − M Λ χ M mp ( u ) = Θ mp λ p − λ m ( with Θ ij = h ˜ F ( u ) φ j , φ i i ) if p 6 = m and λ p 6 = λ m ( otherwise M mp = 0) 19

COLLOQUE EDP-NORMANDIE CAEN 2013 Representation of the “multiplication by F(u)” T ijk { N M Θ ij = h ˜ X F ( u ) φ j , φ i i = γ k h φ k φ j , φ i i k =1 ˙ T ijk = h ∂ t φ k φ j , φ i i + h φ k ∂ t φ j , φ i i + h φ k φ j , ∂ t φ i i N M T ijk = { M, T } (3) X ˙ Θ ij = γ k T ijk ijk k =1 { M, T } (3) X ijk = ( M li T ljk + M lj T ilk + M lk T ijl ) 20

COLLOQUE EDP-NORMANDIE CAEN 2013 Summary Full Order Space Reduced Order Space β m u F ( u ) γ m L χ ( u ) Λ = diag( λ m ) M ( u ) M mp F ( u ) · Θ mp d Λ ( ∂ t L χ + [ L , M ]) ϕ m = ∂ t λ m ϕ m dt + χ Θ = Λ M − M Λ ˙ ∂ t u = F ( u ) β + M β = γ 21

COLLOQUE EDP-NORMANDIE CAEN 2013 Summary • Generic relations between β , γ , T , M , λ , θ ˙ ( i ) β + M β = γ N M X ( ii ) = Θ ij γ k T ijk k =1 ˙ ( iii ) = λ m − χ Θ mm χ ( iv ) M mp ( u ) = Θ mp λ p − λ m { M, T } (3) ˙ ( v ) = T ijk ijk • A relation between β and γ has to be derived for the specific PDE 22

COLLOQUE EDP-NORMANDIE CAEN 2013 Example : Fisher-Kolmogorov • The Fisher-Kolmogorov equation: ∂ u ∂ t − ∆ u = ν u (1 − u ) • Using the approximation of u : N M N M N M ( ˙ X X X β j φ j + β j ∂ t φ j − β j ∆ φ j ) = ν β j φ j − ν β j β k φ i φ k j =1 j =1 j,k =1 • Hence, the relation between β and γ : N M X ( vi ) γ i = ( ν − λ i ) β i − ( ν + χ ) T ijk β j β k j,k =1 23

COLLOQUE EDP-NORMANDIE CAEN 2013 Example 1 : FKPP ∂ u  t=0 = ∆ u + ν u (1 − u ) in Ω ,   ∂ t    ∂ u = 0 on ∂ Ω , ∂ n    e − 50(( x − 0 . 5) 2 +( y − 0 . 25) 2 )  u 0 ( x, y ) =  24

COLLOQUE EDP-NORMANDIE CAEN 2013 Example 1 : FKPP ALP ( 40 modes) FEM ( 5700 DOF) 25

COLLOQUE EDP-NORMANDIE CAEN 2013 Z T k u � u ALP ) k 2 L 2 = 1 L 2 ε 2 dt k u k 2 T 0 L 2 ε T = k u ( T ) � u ALP ( T ) k L 2 k u ( T ) k L 2 26

COLLOQUE EDP-NORMANDIE CAEN 2013 Example 1 : FKPP 11350 vertices 27

COLLOQUE EDP-NORMANDIE CAEN 2013 Example 1 : FKPP FEM (11350 dofs) ALP (30 modes) 28