NON-OVERLAPPING SCHWARZ ALGORITHMS FOR THE INCOMPRESSIBLE NAVIER–STOKES EQUATIONS WITH DDFV DISCRETIZATIONS

We propose and analyze non-overlapping Schwarz algorithms for the domain decomposition of the unsteady incompressible Navier–Stokes problem with Discrete Duality Finite Volume (DDFV) discretization. The design of suitable transmission conditions for the velocity and the pressure is a crucial issue. We establish the well-posedness of the method and the convergence of the iterative process, pointing out how the numerical fluxes influence the asymptotic problem which is intended to be a discretization of the Navier–Stokes equations on the entire computational domain. Finally we numerically illustrate the behavior and performances of the algorithm. We discuss on numerical grounds the impact of the parameters for several mesh geometries and we perform simulations of the flow past an obstacle with several domains. Mathematics Subject Classification. 65M08, 35Q30, 76D05. Received April 20, 2020. Accepted May 22, 2021.


Introduction
The aim of this paper is to develop a non-overlapping iterative Schwarz algorithm for the incompressible Navier-Stokes problem with DDFV schemes. The problem we are interested in reads where Ω is an open connected bounded polygonal domain of R 2 , f ∈ ( 2 (Ω)) 2 is a given force field, u ∈ ( ∞ (Ω)) 2 . The unknowns u : Ω × [0, ] → R 2 and : Ω × [0, ] → R are respectively the velocity and the pressure; (u, ) = 2 Re Du − Id stands for the stress tensor, and Re > 0 is the Reynolds number. Here and below, the strain rate tensor is defined by the symmetric part of the velocity gradient Du = 1 2 (∇u + ∇u). Non-overlapping Schwarz algorithms enter the class of domain decomposition methods, in which a domain is decomposed into smaller subdomains. The main advantage is that, contrarily to direct methods, decomposition methods are naturally parallel; in fact, subdomain problems are connected by some transmission conditions on the interface, but they are uncoupled by an iterative procedure. This makes those methods interesting for high performance computing perspectives. The classical Schwarz algorithm, proposed in 1870 by Schwarz [39] for the Laplace problem, is an iterative method that consists in transmitting the solution, or its normal derivative, from a subdomain to the others, in order to deal with complex domains. This method converges only if the subdomains overlap. Moreover, this convergence becomes slower as the overlap between the subdomains is smaller. In non-overlapping Schwarz algorithms, the subdomains intersect only on their interfaces and in order to obtain convergence, more elaborated transmission conditions should be defined on the interfaces. In 1990, P.-L. Lions [35] showed that, with Fourier-Robin transmission conditions, the Schwarz algorithm for the Laplace operator converges even without overlap. This method has been adapted to the discrete setting for many problems of isotropic diffusion [1,7,19], and for advection-diffusion-reaction problems [17,26]. For the Navier-Stokes problem, different approaches have been proposed, with different design of the interface conditions, that depend also on the discretization framework. In the spirit of Halpern and Schatzman [27], Blayo et al. [3] derives optimal transparent boundary conditions for the Stokes equation; these conditions are tested with finite difference methods. In the finite element setting, it is proposed in [37] a non-overlapping domain decomposition algorithm of Robin-Robin type for the discretized Oseen equations (i.e. linearized Navier-Stokes). In [40], still in the finite element setting, the authors build a Dirichlet-Neumann domain decomposition method for the nonlinear steady Navier-Stokes equations, under the hypothesis that the Reynolds number is sufficiently small and Girault et al. [22] studies a family of discontinuous Galerkin finite element methods for Stokes and Navier-Stokes problems on triangular meshes. The Inf-Sup condition, which is a crucial ingredient of the stability analysis of numerical methods for Incompressible Navier-Stokes equations, has to be adapted to the domain decomposition formulation, in particular to satisfy the incompressibility constraint, and it might depend on the Reynolds number, see [22,37]. Therefore, our objective is to decompose the domain Ω of problem (1.1) into smaller subdomains, to solve the incompressible Navier-Stokes problem on those subdomains by imposing some transmission conditions on the interfaces, and to recover by an iterative Schwarz algorithm the discrete solution of (1.1) on the entire computational domain Ω. Since we are interested in the unsteady problem, we shall apply this iterative algorithm at each time iteration. Moreover, we want the interface conditions to be local and we wish the method to remain free of any restrictive condition on the Reynolds number. We address these issues in the framework of finite volume methods, and more specifically by using Discrete Duality Finite Volume discretizations.
The introduction of the DDFV formalism dates back to [9,12,29,30], in order to approximate anisotropic diffusion problems on general meshes, including non-conformal and distorted meshes. Such schemes require unknowns on both the vertices and centers of primal control volumes and allow us to build two-dimensional discrete gradient and divergence operators that satisfy discrete duality relations analogous to the standard integration by parts formula. The DDFV scheme is extended in [2] to general linear and nonlinear elliptic problems with non homogeneous Dirichlet boundary conditions, including the case of anisotropic elliptic problems. Applying the DDFV method for Stokes and Navier-Stokes problems leads naturally to locate the unknowns of velocity and pressure in different points; the velocity unknowns are associated to the vertices and centers of primal control volumes, while the pressure unknowns are located on the edges of the mesh [5,10,11,24,32,33]. Hence, DDFV enters the class of staggered methods, reminiscent of the MAC scheme [28] constructed on Cartesian meshes for incompressible flows. The DDFV approach has, at least, two important advantages. First of all, it applies to very general meshes. It is useful, for instance, in the domain decomposition setting where the subdomains can be meshed separately and non-conformal edges appear on the interface, or simply if one wants to locally refine the mesh with cells adapted to complex geometries. Second of all, DDFV operators mimic at the discrete level the duality properties of the continuous differential operators, which leads to important properties for the numerical analysis of the schemes.
As a starting point of this study, we refer the reader to [20,26]: they both build a non-overlapping Schwarz algorithm in a finite volume framework with Fourier-like transmission conditions between subdomains, respectively for anisotropic diffusion with a DDFV discretization, see also [4], and for advection-diffusion-reaction in a TPFA discretization. The case of the Navier-Stokes equations (1.1) is more demanding, since it combines further difficulties: the vectorial nature of the unknowns, the non-linear convection terms and the incompressibility constraint. Of course, there are no explicit interface conditions, and one should construct suitable transmission conditions between the subdomains, which have the shape of Fourier-like conditions on the velocity and account for the constraint by involving the divergence of the velocity and the pressure. Let us split the computational domain Ω into two smaller subdomains Ω = Ω 1 ∪ Ω 2 . We denote by Γ the interface Ω 1 ∩ Ω 2 . The Schwarz algorithm defines a sequence of solutions (︀ u )︀ ∈N of the Navier-Stokes problem in Ω , with ∈ {1, 2}, endowed with the following two-fold transmission condition where ̸ = and ⃗ n is the outward normal to Ω . The former condition, which involves the parameters > 0, > 0, is inspired by the classical Fourier condition; it linearly combines the values of the unknown and the values of its derivative; here, also the convection is included. The latter, which depends only on , combines the divergence of the velocity with the pressure; it will be useful to satisfy the incompressibility constraint at the convergence of the algorithm. The first condition is comparable to the transmission conditions in [37]. However they need to justify a modified Inf-Sup condition which induces Reynolds-dependent stability constraints. This can be relaxed by imposing the new condition for the pressure on the interface. Once the transmission condition is fixed -which in practice will be solved in an approximated form -it remains to establish the convergence of the iterative process: as → ∞, one expects to recover the solution of a discrete version of the Navier-Stokes equations on the entire domain Ω. We shall see that the asymptotic scheme depends on the details of the numerical fluxes of the domain decomposition method. To analyze this issue, it is convenient to discuss general discretizations of the convection terms, inspired by Chainais-Hillairet and Droniou [8] and Halpern and Hubert [26].

