Ricardo Costa
email: ricardodpcosta@gmail.com

Stéphane Clain
email: clain@math.uminho.pt

Gaspar J Machado

A sixth-order finite volume scheme for the steady-state incompressible Stokes equations on staggered unstructured meshes

Keywords: Stokes equations, incompressible fluid, finite volume, high-order scheme, preconditioning

We propose a new sixth-order finite volume scheme to solve the bidimensional linear steady-state Stokes problem on staggered unstructured meshes and complex geometries. The method is based on several classes of polynomial reconstructions to accurately evaluate the diffusive fluxes, the pressure gradient, and the velocity divergence. The main difficulty is to handle the div-grad duality to avoid numerical locking and oscillations. A new preconditioning technique based on the construction of a pseudo-inverse matrix is also proposed to dramatically reduce the computational effort. Several numerical simulations are carried out to highlight the performance of the new method.

Introduction

The Stokes problem together with the Darcy system are two classical prototype models of mixed problems where the pressure function derives from the divergence-free velocity constraint. The finite volume method turns to be a natural framework to design built-in conservation schemes since the pioneer book of Patankar [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF] and we refer to the textbook of Ferziger and Peric [START_REF] Ferziger | Computational methods for fluids dynamics[END_REF] for an overview of the finite volume for the Navier-Stokes equations. A large range of schemes has been developed to provide both accurate and stable solutions where one can distinguish different kinds of approach namely staggered or collocated discretizations, structured grids or unstructured meshes for complex geometries, and coupled or segregated velocity and pressure leading to a saddle point problem or a projection method in the divergence-free space (see the introductions of [START_REF] Griffith | An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as preconditioner[END_REF][START_REF] Ramírez | A new higher-order finite volume method based on moving least squares for the resolution of the incompressible Navier-Stokes equations on unstructured grids[END_REF] for a short overview). Another fundamental challenge concerns the preconditioning of the linear system deriving from the space discretization.

Second-order methods are a standard in industry for the computation of incompressible flow and in commercial software development. There exists an important literature and books on the subject using the finite difference [START_REF] Chorin | Numerical method for solving incompressible viscous flow problems[END_REF][START_REF] Peyret | Computational methods for fluid flow[END_REF], the finite element [START_REF] Pironneau | Finite element methods for fluids[END_REF][START_REF] Temam | Navier-Stokes equations. Theory and numerical analysis[END_REF][START_REF] Brezzi | Mixed and hybrid finite element methods[END_REF][START_REF] Guermont | An overview of the projection methods for incompressible flows[END_REF][START_REF] Zienkiewicz | The finite element method for fluid dynamics[END_REF], or the finite volume approach [START_REF] Patankar | Numerical heat transfer and fluid flow[END_REF][START_REF] Gao | The finite volume projection method with hybrid unstructured triangular collocated grids for incompressible flows[END_REF][START_REF] Vidović | A superlinearly convergent finite volume method for the incompressible Navier-Stokes equations on staggered unstructured grids[END_REF][START_REF] Shang | Generalized formulations for the Rhie-Chow interpolation[END_REF][START_REF] Eymard | Convergence of a locally stabilized collocated finite volume scheme for incompressible flows[END_REF][START_REF] Boivin | A finite volume method to solve the Navier-Stokes equations for incompressible flows on unstructured meshes[END_REF]. However, there are fewer papers on higher order approximations (greater than the second-order of accuracy). Very high-order schemes for incompressible fluid flow have been developed using the finite difference framework using the Padé methodology (the so-called compact scheme) [START_REF] Hirsh | Higher order accurate difference solutions of fluid mechanics problems by a compact differencing technique[END_REF] on staggered structured grids (see [START_REF] Kampanis | A staggerid grid, high-order accurate method for the incompressible Navier-Stokes Equations[END_REF][START_REF] Griffith | An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as preconditioner[END_REF] and references herein) for the fourth-order case, and more recently for the sixth-order case [START_REF] Boersma | A 6th order staggered compact finite difference method for the incompressible Navier-Stokes and scalar transport equations[END_REF]. Finite element [START_REF] Guermont | An overview of the projection methods for incompressible flows[END_REF][START_REF] Frochte | A splitting technique of higher order for the Navier-Stokes equations[END_REF] and discontinuous Galerkin methods [START_REF] Montlaur | Discontinuous Galerkin methods forthe Stokes equations using divergence-free approximations[END_REF][START_REF] Ferrer | A high order discontinuous Galerkin finite element solver for the incompressible Navier-Stokes equations[END_REF][START_REF] Nigro | Up to sixth-order accurate A-stable implicit schemes applied to the discontinuous Galerkin discretized Navier-Stokes equations[END_REF] also received important contributions to achieve very high-order approximation both in time and space.

Accurate finite volume approximations are receiving considerable attention to compute approximations for the Euler system and the compressible Navier-Stokes equations in two-and three-dimensional geometries [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Michalak | Unstructured high-order accurate finite-volume solutions of the Navier-Stokes equations[END_REF][START_REF] Ivan | High-order solution-adaptative central essentially non-oscillatory (CENO) method for viscous flows[END_REF][START_REF] Diot | The MOOD method in the three-dimensional case: very-high-order finite volume method for hyperbolic systems[END_REF] but the incompressible case is far from being so well-developed. Pereira et al. [START_REF] Pereira | A fourth-order-accurate Finite volume Ccmpact method for the incompressible Navier-Stokes solutions[END_REF] and Smirnov et al. [START_REF] Smirnov | A finite volume formulation for compact scheme with application to LES[END_REF] proposed a fourth-order finite volume method on structured grids based on the Padé technique and several very high-order schemes have been developed involving compact stencils to provide accurate approximations [START_REF] Piller | Finite volume compact schemes on staggered grids[END_REF][START_REF] Hokpunna | Compact fourth-order finite volume method for numerical solutions of Navier-Stokes equations on staggered grids[END_REF][START_REF] Lacor | A finite volume formulation of compact central schemes on arbitrary structural grids[END_REF][START_REF] Fosso | Curvilinear finite volume schemes unising high-order compact interpolation[END_REF]. The major drawback of the compact technology derives from the restriction to structured grids since it turns to be very complicated for unstructured meshes. Up to our knowledge, the use of very high-order finite volume methods for the (Navier)-Stokes problem with unstructured meshes has only been tackled in [START_REF] Nogueira | High Resolution finite-volume methods on unstructured grids for turbulence and aeroacustics[END_REF][START_REF] Ramírez | A new higher-order finite volume method based on moving least squares for the resolution of the incompressible Navier-Stokes equations on unstructured grids[END_REF] introducing fourth-order and sixth-order schemes based on the application of a mesh-free technique (Moving Least Squares refereed as FV-MLS method) to the finite volume framework.

In the present paper we propose a new sixth-order finite volume method for the steady-state incompressible Stokes equations with unstructured meshes based on the technology initially developed for the convection-diffusion problem in [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF][START_REF] Boularas | A sixth-order finite volume method for diffusion problem with curved boundaries[END_REF]. We use a staggered discretization with a primal unstructured mesh for the pressure and the associated diamond mesh for the velocity to avoid the Rhie-Chow interpolation [START_REF] Rhie | A numerical study of the turbulent flow past an isolated airfoil with trailing edge separation[END_REF][START_REF] Shang | Generalized formulations for the Rhie-Chow interpolation[END_REF]. The coupled velocity-pressure approach is employed to avoid the pressure correction intermediate step to provide the divergence-free velocity [START_REF] Griffith | An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as preconditioner[END_REF]. Moreover, we do not treat the steady-state as the asymptotic limit of an artificial time marching problem but we directly solve the linear system associated to the saddle point problem. The main difficulty is to achieve an efficient approximation of the solution taking into account the divergence-free velocity constraint to determine the pressure. The method is based, on the one hand, in different kinds of polynomial reconstructions to compute the viscous flux, the pressure gradient, and the velocity divergence up to a sixth-order of accuracy and, on the other hand, in a matrix-free formulation using the residual method as in [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF][START_REF] Boularas | A sixth-order finite volume method for diffusion problem with curved boundaries[END_REF]. Unlike the popular second-order methods, the underlying global matrix is not symmetric hence the classical argument to guarantee the existence of divergence-free numerical approximations does not hold any longer. Nevertheless, we show that the method enables to carry out accurate approximations and the algebraic solver (here GMRES) converges to the steady-state solution. A new kind of preconditioning matrix is also given based on the original method proposed in [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF]. The preconditioning technique has been adapted to the specific diamond structure of the dual mesh and we numerically prove the efficiency of our preconditioner.

The document is divided in seven sections. After the introduction, we present in Section 2 all the geometrical ingredients and notations we need to develop the generic finite volume scheme. We detail in Section 3 the different types of polynomial reconstructions to achieve local sixth-order representations of the underlying solution, the numerical fluxes and the numerical scheme and in Section 4 we present the new preconditioning procedure. In Section 5 we address the numerical assessments to show the robustness and accuracy of the proposed techniques and in Section 6 we present the simulation of a 2 polymer extruder apparatus. We end the document with some conclusions.

Finite volume for the Stokes equations

Let Ω be an open bounded polygonal domain of R 2 with boundary BΩ and x " px 1 , x 2 q. We seek functions U " pU 1 , U 2 q " pU 1 pxq, U 2 pxqq, the velocity field, and P " P pxq, the pressure, solutions of the steady-state flow of an incompressible Newtonian fluid governed by the Stokes equations

∇ ¨p´µ∇U `P I 2 q " f, in Ω, (1)
∇ ¨U " 0, in Ω, (2)
where the dynamic viscosity µ " µpxq and the source term f " pf 1 , f 2 q " pf 1 pxq, f 2 pxqq are given regular functions. The tensor ∇U is defined as r∇U s αβ " BUα Bx β , α, β " 1, 2, and I 2 stands for the identity matrix in R 2ˆ2 . The system (1-2) is completed with the Dirichlet boundary condition

U " U D , on BΩ, (3)
where U D " pU 1,D , U 2,D q " pU 1,D pxq, U 2,D pxqq is a given regular function on BΩ which satisfies the compatibility condition

ż BΩ U D ¨n ds " 0,
with n " pn 1 , n 2 q the outward unit normal vector on BΩ. Moreover, uniqueness for the pressure is guaranteed by the additional constraint ş Ω P dx " 0.

Primal and diamond meshes

The primal mesh of Ω, that we denote by M, is a partition of Ω into I non-overlapping convex polygonal cells c i , i P C M " t1, . . . , Iu, and adopt the notations we detail hereafter (see Fig. 1, left):

• for any cell c i , i P C M , we denote by Bc i its boundary and by |c i | its area; the reference cell point is denoted by m i which can be any point in c i (in the present work we shall consider the centroid);

• two cells c i and c j share a common edge e ij whose length is denoted by |e ij | and n ij " pn 1,ij , n 2,ij q is the unit normal vector to e ij outward to c i , i.e. n ij " ´nji ; the reference edge point is m ij which can be any point in e ij (in the present work we consider the midpoint); if an edge of c i belongs to the boundary, the index j is tagged by the letter D;

• for any cell c i , i P C M , we associate the index set νpiq Ă t1, ¨¨¨, Iu Y tDu such that j P νpiq if e ij is a common edge of cells c i and c j or with the boundary if j " D. The diamond mesh of Ω, that we denote by D, derives from the primal mesh M and is constituted of K non-overlapping diamond-shape cell (which degenerate to triangles in the boundary) c k , k P C D " tI `1, . . . , I `Ku. Indeed, for each inner primal edge e ij corresponds a unique cell of the diamond mesh defined by the reference points m i and m j and the vertices of the edge (the dual cell associated to a boundary edge e iD is defined by the reference point m i and the vertices of the edge).

m i m j m ij e ij n ij c i c j m iD n iD n kℓ c ℓ c k m vD n vD m kℓ c v m k e kℓ m ℓ
The notation for the diamond mesh follows the notation introduced for the primal mesh where we substitute the index i P C M by k P C D and the index j P νpiq by P νpkq (see Fig. 1, right). In particular m k is any point in c k (in the present work we shall consider the centroid) and m k is any point in e k (in the present work we consider the midpoint).

To define the association between diamond cells and primal edges, we introduce the correspondence operator Π D such that for given arguments pi, jq, i P C M , j P νpiq, we associate the corresponding diamond cell index k " Π D pi, jq P C D . In the same way, for each diamond edge, we introduce the correspondence operator Π M such that for given arguments pk, q, k P C D , P νpkq, we associate the corresponding primal cell index i " Π M pk, q P C M .

The numerical integrations on the edges are performed with Gaussian quadrature where for the primal edges e ij , i P C M , j P νpiq, we denote by q ij,r , r " 1, . . . , R, their Gauss points and for the diamond edges e k , k P C D , P νpkq, we denote by q k ,r , r " 1, . . . , R, their Gauss points, both sets with weights ζ r , r " 1, . . . , R (see Fig. 2).

Generic finite volume scheme

To provide the generic very high-order finite volume scheme, we first integrate equation (1) over each diamond cell c k , k P C D , and then apply the divergence theorem, yielding ż which can be rewritten in the scalar form as

Bc k p´µ∇U `P I 2 q n ds " ż c k f dx, q ij,2 q ij,1 q ij,3 q kℓ,1 q kℓ,2 q kℓ,3 e ij e kℓ
ż Bc k p´µ∇U β ¨n `P n β q ds " ż c k f β dx, β " 1, 2.
Considering the Gaussian quadrature with R points, i.e. of order 2R, for the line integrals, we get the residual expression

ÿ Pνpkq |e k | |c k | « R ÿ r"1 ζ r `FU β,k ,r `FP β,k ,r ˘ff ´fβ,k " O `h2R k ˘, β " 1, 2, (4)
with the physical fluxes given by F U β,k ,r " ´µpq k ,r q∇U β pq k ,r q ¨nk , F P β,k ,r " P pq k ,r qn β,k , h k " max Pνpkq |e k |, and f β,k an approximation of order 2R of the mean value of f β over cell c k (if cell c k is not triangular, we split it into sub-triangles which share the cell centroid as a common vertex and apply the quadrature rule on each sub-triangle as in [START_REF] Ern | Theory and Practice of Finite Elements[END_REF]). We now integrate equation ( 2) over each primal cell c i and apply again the divergence theorem, yielding ż Bci U ¨n ds " 0.

Considering again Gaussian quadrature with R points for the line integrals, we get the residual expression

ÿ jPνpiq |e ij | |c i | R ÿ r"1 ζ r F ∇ ij,r " Oph 2R i q, (5)
with the physical flux given by

F ∇ ij,r " U pq ij,r q ¨nij
and h i " max jPνpiq |e ij |.

Polynomial reconstructions and high-order finite volume scheme

The polynomial reconstruction is a powerful tool to provide an accurate local representation of the underlying solution and was initially introduced in [START_REF] Barth | Higher order solution of the Euler equations on unstructured grids using quadratic reconstruction[END_REF][START_REF] Barth | Recent developments in high order k-exact reconstruction on unstructured meshes[END_REF] for hyperbolic problems. In [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF] a new methodology was proposed in the context of convection-diffusion problems in order to achieve very accurate approximations of the gradient fluxes and to take into account the boundary conditions. The authors introduced different types of polynomial reconstructions namely the conservative reconstruction in cells and on Dirichlet boundary edges and the non-conservative reconstruction on inner edges, in order to compute approximations of the convective and the diffusive fluxes. We now adapt this technology for the specific Stokes problem where the main difficulty is to handle the two meshes.

Stencil and data

A stencil is a collection of cells situated in the vicinity of a reference geometrical entity, namely an edge or a cell where the number of elements of the stencil shall depend on the degree d of the polynomial function we intend to construct. For each diamond edge e k , k P C D , P νpkq, we associate the stencil S k consisting of the indices of neighbor diamond cells. Analogously, we associate the stencil S k for each diamond cell c k , k P C D , and the stencil S i for each primal cell c i , i P C M , consisting of the indices of neighbor dual and primal cells, respectively. Remark 1. A polynomial reconstruction of degree d requires n d " pd `1qpd `2q{2 coefficients. So, in practice, a stencil consists of the N d closest cells to each geometrical entity (edge or cell) in the respective mesh, with N d ě n d (we consider N d « 1.5n d for the sake of robustness). Now, we want to compute the polynomial reconstructions based on the data of the associated stencil. To this end, we assume that vectors U 1 " pU 1,k q k"I`1,...,I`K , U 2 " pU 2,k q k"I`1,...,I`K , and P " pP i q i"1,...,I gather the approximations of the mean values of U 1 and U 2 over the diamond cells and P over the primal cells, i.e.

U 1,k « 1 |c k | ż c k U 1 dx, U 2,k « 1 |c k | ż c k U 2 dx, P i « 1 |c i | ż ci P dx.

Conservative reconstruction for primal cells

For each primal cell c i , i P C M , the local polynomial approximation of the underlying solution P based on vector P of degree d is defined as

P i pxq " P i `ÿ 1ď|α|ďd R α i rpx ´mi q α ´M α i s ,
where α " pα 1 , α 2 q with |α| " α 1 `α2 and the convention x α " x α1 1 x α2 2 , vector R i " pR α i q 1ď|α|ďd gathers the polynomial coefficients, and M α i " 1 |ci| ş ci px ´mi q α dx in order to guarantee the conservation property

1 |c i | ż ci P i pxq dx " P i .
For a given stencil S i , we consider the quadratic functional

E i pR i q " ÿ qPSi « 1 |c q | ż cq P i pxq dx ´Pq ff 2 . ( 6
)
We denote by p R i the unique vector which minimizes the quadratic functional (6) and we set p P i pxq the polynomial which corresponds to the best approximation in the least squares sense.

Conservative reconstruction for diamond cells

For each diamond cell c k , k P C D , the local polynomial approximation of the underlying functions U 1 and U 2 based on vectors U 1 and U 2 of degree d are defined as

U β,k pxq " U β,k `ÿ 1ď|α|ďd R α β,k rpx ´mk q α ´M α k s , β " 1, 2,
where vector R β,k " pR α β,k q 1ď|α|ďd gathers the polynomial coefficients and

M α k " 1 |c k | ş c k pxḿ k q α dx
in order to guarantee the conservation property

1 |c k | ż c k U β,k pxq dx " U β,k .
For a given stencil S k , we consider the quadratic functional

E β,k pR β,k q " ÿ qPS k « 1 |c q | ż cq U β,k pxq dx ´Uβ,q ff 2 . (7)
We denote by p R β,k the unique vector which minimizes the quadratic functional (7) and we set p U β,k pxq the polynomial which corresponds to the best approximation in the least squares sense.

Non-conservative reconstruction for inner diamond edges

For each inner diamond edge e k , k P C D , P νpkq, the local polynomial approximations of degree d of the underlying functions U 1 and U 2 are defined as

U β,k pxq " ÿ 0ď|α|ďd R α β,k px ´mk q α , β " 1, 2,
where vector R β,k " pR α β,k q 0ď|α|ďd gathers the polynomial coefficients (notice that in this case |α| starts with 0 since no conservation property is required). For a given stencil S k with #S k elements and vector ω β,k " pω β,k ,q q q"1,...,#S k of the positive weights of the reconstruction, we consider the quadratic functional

E β,k pR β,k q " ÿ qPS k ω β,k ,q « 1 |c q | ż cq U β,k pxq dx ´Uβ,q ff 2 . ( 8
)
We denote by r R β,k the unique vector which minimizes the quadratic functional (8) and we set r U β,k pxq the polynomial which corresponds to the best approximation in the least squares sense.

