LINEBURG


<< . .

 4
( 23)



. . >>

Similar to the second-order accurate correction, the third-order ac-
φ(x) = Q0 (x) + Q1 (x) + Q2 (x) + Q3 (x) (3.17) 3 3
curate correction is obtained by comparing |Dk +1/2 φ| and |Dk +3/2 φ|.
that can be di¬erentiated and evaluated at xi to ¬nd (φ’ )i and (φ+ )i . That 3 3 3
If |Dk +1/2 φ| ¤ |Dk +3/2 φ|, we set c = Dk +1/2 φ; otherwise, we set
x x
is, we use 3
c = Dk +3/2 φ. Then we de¬ne
φx (xi ) = Q1 (xi ) + Q2 (xi ) + Q3 (xi ) (3.18) Q3 (x) = c (x ’ xk )(x ’ xk ’ xk
+1 )(x +2 ), (3.23)
(φ’ )i (φ+ )i .
to de¬ne and Note that the constant Q0 (x) term vanishes so that
x x
upon di¬erentiation.
3(i ’ k )2 ’ 6(i ’ k ) + 2 ( x)2
Q3 (xi ) = c (3.24)
To ¬nd φ’ we start with k = i ’ 1, and to ¬nd φ+ we start with k = i.
x x
Then we de¬ne is the third-order accurate correction to the approximation of φx in
equation (3.18).
1
Q1 (x) = (Dk+1/2 φ)(x ’ xi ), (3.19)
so that
3.4 Hamilton-Jacobi WENO
1
Q1 (xi ) = Dk+ 1 φ, (3.20)
2

When calculating (φ’ )i , the third-order accurate HJ ENO scheme uses a
implying that the contribution from Q1 (xi ) in equation (3.18) is the back- x
ward di¬erence in the case of φ’ and the forward di¬erence in the case subset of {φi’3 , φi’2 , φi’1 , φi , φi+1 , φi+2 } that depends on how the stencil
x
34 3. Motion in an Externally Generated Velocity Field 3.4. Hamilton-Jacobi WENO 35