Meshes
The DDFV method requires unknowns on vertices, centers and edges of control volumes; for this reason, it works on (three) staggered meshes. From an initial mesh, called the "primal mesh" (denoted with M ∪ M), we construct the "dual mesh" (denoted with M * ∪ M * ), centered on the vertices of the primal mesh, and the "diamond mesh" (denoted with D), centered on the edges of the primal mesh; see Figure 1. The union of the primal and dual meshes will be denoted by T.
More precisely, we consider a primal mesh M consisting of open disjoints polygons K such that ⋃︀ K∈MK =Ω. We denote M the set of edges of the primal mesh included in Ω, considered as degenerated primal cells. We associate to each K a point K , called center. For the volumes of the boundary, K is situated at the midpoint of the edge. Hp 2.1. All control volumes K are star-shaped with respect to K .
We build two others meshes: the dual mesh M * ∪ M * and the diamond mesh D. The dual mesh M * ∪ M * is made of dual cells K * and the diamond mesh D is made of diamond D. We set * the set of vertices of the primal mesh M and K * is an element of this set * . When K and L are neighboring volumes, that is if the measure of K ∩ L is positive, we define a diamond D as the quadrilateral (see Fig. 2) whose vertices are K , L , K * and L * where K * and L * are common vertices of K and L such that [ K * , L * ] ⊂ K ∩ L. To each diamond, we define its diagonals as a primal edge = K|L = [ K * , L * ] and a dual edge * = K * |L * = [ K , L ]. We denote a diamond D or D , * . (In the framework adopted here this is an interpretation and not a restrictive assumption; in particular, it is possible to consider non-conformal meshes. For instance, this is meaningful for the shadowed cell in Fig. 1-left since it is actually considered as a pentagon, not as a rectangle.) Let ℰ be the set of all primal edges and ℰ int = ℰ ∖ { ∈ ℰ such that ⊂ Ω}. The DDFV framework is free of "admissibility constraint", in particular we do not need to assume the orthogonality of the segment K , L with = K|L. Let ℰ * be the set of all dual edges. The diagonal and * intersect at the center of the diamond, D ∈ D and every diamond is star-shaped with respect to . As a consequence of this setting, all segments * belong to the physical domain Ω. We distinguish the diamonds on the interior and of the boundary: From the primal mesh, we build the associated dual mesh. A dual cell K * is associated to a vertex K * of the primal mesh. The dual cells are obtained by joining the centers of the primal control volumes that have K * as vertex. We distinguish interior dual mesh, for which K * does not belong to Ω, denoted by M * and the boundary dual mesh, for which K * belongs to Ω, denoted by M * . In what follows, we assume: Hp 2.2. All control volumes K * are star-shaped with respect to K * .
There are several possible constructions of the dual mesh, and it can happen that dual cells overlap. To avoid this inconvenience, we can either suppose that the diamonds are convex or consider the barycentric dual mesh, obtained by joining the centers K of the primal control volumes to the middle point of the edges that have K * as a vertex. Thanks to Hypothesis 2.1, barycentric dual cells have disjoint interiors. Throughout the paper, we restrict to the case where all diamond cells are convex.

Notations
The following notation will be used throughout the paper. The reader familiar with DDFV may skip this section.
For a volume V ∈ M ∪ M ∪ M * ∪ M * we define: For a diamond D , * whose vertices are ( K , K * , L , L * ), we denote: -D the center of the diamond D: D = ∩ * , -the length of the edge , -* the length of * , -D the measure of the diamond D , * , -d D the diameter of the diamond D , * , -D the angle between and * .
We introduce for every diamond two orthonormal basis (⃗ K * L * , ⃗ n K ) and (⃗ n * K * , ⃗ KL ), where: -⃗ n K the unit normal to going out from K, -⃗ K * L * the unit tangent vector to oriented from K * to L * , -⃗ n * K * the unit normal vector to * going out from K * , -⃗ KL the unit tangent vector to * oriented from K to L.
We denote for each diamond: -its edges s (e.g. s = [ K , K * ]), such segments are interfaces between diamond cells, and when necessary we will write s = D|D ′ to emphasize that s separates the diamonds D and D ′ , -ℰ D = {s, s ⊂ D and s Ω} the set of all interior sides of the diamond, s the length of s, -⃗ n sD the unit normal to s going out from D, -S = {s ∈ ℰ D , ∀D ∈ D} the set of interior edges of all diamond cells D ∈ D, -S K = {s ∈ S, such that s ⊂ K} and S K * = {s ∈ S, such that s ⊂ K * }.
Let size(T) be the maximum of the diameters of the diamond cells in D. The flattening of the triangles is measured by the angle T ∈]0, 2 ] such that sin( T ) := min D∈D | sin( D )|. We introduce a positive number reg(T) that measures the regularity of the mesh: where and * are the maximal number of edges of each primal cell and the maximal number of edges incident to any vertex.
Hp 2.3. The number reg(T) is uniformly bounded from above and below as size(T) → 0.
Accordingly, there exist two constants 1 and 2 , which both depend on reg(T), such that ∀K ∈ M, ∀K * ∈ M * ∪ M * and ∀D ∈ D such that D ∩ K ̸ = 0 and D ∩ K * ̸ = 0 we have:

Unknowns and meshes
The DDFV method for Navier-Stokes problem uses staggered unknowns. We associate to each primal volume K ∈ M ∪ M an unknown u K ∈ R 2 for the velocity, to every dual volume K * ∈ M * ∪ M * an unknown u K * ∈ R 2 for the velocity and to each diamond D ∈ D an unknown D ∈ R for the pressure. Those unknowns are collected in the families: We define now two discrete average projections, for all functions v in ( 1 (Ω)) 2 : -one on the interior: -one on the boundary: We can collect them in a shorthand notation We introduce also a centered projection on the mesh T: (This projection makes sense owing to the regularity assumption 2 (Ω) ⊂ (Ω).) Next, we define subspaces of (︀ R 2 )︀ T , which take in account Dirichlet boundary conditions. Let Γ Dir ⊂ Ω, the boundary on which homogeneous Dirichlet conditions will be imposed. When Γ Dir ̸ = Ω, we need to distinguish the subsets of the boundary mesh When Γ Dir = Ω, we simply denote E 0 the discrete space satisfying the Dirichlet condition.

