Numerical solution algorithms for the Navier Stokes equations is the rapidly developing field, in which various cost-effective methods are being proposed. This can be domain decomposition methods (Rønquist, 1996; Girault and Wheeler, 2008), methods based on splitting (Henriksen and Holmen, 2002; Viguerie and Veneziani, 2018), asymptotic expansions (Panasenko, 1998; Hoanga and Martinez, 2018), multigrid methods (Griebel et al., 1998; Pernice and Tocci, 2001) etc. In this paper we use a simple Schur complement algorithm. A survey of the Schur complement methods can be found in Loghin and Wathen (2002).

The paper focuses on the generalized Navier–Stokes system in three space dimensions (3D) approximated by the mixed finite element method using the MINI element called also the P1-bubble/P1 pair (Arnold et al., 1984). The main contribution of the paper consists in the development of the vectorized codes for fast assembly finite element matrices in the MATLAB language so that the loop over tetrahedra is avoided. These codes are very fast and enable to perform experiments with relatively large scale problems. Vectorized codes were proposed for different differential operators, e.g. see Rahman and Valdman (2015) and references therein. Our research is inspired by Koko (2019), where the vectorized codes for the generalized Stokes problem are proposed. However, an extension to the Navier–Stokes problem is not trivial. Formally, it consists of adding the nonlinear convective term in the momentum equation. This nonlinearity is typically treated iteratively using the Oseen or the Newton type linearization (Elman et al., 2014). In the first case, the Oseen convection matrix depends on the nodal velocity field from the previous iteration that is presented in computations. The situation for the Newton convective matrix is more involved, since it depends, in addition, on the partial derivatives of the velocity field that are not immediately presented. We approximate them from appropriate directional derivatives so that computed approximations are invariant with respect to local renumbering of nodes and can by easily vectorized. Another ingredient in the assembling process is the bubble component elimination that is performed by the Schur complement reduction on the element level. This elimination requires inverting blocks of the local matrices that is done by the vectorized Cramer rule. The elimination itself uses a vectorized variant of linear combinations of vectors and of sums of vector outer products. Note that the vectorized codes for the Navier–Stokes problem play an important role in the whole solution process, since the finite element matrices for the linearized subproblems are repeatedly assembled in each iterative step that have to be fast.

Numerical experiments with our codes show superconvergence rate of the finite element approximation of the pressure component that is close to O ( h 3 / 2 ) or higher. Similar results were observed for the pure Stokes problem by Cioncolini and Boffi (2022).

The rest of the paper is organized as follows. Section 2 presents the classical and weak formulation of the problem and introduces the basic iterative schemes. In Section 3, we describe mixed finite element approximation based on the MINI element, for which we derive local linear systems that are split on the bubble and non-bubble components. We present also an idea how to approximate partial derivatives from the discrete vector field of the previous iteration. In Section 4, we discuss element matrices in more details so that their final forms may be coded by vectorized operations. Section 5 introduces ideas of the vectorization and refers to our free available codes. In Section 6, we describe the dual implementation of the iterative scheme that we use in computations. Section 7 is devoted to numerical experiments. First, we demonstrate low time requirements of the vectorized codes and than we compute the convercence rates. Finally, we conclude with some remarks and comments in Section 8.

To better understand our presentation we use different font styles. The mathematics bold symbols are used for the vector functions or their vector arguments, e.g. u, p. The “text up” symbols are used for matrices, vectors and scalars on the level of finite elements, e.g. A k k , B l , p. The “bold up” symbols are used for global matrices and vectors, e.g. A, B U , p. Finally, the typewriter style is used for codes, e.g. for p = 1:4.

