A Fixed‑Point Fast Sweeping WENO Method with Inverse Lax‑Wendroff Boundary Treatment for Steady State of Hyperbolic Conservation Laws

Fixed-point fast sweeping WENO methods are a class of efficient high-order numerical methods to solve steady-state solutions of hyperbolic partial differential equations (PDEs). The Gauss-Seidel iterations and alternating sweeping strategy are used to cover characteristics of hyperbolic PDEs in each sweeping order to achieve fast convergence rate to steady-state solutions. A nice property of fixed-point fast sweeping WENO methods which distinguishes them from other fast sweeping methods is that they are explicit and do not require inverse operation of nonlinear local systems. Hence, they are easy to be applied to a general hyperbolic system. To deal with the difficulties associated with numerical boundary treatment when high-order finite difference methods on a Cartesian mesh are used to solve hyperbolic PDEs on complex domains, inverse Lax-Wendroff (ILW) procedures were developed as a very effective approach in the literature. In this paper, we combine a fifth-order fixed-point fast sweeping WENO method with an ILW procedure to solve steady-state solution of hyperbolic conservation laws on complex computing regions. Numerical experiments are performed to test the method in solving various problems including the cases with the physical boundary not aligned with the grids. Numerical results show high-order accuracy and good performance of the method. Furthermore, the method is compared with the popular third-order total variation diminishing Runge-Kutta (TVD-RK3) time-marching method for steady-state computations. Numerical examples show that for most of examples, the fixed-point fast sweeping method saves more than half CPU time costs than TVD-RK3 to converge to steady-state solutions


