1 Introduction

Shallow water equations (SWEs) are commonly used to describe the motion of the “shallow” free-surface flows subjected to gravitational force. They play a critical role in the modelling and simulation of the flows in rivers, channels, ocean tides, etc. Under the assumption of incompressible fluid and hydrostatic pressure distribution, when the vertical acceleration of the water particles is neglected, the SWE may be derived by depth-integrating the Navier-Stokes equations as follows^[1]:

(1)

where t denotes the time. The effect of the bed slope on the flow has been modelled by the inclusion of the source terms at the right-hand side of Eq. (1), which modifies the momentum equations. In the two-dimensional (2D) Cartesian coordinate x=(x, y), the conservative vector U, the flux vector F, and the source term S in Eq. (1) are given by

where u=(u, v) denotes the velocity vector, g is the acceleration due to gravity, I is the identity matrix, and h is the water depth, i.e., the height of water above the riverbed topography z(x). Dropping the time derivatives in Eq. (1) gives the the steady-state SWE as follows:

(2)

If the water is at rest, i.e., u(x)=0, then the momentum parts in the above equation can be reduced to

(3)

which implies the so-called well-balanced property, i.e., “steady state of the water at rest”,

Up to now, there exist various numerical methods for SWEs, e.g., the finite difference scheme based on the flux-difference splitting^[2], the generalized Riemann problem scheme^[3], the highorder weighted essentially non-oscillatory schemes^[4-5], the gas-kinetic schemes^[6-7], the moving mesh method^[8], and the Runge-Kutta discontinuous Galerkin methods^[9-10]. Most of them are explicit for the time-dependent SWEs, and can successfully simulate the evolution of the time-dependent solutions with good accuracy in time. If the explicit time advancing method is used to investigate the steady-state behavior of the flow, then the long time simulation will be needed and become time-consuming due to the small time step size satisfying a CourantFriedrichs-Lewy condition to guarantee the stability. For example, in the sediment transport and morphodynamic change model^[11-12], the time stiffness arising from the characteristic time scales in the flow and sediment transport seriously challenges the long time simulation if the interaction of the water flow with the bed topography is very weak.

For such cases, the implicit or semi-implicit schemes are attractive, such as the multigrid semi-implicit finite difference method^[13], the linearized implicit scheme with a modified Roe flux^[11], the space-time discontinuous residual distribution scheme^[14], the implicit higher-order compact scheme^[15], and the semi-implicit discontinuous Galerkin methods^[16-17]. In fact, the time step size for an implicit scheme is also often constrained by convergence. Even for the same time step size as used in explicit cases, the unsteady solutions of the implicit schemes may be less accurate. An implicit scheme usually requires to solve a nonlinear equation by some iteration method. Therefore, it is also very time-consuming. For the steady-state behavior of the SWE flow, another way is to directly solve the steady-state SWEs^[18]. Therefore, developing a robust and efficient solver for the corresponding nonlinear algebraic system is the key. Unfortunately, there are few such studies for the steady-state SWEs. The steady-state Euler equations and the Navier-Stokes equations have been well solved numerically^[19-23]. For example, a multigrid block lower-upper symmetric Gauss-Seidel (BLU-SGS) algorithm was proposed for the 2D steadystate Euler equations on the unstructured grid^[22]. Unlike the existing methods which added the pseudo-time terms to the steady-state equations, the norm of the local residual in each cell was used to regularize the nonlinear algebraic system arising from the spatial discretization of the steady-state Euler equations. The Newton iteration was then adopted to solve the nonlinear algebraic system, and the multigrid method was used as the inner iteration with the BLU-SGS smoother.

The aim of the paper is to extend the Newton multigrid method^[22] to the steady-state SWEs and investigate the convergence, in which the steady-state SWEs are discretized by a well-balanced hydrostatic reconstruction, the wet/dry transition is numerically handled, the resulting nonlinear algebraic system is iteratively solved by the Newton multigrid iteration, and the convergence behavior of the method for different numerical fluxes is investigated in the numerical experiments through the distribution of the eigenvalues of the iteration matrix. The paper is organized as follows. Section 2 presents the Newton multigrid method, including the well-balanced spatial discretization of the steady-state SWEs in Subsection 2.1, the regularization of the resulting nonlinear system and its Newton iteration linearization in Subsection 2.2, the geometric multigrid method in Subsection 2.3, and the solution procedure of the Newton multigrid method in Subsection 2.4. Section 3 conducts several numerical experiments to demonstrate the efficiency and robustness of the proposed Newton multigrid method, and detailedly discusses the relation between the convergence behavior of the proposed method and the eigenvalue distribution of the block symmetric Gauss-Seidel (SGS) iteration matrix for different numerical fluxes. Section 4 concludes the paper with several remarks.

