The spectral implementation of the nonstationary Stokes problem with nonstandard boundary conditions

*Correspondence: houertani@ksu.edu.sa 2Information Technology Department, College of Computer and Information Sciences, King Saud University, Riyadh, Saudi Arabia Full list of author information is available at the end of the article Abstract This paper deals with the implementation of the spectral discretization of the vorticity–velocity–pressure formulation of the nonstationary Stokes problem. We use the implicit Euler scheme for the discretization in time. We propose a global method for the resolution of the linear system resulting from the discrete problem. This method is implemented, and some numerical results are presented which confirm the good convergence for the three unknowns.


Introduction
The low velocity flow of an incompressible viscous fluid can be modeled by the nonstationary Stokes equations. Such problem is typically provided with a Dirichlet boundary condition on the velocity [1][2][3]. In this paper, we deal with this problem in two-or threedimensional bounded domain using nonstandard boundary conditions. In dimension 2, the problem is supplied with boundary conditions on the normal component of the velocity and the Dirichlet condition on the vorticity (since it is a scalar). However, in dimension 3, it is supplied with boundary conditions on the normal component of the velocity and on the tangential components of the vorticity vector field. This permits us to prove that the above problem is equivalent to a time-dependent variational formulation with unknown vorticity, velocity, and pressure [4][5][6][7]. For the discretization of this problem, we combined the implicit Euler's scheme in time and the spectral method in space. One can see [8] for the analysis of this discretization in the case of the stationary Stokes problem.
In [9], the Generalized Minimal Residual method (GMRES) [10] is used to solve the stationary Stokes problem. However, this method is not easy to implement since the matrix is not symmetric and needs a huge number of iterations to converge. In [11], a direct method was proposed permitting to simplify and to improve the resolution.
For the nonstationary case, which is the subject of this paper, we will use the global resolution based on the good results established in [11]. The discretization in time by the implicit Euler method permits us to stabilize the discrete problem. The resulting global matrix is now symmetric and positive definite, making possible the use of the gradient conjugate method. As a result, we obtained high accuracy with an optimized number of iterations. Some numerical results are presented to show the time and space convergence.
The paper is organized as follows: • In Sect. 2 the continuous problem and some regularity results are presented. • Sect. 3 deals with the discrete problem and error estimates.
• The nonhomogeneous case is handled in Sect. 4.
• The implementation details of the discrete problem are developed in Sect. 5.
• Sect. 6 is dedicated to describe thoroughly the numerical results.

