Effects of three-dimensional surface corrugations on electromagnetohydrodynamic flow through microchannel

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization ∗

Wenbin Li a , Jianliang Qian b , a b

School of Science, Harbin Institute of Technology, Shenzhen, Shenzhen 518055, China Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA

article

info

a b s t r a c t

Article history: Received 15 April 2019 Received in revised form 9 August 2019 Accepted 28 August 2019 Available online xxxx Keywords: Generalized eikonal equation Newton’s method Gauss–Seidel iterations Lax–Friedrichs scheme

We propose a Newton-type Gauss–Seidel Lax–Friedrichs sweeping method to solve the generalized eikonal equation arising from wave propagation in a moving fluid. The Lax– Friedrichs numerical Hamiltonian is used in discretization of the generalized eikonal equation. Different from traditional Lax–Friedrichs sweeping algorithms, we design a novel approach with a line-wise sweeping strategy. In the local solver, the values of traveltime on an entire line are updated simultaneously by Newton’s method. The global solution is then obtained by Gauss–Seidel iterations with line-wise sweepings. We first develop the Newton-based first-order scheme, and on top of that we further develop high-order schemes by applying weighted essentially non-oscillatory (WENO) approximations to derivatives. Extensive 2-D and 3-D numerical examples demonstrate the efficiency and accuracy of the new algorithm. The combination of Newton’s method and Gauss–Seidel iterations improves upon the convergence speed of the original Lax–Friedrichs sweeping algorithm. In addition, the Newton-type sweeping method manipulates data in a vectorized manner so that it can be efficiently implemented in modern programming languages that feature array programming, and the resulting advantages are extremely significant for large-scale 3-D computations. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction We study the generalized eikonal equation with the following form:

∥∇ T (x)∥ −

|1 − v(x) · ∇ T (x)| F (x)

= 0, x ∈ Ω ;

T (x) = g(x), x ∈ Γ .

(1.1)

It models the traveltime T (x) of wave propagation in an isotropic acoustic medium occupied by a moving fluid [1]. F (x) denotes the speed of wave propagation, v(x) is the velocity of the moving fluid, and Γ is a subset of the domain Ω where the boundary condition is imposed. This equation plays an important role in a variety of applications. For example, in [2], a generalized eikonal equation is derived to model the phase function of the geometric-optics ansatz for the wave field of steady wind flow over mountains; in [3], it shows that the Class 1C fold in structural geology can be modeled by a modified generalized eikonal equation; and in [4,5], the generalized eikonal equation appears in the study of sonic boom propagation from cruising aircrafts. ∗ Corresponding author. E-mail addresses: [email protected] (W. Li), [email protected] (J. Qian). https://doi.org/10.1016/j.camwa.2019.08.031 0898-1221/© 2019 Elsevier Ltd. All rights reserved.

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

2

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Eq. (1.1) is a static Hamilton–Jacobi equation, and a classic approach for computing the traveltime T (x) is the Lagrangian ray tracing method [6,7]. Here, a system of ordinary differential equations is solved to integrate the characteristics of the Hamilton–Jacobi equation, resulting in a group of rays along which the traveltime is evaluated. However, such Lagrangian methods yield the solution in the physical space which might not be uniformly distributed, leading to shadow regions with few passing rays. An alternative approach is to work in the Eulerian framework, in which the partial differential equation is solved directly on mesh grids. Osher has developed a level set formulation to solve the static Hamilton–Jacobi equations [8]. One can embed the characteristic system of (1.1) into the evolution of a level set function, and solve a time-dependent Hamilton–Jacobi equation with one dimension higher. On the other hand, recent techniques solve the eikonal equation in physical space as a stationary boundary value problem, and among them two main techniques are fast marching method [9,10] and fast sweeping method [11–17]. In the fast marching method, the update of the solution follows the causality or the characteristics in a sequential way. Once a solution is updated, it means that the wave front has crossed the grid point and the value will not change in further updating. As a result, the fast marching method must enforce a strict causality principle to sort the order of computations. For the generalized eikonal equation (1.1), the front of the traveltime may propagate in both tangential and normal directions, and the updating scheme must be carefully developed to respect the characteristic direction. In [18], the authors propose a characteristic fast marching method for Eq. (1.1). The algorithm exhaustively studies the characteristic direction at each updating, and thus chooses the correct finite difference approximation in the scheme. The deficiency of the algorithm is that a CFL-like condition in the finite difference scheme is unavailable, which can make the computation have unpredictable behaviors [18]. In the fast sweeping method, the solution is updated in an alternating sweeping ordering. The characteristics are divided into a finite number of groups according to their directions, and each sweeping ordering covers a group of characteristics simultaneously [11,19]. With alternating ordering, the fast sweeping method can follow causality along characteristics as well as handle global relaxation and smoothing effects [14]. Even if there is no strict causality, such as the Hamiltonian is not convex, or viscosity terms are involved, the iteration of fast sweeping method may still converge. In [11], the author proposes a fast sweeping method for the classical eikonal equation, where a Godunov upwind difference scheme is used to update the traveltime in local stencil. For the generalized eikonal equation (1.1), however, it can be extremely involved to carry out the min–max optimization resulting from the Godunov scheme. In [13], the authors propose a Lax–Friedrichs sweeping algorithm for the general static Hamilton–Jacobi equations. The Lax–Friedrichs scheme is used to discretize the Hamiltonian, and a simple updating formula is derived for the fast sweeping method. This algorithm is further developed for high-order computation in [14], where the weighted essentially non-oscillatory (WENO) approximation [20] is used to achieve high-order accuracy. The Lax–Friedrichs sweeping algorithm can simply carry out the local solution process for inhomogeneous and non-convex Hamiltonians; therefore one can apply it to the solution of the generalized eikonal equation. In [21], the authors propose various operator-splitting based fast sweeping methods for Eq. (1.1), where the Lax–Friedrichs sweeping algorithm is adopted to solve the equations after splitting. In this work, we propose a Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm for the solution of the generalized eikonal equation. The Lax–Friedrichs scheme is adopted to construct the numerical Hamiltonian. Different from traditional Lax–Friedrichs sweeping algorithms, we propose a novel approach with a line-wise sweeping strategy. In the typical fast sweeping method [11,13–15,22], the solution is updated point by point and Gauss–Seidel iterations with pointwise sweepings are performed to achieve global convergence. In the new algorithm, we take advantage of the simplicity of the Lax–Friedrichs scheme and update the values of traveltime on an entire line simultaneously by Newton’s method. To avoid solving a very large linear system, the traveltime values outside the updating line are treated explicitly when applying Newton’s method. Gauss–Seidel iterations with line-wise sweepings are then performed to search for the global solution. It turns out that the linear system generated by Newton’s method is tridiagonal when computing the values of traveltime along a line, and hence the Thomas algorithm [23] is applicable to solving it quickly. In a nutshell, the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm applies Newton’s method to solve a block of the global solution, and further sweeps the entire computational domain in a block-wise manner to obtain the global solution. As we will see from numerical examples, this strategy improves upon the convergence speed of the sweeping method. In addition, since the solution is updated in a line-wise rather than pointwise manner, the new algorithm is better vectorized and is suitable to modern programming languages that feature array programming. In Section 2, we introduce the Lax–Friedrichs scheme for the numerical Hamiltonian, and briefly review the traditional pointwise sweeping method as introduced [13,14]. In Section 3, we develop the local solver and sweeping strategy for the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. In Section 4, the algorithm is extended to construct a high-order scheme by applying the WENO approximation to derivatives. Numerical examples are shown in Section 5 to illustrate the efficiency and accuracy of the new algorithm. Lastly in Section 6 we draw our conclusion. 2. Lax-Friedrichs schemes Eq. (1.1) is a specific form of the general static Hamilton–Jacobi equation

