Transcript Optimal transport methods for mesh generation Chris Budd (Bath), JF Williams (SFU)

Optimal transport methods
for mesh generation
Chris Budd (Bath), JF Williams (SFU)
Have a PDE with a rapidly evolving solution u(x,t)
How can we generate a mesh which effectively follows the
solution structure?
• 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 an efficient n-dimensional r-refinement
strategy using a-priori estimates
r-refinement
Strategy for generating a mesh by mapping a uniform mesh
from a computational domain  C into a physical domain  P
F
C ( , )
P ( X , Y )
Need a strategy for computing the mesh mapping
function F which is both simple and fast
Equidistribution: 2D
Introduce a positive unit measure M(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
 dd  
A
F ( A)
M ( X , Y , t )dXdY
JF: I wrote this talk using the assumption of
unit measure for M as it simplifies the
presentation. Of course we don’t need this in
general
Differentiate:
( X , Y )
M ( X ,Y , t)
1
 ( , )
Basic, nonlinear, equidistribution
mesh equation
Equidistribution in one-dimension
This is a very well defined process
Let computational and physical domains both be the unit
interval [0,1]. Mesh function X(xi)
dX
M ( X , t)
1
d
X(0) = 0, X(1) = 1
Basic mesh equation: unique solution if M>0
Hard, and unecessary, to solve exactly!
Various approaches to the solution …
1. Geometric conservation
d
dX 
 M ( X , t )
  0
dt 
d 
M t  (MV ) X  0
Solve for V
X
V
t
Some CJB notes on possible issues with GCL
• Have to start with an equidistributed mesh
• Have to know M_t
• Have to calculate V then calculate X
• Potential problems with mesh crossing
• Generalises to higher dimensions
M t   X .(MV )  0,  X V  0
2. Relaxation
Evolve towards an equidistributed mesh
or

X
X t   M ( X , t )


 X t

X
  M ( X , t )








(MMPDE5)
Mesh PDE
[Russell]
(MMPDE6)
Very effective provided that the time-scales for the
mesh evolution are smaller than those for the evolution
of the underlying PDE: MOVCOL Code [Huang, Russell]
Implementation:
Underlying PDE solution:
u ( x, t )
Moving mesh:
X i (t )
Approximate solution:
U i (t )  u( X i (t ), t )
Discretise Underlying PDE (in Lagrangian form) and
Mesh PDE in the computational variable
Solve the resulting ODEs
 [0,1]
Back to two-dimensions
Problem: in two-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
Optimally transported meshes
F
C ( , )
P ( X , Y )
Argue: A good mesh is one which is as close as
possible to a uniform mesh
Monge-Kantorovich optimal transport problem
Minimise
I ( X ,Y ) 

( X , Y )  ( , ) dd
C
Subject to
( X , Y )
M ( X ,Y , t)
1
 ( , )
Also used in image registration,meteorology
2
Intuitively a good approach
A=2
Small I
A=1
A=2
Larger I
Optimal transport helps to prevent small angles
Brenier’s polar factorisation theorem
Key result which makes everything work!!!!!
Theorem: [Brenier]
(a)There exists a unique optimally transported mesh
(b) For such a mesh the map F is the gradient of a
convex function P( , )
P : Scalar mesh potential
Map F is a Legendre Transformation
Some corollaries of the Polar Factorisation Theorem
( X , Y )   P  (P , P )
X  Y
Gradient map
Irrotational mesh
Convexity of P prevents mesh
tangling
Some simple examples
1
P
M
P
 2 2 
  
2
 2
 A
T
2
b 
T
P  A( )  B( )
Uniform enlargement scale
factor 1/M
Linear map. A is symmetric
positive definite
Tensor product mesh
It follows immediately that
 P
( X , Y )
 H ( P)  det
( , )
 P
P 
  P P  P 2
P 
Hence the mesh equidistribustion equation becomes
M (P, t ) H ( P)  1
(MA)
Monge-Ampere equation: fully nonlinear elliptic PDE
Basic idea: Solve (MA) for P with appropriate
(Neumann) 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!
Use relaxation as in the MMPDE equations
Relaxation uses a combination of a rescaled
version of MMPDE5 and MMPDE6 in 2D
 I   Qt  M (Q) H (Q) 
1/ 2
Spatial smoothing
(Invert operator
using a spectral
method)
Averaged
monitor
(PMA)
Ensures right-handside scales like Q in
2D to give global
existence
Parabolic Monge-Ampere equation
Discretise and solve PMA in parallel with the Lagrangian
form of the PDE (possibly using a temporal rescaling)
Useful properties of PMA
Lemma 1: [Budd,Williams 2006]
(a) If M(X,t) = M(X) then PMA admits the solution
Q( , t )  P( )  t
X ( )   Q  P
(b) This solution is locally stable.
Proof: Follows from the convexity of P which ensures
that PMA behaves locally like the heat equation
Lemma 2: [Budd,Williams 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: [Budd,Williams 2006]
The mapping is 1-1 and convex for all times
Additional issues for CJB and JFW
• If M is slowly varying then Q stays epsilon close
to the solution of the Monge-Ampere equation
• Need rigorous proof that Q remains convex and
of global convergence and a maximum principle
• When we solve in a rectangular domain then
there is a mild loss of regularity in the corners.
• There is also an boundary orthogonality of the
grid in the physical domain. This is both good and
bad
Lemma 4: [Budd,Williams 2005]
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
Q( , t )  L(t ) P( )
X  L(t )Y ( )
Extremely useful property when working with
PDEs which have natural scaling laws eg.
ut   X u  u
3
L(t )  (T  t )
1
u4
1
M ( X , t) 

4
2  u dXdY 2
1/ 2
log(T  t )
1/ 2
Examples of applications
1. Prescribe M(X,t) and solve PMA
3
u


u

u
2. Solve in parallel with the PDE
t
X
10
10^5
Solution:
Y
Mesh:
X
Solution in the computational domain
10^5


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
• Still lots of work to be done eg. Finding efficient
ways to couple PMA to the underlying PDE