Let Ω ⊂ R 3 be a bounded domain with a sufficiently smooth boundary Γ = ∂ Ω. We consider the steady Navier–Stokes system with the homogeneous Dirichlet boundary condition: (1) − ν Δ u + u · ∇ u + α u + ∇ p = f in Ω , ∇ · u = 0 in Ω , u = u D on Γ , where ν > 0 is the kinematic viscosity, α > 0 is a constant from a time discretization of the unsteady problem, and f : Ω → R 3 represents the external forces. We are searching for the vector velocity field u : Ω → R 3 and the scalar pressure field p : Ω → R. The existence and uniqueness of a weak solution to (1) are discussed in Girault and Raviart (1986).

We will consider two iterative methods for solving (1) with different linearizations of the convection term u · ∇ u. Let w be an approximation of u and δ w be an increment such that u = w + δ w. The Oseen type linearization is based on the omission of the linear term with δ w: u · ∇ u = ( w + δ w ) · ∇ u = w · ∇ u + δ w · ∇ u ≈ w · ∇ u . In the Newton type linearization we omit the quadratic term with δ w: u · ∇ u = w · ∇ u + δ w · ∇ ( w + δ w ) = w · ∇ u + ( u − w ) · ∇ w + δ w · ∇ δ w ≈ w · ∇ u + u · ∇ w − w · ∇ w . We arrive at the following linearizations of (1): (2) − ν Δ u + w · ∇ u + ρ u · ∇ w + α u + ∇ p = f + ρ w · ∇ w in Ω , ∇ · u = 0 in Ω , u = 0 on Γ , where w : Ω → R 3 is given. For ρ = 0 / 1, we get the Oseen/Newton linearization of (1).

The weak formulation of (2) requires the following spaces: V = ( H 0 1 ( Ω ) ) 3 , Q = { q ∈ L 2 ( Ω ) : ∫ Ω q d x = 0 } and the forms: a ( u , v ) = ν ∑ i = 1 3 ∫ Ω ∇ u i · ∇ v i d x , c ( w ; u , v ) = ∫ Ω ( w · ∇ u ) · v d x , m ( u , v ) = α ∫ Ω u · v d x , b ( v , q ) = − ∫ Ω q ( ∇ · v ) d x , a ρ ( w ; u , v ) = a ( u , v ) + c ( w ; u , v ) + ρ c ( u ; w , v ) + m ( u , v ) , l ρ ( w ; v ) = ∫ Ω f · v d x + ρ c ( w ; w , v ) , where u , v , w ∈ ( H 1 ( Ω ) ) 3 , q ∈ L 2 ( Ω ), and u = ( u 1 , u 2 , u 3 ), v = ( v 1 , v 2 , v 3 ).

The weak formulations of (2) read as follows: (3) Find ( u , p ) ∈ V × Q such that for all ( v , q ) ∈ V × Q , a ρ ( w ; u , v ) + b ( v , p ) = l ρ ( w ; v ) , b ( u , q ) = 0 . The problem (1) can be solved by the following iterative scheme: (4) Given ( u ( 0 ) , p ( 0 ) ) ∈ V × Q . For κ ⩾ 1 , find ( u ( κ ) , p ( κ ) ) ∈ V × Q solving (3) with w = u ( κ − 1 ) . This scheme will be called the Oseen/Newton iteration for ρ = 0 / 1, respectively. A comprehensive discussion of convergence results can be found in Elman et al. (2014), pp. 344–346. The Oseen iteration are also caled the Picard iteration.

We approximate (3) by the mixed finite element method. It requires to choose a finite element pair satisfying the inf-sup condition (Brezzi and Fortin, 1991). Here, we use the P1-bubble/P1 pair called also the MINI element proposed by Arnold et al. (1984).

