F Boyer
email: fboyer@latp.univ-mrs.fr

F Dardalhon
email: fanny.dardalhon@latp.univ-mrs.fr

C Lapuerta
email: celine.lapuerta@irsn.fr

J.-C Latché

Stability of a Crank-Nicolson Pressure Correction Scheme Based on Staggered Discretizations

Keywords: Pressure correction scheme, compressible Navier-Stokes equations, stable convection operator, low-order finite element approximation

In the context of Large Eddy Simulation of turbulent flows, the control of kinetic energy seems to be an essential requirement for a numerical scheme. Designing such an algorithm, i.e. as less dissipative as possible while being simple, for the resolution of variable density Navier-Stokes equations is the aim of the present work. The developed numerical scheme, based on a pressure correction technique, uses a Crank-Nicolson time discretization and a staggered space discretization relying on the Rannacher-Turek finite element. For the inertia term in the momentum balance equation, we propose a finite volume discretization, for which we derive a discrete analogue of the continuous kinetic energy local conservation identity. Contrary to what was obtained for the backward Euler discretization, the dissipation defect term associated to the Crank-Nicolson scheme is second order in time. This behaviour is evidenced by numerical simulations.

INTRODUCTION

We consider the time-dependent variable density Navier-Stokes equations, posed on a finite time interval (0, T ) and in an open, connected, bounded domain Ω of R d (d = 2 or 3), which is supposed to be polygonal (d = 2) or polyhedral (d = 3): ∂ t ρ + div(ρu) = 0, ∂ t (ρu) + div(u ⊗ ρu)div(τ (u)) + ∇p = 0, [START_REF] Ansanay-Alex | An L2-Stable Approximation of the Navier-Stokes Convection Operator for Low-Order Non-Conforming Finite Elements[END_REF] where u, p and ρ are respectively the velocity, the pressure and the density of the flow. The quantity

τ (u) = µ ∇u + ∇ T u - 2 3 (divu) I
is the shear stress tensor, with µ the positive kinematic viscosity, and u ⊗ ρu is the R d × R d -tensor of components (u ⊗ ρu) i,j = ρ u i u j , ∀ i, j ∈ {1, • • • , d}. As such, System (1) is not closed. The simplest assumption yielding a self-contained problem is just supposing that ρ is a given function of space and time; the mass balance equation (i.e. the first relation of ( 1)) thus has to be considered as a constraint on the velocity, as for incompressible flows (and, indeed, choosing a constant value for ρ yields divu = 0). This equation basically plays the same role in the so-called asymptotic model for low Mach number flows [START_REF] Majda | The Derivation and Numerical Solution of the Equations for Zero Mach Number Combustion[END_REF], where ρ is given as a function of an additional unknown θ (usually, the temperature or a concentration), which satisfies a balance equation:

∂ t (ρ θ) + div(ρ θ u) -div(λ ∇θ) = 0, ρ = ̺(θ). (2)
The function θ → ̺(θ) is the given equation of state for the flow under study, and λ is a non-negative diffusion coefficient.

System (1)-( 2) must be supplemented by suitable boundary conditions, for instance Dirichlet conditions for u and θ, or slip and Neumann conditions: u • n = 0, ((τ (u)p I)n) • t = 0, and ∇θ • n = 0 on (0, T ) × ∂Ω, [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF] where n and t are respectively the outward unit normal and tangential vectors to ∂Ω. Initial conditions are given for u and θ, namely u = u 0 and θ = θ 0 in Ω. We suppose that the equation of state satisfies ̺(θ) > 0 for θ > 0, that θ 0 > 0, and that boundary conditions are such that θ remains positive for all t ∈ (0, T ).

The density ρ and the velocity u are known to satisfy the so-called kinetic energy identity. This relation stems from the following formal computation. Property 1.1 (Property K.E.) Let us assume that ∂ t ρ + div(ρβ) = 0, for a regular field β and that the functions ρ and u are smooth. Then, we have:

∂ t (ρu) + div(u ⊗ ρ β) • u = 1 2 ∂ t ρ|u| 2 + 1 2 div ρ|u| 2 β . ( 4
)
Integrating this relation over Ω yields, assuming that β • n = 0 on ∂Ω:

Ω ∂ t (ρu) + div(u ⊗ ρ β) • u dx = 1 2 d dt Ω ρ |u| 2 dx. (5)
Applying identity [START_REF] Ciarlet | Handbook of Numerical Analysis: Finite Elements Methods -Basic Error Estimates for Elliptic Problems[END_REF], with β = u, to the inner product of the momentum balance equation (second relation of ( 1)) with u yields the so-called (local) kinetic energy identity:

1 2 ∂ t (ρ |u| 2 ) + 1 2 div ρ|u| 2 u) + u • ∇p -div τ (u) • u = 0. ( 6
)
The first two terms of this relation are the time derivative and the transport term of the kinetic energy, respectively; the third one corresponds to the power of the pressure forces, and the last term of the left-hand side reads -div τ (u) • u = -div τ (u) u + τ (u) : ∇u, which makes appear the so-called viscous dissipation τ (u) : ∇u. Integrating over Ω, then integrating by parts the viscous term and using the boundary conditions, we obtain the (global) kinetic energy conservation relation:

1 2 d dt Ω ρ |u| 2 dx + Ω u • ∇p dx + Ω τ (u) : ∇u dx = 0. (7)
This relation does not yield, in its present form, a stability estimate, because of the presence of the pressure term. For incompressible flows, a simple integration by parts shows that it vanishes; for barotropic compressible flows, it may be recast as the time derivative of an energy term. For the low Mach number model, the correct way to treat this term is unclear.

Obtaining a scheme satisfying a discrete analogue of ( 6)-( 7) has many advantages.

-First, possibly combined with the above-mentioned additional arguments to control the pressure term, (7) yields a stability estimate, which is classically observed to considerably enhance the robustness of the scheme, especially for the computation of convection dominant flows.

-Second, in the context of Large Eddy Simulation (LES), a subgrid-scale model is introduced to simulate the (kinetic energy) dissipative role of the small structures in the momentum balance equation. It is thus crucial that the scheme respects this "physical" dissipation, i.e. does not overwhelm the viscous dissipation term in (the discrete counterpart of) Equation ( 6) by too large numerical residual terms [START_REF] Mittal | Suitability of Upwind-Biased Finite Difference Schemes for Large-Eddy Simulation of Turbulent Flows[END_REF][START_REF] Beaudan | Numerical Experiments on the Flow Past a Circular Cylinder at Sub-Critical Reynolds Number[END_REF]. Defining these latter as the "dissipation defect", when designing a numerical scheme for LES, one aim should be to minimize the "dissipation defect".

The Marker And Cell (MAC) space discretization, first introduced in [START_REF] Harlow | Numerical Calculations of Time Dependent Viscous Incompressible Flow of Fluid with a Free Surface[END_REF] and now widely used for the simulation of incompressible flows, applies to Cartesian meshes and is of staggered type, with the pressure unknowns located at the cell centers and the velocity normal components unknowns at the face centers. For divergence-free constant density flows, it has been observed since the middle of the sixties that the natural discrete convection operator associated to this discretization conserves the discrete kinetic energy [START_REF] Lilly | On the Computational Stability of Numerical Solutions of Time-Dependent Non-Linear Geophysical Fluid Dynamics Problems[END_REF], i.e. satisfies a discrete analogue of the integral identity (5) (of course, supposing, in this latter, a constant density ρ). Higher order convection operators, still enjoying the same property, and also a discrete analogue of the local kinetic energy transport equation [START_REF] Desjardins | High Order Conservative Finite Difference Scheme for Variable Density Low Mach Number Turbulent Flows[END_REF], have been proposed in [START_REF] Morinishi | Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow[END_REF][START_REF] Vasilyev | High Order Finite Difference Schemes on Non-Uniform Meshes with Good Conservation Properties[END_REF][START_REF] Verstappen | Symmetry-Preserving Discretization of Turbulent Flow[END_REF][START_REF] Subbareddy | A Fully Discrete, Kinetic Energy Consistent Finite-Volume Scheme for Compressible Flows[END_REF]; combined with a suitable discretization of the pressure gradient term and a suitable time-marching algorithm, these results yield (discrete analogues to) the local and global kinetic energy balances ( 6) and [START_REF] Eymard | Discretization of Heterogeneous and Anisotropic Diffusion Problems on General Nonconforming Meshes SUSHI: a Scheme Using Stabilisation and Hybrid Interfaces 1[END_REF].

Works concerning the MAC scheme for the compressible low Mach number model are more recent and scarcer. First, [START_REF] Nicoud | Conservative High-Order Finite-Difference Schemes for Low-Mach Number Flows[END_REF] generalized the finite difference formula of [START_REF] Morinishi | Fully Conservative Higher Order Finite Difference Schemes for Incompressible Flow[END_REF], to obtain a scheme which appears, from numerical experiments, to satisfy a discrete version of the global identity [START_REF] Eymard | Discretization of Heterogeneous and Anisotropic Diffusion Problems on General Nonconforming Meshes SUSHI: a Scheme Using Stabilisation and Hybrid Interfaces 1[END_REF] up to a second order residual. A similar approach, adapted to cope with cylindrical coordinates, is proposed in [START_REF] Desjardins | High Order Conservative Finite Difference Scheme for Variable Density Low Mach Number Turbulent Flows[END_REF]. Recently, a (structured) MAC scheme is presented in [START_REF] Morinishi | Skew-Symmetric Form of Convective Terms and Fully Conservative Finite Difference Schemes for Variable Density Low-Mach Number Flows[END_REF], followed by a discussion about its conservation properties: both global and local kinetic energy conservation are shown, in both semi-discrete and fully-discrete settings. However, the resulting scheme seems to be rather expensive, and more efficient alternatives with nonzero but still high order residuals are also proposed. For related works in the collocated framework, see [START_REF] Ducros | High-Order Fluxes for Conservative Skew-Symmetric-Like Schemes in Structured Meshes: Application to Compressible Flows[END_REF][START_REF] Felten | Kinetic Energy Conservation Issues Associated with the Collocated Mesh Scheme for Incompressible Flow[END_REF][START_REF] By | Accurate and Stable Finite Volume Operators for Unstructured Flow Solvers[END_REF][START_REF] Georges | Improving Shock-Free Compressible RANS Solvers for LES on Unstructured Meshes[END_REF].

In this article, we pursue the development of kinetic energy conserving schemes for low Mach number flows, with the aim to deal with unstructured staggered discretizations. Precisely speaking, we build a discrete convection operator satisfying the control of the kinetic energy and the reduction of the dissipation defect (discrete analogues of ( 4) and ( 5)) with high-order numerical residuals. The approach we adopt is based on a finite volume structure for the convection operator first introduced in [START_REF] Ansanay-Alex | An L2-Stable Approximation of the Navier-Stokes Convection Operator for Low-Order Non-Conforming Finite Elements[END_REF]; the approximation of the convected velocity at the face is centered, and the time-algorithm is obtained by a Crank-Nicolson technique. Then we implement this discretization in a pressure correction scheme, and make the expression of the dissipation defect explicit. This quantity appears to be formally second order in time, but is unsigned (unlike the first order counterpart which would arise from a backward Euler discretization). It is plotted in numerical experiments to assess the performance of the scheme.

The outline of the article is as follows. After introducing the meshes and the discrete spaces in Section 2, we deal in Section 3 with the construction of the convection operator (Section 3.1) and prove its stability (Property K.E., Section 3.2). The pressure correction scheme for low Mach number flows is presented in Section 4, first in the semi-discrete setting, then in the fully discrete one (Sections 4.1 and 4.2 respectively). The kinetic energy identity satisfied by the fully discrete scheme is established (Section 4.3), then the size of the dissipation defect is evaluated (Section 4.4). Numerical tests are presented in Section 5; they confirm theoretical results and illustrate the ability of the scheme to compute quite efficiently turbulent structures for large time steps.