Introduction
Steady-state problems of hyperbolic conservation laws are important mathematical models arising in many applications of compressible fluid dynamics. High-order accuracy WENO schemes are a class of popular numerical methods for spatial discretizations of nonlinear hyperbolic conservation laws, e.g., see [1,7,15,[17][18][19]. After spatial discretizations of a steady-state problem by a high-order WENO scheme, a large nonlinear algebraic system is formed. One of the important factors which determines efficiency of the numerical simulations is how to design fast algorithms for solving the large nonlinear system. For both linear and nonlinear stability of high-order schemes, a common approach is to use the third-order total variation diminishing Runge-Kutta (TVD-RK3) scheme [5,19] to timemarch numerical solution to steady state. To improve computational efficiency, explicit fast sweeping WENO methods (e.g., [24,25,29,30]) were developed to accelerate the convergence to steady-state solutions of high-order schemes. Fast sweeping methods are a class of efficient iterative methods developed to solve steady-state problems of hyperbolic PDEs. The methods use alternating sweeping strategy to cover a family of characteristics in a certain direction simultaneously in each sweeping order. Coupled with the Gauss-Seidel iterations, these methods can achieve a fast convergence speed. Implicit fast sweeping methods are very efficient and often have linear computational complexity (see, e.g., [4,13,14,31]); however, the algorithms are much more complicated than explicit methods when they are coupled with high-order schemes [9,23,28]. To solve steady-state problems of the nonlinear hyperbolic conservation laws by high-order WENO fast sweeping schemes, the fixed-point fast sweeping methods were applied [3,24] which were based on the fixed-point fast sweeping methods for static Hamilton-Jacobi equations [29]. This kind of fast sweeping methods have nice properties such as that they have explicit forms and do not involve inverse operation of nonlinear local systems, and they can be adopted to solve general hyperbolic systems with any monotone numerical fluxes and high-order WENO approximations easily. Numerical experiments show that the methods converge to steady-state solutions much faster than regular time-marching approach. Since it shows that the WENO schemes with unequal-sized substencils [32,33] have very nice convergence property in solving steady-state problems, e.g., the high-order multi-resolution WENO schemes [35], the fixed-point fast sweeping methods were combined with the multi-resolution WENO scheme in [10] to overcome the issue that fast sweeping iteration residue does not fully converge to values of round-off error level for difficult steady-state problems, and a fully convergent fast sweeping WENO method (called absolutely convergent fixed-point fast sweeping WENO method) was obtained.
When the computational domain has complex geometries, the efficient finite difference method is difficult to be applied, because the Cartesian grid points are usually not located on the physical boundary and the boundary may intersect the grids in an arbitrary fashion. There are many excellent approaches to deal with this. For example, an accurate implementation of solid wall boundary conditions in curved geometries was developed in [8], and a Cartesian embedded boundary method based on the finite difference formulation was designed in [20]. In [6,25], the inverse Lax-Wendroff (ILW) numerical boundary condition procedure was developed for solving static Hamilton-Jacobi equations. The effective ILW methods were successfully developed to solve hyperbolic conservation laws on complex geometry domains in [11,21,22].
In this paper, we combine the absolutely convergent fixed-point fast sweeping WENO method [10] with the ILW procedure to solve steady-state solution of hyperbolic 1 3 conservation laws on complex domains. Numerical experiments are performed to test the method in solving various problems including the cases with the physical boundary not aligned with the grids. It is shown that the absolutely convergent fixed-point fast sweeping method with the ILW procedure is more efficient than other popular approach such as the TVD-RK3 time-marching method to converge to steady-state solutions of the fifth-order multi-resolution WENO scheme. Furthermore, it can drive the residue of iterations to converge to machine zero/round off errors for all benchmark problems tested. The rest of the paper is organized as follows. The detailed numerical algorithm is described in Sect. 2. In Sect. 3, we provide extensive numerical experiments to test and study the proposed method, and to perform comparisons of different methods. Concluding remarks are given in Sect. 4.

Numerical Methods
We consider steady-state problems of hyperbolic conservation laws on a bounded domain with appropriate boundary conditions prescribed on . Here, u is the vector of the unknown conservative variables, f(u) and g(u) are the vectors of flux functions, and R is the source term. In the following, we first describe the absolutely convergent fixed-point fast sweeping WENO method [10] to solve (1).

A Fifth-Order Absolutely Convergent Fixed-Point Fast Sweeping WENO Method
A spatial discretization of (1) by a high-order WENO scheme leads to a large nonlinear system of algebraic equations with the size determined by the number of spatial grid points. We partition by a uniform Cartesian mesh {(x i , y j )} with the grid size Δx in the x-direction and Δy in the y-direction. Since high-order WENO schemes with unequal-sized substencils [34,35] show nice convergence property in solving steady-state problems, a fifth-order multi-resolution WENO scheme [33] is used in [10] for the spatial discretization. The procedure is outlined here. For upwinding and linear stability of the scheme, we perform a flux splitting. For example, the , which satisfies the condition that df + (u) du ⩾ 0 and df − (u) du ⩽ 0 for the scalar case, or the corresponding eigenvalues of the Jacobian matrix are positive/negative for the system case with a local characteristic decomposition. Here, a simple Lax-Friedrichs flux splitting f ± (u) = 1 2 (f (u) ± u) is used, where = max u |f � (u)| . Then, each of them is approximated separately by the fifth-order multi-resolution WENO scheme to obtain numerical fluxes f + i+1∕2,j and f − i+1∕2,j at the point (x i+1∕2 , y j ) , respectively, where The nonlinear weights are computed as

3
where is taken as = 10 −6 to avoid the case that the denominator becomes zero. The details on the linear weights l 1 ,3 , computations of smoothness indicators l 1 , the quantity , and reconstruction polynomials p l 1 (x) can be found in [10,33]. The reconstruction of the numerical flux f − i+1∕2,j is mirror-symmetric with respect to x i+1∕2 , and the similar procedure is applied to the y-direction flux g (u).
After the fifth-order multi-resolution WENO finite difference spatial discretization, we obtain a nonlinear algebraic system where N and M are the numbers of grid points in the x and y directions, respectively. f i+1∕2,j and ĝ i,j+1∕2 are the numerical fluxes. The right-hand-side (RHS) of (3) is a nonlinear function of the numerical values at the grid points of the computational stencils. We denote it by L. To solve (3), the forward Euler type fixed-point fast sweeping method [10,24,29] is used. It has the following form: u n ij is the nth step iteration value at the grid point (x i , y j ) . is the Courant-Friedrichs-Lewy (CFL) number and x ∕Δx+ y ∕Δy is corresponding to the time step size Δt n . x = max u |f � (u)| and y = max u |g � (u)| represent the maximum characteristic speeds in each spatial direction. r, s are values which depend on the stencil of the WENO reconstruction. For the fifthorder multi-resolution WENO scheme used here, we have r = s = 3 . In the fast sweeping method, the Gauss-Seidel iterations and alternating direction sweepings are used. The Gauss-Seidel philosophy requires that the newest numerical values of u are used in the fifth-order multi-resolution WENO reconstruction stencils whenever they are available. The iterations do not just proceed in only one direction i = 1 ∶ N, j = 1 ∶ M , but in the following four alternating directions repeatedly: This is marked as " i = i 1 , ⋯ , i N ; j = j 1 , ⋯ , j M " in the scheme (4) to show the alternating iteration directions. The purpose of using alternating direction sweepings is to utilize the characteristics property of hyperbolic PDEs. By combining with the Gauss-Seidel philosophy, a fast convergence to the steady-state solution is achieved, as shown in the following numerical experiments. In the Gauss-Seidel iterations, the newest numerical values on the stencil of the WENO scheme are used if they are available. This fact is shown by the notation u * in the scheme (4) to represent the numerical values in the WENO stencil, with the understanding that u * k,l could be u n k,l or u n+1 k,l , depending on the current sweeping direction.

Remark 1
The similar fast sweeping techniques can also be applied to a TVD-RK scheme as that in [24]. However, it was found in [24] that the scheme resulted by applying these fast sweeping techniques to the TVD-RK type scheme is less efficient than the forward Euler type fixed-point fast sweeping method. Hence, in this paper, we only focus on the forward Euler type fixed-point fast sweeping method coupled with the fifth-order multiresolution WENO scheme.

The ILW Boundary Treatment
In this section, we outline the ILW numerical boundary treatment in both one-dimensional (1D) and two-dimensional (2D) problems. The ILW methods have been developed and studied extensively in, e.g., [21,22], so here we just give a very brief description of the ILW methods, and the details are omitted and referred to [21,22]. The only difference from [21,22] is that in this paper, the fifth-order multi-resolution WENO method is used for the extrapolation step in the ILW methods, and we provide detailed description about that.

1D Case
First we consider 1D compressible steady-state Euler equations where the conservative variables U ≜ (U 1 , U 2 , U 3 ) T = ( , u, E) T , and the flux F(U) = . , u, p, and E denote the density, velocity, pressure, and total energy, respectively. R(U, x, y) = (R 1 (U, x, y), R 2 (U, x, y), R 3 (U, x, y)) T are the source terms. The equation of state has the form E = p � −1 + 1 2 u 2 , where ′ is the specific heat ratio and � = 1.4 for the case of air. The sound speed is c = √ � p∕ . Without loss of generality, the computational domain is taken as x ∈ (−1, 1).
The computational domain is partitioned by a uniform mesh where h is the size of the grid. We proceed as in [21] to compute numerical values on ghost points. Here, the computation of U(x N+1 ) is taken as an example. Other values such as U(x −1 ), U(x −2 ), ⋯ , and U(x N+2 ), U(x N+3 ), ⋯ are treated in a similar way. By a fifth-order Taylor approximation where U * (k) m is a (5 − k)th-order approximation of the spatial derivative k U m x k | x=1 . Then, a local characteristic decomposition is performed to take care of the inflow and outflow boundary conditions. At the boundary, the Jacobian matrix of the flux is denoted by , and A ⟂ (U(x N )) has three eigenvalues and a complete set of left eigenvectors l 1 (U(x N )), l 2 (U(x N )), l 3 (U(x N )) . The local characteristic variables V m at interior grid points near the boundary are defined by Using the numerical values on these interior points, we can find fourth-order polynomials to compute extrapolation values for both functions and their derivatives at the bound- For a more robust algorithm to deal with possible discontinuities near the boundary, we use a multi-resolution WENO type extrapolation which is described in details below. Then, following the procedure in [21], these extrapolation values V * (k) m are coupled with the inflow boundary conditions and the PDEs itself to compute these values U * (k) m which are needed in (6). The procedure is summarized as follows. (i) Compute the eigenvalues m and left eigenvectors l m of the Jacobian matrix A ⟂ (U(x N )) described above, for m = 1, 2, 3 . Determine the inflow boundary conditions by the signs of the eigenvalues m . The inflow boundary conditions are prescribed and given in the problem. (ii) Compute the outgoing local characteristic variables

using the prescribed inflow boundary conditions and the extrapolated values
, use the ILW procedure of the PDEs to write the kth-order derivatives of the prescribed inflow boundary conditions as a linear combination of kth-order spatial derivatives and other terms involving lower order spatial derivatives, and together with the extrapolated values, we form a linear system for U * (k) m , m = 1, 2, 3 . Solve the linear system to obtain U * (k) m , m = 1, 2, 3 . (v) Use the Taylor approximation (6) to find the values of the ghost points.

1D Multi-resolution WENO Extrapolation
Here, we describe the multi-resolution WENO extrapolations for the 1D case, i.e., construct a (5 − k)th-order approximation V * (k) m at the boundary x = 1 , using the numerical values at the interior grid points of the domain V m (x j ), j = N − 4, ⋯ , N . For the simplicity of notation, we omit the subscript m and use V j to denote V m (x j ) for a specific m. Three substencils are used to construct WENO extrapolations, i.e., the substencils The extrapolation procedure is given in the following.
Step 3 Compute the smoothness indicators l which measure how smooth the functions p l (x) for l = 1, 2, 3 are in the interval [x N , x N+1 ] . We use the same recipe for the smoothness indicators as in [7,17], Notice here 1 = 0 which is a simplification of that in [33]. Numerical experiments in this paper show that high-order accuracy and resolution, and nonlinear stability of the computations are still maintained by taking this simplified value of 1 .
Step 4 Compute the nonlinear weights based on the linear weights and the smoothness indicators. The WENO-Z type nonlinear weights as in [1,2] are used. First, a quantity which depends on the absolute differences between the smoothness indicators is calculated . The nonlinear weights are Here, = 10 −6 is used to avoid the case that the denominator becomes zero. Then, the final (5 − k)th-order WENO approximation V * (k) at the boundary x = 1 is

2D Case
We consider 2D compressible steady-state Euler equations with appropriate boundary conditions, where U = (U 1 , T . , u, v, p, and E denote the density, x-velocity, y-velocity, pressure, and total energy, respectively. The equation of state has the form E = p The stencil S for 2D fifthorder multi-resolution WENO extrapolation We follow the procedure of the ILW boundary treatment in [21]. For example, a ghost point P = (x i , y j ) , we first identify the point P 0 = (x 0 , y 0 ) on the boundary = , such that the outward normal n(P 0 ) at P 0 goes through P. See Fig. 1 for an illustration. The coordinates of the point P 0 and the outward normal n(P 0 ) are obtained from the given geometry of . We set up a local coordinate system (x,ŷ) at P 0 by where is the angle between the outward normal n(P 0 ) and the x-axis and T is the corresponding rotation matrix. The x-axis has the same direction as n(P 0 ) and the ŷ-axis is the tangential direction. Under this local coordinate system, the Euler equations (8) For the fifth-order boundary treatment, the numerical value at the ghost point P is approximated by the Taylor expansion To decide the inflow or outflow boundary conditions to apply in the ILW procedure, a local characteristic decomposition is performed. We use the numerical values Û (0) of a grid point which is the closest to P 0 among all the grid points inside the domain , to form the Jacobian matrix of the normal flux by where S is the stencil for the fifth-order multi-resolution WENO extrapolation, which is used to compute V * (k) m , the (5 − k)th-order approximation of The fifth-order multiresolution WENO extrapolation will be described in the following. Then, following the procedure in [21,22], these extrapolation values V * (k) m are coupled with the inflow boundary conditions and the PDEs itself to compute these values Û * (k) m which are needed in the Taylor expansion (9).

2D Multi-resolution WENO Extrapolation
Here, we describe the 2D fifth-order multi-resolution WENO extrapolation. Following [11,22], interior grid points in the neighborhood of the point P 0 which are distributed around the normal n(P 0 ) ( x-axis) are selected as the stencil S of the WENO extrapolation. For the fifthorder multi-resolution WENO extrapolation used here, S contains three substencils S 1 , S 2 , S 3 , where S 3 is the same as S, and S 1 has just one point, while S 2 has nine points and S 3 is composed of 25 points to guarantee the fifth-order accuracy for the extrapolation function. See Fig. 1 for an illustration of the stencil, where S 1 contains the point 1, and S 2 contains the points 1-9, while S 3 includes the points 1-25. Given the numerical values of the variables V m on the stencil S, we construct a constant polynomial q 1 (x, y) based on S 1 , a degree 2 polynomial q 2 (x, y) based on S 2 , and a degree 4 polynomial q 3 (x, y) based on S 3 . Since there are more given values than the degrees of freedom, the least square method is used to find the (V m ) , = l m (Û (0) )Û , , m = 1, ⋯ , 4, (x , y ) ∈ S, polynomials q 2 and q 3 . Following the multi-resolution WENO reconstruction idea in [33], we obtain equivalent expressions for these polynomials, that is p 1 (x, y) = q 1 (x, y) , y) . The linear weights are taken as 1,2 = 1∕11 , 2,2 = 10∕11 , and 1,3 = 1∕111, 2,3 = 10∕111, 3,3 = 100∕111 . Then, the multi-resolution WENO extrapolation gives Here, the nonlinear weights are computed as follows: The quantity =

, and the smoothness indicators
Here, is the multi-index notation for partial derivative operator D , and = 10 −6 is used to avoid the case that the denominator becomes zero. As that in the 1D case, here, the smoothness indicator 1 = 0 which is a simplification of that in [33]. Numerical experiments in the following section verify that high-order accuracy and resolution, and nonlinear stability of the simulations are maintained well by taking this simplified value of 1 .

Numerical Experiments
In this section, we perform numerical experiments to test the performance of the fixedpoint fast sweeping multi-resolution WENO method with the ILW boundary treatment for steady-state computations. A popular approach to solve steady-state problems is to apply time-marching approach to the system (3). For example, the TVD-RK3 scheme in [5,19] has been broadly used for time-marching high-order WENO schemes to steady states, and both linear and nonlinear stability of high-order WENO schemes are maintained. To show the advantage of the fixed-point fast sweeping multi-resolution WENO method with the ILW boundary treatment, the computational efficiency of the TVD-RK3 scheme and this fast sweeping scheme is compared. Time-marching methods are actually Jacobi type fixedpoint iterations (see, e.g., [10,24,29]). For the presentation purpose, the fixed-point fast sweeping method is called "FE fast sweeping scheme", and the TVD-RK3 time-marching scheme "RK Jacobi scheme" in this section. The convergence of the iterations is measured by the average residue which is defined as where the local residue at the grid point i is N is the total number of grid points and n is the iteration step. Δt n = x ∕Δx+ y ∕Δy . For the 2D Euler system, the average residue is defined as . Here, R * i 's are the local residuals of the conservative variables, namely, The convergence criterion is set to be Res A ⩽ , where the threshold value is taken to be at round-off error level 10 −11 − 10 −13 . In the numerical examples, the CPU times which are costed by the TVD-RK3 scheme and the fast sweeping scheme to reach their numerical steady states are recorded and compared. Also, we record the final time which is the accumulation of time step sizes Δt n , i.e., the total time that numerical schemes are time-marched from the initial condition at t = 0 to the time step that the numerical steady states are reached. In this paper, the number of iterations reported in the numerical simulations counts a complete update of numerical values in all grid points once as one iteration. Hence, these four alternating directions in the FE fast sweeping scheme are counted as four iterations. In the numerical experiments, for various examples, we compare the computational efficiency of the RK Jacobi scheme and the FE fast sweeping scheme using the largest possible CFL numbers which lead to iteration convergence with the fastest speed for each method. To identify the largest possible CFL number for a problem, we gradually increase/decrease the values of from an initial value.

Example 1 A 1D Burgers' Equation
We solve the following 1D steady-state Burgers' equation with a source term: The initial condition function u 0 (x) = 2 sin x is used as the initial guess in the iterations, in which case the unique steady-state solution for this problem is u(x) = sin x . An inflow boundary condition is imposed at the left boundary x = (1∕4)π with u((1∕4)π) = √ 2∕2 . The right boundary x = (3∕4)π is the outflow boundary. Both iterative schemes with the fifth-order ILW procedure and fifth-order multi-resolution WENO extrapolation are used to compute the steady-state solution. 10 −13 is taken as the iteration convergence criterion threshold value. Mesh refinement study is performed. L 1 and L ∞ numerical errors and accuracy orders, iteration numbers, and CPU costs of each iterative scheme to converge are presented in Table 1. The results show that both schemes achieve comparable numerical errors and the fifth-order accuracy when they converge. However, the FE fast sweeping scheme is more efficient than the TVD-RK3 scheme in terms of both the iteration numbers and CPU times. In this example, we can see that the FE fast sweeping scheme saves about 60% CPU time cost of that of the RK Jacobi (TVD-RK3) scheme on the most refined grid.

Example 2 A 2D Euler System
The second accuracy test example we consider is the following 2D steady-state Euler system: The exact steady-state solution is (x, y) = 1 + 0.2 sin(x − y), u(x, y) = 1, v(x, y) = 1, and p(x, y) = 1 . The computational domain is (x, y) ∈ [−π, π] × [−π, π] . To start the iterations in the schemes, the exact steady-state solution is used as the initial guess. Since the exact steady-state solution does not satisfy the numerical schemes, it will be driven by the iterative schemes to their numerical steady states. Both iterative schemes with the fifth-order ILW numerical boundary treatment and fifth-order multi-resolution WENO extrapolation are used to compute the steady-state solution. We take 10 −12 as the iteration convergence criterion threshold value and perform mesh refinement study. L 1 and L ∞ numerical errors and accuracy orders, iteration numbers, and CPU costs of each iterative scheme to converge are reported in Table 2. Similar as the last example, we observe that both iterative schemes achieve comparable numerical errors and the fifth-order accuracy when they converge. Again, the FE fast sweeping scheme is more efficient than the RK Jacobi (TVD-RK3) scheme in both the iteration numbers and CPU times, and it saves about 60% CPU time cost of that of the TVD-RK3 scheme in this 2D Euler system example.

Example 3 Regular Shock Reflection
We consider the 2D regular shock reflection problem. Because the residues of highorder WENO schemes in solving this problem are often difficult to settle down to machine zero or the round-off error level (see, e.g., [24,26,27]), this is a typical benchmark problem to test the convergence of high-order WENO methods for solving steady-flow problems. In recent work [10,35], it is found that the iterative schemes based on multi-resolution WENO schemes can drive the iteration residues to the round-off error level. Here, we use this example to test the absolutely convergent fixed-point fast sweeping method with ILW boundary treatment proposed here and verify its efficiency. The computational domain is a rectangle with length 4 and height 1. The flow has the following Dirichlet conditions along the left and the top boundaries: The initial values in the entire domain to start the iterations are set to be those values of the left boundary. The right boundary is the outflow boundary. In the literature (e.g., [10,24,26,27,35]), a reflection condition is applied along the bottom boundary. Here, the fifth-order ILW boundary treatment is applied to the bottom boundary to test the method. Both the fast sweeping method and the RK Jacobi scheme are used to compute the steadystate solution, and 10 −12 is taken as the iteration convergence criterion threshold value. The computational grid is 120 × 30 . Number of iterations, the final time, and the total CPU time when the schemes converge are reported for the RK Jacobi scheme and the FE fast sweeping scheme with the fifth-order ILW boundary treatment in Table 3. It is observed that it takes the FE fast sweeping scheme a much smaller number of iterations to converge to steady state than the RK Jacobi (TVD-RK3) scheme, and about 50% CPU time cost is saved. Figure 2 presents the contour plots of the density variable of the converged steadystate solutions, and the comparable results are obtained for both schemes. Figure 3 shows residue history in terms of iterations for these schemes. We see that for this difficult benchmark problem, the iteration residues converge to small values of round-off error level for both schemes with the fifth-order ILW boundary treatment.

Example 4 Supersonic Flow Past a Cylinder
In this example, we apply the methods to a problem with the physical boundary not aligned with the computational grids, i.e., supersonic flow past a cylinder problem [11,21]. A curved solid wall that is a circular cylinder of unit diameter is positioned at the  . We follow the setup as in [21]. The problem is initialized by a Mach 3 flow moving toward the cylinder wall from the left. At the left inflow boundary x = − 5.0 , the uniform far-field data are imposed. At the top boundary y = 5.0 and the right boundary x = 0 , we use constant extrapolation as in [21]. The simulation is carried out at the upper half of the whole domain by the symmetry of the problem. Hence, at y = 0 , the reflection technique is applied. The fifth-order ILW procedure is used to deal with the surface of the solid wall cylinder. The FE fast sweeping scheme and the RK Jacobi scheme are used to compute the steady-state solution, and 10 −12 is taken as the iteration convergence criterion threshold value. The computational grid sizes are Δx = Δy = 1 32 . Number of iterations, the final time, and the total CPU time when the iterations converge are reported for both schemes with the fifth-order ILW boundary treatment in Table 4. Again, it is observed that the FE fast sweeping scheme is much more efficient than the RK Jacobi scheme to converge to steady state, in terms of both the number of iterations and CPU time costs. About 60% CPU time cost is saved using the FE fast sweeping scheme with the fifth-order ILW boundary treatment. The contour plots of the pressure variable of the converged steady-state solutions are shown in Fig. 4. We observe similar results for both iterative schemes and the well-captured bow shock. Residue history in terms of iterations for these two schemes is presented in Fig. 5. As in the last example, the iteration residues settle down to the 10 −12 level, the round-off error level, for both iterative schemes with the fifth-order ILW boundary method.

Example 5 Supersonic Flow Past a NACA0012 Airfoil
In this example, we consider the problem of an inviscid Euler supersonic flow past a single NACA0012 airfoil configuration [16,35].  schemes with the largest possible CFL numbers are reported in Table 5 for the Ma ∞ = 2 case and Table 6 for the Ma ∞ = 3 case. Contour plots of the pressure variable of the converged steady-state solutions for these schemes are shown in Figs. 6 and 8 for two different Mach number cases. It is observed that the shock waves on the complex surface boundary are well captured. Residue history in terms of iterations is reported in Figs. 7 and 9, which show that the iteration residues settle down to the round-off error 10 −11 − 10 −12 level for both iterative schemes with the fifth-order ILW boundary method. As shown by the results of this example, if the largest CFL number allowed by each scheme to converge to steady states is used, the FE fast sweeping scheme saves about 50% − 70% CPU time of that by the RK Jacobi scheme, while similar numerical steady-state solutions are obtained by them. , and the computational grid is 200 × 200 . As in the supersonic cases of the last example, the fifthorder ILW procedure is used to deal with the complex surface boundary in this NACA0012 airfoil problem. Number of iterations required to reach the convergence threshold value, the final time, and the total CPU time when convergence is obtained for the RK Jacobi and the FE fast sweeping schemes with the largest possible CFL numbers are reported in Table 7 for the Ma ∞ = 0.8 case and Table 8 for the Ma ∞ = 0.5 case. Contour plots of the pressure variable of the converged steady-state solutions for these schemes are shown in Figs. 10 and 12  for these different Mach number cases. As in the previous example, we observe that the shock waves on the complex surface boundary are well captured. In Figs. 11 and 13, residue history in terms of iterations is presented. Again, the iteration residues settle down to the round-off error 10 −11 −10 −12 level for both iterative schemes with the fifth-order ILW boundary procedure. From the simulation results of this example, if the largest CFL number allowed by each scheme to converge to steady states is used, the FE fast sweeping scheme saves about 46% CPU time of that by the RK Jacobi scheme for the Ma ∞ = 0.5 case to obtain similar numerical steady-state solutions, while about 23% CPU time saving is obtained for the Ma ∞ = 0.8 case. Comparing with the supersonic cases of the last example, we observe less CPU time saving for the subsonic cases here. This will be further verified in the next example.

Example 7 Subsonic Flow Past a Cylinder
In this example, we solve the problem of a subsonic flow past a circular cylinder with the Mach number Ma ∞ = 0.38 in [12,34]. The radius of the circular cylinder is 0.5, which is the curved solid wall and located at the center of the computational domain [− 5.0, 5.0] × [− 5.0, 5.0] . The computational grid is 320 × 320 . To deal with the physical boundary which is not aligned with the computational grids, we apply the fifth-order ILW procedure to the complex surface boundary. Number of iterations required to reach the convergence threshold value, the final time, and the total CPU time when convergence is obtained for the RK Jacobi and the FE fast sweeping schemes with the largest possible  Table 9. Contour plots of the Mach number variable of the converged steady-state solutions for both schemes are shown in Fig. 14. We observe that comparable numerical steady-state solutions are obtained by these two schemes. However, it is noticeable that the FE fast sweeping scheme gives a better symmetric numerical steady-state solution around the cylinder than the RK Jacobi scheme. Residue history in terms of iterations provided in Fig. 15 shows that the iteration residues settle down to the round-off error 10 −13 −10 −14 level for both iterative schemes with the fifth-order ILW boundary procedure. From the simulation results of this example with a low Mach number flow, we observe less CPU time saving using the FE fast sweeping scheme than the supersonic cases of the previous examples. Table 9 shows that the number of iterations  required to reach the convergence of the FE fast sweeping scheme is much less than the RK Jacobi scheme. However, since in the fast sweeping iterations, the WENO reconstructions for numerical fluxes at the cell interfaces need to be performed individually for grid points sharing the same cell interface due to the Gauss-Seidel procedure [10], one iteration of the FE fast sweeping scheme has more operations than that of the RK Jacobi scheme. Hence, the CPU time cost of the FE fast sweeping scheme in this example is just slightly less than that of the RK Jacobi scheme. It will be interesting to study the possible reason which causes the FE fast sweeping scheme less efficient for these low Mach number flow problems, for example, by investigating the key factors that fast sweeping methods can accelerate convergence of numerical simulations to steady-state solutions [31]. This will be among our future work on this topic.

Concluding Remarks
In this paper, the absolutely convergent fixed-point fast sweeping WENO method in [10] is combined with the ILW boundary treatment procedure to solve the steady-state solution of hyperbolic conservation laws, especially for problems with the physical boundary not aligned with the computational grids. By applying WENO reconstructions with unequalsized substencils such as the fifth-order multi-resolution WENO scheme, iteration residues of the fast sweeping method with the ILW procedure converge to tiny values at the roundoff error level, for all test problems including difficult cases. Via numerical experiments on various examples, we also verify that the absolutely convergent fixed-point fast sweeping   method with the ILW boundary treatment presented in this paper is more efficient than the popular TVD-RK3 time-marching method to converge to steady-state solutions of a highorder WENO scheme, and more than 50% CPU time costs are saved for most of examples in this paper. In the future work, we will study the approach about how to further improve the computational efficiency of the absolutely convergent fixed-point fast sweeping WENO method with the ILW boundary treatment procedure, especially for some subsonic flow problems.

Compliance with Ethical Standards
Conflict of Interest The authors declare that there is no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.