Discrete operators
In this section we define the discrete operators of the DDFV scheme.
Definition 2.4. We define the discrete gradient of a vector field of (︀ R 2 )︀ T as the operator with ℳ 2 (R) the space of 2 × 2 matrices with real entries, such that for D ∈ D: where ⊗ represents the tensor product. It can equivalently be written as The discrete strain rate tensor D D : We define the discrete divergence of a vector field of (︀ R 2 )︀ T as the operator Definition 2.5. We define the discrete divergence of a tensor field of (ℳ 2 (R)) D as the operator For the boundary dual cells K * ∈ M * , we have in the definition two sums corresponding to two differents parts of the boundary of K * . Indeed, we can express K * (see Fig. 2) as an union over * = [ K , L ] and an union over a part [ K * , L ] of = [ K * , L * ]: Since we only have [ K * , L ], we have the factor 1/2 when considering the edges .
On the diamond mesh we set D : which is the operator of restriction to the boundary diamonds.
The discrete gradient and divergence operators are linked by a discrete Stokes formula. This is precisely the duality property that gives its name to the method [12], see for instance Theorem IV.9 of [32]. Theorem 2.6 (Discrete Green's formula). For all D ∈ (ℳ 2 (R)) D , u T ∈ (︀ R 2 )︀ T , we have: → n is the unitary outward normal.

Brezzi-Pitkäranta stabilization
The Inf-Sup condition is a crucial structure property for the stability of a scheme for the simulation of incompressible viscous flows. A stabilization term involving the pressure can be added to enforce this condition. This idea dates back to [6] for finite element methods. It has been adapted to the finite volume framework too [15,16] and we refer the reader to [33] for the specific case of DDFV schemes. Note that the Inf-Sup condition actually holds for a large class of meshes, which do not require any stabilization [5].
The stabilization term involves the second order discrete operator, denoted by Δ D : It resembles an approximation of the Laplace's operator (endowed with the homogeneous Neumann boundary condition), however it is consistent only under orthogonality condition (as in the case of admissible meshes, see [13,14]); that is not true in general for diamond meshes obtained from M. In relation with this operator we define a semi-norm | · | on R D that depends on the mesh: Observe that This operator is just needed to introduce a stabilisation, which is necessary in a few very specific situations, like when working with Cartesian meshes where the stabilisation prevents oscillations due to spurious modes associated to the chessboard pressure lying in the kernel of the Stokes operator. Other constructions can be considered as well, with similar purposes.

DDFV scheme for the Navier-Stokes problem on Ω
This section is concerned by the analysis of DDFV schemes for the Navier-Stokes problem with Dirichlet boundary conditions on the entire domain Ω. The choice of the Dirichlet boundary conditions is adopted on Ω for the sake of simplicity; the discussion can be readlily adapted to handle more general boundary conditions on Ω, see [23,36]. As far as the convection is treated by upwind discretization, the analysis has been performed in [32]. As mentioned above, it is convenient to extend this analysis to general -schemes where the convection term is approximated by a centered discretization plus a diffusive perturbation, which depends on a certain function , see [8,26]. In what follows, D , * will be denoted by D, to simplify the notations.

Let
∈ N * and 0 < < ∞. We note = and = for ∈ {0, . . . , }. We use an implicit Euler time-discretization, except for the nonlinear convection term, which is linearized by using a semi-implicit approximation. Here and below, the time step is supposed to be constant; of course the discussion can be directly adapted to handle variable time steps. At each time step, we shall enforce the equality which takes into account the Brezzi-Pitkäranta stabilization, with a parameter > 0.
We look for u T,[0, ] = (u ) ∈{0,..., } ∈ (︀ E 0 )︀ +1 , where, as stated in Section 2.3, E 0 is the space of discrete velocities satisfying the homogeneous Dirichlet condition, and D,[0, ] = ( ) ∈{0,..., } ∈ (R D ) +1 . The sequence is initialized with: The vector 0 is well defined since it is solution of a square system, whose matrix is invertible (owing to the fact that u 0 is supposed to belong to E 0 , and because in the formulation of the problem we take into account the constraint of null average on the pressure). With those choices of (u 0 , 0 ) we guarantee the property (3.1) at the initial time step and it will be propagated at each time step. The discrete force term is also defined by a projection over T, with P M f and P M * f . From now on, to simplify the notations we will denote (u +1 , +1 ) with (u T , D ) and (u , ) with (ū T ,¯D).
Given (ū T ,¯D) satisfying (3.1) the update u T ∈ E 0 and D ∈ R D is such that: ( ) The fluxes are defined as a sum of a "diffusion" term and a "convection" term: The diffusion fluxes are defined as: The convection fluxes are expressed as the sum of a centered discretization and a diffusive perturbation The diffusive part depends on the function , which describes the different schemes that we can work with. The centered scheme corresponds to ( ) = 0 and the upwind scheme corresponds to ( ) = 1 2 | |. However, for further purposes and the analysis of the domain decomposition method, it is relevant to consider a quite general framework where can be matrix-valued. In what follows, we denote The total fluxes then become: The definition of K , * K * comes from [32,34], up to the boundary terms. They are approximations of the fluxes: ∫︀ (u · ⃗ n K ) K (u T ) and ∫︀ * (u · ⃗ n * K * ) * * K * (u T ). Note that this part of the scheme is explicit: the velocity-pressure pair is updated by solving a linear system corresponding to a semi-explicit time discretisation.
For D ∈ D int , we can rewrite the discrete divergence div D as follows and deduce that We introduce the flux s s,D to approximate ∫︀ sū T · ⃗ n sD d for s = [ K , K * ] = D|D ′ ∈ ℰ D , noting that this flux must be perturbed by the stabilisation term imposed in (3.1), so that we have s s,D = sū Thanks to (3.1), it implies For D ∈ D ext , (see Fig. 3), a similar reasoning leads to Based on these considerations, we define the convections fluxes as follows: -For the primal edges that is for a primal cell K ∈ M and D ∈ D K : This sum contains two terms s = [ K , K * ] and s = [ K , L * ]. -For the dual edges, we have two differents cases (see Fig. 3 This sum contains two terms s = [ K , K * ] and s = [ L , K * ].
• for K * ∈ M * and D ∈ D K * ∩ D ext where s = [ K , K * ] and Ω∩ K * indicates the measure of the intersection between K * ∩ Ω and That the scheme preserves the conservation laws of the continuous problem is a remarkable property of the construction. It means that mass exchanges between the cells are well reproduced by the scheme. Proposition 3.1. Let (3.1) be satisfied. Then the fluxes K and * K * are conservative, that is to say Proof. For the interior mesh, we proceed as in [32]. If K ∈ M, by reorganizing the sum on the sides s ∈ G K belonging to the primal cell K, we obtain: since ⃗ n sD = −⃗ n sD ′ , where D and D ′ denote the two neighboring diamonds which share the edge s, of vertices K , K * . In the same way, We deduce that ∑︀ where the first sum vanishes thanks to (3.3), (3.4), and for the second term we use the fact that each vertex K * is shared by two boundary diamonds.

Well-posedness of problem ( )
The well-posedness of the scheme ( ), which is known when the centered or upwind discretization is used, generalizes to a wide class of functions . In what follows, for a × matrix , we write ≥ 0 when the symmetric part of is positive semi-definite, which means that · ≥ 0 holds for any vector ∈ R .

the (possibly matrix valued) coefficients arising in the definition of the fluxes (3.2), as the diffusive correction with respect to the centered approximation. Assume that
Then the problem ( ) is well-posed.
) is a linear system = with a rectangular matrix ∈ ℳ +1, (R), ∈ R and ∈ R +1 . Let be the following set: We have dim( ) = , (f M , f M * , 0, 0, 0) belongs to X and Im( ) ⊂ as a consequence of the Green formula in Theorem 2.6. If we show that the matrix is injective, we conclude that dim(Im( )) = and that Im( ) = . We are going to show that if f M = f M * = 0, then u T = 0 and D = 0.
We multiply the equations on the primal and dual mesh of ( ) by u T and we sum over all the control volumes: By definition of the scalar products we have 1 and, by replacing the definition of the fluxes, we get We can consider separately the terms. By replacing the definition of the divergence operator and then by appying Green's formula (Thm. 2.6) for u T ∈ E 0 , we obtain where for the last equality we use that div D (u T ) − d 2 D Δ D D = 0 and we apply (2.1). For the convection terms, we sum over diamonds recalling that u T ∈ E 0 , so we do not have boundary terms. For the centered part, we apply Propositions 3.2 and 3.1, to conclude that For the diffusive perturbation, (ℋ ) implies Putting all together, equation (3.5) becomes: from which we deduce that u T = 0 and D is a constant (we recall that > 0). Since D verifies ∑︀ D∈D D D = 0, we have D = 0.

DDFV domain decomposition
We start by defining a discretization for the problem set on the subdomain Ω . As in Section 3, the nonlinear convection term will be approximated through -schemes; we will see that the coefficients K , * K * play an important role in the convergence of the Schwarz algorithm. We start by defining the meshes, and we analyse the scheme on each subdomain, denoted by ( ), and we introduce the Schwarz algorithm for the domain decomposition. We present the study for two subdomains for the sake of simplicity, but it could be extended to an arbitrary number of adjacent subdomains. (Difficulties arise when more than two domains have common points on their interfaces, a situation which deserves a specific analysis see e.g. [18].)

DDFV on composite meshes
For each subdomain Ω of Ω, = 1, 2, we consider a DDFV mesh T = (M ∪ M , M * ∪ M * ) and the associated diamond mesh D . Note that the DDFV approach allows us to work with non conformal meshes,  (1) The two meshes share the same vertices on Γ. This, in particular, implies that the two meshes have the same degenerate volumes on Γ, i.e. M 1,Γ = M 2,Γ . (2) The center L of the degenerate volumes of the interface L = [ K * , L * ] ∈ M 1,Γ = M 2,Γ is the intersection between ( K * , L * ) and ( K1 , K2 ), where K 1 ∈ M 1 and K 2 ∈ M 2 are the two primal cells such that L ∈ K 1 and L ∈ K 2 (see Fig. 4).
In order to build a composite mesh, it is equivalent to either build a global mesh by gluing two subdomain meshes, or to build the subdomain meshes starting from a global mesh (which has edges along Γ).
Consider the composite mesh of Figure 4; remark that: K2 that intersects Γ in the domain Ω can be written as the union of diamonds D 1 , of vertices K1 , K * , L * , L , and D 2 , of vertices K2 , K * , L * , L , respectively in Ω 1 , Ω 2 . Moreover, on the subdomain meshes we have additional unknowns on L on Γ with respect to the mesh on Ω; -equivalently, a volume K * that intersects Γ in Ω is the union of K * 1 , K * 2 in Ω 1 , Ω 2 . In particular, an edge ; -an edge = [ K * , L * ] on the interface Γ is shared by all the meshes.
Due to the fact that each dual cell on the global mesh that intersects Γ is split in two between the subdomains, it is necessary to introduce some additional unknowns fluxes Ψ K * j , for all K * j ∈ M * ,Γ , as in [20]. Those unknowns are intended to approximate the dual fluxes ℱ * K * on the interface. For a diamond D ∈ D Γ , the unknowns are illustrated in Figure 5.

The subdomain problem: transmission conditions
On each subdomain Ω of Ω, we want to solve a Navier-Stokes system with mixed boundary conditions. On the fraction of the boundary that intersects Ω, we impose Dirichlet boundary conditions. On the the interface Γ between the two subdomains, we impose the discretized version of the transmission conditions (1.2).
To Precisely, the solenoidal constraint is approximated on the diamond mesh D and for the diamonds in D Γ a transmission term is added, controlled by the parameter . We give now, formally, an hint of why it is necessary to add this condition on the interface diamonds D Γ : our goal is to recover, at convergence of the Schwarz algorithm, the solenoidal constraint div D (u T ) = 0 for all D that intersect Γ in Ω (that we write here for sake of simplicity without the stabilization term). As described in Definition 4.1, a diamond D in Ω that intersects Γ can be written as the union of diamonds D 1 , D 2 in Ω 1 , Ω 2 . By definition of the discrete divergence, see Section 2.4, we wish to decompose D div D u T = D1 div D1 u T 1 + D2 div D2 u T 2 . Therefore, we expect on Γ an expression which would look like However, equation (4.1) does not make sense and a detailed construction of the . This can be understood by a dimensional argument. Of course, we naturally identify the values of u T and u T when they are evaluated on common points of the grids T and T . But, we should bear Moreover, imposing a condition of this kind along the iterations of the Schwarz algorithm is not sufficient to prove convergence of the algorithm, as we will show later in Theorem 5.8; in order to apply the analytical tools of the proof, it is necessary to add to (4.1) a Fourier-like term for the pressure, controlled by .
The DDFV discretization leads to the following system on Ω : Here, we denote f T = P T f . We will refer to the system ( ) in the shorthand form:

Remark 4.2.
When we impose transmission conditions in Schwarz' algorithm, we are led to approximate the boundary term , which keeps track of the anti-symmetrization of the convection term. Formally, at the continuous level, if is a test function in = { ∈ ( 1 (Ω)) 2 , Ψ| Γ Dir = 0, div( ) = 0}, the variational formulation of (1.1) reads: The convection term can be written as by integration by parts, since u is divergence free. Coming back to (4.2), we integrate by parts also the diffusion terms, and we end up with: This is the reason why, when working with transmission conditions, we impose a condition on (u, )⃗ n− 1 2 (u·⃗ n)u, that contains just "half of the convection". Besides, the numerical flux ℱ K is constructed to approximate the term This is why in the approximation it gives: The proof relies on the following energy estimate (where we bear in mind that K and * K * can be matrixvalued).
Theorem 4.4 (Energy estimate on ( )). Under the hypothesis (ℋ ), the scheme ( ) satisfies the following relation -the conditions on the equation of mass conservation: This implies: If now we impose that all data f T , h T , D and u T in (4.4) vanish, we have: from which we deduce that u T = 0 and D is a constant (since > 0). Thanks to the transmission conditions on D Γ , since > 0 and u T = 0, we obtain D = 0. Finally, thanks to the transmission condition on M * ,Γ and u T = 0, we also have Ψ T = 0.
Proof of Theorem 4.4. We multiply the equations on the primal and dual mesh of ( ) by u T and we sum over all the control volumes: By definition of the scalar products we have and, by rewriting the fluxes as a sum of the diffusive and convective contribution we obtain: We consider separately the two contributions. For the diffusion terms, we have, by the definition of the divergence operator We can now apply Green's formula to the RHS, and remark that We thus find: By the definition of the trace operator, we obtain: For the convection terms, we get We estimate the term 1 ; we first integrate by parts thanks to Proposition 3.1 and (ℋ ): We replace the definition of ℱ K for all D ∈ D : Passing to the sum over primal cells K for the first term and applying Proposition 3.2 we get: It can be rewritten as: We estimate the term 2 ; we first integrate by parts, by using Proposition 3.1 and (ℋ ): We replace the definition of ℱ K for all D ∈ D : Passing to the sum over dual cells K * for the first term we get: From the definition of * K * and by Proposition 3.2 we have that ∑︀ ,Γ , that gives: Gathering (4.5)-(4.7) together, we find: Since ℱ K + ℱ K = ℱ K , it leads to (4.3).

Convergence analysis of the DDFV Schwarz algorithm
Bearing in mind the properties of the mesh discussed after Definition 4.1, we infer that the asymptotic fluxes as → ∞ should satisfy In order to obtain these relations, it will become necessary to modify the fluxes on the interface, either for the limit or for the subdomain problem. For this reason, the convergence will be studied in two steps. In this Section, we shall identify the limit of the Schwarz algorithm defined in Section 4.3. We focus here on the natural situation where K , * K * take scalar values, like with the upwind and centered discretizations. We will show that this limit is still a DDFV scheme for the problem (1.1), but with modified fluxes on Γ. We will then prove convergence to this limit scheme, to which we will refer to as (̃︀). In the next Section 6, we will show that it is possible to obtain ( ) asymptotically, at the price of modifying the fluxes of the Schwarz algorithm ( 1 ), dealing with matrix coefficients K , * K * . That ( ) provides a consistent approximation of the Navier-Stokes equation is justified, with error estimates, in [33,34]; the arguments can be adapted for (̃︀) and numerical tests of convergence are presented in [25]. Here, we focus on the convergence of the Schwarz iterations.

The limit problem (̃︀)
We consider the following DDFV scheme for (1.1), on the domain Ω: given (ū T ,¯D), satisfying (3.1), we look for u T ∈ E 0 and D ∈ R D such that: In the interior of the domain, the fluxes coincide with the fluxes in ( ), see (3.2). On the interface, they are defined as:̃︀ wherẽ︀ K and̃︀ * K * are matrix-valued quantities that come from the transmission condition of the iterative process. Their expressions are established in Propositions 5.4 and 5.5.

Definition of̃︁ K and̃︁ * K *
Let us start with some preliminary definition, bearing in mind that K , * K * are supposed to be scalars. Definition 5.1. For = 1, 2, and ∈ M ,Γ , we set = Id + ⃗ n K ⊗ ⃗ n K and Remark 5.2. The matrix = 1 + 2 is symmetric and definite positive, thus invertible, since it is the sum of two symmetric and definite positive matrices. In fact, with ⃗ n K = (︂ )︂ , we have: which is symmetric and for any = owing to Hypothesis (ℋ ). For , = 1, 2, ̸ = , since and are polynomial in , the following properties hold: The fluxes̃︀ ℱ K ,̃︀ ℱ * K * are constructed in order to satisfy the properties (5.1) and (5.2). The system (̃︀) is a scheme defined on the mesh T on Ω; in particular, this means that there are no additional unknowns u L on the interface Γ, see Figure 4. The following results apply for a general diamond: Proposition 5.3. Let D ∈ D Γ be a diamond and let D 1 , D 2 be the two semi-diamonds such that D = D 1 ∪ D 2 , see Figure 6. We denote by ( K , K * , L * , L ) the vertices of D and by ( K1 , K * , L * , ), ( K2 , K * , L * , ) the vertices of D 1 and D 2 . Let = K 1 |K 2 , and let , 1 , 2 be as in Definition 5.1. Then, there exists a unique u , given by  The strain rate tensors can be written by using the matrix as: , the contributions of the pressure D and of the velocity u K * , u L * on the vertices cancel out. So (5.5) becomes: We group the terms in u thanks to K1 = − K2 , and we obtain: By Definition 5.1, equation (5.8) becomes: It is sufficient to show that this expression is injective; if (u T , D ) is equal to zero, we are going to show that u is zero. This is true because, if (u T , D ) vanishes, this means in particular u K1 = u K2 = 0 and (5.9) becomes u = 0. Since is definite positive, see Remark 5.2, we deduce u = 0.
It is possible to obtain property (5.1), by adapting the fluxes on the interface.

Proposition 5.4.
Let D be a diamond and let D 1 , D 2 be the two semi-diamonds such that D = D 1 ∪D 2 , see Figure 6. Then there exists a unique flux̃︀ ℱ K on = K 1 |K 2 such that︀ given bỹ︀ Proof. We consider ℱ K1 and we refer the reader to Figure 6: we recall that it is a flux on the semi-diamond D 1 of vertices K1 , K * , L * , . Thanks to (5.6), it can be written as: . By grouping the terms in u K1 and u in ℱ K1 , we are thus led to that can be written as: According to Remark 5.2, the matrices and commute, for = 1, 2. Hence, we can write We develop the computations and we find: Let̃︀ K be the matrix defined in (5.11). We get: Fig. 6), we end up with: We remark that now the expression of ℱ K1 depends only on the unknowns u K1 , u K2 , u K * , u L * ; so it is a flux defined on the entire diamond D (see Fig. 6). It can be rewritten as: So that we find (5.10).
We proceed similarly to obtain (5.2): Proposition 5.5. Let D be a diamond and let D 1 , D 2 be the two semi-diamonds such that D = D 1 ∪D 2 , see Figure 6.
given bỹ︀ Proof. This is a direct consequence of the computation of (5.12). By definition, So if we take the sum of the two fluxes, we get: which ends the proof.

