On the numerical simulation of nonlinear convection-aggregation - PowerPoint PPT Presentation

On the numerical simulation of nonlinear convection-aggregation equations by evolving diffeomorphisms J. Carrillo, H. Ranetbauer and M.T. Wolfram RICAM Workshop on ’Optimal Transport in the Applied Sciences’ Linz, December 2014

Special thanks to Jose Carrillo Jean David Benamou Michael Neilan 2 / 26

Outline Introduction 1 Gradient flows Evolving diffeomorphisms Numerical discretization by evolving diffeomorphisms 2 Implicit-in-time discretization Calculating the initial diffeomorphism 3 The Monge Ampere equation Diffusion based algorithms for density equalizing maps Numerical results 4 3 / 26

Nonlinear continuity equations Let us consider a time dependent unknown probability density ρ (⋅ , t ) on a domain Ω ⊂ R d , which satisfies the nonlinear continuity equation ∂ t ρ = −∇ ⋅ ( ρ v ) ∶ = ∇ ⋅ ( ρ ∇ [ U ′ ( ρ ) + V + W ∗ ρ ]) . U ∶ R + → R denotes the internal energy. V ∶ R d → R is the confining potential. W ∶ R d → R corresponds to an interaction potential. Nonlinear velocity is given by v = −∇ δ F δρ , where F denotes the free energy or entropy functional F( ρ ) = ∫ R d U ( ρ ) dx + ∫ R d V ( x ) ρ ( x ) dx + 1 2 ∫ R d × R d W ( x − y ) ρ ( x ) ρ ( y ) dxdy . Free energy is decreasing along trajectories dt F( ρ )( t ) = − ∫ R d ρ ( x , t )∣ v ( x , t )∣ 2 dx . d 4 / 26

Nonlinear continuity equations Included models: U ( s ) = s log s , V = 0 , W = 0 heat equation. U ( s ) = m − 1 s m , V = 0 , W = 0 porous medium ( m > 1 ) or fast diffusion ( m < 1 ) 1 equation. U ( s ) = s log s , V given, W = 0 Fokker Planck equations or Patlak-Keller-Segel model. 2 ∣ x ∣ 2 − 1 U = 0 , V = 0 , W = log (−∣ x ∣) or W = 1 4 ∣ x ∣ 4 correspond to attraction-(repulsion) potentials in swarming, herding and aggregation models. (a) Dictyostelium discoideum (b) Fish school 5 / 26

Gradient flow formalism 1 Solutions ρ can be constructed by the following variational scheme: ∈ arg inf ρ ∈K { W ( ρ n ∆ t , ρ ) + F( ρ )} , 1 ρ n + 1 2∆ t d 2 ∆ t for a fixed ∆ t > 0 and K = { ρ ∈ L 1 + ( R d ) ∶ ∫ R d ρ ( x ) dx = M , ∣ x ∣ 2 ρ ∈ L 1 ( R d )} . Quadratic Euclidean Wasserstein distance d W between two probability measures µ and ν , W ( µ, ν ) ∶= T ∶ ν = T # µ ∫ R d ∣ x − T ( x )∣ 2 d µ ( x ) . d 2 inf Variational scheme corresponds to the time discretization of an abstract gradient flow in the space of probability measures. Solutions can be constructed by this variational scheme; naturally preserve positivity and the free-energy decreasing property. 1 Jordan, Kinderlehrer and Otto (1999); Otto (1996, 2001); Ambrosio, Gigli and Savare (2005); Villani(2003)..... 6 / 26

Numerical schemes based on the gradient flow formalism Lagrangian schemes: discrete solution is defined by the discrete minimizer of the energy functional E : ρ n + 1 = arg inf E ( ρ, ρ n − 1 ) for n = 1 , 2 , . . . Carrillo and Moll (2009), Westdickenberg and Wilkening (2010), Matthes and Osberger (2014),.... Moving mesh methods: ’monitor’ evolution of the density by a time-dependent spatial mesh. Fixed mesh is mapped onto a moving mesh using the parabolic Monge-Ampere equation. Budd and Williams (2006, 2009); Budd, Cullen and Walsh (2012), ... Finite element and finite volume methods Burger, Carrillo and W. (2010), Bessemoulin-Chatard and Filbet (2012), Carrillo, Chertock and Huang (2014)... 7 / 26