{

H(x, ∇ T (x)) = f (x), T (x) = g(x),

x∈Ω x ∈ Γ.

(2.1)

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

3

The Lax–Friedrichs scheme for the numerical Hamiltonian reads as follows in the 2-D space [24],

ˆ LF (u+ , u− ; v + , v − ) = H H

(

u+ + u− v + + v −

,

2

)

2

1

1

2

2

− α x (u+ − u− ) − α y (v + − v − )

(2.2)

where u := ∂ T /∂ x and v := ∂ T /∂ y with u± and v ± their corresponding forward and backward finite-difference approximations, and α x and α y are artificial viscosities satisfying

⏐ ⏐ ⏐ ∂H ⏐ α x ≥ max ⏐⏐ ⏐⏐ , ∂u

⏐ ⏐ ⏐ ∂H ⏐ α y ≥ max ⏐⏐ ⏐⏐ . ∂v

(2.3)

Based on this numerical Hamiltonian, Kao, Osher and Qian proposed a first-order Lax–Friedrichs sweeping algorithm for the static Hamilton–Jacobi equation [13],

( H

Ti+1,j − Ti−1,j Ti,j+1 − Ti,j−1

,

2∆x 2∆y Ti+1,j − 2Ti,j + Ti−1,j

−α · x

2∆x

)

− αy ·

Ti,j+1 − 2Ti,j + Ti,j−1 2∆y

= fi,j .

(2.4)

Consequently, T is solved pointwise with the following updating formula:

( Ti,j =

+ αx ·

)[

1 αx ∆x

+

(

fi,j − H

αy ∆y

Ti+1,j + Ti−1,j 2∆x

+ αy ·

Ti+1,j − Ti−1,j Ti,j+1 − Ti,j−1

,

2∆x Ti,j+1 + Ti,j−1

]

2∆y

)

2∆y

,

(2.5)

and Gauss–Seidel iterations with alternating-direction sweepings are applied to obtain the global solution. In [14], Zhang, Zhao and Qian applied WENO schemes to approximate derivatives and extend the algorithm to high-order computations. Since the traveltime T is updated point by point in Gauss–Seidel iterations, we call this type of algorithm the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm or simply pointwise sweeping in the context of this paper. 3. Newton-type Gauss–Seidel sweeping methods We propose a Newton-type Gauss–Seidel sweeping method for the Lax–Friedrichs based numerical discretization system (2.4), and apply it to solve the generalized eikonal equation. Instead of updating Ti,j pointwise, the new method updates the solution in a line-wise manner. The algorithm is composed of two parts: the Newton-type local solver and the Gauss–Seidel iterations with line-wise sweepings for the global solution. For simplicity, we describe the algorithm in the two-dimensional space with mesh grids {xi,j , i, j = 1, 2, . . . , N }. We mention that the formulations can be directly extended to high-dimensional computations. 3.1. Column-wise Newton solver Instead of solving Ti,j pointwise, the column-wise solver looks for the jth column of the solution {Ti,j , i = 1, 2, . . . , N } simultaneously from the Lax–Friedrichs scheme. The discretization in (2.4) is then rewritten in the following form, Fi,j Ti∗−1,j , Ti∗,j , Ti∗+1,j

(

− αx ·

Ti∗+1,j

+

)

Ti∗−1,j

2∆x

( =H − αy ·

Ti∗+1,j − Ti∗−1,j Ti,j+1 − Ti,j−1 2∆x Ti,j+1 + Ti,j−1 2∆y

,

2∆y

( +

)

αx αy + ∆x ∆y

)

· Ti∗,j − fi,j = 0,

(3.1)

where we use the superscript ‘∗’ to indicate the variables to be solved at this stage. In (3.1), the other variables without ‘∗’ are treated explicitly so that the jth column of the solution can be computed at a reasonable cost. The entire solution of Ti,j at all grid points will be achieved by a line-wise sweeping algorithm, which will be discussed in the next subsection. We mention that although we propose the line-wise solver in a column-wise manner here, one can also consider a rowwise solver and solve the ith row of the solution, which can be equivalent to the column-wise solver in the sense of interchanging i and j. We propose to apply Newton’s method to solve Eq. (3.1). The interior points of the jth column, {Ti,j , i = 2, . . . , N − 1}, are governed by a system of (N − 2) equations,

( ) ⎧ F2,j T1,j , T2∗,j , T3∗,j ⎪ ( ) ⎪ ⎪ ∗ ∗ ∗ ⎪ ⎪ ⎨ F3,j T2,j , T3,j , T4,j .. . ( ⎪ ) ⎪ ⎪ ⎪ FN −2,j (TN∗−3,j , TN∗−2,j , TN∗−)1,j ⎪ ⎩ FN −1,j TN∗−2,j , TN∗−1,j , TN ,j

= 0 = 0 (3.2)

= 0 = 0.

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

4

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Denoting Tj = (T2,j , T3,j , . . . , TN −1,j )⊺ and Fj = (F2,j , F3,j , . . . , FN −1,j )⊺ , Newton’s method for the nonlinear system Fj (Tj ) = 0 has the following updating formula:

]−1

· Fj , [ ] where ∇ Fj denotes the Jacobian matrix with the following form, ⎛ ∂F ⎞ ∂F − ∇ Fj = Told Tnew j j