Well-posedness of the limit problem (̃︀)
The expression of the new fluxes̃︀ ℱ K ,̃︀ ℱ * K * permits us to justify the well-posedness of (̃︀).
Theorem 5.6. Under Hypothesis (ℋ ) for K , * K * , problem (̃︀) is well-posed. Proof. By Theorem 3.3, we need to verify that Hypothesis (ℋ ) holds. Since we are supposing it for K , * K * , we just need to check it for̃︀ K ,̃︀ * K * , the modified fluxes on the interface. As a direct consequence of (5.11) and (5.14), we have:̃︀ In fact, if we consider a diamond on the interface between the two subdomains Ω 1 , Ω 2 , it can be seen as the one in Figure 6. For K = K 1 and L = K 2 , equation (5.11) reads︀ Observe that , do not depend on the index of the subdomain; moreover, we have K1 = − K2 , so that ( K1 ) 2 = ( K2 ) 2 and 1 2 = 2 1 from Remark 5.2. We conclude that̃︀ K =̃︀ L . For the dual flux, equation (5.14) becomes︀ Thanks to (ℋ ), we have )︂ ; then: thanks to Hypothesis (ℋ ) on K , that ensures ≥ 0 and den > 0. Sõ︀ K is semi-definite positive. For what concerns the dual flux, by (5.14) we obtain directly that̃︀ * K * is semi-definite positive since it is the sum of the two semi-definite positive matrices Further comments on problem (̃︀) can be found in [25].

