Enhanced algorithms for solving the spectral discretization of the vorticity–velocity–pressure formulation of the Navier–Stokes problem

The objective of the article is to improve the algorithms for the resolution of the spectral discretization of the vorticity–velocity–pressure formulation of the Navier–Stokes problem in two and three domains. Two algorithms are proposed. The first one is based on the Uzawa method. In the second one we consider a modified global resolution. The two algorithms are implemented and their results are compared.


Introduction
The Navier-Stokes system models the flow of a fluid, for example, the movements of air in the atmosphere, ocean currents, the flow of water in a pipeline and many other fluidflow phenomena. The change of one of the parameters associated with the Navier-Stokes equations (the domain, the boundary conditions, the data, the variational formulation, the approximation method, . . . ) gives us a new research problem. In the reference article [1], the authors handle the Navier-Stokes equations with no standard boundary conditions on the velocity and the pressure in a convex domain.
In this work, is assumed to be a bounded simply connected domain of R d (d = 2, 3) and ∂ is its continuous Lipschitz boundary. We focus on the following Navier-Stokes system: where ϕ is the velocity, P the pressure, f is the density of forces, ϕ is the fluid viscosity and n is an outer unit vector normal to ∂ . System (1) can be written as (see [2,3]): where • μ = curl ϕ is the vorticity verifying ϕ · ∇ϕ = μ × ϕ + 1 2 grad |ϕ| 2 , • p = P + 1 2 |ϕ| 2 is the dynamic pressure, • η is an operator such that η(curl ϕ) is the normal (resp. the tangential) boundary component of curl ϕ when d = 2 (resp. d = 3). The above problem has been solved in several previous works using finite-element discretization, see [2][3][4][5][6], while the extension to the discretization by the spectral method has been handled in [7,8] for solving stationary and nonstationary Navier-Stokes equations. In a previous work [8], we studied (2) using the spectral discretization method. The resulting linear system has been solved using the GMRES iterative method (see [9]) since the matrix of the resulting linear system is not symmetric.
To improve the resolution of this problem, we propose in this work two new methods in order to reduce the number of iterations and execution time. The first method consists in applying the Uzawa algorithm (see [10,11]) after reducing the unknowns of the system by eliminating the vorticity. The second method consists of modifying the iterative algorithm used to linearize the resolution of the Navier-Stokes problem in order to obtain a linear system with a symmetric positive-definite matrix. Then, the conjugate gradient method is used.
The paper is organized as follows: • The continuous and discrete weak formulation are described in Sect. 2.
• Sect. 3 provides the details of the obtained linear system and its implementation using the two proposed algorithms.
• In Sect. 4, we present and compare the two algorithms for the resolution of the discrete problem based on some numerical tests.
Hereinafter, we suppose that the domain is a square in dimension two and a cube in dimension three. Using the Nédélec's finite-element spaces on cubic three-dimensional meshes (see [14,Sect. 2]), we construct our discrete spectral spaces as follows.
Let an integer N ≥ 2; we consider X N the discrete space of the velocity and M N the discrete space of the pressure In the following, we use the Gauss-Lobatto quadrature formula to construct the discrete problem where ξ j , 0 ≤ j ≤ N , are the roots of the polynomial (1x 2 )χ N , χ N is the Legendre polynomial of degree N defined on [-1, 1], and ρ j , 0 ≤ j ≤ N , is the associated set of the positive weights. The following inequality enables us to show that the continuous and discrete norms are equivalent (see [15]). The discrete scalar product is defined such that for the continuous functions ϕ and ψ on¯ , we have: Then, using the Galerkin method and the numerical integration based on the Gauss-Lobatto quadrature formula (5), we deduce for a continuous data f on the following discrete variational formulation: The bilinear forms a N (·, ·; ·), b N (·, ·), and c N (·, ·; ·) are continuous and are formulated as follows: However, the discrete trilinear form T N (·, ·; ·) is defined as: Using Brouwer's fixed-point theorem ( [10], Chap. IV, Cor. 1.1), problem (6) has a unique solution, since the discrete bilinear form b N (·, ·) coincides with the continuous one on X N × M N due to the exactness of the quadrature formula and verifies the inf-sup condition (see [8], Lem. 3.9) In the case of d = 2 the error estimate between the continuous and discrete problems is given by (see [7] for the proof ) where c is a positive constant independent of N , In the case d = 3, the error estimate is difficult to prove and remains an open problem.

