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

1

Introduction Gradient flows Evolving diffeomorphisms

2

Numerical discretization by evolving diffeomorphisms Implicit-in-time discretization

3

Calculating the initial diffeomorphism The Monge Ampere equation Diffusion based algorithms for density equalizing maps

4

Numerical results

3 / 26

Nonlinear continuity equations

Let us consider a time dependent unknown probability density ρ(⋅, t) on a domain Ω ⊂ Rd, which satisfies the nonlinear continuity equation ∂tρ = −∇ ⋅ (ρv) ∶= ∇ ⋅ (ρ∇ [U ′(ρ) + V + W ∗ ρ]) . U ∶ R+ → R denotes the internal energy. V ∶ Rd → R is the confining potential. W ∶ Rd → R corresponds to an interaction potential. Nonlinear velocity is given by v = −∇ δF

δρ , where F denotes the free energy or entropy

functional F(ρ) = ∫Rd U (ρ)dx + ∫Rd V (x)ρ(x)dx + 1 2 ∫Rd ×Rd W (x − y)ρ(x)ρ(y)dxdy. Free energy is decreasing along trajectories d dt F(ρ)(t) = − ∫Rd ρ(x, t)∣v(x, t)∣2dx.

4 / 26

Nonlinear continuity equations

Included models: U (s) = s log s, V = 0, W = 0 heat equation. U (s) =

1 m−1 sm, V = 0, W = 0 porous medium (m > 1) or fast diffusion (m < 1)

equation. U (s) = s log s, V given, W = 0 Fokker Planck equations or Patlak-Keller-Segel model. U = 0, V = 0, W = log(−∣x∣) or W = 1

2∣x∣2 − 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 formalism1

Solutions ρ can be constructed by the following variational scheme: ρn+1

∆t

∈ arg inf

ρ∈K{

1 2∆t d2

W (ρn ∆t, ρ) + F(ρ)},

for a fixed ∆t > 0 and K = {ρ ∈ L1

+(Rd) ∶ ∫Rd ρ(x)dx = M , ∣x∣2ρ ∈ L1(Rd)}.

Quadratic Euclidean Wasserstein distance dW between two probability measures µ and ν, d2

W (µ, ν) ∶=

inf

T∶ν=T#µ ∫Rd ∣x − T(x)∣2dµ(x).

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.

1Jordan, 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 formalism2

Let Ω and ˜ Ω be smooth, open, bounded and connected subsets of Rd. We denote by D the set of diffeomorphisms from ¯ Ω to ¯ ˜ Ω (mapping ∂Ω onto ∂ ˜ Ω). Classic L2-gradient flow: Φn+1

∆t

∈ arg inf

Φ∈D{

1 2∆t ∥Φn

∆t − Φ∥2 L2(Ω) + I(Φ)}

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.

2Evans, 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 ρ = Φ#LN ⌞Ω 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 simulations3

Heat equation ∂ρ ∂t = ρxx ⇒ ∂Φ ∂t = − ∂ ∂x ( 1 Φx ) = Φxx (Φx )2 , Porous medium equation ∂ρ ∂t = ∂xx (ρm) ⇒ ∂Φ ∂t = − ∂ ∂x ( 1 (Φx )m ) = m Φxx (Φx )m+1 . 1D simulation of the PME for m = 2, where ∫

Φ(xi )

ρ0(y)dy = xi.

(c) Diffeomorphism Φ (d) Density ρ

3Carrillo 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 simulations3

Heat equation ∂ρ ∂t = ρxx ⇒ ∂Φ ∂t = − ∂ ∂x ( 1 Φx ) = Φxx (Φx )2 , Porous medium equation ∂ρ ∂t = ∂xx (ρm) ⇒ ∂Φ ∂t = − ∂ ∂x ( 1 (Φx )m ) = m Φ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.

3Carrillo 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 nT (cof DΦ)T ∂Φ ∂t = (cof DΦ)n ⋅ ∂Φ ∂t = 0. On the square Ω = [−1, 1]2 these boundary conditions translate to ∂Φ1 ∂t ∂Φ2 ∂x1 − ∂Φ2 ∂t ∂Φ1 ∂x1 = 0 for x1 = −1, x1 = 1 −∂Φ1 ∂t ∂Φ1 ∂x2 + ∂Φ2 ∂t ∂Φ1 ∂2 = 0 for x2 = −1, x2 = 1. Consider only diffeomorphisms which map each edge of ∂Ω to the corresponding

• ne of ∂ ˜

Ω without rotation.

11 / 26

Implicit Euler scheme

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

k

(Φ1

k

Φ2

k

) ϕk(x1, x2).

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 tn+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