Identification of the limit
In order to prove the convergence of the Schwarz algorithm towards the solution of (̃︀), it is necessary to project this solution, that is defined on Ω, on the subdomains Ω , = 1, 2.
, such that: Proof -for all L = L ∈ M ,Γ and for all K * ∈ M * such that K * ∈ Γ,

Consequence on the equations
We now show that from a solution (u T , D ) of the DDFV scheme (̃︀) we built a solution to (̃︀ ∞ ): -∀K ∈ M, (u T , D ) satisfies: If we look at the composite mesh (see Fig. 4), we remark that the primal cells K ∈ M correspond to Kj ∈ M (or to Ki ∈ M ). This implies that Moreover, for a diamond D ∈ D K ∖D Γ K , remark that the limit unknowns we have: ∑︁ thanks to the choice (5.15) of u ∞ L for all L ∈ M ,Γ and thanks to Proposition 5.4, we havẽ︀ -∀K * ∈ M * , (u T , D ) satisfies: We need to distinguish two cases.
(1) If K * ∩ Γ = ∅, equation (5.17) reduces to: and the cells K * ∈ M * correspond to K * j ∈ M * (or to Ki ∈ M * ). This implies that j and * = * . Moreover, for a diamond D ∈ D K * ∖ D Γ K * , (that is the case here since we are supposing K * ∩ Γ = ∅) remark that the limit unknowns u ∞ we have: ∑︁ satisfies on the interior dual mesh: (2) If K * ∩ Γ ̸ = ∅, the cell K * can be written as the union of Kj ∈ M * ,Γ and Ki ∈ M * ,Γ . This implies we have: ∑︁ For a diamond D ∈ D Γ K * , thanks to (5.12), we have We deduce from (5.17): satisfies on the boundary dual mesh: -∀K ∈ M, with K ∩ Γ ̸ = ∅, if we look at the composite mesh, the diamond D ∈ D Γ K can be written as the union of D ∈ D Γ K and D ∈ D Γ K . By definition, we have K = − K ; moreover, thanks to the choice (5.15) of u ∞ L for all L ∈ M ,Γ and thanks to Proposition 5.4, we have T , we get the relation: -∀K * ∈ M * , with K * ∩ Γ ̸ = ∅, the cell K * can be written as the union of Kj ∈ M * ,Γ and Ki ∈ M * ,Γ . By definition, we have K * . This leads, from the definition of h ∞ T , to the relation: -for all D ∈ D, (u T , D ) satisfies: We need to distinguish two cases: (1) If D ∩Γ = ∅, the diamond D coincides with a diamond D ∈ D (or with a diamond D ∈ D ). For a diamond D ∈ D ∖ D Γ , remark that the limit unknowns u ∞ K , u ∞ K * j , ∞ D on T for = 1, 2 coincide with u K , u K * , D on T. Thus we can directly deduce that (2) If D ∩ Γ ̸ = ∅, the diamond D can be written as the union of D ∈ D Γ and D ∈ D Γ . This implies that the divergence can be split as: . From (5.22), the choice of unknowns ∞ D and from the definition of ∞ D we obtain:

Convergence of the DDFV Schwarz algorithm towards (̃︀)
Theorem 5.8 (Convergence of the discrete Schwarz algorithm). Under the hypothesis that * = 2 * = 2 * for , = 1, 2, ̸ = , the iterates of the Schwarz algorithm ( 1 ) and ( 2 ) converge as tends to infinity to the solution of the DDFV scheme (̃︀) (up to a constant for the pressure).
The relation imposed on the * 's -which also appears in Theorem 6.1 below -is not a restriction in practice. Given a mesh of the entire domain and the interface Γ, we can adjust the centers K neighboring the interface so that the condition is fulfilled. It ensures that the diamonds of the global mesh are split into two half-diamonds with equal area, see Figure 4. All meshes used for the simulations satisfy this condition.
Proof. The iterates of ( 1 ) and ( 2 ) satisfy: , constructed from the solution of (̃︀) is solution to: We define the errors By linearity, they satisfy: To prove the convergence of the iterates of Schwarz algorithm, it is sufficient to prove the convergence to 0 of the solution of (5.26). In the expanded form, equation (5.26) is written as: Thanks to the hypothesis * = 2 * = 2 * , we have D = D ; so, in the equation on D ∈ D Γ , we can simplify the measures and it becomes: We multiply the equations by e T and we sum over all the control volumes, as in the proof of Theorem 4.4. We obtain, analogously to (4.3), the following: By the equations on D Γ , we can split the scalar product into interior diamonds D ∖ D Γ and boundary diamonds D Γ : for the diamonds D ∈ D ∖ D Γ we apply the equation of conservation of mass, for the diamonds D ∈ D Γ we add and subtract the term ∑︀ We apply (2.1) to the term − ∑︀ ; we then multiply and divide ∑︀ )︁ by to finally obtain: So (5.27) becomes: where we multiplied and divided by > 0 the terms on the second line.
We start by considering 1 ︁ · e L j . By applying now the equality − = 1 4 (︀ (− + ) 2 − ( + ) 2 )︀ we can write: Owing to the transmission conditions it becomes: Equivalently for 1 , we obtain: The hypothesis * = 2 * = 2 * implies D = D , so that this expression becomes: Replacing those results into (5.28), we have: Summing over = 0, . . . , max and = 1, 2 we obtain: that shows how the total energy stays bounded as the iteration index max goes to infinity. The series ∑︀ =1,2 |Π D | 2 converge, so their general term tends to zero, that implies the convergence to zero of the errors ‖e T ‖ 2 2 , |Π D | 2 , defined in (5.25). Thus the algorithm converges. The limit is the solution of problem (̃︀), that is problem ( ) with an appropriate choice of the flux on Γ; in fact, we can deduce that, as max goes to infinity: -‖e T ‖ 2 2 tends to zero implies u T → u ∞ T for = 1, 2. -|Π D | 2 tends to zero implies (since | · | is a semi-norm): D + const (Ω ) → ∞ D for = 1, 2. Thus the pressure converges up to a constant that depends on the subdomain. In some cases we are able to determine const (Ω ).
Remark 5.9. We can determine the constant const (Ω ) if we suppose that the mesh satisfies the Inf-Sup inequality [5]. In fact, this implies that the norm

A modified DDFV Schwarz algorithm
We now investigate whether it is possible to construct a discrete Schwarz algorithm with modified fluxes that converges to the solution of ( ). We show that this is possible if we suppose an asymmetric discretization of (1.1), in the sense that we need to consider an upwind discretization of the convection term on the primal mesh and a centered scheme on the dual mesh, that corresponds to the choice in ( ). We remind the reader that the convergence to (̃︀) holds if and only if both (5.11) and (5.14) hold, which can be seen as a definition of̃︀ K (resp.̃︀ * K * ) as a function of K1 , K2 (resp. * 1 K * , * 2 K * ). The idea is to modify the Schwarz algorithm, so that it converges to the solution of ( ); this relies on the ability to invert these relations. Accordingly, the fluxes of the limit equation depend only on K , * K * and a different definition of the fluxes is not required on the interface Γ.
with:¯K Proof. The assumption * = 2 * = 2 * implies that D1 = D2 = 1 2 D and K1 = K2 = K . This means that Moreover, −1 = Re D 2 2 ( + K Id) −1 . Therefore, we get︀ Expanding this expression, we get:︀ by using the definition of and −1 . Let us set = Re D K , We have 2 , so we end up with:̃︀ If we make explicit the dependences of̃︀ K , K as a function of , since K is a function of D Re K and̃︀ K a function of 2 DRe K , we are led tõ︀ We can rewrite this condition as:̃︀ . This relation implies that the Schwarz algorithm ( 1 ) and ( 2 ), whose convection fluxes depend on K , converges towards the solution of (̃︀), whose convection fluxes depend oñ︀ K for ∈ ℰ Γ .
We want to build a new Schwarz algorithm (¯) that converges toward ( ), whose fluxes are defined by K ; so we need to build¯K such that: K can be a full matrix. Since our goal is to converge towards the fluxes that define an upwind scheme, i.e. defined by ( ) = 1 2 | |, K is actually a diagonal matrix, that will be denoted by K Id to make its matrix nature clear.
Thus we need to invert the function ϒ defined above to find the new coefficients¯K. The inverse of ϒ does not exist for every K . Given and K (2 ), we have a second-degree equation for¯K( ): that is:¯K ( ) 2 +¯K( ) + = 0.
Since the matrices , are symmetric and they commute (because they are polynomials in ), they can be diagonalized using the same basis of eigenvectors. The matrix =