Conservative reconstruction for diamond boundary edges

We treat the boundary diamond edges in a particular way in order to take into account the Dirichlet boundary conditions prescribed for the velocity. For each boundary diamond edge e kD , k P C D , the local polynomial approximations of degree d of the underlying functions U 1 and U 2 are defined as

U β,kD pxq " U β,kD `ÿ 1ď|α|ďd R α β,kD rpx ´mkD q α ´M α kD s , β " 1, 2,
where vector R β,kD " pR α β,kD q 1ď|α|ďd gathers the polynomial coefficients, U β,kD is an approximation of the mean value U β,D of order 2R over the diamond boundary edge e kD , and M α kD " 1 |e kD | ş e kD px ´mkD q α dx in order to guarantee the conservation property

1 |e kD | ż e kD U β,kD pxq ds " U β,kD .
For a given stencil S kD with #S kD elements and vector ω β,kD " pω β,kD,q q q"1,...,#S kD of the positive weights of the reconstruction, we consider the quadratic functional

E β,kD pR β,kD q " ÿ qPS kD ω β,kD,q « 1 |c q | ż cq U β,kD pxq dx ´Uβ,q ff 2 . (9)
We denote by p R β,kD the unique vector which minimizes the quadratic functional (9) and we set p U β,kD pxq the polynomial which corresponds to the best approximation in the least squares sense.

Remark 2. The motivation for introducing the weights in the case of a non-conservative polynomial reconstruction and in the case of a conservative polynomial reconstruction for Dirichlet boundary edges, is presented in [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF] as well as the importance to set larger values for the adjacent cells.

High-order finite volume scheme

This subsection is dedicated to design high-order numerical flux approximations based on the polynomial reconstructions presented in the previous subsections to provide the global residual operator.

Numerical fluxes

For a given polynomial degree d and the associated stencils which guarantee the d-consistency property (see [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF]), four numerical fluxes situations arise:

• for an inner diamond edge e k , the fluxes at the quadrature point q k ,r write

F U β,k ,r " ´µpq k ,r q∇ r U β,k pq k ,r q ¨nk and F P β,k ,r " p P i pq k ,r qn β,k , β " 1, 2,
with the correspondence i " Π M pk, q;

• for a boundary diamond edge e kD , the fluxes at the quadrature point q kD,r write F U β,kD,r " ´µpq kD,r q∇ p U β,kD pq kD,r q¨n kD and F P β,kD,r " p P i pq kD,r qn β,kD , β " 1, 2, with the correspondence i " Π M pk, Dq;

• for an inner primal edge e ij , the flux at the quadrature point q ij,r writes

F ∇ ij,r " p U 1,k pq ij,r qn 1,ij `p U 2,k pq ij,r qn 2,ij ,
with the correspondence k " Π D pi, jq;

• for a boundary primal edge e iD , the flux at the quadrature point q iD,r writes

F ∇ iD,r " U 1,D pq iD,r qn 1,iD `U2,D pq iD,r qn 2,iD .

Residual operators

For any vector Φ " pU 1 , U 2 , Pq in R 2K`I , we define the residual operators for each diamond cell c k , k P C D , as

G β k pΦq " ÿ Pνpkq |e k | |c k | « R ÿ r"1 ζ r `FU β,k ,r `FP β,k ,r ˘ff ´fβ,k , β " 1, 2,
and for each primal cell c i , i P C P , as

