Ryadh Haferssas
email: rhaferss@mit.edu

Pierre Jolivet
email: pierre.jolivet@enseeiht.fr

Samuele Rubino
email: samuele@us.es

Efficient and scalable discretization of the Navier-Stokes equations with LPS modeling

Keywords: Navier-Stokes equations, LPS by interpolation, Pressure-correction methods, Large eddy simulation, Turbulent incompressible flows, Domain decomposition

In this work, we address the solution of the Navier-Stokes equations (NSE) by a Finite Element (FE) Local Projection Stabilization (LPS) method. The focus is on a LPS method that has one level, in the sense that it is defined on a single mesh, and in which the projection-stabilized structure of standard LPS methods is replaced by an interpolation-stabilized structure, which only acts on the high frequency components of the flow. As a main contribution, we propose and test an efficient discretization of the model via a stable velocity-pressure segregation, using semi-implicit Backward Differentiation Formulas (BDF) in time. On the one hand, numerical studies illustrate that the solver accurately reproduces first and second-order statistics of benchmark turbulent flows for relatively coarse meshes. On the other hand, they show that the solver works in an efficient (i.e., robust and fast) way, especially when interfaced with scalable domain decomposition methods. Such scalability results are obtained on up to 16,384 cores with a near-ideal speedup.

Introduction

In the present paper, we propose an efficient space-time discretization of the Navier-Stokes Equations (NSE) for the simulation of laminar and turbulent incompressible flows, with a special emphasis on their numerical solution in a parallel setting. Turbulent flows occur in many physical contexts (e.g., external flows in civil engineering, hydrodynamics, aeronautical applications; internal flows in hemodynamics).

The major problem in treating turbulent flows is due to the wide range of scales involved. For laminar flows already, a substantial range of scales may be encountered. Dealing with turbulent flows, however, implies an even broader range of scales in comparison to laminar flows, which are also in nonlinear interactions with each other. This makes their accurate and efficient simulation a really challenging task, even in a High Performance Computing (HPC) framework. Indeed, Direct Numerical Simulation (DNS) demands high computational effort (beyond the limits of the currently available computer power in most cases) to accurately solve with extremely fine grids the broad range of scales involved.

Conversely, Large Eddy Simulation (LES) [START_REF] Sagaut | Large Eddy Simulation for Incompressible Flows[END_REF] is an intermediate approach in its requirement of computational effort and degree of modeling. The strategy of LES consists in solving the largest flow structures and modeling the effect of the smallest flow structures on the largest ones. The traditional LES model relies on a filter to separate resolved and unresolved scales at the continuous level. On the one hand, a coarser discretization, which is substantially coarser than a comparable DNS discretization, is sufficient for resolving the large scales and, on the other hand, the universal character of the statistical behavior for the small scales justifies the modeling process, cf. [START_REF] Kolmogorov | The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers[END_REF]. These models are based on a physical approach defined considering the physical phenomena that take place on the smallest scales. Nevertheless, another way to model the fine scales in a LES method can be developed by a purely numerical approach that does not introduce any modification of the governing equations at the continuous level. This last numerical approach, which hence relies on purely numerical artifacts without any modification of the continuous problem, was seldom followed, the MILES (Monotone Integrated LES) approach [START_REF] Boris | New insights into large eddy simulation[END_REF] being the main exception, until the residual-based Variational MultiScale (VMS) models were introduced in the seminal papers [START_REF] Hughes | Multiscale phenomena: Green's functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods[END_REF][START_REF] Hughes | The variational multiscale method-a paradigm for computational mechanics[END_REF] and subsequently proposed as implicit LES techniques (ILES) for turbulent flows in [START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF][START_REF] Codina | Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods[END_REF]. These models provide a unified framework for the definition of spatial approximation schemes capable of preventing numerical instabilities that arise when the standard Galerkin Finite Element (FE) method is used, and adequate to represent the turbulence LES modeling. The basic concept consists in differentiating two scale groups (resolved and unresolved scales) using a residual-based model of the unresolved scales to account for their influence into the resolved ones. In contrast to the use of a filter in the aforementioned traditional LES, a variational projection between function spaces separates scale ranges within the VMS method. Thus, residual-based VMS are intrinsically discrete models, and no approximation of an intermediate averaged model is needed. The residual-based VMS procedure does not make use of the statistical theory of equilibrium turbulence, and thus no ad-hoc eddy viscosity is required. Furthermore, it strongly retains numerical consistency. However, the subgrid terms have a rather complex structure, since they involve the full residual with convective interactions between resolved and unresolved scales, thus increasing computational complexity and setting serious numerical difficulty just to prove stability.

In this work, we focus on an alternative strategy starting from a high-order term-by-term stabilization, cf. [START_REF] Chacón Rebollo | A high order term-by-term stabilization solver for incompressible flow problems[END_REF], which does not involve the full residual, and presents a simpler and less expensive structure for practical implementations. This method is a particular type of Local Projection Stabilization (LPS) scheme that may be cast in the VMS framework [START_REF] Ahmed | A review of variational multiscale methods for the simulation of turbulent incompressible flows[END_REF], and constitutes a low-cost, accurate solver (of optimal order) for incompressible flows, despite being only weakly consistent. It presents the same structure of the Streamline Derivative-based (SD-based) LPS model [START_REF] Braack | Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method[END_REF][START_REF] Knobloch | Local projection stabilization for advection-diffusion-reaction problems: one-level vs. two-level approach[END_REF], but it differs from it because at the same time it uses continuous buffer functions, it does not need enriched FE spaces, it does not need element-wise projections satisfying suitable orthogonality properties, and it does not need different nested meshes. An interpolant-stabilized structure of Scott-Zhang type replaces the projectionstabilized structure of standard LPS methods. The interpolation operator takes its values in a continuous buffer space, different from the discrete velocity space, but defined on the same mesh, constituted by standard polynomials with one degree less than the FE space for the velocity. This approach gives rise to a method with reduced computational cost for some choices of the interpolation operator. This method has been recently supported by a thorough numerical analysis (existence and uniqueness, stability, convergence, error estimates, asymptotic energy balance) for the nonlinear problem related to the evolution NSE, cf. [START_REF] Ahmed | Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method[END_REF][START_REF] Chacón Rebollo | Numerical analysis of penalty stabilized finite element discretizations of evolution Navier-Stokes equation[END_REF], using a semi-implicit Euler scheme for the monolithic discretization in time. In particular, the error analysis reveals a self-adapting high spatial accuracy in laminar regions of a turbulent flow that turns to be of overall optimal high accuracy if the flow is fully laminar. Numerical simulations of 3D Beltrami flow in laminar regimes [START_REF] Ahmed | Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method[END_REF] confirm this fact. This also allows to obtain an asymptotic energy balance for smooth flows.

The main contribution of this work is to propose an efficient space-time discretization of the incompressible NSE with LPS modeling of subgrid interactions by aiming at the simulation of laminar and, especially, turbulent flows in a parallel framework. In particular, the problem is first discretized in space by using a FE LPS by interpolation, and then in time by a two-step pressure-correction projection algorithm based on semi-implicit Backward Differentiation Formulas (BDF) [START_REF] Gervasio | Algebraic fractional-step schemes with spectral methods for the incompressible Navier-Stokes equations[END_REF]. Since we aim at solving large scale problems at high Reynolds numbers, the use of parallel architectures is necessary. In such parallel framework, in addition to suitable choices of the spatial and time discretization schemes, we stress the fact that the use of domain decomposition methods and efficient linear solvers and preconditioning strategies is often mandatory to make the numerical simulations computationally feasible.