2 Numerical method

This section is devoted to present the Newton multigrid method for Eq. (2).

Let be a partition of the spatial domain Ω_p, and _i ∈ be the ith cell whose centroid is x_i. Let denote the edge set of _i, and e_ij be the edge of _i sharing with the neighboring cell _j, i.e., e_ij= ∩ . Let n_ij=(n_ij^x, n_ij^y)^T be the unit normal vector of e_ij pointing from _i to _j, and |e_ij| denote the length of the edge e_ij.

2.1 Well-balanced finite volume discretization

This section presents the well-balanced finite volume discretization of Eq. (2).

Integrating Eq. (2) over the cell _i and using the divergence theorem give

(4)

where F_nij(U) :=F (U) · n_ij. If u=0, then Eq. (4) can be reduced to

(5)

where ñ_ij=(0, n_ij^x, n_ij^y)^T. It is a starting point of the hydrostatic reconstruction method^[24] to derive a well-balanced method.

Reconstruct a piecewise polynomial U_h(x) to approximate the solution U(x), e.g.,

(6)

where U_i is the cell average approximation of U(x) over _i. Replacing U(x) in Eq. (4) with U_h(x), using the midpoint quadrature to evaluate the integral, and approximating the values of F_nij(U) and S_ij(h) at the midpoint x_eij of the edge e_ij by numerical fluxes give the finite volume discretization of Eq. (2) as follows:

(7)

where

(8)

In the above equation, U_eij, ±=: (h_eij, ±, h_eij, ±u_eij, ±)^T denote the left and right limit values of U_h(x) at the point x_eij in the direction n_ij, _nij(U_{eij, −}, U_{eij, +}) is any given numerical flux satisfying

and

where

(9)

Such locally reconstructed heights can roughly capture the dry regions, where h=0^[24].

Remark 1  For one-dimensional (1D) steady-state SWEs, the first-order accurate wellbalanced scheme may be described as follows:

(10)

where the numerical flux and the source term S_i are

(11)

respectively, and the left and right limits of the approximate solution at are

(12)

Remark 2  Three numerical fluxes will be considered in this work. The first flux is the Harten-Lax-van Leer contact (HLLC) (resp.Harten-Lax-van Leer (HLL)) flux for 2D (resp.1D) SWEs^[25-26]. The 2D HLL flux is defined by

(13)

where

In the above equations,

Based on the rotational invariance property of 2D SWEs^[27], we have

where

The HLLC flux may be given by

where

The second flux is the local Lax-Friedrichs (LLF) flux

(14)

where s_max denotes an estimation of the fastest wave speed in the local 1D Riemann problem solution, and is usually taken as an upper bound for the absolute value of the eigenvalues of the Jacobian , i.e.,

If there exists wet/dry transition, then because the speed of a wet/dry front is of the form S_*+=u_{nij, +} − for a left dry state and S_*−=u_{nij, −} + for a right dry state^[28], s_max should be appropriately larger in the presence of wet/dry fronts to avoid the numerical instabilities and taken as follows:

where ε_r=0.03 in our calculation.

The third flux is the Roe flux with Harten’s entropy fix

(15)

where m is equal to 2 (resp. 3) for the 1D (resp. 2D) case, =(U_{eij, −}, U_{eij, +}) and ȓ_k= ȓ_k(U_{eij, −}, U_{eij, +}) are the eigenvalues and the corresponding right eigenvectors of the Roe matrix, respectively, can be obtained by

and

in which ε_f is a small positive constant, e.g., 0.4 in the later computation.