G ∇ i pΦq " ÿ jPνpiq |e ij | |c i | R ÿ r"1 ζ r F ∇ ij,r ,
which correspond to the finite volume scheme (4-5) cast in residual form. Gathering all the components of the residuals in vectors G β pΦq " ´Gβ k pΦq ¯k"I`1,...,I`K and G ∇ pΦq " `G∇ i pΦq ˘i"1,...,I , we introduce the global affine operator from R 2K`I into R 2K`I , given by HpΦq " `G1 pΦq, G 2 pΦq, G ∇ pΦq ˘T , such that vector Φ ‹ " pU ‹ 1 , U ‹ 2 , P ‹ q T P R 2K`I , solution of the problem HpΦq " 0, provides a constant piecewise approximation of U 1 , U 2 , and P .

A new preconditioning technique

This section is dedicated to the development of a new efficient preconditioning method to solve the underlying affine problem.

Matrix-free problem and preconditioning

The underlying operator corresponding to the Stokes problem takes the form HpΦq " AΦ ´B "

" A B T C 0  " U P  ´"F V  " 0 (10)
where U " pU 1 , U 2 q T P R 2K gathers the two unknown velocity components, F " pF 1 , F 2 q T P R 2K with F β " pf β,k q k"I`1,...,I`K , β " 1, 2, gathers the source term, and V P R I gathers the prescribed normal velocities on the boundary. Since the matrix A is unknown (and we do not want to explicitly assemble it), we use linear solvers which only require the residual term such as the GMRES method [START_REF] Saad | Iterative methods for spares linear systems[END_REF][START_REF] Sadd | GMRES: a general minimal residual algorithm for solving nonsymmetric linear systems[END_REF].

The preconditioning of the matrix-free problem is a critical stage to dramatically reduce the computational effort and achieve accurate approximations for the affine system. Numerous techniques have been developed, namely the multigrid preconditioner [START_REF] Lonsdale | An algebraic multigrid scheme for solving the NavierStokes equations on unstructured meshes[END_REF][START_REF] Hutchinson | A multigrid method based on the additive correction strategy[END_REF] and the SIMPLE-type preconditioner [START_REF] Klaij | SIMPLE-type preconditioners for cell-centered, colocated finite volume discretization of incompressible Reynolds-averaged Navier-Stokes equations[END_REF][START_REF] Rehman | Preconditioners for the steady incompressible Navier-Stokes problem[END_REF]. Also a very popular method for preconditioning the system uses matrix

P " " A 0 0 S  , (11)
where S " CA ´1B T is the Schur complement of A with respect to A [START_REF] Griffith | An accurate and efficient method for the incompressible Navier-Stokes equations using the projection method as preconditioner[END_REF]. Indeed, using P as a preconditioning matrix theoretically provides the solution with at most four GMRES iterations as proven in [START_REF] Murphy | A Note on preconditioning for indefinite linear systems[END_REF]. Most of the preconditioning techniques have been proposed for the saddle point problem when C " B corresponding to a "symmetric" discretization which preserves the duality between the divergence and the gradient operators. Such a situation arises in many first-and second-order finite volume discretizations such as the DDFV method [START_REF] Delcourte | Preconditioning NavierStokes problem discretized by discrete duality finite volume schemes[END_REF][START_REF] Delcourte | Saddle point preconditioners for linearized Navier-Stokes equations discretized by a finite volume method[END_REF] or in the context of the finite element method [START_REF] Silvester | Efficient preconditioning of the linearized Navier-Stokes equations for incompressible flow[END_REF]. Unfortunately, such symmetry does not hold any longer in the case of very high-order finite volume discretizations and there exists literature which provides compatibility conditions for system [START_REF] Delcourte | Preconditioning NavierStokes problem discretized by discrete duality finite volume schemes[END_REF] to provide existence of a solution [START_REF] Benzi | Numerical solution of saddle point problems[END_REF][START_REF] Haslinger | Projected Schur complement method for solving non-symmetric systems arising from a smooth fictitious domain approach[END_REF].

In the present work, two major difficulties arise for applying such preconditioning technique. First, the matrices A, B, and C are not available since we consider a matrixfree problem. Second, the computation of A ´1 and S ´1 (or an incomplete version such as ILU) requires an important computational effort which turns unsustainable for large linear systems. To overcome these problems, we first introduce a simple and computationally fast approximation of matrices A, B, and C and then propose an efficient preconditioning technique as an extension of the one proposed in [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF].

Approximations for unknown matrices A, B, and C

We introduce the approximation r A P R 2Kˆ2K of A using a Patankar-like discretization (also very similar to the FV4 scheme as in [START_REF] Eymard | Finite volume methods, Handbook of numerical analysis[END_REF]) given by

r Apk ´I, k ´Iq " r Apk `K ´I, k `K ´Iq " ÿ Pνpkq |e k | |c k | µpm k q |m k m | , k P C D ,
for the diagonal entries and by

r Apk ´I, ´Iq " r Apk `K ´I, `K ´Iq " ´|e k | |c k | µpm k q |m k m | , k P C D , P νpkq,
for the extra-diagonal non-null entries.

In the same way, we introduce the approximation r B P R Iˆ2K of B given by (using local indices)

r Bpi, k ´Iq " ÿ Pνpkq with Π M pk, q"i |e k | |c k | n 1,k , i P C M , k P C D , r Bpi, k `K ´Iq " ÿ Pνpkq with Π M pk, q"i |e k | |c k | n 2,k , i P C M , k P C D ,
and the approximation r C P R Iˆ2K of C given by (using local indices)

r Cpi, k ´Iq " |e ij | |c i | n 1,ij , i P C M , j P νpiq, k " Π D pi, jq, r Cpi, k `K ´Iq " |e ij | |c i | n 2,ij , i P C M , j P νpiq, k " Π D pi, jq.
Note that we do not have, a priori, B " C since B depends on the diamond mesh while C depends on the primal mesh.

Incomplete inverse matrix

In [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF], we have proposed a new preconditioning technique based on the evaluation of an incomplete inverse sparse matrix. We here propose an important extension to generalize that idea and introduce a simple way to compute an approximation of A ´1. Let us denote by M a square matrix in R nˆn and denote by E M piq, i " 1, . . . , n, the index set of the non-null entries of row i, that is, j P E M piq if and only if M pi, jq ‰ 0. Definition 1. We say that matrix M has a symmetric structure if E M piq " E M T piq, i " 1, . . . , n, and two matrices M , N have the same structures if E M piq " E N piq, i " 1, . . . , n. Definition 2. For a given square matrix M P R nˆn , the admissible pair pi, jq, i P t1, . . . , nu, j P E M piq, is of type

T 0 if E M piq X E M pjq " H,
and of type T 1 if there exists a unique k ‰ i, j such that

E M piq X E M pjq " k, E M piq X E M pkq " j, E M pjq X E M pkq " i.
We say that matrix M enjoys the T 0 property (resp. T 1 ) if all the admissible pairs are of type T 0 (resp. T 1 ).

Let M and N be two square matrices in R nˆn with the same symmetric structure, i.e.

Epiq " E M piq " E N piq, i " 1, . . . , n, and let R " N M . One can verify the following properties:

• each diagonal element of R is given by

R ii " N ii M ii `ÿ jPEpiq N ij M ji , i " 1, . . . , n; (12)
• @i " 1, . . . , n and @j P Epiq, if pi, jq if of type T 0 then

R ij " N ij M jj `Nii M ij ; (13)
• @i " 1, . . . , n and @j P Epiq, if pi, jq if of type T 1 then

R ij " N ij M jj `Nik M kj `Nii M ij , (14a) R ik " N ik M kk `Nij M jk `Nii M ik , (14b)
with k " Epiq X Epjq.

Of course, matrix R may contain other non-null entries than the ones associated to the admissible pairs pi, jq but the idea is to only consider such pairs to design the incomplete inverse. To this end, we have the following definition.

Definition 3. Assume that M has a symmetric structure such that for each non-null entry M ij , the pair pi, jq is of type T 0 or T 1 . We say that M admits an incomplete left inverse matrix M : if

• the matrix M : has the same structure than M and

• let R be defined by R ii " 1 and R ij " 0, i " 1, . . . , n, j P E M piq, then M : M : " R where relation : " means that we require the equality only for the entries R ij , i " 1, . . . , n, j P tiu Y Epiq.

In other words, we seek an incomplete left inverse matrix M : which have the same structure than M such that the product gives the identity matrix for indices corresponding to the non-null entries of M . Such a specific structure enables to easily compute the incomplete left inverse as we present in the following theorem. Theorem 1. Let M : be the left incomplete inverse of M P R nˆn . Then @i " 1, . . . , n, @j P E M piq:

• if pi, jq is of type T 0 , there exists a coefficient χ ij such that

M : ij " χ ij M : ii with χ ij " ´Mij M jj ; ( 15
) • if pi, jq is of type T 1 , there exist coefficients χ ij , χ ik , with k " E M piq X E M pjq, solution of the linear system " M jj M kj M jk M kk  " χ ij χ ik  " ´"M ij M ik  (16
)

such that M : ij " χ ij M : ii , M : ik " χ ik M : ii . (17)
Moreover, we have

M : ii " 1 M ii `řjPE M piq χ ij M ji . Proof. If pi, jq is of type T 0 , relation (13) yields M : ij M jj `M : ii M ij " 0 then M : ij " ´Mij Mjj M : ii hence M : ij " χ ij M :
ii with χ ij " ´Mij Mjj . Assume now that pi, jq is of type T 1 . Then, relation [START_REF] Eymard | Finite volume methods, Handbook of numerical analysis[END_REF] yields

M : ij M jj `M : ik M kj `M : ii M ij " 0, M : ik M kk `M : ij M jk `M : ii M ik " 0, that we rewrite as " M jj M kj M jk M kk  « M : ij M : ik ff " ´M : ii " M ij M ik  .
Solving the linear system ( 16), we provide coefficients χ ij , χ ik and relations (17) hold.

From relation [START_REF] Diot | The MOOD method in the three-dimensional case: very-high-order finite volume method for hyperbolic systems[END_REF], we have

1 " M : ii M ii `ÿ jPE M piq M : ij M ji .
Having the coefficients χ ij in hand, we compute the diagonal entries M : ii with

1 " M : ii M ii `ÿ jPE M piq χ ij M ji M : ii .
Hence we deduce M :

ii "

1 M ii `řjPE M piq χ ij M ji .
At last we compute M : ij , i " 1, . . . , n, j P E M piq with equations ( 15) and [START_REF] Ferziger | Computational methods for fluids dynamics[END_REF]. l Definition 4. Let N be a R nˆn matrix and M a matrix of symmetric structure with entries of type T 0 or T 1 . We say that M : is a left incomplete inverse preconditioning matrix of N when we solve M : pN x ´bq " 0 in place of solving N x ´b " 0.

Remark 3. Notice that the usual expression "preconditioning matrix" given for example in [START_REF] Saad | Iterative methods for spares linear systems[END_REF] refers to M and not to M ´1 since one solves N x " b using preconditioning technique such as ILU or diagonal matrix where M is simpler than N . In particular, we do not compute M ´1 explicitly, except for some very elementary cases (diagonal matrices).

On the other hand, we use the expression "incomplete inverse preconditioning matrix" to underline that explicitly compute M : and perform the preconditioning multiplying M : with N to reduce the conditioning number of the global linear system. This yields great advantages since the preconditioning procedure just requires, in practice, the product of a sparse matrix with a vector that could be efficiently parallelized.

c o c k c ℓ c k c ℓ c o c p

Projection operators T M and T D

In the numerical section, we need to introduce two projection operators for filtering the matrices in the sense that we shall cancel some extra-diagonal entries to suit to a specific pattern deriving from the mesh connectivities. For a given primal mesh M, we define operator T M on matrices M P R IˆI setting rT M pM qs ii " M ii ; rT M pM qs ij " M ij , j P νpiq; rT M pM qs ij " 0, otherwise.

In short, the operator keeps the extra-diagonal entries belonging to the adjacent cells and cut the other entries. Notice that T M pT M pM qq " T M pM q which justifies the term projection. Moreover, we highlight that νpiq plays the role of Epiq.

Theorem 2. Assume that the primal mesh is a Delaunay triangulation with I cells such that there do not exist three triangles where each triangle is adjacent with the two others as shown in Fig. 3, right panel. Then for any matrix M P R IˆI , associated to the triangular mesh M, the projection result of T M pM q enjoys the T 0 property.

Proof. Let c i be a cell and j P νpiq. If νpiq X νpjq ‰ H then there exists a cell c q q P νpiq X νpjq which share an edge with the two other triangles, i.e. triangles c i , c j , c q corresponds to the configuration of Fig. 3, left panel. Since we assume that such situation is excluded hence q P νpiq X νpjq " H and T M pM q enjoys the T 0 property. l In the same way, for the associated diamond mesh D and any matrix M P R KˆK , we define the operator T D by rT D pM qs kk " M kk ; rT D pM qs k " M k , P νpkq; rT D pM qs k " 0, otherwise.

To deal with the vectorial case (velocity) with matrix M P R 2Kˆ2K of the form

M " " M 11 M 12 M 21 M 22  ,
we extend operator T D pM q setting

T D pM q " " T D pM 11 q 0 0 T D pM 22 q  .
Theorem 3. Assume that a triangular primal mesh satisfies the condition of Theorem 2. Then for any matrix M P R KˆK , the projection result of T D pM q enjoys the T 1 property.

Proof. For any triangular primal cell we have three diamond cells, lets say c k , c , and c p which share two sides between themselves, i.e. c k shares an edge with c and another edge with c p , hence tc p u Ă νpkq X νp q. We claim that there is not another cell (lets say c o ) such that cells c k and c share a common edge with cell c o , as shown in Fig. 3, right panel. Indeed, if there exists such a cell c o different of c p , then c o contains a common edge of c k and c . In the present work we assume that such undesirable configuration is excluded. In conclusion, νpkq X νp q " tc p u such that c k and c are of type T 1 . l Remark 4. Notice that the projection result T D pM q of a matrix M P R KˆK is of type T 1 only if the primal mesh is a triangular mesh. On the other cases, p.e. quadrilateral meshes, matrix T D pM q is of type T 0 .

In the next section, we detail how we compute the incomplete inverse preconditioning matrix r P : based on the approximations r A, r B, and r C and on the procedure we have described before.

Numerical results

This section is dedicated to quantitatively and qualitatively assess the robustness and accuracy of the proposed method. Time consumption and computational scalability assessments are also provided in order to prove the effectiveness of the method and its capacity to be highly parallelizable.

To perform the numerical tests, we consider a fluid with viscosity µ " 1 flowing in a unit square domain Ω " s0, 1r

2 . In order to check the implementation of the method and assess the convergence rates, we manufacture an analytical solution for the given problem setting

U 1 pxq " 1 2
p1 ´cospπx 1 qq sinpπx 2 q, U 2 pxq " 1 2 sinpπx 1 q pcospπx 2 q ´1q , P pxq "

1 2 cos ´π 2 px 1 `x2 q ¯.
Then, the source terms are given by

