Optimal transport methods for mesh generation, with applications to meteorology Chris Budd,

1. Optimal transport methods for mesh generation, with applications to meteorology Chris Budd, Emily Walsh JF Williams (SFU)

2. Have a PDE with a rapidly evolving solution u(x,t) eg. A storm How can we generate a mesh which effectively follows the solution structure on many scales? • h-refinement • p-refinement • r-refinement Also need some estimate of the solution/error structure which may be a-priori or a-posteriori Will describe the Parabolic Monge-Ampere method: an efficient n-dimensional r-refinement strategy using a-priori estimates, based on optimal transport ideas C.J. Budd and J.F. Williams, Journal of Physics A, 39, (2006), 5425--5463. C.J. Budd and J.F. Williams, SIAM J. Sci. Comp (2009) C.J. Budd, W-Z Huang and R.D. Russell, ACTA Numerica (2009)

3. r-refinement Strategy for generating a mesh by mapping a uniform mesh from a computational domain into a physicaldomain F Need a strategy for computing the mesh mapping function F which is both simple and fast

4. Equidistribution: 2D Introduce a positive unit measureM(x,y,t) in the physical domain which controls the mesh density A : set in computational domain F(A) : image set Equidistribute image with respect to the measure

5. Differentiate to give: Basic, nonlinear, equidistribution mesh equation Choose M to concentrate points where needed without depleting points elsewhere: error/physics/scaling Note: All meshes equidistribute some function M [Radon-Nicodym]

6. Choice of the monitor function M(X) • Physical reasoning • eg. Vorticity, arc-length, curvature • A-priori mathematical arguments • eg. Scaling, symmetry, simple error estimates • A-posteriori error estimates • eg. Residuals, super-convergence

7. Mesh construction Problem: in two/three -dimensions equidistribution does NOT uniquely define a mesh! All have the same area Need additional conditions to define the mesh Also want to avoid mesh tangling and long thin regions

8. Optimally transported meshes F Argue: A good mesh for solving a pde is often one which is as close as possible to a uniform mesh

9. Monge-Kantorovich optimal transport problem Minimise Subject to Also used in image registration,meteorology ….. Optimal transport helps to prevent small angles, reduce mesh skewness and prevent mesh tangling.

10. Key result which makes everything work!!!!! • Theorem: [Brenier] • There exists a unique optimally transported mesh • (b) For such a mesh the map F is the gradient of a convex function P : Scalar mesh potential Map F is a Legendre Transformation

11. Some corollaries of the Polar Factorisation Theorem Gradient map Irrotational mesh Same construction works in all dimensions

12. It follows immediately that Hence the mesh equidistribustion equation becomes (MA) Monge-Ampere equation: fully nonlinear elliptic PDE

13. Basic idea: Solve (MA) for P with appropriate (Neumann or Periodic) boundary conditions Good news: Equation has a unique solution Bad news: Equation is very hard to solve Good news: We don’t need to solve it exactly, and can instead use parabolic relaxation Q P Alternatively: Use Newton [Chacon et. al.] Use a variational approach [van Lent]

14. Relaxation uses a combination of a rescaled version of MMPDE5 and MMPDE6 in 2D (PMA) Spatial smoothing [Hou] (Invert operator using a spectral method) Ensures right-hand-side scales like Q in 2D to give global existence Averaged monitor Parabolic Monge-Ampere equation

15. Basic approach • Carefully discretise PDE & PMA in computational domain

16. Solve the coupled mesh and PDE system either (i) As one large system(stiff!) Velocity based Lagrangian approach. Works well for paraboloic blow-up type problems (JFW) or (ii) By alternating between PDE and mesh Interpolation based - rezoning - approach Need to be careful with conservation laws here!!! But this method works very well for advective meteorological problems. (EJW)

17. Because PMA is based on a geometric approach, it has a set of useful regularity properties 1. System invariant under translations, rotations, periodicity

18. 2. Convergence properties of PMA • Lemma 1: [Budd,Williams 2006] • If M(x,t) = M(x) then PMA admits the solution • (b) This solution is locally stable/convergent and the mesh evolves to an equidistributed state Proof: Follows from the convexity of P which ensures that PMA behaves locally like the heat equation

19. Lemma 2: [B,W 2006] If M(x,t) is slowly varying then the grid given by PMA is epsilon close to that given by solving the Monge Ampere equation. Lemma 3: [B,W 2006] The mapping is 1-1 and convex for all times: No mesh tangling or points crossing the boundary

20. 4. For appropriate choices of M the coupled system is scale-invariant Lemma 4: [B,W 2005] Multi-scale property If there is a natural length scale L(t) then for careful choices of M the PMA inherits this scaling and admits solutions of the form Lemma 5: [B,W 2009] Multi-scale regularity The resulting scaled meshes have the same local regularity as a uniform mesh Extremely useful properties when working with PDEs which have natural scaling laws

21. Example 1: Parabolic blow-up M is locally scale-invariant, concentrates points in the peak and keeps 50% of the points away from the peak

22. Solve in parallel with the PDE 10 10^5 Solution: y x Mesh:

23. Solution in the computational domain 10^5

24. Example 2: Tropical storm formation (Eady problem)

27. Monitor function: Maximum eigen-value of R: • Solve using rezoning and interpolation • Update solution every 10 mins • Update mesh every hour • Advection and pressure correction on adaptive mesh • Discontinuity singularity after 6.3 days

29. Conclusions • Optimal transport is a natural way to determine meshes in dimensions greater than one • It can be implemented using a relaxation process by using the PMA algorithm • Method works well for a variety of problems, and there are rigorous estimates about its behaviour • Looking good on meteorological problems • Still lots of work to be done eg. Finding efficient ways to couple PMA to the underlying PDE