Remark 3  Let us discuss the boundary conditions for the discrete problem (7). If the boundary is open, according to the local Froude number Fr :=|u|/, the numerical boundary conditions can be specified as (i) the subcritical inflow boundary (Fr < 1, u_n < 0): _g⁻=_I⁻, and h_gu_n, _g and u_{τ, g} are prescribed; (ii) the subcritical outflow boundary (Fr < 1, u_n > 0): _g⁻=_I⁻, u_{τ, g}=u_{τ, I}, and h_g is prescribed; (iii) the supercritical inflow boundary (Fr > 1, u_n < 0): h_g, u_{n, g}, and u_{τ, g} are prescribed; (iv) the supercritical outflow boundary (Fr > 1, u_n > 0): h_g=h_I, u_{n, g}=u_{n, I}, and u_{τ, g}=u_{τ, I}, where ^∓ :=u_n± denote 1D Riemann invariants associated with the eigenvalues u_n∓, and u_n and u_τ are the normal and tangential velocity components to the cell interface located on the domain boundary, respectively. The quantities with the subscripts g and I denote the values from the ghost and interior cells adjacent to the domain boundary, respectively. The slip boundary conditions (h_g=h_I, u_{n, g}=0, and u_{τ, g}=u_{τ, I}) or the reflective boundary conditions (h_g=h_I, u_{n, g}=−u_{n, I}, and u _{τ, g}=u_{τ, I}) may be specified on the wall.

2.2 Newton’s iterative method

As soon as the boundary conditions are specified, the approximate solutions of Eq. (2) may be obtained by iteratively solving the nonlinear algebraic system (7) with respect to the unknown variables U_i. Here, Newton’s iteration method is employed to solve Eq. (7) with the following formula:

(16)

where the unknown δU_j⁽ⁿ⁾ is the increment, U_i⁽⁰⁾ is the given initial guess, and R_i⁽ⁿ⁾ is the local residual at the nth Newtonian iterative step defined by

(17)

The partial derivatives are

(18)

(19)

contributing to the Jacobian matrix in the Newtonian method approximately calculated by the following numerical differentiation:

(20)

which denotes the derivative of the sth component of the vector , denoted by , with respect to the kth component of the vector U_i, where e_k is the kth column vector of the identity matrix of the same size as . If an edge of the interior cell _i is located on the boundary of Ω_p, then Eq. (20) should be replaced with

where U_g denotes the approximate cell-average value of U in the ghost cell, and is considered as a function of U_i⁽ⁿ⁾ according to the boundary conditions given in Remark 3.

Although the calculation of can also be obtained by using the chain rule, it may be seriously tedious. Moreover, their analytical expressions are difficultly derived near the domain boundary.

The linear system (16) is generally singular, and should be regularized. One way is to add an artificial time derivative term into Eq. (16).

An alternative approach is to use the ℯ¹-norm of the local residual to regularize Eq. (16) as follows:

(21)

where α is the positive regularization parameter. Such a regularization technique has been used in solving steady-state Euler equations^[22].

Solving the linear system (21) for the unknown δU_j⁽ⁿ⁾ by the multigrid method will be discussed in Subsection 2.3. If the solution of the linear system (21) is gotten, then the approximate solution of Eq. (7) can be updated by

(22)

where τ_i is a relaxation parameter on the cell _i.

2.3 Geometric multigrid solver

This section extends the geometric multigrid method^[22] to the linear system (21).

The geometric multigrid methods described so far need a hierarchy of geometric grids or meshes {_l : l=0, 1, …, N_L} for the spatial domain Ω_p, from the coarsest one (l=N_L) to the finest one (l=0). On all levels but the coarsest one, the smoother will be applied. On the coarsest level, the system is usually solved exactly. Assume that the ratio of the grid points on the “neighboring” grids is constant throughout the grid hierarchy, and each coarser cell _{i, l+1} ∈ _l+1 is a union of several neighboring finer cells (two cells for 1D case and four cells in 2D case) in the mesh _l, i.e.,

where I_{i, l+1} is the corresponding index set of those finer cells _{j, l}.

Under the above assumptions, Eq. (21) is discretized and solved by the Newton iteration on the finest mesh ₀ (see Subsections 2.1 and 2.2). The linear system (21) on ₀ is reformulated as follows:

(23)

where the subscript l marks the mesh level, the superscript (n) in Eq. (21) has been omitted for convenience,

when j=i, and

when j ≠ i.

On the coarser mesh, the coarse mesh matrices A_{ij, l+1} are defined by the Galerkin projection, and the source term R_{i, l+1} is derived by the restriction operator I_l^l+1 that restricts those on the fine mesh _l to the coarse mesh _l+1 (l≥0). In this paper, specifically, they are

(24)

where δU_{j, l} is the (approximate) solution of Eq. (23) if l=0, otherwise, the solution can be obtained by solving the following linear system:

(25)