In this work, we interface the proposed fully discrete scheme for the NSE with LPS modeling with HPDDM [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF], a high performance unified framework for domain decomposition methods. In particular, we use a parallel iterative linear solver based on an optimized Schwarz domain decomposition method as preconditioner [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF][START_REF] Haferssas | An additive Schwarz method type theory for Lions's algorithm and a symmetrized optimized restricted additive Schwarz method[END_REF]. In this manner, we obtain an efficient, i.e., robust and fast, solver for the HPC of laminar and turbulent flows in the open-source FE software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF] interfaced with HPDDM. A similar recent study has been performed in [START_REF] Forti | Semi-implicit BDF time discretization of the Navier-Stokes equations with VMS-LES modeling in a high performance computing framework[END_REF], where a semi-implicit BDF time discretization scheme for the NSE with residual-based VMS-LES modeling is combined together with a parallel multigrid preconditioner applied on the right and the GMRES iterative method. However, the cited study differs also because it is applied to the monolithic (coupled velocity-pressure) form of the linear system associated to the problem. Other similar studies can be found in the literature, but they always differ in some aspects with respect to the present work. In [START_REF] Colomés | Assessment of variational multiscale models for the large eddy simulation of turbulent incompressible flows[END_REF], a balancing Neumann-Neumann domain decomposition method is used for preconditioning the GMRES iterative method applied to the fully implicit monolithic system associated to a residual-based Orthogonal Subscales (OSS) modeling of the NSE. In [START_REF] Gravemeier | An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow[END_REF], algebraic multigrid strategies for VMS-LES modeling are discussed, where a Smagorinsky-type eddy viscosity model is required. In our parallel framework, we aim at further reducing the computational cost by using a semi-implicit fractional-step (FS) projection algorithm, and the computational complexity by using a LPS modeling, which is not fully residual-based, but allows similarly to obtain a high-order accuracy without requiring any ad-hoc eddy viscosity. The parallel framework set aside, a fully discrete LPS method that uses a pressure-correction scheme based on semi-implicit BDF and inf-sup stable FE has been analyzed in [START_REF] Arndt | Quasi-optimal error estimates for the incompressible Navier-Stokes problem discretized by finite element methods and pressure-correction projection with velocity stabilization[END_REF]. Here, we aim to test the practical performances of a similar method based on inf-sup stable FE in conjunction with grad-div and streamline derivative-based LPS in our parallel framework, similarly to the numerical investigation performed in [START_REF] Colomés | Mixed finite element methods with convection stabilization for the large eddy simulation of incompressible turbulent flows[END_REF] for a monolithic implicit time discretization with a recursive block preconditioning.

The proposed fully discrete scheme is tested towards the benchmark problem of recirculating flow in a lid-driven cavity. Firstly, we validate the proposed numerical scheme for the two-dimensional lid-driven cavity flow at high Reynolds number (up to Re = 10,000). Then, we address simulations of the three-dimensional lid-driven cavity flow, for which there exist experimental measurements and numerical results in the literature, by considering three significant Reynolds numbers (Re = 3200, 7500, and 10,000) to cover the spectrum from laminar to turbulent regime. A comparison of first and second-order statistics with experimental data so as to other numerical results justifies the interest of our approach: the proposed method exhibits a high-order accuracy in predicting these sensitive measures already for relatively coarse meshes, and also computational efficiency and strong scalability results of the solver in a HPC framework are showcased.

The outline of the paper is as follows. In Section 2, we introduce the model problem and its continuous variational formulation for time-dependent incompressible NSE. In Section 3, we introduce the fully discrete problem. We initially describe the proposed LPS spatial approximation of the incompressible evolution NSE, and we state its main properties. The rest of the section is devoted to the time discretization of the LPS model by means of an incremental pressure-correction algorithm with semi-implicit BDF. In Section 4, we describe the parallel solver developed for the fully discrete problem, as well as the preconditioning technique used in the framework of HPDDM. The proposed strategy is tested for the recirculating flow in a lid-driven cavity in Section 5. Firstly, we report and discuss the numerical results for the two-dimensional case at high Reynolds number. Then, we show the potential of the proposed method for simulating turbulent recirculating flow in a three-dimensional lid-driven cavity, pointing out also the parallel performances of the solver. Section 6 states the main conclusions of the paper.

Time-dependent incompressible Navier-Stokes equations: problem statement and variational formulation

We introduce an Initial-Boundary Value Problem (IBVP) for the incompressible evolution NSE. For the sake of simplicity, we just impose homogeneous Dirichlet boundary condition on the whole boundary.

Let [0, T ] be the time interval, and Ω a bounded polyhedral domain in R d , d = 2 or 3, with a Lipschitz-continuous boundary Γ = ∂Ω . The transient NSE in strong form for an incompressible fluid are given by:

find u : Ω × (0, T ) -→ R d and p : Ω × (0, T ) -→ R such that:        ∂ t u + ∇ • (u ⊗ u) -2ν∇ • D(u) + ∇ p = f in Ω × (0, T ), ∇ • u = 0 in Ω × (0, T ), u = 0 on Γ × (0, T ), u(x, 0) = u 0 (x) in Ω , (2.1)
where u ⊗ u is the tensor function of components u i u j , and D(u) is the symmetric deformation tensor given by D(u) = (1/2)(∇u + ∇u T ). The unknowns are the velocity u and the pressure p of the incompressible fluid. The data are the source term f , which represents a body force per mass unit (typically the gravity), the kinematic viscosity ν of the fluid, which is a positive constant, and the initial velocity u 0 .

To define the weak formulation of problem (2.1), we need to introduce some useful notations for spaces. We consider the Sobolev spaces H s (Ω ), s ∈ R, and L p (Ω ), 1 p ∞. We use the following notation for vector-valued Sobolev spaces: H s and L p respectively shall denote [H s (Ω )] d and [L p (Ω )] d (similarly for tensor spaces of dimension d × d). Also, the parabolic Bochner function space L p (0, T ; X ) (resp. L p (0, T ; X)), where X (resp. X) stands for a scalar (resp. vector-valued) Sobolev space shall be denoted by L p (X ) (resp. L p (X)). In order to give a variational formulation of problem (2.1), let us consider the velocity space:

H 1 0 = [H 1 0 (Ω )] d = w ∈ [H 1 (Ω )
] d : w = 0 on Γ . This is a closed linear subspace of H 1 , and thus a Hilbert space endowed with the H 1 -norm. Thanks to Korn's inequality, cf. [START_REF] Horgan | Korn's inequalities and their applications in continuum mechanics[END_REF], the H 1 -norm is equivalent on H 1 0 to the norm w H 1 0 = D(w) L 2 . Also, let us introduce the space of divergence-free functions:

H 1 0,div = w ∈ H 1 0 : ∇ • w = 0 a.e.
in Ω . The space H 1 0,div is a closed linear subspace of H 1 0 , and thus a Hilbert space endowed with the H 1 -norm. We shall consider the following variational formulation of (2.1):

given f ∈ L 2 (H -1 ) and u 0 ∈ H -1 , find u ∈ L ∞ (L 2 ) ∩ L 2 (H 1 0,div ), P ∈ L 2 (L 2 0 ) such that:                  - T 0 (u(t), v) Ω ϕ ′ (t) dt -u 0 , v ϕ(0) + T 0 [b(u(t), u(t), v) + a(u(t), v)]ϕ(t) dt + T 0 (P(t), ∇ • v) Ω ϕ ′ (t) dt = T 0 f (t), v ϕ(t) dt, (2.2) for any v ∈ H 1 0 , ϕ ∈ D([0, T ]) such that ϕ(T ) = 0, •,
• stands for the duality pairing between H 1 0 and its dual H -1 , and L 2 0 consists of L 2 -functions with zero mean in Ω . The forms b and a are given by:

b(w, u, v) = 1 2 [(w • ∇ u, v) Ω -(w • ∇ v, u) Ω ] , (2.3
)