We suppose that Ω is a polyhedral domain. Let T h be a regular partition of Ω ‾ given by tetrahedra T j ∈ T h , 1 ⩽ j ⩽ n t , Let x i ∈ Ω ‾, 1 ⩽ i ⩽ n p , be finite element nodes. Each tetrahedron T = T j has 4 vertices x i 1 , x i 2 , x i 3 , x i 4 ∈ T. Each vertex of T is associated with the local linear basis function ϕ j = ϕ j ( x ) such that ϕ j ( x i k ) = δ j k , 1 ⩽ j , k ⩽ 4. The bubble function is defined on T by the product: ϕ b = ϕ b ( x ) = 4 4 ϕ 1 ( x ) ϕ 2 ( x ) ϕ 3 ( x ) ϕ 4 ( x ), x ∈ T. We denote also ϕ 5 = ϕ b . Integrals over T can be evaluated by the formula: ∫ T ϕ 1 α 1 ( x ) ϕ 2 α 2 ( x ) ϕ 3 α 3 ( x ) ϕ 4 α 4 ( x ) d x = 6 | T | α 1 ! α 2 ! α 3 ! α 4 ! ( α 1 + α 2 + α 3 + α 4 ) ! , where | T | is the tetrahedron volume.

On T h we introduce the space of the bubble functions and the piecewise linear functions: B h = { v h ∈ C ( Ω ¯ ) : v h | T = c T ϕ b , c T ∈ R ∀ T ∈ T h } , W h = { v h ∈ C ( Ω ¯ ) : v h | T ∈ P 1 ( T ) ∀ T ∈ T h } , respectively, where P 1 ( T ) is the space of all polynomials on T of the degree at least one. We approximate V and Q as follows: V h = { v h ∈ ( W h ⊕ B h ) 3 : v h ( x i ) = 0 ∀ x i ∈ Γ } , Q h = { q h ∈ W h : ∫ Ω q h d x = 0 } , respectively. The finite element approximation of (3) reads as follows: (5) Find ( u h , p h ) ∈ V h × Q h such that for all ( v h , p h ) ∈ V 0 h × Q h , a ρ ( w h ; u h , v h ) + b ( v h , p h ) = l ρ ( w h ; v h ) , b ( u h , q h ) = 0 , where w h ∈ V h is an approximation of w from (3).

On T ∈ T h the components of u h = ( u h 1 , u h 2 , u h 3 ) are the linear combinations of the basis functions: u h k ( x ) = ∑ j = 1 4 u k j ϕ j ( x ) + u k b ϕ b ( x ) , 1 ⩽ k ⩽ 3 , and p h ( x ) = ∑ j = 1 4 p j ϕ j ( x ) . The linear systems arising from (5) over one element T read as: (6) A ¯ 11 A ¯ 12 A ¯ 13 B ¯ 1 ⊤ A ¯ 21 A ¯ 22 A ¯ 23 B ¯ 2 ⊤ A ¯ 31 A ¯ 32 A ¯ 33 B ¯ 3 ⊤ B ¯ 1 B ¯ 2 B ¯ 3 0 u ¯ 1 u ¯ 2 u ¯ 3 p = f ¯ 1 f ¯ 2 f ¯ 3 0 , where u ¯ k = ( u k 1 , u k 2 , u k 3 , u k 4 , u k b ) ⊤ , 1 ⩽ k ⩽ 3, and p = ( p 1 , p 2 , p 3 , p 4 ) ⊤ . The matrices A ¯ k l ∈ R 5 × 5 are algebraic counterparts of the form a ρ : A ¯ k k = R ¯ + C ¯ + ρ N ¯ k k + M ¯ , A ¯ k l = ρ N ¯ k l for k ≠ l , where R ¯ is the diffusion matrix, C ¯ is the Oseen convection matrix, N ¯ k l are the Newton convection matrices, 1 ⩽ k , l ⩽ 3, and, M ¯ is the mass matrix. Further, B ¯ k ∈ R 4 × 5 are the components of the divergence matrix representing the form b and the right-hand sides f ¯ k ∈ R 5 , 1 ⩽ k ⩽ 3, correspond to the form l ρ . The entries are given as follows: (7) ( R ¯ ) i j = ν ∫ T ∇ ϕ j ( x ) · ∇ ϕ i ( x ) d x , ( C ¯ ) i j = ∑ k = 1 3 ∫ T w h k ( x ) ∂ k ϕ j ( x ) ϕ i ( x ) d x , ( N ¯ k l ) i j = ∫ T ∂ l w h k ( x ) ϕ j ( x ) ϕ i ( x ) d x , ( M ¯ ) i j = α ∫ T ϕ j ( x ) ϕ i ( x ) d x , ( f ¯ k ) i = ∫ T f k ( x ) ϕ i ( x ) d x + ρ ∑ l = 1 3 ∫ T w h l ( x ) ∂ l w h k ( x ) ϕ i ( x ) d x for 1 ⩽ i , j ⩽ 5, where w h k are the components of w h , and (8) ( B ¯ k ) i j = − ∫ T ∂ k ϕ j ( x ) ϕ i ( x ) d x for 1 ⩽ i ⩽ 4 and 1 ⩽ j ⩽ 5.

