LINEBURG

φ(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

Copyright Design by: Sunlight webdesign