a(u, v) = 2ν (D(u), D(v)) Ω , (2.4)
for u, v, w ∈ H 1 0 . The skew-symmetric form of the convective term b is chosen for conservation purposes: note that b(w, v, v) = 0 for all w, v ∈ H 1 0 . The physical pressure is the time derivative of the unknown P :

p = ∂ t P ∈ H -1 (L 2 0 ) = H 1 0 (0, T ; L 2 0 ) ′ .
The interest of considering P as unknown instead of p is that there are high technical difficulties to obtain uniform bounds for the discrete pressures in a Banach space of space-time functions, see [24, remark 10.2], while one obtains uniform bounds in the Banach space L ∞ (L 2 ) for the numerical approximation of P, see [12, theorem 4.3]. It is known that for domains which satisfy the cone condition, as bounded polyhedral domains, P ∈ L ∞ (L 2 ), e.g., see [25, remark 2.5]. We notice, however, that for practical computations one would approximate the physical pressure p, and P is introduced just for the numerical analysis.

Fully discrete problem

In this section, we describe both the spatial approximation and the time discretization proposed for the model problem.

Space approximation: a finite element local projection stabilization model

We consider a FE LPS method applied to the weak form of the NSE (2.2). LPS schemes were originally proposed for the Stokes problem [START_REF] Becker | A finite element pressure gradient stabilization for the Stokes equations based on local projections[END_REF], and then successfully extended to transport problems [START_REF] Ahmed | Discontinuous Galerkin time stepping with local projection stabilization for transient convectiondiffusion-reaction problems[END_REF][START_REF] Barrenechea | A local projection stabilization finite element method with nonlinear crosswind diffusion for convection-diffusion-reaction equations[END_REF][START_REF] Becker | A two-level stabilization scheme for the Navier-Stokes equations[END_REF][START_REF] Knobloch | A generalization of the local projection stabilization for convection-diffusion-reaction equations[END_REF][START_REF] Matthies | Stabilization of local projection type applied to convection-diffusion problems with mixed boundary conditions[END_REF][START_REF] Roos | Robust Numerical Methods for Singularly Perturbed Differential Equations. Convection-Diffusion-Reaction and Flow Problems[END_REF]. As classical stabilization procedures, these numerical discretizations are based upon an "augmented" variational formulation of the flow equations, which includes additional terms to the standard Galerkin discretization to provide specific stabilization of any single operator term that could be a source of instability. In particular, LPS schemes allow to circumvent the standard discrete inf-sup condition and to use equal order interpolation for velocity and pressure, and they also provide local stabilization of convection-dominant effects and improvement of local mass conservation [START_REF] Linke | Collision in a cross-shaped domain-a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD[END_REF]. Different variants of LPS methods have been investigated during the recent years for incompressible flow problems. The main common feature is that, thanks to local projection, the symmetric stabilization terms only act on the small scales of the flow, thus ensuring a higher accuracy with respect to more classical stabilization procedures, such as penaltystabilized methods, cf. [START_REF] Rebollo | A term by term stabilization algorithm for finite element solution of incompressible flow problems[END_REF]. This also guarantees a self-adapting high accuracy in laminar regions of a turbulent flow, which turns out to be of overall optimal high accuracy if the flow is fully laminar, and allows to obtain an asymptotic energy balance for smooth flows [START_REF] Ahmed | Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method[END_REF]. Thus, the effect of LPS is on the one hand to improve the convergence to smooth solutions. On the other hand, for rough solutions, LPS limits the propagation of perturbations generated in the vicinity of sharp gradients, potentially maintaining these schemes as suitable and useful tools for the simulation of turbulent flows. Moreover, an important advantage of their term-by-term structure is that the projection can be easily treated as implicit, without having all the residual terms coupled, as for more complex residual-based VMS methods, cf. [START_REF] Bazilevs | Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows[END_REF][START_REF] Codina | Time dependent subscales in the stabilized finite element approximation of incompressible flow problems[END_REF]. Actually, LPS schemes can be viewed as simplifications of residual-based VMS methods [START_REF] Ahmed | A review of variational multiscale methods for the simulation of turbulent incompressible flows[END_REF], since they are not fully consistent (only specific dissipative interactions are retained), but are of optimal order with respect to the FE interpolation. For a detailed description of different variants of LPS schemes, we refer to [START_REF] Knobloch | Local projection stabilization for advection-diffusion-reaction problems: one-level vs. two-level approach[END_REF][START_REF] He | The two-level local projection type stabilization as an enriched one-level approach[END_REF][START_REF] Tobiska | The two-level local projection type stabilization as an enriched one-level approach. A one-dimensional study[END_REF].

In order to describe in detail the spatial approximation of the model problem (2.2), let {T h } h>0 be a family of affine-equivalent, conforming (i.e., without hanging nodes) and regular triangulations of Ω , formed by triangles or quadrilaterals (d = 2), tetrahedra or hexahedra (d = 3). For any mesh cell K ∈ T h , its diameter will be denoted by

h K and h = max K ∈T h h K .
Given a positive integer l and a mesh cell K ∈ T h , denote by R l (K ) either P l (K ) (i.e., the space of Lagrange polynomials of degree l, defined on K ), if the grids are formed by triangles (d = 2) or tetrahedra (d = 3), or Q l (K ) (i.e., the space of Lagrange polynomials of degree l on each variable, defined on K ), if the family of triangulations is formed by quadrilaterals (d = 2) or hexahedra (d = 3). We consider the following FE spaces for the velocity:

   Y l h = V l h (Ω ) = {v h ∈ C 0 (Ω ) : v h | K ∈ R l (K ), ∀K ∈ T h }, Y l h = [Y l h ] d = {v h ∈ [C 0 (Ω )] d : v h | K ∈ [R l (K )] d , ∀K ∈ T h }, X h = Y l h ∩ H 1 0 . (3.1)
Hereafter, Y l h (resp., Y l h ) will constitute the discrete foreground vector-valued (resp. scalar) spaces in which we will work on.

We initially approximate the weak formulation (2.2) of the initial-boundary value problem (2.1) for the incompressible evolution NSE by a high-order term-by-term stabilization procedure in space [START_REF] Chacón Rebollo | A high order term-by-term stabilization solver for incompressible flow problems[END_REF]. The stabilization effect is achieved by adding least-squares terms that give a weighted control on the fluctuations of the quantity of interest, based upon a specific locally stable projection or interpolation operator on a continuous buffer space. This provides an efficient discretization with a reduced computational cost that keeps the same high-order accuracy with respect to standard projection-stabilized methods. We initially state the proposed LPS discretization as:

find (u h , p h ) : (0, T ) → X h × M h such that:    (∂ t u h , v h ) Ω + b(u h , u h , v h ) + a(u h , v h ) -( p h , ∇ • v h ) Ω + s conv (u h , u h , v h ) + s div (u h , v h ) = f , v h , (∇ • u h , q h ) Ω + s pres ( p h , q h ) = 0, (3.2) for any (v h , q h ) ∈ X h × M h , where M h = Y l h ∩ L 2 0 .
The forms s conv , s div and s pres in (3.2) correspond to a high-order term-by-term stabilized method, cf. [START_REF] Chacón Rebollo | A high order term-by-term stabilization solver for incompressible flow problems[END_REF], and are given by:

s conv (u h , w h , v h ) = K ∈T h τ m,K (σ * h (u h • ∇w h ), σ * h (u h • ∇v h )) K , (3.3
)

s div (u h , v h ) = K ∈T h τ d,K (σ * h (∇ • u h ), σ * h (∇ • v h )) K , (3.4
)