f 1 pxq " ´π2 2 p2 cospπx 1 q ´1q sinpπx 2 q ´π 4 sin ´π 2 px 1 `x2 q ¯, f 2 pxq " π 2 2 sinpπx 1 q p2 cospπx 2 q ´1q ´π 4 sin ´π 2 px 1 `x2 q ¯.
Boundary conditions derive from the exact solution, namely on the top side we prescribe U D px 1 , 1q " p0, ´sinpπx 1 qq, x 1 P s0, 1r, while we set U D p1, x 2 q " psinpπx 2 q, 0q, x 2 P s0, 1r on the right side. For the other sides, the homogeneous Dirichlet boundary condition U D pxq " p0, 0q, on tBΩ : x 1 ‰ 1, x 2 ‰ 1u is prescribed. We plot in Fig. 4 the isocontours of the x-component of the velocity, U 1 , and of the y-component of the velocity, U 2 , the isocontours of the magnitude of the velocity, ||U ||, and the isocontours of the pressure, P . In all the simulations we have carried out, the weights in functional (8) are set ω β,k ,q " 3, k P C D , P νpkq, q P S k , β " 1, 2, if e k is an edge of c q and ω β,k ,q " 1, otherwise, following [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF].

Accuracy and convergence rates assessment

We first assess the accuracy and the convergence rates of the numerical approximation using the manufactured solution. Vectors U ‹ β " pU ‹ β,k q kPC D , β " 1, 2, and P ‹ " pP ‹ i q iPC M gather the approximate mean values while vectors U β " pU β,k q kPC D , β " 1, 2, and P " pP i q iPC M gather the exact mean values of the solution given by

U β,k " 1 |c k | ż c k U β dx, β " 1, 2, and P i " 1 |c i | ż ci P dx.
The L 1 -norm errors are given by

E β 1 pDq " ÿ kPC D |U ‹ β,k ´U β,k ||c k | ÿ kPC D |c k | , β " 1, 2
, and E P 1 pMq "

ÿ iPC M |P ‹ i ´P ‹ ´P i ||c i | ÿ iPC M |c i |
, and the L 8 -norm errors are given by

E β 8 pDq " max kPC D |U ‹ β,k ´U β,k |, β " 1, 2
, and

E P 8 pMq " max iPC M |P ‹ i ´P ‹ ´P i |.
where P ‹ is the mean value of the values gather in vector P ‹ ,

P ‹ " ÿ iPC M P ‹ i |c i | ÿ iPC M |c i |
, since the GMRES procedure does not guarantee a solution of null mean value.

Remark 5. One can verify that the problem given by equations (1-2) is singular for the pressure since if P a solution, then P `C, C P R, is also a solution. We have consider no Dirichlet condition for the pressure (or other procedure to fix the singularity) since we have noticed that the GMRES procedure always finds the null mean value solution. However, such property does not hold for a preconditioned GMRES and therefore no convergence is achieved (the C value does not converge when we consider a finner mesh although the gradient of P is well evaluated). To overcome the problem, we simply fix P ‹ i subtracting the mean value of P ‹ such that one gets a unique piecewise approximation of P with a null mean value.

We evaluate the convergence rate of the L 1 -norm (and L 8 -norm error) between two different and successive finer primal meshes M 1 and M 2 , with I 1 and I 2 cells, respectively, as

O P 1 pM 1 , M 2 q " 2 | logpE P 1 pM 1 q{E P 1 pM 2 qq| | logpI 1 {I 2 q| .
In the same way, we define the convergence order between two different and successive finer diamond meshes D 1 and D 2 with K 1 and K 2 cells, respectively, as As a first test, we carry out simulations with successive finer regular triangular Delaunay meshes (see Fig. 5, left) and the associated diamond meshes (see Fig. 5, right). We report in Tables 1,2, and 3 the L 1 -and L 8 -norm errors and the convergence rates using the P 1 , P 3 , and P 5 polynomial reconstructions, respectively, where the number of unknowns (the same as degrees of freedom) is DOF " K for U 1 and U 2 and DOF " I for P . The notation E 1 is a generalization which stands for E β 1 or E P 1 depending on the variable we are dealing with (U β or P , respectively). The same convention is valid for E 8 , O 1 , and O 8 .

O β 1 pD 1 , D 2 q " 2 | logpE β 1 pD 1 q{E β 1 pD 2 qq| | logpK 1 {K 2 q| .
The P 1 polynomial reconstruction provides a second-order approximation for the velocity and a first-order approximation for the pressure. The scheme based on the P 3 reconstruction achieved an effective fourth-order approximation for the velocity and a third-order (slightly better) approximation for the pressure and scheme based on the P 5 reconstruction achieved a sixth-and fifth-order approximations for the velocity and the pressure, respectively. We also mention that no oscillations or numerical locking are reported in all the experiences.

Scheme robustness and accuracy assessments with deformed meshes are also of crucial importance in order to check the method capacity to handle complex meshes still preserving high-order convergence rates. To this end, we consider successive finer deformed quadrilateral meshes (see Fig. 6) applying a random displacement of each inner vertex of structured meshes controlled by a deformation factor (see the detailed procedure in [START_REF] Clain | A sixth-order finite volume method for the convection-diffusion problem with discontinuous coefficients[END_REF]). In the present experience, we choose a test case with 30% of deformation. We report in Tables 4,5, and 6 the L 1 -and L 8 -norm errors and convergence rates using the P 1 , P 3 , and P 5 polynomial reconstructions, respectively (the meaning of DOF is the one just presented).

We observe that the scheme correctly handles complex meshes and the convergence rates are optimal both for the velocity and for the pressure. We achieved up to a sixthorder convergence rate for the velocity using the P 5 polynomial reconstruction and no oscillations were noticed in all the tests. 20

Assessment of the preconditioning technique

We devote this subsection to assess the quality of the preconditioning technique proposed in Section 4. Two topics are here considered: the use of approximate matrices r A, r B, and r C deriving from a Patankar-like discretization to compute an approximation of P and the use of the incomplete inverse matrix for preconditioning. To do so, we consider four situations.

• Accordingly to Section 4, matrix P defined by relation [START_REF] Delcourte | Saddle point preconditioners for linearized Navier-Stokes equations discretized by a finite volume method[END_REF] stands for the ideal preconditioning matrix as shown in [START_REF] Murphy | A Note on preconditioning for indefinite linear systems[END_REF]. Nevertheless, computation P ´1 is not possible since the Schur complement S is not invertible. To overcome such a difficulty, we add a slight perturbation setting S ε " CA ´1B T `εI with ε ą 0 small enough (we take ε " 10 ´8 in the numerical tests). Since S ε is now invertible, we define

P ´1 ε " " A ´1 0 0 S ´1 ε  .
• To evaluate the impact of the approximate matrices deriving from the Patankarlike discretization, we shall consider a preconditioning matrix r P based on r A, r B, and r C. One more time, the Schur complement r S " r C r A ´1 r B T is not invertible and we slightly perturbe the matrix setting r S ε " r S `εI such that r P ε is invertible and we set

r P ´1 ε " « r A ´1 0 0 r S ´1 ε ff .
• To assess the efficiency of the incomplete inverse procedure, we substitute the inverse operator ¨´1 by the ¨: one. To deal with the Schur matrix, we use projector operator T M defined in subsection 4.4 and we set p r S " T M ´r C r A : r B T ¯. Indeed, matrix r C r A : r B T has non-null entries which do not correspond to adjacent cells hence we apply the T M projector to guarantee that matrix p r S enjoys the T 0 property given by definition 2. Since p r S now satisfies the T 0 property and r A satisfies the T 0 or T 1 properties, we can apply the ¨: operator and we define

r P : " « r A : 0 0 p r S :
ff .

• For the last issue, we want to substitute the inverse operator ¨´1 by the ¨: when one employs the exact matrices A, B, and C of the linear system [START_REF] Delcourte | Preconditioning NavierStokes problem discretized by discrete duality finite volume schemes[END_REF]. Since S does not enjoy the T 0 property and A does not enjoy the T 0 or T 1 properties for quadrilateral or triangular primal meshes, respectively, a specific treatment is required to apply the ¨: operator. For matrix A we simply apply the T D operator setting p A " T D pAq. To build an approximation of the Schur matrix, we define p S " T 0 `C p A : B T ˘which turns to be a matrix of type T 0 . We then define the incomplete inverse preconditioning matrix by

P : " « p A : 0 0 p S :
ff .

Remark 6. We do not explicitly build the matrices A, B, and C since the method only deals with operator H. Nevertheless it is possible to recover the matrix by computing HpΦq where Φ spans the canonical basis. In that way, we obtain all the columns of the matrices. Of course, such a technique is not computationally interesting but it enables to evaluate explicitly all the matrices to carry out the numerical tests of this subsection.

Efficiency assessments of the preconditioning matrices

To perform the simulations, we consider a uniform triangular Delaunay mesh with 944 primal cells and its associated diamond mesh with 1454 cells corresponding to a linear system with 3852 unknowns. We plot in Fig. 7 the residual curves of the GMRES method for P 1 , P 3 , and P 5 polynomial reconstructions, where R iter is the L 2 -norm residual at iteration iter, with no preconditioning and with preconditioning, where for the latter we consider the four preconditioning matrices just presented. The main observations are the following.

• The exact inverse matrix P ´1 ε enables to determine the exact solution in 4 iterations as expected and the small perturbation controlled by ε does not perturb the convergence. The preconditioner is efficient for all the reconstructions and suggest that approximation of P ´1 ε will provide a very good matrix for reducing the conditioning number of matrix A.

• The use of any preconditioning matrix dramatically reduces the number of iterations with respect to the identity matrix (no preconditioning). We also remark that the number of iterations is mainly the same, independent of the polynomial degree reconstruction since the number of unknowns is preserved (only the stencil size changes, hence the matrix connectivity) with the exception of the nonpreconditioning case where the computational effort increases with the polynomial degree.

• Substituting the exact matrices A, B, and C with their respective approximations and using the exact inverse matrix r P ´1 ε provides excellent results, even with higher degree. For the global system with 3832 unknowns we just need about 80 iterations to provide the solution. In other word, the use of Patankar-like matrices to perform the preconditioning of the matrix-free problem is very efficient.

• We now deal with the incomplete inverse preconditioning matrix using the exact matrices or the approximated matrices. We observe that the two incomplete inverse matrices provide approximatively the same iterations numbers which means that the loose of performance may be essentially imputed to the inverse procedure and not the matrix approximation. In all the situations, we cut by 3 or 4 the number of iterations but we recall that the preconditioning procedure is highly parallelizable since it just corresponds to a sparse matrix product with a vector.

As an overall conclusion, the use of the Patankar-like matrices in substitution of the real matrices is very efficient and the incomplete inverse preconditioning enables to substantively reduce the number of iterations of the GMRES procedure.

Not Preconditioned Preconditioned

Figure 8: Comparison of the residual curves of the GMRES method with and without preconditioning matrix for P 1 polynomial reconstructions.

Convergence and saturation for fine meshes

