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

HAL Id: hal-01613834 https://hal.archives-ouvertes.fr/hal-01613834v2 Submitted on 9 Jan 2018 HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers. L’archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d’enseignement et de recherche français ou étrangers, des laboratoires publics ou privés. Efficient and scalable discretization of the Navier–Stokes equations with LPS modeling Ryadh Haferssas, Pierre Jolivet, Samuele Rubino To cite this version: Ryadh Haferssas, Pierre Jolivet, Samuele Rubino. Efficient and scalable discretization of the Navier– Stokes equations with LPS modeling . 2018. <hal-01613834v2> Efficient and scalable discretization of the Navier–Stokes equations with LPS modeling Ryadh Haferssas ∗ , Pierre Jolivet † , Samuele Rubino ‡ January 9, 2018 Abstract 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. Keywords: Navier–Stokes equations, LPS by interpolation, pressure-correction methods, large eddy simulation, turbulent incompressible flows, domain decomposition 2010 Mathematics Subject Classification: 65M55, 65M60, 76D05. 1 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 ∗Department of Aeronautics and Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, United States of America [email protected] †CNRS, IRIT–ENSEEIHT, 2 rue Charles Camichel, 31071 Toulouse Cedex 7, France [email protected] ‡Departamento EDAN & IMUS, Universidad de Sevilla, Avda. Reina Mercedes s/n, 41012 Sevilla, Spain [email protected] 1 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) [65] 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. [55]. 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 [12] being the main exception, until the residual-based Variational MultiScale (VMS) models were introduced in the seminal papers [47,48] and subsequently proposed as implicit LES techniques (ILES) for turbulent flows in [9,25]. 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. [17], 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 [2], 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 [13,54], 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 projection-stabilized 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 2 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. [1,19], 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 [1] 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) [38]. 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 [32], 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 [32,43]. 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++ [45] interfaced with HPDDM. A similar recent study has been performed in [34], 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 [29], 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 [39], 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 [4]. 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 [30] for a monolithic implicit time discretization with a recursive block preconditioning. The proposed fully discrete scheme is tested towards the benchmark problem of recircu- 3 lating 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 = 3,200, 7,500, 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 6states the main conclusions of the paper. 2 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 Rd,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)−→ Rdand p: Ω ×(0, T)−→ Rsuch that: (2.1)        ∂tu+∇·(u⊗u)−2ν∇·D(u) + ∇p=fin Ω ×(0, T), ∇·u= 0 in Ω ×(0, T), u=0on Γ ×(0, T), u(x,0) = u0(x) in Ω, where u⊗uis the tensor function of components uiuj, and D(u) is the symmetric deformation tensor given by D(u) = (1/2)(∇u+∇uT). The unknowns are the velocity u and the pressure pof 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 u0. To define the weak formulation of problem (2.1), we need to introduce some useful notations for spaces. We consider the Sobolev spaces Hs(Ω), s∈R, and Lp(Ω), 1 ⩽p⩽∞. We use the following notation for vector-valued Sobolev spaces: Hsand Lprespectively 4 shall denote [Hs(Ω)]dand [Lp(Ω)]d(similarly for tensor spaces of dimension d×d). Also, the parabolic Bochner function space Lp(0, T;X) (resp. Lp(0, T;X)), where X(resp. X) stands for a scalar (resp. vector-valued) Sobolev space shall be denoted by Lp(X) (resp. Lp(X)). In order to give a variational formulation of problem (2.1), let us consider the velocity space: H1 0= [H1 0(Ω)]d=nw∈[H1(Ω)]d:w=0on Γo. This is a closed linear subspace of H1, and thus a Hilbert space endowed with the H1- norm. Thanks to Korn’s inequality, cf. [46], the H1-norm is equivalent on H1 0to the norm kwkH1 0=kD(w)kL2. Also, let us introduce the space of divergence-free functions: H1 0,div =w∈H1 0:∇·w= 0 a.e. in Ω. The space H1 0,div is a closed linear subspace of H1 0, and thus a Hilbert space endowed with the H1-norm. We shall consider the following variational formulation of (2.1): given f∈L2(H−1)and u0∈H−1, find u∈L∞(L2)∩L2(H1 0,div),P∈L2(L2 0)such that: (2.2)                        −ZT 0 (u(t),v)Ωϕ0(t)dt −hu0,viϕ(0) +ZT 0 [b(u(t),u(t),v) + a(u(t),v)]ϕ(t)dt +ZT 0 (P(t),∇·v)Ωϕ0(t)dt =ZT 0hf(t),viϕ(t)dt, for any v∈H1 0,ϕ∈ D([0, T]) such that ϕ(T) = 0, h·,·i stands for the duality pairing between H1 0and its dual H−1, and L2 0consists of L2-functions with zero mean in Ω. The forms band aare 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∈H1 0. The skew-symmetric form of the convective term bis chosen for conservation purposes: note that b(w,v,v) = 0 for all w,v∈H1 0. The physical pressure is the time derivative of the unknown P:p=∂tP∈H−1(L2 0) = H1 0(0, T;L2 0)0. The interest of considering Pas unknown instead of pis that there are high technical difficulties to obtain uniform bounds for the discrete pressures in a Banach space of space-time functions, see [23, remark 10.2], while one obtains uniform bounds in the Banach space L∞(L2) for the numerical approximation of P, see [1, theorem 4.3]. It is known that for domains which satisfy the cone condition, as bounded polyhedral domains, P∈L∞(L2), e.g., see [35, remark 2.5]. We notice, however, that for practical computations one would approximate the physical pressure p, and Pis introduced just for the numerical analysis. 3 Fully discrete problem In this section, we describe both the spatial approximation and the time discretization proposed for the model problem. 5 3.1 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 [10], and then successfully extended to transport problems [3,8,11,53,58,63]. 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 [56]. 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 penalty-stabilized methods, cf. [16]. 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 [1]. 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. [9,28]. Actually, LPS schemes can be viewed as simplifications of residual-based VMS methods [2], 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 [44,54,71]. In order to describe in detail the spatial approximation of the model problem (2.2), let {Th}h>0be 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∈ Th, its diameter will be denoted by hKand h= maxK∈ThhK. Given a positive integer land a mesh cell K∈ Th, denote by Rl(K) either Pl(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 Ql(K) (i.e., the space of Lagrange polynomials of degree ⩽lon 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: (3.1)            Yl h=Vl h(Ω) = {vh∈C0(Ω) : vh|K∈Rl(K),∀K∈ Th}, Yl h= [Yl h]d={vh∈[C0(Ω)]d:vh|K∈[Rl(K)]d,∀K∈ Th}, Xh=Yl h∩H1 0. Hereafter, Yl h(resp., Yl h) will constitute the discrete foreground vector-valued (resp. scalar) spaces in which we will work on. 6 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 [17]. 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 (uh, ph) : (0, T)→Xh×Mhsuch that: (3.2)        (∂tuh,vh)Ω+b(uh,uh,vh) + a(uh,vh)−(ph,∇·vh)Ω +sconv(uh,uh,vh) + sdiv(uh,vh) = hf,vhi, (∇·uh, qh)Ω+spres(ph, qh)=0, for any (vh, qh)∈Xh×Mh, where Mh=Yl h∩L2 0. The forms sconv,sdiv and spres in (3.2) correspond to a high-order term-by-term stabilized method, cf. [17], and are given by: sconv(uh,wh,vh) = X K∈Th τm,K(σ∗ h(uh·∇wh), σ∗ h(uh·∇vh))K,(3.3) sdiv(uh,vh) = X K∈Th τd,K(σ∗ h(∇·uh), σ∗ h(∇·vh))K,(3.4) spres(ph, qh) = X K∈Th τm,K(σ∗ h(∇ph), σ∗ h(∇qh))K.(3.5) The aims of the stabilization terms is to prevent spurious instabilities due to dominant convection (sconv), achieve additional control on the incompressibility condition by means of pressure subscales effect (sdiv), and ensure inf-sup stability when using equal order interpolation for velocity and pressure (spres). Here, τm,K,τd,K are stabilization coefficients for convection-pressure gradient and divergence, respectively, and σ∗ h=Id −σh, where Id is the identity operator, and σhis some locally stable projection or interpolation operator from L2on the foreground vector-valued space Yl−1 h(also called “buffer space” in this context), satisfying optimal error estimates. In practical implementations, we choose σh as a Scott–Zhang-like [68] linear interpolation operator in the space Yl−1 h, implemented in the software FreeFem++ [45]. This interpolant may be defined as: ∀x∈Ω, σh(v)(x) = X a∈N Πh(v)(a)ϕa(x), where Nis the set of Lagrange interpolation nodes of Yl−1 h,ϕaare the Lagrange basis functions associated to N, and Πhis the interpolation operator by local averaging of Scott– Zhang kind, which coincides with the standard nodal Lagrange interpolant when acting on continuous functions (cf. [17], 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 [6] 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 sconv, the following operators have to be composed: 7 •bilinear form τdiscretized as C: τ:dcYl−1 h×dcYl−1 h→R (uh,vh)7→ X K∈Th τm,K(uh,vh)K, •linear operator σdiscretized as P: σ:dcYl−1 h→dcYl−1 h vh7→ σ∗ h(vh), •linear operator Ddiscretized as D: D:Xh→dcYl−1 h vh7→ uh·∇vh, where dcYl−1 hstands for the discontinuous version of the buffer space Yl−1 h. In practice, we express Das the sum of two (resp. three) operators in 2D (resp. 3D). This is done by considering vhcomponent-by-component, i.e., D=Dx+Dy(resp. D=Dx+Dy+Dz). The discretization of (3.3) is thus the symmetric local matrix Sdefined as: (3.6) S= (P·D)TC(P·D). In formula (3.4), σhdenotes an operator between the scalar spaces L2and Yl−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 spres, so this term is neglected in this case. However, the term sdiv 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 [60]. For this term, just a pure grad-div penalty stabilization is used in practical implementations, i.e., σ∗ h=Id 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 [1]. The working expressions of the stabilization coefficients are: (3.7) τm,K =d c1 ν (hK/l)2+c2 Un K (hK/l)−1 , (3.8) τd,K =(hK/l)2 d c1τm,K , by following the form proposed in [27,28], 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), dis the dimension of the problem, c1and c2are user-chosen positive constants, lis the polynomial degree of the velocity FE approximation, and UKis some local speed on the mesh cell K. The values of the constants c1and c2are chosen to be c1= 4, c2=√c1= 2, cf. [26]. If UK∈L∞(K), the following technical hypothesis on the stabilization coefficients required to perform the numerical analysis is ensured: 8 The mean velocities hu1iand hu3iare computed as a discrete time average according to: huii(x) = 1 N/2 N−1 X n=N/2 ui(x, tn), i = 1,3, N = # time steps = 1,000. Figure 4shows the mean velocity hu1ion the centerline z= 0.5 of the longitudinal midplane y= 0.5, for the various Reynolds numbers under consideration. The proposed method is in agreement with the experimental data of Prasad and Koseff [62], even with the coarse basic discretization at hand, and performs better than the VMS-3L method [40]. A similar accuracy is found for the mean velocity hu3ion the centerline x= 0.5 of the longitudinal mid-plane y= 0.5, see fig. 5. Also, higher-order moments h˜ uni, with n > 1 and ˜ udenoting 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 h˜u2 ii=hu2 ii−huii2, with the standard deviation (root mean square, RMS) defined as qh˜u2 ii(i= 1,3). Also, the offdiagonal component h˜u1˜u3i=hu1u3i−hu1ihu3iof the Reynolds stress tensor is depicted, and finally the turbulent kinetic energy (TKE) defined as 1 2 3 X i=1h˜u2 iiis reported, together with the corresponding energy spectrum. As in Prasad and Koseff [62], 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. 6and 7. Mispredictions of various peaks of these curves may also be found in the numerical results of Zang et al. [72], 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. [40] 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 kis 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 threedimensional 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 = 3,200, in accordance to Prasad and Koseff [62], 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 highfrequency turbulent fluctuations become dominant, and they partially destroy the integrity (or coherence) of the TGL vortices. 5.2.1 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 15 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 [52]. The operations related to the discretization of continuous operators are performed using FreeFem++ [45]. The linear solvers and preconditioners are implemented in HPDDM [51]. Internally, the underlying linear solver for all local subdomain solves is Intel MKL PARDISO [50,66]. 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 Sto 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 [7] and its MatPtAPNumeric1routine (where the symbolic product may be bypassed). The matrix–matrix product P·D, with Dexpressed 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 Avu=fvwith 24.7×106 unknowns, while the single discretization of eq. (3.17) yields a linear system Apq=fpwith 1.1×106unknowns. Both systems are left-preconditioned by either Mv−1 ORAS or Mp−1 GenEO-2 and the GMRES method is stopped when the relative preconditioned residual is lower than 10−8for the velocity unknowns and 10−6for 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 Mp−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 Mv−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 decompositions for subdomain solves, we can see that the fifth column (setup Mv−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 1www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatPtAPNumeric.html 16 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 domain decomposition methods with much less communication. 512 1,024 2,048 4,096 8,192 16,384 5 10 20 50 100 200 # of processes Runtime (seconds) Linear speedup BDF2 FS-LPS Figure 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. 512 1,024 2,048 4,096 8,192 16,384 0 5 10 15 20 # of subdomains # of iterations Velocity systems Pressure systems Figure 12: Number of iterations for the GMRES method to reach convergence for a 3D problem of 24 million velocity unknowns and 1 million pressure unknowns. NStabilization Optimized operators BDF2 & aux. tasks Setup Mv−1 ORAS Solve Avu=fv Solve Apq=fpTotal 512 16.8 s 39.8 s 2.0 s 140.8 s 7.8 s 0.1 s 207.4 s 1,024 7.5 s 22.8 s 1.2 s 40.4 s 6.0 s 0.2 s 78.0 s 2,048 4.6 s 11.6 s 0.7 s 13.3 s 2.7 s 0.2 s 33.1 s 4,096 2.1 s 7.0 s 0.4 s 4.4 s 1.5 s 0.3 s 15.9 s 8,192 1.5 s 4.4 s 0.3 s 2.4 s 0.9 s 0.7 s 10.2 s 16,384 0.9 s 2.8 s 0.2 s 1.1 s 0.4 s 0.8 s 6.1 s Table 1: Breakdown of the time spent in the various subroutines involved in the scalability analysis of fig. 11 for a single time step. Avand Apare of order 24 millions and 1 million, respectively. In fig. 13, we display the relative time (in 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. 17 512 1,024 2,048 4,096 8,192 16,384 0% 20% 40% 60% 80% 100% # of processes Ratio Stabilization Optimized operators Setup Mv−1 ORAS BDF2 & aux. tasks Solve Avu=fv Solve Apq=fp Figure 13: Comparison of the time spent in the various subroutines involved in the scalability analysis of fig. 11 for a single time step. 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., [34]. 6 Summary and conclusions In this work, we have proposed a FE LPS spatial approximation of the NSE, combined with an efficient velocity-pressure 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 [40]. 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. [42], to an alternative strategy proposed in [61], 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. [57]. This will eventually allow to remove grad-div stabilization ; 18 •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 [21] 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 [18], 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. 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´e de Bordeaux) and FSMP (Fondation Sciences Math´ematiques 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. References [1] N. Ahmed, T. Chac´on Rebollo, V. John, and S. Rubino. Analysis of a full spacetime discretization of the Navier–Stokes equations by a local projection stabilization method. IMA J. Numer. Anal., 37(3):1437–1467, 2017. [2] N. Ahmed, T. Chac´on Rebollo, V. John, and S. Rubino. A review of variational multiscale methods for the simulation of turbulent incompressible flows. Arch. Comput. Methods Engrg., 24(1):115–164, 2017. [3] N. Ahmed, G. Matthies, L. Tobiska, and H. Xie. Discontinuous Galerkin time stepping with local projection stabilization for transient convection-diffusion-reaction problems. Comput. Methods Appl. Mech. Engrg., 200(21-22):1747–1756, 2011. [4] D. Arndt, H. Dallmann, and G. Lube. Quasi-optimal error estimates for the incompressible Navier–Stokes problem discretized by finite element methods and pressurecorrection projection with velocity stabilization. NAM Preprint 2015-8, Univeristy of G¨ottingen. Submitted to ESAIM Math. Model. Numer. Anal., 2016. [5] F. Auteri, N. Parolini, and L. Quartapelle. Numerical investigation on the stability of singular driven cavity flow. J. Comput. Phys., 183(1):1–25, 2002. [6] S. Badia. On stabilized finite element methods based on the Scott-Zhang projector. Circumventing the inf-sup condition for the Stokes problem. Comput. Methods Appl. Mech. Engrg., 247/248:65–72, 2012. [7] S. Balay, S. Abhyankar, M. F. Adams, J. Brown, P. Brune, K. Buschelman, L. Dalcin, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, K. Rupp, B. F. Smith, S. Zampini, and H. Zhang. PETSc web page. http://www.mcs.anl. gov/petsc, 2017. 19 [8] G. R. Barrenechea, V. John, and P. Knobloch. A local projection stabilization finite element method with nonlinear crosswind diffusion for convection-diffusion-reaction equations. ESAIM Math. Model. Numer. Anal., 47(5):1335–1366, 2013. [9] Y. Bazilevs, V. M. Calo, J. A. Cottrell, T. J. R. Hughes, A. Reali, and G. Scovazzi. Variational multiscale residual-based turbulence modeling for large eddy simulation of incompressible flows. Comput. Methods Appl. Mech. Engrg., 197(1-4):173–201, 2007. [10] R. Becker and M. Braack. A finite element pressure gradient stabilization for the Stokes equations based on local projections. Calcolo, 38(4):173–199, 2001. [11] R. Becker and M. Braack. A two-level stabilization scheme for the Navier-Stokes equations. In Numerical mathematics and advanced applications, pages 123–130. Springer- Verlag, 2004. [12] J. P. Boris, F. F. Grinstein, E. S. Oran, and R. L. Kolbe. New insights into large eddy simulation. Fluid Dyn. Res., 10(4-6):199–228, 1992. [13] M. Braack and E. Burman. Local projection stabilization for the Oseen problem and its interpretation as a variational multiscale method. SIAM J. Numer. Anal., 43(6):2544–2566, 2006. [14] X.-C. Cai and M. Sarkis. A restricted additive Schwarz preconditioner for general sparse linear systems. SIAM J. Sci. Comput., 21(2):792–797, 1999. [15] F. E. Cellier. Continuous system modeling. Springer-Verlag, New York, 1991. [16] T. Chac´on Rebollo. A term by term stabilization algorithm for finite element solution of incompressible flow problems. Numer. Math., 79(2):283–319, 1998. [17] T. Chac´on Rebollo, M. G´omez M´armol, V. Girault, and I. S´anchez Mu˜noz. A high order term-by-term stabilization solver for incompressible flow problems. IMA J. Numer. Anal., 33(3):974–1007, 2013. [18] T. Chac´on Rebollo, M. G´omez M´armol, F. Hecht, S. Rubino, and I. S´anchez Mu˜noz. A high-order local projection stabilization method for natural convection problems. J. Sci. Comput., Published online:1–26, 2017. [19] T. Chac´on Rebollo, M. G´omez M´armol, and M. Restelli. Numerical analysis of penalty stabilized finite element discretizations of evolution Navier-Stokes equation. J. Sci. Comput., 61(1):1–28, 2014. [20] T. Chac´on Rebollo, M. G´omez M´armol, and S. Rubino. Finite element approximation of an unsteady projection-based VMS turbulence model with wall laws. In P. Knobloch, editor, Boundary and Interior Layers, Computational and Asymptotic Methods - BAIL 2014, volume 108 of Lecture Notes in Computational Science and Engineering, pages 47–73. Springer, 2015. [21] T. Chac´on Rebollo, M. G´omez M´armol, and S. Rubino. Numerical analysis of a finite element projection-based VMS turbulence model with wall laws. Comput. Methods Appl. Mech. Engrg., 285:379–405, 2015. [22] T. Chac´on Rebollo, F. Hecht, M. G´omez M´armol, G. Orzetti, and S. Rubino. Numerical approximation of the Smagorinsky turbulence model applied to the primitive equations of the ocean. Math. Comput. Simulation, 99:54–70, 2014. 20 [23] T. Chac´on Rebollo and R. Lewandowski. Mathematical and numerical foundations of turbulence models and applications. Birkh¨auser, 2014. [24] R. A. Clark, J. H. Ferziger, and W. C. Reynolds. Evaluation of subgrid-scale models using an accurately simulated turbulent flow. J. Fluid Mech., 91:1–16, 1979. [25] R. Codina. Stabilization of incompressibility and convection through orthogonal sub-scales in finite element methods. Comput. Methods Appl. Mech. Engrg., 190(13- 14):1579–1599, 2000. [26] R. Codina. A stabilized finite element method for generalized stationary incompressible flows. Comput. Methods Appl. Mech. Engrg., 190(20-21):2681–2706, 2001. [27] R. Codina and J. Blasco. Analysis of a stabilized finite element approximation of the transient convection-diffusion-reaction equation using orthogonal subscales. Comput. Vis. Sci., 4(3):167–174, 2002. [28] R. Codina, J. Principe, O. Guasch, and S. Badia. Time dependent subscales in the stabilized finite element approximation of incompressible flow problems. Comput. Methods Appl. Mech. Engrg., 196(21-24):2413–2430, 2007. [29] O. Colom´es, S. Badia, R. Codina, and J. Principe. Assessment of variational multiscale models for the large eddy simulation of turbulent incompressible flows. Comput. Methods Appl. Mech. Engrg., 285:32–63, 2015. [30] O. Colom´es, S. Badia, and J. Principe. Mixed finite element methods with convection stabilization for the large eddy simulation of incompressible turbulent flows. Comput. Methods Appl. Mech. Engrg., 304:294–318, 2016. [31] V. Dolean, M. J. Gander, and L. Gerardo-Giorda. Optimized Schwarz methods for Maxwell’s equations. SIAM J. Sci. Comput., 31(3):2193–2213, 2009. [32] V. Dolean, P. Jolivet, and F. Nataf. An introduction to domain decomposition methods: algorithms, theory, and parallel implementation. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, 2015. [33] E. Erturk, T. C. Corke, and C. G¨ok¸c¨ol. Numerical solutions of 2-D steady incompressible driven cavity flow at high Reynolds numbers. Internat. J. Numer. Methods Fluids, 48:747–774, 2005. [34] D. Forti and L. Ded`e. Semi-implicit BDF time discretization of the Navier-Stokes equations with VMS-LES modeling in a high performance computing framework. Comput. & Fluids, 117:168–182, 2015. [35] G. P. Galdi. An introduction to the Navier-Stokes initial-boundary value problem. In Fundamental directions in mathematical fluid mechanics, Adv. Math. Fluid Mech., pages 1–70. Birkh¨auser, Basel, 2000. [36] G. P. Galdi and W. J. Layton. Approximation of the larger eddies in fluid motions. II. A model for space-filtered flow. Math. Models Methods Appl. Sci., 10(3):343–350, 2000. [37] M. Germano, U. Piomelli, P. Moin, and W. H. Cabot. A dynamic subgrid-scale eddy viscosity model. Phys. Fluids, 3:1760–1765, 1991. 21 [38] P. Gervasio, F. Saleri, and A. Veneziani. Algebraic fractional-step schemes with spectral methods for the incompressible Navier-Stokes equations. J. Comput. Phys., 214(1):347–365, 2006. [39] V. Gravemeier, M. W. Gee, M. Kronbichler, and W. A. Wall. An algebraic variational multiscale-multigrid method for large eddy simulation of turbulent flow. Comput. Methods Appl. Mech. Engrg., 199(13-16):853–864, 2010. [40] V. Gravemeier, W. A. Wall, and E. Ramm. Large eddy simulation of turbulent incompressible flows by a three-level finite element method. Internat. J. Numer. Methods Fluids, 48(10):1067–1099, 2005. [41] J.-L. Guermond. Un r´esultat de convergence d’ordre deux en temps pour l’approximation des ´equations de Navier-Stokes par une technique de projection incr´ementale. ESAIM Math. Model. Numer. Anal., 33(1):169–189, 1999. [42] J.-L. Guermond, P. Minev, and J. Shen. Error analysis of pressure-correction schemes for the time-dependent Stokes equations with open boundary conditions. SIAM J. Numer. Anal., 43(1):239–258, 2005. [43] R. Haferssas, P. Jolivet, and F. Nataf. An additive Schwarz method type theory for Lions’s algorithm and a symmetrized optimized restricted additive Schwarz method. SIAM J. Sci. Comput., 39(4):A1345–A1365, 2017. [44] L. He and L. Tobiska. The two-level local projection type stabilization as an enriched one-level approach. Adv. Comput. Math., 36(4):503–523, 2012. [45] F. Hecht. New development in FreeFem++. J. Numer. Math., 20(3-4):251–265, 2012. [46] C. O. Horgan. Korn’s inequalities and their applications in continuum mechanics. SIAM Rev., 37(4):491–511, 1995. [47] T. J. R. Hughes. Multiscale phenomena: Green’s functions, the Dirichlet-to-Neumann formulation, subgrid scale models, bubbles and the origins of stabilized methods. Comput. Methods Appl. Mech. Engrg., 127(1-4):387–401, 1995. [48] T. J. R. Hughes, G. R. Feij´oo, L. Mazzei, and J.-B. Quincy. The variational multiscale method—a paradigm for computational mechanics. Comput. Methods Appl. Mech. Engrg., 166(1-2):3–24, 1998. [49] T. Iliescu, V. John, W. J. Layton, G. Matthies, and L. Tobiska. A numerical study of a class of LES models. Int. J. Comput. Fluid Dyn., 17(1):75–85, 2003. [50] Intel. MKL web page. https://software.intel.com/en-us/intel-mkl, 2017. [51] P. Jolivet, F. Hecht, F. Nataf, and C. Prud’homme. Scalable domain decomposition preconditioners for heterogeneous elliptic problems. In Proceedings of the 2013 ACM/IEEE conference on Supercomputing, SC13, pages 80:1–80:11. ACM, 2013. [52] G. Karypis and V. Kumar. Multilevel k-way partitioning scheme for irregular graphs. Journal of Parallel and Distributed computing, 48(1):96–129, 1998. [53] P. Knobloch. A generalization of the local projection stabilization for convectiondiffusion-reaction equations. SIAM J. Numer. Anal., 48(2):659–680, 2010. 22 [54] P. Knobloch and G. Lube. Local projection stabilization for advection-diffusion- reaction problems: one-level vs. two-level approach. Appl. Numer. Math., 59(12):2891–2907, 2009. [55] A. N. Kolmogorov. The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers. Proc. Roy. Soc. London Ser. A, 434(1890):9–13, 1991. [56] A. Linke. Collision in a cross-shaped domain—a steady 2d Navier-Stokes example demonstrating the importance of mass conservation in CFD. Comput. Methods Appl. Mech. Engrg., 198(41-44):3278–3286, 2009. [57] A. Linke, G. Matthies, and L. Tobiska. Robust arbitrary order mixed finite element methods for the incompressible Stokes equations with pressure independent velocity errors. ESAIM Math. Model. Numer. Anal., 50(1):289–309, 2016. [58] G. Matthies, P. Skrzypacz, and L. Tobiska. Stabilization of local projection type applied to convection-diffusion problems with mixed boundary conditions. Electron. Trans. Numer. Anal., 32:90–105, 2008. [59] P. Moin. Advances in large eddy simulation methodology for complex flows. Internat. J. Numer. Methods Heat Fluid Flow, 23:710–720, 2002. [60] M. Olshanskii, G. Lube, T. Heister, and J. L¨owe. Grad-div stabilization and subgrid pressure models for the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 198(49-52):3975–3988, 2009. [61] A. Poux, S. Glockner, and M. Aza¨ıez. Improvements on open and traction boundary conditions for Navier-Stokes time-splitting methods. J. Comput. Phys., 230(10):4011– 4027, 2011. [62] A. K. Prasad and J. R. Koseff. Reynolds number and end-wall effects on a lid-driven cavity flow. Phys. Fluids, 1:208–218, 1989. [63] H.-G. Roos, M. Stynes, and L. Tobiska. Robust numerical methods for singularly perturbed differential equations. Convection-diffusion-reaction and flow problems, volume 24 of Springer Series in Computational Mathematics. Springer-Verlag, 2nd edition, 2008. [64] S. Rubino. Numerical modeling of turbulence by Richardson number-based and VMS models. PhD thesis, Univeristy of Seville, 2014. [65] P. Sagaut. Large eddy simulation for incompressible flows. Springer-Verlag, Berlin, third edition, 2006. [66] O. Schenk and K. G¨artner. Solving unsymmetric sparse systems of linear equations with PARDISO. Future Generation Computer Systems, 20(3):475–487, 2004. [67] H. Schlichting and K. Gersten. Boundary-layer theory. Springer-Verlag, Berlin, 2000. [68] R. L. Scott and S. Zhang. Finite element interpolation of non-smooth functions satisfying boundary conditions. Math. Comput., 54(190):483–493, 1990. [69] J. Smagorinsky. General circulation experiment with the primitive equations: I. The basic experiment. Mon. Weather Rev., 91(3):99–164, 1963. 23 [70] G. Tiesinga, F. W. Wubs, and A. E. P. Veldman. Bifurcation analysis of incompressible flow in a driven cavity by the Newton-Picard method. In Proceedings of the 9th International Congress on Computational and Applied Mathematics (Leuven, 2000), volume 140, pages 751–772, 2002. [71] L. Tobiska and C. Winkel. The two-level local projection type stabilization as an enriched one-level approach. A one-dimensional study. Int. J. Numer. Anal. Model., 7(3):520–534, 2010. [72] Y. Zang, R. L. Street, and J. R. Koseff. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows. Phys. Fluids, 5:3186–3196, 1993. 24 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 −1.0 −0.5 0.0 z Reynolds stress tensor h˜u1˜u3i Re = 10,000 BDF2 FS-LPS VMS-3L, Gravemeier et al. Experimental data, Prasad and Koseff 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.0 0.5 1.0 x Reynolds stress tensor h˜u1˜u3i Re = 10,000 BDF2 FS-LPS VMS-3L, Gravemeier et al. Experimental data, Prasad and Koseff Figure 7: h˜u1˜u3ion the centerlines of the mid-plane y= 0.5 for Re = 10,000 (factor 500). 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0.0 0.1 0.2 0.3 Distance TKE Re = 10,000 x-direction, BDF2 FS-LPS z-direction, BDF2 FS-LPS Figure 8: TKE on the centerlines of the mid-plane y= 0.5 for Re = 10,000 (factor 100). 31 1 10 10−3 10−2 10−1 100 k E(k) Re = 10,000 x-direction, BDF2 FS-LPS z-direction, BDF2 FS-LPS k−5/3 Figure 9: Energy spectrum in logarithmic scale on the centerlines of the mid-plane y= 0.5 for Re = 10,000. 32 Figure 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. 33