s pres ( p h , q h ) = K ∈T h τ m,K (σ * h (∇ p h ), σ * h (∇q h )) K . (3.5)
The aims of the stabilization terms is to prevent spurious instabilities due to dominant convection (s conv ), achieve additional control on the incompressibility condition by means of pressure subscales effect (s div ), and ensure infsup stability when using equal order interpolation for velocity and pressure (s pres ). Here, τ m,K , τ d,K are stabilization coefficients for convection-pressure gradient and divergence, respectively, and σ * h = I dσ h , where I d is the identity operator, and σ h is some locally stable projection or interpolation operator from L 2 on the foreground vector-valued space Y l-1 h (also called "buffer space" in this context), satisfying optimal error estimates. In practical implementations, we choose σ h as a Scott-Zhang-like [START_REF] Scott | Finite element interpolation of non-smooth functions satisfying boundary conditions[END_REF] linear interpolation operator in the space Y l-1 h , implemented in the software FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. This interpolant may be defined as:

∀x ∈ Ω , σ h (v)(x) = a∈N Π h (v)(a)ϕ a (x),
where N is the set of Lagrange interpolation nodes of Y l-1 h , ϕ a are the Lagrange basis functions associated to N , and Π h is the interpolation operator by local averaging of Scott-Zhang kind, which coincides with the standard nodal Lagrange interpolant when acting on continuous functions (cf. [START_REF] Chacón Rebollo | A high order term-by-term stabilization solver for incompressible flow problems[END_REF], section 4). This is an interpolant that just uses nodal values, and so is simpler to work out and more computationally efficient than the variant of the Scott-Zhang operator introduced in [START_REF] Badia | On stabilized finite element methods based on the Scott-Zhang projector. Circumventing the inf-sup condition for the Stokes problem[END_REF] for the Stokes problem, which is instead an operator defined from a node-to-element map and requires integration on mesh elements. According to the structure of s conv , the following operators have to be composed:

• bilinear form τ discretized as C:

τ : dcY l-1 h × dcY l-1 h → R (u h , v h ) → K ∈T h τ m,K (u h , v h ) K ,
• linear operator σ discretized as P:

σ : dcY l-1 h → dcY l-1 h v h → σ * h (v h ), • linear operator D discretized as D: D : X h → dcY l-1 h v h → u h • ∇v h ,
where dcY l-1 h stands for the discontinuous version of the buffer space Y l-1 h . In practice, we express D as the sum of two (resp. three) operators in 2D (resp. 3D). This is done by considering v h component-by-component, i.e., D = D x + D y (resp. D = D x + D y + D z ). The discretization of (3.3) is thus the symmetric local matrix S defined as:

S = (P • D) T C(P • D). (3.6)
In formula (3.4), σ h denotes an operator between the scalar spaces L 2 and Y l-1 h , but we use the same notation for the sake of simplicity. Actually, if needed, specific stabilizations for convection, divergence and pressure gradient may be used, through different approximation operators. When using inf-sup stable FE (as in the case of the computations performed in this work), there is no need to consider the pressure stabilization term s pr es , so this term is neglected in this case. However, the term s div seems to be crucial in this case, due to the poor resolution of the pressure typical of mixed interpolations that satisfy the inf-sup condition [START_REF] Olshanskii | Grad-div stabilization and subgrid pressure models for the incompressible Navier-Stokes equations[END_REF]. For this term, just a pure grad-div penalty stabilization is used in practical implementations, i.e., σ * h = I d in formula (3.4), to reduce the computational cost. However, it can be seen from numerical analysis that this approximation does not introduce any consistency error, that is it does not affect the optimal accuracy of the method [START_REF] Ahmed | Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method[END_REF]. The working expressions of the stabilization coefficients are:

τ m,K = d c 1 ν (h K /l) 2 + c 2 U n K (h K /l) -1 , (3.7)
τ d,K = (h K /l) 2 d c 1 τ m,K , (3.8)
by following the form proposed in [START_REF] Codina | Time dependent subscales in the stabilized finite element approximation of incompressible flow problems[END_REF][START_REF] Codina | Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales[END_REF], designed by asymptotic scaling arguments applied in the framework of stabilized methods aimed at taking into account the local balance between convection and diffusion. In expressions (3.7) and (3.8), d is the dimension of the problem, c 1 and c 2 are user-chosen positive constants, l is the polynomial degree of the velocity FE approximation, and U K is some local speed on the mesh cell K . The values of the constants c 1 and c 2 are chosen to be

c 1 = 4, c 2 = √ c 1 = 2, cf. [42]. If U K ∈ L ∞ (K )
, the following technical hypothesis on the stabilization coefficients required to perform the numerical analysis is ensured:

Hypothesis 3.1. The stabilization coefficients τ m,K , τ d,K satisfy the following conditions:

α 1 h 2 K τ m,K α 2 h 2 K , 0 < τ d,K β, (3.9)
for all K ∈ T h , and some positive constants α 1 , α 2 , β independent of h.

In practical implementations,

U n K = u n h L 2 (K ) /|K | 1/2
, with |K | denoting the surface (d = 2) or volume (d = 3) of element K , considering that in practice u h solves the NSE. Also, in three dimensions, we define the element size

h K = 3 √ |K |
, for all K ∈ T h , for simplicity of implementation. Note that when dealing with equal-order interpolations, a specific discrete inf-sup condition that is essential for the stability of the proposed method can be derived: Lemma 3.2. Assume that Hypothesis 3.1 holds. Then, for a uniformly regular family of triangulations {T h } h>0 , we have the following inf-sup condition:

∀q h ∈ M h , q h L 2 C sup v h ∈X h (∇ • v h , q h ) Ω D(v h ) L 2 + σ * h (∇q h ) τ m , (3.10)
for some positive constant C independent of h, where τ m denotes here the weighted L 2 -norm with stabilization coefficient τ m,K .

The proof of this lemma can be derived from [START_REF] Chacón Rebollo | A high order term-by-term stabilization solver for incompressible flow problems[END_REF], where it is also shown that the discrete inf-sup condition (3.10) can be extended to a more complex condition that holds for a simply regular family of triangulations.

Remark 3.3. The presented high-order term-by-term stabilization procedure in space by using a Scott-Zhang-like interpolation operator has been extended to geophysical flows governed by the primitive equations of the ocean [START_REF] Chacón Rebollo | Numerical approximation of the Smagorinsky turbulence model applied to the primitive equations of the ocean[END_REF] and buoyant flows governed by the Boussinesq equations [START_REF] Chacón Rebollo | A high-order local projection stabilization method for natural convection problems[END_REF], giving high-order accuracy for relatively coarse grids. Also, it has been combined with a VMS-Smagorinsky term and wall laws for the accurate simulation of turbulent boundary layers in [START_REF] Chacón Rebollo | Finite element approximation of an unsteady projection-based VMS turbulence model with wall laws[END_REF][START_REF] Chacón Rebollo | Numerical analysis of a finite element projection-based VMS turbulence model with wall laws[END_REF][START_REF] Rubino | Numerical Modeling of Turbulence by Richardson Number-Based and VMS Models[END_REF].

Time discretization: incremental pressure-correction algorithm with semi-implicit backward differentiation formulas

To state the proposed time discretization of the unsteady LPS model described by system (3.2), consider a positive integer number N and define ∆t = T /N , t n = n∆t, n = 0, 1, . . . , N . We compute the approximations u n h and p n h to u n = u(•, t n ) and p n = p(•, t n ), respectively, by using an incremental pressure-correction scheme based on semi-implicit BDF, for which the nonlinear terms are extrapolated by means of Newton-Gregory backward polynomials [START_REF] Cellier | Continuous System Modeling[END_REF]. Let us denote by u n h an intermediate approximate velocity at time t n . In order to abbreviate its discrete time derivative, we define the operator D t by:

D r t u h n = α r u h n+1 -u h,r n ∆t , (3.11)
where for BDF schemes of orders r = 1, 2 we have:

u n h,r =    u n h if n 0, for r = 1 (BDF1) 2 u n h - 1 2 u n-1 h if n 1, for r = 2 (BDF2) (3.12)
and

α r =    1, for r = 1 (BDF1) 3 2 , for r = 2 (BDF2). (3.13)
We consider the following extrapolations [START_REF] Cellier | Continuous System Modeling[END_REF] of order r = 1, 2 for the intermediate convection velocity:

u n, * h,r = u n h if n 0, for r = 1 (BDF1) 2 u n h -u n-1 h if n 1, for r = 2 (BDF2) (3.14)
and the pressure:

p n, * h,r =    p n h if n 0, for r = 1 (BDF1) 1 3 (7 p n h -5 p n-1 h + p n-2 h ) if n 1, for r = 2 (BDF2) (3.15)
using in the last case the convention

p n-1 h = p n-2 h
for n = 1. In this way, after applying a standard incremental pressure-correction approach (cf. [START_REF] Codina | A stabilized finite element method for generalized stationary incompressible flows[END_REF][START_REF] Guermond | Un résultat de convergence d'ordre deux en temps pour l'approximation des équations de Navier-Stokes par une technique de projection incrémentale[END_REF] for derivation in the case of equal-order and inf-sup stable FE, respectively) to system (3.2), the fully discrete semi-implicit formulation consists in solving, for n = 0, . . . , N -1, the two-step algorithm:

find u n+1 h ∈ X h such that: D r t u n h , v h Ω + b( u n, * h,r , u n+1 h , v h ) + a( u n+1 h , v h ) + s conv ( u n, * h,r , u n+1 h , v h ) + s div ( u n+1 h , v h ) = f n+1 , v h + p n, * h,r , ∇ • v h Ω , (3.16)
for any v h ∈ X h , where f n+1 is the average value of f in [t n , t n+1 ], and:

find p n+1 h ∈ M h such that: ∇( p n+1 h -p n h ), ∇q h Ω + s pres ( p n+1 h , q h ) = - α r ∆t ∇ • u n+1 h , q h Ω n • ∇( p n+1 h -p n h ) | Γ = 0, (3.17)
for any q h ∈ M h , where n is the outer normal to Γ . The final velocity can then be recovered according to:

u n+1 h = u n+1 h - ∆t α r ∇( p n+1 h -p n h ). (3.18)
This is the temporal approach used in practical implementations. In particular, for the first time step (n = 0) we use a BDF1 scheme (r = 1) to initialize the algorithm with u 0 h = u 0 h and p 0 h some stable approximations to u 0 and p 0 , respectively. Note that this scheme coincides with a Backward Euler method. Then, a BDF2 scheme (r = 2) is applied for n 1.

In order to achieve a global (space-time) second-order accuracy for velocities, P 2 FE are used to approximate them in space. Thus, in the expressions (3.7) and (3.8) of the stabilization coefficients, l = 2. Also,

U n K = u n, * h,r L 2 (K ) /|K | 1/2 .
Note that the term s pres in Eq. (3.17) is neglected when dealing with inf-sup stable interpolations. This obviously leads to cheaper (amortized setup) Poisson solves for the pressure equation, maintaining a reasonable balance between accuracy and computational complexity. This is why, in the practical computations of this article, we choose mixed FE for which the pressure is approximated with P 1 FE (i.e., Taylor-Hood FE pairs are considered).

Remark 3.4. The semi-implicit discretization in time segregating velocity and pressure through a standard incremental time-splitting helps to construct an efficient linear solver to the NSE system for the LES of turbulent flows. In the first step Eq. (3.16), a convection-dominated convection-diffusion-reaction problem for the intermediate velocity must be solved. The second step Eq. (3.17) consists of a (stabilized) pressure-Poisson problem. Both steps are solved by using a domain decomposition method preconditioner with the GMRES iterative method applied to the associated system in the parallel framework proposed in Section 4, and a convincing strong scaling analysis of the used algorithm is showcased in Section 5.2.1.

Parallel strategy for the fully discrete problem: interface with HPDDM

The fully discrete problem described in Section 3 gives rise to two large linear systems of the form A v u = f v for subproblem (3.16) and A p q = f p for subproblem (3.17), respectively. Solving such large linear systems could become extremely expensive from the computational point of view, that is why we adopt in the numerical implementation a highly parallel strategy based on domain decomposition methods.

To describe it, we first divide the mesh {T h } h>0 in N non-overlapping meshes {T i } 1 i N using standard graph partitioners, e.g., ParMETIS [START_REF] Karypis | Multilevel k-way partitioning scheme for irregular graphs[END_REF]. If δ is a positive integer, the overlapping decomposition T δ i 1 i N is defined recursively as follows: T δ i is obtained by including all elements of T δ-1 i plus all adjacent elements of T δ-1 i . For δ = 0,

T δ i = T i . Let X δ i 1 i N (resp. M δ i 1 i N ) be the local velocity (resp. pressure) FE spaces defined on T δ i 1 i N . Now, consider the restrictions R v i 1 i N (resp. R p i 1 i N ) from X h (resp. M h ) to X δ i 1 i N (resp. M δ i 1 i N ),

and two local partitions of unity D

v i 1 i N and D p i 1 i N such that: N j=1 R v T j D v j R v j = I d n v ×n v N j=1 R p T j D p j R p j = I d n p ×n p ,
where I d denotes the identity matrix and n v (resp. n p ) is the global number of unknowns in the velocity (resp. pressure) space. Algebraically speaking, if n v is the global number of velocity unknowns and {n v i } 1 i N are the numbers of degrees of freedom in each local velocity FE space, then R v i is a Boolean matrix of size n v i × n, and

D v i is a diagonal matrix of size n v i × n v i , for all 1 i N .
Of course, this also holds for the operators related to the pressure spaces.

Using the partition of unity, one can define the following one-level preconditioner as an extension of the Restricted Additive Schwarz (RAS) method proposed by Cai and Sarkis [START_REF] Cai | A restricted additive Schwarz preconditioner for general sparse linear systems[END_REF]:

M v -1 ORAS = N i=1 R v T i D v i B v -1 i R v i , (4.1)
where the B v i 1 i N are local operators that resemble the submatrices

R v i A v R v T i 1 i N
, but with more efficient transmission conditions between subdomains, e.g., see [START_REF] Dolean | Optimized Schwarz methods for Maxwell's equations[END_REF]. In practice, we use simple Robin boundary conditions. This produces a rather efficient preconditioner for the successive linear systems related to the velocity unknowns even at large scale, when considering thousands of subdomains and millions of degrees of freedom, cf. Section 5.2.1. Note that the linear systems related to the velocity Eq. (3.16) are non-symmetric. Their sparsity pattern is however the same throughout the simulation as soon as the time discretization is properly initialized, as proposed in Section 3.2. Since we use in all simulations direct solvers to compute the action of the inverse of the local B v i 1 i N , needed in Eq. ( 4.1), we can perform a single symbolic factorization after startup, followed by as many numeric factorizations as needed.

The same formalism cannot exactly be used as it is for solving the linear systems related to the pressure Eq. (3.17). First, note that for these systems only the right-hand sides are changing when using inf-sup stable FE, as it is the case in the present computations, for which the pressure stabilization term s pres in Eq. (3.17) is neglected. We thus need to assemble a single matrix and setup a preconditioner once. These costs will then be amortized over the many successive time steps. Even with more complex transmission conditions between subdomains, a so-called "one-level" method such as M v -1