As soon as the correction δU_{j, l+1} on the coarse mesh _l+1 is obtained, the correction δU_{j, l} on the fine mesh _l will be improved as follows:

(26)

where I_l+1^l denotes the prolongation or coarse-to-fine operator that prolongates or interpolates the correction to the fine mesh from the coarse mesh.

A multigrid cycle can be defined as a recursive procedure applied at each mesh level as it moves through the grid hierarchy. For example, multigrid methods with γ-cycle have the following compact recursive definition.

Algorithm 0 δU_{j, l} ← MG_l^γ (δU_{j, l}, R_{j, l})

(i) Pre-smoothing: apply the smoother ν₁ times to Eq. (23) or Eq. (25) with the initial guess δU_{j, l}.

(ii) If _l is the coarsest grid, i.e., l=N_L, solve the problem (25) with l=N_L. Else, restrict to the next coarser grid _l+1 by Eq. (24), and set the initial increment on the next coarser grid as follows: δU_{j, l+1}=0. If _l is the finest grid, set γ=1, and call the γ-cycle scheme γ times for the next coarser grid _l+1 as follows:

(iii) Correct with the prolongated update (26).

(iv) Post-smoothing: apply the smoother ν2 times to Eq.(23) or Eq.(25) with the initial guess δU_{j, l}.

This paper only focuses on two types of multigrid cycles, the V-cycle (γ=1) and the Wcycle (γ=2). Figure 1 shows their schematic descriptions with N_L=3, where symbols “◦” and “•” denote pre-smoothing and post-smoothing, respectively, while the oblique lines between two symbols “◦” (resp.“•”) correspond to the restriction I_l^l+1 (resp.the prolongation I_l+1^l) steps. The smoother is taken as the block SGS iteration.

2.4 Solution procedure

This section summarizes the solution procedure of our Newton multigrid method for Eq.(2). It is well-known that the convergence of Newton’s iteration seriously depends on the choice of the initial guess. To overcome this difficulty, the initial guess is obtained by using the improved BLU-SGS method^[19] so as to solve Eq.(25) from the coarsest mesh _NL to the mesh ₁ successively. The idea of the BLU-SGS method is to retain the block diagonal matrices but employ BLU-SGS-like backward and forward Gauss-Seidel iterations to include the implicit contributions from the off diagonal blocks. This method may be regarded as a “nonlinear extension” of the block SGS method for solving the nonlinear algebraic system (7).

The detailed solution procedure is illustrated by the following flowchart.

Algorithm 1   Newton multigrid method (NMGM)

Step 1   Initialization

Give the “initial data” U(x) and successively refined partitions {_l : l=N_L, …, 1, 0} of the spatial domain Ω_p.

(i) Compute U_{i, NL}, the cell average value of U(x) over the coarsest cell _{i, NL}.

(ii) Let l=N_L, …, 1. Then, perform the BLU-SGS iteration on the mesh _l by

until

where I denotes the identity matrix, and prolongate the solution to the fine mesh from the coarse mesh by U_{j, l−1}=U_{i, l} for all j ∈ I_{i, l}.

(iii) Set n=0 and the initial guess for Newton’s iteration by U_{j, 0}⁽⁰⁾=U_{i, 0} for all j ∈ I_{i, 1}.

Step 2   Newton multigrid iteration

Let n=1, 2, …, N_step. Then, do as follows:

(i) Perform the BLU-SGS iteration on the finest mesh ₀ by

(ii) Solve Eq.(23) by calling Algorithm 0 N_mg times.

(iii) Update the solution U_{i, 0}⁽ⁿ⁺¹⁾=U_{i, 0}⁽ⁿ⁾ + τ_iδU_{i, 0}⁽ⁿ⁾.

Step 3   Check < ε. If yes, output the results and stop; otherwise, set n ← n+1, and go to Step 2.

Before ending this section, several remarks are given below.

Remark 4  The parameter ∊ _p in Step 1 of Algorithm 1 should be chosen appropriately. If ∊ _p is very small, the cost of Step 1 becomes huge so that the steady-state SWE solver is inefficient. According to the numerical experiences, Algorithm 1 works satisfactory if ∊ _p is chosen about one percent of the residual given by the initial data.

Remark 5  If the approximate solution U_{i, 0}(n) given in the Newton multigrid iteration satisfies h_{i, 0}(n) < h_∊, the cell _{i, 0} can be temporarily regarded as a dry cell, and the value of U_{i, 0}(n) can be reset as zero. In our calculations, h_∊=10⁻⁶.

