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 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 PTDC/MAT/121185/2010.


Nicolas Ribière-Tharaud

Aronne Casagranda

François Marc Lambert

François Jouvie

F Jouvie Outline

Numerical and experimental assessment of a phase retrieval technique applied to planar near field distribution for wide band application

! Sampling parameter : δ x ! Far field from using the reconstructed phase

• δ x = λ 0 /6 Δ n cmplx = 7.94 % • δ x = λ 0 /3 Δ n cmplx = 6.07 % • δ x = λ 0 /2 Δ n cmplx = 29.53 % θ(°) (ϕ = 0 °)# ϕ(°) (θ = 90 °)# |E θ | (dB) |E θ | (

•=•

  Experimental antenna far field pattern characterization • Wide frequency range • Pulsed sources ! Measurement constraints • Planar near field test setup • Magnitude only measurements ! APPROACH ! Implementation of a phase retrieval algorithm ! Assumptions on sources • CW • High directivity Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) 2 THE PHASE RETRIVAL TECHNIQUE ! HISTORICAL BACKGROUND ! Gerschberg-Saxton [1972] : electron microscopy (1 plane) ! « Misell » variant [1973] : electron microscopy (2 defocused planes) ! Anderson & Ali [1984] : microwave applications ! PRINCIPLE ! Algorithm initialization : • Electric field magnitude known in two planes P1 & P2 in front of the source • Phase is set to zero on the first plane at iteration #Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) Fields calculated in magnitude and phase using Feko • distances : 11 values [x 1 = 2λ, …, x 11 = 1 m] • en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) 5 NUMERICAL ASSESMENTS (3/8) ! Planes positions parameter : x 1 and x 2 ! Far field from using the reconstructed phase • Worst case : Δ n cmplx = 222.54 % • Best case : Δ n cmplx = 7.94 % |E θ | (dB) θ(°) (ϕ = 0 °)# |E θ | (dB) ϕ(°) (θ = 90 °)# Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) 6 NUMERICAL ASSESMENTS (4/8)

!

  dB) Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) ) (ϕ = 0 °)# ϕ(°) (θ = 90 °)# |E θ | (dB) |E θ | (dB) Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) Experimental validation : planar near field measurement setup ! AUT = Horn, Probe = Dipole, open end waveguide ! Frequencies = 2 GHz, 8 GHz and 18 GHz ! Distance AUT/Probe : 2λ, 3λ, … 1m Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) Département de Recherche en Électromagnétisme de Supélec et du L2S (UMR8506 / CNRS -Supélec -UPS) 14 CONCLUSION & FUTURE WORKS ! Conclusions on the phase reconstruction algorithm ! Validation on numerical and experimental data ! Parametric study for optimal use ! Wide band validation ! Future tasks ! Adding information to increase performances ! Expanding the parametric study ! Towards pulsed sources …


Mena Lafkioui

DIALECTOMETRY ANALYSES OF BERBER LEXIS 1

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

phonological and therefore lack an interpretation basis that is connected to other linguistic dimensions. 5The computational dialectometry methods are numerous and are currently considered most adequate for reasons I will explain later on in this article. The foundations of digital automated dialectometry were established by [START_REF] Séguy | La dialectométrie dans l'Atlas linguistique de la Gascogne[END_REF] with his analytical method to calculate the linguistic differences between varieties of Gascogne. The comparison is based on an algorithm which classifies data as identical or non-identical. The sum of the measured distances between two varieties matches their linguistic distance. The visualisation of the classification analysis is conducted through lines of various types (bold/non bold, dotted/non dotted, etc.), which divide the region according to the linguistic differences of the varieties. As a counterpart, Goeble (1982Goeble ( , 1993) ) has calculated the similarities between varieties from Italy and Southern Switzerland. Even though the results of the calculation of Séguy and Goeble have the merit of being objective, they lack refinement because their technique excludes distance graduation.