ORAS defined in Eq. (4.1) would not scale properly in terms of number of iterations for solving a pressure-Poisson problem. In this case, we thus chose to use a more sophisticated multilevel domain decomposition method using the GenEO-2 approach [START_REF] Haferssas | An additive Schwarz method type theory for Lions's algorithm and a symmetrized optimized restricted additive Schwarz method[END_REF]. This preconditioner, M p -1 GenEO-2 , uses a spectral coarse grid to better couple all subdomains. Thus, as clearly shown in Section 5.2.1, the number of iterations for solving the linear systems related to the pressure unknowns remains really low even when considering large number of subdomains.

The applied domain decomposition methods are implemented in the library HPDDM [START_REF] Dolean | An Introduction to Domain Decomposition Methods: Algorithms, Theory, and Parallel Implementation[END_REF][START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF] (more details about the implementation of the cited spectral preconditioners are given in [START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF]). HPDDM can be interfaced with various programming languages and open-source FE libraries such as FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF], the one used in the present simulations.

Numerical studies

In this section, we present two numerical experiments respectively based on a two-and a three-dimensional lid-driven cavity. The former has the aim of validating the proposed method, the latter of evaluating its numerical performances at laminar, transient, and turbulent regimes, also on massive parallel settings.

2D lid-driven cavity flow

The two-dimensional lid-driven cavity test is one of the most popular validation problems for fluid flow simulations. In this test, the fluid is contained in a unitary squared domain and it has Dirichlet boundary conditions on all sides: three stationary sides and one moving side at the top, characterized by a unitary tangent velocity. The right-hand side of the momentum equation vanishes in Ω , i.e., f = 0.

Three different flow conditions have been simulated, respectively for Reynolds numbers Re = 1,000, 5,000, 10,000 (Re = 1/ν). In all cases, an impulsive start is performed, i.e., the initial condition is a zero velocity field, and the time step is ∆t = 0.01. All computations are carried out with inf-sup stable FE of Taylor-Hood-type, i.e., P 2 -P 1 FE are used to approximate the velocity-pressure pair on grids refined towards the walls in both spatial directions using the hyperbolic tangent function:

f (x) = 0.5 1 + tanh[2(2x -1)] tanh (2)
.

(5.1)

The partition of the cavity is 32 2 , 50 2 , 64 2 respectively for the Reynolds numbers Re = 1,000, 5,000, 10,000. Numerically, a steady solution seems to exist up to Re = 7,500 approximately. Bifurcations to periodic unsteady solutions, obtained by DNS simulations, are reported by many authors, e.g., [START_REF] Auteri | Numerical investigation on the stability of singular driven cavity flow[END_REF][START_REF] Tiesinga | Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method[END_REF]. So, for Re = 1,000 and 5,000, we assume that a stationary solution is reached if the velocity deviation between two consecutive time steps is lower than a chosen tolerance, fixed as:

u n+1 h -u n h L 2 10 -6 .
For Re = 10,000, we just let run the simulation till a final time T = 1,000. To compare our results, we consider those of Erturk et al. [START_REF] Erturk | Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers[END_REF] obtained using the streamfunction-vorticity formulation on a very fine uniform grid mesh of 601 2 . The most interesting feature is the genesis of a varying number of vortices inside the cavity. Fig. 1 shows the streamlines on a colored velocity distribution for all the Reynolds numbers considered, Re = 1,000, 5,000, 10,000, at the final computation time. We observe that secondary and tertiary vortices are well-resolved. In Fig. 2, we show cross-sections of the horizontal and vertical final velocities, respectively at y = 0.5 and x = 0.5, for all the Reynolds numbers considered. Hereafter, in the legends of figures, we denote by BDF2 FS-LPS our proposed method. It can be observed that the results for the used coarse grids are in agreement with the ones obtained by Erturk et al. [START_REF] Erturk | Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers[END_REF], also in the boundary layers, thus validating the method.

3D lid-driven cavity flow

In this section, the three-dimensional lid-driven cavity test is performed to investigate the numerical performances of the proposed solver at laminar, transient, and turbulent regimes, also on massive parallel settings. Similarly to the two-dimensional problem, the three-dimensional case is characterized by a fluid flow in a cubic domain driven by a tangential unitary velocity along one of the six boundary surfaces. Homogeneous Dirichlet conditions are adopted on all the other boundaries.

The recirculating flow in a 3D lid-driven cavity presents the occurrence of some considerable 3D features, even at relatively low Reynolds numbers. One of the most remarkable is the formation of Taylor-Görtler-like (TGL) vortices at the corners of the bottom of the cavity. Small counter-rotating vortices are formed as a result of the curvature of the streamlines due to the main vortex in the middle of the cavity [START_REF] Schlichting | Boundary-Layer Theory[END_REF], similarly to the 2D case, but with a genuine 3D structure. Zang et al. [START_REF] Zang | A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows[END_REF] report the numerical simulation of the 3D cavity flow by a LES model using the Finite Volume Method (FVM), using the dynamic procedure of Germano et al. [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF]. Based on experimental experiences performed in Prasad and Koseff [START_REF] Prasad | Reynolds number and end-wall effects on a lid-driven cavity flow[END_REF], they describe the flow at Reynolds number 3,200 to be essentially laminar, although an inherent unsteadiness may occur. For Reynolds number 7,500, a transitional stage is reached, since the flow becomes unstable near the downstream eddies at Reynolds numbers higher than about 6,000. With even higher Reynolds numbers at about 10,000, the flow becomes fully turbulent. Thus, laminar, transient, and turbulent regimes are covered by the choice of these three cases.

The primary goal of the simulation of the 3D cavity flow is to obtain a bounded kinetic energy as time increases, during the complete simulation time needed to reach a stable equilibrium [START_REF] Moin | Advances in large eddy simulation methodology for complex flows[END_REF]. This may look as a simple requirement, but some turbulence models violate it. Indeed, Iliescu et al. [START_REF] Iliescu | A numerical study of a class of LES models[END_REF] reported the numerical results obtained with three subgrid scale models for Re = 10,000: the Smagorinsky model [START_REF] Smagorinsky | General circulation experiment with the primitive equations: I. The basic experiment[END_REF], a traditional Taylor LES model of Clark et al., [START_REF] Clark | Evaluation of subgrid-scale models using an accurately simulated turbulent flow[END_REF], and two variants of a new rational LES model developed in Galdi and Layton [START_REF] Galdi | Approximation of the larger eddies in fluid motions. II. A model for space-filtered flow[END_REF]. It was shown that the Taylor LES model produces an energy blow-up in finite time. The two rational LES models did not cause an energy blow-up, but exhibited important oscillations. The standard Smagorinsky model (with the Smagorinsky constant C S = 0.1) turned out to be notably more diffusive, as expected.

Following the works of Gravemeier et al. [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF] and Zang et al. [START_REF] Zang | A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows[END_REF], we simulate the 3D cavity flow at Reynolds numbers Re = 3,200, 7,500, 10,000. In [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF], Gravemeier et al. analyze the performances of a VMS model based on the Residual Free Bubbles (RFB) method and the use of several nested meshes. This is called the three-level FE method (VMS-3L), as it includes three grid levels, and it takes into account the effect of small unresolved scales onto small resolved scales by a subgrid eddy viscosity approach in Smagorinsky's form.