(︂ −
)︂ is an orthogonal matrix, and we can write: with̃︀ and̃︀ diagonal matrices, whose expressions are:︀ We then look for¯K( ) of the form¯K( ) =̃︁ −1 , with̃︁ a diagonal matrix such that︁ Since we are supposing K ( ) = 1 2 | |, the solution exists and is given bỹ︁ that leads to our result (6.1).
For what concerns property (5.14), we would like to define a unique * K * ( * ) for * ∈ ℰ * in the limit scheme ( ). With the assumption * = 2 * 1 = 2 * 2 , we can define * = Re D * * K * and * = Re D * * K * for = 1, 2: remark that there is no relation between the * . We have * = * + * , since * * K * + * * K * = * * K * . This leads to the new expression for (5.14):︀ This is true only if * K * =̃︀ * K * = 0; in this way, property (5.13) is verified. So the dual flux for the algorithm (¯) and for the limit ( ) correspond to a centered discretization of the convection flux on the dual mesh.

Numerical results
In this section, the objectives are the following: -to show and compare the convergence properties of the Schwarz algorithms ( 1 ) and ( 2 ) (presented in Sect. 4.3)) and (¯) (presented in Sect. 6); -to study on numerical grounds the influence of the parameters , , of (1.2) on the convergence; -to further validate the method with the simulation of a benchmark of a flow past an obstacle.
We will refer to ( 1 ) and ( 2 ) as "first Schwarz algorithm", and to (¯) as "second Schwarz algorithm". We recall that the difference between the two algorithms relies in the definition of the fluxes at the interface; the former converges towards the solution of (̃︀) (see Thm. 5.8), the latter towards the solution of ( ) (see Thm. 6.1). For the first Schwarz algorithm, in all the following test cases, we will consider an upwind discretization of the convection flux, i.e. we set ( ) = 1 2 | |. The scheme needs the resolution of large linear systems; for the simulations discussed below, the linear systems -possibly non symmetric due to the interface -are treated by a direct method, appealing to Umfpack libraries.

Numerical study of the convergence
We recall that the domain decomposition algorithm is an iterative algorithm that is employed at each time step; this, in particular, implies that at each iteration of the Schwarz algorithm we solve a steady problem. In the following tests, we fix the time step ( = 10 −4 ), and we apply the iterative method on the time interval [0, ]. The time step is voluntarily picked quite small here in order to focus the discussion on the spacial error and the effects of the interface; its influence is also discussed below (see Fig. 14). In all the test cases, the domain Ω = [−1, 1]×[0, 1] will be divided into two subdomains Ω = Ω 1 ∪Ω 2 . The meshes we will consider are illustrated, in their first level of refinement, in Figure 7. The subscript in the name of the mesh (see Fig. 7) denotes the level of refinement, i.e. Mesh 1 represents the coarse mesh of a family of refined meshes (Mesh ) . More precisely, Mesh is obtained by dividing by two all the edges of Mesh −1 . We consider the following exact solution to (1.1): The algorithms, in all the following simulations, are initialized with initial random guesses h 0 T and 0 D for = 1, 2. As a stopping criterion, we impose: where the errors are defined in (5.25).

Error on the interface
In this first test case, we consider the first Schwarz algorithm; our goal is to point out that the error computed with respect to the solution of (̃︀), along the iterations of the algorithm, stays localized at the interface between the two subdomains.
The domain Ω is meshed with Mesh 3 5 , we fix the parameters = 100, = 1, = 10 −2 . Since the initialization assigns random values, the initial error is ‖u 0 T − u T ‖ ∞ = 100. for both primal and dual mesh. As we pass to the 1st iteration, we observe in Figure 8 how it immediately locates on the interface between the subdomains; it decreases, passing from 100 to 1.9 on the primal mesh and to 6.9 on the dual mesh. Already at the 10th iteration we see in Figure 9 how it has diminished, staying localized on the interface, passing from 1.9 to 0.52 on the primal mesh and from to 6.9 to 0.05 on the dual mesh.

Study of the parameters
In this section our goal is to study the influence of the parameters , , and of the mesh on the convergence of the first and second Schwarz algorithms. We recall that is associated to the Brezzi-Pitkäranta stabilization   (see Sect. 2.6) while the parameters and are associated to the transmission conditions between subdomains. In each test case we fix all parameters, but one. First, the value of associated to the stabilization is set to 10 −2 . In Figures 10 and 11 we represent on the -axis the number of iterations, on the -axis the error.
We start by considering the first Schwarz algorithm; we can observe in Figure 10 the convergence of the algorithm to the solution of Test 1 on Mesh 1 1 .  In particular, on the left of Figure 10, is fixed to 1, and we observe how, as increases, the number of iterations necessary to converge decreases until = 200; past this critical value, the number of iterations starts to increase again. This suggests that on Mesh 1 1 , for = 1 and = 10 −2 , = 200 is a good choice to have a better convergence. On the right of Figure 10, we set = 100 and we let vary: we observe the same kind of behavior as the one of . We consider now the second Schwarz algorithm on the same test case, i.e. Test 1 on Mesh 1 1 . We show its convergence in Figure 11. This indicates that on Mesh 1 1 , for = 100 and = 10 −2 , = 0.25 is a good choice to have a better convergence. We remark that the second Schwarz algorithm behaves similarly to the first one, if we compare Figures 10  and 11; thus, both algorithms converge and the speed of convergence is influenced by the choice of or , once fixed the value of and the mesh. Since the parameters have the same behavior and the number of iterations necessary to the convergence is almost identical between the two algorithms, from now on we will only focus on the first one.  In the following, in Figures 12-16 we represent on the -axis the value of the parameter under study, and on the -axis the number of iterations.
In the first test case of Figure 12, our goal is to show how the level of refinement of the mesh can influence the choice of the optimal parameter; we consider Test 1 on the family (Mesh 1 ) , = 1, 2, 3, 4. As before, we set the value of , we fix one of the two parameters between and and we let the other vary; we represent on the -axis the value of the parameter that changes, on the -axis the number of iterations required to obtain an error of order 10 −5 .
As illustrated in Figure 12 and summarized in Table 1, we observe different results for the two parameters; the mesh refinement has an impact on but not really on . The mesh size ℎ is divided by two at each level of refinement, and we see that it has an influence on the value of ; unfortunately, we can not conclude by defining a relation between the two.
In Figure 13 (left) and Table 2 we want to confirm the results obtained for on Figure 12 (left) and Table 1, by considering the same test case (Test 1) on a different family of meshes, (Mesh 3 ) , = 1, 2, 3, 4. As before,   is influenced by the mesh discretization step but we can not conclude by defining a relation between the two; moreover, we remark that its optimal values change with respect to Table 1, due to the different meshes.
In Figure 14 (left) and Table 3 we want to point out that also the choice of the time step influences the optimal . In fact, we can see that for a bigger (such as = 10 −3 ), the optimal is smaller ( = 21.18), and the more the decreases, the bigger becomes the value of .   In Figure 15 and in Table 4 we study the influence of the parameter , associated to the Brezzi-Pitkäranta stabilization. We see how the choice of this parameter affects the convergence of the algorithm and how it affects the optimal value of : we pass from 818 iterations with = 436.81 (for = 10 −4 ) to 40 iterations with = 122 (for = 10 −1 ). There is then an optimal choice even for this parameter.
As last simulation, on Figure 16 and Table 5 we compare the optimal values of for Test 1 on different meshes. We see that even the choice of the mesh influences the optimal choice of the parameter: for a Cartesian mesh, = 105.91 while for Mesh 2 1 = 154.3. For every test case, we have observed that the following four parameters -the parameters and associated to the transmission conditions, the parameter of the Brezzi-Pitkäranta stabilization and the mesh choice (its geometry and its level of refinement) -impact the convergence of the algorithm. Considering three of these parameters as fixed, it is possible to optimize the remaining one in order to reach a faster convergence. This preliminary study also confirms the high interdependence between the parameters: a conclusion can be substantially changed by modifying the conditions, like the geometry of the mesh. The understanding of the elliptic equation, which is still in its infancy, can give relevant hints [20,21], as well as a further analysis of the continuous system.