The continuous problem
Let Ω be a bounded simply connected domain of R d (d = 2 or d = 3), and Γ be its boundary that we suppose connected and Lipschitz; x = (x, y) in d = 2 or x = (x, y, z) in d = 3 is the generic point in Ω. We introduce [0, T] an interval of R where T is a positive constant and n is the unit outward vector to Ω on Γ . The nonstationary Stokes problem is presented as follows: where f is the data function and ν is the positive constant viscosity. The unknowns are the velocity u and pressure p. The boundary operator ζ is defined according to the dimension as follows: • In dimension 2 (d = 2), ζ represents the trace operator on the boundary Γ . • In dimension 3 (d = 3), if the vector field w = (w x , w y , w z ), then ζ (w) = curl w × n on the boundary Γ . Using the formula -w = curl(curl w) -∇(div w), and the vorticity τ = curl u, we prove that (1) is equivalent to the following system: We suppose that the initial velocity and vorticity satisfy div u 0 = 0 and τ (x, 0) = τ 0 = curl u 0 in Ω.
We need to introduce the following Sobolev spaces: associated with the following norm and seminorm: Also H m (Ω) = W m,2 (Ω) is a Hilbert space provided with the scalar product where (·, ·) is the L 2 (Ω) scalar product. We denote by L 2 0 (Ω) the space of functions in L 2 (Ω) which have a null integral on Ω, and D(Ω) the space of infinitely many times differentiable functions with a compact support in Ω. Then we introduce the space H(div, Ω) as 1 2 and H 0 (div, Ω) = ϕ ∈ H(div, Ω) : ϕ · n = 0 on ∂Ω .
In the same way, we introduce the space H(curl, Ω) by Remark 1 Notice that the spaces H(curl, Ω) and H 0 (curl, Ω) coincide with the spaces H 1 (Ω) and H 1 0 (Ω) in dimension d = 2.
Now we introduce the following spaces to deal with time issues. Let B be a separable Banach space. We define C m (0, T; B) as the set of C m -class time functions with values in B.
Then C m (0, T; B) is a Banach space with the norm where ∂ k t ϕ is the partial derivative of order k in time of function ϕ. Consider also the spaces .
We also introduce the kernel of the bilinear form c(·, ·; ·): Then (τ , u) is a solution of the following reduced problem: We refer to [

Theorem 1
For a function f ∈ L 2 (0, T; (H 0 (div, Ω)) ), and u 0 ∈ L 2 (Ω) d which satis- . This solution satisfies where c is a positive constant that depends on Ω and T.

Theorem 2
For a data function f belonging to L 2 (0, T; H r-1 (Ω) d ), the solution (τ , u, p) of problem (4) belongs to L 2 (0, T; H r (Ω) ) for any r > 0 such that: when Ω is a polygon with largest angle equal to α.

The time discrete problem
In this section, we start by the discretization in time of problem (4) using implicit Euler method. We propose a partition of the interval For any function f ∈ L 2 (0, T; (H 0 (div, Ω)) ) and u 0 ∈ H 0 (div, Ω) satisfying condition (3), we are considering the following scheme: where f k = f(·, t k ).
For any k, The functional L is linear and continuous on V . Thus if (τ k , u k , p k ) is the solution of problem (9) then (τ k , u k ) belongs to U and is the solution of the following reduced problem: Based on the property of the bilinear formâ(·, ·, ·), proved in [8, Sect. 2], we state the following proposition.

Proposition 2
Suppose that the data function f ∈ L 2 (0, t; H 0 (div, Ω) ) and that the initial velocity u 0 ∈ H 0 (div, Ω) satisfies condition (3). Then, for any k, where c is a positive constant independent of k.
Then, according to the inf-sup condition (6), we can state the following theorem (see [8,Sect. 2] for its proof).

The spectral discrete problem
Hereinafter, we assume that the domain Ω is a square in dimension d = 2 or a cube in dimension d = 3. We adopt the Nédélec's finite element method analog for the spectral discretization [16,Sect. 2].
To define the discrete spaces, we introduce P n,m (Ω) as the space of polynomials with degree ≤ n with respect to x and ≤ m with respect to y, and P n,m,p (Ω) as the space of polynomials with degree ≤ n with respect to x, ≤ m with respect to y, and ≤ p with respect to z. When n = m = p, these spaces are equal to P n (Ω).
Let N be an integer ≥ 2. The space D N which approximates H 0 (div, Ω) is defined by The space C N approximating H 0 (curl, Ω) depends on the dimension (see Remark 1) and is defined by Finally, we consider the space M N for the approximation of L 2 0 (Ω) given by If ξ 0 = -1 and ξ N = 1, we consider the nodes ξ i , 1 ≤ i ≤ N -1, and the set of N + 1 weights ρ i , 0 ≤ i ≤ N , of the Gauss-Lobatto quadrature formula. We denote P n (-1, 1) the space of restrictions to [-1, 1] of polynomials with degree ≤ n, and then We also recall a property from [17, (13.20)], which is useful in what follows: According to (16), we define the discrete product on P N (Ω), for two continuous functions u and v by Finally, let I N denote the Lagrange interpolation operator on the nodes (ξ i , ξ j ), 0 ≤ i, j ≤ N , in dimension d = 2, and on the nodes (ξ i , ξ j , ξ k ), 0 ≤ i, j, k ≤ N , in dimension d = 3, with values in P N (Ω). Henceforth, we suppose that the data f is continuous on Ω × [0, T]. The spectral discrete problem is constructed from the time semidiscrete problem (8), (9) using the Galerkin method combined with numerical integration. It reads as follows: where the bilinear forms a N (·, ·; ·), b N (·, ·), and c N (·, ·; ·) are defined by Given that using (16) combined with Cauchy-Schwarz inequality, we prove that the bilinear formŝ a N (·, ·; ·), b N (·, ·), and c N (·, ·; ·) are continuous on (C N × D N ) × D N , D N × M N , and (C N × D N ) × C N , respectively, with norms bounded independently of N . Moreover, as a consequence of the exactness property (15), the forms b(·, ·) and b N (·, ·) coincide on D N × M N . We remark also that the functional L N is linear and continuous on D N . Let be the kernel of the bilinear form b N (·, ·). It corresponds to the space of divergence-free polynomials in D N . We consider also the kernel of the bilinear form c N (·, ·; ·).
We consider now the following reduced discrete problem: The well-posedness of the reduced discrete problem (20) is based on the positivity condition and the inf-sup condition satisfied by the bilinear formâ N (·, ·; ·) proved in [8, Lemmas 1 and 2] (see [18,Theorem 2.1] and [19] for similar results in finite element methods). We can presently write the following proposition.
where c is a positive constant independent of N and k.
We recall the inf-sup condition on the bilinear form b N (·, ·). We refer to [20, Lemma 3.1] for details of its proof.

Lemma 1 There exists a positive constant β independent of N such that the bilinear form
Using above results, we obtain the following theorem.

The error estimate
For the estimation of the error between the solution (τ k , u k , p k ) of problem (8)-(9) and the solution (τ k N , u k N , p k N ) of problem (18), we need to define the following space for r ≥ 0: We notice that in dimension d = 2 this space coincides with H r+1 (Ω). We have the following error estimate (see [8,Sect. 5] for the proof).

Theorem 5
Suppose that function f ∈ L 2 (0, T; H σ (Ω) d ) for a real number σ > d 2 and u 0 ∈ H μ (Ω) d for a real number μ > d 2 . For any k, 1 ≤ k ≤ K , if the solution (τ k , u k , p k ) of problem (8)-(9) belongs to H r (curl, Ω) × H r (Ω) d × H r (Ω) for a real number r ≥ d -1, then the following error estimate holds between this solution and the solution (τ k N , u k N , p k N ) of problem (18): Estimate (21) is optimal for the three unknowns, while the estimate of the pressure is not obtained for most spectral discretizations of the Stokes problem.

Corollary 1
Assume that the data f belongs to L 2 (0, T; H σ (Ω) d ) for a real number σ > d 2 and u 0 ∈ H μ (Ω) d for a real number μ > d 2 . Then for any k, 1 ≤ k ≤ K , the following error estimate holds between the solution (τ k , u k , p k ) of problem (8)-(9) and the solution (τ k N , u k N , p k N ) of problem (18): where σ Ω is a real number ≥ 1 depending only on Ω.
We obtain the following result.
Proof For any k, is a solution of problem (28)-(29) if and only if (τ k , u k h , p k ) is a solution of problem (8)- (9). Theorem 3 permits us to deduce the well-posedness of problem (28)-(29), and the estimate (31) is derived by combining (11) and (30).
Since we handle the nonhomogeneous boundary condition on the velocity, we define the following new discrete space for the discrete velocity: Suppose that for any k, 1 ≤ k ≤ K , the function g k ∈ L 2 (Γ ) and g k N is the approximation of g k defined as follows: On each edge (d = 2) or face (d = 3) Γ r of Ω, 1 ≤ r ≤ 2d, g k N|Γ r is equal to the image of g k |Γ r by the orthogonal projection operator from L 2 (Γ r ) onto P N-1 (Γ r ). Thus we write the following spectral discrete problem: and Theorem 8 If for any k, 1 ≤ k ≤ K , the data function f k is continuous on Ω and g k is in L 2 (Γ ) satisfying (23), the problem (32)-(33) has a unique solution Proof Writing problem (32)-(33) as a square linear system, we deduce from Theorem 4 that the unique solution of this problem when the data f k and g k N are zero is (0, 0, 0). This yields the existence and uniqueness property.
We derive the following error estimate from those proved in the homogeneous case (see [8,Sect. 5]) using slight modifications.

Theorem 9
If the assumptions of Theorem 5 hold and for any k, 1 ≤ k ≤ K , the data g k satisfies condition (23) and is such that each g k |Γ r , 1 ≤ r ≤ 2d, belongs to H τ (Γ r ) for a nonnegative real number τ , then the following error estimate holds between the solution (τ k , u k , p k ) of problem (24)-(25) and the solution (τ k N , u k N , p k N ) of problem (29)-(30): Corollary 2 Assume that for any k, for a real number σ > d 2 and that condition (23) is satisfied. Then the following error estimate holds between the solution (τ k , u k , p k ) of problem (24)-(25) and the solution where the real number σ Ω is same as in Corollary 1.

The implementation of the discrete problem
In this section, we propose a global method for the resolution of the discrete problem (18). This method was used to solve the stationary Stokes problem for the same formulation (see [11]). This new algorithm enhanced the performance of the previous resolution of Stokes problem (2D and 3D bounded domain) (see [9]) and optimized the execution time. In the following, we start by presenting the linear system. Afterwards, we describe the resolution algorithm.
For the discrete problem (18), as a linear system, we have to choose a basis of the discrete spaces C N , D N , and M N .
The Lagrange polynomials in P N (-1, 1) linked with the nodes ξ j are denoted as ϕ j , 0 ≤ j ≤ N . We define where J * is the set {0, . . . , N} \ {j * } and j * is the integer part of N 2 . For any k, 1 ≤ k ≤ K , the unknowns τ k N , u k N = (u k Nx , u k Ny ), and a pseudopressurep k N admit the expansions, in dimension d = 2 for simplicity, The functionp k N vanishes in (-1, -1) but no longer belongs to L 2 0 (Ω); however, the real pressure p k N can easily be recovered in a postprocessing step, thanks to the formula For any k, 1 ≤ k ≤ K , we denote by Φ k , U k , and P k the vectors made of these coefficients. Their respective dimensions are equal to d(d-1) 2 N d-2 (N -1) 2 , dN d-1 (N -1), and N d -1. Hereinafter, by supposing the viscosity ν = 1, problem (18) is equivalent to the following square linear system: If U 0 = (U 0 1 , U 0 2 ), the components of the vectors U 0 1 and U 0 2 are respectively u 1 0 (ξ i , ξ j ) and where A k T and B k T denote the transposed matrices of A k and B k , respectively. Matrix A k : Matrix A k is written as For any k, 1 ≤ k ≤ K , curl(τ k N ) = (∂ y τ k N , -∂ x τ k N ), and then the coefficients of the matrices A k 1 and A k 2 are deduced respectively from the two terms (∂ y τ k N , u k Nx ) N and (∂ x τ k N , u k Ny ) N . Thus, The coefficients of the matrices B k 1 and B k 2 are found respectively from the terms (∂ x u k Nx , p k N ) N and (∂ y u k Ny , p k N ) N . Let then We also note that matrices B k 1 and B k 2 are not square, having N d -1 lines and dN d-1 (N -1) columns.
Matrix D k : Matrix D k is a diagonal matrix with d(d-1) 2 N d-2 (N -1) 2 lines and columns. Its coefficients are equal to Matrix I k : Matrix I k is written as For any k, 1 ≤ k ≤ K , the coefficients of the matrices I k 1 and I k 2 are respectively equal to Lastly, for any k, 1 ≤ k ≤ K , we denote The components of the two entries F k 1 and F k 2 are respectively (u k-1 Nx , ϕ r ϕ * s ) N + h k (f k 1 , ϕ r ϕ * s ) N , 1 ≤ r ≤ N -1; s ∈ I * and (u k-1 Ny , ϕ * r ϕ s ) N + h k (f k 2 , ϕ * r ϕ s ) N , 1 ≤ s ≤ N -1; r ∈ I * , where the data function is f = (f 1 , f 2 ).
We solve the linear system (38) using the gradient conjugate method preconditioned by LU incomplete factorization.

Two-dimensional experiments
In this section, at first, we focus on the time convergence. We consider the square Ω = ] -1, 1[ 2 . We look at a given solution obtained from the formulas u = curl ψ and τ = curl u in the two cases: Case (i). Regular functions ψ and p, defined by ψ(x, y) = t sin(πx) sin(πy), p(x, y) = e -t xy.
The velocity is a Gaussian which is null for t = 0. We consider the spectral discrete param-  Hereinafter, we consider h = 0.0001 and T = 1, Figure 2 presents the spectral error curves on the vorticity, velocity and the pressure (log(|error|) as a function of log(N) for N varying from 5 to 30).
In Fig. 2(a), the discrete solution is computed from (40). We obtain a very good (exponential) convergence due to spectral discretization.
In Fig. 2(b), the discrete solution is calculated from (41). We get lower convergence due to the irregularity of the solution.
We notice that the convergence of the pressure is of the same order as for the vorticity and velocity.
homogeneous boundary conditions g = 0 and N = 30. Figure 4 corresponds from top to bottom and left to right to the discrete vorticity, the two components of the discrete velocity, and the discrete pressure for the data f and u 0 defined in (42), with g given by g(x, -1) = -t 1x 2 3 2 , g(x, 1) = t 1x 2 3 2 , g(±1, y) = 0, and N = 30. This estimation depends only on the regularity of the solution. Our forthcoming work concerns the nonlinear nonstationary vorticity-velocity-pressure formulation of the Navier-Stokes problem.