Remark 6  The Newton multigrid scheme in Algorithm 1 can also be extended to solving the nonlinear system arising from an implicit or semi-implicit scheme for the time-dependent SWE (1), e.g.,

where Δt_n denotes the time step size, and the weight β ∈ [0, 1).

3 Numerical experiments

The section presents several numerical examples to demonstrate the robustness and efficiency of the NMGM for 1D and 2D steady-state SWEs (2), and investigates the relation between the convergence behavior of the NMGM and the distribution of the eigenvalues of the iteration matrix detailedly. Unless specifically stated, the parameter ∊ _p in the initialization step of Algorithm 1 is taken as 0.2, and the multigrid iteration number N_mg is set to be 2 (resp.3) for 1D (resp.2D) problems. Moreover, the parameters ∊ in Eq.(20), α in Eq.(21), and τ_i in Eq.(22) are always taken as 10⁻⁸, 3, and 1, respectively. All calculations are carried out on the Linux environment of a personal computer of Lenovo (Intel(R) Core(TM) i5 CPU 3.2GHZ 4GB RAM).

3.1 1D case

Example 1 (Smooth subcritical flow) This problem has been studied^[7] to check the dissipative and dispersive errors in the kinetic schemes. The bottom shape of the river is

and the boundary conditions at x=±10 are specified as h=1 and hu=1.

Figure 2 shows the numerical steady-state solutions obtained by the NMGM on the mesh of 512 uniform cells, in comparison with the exact solutions obtained by solving the algebraic system

Figure 3 displays the numerical residuals obtained by the NMGM versus the NMGM iteration number N_step and CPU time for three meshes of 512, 1 024, and 2 048 uniform cells, respectively. The results show that the NMGM is very efficient and fast to get the correct steadystate solutions. Moreover, the convergence behaviors are similar on those different meshes, and the NMGM iteration number does not increase with refining the mesh. In those calculations, the HLL flux is used, and the grid level number in the V-cycle multigrid is set to be 4, i.e., N_L=3.

Table 1 investigates the effects of the HLL, LLF, and Roe fluxes shown in Remark 2 on the convergence behavior of the NMGM, in comparison with the BLU-SGS iteration in solving Eq. (21), where N, N_L, N_step, and T_cpu denote the cell number, the coarse mesh level number of the geometric multigrid, the total Newton iteration number, and the CPU time T_cpu, respectively, ρ denotes the spectral radius of the iteration matrix of the block SGS method around the steady-state solution, and R_∞ :=− ln ρ corresponds to the asymptotic convergence rate.

The results show that the NMGM is a little more efficient than the BLU-SGS iteration for the HLL and Roe fluxes, and exhibits great advantages for the LLF flux, specially, in the case of the W-cycle multigrid with N_L=5. The W-cycle multigrid is also tested for the HLL and Roe fluxes, and its performance is almost the same as the V-cycle multigrid.

The BLU-SGS iteration with the HLL and Roe fluxes are more efficient than the LLF flux. For the former, N_step increases very slowly with refining the mesh, but the latter does nearly linearly increase in terms of the cell number N. Such phenomenon can be explained by comparing the spectral radius ρ of the iteration matrix, whose values are around 0.5 or larger than 0.9.

Figure 4(a) displays the distribution (in the complex plane) of the eigenvalues of the block SGS iteration matrix on the mesh with 256 uniform cells for the HLL, LLF, and Roe fluxes. We see that the eigenvalues are almost around 0 for the HLL flux (except for two eigenvalues with relatively large imaginary parts of±0.4, respectively) and the Roe flux (except for only one eigenvalue located around −0.5), while those are widely distributed for the LLF flux (some of them are near 1). Thus, the “high frequency” eigenvalues for the LLF flux are much more than those for the HLL and Roe fluxes. Figure 4(b) plots an asymptotic relation of the spectral radius ρ of the iteration matrix with respect to the spatial step size Δx for the LLF flux as follows:

(27)

where C is about 0.092.

Example 2 (Two transcritical flows over a bump) The example simulates two transcritical flows over a bump, which have been widely used to test the SWE solvers^{[3, 7, 29]}. The bed topography is

(28)

for a channel of length 25.

(Ⅰ) Transcritical flow without a shock

