Novel Approach for solving the discrete Stokes problems based on Augmented Lagrangian and Global Techniques: Application to Saddle-Point Linear Systems from Incompressible flow
A. Badahmane A. Ratnani H. Sadok






1]A. Badahmane
1]A. Ratnani
2]H. Sadok
[1]The UM6P Vanguard Center, Mohammed VI Polytechnic University, Benguerir 43150, Lot 660, Hay Moulay Rachid, Morocco.
[2]LMPA, Université du Littoral Côte d'Opale, 50 Rue F. Buisson, BP 699 - 62228 Calais cedex, France.














































=msbm10


Z-5.5muZ
exoExercice
ideaIdea

prePreuve
proPropriété
exeExample
theoremTheorem[section]
propositionProposition
definitionDefinition[section]
remarkRemark[section]
lemLemma[section]


§ ABSTRACT 
In this paper, a novel   augmented Lagrangian preconditioner based on global Arnoldi for  accelerating the convergence of Krylov subspace methods applied to linear systems of equations with a block three-by-three structure, these systems typically arise from   discretizing  the Stokes equations using mixed-finite element methods. In practice, the components of velocity are always approximated using a single finite element space. More precisely, in two dimensions, our new approach based on  standard space of scalar finite element basis functions to discretize the velocity space. This componentwise splitting can be shown to induce  a natural block  three-by-three structure.
Spectral analyses is established for the exact versions of these preconditioners. Finally, the obtained numerical results claim that our  novel approach is more efficient and robust  for solving the  discrete Stokes problems. The efficiency of our new approach is evaluated by measuring computational time.


Stokes equation, saddle point problem, Krylov subspace method, global Krylov subspace method, augmented Lagrangian-based preconditioning.





§ INTRODUCTION

 The Stokes problem is discretized using conforming finite element spaces X^h⊂ Q_2 and Q^h_1⊂ Q_1 that satisfy the inf-sup condition for the Stokes velocity and pressure such as  Taylor–Hood elements <cit.>.
The discrete form of the weak formulation can be cast as a block linear system of the form:

    𝒜_3× 3𝐮=[     A     O B_x^T;     O     A B_y^T;   B_x   B_y     O ][ u_x; u_y;   p ]
    =
    [ f_x; f_y;   g ]_b,

assuming that  n_u=2n and n_p are  respectively,  dimensions of velocity solution and  pressure finite-dimensional spaces with (n_u+n_p=N). 
Where  A∈ℝ^n× n is the scalar-Laplacian matrix, it is worth nothing  that A  is  symmetric positive definite (SPD) matrix, the n_p × n matrices B_x and B_y represent weak derivatives in the x and y directions,
f_x, f_y and  g are given vectors. Besides, we assume, as is typically the case in most applications of Stokes problem, n_u>>n_p. The increasing popularity of mixed finite element methods for Stokes and Navier-Stokes flows has been a significant cause of saddle-point systems, such as the one in (<ref>). A major source of applications for saddle-point problems, can be found in <cit.>. 
In general, inasmuch as the large dimension and sparsity of the matrices A and  B,
it is sensible for systems (<ref>) to be solved by iterative methods. Additionally, since the
coefficient matrix A is nonsingular, numerous effective methods have been put forward
with great popularity. Such as the successive overrelaxation (SOR)-like methods <cit.>, variants of the Uzawa-type methods <cit.>,
Hermitian and skew-Hermitian (HSS) method, which was initially  introduced by Bai, Golub, and Ng in <cit.>. Additionally, the  PHSS iteration method has been presented in <cit.>. For a more in-depth understanding of the works related to the stationary iterative
methods, please refer to references <cit.>.
Generally speaking, iteration methods become more attractive than direct methods from
two aspects of storage requirements and computing time. In order to solve the linear system (<ref>) in an efficient manner, we often use valid preconditioning techniques to accelerate
Krylov subspace methods, such as GMRES  method <cit.>. As is well known, a clustered spectrum of preconditioned matrix often results in rapid rate of convergence for Krylov
subspace methods. Therefore, to achieve rapid convergence rate and improve computational
efficiency, a large number of efficient iteration methods and preconditioning techniques have
been presented in recent years, such as block triangular preconditioner applied  to the augmented linear system <cit.>, augmented Lagrangian-based preconditioning technique for a class of block three-by-three linear systems <cit.>, and so forth. We make organizations of this paper as follows. An example of modelling that leads to this
type of system is outlined in Section 1.
Section 2 introduces the 3×3 strategy. In Section 3, we recall and define the 2× 2 strategy.
Some numerical tests are implemented to show the eﬀectiveness of  the proposed  preconditioners, in particular in the presence of inexact solvers. At the end, we conclude
with a brief summary in Section 5.


 §.§ The Problem Setting
 