Gradient flow formalism 2 Let Ω and ˜ Ω be smooth, open, bounded and connected subsets of R d . We denote by D the set of diffeomorphisms from ¯ Ω to ¯ Ω (mapping ∂ Ω onto ∂ ˜ ˜ Ω ). Classic L 2 -gradient flow: ∈ arg inf Φ ∈D { 2∆ t ∥ Φ n ∆ t − Φ ∥ 2 L 2 ( Ω ) + I( Φ )} 1 Φ n + 1 ∆ t converges to solutions of the PDE ∂ t ∶= u ( t ) ⋆ Φ ∂ Φ = ∇ ⋅ [ Ψ ′ ( det D Φ )( cof D Φ ) T ] − ∇ V ○ Φ − ∫ Ω ∇ W ( Φ ( x ) − Φ ( y )) dy , where I ( Φ ) = ∫ Ω Ψ ( det d Φ ) dx + ∫ Ω V ( Φ ( x )) dx + 1 2 ∫ Ω ∫ Ω W ( Φ ( x ) − Φ ( y )) dxdy . 2 Evans, Savin and Gangbo, Diffeomorphisms and nonlinear heat flows , SIAM J. Math. Anal. (2004) 8 / 26

Gradient flow formalism For a given diffeomorphism Φ ∈ D the corresponding density ρ ∈ K is given by ρ = Φ# L N ⌞ Ω which is equivalent to ρ ( Φ ( x )) det ( D Φ ) = 1 on Ω for sufficiently smooth functions. PDE for the evolving diffeomorphisms Φ = Φ ( x , t ) is the Lagrangian coordinate representation of the original Eulerian formulation for ρ = ρ ( x , t ) . Idea can be generalized for different transportation costs, e.g. relativistic cost instead of the Euclidean. 9 / 26

Examples and first simulations 3 Heat equation ∂ t = ρ xx ⇒ ∂ t = − ∂ ∂ x ( 1 ) = ∂ρ ∂ Φ Φ xx ( Φ x ) 2 , Φ x Porous medium equation ∂ t = ∂ xx ( ρ m ) ⇒ ∂ t = − ∂ ∂ x ( ( Φ x ) m ) = m ∂ρ ∂ Φ 1 Φ xx ( Φ x ) m + 1 . 1D simulation of the PME for m = 2 , where ∫ ρ 0 ( y ) dy = x i . Φ ( x i ) 0 (c) Diffeomorphism Φ (d) Density ρ 3 Carrillo and Moll, Numerical simulation of diffusive and aggregation phenomena in nonlinear continuity equations by evolving diffeomorphism , SIAM J. Sci. Comput. 31 (2009) 10 / 26

Examples and first simulations 3 Heat equation ∂ t = ρ xx ⇒ ∂ t = − ∂ ∂ x ( 1 ) = ∂ρ ∂ Φ Φ xx ( Φ x ) 2 , Φ x Porous medium equation ∂ t = ∂ xx ( ρ m ) ⇒ ∂ t = − ∂ ∂ x ( ( Φ x ) m ) = m ∂ρ ∂ Φ 1 Φ xx ( Φ x ) m + 1 . Why do we make our life so much harder ? Automatic mesh adaptation in regions of high density. Delta Dirac corresponds to a degeneration of the transport map - numerically more tractable than blow up maps. Energy dissipation. 3 Carrillo and Moll, Numerical simulation of diffusive and aggregation phenomena in nonlinear continuity equations by evolving diffeomorphism , SIAM J. Sci. Comput. 31 (2009) 10 / 26

Boundary conditions We consider with no flux boundary conditions, i.e. ∇ ρ ⋅ n = 0 on ∂ Ω , The corresponding natural boundary conditions for the equation in Lagrangian formulation are n T ( cof D Φ ) T ∂ Φ ∂ t = ( cof D Φ ) n ⋅ ∂ Φ ∂ t = 0 . On the square Ω = [ − 1 , 1 ] 2 these boundary conditions translate to ∂ Φ 1 ∂ Φ 2 − ∂ Φ 2 ∂ Φ 1 = 0 for x 1 = − 1 , x 1 = 1 ∂ t ∂ x 1 ∂ t ∂ x 1 − ∂ Φ 1 + ∂ Φ 2 = 0 for x 2 = − 1 , x 2 = 1 . ∂ Φ 1 ∂ Φ 1 ∂ t ∂ x 2 ∂ t ∂ 2 Consider only diffeomorphisms which map each edge of ∂ Ω to the corresponding one of ∂ ˜ Ω without rotation. 11 / 26