The discharge hu=1.53 is imposed at the upstream boundary condition x=0, while the water height h=0.66 is imposed at the downstream end of the channel x=25 when the flow is sub-critical.

Figure 5 shows the numerical steady-state solution h + z obtained by the HLL flux and the V-cycle multigrid with N_L=3 on the mesh of 512 uniform cells in comparison with the exact solution in Ref. ^[30]. Figure 6 gives the convergence history of the NMGM in terms of the NMGM iteration number N_step and CPU time on three meshes of 512, 1 024, and 2 048 uniform Fig. 5 Example 2(I): steady-state solution h + z obtained by NMGM on mesh of 512 uniform cells cells, respectively. The results show that the steady-state solutions can be correctly and fast obtained by the NMGM, the convergence behaviors are similar to Example 1, and the NMGM iteration number does not increase with the mesh refinement.

Similar to Table 1, Table 2 investigates the effects of the HLL, LLF, and Roe fluxes on the convergence behavior of the NMGM, in comparison with the BLU-SGS iteration. The convergence behaviors of the NMGM and the BLU-SGS iteration are almost the same as those in Example 1.

The distribution (in the complex plane) of the eigenvalues of the block SGS iteration matrix in Fig. 7(a) shows that the “high frequency” eigenvalues for the LLF flux are much more than those for the HLL or Roe flux. Figure 7(b) shows that the spectral radius ρ of the iteration matrix has an asymptotic relation (27) with C~0.3 for the LLF flux in terms of the spatial step size Δx.

(Ⅱ) Transcritical flow with a shock

The discharge hu is taken as 0.18 on the upstream boundary, and h=0.33 is specified on the downstream boundary condition. In this case, the Froude number Fr= increases to a value larger than 1 above the bump, and then decreases to less than 1.

Figure 8 shows the numerical steady-state solution on the mesh of 512 uniform cells in comparison with the exact solution^[30], where a stationary shock appears and is well resolved.

Figure 9 gives the convergence history of the NMGM versus the NMGM iteration number N_step and CPU time on three meshes of 512, 1 024, and 2 048 uniform cells, respectively. The results show that the NMGM is very efficient and fast for achieving the correct steady-state solution with a shock, the convergence behaviors are almost similar on the three uniform meshes, and the iteration number N_step scarcely changes with the mesh refinement. The HLL flux and V-cycle multigrid with N_L=3 have been adopted in those calculations.

The effects of the HLL, LLF, and Roe fluxes on the convergence behavior of the NMGM and the BLU-SGS iteration are investigated in Table 3, showing that the NMGM performs much more efficiently than the BLU-SGS iteration, and it becomes obvious as N_L increases in the multigrid and the W-cycle multigrid is adopted with the LLF flux. The convergence behaviors depend on the numerical flux, and the block LU-SGS method with the HLL or Roe flux is more efficient than the LLF flux. Such observation is consistent with the previous study. However, different from the results of Examples 1 and 2(Ⅰ), the results in Table 3 show that N_step increases almost linearly with the mesh refinement for the BLU-SGS iteration with the HLL, LLF, and Roe fluxes, and the spectral radius ρ of the iteration matrix increases asymptotically to 1 as N increases.

The distribution in the complex plane of the eigenvalues of the block SGS iteration matrix in Fig. 10(a) still shows that the “high frequency” eigenvalues for the LLF flux are much more than those for the HLL or Roe flux. Figures 10(b), 10(c), and 10(d) show the asymptotic relation (27) of ρ with respect to the spatial step size Δx, where C is about 2.58, 0.15, and 4.53 for the HLL, LLF, and Roe fluxes, respectively.

Example 3 (Wet-dry boundary problem) The bed topography for this problem is the same as that in Example 2, but for a channel of length 20. Initially, the flow with the water height h⁽⁰⁾(x)=0.22 − z(x) is static in the channel, i.e., u(x)=0. A steady state with the dry bed h(x)=max(0.1 − z(x), 0) and u(x)=0 will be reached if imposing h=0.1 at the interval ends x=0 and 20. Since the steady solution involves wet/dry transition, the Jacobian matrix of these numerical fluxes near the wet-dry boundary becomes singular so that it is much more difficult and challenging.

Figure 11 shows the steady solutions h + z and the error in the rest water surface obtained on the mesh of 512 uniform cells, where the LLF flux and W-cycle multigrid with N_L=3 are used. It is seen that the correct steady solutions are efficiently obtained by the NMGM, and the rest water surface is preserved exactly up to the machine precision.