MESHES AND DISCRETE SPACES

Let M be a decomposition of the domain Ω into quadrangles (d = 2) or hexahedra (d = 3), supposed to be regular in the usual sense of the finite element literature [4, conditions T h 1 to T h 5, p. 61 and 71]. We denote by E the set of all faces σ of the mesh, by E ext the set of faces included in the boundary of Ω, by E int the set of internal faces (i.e. E \ E ext ) and by E(K) the faces of a particular cell K ∈ M. The internal face separating the neighbouring cells K and L is denoted by σ = K|L. Moreover, we also use the notation σ = K|ext if the face σ is included in the boundary of the control volume K and σ belongs to E ext . For each cell K ∈ M and each face σ ∈ E(K), n K,σ stands for the unit normal vector to σ outward K. By |K| and |σ|, we respectively denote the measures of the control volume K and of the face σ.

The velocity u and the pressure p are discretized using the so-called Rannacher-Turek element (also called rotated bilinear finite element [START_REF] Rannacher | Simple Nonconforming Quadrilateral Stokes Element[END_REF]). The reference element K is the unit d-cube (with faces orthogonal to the coordinate axes); the discrete functional space on K is Q 1 ( K), where Q 1 ( K) is defined as follows:

Q 1 ( K) = span 1, (x i ) i=1,...,d , (x 2 i -x 2 i+1 ) i=1,••• ,d-1 .
The mapping T K from the reference element K to the current one K is the standard Q 1 mapping. Besides, we choose the version of the element where the nodal function N K,σ on the cell K is the average of the velocity through the face σ:

∀K ∈ M, ∀σ ∈ E(K), ∀v ∈ C ∞ (K), N K,σ (v) = 1 |σ| σ v |K dσ.
Thus the discrete space X h is defined as:

X h = v ∈ (L 2 (Ω)) d ; ∀K ∈ M, v • T K ∈ Q 1 ( K) d , and, ∀σ = K|L, N K,σ (v) = N L,σ (v) .
The approximation for the velocity is thus non-conforming in H 1 (Ω) d because the space X h is composed of discrete functions which are discontinuous through a face, and thus

X h H 1 (Ω) d .
The set of degrees of freedom for the velocity is:

{u σ,i = N σ (u i ), σ ∈ E, 1 ≤ i ≤ d}, where N σ (u i ) = N K,σ (u i ) = N L,σ (u i ) if σ = K|L. We denote by ϕ (i)
σ the vector shape function associated to u σ,i , which, by definition, reads ϕ

(i) σ = ϕ σ e (i)
, where ϕ σ is the scalar shape function and e (i) is the i th vector of the canonical basis of R d , and we define u σ by u σ = d i=1 u σ,i e (i) . With these definitions, we have the identity:

u(x) = σ∈E d i=1 u σ,i ϕ (i) σ (x) = σ∈E u σ ϕ σ (x), a.e. in Ω.
We respectively denote by ∇v and div v the broken gradient and divergence of any function v ∈ X h . The pressure belongs to the space M h composed of piecewise constant functions:

D σ K L |σ | σ = K |L D σ ′ M σ ′ = K |M ǫ = D σ |D σ ′
M h = p ∈ L 2 (Ω); ∀K ∈ M, p is constant on K .
The degrees of freedom for the pressure are consequently {p K ; K ∈ M}. This pair of finite elements is inf-sup stable [START_REF] Brezzi | Mixed and Hybrid Finite Element Methods[END_REF]. The unknown θ and the density are also approximated by piecewise constant functions, and their degrees of freedom are denoted by {θ K ; K ∈ M} and {ρ K ; K ∈ M}.

For the definition of the schemes, we need a dual mesh which is defined as follows. When K ∈ M is a rectangle or a cuboid, for σ ∈ E(K), we define D K,σ as the cone with basis σ and with vertex the mass center of K. We thus obtain a partition of K in m sub-volumes, where m is the number of faces of the cell (i.e. m = 4 for d = 2 and m = 6 for d = 3), each sub-volume having the same measure |D K,σ | = |K|/m. We extend this definition to general quadrangles and hexahedra, by supposing that we still have built a partition of equal-volume sub-cells, and with the same connectivities; note that this is of course always possible, but that such a volume D K,σ may be no longer a cone, since, if K is far from a parallelogram, it may not be possible to build a cone having σ as basis, the opposite vertex lying in K and a volume equal to |K|/m. The volume D K,σ is referred to as the half-diamond cell associated to K and σ. For σ ∈ E int , σ = K|L, we now define the diamond cell D σ associated to σ by D σ = D K,σ ∪ D L,σ (see Figure 1). For an external face σ ∈ E ext ∩ E(K), D σ is just the same volume as D K,σ . We use the following notation for the set of dual faces of the dual mesh thus built: the set of internal dual faces is referred to as Ēint and the set of faces of a particular dual cell D σ is denoted by Ē(D σ ); the external faces of the dual mesh coincide with primal faces, and the set of dual external faces is thus denoted by E ext . A dual face ǫ separating the diamond cells D σ and D σ ′ is denoted by ǫ = D σ |D σ ′ and a dual face of D σ included in the boundary of the computational domain is denoted by D σ |ext.

CONSTRUCTION OF A STABLE CONVECTION OPERATOR

The aim of this section is to establish an analogue of Property K.E. for the velocity convection operator discretized by the Crank-Nicolson time-algorithm.

Construction of the operator

Let us give the main ideas of this construction, details of which are given in [START_REF] Ansanay-Alex | An L2-Stable Approximation of the Navier-Stokes Convection Operator for Low-Order Non-Conforming Finite Elements[END_REF][START_REF] Gastaldo | Staggered Discretizations, Pressure Correction Schemes and all Speed Barotropic Schemes[END_REF], by considering the problem at a fixed time step. We start from the discrete mass balance equation, supposed to be obtained by a backward Euler finite volume technique based on the primal mesh (see Section 4 for the description of the whole algorithm):

∀K ∈ M, |K| δt (ρ K -ρ * K ) + σ∈E(K) F K,σ = 0. (8)
In this relation, ρ K (resp. ρ * K ) stands for the approximation of the density ρ on the cell K at the end (resp. at the beginning) of the time step, and the quantity F K,σ is the discrete mass flux through the face σ outward K, i.e.

F K,σ = |σ| ρσ u σ • n K,σ ,
ρσ being some approximation of the density at the face σ.

Let D σ be a given dual cell. Our aim is to build approximations of the density ρ σ and ρ * σ and of the mass fluxes through the dual faces of D σ , (F σ,ǫ ) ǫ∈ Ē(Dσ ) , such that the following mass balance equation holds:

|D σ | δt (ρ σ -ρ * σ ) + ǫ∈ Ē(Dσ ) F σ,ǫ = 0. (9)
We first define

ρ σ by    ρ σ |D σ | = |D K,σ | ρ K + |D L,σ | ρ L , if σ = K|L ∈ E int , ρ σ = ρ K , if σ = K|ext ∈ E ext . (10)
Then, we need to build the fluxes F σ,ǫ . The essential idea is to impose the following mass balance over the half dual cell

F K,σ + ǫ∈ Ē(Dσ) ǫ⊂K F σ,ǫ = - |D K,σ | δt (ρ K -ρ * K ). (11)
If we manage to obtain [START_REF] Guermond | A Projection FEM for Variable Density Incompressible Flows[END_REF] for every K and σ, then (9) will hold thanks to the following computations:

-For an external face σ, using that ρ σ = ρ K , ( 11) is exactly [START_REF] Gastaldo | Staggered Discretizations, Pressure Correction Schemes and all Speed Barotropic Schemes[END_REF], provided that we choose the primal mass flux F K,σ as the dual mass flux through σ (remember that σ is both a primal and dual face), which we indeed do.

-For an internal face σ = K|L, combining [START_REF] Guermond | A Projection FEM for Variable Density Incompressible Flows[END_REF] written for K and L with the definition [START_REF] Georges | Improving Shock-Free Compressible RANS Solvers for LES on Unstructured Meshes[END_REF] of the density and using the conservativity of the primal mass fluxes (i.e. F K,σ + F L,σ = 0) yields (9) once again. Let us now use the primal mass balance [START_REF] Felten | Kinetic Energy Conservation Issues Associated with the Collocated Mesh Scheme for Incompressible Flow[END_REF] to rewrite [START_REF] Guermond | A Projection FEM for Variable Density Incompressible Flows[END_REF] as

F K,σ + ǫ∈ Ē(Dσ) ǫ⊂K F σ,ǫ = |D K,σ | |K| σ ′ ∈E(K) F K,σ ′ = 1 m σ ′ ∈E(K) F K,σ ′ . ( 12
)
Writing this relation for all the faces σ of K, we obtain a linear system for the dual fluxes (F σ,ǫ ) ǫ∈ Ē(Dσ) , which is singular (this fact is easily checked by summing these m relations and using the conservativity of the dual mass fluxes). However, thanks to the particular form of the right-hand side of [START_REF] By | Accurate and Stable Finite Volume Operators for Unstructured Flow Solvers[END_REF], this system has a solution and we can select the one which is orthogonal to the kernel of the system. We obtain an expression of the form

F σ,ǫ = σ ′ ∈E(K) α ǫ,σ ′ F K,σ ′ ,
where the coefficients (α ǫ,σ ′ ) are real numbers which appear to be independent of the cell K. For instance, in two dimensions and for a cell K, if we denote by F N , F S , F E and F W the four primal mass fluxes and by F W|S , F S|E , F E|N and F N|W the four dual mass fluxes, as sketched on Figure 2, we obtain expressions of these latter of the form

F σ,ǫ = α W F W + α E F E + α S F S + α N F N , (13)
with the coefficients gathered in Table I.

We are now in position to define a finite volume convection operator C h for any discrete fields u and u * in X h as follows: Table I. Coefficients for the dual mass fluxes in Equation [START_REF] Harlow | Numerical Calculations of Time Dependent Viscous Incompressible Flow of Fluid with a Free Surface[END_REF].

∀σ ∈ E, C h (u, u * ) σ = 1 δt (ρ σ u σ -ρ * σ u * σ ) + 1 |D σ | ǫ∈ Ē(Dσ ) F σ,ǫ u ǫ + u * ǫ 2 , ( 14
)
F S F N F W F E F W | S F S | E F E | N F N | W
F σ,ǫ α W α E α S α N F W|S -3/8 1/8 3/8 -1/8 F S|E -1/8 3/8 -3/8 1/8 F E|N 1/8 -3/8 -1/8 3/8 F N|W 3/8 -1/8 1/8 -3/8
with, for ǫ ∈ Ē, u ǫ (resp. u * ǫ ) given as a function of (u) σ∈E (resp. (u * ) σ∈E ) by:

u ǫ =      1 2 (u σ + u σ ′ ), if ǫ = D σ |D σ ′ ∈ Ēint , u σ , if ǫ = D σ |ext ∈ E ext . ( 15
)
Remark 3.1 With the same arguments, the construction of the convective operator and the dual mass fluxes can be generalized in three dimensions. For instance, 3D computations are performed in [START_REF] Ansanay-Alex | An L2-Stable Approximation of the Navier-Stokes Convection Operator for Low-Order Non-Conforming Finite Elements[END_REF].

Stability analysis

In this section, we prove the discrete counterpart of Property K.E. for the convection operator given by ( 14)-( 15).

Proposition 3.2