Reference [89] pointed out that setting ω1 = 0.1 + O(( x)2 ), ω2 = 0.6 +
is chosen. In fact, there are exactly three possible HJ ENO approximations
to (φ’ )i . De¬ning v1 = D’ φi’2 , v2 = D’ φi’1 , v3 = D’ φi , v4 = D’ φi+1 , O(( x)2 ), and ω3 = 0.3 + O(( x)2 ) still gives the optimal ¬fth-order
x
and v5 = D’ φi+2 allows us to write accuracy in smooth regions. In order to see this, we rewrite these as ω1 =
0.1 + C1 ( x)2 , ω2 = 0.6 + C2 ( x)2 and ω3 = 0.3 + C3 ( x)2 and plug them
v1 7v2 11v3
φ1 = into equation (3.28) to obtain
’ + , (3.25)
x
3 6 6
0.1φ1 + 0.6φ2 + 0.3φ3
v2 5v3 v4 (3.29)
φ2 = ’ + x x x
+, (3.26)
x
6 6 3
and
and
C1 ( x)2 φ1 + C2 ( x)2 φ2 + C3 ( x)2 φ3 (3.30)
x x x
v3 5v4 v5
φ3 = ’
+ (3.27) as the two terms that are added to give the HJ WENO approximation
x
3 6 6
to φx . The term given by equation (3.29) is the optimal approximation
as the three potential HJ ENO approximations to φ’ . The goal of HJ ENO that gives the exact value of φx plus an O(( x)5 ) error term. Thus, if
x
the term given by equation (3.30) is O(( x)5 ), then the entire HJ WENO
is to choose the single approximation with the least error by choosing the
approximation is O(( x)5 ) in smooth regions. To see that this is the case,
smoothest possible polynomial interpolation of φ.
¬rst note that each of the HJ ENO φk approximations gives the exact value
In [107], Liu et al. pointed out that the ENO philosophy of picking exactly x
of φx , denoted by φx , plus an O(( x)3 ) error term (in smooth regions).
E
one of three candidate stencils is overkill in smooth regions where the data
are well behaved. They proposed a weighted ENO (WENO) method that Thus, the term in equation (3.30) is
takes a convex combination of the three ENO approximations. Of course,
C1 ( x)2 φE + C2 ( x)2 φE + C3 ( x)2 φE (3.31)
if any of the three approximations interpolates across a discontinuity, it is x x x
given minimal weight in the convex combination in order to minimize its plus an O(( x)2 )O(( x)3 ) = O(( x)5 ) term. Since, each of the Ck is O(1),
contribution and the resulting errors. Otherwise, in smooth regions of the as is φE , this appears to be an O(( x)2 ) term at ¬rst glance. However,
x
¬‚ow, all three approximations are allowed to make a signi¬cant contribu- since ω1 + ω2 + ω3 = 1, we have C1 + C2 + C3 = 0, implying that the term
tion in a way that improves the local accuracy from third order to fourth in equation (3.31) is identically zero. Thus, the HJ WENO approximation
order. Later, Jiang and Shu [89] improved the WENO method by choosing is O(( x)5 ) in smooth regions. Note that [107] obtained only fourth-order
the convex combination weights in order to obtain the optimal ¬fth-order accuracy, since they chose ω1 = 0.1 + O( x), ω2 = 0.6 + O( x), and
accuracy in smooth regions of the ¬‚ow. In [88], following the work on HJ ω3 = 0.3 + O( x).
ENO in [127], Jiang and Peng extended WENO to the Hamilton-Jacobi In order to de¬ne the weights, ωk , we follow [88] and estimate the
framework. This Hamilton-Jacobi WENO, or HJ WENO, scheme turns smoothness of the stencils in equations (3.25), (3.26), and (3.27) as
out to be very useful for solving equation (3.2), since it reduces the errors
13 1
by more than an order of magnitude over the third-order accurate HJ ENO (v1 ’ 2v2 + v3 )2 + (v1 ’ 4v2 + 3v3 )2 ,
S1 = (3.32)
12 4
scheme for typical applications.
13 1
The HJ WENO approximation of (φ’ )i is a convex combination of the (v2 ’ 2v3 + v4 )2 + (v2 ’ v4 )2 ,
S2 = (3.33)
x
12 4
approximations in equations (3.25), (3.26), and (3.27) given by
and
φ x = ω1 φ 1 + ω2 φ 2 + ω3 φ 3 , (3.28)
x x x
13 1
(v3 ’ 2v4 + v5 )2 + (3v3 ’ 4v4 + v5 )2 ,
where the 0 ¤ ωk ¤ 1 are the weights with ω1 + ω2 + ω3 = 1. The key S3 = (3.34)
12 4
observation for obtaining high-order accuracy in smooth regions is that
weights of ω1 = 0.1, ω2 = 0.6 and ω3 = 0.3 give the optimal ¬fth-order respectively. Using these smoothness estimates, we de¬ne
accurate approximation to φx . While this is the optimal approximation, it is
0.1
valid only in smooth regions. In nonsmooth regions this optimal weighting ±1 = , (3.35)
(S1 + )2
can be very inaccurate, and we are better o¬ with digital (ωk = 0 or
0.6
ωk = 1) weights that choose a single approximation to φx , i.e., the HJ ±2 = , (3.36)
(S2 + )2
ENO approximation.
36 3. Motion in an Externally Generated Velocity Field 3.5. TVD Runge-Kutta 37