Figure 12 gives the convergence history of the NMGM in terms of the NMGM iteration number N_step and CPU time on four meshes of 512, 1 024, 2 048, and 4 096 uniform cells, respectively. The convergence behaviors of the NMGM and Nstep do not depend on the uniform mesh number N.

Because the NMGM and the BLU-SGS iteration with the HLL or Roe flux fail to work now, Table 4 only investigates the convergence behaviors of the NMGM and the BLU-SGS iteration with the LLF flux. We see that the NMGM is much more efficient than the BLU-SGS iteration and it is obvious as N_L increases and the W-cycle multigrid is adopted. In this problem, the big solution error is mainly introduced around the wet-dry boundary, and can be fast reduced by using the W-cycle multigrid with properly increasing the coarse level number N_L.

3.2 2D case

Unless specifically stated, in the following, the HLLC flux and V-cycle multigrid are adopted, and the coarse mesh level number N_L is set to be 3.

Example 4 The first 2D example is to simulate the fluid flows in two symmetric channels with a flat bottom constricted from both sides in the y-direction^[8]. Specifically, Channel I is with a constriction angle α=5◦ started at x=10, and Channel Ⅱ is with a constriction angle α=15° started at x=10, ended at x=30, and then followed by a straight narrower channel (see Fig. 13 for their geometries and corresponding structured mesh of 72×40 cells). The inflow conditions h=1 and v=0 and the Froude number Fr==2.5 are imposed at x=0, the outflow boundary condition is specified at x=90, while the slip boundary conditions are employed on the top and bottom boundaries.

Figure 14 displays the contour plots of the steady solutions for two channels obtained by the NMGM on the mesh of 576 × 320 cells. For the steady state flow in Channel Ⅰ, two bore waves starting at the points (10, ±20) interact at the point (45, 0), and then two regular reflections happen around x=74 due to the constriction. The present result is well comparable to the numerical result given by the adaptive moving mesh method to solve the time-dependent SWEs^[8] and the analytical result in Ref. ^[31], especially, the numerical values of the water heights of the first and second plateaus are 1.25 and 1.527 1, respectively, which agree well with those results in Refs. [8] and [31]. For the steady flows in Channel Ⅱ, two shock waves started at (10, ±20) meet around (33, 0) interact with two expansion waves generated at the points (30, ±15), and then two regular reflections happen at (50, ±15). Finally, two reflected shock waves meet around (69, 0). From the contour plots in Fig. 14, we see that the NMGM can well capture those waves and their interactions with high resolution and non-oscillation.

Figures 15 and 16 show the convergence history of the NMGM versus the NMGM iteration number N_step and CPU time on three successively refined meshes for Channel Ⅰ and Channel Ⅱ, respectively. It is seen that the NMGM works successfully on those meshes, and the convergence behaviors of the NMGM are independent of the cell number N. Compared with Channel Ⅰ, the computation of the steady solution for Channel Ⅱ seems more difficult for the NMGM, and takes more iteration steps. In those calculations, ∊ _p=2 × 10².

Example 5 This example is to solve another channel constriction problem^[32], in which a channel of length 3 units is with a symmetric constriction of length 1 unit at its center x=1.5, and the variable width

where W_min is the minimum channel width (see Fig. 17(a) for its geometry and corresponding structured mesh of 96 × 32 cells). The boundary conditions are specified as follows:

The supercritical outflow boundary condition is at x=3, and the slip boundary conditions are on the top and bottom boundaries.

Our calculations, where h₀=1, W_min=0.9, and Fr_in=0.5, 0.67, 1.2, 1.7, or 2, investigate the effect of the Froude number Fr_in on the convergence of the NMGM, and demonstrate the robustness of the NMGM. Figures 17(b)-17(f) display the contours of the steady water depth h obtained by the NMGM on the mesh of 384 × 128 cells, where the subcritical region (in which the Froude number Fr(x) < 1) is shaded. Those plots show that the steady solutions resolved by the NMGM are non-oscillatory, and the symmetry of the flow is well preserved. The results are very similar to those in Ref.[32], and the following five different types of steady flow patterns can be observed:

Case 1 (Fr_in=0.5) The flow is smooth and purely subcritical.

Case 2 (Fr_in=0.67) The flow is subcritical at the inflow and outflow, critical at the channel throat, and with a steady discontinuity in the divergent region of the channel.