Let us assume that the discrete mass balance (9) holds, for any σ ∈ E. Provided that F σ,ǫ = 0 for any ǫ ∈ E ext , we have:

∀σ ∈ E, ∀u, u * ∈ X h , |D σ | C h (u, u * ) σ • u σ + u * σ 2 = 1 2δt |D σ | ρ σ |u σ | 2 -ρ * σ |u * σ | 2 + 1 2 ǫ=Dσ|D σ ′ F σ,ǫ u σ + u * σ 2 • u σ ′ + u * σ ′ 2 + R σ ,
where R σ is a remainder term given by:

R σ = - 1 8δt |D σ | (ρ σ -ρ * σ ) |u σ -u * σ | 2 .
We recognize at the right-hand side a discrete kinetic energy balance with a conservative finite volume discretization of the kinetic energy convection terms.

Proof

We have:

|D σ | C h (u, u * ) σ • u σ + u * σ 2 = T 1 + T 2 , with T 1 = 1 δt |D σ | (ρ σ u σ -ρ * σ u * σ ) • u σ + u * σ 2 , T 2 = u * σ + u σ 2 • ǫ∈ Ē(Dσ ) F σ,ǫ u ǫ + u * ǫ 2 .
The term T 1 reads:

T 1 = 1 2δt |D σ | (ρ σ |u σ | 2 -ρ * σ |u * σ | 2 ) T1,1 + 1 2δt |D σ | (ρ σ -ρ * σ ) u σ • u * σ T1,2 . ( 16
)
We now turn to T 2 . Let us introduce the following notation:

ūσ = u σ + u * σ 2 , ∀σ ∈ E and ūǫ = u ǫ + u * ǫ 2 , ∀ǫ ∈ Ē. (17)
We get:

T 2 = ūσ • ǫ∈ Ē(Dσ)
F σ,ǫ ūǫ .

Using ūσ • ūǫ = |ū σ | 2 + ūσ • (ū ǫ -ūσ ),
we obtain:

T 2 = |ū σ | 2 ǫ∈ Ē(Dσ) F σ,ǫ T2,1 + ǫ∈ Ē(Dσ ) F σ,ǫ ūσ • (ū ǫ -ūσ ) T2,2
.

Let us consider the second term T 2,2 . Using the identity 2a

• (a -b) = |a| 2 + |a -b| 2 -|b| 2 , valid
for any real vectors a and b, we get:

T 2,2 = - 1 2 |ū σ | 2 ǫ∈ Ē(Dσ ) F σ,ǫ + 1 2 ǫ∈ Ē(Dσ ) F σ,ǫ (|ū ǫ | 2 -|ū σ -ūǫ | 2 ).
Consequently,

T 2 = 1 2 |ū σ | 2 ǫ∈ Ē(Dσ) F σ,ǫ + 1 2 ǫ∈ Ē(Dσ) F σ,ǫ (|ū ǫ | 2 -|ū σ -ūǫ | 2 ).
Using the discrete mass balance equation ( 9), we get:

T 2 = - 1 2 δt |D σ | (ρ σ -ρ * σ ) |ū σ | 2 T2,2,1 + 1 2 ǫ∈ Ē(Dσ) F σ,ǫ (|ū ǫ | 2 -|ū σ -ūǫ | 2 ) T2,2,2
.

First, the term T 2,2,1 may be combined with the term T 1,2 in ( 16), and using the definition (17) of ūσ , this gives:

T 2,2,1 + T 1,2 = - 1 8δt |D σ | (ρ σ -ρ * σ ) |u σ -u * σ | 2 = R σ .
Second, by the centered definition [START_REF]A CFD computer code for the simulation of reactive turbulent flows[END_REF] of the approximation of the convected velocity at the face, T 2,2,2 reads also:

T 2,2,2 = 1 2 ǫ=Dσ|D σ ′ F σ,ǫ ūσ • ūσ ′ + 1 2 ǫ=Dσ |ext F σ,ǫ |ū σ | 2 .
Finally, with Equation ( 17), we obtain:

T 2,2,2 = 1 2 ǫ=Dσ |D σ ′ F σ,ǫ (u σ + u * σ ) 2 • (u σ ′ + u * σ ′ ) 2 + 1 8 |u σ + u * σ | 2 ǫ=Dσ|ext F σ,ǫ ,
and we conclude the proof invoking the assumption F σ,ǫ = 0 for any ǫ = D σ |ext.

A CRANK-NICOLSON-LIKE PRESSURE CORRECTION SCHEME

We now present a pressure correction scheme for the resolution of the low Mach number Navier-Stokes system (1). We begin by presenting the algorithm in the time semi-discrete setting (Section 4.1), then give the space discretization (Section 4.2). Using the property of the convection operator proved in Section 3, we then establish the kinetic energy identity satisfied by the scheme (Section 4.3). In this relation, a remainder term appears when the density is not constant, the time order of which is discussed (Section 4.4).

The time-marching algorithm

We consider a uniform subdivision of the interval (0, T ) denoted by

0 = t 0 < t 1 < • • • < t N = T . Let δt be the constant time step, namely δt = t n+1 -t n for any integer n in {0, • • • , N -1}.
Let n in {0, • • • , N -1} be fixed and let p n , θ n and ρ n-1 , ρ n , u n be given such that

ρ n -ρ n-1 δt + div(ρ n u n ) = 0.
The Crank-Nicolson-like time splitting algorithm consists in four steps performed successively:

1-Balance equation for θ -Solve for θ n+1 :

1 δt (ρ n θ n+1 -ρ n-1 θ n ) + div(ρ n θ n+1 u n ) -div(λ ∇θ n+1 ) = 0 in Ω, ∇θ n+1 • n = 0 on ∂Ω.
2-Update of the density -Update the density with the equation of state:

ρ n+1 = ̺(θ n+1 ).
3-Velocity prediction -Solve for u n+1 :

1 δt (ρ n u n+1 -ρ n-1 u n ) + div( u n+1/2 ⊗ ρ n u n ) -div τ ( u n+1/2 ) +∇p n = 0 in Ω, u n+1/2 • n = 0 and τ ( u n+1/2 ) • n • t = 0 on ∂Ω, (18)
with

u n+1/2 = ( u n+1 + u n )/2.
4-Velocity and pressure correction -Solve for u n+1 and p n+1 :

ρ n δt (u n+1 -u n+1 ) + 1 2 ∇(p n+1 -p n ) = 0 in Ω, 1 δt (ρ n+1 -ρ n ) + div(ρ n+1 u n+1 ) = 0 in Ω, u n+1 • n = 0 on ∂Ω. (19)
The first step of the algorithm corresponds to the resolution of the conservation equation for θ and yields θ n+1 ; then, thanks to the equation of state, the second step gives the density ρ n+1 at time t n+1 = (n + 1) δt.

The aim of the third step, called velocity prediction step, is to compute an intermediate velocity field denoted by u n+1 , using the pressure p n and the density fields ρ n and ρ n-1 already evaluated at the previous time steps. Since we have in view to use the theory developed in Section 3, we need the densities and convection field involved in the convection operator to satisfy a (semi-)discrete mass balance; this latter equation being solved only in Step 4 of the algorithm, we have to rely for this on the mass balance at the previous time step, and we thus are lead to perform a backward time-shift of the densities in the time derivative term.

Finally, the fourth step, referred to as the velocity and pressure correction step, and formulated as a Darcy problem, may seem numerically expensive. In fact, by taking the divergence of the first equation multiplied by ρ n+1 /ρ n and using the second one, we are able to eliminate the end-of-step velocity and reformulate this step as an elliptic problem for the pressure increment as follows:

- δt 2 div ρ n+1 ρ n ∇(p n+1 -p n ) = - 1 δt (ρ n+1 -ρ n ) -div(ρ n+1 u n+1 ) in Ω, ∇(p n+1 -p n ) • n = 0 on ∂Ω. ( 20
)
The underlying operator is the Schur complement of the Darcy problem. This computation is actually done at the algebraic level (i.e. equations are first discretized in space before making this computation), which necessitates a priori to compute the inverse of the velocity mass matrix. However, since we adopt here a finite volume discretization of the convection operator (see Section 3), the corresponding velocity mass matrix (which is not the usual Rannacher-Turek mass matrix) is diagonal, see [START_REF] Morinishi | Skew-Symmetric Form of Convective Terms and Fully Conservative Finite Difference Schemes for Variable Density Low-Mach Number Flows[END_REF]. Once the pressure increment is known, we deduce the values of the pressure and the end-of-step velocity.

The time accuracy of this algorithm is not analyzed in details here. However, we can see that the system (20) can be understood as a consistent discretization of the following approximate mass balance equation:

∂ρ ∂t + div(ρu) - δt 2 2 ∆ ∂p ∂t = 0.
Following the approach proposed in [START_REF] Rannacher | On Chorin's Projection Method for the Incompressible Navier-Stokes Equations[END_REF] for studying Chorin's projection methods, this remark can certainly lead to a suitable convergence analysis that we may develop in future works. Moreover, we give in Section 5 some numerical evidences that, at least for analytic smooth solutions of the system, our scheme is actually convergent.

Remark 4.1 (Link with a Crank-Nicolson time approximation) By summing up the velocity prediction [START_REF] Lilly | On the Computational Stability of Numerical Solutions of Time-Dependent Non-Linear Geophysical Fluid Dynamics Problems[END_REF] and correction ( 19) equations, we obtain a discrete momentum balance equation with a Crank-Nicolson-like discretization (i.e. with a middle-of-step approximation of the unknowns in the convection, diffusion and pressure gradient terms):

1 δt (ρ n u n+1 -ρ n-1 u n ) + div( u n+1/2 ⊗ ρ n u n ) -div τ ( u n+1/2 ) + ∇p n+1/2 = 0, with p n+1/2 = (p n+1 + p n )/2.

Remark 4.2 (Splitting error in time)

This algorithm is formally first order in time when the density varies with time, thus, even if we have derived the prediction and correction steps by using a Crank-Nicolson-like technique, we cannot expect a second order convergence in this case. This would require a more expensive second order splitting technique. Note that our aim here is only to limit the dissipation defect of the scheme.

Remark 4.3 (Convection operator for θ)

At the continuous level, θ satisfies a maximum principle. The essential argument which yields this conclusion is that the convection operator in the conservation equation for θ may be recast, using the mass balance equation, as a transport operator:

∂ t (ρθ) + div(ρθu) = ρ ∂ t θ + u • ∇θ .
To obtain this property at the discrete level, we thus have to build a convection operator consistent with the mass balance [START_REF] Larrouturou | How to Preserve the Mass Fractions Positivity when Computing Compressible Multi-Component Flows[END_REF] (in other words, which is seen to vanish, invoking the mass balance equation, for constant θ fields). To this purpose, once again, we have to shift in time the density fields in the time derivative term. Note however that, at the fully discrete level, we fall short to obtain a discrete maximum principle in the general case, since the Laplace operator we use enjoys this property only for specific meshes (see Remark 4.5 below).

The fully discrete scheme

We now give the fully discrete form of the scheme. The construction of the velocity convection operator in the discrete momentum balance equation and the form of the discrete mass balance equation are conform to the developments of Section 3.

Let us assume that (ρ n-1 K ) K∈M , (ρ n K ) K∈M , (θ n K ) K∈M , (u n σ ) σ∈E and (p n K ) K∈M are known, and let us suppose that:

∀K ∈ M, |K| δt (ρ n K -ρ n-1 K ) + σ∈E(K) F n K,σ = 0.
This relation is the mass balance equation at the previous time step, and the expression of the mass fluxes (F K,σ ) σ∈E(K) as a function of the density and the velocity is given below (Step 4 of the algorithm).