The function (φ+ )i is constructed with a subset of {φi’2 , φi’1 , φi ,
and x
φi+1 , φi+2 , φi+3 }. De¬ning v1 = D+ φi+2 , v2 = D+ φi+1 , v3 = D+ φi ,
.3
±3 = (3.37) v4 = D+ φi’1 , and v5 = D+ φi’2 allows us to use equations (3.25), (3.26),
(S3 + )2
and (3.27) as the three HJ ENO approximations to (φ+ )i . Then the HJ
x
with WENO convex combination is given by equation (3.28) with the weights
given by equations (3.39), (3.40), and (3.41).
= 10’6 max v1 , v2 , v3 , v4 , v5 + 10’99 ,
22222
(3.38)
where the 10’99 term is set to avoid division by zero in the de¬nition of
the ±k . This value for epsilon was ¬rst proposed by Fedkiw et al. [69],
3.5 TVD Runge-Kutta
where the ¬rst term is a scaling term that aids in the balance between
the optimal ¬fth-order accurate stencil and the digital HJ ENO weights.
HJ ENO and HJ WENO allow us to discretize the spatial terms in
In the case that φ is an approximate signed distance function, the vk that
equation (3.2) to ¬fth-order accuracy, while the forward Euler time dis-
approximate φx are approximately equal to one, so that the ¬rst term in
equation (3.38) can be set to 10’6 . This ¬rst term can then absorb the cretization in equation (3.4) is only ¬rst-order accurate in time. Practical
second term, yielding = 10’6 in place of equation (3.38). Since the ¬rst experience suggests that level set methods are sensitive to spatial accu-
racy, implying that the ¬fth-order accurate HJ WENO method is desirable.
term in equation (3.38) is only a scaling term, it is valid to make this vk ≈ 1
On the other hand, temporal truncation errors seem to produce signi¬-
estimate in higher dimensions as well.
cantly less deterioration of the numerical solution, so one can often use the
A smooth solution has small variation leading to small Sk . If the Sk
low-order accurate forward Euler method for discretization in time.
are small enough compared to , then equations (3.35), (3.36), and (3.37)
become ±1 ≈ 0.1 ’2 , ±2 ≈ 0.6 ’2 , and ±3 ≈ 0.3 ’2 , exhibiting the proper There are times when a higher-order temporal discretization is necessary
in order to obtain accurate numerical solutions. In [150], Shu and Osher
ratios for the optimal ¬fth-order accuracy. That is, normalizing the ±k to
proposed total variation diminishing (TVD) Runge-Kutta (RK) methods to
obtain the weights
increase the accuracy for a method of lines approach to temporal discretiza-
±1
ω1 = , (3.39) tion. The method of lines approach assumes that the spatial discretization
±1 + ±2 + ±3
can be separated from the temporal discretization in a semidiscrete manner
±2
ω2 = , (3.40) that allows the temporal discretization of the PDE to be treated indepen-
±1 + ±2 + ±3
dently as an ODE. While there are numerous RK schemes, these TVD RK
schemes guarantee that no spurious oscillations are produced as a conse-
and
quence of the higher-order accurate temporal discretization as long as no
±3
ω3 = (3.41) spurious oscillations are produced with the forward Euler building block.
±1 + ±2 + ±3
The basic ¬rst-order accurate TVD RK scheme is just the forward Euler
gives (approximately) the optimal weights of ω1 = 0.1, ω2 = 0.6 and ω3 = method. As mentioned above, we assume that the forward Euler method
0.3 when the Sk are small enough to be dominated by . Nearly optimal is TVD in conjunction with the spatial discretization of the PDE. Then
weights are also obtained when the Sk are larger than , as long as all the Sk higher-order accurate methods are obtained by sequentially taking Euler
are approximately the same size, as is the case for su¬ciently smooth data. steps and combining the results with the initial data using a convex com-
On the other hand, if the data are not smooth as indicated by large Sk , bination. Since the Euler steps are TVD (by assumption) and the convex
then the corresponding ±k will be small compared to the other ±k ™s, giving combination operation is TVD as long as the coe¬cients are positive, the
that particular stencil limited in¬‚uence. If two of the Sk are relatively large, resulting higher-order accurate TVD RK method is TVD. Unfortunately, in
then their corresponding ±k ™s will both be small, and the scheme will rely our speci¬c case, the HJ ENO and HJ WENO schemes are not TVD when
most heavily on a single stencil similar to the digital behavior of HJ ENO. used in conjunction with upwinding to approximate equation (3.4). How-
In the unfortunate instance that all three of the Sk are large, the data ever, practical numerical experience has shown that the HJ ENO and HJ
are poorly conditioned, and none of the stencils are particularly useful. WENO schemes are most likely total variation bounded (TVB), implying
This case is problematic for the HJ ENO method as well, but fortunately that the overall method is also TVB using the TVD RK schemes.
it usually occurs only locally in space and time, allowing the methods to The second-order accurate TVD RK scheme is identical to the standard
repair themselves after the situation subsides. second-order accurate RK scheme. It is also known as the midpoint rule,
38 3. Motion in an Externally Generated Velocity Field 3.5. TVD Runge-Kutta 39