Recall that w h represents the velocity field from the previous iteration of discretized analogy of (4). In the discrete case we know only nodal values of w h at the vertices x i of T, 1 ⩽ i ⩽ 4. In ( C ¯ ) i j , ( N ¯ k l ) i j , and ( f ¯ k ) i we will use constant approximations of w h k and ∂ l w h k on T. The approximation of w h k is defined by the arithmetic mean: w k = 1 4 ( w k 1 + w k 2 + w k 3 + w k 4 ) , 1 ⩽ k ⩽ 3 , where w k i = w h k ( x i ), 1 ⩽ i ⩽ 4. The approximation of ∂ l w h k denoted by δ l w k will be taken by their value at x 1 : δ l w k = ∂ l w h k ( x 1 ). The relations between the gradient ∇ w h k ( x 1 ) = ( δ 1 w k , δ 2 w k , δ 3 w k ) and the derivatives in the directions x i + 1 − x 1 lead to three linear systems of the form: (9) ∇ w h k ( x 1 ) · ( x i + 1 − x 1 ) = w h k ( x i + 1 ) − w h k ( x 1 ) , 1 ⩽ i ⩽ 3 , from which δ l w k , 1 ⩽ k , l ⩽ 3, can be computed. Note that the solutions are invariant with respect to (local) renumbering of the vertices x i of T. Using the same symbols for the respective approximations of ( C ¯ ) i j , ( N ¯ k l ) i j , and ( f ¯ k ) i , we get: (10) ( C ¯ ) i j = ∑ k = 1 3 w k ∫ T ∂ k ϕ j ( x ) ϕ i ( x ) d x , ( N ¯ k l ) i j = δ l w k ∫ T ϕ j ( x ) ϕ i ( x ) d x , (11) ( f ¯ k ) i = ∫ T f k ( x ) ϕ i ( x ) d x + ρ ∑ l = 1 3 w l δ l w k ∫ T ϕ i ( x ) d x . Hence, C ¯ is the linear combination of B ¯ k and N ¯ k l are multiples of M ¯ so that the Ossen and Newton matrices may be assembled from the matrices of the Stokes problem.

The element matrices and vectors will be divided on the non-bubble and bubble components: R ¯ = R r r ⊤ ω R , C ¯ = C c U c L ⊤ ω C , N ¯ k l = N k l n k l n k l ⊤ ω N k l , M ¯ = M ( α ) m ( α ) m ( α ) ⊤ ω M ( α ) , B ¯ k = ( B k , B k b ) , f ¯ k = f k f k b , u ¯ k = u k u k b , where R , M ( α ) , B k , C , N k l ∈ R 4 × 4 , r , m ( α ) , B k b , c U , c L , n k l , f k , u k ∈ R 4 , and ω R , ω M ( α ) , ω C , ω N k l , f k b , u k b ∈ R for 1 ⩽ k , l ⩽ 3.

The tetrahedron T ∈ T h is represented by four vertices x i = ( x i , y i , z i ), 1 ⩽ i ⩽ 4. We denote the entries of x i − x j by x [ i j ] = x i − x j , y [ i j ] = y i − y j , z [ i j ] = z i − z j and w k [ i j ] = w k i − w k j , i ≠ j.