Numerical algorithms
We present in this section two numerical algorithms to improve the numerical resolution of discrete problem (6). In the following we consider =] -1, 1[ 2 .

The iterative algorithm
To linearize the Navier-Stokes problem, we consider the following iterative algorithm.
• Step 1: Solve the following linear Stokes problem: • Step 2: Solve the following problem, by assuming that the i -1 iterative solution We stop the iterations when the following inequality is satisfied where a positive small real number.
Since the discrete velocity function does not have the same degree in the different directions, we consider j ∈ 0, . . . , N the fixed integer and we define the following polynomial of degree N -1 where J is the set {0, . . . , N} \ {j }. Then, we write the unknowns (μ i N , ϕ i N ; p i N ) in two dimensions as: We consider the real pressure since the pseudopressurep is not in L 2 0 ( ). We denote by (μ i N , ϕ i N ; p i N ) the components of the solution (W , U, P) on grid (ξ j , ξ k ) 0 ≤ j, k ≤ N . Problem (10) can be written as the following linear system: where • A T and B T are, respectively, the transposed matrix of A and B.
• The matrix A is defined as: Since curl(μ i N ) = (∂ y μ i N , -∂ x μ i N ), the coefficient of the matrices A 1 and A 2 correspond, respectively, to the terms In two dimensions, the matrices A 1 and A 2 are square, their dimension is (N -1) 2 . • The matrix B = [B 1 , B 2 ] where B 1 and B 2 are two matrices, the coefficients of which are deduced, respectively, from the terms (∂ x ϕ i Nx , p N ) N and (∂ y ϕ i Ny , p N ) N . We note that the matrices B 1 and B 2 are not square, their dimension is N 2 -1 lines and 2N(N -1) columns.
• M is a diagonal matrix, the coefficients of which are found from the values of (μ i N , μ i N ) N . Its dimension is N(N -1) 2 . • The matrices D and N are, respectively, written as D = (D 1 , D 2 ) and N 1 = (N 1 , N 2 ) and are, respectively, calculated from the terms T N (μ i-1 N , ϕ i N ; v N ) and T N (μ i-1 N , ϕ i-1 N ; v N ).

Algorithm 1
To solve the linear system (11), we propose to remove the vorticity from the first equation and replace it in the second one.
This modification permits us to use the Uzawa algorithm (see [10,11]) with only the velocity and the pressure as unknowns.
LetÃ = AM -1 A T + 2D. The second and third equations of (12) can be decoupled as and Since BÃ -1 B T is symmetric and positive defined, we use the conjugate gradient method to solve system (14) at each iteration. Then, we deduce the velocity U and vorticity W , respectively, from (13) and the first equation of (12). The Uzawa algorithm used is written as: Uzawa algorithm • Let P 0 = 0.
• Initialization step: • Iterations: n ≥ 0 From U n and P n :

Algorithm 2
We propose a second algorithm to solve problem (11). We begin first by inverting the matrix G given by and then find U, W and P.
Since G is symmetric and positive defined, we use the preconditioned gradient conjugate method.
In Fig. 1 we present a comparison of the convergence order of the two presented algorithms based on the error curves on the solution (μ, ϕ, p) when N varies between 5 and 20. Figure 1(a) and (b) correspond, respectively, to the resolution using Algorithms 1 and 2. This shows clearly that the obtained error is better when using Algorithm 2.
κ(Ã) and κ(G) are, respectively, the condition number of the matrixÃ and that of the matrix G. Tables 1 and 2 show the number of convergence iterations and the corresponding average error for the unknowns (vorticity, velocity and pressure). While comparing the two results, it is clear once again that the second method has a better accuracy and converges within a few iterations.    We conclude that: • Although Algorithm 1 based on the Uzawa method it is easier to implement and requires much less memory space, it converges after a large number of iterations, which is due to the fact that the conditioning κ(Ã) of the matrix A is high if we compare it with the conditioning κ(G) of the matrix G, see Table 3 (hundreds of iterations for an error of 10 -5 order, see Fig. 1 and Table 1). • The convergence order of Algorithm 2 is much better than the order produced by Algorithm 1 (22 iterations for a relative error approximatively equal to 10 -25 , see Fig. 1 and Table 2).

Solution using Algorithm 2
Consider the numerical resolution of the problem (6) with f = (0, x 2 y) using Algorithm 2. Figure 2 shows, from top to bottom and left to right, the obtained values of the vorticity μ N , the two components of the velocity (ϕ Nx , ϕ Ny ) and the pressure p N , for N = 30.

Figure 3
The vector field of the velocity for f = (0, 0) and g defined in (16) In Fig. 3 we present the vector field of the velocity corresponding to the data f = (0, 0), and a boundary condition g given by (see [8]), g(-1, y) = -1y 2 , g(1, y) = 1y 2 , g(x, ±1) = 0, (16) where N = 40. We note that the vorticity μ N and pressure p N are null, since we handle a Poiseuille linear flow.

Conclusion and future work
In this work, we show the efficiency of the global resolution compared to the Uzawa algorithm adapted to the resolution of the discrete problem issued from the spectral discretization of the vorticity, velocity and pressure formulation of the Navier-Stokes problem. We achieved a good convergence with the global resolution through the transformation of the matrix to a symmetric and positive defined one. As future work, we are looking at applying the global method to the discretization by the spectral element methods for handling more complex domains.