as the modi¬ed Euler method, and as Heun™s predictor-corrector method. in practical calculations, especially since the HJ WENO scheme usually
First, an Euler step is taken to advance the solution to time tn + t, loses accuracy and looks a lot like the third-order accurate HJ ENO scheme
in many interesting areas of the ¬‚ow. Also, the fourth-order accurate (and
φn+1 ’ φn
+Vn· φn = 0, higher) TVD RK methods require both upwind and downwind di¬erencing
(3.42)
t approximations, doubling the computational cost of evaluating the spatial
followed by a second Euler step to advance the solution to time tn + 2 t, operators. See [150] for fourth- and ¬fth-order accurate TVD RK schemes.
Finally, we note that a rather interesting approach to TVD RK schemes
φn+2 ’ φn+1
+ V n+1 · φn+1 = 0, (3.43) has recently been carried out by Spiteri and Ruuth [154], who proposed
t
increasing the number of internal stages so that this number exceeds the
followed by an averaging step order of the method.
1 n 1 n+2
φn+1 = φ+φ (3.44)
2 2
that takes a convex combination of the initial data and the result of two
Euler steps. The ¬nal averaging step produces the second-order accurate
TVD (or TVB for HJ ENO and HJ WENO) approximation to φ at time
tn + t.
The third-order accurate TVD RK scheme proposed in [150] is as follows.
First, an Euler step is taken to advance the solution to time tn + t,
φn+1 ’ φn
+Vn· φn = 0, (3.45)
t
followed by a second Euler step to advance the solution to time tn + 2 t,
φn+2 ’ φn+1
+ V n+1 · φn+1 = 0, (3.46)
t
followed by an averaging step
3 n 1 n+2
1
φn+ 2 = φ+φ (3.47)
4 4
that produces an approximation to φ at time tn + 1 t. Then another Euler
2
step is taken to advance the solution to time tn + 3 t,
2
3 1
φn+ 2 ’ φn+ 2 1 1
+ V n+ 2 · φn+ 2 = 0, (3.48)
t
followed by a second averaging step
1 n 2 n+ 3
φn+1 = φ+φ 2 (3.49)
3 3
that produces a third-order accurate approximation to φ at time tn + t.
This third-order accurate TVD RK method has a stability region that
includes part of the imaginary axis. Thus, a stable (although ill-advised)
numerical method results from combining third-order accurate TVD RK
with central di¬erencing for the spatial discretization.
While fourth-order accurate (and higher) TVD RK schemes exist, this
improved temporal accuracy does not seem to make a signi¬cant di¬erence
42 4. Motion Involving Mean Curvature

1 1 1

0.5
0.5 0.5

0
0 0

4 ’0.5 ’0.5 ’0.5

’1 ’1 ’1
1
0
’1 ’1 1
0 ’1 0 1
Motion Involving Mean Curvature 1 1 1

0.5 0.5 0.5

0 0 0

’0.5 ’0.5 ’0.5

’1 ’1 ’1
’1 1
0 ’1 0 1 0
’1 1

1 1 1

0.5 0.5 0.5

0 0 0

’0.5 ’0.5 ’0.5

’1 ’1 ’1
0 1
’1 1
’1 0 0
’1 1


4.1 Equation of Motion Figure 4.1. Evolution of a wound spiral in a curvature-driven ¬‚ow. The
high-curvature ends of the spiral move signi¬cantly faster than the elongated
body section.
In the last chapter we discussed the motion of an interface in an externally
generated velocity ¬eld V (x, t). In this chapter we discuss interface motion
for a self-generated velocity ¬eld V that depends directly on the level set T · φ = 0 for any tangent vector T , implying that the tangential velocity
function φ. As an example, we consider motion by mean curvature where components vanish when plugged into the level set equation. For example,
the interface moves in the normal direction with a velocity proportional in two spatial dimensions with V = Vn N + Vt T , the level set equation
to its curvature; i.e., V = ’bκN , where b > 0 is a constant and κ is the
curvature. When b > 0, the interface moves in the direction of concavity, φ t + Vn N + Vt T · φ=0 (4.1)
so that circles (in two dimensions) shrink to a single point and disappear.
is equivalent to
When b < 0, the interface moves in the direction of convexity, so that
circles grow instead of shrink. This growing-circle e¬ect leads to the growth
φ t + Vn N · φ = 0, (4.2)
of small perturbations in the front including those due to round-o¬ errors.
Because b < 0 allows small erroneous perturbations to incorrectly grow since T · φ = 0. Furthermore, since
into O(1) features, the b < 0 case is ill-posed, and we do not consider it
| φ|2
φ
here. Figure 4.1 shows the motion of a wound spiral in a curvature-driven = | φ|,
N· · (4.3)
φ= φ=
| φ|
| φ|
¬‚ow. The high-curvature ends of the spiral move signi¬cantly faster than
the relatively low curvature elongated body section. Figure 4.2 shows the
we can rewrite equation (4.2) as
evolution of a star-shaped interface in a curvature-driven ¬‚ow. The tips of
the star move inward, while the gaps in between the tips move outward. φt + Vn | φ| = 0 (4.4)
The velocity ¬eld for motion by mean curvature contains a component
where Vn is the component of velocity in the normal direction, other-
in the normal direction only, i.e., the tangential component is identically
wise known as the normal velocity. Thus, motion by mean curvature is
zero. In general, one does not need to specify tangential components when
characterized by Vn = ’bκ.
devising a velocity ¬eld. Since N and φ point in the same direction,
4.1. Equation of Motion 43 44 4. Motion Involving Mean Curvature

