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

Abstract

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.

1 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 time-march 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 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.

2 Numerical Methods

We consider steady-state problems of hyperbolic conservation laws

$$\begin{aligned} f(u)_x+g(u)_y=R, \qquad (x,y)\in \varOmega \end{aligned}$$

(1)

on a bounded domain $\varOmega$ with appropriate boundary conditions prescribed on $\partial \varOmega$. 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).

2.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 $\varOmega$ by a uniform Cartesian mesh $\{(x_i, y_j)\}$ with the grid size $\Delta {x}$ in the x-direction and $\Delta {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 x-direction flux f(u) is split into $f(u)=f^{+}(u)+f^{-}(u)$, which satisfies the condition that $\frac{{\text{d}}f^{+}(u)}{{\text{d}}u}\geqslant {0}$ and $\frac{{\text{d}}f^{-}(u)}{{\text{d}}u}\leqslant {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^{\pm }(u)=\frac{1}{2}(f(u)\pm \alpha u)$ is used, where $\alpha =\max _{u}|f^{'}(u)|$. Then, each of them is approximated separately by the fifth-order multi-resolution WENO scheme to obtain numerical fluxes ${\hat{f}}^+_{i+1/2,j}$ and ${\hat{f}}^-_{i+1/2,j}$ at the point $(x_{i+1/2}, y_j)$, respectively, where $x_{i+1/2}=(x_i+x_{i+1})/2$. The final numerical flux is $\hat{f}_{i+1/2,j}=\hat{f}_{i+1/2,j}^{+}+\hat{f}_{i+1/2,j}^{-}$. Specifically, the reconstructed numerical flux $\hat{f}_{i+1/2,j}^{+}$ is

$$\begin{aligned} \hat{f}_{i+1/2,j}^{+}=\sum _{l_1=1}^{3}\omega _{l_1}p_{l_1}(x_{i+1/2}). \end{aligned}$$

(2)

The nonlinear weights are computed as

$$\begin{aligned} \omega _{l_1}=\frac{\bar{\omega }_{l_1}}{\sum _{l_2=1}^{3}\bar{\omega }_{l_2}}, \quad \bar{\omega }_{l_1}=\gamma _{l_1,3}\left(1+\frac{\tau }{\varepsilon +\beta _{l_1}}\right), \qquad l_1=1,2,3, \end{aligned}$$

where $\varepsilon$ is taken as $\varepsilon = 10^{-6}$ to avoid the case that the denominator becomes zero. The details on the linear weights $\gamma _{l_1,3}$, computations of smoothness indicators $\beta _{l_1}$, the quantity $\tau$, and reconstruction polynomials $p_{l_1}(x)$ can be found in [10, 33]. The reconstruction of the numerical flux $\hat{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

$$\begin{aligned}&0=-(\hat{f}_{i+1/2,j}-\hat{f}_{i-1/2,j})/ \Delta x-(\hat{g}_{i,j+1/2}-\hat{g}_{i,j-1/2})/\Delta y+R(u_{ij},x_i,y_j),\nonumber \\&\quad i=1,\cdots ,N; \,j=1,\cdots , M, \end{aligned}$$

(3)

where N and M are the numbers of grid points in the x and y directions, respectively. $\hat{f}_{i+1/2,j}$ and $\hat{g}_{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:

$$\begin{aligned} u_{ij}^{n+1}=u_{ij}^{n}+\frac{\gamma }{\alpha _{x}/\Delta x +\alpha _{y}/\Delta y}L(u_{i-r,j}^{*},\cdots ,u_{i+s,j}^{*};u_{i,j}^{n};u_{i,j-r}^{*},\cdots , u_{i,j+s}^{*}), \end{aligned}$$

$$\begin{aligned} i=i_{1},\cdots ,i_{N};\,j=j_{1},\cdots ,j_{M}. \end{aligned}$$

(4)

$u_{ij}^{n}$ is the nth step iteration value at the grid point $(x_i, y_j)$. $\gamma$ is the Courant-Friedrichs-Lewy (CFL) number and $\frac{\gamma }{\alpha _{x}/\Delta {x}+\alpha _{y}/\Delta {y}}$ is corresponding to the time step size $\Delta {t_{n}}$. $\alpha _{x}=\max _{u}{|f'(u)|}$ and $\alpha _{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 fifth-order 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:

$$\begin{aligned} \text{(i) } i=1:N, j=1:M; \\ \text{(ii) } i=N:1, j=1:M; \\ \text{(iii) } i=N:1, j=M:1; \\ \text{(iv) } i=1:N, j=M:1. \end{aligned}$$

This is marked as “$i=i_{1},\cdots ,i_{N};\,j=j_{1},\cdots ,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_{k,l}^{n}$ or $u_{k,l}^{n+1}$, 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 multi-resolution WENO scheme.

2.2 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.

2.2.1 1D Case

First we consider 1D compressible steady-state Euler equations

$$\begin{aligned} F(U)_{x}=R(U,x,y), \end{aligned}$$

(5)

where the conservative variables $U \triangleq (U_1, U_2, U_3)^{\text{T}} = (\rho , \rho u, E)^{\text{T}}$, and the flux $F(U) = (\rho u, \rho u^2 + p, u (E + p))^{\text{T}} = \left(U_{2}, (\gamma '-1)U_{3}+\frac{3-\gamma '}{2}\frac{U_{2}^{2}}{U_{1}}, \left(\gamma ' U_{3}-\frac{\gamma '-1}{2}\frac{U_{2}^{2}}{U_{1}}\right)\frac{U_{2}}{U_{1}}\right)^{\text{T}}$. $\rho ,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))^{\text{T}}$ are the source terms. The equation of state has the form $E=\frac{p}{\gamma '-1}+\frac{1}{2}\rho u^{2},$ where $\gamma '$ is the specific heat ratio and $\gamma '=1.4$ for the case of air. The sound speed is $c=\sqrt{\gamma ' p /\rho }$. Without loss of generality, the computational domain is taken as $x\in (-1,1)$.

The computational domain is partitioned by a uniform mesh

$$\begin{aligned} -1+h/2=x_{0}<x_{1}<\cdots<x_{N}=1-h/2, \end{aligned}$$

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}), \cdots,$ and $U(x_{N+2}),U(x_{N+3}), \cdots$ are treated in a similar way. By a fifth-order Taylor approximation

$$\begin{aligned} U_{m}(x_{N+1})=\sum _{k=0}^{4}\frac{(x_{N+1}-1)^{k}}{k!}U^{*(k)}_{m}, \qquad m=1,2,3, \end{aligned}$$

(6)

where $U^{*(k)}_{m}$ is a $(5-k)$th-order approximation of the spatial derivative $\frac{\partial ^{k}{U_{m}}}{\partial {x^{k}}}\mid _{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 $A_{\perp }(U(x_{N}))=\frac{\partial {F(U)}}{\partial U}\Big|_{U=U(x_{N})}$, and $A_{\perp }(U(x_{N}))$ has three eigenvalues

$$\lambda _{1}=u(x_{N})-c(x_{N}), \lambda _{2}=u(x_{N}), \lambda _{3}=u(x_{N})+c(x_{N})$$

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

$$\begin{aligned} V_{m}(x_{j})=l^{\text{T}}_{m}(U(x_{N}))U(x_{j}), \quad j=N,N-1, \cdots , N-4; \quad m=1,2,3. \end{aligned}$$

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 boundary $x=1$, i.e., $V_{m}^{*(k)}, m=1,2,3;\, k=0,1,2,3,4$. 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_{m}^{*(k)}$ are coupled with the inflow boundary conditions and the PDEs itself to compute these values $U_{m}^{*(k)}$ which are needed in (6). The procedure is summarized as follows. (i) Compute the eigenvalues $\lambda _m$ and left eigenvectors $l_{m}$ of the Jacobian matrix $A_{\perp }(U(x_{N}))$ described above, for $m=1,2,3$. Determine the inflow boundary conditions by the signs of the eigenvalues $\lambda _m$. The inflow boundary conditions are prescribed and given in the problem. (ii) Compute the outgoing local characteristic variables $V_{m}(x_j), j=N-4, \cdots , N$ as described above. Extrapolate $V_{m}(x_j)$ to the boundary to find $V_{m}^{*(k)}, k=0, \cdots , 4$ with the multi-resolution WENO extrapolation described below. (iii) Solve for $U_{m}^{*(0)}, m=1,2,3$, using the prescribed inflow boundary conditions and the extrapolated values $V_{m}^{*(0)}$. (iv) For $k=1, \cdots , 4$, 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_{m}^{*(k)}, m=1,2,3$. Solve the linear system to obtain $U_{m}^{*(k)}, 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_{m}^{*(k)}$ at the boundary $x=1$, using the numerical values at the interior grid points of the domain $V_{m}(x_{j}), j=N-4, \cdots , 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 $S_{1}=\{x_{N}\}, S_{2}=\{x_{N-2}, x_{N-1}, x_{N}\},$ and $S_{3}=\{x_{N-4}, \cdots , x_{N}\}$. The extrapolation procedure is given in the following.

Step 1 On each substencil, a Lagrange polynomial $q_{r}(x)$ of degree $2r-2$ is constructed, such that $q_{r}(x_{j})=V_{j}, j=N-2r+2, \cdots , N;\, r=1,2,3$.

Step 2 Obtain equivalent expressions for these polynomials as in multi-resolution WENO reconstruction procedure [33]. Namely, take $p_{1}(x)=q_{1}(x)$, and construct

$$\begin{aligned} p_{2}(x)=\frac{1}{\gamma _{2,2}}q_{2}(x)-\frac{\gamma _{1,2}}{\gamma _{2,2}}p_{1}(x), \\ p_{3}(x)=\frac{1}{\gamma _{3,3}}q_{3}(x)-\frac{\gamma _{2,3}}{\gamma _{3,3}}p_{2}(x)-\frac{\gamma _{1,3}}{\gamma _{3,3}}p_{1}(x) \end{aligned}$$

with $\sum _{l=1}^{2}\gamma _{l,2}=1$, $\sum _{l=1}^{3}\gamma _{l,3}=1,$ and $\gamma _{2,2}\ne 0,\gamma _{3,3}\ne 0.$ We take $\gamma _{1,2}=1/11$, $\gamma _{2,2}=10/11$, $\gamma _{1,3}=1/111$, $\gamma _{2,3}=10/111$, $\gamma _{3,3}=100/111$ as in [33].

Step 3 Compute the smoothness indicators $\beta _{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],

$$\begin{aligned} \beta _{l}=\sum _{\alpha =1}^{2l-2}\int _{x_{N}}^{x_{N+1}}h^{2\alpha -1}\left(\frac{{\text{d}}^{\alpha }p_{l}(x)}{{\text{d}}x^{\alpha }}\right)^{2}{\text{d}}x, \qquad l=1,2,3. \end{aligned}$$

Notice here $\beta _{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 $\beta _{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 $\tau$ which depends on the absolute differences between the smoothness indicators is calculated as $\tau =\left(\frac{|\beta _{3}-\beta _{1}|+|\beta _{3}-\beta _{2}|}{2}\right)^2$. The nonlinear weights are

$$\begin{aligned} \omega _{l_1}=\frac{\bar{\omega }_{l_1}}{\sum _{l_2=1}^{3}\bar{\omega }_{l_2}}, \quad \bar{\omega }_{l_1}=\gamma _{l_1,3}\left(1+\frac{\tau }{\varepsilon +\beta _{l_1}}\right), \qquad l_1=1,2,3. \end{aligned}$$

Here, $\varepsilon = 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

$$\begin{aligned} V^{*(k)}=\sum _{r=1}^{3}\omega _{r}\frac{{\text{d}}^{k}p_{r}(x)}{{\text{d}}x^{k}}\bigg |_{x=1}. \end{aligned}$$

(7)

2.2.2 2D Case

We consider 2D compressible steady-state Euler equations

$$\begin{aligned} F(U)_{x}+G(U)_{y}=0, \qquad (x,y)\in \varOmega \end{aligned}$$

(8)

with appropriate boundary conditions, where $U=(U_{1},U_{2},U_{3},U_{4})^{\text{T}}=(\rho ,\rho {u},\rho {v},E)^{\text{T}}, F(U)=(\rho {u},\rho {u^{2}}+p,\rho {u}v,u(E+p))^{\text{T}}, G(U)=(\rho {v},\rho {u}v,\rho {v^{2}}+p,v(E+p))^{\text{T}}$. $\rho , 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=\frac{p}{\gamma '-1}+\frac{1}{2}\rho (u^{2}+v^{2})$, where $\gamma '$=1.4 for the case of air.

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 $\varGamma =\partial \varOmega$, 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 $\varGamma$. We set up a local coordinate system $(\hat{x}, \hat{y})$ at $P_{0}$ by

$$\begin{aligned} \left( \begin{array}{cc} \hat{x} \\ \hat{y} \end{array}\right) = \left( \begin{array}{cc} \ \cos \theta \ \ \ \sin \theta \\ -\sin \theta \ \ \cos \theta \end{array}\right) \left( \begin{array}{cc} x\\ y \end{array}\right) =T\left( \begin{array}{cc} x\\ y \end{array}\right) , \end{aligned}$$

where $\theta$ is the angle between the outward normal $n(P_{0})$ and the x-axis and T is the corresponding rotation matrix. The $\hat{x}$-axis has the same direction as $n(P_{0})$ and the $\hat{y}$-axis is the tangential direction. Under this local coordinate system, the Euler equations (8) are transformed to $F(\hat{U})_{\hat{x}}+G(\hat{U})_{\hat{y}}=0$, where $\hat{U}=(\hat{U}_{1},\hat{U}_{2},\hat{U}_{3},\hat{U}_{4})^{\text{T}}$. For the fifth-order boundary treatment, the numerical value at the ghost point P is approximated by the Taylor expansion

$$\begin{aligned} (\hat{U}_{m})_{P}=\sum _{k=0}^{4}\frac{\Delta ^{k}}{k!}\hat{U}_{m}^{*(k)}, \qquad m=1, \cdots , 4. \end{aligned}$$

(9)

Here, $\Delta$ is the $\hat{x}$-coordinate of the ghost point P and $\hat{U}_{m}^{*(k)}$ is a $(5-k)$th-order approximation of the normal derivative $\frac{\partial ^{k}\hat{U}_{m}}{\partial \hat{x}^{k}}\big |_{(x,y)=P_{0}}$. 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 $\hat{U}_{(0)}$ of a grid point which is the closest to $P_{0}$ among all the grid points inside the domain $\varOmega$, to form the Jacobian matrix of the normal flux by $A_{\bot }(\hat{U}_{(0)})=\frac{\partial F(\hat{U})}{\partial (\hat{U})}\big |_{\hat{U}=\hat{U}_{(0)}}$. $A_{\bot }(\hat{U}_{(0)})$ has four eigenvalues $\lambda _{1}=\hat{u}_{0}-c_{0}$, $\lambda _{2}=\lambda _{3}=\hat{u}_{0}$, $\lambda _{4}=\hat{u}_{0}+c_{0}$ and a complete set of left eigenvectors $l_{m}(\hat{U}_{(0)})$, $m=1,\cdots ,4$, where $\hat{u}_{0}$ and $c_0$ are the corresponding x-velocity and the sound speed. The local characteristic variables $V_{m}$ at interior grid points near $P_{0}$ are given by

$$\begin{aligned} (V_{m})_{\mu ,\nu }=l_{m}(\hat{U}_{(0)})\hat{U}_{\mu ,\nu }, \quad m=1, \cdots , 4, \quad (x_{\mu },y_{\nu })\in S, \end{aligned}$$

where S is the stencil for the fifth-order multi-resolution WENO extrapolation, which is used to compute ${V}_{m}^{*(k)}$, the $(5-k)$th-order approximation of $\frac{\partial ^{k}V_{m}}{\partial \hat{x}^{k}}\big |_{P_{0}}$. The fifth-order multi-resolution WENO extrapolation will be described in the following. Then, following the procedure in [21, 22], these extrapolation values $V_{m}^{*(k)}$ are coupled with the inflow boundary conditions and the PDEs itself to compute these values $\hat{U}_{m}^{*(k)}$ which are needed in the Taylor expansion (9).

Fig. 1

The stencil S for 2D fifth-order multi-resolution WENO extrapolation

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})$ ($\hat{x}$-axis) are selected as the stencil S of the WENO extrapolation. For the fifth-order 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 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)$, $p_{2}(x,y)=\frac{1}{\gamma _{2,2}}q_{2}(x,y)-\frac{\gamma _{1,2}}{\gamma _{2,2}}p_{1}(x,y)$, and $p_{3}(x,y)=\frac{1}{\gamma _{3,3}}q_{3}(x,y)-\frac{\gamma _{1,3}}{\gamma _{3,3}}p_{1}(x,y)-\frac{\gamma _{2,3}}{\gamma _{3,3}}p_{2}(x,y)$. The linear weights are taken as $\gamma _{1,2}=1/11$, $\gamma _{2,2}=10/11$, and $\gamma _{1,3}=1/111,\gamma _{2,3}=10/111,\gamma _{3,3}=100/111$. Then, the multi-resolution WENO extrapolation gives

$$\begin{aligned} V_{m}^{*(k)}=\sum ^{3}_{r=1}\omega _{r}\frac{\partial ^{k}p_{r}(x,y)}{\partial {\hat{x}^{k}}}\bigg |_{(x,y)=(x_{0},y_{0})}, \qquad k=0, \cdots , 4. \end{aligned}$$

(10)

Here, the nonlinear weights are computed as follows:

$$\begin{aligned} \omega _{l_1}=\frac{\bar{\omega }_{l_1}}{\sum _{l_2=1}^{3}\bar{\omega }_{l_2}}, \quad \bar{\omega }_{l_1}=\gamma _{l_1,3}\left(1+\frac{\tau }{\varepsilon +\beta _{l_1}}\right), \qquad l_1=1,2,3. \end{aligned}$$

The quantity $\tau =\left(\frac{|\beta _{3}-\beta _{1}|+|\beta _{3}-\beta _{2}|}{2}\right)^{2}$, and the smoothness indicators

$$\begin{aligned} \beta _{l}=\sum _{\mid \alpha \mid =1}^{2l-2}\int _{x_{0}-\frac{\Delta {x}}{2}}^{x_{0}+\frac{\Delta {x}}{2}}\int _{y_{0}-\frac{\Delta {y}}{2}}^{y_{0}+\frac{\Delta {y}}{2}}(\Delta {x}\Delta {y})^{\mid \alpha \mid -1}({\text{D}}^{\alpha }p_{l}(x,y))^{2}{\text{d}}x{\text{d}}y, \qquad l=1,2,3. \end{aligned}$$

Here, $\alpha$ is the multi-index notation for partial derivative operator ${\text{D}}^{\alpha }$, and $\varepsilon = 10^{-6}$ is used to avoid the case that the denominator becomes zero. As that in the 1D case, here, the smoothness indicator $\beta _{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 $\beta _{1}$.

3 Numerical Experiments

In this section, we perform numerical experiments to test the performance of the fixed-point 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 fixed-point 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

$$\begin{aligned} {\text {Res}}\,A=\sum _{i=1}^{{\mathcal {N}}}\frac{|R_{i}|}{{\mathcal {N}}}, \end{aligned}$$

(11)

where the local residue at the grid point i is

$$\begin{aligned} R_{i}=\frac{u_{i}^{n+1}-u_{i}^{n}}{\Delta t_{n}}. \end{aligned}$$

${\mathcal {N}}$ is the total number of grid points and n is the iteration step. $\Delta t_n=\frac{\gamma }{\alpha _x/\Delta x+\alpha _y/\Delta y}$. For the 2D Euler system, the average residue is defined as $\text {Res}\,A=\sum _{i=1}^{{\mathcal {N}}}\frac{|R1_{i}|+|R2_{i}|+|R3_{i}|+|R4_{i}|}{4 {\mathcal {N}}}$. Here, $R*_{i}$’s are the local residuals of the conservative variables, namely, $R1_{i}=\frac{\rho _{i}^{n+1}-\rho _{i}^{n}}{\Delta t_n}, R2_{i}=\frac{(\rho u)_{i}^{n+1}-(\rho u)_{i}^{n}}{\Delta t_n}, R3_{i}=\frac{(\rho v)_{i}^{n+1}-(\rho v)_{i}^{n}}{\Delta t_n}, R4_{i}=\frac{E_{i}^{n+1}-E_{i}^{n}}{\Delta t_n}$. The convergence criterion is set to be ${\text{Res}}\,A \leqslant \delta$, where the threshold value $\delta$ 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 $\Delta 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 $\gamma$ 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 $\gamma$ from an initial value.

Example 1

A 1D Burgers’ Equation

We solve the following 1D steady-state Burgers’ equation with a source term:

$$\begin{aligned} \left( \frac{u^{2}}{2}\right) _{x}=\sin {x}\cos {x} ,\qquad x\in \left[ \frac{1}{4}{\uppi} , \frac{3}{4}{\uppi} \right] . \end{aligned}$$

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){\uppi}$ with $u((1/4){\uppi} )=\sqrt{2}/2$. The right boundary $x=(3/4){\uppi}$ 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_{\infty }$ 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.

Table 1 Example 1: a 1D Burgers’ equation

Example 2

A 2D Euler System

The second accuracy test example we consider is the following 2D steady-state Euler system:

$$\begin{aligned} \frac{\partial }{\partial {x}} \left( \begin{array}{cccc} \rho u\\ \rho u^{2}+p\\ \rho uv\\ u(E+p) \end{array}\right) +\frac{\partial }{\partial {y}} \left( \begin{array}{cccc} \rho v\\ \rho uv\\ \rho v^{2}+p\\ v(E+p) \end{array}\right) = 0. \end{aligned}$$

The exact steady-state solution is $\rho (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) \in [-{\uppi} ,{\uppi} ] \times [-{\uppi} ,{\uppi}]$. 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_{\infty }$ 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.

Table 2 Example 2: a 2D Euler system

Example 3

Regular Shock Reflection

We consider the 2D regular shock reflection problem. Because the residues of high-order 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:

$$\begin{aligned} (\rho ,u,v,p)^{\text{T}}= {\left\{ \begin{array}{ll} (1.0,2.9,0,5/7)^{\text{T}}\mid _{(0,y)},\\ (1.699\,97,2.619\,34,-0.506\,32,1.528\,19)^{\text{T}}\mid _{(x,1)}.\\ \end{array}\right. } \end{aligned}$$

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 steady-state solution, and $10^{-12}$ is taken as the iteration convergence criterion threshold value. The computational grid is $120 \times 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 steady-state 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.

Table 3 Example 3: regular shock reflection

Fig. 2

Example 3: regular shock reflection. Thirty equally spaced density contours from 1.1 to 2.6 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 3

Example 3: regular shock reflection. The convergence history of the residue as a function of number of iterations for two schemes

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 origin (0, 0) of the computational domain $[-\,5.0, 0]\times [-\,5.0, 5.0]$. 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 $\Delta x=\Delta y=\frac{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.

Table 4 Example 4: supersonic cylinder problem

Fig. 4

Example 4: supersonic cylinder problem. Thirty equally spaced pressure contours from 0 to 12 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 5

Example 4: supersonic cylinder problem. The convergence history of the residue as a function of number of iterations for two schemes

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]. Two cases of supersonic flows are tested, i.e., a flow with Mach number $Ma_{\infty }=2$ and angle of attack $\alpha =1^{\circ }$; and a flow with Mach number $Ma_{\infty }=3$ and angle of attack $\alpha =3^{\circ }$. The computational domain is $[-\,1.5,1.5]\times [-\,1.5,1.5]$, and the computational grid is $200\times 200$. The fifth-order ILW procedure is used to deal with the complex surface boundary in this 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 $\gamma$ are reported in Table 5 for the $Ma_{\infty }=2$ case and Table 6 for the $Ma_{\infty }=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.

Table 5 Example 5: supersonic flow past a single NACA0012 airfoil

Fig. 6

Example 5: supersonic flow past a single NACA0012 airfoil. $Ma_{\infty }=2$. Thirty equally spaced pressure contours from 0.4 to 4.2 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 7

Example 5: supersonic flow past a single NACA0012 airfoil. $Ma_{\infty }=2$. The convergence history of the residue as a function of number of iterations for two schemes

Table 6 Example 5: supersonic flow past a single NACA0012 airfoil

Fig. 8

Example 5: supersonic flow past a single NACA0012 airfoil. $Ma_{\infty }=3$. Thirty equally spaced pressure contours from 0.4 to 8.0 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 9

Example 5: supersonic flow past a single NACA0012 airfoil. $Ma_{\infty }=3$. The convergence history of the residue as a function of number of iterations for two schemes

Example 6

Subsonic Flow Past a NACA0012 Airfoil

Following the last example, here, we test the cases of subsonic flow past a single NACA0012 airfoil configuration. Two cases of subsonic flows are considered, i.e., a flow with Mach number $Ma_{\infty }=0.8$ and angle of attack $\alpha =1.25^{\circ }$; and a flow with Mach number $Ma_{\infty }=0.5$ and angle of attack $\alpha =1.25^{\circ }$. The computational domain is $[-\,1.0,2.0]\times [-\,1.5,1.5]$, and the computational grid is $200\times 200$. As in the supersonic cases of the last example, the fifth-order 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 $\gamma$ are reported in Table 7 for the $Ma_{\infty }=0.8$ case and Table 8 for the $Ma_{\infty }=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_{\infty }=0.5$ case to obtain similar numerical steady-state solutions, while about $23\%$ CPU time saving is obtained for the $Ma_{\infty }=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.

Table 7 Example 6: subsonic flow past a single NACA0012 airfoil

Fig. 10

Example 6: subsonic flow past a single NACA0012 airfoil. $Ma_{\infty }=0.8$. Thirty equally spaced pressure contours from 0.3 to 1.35 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 11

Example 6: subsonic flow past a single NACA0012 airfoil. $Ma_{\infty }=0.8$. The convergence history of the residue as a function of number of iterations for two schemes

Table 8 Example 6: subsonic flow past a single NACA0012 airfoil

Fig. 12

Example 6: subsonic flow past a single NACA0012 airfoil. $Ma_{\infty }=0.5$. Sixty equally spaced pressure contours from 0.5 to 1.1 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 13

Example 6: subsonic flow past a single NACA0012 airfoil. $Ma_{\infty }=0.5$. The convergence history of the residue as a function of number of iterations for two schemes

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_{\infty } = 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]\times [-\,5.0,5.0]$. The computational grid is $320\times 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 CFL numbers $\gamma$ are reported in 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.

Table 9 Example 7: subsonic flow past a circular cylinder

Fig. 14

Example 7: subsonic flow past a circular cylinder. Forty-four equally spaced Mach number contours from 0.0 to 0.92 of the converged steady states of numerical solutions by two different iterative schemes

Fig. 15

Example 7: subsonic flow past a circular cylinder. The convergence history of the residue as a function of number of iterations for two schemes

4 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 unequal-sized 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 round-off 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 high-order 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.