We now turn to assess the preconditioning using r P : when dealing with bigger linear systems. For that purpose, we consider a uniform triangular Delaunay mesh with 16374 primal cells and 24719 diamond cells corresponding to 65812 unknowns.

We plot in Fig. 8 the residual curves of the GMRES method with no preconditioning and with the preconditioning matrix r P : for the P 1 reconstruction while we present the results for the other reconstruction in Table 7 where N ITER is the number of iterations to achieve a residual of R iter ă 10 ´14 , T E is the execution time (normalized for the sake of comparison), and T PREC is the percentage of T E in the preconditioning of the residual vectors during the GMRES procedure.

We observe that the non-preconditioning version converges very slowly and the dependance on the polynomial degree is noticeable. On the contrary, the new preconditioning procedure is low sensitive to the polynomial degree. We also note the very low computational cost of the preconditioning procedure, the huge differences of running time between with and without preconditioning, and the increase of the computational effort just by 20% and 60% fot the P 3 and P 5 versions, respectively, which demonstrates the great advantage to perform high-order approximation on coarse meshes when comparing with second-order approximations on fine meshes. Of course fine meshes could be necessary to capture specific structures of the solution (vortex, obstacle, discontinuities) but the very high-order approximation should be privileged when possible. In spite of that positive conclusion, a more representative preconditioning matrix based on the second-order discretization should be considered in the future. Indeed, the matriz approximations r A, r B, and r C derive from a Patankar-like discretization which is non-consistent for unstructured meshes.

We carried out the same simulation using deformed quadrilateral meshes and we noticed the same effectiveness of the proposed preconditioning as for the presented tests.

Scalability and time consumption

Multithreading/multicore implementation is of crucial importance when dealing with a large number of unknowns and long running times (larger than several hours for instance). Computers are equipped with multicore processors and GPU cards allowing efficient time simulations cut by running specific portions of the code in parallel. High Performance Computing aims to take advantage of algorithms supporting the multithreading so we have developed and implemented a parallel version of the method using the openMP framework for a small multicore machine to assess the scalability and the potentiality to be parallelized for large multicore machines.

Speed-up assessment

For each simulation, we evaluate the following parameters:

• T PREC -execution time to build the preconditioning matrix;

• T GMRES -execution time to solve the linear system using the GMRES method (including the preconditioning, the Krylov basis construction, and the residual vectors computation);

• T E -execution time of the entire simulation;

• S n -speedup with n cores given by

S n " T E,1 T E,n , (18)
where T E,1 is the execution time with one CPU core and T E,n is the execution time with n CPU cores;

• E n -strong efficiency with n cores (also given in percentage) is defined as

E n " S n n . (19)
All the time parameters correspond to the wall-clock time and the measurements are given in seconds. The efficiency of an algorithm is related to the relative usage of each allocated core. Very scalable algorithms have high values of efficiency and the Amdahl's Law predicts the speed-up with n cores as where P n is the fraction of the algorithm running in parallel. Notice that, according to Amdahl, a linear speed-up is not usually achieved since it implies a full parallel algorihtm, which is not always possible. Moreover, the Amdahl's Law does not take into account some effects such as syncronization, overhead, or cache saturation, which can reduce the speed-up or even increase it. In fact, a superlinear speed-up (S n ą n) is, in practice, possible since more CPUs implies a larger available cache memory, avoiding cache saturation effects. Since we do not know, a priori, the value of P n , we cannot predict the speed-up using the Amdahl's Law. Instead, we obtain P n using equation [START_REF] Gao | The finite volume projection method with hybrid unstructured triangular collocated grids for incompressible flows[END_REF] for a given speed-up and number of cores, as a indicator of the quality of the scalability. Simulations were carried out with uniform Delaunay meshes and we perform the computation with n " 1, 2, 4, 8 cores. Having T E,1 and T E,n in hand, we compute the speed-ups with relation [START_REF] Fosso | Curvilinear finite volume schemes unising high-order compact interpolation[END_REF], the efficiency with equation [START_REF] Frochte | A splitting technique of higher order for the Navier-Stokes equations[END_REF], and finally deduce P n from the Amdahl's Law [START_REF] Gao | The finite volume projection method with hybrid unstructured triangular collocated grids for incompressible flows[END_REF]. The algorithm was implemented in C``and parallelized using the openMP framework. The machine has 8 Intel Xeon processors with 2.2 GHz of clock rate and 2 MB of cache memory for each core. We display the speed-up and the efficiency for the P 1 , P 3 , and P 5 polynomial reconstructions in Fig.s 9, 10, and 11, respectively, where DOF is the total number of unknowns, i.e. DOF " 2K `I.

S n " 1 1 ´Pn `Pn n . ( 20
We report a very good scalability using 2 CPU cores with speed-ups close to 2 for the three kinds of reconstruction. We obtain larger speed-ups allocating more CPU cores until reaching a speed-up around 4.5 with 8 cores with small variations with respect to the polynomial degree and the mesh size (except for very small meshes where all the data is contained in the cache memory). Such a situation is expected since more 0 20 threads lead to larger overhead delays, that is, task start-up time, synchronizations, data communications, task termination time, among other factors, slowing down the simulations. Moreover, the bottleneck between the main and the cache memory slows down the data transfer when dealing with more than two cores. We notice some efficiency rates above 100% because some cache memory effects can lead to a superlinear speedup in some scenarios. Moreover, the smallest mesh has the worst efficiency (36% with speed-up of 2.9) since the overhead affects very much the simulation.

Considering the Amdahl's Law, we predict a fraction of parallelized code around P " 90% which is a indicator of a good scalability. Of course, as stated before, the Amdahl's prediction does not take into account some factors which can slow down the algorithm.

To sum up, the tests we have carried out show that the method we propose has a great potential to be parallelized for large multicore machines to drastically reduce the execution time.

Time consumption versus mesh size