Cylinder test case
In this section, we test the first Schwarz algorithm ( 1 ) and ( 2 ), on a test case inspired by the benchmark of [38] (we precisely use the detailed results in [31]). In both [31,38], the drag and lift coefficients of the flow past an obstacle are computed from simulations on a domain Ω, with Dirichlet boundary conditions. Our goal is to measure the quality of the DDFV solution obtained on Ω with a domain decomposition algorithm.  The benchmark is defined with dimensional equations, so we adopt the same framework, see Figure 17. References [31,38] consider a long channel Ω = [0, 2.2 m] ×[0, 0.41 m] with a cylindrical obstacle whose center is at (0.2 m, 0.2 m); we decompose the domain Ω into two subdomains Ω 1 , Ω 2 and we place the interface Γ at = 0.56 m. On Ω we impose Dirichlet boundary conditions, as in [31].
The mesh that we consider on Ω is represented in Figure 18; it is obtained with GMSH, it has 34 634 cells and it is locally refined around the cylinder. Remark that on the left domain Ω 1 (the one with triangles) there are 18 250 cells and on the right domain Ω 2 (the one with rectangles) there are 16 384 cells.
The initial condition is u ( , ) = (0, 0). The density of the fluid is given by = 1 kg m −3 , and the reference velocity is¯= 1 m s −1 (note that the maximum velocity is 3 2¯) . The diameter of the cylinder is = 0.1 m, so that the Reynold's number is 0 ≤ Re( ) ≤ 100. The time step is = 0.00166 s.
We use the limit scheme (̃︀), but at some fixed times, we use instead the iterative Schwarz algorithm ( 1 ) and ( 2 ), with the initial guesses h 0 T and 0 D , for = 1, 2, given by (̃︀) at the previous time step. The stopping criterion is max (︁ ‖e T ‖ 2 , ‖Π D ‖ 2 )︁ < 10 −3 (7.1) and the values of the parameters are set to = 200, = 1 and = 0.01. To discretize the convection term, we choose = 0 which gives a second order accurate method.
To start with, we consider the profile of the first component of the velocity. The iterative algorithm is applied at each time step. We compare the solution of the limit problem (̃︀) (Figs. 19 and 21) to the solution obtained by the iterative algorithm ( 1 ) and ( 2 ) (Figs. 20 and 22) at times = 2 s and = 6 s. As we can see, the profile is the same and the domain decomposition does not introduce any spurious modes to the solution; the    Drag and lift coefficients. In order to measure the quality of the obtained results, we focus on the computation of the drag and lift coefficients for the limit problem (̃︀) and for the solution of ( 1 ) and ( 2 ). We define the drag coefficient ( ) and the lift coefficient ( ) as where stands for the boundary of the obstacle, ⃗ n = ( , ) is the normal vector on pointing to Ω, t = ( , − ) the tangential vector and u the tangential velocity. In the DDFV setting, We study the evolution of the coefficients in Figure 23 and their maximum value in Table 6, defined as: The results shown in Table 6 and in Figure 23 prove that the approximation given by the limit problem (̃︀) and the results obtained with the Schwarz algorithm ( 1 ) and ( 2 ) are robust and quantitatively correct. The behavior of the drag and lift coefficients of (̃︀) is coherent with the reference values from [31], and the extreme values of both coefficients are similar, see Table 6. The slight discrepancy in the maximum value of the coefficients is due to level of refinement of the mesh and to the order of the scheme: we work with approximately 90 000 unknowns, for all velocity components and pressure, compared to the approximately 500 000 unknowns used in [31]. Figure 23 shows that the lift coefficient is sensitive to the choice of the time discretization: the time step in [31] is = 0.00125 s with a second order scheme in time. Our scheme is first order in time, and we work with = 0.00166 s. We have implemented a second order backward difference formula in time, as in [23]: the first iteration in time remains unchanged, while for ∈ {1, . . . , } the term u is discretized by 1 (︀ 3 2 u +1 − 2u + 1 2 u −1 )︀ instead of 1 (︀ u +1 − u )︀ and the convection fluxes K depend on (︀ 2u − u −1 )︀ instead of u . This approach indeed improves the quality of the approximation of the lift coefficient, see Figure 23. The drag and lift coefficients associated to the domain decomposition method ( 1 ) and ( 2 ) have been computed with the second order scheme. The iterative process is applied at each time step; we then compute the coefficients associated to the solution given by the algorithm. The results are illustrated in Figure 23, where we can observe that the values of the coefficients associated to the Schwarz algorithm reproduce the curves given by the solution of (̃︀); since we established that they are a coherent reproduction of the reference values in [31], we can conclude that the algorithm produces a good approximation of the solution of the Navier-Stokes problem on the entire domain. Multi-domains. We study now the convergence of the Schwarz algorithm in the case of more than two subdomains; in particular, we decompose the domain Ω into four and five subdomains.
Then, we decompose Ω into five subdomains Ω =      and the values of the parameters are set, as for the two-subdomains case, to = 200, = 1, = 0.01, = 0. We compare the results at = 6 s, when the flow is more turbulent and with sensitive variations of the flow in all subdomains; we take as a reference the solution shown in Figure 21. As we can see in Figures 24-26, the profiles are similar and the introduction of more subdomains does not affect the solution. As resumed in Table 7, we see that by increasing the number of subdomains, we increase the number of iterations necessary to the convergence; this is due to the fact that the subdomains have to share the information between an increasing number of interfaces. Nevertheless, we gain computational time since, as we decompose the domain, we have to solve smaller linear systems.

Conclusion
This paper establishes the well-posedness of DDFV schemes for solving the incompressible Navier-Stokes system on the entire domain Ω with general convection fluxes defined by means of -schemes, and it proposes two non-overlapping DDFV Schwarz algorithms. DDFV discretizations are constructed with suitable transmission conditions, which are equally well-posed. When using standard convection fluxes in the domain decomposition method, the iterative process converges to a system with modified fluxes at the interface. However, it is possible to modify the fluxes of the domain decomposition algorithm so that it converges to the reference scheme on the entire domain. The algorithms are numerically tested on classical benchmarks, and the numerical experiments also shed some light on the role of the parameters of the method.