Then we perform successively:

1-Balance equation for θ -Solve for (θ n+1 K ) K∈M : ∀K ∈ M, |K| δt (ρ n K θ n+1 K -ρ n-1 K θ n K ) + σ∈E(K) F n K,σ θ n+1 σ -div(λ∇θ n+1 ) K = 0,
where θ σ is an upwind approximation for the unknown θ at the face σ and div(λ∇•) K is a finite volume approximation of the (variable coefficient) Laplace operator: for admissible meshes (i.e. here, for structured grids), we use the usual two-point flux approximation, while the variant of the SUSHI method using only cell centered unknowns [START_REF] Eymard | Discretization of Heterogeneous and Anisotropic Diffusion Problems on General Nonconforming Meshes SUSHI: a Scheme Using Stabilisation and Hybrid Interfaces 1[END_REF][START_REF] Piar | A Formally Second-Order Cell Centered Scheme for Convection-Diffusion Equations on General Grids[END_REF] is implemented when the mesh is unstructured.

2-Update of the density -Compute the end-of-step density with the equation of state:

∀K ∈ M, ρ n+1 K = ̺(θ n+1 K ).
3-Velocity prediction -Disregarding for a while the boundary conditions, this step consists in solving for the predicted velocity unknowns ( u n+1 σ ) σ∈E , the following system:

∀σ ∈ E, for 1 ≤ i ≤ d, e (i) • |D σ | δt (ρ n σ u n+1 σ -ρ n-1 σ u n σ ) + ǫ∈ Ē(Dσ ) F n σ,ǫ u n+1/2 ǫ + |D σ | (∇p n ) σ + K∈M K τ ( u n+1/2 ) : ∇ϕ (i) σ dx = 0, (21)
where we recall that u n+1/2 = ( u n+1 + u n )/2. The convection operator (i.e. the computation of the density at the faces ρ σ , of the mass fluxes at the dual faces (F σ,ǫ ) and the (centered) approximation u n+1/2 ǫ of the velocity at the dual faces) is described in Section 3.1. The pressure gradient term reads:

|D σ | (∇p n ) σ • e (i) = - K∈M K p n divϕ (i) σ dx,
and thus, for internal faces,

|D σ | (∇p n ) σ = |σ| (p n L -p n K ) n K,σ , ∀σ ∈ E int , σ = K|L.
On an external face, the impermeability boundary condition must be taken into account by (possibly, i.e. if the normal of the face is not collinear to a coordinate axis) making a change of unknowns for the velocity to make appear the component(s) of the velocity tangent and normal to the boundary, and prescribing this latter to zero; at each σ ∈ E ext , an equation is thus suppressed from the system. Since the pressure gradient at a face σ is normal to σ, it is thus useless to define it for σ ∈ E ext . Note that, since the Neumann boundary condition associated to the tangential component(s) is homogeneous, no extra surface term appears from the integration by parts of the viscous term.

4-Velocity and pressure correction step -Solve for the end-of-step velocity unknowns (u n+1 σ

) σ∈E and the pressure unknowns (p n+1 K ) K∈M the following system:

∀σ ∈ E int , |D σ | δt ρ n σ (u n+1 σ -u n+1 σ ) + 1 2 |D σ | ∇(p n+1 -p n ) σ = 0, ∀K ∈ M, |K| δt (ρ n+1 K -ρ n K ) + σ∈E(K) F n+1 K,σ = 0. (22)
Note that, since the pressure gradient is normal to a face and the velocity normal to the boundary is prescribed to zero, the velocity on the boundary is left unchanged, at this step. This explains why the system associated to the velocity correction (first equation of ( 22)) may be restricted to internal faces.

Since the density is given by the equation of state, there is no need to use upwinding in the mass balance equation (by opposition with compressible flows problem [START_REF] Herbin | Staggered Schemes for All Speed Flows[END_REF], for instance, where upwinding ensures the positivity of ρ), and we thus use a centered approach:

∀σ ∈ E int , σ = K|L, F n+1 K,σ = |σ| (ρ n+1 K + ρ n+1 L ) 2 u n+1 σ • n K,σ , ∀σ ∈ E ext , σ = K|ext, F n+1 K,σ = 0.
Remark 4.4 (Spatial order of accuracy) Despite the fact that we use a finite volume approximation of the convective operator instead of its finite element form, the overall accuracy of the method is preserved, see Section 5.1.

Remark 4.5 (Maximum principle for θ and positivity of the density) For structured meshes, the upwind choice for θ in the discrete convection operator of the associated balance conservation equation, together with the particular form of the Laplace operator, ensure a discrete maximum principle, and thus the positivity of the density. When the mesh is unstructured, no such property is known for the SUSHI scheme. However, we did not observe under-or overshoots in practical computations.

The discrete kinetic energy identity

We first define the following density-dependent discrete semi-norm for the pressure p in M h :

|p| 2 1,M,ρ = σ∈Eint 1 ρ σ |D σ | |(∇p) σ | 2 = σ∈Eint, σ=K|L |σ| 2 ρ σ |D σ | (p K -p L ) 2 . ( 23
)
This semi-norm is well defined as soon as the density is positive, which is always the case in practice (see Remark 4.5).

The application of the discrete counterpart of the continuous procedure described in introduction yields the following result.

Proposition 4.6 (Global kinetic energy balance)

Let us assume that F σ,ǫ = 0 for any ǫ ∈ E ext . Provided that the discrete mass balance ( 9) holds, we have, for 0 ≤ n ≤ N -1:

1 2 σ∈E |D σ | (ρ n σ |u n+1 σ | 2 -ρ n-1 σ |u n σ | 2 ) + δt K∈M K τ ( u n+1/2 ) : ∇ u n+1/2 dx - δt 2 K∈M p n K σ∈E(K) |σ| u n σ • n K,σ + p n+1 K σ∈E(K) |σ| u n+1 σ • n K,σ + D n+1 = 0, ( 24
)
where D n+1 stands for a numerical remainder term which reads

D n+1 = δt 2 8 ( |p n+1 | 2 1,M,ρ n -|p n | 2 1,M,ρ n ) - 1 8 σ∈E |D σ | (ρ n σ -ρ n-1 σ ) | u n+1 σ -u n σ | 2 . ( 25
)

Proof

Let n ∈ {0, . . . , N -1} and σ ∈ E. Let us multiply the discrete velocity prediction step ( 21) associated to the face σ and the component i by ( u n+1/2 σ ) i for 1 ≤ i ≤ d and then sum over i. Invoking Proposition 3.2, we get

|D σ | 2δt ρ n σ | u n+1 σ | 2 -ρ n-1 σ |u n σ | 2 + 1 2 ǫ=Dσ |D σ ′ F n σ,ǫ u n+1/2 σ • u n+1/2 σ ′ + |D σ | (∇p n ) σ • u n+1/2 σ + |D σ | -divτ ( u n+1/2 ) σ • u n+1/2 σ + R n+1 σ = 0, ( 26
)
where the discrete approximation of the viscous diffusion operator applied to u n+1/2 reads

for 1 ≤ i ≤ d, for σ ∈ E, |D σ | -divτ ( u n+1/2 ) σ,i = K∈M K τ ( u n+1/2 ) : ∇ϕ (i) σ dx,
and the remainder term

R n+1 σ reads R n+1 σ = - |D σ | 8 δt (ρ n σ -ρ n-1 σ ) | u n+1 σ -u n σ | 2 .
We now recast the first equation of the correction step [START_REF] Morinishi | Skew-Symmetric Form of Convective Terms and Fully Conservative Finite Difference Schemes for Variable Density Low-Mach Number Flows[END_REF] as follows:

|D σ | ρ n σ 2δt 1/2 u n+1 σ + |D σ | δt 8ρ n σ 1/2 ∇p n+1 σ = |D σ | ρ n σ 2δt 1/2 u n+1 σ + |D σ | δt 8ρ n σ 1/2 ∇p n σ .
Taking the square of the norm of both sides of this relation, we get

|D σ | 2δt ρ n σ |u n+1 σ | 2 + 1 2 |D σ | ∇p n+1 σ • u n+1 σ + δt 8 |D σ | ρ n σ |(∇p n+1 ) σ | 2 = |D σ | 2δt ρ n σ | u n+1 σ | 2 + 1 2 |D σ | ∇p n σ • u n+1 σ + δt 8 |D σ | ρ n σ |(∇p n ) σ | 2 .
Summing with [START_REF] Piar | A Formally Second-Order Cell Centered Scheme for Convection-Diffusion Equations on General Grids[END_REF] and multiplying by δt, we obtain

|D σ | 2 ρ n σ |u n+1 σ | 2 -ρ n-1 σ |u n σ | 2 + δt 2 ǫ=Dσ |D σ ′ F n σ,ǫ u n+1/2 σ • u n+1/2 σ ′ + δt |D σ | 2 (∇p n ) σ • u n σ + (∇p n+1 ) σ • u n+1 σ + |D σ | δt -divτ ( u n+1/2 ) σ • u n+1/2 σ + D n+1 σ = 0, ( 27
)
with

D n+1 σ = - |D σ | 8 (ρ n σ -ρ n-1 σ ) | u n+1 σ -u n σ | 2 + δt 2 8 |D σ | ρ n σ |(∇p n+1 ) σ | 2 - |D σ | ρ n σ |(∇p n ) σ | 2 .
We then sum over all the faces. The kinetic energy convection term (i.e. the second term at the left-hand side) vanishes by conservativity (for

ǫ = D σ |D σ ′ , F n σ,ǫ = -F n σ ′ ,ǫ
) and by the boundary conditions (for ǫ = D σ |ext, F n σ,ǫ = 0). We make a discrete integration by parts in the gradient terms, i.e. we use the fact that, thanks to the definition of the pressure gradient and, once again, to the boundary conditions, for any discrete functions u and p, σ∈Eint, σ=K|L

|D σ | (∇p) σ • u σ = σ∈Eint, σ=K|L |σ| (p L -p K ) n K,σ • u σ = - K∈M p K σ∈E(K) |σ| u σ • n K,σ .
Finally, we easily get that, by definition,