The main computational methods based on the frequency of linguistic variants are the "Corpus Frequency Method" [START_REF] Hoppenbrouwers | De feature frequentiemethode en de classificatie van Nederlandse dialecten. -TABU[END_REF], 2001) and the "Frequency per Word Method" [START_REF] Nerbonne | Computationale vergelijking en classificatie van dialecten[END_REF], 2001). The basic principle of the first approach is that the degree of difference/similarity between two varieties is derived from comparing the frequency of the marked linguistic features of their variants. The problem in this approach is that the entity "word" is not considered as a linguistic unit. However, this obstacle is removed by the second approach which assigns to words the status of "units" functioning as such. Nevertheless, the two classification tools do not take into account the order of the phonic units in the sequence.

The "Levenshtein distance" (Lv), on the other hand, allows incorporating the parameter of sequential ordering of phonic units in the classification, which makes it more appropriate than other digital/numerical methods. This tool has been introduced in dialectometry by [START_REF] Kessler | Computational dialectology in Irish Gaelic[END_REF], who has applied it to a corpus of Irish Gaelic. The Levenshtein distance measure corresponds to the numerical value of the lowest cost of operations (insertions, deletions and substitutions) needed to convert a string of characters into another [START_REF] Kruskal | An overview of sequence comparison[END_REF]. One of the most employed techniques of comparison is the "phone string comparison" in which all operations have the same cost, regardless of the degree of affinity between the phonic units: the pair [t, d] has the same cost as the pair [u, t] and [u, u:]. Yet, with the technique of "feature string comparison" phonetic features of phonic units can be compared: the cost of the pairs [u, t] and [u, u:] is not equal because the phonetic affinity between the phonic units of [u, u:] is greater than that of [u, t].

Dialectometry analyses of Berber lexis

Among the different existing dialectometry approaches, I prefer the computational methods because they allow handling large data corpora with certain ease, while ensuring the accuracy and consistency of the analyses. These aims can be achieved thanks to the fact that -Distances and frequencies are measured automatically.

-Data are classified digitally.

-Mapping can be assisted by the computer.

-Statistical analyses can be made and displayed automatically.

The dialectometry analyses that I present in this article were performed with the free software of Kleiweg (RuG/L04). 6 In order to complete a displayed dialectometry analysis, all the procedural steps summarised below are indispensable (Lafkioui forthcoming 1): The lexical data which are compared and classified in this study are collected from the Atlas linguistique des variétés berbères du Rif [START_REF] Lafkioui | Atlas linguistique des variétés berbères du Rif[END_REF] or ALR. The digital data corpus consists of sixty-two lexemes regarding the human body (maps 295 to 315), kinship (maps 316 to 321 cards), animals (maps 322 to 327), colours (maps 328 and 329), numbers (maps 330 to 332), besides a subset of various nouns and verbs (maps 333 to 356). Of these lexemes, eleven have only one variant per variety; all fiftyone other lexemes display the co-occurrence of multiple variants for each lexeme.

Due to the completion of the automated ALR, the data obtained from it are already in digital format, which has avoided a great task of digitising. However, an adaptive conversion to the software RuG/L04 (Kleiweg) was necessary. The ALR also offers a precise digital map of the Rif region (see Figure 1), which is essential to the visualisation of the dialectometry analyses, except for the dendrogram. One hundred forty-one georeferenced points -belonging to thirty-two Rif tribes -were selected from a group of four hundred fifty-two localities in the Rif according to their degree of linguistic variation [START_REF] Lafkioui | Atlas linguistique des variétés berbères du Rif[END_REF]). 7 7 The survey points were selected on the basis of the principle of equidistance dividing the inquiry field into several grids to which were assigned points that could match with localities on the field. The greater the variation was, the more the grids were reduced. The four hundred fifty-two locations selected for this research were for the most part chosen so that they could, a priori, indicate linguistic borders. This selection mainly stemmed from the scientific and empirical knowledge of the investigator on the different varieties spoken in the Rif area.

