On Efficient Numerical Methods for Phase Field equations
Download
Report
Transcript On Efficient Numerical Methods for Phase Field equations
Efficient Numerical Methods
for Phase-Field Equations
Tao Tang
Hong Kong Baptist University
September 11-13, 2013
Russian-Chinese Workshop on Numer
Analysis and Scientific Computing
John W. Cahn (1928 -- )
John E. Hilliard (1826-1987)
The Cahn-Hilliard equation: describes the process of phase separation:
u
u f (u )
t
Microstructural evolution under the Cahn–Hilliard equation, demonstrating
distinctive coarsening and phase separation.
PHYSICS
Cahn-Hilliard equation
u
u f (u )
t
Phase separation in a binary alloy (metal, liquid,
…)
Spinodal decomposition
Mass conservation
Interface minimization
MATHEMATICS
2
E (u ) u F (u ) dx
2
d
u (t ) H 1 E (u (t ))
dt
d
d
E (u (t )) u (t )
dt
dt
2
0
H 1
High-order nonlinear diffusion equations:
How
to do time integration?
If both dynamics and steady state are
required,how to do efficient time
discretization?
Higher order methods vs. efficiency;
Adaptivity
Examples of high order nonlinear diffusion equation
Molecular Beam Epitaxy (MBE) Model [T., Xu; WB Chen et ]
2
ht 2 h (1 h )h
Inpaiting with Cahn-Hilliard Equation
[A.L. Bertozzi, etc]
ut u 1 W ' (u ) ( x )( f u )
Phase field crystal equation [Lowengrub, Wang, Wise etc]
tt t (3 2 2)
Thin Film epitaxy [J. Shen, X.M. Wang, Wise, etc]
t ( )- - 2 2
2
EXAMPLE
Cahn-Hilliard impainting
1 '
ut u F (u ) ( f u )
[Bertozzi etc. IEEE Tran. Imag. Proc. 2007, Commun. Math. Sci, 2011]
Monotonic decrease of the energy functional during the coarsening process
NUMERICAL CHALLENGES
Interior layers (i.e. thin interface)
see a figure
Time discretization:
Lower order (good for stability)?
“h”-adaptivity; “p”-adaptivity?
Thin-film epitaxy without
slope selection
Energy
Cahn-Hilliard equation
Energy
Energy
Allen-Cahn equation
Energy curves of three different models
EXAMPLE
Consider the IBV problem for the Cahn-Hilliard Eq.
t u (u u 3 u ) 0,
u ( x,0) u0 ( x),
( x, t ) R
x
Explicit Euler’s scheme (10-7)
Semi-implicit Euler’s scheme (10-4)
Implicit Euler’s scheme (10-6)!
Non-linearly stabilized scheme (implicit for the biharmonic and
non-linear terms : -- 10-3)
Linearly stabilized splitting scheme: introducing two splitting
functionals; one is contractive and the other is expansive: (10-2 ~
10-3)
A stable first-order method:
t 0, in
ut E (u )
u (t 0) u0 in
E (u) 0, u
E(u) , as u
J (E)(u)u, u , u,
J (E )(u)
dE (u )
E (u )
dt
E
2
Eyre’s method
Convexity splitting
E(u) Ec (u) Ee (u)
where Ec , Ee C 2 and are strictly convex. The
semi-implicit discretization is given by
U k 1 U k t (Ec (U k 1 ) Ee (U k ))
Various Eyre’s type or various extension:
Inpaiting problem [Schönlieb & Bertozzi]
Coarsening simulations [Vollmayr-Lee & Rutenberg]
Second-order convex splitting [J. Shen, C. Wang, Wise]
Question:
Given Eyre’s GS scheme, can we use some
iterative ideas to obtain a higher order semistable method?
Spectral deferred Correction (SDC) for y’=f(t,y)
The method is introduced by Dutt, Greengard and Rokhlin
(BIT, 2000)
Multi-implicit SDC method (Layton and Minion, 2004)
SDC with high-order RK schemes [Christliek, Qiu and Ong,
2010]
t
y(t ) y(tn ) f (s, y(s))ds, t [tn , tn1 ].
tn
Collocation Method:
Y (t ) U n K mY
tn ,i
m
m f (tn, j , j ) ln , j ( s)ds
tn
j 1
1i m
Algorithm (SDC method)
1 (Prediction).
[0]
[0]
[0] T
Use a k0-th order numerical method to compute [1 , ,m ]
2 (Correction). For j=1,…,J
Compute the residual for [ j 1]
[ j 1] un Km [ j 1] [ j 1]
Define the error function for [ j 1]
e[ j 1] (t ) ym (t ) Lm ( [ j 1] )
Form the error equation
e(t ) Km (e(t ) Lm ( [ j 1] )) Km [ j 1] Lm ( [ j 1] )
Use an k-th order method to compute i[ j 1] e(tn,i )
at the grid points tn,i on [tn,tn+1].
Define a new approximation solution [ j ] [ j 1] [ j 1]
Convergence analysis
[T., H.-H. Xie, X.-B. Yin, JSC 2012]
Theorem 1: Let [ J ] be computed in the Correction
step of the SDC Algorithm. If the step-size h is
sufficiently small, then the following error estimate
holds:
ym
[J ]
Ch
k0 J
y
k0 J
Ch
m
y
m 1
Outline
Ym [ j ] Km [ j 1] Km [ j ] KmYm Km [ j 1]
Ch Y
Ym [ j ] Ch [ j 1] [ j ] Ym [ j 1]
Ym [ j ]
[ j]
[ j 1]
Y
m
m
Ym [ j ] Ch Ym [ j 1]
Ym [ J ] Ch J Ym [0]
Note
K m1 K m 2 Ch 1 2
Algorithm (Modified SDC method) :
Same as Algorithm 1, except that after
Picard smoothing
[ ] un K m [ 1] , 1,
[ j 1] is computed we use some
, k 1,
( S1)
[0] [ j 1] .
After the above k-1 iterations, let
( S 2)
[ j 1] [k 1].
Theorem 2
Let
be computed in the Correction step of the modified SDC
Algorithm. If the step size h is sufficiently small, and if special
collocation points (e.g. Gauss) are used, then
[J ]
Ym [ J ] Ch Jk k0 y
k0
2 m 1
Ch
y
1,
m 1,
.
Proof. It follows from (S1) that
[ ] Ym KmYm Km [ 1]
[ k 1] Ym Chk 1 [0] Ym Chk 1 [ j 1] Ym
Ym [ j ] Chk [ j 1] Ym
Ym [ J ] Ch Jk Ym [0]
Efficiency comparison using Legendre-Gauss quadrature nodes
Energy evolutions with different time steps
FIRST-ORDER LINEAR SCHEME
Simple, linear discretization in time;
First-order with energy decreasing;
In space, central differencing
FFT is used
STABILTY+ACCURACY VIA P-
ADAPTIVITY
Use Energy difference at t n and tn 1 steps
Eh (u n ) Eh (u n 1 )
If the difference is small, no correction;
If the difference is large, judge how many SDC
corrections are needed.
Note most of time regimes, no corrections are
needed
coarse mesh
with
correction
u0 ( x, y) 0.05sin x sin y 0.001
Energy evolutions with different time steps and different numbers
of corrections for the Cahn-Hilliard equation
Blowing up phenomenon of semi-implicit spectral deferred correction
with uniform number of corrections
HOW MANY ITERATIONS NEEDED
N max
n
n 1
0,
if Eh (u ) Eh (u )
N max k
N max k 1
n
n 1
N p k ,
if
Eh (u ) Eh (u )
,
n
n 1
1
N max , if Eh (u ) Eh (u )
Energy curves of the thin film model without
slope selection and number of corrections
CPU time comparison
Adaptive Time Stepping:
Energy is an important physical quantity to reflect the
structure evolution.
Adaptive time step is determined by
t max(tmin ,
tmax
1 | E(t ) |2
)
∆tmin corresponds to quick evolution of the solution,
while ∆tmax to slow evolution.
[Qiao, T., Zhang, SISC, 2010]
Time adaptivity via energy variation
(Xie; T., Luo)
Numerical Scheme for C-H eqn:
U nj ,k1 U nj ,k
t
n 12
j ,k
h
U nj ,k1 U nj ,k
Eh (U )
2
2
hU
n 12
j ,k
,
U nj ,k1 U nj ,k (U nj ,k1 )2 (U nj ,k )2
2
2
2
1
U
1
4
h
h
2
h.
Stability: Eh (U n1 ) Eh (U n ).
Discrete energy identity:
Eh (U n 1 ) Eh (U n )
n 12
h
t
2
h
0
U nj ,k1 U nj ,k
2
Equi-energy:
• It follows from the numerical scheme and the energy identity that
n 1
h 2
U n 1 U n
n 1
n 12
Eh
h
2
h
• The prescribed energy decrease ( E ) equation
n 1
h 2
U n 1 U n
n 12
E
h
2
h
•Time stepping formula
t
E
h
n 12
2
h
•One step fix-point iteration to solve the prescribed energy
decrease equation
The initial condition is random in [-0.1,0.1], with periodic boundary condition and
0.001
Example [artificial dissippation]
Molecular Bean Epitaxy (MBE) Model:
Model eqn:
ht = -2h - [ (1 - |h|2)h ]
Energy identity:
d
E ( h ) ht 0
dt
where
2
1
2
E (h) h 1 h
4
2
2
ARTIFICIAL DISSIPATION:
Remedy:
h n1 h n
2 h n1 Ah n1 [(1 | h n |2 A)h n ]
t
i.e. an O(t) is added, where A > 0 is an O(1) constant.
Property: If the constant A is sufficiently large, then
E(hn+1) E(hn)?
If the numerical solution is convergent, then the
condition for A is
3
1
2
A h ,
2
2
T & C.Xu: [SINUM, 2006]
a.e. in (0, T ]
MORE ON REGULARITY:
Consider the nonlinear 2-D model for epitaxial growth
of thin films:
ht [(| h |2 1)h] 2h,
0,
Here, we prove an a-priori bound on the L-norm of h
in the 2-D case with =0.
The proof heavily relies on the maximum principle.
It is hard to see how it can be extended to the case of
(small) positive .
However, it is intuitively clear that in the case of
positive , the solution should be more regular, and
one may expect that the similar bound on h still holds.
Conclusions/Remarks
High-order time discretization is needed for highorder nonlinear diffusion equations.
The use of the SDC method seems a useful way.
Analysis of nonlinear stability and convergence
require deep understanding of the relevant PDEs
and numerical methods. [local estimates … T. & Xu
SINUM 2006, Bertozzi etc]
The analysis for adaptive schemes is highly
nontrivial. Most of the existing numerical
methods are lack of rigorous mathematical
justification.
Thanks!
http://www.math.hkbu.edu.hk/~ttang