σ∈E |D σ | -divτ ( u n+1/2 σ • u n+1/2 σ = K∈M K τ ( u n+1/2 ) : ∇ u n+1/2 dx,
which concludes the proof.

Remark 4.7 (Local kinetic energy identity)

Note that Equation ( 27) is the discrete analogue of the (local) kinetic energy identity ( 6), to which it may be identified term by term. Both the kinetic energy flux (second term) and the pressure gradient term (third one) are consistent with centered in time and space discretizations.

Such a relation may be used for different purposes. For instance, one may find in [START_REF] Herbin | Staggered Schemes for All Speed Flows[END_REF] a way to exploit such an identity to switch in Euler equations from the total energy to the internal energy balance without loosing the consistency of the scheme.

Order in time of the dissipation defect

In the previous subsection, it has been shown (Proposition 4.6) that we are able to write a discrete equation governing the variation of kinetic energy for the Crank-Nicolson-like pressure correction scheme presented before. Adding up this relation (i.e. Equation ( 24)) for n ranging from 0 to N -1, we get

1 2 σ∈E |D σ | (ρ N -1 σ |u N σ | 2 -ρ -1 σ |u 0 σ | 2 ) + δt N -1 n=0 K∈M K τ ( u n+1/2 ) : ∇ u n+1/2 dx - δt 2 N -1 n=0 K∈M σ∈E(K) |σ| (p n K u n σ + p n+1 K u n+1 σ ) • n K,σ + D = 0, ( 28
)
with

D = δt 2 8 N -1 n=0 |p n+1 | 2 1,M,ρ n -|p n | 2 1,M,ρ n - 1 8 N -1 n=0 σ∈E |D σ | (ρ n σ -ρ n-1 σ ) | u n+1 σ -u n σ | 2 .
Relation [START_REF] Rannacher | On Chorin's Projection Method for the Incompressible Navier-Stokes Equations[END_REF] does not provide a stability estimate, for three reasons 1 -The third term at the left-hand side is the discrete counterpart of the continuous term

T 0 Ω
p divu dx dt, the control of which is not possible in the low Mach number model;

for an incompressible flow, this term would simply vanish.

-

The first term in the remainder D is not non-negative, because of the presence of different density weights in pressure semi-norms; this problem would be cured if the quantity |p n | 1,M,ρ n was replaced by |p n | 1,M,ρ n-1 , and this can be done by adding to the algorithm a pressure renormalization step, following the study presented in [START_REF] Guermond | A Projection FEM for Variable Density Incompressible Flows[END_REF]. This step is not implemented here, because it does not seem to be crucial for the robustness of our scheme.

-

The last term in D can take any sign. Note that none of these issues occurs when the density is constant (divu = 0 and D = 0), so that the scheme is unconditionally stable in this case.

The quantity D is a numerical remainder term, which is referred to hereafter as the dissipation defect of the scheme. We now show that D is formally second order in time. Indeed, we have, using Equation ( 23) and reordering the sum in time

δt 2 8 N -1 n=0 |p n+1 | 2 1,M,ρ n -|p n | 2 1,M,ρ n = δt 2 8 |p N | 2 1,M,ρ N -1 -|p 0 | 2 1,M,ρ 0 + δt 2 8 N -1 n=1 σ∈Eint |D σ | |(∇p n ) σ | 2 1 ρ n-1 σ - 1 ρ n σ .
Therefore, each sum in D is a product of δt 2 multiplied by a sum that behaves like O(1) in time.

Remark 4.8 (First order backward Euler time discretization) When approximating the velocity time derivative with first order backward Euler scheme [START_REF] Ansanay-Alex | An L2-Stable Approximation of the Navier-Stokes Convection Operator for Low-Order Non-Conforming Finite Elements[END_REF], a similar analysis yields, for

0 ≤ n ≤ N -1 1 2 σ∈E |D σ | (ρ n σ |u n+1 σ | 2 -ρ n-1 σ |u n σ | 2 ) + δt K∈M K τ ( u n+1 ) : ∇ u n+1 dx -δt K∈M p n+1 K σ∈E(K) |σ| u n+1 σ • n K,σ + D n+1 Euler = 0, (29)
where the numerical remainder term D n+1 Euler now reads

D n+1 Euler = δt 2 2 ( |p n+1 | 2 1,M,ρ n -|p n | 2 1,M,ρ n ) + 1 2 σ∈E |D σ | ρ n σ | u n+1 σ -u n σ | 2 . ( 30
)
Summing up Equalities [START_REF] Roig | Zone de Mélange d'Ecoulements Diphasiques à Bulles[END_REF] for n ranging from 0 to N -1, we get the following expression for the dissipation defect of the scheme

D Euler = δt 2 2 N -1 n=0 |p n+1 | 2 1,M,ρ n -|p n | 2 1,M,ρ n + 1 2 N -1 n=0 σ∈E |D σ | ρ n σ | u n+1 σ -u n σ | 2 ,
and, if, as in the Crank-Nicolson-like pressure correction scheme, the first sum behaves like O(δt 2 ), the second one behaves as O(δt) and the dissipation defect is thus first order in time.

NUMERICAL EXPERIMENTS

In this section, we report some numerical experiments performed to assess the behaviour of the Crank-Nicolson-like pressure correction scheme proposed in this paper. To this purpose, we compare it to a classical pressure correction scheme [START_REF] Ansanay-Alex | An L2-Stable Approximation of the Navier-Stokes Convection Operator for Low-Order Non-Conforming Finite Elements[END_REF], based on the same space discretization and on the first order backward Euler approximation in time. The simulations are performed with the ISIS software [START_REF]A CFD computer code for the simulation of reactive turbulent flows[END_REF] based on the development platform and software component library PELICANS [START_REF] Pelicans | Collaborative development environment[END_REF], both developed at IRSN (Institut de Radioprotection et de Sûreté Nucléaire).

Convergence and dissipation defect

We address here a test case with an analytical solution, in order to assess the convergence properties of the scheme under study. We also check that the behaviour of the dissipation defect with respect to the time step is consistent with the theoretical results of Section 4.4, for both Euler and Crank-Nicolson time discretizations. θ, in the following discrete L 2 norm:

||θ|| L 2 = K∈M |K| θ 2 K 1 2
, which is usual in the finite volume framework. For each unknown, we observe a first order in time convergence, until reaching a plateau, the value of which corresponds to the space accuracy of the method. These figures confirm in particular the convergence properties of the overall algorithm even though we do not provide any convergence result in this paper.

We now investigate the behaviour of the dissipation defect for both the Euler and the Crank-Nicolson schemes. For this study, the mesh is fixed and chosen to be the 128 × 128 regular grid. In the definition of the dissipation defect given by [START_REF] Pelicans | Collaborative development environment[END_REF] for Crank-Nicolson and by [START_REF] Schäfer | Benchmark Computations of Laminar Flow Around a Cylinder[END_REF] for Euler, we distinguish the pressure gradient part, noted D n P , from the kinetic energy part, noted D n E , and we define D P and D E as the the maximum values in time of each part:

D P = max 0≤n≤N -1 n k=0 D k+1 P , D E = max 0≤n≤N -1 n k=0 D k+1 E .
We draw on Figure 6 the obtained results as a function of the time step for both algorithms. As predicted by the analysis of Section 4.4, the pressure gradient part is of order two for both schemes and the kinetic energy part is of order two for the Crank-Nicolson scheme while the order is one for the Euler discretization. Moreover, it is also observed that the Crank-Nicolson scheme leads to a much smaller kinetic energy part (note that, for this scheme, the dissipation defect is multiplied by 10 2 in the figure).

2D laminar flow around a cylinder

We address in this section some test cases which are derived from a benchmark referred to as the 2D-2 case in [START_REF] Schäfer | Benchmark Computations of Laminar Flow Around a Cylinder[END_REF]. The geometry is sketched on Figure 7. A disk of diameter D = 0.1 is located inside a 2D channel (the axis of the cylinder lying slightly below the centerline of the channel) and the fluid enters the domain on the left boundary, with an imposed parabolic velocity profile:

u x (0, y) = 4u m y H -y H 2 , u y (0, y) = 0,
where H = 0.41m is the height of the channel and u m = 1.5. The velocity is prescribed to zero at the other boundaries except for the right-hand side, where the flow leaves the domain; in our computation, we impose on this latter boundary an homogeneous Neumann boundary condition.

Our study is split into two parts. In the first step, the flow is supposed to be isotherm; the density ρ is then constant, the flow is incompressible, and we can compare our results to the reference ones provided in [START_REF] Schäfer | Benchmark Computations of Laminar Flow Around a Cylinder[END_REF]. In the second case, density variations are created in the flow by prescribing at the cylinder boundary a temperature greater than in the bulk.

The coarse mesh (mesh#1) used for the presented computation is given on Figure 8; other meshes are refined with respect to this one, by setting to smaller and smaller values the discretization step along some characteristic lines (the boundaries and concentric circles around the cylinder). The numbers of cells and unknowns for each mesh are reported in Table II.

The flow is unsteady (see Figure 9 for a visualization at a given time). The characteristic flow quantities quoted in [START_REF] Schäfer | Benchmark Computations of Laminar Flow Around a Cylinder[END_REF] are the maximum drag coefficient c D max , the maximum lift coefficient c Lmax , the Strouhal number St and an instantaneous pressure difference ∆P between the front and end points of the cylinder, i.e. the points (0.15, 0.20) and (0.25, 0.20) (see [30, section 2.2]). Drag and lift coefficients are defined as follows. First, the drag and lift forces are expressed as: where S stands for the boundary of the cylinder and u t for the tangential component of the velocity (namely u t = u • t with t the tangential vector to S). Then, c D and c L are obtained by scaling F D and F L by 2/(ρ ū2 D). The Strouhal number is St = Df /ū, f being the frequency of separation, evaluated from the temporal evolution of the lift coefficient.

F D = S (ρν ∂u t ∂n n y -pn x ) dS, F L = - S (ρν ∂u t ∂n n x + pn y ) dS,

Isothermal case

The density is ρ = 1 and the viscosity is µ = 0.001, so that the Reynolds number, defined as Re = ρūD/µ, where ū = 2 u x (0, H/2)/3 = 1, is equal to 100. The time step varies from 10 -3 for the coarsest mesh to 2 10 -4 for the finest mesh, in order to keep a constant CFL number.

The computed values (c D max , c Lmax , St and ∆P ) are gathered in Tables III (Euler scheme) and IV (Crank-Nicolson scheme), together with the plausible range for the results derived from the set of contributions to the benchmark. Values entering this reference interval are typeset in bold. For the finest mesh, both schemes give results included in the benchmark reference range, except for the Strouhal number which is underestimated by the Euler scheme; in addition, with coarser meshes, the Crank-Nicolson scheme is clearly more accurate.

Anisothermal case

The cylinder is now heated at a temperature θ h . The fluid enters in the domain with a lower temperature θ c and the temperature is also fixed at θ c at the top and bottom boundaries. We characterize the deviation from the isothermal case by the parameter ε defined by

ε = θ h -θ c θ h + θ c .
The Prandtl number is equal to 0.7, the viscosity is fixed to 0.001 and the specific heat to 1000. The density is given by the following law:

ρ = ρ 0 θ c θ ,
with ρ 0 = 1. The parameter ε varies in our tests from ε = 0 (isothermal case) to ε = 0.8, and thus the ratio ρ 0 /̺(θ h ) varies from 1 (isothermal case) to 9 (ε = 0.8), since this ratio is equal to (ε + 1)/(ε -1).

Since the mesh here is unstructured, the SUSHI scheme is used for the approximation of the Laplace operator (see Remark 4.5) and the MUSCL scheme for the convective term [START_REF] Piar | A Formally Second-Order Cell Centered Scheme for Convection-Diffusion Equations on General Grids[END_REF] in the balance equation for θ. The (constant) time step satisfies a CFL condition required by the MUSCL scheme and so varies between 10 -4 and 2 10 -5 in the computations.

On Figure 10, the density field is represented for ε = 0.4, 0.6, 0.8. For small values of ε, as for the isothermal case, the flow is unsteady and oscillations appear in the wake of the cylinder. For ε = 0.4, the amplitudes seem to be greater than for ε = 0.6. This behaviour is coherent with an increase of Strouhal number when ε increases (Table V). For ε = 0.8, contrary to the other computations, the flow is steady and no oscillation appears. Heating the cylinder leads to a depressurization behind the cylinder which increases with ε (Table V). The drag coefficient also increases with ε, whereas the lift coefficient decreases. For the steady case (ε = 0.8), the lift coefficient is very small and its value is similar to the reference value given in the benchmark [START_REF] Schäfer | Benchmark Computations of Laminar Flow Around a Cylinder[END_REF] for the 2D steady case.

3D turbulent mixing layer

In this section, we address the simulation of an experimental study of a turbulent mixing layer in a vertical channel, described in [START_REF] Roig | Zone de Mélange d'Ecoulements Diphasiques à Bulles[END_REF]. The channel is divided at its bottom in two equal parts by a vertical plate, and water is injected in both parts with a different velocity. The mean velocity in the right side of the channel, u 1 = 0.61m.s -1 , is greater than in the left one u 2 = 0.26m.s -1 ; a mixing layer, the thickness of which increases along the flow, is thus generated at the center of the channel (see Figure 11). The flow is incompressible and the density and the viscosity are respectively fixed at ρ = 10 3 kg.m -3 and µ = 10 -3 Pa.s.

Large Eddy Simulation is used to describe this experiment, and the subgrid scale eddy viscosity is given by the Smagorinsky model, with C s = 0.12 [START_REF] Milelli | A Numerical Analysis of Confined Turbulent Bubble Plumes[END_REF]. The computational domain is Ω = [-0.2; 0.2] × [-0.01; 1.20] × [-0.2; 0.2] (lengthes in meters), and the top of the plate is taken at y = 0. Slip boundary conditions are prescribed on the lateral planes of Ω and outlet boundary conditions on the upper plane y = 1.20. The separation plate is not taken into account in the simulation and mean experimental velocity profiles are imposed on the bottom plane y = -0.01. To generate the turbulent inflow conditions, the vortex method [START_REF] Jarrin | A Synthetic-Eddy-Method for Generating Inflow Conditions for Large-Eddy Simulation[END_REF] is used.

The mesh is an uniform grid of 90 × 90 × 270 cells and simulations are carried out with two time steps, δt = 10 -2 s and δt = 10 -3 s, over the time interval [0s, 20s].

We compare the solutions given by the Crank-Nicolson-like pressure-correction scheme introduced here to the usual first order Euler pressure-correction scheme. Representations of the vorticity norm isosurfaces for the isovalue 6 are given on Figure 12 for several time steps. For the time step δt = 10 -2 s, Figure 12 shows that vortices present more complex structures with the Crank-Nicolson-like scheme (in other words, turbulent vortices are over-damped by the numerical diffusion generated by the first order Euler discretization). For δt = 10 -3 s, both schemes give similar results.

In the experiment, the (time-averaged) mean and root mean square (rms) axial velocity profiles are measured along the x-direction at several heights [START_REF] Roig | Zone de Mélange d'Ecoulements Diphasiques à Bulles[END_REF]. In our computations, we evaluate these quantities by averaging over a time-interval [T = 8 s; T ′ = 20 s] for which the fully turbulent state is reached. Let L be the integer number such that T ′ = T + L δt. On a cell K, the mean value of a the mean axial velocity is given by < u y > |K and the root mean square (rms) axial velocity by:

(u y ) rms =< u 2 y > |K -< u y > 2 |K .
The comparison of experimental and numerical mean velocity profiles at y = 0.5 (Figure 13) shows that the time convergence is achieved for δt = 10 -2 s for both schemes and a good agreement is observed with experimental results. Figure 14 shows the same profiles for the rms axial velocity. For this quantity, a difference is observed between both schemes at δt = 10 -2 s, the Crank-Nicolsonlike algorithm reproducing correctly the experimental peak, contrary to Euler. For δt = 10 -3 s, both schemes seem to have achieved convergence in time and, once again, results are in reasonable agreement with the experimental data.

Figure 1 .

 1 Figure 1. Notations for control volumes and diamond cells.

Figure 2 .

 2 Figure 2. Notation for the mass fluxes.

  Crank-Nicolson-like.

Figure 3 .

 3 Figure 3. Test with a known solution: time convergence in L 2 norm for the velocity, for various meshes.

Figure 4 .

 4 Figure 4. Test with a known solution: time convergence in L 2 norm for the pressure, for various meshes.

Figure 5 .

 5 Figure 5. Test with a known solution: time convergence in discrete L 2 norm for θ, for various meshes.

Figure 6 .

 6 Figure 6. Test with a known solution: magnitude of the two parts of the dissipation defect as a function of the time step for the Crank-Nicolson and the Euler time discretizations.

uxFigure 7 .

 7 Figure 7. Laminar 2D flow around a cylinder -Geometry.

Figure 8 .

 8 Figure 8. The coarse mesh (mesh#1).

Figure 9 .

 9 Figure 9. Laminar 2D flow around a cylinder -First component (x-component) of the velocity.

8 Figure 10 .

 810 Figure 10. Laminar flow around a heated cylinder -Density field isolines for several ε at time t = 10.

2 Figure 11 .

 211 Figure 11. Turbulent mixing layer -Numerical set-up (front view).

  (a) Euler, δt = 10 -2 s. (b) Crank-Nicolson-like, δt = 10 -2 s. (c) Euler , δt = 10 -3 s. (d) Crank-Nicolson-like, δt = 10 -3 s.

Figure 12 .

 12 Figure 12. Turbulent mixing layer -Vorticity norm isosurface for the isovalue 6 at time 20s.

  field v, denoted by < v > |K , is thus < v > |K = 1 T ′ -T L j=0 δt |K| K v(x, t + j δt) dx,

  Crank-Nicolson-like.

Figure 13 .

 13 Figure 13. Turbulent mixing layer -Time averaged axial velocity at y = 0.5 and z = 0.

Figure 14 .

 14 Figure 14. Turbulent mixing layer -Time rms axial velocity at y = 0.5 and z = 0.

Table II .

 II Laminar 2D flow around a cylinder -Description of the meshes used.

	Computation	mesh#1 mesh#2 mesh#3 mesh#4 mesh#5
	Number of cells	4033	12913	43009	76091	106101
	Space unknowns 12256	39014	129527 228937 409099

Table III .

 III Laminar 2D flow around a cylinder, isothermal case -Euler time discretization.

	Computation	c Dmax	c Lmax	St	∆P
	mesh#1	3.62	0.75	0.269	2.36
	mesh#2	3.40	0.92	0.288	2.45
	mesh#3	3.25	0.97	0.294	2.49
	mesh#4	3.23	0.98	0.294	2.484
	mesh#5	3.22	1.00	0.294	2.478
	Reference range [30] 3.22 -3.24 0.99 -1.01 0.295 -0.305 2.46 -2.50

Table IV .

 IV Laminar 2D flow around a cylinder, isothermal case -Crank-Nicolson time discretization.

	Computation	c Dmax	c Lmax	St	∆P
	mesh#1	3.66	0.79	0.270	2.30
	mesh#2	3.41	0.95	0.294	2.45
	mesh#3	3.25	0.98	0.303	2.50
	mesh#4	3.23	1.00	0.303	2.49
	mesh#5	3.22	1.01	0.303	2.48
	Reference range [30] 3.22 -3.24 0.99 -1.01 0.295 -0.305 2.46 -2.50

Table V .

 V Laminar flow around a heated cylinder -Drag and lift coefficients, Strouhal number and ∆P for several ε.

	Computation c Dmax c Lmax	St	∆P
	isotherm	3.22	1	0.303 2.48
	ε = 0.4	3.54	0.56	0.316 2.59
	ε = 0.6	3.84	0.236	0.342 2.71
	ε = 0.8	4.62	0.0055 steady 3.04

We consider the system (1)- [START_REF] Beaudan | Numerical Experiments on the Flow Past a Circular Cylinder at Sub-Critical Reynolds Number[END_REF], with the following expression for the density ρ as a function of the unknown θ:

With this closure laws, the system at hand may be seen as the set of equations governing a twocomponent flow, of component densities ρ 1 and ρ 2 (here supposed to be two constant positive real numbers); the unknown θ stands for the mass fraction of the first component.

Let us suppose that the normal velocity is prescribed to zero on the whole boundary of the computational domain ∂Ω. In this condition, integrating the mass balance over Ω yields the total mass balance:

On the other side, the density ρ is given as a function of θ, itself solution of Equation ( 2) for λ = 10 -4 . Integrating this relation over Ω and supposing that the mass diffusion flux of the first component vanishes at the boundary, i.e. λ∇θ • n = 0 on ∂Ω, we obtain:

By an easy manipulation of the equation of state [START_REF] Subbareddy | A Fully Discrete, Kinetic Energy Consistent Finite-Volume Scheme for Compressible Flows[END_REF], we get:

and thus, since ρ appears as an affine function of ρθ, the relations ( 32) and ( 33) are fortunately compatible. With the proposed time-stepping algorithm, this property does not hold anymore, because of the time shift of the density in the solution of the balance equation for θ. Consequently, a renormalisation of the density is necessary to ensure the existence of a solution to the projection step:

Note that this relation is reminiscent of the scaling of the density obtained through its dependency versus the so-called thermodynamical pressure in the asymptotic model for low Mach number flows [START_REF] Majda | The Derivation and Numerical Solution of the Equations for Zero Mach Number Combustion[END_REF].

We consider now the following solution to System (1)-( 2):

We suppose that ρ 1 = 0.25 and ρ 2 = 2, so relation (31) yields θ = (2ρ)/(7ρ). The viscosity is supposed to be µ = 0.001. The computational domain is Ω = (0, 1) × (0, 1), so the normal velocity is always zero at the boundary, and the integral of the density over Ω does not vary with time. With this choice for ρ and ρu, the mass balance is satisfied. The right-hand side in the momentum balance and in the θ transport equation, the initial conditions and the boundary conditions (prescribed value for u) are chosen to match the analytical solution. The non zero right-hand side for this system is numerically treated implicitly.

The domain is meshed by n × n regular grids, with n = 16, 32, 64 and 128. The discretization of the convection term for θ is performed by an upwind finite volume scheme. Errors obtained at t = 0.5 are displayed on figures 3, 4 and 5, in the L 2 norm for the velocity and the pressure, and, for


Ilya Loshchilov
email: ilya.loshchilov@epfl.ch

Marc Schoenauer
email: marc.schoenauer@inria.fr

Michèle Sebag
email: michele.sebag@lri.fr

BI-population CMA-ES Algorithms with Surrogate Models and Line Searches

Keywords: G.1.6 [Numerical Analysis]: Optimization-global optimization, unconstrained optimization; F.2.1 [Analysis of Algorithms and Problem Complexity]: Numerical Algorithms and Problems Benchmarking, black-box optimization, evolution strategy, CMA-ES, self-adaptation, surrogate models, surrogate-assisted optimization, separable optimization

In this paper, three extensions of the BI-population Covariance Matrix Adaptation Evolution Strategy with weighted active covariance matrix update (BIPOP-aCMA-ES) are investigated. First, to address expensive optimization, we benchmark a recently proposed extension of the self-adaptive surrogate-assisted CMA-ES which benefits from more intensive surrogate model exploitation (BIPOP-saACM-k). Second, to address separable optimization, we propose a hybrid of BIPOP-aCMA-ES and STEP algorithm with coordinatewise line search (BIPOP-aCMA-STEP). Third, we propose HCMA, a hybrid of BIPOP-saACM-k, STEP and NEWUOA to benefit both from surrogate models and line searches. All algorithms were tested on the noiseless BBOB testbed using restarts till a total number of function evaluations of 10 6 n was reached, where n is the dimension of the function search space.

The comparison shows that BIPOP-saACM-k outperforms its predecessor BIPOP-saACM up to a factor of 2 on illconditioned problems, while BIPOP-aCMA-STEP outperforms the original BIPOP-based algorithms on separable functions. The hybrid HCMA algorithm demonstrates the best overall performance compared to the best algorithms of the BBOB-2009, BBOB-2010 and BBOB-2012 when running for more than 100n function evaluations.

INTRODUCTION

The Covariance Matrix Adaptation Evolution Strategy (CMA-ES) [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] has become a popular tool for continuous black-box optimization. However, CMA-ES is a rather local optimization algorithm and can be ineffective when dealing with multi-modal optimization problems. In order to avoid premature convergence to local optima, several restart strategies for CMA-ES have been proposed [START_REF] Auger | A Restart CMA Evolution Strategy With Increasing Population Size[END_REF][START_REF] Hansen | Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed[END_REF][START_REF] Loshchilov | Alternative Restart Strategies for CMA-ES[END_REF]. These strategies, typically based on restarting and increasing (more generally, on varying) the population size and initial mutation step-size each time at least one stopping criterion of CMA-ES is met, implement the simplest niching approach but rather in the time than in space [START_REF] Preuss | Niching the cma-es via nearest-better clustering[END_REF]. Thanks to its population-based nature, the CMA-ES benefits from using larger population sizes on multi-modal functions with some global structure and, thus, typically outperforms quasi-Newton methods such as BFGS [START_REF] Shanno | Conditioning of Quasi-Newton Methods for Function Minimization[END_REF]. The latter, however, often outperforms CMA-ES on a special but quite popular class of noiseless unimodal functions, when gradient information is useful. This drawback of CMA-ES becomes especially undesirable when dealing with expensive optimization, where each function evaluation takes some time and/or costs some money.

A common way to reduce the cost of optimization within Evolutionary Algorithms is to periodically learn a surrogate model of the expensive function and exploit it during the search by pre-filtering promising solutions and/or by optimizing the surrogate directly instead of the expensive function. These approaches have been adopted within Evolution Strategies and CMA-ES by using different surrogate learning approaches: Radial Basis Functions network [START_REF] Hoffmann | Controlled Model Assisted Evolution Strategy with Adaptive Preselection[END_REF], Gaussian Processes [START_REF] Ulmer | Evolution Strategies assisted by Gaussian Processes with Improved Pre-Selection Criterion[END_REF], Artificial Neural Network [START_REF] Jin | A Comprehensive Survey of Fitness Approximation in Evolutionary Computation[END_REF], Support Vector Regression [START_REF] Kramer | Covariance Matrix Self-Adaptation and Kernel Regression -Perspectives of Evolutionary Optimization in Kernel Machines[END_REF], Local-Weighted Regression [START_REF] Kern | Local Meta-Models for Optimization Using Evolution Strategies[END_REF][START_REF] Bouzarkouna | Investigating the Local-Meta-Model CMA-ES for Large Population Sizes[END_REF], Ranking Support Vector Machine (Ranking SVM) [START_REF] Kern | Local Meta-Models for Optimization Using Evolution Strategies[END_REF][START_REF] Loshchilov | Comparison-Based Optimizers Need Comparison-Based Surrogates[END_REF][START_REF] Ingimundardottir | Sampling Strategies in Ordinal Regression for Surrogate Assisted Evolutionary Optimization[END_REF]. More recently it has been shown [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] that key invariance properties of CMA-ES can be preserved even in surrogateassisted scenario by using Ranking SVM approach. Moreover, hyper-parameters used to build surrogate models can be also adapted during the search allowing the user to define only the range of these hyper-parameters, while their actual values will be adapted by the s * aACM-ES algorithm [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF].

In this paper, we benchmark the recently proposed extension of s * aACM-ES, referred to as s * aACM-ES-k, where the surrogate model is more intensively exploited by increasing the population size used for its optimization. This extension further improves the performance of the surrogate-assisted CMA-ES on uni-modal functions, however, almost does not improve the performance of multi-modal functions. This is especially true for separable multi-modal functions, where CMA-ES is outperformed by Evolutionary Algorithms which explicitly or implicitly exploit separability. In order to investigate that would be the effect of such an exploitation by CMA-ES, we present a hybrid version of CMA-ES algorithm and select the easiest point (STEP) [START_REF] Swarzberg | Step: The easiest way to optimize a function[END_REF][START_REF] Pošík | Restarted local search algorithms for continuous black box optimization[END_REF] algorithm, which performs coordinate-wise search and is quite efficient on separable functions [START_REF] Swarzberg | Step: The easiest way to optimize a function[END_REF][START_REF] Pošík | Restarted local search algorithms for continuous black box optimization[END_REF]. Finally, to benefit both from surrogate models and line searches, we propose a hybrid CMA (HCMA) which combines s * aACM-ES-k, STEP and NEWUOA algorithm [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF].

THE ALGORITHMS

The (µ/µw, λ)-CMA-ES

In each iteration t, (µ/µw, λ)-CMA-ES [START_REF] Hansen | Reducing the time complexity of the derandomized evolution strategy with covariance matrix adaptation (CMA-ES)[END_REF] samples λ new solutions xi ∈ R n , where i = 1, . . . , λ, and selects the best µ among them. These µ points update the distribution of parameters of the algorithm to increase the probability of successful steps in iteration t + 1. The sampling is defined by a multi-variate normal distribution, N (m t , σ t 2 C t ), with current mean of distribution m t , n × n covariance matrix C t and step-size σ t . Most of algorithms which will be benchmarked in this paper are based on an active version of the CMA-ES (in short, aCMA-ES) with negative covariance matrix update [START_REF] Jastrebski | Improving Evolution Strategies through Active Covariance Matrix Adaptation[END_REF][START_REF] Hansen | Benchmarking a weighted negative covariance matrix update on the BBOB-2010 noiseless testbed[END_REF].

The BIPOP-aCMA-ES

In BIPOP-CMA-ES [START_REF] Hansen | Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed[END_REF] after the first single run (when at least one of stopping criteria is met) with default population size, the CMA-ES is restarted in one of two possible regimes by taking into account the budget of function evaluations spent in the corresponding regime. Each time the algorithm is restarted, the regime with smallest budget used so far is used.

Under the first regime the population size is doubled as λ large = 2 irestart λ def ault in each restart irestart and uses some fixed initial step-size σ 0 large = σ 0 def ault . Under the second regime the CMA-ES is restarted with some small population size λ small and step-size σ 0 small , where λ small is set to

λ small = λ def ault 1 2 λ large λ def ault U [0,1] 2 , (1)
Here U [0, 1] denotes independent uniformly distributed numbers in [0, 1] and 1] . The active version of BIPOP-CMA-ES (BIPOP-aCMA-ES) has been proposed in [START_REF] Loshchilov | Alternative Restart Strategies for CMA-ES[END_REF].

λ small ∈ [λ def ault , λ/2]. The initial step-size is set to σ 0 small = σ 0 def ault × 10 -2U [0,

The BIPOPs * ACM-ES

The s * ACM-ES [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF] is the surrogate-assisted version of the (µ/µw, λ)-CMA-ES, where Ranking SVM-based surrogate model is used periodically instead of the expensive function for direct optimization. The use of Ranking SVM allows to preserve invariance properties of CMA-ES w.r.t. rankpreserving transformations of the objective function and orthogonal transformations of the search space. The main loop of s * ACM-ES: the surrogate model f is optimized for n generations by the CMA-ES, then the expensive function f is optimized for one generation. To adjust the number of generations n for the next time, the model error is computed as a fraction of incorrectly predicted comparison relations that is observed after ranking the last λ evaluated points according to f and f . The s * ACM-ES performs an online optimization of the surrogate model hyper-parameters during the optimization of the objective function by performing a search in a space of model hyper-parameters, and by generating λ hyp surrogate models in each iteration. The information about the fittest models with the smaller prediction error on the last λ evaluated solutions is used to adjust the surrogate hyper-parameters for the next iteration. The detailed description of s * ACM-ES is given in [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF].

The BIPOPs * ACM-ES-k

A more intensive exploitation of a surrogate model can be useful if the model is sufficiently precise. An extension of s * ACM-ES, based on this idea, referred to as s * ACM-ES-k, has been proposed in [START_REF] Loshchilov | Intensive Surrogate Model Exploitation in Self-adaptive Surrogate-assisted CMA-ES (saACM-ES)[END_REF]. In s * ACM-ES-k, f is optimized for n generations by CMA-ES with population size λ = k λ λ def ault and number of parents µ = kµµ def ault , where k λ ≥ 1 and kµ ≥ 1. Then, f is optimized for 1 generation by CMA-ES with population size λ = λ def ault and number of parents µ = µ def ault . The original s * ACM-ES corresponds to s * ACM-ES-k with k λ = 1 and kµ = 1. Thus, larger values of k λ and kµ lead to a more intensive exploitation of f . To prevent the algorithm from a potential divergence when n oscillates around 0 and 1, k λ > 1 is used only if n ≥ nk λ , where nk λ is the number of generations which corresponds to an "accurate enough" model to be intensively exploited.

The BIPOP-aCMA-STEP

The select the easiest point (STEP) [START_REF] Swarzberg | Step: The easiest way to optimize a function[END_REF][START_REF] Pošík | Restarted local search algorithms for continuous black box optimization[END_REF] is a simple line search method which is based on iterative interval division. The STEP determines the next point to evaluate by analyzing the usefulness of evaluation the function at the middle of an interval [xi, xi-1] with respect to the interval difficulty computed as [START_REF] Swarzberg | Step: The easiest way to optimize a function[END_REF]

D = 4ŷ -2∆y + 4 ŷ2 -ŷ∆y ∆x 2 , ( 2
)
where ∆x = xi -xi-1, ∆y = f (xi) -f (xi-1), ŷ = f * - f (xi-1) + l, f *
is the best solution found so far by the STEP and l is the tolerance (i.e., the precision required for the value of the optimum).

The STEP can be applied to n-dimensional optimization problem by performing n coordinate-wise iterative or parallel searches. A variant of iterative search: optimize f along the coordinate i for NST EP function evaluations or until a stopping criterion is met, then fix the best xi, increment i and repeat the procedure again while i ≤ n. This approach is quite straightforward but requires NST EP to be set a priori. However, NST EP can be initialized by a small value and then incremented (e.g., doubled) in a similar way as the population size in BIPOP-aCMA-ES each time STEP finishes the search on n coordinates. This approach might lead to a loss of a factor of about 2 in terms of number of function evaluations compared to an optimal setting of NST EP . Moreover, some step function-like improvement of the objective function can be observed due to switch between variables, this also means that objective values of tested solutions are far away in the objective space during the coordinate-wise search on all but n-th coordinate.

We propose a parallel STEP search where at each iteration one step per coordinate is performed, then x * i which corresponds to the best solution so far along the coordinate i is used to construct a recommended solution x * combined of x * i , i = 1, . . . , n. The objective value of this solution is estimated and is used for two purposes. First, if the value of this solution is worse (up to some precision p) than the one of a solution recommended in the previous iteration, then f is non-separable or/and noisy. In this case, it is not recommended to use the STEP. Second, the new recommended solution can be of the same quality as the previous one even if the problem is separable. In this case, the STEP still can run while its best recommended solution is better than the one found by its competitor algorithm running in parallel with, e.g., the same budget of function evaluations.

In BIPOP-aCMA-STEP algorithm, a fraction ρST EP of function evaluations is allocated to the STEP, while 1-ρST EP is allocated to the BIPOP-aCMA-ES. At each iteration after the first nMinIterST EP iterations, the best solutions of two algorithms are compared and the STEP is stopped if its best solution is worse or if its current solution is worse than its previously recommended solution.

The HCMA

The HCMA algorithm is a hybrid of BIPOP-s * aACM-ESk and STEP algorithms coupled as described in the previous section. Additionally, the NEWUOA algorithm [START_REF] Powell | The NEWUOA software for unconstrained optimization without derivatives[END_REF] is used for the first 10n function evaluations to insure good convergence of the HCMA on simple functions which can be efficiently approximated by quadratic surrogate models. All solutions generated by STEP and NEWUOA are not used by BIPOP-s * aACM-ES-k for surrogate learning since the sampling/test distribution of CMA-ES might not correspond to the training distribution defined by these solutions.

The Benchmarked Algorithms

For benchmarking we consider five CMA-ES algorithms in BIPOP restart scenario [START_REF] Hansen | Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed[END_REF]: BIPOP-aCMA-ES [START_REF] Loshchilov | Black-box Optimization Benchmarking of NIPOP-aCMA-ES and NBIPOP-aCMA-ES on the BBOB-2012 Noiseless Testbed[END_REF][START_REF] Hansen | Benchmarking a BI-population CMA-ES on the BBOB-2009 function testbed[END_REF], BIPOP-s * aACM-ES [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF], BIPOP-s * aACM-ES-k [START_REF] Loshchilov | Intensive Surrogate Model Exploitation in Self-adaptive Surrogate-assisted CMA-ES (saACM-ES)[END_REF], BIPOP-aCMA-STEP and HCMA 1 . For all but BIPOP-s * aACM-ESk [START_REF] Loshchilov | Intensive Surrogate Model Exploitation in Self-adaptive Surrogate-assisted CMA-ES (saACM-ES)[END_REF] and BIPOP-aCMA-STEP algorithms the BBOB results are taken from the literature. The setting of parameters of BIPOP-s * aACM-ES-k is given in [START_REF] Loshchilov | Intensive Surrogate Model Exploitation in Self-adaptive Surrogate-assisted CMA-ES (saACM-ES)[END_REF], the maximum number of function evaluations is set to 10 6 n. However, surrogate models are used only for the first 2 • 10 5 (10 4 n for n = 20) function evaluations for n = 2, 3, 5, 10, 20 and 2 • 10 3 n evaluations for n = 40 in order to reduce the overall running time and to fit to the BBOB-2013 focus on expensive optimization with the budget of up to 10 3 n function evaluations. The value of k λ is set to 1, 1, 1, 10, 100, 1000 for n = 2, 3, 5, 10, 20, 40, respectively; kµ = 1 for all dimensions; nk λ = 4. The value of l (respectively, p) of the STEP is set to 10 -10 (respectively, 10 -8 ), nMinIterST EP = 10 and ρST EP = 0.5. The following parameters of NEWUOA are used [START_REF] Ros | Benchmarking the newuoa on the bbob-2009 function testbed[END_REF]: m = 2n + 1 interpolation points, the initial radius ρ beg of the search region is set to 10, the final radius ρ end is set to 10 -16 , the maximum number of function evaluations is set to 10n.

RESULTS

Results from experiments according to [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF] on the benchmark functions given in [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF] are presented in Figures 1 and2 and in Table 1. The expected running time (ERT), used in the figures and table, depends on a given target function value, ft = fopt + ∆f , and is computed over all relevant trials (on the first 15 instances) as the number of function evaluations executed during each trial while the best function value did not reach ft, summed over all trials and divided by the number of trials that actually reached ft [START_REF] Hansen | Real-Parameter Black-Box Optimization Benchmarking[END_REF]. Statistical significance is tested with the rank-sum test for a given target ∆ft (10 -8 as in Figure 1) using, for each trial, either the number of needed function evaluations to reach ∆ft (inverted and multiplied by -1), or, if the target was not reached, the best ∆f -value achieved, measured only up to the smallest number of overall function evaluations for any unsuccessful trial under consideration.

The BIPOP-s * aACM-ES and BIPOP-s * aACM-ES-k perform similarly for n = 2, 3, 5 and before 10 4 n function evaluations because k λ = 1 is used, after 10 4 n evaluations BIPOPs * aACM-ES-k represents BIPOP-aCMA-ES. For larger n, BIPOP-s * aACM-ES-k outperforms BIPOP-s * aACM-ES on unimodal ill-conditioned problems, e.g., for 20-dimensional ones the target function values 10 -8 can be reached faster (see Figure 1, Table 1): by a factor of 1.8 on f2 Separable Ellipsoid and f10 Rotated Ellipsoid functions, by a factor of 1.5 on f11 Discus, by a factor of 1.1 on f12 Bent cigar, by a factor of 1.8 on f13 Sharp ridge and by a factor of 1.2 on f14 Sum of different powers. On f8 and f9 Rosenbrock's function the speedup of a factor of about 1.3 is obtained.

A relative loss of performance (compared to BIPOP-s * aACM-ES) on multi-modal functions can be explained by not using surrogate models after 10 4 n and thus using the original BIPOP-aCMA-ES which usually is outperformed by BIPOP-s * aACM-ES. A loss of performance is observed on f6 Attractive sector (by a factor of 1.2), f17 Schaffer (by a factor of 1.2) and f23 Katsuuras (by a factor of 1.5). However, the reason of this loss remains unclear for f17 and f23 functions because a similar loss of performance is also observed on the same problems in dimension 5 where k λ = 1 for both algorithms (i.e., there is no intensive surrogate model exploitation). This difference might be explained by the stochasticity of multi-modal optimization or by fact of adapting an additional hyper-parameter -a stopping criterion of SVM surrogate learning procedure in BIPOP-s * aACM-ES-k [START_REF] Loshchilov | Intensive Surrogate Model Exploitation in Self-adaptive Surrogate-assisted CMA-ES (saACM-ES)[END_REF], whose values were fixed in BIPOP-s * aACM-ES [START_REF] Loshchilov | Self-Adaptive Surrogate-Assisted Covariance Matrix Adaptation Evolution Strategy[END_REF]. The adaptation of this hyper-parameter probably is not efficient enough compared to the offline tuned value, but this adaptation makes the algorithm more parameter-free. While there is no data to assess the effect of more intensive surrogate model exploitation for budgets larger than 10 4 n, the loss of performance on f6 is clearly due to k λ > 1. A simple way to reduce this loss would be to increment nk λ , but this might decrease the speedup on ill-conditioned problems back toward the level of BIPOP-s * aACM-ES.

The results shown in Figure 2 suggest that separable multimodal functions such as f3 Separable Rastrigin and f4 Separable Skew Rastrigin-Bueche can be easily solved by BIPOP-aCMA-STEP and HCMA thanks to its STEP part. Thus, this kind of hybridization is beneficial if the main target is to solve the maximum number of problems within a relatively large budget of function evaluations.

The HCMA demonstrates the best performance among the tested algorithms after 100n function evaluations. The hybridization with STEP allows to exploit the separability when it is detected, i.e., when STEP generates better solutions than BIPOP-s * aACM-ES-k. In a few cases when separability was not detected successfully 3-dimensional f3 Rastigin and 3-, 5-dimensional f4 Skew Rastrigin-Bueche separable functions, it seems that 10n function evaluations (nMinIterST EP = 10) sometimes is not sufficient to see the difference between STEP and BIPOP-s * aACM-ES-k on these problems.

CPU TIMING EXPERIMENT

The time complexity of BIPOP-aCMA-STEP is essentially the same as the one of BIPOP-aCMA-ES since the STEP has linear time complexity and its relatively small "constant factor" (see Eq.( 2)) depends only on a particular implementation. The time complexity of BIPOP-s * aACM-ES-k was measured after 15 runs on 20-dimensional f8, the cost per function evaluation is about 0.45 seconds (runs were performed on different machines with 2.2 and 2.4 GHz cores under Ubuntu 10.04 using Octave 3.2.3). These results are similar to the ones of [START_REF] Loshchilov | Black-box Optimization Benchmarking of IPOP-saACM-ES and BIPOP-saACM-ES on the BBOB-2012 Noiseless Testbed[END_REF], the time complexity can be reduced by a factor of λ hyp = 20 by switching off the adaptation of surrogate hyper-parameters.

CONCLUSION

In this paper, we have compared the recently proposed self-adaptive surrogate-assisted BIPOP-s * aACM-ES-k with intensive surrogate model exploitation and its predecessor BIPOP-s * aACM-ES with the BIPOP-aCMA-ES. The comparison shows that the performance of surrogate-assisted CMA-ES can be further improved on ill-conditioned functions, e.g., on 20-dimensional f10 Rotated Ellipsoid the BIPOPs * aACM-ES-k algorithm is about 6 times faster than BIPOP-aCMA-ES and about 1.8 times faster than BIPOP-s * aACM-ES. However, the performance on multi-modal functions should be further analyzed for a larger maximum number of function evaluations.

The proposed BIPOP-aCMA-STEP as expected shows good results on separable multi-modal functions by exploitation when it is detected (when STEP performs better than BIPOP-aCMA-ES). The proposed hybrid HCMA algorithm demonstrates the best overall performance when running for more than 100n function evaluations on 10-, 20-and 40dimensional problems, its overall performance for this budget of function evaluations and for these dimensions is better than of all algorithms tested during the BBOB-2009, BBOB-2010 and BBOB-2012.

It would be interesting to further develop HCMA and BIPOP-aCMA-STEP and make these algorithms better suited to partially-separable functions and/or adjust the algorithm to work in a sampling space of CMA-ES defined by its covariance matrix. A natural perspective is to improve the results of HCMA on multi-modal function by applying the recently proposed alternative restart strategies for CMA-ES [START_REF] Loshchilov | Black-box Optimization Benchmarking of NIPOP-aCMA-ES and NBIPOP-aCMA-ES on the BBOB-2012 Noiseless Testbed[END_REF].

.7e4(2e4) ∞ ∞ ∞ ∞ 2e7 0/15 BIPOP-s 10(7) ∞ ∞ ∞ ∞ ∞ 2e7 0/5 BIPOP-s 20(23) ∞ ∞ ∞ ∞ ∞ 2e7 0/
.0e4(3e4) ∞ ∞ ∞ ∞ ∞ 2e7 0/15 BIPOP-s ∞ ∞ ∞ ∞ ∞ ∞ 2e7 0/5 BIPOP-s 1.9e4(2e4) ∞ ∞ ∞ ∞ ∞ 2e7 0/
BIPOP-a 1.6(0.

BIPOP-s 1.0(0.9) 1.6(0.6) 0.84(0.3) 0.61(0.1) ↓3 0.61(0.1) ↓3 0.60(0.1) ↓3 15/15 BIPOP-s 0.52(0.2) ↓2 1.2(0.9) 0.65(0.3) 0.56(0.2) 0.56(0.2) 0.55(0.2) 15/15 BIPOP-a 1.5(

Figure 1 :

 1 Figure 1: Expected running time (ERT in number of f -evaluations) divided by dimension for target function value 10 -8 as log 10 values versus dimension. Different symbols correspond to different algorithms given in the legend of f1 and f24. Light symbols give the maximum number of function evaluations from the longest trial divided by dimension. Horizontal lines give linear scaling, slanted dotted lines give quadratic scaling. Black stars indicate statistically better result compared to all other algorithms with p < 0.01 and Bonferroni correction number of dimensions (six). Legend: •:BIPOP-aCMA, ▽:BIPOP-saACM, ⋆:BIPOP-saACM-k, 2:BIPOP-aCMA-STEP, △:HCMA.

Figure 2 :

 2 Figure 2: Bootstrapped empirical cumulative distribution of the number of objective function evaluations divided by dimension (FEvals/D) for 50 targets in 10 [-8..2] for all functions and subgroups in 20-D. The "best 2009" line corresponds to the best ERT observed during BBOB 2009 for each single target.

  1.0(0.1) ⋆ 0.97(0.1) ⋆2 1.0(0.1) ⋆2 1.1(0.1) ⋆2 1.1(0.1) ⋆2 1.1(0.1) ⋆2 15/15 BIPOP-a 3.7

Table 1 :

 1 Expected running time (ERT in number of function evaluations) divided by the respective best ERT measured during BBOB-2009 (given in the respective first row) for different ∆f values in dimension 20. The central 80% range divided by two is given in braces. The median number of conducted function evaluations is additionally given in italics, if ERT(10 -7 ) = ∞. #succ is the number of trials that reached the final target fopt + 10 -8 . Best results are printed in bold.

For the sake of reproducibility the source code is available at http://sites.google.com/site/bbobgecco2013/