In this section, we will describe how to handle |∇φ| in Eq.(6). If we just use central difference to compute |∇φ| , then we will get terrible result.Before handling this term, let us look at 1D example :
Example.2 (Burger’s equation)
ut+ uux = 0 (10)
We discretize ux in Eq.(10) by upwind difference.We denote : Di+u = ui+1− ui
∆x and Di−u = ui− ui−1
∆x
Note that the coefficient u of ux decides the information from which direc-tion. That is, if u > 0 we know that the characteristics come from left side so we use D−i u to compute (ui)x. On the other hand, if u < 0,it means that the information and characteristics come from right side so we use D+i u to compute (ui)x.
The above equations explain what is upwind difference for 1D case.Note that Eq.(6) and Eq.(10) are similar for some sense,but Eq.(6) is a 2D PDE and is not as simple as Eq.(10) when handling first derivative term.
3.4.1 Godunov’s Method
Since Eq.(6) is a 2D PDE which is not as simple as Eq.(10), we introduce Godunov’s Method to handle |∇φ| in Eq.(6). Actually, Godunov’s method uses the same concept which is characteristic coming from which direction.
But in 2D domain you may doubt that the characteristic can come from infinite directions and how to decide? To overcome this problem, we just simplify and separate x-direction and y-direction to compute respectively that is just compute φx and φy respectively. Since we know the information propagated from interface, so when computing (φi,j)x, we choose the point of φi+1,j and φi−1,j by the direction of information propagating and then compute (φi,j)x. For example,if (i, j) is outside interface and if φi−1,j <
φi,j < φi+1,j which means (i − 1, j) is nearest interface within these three points and the information is propagated from φi−1,j to φi+1,j, so we use D−i,jφ = (φi,j−φ∆xi−1,j). On the other hand, if φi−1,j > φi,j > φi+1,j which means information comes from φi+1,j to φi−1,j, then we use D+i,jφ = (φi+1,j∆x−φi,j).
Having the basic concept of Godunov’s mothod, we find that in some sense this method is as upwind method. Here we describe Godunov’s method for solving Eq.(6) in the following :
Details of Godunov’s Method
1. If φi,j > 0 which means the point is outside interface, then
(a) If Di,j+φ > 0 and Di,j−φ > 0 then it means φi−1,j < φi,j < φi+1,j and (φi,j)x = Di,j−φ
(b) If Di,j+φ < 0 and Di,j−φ < 0 then it means φi−1,j > φi,j > φi+1,j and (φi,j)x = Di,j+φ
(c) If Di,j+φ ≤ 0 and Di,j−φ ≥ 0 then it means φi−1,j ≤ φi,j and φi,j ≥ φi+1,j. In this case, the information comes form both two sides and we treat the shock by the larger influence between D−i,j and Di,j+, that is (φi,j)x = max(|Di,j−|, |Di,j+|)
(d) If Di,j+φ ≥ 0 and Di,j−φ ≤ 0 then it means φi−1,j ≥ φi,j and φi,j ≤ φi+1,j and the information flow to both sides from (i, j). This does not make sense since interface do not pass grid (i, j). So in this case we just set (φi,j)x = 0.
2. If φi,j < 0 which means the point is inside interface.
(Note that inside interface the φ-value is negative and if φ is larger, it is much near the interface.)
(a) If Di,j+φ > 0 and Di,j−φ > 0 then it means φi−1,j < φi,j < φi+1,j and (φi,j)x = Di,j+φ.
(b) If Di,j+φ < 0 and Di,j−φ < 0 then it means φi−1,j > φi,j > φi+1,j and (φi,j)x = Di,j−φ
(c) If Di,j+φ ≥ 0 and Di,j−φ ≤ 0 then it means φi−1,j ≥ φi,j and φi,j ≤ φi+1,j. In this case, the information comes form both two sides and we treat the shock by the larger influence between D−i,j and Di,j+, that is (φi,j)x = max(|Di,j−|, |Di,j+|)
(d) If Di,j+φ ≤ 0 and Di,j−φ ≥ 0 then it means φi−1,j ≤ φi,j and φi,j ≥ φi+1,j and the information flow to both sides from (i, j). This does not make sense since interface do not pass grid (i, j). So in this case we just set (φi,j)x = 0.
After simplifying the Godunov’s method as above, the scheme is obvious and easy to implement by numerical computing. Here we give three remarks as follows:
Remark
1. In Godunov’s method, we discetize φx by one-sided difference, so the scheme is first order accurate.
2. The cases φ > 0 and φ < 0 are opposite since the signed function S(φ0) is 1 outside interface and −1 inside interface.
3. The y-direction case φy is as same as computing φx. So after solving φx and φy, the spatial discretization of |∇φ| in Eq.(6) is complete.
3.4.2 Second-Order Hamilton-Jacobi ENO Scheme
Since the Godunov’s method is first order, if we want the scheme more pre-cise, we introduce second-order ENO (essential non-oscillatory) scheme to discretize. The ENO scheme is based on Godunov’s method for first-order term which described above when solving Eq.(6). We write the Taylor ex-pansion of φ:
φ(x) = φ(˜x) + φ0(˜x)(x − ˜x) + φ00(˜x)
2! (x − ˜x)2+ · · · (11) If we use D−φi, then it means we plug xi and xi−1 into Eq.(11) and get (φi)x. Now the ENO scheme is to add a adjacent point plug into Eq.(11).
The point is chosen to be much smooth between two sides. Let us take a example to explain :
Example.3 If Godunov’s method tell us that we choose the point xi−1 when computing (φi)x then we need to choose another point between xi+1
and xi−2 to plug into Eq.(11) and get second-order result. We want choose much smooth side point so we need some criterion for choosing point. Let us denote: and we choose the point xi−2to plug into Eq.(11). Otherwise we choose xi+1
to plug into Eq.(11). No matter what point is chosen, the scheme is always second-order.
Now we describe the whole process for 1D case for ENO scheme based on Godunov’s method :
1. If we use D−φ for first-order discretization to compute (φi)x
Figure 3: (a).left-biased (b).right-biased
To handle 2D case, it is as the same method and we split the domain into x-direction and y-direction to handle respectively. Here we have intro-duced the second-order ENO scheme for solving Hamilton-Jacobi equation like Eq.(6). Next we will introduce more higher order scheme.
3.4.3 Third-Order ENO Scheme
The third-order ENO scheme for spatial is based on second-order ENO scheme which is based on Godunov’s method, so we realize that the concept of Godunov method is very important when solving Eq.(6). In second-order ENO, we use three points to calculate φx and get second-order result, the third-order is as the same concept as ENO-2. While deciding which three points are needed, we still demand one point to reach third order ENO(ENO-3) which is either left side or right side adjacent point. The fourth point is also chosen in smoother side and the criterion is similar to above ENO-2.
Since the concept is similar, we do not take example again and we just give the details of algorithm as below:
1. If we use D−φ for first-order discretization to compute (φi)x
2. If we use D+φ for first-order discretization to compute (φi)x
It seems that there are eight cases for ENO-3, but actually it can be reduced to just three cases. That is using {xk|k = i − 3, i − 2, i − 1, i} or {xk|k = i − 2, i − 1, i, i + 1} {xk|k = i − 1, i, i + 1, i + 2} for D− case to plug into Eq.(11). Similarly the right-biased stencil case for D+ is {xk|k = i, i+1, i+2, i+3} or {xk|k = i−1, i, i+1, i+2} or {xk|k = i−2, i−1, i, i+1}.
Finally we get just three discretization forms.
3.4.4 Third-Order WENO
In this section we will discuss the third-order WENO(weighted essential non-oscillatory) scheme to discretize φx. Again W3 is based on ENO-2.ENO-2 has two possibilities to compute φx under D− (See Eq.(13) and
Eq.(15))and we combine this two by weights instead using only one to com-pute φx. That is why we call this method Weighted ENO.
The scheme of WENO-3 is more easy than ENO-3 since we only need to judge D− or D+. If we need D− for first-order then we use left-biased points {xk|k = i − 2, i − 1, i, i + 1} to compute φx. On the other case we use right-biased points {xk|k = i − 1, i, i + 1, i + 2} to compute φx. The process of WENO is more simple since we do not need to judge which points are needed to add.
We write down the details of WENO-3 which we refer [12]:
1. If D−φ is needed then
Note that the form is more simple and is easy to implement, so we rec-ommend the WENO-3 scheme to discretize the φx and φy of reinitialization PDE.
Remark
Although the sign function of Eq.(6) is just simple 0,1 or −1, [12] suggested that
S(φ) = φ
pφ2+ |∇φ|2(∆x)2 (20)
for a better result and reduce the time steps when solving Eq.(6) to steady state.
3.5 Temperal Discretization for Reinitialization Pro-cess
Since we prefer to use WENO-3 to discretize |∇φ| for spatial, if we want the scheme is higher order, we also need third-order method in time. Here we use TVD Runge-Kutta method which we refer [1].
3.5.1 Second-Order TVD Runge-Kutta Method
We let current time step be n and first we move time to n + 2,then take average as below
φn+1= 1
2φn+ 1 2φn+2 3.5.2 Third-Order TVD Runge-Kutta Method
For third-order TVD-RK scheme we first need φn+12 which is weighted by φn and φn+2 and so that
φn+12 = 3
4φn+ 1 4φn+2
Next we move φn+12 to φn+32 by iteration and finally we get φn+1 as follows:
φn+1 = 1
3φn+2 3φn+32
Note that if ENO-2 is used for spatial, then TVD-RK2 have better been used for time and both ENO-3 and WENO-3 match TVD-RK3.
3.6 Numerical Test for Reinitialization Process
In this section, we will test some examples for reinitialization. We choose above example φ(x, y) = p(0.6x )2+ (0.3y )2 − 1 and a concave case which we refer [12] with polar coordinate φ(r, θ) = (r − 0.5 + 0.1r sin(7θ))3.We contour [−10∆x : 2∆x : 10∆x] for each case. Note that we use WENO-3 and TVD-RK3 method for these examples.
Figure 4: Reinitialization of ellipse
Figure 5: Reinitialization of plum blossom shape
Remark
1. From the figures of two examples, we find that the level function be-comes signed distance function from interface to Ω+ and Ω− within a neighborhood of interface. As the time steps increase, the width of neighborhood also increase. The reason is that the information flows from interface to Ω+ and Ω−.
2. After reinitialization process, the level function becomes signed distance function in a neighborhood of interface and we can not write the level function by mathematical formulation exactly.
3. There is no clear condition of stopping criterion for reinitialization pro-cess. The iteration steps is determined by the width of signed distance function neighborhood of interface. If you want the width wider, the iteration steps needed more.
4. Eq.(6) is solved in whole domain and the above two examples are lo-cated at center region of computation domain, so if we just want a signed distance function of just a neighborhood of interface, we do not need the boundary condition.But if the interface contacts with the boundary, how do we set the boundary condition for Eq.(6)? We will give the idea and algorithm in next section.