[

2,j

2,j

⎜ ⎜ ⎜ [ ] ⎜ ⎜ ∇ Fj = ⎜ ⎜ ⎜ ⎜ ⎝

(3.3)

∂ T2,j

∂ T3,j

∂ F3,j ∂ T2,j

∂ F3,j ∂ T3,j

∂ F3,j ∂ T4,j

.

.

..

..

..

∂ FN −2,j ∂ TN −3,j

.

∂ FN −2,j ∂ TN −2,j

∂ FN −2,j ∂ TN −1,j

∂ FN −1,j ∂ TN −2,j

∂ FN −1,j ∂ TN −1,j

⎟ ⎟ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎟ ⎠

(3.4)

The Jacobian matrix is tridiagonal due to the specific form of Fi,j , and it can be evaluated explicitly from Eq. (3.1), a2 ⎜ b3

c2 a3

⎛

..

[ ] ⎜ ∇ Fj = ⎜ ⎜ ⎝

.

⎞ c3

..

..

.

bN −2

.

aN −2 bN −1

cN −2 aN −1

⎟ ⎟ ⎟ ⎟ ⎠

(3.5)

with

⎧ y ∂F αx ⎪ + ∆α y := A, ai = ∂ Ti,j = ∆ ⎪ x i ,j ⎪ ⎨ ( ) x ∂F bi = ∂ T i,j = ∂∂Hu i,j · 2−∆1x − 2α∆x , i − 1 , j ⎪ ⎪ ⎪ ⎩ c = ∂ Fi,j = ( ∂ H ) · 1 − αx , i ∂ Ti+1,j ∂ u i,j 2∆x 2∆x ( ∂H )

i = 2, . . . , N − 1 i = 3, . . . , N − 1

(3.6)

i = 2, . . . , N − 2.

In Eq. (3.6), ∂ u i,j is a function of ui,j and vi,j , where ui,j and vi,j are finite-difference approximations of (∂ T /∂ x)i,j and (∂ T /∂ y)i,j ; for example, we can take ui,j =

Ti+1,j − Ti−1,j 2∆x

,

vi,j =

Ti,j+1 − Ti,j−1 2∆y

[

.

(3.7)

]

Proposition 3.1. The Jacobian matrix ∇ Fj is strictly diagonally dominant thus invertible. Proof. The viscosity parameters α x and α y are chosen according to formula (2.3), therefore

) ) ( ( ) (( ) 1 ∂H ∂H x x − α ≤ 0 and ci = − α ≤ 0, − bi = 2∆x ∂ u i ,j 2∆x ∂ u i,j 1

α which implies that |bi | + |ci | = −(bi + ci ) = ∆ < A = ∆α x + ∆α y = |ai |. The matrix ∇ Fj is strictly diagonally dominant x and thus invertible according to the Levy–Desplanques theorem. □ x

x

y

[

]

The updating formula in (3.3) is implemented with two steps: 1. solve the perturbation term δT j from

[ ] ∇ Fj · δT j = Fj ;

(3.8)

2. update Tj by Tnew = Told − δT j . j j

(3.9)

Since Eq. (3.8) is a tridiagonal linear system and the coefficient matrix is strictly diagonally dominant, one can solve it quickly using the Thomas algorithm [23]. On the computational boundary, we adopt the boundary condition proposed by Kao et al. in [13], which performs linear extrapolation if there is no source point, and imposes a zero Neumann boundary condition if a source point appears on the computational boundary. The detailed formulas in the 2-D case are given below:

⎧ T1,j = ⎪ ⎪ ⎨ T N ,j = Ti,1 = ⎪ ⎪ ⎩ Ti,N =

max (2T2,j − T3,j , T3,j ) max (2TN −1,j − TN −2,)j , TN −2,j max (2Ti,2 − Ti,3 , Ti,3 ) max 2Ti,N −1 − Ti,N −2 , Ti,N −2 .

(

)

(3.10)

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

5

Eqs. (3.8), (3.9) and (3.10) complete the Newton-type local solver for the Lax–Friedrichs scheme (3.1). When it is applied to the generalized eikonal equation (1.1), the Hamiltonian has the specific form, H(u, v ) =

√

u2 + v 2 −

|1 − v1 u − v2 v| F (x)

(3.11)

where (u, v ) = ∇ T and (v1 , v2 ) = v(x). It implies that v1 ∂H u + = √ · sgn(1 − v1 u − v2 v ), (3.12) 2 2 ∂u F (x) u +v v v2 ∂H = √ + · sgn(1 − v1 u − v2 v ), (3.13) ∂v F (x) u2 + v 2 where sgn(·) denotes the signum function. The viscosity parameters α x and α y satisfying (2.3) can be set as follows, ⏐v ⏐ ⏐v ⏐ ⏐ 1⏐ ⏐ 2⏐ α x = 1 + max ⏐ ⏐ , α y = 1 + max ⏐ ⏐ . (3.14) F F √ In the numerical discretization one may add a small positive parameter to u2 + v 2 to avoid zero denominator in Eqs. (3.12) and (3.13). 3.2. Gauss–Seidel iterations with line-wise sweepings In the Newton-type local solver, the jth column of the solution {Ti,j , i = 1, 2, . . . , N } is solved from Eq. (3.1) using its neighboring values Ti,j+1 and Ti,j−1 , and in this procedure the latest values of T are utilized if available. This is similar to the strategy of the Gauss–Seidel method for linear systems, and moreover, we can use an iterative approach to search for the global solution of this nonlinear system (3.1). The basic procedure is to sweep the whole domain with two alternating orderings repeatedly: from left to right and from right to left; at every sweeping the solution is updated column by column according to Eqs. (3.8)–(3.10). The correct values of Ti,j are obtained at all grid points as the iteration converges. This iterative approach is similar to those proposed in [13,14], and the difference is that in the alternating sweepings the values of T are updated in a column-wise rather than point-wise manner. In the following, we list the detailed steps of this algorithm. 1. Initialization. To enforce the source condition, T (x) = g(x) for x ∈ Γ ⊂ Ω , assign exact values or interpolated values at grid points on or near Γ , and these values are fixed in later computations. Assign large positive values M at all other grid points, where M should be larger than the maximum of the true solutions, and these values will be updated in the iterations. 2. Newton solver with alternating sweeping orderings. At each column, update solutions by the column-wise Newton solver; sweep the whole domain with two alternating directions repeatedly: (1) j = 1 : N (left to right),

(2) j = N : 1 (right to left).

We use a pseudo code to illustrate the detailed operations in one sweeping. 1: 2: 3: 4: 5: 6: 7: 8: 9: 10: 11:

for j = 2 to N − 1 do New tonError = 1; New tonIter = 0; while New tonError > ϵ0 and New tonIter < MaxIter do compute Tnew = (T2,j , T3,j , . . . , TN −1,j )⊺ according to Eqs. (3.8), (3.9); j New tonError = ∥δ Tj ∥l∞ ; compute T1,j and TN ,j according to Eq. (3.10); New tonIter = New tonIter + 1; end while ( ) {( ) ( n )⊺ } 1 ⊺ new ⊺ update T1n,+j 1 , . . . , TNn+ = min T1new , T1,j , . . . , TNn,j ; ,j ,j , . . . , TN ,j end for ( ) ( ) 1 ⊺ new ⊺ j = 1 and j = N, compute T1new according to Eq. (3.10) and update T1n,+j 1 , . . . , TNn+ = ,j ,j ), .}. . , TN ,j {( new ) ( ⊺ n n ⊺ min T1,j , . . . , TNnew , T , . . . , T . ,j N ,j 1,j

In line 3, ϵ0 denotes the tolerance error of the Newton solver, and MaxIter denotes the maximum number of iterations in the local solver. In line 5, ∥ · ∥l∞ denotes the vector l∞ -norm. It is unnecessary to achieve convergence in the Newton solution of the jth column in every sweeping, since the local solution is non-essential before the correct value is propagated to its neighbors. In our application, we set MaxIter = 1 and only perform one iteration in the Newton-type local solver. Numerical tests suggest that such strategy can significantly improve computational efficiency of this algorithm. 3. Convergence test. Given convergence criterion ϵ > 0, check whether

∥Tn+1 − Tn ∥l∞ ≤ ϵ. Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

6

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

3.3. Comparing with the pointwise sweeping algorithm As mentioned in Section 2, we use the name ‘pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm’ to refer to the first-order and high-order Lax–Friedrichs schemes proposed in [13,14], since both of them update Ti,j in a pointwise manner. The newly designed Newton-type sweeping algorithm inherits the Lax–Friedrichs discretization and the strategy of Gauss–Seidel iterations except that the values of T are updated in a line-wise manner, and this difference changes the nature of the algorithm. The pointwise sweeping algorithm can be viewed as an explicit scheme because the value of Ti,j is evaluated explicitly through Eq. (2.5). On the other hand, the Newton-type sweeping algorithm has a flavor of semi-implicit scheme. When updating the jth column of T , the column values are treated implicitly and solved by a Newton approach, while the neighboring values of T are treated explicitly. The algorithm thus avoids a direct application of Newton’s method to the entire system of the discretized generalized eikonal equation, which will generate a very large system with a non-sparse coefficient matrix. This semi-implicit strategy improves upon the convergence speed of the original pointwise sweeping algorithm and reduces the size of the to-be-solved linear system when applying Newton’s method. As we will see from numerical examples, the Newton-type sweeping algorithm usually takes much smaller numbers of sweepings to achieve convergence comparing to the pointwise sweeping algorithm. As discussed in the preceding subsection, we perform only one step of Newton’s iteration when updating the column values of T . Since the complexity of Thomas algorithm is O(N), the computational cost of updating one column of T is also O(N), and the operations in one full sweeping are O(N 2 ) for 2-D cases. This complexity is the same as that of pointwise Gauss–Seidel Lax–Friedrichs sweeping method, although the proportionality constant in front of N 2 can be different. Furthermore, the Newton-type sweeping algorithm manipulates data in a vectorized manner; namely, the values of T on an entire column are updated simultaneously by vector operations. The vectorized manner of the new algorithm makes it more suitable for modern programming languages that are designed to feature array programming, such as Fortran 90, Matlab, Octave, R and Python. The feature of vectorization reduces the computational time as the algorithm is implemented by array programming. Such advantages turn out to be significant when the computational cost of each updating becomes large, in which situation the effect of vectorization in time-saving becomes substantial. In summary, the Newton-type sweeping algorithm has two advantages comparing to the traditional pointwise sweeping algorithm. Firstly, the newly designed algorithm takes a smaller number of sweepings to achieve convergence; secondly, it manipulates the solution in a vectorized manner, which is amenable to array programming and helps to save computational time. As we will see in the numerical examples, the Newton-type sweeping algorithm outperforms the traditional pointwise sweeping algorithm in computational time when the mesh size being very fine and/or when involving 3-D computations. This feature suggests that the Newton-type sweeping algorithm is appropriate for the application of high-order schemes, where the computational cost of every updating step increases substantially. 4. High-order schemes In this subsection we illustrate the construction of high-order schemes by using the third-order WENO approximation [20] to derivatives. The skeleton of the Newton-type sweeping algorithm does not change, and it is still composed of the Newton-type local solver and the Gauss–Seidel iterations with line-wise sweepings. 4.1. Newton-type local solver

)⊺

The Newton-type local solver for the jth column T1,j , . . . , TN ,j is still completed by Eqs. (3.8)–(3.10), while the difference is that the WENO scheme should be incorporated into the discretization of derivatives. Firstly, in Eq. (3.8), the right-hand side Fj = (F2,j , F3,j , . . . , FN −1,j )⊺ is evaluated according to the following formula,

(

ˆ LF (u+ , u− ; v + , v − ) − fi,j Fi,j = H i,j i,j i,j i,j ˆ LF is the numerical Hamiltonian as shown in Eq. (2.2); therefore where H

( Fi,j = H

− + − u+ i,j + ui,j vi,j + vi,j

2

,

2

)

1

1

2

2

− − α x (u+ α y (vi+,j − vi−,j ) − fi,j i , j − ui , j ) −

(4.1)

+ − − + + − − with u+ i,j = (∂ T /∂ x)i,j , ui,j = (∂ T /∂ x)i,j , vi,j = (∂ T /∂ y)i,j , and vi,j = (∂ T /∂ y)i,j . The derivatives are approximated by third-order WENO schemes, and the details are as follows:

⎧ ( )+ ( ) ( ) ∂T Ti+1,j − Ti−1,j −Ti+2,j + 4Ti+1,j − 3Ti,j ⎪ + x x ⎪ ⎪ u = = (1 − w+ ) + w+ ⎪ ⎨ i,j ∂ x i,j 2∆x 2∆x ( )2 ⎪ ϵ + Ti+2,j − 2Ti+1,j + Ti,j 1 ⎪ x ⎪ wx = ⎪ ( x ) 2 , r+ = ( )2 ⎩ + 1 + 2 r+ ϵ + Ti+1,j − 2Ti,j + Ti−1,j

(4.2)

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

7

⎧ ( ( ) ) ( )− Ti+1,j − Ti−1,j 3Ti,j − 4Ti−1,j + Ti−2,j ∂T ⎪ − x x ⎪ ⎪ ) = (1 − w + w u = − − ⎪ ⎨ i,j ∂ x i,j 2∆x 2∆x ( )2 ⎪ ϵ + Ti,j − 2Ti−1,j + Ti−2,j 1 ⎪ x x ⎪ ⎪ ( x )2 , r− = ( )2 ⎩ w− = 1 + 2 r− ϵ + Ti+1,j − 2Ti,j + Ti−1,j

(4.3)

⎧ ( ) ) ( )+ ( Ti,j+1 − Ti,j−1 ∂T −Ti,j+2 + 4Ti,j+1 − 3Ti,j ⎪ y y + ⎪ ⎪ = (1 − w ) + w v = + + ⎪ ⎨ i,j ∂ y i,j 2∆y 2∆y ( )2 ⎪ ϵ + Ti,j+2 − 2Ti,j+1 + Ti,j 1 ⎪ y y ⎪ ⎪ w+ = ( y )2 , r+ = ( )2 ⎩ 1 + 2 r+ ϵ + Ti,j+1 − 2Ti,j + Ti,j−1

(4.4)

⎧ ( )− ( ) ) ( ∂T Ti,j+1 − Ti,j−1 3Ti,j − 4Ti,j−1 + Ti,j−2 ⎪ y y − ⎪ ⎪ v = = (1 − w− ) + w− ⎪ ⎨ i,j ∂ y i,j 2∆y 2∆y )2 ( ⎪ ϵ + Ti,j − 2Ti,j−1 + Ti,j−2 1 ⎪ y y ⎪ w− = ⎪ )2 ( y )2 , r− = ( ⎩ 1 + 2 r− ϵ + Ti,j+1 − 2Ti,j + Ti,j−1

(4.5)

At the grid points near the computational boundary, there are no enough grid points for constructing the WENO approximation; consequently, we use first-order finite differences therein:

⎧ TN ,j − TN −1,j ⎪ , ⎨ u+ N −1,j = ∆x Ti,N − Ti,N −1 ⎪ ⎩ v+ , i,N −1 = ∆y

vi−,2

T2,j − T1,j

, ∆x Ti,2 − Ti,1 = . ∆y

u− 2,j =

(4.6)

Eqs. (4.1)–(4.6) complete the computation of Fj in Eq. (3.8). Secondly, the coefficient matrix [∇ Fj ] in Eq. (3.8) should also be constructed ( ) by the WENO approximation. Recall that the matrix [∇ Fj ] is determined by Eqs. (3.5) and (3.6), where one needs ∂∂Hu i,j to evaluate entries of the matrix. Since

( ∂H )

∂ u i,j

is a function of ui,j and vi,j , we propose a high-order scheme to evaluate their values,

ui,j =

− u+ i,j + ui,j

2

,

vi,j =

vi+,j + vi−,j 2

,

(4.7)

± where u± i,j and vi,j are finite difference approximations as shown in Eqs. (4.2)–(4.6). The high-order scheme does not change the structure of the coefficient matrix; namely, [∇ Fj ] is still tridiagonal and strictly diagonally dominant so that Eq. (3.8) can still be solved by the Thomas algorithm.

4.2. Line-wise sweepings The application of high-order schemes makes slight differences in the Gauss–Seidel iterations with line-wise sweepings. Firstly, in the initialization, one should assign exact values or interpolated values in a neighborhood of the source location Γ , since the high-order scheme needs more stencils for discretization. In the third-order WENO approximation, we assign exact values or interpolated values at grid points whose distances to Γ are less than or equal to 2h, where h = ∆x = ∆y is the mesh size. Secondly, when updating the jth column of the solution, we do not impose monotonicity in the high-order scheme; that is, we do not take Tn+1 = min {Tnew , Tn } but simply Tn+1 = Tnew . Moreover, since the high-order scheme is not monotonic and can be less stable, we use the solution of the first-order scheme to initialize the grid points outside the neighborhood of the source location, and this process leads to a better convergence. 5. Numerical examples Although there is no rigorous proof for the global convergence of the Newton-type Lax–Friedrichs sweeping algorithm at this moment, our numerical experiments show that this algorithm works very well for solving generalized eikonal equations at large-scale discretization. In this section, we provide numerical examples to demonstrate the effectiveness of first-order and high-order schemes in both 2-D and 3-D computations. We program the algorithm in Matlab and the computation is performed on an equivalent 2.5 GHz processor. Since the Newton-type Lax–Friedrichs sweeping algorithm updates solution in a vectorized manner, it is natural to implement it by array programming that manipulates data by operations on vectors. In the Matlab programming, we use vector operations Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

8

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. (Example 1) Convergence of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping method. Both the x-axis and y-axis are in the logarithmic scale. (a) First-order scheme. (b) High-order scheme with third-order WENO approximations.

instead of ‘for’ iterations whenever possible. For the prospective users of the Newton-type sweeping algorithm, we also suggest that the algorithm should be programmed in such a vectorized way; otherwise, the coding of the program can be fussy and the implementation of the algorithm will be inefficient. The first step of the sweeping algorithm is to initialize the values of traveltimes T (x). As we have discussed, the boundary condition is imposed by assigning exact or interpolated values of T at grid points on or near Γ , but the initialization strategy away from Γ for the first-order scheme is different from that for the high-order scheme. For the first-order scheme, we initialize the traveltime T to be T = 108 away from the source location Γ . One can certainly apply a better value of initialization which is closer to the exact solution and improve the convergence speed. Here we use the very large initial value T = 108 in order to demonstrate the robustness of the Newton-type sweeping algorithm. For the high-order scheme, the first-order solution is used to initialize T away from the source location. The sweeping algorithm is stopped when the convergence criterion is achieved: ∥Tn+1 − Tn ∥l∞ ≤ ϵ . We set ϵ = 10−8 in the first-order scheme and ϵ = 10−9 in the high-order scheme. Moreover, the small parameter ϵ in the third-order WENO discretizations (4.2)–(4.5) is also set to be 10−9 . We will study the performance of the Newton-type sweeping algorithm by carrying out the comparison with the point-wise sweeping algorithm in terms of numbers of sweepings and the overall computational time. For the details of the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm, we refer readers to [13] for the first-order scheme and [14] for the high-order scheme. 5.1. 2-D examples 5.1.1. Example 1 In the first example we use a smooth model to illustrate the convergence property of the proposed algorithm. The exact solution is set as T (x, y) = sin(x) sin(y), which is smooth everywhere. The speed function F (x) is computed analytically from Eq. (1.1) by the following formula: F (x, y) =

|1 − v(x, y) · ∇ T (x, y)| , ∥∇ T (x, y)∥

(5.1)

where the velocity of the moving fluid is set to be v = (0.1, 0.2). The computational domain is [0, 1] × [0, 1], and the boundary condition T = 0 is imposed on two sides of the domain along x = 0 and y = 0. To study the order of convergence of the Newton-type sweeping algorithm, we carry out computations on a series of meshes with mesh size h = 1/40, 1/80, 1/160, 1/320, 1/640 and 1/1280, respectively. Since the exact solution is given in this synthetic model, we can evaluate errors of the numerical solutions exactly. Fig. 1 shows the plot of convergence for both first-order and high-order schemes. The errors are measured at grid points using L1 -, L2 -, and L∞ norms, respectively, and we have excluded the computational boundary when evaluating the errors. It shows that the algorithm enjoys good convergence for this example. The solutions are approximately first-order accurate for the firstorder scheme and approximately third-order accurate for the high-order scheme with third-order WENO approximations. The solutions on the 1280 × 1280 mesh are shown in Fig. 2, where Figs. 2 (b) and 2 (d) show the corresponding errors. To illustrate the performance of the algorithm, we record the numbers of sweepings and the computational time in computations, which are shown in Table 1. In Table 2, we provide the similar records of the pointwise Gauss–Seidel Lax– Friedrichs sweeping algorithm proposed in [13,14]. Both the Newton-type sweeping and the pointwise sweeping have the same initialization setup and convergence criterion so that we can compare their performances on equal footing. Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

9

Fig. 2. (Example 1) Solutions on 1280 × 1280 mesh. (a) First-order solution. (b) Error of first-order solution. (c) High-order solution with third-order WENO approximations. (d) Error of high-order solution.

Table 1 (Example 1) Performance of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. The computational time is in seconds. Mesh Sweeping Sweeping Sweeping Sweeping

number (first-order scheme) time (first-order scheme) number (high-order scheme) time (high-order scheme)

1/40

1/80

1/160

1/320

1/640

1/1280

99 0.5631 55 0.6339

129 1.4001 95 1.3243

183 4.5878 167 5.4158

281 17.004 299 23.887

455 71.964 551 119.58

777 371.58 1049 694.50

It shows that the Newton-type algorithm takes a smaller number of sweepings to converge than the pointwise sweeping algorithm. In all the mesh sizes we have used and in both the first-order and high-order scheme, the sweeping number of the Newton-type algorithm is about half less than that of the pointwise sweeping. It demonstrates that the semi-implicit feature of combining Newton’s method with the Gauss–Seidel iterations does improve upon the convergence property of the eikonal solver. On the other hand, the performance of the Newton-type algorithm in computational time is more sophisticated. For the first-order scheme, the Newton-type algorithm outperforms the traditional pointwise sweeping only on a very fine mesh such as h = 1/1280, and for the high-order scheme the new algorithm prevails at h ≤ 1/160 or smaller. The performance comparison suggests that the Newton-type algorithm can incur some overhead in setting up the overall computation, but such overhead eventually pays for large-scale computation in terms of very fine meshes and/or high-order schemes. Specifically, we can observe that on relatively coarser meshes, the Newton-type algorithm may spend more time in computation although it takes fewer sweepings to converge; but on finer meshes or for high-order schemes, the Newton-type algorithms outperform the pointwise sweeping algorithms, and the resulting efficiency and advantages are clear and significant. This phenomenon is consistent with the feature of the algorithm, as we have discussed in Section 3.3. The effect of vectorization in time-saving becomes substantial when the computational cost of each updating becomes large. We will see this phenomenon repeatedly in the following examples, and in 3-D computations the Newtontype algorithm simply outperforms the pointwise sweeping since the computational cost of each updating is sufficiently high. Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

10

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Table 2 (Example 1) Performance of the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm. The computational time is in seconds. The initialization and convergence criterion are the same as those of Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm, and so one can compare their performances. Mesh Sweeping Sweeping Sweeping Sweeping

number (first-order scheme) time (first-order scheme) number (high-order scheme) time (high-order scheme)

1/40

1/80

1/160

1/320

1/640

1/1280

154 0.3162 89 0.1626

218 0.2748 155 0.8182

333 1.5264 291 6.0205

541 9.5402 563 46.244

923 60.953 1095 358.65

1631 413.89 2115 2788.2

Table 3 (Example 2) Performance of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. The computational time is in seconds. Mesh Sweeping Sweeping Sweeping Sweeping

number (first-order scheme) time (first-order scheme) number (high-order scheme) time (high-order scheme)

2/40

2/80

2/160

2/320

2/640

2/1280

88 0.5370 68 0.7594

102 1.1670 89 1.3460

130 3.2933 136 4.4953

181 11.071 231 18.933

277 45.330 399 88.965

457 224.02 590 401.06

Fig. 3. (Example 2) Convergence of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping method. Both the x-axis and y-axis are in the logarithmic scale. (a) First-order scheme. (b) High-order scheme with third-order WENO approximations.

5.1.2. Example 2 In this example we study the convergence of the proposed algorithm in a model with a point-source boundary condition. The exact solution of Eq. (1.1) is taken as T (x, y) = (1 + x2 )(1 + y2 ). The speed function F (x) is again computed by Eq. (5.1) with v = (0, 0). The computational domain is [−1, 1] × [−1, 1], and the boundary condition T (xs , ys ) = 1 is imposed at the point source (xs , ys ) = (0, 0). Fig. 3(a) shows the convergence property of the first-order scheme, and Fig. 3(b) shows the convergence property of the high-order scheme with third-order WENO approximations. In Fig. 4, we provide the numerical solutions and the corresponding errors on the 1280 × 1280 mesh. In Table 3, we provide sweeping numbers and computational times of the Newton-type sweeping algorithm; this table can be compared with Table 4, which shows the corresponding results of the pointwise sweeping. We conclude that the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm enjoys good convergence behavior in the example with a point-source boundary condition. The first-order scheme achieves a first-order of convergence, and the high-order scheme with third-order WENO approximations has an order of convergence greater than two. In comparison to the traditional pointwise sweeping, the Newton-type algorithm takes much fewer sweepings to converge in all the mesh sizes for both first-order and high-order schemes. Similar to the situation of Example 1, the efficiency of the Newton-type algorithm in computational time becomes significant in the application of high-order schemes and in computations on fine meshes. 5.1.3. Example 3 We consider an example of two incident planar wavefronts from [18], where the wave speed F is given by F (x, y) =

{

2− 2

1 2

cos2 π y −

( (

1 2

))

if 0 ≤ y ≤ 1 otherwise,

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Fig. 4. (Example 2) Solutions on 1280 × 1280 mesh. (a) First-order solution. (b) Error of first-order solution. (c) High-order solution with third-order WENO approximations. (d) Error of high-order solution.

Table 4 (Example 2) Performance of the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm. The computational time is in seconds. Mesh Sweeping Sweeping Sweeping Sweeping

number (first-order scheme) time (first-order scheme) number (high-order scheme) time (high-order scheme)

2/40

2/80

2/160

2/320

2/640

2/1280

123 0.3584 91 0.1782

154 0.2208 135 0.7567

211 1.1049 236 5.0059

319 6.2269 438 36.743

517 37.899 778 263.00

884 246.09 1235 1667.4

Table 5 (Example 3) Performances of the Newton-type and pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithms. The computational time is in seconds. Newton-type pointwise

Sweeping number (first-order)

Sweeping time (first-order)

Sweeping number (high-order)

Sweeping time (high-order)

1584 1753

611.17 337.68

636 1357

339.53 1303.8

and the ambient velocity is set to be v = (0, 0.8). The boundary condition T = 0 is imposed along x = 0 and x = 1 such that there are two planar wavefronts normal to the x-axis entering into the computational domain. We take the computational domain to be [0, 1] × [−0.1, 1.1] and the mesh size h = 0.001. For the implementation of high-order schemes with third-order WENO discretizations, one needs exact or interpolated values of the solution T (x) at grid points whose distance to the boundary is less than or equal to 2h. In this example, however, the exact solution is unavailable, and we simply apply the first-order solution to those grid points. Fig. 5 shows the numerical solutions by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm, where Fig. 5 (a) shows the first-order solution and Fig. 5 (b) shows the high-order solution. The two wavefronts entering into the computational domain merge to form a single wavefront, and the wavefronts are symmetric with respect to x = 0.5 but not with respect to y = 0.5 due to the ambient velocity v = (0, 0.8). In Table 5, we show the performance of the algorithm, and compare it with that of the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm. Again, the Newton-type sweeping algorithm needs a smaller number of sweepings to achieve convergence, and it outperforms the pointwise sweeping algorithm in the application of high-order schemes. Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

12

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 5. (Example 3) Solutions of two incident planar waves with mesh size h = 0.001. (a) First-order solution. (b) High-order solution with third-order WENO approximations.

Fig. 6. (Example 4) Marmousi model. Distribution of wave speed. Table 6 (Example 4 (1): v = (0, 0)) Performances of the sweeping algorithms. The computational time is in seconds. Newton-type pointwise

Sweeping number (first-order)

Sweeping time (first-order)

Sweeping number (high-order)

Sweeping time (high-order)

253 421

6.1721 4.4162

418 1163

13.926 43.540

Table 7 (Example 4 (2): v = (0.4, 0.2)) Performances of the sweeping algorithms. The computational time is in seconds. Newton-type pointwise

Sweeping number (first-order)

Sweeping time (first-order)

Sweeping number (high-order)

Sweeping time (high-order)

392 720

9.5159 6.3898

333 608

10.926 22.384

5.1.4. Example 4 In this example, we use a Marmousi velocity model to demonstrate the robustness of the Newton-type sweeping algorithm. The Marmousi model was built to resemble an overall continental drift geological setting, and it contains plenty of complex velocity settings. Fig. 6 shows the profile of the wave speed function F (x). The computational domain is [0, 9.192] km ×[0, 2.904] km, and the boundary condition T (xs , ys ) = 0 is imposed at the point source (xs , ys ) = (4, 0). We have considered three configurations of the external velocity: (1) v = (0, 0); (2) v = (0.4, 0.2); and (3) v = (0.4 sin(x), 0.2 cos(x)). The Newton-type sweeping algorithm is used to solve the generalized eikonal equation (1.1) with the mesh size h = 0.024 km. The solutions of traveltimes T are shown in Figs. 7, 8 and 9, respectively, where both first-order and high-order solutions are displayed. The contours of the traveltime solutions depict the shape of wavefronts in the wave propagation, which appear irregular in this example due to the complex settings of wave speeds and ambient velocities. We mention that the solutions by the Newton-type sweeping algorithm are consistent with those obtained by the traditional pointwise sweeping algorithm. It demonstrates that the newly designed Newton-type algorithm is effective in solving the generalized eikonal equation. In Tables 6, 7 and 8, respectively, we compare the Newton-type algorithm with the original pointwise sweeping method. Similar to the situation of the preceding examples, the Newtontype algorithm needs a smaller number of sweepings to converge in both first-order and high-order schemes, and in terms of computational time the efficiency of the Newton-type algorithm becomes evident in the application of high-order schemes. Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Fig. 7. (Example 4 (1): v = (0, 0)) Solutions of traveltimes by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. (a) First-order solution. (b) High-order solution with third-order WENO approximations.

Fig. 8. (Example 4 (2): v = (0.4, 0.2)) Solutions of traveltimes by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. (a) First-order solution. (b) High-order solution with third-order WENO approximations.

Fig. 9. (Example 4 (3): v = (0.4 sin(x), 0.2 cos(x))) Solutions of traveltimes by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. (a) First-order solution. (b) High-order solution with third-order WENO approximations. Table 8 (Example 4 (3): v = (0.4 sin(x), 0.2 cos(x))) Performances of the sweeping algorithms. The computational time is in seconds. Newton-type pointwise

Sweeping number (first-order)

Sweeping time (first-order)

Sweeping number (high-order)

Sweeping time (high-order)

355 647

9.1302 6.2926

290 614

10.218 22.620

5.2. 3-D examples The proposed Newton-solver based linewise sweeping algorithm can be directly extended to 3-D computations. In the computational domain [xmin , xmax ] × [ymin , ymax ] × [zmin , zmax ] with mesh size h = ∆x = ∆y = ∆z, the grid point xi,j,k = (xi , yj , zk ) has the spatial coordinate xi = xmin + (i − 1)h, yj = ymin + (j − 1)h and zk = zmin + (k − 1)h. At each update, the column-wise Newton solver updates the jth column of T (x) on the kth page corresponding to the coordinates {(x, y, z) | y = yj , z = zk }. The global solver sweeps the 3-D computational domain in four directions: (1) k = 1 : N , j = 1 : N ;

(2) k = 1 : N , j = N : 1;

(3) k = N : 1, j = 1 : N ;

(4) k = N : 1, j = N : 1.

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

14

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Table 9 (Example 5) Performance of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. The computational time is in seconds. Mesh Sweeping Sweeping Sweeping Sweeping

number (first-order scheme) time (first-order scheme) number (high-order scheme) time (high-order scheme)

1/20

1/40

1/80

1/160

1/320

166 8.6515 58 3.8719

217 47.438 106 27.954

299 284.32 202 237.38

441 1973.6 385 2158.0

697 14815 739 20919

Table 10 (Example 5) Performance of the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm, where the computational time is in seconds. We did not do computation in the mesh size h = 1/320 because it can be too costly. Mesh Sweeping Sweeping Sweeping Sweeping

number (first-order scheme) time (first-order scheme) number (high-order scheme) time (high-order scheme)

1/20

1/40

1/80

1/160

1/320

201 13.681 71 6.2981

274 126.17 138 102.28

390 1245.5 277 1750.9

589 13176 547 28093

∗ ∗ ∗ ∗

Table 11 (Example 6) Performance of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm in 3-D computation. The computational time is in seconds. Sweeping number (first-order)

Sweeping time (first-order)

Sweeping number (high-order)

Sweeping time (high-order)

3084

24518

1556

15722

Fig. 10. (Example 5) Convergence of the Newton-type Gauss–Seidel Lax–Friedrichs sweeping method. Both the x-axis and y-axis are in the logarithmic scale. (a) First-order scheme. (b) High-order scheme with third-order WENO approximations.

The setups of initialization and convergence criterion are the same as those of 2-D computations. In the following, we use numerical examples to demonstrate the Newton-type Gauss–Seidel Lax–Friedrichs algorithm in 3-D spaces. We will see that the advantages of the Newton-type algorithm are more significant in 3-D applications. 5.2.1. Example 5 To illustrate the convergence property, we consider a synthetic smooth model which is similar to the 2-D model in Example 1. The exact solution of Eq. (1.1) is given by T (x, y, z) = sin(x) sin(y) sin(z). The external velocity is set to be v = (0.1, 0.2, 0.3), and the speed function F is analytically computed by F (x, y, z) =

|1 − v · ∇ T (x, y, z)| . ∥∇ T (x, y, z)∥

The computational domain is [0, 1]3 , where the boundary condition T = 0 is imposed along three sides: x = 0, y = 0 and z = 0. We have performed computations in a set of meshes with mesh size h = 1/20, 1/40, 1/80, 1/160 and 1/320, respectively. Fig. 10 shows the plot of convergence, where the errors are measured at the grid points using the L1 -, L2 -, and L∞ - norm, respectively, and we have excluded the computational boundary when evaluating the errors. It shows that both the first-order and high-order schemes have achieved their expected accuracies in 3-D computations. In Fig. 11 and Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

15

Fig. 11. (Example 5) First-order solution by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm; the mesh size h = 1/320. (a) Isosurfaces of the solution from T = 0.05 to T = 0.55 with an increment of 0.05. (b) Cross-section of the solution along x = 0.4. (c) Cross-section along y = 0.5. (d) Cross-section along z = 0.6.

Fig. 12, we provide the numerical solutions on the 320 × 320 × 320 mesh. To visualize the 3-D solutions, we plot the isosurfaces of the traveltimes and 2-D cross-sections of the solutions. To study the performance of the Newton-type algorithm in the 3-D computation, we record the numbers of sweepings and computational times, as shown in Table 9, and compare them with those of the pointwise Gauss–Seidel Lax–Friedrichs sweeping algorithm, which are shown in Table 10. It shows that the Newton-type algorithm outperforms the pointwise sweeping algorithm in both sweeping numbers and computational times. In all the mesh sizes we have used and in both first-order and high-order schemes, the Newton-type algorithm takes fewer sweepings to converge and uses less time to finish computations. At the mesh size h = 1/320, the expected computational time of the pointwise sweeping algorithm will go beyond three days, and we choose not to carry out such computations. It demonstrates that the advantages of the Newton-type algorithm become significant in 3-D applications. As the numerical complexity of each updating iteration increases substantially, the vectorized manner of the Newton-type algorithm in manipulating the traveltime solution takes full advantage of array programming to save significant computational times. 5.2.2. Example 6 In this example we consider a Gaussian wave speed, F (x, y, z) = 2.7 − 0.5y − 0.6 exp{−10(x − 0.3)2 − 15(y − 0.4)2 − 20(z − 0.3)2 }, and a constant ambient velocity v = (0.4, 0.8, 1.2); the source condition T = 0 is imposed at the point source: (xs , ys , zs ) = (0.5, 0.5, 0.5). We perform computations in Ω = [0, 1]3 with mesh size h = 1/200. Fig. 13 shows the first-order solution by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm, where Fig. 13 (a) plots isosurfaces of the traveltime, and Fig. 13 (b)–(d) display cross-sections of the solution. Correspondingly, the high-order solution with third-order WENO approximations are shown in Fig. 14. It is clear that the propagation of wavefronts speeds up in the top-right direction due to the ambient velocity v = (0.4, 0.8, 1.2). To illustrate the performance of the Newton-type algorithm, we show sweeping numbers and computational times in Table 11. We conclude that the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm is successfully and efficiently implemented in 3-D computations. 6. Conclusion We have developed a Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm to solve the generalized eikonal equation arising from isotropic wave propagation in a moving fluid. The Lax–Friedrichs scheme is used for discretization, Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

16

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 12. (Example 5) High-order solution by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm; third-order WENO approximations are implemented, and the mesh size h = 1/320. (a) Isosurfaces of the solution from T = 0.05 to T = 0.55 with an increment of 0.05. (b) Cross-section of the solution along x = 0.4. (c) Cross-section along y = 0.5. (d) Cross-section along z = 0.6.

Fig. 13. (Example 6) First-order solution by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm. (a) Isosurfaces of the solution from T = 0.1 to T = 1.1 with an increment of 0.1. (b) Cross-section of the solution along x = 0.5. (c) Cross-section along y = 0.5. (d) Cross-section along z = 0.5.

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

17

Fig. 14. (Example 6) High-order solution by the Newton-type Gauss–Seidel Lax–Friedrichs sweeping algorithm; third-order WENO approximations are implemented. (a) Isosurfaces of the solution from T = 0.1 to T = 1.1 with an increment of 0.1. (b) Cross-section of the solution along x = 0.5. (c) Cross-section along y = 0.5. (d) Cross-section along z = 0.5.

and different from the traditional pointwise Lax–Friedrichs sweeping algorithm, the Newton-type algorithm solves the discretized system in a line-wise manner. The values of T on an entire line are updated simultaneously by Newton’s method, and the Gauss–Seidel iterations with line-wise sweepings are proposed to achieve a global solution. We have developed the first-order scheme and a high-order scheme by applying third-order WENO approximations to derivatives. Extensive numerical examples have been presented to illustrate the performance of the algorithm in both 2-D and 3-D computations. Numerical examples demonstrate that the semi-implicit feature of combining Newton’s method with Gauss–Seidel iterations does improve upon the convergence speed of the original sweeping algorithms, and the Newton-type sweeping algorithm takes much fewer sweepings to converge in comparison to the pointwise sweeping algorithm. In addition, the new algorithm significantly outperforms the pointwise sweeping algorithm in computational times when the computational cost of updating T becomes high such as in 3-D computations, in high-order schemes, and in computations with fine meshes. These examples also validate that the vectorized manner of the newly designed algorithm is naturally amenable to array programming so as to save computational times. To conclude, the Newton-type Gauss–Seidel Lax– Friedrichs sweeping algorithm turns out to be a robust and efficient method for the solution of generalized eikonal equations at large-scale discretization. Acknowledgments The authors thank Dr. Wangtao Lu for discussions related to this project. Wenbin Li is supported by NSFC grant 41804096 and Natural Science Foundation of Guangdong Province grant 2018A030313341. Qian’s research is partially supported by NSF. References [1] E.T. Kornhauser, Ray theory for moving fluids, J. Acous. Soc. Amer. 25 (1953) 945–949. [2] N. Tanushev, J. Qian, J. Ralston, Mountain waves and Gaussian beams, SIAM J. Multiscale Model. Simul. 6 (2007) 688–709. [3] Ø. Hjelle, S. Petersen, A Hamilton–Jacobi framework for modeling folds in structural geology, Math. Geosci. 43 (7) (2011) 741–761.

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.

18

W. Li and J. Qian / Computers and Mathematics with Applications xxx (xxxx) xxx

[4] R. Blumrich, F. Coulouvrat, D. Heimann, Meteorologically induced variability of sonic-boom characteristics of supersonic aircraft in cruising flight, J. Acous. Soc. Amer. 118 (2005) 707–722. [5] A. Loubeau, F. Coulouvrat, Effects of meteorological variability on sonic boom propagation from hypersonic aircraft, AIAA J. 47 (11) (2009) 2632–2641. [6] P. Bois, M. La Porte, M. Lavergne, G. Thomas, Well-to-well seismic measurements, Geophysics 37 (3) (1972) 471–480. [7] V. Cerveny, Seismic Ray Theory, Cambridge University Press, Cambridge, 2001. [8] S. Osher, A level set formulation for the solution of the Dirichlet problem for Hamilton-Jacobi equations, SIAM J. Math. Anal. 24 (1993) 1145–1152. [9] J. Tsitsiklis, Efficient algorithms for globally optimal trajectories, IEEE Tran. Autom. Control 40 (1995) 1528–1538. [10] J. Sethian, A fast marching level set method for monotonically advancing fronts, Proc. Natl. Acad. Sci. 93 (4) (1996) 1591–1595. [11] H.K. Zhao, Fast sweeping method for eikonal equations, Math. Comp. 74 (2005) 603–627. [12] R. Tsai, L.-T. Cheng, S.J. Osher, H.K. Zhao, Fast sweeping method for a class of Hamilton-Jacobi equations, SIAM J. Numer. Anal. 41 (2003) 673–694. [13] C.Y. Kao, S.J. Osher, J. Qian, Lax-Friedrichs sweeping schemes for static Hamilton-Jacobi equations, J. Comput. Phys. 196 (2004) 367–391. [14] Y.T. Zhang, H.K. Zhao, J. Qian, High order fast sweeping methods for static Hamilton-Jacobi equations, J. Sci. Comput. 29 (2006) 25–56. [15] J. Qian, Y.-T. Zhang, H.-K. Zhao, A fast sweeping method for static convex Hamilton–Jacobi equations, J. Sci. Comput. 31 (1) (2007) 237–271. [16] C.Y. Kao, S. Osher, J. Qian, Legendre Transform based fast sweeping methods for static Hamilton-Jacobi equations on triangulated meshes, J. Comput. Phys. 227 (2008) 10209–10225. [17] S. Serna, J. Qian, A stopping criterion for higher-order sweeping schemes for static Hamilton-Jacobi equations, J. Comput. Math. 28 (2010) 552–568. [18] D. Dahiya, S. Baskar, F. Coulouvrat, Characteristic fast marching method for monotonically propagating fronts in a moving medium, SIAM J. Sci. Comput. 35 (2013) A1880–A1902. [19] H. Zhao, S. Osher, B. Merriman, M. Kang, Implicit and non-parametric shape reconstruction from unorganized points using variational level set method, Comput. Vis. Image Underst. 80 (2000) 295–319. [20] G.S. Jiang, D. Peng, Weighted ENO schemes for Hamilton-Jacobi equations, SIAM J. Sci. Comput. 21 (2000) 2126–2143. [21] R. Glowinski, S. Leung, J. Qian, Operator-splitting based fast sweeping methods for isotropic wave propagation in a moving fluid, SIAM J. Sci. Comput. 38 (2016) A1195–A1223. [22] S. Bak, J.R. McLaughlin, D. Renzi, Some improvements for the fast sweeping method, SIAM J. Sci. Comput. 32 (2010) 2853–2874. [23] B. Datta, Numerical linear algebra and applications, SIAM, 2010. [24] S.J. Osher, C.W. Shu, High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations, SIAM J. Numer. Anal. 28 (1991) 907–922.

Please cite this article as: W. Li and J. Qian, Newton-type Gauss–Seidel Lax–Friedrichs high-order fast sweeping methods for solving generalized eikonal equations at large-scale discretization, Computers and Mathematics with Applications (2019), https://doi.org/10.1016/j.camwa.2019.08.031.