The Stokes equation describes the flow of a viscous fluid and is used in various fields, including aerodynamics, propulsion, and biomedical fluid analysis. In many cases, finding an exact solution to the Stokes equation can be challenging, so we often use numerical methods to approximate the solution <cit.>. Their discretization results in a linear system, as shown in Eq. (<ref>). In the incompressible case, the Stokes equation can be written as follows :

    {[ -∇⃗^2u⃗ +∇⃗ p             =     0⃗ in  Ω,;         ∇⃗·u⃗             =      0  in Ω. ].

The variable u⃗ is the unknown velocity field, the scalar function  p is the unknown pressure field. It is important to acknowledge that the Laplacian  and divergence operators are defined in <cit.>.
The first equation in  Eq. (<ref>) represents conservation of the momentum of the fluid (and so is the momentum equation), and the second equation enforces conservation of mass. We consider the problem posed on a domain Ω of dimension d=2 with boundary conditions ∂Ω=∂Ω_D∪∂Ω_N defined by

    [                                  u⃗                                   = w⃗ on ∂Ω_D,∂u⃗/∂ n-n⃗p= s⃗ on ∂Ω_N, ]

where :

    
  * w⃗: is the vorticity variable, given by:

    w⃗   =   ∇⃗×u⃗,

where × is the curl operator,


  * s⃗: function depends on  the outflow boundary to ensure that mass is conserved,  

  * n⃗: the outward-pointing normal to the boundary,

  * ∂u⃗/∂ n: denotes
the directional derivative in the normal direction.

In practice, the d components of velocity are always approximated using
a single finite element space <cit.>, then the discrete formulation of Eq. 
 (<ref>) can be expressed as a two-by-two partitioning of the discrete Stokes system, which the matrix of the system is a saddle point matrix  defined as follows :

    𝒜_2×2x=([   A B^T;   B   0 ])([ u; p ])=([ f; g ])_h,

where  A∈ℝ^n_u× n_u is the vector-Laplacian matrix, it is worth nothing  that A  is  symmetric positive definite (SPD) matrix, B∈ℝ^n_p× n_u is divergence matrix  with rank(B^T)=n_p, f∈ℝ^n_u and g∈ℝ^n_p are given vectors.



§ MOTIVATION:

The main motivation of this work, instead of using a single finite element space to discretize the velocity space and to obtain the two-by-two partitioning (<ref>), we use a standard space of scalar finite element basis functions {ϕ_j}_j=1^n, we set n_u = 2n and define the velocity basis set

    {ϕ⃗_1, …, ϕ⃗_2n} := {(ϕ_1, 0)^T, …, (ϕ_n, 0)^T, (0, ϕ_1)^T, …, (0, ϕ_n)^T}. 5.38

This component-wise splitting can be shown to induce  a natural block three-by-three partitioning of the discrete Stokes system (<ref>), for more details, we refer to <cit.>.
Specifically, with

    u:= ([u_x]_1, …, [u_x]_n, [u_y]_1, …, [u_y]_n),

(<ref>) can be rewritten as :

    [     A     O B_x^T;     O     A B_y^T;   B_x   B_y     O ][ u_x; u_y;   p ]
    =
    [ f_x; f_y;   g ],

where the n × n matrix A is the scalar Laplacian matrix (discussed in detail in <cit.>), and the n_p × n matrices B_x and B_y represent weak derivatives in the x and y directions, 
where
 
    A = [a_ij], a_ij=∫_Ω∇ϕ_i·∇ϕ_j,

 
    B_x = [b_x,kj], b_x,ki=-∫_Ωψ_k∂ϕ_i/∂ x,

 
    B_y = [b_y,kj], b_y,kj=-∫_Ωψ_k∂ϕ_j/∂ y.

 
 
 


§ MATHEMATICAL BACKGROUND:

Given a square matrix A, the set of all eigenvalues (spectrum) of A is denoted by σ(A). When the spectrum of A is real, we use λ_min(A) and λ_max(A) to respectively denote its minimum and maximum eigenvalues. When A is symmetric positive (semi)definite, we write A ≻ 0 (A ≽ 0). In addition, for two given matrices A and B, the relation A ≻ B (A ≽ B) means A - B ≻ 0 (A - B ≽ 0). Finally, for vectors x, y, and z of dimensions n, m, and p, (x; y; z) will denote a column vector of dimension n+m+p.
In this paper, I will denote the identity matrix, specifying its size as appropriate to the context.



§ 3×3  STRATEGY FOR SOLVING THREE-BY-THREE LINEAR SYSTEM (<REF>)

The 3×3 strategy, based on the motivation outlined in Section 1, is designed to solve three-by-three saddle-point problem (<ref>).
 The 3×3  strategy can significantly reduce the computational cost compared  using 2× 2  strategy for solving the classical structure of saddle-point problem (<ref>). The preconditioning technique helps to improve the convergence rate of the Krylov subspace methods. 
This strategy is motivated by the use of a set of standard scalar finite element basis functions within a defined space, aimed at obtaining the three-by-three partitions of the saddle-point matrix (<ref>).




 §.§ Novel Augmented Lagrangian-based preconditioning and global
techniques:

Krylov subspace methods (such as GMRES) in conjunction with suitable preconditioners are frequently the method of choice for computing approximate solutions of such linear systems of equations.
First, problem (<ref>)  is reformulated as the equivalent augmented system 𝒜̅_3× 3𝐮̅= 𝐛̅, where

    𝒜̅_3× 3= 
    [ A + γ B^T_x Q^-1B_x                   0               B_x^T;                   0 A + γ B^T_y Q^-1B_y               B^T_y;                 B_x                 B_y                   0 ],

and 𝐛̅ = (f_x+B^T_x Q^-1g;f_y + γ B^T_y Q^-1g; g), with Q being an arbitrary SPD matrix and γ > 0 a user-defined parameter. Evidently, the linear system of equations 
𝒜̅_3× 3𝐮̅= 𝐛̅ is equivalent to 𝒜_3× 3𝐮 = 𝐛. This approach is inspired by the effectiveness of employing grad-div stabilization and augmented Lagrangian techniques to solve saddle-point problems.




  §.§.§  Preconditioning:

In this section, we investigate a new augmented Lagrangian-based preconditioning and global approach for solving (<ref>).  Left preconditioning of  (<ref>) gives the following new linear system:

    𝒫^-1𝒜̅_3× 3𝐮̅ = 𝒫^-1𝐛̅,

where  𝒫 is one of the preconditioners  below:

    
  * 𝒫_γ, α, x : is the  augmented Lagrangian preconditioner in the x direction.
 
  * 𝒫_γ, α, y : is the  augmented Lagrangian preconditioner in the y direction.

 The following two constraint-type preconditioners were proposed for accelerating the convergence of Krylov subspace methods, given as follows:

    𝒫_γ, α,x = 
    [ A + γ B^T_x Q^-1B_x                   0               B_x^T;                   0 A + γ B^T_x Q^-1B_x   (1 - γα^-1) B^T_y;                   0                   0              -α^-1Q ],



    𝒫_γ, α,y = 
    [ A + γ B^T_y Q^-1B_y                   0               B_x^T;                   0 A + γ B^T_y Q^-1B_y   (1 - γα^-1) B^T_y;                   0                   0              -α^-1Q ],

where α and γ are prescribed positive parameters.


  §.§.§ Algorithmic implementation  of the augmented Lagrangian preconditioners  𝒫_γ, α,x and 𝒫_γ, α,y.

In this part, we display the algorithmic implementation of 𝒫_γ, α,x and 𝒫_γ, α,y,
in which, inside Krylov subspace methods, the SPD subsystems were solved inexactly by the preconditioned conjugate gradient (PCG) method using loose tolerances. More precisely, the inner PCG solver for linear systems with coefficient matrix A, A + γ B^T_xQ^-1 B_x and A + γ B^T_yQ^-1 B_y was terminated when the relative residual norm was below 10^-6, with the maximum number of 100 iterations was reached. The preconditioner for PCG is incomplete Cholesky factorizations constructed using the  function  where opts.type = 'ict' with drop tolerance 10^-2. 
In the following parts, we will work on some specific problems. Every step of the Krylov subspace
method such as GMRES method is used in combination with the augmented Lagrangian  preconditioner to solve the  saddle-point
problem (<ref>).
We summarize the implementation of preconditioners 
𝒫_γ, α,x and 𝒫_γ, α,y in
Algorithms 1 and 2.
For the linear systems corresponding to A+ γ B^T_x Q^-1 B_x and A+ γ B^T_y Q^-1 B_y, we distinguish between two approaches:

    
  * Approach I. Since A+ γ B^T_x Q^-1 B_x 
    is  SPD matrix, we solve the linear systems corresponding to this  matrix independently 
by the preconditioned conjugate gradient method (PCG), the matrix is  formed  inside PCG  with incomplete Cholesky preconditioning, ichol(A).

  
 To implement the preconditioners 𝒫_γ, α,x and 𝒫_γ, α,y, we
use the following algorithms:
 


  * The subsystems corresponding to (A+ γ B^T_x Q^-1 B_x) are solved by PCG method. Within the PCG process, we perform sequence of matrix-vector product, first multiplying vectors by B_x, Q^-1
and then by B^T_x. 

 
We use  the steps described  in Algorithm 1 to implement Algorithm 2.

    
  * Approach II.
    In step 2 and 3 of Algorithms 1 and 2,
the secondary objective of this work is not to solve it independently, but instead to utilize  𝒫GCG method <cit.> for  solving  linear system with several right-hand sides of the following  form :

    (A+ γ B^T_y Q^-1 B_y)𝒳=ℋ,

where: 𝒳 and ℋ  are both  an n×2 matrices.  Each column of matrix 𝒳 is denoted as 𝒳^(1)=x  and 
𝒳^(2)=y, each column of matrix ℋ is denoted as  ℋ^(1)=r_1-B_x^Tz and ℋ^(2)=r_2-(1-γα^-1)B_y^Tz,
𝒳_0 is the initial guess of  solution (<ref>) and R_0=ℋ-(A+ γ B^T_y Q^-1 B_y)𝒳_0 is the initial residual.

By leveraging the structure of the augmented based-Lagrangian preconditioners 𝒫_γ, α, x,
𝒫_γ, α, y and the approach II, in the rest of the paper, we refer
to the new preconditioners as  𝒫_γ, α,x, G and 𝒫_γ, α, y, G where

    
  * 𝒫_γ, α,G,x: denotes the global augmented Lagrangian preconditioner in the x direction.
 
  * 𝒫_γ, α,G,y: denotes  global augmented Lagrangian preconditioner in the y direction.

To implement the preconditioners 𝒫_γ, α,x,G and 𝒫_γ, α,y,G, we
use the following algorithms:
 


  * The subsystem with several right-hand sides corresponding to (A+ γ B^T_x Q^-1 B_x) is solved by 𝒫GCG method. Within the 𝒫GCG process, we perform sequence of matrix-vector product, first multiplying vectors by B_x, Q^-1
and then by B^T_x. 



 
We  apply a similar approach as in Algorithm 3 to implement Algorithm 4.


§ 2×2  STRATEGY FOR SOLVING TWO-BY-TWO LINEAR SYSTEM (<REF>)

In this strategy, we employ a single finite element space to discretize the velocity field and achieve the two-by-two partitioning (<ref>).


 §.§ Novel Augmented Lagrangian-based preconditioning and global
techniques:

The iterative solution of the discrete Stokes equations has attracted considerable attention in recent years. Here we limit ourselves to discussing solution algorithms based on preconditioned Krylov subspace methods. In [13], the following two constraint-type preconditioners were proposed for accelerating the convergence of Krylov subspace methods.
First, problem (<ref>)  is reformulated as the equivalent augmented system 𝒜̅𝐮 = 𝐛̅, where

    𝒜̅_2× 2= 
    [ A + γ B^T Q^-1B             B^T;               B               0 ],

and 𝐛̅ = (f+B^T Q^-1g; g), with Q being an arbitrary SPD matrix and γ > 0 a user-defined parameter. Evidently, the linear system of equations 
𝒜̅_2× 2𝐮 = 𝐛̅ is equivalent to 𝒜𝐮 = 𝐛 . This approach is motivated by the success of the use of grad-div stabilization and augmented Lagrangian techniques for solving saddle point problems.


  §.§.§  Preconditioning:

In this section, we investigate a new augmented Lagrangian-based preconditioning and global approach for solving (<ref>). The idea of preconditioning is to transform the linear system (<ref>) into another one that is easier to solve. Left preconditioning of  (<ref>) gives the following new linear system:

    𝒫^-1_γ,α𝒜̅_2× 2𝐮 = 𝒫^-1𝐛̅,

where  𝒫_γ,α is given as follows

    𝒫_γ, α= 
    [ A + γ B^T Q^-1B    (1-γα^-1)B^T;               0           α^-1Q ],

To apply the preconditioner, we need to solve systems of the following form:
 


  * We compute y,

  * The matrix  A+ γ B^T Q^-1 B is SPD, we solve it iteratively by PCG method. We address this by employing the PCG method for iterative solution. Within the PCG process, we
perform sequence of matrix-vector product, first multiplying vectors  by B, by Q^-1  and then by B^T.



 §.§ Spectral analysis

The distribution of eigenvalues and eigenvectors of a preconditioned matrix has a significant connection to how quickly Krylov subspace methods converge. Hence, it's valuable to analyze the spectral characteristics of the preconditioned matrix, denoted as 𝒫_γ,α^-1𝒜̅_2× 2. In the upcoming theorem, we  will estimate the 
 lower and upper bounds for the eigenvalues of preconditioned
matrix 𝒫_γ,α^-1𝒜̅_2× 2. 


Let the preconditioner 𝒫_γ,α be defined as in (<ref>). Then the eigenvalues of 𝒫_γ,α^-1𝒜̅_2× 2 are all real, positive and bounded. Furthermore
 the matrix 𝒫_γ,α^-1𝒜̅_2× 2 is diagonalizable and has n_p+1 distinct eigenvalues 
 {1,λ_1,...,λ_n_p}. 


Assume that λ represents an eigenvalue of the preconditioned matrix  and u̅=(u;p) is the associated eigenvector. In order to deduce the distribution of eigenvalues, we analyze the following generalized eigenvalue problem
	
    𝒜̅_2× 2u̅ =λ𝒫_γ,αu̅.

	(<ref>) can be reformulated  as follows
 
    {[ (1-λ)(A+γ B^TQ^-1B)u+(1+λ(γα^-1-1))B^Tp                                     =0,;                                      Bu                             =-λα^-1 Qp. ].

In the case where λ=1, equation (<ref>) is always true for u∈Null(B), consequently, there exist n_u-n_p linearly independent eigenvectors ([ u^(i);
 0 ]), i=1,..,n_u-n_p, corresponding to the eigenvalue 1, where u^(i)∈Null(B).
     If λ=1 and u=0, from the second equation of (<ref>), it can be deduced that p=0. This conflicts with the initial assumption that the column vector (
     u;
p) is an eigenvector of the preconditioned matrix 𝒫_γ,α^-1𝒜̅_2× 2. If λ≠ 1 and p=0, from the first equation of (15), it can be deduced that  u must be 0. This contradicts
the initial assumption that (u; p) is the eigenvector of the preconditioned matrix and therefore u≠ 0 and
p≠ 0. 
Since λ≠ 1, from (<ref>) we further obtain :

    p =-α/λQ^-1 B u.

Substituting p from the above relation  in the first equation of (<ref>), we get :

    λ^2(A+γ B^T Q^-1 B)u-λ(A+α B^T Q^-1 B )u +α B^T Q^-1 Bu  = 0.

Premultiplying (<ref>) with u^T/u^Tu 
(<ref>) gives:

    (a+γ q)λ^2-(a+α q )λ +α q=0,

which can be written 

    λ^2-bλ+c=0,

where a, q, b and c are given  as follows:

    a = u^T A u/u^Tu,
    q=u^T B^T Q^-1 B u/u^Tu, b=a+α q/a+γ qand c=α q/a+γ q.

As a result, it is immediate to see that the roots of  (<ref>)  are real and positive, given by
n_p eigenvalues 
λ_1 = b - √(b^2 - 4c)/2 and n_p eigenvalues 
λ_2 = b + √(b^2 - 4c)/2 of the preconditioned matrix.
After some manipulations, λ_1 and λ_2 must hold the
following inequalities:

    λ_1  ≥2λ_min(B^T Q^-1 B) /λ_max(A)+(1+α-γ)λ_max(B^T Q^-1 B),λ_2  ≤2αλ_max(B^T Q^-1 B) /λ_min(A)+(α-γ)λ_min(B^T Q^-1 B).







§ NUMERICAL RESULTS

In this section, we report on the performance of inexact variants of the proposed block preconditioner using a test problem taken from <cit.>, which corresponds to a 2D Stokes flow problem. The programs are performed on a computer with an Intel Core i7-10750H CPU @ 2.60 GHz processor and 16.0 GB RAM using MATLAB R2020b. In Tables 1 and 2, we report the total required number of outer GMRES iterations and elapsed CPU time (in seconds) under “Iter” and "CPU", respectively. The total number of inner GMRES (PCG) iterations to solve subsystems with coefficient matrices (A+ γ B^T_x Q^-1 B_x)  and (A+ γ B^T_x Q^-1 B_x)   are reported under Iter (Iter_pcg). No restart is used for either  GMRES iteration. The initial guess is taken to be the zero vector and the iterations are stopped as soon as

    𝒜x_k - b_2 ≤ 10^-7b_2,

where x_k is the computed k-th approximate solution. In the tables, we also include the relative error and relative residual

    Err :=  x_k - x^* _2/ x^* _2,

and 

    Res := 
    𝒜x_k - b_2/b_2
    ,

where x^* and x_k are respectively, the exact solution and its approximation obtained in the k-th iterate. In addition, we have used right-hand sides corresponding to random solution vectors.


L-shaped two dimensional domain Ω_∟, parabolic inflow boundary condition, natural outflow boundary condition. Consider the  Stokes  equation system (<ref>) posed in  Ω_∟=( -1,5)×( -1,1). In this scenario, we have a situation where there is a slow flow occurring in a rectangular duct with a sudden expansion. This configuration is often referred to as "flow over a backward facing step". Dirichlet no-flow (zero velocity) boundary conditions on uniform streamline imposed in the inflow boundary (x=-1; 0≤ y ≤ 1), the Neumann condition (<ref>) is again applied at the outflow boundary  (x = 5;-1<y< 1).
We use Q_2-P_1  mixed finite element approximation from  IFISS library <cit.> to discretize this problem in  Ω_∟, where:

    
  * Q_2: biquadratic finite element approximation on rectangles for the velocity,
    
  * P_1: triangular finite element approximation on triangle for the pressure,

the nodal positions of this mixed finite element  are illustrated in the following Fig. <ref>:

Then we obtain the  nonsingular saddle point problem (<ref>).

The numerical results of strategy 3×3 with approaches I and II for the tested example are listed in Tables 1 and 3.  
In Tables 2 and  4, we list numerical results with respect to Iter, CPU and Res in the case of 2×2  and 3×3 strategies.


§ IN THE CASE Γ=1E-04 AND Α=1E+01.


 It can be seen numerically that the Approach II incorporated with  𝒫_γ, α,x,G and 𝒫_γ, α,y,G
preconditioners 
is  more convenient than the Approach I  
incorporated with  𝒫_γ, α a
in terms of
both iteration number and CPU time.

Table 2 reports the corresponding results of the two strategies with the proposed preconditioners, which show that the  3×3 strategy with 𝒫_γ, α,x,G perform much better than the 2×2 strategy with 𝒫_γ, α, especially for the large
problems. Numerical results are reported in Tables 3 for the tested methods
with respect to the number of outer iteration steps, inner iteration steps and elapsed CPU time in seconds, denoted
as "Iter", "Iter_pcg" and "CPU", respectively.

The Approach II incorporated with  𝒫_γ, α, y, G outperforms the Approach I with  𝒫_γ, α, y on efficiency and performance
concerning both iteration steps and CPU times. Moreover, the Approach II incorporated with  𝒫_γ, α, y
preconditioner is more economical and it is superior to the other two preconditioners regarding
execution time, especially for relatively large size problems.

In the case γ=1e-02 and α=1e+01.








It was observed in all the Tables  that 2× 2 and 3× 3  strategies with the inexact augmented Lagrangian-based preconditioner  exhibits faster convergence for smaller  values of γ. However, for large γ the total timings increase due to the fact that the condition number of the blocks (A+ γ B^T_y Q^-1 B_y) and (A+ γ B^T_x Q^-1 B_x) goes up as increase. The 3× 3  strategy  incorporated with 𝒫_γ, α,x,G and 𝒫_γ, α,y,G  preconditioners
 demonstrates significantly better performance. 
This superiority is observed across various comparisons with 2× 2  strategy incorporated with 𝒫_γ, α. Moreover, 3× 3  strategy consistently requires less CPU time for convergence. Therefore, it can be concluded that the convergence behavior of 3× 3 strategy with 𝒫_γ, α,x,G and  𝒫_γ, α,y,G outperforms that of other methods.

To discretize problem (<ref>)  using  Taylor-Hood Q_2-Q_1  mixed-finite element approximation in  Ω_∟, we utilize the nodal positions of Q_2-Q_1  from  IFISS library <cit.>, where:

    
  * Q_1: denotes a quadratic finite element approximation on rectangle,

and  the nodal positions of Q_2-Q_1  are given below in the following Fig. <ref>:

Then we derive the nonsingular saddle point problem (<ref>). 

To further confirm the effectiveness of the 3× 3 strategy incorporated with  𝒫_γ, α, x, G or  𝒫_γ, α, y, G preconditioners, numerical results of the 2× 2 and  3× 3 strategies incorporated with  various preconditioners, 
with respect to Iter, Iter_pcg, CPU, Res and  Err for saddle point problems with different values of l,
are reported in the following Tables. 


§ IN THE CASE Γ=1E-04 AND Α=1E+01.









It can be observed from Tables 10 and  11, that the  𝒫_γ, α, x, GGMRES and 𝒫_γ, α, y, GGMRES 
 methods  have great advantage in the CPU compared with 𝒫_γ, α, xGMRES and 𝒫_γ, α, yGMRES  methods, which shows with Approach II  the total timings are much smaller than in the case  of Approach I.  Although the results in Tables 2 and 4 indicates 
  applying 3× 3 strategy  and Algorithms 3 and 4 to solve the problem with several right-hand sides (A+ γ B^T_x Q^-1 B_x)𝒳=ℋ or  (A+ γ B^T_y Q^-1 B_y)𝒳=ℋ,  need less computing time than using 2× 2 strategy with Algorithms 1 and 2. 
 

§ IN THE CASE Γ=1E-02 AND Α=1E+01.











By comparing the results in Tables 13, 14, 15 and 16 it can be seen that our proposed strategy  incorporated with the 
preconditioned 𝒫_γ, α, x, G and 𝒫_γ, α, y, G
GMRES methods succeed in producing high-quality approximate solutions in
all cases, while the
 3× 3 strategy incorporated with 
preconditioned 𝒫_γ, α,x,GGMRES, 𝒫_γ, α,y,GGMRES 
methods, outperforms the classical 2× 2 strategy incorporated with 
preconditioned 𝒫_γ, αGMRES method, 
in terms of Iter and CPU times. Besides, numerical results in Tables  above show that the 3× 3 strategy incorporated with 
preconditioned 𝒫_γ, α,x,GGMRES and 
𝒫_γ, α,y,GGMRES 
methods with
proper α and γ is still very efficient  even for larger  size of problems.


§ CONCLUSION

In this paper, we  introduce a new class of augmented Lagrangian-preconditioners  based on global conjugate gradient (GCG) method for solving three-by-three linear systems, focusing on systems arising from finite element discretizations of the Stokes flow problem. Numerical experiments on a challenging 2D model problem indicate that the corresponding inexact preconditioner with 3×3 strategy can achieve significantly faster convergence compared to previous versions of the augmented Lagrangian-based preconditioner.
Future work will concentrate on replacing the incomplete Cholesky inner preconditioners with multilevel preconditioners to enhance the scalability of the global conjugate gradient.
elsarticle-num-names

9
Arno 
W.E. Arnoldi, The principle of minimized iterations in the solution of the matrix eigenvalue problem, Quart. Appl. Math., 9 (1951), pp. 17-29.
badahmane
A. Badahmane, A.H. Bentbib and H. Sadok, Preconditioned  global  Krylov subspace methods for solving saddle point problems with multiple right-hand sides, Electron. Trans. Numer. Anal., 51 (2019), pp. 495-511.
AHSS
Z.-Z. Bai, G.H. Golub, Accelerated Hermitian and skew-Hermitian splitting iteration methods for saddle-point problems, IMA J. Numer. Anal., 27 (2007), pp. 1-23.
PHSS
Z.-Z. Bai, G.H. Golub and J.-Y. Pan, Preconditioned Hermitian and skew-Hermitian splitting methods for non-Hermitian positive semidefinite linear systems, Numer. Math., 98 (2004), pp. 1-32.
HSS
Z.-Z. Bai, G.H. Golub and M.K. Ng, Hermitian and skew-Hermitian splitting methods for non-Hermitian
positive definite linear systems, SIAM J. Matrix Anal. Appl., 24 (2003), pp. 603-626.


Wang
Z.-Z. Bai, Z.-Q. Wang, On parameterized inexact Uzawa methods for generalized saddle point problems, Linear Algebra Appl., 428 (2008), pp. 2900–2932.
Parlett
Z.-Z. Bai, B.N. Parlett,  Z.-Q.Wang, On generalized successive overrelaxation methods for
augmented linear systems. Numer. Math., 102 (2005), pp. 1–38.
Benzi
M. Benzi, G. H. Golub and J. Liesen, Numerical solution of saddle point problems, Acta
Numerica, 14 (2005), pp. 1-137.
Bramble
 J.H. Bramble, J.E. Pasciak, A.T. Vassilev, Analysis of the inexact Uzawa algorithm for saddle point
problems. SIAM J. Numer Anal., 34 (1997), pp. 1072–1092
Vuik
 G. Ebadi, N. Alipour, and C. Vuik, Deflated and augmented global Krylov subspace methods for the
matrix equations, Appl. Numer. Math., 99 (2016), pp. 137–150
benzi2024
F. P. A. Beik, M. Benzi, An augmented Lagrangian-based preconditioning technique for a class of block three-by-three linear systems, Applied Mathematics Letters., 149 (2024).
beik2022
F. P. A. Beik, M. Benzi,
Preconditioning techniques for the coupled Stokes-Darcy problem: Spectral and field-of-values analysis, Numerische Mathematik, 150 (2022), pp. 257-298.
Benzi2
M. Benzi, G. H. Golub and J. Liesen, Numerical solution of saddle point problems, Acta
Numerica, 14 (2005), pp. 1-137.


Elman
H.C. Elman, D.J. Silvester and  A.J. Wathen, Finite Elements and Fast Iterative Solvers with Applications in Incompressible Fluid Dynamics, Oxford University Press,  New York., (2005).


Golub2
G.H. Golub, X. Wu, J.-Y. Yuan, SOR-like methods for augmented systems, BIT Numer. Math., (2001), pp. 71-85.


Guo
P. Guo, C.-X. Li, S.-L. Wu, A modified SOR-like method for the augmented systems, J. Comput.
Appl. Math., 274 (2015), pp, 58–69.
GBICG
K. Jbilou, H. Sadok and A. Tinzefte, Oblique projection methods for linear systems with multiple right-hand sides, Electron. Trans. Numer. Anal., 20 (2005), pp. 119-138.
Ichol:Manteuffel
T.A. Manteuffel, An incomplete factorization technique for positive definite linear systems, Math. Comput., 34 (1980), pp. 473-497.
Wang
N.-N. Wang, J.-C. Li, A class of new extended shift-splitting preconditioners for saddle
point problems, J. Comput. Appl. Math., 357 (2019), pp. 123-145.
saad:2003
Y. Saad, Iterative Methods for Sparse Linear Systems, SIAM, Philadelphia, PA., (2003).
Saad
Y. Saad, M.H. Schultz, A generalized minimal residual algorithm for solving nonsymmetric linear systems, SIAM J. Sci. Stat. Comp., 7 (1986), pp. 856-869.
saad1993flexible
Y. Saad, A flexible inner-outer preconditioned GMRES algorithm, SIAM Journal on Scientific Computing, 14 (1993), pp. 461–469.
Zhang
J.-J. Zhang, J.-J. Shang, A class of Uzawa-SOR methods for saddle point problems, Appl. Math
Comput., 216  (2010), pp. 2163–2168.


