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)
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  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 , 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  improved the WENO method by choosing is O(( x)5 ) in smooth regions. Note that  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 , following the work on HJ П‰3 = 0.3 + O( x).
ENO in , Jiang and Peng extended WENO to the Hamilton-Jacobi In order to deп¬Ѓne the weights, П‰k , we follow  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. ,
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 , 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  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 , 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  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)ОГЛАВЛЕНИЕ След. стр. >>