Data matrix of Rif-Berber lexis

The data matrix is composed of digital lexical excerpts from the ALR [START_REF] Lafkioui | Atlas linguistique des variétés berbères du Rif[END_REF]) converted following the format of the software RuG/L04 (Kleiweg). Here below, a small sample in digital format of the ALR (Mapinfo Professional format; Table 2) and in text format of RuG/L04 (Table 3) are given:

Distance matrix of Rif-Berber lexis

This section contrasts the three most employed digital comparison techniques: the Binary distance (Hamming algorithm), the Gewichteter Identitätswert distance (Weighted identity value), and the Levenshtein distance. I will apply these techniques on the Rif-Berber lexical corpus to test their validity and to select the most appropriate to Berber. Each distance measuring allows acquiring precise numerical values derived from the linguistic comparison between the varieties of the Rif area. These values make up the distance matrices (symmetric matrices N x N, N= sum of varieties), whose configuration differs depending on the adopted algorithm.

Binary distance

The Binary distance (Bin) is used to classify lexical units as being identical or non-identical: comparison of type 0-1; 0= resemblance and 1= difference. Table 4 presents an excerpt from the Binary distance matrix of the lexeme "heel" (ALR, map 312):

Gewichteter Identitätswert distance

The Gewichteter Identitätswert distance (GIW) deviates from the Binary distance in that the frequency of the lexical variants is considered in the comparison: low-frequency variants weigh heavier than high-frequency variants. The distance obtained by this technique varies between 0 and 1; {0 ≤ d ≤ 1}. Table 5 presents an extract from the distance matrix of the lexeme "heel":

Levenshtein distance

The distance values resulting from a Levenshtein-based comparison -an algorithm taking into account the sequential order of phonic units composing lexemes -fluctuate between 0 and 1, {0 ≤ d ≤ 1}, as shown in the following excerpt: These values result from the selection of the least costly calculation to transform a lexical unit -as a string of phonic units -into another. Table 7 depicts the lowest cost of operations which allow modifying the string awrez (heel) into inerz (heel): The lowest cost of operations amending awrez into inerz is 3, which implies that the distance between these two lexemes is 3/5 (5 being the total number of features); consequently, the Levenshtein distance is 60%. These calculations are based on operations that cost 0.5 for an insertion or deletion and 1 for a substitution. Table 8 illustrates this calculation technique: "[…] a technique that, using a table of differences, tries to position a set of elements into some space, such that the relative distances in that space between all elements corresponds as close as possible to those in the table of differences." (Kleiweg, RuG/L04).

Visualisation of dialectometry analyses of Rif-Berber lexis

Classification by clustering (CA) necessarily uses a dendrogram for its display. A dendrogram is a complex ranking structure, usually in colour, whose branches represent the linguistic varieties. It can be matched with a digital map, resulting in a geolinguistic map that shows the distribution of linguistic varieties depending on the linguistic differences and the selected classification criteria. In contrast, analyses by Multidimensional Scaling (MDS) directly offer maps on which the relative linguistic variation is gradually represented by different colours.

Visualisation and interpretation of the CA analyses

The hierarchical structure of the dendrogram and the corresponding varieties' distribution on the Rif map vary significantly depending on the distance algorithm (Bin, GIW or Lv) applied. Thus, the outcome of the Binary distance comparison (Figure 2) is a configuration consisting of seven main groups, clustered into two subgroups: the minor subgroup containing groups 6 and 7 and the major subgroup containing groups 1 to 5; the distance between these two subgroups is 16.17. This distance value indicates a relatively high linguistic boundary after group 7, which is delimited on the right by the varieties of Ayt Weryaghel and Ayt 'Ammart.