Implicit Euler scheme Implicit in time scheme: Φ n + 1 − Φ n = [∇ V ○ Φ n + 1 + ∫ ∇ W ( Φ n + 1 ( x ) − Φ n + 1 ( y )) dy ] . ∆ t Conforming finite element discretization ∆ t ∫ Ω ( Φ n + 1 − Φ n ) ϕ ( x ) dx − ∫ ∇ V ( Φ n + 1 ) ϕ ( x ) dx F ( Φ n + 1 , ϕ ) = 1 − ∫ Ω [∫ Ω ∇ W ( Φ n + 1 ( x ) − Φ n + 1 ( y )) dy ] ϕ ( x ) dx . We use lowest order H 1 conforming finite elements, i.e. Φ ( x 1 , x 2 ) = ∑ ( Φ 1 ) ϕ k ( x 1 , x 2 ) . k Φ 2 k k 12 / 26

Newton-Raphson method Solve the nonlinear operator equation F ( Φ , ϕ ) = 0 using a Newton Raphson method in every time step. Calculate the Jacobi Matrix JF and Newton update Υ k + 1 by solving JF ( Φ n + 1 , k ) Υ n + 1 , k + 1 = − F ( Φ n + 1 , k ) . Jacobi Matrix Υ n + 1 , k + 1 is a full matrix due to the convolution operator. Newton updates are calculated via Φ n + 1 , k + 1 = Φ n + 1 , k + α Υ n + 1 , k + 1 , where α is a suitable damping parameter. Newton iteration in every time step t n + 1 = ( n + 1 ) ∆ t is terminated if ∣ F ( Φ n + 1 , k + 1 , ϕ )∣ ≤ ǫ 1 or ∥ Φ n + 1 , k + 1 − Φ n + 1 , k ∥ ≤ ǫ 2 , with given error bounds ǫ 1 and ǫ 2 . No CFL type condition on the time steps ∆ t necessary. 13 / 26

Calculating the initial diffeomorphism Rectangular mesh: diffeomorphism can be constructed by subsequently solving two one-dimensional Monge Kantorovic problems in x 1 and x 2 direction. Triangular mesh: No natural ordering Monge Ampere equation gives the optimal transportation plan T = T ( x ) ∶ Ω → ˜ Ω , which maps a given probability density ρ X on Ω to another probability density ρ Y on ˜ Ω with respect to the quadratic cost ∫ Ω ∥ x − T ( x )∥ 2 ρ X ( x ) dx . Unique minimizing map is the gradient of a convex function u , i.e. T = ∇ u , which satisfies the MA equation ρ X ( x ) det ( D 2 u ( x )) = ρ Y (∇ u ( x )) Initial diffeomorphism Φ 0 maps the constant density one to ρ I , hence det ( D 2 u ( x )) = 1 ρ I (∇ u ( x )) . 14 / 26

Vanishing viscosity approach for the MA equation Vanishing viscosity approach: approximate a given fully nonlinear second order PDE by a sequence of well chosen higher order quasi-linear PDEs. Quasi-linear fourth order problem: − ε ∆ 2 u ε + det ( D 2 u ε ) = 0 < ε ≪ 1 , 1 ρ 0 (∇ u ε ) Abstract discrete formulation: h + B h ( u ε h ) = C h ( u ε h ) , ε A h u ε where B ( u ) = det ( D 2 u ) and C ( u ) = 1 ρ i (∇ u ) . Corresponding Newton iteration creates a sequence { u ε k } satisfying k + 1 + ( DB h [ u ε k ] − DC h [ u ε ] k ])( u ε k + 1 − u ε k ) = − B h ( u ε k ) + C h ( u ε ε A h u ε k ) . 15 / 26