• ne-dimensional Monge Kantorovic problems in x1 and x2 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 det(D2u(x)) = ρX (x) ρY (∇u(x)) Initial diffeomorphism Φ0 maps the constant density one to ρI , hence det(D2u(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: −ε∆2uε + det(D2uε) = 1 ρ0(∇uε) 0 < ε ≪ 1, Abstract discrete formulation: εAhuε

h + Bh(uε h) = Ch(uε h),

where B(u) = det(D2u) and C(u) =

1 ρi (∇u).

Corresponding Newton iteration creates a sequence {uε

k } satisfying

εAhuε

k+1 + (DBh[uε k ] − DCh[uε]k])(uε k+1 − uε k ) = −Bh(uε k ) + Ch(uε k ).

15 / 26

Vanishing viscosity approach for the MA equation

Spatial discretization: Hybrid C 0 conforming finite elements. Discrete fourth order operator: find u ∈ V = H 1 × H (div) such that ⟨Ahu, v⟩ = ∫Ω ∆u∆vdx − ∑

e∈Ei

∫e (∆u [ ∂v ∂n ] + ∆v [ ∂u ∂n ] + α he [ ∂u ∂n ] [ ∂v ∂n ]) ds, for all test functions v ∈ V , where [ ∂u

∂n ] ∶= ∂u ∂n − ̂

∂nu and α ∈ R. Nonlinear discrete operator Bh is calculated from the discrete linearization Lh ∶= DBh, which is consistent with L given by L(w) ∶= lim

t→0

B(u + tw) − B(u) t = − cof(D2u) ∶ D2(w) = −∇ ⋅ (cof(D2u)∇w). Hence we obtain ⟨Lhw, v⟩ = ∫Ω cof(D2u)∇w∇vdx + βhe ∑

e∈Ei

∫e [∂w ∂n ] [ ∂v ∂n ] ds. with β ∈ R.

16 / 26

Using the Monge Ampere equation is like ....

17 / 26

Alternative approach - density equalizing maps

Main idea: construct cartograms in which the sizes of the geographic regions appear in proportion to their population or an analogous property.

Figure : Gastner, Newman, Diffusion based maps for producing density equalizing maps, PNAS 101(20), 7499–7504, 2004

18 / 26

Density equalizing maps

Mathematical formulation: Given a domain Ω and a positive density ρ with ∫Ω ρ(x)dx = ∣Ω∣, find a map Φ ∶ Ω → Ω such that det DΦ(x) = ρ(x). for all x ∈ Ω. Consider the heat equation: ∂u ∂t (x, t) = ∆u(x, t) u(x, 0) = ρ(x), with homogeneous Neumann boundary conditions. Then u(x, t) → 1 as t → ∞. Map is constructed from the evolution: ˙ x(t) = −∇u(x, t) u(x, t) , lim

t→0 x(t) = x0

by Φ(x0) = limt→∞ x(t; x0).

19 / 26

Calculating the initial diffeomorphism

Initial diffeomorphism corresponds to the inverse map Φ−1, which maps the constant density one to ρ = ρ(x, t). Hence we determine the solution of , then solve the ODE system backward in time to obtain the inverse mapping.

(a) Density ρ (b) Initial mesh (c) Mapped mesh

20 / 26

Numerical simulations

Numerical solver:

1

Given an initial density ρ0 = ρ0(x) calculate the corresponding initial diffeomorphism by determining the solution uε = uε(x) the Monge Ampere equation for a sequence of decreasing ε.

2

Map the optimal transportation plan to the initial diffeomorphism Φ0(x) = ∇uε.

3

Solve the implicit-in-time discretization for Φ = Φ(x, t) via Newton’s method in every time step.

4

Calculate the corresponding density ρ = ρ(x, t) via ρ(Φ(x)) det(DΦ(x)) = 1. Simulation parameters: Computational domain Ω = [−1, 1]2 discretized into 5732 triangles. Monge Ampere solver: H 1 conforming finite elements of order 3. Vanishing viscosity approach: ε = 1, 0.1, 0.01, 0.001. Max error for the Newton solver: 10−6. Diffeomorphism solver: Time steps ∆t = 0.2. Max error for the Newton solver: 10−4.

21 / 26

Attraction: W (x) = 1

2∣x∣2

(d) t = 0.2 (e) t = 0.8 (f) Entropy (g) t = 0.2 (h) t = 0.8

22 / 26

Attraction-repulsion:W (x) = 1

2∣x∣2 − 1 4∣x∣4

(i) t = 0.8 (j) t = 1.4 (k) Entropy (l) t = 0.8 (m) t = 1.4

23 / 26

Attraction-repulsion: W (x) = 1

2∣x∣2 − ln(∣x∣),

V (x) = −1

4 ln(∣x∣)

(n) t = 0.4 (o) t = 1 (p) Entropy (q) t = 0.4 (r) t = 1

24 / 26

Conclusion and future work

Conclusions: Presented an numerical algorithm for solving nonlinear convection-aggregation equations, which automatically adapts the mesh in case of mass aggregation. It is based on the Lagrangian formulation and preserves energy dissipation. Can be used for general geometries due to triangular mesh. Future work: Include nonlinear diffusion. Speed up simulations (operator splitting, assembling of the convolution matrix). Higher order basis functions for the diffeomorphism solver. Boundary conditions for more general geometries; better ’treatment’ in the numerical solver.

25 / 26

Conclusion and future work

Conclusions: Presented an numerical algorithm for solving nonlinear convection-aggregation equations, which automatically adapts the mesh in case of mass aggregation. It is based on the Lagrangian formulation and preserves energy dissipation. Can be used for general geometries due to triangular mesh. Future work: Include nonlinear diffusion. Speed up simulations (operator splitting, assembling of the convolution matrix). Higher order basis functions for the diffeomorphism solver. Boundary conditions for more general geometries; better ’treatment’ in the numerical solver. Thank you very much for your attention !

25 / 26

Bibliography

• L. Evans, O. Savin, and W. Gangbo,

Diffeomorphisms and nonlinear heat flows, SIAM J. Math. Anal., 37 (2005).

• J. A. Carrillo and J. S. Moll,

Numerical simulation of diffusive and aggregation phenomena in nonlinear continuity equations by evolving diffeomorphisms, SIAM J. Sci. Comput., 31 (2009/10).

• S. Brenner, T. Gudi, M. Neilan, and L.-Y. Sung

C 0 penalty method for the fully nonlinear Monge-Ampere equation, Mathematics of Computation, 80 (2011).

• A. Avinyo, J. Sola-Morales and M. Valencia

On the diffusion algorithm for density equalizing maps with piecewise constant initial data, Mathematical Methods in the Applied Sciences, 35, (2012) M.T. Gastner and M.E.J. Newman Diffusion-based method for producing density-equalizing maps, PNAS, 101(20), 2004.

26 / 26