The major subgroup shows a rather balanced subdivision (d = 9.34) between groups 4 and 5 (variety of Targist included) and groups 1 to 3, which have also been subdivided. The second important linguistic border therefore coincides with the bordering varieties of groups 4 (Igzennayen) and 5 (Ayt S'id and Ayt Tuzin). The classification based on the GIW algorithm diverges considerably from that based on the Bin algorithm, because it leads to a set of five clusters (Figure 3), of which cluster 1 includes the sub-clusters 1, 2 and 4 of the Bin classification (Figure 2). However, the main linguistic boundary detected through GIW -boundary drawn after the varieties of cluster 5 -is identical to the one that emerged from the Bin dendrogram. Although, the distance between the two major sub-clusters is lower for GIW (d GIW = 10.87) than for Bin (d Bin = 16.17). This difference can be explained by the integration of the frequency parameter in the comparison.

Visualisation and interpretation of the MDS analyses

The MDS technique has the major advantage of ensuring objectivity and accuracy during the analysis stage of the materials because it excludes any external parametering. For example, the number of clusters cannot be changed because the analysis system provides it automatically. Each variety has its own colour. The colour contrasts are used to interpret the compared linguistic data: a colour continuity points to a perfect correlation between lexemes, while a colour mosaic reveals a low correlation between them.

The Rif region is divided into 7 major areas, regardless of the distance measuring applied (Figure 6). The distribution of the varieties on the MDS maps is almost similar to Bin and GIW; only a few minor differences in shades of certain colours were observed. The MDS Lv map closely resembles the other two; the only significant distinction observed is the emergence of a small subdivision inside the group of Western varieties.

Contrastive results

Because of its accuracy, the MDS method is most appropriate for dialectometry analysis of Berber lexis. Accordingly, it forms a yardstick against which other dialectometry methods can be contrasted. Among the Cluster Analyses classifications (CA), CA Bin and CA Lv join best the distribution maps displayed by MDS (7 groups). Moreover, the CA Lv classification shows a further refinement because it takes into account the phonic variation of the units as much as their arrangement in the lexemes. However, any analysis based on Lv distance (CA as well as MDS) ignores the existence of the hierarchy between the phonic units (phonetic units= phonemic units), unless various weights are granted to them through a specific parametring. This method implies the construction of a phonological system within the software, involving a time and energyconsuming effort that is much too expensive compared to its profits. The Cluster Analysis classification has the benefit of precisely indicating significant linguistic boundaries. The CAC maps (Composite Cluster map; Figures 7 and8) designate these boundaries by dark lines. Compared to the distinctive boundaries drawn by the dendrograms and CA maps of Figures 2 to 5, the principal linguistic delimitation of the CAC map of Figure 7 is drawn further to the West. It is important, nevertheless, to note that the CAC maps do not seem suitable to display the classification of Rif-Berber lexis because of the difficulty of interpreting the data, due to their rather chaotic representation (Figure 8).8

Figure 1 :

 1 Figure 1: Map of the georeferenced survey points of the Rif (Lafkioui 2007)

Figure 2 :

 2 Figure 2: Dendrogram vs. CA Map -Bin -All lexis

Figure 6 :

 6 Figure 6: MDS Lv map -All lexis

Figure

  Figure 7: CA Lv vs. MDS Lv vs. CAC Lv maps

Figure 8 :

 8 Figure 8: CAC -Bin vs. GIW vs. Lv maps -All lexis

Table 1 : General procedure of computational dialectometry analysis Step 1 Linguistic Atlas = georeferenced data source Step 2 Data Matrix Step 3 Distance Matrix Step 4 Analysis Step 5 Visualisation 2.1. The Linguistic Atlas of the Rif as a data source

 1

	The Rif is that region of North Morocco stretching from the Strait of
	Gibraltar in the West to the Algerian frontier in the East. The Rif-Berber
	varieties (Tarifit) belong to the northern Berber languages and thus are part
	of the large Afro-Asiatic language phylum. The Berber-speaking area of the
	Rif is delimited:
	-In the West, by the varieties of the Ktama tribe, (the so-called
	Senhaja varieties).
	-In the South, by the koinè of Gersif, which is the ultimate
	geographic point where Rif-Berber (Tarifit) is spoken before
	reaching the corridor of Taza.
	-In the East, by the varieties of Iznasen, which have spread to the
	regions of Arabic speaking varieties to the Morocco-Algerian
	border.

Table 2 : Data excerpt in digital format of ALR

 2

Table 3 : Data excerpt in text format of RuG/L04

 3

	: Asammer	: Wersan	: Bu'di
	-aqnin	-aqnin	-aqenni
	: Lmexzen	: A'raben	: Aghir Hmed
	-aqnin	-aqnin	-aqenni
	: Ssahel	: Azila	: Asammar
	-aqnin	-aqnin	-aqenni
	: Lqel'a	: Tamadda	: Ayt Hmed
	-aqnin	-aqnin	-aqenni
	: Ssaqya	: Azru n tili	: Sidi Bucetta
	-aqnin	-aqenni	-aqenni
	: Tarya	: Tizi	: Tazrut
	-aqnin	-aqenni	-aqenni
	: Mazuz	: Igzennayen	: Ufis
	-aqnin	-aqenni	-aqenni
	: Luta	: Iharunen	: Wad Mahkim
	-aqnin	-aqenni	-aqenni
	: Tizirt	: Ayt 'Azza	: L'ars
	-aqnin	-aqnenniy	-aqenni
	: Aghennuy	: Bni Budjay	: Tufist-Imuruten
	-aqnin	-aqenni	-aqenni

Table 4 : Excerpt from the Binary distance matrix of the lexeme "heel"

 4

Table 5 : Excerpt from the GIW distance matrix of the lexeme "heel"

 5

Table 6 : Excerpt from the Lv distance matrix of the lexeme "heel"

 6

Table 7 : Cost of operations allowing modification of awrez into inerz (heel)

 7

			a	w	r	e	z
		0	0.5	1	1.5	2	2.5
	i	0.5	1	1.5	2	2.5	3
	n	1	1.5	2	2.5	3	3.5
	e	1.5	2	2.5	3	2.5	3
	r	2	2.5	3	2.5	3	3.5
	z	2.5	3	3.5	3	3.5	3

Table 8 : Example of calculation of Lv distance for modifying awrez into inerz (heel)

 8

	Tamadda	a	w	r	e		z
	Tizi	i	n		e	r	z
	Lv Distance 1	1	0,5	0	0,5 0	3/5 * 100 = 60 %
	2.4.

Numerical dialectometry analyses of Rif-Berber lexis

  From the distance matrices, numerical comparative analyses of Berber lexis can be accomplished through two techniques: Cluster Analysis and Multidimensional scaling. The technique of Cluster Analysis (CA) consists of regrouping data by reducing the distance matrix by means of various algorithms. According to Kleiweg (RuG/L04), I have implemented the Ward algorithm (minimum variance), which is generally regarded as one of the most appropriate algorithms for this type of analysis. On the other hand, multidimensional scaling (MDS) is:

Although the frequency of the compared variants is taken into account (Kocks 1970), this approach does not seem to be the most appropriate(Heeringa 2004: 24-25).

Http://odur.let.rug.nl/~kleiweg/L04.

Kleiweg (RuG/L04) offers some alternatives to the Ward algorithm which seems to be causing visual disorder by which CA GIW maps are mostly affected.

The lexis classification obtained through Lv distance corresponds with an asymmetrical configuration of 7 clusters which are structured into 2 major clusters distanced from one another by 8.08 (Figure 4). The matching dendrogram shares the same linguistic main delimitation (between groups 6 and 3-4) with the other dendrograms. This observation is corroborated by the CA Lv maps displayed in Figure 5, of which the 2-cluster map clearly indicates the most distinctive linguistic boundary. It is important to note that the CA Lv map (Figure 4) shows a distribution of the varieties similar to the CA Bin distribution, even though the composition of their respective dendrogram is divergent.