A first difficulty we face in the numerical simulations is to obtain a high-order accuracy with a relatively coarse basic discretization (i.e., low computational cost), for all flow situations. The computational grid consists of a 32 3 partition of the unit cube, uniform in the y-direction, and refined towards the walls in both the x-and z-directions using the hyperbolic tangent function (5.1), in order to handle large velocity gradients. This already provides a large improvement in the accuracy of the numerical results. On this mesh, we consider inf-sup stable three-dimensional FE. In particular, Taylor-Hood P 2 -P 1 FE are used to approximate the velocity-pressure pair (no need to consider pressure stabilization). In any case, we are considering a number of degrees of freedom comparable to the one of the VMS-3L numerical simulation of Gravemeier et al. [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF]. For all Reynolds numbers considered, an impulsive start is performed, i.e., the initial condition is a zero velocity field, and ∆t = 0.1, as in [START_REF] Iliescu | A numerical study of a class of LES models[END_REF][START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF]. A characteristic time scale T cav is defined in Zang et al. [START_REF] Zang | A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows[END_REF] to be the estimated time for a fluid particle at the edge of the top boundary layer to turn and (approximately) reach its starting position in the cavity. This time scale is roughly estimated to be about 10 time units for the current computations. Initially, the simulation is run for five time scales T cav , i.e., 50 time units or 500 time steps. Within this time period, the flow is expected to develop to full extent [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF], including a subsequent relaxation time. Afterwards, statistics (quasi-steady results) are collected for another five time scales T cav .

When available, results are graphically compared to the experimental data of Prasad and Koseff [START_REF] Prasad | Reynolds number and end-wall effects on a lid-driven cavity flow[END_REF], and numerical results of Gravemeier et al. [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF]. The experimental data for the flow at Reynolds number Re = 7,500 have only been evaluated for half of the cavity.

The temporal evolution of the total kinetic energy subject to:

E kin (u n h ) = 1 2 Ω u n h • u n h dx,
for the three cases Re = 3,200, 7,500, and 10,000 is displayed in Fig. 3. The flows become roughly stationary at t ≃ 5 T cav (i.e., at about 50 time units), as expected. These results are almost comparable with the ones obtained by the VMS-3L method in Gravemeier et al. [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF].

The mean velocities u 1 and u 3 are computed as a discrete time average according to:

u i (x) = 1 N /2 N -1 n=N /2 u i (x, t n ), i = 1, 3, N = # time steps = 1,000.
Fig. 4 shows the mean velocity u 1 on the centerline z = 0.5 of the longitudinal mid-plane y = 0.5, for the various Reynolds numbers under consideration. The proposed method is in agreement with the experimental data of Prasad and Koseff [START_REF] Prasad | Reynolds number and end-wall effects on a lid-driven cavity flow[END_REF], even with the coarse basic discretization at hand, and performs better than the VMS-3L method [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF]. A similar accuracy is found for the mean velocity u 3 on the centerline x = 0.5 of the longitudinal mid-plane y = 0.5, see Fig. 5. Also, higher-order moments ũn , with n > 1 and ũ denoting the fluctuating part of u, are achieved by collecting values in the sense of a discrete time average, which is an appropriate procedure for stationary turbulence. In particular, we have considered the variance (n = 2) for the first and third component of the velocity, that reads

ũ2 i = u 2 i -u i 2 ,
with the standard deviation (root mean square, RMS) defined as ũ2 i (i = 1, 3). Also, the off-diagonal component ũ1 ũ3 = u 1 u 3u 1 u 3 of the Reynolds stress tensor is depicted, and finally the turbulent kinetic energy (TKE) defined as 1 i is reported, together with the corresponding energy spectrum. As in Prasad and Koseff [START_REF] Prasad | Reynolds number and end-wall effects on a lid-driven cavity flow[END_REF], the RMS values and the off-diagonal Reynolds stress component are multiplied by the amplification factors 10 and 500, respectively, in order to ensure a reasonable visual impression of these values within the respective graphs, and analogously the TKE is multiplied by the amplification factor 100. With respect to the experimental data, larger errors appear for the RMS values and the crossed component of the Reynolds stress tensor. These deviations are shown for Re = 10,000 in Figs. 6 and7. Mispredictions of various peaks of these curves may also be found in the numerical results of Zang et al. [START_REF] Zang | A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows[END_REF], achieved with a two-time finer discretization in every coordinate directions. This underlines the difficulty in predicting these sensitive measures. For completeness, we also report the TKE obtained with our scheme on the centerlines of the mid-plane y = 0.5 for Re = 10,000 in Fig. 8, although experimental data and numerical results from Gravemeier et al. [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF] are not available for this metric. From the TKE, we derive the corresponding energy spectrum on the centerlines of the mid-plane y = 0.5 for Re = 10,000 in Fig. 9, where we denote by E(k) the energy spectrum, and k is the wavenumber. The slope of the energy spectrum in the inertial subrange is in line with the one of the theoretical Kolmogorov curve k -5/3 . Qualitatively, we have observed that the flow exhibits effectively the formation of three-dimensional TGL corner vortices at the cavity end walls, that interact with the primary circulation vortex, thus influencing the distribution of momentum within the entire cavity, see Fig. 10. In the case Re = 3200, in accordance to Prasad and Koseff [START_REF] Prasad | Reynolds number and end-wall effects on a lid-driven cavity flow[END_REF], it is possible to discern these vortices as organized structures, while for higher Re, increasing turbulent effects cause the breakdown of these organized structures, resulting in a "weaker" flow when compared with the pure two-dimensional flow. This suggests that the high-frequency turbulent fluctuations become dominant, and they partially destroy the integrity (or coherence) of the TGL vortices.

Parallel performances of the solver

The goal of this section is to assess the parallel efficiency of the proposed method. In particular, we are interested in the strong scalability of the described algorithm. In all that follows, we start by generating a global mesh using 100 grid points in each direction, distributed in the x-and z-directions using the hyperbolic tangent function (5.1). The mesh is then decomposed using ParMETIS [START_REF] Karypis | Multilevel k-way partitioning scheme for irregular graphs[END_REF]. The operations related to the discretization of continuous operators are performed using FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]. The linear solvers and preconditioners are implemented in HPDDM [START_REF] Jolivet | Scalable domain decomposition preconditioners for heterogeneous elliptic problems[END_REF]. Internally, the underlying linear solver for all local subdomain solves is Intel MKL PARDISO [START_REF]MKL web page[END_REF][START_REF] Schenk | Solving unsymmetric sparse systems of linear equations with PARDISO[END_REF]. The stabilization terms from Eq. (3.6) have to be evaluated at each time step but involve purely local computations. In our experiments, the velocity field was initially set to zero. Thus, it takes one time step for the sparsity pattern of S to remain the same. We exploit this property by computing the matrix-matrix-matrix product first by using the domain-specific language of FreeFem++, and then by using PETSc [START_REF] Balay | PETSc web page[END_REF] and its MatPtAPNumeric1 routine (where the symbolic product may be bypassed). The matrix-matrix product P • D, with D expressed as a sum of two or three matrices, is computed using a dedicated routine. Results were obtained on Curie, a system composed of 5,040 nodes with two eight-core Intel Sandy Bridge clocked at 2.7 GHz. The interconnect is an InfiniBand QDR full fat tree and the MPI implementation exploited was IntelMPI version 2017.0.2.174. All binaries and shared libraries were compiled with Intel compilers and Math Kernel Library support (for dense linear algebra computations). We used from 512 up to 16,384 MPI processes with a single OpenMP thread per process (flat MPI parallelism). Preconditioners are thus defined with as many subdomains as the number of MPI processes.

We used a test case in the turbulent regime with Re = 15,000 and ∆t = 0.05. The successive discretizations of Eq. (3.16) yield linear systems A v u = f v with 24.7×10 6 unknowns, while the single discretization of Eq. (3.17) yields a linear system A p q = f p with 1.1×10 6 unknowns. Both systems are left-preconditioned by either M v -1

ORAS or M p -1

GenEO-2 and the GMRES method is stopped when the relative preconditioned residual is lower than 10 -8 for the velocity unknowns and 10 -6 for the pressure unknowns. Since we have observed in the two previous sections that the proposed solver is stable throughout the time steps, we will consider in the present scalability analysis only the 10th time step. This means that the startup phase is over. Since the single construction of M p -1 GenEO-2 is amortized over time, it will not be included in the present analysis.