step (or a forward Euler substep in the case of RK), the new value of φ
1 1 1
is not a signed distance function, and equations (4.5) and (4.6) can no
0.5
0.5 0.5
longer be interchanged. If this new value of φ is reinitialized to a signed
0
0 0
distance function (methods for doing this are outlined in Chapter 7), then
’0.5 b φ can be used in place of bκ| φ| in the next time step as well. In
’0.5 ’0.5

summary, equations (4.5) and (4.6) have the same e¬ect on the interface
’1 ’1 ’1
1
0
’1 ’1 1
0 ’1 0 1
location as long as one keeps φ equal to the signed distance function o¬ the
interface. Note that keeping φ equal to signed distance o¬ the interface does
1 1 1
not change the interface location. It only changes the implicit embedding
0.5 0.5 0.5
function used to identify the interface location.
0 0 0

’0.5 ’0.5 ’0.5

’1 ’1 ’1
’1 1
0 ’1 0 1
4.2 Numerical Discretization
0
’1 1

1 1 1
Parabolic equations such as the heat equation need to be discretized using
0.5 0.5 0.5
central di¬erencing since the domain of dependence includes information
0 0 0
from all spatial directions, as opposed to hyperbolic equations like equa-
’0.5 ’0.5 ’0.5 tion (3.2), where information ¬‚ows in the direction of characteristics only.
Thus, the φ term in equation (4.6) is discretized using the second-order
’1 ’1 ’1
0 1
’1 1
’1 0 ’1 0 1
accurate formula in equation (1.9) in each spatial dimension (see equa-
tion (2.7)). A similar approach should therefore be taken in discretizing
Figure 4.2. Evolution of a star-shaped interface in a curvature-driven ¬‚ow. The
equation (4.5). The curvature κ is discretized using second-order accurate
tips of the star move inward, while the gaps in between the tips move outward.
central di¬erencing as outlined in equation (1.8) and the discussion follow-
ing that equation. Likewise, the φ term is discretized using the second
Equation (4.4) is also known as the equation of the level set equation. Like order accurate central di¬erencing in equation (1.5) applied independently
equation (3.2), equation (3.2) is used for externally generated velocity ¬elds, in each spatial dimension. While these discretizations are only second-order
while equation (4.4) is used for (internally) self-generated velocity ¬elds. accurate in space, the dissipative nature of the equations usually makes this
As we shall see shortly, this is more than a notational di¬erence. In fact, second-order accuracy su¬cient.
slightly more complicated numerical methods are needed for equation (4.4) Central di¬erencing of φ in equation (4.6) combined with a forward
than were proposed in the last chapter for equation (3.2). Euler time discretization requires a time-step restriction of
Plugging Vn = ’bκ into the level set equation (4.4) gives
2b 2b 2b
φt = bκ| φ|, (4.5) t + + <1 (4.7)
2 2 2
( x) ( y) ( z)
where we have moved the spatial term to the right-hand side. We note
to maintain stability of the numerical algorithm. Here t is O(( x)2 ),
that bκ| φ| is a parabolic term that cannot be discretized with an upwind
approach. When φ is a signed distance function, equation (4.5) becomes which is signi¬cantly more stringent than in the hyperbolic case, where
t is only O( x). Enforcing t = O(( x)2 ) gives an overall O(( x)2 )
the heat equation
accurate discretization, even though forward Euler is used for the time
φt = b φ, (4.6)
di¬erencing (i.e., since the ¬rst-order accurate O( t) time discretization
is O(( x)2 )). Equation (4.5) can be discretized using forward Euler time

<< . .

 4
( 23)



. . >>

Copyright Design by: Sunlight webdesign