Case 3 (Fr_in=1.2) The flow shows the cross-wave pattern with the oblique downstream jumps, and is subcritical at the inflow, critical at the throat, and supercritical at the outflow.

Case 4 (Fr_in=1.7) The flow is supercritical at the inflow and outflow, and with oblique jumps and a subcritical pocket within the constriction.

Case 5 (Fr_in=2) The flow is purely supercritical everywhere throughout the channel.

Figures 18(a)-18(e) present the convergence history of the NMGM in terms of the NMGM iteration number N_step on three meshes. It is seen that the NMGM exhibits good robustness and works well on those meshes for the above steady flows. Case 2 requires more iteration steps than the other cases. Figure 18(f) gives the histogram of the NMGM iteration number N_step of the NMGM in terms of the Froude number Fr_in on the mesh of 384 × 128 cells. We further see that the NMGM requires more iteration steps when the Froude number Fr_in is less than 0.3 or equal to 0.7.

Example 6 The test describes a transcritical flow over a 2D bump, which is a 2D extension of Example 2. The bottom topography is defined by

in the channel of [0, 25] × [−5, 5]. The discharge hu=1.53, hv=0 is imposed at x=0, the water height h=0.52 is imposed at the downstream x=25, and the reflective boundaries are specified at y=±5 in the y-direction.

The contours of the steady water surface h + z obtained by the NMGM are displayed in Fig. 19 on the mesh of 640 × 320 uniform cells, where the subcritical region is shaded. In the steady flow, the cross-wave pattern with the oblique downstream jumps is observed. Moreover, the flow varies from the subcritical region at the inflow to the supercritical region at the outflow, similar to Example 2(Ⅰ). The efficiency of the NMGM is demonstrated in Fig. 20 by displaying its convergence history versus iterations and CPU time on three meshes of 160 × 80 uniform cells, 320 × 160 uniform cells, and 640 × 320 uniform cells, respectively. In those calculations, ∊ _p=2.

Example 7 The last example simulates a steady channel flow around a hill (a relatively high bump). It is defined by

in the channel of [0, 25] × [−5, 5] with the following boundary conditions:

and the reflective boundary conditions at y=±5 in the y-direction. Initially, the channel is full of static flow with the water height h(x, y)=max(0.2 − z(x, y), 0). Since the steady solution involves wet/dry transition, it becomes much more difficult to be obtained. Therefore, the steady solution is obtained by solving the time-dependent SWEs with numerical schemes. The NMGM and the BLU-SGS iteration with the HLLC or Roe flux fail to work.

Figure 21 displays the contours of the steady water surface h + z obtained by the NMGM with the LLF flux on the mesh of 512 × 256 uniform cells, in comparison with the result given by using a first-order accurate explicit finite volume scheme with the LLF flux to solve the time-dependent SWEs (1) on the same mesh. The latter takes 412 652 time steps and 13 804 s CPU time to reduce the residual to 1.5 × 10⁻¹² at the physical time t=1 701.96. The results show that the steady solution given by the NMGM agrees well with that given by the explicit scheme, and the NMGM with the LLF flux works efficiently and robustly for the steady solution with wet/dry transition.

Figure 22 gives the convergence history of the NMGM with the LLF flux in terms of the NMGM iteration number N_step and CPU time on three meshes of 128 × 64, 256 × 128, and 512 × 256 uniform cells, respectively. It is seen that the convergence behaviors of the NMGM are similar to those meshes, and N_step scarcely changes with the mesh refinement.

4 Conclusions

A Newton multigrid method is developed for 1D and 2D SWEs with topography and dry areas. The steady-state SWEs are first approximated by the well-balanced finite volume discretization based on the hydrostatic reconstruction technique. The resulting nonlinear system is linearized by the Newton method, and the geometric multigrid method with the block SGS smoother is used to solve the linear system. The proposed Newton multigrid method makes use of the local residual to regularize the Jacobian matrix of the Newton iteration, and can handle the steady-state problem with wet/dry transition. Several numerical experiments are conducted to demonstrate the efficiency, robustness, and well-balanced properties of the proposed method. Moreover, the relation between the convergence behavior of the proposed method with the HLL, LLF, or Roe flux, and the distribution of the eigenvalues of the iteration matrix is detailedly discussed, in comparison with the block LU-SGS method. The numerical results show that the convergence behaviors of the NMGM depend on the numerical flux and the Froude number of the flow.