In order to give a complete overview of the performances for the proposed method, we first represent in Fig. 11 the total time to complete the 10th time step. Clearly, the implementation scales very well. This may be explained by multiple facts. As displayed in Fig. 12, the preconditioners are indeed both numerically extremely stable, with numbers of iterations remaining in the same low range. In Table 1, we report the time spent in all subroutines of the 10th time step. The first column (N ) represents the number of subdomains (or MPI processes), to which most of subroutines scale almost linearly. The second column (stabilization) is the time needed by FreeFem++ and PETSc to assemble the stabilization terms. The third column (optimized operators) is the sum of the time spent assembling the velocity linear system and adding the optimized boundary conditions used for the preconditioner M v -1

ORAS . The fourth column (BDF2 & aux. tasks) represents the amount of time spent incrementing the time step, mostly vector summations, as well as auxiliary tasks involving some communications, e.g., updating ghost elements. Since we are using exact LU

Table 1

Breakdown of the time spent in the various subroutines involved in the scalability analysis of Fig. 11 for a single time step. A v and A p are of order 24 millions and 1 million, respectively.

N

Stabilization

Optimized operators

BDF2 & aux. tasks Setup M v -1 ORAS Solve A v u = f v

Solve

A p q = f p Total 512 decompositions for subdomain solves, we can see that the fifth column (setup M v -1 ORAS ) exhibits a super-linear speedup. It is also worthwhile to note that the solution phase for the pressure unknowns (seventh column) does not scale as the solution phase for the velocity unknowns (sixth column). Indeed, because of the use of a multilevel preconditioner, the volume of communication tends to outgrow the number of local computations, which decreases in a strong scaling experiment. There is no such phenomenon for the velocity unknowns since we are using a simpler one-level decomposition with much less communication.

Fig. 13, we display the relative time percent) spent in each subroutines of Table 1. The first three colors represent tasks that are purely concurrent and that do not involve any kind of communication. The last three colors represent tasks with some communication between processes. For example, the GMRES method requires global reductions at each inner iteration for computing scalar products. Overall, we can see that even for 16,384 processes, most of the time (around 80% of the total) of a single time step is spent on local computations. This explains why the overall scalability of the proposed method is satisfactory on a wide range of process counts.

To sum up, the parallel performances of the proposed method are rather satisfactory, and seem to be in accordance with the current state-of-the-art, e.g., [START_REF] Forti | Semi-implicit BDF time discretization of the Navier-Stokes equations with VMS-LES modeling in a high performance computing framework[END_REF].

Summary and conclusions

In this work, we have proposed a FE LPS spatial approximation of the NSE, combined with an efficient velocitypressure segregation, using semi-implicit BDF in time. In order to face the computational complexity of the problem, we have developed an ad-hoc parallel solver for the proposed method, based on HPDDM.

We have tested the proposed numerical scheme by solving the benchmark problem of the recirculating flow in a lid-driven cavity at high Reynolds numbers. This numerical study shows that the solver is able to reproduce first and second-order statistics up to a turbulent regime for relatively coarse meshes, with a similar (or even higher) accuracy than a more complex VMS-LES method [START_REF] Gravemeier | Large eddy simulation of turbulent incompressible flows by a three-level finite element method[END_REF]. We studied the practical performances of the solver implemented in a HPC framework, showing strong scalability results up to thousands of cores. This suggests that the present method is efficient, and also provide a reasonable compromise between accuracy and computational complexity, thus proposing it as a suitable and useful tool in the challenging simulation of turbulent flows.

Finally, an extension of the present approach is possible in the following directions:

• problems with outflow boundary conditions. In particular, we aim to compare in the framework proposed in this paper the rotational incremental pressure-correction scheme, cf. [START_REF] Guermond | Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions[END_REF], to an alternative strategy proposed in [START_REF] Poux | Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods[END_REF], which has demonstrated to improve the accuracy for the standard incremental version while remaining compatible with the rotational one; • enhancement of inf-sup stable FE to exactly divergence-preserving schemes, cf. [START_REF] Linke | Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors[END_REF]. This will eventually allow to remove grad-div stabilization; • further improvement of the velocity subgrid model based on the local projection of the streamline derivative in the proposed parallel framework, maybe combined with a VMS-eddy viscosity term and wall laws [START_REF] Chacón Rebollo | Numerical analysis of a finite element projection-based VMS turbulence model with wall laws[END_REF] for a more accurate simulation of turbulent boundary layers in particular, and more complex turbulent flows in general; • coupled flow models like non-isothermal incompressible flows. A step forward this direction has been recently done in [START_REF] Chacón Rebollo | A high-order local projection stabilization method for natural convection problems[END_REF], and the construction of an efficient solver for the method introduced in that paper is today in preparation, following the guidelines given in the present work.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Streamlines on colored velocity distribution at Re = 1,000, 5,000, 10,000 (from top to bottom). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

  ũ2

Fig. 3 .

 3 Fig. 3. Temporal evolution of the total kinetic energy for Re = 3,200, 7,500, 10,000 (from top to bottom).

Fig. 4 .

 4 Fig. 4. u 1 on the vertical centerline (z = 0.5) of the mid-plane y = 0.5 for Re = 3,200, 7,500, 10,000.

Fig. 5 .

 5 Fig. 5. u 3 on the horizontal centerline (x = 0.5) of the mid-plane y = 0.5 for Re = 3,200, 7,500, 10,000.

Fig. 6 .

 6 Fig. 6.

ũ2 1 and ũ2 3

 3 on the centerlines of the mid-plane y = 0.5 for Re = 10,000 (factor 10).

Fig. 7 .

 7 Fig. 7. ũ1 ũ3 on the centerlines of the mid-plane y = 0.5 for Re = 10,000 (factor 500).

Fig. 8 .

 8 Fig.8. TKE on the centerlines of the mid-plane y = 0.5 for Re = 10,000 (factor 100).

Fig. 9 .

 9 Fig. 9. Energy spectrum in logarithmic scale on the centerlines of the mid-plane y = 0.5 for Re = 10,000.

Fig. 10 .

 10 Fig. 10. Flow streamlines at Re = 3,200 (top) and Re = 7,500 (bottom); results for the proposed BDF2 FS-LPS method at final simulation time T = 100.

Fig. 11 .

 11 Fig. 11. Strong scaling analysis of the BDF2 FS-LPS implementation in 3D for a problem of 24 million velocity unknowns and 1 million pressure unknowns.

Fig. 12 .

 12 Fig.[START_REF] Ahmed | Analysis of a full space-time discretization of the Navier-Stokes equations by a local projection stabilization method[END_REF]. Number of iterations for the GMRES method to reach convergence for a 3D problem of 24 million velocity unknowns and 1 million pressure unknowns.

Fig. 13 .

 13 Fig. 13. Comparison of the time spent in the various subroutines involved in the scalability analysis of Fig. 11 for a single time step. (For of the references to color in this figure legend, the reader is referred to the web version of this article.)

www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPtAPNumeric.html.

Acknowledgments

The authors acknowledge Prof. F. Hecht and F. Nataf for the fruitful discussions, suggestions, and advice. S. Rubino gratefully acknowledge the financial support received from IdEx Bordeaux (Excellence Initiative of Université de Bordeaux) and FSMP (Fondation Sciences Mathématiques de Paris) during his postdoctoral research involved in this article. R. Haferssas partially carried out this work while he was funded as a Ph.D. student at Université Pierre et Marie Curie in Laboratoire Jacques-Louis Lions. This work was granted access to the HPC resources of TGCC@CEA under the allocations 2017-A0010607519 made by GENCI.