Tables 8,9, and 10 provide the execution time per task and the total execution time, where DOF " 2K `I. The time computation of the preconditioning matrices and the incomplete inverse r P : during the pre-processing stage is only given in Table 8 since it does not depend on the polynomial reconstruction degree. The pre-processing stage is very low consuming with respect to the whole computational process and the essential cost derives from the GMRES routine. We observe that the P 3 situation generates an overcost of about 30% whatever the number of core and the mesh size are. For the P 5 reconstruction, the situation is less clear: for small mesh size the ratio between the time consumption with P 5 and P 1 ranges between 3 and 4 in function of the cores number. For the two finest meshes, the ratio is stabilized around a factor 2 independently of the core numbers. At last we observe that the time consumption roughly increases as I 3{4 or, taking into account the space dimension, the time consumption increases as h ´3{2 where h " 1{ ? I is the mesh parameter and we observe low sensitivity regarded to the number of cores and the polynomial degree. Of course the time is reduced when employing more cores or lower degree but the time is still a function of order h ´3{2 .

Time consumption versus the error of approximation

We now compare the computational cost of the proposed method in terms of accuracy. For a given tolerance, we want to determine the mesh and the associated execution time one has to use to provide an accurate approximation up to the prescribed error for several type of reconstructions. To do so, we compute the L 2 -norm errors of the velocity and the pressure, given by

E β 2 pDq " ¨ÿ kPC D |c k |pU ‹ β,k ´U β,k q 2 ÿ kPC D |c k |U 2 β ‹ ‹ ' 1 2
and E P 2 pMq " using successive finer uniform triangular Delaunay meshes and the execution time of the simulation, T E . We plot in Fig. 12 the L 2 -norm error as a function of the execution time with P 1 , P 3 , and P 5 polynomial reconstructions for four meshes. Clearly the P 5 reconstruction has the lower computational cost to achieve an approximate solution up to a given tolerance. As an example, to provide a numerical solution with one core such that the L 2 -error of the velocity is lower than 10 ´6, the P 5 version requires 10s with a mesh of 2500 cells, the P 3 reconstruction needs 100s with a mesh of 12000 cells, and the P 1 second-order method performs the computation in 10 8 s (about 3 years and 2 months) with a mesh of about 10 7 cells (estimated value).

¨ÿ iPC M |c i |pP ‹ i ´P i q 2 ÿ iPC M |c i |P 2 i ‹ ‹ ' 1 2 , 1E0 1E1 1E2 1E3 1E4 1E0

Simulation of a polymer extruder apparatus

We dedicate this section to the polymer extrusion simulation in order to highlight the method capacity to handle complex geometries with unstructured meshes within the very high-order finite volume context. We also intend to evaluate the preconditioning technique efficiency in the case of a concrete application. For this purpose, we consider a geometry which corresponds to a standard extruder machine [START_REF] Gonçalves | Design of complex profile extrusion dies through numerical modeling[END_REF] displayed in Fig. 13 where a molten polymer material flows out to the extruder head with a specific constant cross section shape (plastic pipes, window frames, etc.). Due to the high viscosity and the low velocity, the nonlinear contribution of the Navier-Stokes equation is neglected and we assume that the process is governed by the steady-state Stokes equation for incompressible fluids. To prescribe the boundary conditions, we consider an inlet velocity U D,in defined on x P Γ in " tpx 1 , x 2 q P BΩ : x 1 " 0 ^p0 ă x 2 ă 0.2 _ 0.8 ă x 2 ă 1qu given by U D,in p0, x 2 q " # p4x 2 p0.2 ´x2 q, 0q , if 0 ă x 2 ă 0.2, p6p1 ´x2 qpx 2 ´0.8q, 0q , if 0.8 ă x 2 ă 1, and an outlet velocity U D,out defined on x P Γ out " tpx 1 , x 2 q P BΩ : x 1 " 1.4 ^p0.4 ă x 2 ă 0.45 _ 0.55 ă x 2 ă 0.6qu given by U D,out p1.4, x 2 q " # p320p0.45 ´x2 qpx 2 ´0.4q, 0q , if 0.4 ă x 2 ă 0.45, p320p0.6 ´x2 qpx 2 ´0.55q, 0q , if 0.55 ă x 2 ă 0.6.

The other portions of the boundary are walls with null velocity U D pxq " p0, 0q. We consider a null source term, that is f pxq " p0, 0q, and we assume the viscosity to be unitary, that is µ " 1, for the sake of simplicity since the problem is linear. We carried out three simulations with triangular Delaunay meshes as primal meshes (with 2406, 9354, and 26572 cells) and deduce the respective diamond meshes (with 3750, 14334, and 40380 cells, respectively). We display a coarse version of the primal and dual meshes in Fig. 14. Fig. 15 shows the isocontours of the horizontal (top panel) and vertical (middle panel) components of velocity together with the isocontours of the pressure (bottom panel) for the finer mesh and with P 5 polynomial reconstructions. An enlargement near the extruder outlet is plotted on the right side to detail the flow structures. The isolines are slightly asymmetric since the inflow conditions are not the same in the upper and lower inflow pipes. To assess the accuracy of the different polynomial reconstruction, we report in Table 11 the minimum and the maximum mean values of the three scalar quantities, where DOF " K for U 1 and U 2 and DOF " I for P . The minimum of the horizontal and vertical velocities presents the main variations. We notice that the fourth-and sixth-order schemes provide very similar values whereas we observe larger differences between the second-and fourth-order ones. The most noteworthy difference occurs with the pressure where the P 5 reconstruction with the smaller mesh provides the same approximation than the P 1 reconstruction with the middle size mesh. The fourth-and sixth-order methods provide very similar pressure results with the finer mesh (difference of about 7 units for the maximum) whereas the second-order scheme is 80 units far from the maximum pressure using the P 5 reconstruction. To assess the efficiency of the preconditioning method, we have carried out the previous simulations with and without the inverse preconditioning matrix and we plot the L 2 -norm residuals curves for the P 1 polynomial reconstruction in Fig. 16 (we do not represent the curves with the other reconstructions since they are very similar). We report 12 a set of output parameters for the simulations we have done.

We observe that the non-preconditioning versions do not converge and the residual norm saturate above the prescribe tolerance threshold. We report the minimum residual we achieved and notice the low sensitivity to the polynomial degree and the high dependency to the mesh size. When the preconditioning technique is applied, the GMRES method always converge up to the prescribed tolerance and we report the number of iterations that took to obtain convergence and the associated running time.

Conclusion

We have presented an high-order finite volume scheme to solve the bidimensional incompressible Stokes problem based on a new class of polynomial reconstructions. The scheme achieves an effective sixth-order accuracy for the velocity and the fifth-order accuracy for the pressure. We have also presented a new preconditioning technique based on the Schur complement. A large number of numerical tests were carried out to prove its high-efficiency in reducing the computational effort.

Figure 1 :

 1 Figure 1: Notation for the primal mesh (left) and for the diamond mesh (right).

Figure 2 :

 2 Figure 2: Gauss points on the edges of the primal cells (dashes lines) and on the edges of the diamond cells (solid lines).

Figure 3 :

 3 Figure 3: Correct configuration (left) for a triangular mesh where the corresponding diamond cells c k and c share a unique cell co, and undesirable configuration (right) where the diamond cells c k and c share two cells, co and cp.

Figure 4 :

 4 Figure 4: Isocontours of the x-component (top, left) and the y-component of the velocity (top, right), isocontours of the velocity magnitude (bottom, left), and isocontours of the pressure (bottom, right).

Figure 5 :

 5 Figure 5: A coarse uniform triangular Delaunay mesh (left) and the associated diamond mesh (right).

Figure 6 :

 6 Figure 6: A coarse deformed quadrilateral mesh (left) and the associated diamond mesh (right).

Figure 7 :

 7 Figure 7: Comparison of the residual curves of the GMRES method with different preconditioning matrices for P 1 (top), P 3 (middle), and P 5 (bottom) polynomial reconstructions.

Figure 9 :

 9 Figure 9: Speed-ups and efficiency using uniform triangular Delaunay meshes with P 1 polynomial reconstructions.

Figure 11 :

 11 Figure 11: Speed-ups and efficiency using uniform triangular Delaunay meshes with P 5 polynomial reconstructions.

Figure 13 :

 13 Figure 13: Polymer extruder machine geometry where a rotation screw (1) forces the molten polymer through a die (3) after passing through the breaker plate (2). The breaker plate creates back pressure in the barrel required for uniform melting and proper mixing of the polymer.

Figure 14 :

 14 Figure 14: A coarse non-uniform triangular Delaunay mesh (left) and its associated diamond mesh (right).

Figure 15 :Figure 16 :

 1516 Figure 15: Isocontours of the x 1 -component of the velocity (top), isocontours of the x 2 -component of the velocity (middle), and isocontours of the pressure (bottom), and respective exit zooms.

Table 1 :

 1 Errors and convergence rates for the P 1 scheme with uniform triangular Delaunay primal meshes.

		DOF	E 1	O 1	E 8	O 8
		363	1.23E´03	-	9.54E´03	-
	U 1	1454 3.34E´04 1.88 2.22E´03 2.10 6135 7.96E´05 1.99 6.80E´04 1.64
		24719 1.90E´05 2.06 1.35E´04 2.33
		363	1.24E´03	-	5.31E´03	-
	U 2	1454 3.42E´04 1.86 1.88E´03 1.50 6135 7.78E´05 2.06 5.10E´04 1.81
		24719 1.91E´05 2.01 1.22E´04 2.05
		230	5.23E´02	-	2.32E´01	-
	P	944 4038 9.32E´03 1.20 6.32E´02 1.05 2.23E´02 1.20 1.36E´01 0.76
		16374 4.21E´03 1.14 3.43E´02 0.87

Table 2 :

 2 Errors and convergence rates for the P 3 scheme with uniform triangular Delaunay primal meshes.

		DOF	E 1	O 1	E 8	O 8
		363	3.83E´05	-	1.56E´04	-
	U 1	1454 1.98E´06 4.27 1.34E´05 3.54 6135 1.08E´07 4.04 8.81E´07 3.78
		24719 6.42E´09 4.05 4.58E´08 4.24
		363	3.89E´05	-	1.88E´04	-
	U 2	1454 2.04E´06 4.25 1.43E´05 3.71 6135 1.07E´07 4.09 6.78E´07 4.24
		24719 6.41E´09 4.04 4.75E´08 3.81
		230	1.24E´03	-	6.86E´03	-
	P	944 4038 1.17E´05 3.36 1.23E´04 3.29 1.34E´04 3.15 1.34E´03 2.32
		16374 1.29E´06 3.15 1.82E´05 2.73

Table 3 :

 3 Errors and convergence rates for the P 5 scheme with uniform triangular Delaunay primal meshes.

		DOF	E 1	O 1	E 8	O 8
		363	1.34E´06	-	1.59E´05	-
	U 1	1454 1.40E´08 6.57 1.57E´07 6.66 6135 1.48E´10 6.32 1.01E´09 7.01
		24719 2.66E´12 5.77 1.70E´11 5.87
		363	1.34E´06	-	1.69E´05	-
	U 2	1454 1.44E´08 6.53 1.58E´07 6.74 6135 1.49E´10 6.35 1.31E´09 6.66
		24719 2.76E´12 5.73 1.84E´11 6.12
		230	4.39E´05	-	4.11E´04	-
	P	944 4038 1.32E´08 5.45 1.12E´07 6.07 6.91E´07 5.88 9.19E´06 5.38
		16374 3.65E´10 5.13 4.42E´09 4.62

Table 4 :

 4 Errors and convergence rates for the P 1 scheme with deformed quadrilateral primal meshes.

		DOF	E 1	O 1	E 8	O 8
		220	2.06E´03	-	6.78E´03	-
	U 1	840 3280 1.08E´04 2.12 5.65E´04 1.80 4.58E´04 2.24 1.93E´03 1.88
		12960 2.62E´05 2.06 1.48E´04 1.95
		220	1.94E´03	-	8.35E´03	-
	U 2	840 3280 1.08E´04 2.05 7.74E´04 2.22 4.36E´04 2.23 3.52E´03 1.29
		12960 2.63E´05 2.06 2.44E´04 1.68
		100	6.75E´02	-	2.23E´01	-
	P	400 1600 9.07E´03 1.41 6.83E´02 1.18 2.41E´02 1.48 1.55E´01 0.53
		6400 4.02E´03 1.17 3.61E´02 0.93

Table 5 :

 5 Errors and convergence rates for the P 3 scheme with deformed quadrilateral primal meshes.

		DOF	E 1	O 1	E 8	O 8
		220	6.05E´05	-	3.64E´04	-
	U 1	840 3280 1.83E´07 4.22 1.88E´06 3.46 3.25E´06 4.36 1.98E´05 4.35
		12960 1.11E´08 4.08 9.87E´08 4.29
		220	5.65E´05	-	2.97E´04	-
	U 2	840 3280 1.80E´07 4.16 1.50E´06 3.89 3.06E´06 4.35 2.12E´05 3.94
		12960 1.12E´08 4.04 1.25E´07 3.61
		100	2.10E´03	-	1.03E´02	-
	P	400 1600 1.60E´05 3.56 2.22E´04 3.03 1.89E´04 3.47 1.81E´03 2.51
		6400 1.65E´06 3.28 3.09E´05 2.85

Table 6 :

 6 Errors and convergence rates for the P 5 scheme with deformed quadrilateral primal meshes.

		DOF	E 1	O 1	E 8	O 8
		220	3.08E´06	-	2.26E´05	-
	U 1	840 3280 4.37E´10 6.38 5.28E´09 5.31 3.38E´08 6.73 1.97E´07 7.08
		12960 6.35E´12 6.16 5.62E´11 6.61
		220	2.72E´06	-	1.46E´05	-
	U 2	840 3280 4.43E´10 6.26 3.37E´09 6.68 3.14E´08 6.66 3.20E´07 5.70
		12960 6.36E´12 6.18 6.11E´11 5.84
		100	5.84E´05	-	3.18E´04	-
	P	400 1600 2.92E´08 5.53 4.73E´07 5.48 1.35E´06 5.44 2.12E´05 3.91
		6400 8.03E´10 5.18 1.63E´08 4.86

Table 7 :

 7 Number of GMRES iterations and execution time using uniform triangular Delaunay meshes.

			N ITER	T E	T PREC [%]
	P 1	I d r P :	8792 4700	1.36 1	-3.4
	P 3	I d 25000* 3.63* r P : 4651 1.18	-2.9
	P 5	I d 43000* 9.15* r P : 4662 1.61	-2.3

* Estimated value

Table 8 :

 8 Runtime per task and multithreading performance of the P 1 scheme using uniform triangular Delaunay meshes.

	DOF n T PREC [s] T GMRES [s]	T E [s]
		1	0.04	0.47	0.57
	956	2 4	0.02 0.01	0.26 0.06	0.33 0.26
		8	0.01	0.11	0.20
		1	0.53	10.10	10.94
	3852	2 4	0.28 0.15	5.44 2.97	5.93 3.26
		8	0.09	1.84	2.01
		1	9.31	238.22	249.65
	16308	2 4	4.81 2.47	128.10 81.39	134.17 84.58
		8	1.30	61.53	63.28
		1	169.73	5319.71	5510.26
	65812	2 4	87.87 46.12	2743.64 1741.87	2843.42 1794.20
		8	22.57	1145.93	1172.36

Table 9 :

 9 Runtime per task and multithreading performance of the P 3 scheme using uniform triangular Delaunay meshes.

	DOF n T GMRES [s]	T E [s]
		1	0.95	1.39
	956	2 4	0.44 0.26	0.68 0.39
		8	0.16	0.24
		1	14.23	16.43
	3852	2 4	8.10 4.06	9.28 4.68
		8	3.39	3.74
		1	302.00	319.30
	16308	2 4	161.46 101.50	170.55 106.19
		8	75.66	78.20
		1	6541.88	6759.10
	65812	2 4	3228.46 2123.95	3342.11 2183.35
		8	1477.90	1508.14

Table 10 :

 10 Runtime per task and multithreading performance of the P 5 scheme using uniform triangular Delaunay meshes.

	DOF n T GMRES [s]	T E [s]
		1	1.95	4.06
	956	2 4	0.94 0.47	2.00 1.02
		8	0.28	0.59
		1	23.69	32.24
	3852	2 4	13.11 8.27	17.59 10.53
		8	6.92	8.11
		1	477.74	523.93
	16308	2 4	252.88 159.29	276.73 171.51
		8	115.87	122.14
		1	10673.78	11058.73
	65812	2 4	5234.73 3108.89	5433.43 3215.69
		8	2474.25	2526.57

Table 11 :

 11 Maximum and minimum mean values for U 1 , U 2 , and P . 3750 ´7.05E´2 7.45E´2 ´7.85E´2 7.51E´2 ´7.65E´2 7.43E´2 14334 ´7.59E´2 7.63E´2 ´7.74E´2 7.73E´2 ´7.75E´2 7.71E´2 40380 ´7.78E´2 7.76E´2 ´7.82E´2 7.82E´2 ´7.81E´2 7.81E´2 P 2406 ´80948.13 9868.84 ´88824.25 10707.79 ´88096.95 10711.55 9354 ´88784.94 10699.58 ´90888.84 10938.38 ´90796.49 10932.28 26572 ´90753.83 10907.28 ´91526.43 10993.90 ´91472.65 10987.23

			P 1		P 3		P 5
		DOF	Min	Max	Min	Max	Min	Max
		3750 ´1.17E´3 1.98E´1 ´1.64E´3 1.99E´1 ´2.35E´3 2.00E´1
	U 1	14334 ´6.37E´4 2.00E´1 ´2.31E´4 2.00E´1 ´2.53E´4 2.00E´1
		40380 ´7.36E´4 2.00E´1 ´1.10E´4 2.00E´1 ´1.83E´4 2.00E´1
	U 2

Table 12 :

 12 Number of GMRES iterations (N ITER ), final L 2 -norm residual (R ITER ), and execution time in seconds (T E ).

	DOF		9906			38022		107332
		N ITER	R ITER	T E [s] N ITER	R ITER	T E [s] N ITER	R ITER	T E [s]
	P 1	I d r P : 1220 NC	2.39E´07 9.83E´15 10.30 -	NC 2805	8.06E´08 9.86E´15 161.55 5184 -NC	1.19E´07 9.97E´15 1621.02 -
	P 3	I d r P : 1417 NC	1.15E´07 9.38E´15 17.01 -	NC 2555	3.23E´08 9.95E´15 191.67 5043 -NC	1.41E´08 9.93E´15 1691.06 -
	P 5	I d r P : 1599 NC	7.18E´08 9.59E´07 29.49 -	NC 3109	2.02E´08 9.71E´15 299.33 6513 -NC	9.15E´09 1.00E´14 2180.70 -

NC: No Convergence

in Table

Acknowledgements

This research was financed by FEDER Funds through Programa Operational Fatores de Competitividade -COMPETE and by Portuguese Funds FCT -Fundação para a Ciência e a Tecnologia, within the Projects PEst-C/MAT/UI0013/2014 and FCT-ANR/MAT-NAN/0122/2012.


Reiffers Alexandre

Yezekael Masson

Eitan Hayel

Altman

Game theory approach for modeling competition over visibility on social networks

In Online Social Networks, such as Facebook, Linkedin, Twitter, Google+ and others, many members post messages to walls or to timelines of their friends. There is a permanent competition over content visibility since timelines have finite capacity. As new content arrives, older content gets pushed away from the timeline. A selfish source that wishes to be visible has to send from time to time new content thus pre-empting other content from the timeline. We assume that sending more content comes with some extra cost. We study the problem of selecting the rates of content creation as a non-cooperative game between several sources that share some common destination to which they send content. We identify conditions under which the problem can be reduced to the Kelly mechanism for which we compute explicitly the equilibrium. This is done in a very general probabilistic framework where time between arrivals of content is only required to be stationary ergodic.

Introduction

The visibility of contents in timelines is limited due to capacity limits of the timeline. It is also impacted by visibility restriction due to a finite screen size and by the scrolling habits of users. Sources of content that wish to be visible on a given timeline have to take into account other flows of content that compete over space in the same timeline. Thus they have incentives to resend copy of their contents in order to make their opinions visible. In absence of such actions, the average number of their posts seen in the Timeline may be very small. The choice of rates of content can be formulated as a noncooperative game where a source's utility is the difference between some visibility utility and some cost for sending and creating contents.

In a previous paper, the authors [START_REF] Altman | Competition Over Timeline in Social Networks[END_REF] studied a similar game theory model under the following assumptions:

• The arrival process of contents is Poisson distributed;

• The visibility utility of a source is assumed to be the indicator of having at least one content of that source in the timeline.

In this paper we abandon the restrictive Poisson assumption. We allow general dependence of times between contents arrival. We shall only require that the arrival time of content for each different source i is a stationary ergodic point process, and that their superposition is a simple process (at a given time there may only be a single content arrival).

The assumption of stationary behaviour is well suited to a subscriber's wall. Indeed, there is no reason to expect that the user will access to his timeline at times with a distribution, which is not time-uniform. We assume that a user does not observe the system at a stationary regime but rather at instants of arrivals of new contents. This hypothesis is coming from the fact that we assume that users receive notifications whenever a new content arrives. In the theory of point process, the system observed at arrival instants of points in a stationary ergodic point process corresponds to the so-called Palm probability.

Another difference with respect to [START_REF] Altman | Competition Over Timeline in Social Networks[END_REF] is that we shall consider in this paper the total number of contents at arrival times of contents as the visibility payoff.

Our first result of this paper is to show that the problem can be reduced to the well known generalized Kelly mechanism, which allows us to use many existing results to characterize the equilibrium. Moreover we derive the equilibrium in explicit form and thus obtain new results for the generalized Kelly mechanism. We then extend the model to study rerouting of posts, which corresponds to sharing content.

Related Work

We briefly discuss the generalised Kelly mechanism [START_REF] Vazirani | Algorithmic game theory[END_REF]. Consider a resource K that is to be shared among J strategic users. A user proposes a bid λ j ∈ [0, λ]. Then the resource is shared proportionally to the bids, so that the amount of resource that player j receives is K λ j i λ i and each player j has a linear cost γ for his bid. Finally his objective function is:

U j (K λ j i λ i ) -γλ j ,
where U j (.) is a function that measures the satisfaction that receive the player j when he has K λ j i λ i of the resource. A study of this game was proposed in [START_REF] Tullock | Efficient rent seeking[END_REF]. Some networking applications of this game can be found in [START_REF] Ma | On resource management for cloud users: A generalized kelly mechanism approach[END_REF] and [START_REF] Ramesh | A game theoretic view of efficiency loss in resource allocation[END_REF].

A similar analytical study works like ours has been recently proposed in [START_REF] Altman | Competition Over Timeline in Social Networks[END_REF]. In this paper, the authors propose to study the competition over a timeline between several information providers (sources). There are several differences with our paper. Firstly, they consider Poisson arrival rates whereas, in our model, we consider a more general point process. Secondly, they define the visibility over a news feed as the proportion of time one message is visible on the timeline. Then, they consider in fact the case that only one message of each source can appear in the timeline. In our analysis we also generalize this point by considering that several messages, from the same source, can be visible on the timeline. Finally, we consider the propagation effect between news feed, which is not taking into account in [START_REF] Altman | Competition Over Timeline in Social Networks[END_REF].

The game formulation

Consider a single timeline of size K. It can contain up to K posts. Consider an arrival of content to that timeline which is a superposition of J stationary ergodic point processes N j (j = 1, . . . , N ) assumed to be compatible with the flow θ t , each of which has a finite and non-null intensity λ j . We call by N their superposition. Each source j is associated to the point process of intensity λ j . We use the definition of [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF] for the Palm probability which is the probability that a message is seen at an arrival instant. Let P o N (A) be the Palm probability corresponding to the point process N , and we define similarly the Palm process that corresponds to point process N j . It can be shown in page 37 or [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF] that

P o N (N j ({0} = 1)) = λ j i λ i (1)
This formula states that when an arrival occurs, the probability that it is from source j is

λ j i λ i .
We consider the dissemination utility for player i to be the sum of messages originating from it in the timeline at an arbitrary arrival instant of a content. This event can be written as

K-1 j=0 1(N o j T -j = 1)
whose expectation with respect to the Palm probability measure is given by

K-1 j=0 P o N (N j T -j = 1) = K λ j i λ i
due to (1.3.17) p 25 in [START_REF] Baccelli | Elements of queueing theory: Palm Martingale calculus and stochastic recurrences[END_REF] and eq [START_REF] Altman | Game theoretic approaches for studying competition over popularity and over advertisement space in social networks[END_REF]. We assume that each source j has a cost for sending messages, which depends on his rate λ j . This cost is expressed by γ j λ j . Finally, source's j objective function is defined by the difference between the average number of messages and the cost associated to the creation of contents:

U j (λ j , λ -j ) = K λ j i λ i -γ j λ j , (2)
where λ -j = {λ 1 , . . . , λ j-1 , λ j+1 , . . . , λ J } is the strategy vector of all the other sources.

Competition between sources

We are studying the visibility competition between the sources. We are studying a normal form non-cooperative game with J sources. The utility function of each source depends of its own rate and the one of other sources.

The utility function of source j is given in equation ( 2).

Definition 1. The decision vector λ * = (λ * 1 , . . . , λ * J ) is a Nash equilibrium if for all j ∈ {1, . . . , J},

U j (λ * j , λ * -j ) = max λ j ∈[0,λ] U j (λ j , λ * -j ),
where λ is the maximum source rate.

In other words, a source has not an interest in changing his rate unilaterally at the Nash equilibrium. This game was already originally proposed in [START_REF] Tullock | Efficient rent seeking[END_REF] but never applied in the context of Online Social Networks. It has been also used to study competition over popularity in [START_REF] Altman | Game theoretic approaches for studying competition over popularity and over advertisement space in social networks[END_REF], but without considering a finite number of messages of interest in a Timeline. This game has several important properties. First, at the Nash equilibrium, at least one source send messages in the timeline as it is stated in the next proposition.

Proposition 2. The decision vector (0, . . . , 0) is not a Nash equilibrium.

Next in the following theorem, we generalize the results obtained in [START_REF] Altman | Game theoretic approaches for studying competition over popularity and over advertisement space in social networks[END_REF]. In this theorem, we obtain an explicit solution of the unique equilibrium of the game, when each source sends messages in the timeline.

Theorem 3. If

i γ i -max j {γ j }(J -1) > 0, the Nash equilibrium is uniquely defined by:

∀j, λ * j = (J -1)K i γ i (1 - φγ j K ) > 0.
with φ = (J-1)K i γ i . In the previous theorem, there is a condition over the cost γ j for all j such that λ j > 0 for all j. If this condition is not satisfied, some rates can be equal to 0 at equilibrium for some sources. This result is described in the next theorem. Given the explicit expression of the Nash equilibrium of the game, we are able to understand better the behaviour of the equilibrium when the parameters of the game change. Thus from theorem 3, it can be noticed that for all source j, the rate λ * j is strictly decreasing in γ j and linear increasing in K. Thus we can now compute a general form of the equilibrium and not just restrict our self to an equilibrium where i γ i -max i {γ j }(J -1) > 0, in other word when all λ * j are positive. Assume we rearrange the sources, from smallest γ j to largest, then

γ 1 > γ 2 > . . . > γ J . Theorem 4. Assume that j is such that j = max{j | (J-1)K J i=j γ i (1 - (J-1)γ j J i=j γ i ) ≤ 0}
i.e j is the source with largest γ j who doesn't send any messages in the timeline. In this case the Nash Equilibrium is uniquely defined by:

λ * j = (J -1)K i>j γ i (1 - (J -1)γ j J i>j γ i ) ∀j > j λ * j = 0 ∀j ≤ j where J = J-| {1, . . . , j } |.
In order to study the inefficiency of the Nash Equilibrium we study the price of anarchy [START_REF] Nisan | Algorithmic game theory[END_REF]. We adopt the following assumption, for each j, λ j ∈ [λ, λ], with λ > 0 . First we define the social welfare of the game:

V (λ 1 , . . . , λ J ) = i [K λ i i λ i -γ j λ i ] = K - i γ j λ i
We maximize the social welfare:

(λ 1 , . . . , λ J ) = argmax V (λ 1 , . . . , λ J ) = (λ, . . . , λ)
. The price of anarchy is define by:

P oA = V (λ 1 , . . . , λ J ) j U i (λ * j , λ -j * )
Theorem 5. If i γ i -max j {γ j }(J -1) > 0 then the price of anarchy is:

P oA = ( i γ i )(K -λ i γ i ) j ( i γ i -(J -1)γ j )( i γ i -γ j (J -1)K)
.

Propagation effect within 2 timelines

We consider a generalization of our model to take into account the propagation effect between several timelines. In fact, some messages posted on one timeline can be relayed (propagated) to another timeline. A message which is of interest to one community can be of interest to another community for which the center of interest is closed. We are interested in one particular question for the sources:

Knowing the propagation effect between timelines , how will each source decide to send his information flow to the timelines, directly or indirectly using the propagation effect ? Indeed, we have assumed that posting a message on a timeline has a cost for the source, whereas if a message is relayed, it is free for the source. Thus, in this section, we study this topological effect on the non-cooperative game between the sources. In order to get interesting closed-form results, we reduce the complexity of the game analysis by considering two timelines, L = 2, but our results will be extended to more complex topologies in future works. First, in the next sub-section, we define the point process model the characterizes the arrival of the messages for the different timelines.

Point process model

Now assume that you have 2 timelines l ∈ {1, 2}. Each source j ∈ {1, . . . , J} can send messages to each timeline 1 and 2. It means that each source j controls 2 point processes, which represent flow of messages of source j in each timeline, with intensity λ l j = E[N l j ((0, 1])] for l = 1, 2. For each timeline l, any new message will be copied and posted, with a probability p ll , to the other timeline l = l. This process defines a new point process in each timeline l with intensity p l l j λ l j . If we add all the point processes for a timeline l (direct and indirect messages), this new point process N l has an intensity i (λ l i + p l l λ l i ). We assume that K l is the number of messages that subscribers of timeline l take care of it. We denote by T l n the arrival time of the n th message in the timeline l. Let N l j ((T n , T n+K l ]) the number of messages from source j in timeline l, as long as the n th message is visible on the timeline, i.e. is part of the first K l messages. Source j wants to maximize Nj that is the expected number of messages from him that are visible on all the timelines:

Nj = lim N →+∞ l 1 N N n=1 N l j ((T l n , T l n+K l ]).
By applying the similar analysis as in section for IV, we obtain the following closed-form expression:

Nj = K 1 λ 1 j + p 21 λ 2 j i (λ 1 i + p 21 λ 2 i ) + K 2 λ 2 j + p 12 λ 1 j i (λ 2 i + p 12 λ 1 i )
.

Let γ l j be the cost from source j for sending one message to timeline l. We denote the control vector of source j by λ j = {λ 1 j , λ 2 j }. The objective function of each source j is:

U j (λ j , λ -j ) = K 1 λ 1 j + p 21 λ 2 j i (λ 1 i + p 21 λ 2 i ) + K 2 λ 2 j + p 12 λ 1 j i (λ 2 i + p 12 λ 1 i ) -(γ 1 j λ 1 j + γ 2 j λ 2 j ). (3)

Game analysis

For the analysis of (3), we consider another variable for the controls of the players, and then we define a new game where the objectif function for each player j is: max

Λ 1 j ,Λ 2 j U j (Λ j , Λ -j ) = K 1 Λ 1 j i Λ 1 i + K 2 Λ 2 j i Λ 2 i -Γ 1 j Λ 1 j -Γ 2 j Λ 2 j (4) with (Γ 1 j , Γ 2 j ) = ( γ 1 j -γ 2 j p 21 1-p 21 p 12 , γ 2 j -γ 1 j p 12
1-p 21 p 12 ) for all j. We denote by (Λ 1 * j , Λ 2 * j ) the flow sent by source j at equilibrium. We define a set A:

A := (a, b) ∈ [0, 1) 2 | ∀j, Λ 2 * j b -Λ 1 * j > 0 or Λ 1 * j a -Λ 2 * j > 0
In the first proposition we give an other form of (3) in order to compute the equilibrium.

Proposition 6. If (p 12 , p 21 ) ∈ A (H1), then the game (3) is equivalent to the game (4) with

(Λ 1 j , Λ 2 j ) = (λ 1 j + p 21 λ 2 j , λ 2 j + p 12 λ 1 j ).
Here, we define some new hypothesis.

(H2)

{Γ 1 i } i are such that i Γ 1 i -max i {Γ 1 i }(J -1) > 0 {Γ 2 i } i are such that i Γ 2 i -max i {Γ 2 i }(J -1)
> 0 The next theorem gives us an explicit form of the equilibrium of the game when all sources send messages.

Theorem 7. Assume that (H1) and (H2) are satisfied. In this case the Nash Equilibrium is uniquely defined by:

(λ 1 * j , λ 2 * j ) = (a 1 j , 0) if Λ 1 * j < Λ 2 * j (0, a 2 j ) if Λ 1 * j > Λ 2 * j with a 1 j = p 21 (J-1)K 2 i γ 2 i -γ 1 i p 12 (1 - (J-1)(γ 2 j -γ 1 j p 12 ) i γ 2 i -γ 1 i p 12 ) -(J-1)K 1 i γ 1 i -γ 2 i p 21 (1 - (J-1)(γ 1 j -γ 2 j p 21 ) i γ 1 i -γ 2 i p 21 ) a 2 j = p 12 (J-1)K 1 i γ 1 i -γ 2 i p 21 (1 - (J-1)(γ 1 j -γ 2 j p 21 ) i γ 1 i -γ 2 i p 21 ) -(J-1)K 2 i γ 2 i -γ 1 i p 12 (1 - (J-1)(γ 2 j -γ 1 j p 12 ) i γ 2 i -γ 1 i p 12 ) Λ 1 * j = (J-1)K 2 (p 12 p 21 -1) i γ 2 i -γ 1 i p 12
(1 -

(J-1)(γ 2 j -γ 1 j p 12 ) i γ 2 i -γ 1 i p 12 ) Λ 2 * j = (J-1)K 1 (p 12 p 21 -1) i γ 1 i -γ 2 i p 21 ( 1
- (J-1)(γ 1 j -γ 2 j p 21 ) i γ 1 i -γ 2 i p 21 )
The previous theorem shows that under condition (H1) and (H2) all the sources send messages. And another important result is that each source just sends messages on one timeline. In the next section we study the effect of the parameter on the equilibrium.

Numerical Illustrations

We make our numerical illustrations with γ 1 j = γ 2 j = γ for all j. We start fig. 1. at p 12 = 0.75 because with our parameters, when p 12 < 0.75 the condition (H1) is not verified. We see that λ 1 * j increasing in p 12 and p 21 . In fig. 2. we can observe that sources (because symmetric sources have a symmetric behavior) prefer to send messages in timeline 2 when p 12 < 0.4 and p 21 = 0.45 or when p 12 < 0.7 and p 21 = 0.7. And they prefer to send in

Conclusions

In this paper, we model a competition between information sources that use a social network timeline in order to disseminate and to maximize its visibility, particularly the visibility of its messages inside a timeline. We take into consideration that the messages in a timeline are of interest to the subscribers only if the messages are in the first K positions (more K recent). We have modeled the competition between sources using a noncooperative game and we have proved that this non-cooperative game is a standard game well known in telecommunication which has been studied in different contexts (resource sharing, communication, etc). We have also considered an interesting feature taking into account several timelines and the possibility for a message to be copied and relayed from one timeline to another. Then we have generalized our result and prove the existence of the equilibrium. In future works we plan to analyze more complex networks of timelines in which messages can be propagated between several timelines depending on the behavior of the timeline subscribers.

give us the best reply of the player j.

∂U j λ j (λ j , λ -j ) = 0 ⇔ ( 1 i λ i - λ j ( i λ i ) 2 ) = γ j K ⇔ ( i λ i -λ j ( i λ i ) 2 ) = γ j K . ( * )
In the last equation, by taking the sum over j, we have:

J( i λ i ) -1 K i γ i ( i λ i ) 2 = ( i λ i ) (5)
⇔ (J -1)( i λ i ) -1 K i γ i ( i λ i ) 2 = 0 (6)
⇔ ( i λ i )((J -1) -1 K i γ i ( i λ i )) = 0 (7)
⇔ i λ i = (J-1)K i γ i . (8)
The closed form expression is obtained by using ( 8) and ( * ). The unicity comes from that λ * j is strictly defined by i λ * i and i λ * i is unique.

Proof of theorem 4

We only have to prove that if j plays 0 at equilibrium, all the other sources j < j play also 0. This fact is clear because λ * j is decreasing in γ j .

Proof of theorem 5

We have to compute U j (λ * j , λ -j * ).

U j (λ * j , λ -j * ) = (1 - (J-1)γ j i γ i ) -γ j ( (J-1)K i γ i )(1 - (J-1)γ j i γ i ) = (1 - (J-1)γ j i γ i )(1 -γ j (J-1)K i γ i ) = 1 i γ i ( i γ i -(J -1)γ j )( i γ i -γ j (J -
1)K) And we can deduce the form of the price of anarchy which is:

P oA = ( i γ i )(K -λ i γ i ) j ( i γ i -(J -1)γ j )( i γ i -γ j (J -1)K)
.

Proof of proposition 6

We are going to prove this theorem in two steps.

Step 1: let us look at the limit of the average number of messages in the two timelines for player j (i.e the first part of the objective function of the player j). We can notice that

K 1 λ 1 j + p 21 λ 2 j i (λ 1 i + p 21 λ 2 i ) + K 2 λ 2 j + p 12 λ 1 j i (λ 2 i + p 12 λ 1 i ) = K 1 Λ 1 j i Λ 1 i + K 2 Λ 2 j i Λ 2 i
with Λ 1 j = λ 1 j + p 21 λ 2 j and Λ 2 j = λ 2 j + p 12 λ 1 j for all j.

Step 2: We want to find Γ 1 j and Γ 2 j such that Γ

1 j Λ 1 j + Γ 2 j Λ 2 j = γ 1 j λ 1 j + γ 2 j λ 2 j . Firstly (λ 1 j , λ 2 j ) = ( Λ 2 j p 21 -Λ 1 j p 12 p 21 -1 , Λ 1 j p 12 -Λ 2 j p 12 p 21 -1 ) is solution of Λ 1 j = λ 1 j + p 21 λ 2 j Λ 2 j = λ 2 j + p 12 λ 1 j .
Secondly, with the above result we can calculate Γ 1 j and Γ 2 j . In fact,

γ 1 j λ 1 j + γ 2 j λ 2 j = 1 p 12 p 21 -1 Λ 1 j (γ 2 j p 12 -γ 1 j ) + Λ 2 j (γ 1 j p 21 -γ 2 j ) . Thus we can conclude that (Γ 1 j , Γ 2 j ) = ( γ 2 j p 12 -γ 1 j p 12 p 21 -1 , γ 1 j p 21 -γ 2 j p 12 p 21 - 1
).

Proof of theorem 7

The proof follows two steps. Firstly, we calculate the equilibrium of ( 4) and secondly we make the change of variables in order to have the value of the direct flow of each sources at equilibrium.

Step 1: In order to find the equilibrium of (4) we have to compute the best reply. We do it just for one Λ 1 j . It works the same way as Λ 1 j for all the other j and timeline 2.

∂U j ∂Λ 1 j (Λ j , Λ -j ) = 0 K 1 ( 1 i Λ 1 i - Λ 1 j ( i Λ 1 i ) 2 ) -Γ 1 j = 0 (9)
We remark that the best reply just depend of {λ 1 j } j and it is equivalent to the best reply of the game (2) where the objectif function for each player j is: max

Λ 1 j U j (Λ j , Λ -j ) = K 1 Λ 1 j i Λ 1 i -Γ 1 j Λ 1 j
This is why by using theorem 3 we have ∀j, (Λ 1 * j , Λ 2 * j ) = (

(J -1)K 1 i Γ 1 i (1 - Φ 1 Γ 1 j K 1 ), (J -1)K 2 i Γ 2 i (1 - Φ 2 Γ 2 j K 2
)) [START_REF] Tucker | Social advertising[END_REF] with (Φ 1 , Φ 2 ) = ( (J-1)

K 1 i Γ 1 i , (J-1)K 2 i Γ 2 i
).

Step 2: In the proof of proposition 10 we have proved that (λ 1 j , λ 2 j ) = ( ). This is why by using this fact we have:

λ 1 j = p 21 (J-1)K 2 i γ 2 i -γ 1 i p 12 (1 - (J-1)(γ 2 j -γ 1 j p 12 ) i γ 2 i -γ 1 i p 12 ) -(J-1)K 1 i γ 1 i -γ 2 i p 21 (1 - (J-1)(γ 1 j -γ 2 j p 21 ) i γ 1 i -γ 2 i p 21 ) λ 2 j = p 12 (J-1)K 1 i γ 1 i -γ 2 i p 21 (1 - (J-1)(γ 1 j -γ 2 j p 21 ) i γ 1 i -γ 2 i p 21 ) -(J-1)K 2 i γ 2 i -γ 1 i p 12 (1 - (J-1)(γ 2 j -γ 1 j p 12 ) i γ 2 i -γ 1 i p 12 )
Now in (H1) we have assume that Λ 2 * j p 21 -Λ 1 * j > 0 or Λ 1 * j p 12 -Λ 2 * j > 0. We notice that if Λ 2 * j > Λ 1 * j then Λ 1 * j p 12 -Λ 2 * j < 0 which implies that λ 2 * j = 0. And if Λ 2 * j < Λ 1 * j then Λ 2 * j p 21 -Λ 1 * j < 0 which implies that λ 1 * j = 0. This is why a player sends messages just on one wall.

Figure 1 :

 1 Figure 1: Behavior of λ 1 * j under p 12 with γ = 10, K 1 = K 2 = 10, J = 100, p 21 = 0.45 for -. and p 21 = 0.7 for . . .

Figure 2 :

 2 Figure 2: Behavior of Λ 1 * j -Λ 2 * j under p 12 with γ = 10, K 1 = K 2 = 10, J = 100, p 21 = 0.45 for -. and p 21 = 0.7 for . . .

Λ 2 j p 21 -Λ 1 jp 1 j p 12 -Λ 2 jp

 11122 12 p 21 -1 , Λ 12 p 21 -1

APPENDIX Proof of proposition 2

If all sources {1, • • • , J} -{j} play 0, j can always find a very small λ * j such that λ * j < K γ j . Then, we have that U j (λ * j , 0) > 0 and the decision vector (0, . . . , 0) is not a Nash equilibrium.

Proof of theorem 3

The first order condition:

∂U j λ j (λ j , λ -j ) = 0