Gautier Brèthes

Alain Dervieux
email: alain.dervieux@inria.fr

Anisotropic Norm-Oriented Mesh Adaptation for a Poisson problem

Keywords: goal-oriented mesh adaptation, anisotropic mesh adaptation, adjoint, metric, Poisson problem, finite elements

We present a novel formulation for the mesh adaptation of the approximation of a Partial Differential Equation (PDE). The discussion is restricted to a Poisson problem. The proposed norm-oriented formulation extends the goal-oriented formulation since it is equation-based and uses an adjoint. At the same time, the norm-oriented formulation somewhat supersedes the goal-oriented one since it is basically a solution-convergent method. Indeed, goal-oriented methods rely on the reduction of the error in evaluating a chosen scalar output with the consequence that, as mesh size is increased (more degrees of freedom), only this output is proven to tend to its continuous analog while the solution field itself may not converge. A remarkable quality of goal-oriented metric-based adaptation is the mathematical formulation of the mesh adaptation problem under the form of the optimization, in the well-identified set of metrics, of a well-defined functional. In the new proposed formulation, we amplify this advantage. We search, in the same well-identified set of metrics, the minimum of a norm of the approximation error. The norm is prescribed by the user and the method allows addressing the case of multi-objective adaptation like, for example in aerodynamics, adaptating the mesh for drag, lift and moment in one shot. In this work, we consider the basic linear finite-element approximation and restrict our study to L 2 norm in order to enjoy second-order convergence. Numerical examples for the Poisson problem are computed.

Introduction

This paper addresses anisotropic mesh adaptation. We focus on methods which build a somewhat optimal mesh defined by a parametrization using a Riemannian metric. A typical family of optimal metric-based methods for CFD is the family of Interpolation-based/Hessian-based methods. An attractive property of these methods is that they are based on a mathematical optimization principle.

Iso-distribution /equi-repartition Hessian-based methods tend to minimize a Sup or L ∞ norm of the (main term of) interpolation error with respect to a metric considered in a subset of metrics with a prescribed number of vertices. We refer to the two pioneering works [START_REF] Castro-Díaz | Anisotropic unstructured mesh adaptation for flow simulations[END_REF][START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent CFD[END_REF] for the methods, the two pioneering works [START_REF] Agouzal | Adaptive generation of quasi-optimal tetrahedral meshes[END_REF][START_REF] Vasilevski | An adaptive algorithm for quasi-optimal mesh generation[END_REF] for the analysis, and to [START_REF] Vasilevski | Error bounds for controllable adaptive algorithms based on a Hessian recovery[END_REF][START_REF] Alauzet | Adaptation de maillage anisotrope en trois dimensions[END_REF]. The methods minimizing the L p norm (p < ∞) of the interpolation error of one or several sensors depending on the CFD solution allows to better capture features of different scales in the solution. Cf. [START_REF] Yu | An unified asymptotical analysis of interpolation errors for optimal meshes[END_REF][START_REF] Huang | Metric tensors for anisotropic mesh generation[END_REF][START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in L p -norm[END_REF][START_REF] Loseille | Adaptation de maillage anisotrope 3D multi-échelles et ciblée à une fonctionnelle pour la mécanique des fluides. Application à la prédiction haute-fidélité du bang sonique[END_REF][START_REF] Alauzet | Multi-model and multi-scale optimization strategies. application to sonic boom reduction[END_REF]. Sensors are solutiondependant fields chosen by the user according to their ability to take into account mesh-resolution difficulties of the flow to compute.

The Hessian-based methods are particularly well adapted to the finite-element approximation of second-order elliptic PDEs. It is true, by the projection theorem, that a norm of the approximation error is bounded by the analog norm of the interpolation error but this concerns the H 1 norm while Hessian-based method concentrate on L ∞ or L p norms. More generally, while taking into account the features of the PDE solution, these methods do not take into account the features of the PDE itself. This is penalizing in the case of systems, for which several sensors need to be chosen and weighted by the user. However, if sensors are wisely chosen, a good convergence of the whole approximate solution field to the exact solution field is usually observed.

Taking into account the influence of the PDE on the error through an equation-based estimate has also been an important research topic. The formulation of goal-oriented methods is an important step for a more justified error evaluation. It has been introduced in [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF]. It relies on an a posteriori estimate. A good synthesis concerning a posteriori estimates is [START_REF] Verfürth | A Posteriori Error Estimation for Finite Element Methods[END_REF], see also [START_REF] Ern | Polynomial-degree-robust a posteriori estimates in a unified setting for conforming, nonconforming, discontinuous Galerkin, and mixed discretizations[END_REF]. An interest of a posteriori estimate is that it is expressed in terms of the approximate solution, assumed to be available in a mesh adaptation loop. A second interest is that it does not require the use of higherorder (approximate) derivatives, in contrast to truncation analyses. These estimates show accurately where the mesh should be refined. A method for deducing a better anisotropic mesh from an a posteriori estimate is proposed in [START_REF] Formaggia | Anisotropic a priori error estimates for elliptic problems[END_REF], while a theory for H p norms in [START_REF] Agouzal | Minimization of gradient errors of piecewise linear interpolation on simplicial meshes[END_REF] and a joint analysis of H p and L p norms of the error are presented in [START_REF] Agouzal | Hessian free metric based mesh adaptation via geometry of interpolation error[END_REF]. These methods cannot focus on an arbitrary userspecified error norm but relies on a particular one, specified by the variational formulation of the PDE. A more popular option is to choose, as accuracy target, a particular scalar output depending on the PDE solution. Any scalar output can be considered, except that difficulties can arise for the so-called non-admissible ones, according to [START_REF] Arian | Admitting the inadmissible: Adjoint formulation for incomplete cost functionals in aerodynamic optimization[END_REF]. An a posteriori estimate also allows for building correctors for goal-oriented analyses [START_REF] Giles | Improved lift and drag estimates using adjoint Euler equations[END_REF][START_REF] Pierce | Adjoint recovery of superconvergent functionals from PDE approximations[END_REF]. In [START_REF] Venditti | Adjoint error estimation and grid adaptation for functional outputs: application to quasi-one-dimensional flow[END_REF], the goal-oriented approach is cleverly combined with the correction strategy of [START_REF] Pierce | Adjoint recovery of superconvergent functionals from PDE approximations[END_REF] and with the Hessian-based metric approach, still minimizing the interpolation error of a user-prescribed sensor.

A priori estimates generally rely on Taylor series, either through divided differences or through polynomial approximation of functions. Then, approximations of higher-order derivatives of solution need to be built from the approximate solution. This is a delicate job since there are not many proofs ensuring that a good approximation of a higher-order derivative of the exact solution can be recovered from the low-order approximate solution. However, since the development of the first recovery methods (see for example [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF]), many numerical experiments tend to show that the method is useful and rather reliable.

A remarkable feature of the goal-oriented metric-based adaptation of [START_REF] Loseille | Adaptation de maillage anisotrope 3D multi-échelles et ciblée à une fonctionnelle pour la mécanique des fluides. Application à la prédiction haute-fidélité du bang sonique[END_REF][START_REF] Belme | A fully anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] is the complete and coherent mathematical formulation of the mesh adaptation problem. Indeed, it takes the form of the optimization of a well-defined functional, namely the error for a prescribed scalar output, to be minimized with respect to a parameter, the metric, belonging to a well-identified and compact set. This strategy is applied to the discrete case in [START_REF] Yano | An optimization framework for anisotropic simplex mesh adaptation: application to aerodynamics flows[END_REF]. In [START_REF] Loseille | Adaptation de maillage anisotrope 3D multi-échelles et ciblée à une fonctionnelle pour la mécanique des fluides. Application à la prédiction haute-fidélité du bang sonique[END_REF][START_REF] Belme | A fully anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF], in order to analytically solve the optimum, an a priori analysis is developed. It restricts to the main asymptotic term of the local error in order to exhibit more easily the dependance with respect to metric.

Goal-oriented methods have strongly impacted the applications but, due to its formulation, a goal-oriented method has two inherent limitations. First, it does not naturally extend to several scalar outputs. This "multi-target" issue is well-known and a proposition for addressing it is made in [START_REF] Hartman | Multitarget error estimation and adaptivity in aerodynamic flow simulations[END_REF]. Second, because they are specialized to a given scalar output, the features of the solution field which are not influencing this output might be neglected by the automatic mesh improvement. A goal-oriented method provides the convergence of the approximate prescribed scalar output to its continuous analog. But generally convergence does not hold for the whole flow field itself. To clarify this point, let us consider the mesh adaptative computation of a sonic boom footprint at the ground. The functional depends only on the pressure at the ground level. Now, many details of the flow on upper part of the aircraft do not influence the pressure at the ground level. This vanishing influence is taken into account by the adjoint state which also vanishes on these upper regions. Then, in these regions, the adapted mesh is not refined and the approximation of the flow field does not converge. See an illustration in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. As already noted, another limitation of a goal-oriented method is the scalar character of the error to reduce. It leads to use integrals of solution fields as in (u -u h , g), u being the exact solution, u h its approximation and g a field prescribed by user. Now, these integrals are generally not sufficiently sensitive to oscillating deviations between u and u h .

In the new norm-oriented formulation proposed in this paper, the user can prescribe a norm or a semi-norm |u-u h | of the error in order to minimize it with respect to the mesh. As a typical example of semi-norm, this can be the sum of square deviations on particular outputs. Let us take an example in aerodynamics. The semi-norm |u -

u h | ≡ |C l (u) -C l (u h )| 2 + |C d (u) -C d (u h )| 2 + |C m (u) -C m (u h )| 2
will account for minimizing the errors on lift, drag and moment measured from flow solution u h with respect to mesh. The proposed method will ultimately address this kind of semi-norm, assuming that, as for the goal-oriented method, the possible issue of a non-admisssible norm according to [START_REF] Arian | Admitting the inadmissible: Adjoint formulation for incomplete cost functionals in aerodynamic optimization[END_REF] is solved. As for the goal-oriented method, the norm-oriented method takes into account the PDE features and, in the case where a norm is prescribed, it produces an approximate solution field which does converge to the exact one in this norm.

Although the proposed method is a rather general method extending to more complex CFD models, see for example [START_REF] Loseille | Anisotropic norm-oriented mesh adaptation for compressible flows[END_REF], we consider in this paper a 2D Poisson problem discretized by the usual linear finite-element method. This choice is motivated first by the rather complete set of theoretical works available for the finite-element approximation of a Poisson problem. This amount of theoretical background reduces as much as possible (although far from completely) the heuristics to introduce in building the mesh adaptation analysis. In contrast, convergence results (in which functional spaces?) are not available for Euler equations, for example. A second motivation is the easy availability of exact solutions defined in a simple way. This allows to build a kind of benchmark allowing to compare mesh adaptation methods. Let us mention also that the Poisson problem (with variable coefficient) is a central equation in two-fluid models. The proposed benchmark will be inspired by a few particular two-fluid configurations.

Although the FEM is designed for minimizing the H 1 norm (with first-order accuracy), the user may wish to enjoy a convergence with a different norm. The proposed method, demonstrated for L 2 norm (providing second-order accuracy), is defined in order to many types of norms (or semi-norms). The method relies on the use of a corrector field and on an a priori error estimate from which is extracted the asymptotically largest terms of the local error. After we have the best mesh for the prescribed norm, are we in the best of all possible worlds? Not really. We need that this norm be less than a prescribed level. The central difficulty, in our opinion, is the accurate evaluation of this error norm. Most estimates are conditioned by mesh convergence. Therefore, we compare our estimates with mesh convergence analysis and try to propose a least uncertain error evaluation. In our test cases, this evaluation is compared with the information produced by the analytic exact solution.

After a formulation of the problem, the derivation of two correctors is proposed in Sec.2. Next three sections are devoted to the three identified mesh adaptation formulations: Hessian-based, minimizing an interpolation error in Sec.3, goal-oriented formulation in Sec.4 and our proposal for a norm-oriented in Sec.5. Sec.6 is devoted to a numerical comparison between the two field-convergent formulations, viz. Hessian-based and norm-oriented and the paper is completed by a discussion of methods and numerical examples.

Two correctors for the Poisson problem

On the way for minimizing the error norm ||u -u h || L 2 , where u is the exact PDE solution and u h its approximation, we have to replace u -u h by a so-called corrector which would be rather easy to evaluate.

Notations

Let V = H 1 0 (Ω), Ω being a sufficiently smooth computational domain of R 2 . The continuous PDE system is written in short:

u ∈ V | Au = f or u ∈ V | ∀ ϕ ∈ V, a(u, ϕ) = (f, ϕ) (1)
where

A = -div( 1 ρ ∇) ; a(u, ϕ) = Ω 1 ρ ∇u • ∇ϕ dxdy
and where 1 ρ is a positive, possibly discontinuous, scalar field in L ∞ (Ω). Further, we assume that the bilinear form a is coercive in space V , i.e. there exists a positive α such that:

a(v, v) ≥ α|v| 2 V .
This model exemplifies the pressure equation in some multi-phase incompressible flow formulation for which mesh adaptation is useful, see e.g. [START_REF] Guégan | An L ∞ -L p mesh adaptive method for computing unsteady bi-fluid flows[END_REF]. Let Ω h = Ω for simplicity, τ h a triangulation of Ω h and V h be the usual P 1 -continuous finite-element approximation space related to τ h :

V h = {ϕ h ∈ C 0 ( Ω) ∩ V, ϕ h | T is affine ∀T ∈ τ h }.
The finite-element discretisation of (1) is written in variational and operational form:

u h ∈ V h and ∀ ϕ h ∈ V h , a(u h , ϕ h ) = (f h , ϕ h ) (2)
in such a way that u h is a linear function of f h which we denote:

u h = A -1 h f h .
We denote by Π h the usual interpolation operator:

∀v ∈ C 0 ( Ω) ∩ V, Π h v ∈ V h and, ∀ x i vertex of Ω h , Π h v(x i ) = v(x i ).
Scalar correctors, i.e. correctors for scalar outputs j(u h ) depending on the solution, e.g. j(u h ) = (g, u h ) with g prescribed, have been defined by Giles and Pierce, [START_REF] Giles | Improved lift and drag estimates using adjoint Euler equations[END_REF]. Our interest concerns the correction of the unknown field itself. Two options are now proposed.

A priori corrector for the PDE solution

A first rather simple a priori corrector can be derived from an analysis of the error RHS. We observe that:

a(u -u h , ϕ h ) = (f -f h , ϕ h ) ∀ϕ h ∈ V h .
Assuming that the solution u is continuous, we get:

a(Π h u -u h , ϕ h ) = a(Π h u -u, ϕ h ) + (f -f h , ϕ h ) ∀ϕ h ∈ V h . (3)
We call Π h u -u h the implicit error. By implicit, we mean that it can be obtained through solving a discrete system. It differs from the approximation error by an interpolation error:

u -u h = Π h u -u h + u -Π h u. (4)
In order to find an approximate of the implicit error, we need to evaluate the RHS of (3) for any test function ϕ h . The second term of RHS of (3) is easy to evaluate (we know f and f h ). The first term of RHS of (3) can be transformed as follows (assuming that 1 ρ is constant on each element):

a(Π h u -u, ϕ h ) = T T 1 ρ ∇ϕ h • ∇(Π h u -u) dxdy = T ∂T (Π h u -u) 1 ρ ∇ϕ h • n dσ.
Then, we get:

a(Π h u -u h , ϕ h ) = ∂T ij 1 ρ ∇(ϕ h | T i -ϕ h | T j ) • n ij ∂T ij (Π h u -u) dσ + (f -f h , ϕ h ) (5)
where the sum is taken for all edges ∂T ij separating triangles T i and T j of the triangulation. The unit vector n ij normal to ∂T ij is pointing outward T i . Now, we do not know u but u h . In order to evaluate the interpolation error, we first approximate the Hessian of u by an approximation H h (u h ) in V h of the Hessian of u h . This is done with a Zienkiewicz-Zhu-type ( [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF]) recovery method defined in the seventh chapter of [START_REF] Magoules | Computational Fluid Dynamics[END_REF]. Then, the evaluation of Π h u -u is built on the edge e ij as a quadratic function vanishing at both extremities of e ij and of second-derivative in direction e ij defined from H h (u h ). We replace Π h u -u by π h u h -u h where π h u h -u h is defined on edge ij as follows:

∀ x ∈ e ij , (π h u h -u h )(x) = 1 2 (H h (u h )(x i ) + H h (u h )(x j ))(x -x i )(x -x j )
which allows a mid-edge integration on each triangle of Ω h . This replacement is justified in [START_REF] Yu | Hessian-based anisotropic mesh adaptation in domains with discrete boundaries[END_REF]. We shall see, in the sequel, that it is useful to apply a similar estimate for the f -

f h term, f -f h ≈ -(π h f h -f h ) .
In order to approximate the implicit error Π h u -u h , solution of (5), we define our a priori implicit corrector by:

u prio ∈ V h , and ∀ ϕ h ∈ V h , a(ū prio , ϕ h ) = ∂T ij ( 1 ρ ∇ϕ h | T i -∇ϕ h | T j ) • n ij ∂T ij (π h u h -u h ) dσ -(ϕ h , π h f h -f h ). ( 6
)
This corrector is an approximation of Π h u -u h . Replacing again the unknown interpolation error by its evaluation on the discrete approximation, we define our a priori corrector by:

u prio = ū prio -(π h u h -u h ). (7)
Assuming the approximations made between ( 5) and ( 6) are small, u prio is built in such a way that:

u prio ≈ u -u h .
This corrector is easy to compute but of a priori low accuracy.

Finer-grid defect correction corrector for the PDE solution

A second option for a corrector relies on using a fictitious finer grid. Let us assume that the approximation is in its asymptotic mesh convergence phase for the mesh Ω h under study, of size h. Then, this will be also true for a strictly two-times finer embedding mesh Ω h/2 and, for our secondorder accurate scheme applied to a sufficiently smooth problem, the heuristics of Richardson analysis writes:

u h = A -1 h f h , u h/2 = A -1 h/2 f h/2 ⇒ u -u h/2 ≈ 1 4 (u -u h ) (8)
where u h and u h/2 are respectively the solutions on Ω h and Ω h/2 . We have formally:

Π h u -Π h u h/2 ≈ 1 4 (Π h u -u h )
which implies:

Π h u -u h ≈ 4 3 (Π h u h/2 -u h ). Now, A h/2 (u h/2 -P h→h/2 u h ) = f h/2 -A h/2 P h→h/2 u h .
where P h→h/2 linearly interpolates coarse values on fine mesh. This allows to evaluate Π h u -u h but needs to solve a finer grid system. Let us introduce the residual transfer R h/2→h which accumulates on coarser grid vertices the values at finer vertices in neighboring coarse elements multiplied with barycentric weights. In order to reduce the computational cost to solving a coarser grid system, we approximate

Π h A -1 h/2 by A -1 h R h/2→h : Π h (u h/2 -P h→h/2 u h ) ≈ A -1 h R h/2→h (f h/2 -A h/2 P h→h/2 u h ).
This motivates the definition of a finer-grid Defect-Correction (DC) corrector as follows:

A h ū DC = 4 3 R h/2→h (f h/2 -A h/2 P h→h/2 u h ) (9)
In the case of local singularities, statement (8) is not true for uniform meshes but we have some hints that it holds almost everywhere for a sequence of adapted meshes, according to [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF]. The DC corrector ū DC approximates Π h u -u h instead of u -u h and can be corrected as the previous one:

u DC = ū DC -(π h u h -u h ). ( 10
)

Interpolation error optimization

Main notions and notations are stated by recalling the basics of the Hessian-based approach.

Mesh parametrization

We shall work inside the framework proposed in [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF][START_REF] Loseille | Continuous mesh framework. Part II: validations and applications[END_REF]. The main idea of this framework is to model discrete meshes by continuous Riemannian metric fields. It allows us to define the adaptation problem as a differentiable optimization problem , i.e., to apply, on the class continuous metrics, a calculus of variations which cannot be applied on the class of discrete meshes. This framework lies in the class of metric-based methods. A continuous metric M of the computational domain Ω is a Riemannian metric field [START_REF] Berger | A panoramic view of Riemannian geometry[END_REF] M = (M(x)) x∈Ω where M(x) is a symmetric 2 × 2 matrix. We define the total number of vertices of M as:

C(M) = Ω det(M(x)) dx.
Given a continuous metric M, we shall say that a discrete mesh H of the same domain Ω is a unit mesh with respect to M if each triangle K ∈ H, defined by its list of edges (a i b i ) i=1...3 , verifies:

∀i ∈ [1, 3], 1 0 t a i b i M(a i + t a i b i ) a i b i dt ∈ 1 √ 2 , √ 2 .
The rest of the paper will try to find the best metric M from an error analysis which is asymptotic with respect to mesh size.

Interpolation-based optimal metric

Let u be any sufficiently smooth function defined on Ω. Let M be a mesh/metric of Ω. In the rest of the paper, M replaces h as the discretization index. We consider only meshes M involving enough nodes for justifying the replacement of the complete error by its main asymptotic part. The P 1 interpolation error |Π M u -u| can be approximated in terms of second derivatives of u and of the metric M by the continuous interpolation error defined in [START_REF] Loseille | Continuous mesh framework. Part I: well-posed continuous interpolation error[END_REF]:

|Π M u -u| ≈ |u -π M u| with: |u -π M u|(x) = 1 8 trace(M -1 2 (x) |H u (x)| M -1 2 (x)) (11)
where |H u | is deduced from H u by taking the absolute values of its eigenvalues. Starting from:

u -π M u L p (Ω h ) = Ω trace M -1 2 (x)|H u (x)|M -1 2 (x) p dx 1 p , (12)
we define as optimal metric the one which minimizes the right hand side under the constraint of a total number of vertices equal to a parameter N . After solving analytically this optimization problem (see e.g. [START_REF] Yu | An unified asymptotical analysis of interpolation errors for optimal meshes[END_REF][START_REF] Huang | Metric tensors for anisotropic mesh generation[END_REF][START_REF] Chen | Optimal anisotropic meshes for minimizing interpolation errors in L p -norm[END_REF][START_REF] Loseille | Adaptation de maillage anisotrope 3D multi-échelles et ciblée à une fonctionnelle pour la mécanique des fluides. Application à la prédiction haute-fidélité du bang sonique[END_REF][START_REF] Alauzet | Multi-model and multi-scale optimization strategies. application to sonic boom reduction[END_REF], we get the unique optimal (M L p (x)) x∈Ω as:

M L p = K p (H u ) with K p (H u ) = D L p (det |H u |) -1 2p+2 |H u | and D L p = N Ω (det |H u | dx) p 2p+2 -1 , (13)
where D L p is a global normalization term set to obtain a continuous metric with complexity N and (det

|H u |) -1
2p+2 is a local normalization term accounting for the sensitivity of the L p norm. In the case of an adaptation loop for solving a Partial Differential Equation, a continuous function u is not available but an approximate solution u M is. In that case, the continuous interpolation error [START_REF] Berger | A panoramic view of Riemannian geometry[END_REF] is replaced by:

|u M -π M u M |(x) = 1 8 trace(M -1 2 (x) |H u M (x)| M -1 2 (x)) (14)
where H u M is an approximate Hessian evaluated with a recovery method.

According to the continuous metric framework, statement [START_REF] Castro-Díaz | Anisotropic unstructured mesh adaptation for flow simulations[END_REF] defines directly a continuous optimal metric. In practice, solving (13) is done by approximation, i.e. in a discrete context with a couple (mesh, solution) denoted (H M , u M ) and iteratively obtained through the following fixed point:

Step 1: compute the discrete state u M on mesh H M , Step 2: compute sensor s M = s(u M ) and optimal metric M opt inter = K p (H M (s M )), Step 3: put M = M opt inter , generate a new mesh H M = H M opt inter and go to 1, until convergence. In the above algorithm, the continuous Hessian of s is replaced by an approximate Hessian H M (s M ), evaluated by the patch-recovery approximation defined in [START_REF] Magoules | Computational Fluid Dynamics[END_REF]. In our Hessian-based numerical examples, the L 2 case, p = 2, has been considered. The above notation K p will also be used in the next sections for p = 1.

Implicit a priori error estimate

In contrast to a corrector as defined in Section 2, an asymptotic upper bound of the approximation error should allow an easier error reduction by minimisation of its norm with respect to the metric. In our PDE discretisation notations, we henceforward replace in the discretization index h by the index M which holds for any unit mesh of the metric M. The implicit a priori error system (5) then writes:

∀ ϕ ∈ V M , a(Π M u -u M , ϕ) = ∂T ij 1 ρ (∇ϕ| T i -∇ϕ| T j ) • n ij ∂T ij (Π M u -u) dσ -(ϕ, Π M f -f ) (15)
where T ij are the triangles of a unit mesh for M and the proposed corrector is expressed with the discrete solution:

ū prio ∈ V M , a(ū prio , ϕ) = K(M, ϕ, u M ) ∀ ϕ ∈ V M with K(M, ϕ, u M ) = ∂T ij 1 ρ (∇ϕ| T i -∇ϕ| T j ) • n ij ∂T ij (π M u M -u M ) dσ -(ϕ, π M f -f ), u prio = ū prio -(π M u M -u M ). ( 16
)
We now restart from [START_REF] Alauzet | Multi-model and multi-scale optimization strategies. application to sonic boom reduction[END_REF]. The following result is proven in [START_REF] Belme | Aérodynamique instationnaire et méthode adjointe[END_REF]:

Lemma 4.1. We assume that the metric anisotropy is bounded by a positive number. For any smooth couple of functions (u, ϕ), where u is not necessarily a solution of (1), we have the following bound:

| Ω 1 ρ ∇(u -Π M u)∇Π M ϕdx| K Ω 1 ρ ρ S (H(ϕ)) |u -Π M u|dx ( 17
)
where A B holds for a majoration asymptotically valid, i.e. A ≤ B + o(A). Expression ρ S (H(ϕ)) holds for the largest (in absolute value) eigenvalue of the Hessian H(ϕ) of ϕ.

The next section shows how to use this estimate.

Equation-based adaptation

Scalar output "goal-oriented" analysis

The goal-oriented analysis relies on the minimization of the error δj goal (M) done in the evaluation of the scalar output j = (g, u) , error which we simplify as follows:

δj goal (M) = |(g, u -u M )| = |(g, Π M u -u M + u -Π M u)|. ( 18
)
Let us define the discrete adjoint state u * g,M :

∀ψ M ∈ V M , a(ψ M , u * g,M ) = (ψ M , g). (19)
Then:

(g, Π M u -u M + u -Π M u) = a(Π M u -u M , u * g,M ) + (g, u -Π M u) and, using (3), (g, Π M u -u M + u -Π M u) = a(Π M u -u, u * g,M ) + (f -Π M f, u * g,M ) + (g, u -Π M u) thus δj goal (M) ≈ |a(Π M u -u, u * g,M ) + (f -Π M f, u * g,M ) + (g, u -Π M u)|
Recall that u is unknown. The second term Π M u -u, similar to the main term of the Hessianbased adaptation in Section 3.2, can be explicitly approached in the same way, i.e. introducing the continuous interpolation error ( 14):

δj goal (M) |a(Π M u -u, u * g,M )| + |(f -Π M f, u * g,M )| + |g||π M u M -u M |
For the second Π M u -u term, we first apply Lemma 4.1 and then also introduce the continuous interpolation error. We get

δj goal (M) Ω 1 ρ ρ S (H(u * g,M )) + |g| |π M u M -u M | + |u * g,M | |π M f -f | dx.
It is then reasonable to try to minimize the RHS of this inequality instead of the LHS. But this still involves some difficulty due to the dependancy of adjoint state u * g,M with respect to M. We shall further simplify our functional by freezing, during a part of the algorithm, the adjoint state. The idea is that, when we change the parameter M, the discrete adjoint u * g,M is close to its (non-zero) continuous limit and is thus not much affected, in contrast to the interpolation errors |π M u M -u M | and |π M f -f |. We then consider, for a given M 0 , the following optimum problem:

min M Ω 1 ρ ρ S (H(u * g,M 0 )) + |g| |π M u M -u M | + |u * g,M 0 | |π M f -f | dx.
This will produce an optimum:

M opt,M 0 = arg min M |tr(M -1/2 1 ρ ρ S H(u * g,M 0 ) + |g| |H u | + |u * g,M 0 ||H f | M -1/2 )|.
Observing that, in the integrand,

H goal,0 = [ 1 ρ ρ S (H(u * g,M 0 )) + |g|] |H u | + |u * g,M 0 | |H f |
is a positive symmetric matrix, we can apply the above calculus of variation and get:

M opt,M 0 = K 1 ( [ 1 ρ ρ S (H(u * g,M 0 )) + |g|] |H u | + |u * g,M 0 | |H f |)
where K 1 is defined in [START_REF] Castro-Díaz | Anisotropic unstructured mesh adaptation for flow simulations[END_REF]. This solution can then be introduced in a fixed-point loop and will produce the solution of:

M opt,goal = K 1 ( [ 1 ρ ρ S (H(u * g,M opt,goal )) + |g|] |H u | + |u * g,M opt,goal | |H f |).
Let us precise how the discrete algorithm is organised:

Step 1: compute the discrete state u M on mesh H M ,

Step 2: compute the discrete adjoint state u * M , Step 3: compute optimal metric M opt goal (u M ), Step 4: put M = M opt goal (u M ), generate a new mesh H M = H M opt goal (u M ) and go to 1, until convergence.

The adaptation of this process to the Euler model of Gas Dynamics is studied in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] for the steady case and in [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] for the unsteady case.

Norm-based functional

We are now interested by the minimization of the following expression with respect to the mesh M:

δj(M) = ||u -u M || 2 L 2 (Ω) . (20)
Introducing g M = u -u M , we get a formulation similar to the goal-oriented formulation:

δj(M) = (g M , u -u M ). (21)
Let us define the discrete adjoint state u * M :

∀ψ M ∈ V M , a(ψ M , u * M ) = (ψ M , g M ). ( 22
)
Then, similarly to Section 5.1, we have to solve the following optimum problem.

min

M Ω 1 ρ ρ S (H(u * M )) + |g M | |π M u M -u M | + |u * M | |π M f -f | dx.
Exactly as for Section 5.1, we freeze the dependancy of the adjoint state.

min

M Ω 1 ρ ρ S (H(u * M 0 )) + |g M | |π M u M -u M | + |u * M 0 | |π M f -f | dx. M opt,M 0 = K 1 ( [ 1 ρ ρ S (H(u * M 0 )) + |g M |] |H u | + |u * M 0 | |H f |).
In practice, the g M RHS in [START_REF] Hartman | Multitarget error estimation and adaptivity in aerodynamic flow simulations[END_REF][START_REF] Huang | Metric tensors for anisotropic mesh generation[END_REF] is replaced by the corrector u prio,M . In order to get the final norm-oriented optimum M opt,norm , we apply:

Step 1: compute the discrete state u M (α) on mesh H M (α) ,

Step 2: solve the linearised corrector system:

a(ū prio,M (α) , ϕ) = ∂T ij (∇ϕ| T i -∇ϕ| T j ) • n ij ∂T ij (π M (α) u M (α) -u M (α) ) dσ -(ϕ, π M (α) f M (α) -f M (α) ). ( 23
)
where

π M (α) u M (α) -u M (α)
is expressed in terms of metric and Hessian, as in [START_REF] Arian | Admitting the inadmissible: Adjoint formulation for incomplete cost functionals in aerodynamic optimization[END_REF].

Put u prio,M (α) = ū prio,M (α) -(π M (α) u M (α) -u M (α) ).
Step 3: then, solve the adjoint system:

a(ψ, u * prio,M (α) ) = (u prio,M (α) , ψ) (24)
Step 4: put:

M (α+1) = K 1 ([|u prio | + 1 ρ ρ S H(u * prio )] |H u M (α) | + |u * prio ||H f |) (25)
Step 5: generate a new mesh H M (α) and go to 1, until convergence to a fixed point M opt,norm = M (∞) .

The new algorithm involves the solution of two extra linear systems with respect to a basic Hessian-based algorithm and one extra linear system with respect to the goal-oriented algorithm.

Remark: In contrast to the adjoint of the goal-oriented algorithm, these auxiliary variables are not consistent with a continuous adjoint (and do not converge towards it when mesh is refined). They are correctors and converge to zero. Then it seems that we cannot get a good initial condition to an FMG phase (from a coarse mesh to a finer one) on these variables. For doing it anyway, we rely again on the Richardson heuristics. We assume that our data and solution are sufficiently regular for applying the Aubin-Nitsche L 2 analysis. This implies the second-order L 2 convergence of approximate solution to exact one, i.e. ||u M -u|| L 2 = O(h 2 ). This tends to indicate that our corrector norm would be also O(h 2 ), i.e., formally:

||u prio,M || L 2 = O(h 2 ).
We write this in terms of the number N of vertices for metric M as follows (again formally):

||u prio,M || L 2 = O(N -1 ).
Further, following a (formal) Richardson expansion:

u prio,M ≈ N -1 u prio (26)
where u prio is the limit (assumed to exist) of N u prio,M when mesh is refined. Similarly, we can put:

u * prio,M ≈ N -1 u * prio . (27)
In particular, when passing from a mesh M 1 with N 1 vertices, M 2 with N 2 vertices, approximations of u prio,M 2 and u * prio,M 2 can be defined as

u prio,M 2 ≈ N 1 N 2 u prio,M 1 ; u * prio,M 2 ≈ N 1 N 2 u * prio,M 1 (28)
and used as initial conditions inside the FMG method.

Numerical examples

We restrict our study to a benchmark of two-dimensional Poisson problems. We conjecture that the two following mesh adaptation methods produce L 2 convergent solutions to continuous. The first method, the Hessian-based method (with p = 2), is just heuristically relying on usual finite-element estimates. The second method, our novel norm-oriented method, is directly built on the minimisation of the L 2 error norm. We do not consider goal-oriented applications for which examples of computations can be found in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] and [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF]. As already remarked, the convergence of goal-oriented solutions to continuous is definitively questionable.

Numerical features

In [START_REF] Brèthes | A mesh-adaptative metric-based Full-Multigrid for the Poisson problem[END_REF], a mesh-adaptative full-multigrid (FMG) algorithm relying on the Hessian-based adaptation criterion is designed. We first describe in short this algorithm for the Hessian-based option. A sequence of numbers N k of vertices is specified from a coarse mesh to finer one N 0 = N, N 1 = 4N, N 2 = 16N, N 3 = 64N, .... For each mesh size N k , a sequence of adapted meshes of size N k is built by iterating the following loop: (a) computing a solution, (b) computing the optimal metric, (c) building the adapted mesh. In (a), a multi-grid V-cycle is applied to a sufficient convergence. In (b), approximations of the Hessians are performed as in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. When changing of mesh, an interpolation is applied in order to enjoy a good initial condition. A prescribed number of 4 adaptation iterations is applied at each mesh fineness N k .

The extension of the above loop to norm-oriented adaptation consists of replacing the single Hessian evaluation by: -the computation of the corrector, using MG and, as initial solution, the previous evaluation interpolated to current mesh and corrected according to [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF], -the computation of the adjoint, using MG and, as initial solution, the previous evaluation interpolated to current mesh and corrected according to [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF], -the evaluation of (25). Let us discuss computer efficiency. In the demonstrator of [START_REF] Brèthes | A mesh-adaptative metric-based Full-Multigrid for the Poisson problem[END_REF], a particular feature is the stopping criterion of FMG which applies to the convergence of the solution of the unique solved system, i.e.

the system under study, u = A -1 f . It is then possible to enjoy a better and better initial condition and control the iterative and approximation errors convergence. Consequently, it was possible, in [START_REF] Brèthes | A mesh-adaptative metric-based Full-Multigrid for the Poisson problem[END_REF], to show that mesh adaptation carries large improvement not only in terms of accuracy for a given number of vertices but also in terms of accuracy for a given computational time.

In contrast, the method proposed in this paper involves three systems to solve: (1) the system under study, u = A -1 f , (2) the corrector system, (3) the adjoint system. We shall give an idea of the performances for one test case.

2D Boundary layer

This test case is taken from [START_REF] Formaggia | Anisotropic a priori error estimates for elliptic problems[END_REF]. We solve the Poisson problem -∆u = f in ]0, 1[×]0, 1[ with Dirichlet boundary conditions and a right-hand side f chosen for having:

u(x, y) = [1 -e -αx -(1 -e -α )x]4y(1 -y).
The coefficient α is chosen equal to 100. The graph of the solution is depicted in Figure 1. Before applying our mesh adaptative algorithm, it is interesting to evaluate the accuracy of our correctors. We choose a 161 × 161 uniform mesh and show, in Figure 2 and Figure 3, the cut of u -u h compared with the cut of u . We observe that the a priori corrector does its job in a correct but inaccurate manner while the DC one is rather accurate. We have also observed that the DC corrector does not consume notably more CPU than the a priori one. Therefore, we keep this option for the rest of the test case. In Figure 4, we show a set of FMG calculations for the considered test case. The numbers of vertices of the successive meshes are supported by the horizontal axis, from 120 vertices to 30,000 vertices. The vertical axis gives the L 2 -norm of the approximation error |u -u h | L 2 obtained on the mesh. Its variation with respect to number of vertices is compared in Figure 4 for the three following algorithms: (a) the uniform-mesh FMG, (b) the Hessian-based adaptative FMG and (c) the norm-oriented adaptative FMG. We observe that both adaptation methods carry an important improvement with respect to uniform-grid FMG (25921 vertices on finest mesh). For essentially the same number of vertices (32318), the Hessian option gives an error divided by 47. The norm-oriented option appears as better with an error divided by 208 with 29485 vertices. Since the exact solution u is analytically available, we also propose an a posteriori measure of the correctors efficiency by comparing the convergence of our norm-oriented adaptation (equipped with either corrector) with the same method except that the corrector is replaced by u -u h . Of course, that latter algorithm is not a mesh adaptation method since we assume that we already know the exact solution. In Figure 5, we observe that the error convergence for the three computations are very close to each other. This confirms the interest of the two proposed correctors.

Lastly, when using the proposed method in usual conditions, i.e. without knowing the analytic solution, it can be interesting to have an evaluation of the final error u-u h by using the DC corrector which appeared in experiments as more accurate. Figure 6 shows a comparison of the evaluated error computed from the corrector with the actual error. We observe that this evaluation is useful, giving the good order of magnitude. For the finer mesh of 30, 000 vertices, the evaluated error is 2.14 10 -5 while the actual error is 4.31 10 -5 . However, for most meshes, the evaluation remains somewhat optimistic (smaller than the actual error).

Bubble-like test case with thick interface

We are interested by a Poisson problem the solution of which is a function u equal to 1 on a disk and to 0 in the rest of the domain. This function is the prototype of the pressure in a multi-fluid flow involving capillary forces. The source term is a Dirac derivative. We smooth this computation by defining a thickness ε for defining an annular region separating the two subdomains (outside the disk, inside the disk) and in which u is smoothly varying from 0 to 1: if (x, y) is located inside the annular region, u(x, y) is given by the formula: u(x, y) = 1 2 + 1 2 sin( πψ ε ) with ψ = 0.25 -(x C -x) 2 + (y C -y) 2 . From this solution, a right-hand side f is computed. Given a mesh, vertex values of f h are interlopated from the analytic f . As a result, for rather coarse meshes, the zone where f is not zero can be spuriously missed and f h can be zero even in the neighborhood of the high values of f . We consider first a quite large thickness of ε = 0.1. An approximate solution u h is shown in Figure 7. As for the previous test case, we first evaluate the accuracy of the corrector. We choose a uniform mesh 161 × 161 and show, in Figure 8 and Figure 9, the cut of u -u h compared with the cut after correction, that is u-u h -u h . We observe that both a priori and Defect-Correction correctors do an accurate job.

The three methods are again compared in Figure 10: standard FMG, Hessian-based adaptative FMG and norm-oriented adaptative FMG. For the Hessian-based calculation, we observe a tendancy for a slower convergence for finer meshes, finishing with an error which is worse than the uniform refinement. The proposed norm-oriented adaptative method behaves in a better way with a five times smaller error than for the uniform refinement.

Bubble-like test case with thin interface

In order to evaluate the robustness of the methods with respect to steeper gradients, we consider the same test case with a thinner transition: ε = 0.02. The convergence of the three methods is shown in Figure 14. Due to the very thin support of the right-hand side f , the three methods start with a zero f h . Then, either with adaptation or refinement, the error increases to several units. The convergence of the uniform FMG shows an acceptable slope but the error values remain relatively huge. By comparison with the thick-bubble convergence, we may infer that the slope of uniform FMG will be second-order with even higher number of vertices and a 0.1 % error may be not attained for meshes of less than 10 millions nodes. The Hessian-based final result is a little better but globally disappointing. The norm-oriented convergence starts chaoticly before being monotone and second-order for meshes finer than 2000 vertices. The final L 2 errors produced by the three methods are 0.17589 with 25921 vertices for the uniform FMG, 0.03773 with 32127 vertices for the Hessian-based adaptation and 0.000585 for 29742 vertices for the norm-oriented calculation, 300 times smaller than the first result. In order to show the behavior of the algorithms for more complex solutions, we have computed the case of a solution involving three bubbles of different sizes. A comparison, in Figure 15, of the solutions with and without adaptation demonstrate that the uniform FMG (results of top range) is unable to produce even a rough approximation with 6561 vertices. Only the solution for 25921 vertices is rather good, with an error L 1 norm of 0.50 but with mean values in bubbles showing still 10% deviations. In contrast, the norm-oriented mesh-adapted solutions (bottom range of Figure 15) show a fast convergence. For 1681 vertices, the solution has the same quality as the uniform 25921-vertices one (error is slightly larger, 0.54). For 6561 and 25921 vertices, the accuracy seems good. The apparent convergence order is greater than two. In contrast, we observe that the Hessian-based calculation is rather disappointing.

Poisson problem with discontinuous coefficient

This test case exemplifies the singularity which is met in the simulation of multi-fluid flows with a large deviation between the densities ρ 1 and ρ 2 of each phase. In the case where a projection algorithm is applied to solve the Navier-Stokes equations for incompressible flow, a Poisson problem with discontinuous coefficients has to be solved. An example can be found in [START_REF] Guégan | An L ∞ -L p mesh adaptive method for computing unsteady bi-fluid flows[END_REF]. The present case does not satisfy the smoothness assumptions introduced for deriving our method. However, a usual expectation in mesh adaptation is that the methods should also apply well on non-smooth contexts. We consider the equation of Poisson -div( 1ρ ∇u) = rhs with a discontinuous coefficient taking two different values 1/ρ 1 and 1/ρ 2 on two sub-domains Ω 1 and Ω 2 separated by an interface which is a sufficiently smooth curve for having a normal vector. This PDE is mathematically referred as a transmission problem and the solution is continuous across the interface but of discontinuous normal derivatives since:

1/ρ 1 ∇u 1 • n = 1/ρ 2 ∇u 2 • n
where u 1 and u 2 are the restrictions of the solution u on Ω 1 and Ω 2 . In our example, we define them as follows

u| Ω i = u i = α i + β i (x 2 + y 2 ) i = 1, 2.
Further, Ω 2 is the disk of center (0.5, 0.5) and of radius 0.2 in the computational domain ]0, 1[×]0, 1[ and we have:

1/ρ 1 = 1000. ; α 1 = 1.23579... ; β 1 = -2.47158... 1/ρ 2 = 1. ; α 2 = 100. ; β 2 = -2471.58... (29)
This is sketched in Figure 17. In the discrete model, the interface appears only as values of 1/ρ evaluated on the vertices of each grid. From the examination of cuts of the correctors in Figures 18 and19, we observe that their pointwise accuracy is somewhat a disaster. Our analysis used intensively smoothness of functions and assumes, for the Defect Correction option that second-order convergence applies, which is not true for this discontinuous case. However, the injection of the correctors in the norm-oriented functional and adjoint performs adequately. The convergence of our process is not really affected by the poor accuracy of the correctors, as can be seen from the comparison of convergence with the correctors or with the exact error, Figure 21. The anisotropy of the finer mesh is illustrated in Figure 20. A convergence in terms of number of vertices with non-adaptative and Hessian-based adaptive is given in Figure 22. We observe that, without mesh adaptation (crosses +), the convergence order is around 1. This behavior can be explained by the singularity of the solution. In contrast, the overall convergence order of the adaptative process is about two. Note that this convergence will finally deteriorate when we attain the limits of the stretching capabilities of the mesh generator. However, the second-order numerical convergence observed is a usual bonus obtained by anisotropic mesh adaptation which has already been noted in [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF]. In [START_REF] Brèthes | A mesh-adaptative metric-based Full-Multigrid for the Poisson problem[END_REF], the convergence of an anisotropic adaptation has been compared with its isotropic analogous for the same test case. The anisotropic calculation was converging at order two while the isotropic one was converging at order 3/2. A short analysis of these behaviors is proposed in [START_REF] Courty | Continuous metrics and mesh adaptation[END_REF].

Let us now compare with the Hessian-based method. This is the only one of our test cases for which the convergence of error in terms of number of nodes of the norm-oriented formulation is not neatly better (but it is as good) than the analog convergence of the Hessian-based formulation. The error plays a role of weighting in the new formulation. Then a possible explanation of this -relativelydisappointing behavior is that, with the realistic ratio of densities we have chosen, the error is very small close to the discontinuity.

Although not optimised for CPU, the novel algorithm may give a good increase of efficiency. We present, in Figure 23, the behavior of the error as a function of CPU. Three curves are compared. The upper curve corresponds to the non-adaptative FMG calculation. An error of about 0.06 is obtained on 100000 vertices after 25000 seconds of CPU. We present also the Hessian-based adaptative convergence obtained in [START_REF] Brèthes | A mesh-adaptative metric-based Full-Multigrid for the Poisson problem[END_REF] in which the convergence of MG is adequately controlled by a stopping criterion. Despite the remeshing and re-computing for adaptation, the total time is about twice the non-adaptative FMG time, for an error of 0.001 -0.002, at least 30 times smaller. The third curve shows the CPU behavior of the norm-oriented algorithm. The total time is about ten times the non-adaptative FMG time, for an error of 0.001. After 25000 seconds of CPU, the error is about 10 times smaller than with the non-adaptative FMG at same time.

A 1D Boundary layer

The new method has shown a good behavior for all test cases we tried except one which we describe now. It is a boundary layer case with a 1D solution: u(x, y) = u(x) of a Dirichlet-Neumann We check first the correctors. Both seem adequate on a uniform grid, as shown with a horizontal cut depicted on Figure 24 and Figure 25.

In Figure 26, the non-adaptative FMG produces an approximation error of 0.003 (30000 vertices). This convergence is relatively satisfactory, being a second order convergence. However, in order to reach a 10 -7 error level, several hundred millions vertices will be necessary with this sequence of uniform meshes. A second curve is obtained with the adaptative FMG with the Hessian-based criterion. Final convergence is disappointing since the slope is first-order. The same problem appears with our new algorithm. A deeper examination of adaptation criteria has shown that the high derivatives of the right-hand side f are very close to boundary x = 1. It could not be seen by the algorithm, because of the weighting by the adjoint u * , which is zero at this boundary. We have replaced the norm-oriented optimum metric by its intersection with the metric based on the Hessian of f . Then, the convergence improved a lot. In contrast, introducing the same metric intersection in the other cases did not produce second-order convergence.

Conclusion

The norm-oriented mesh adaptation method is an answer to a well formulated problem. Considering a numerical scheme, here the most used FEM, and prescribing a norm, we want to find the mesh giving the smallest approximation error in that norm for a given number of vertices. The norm-oriented mesh adaptation method transforms the problem into an optimization problem which is mathematically well-posed. It relies on the following other features. A corrector represents the approximation error. We give two examples of correctors. An a priori corrector is built from the variational discrete statement. A Defect-Correction corrector is built from a finer-mesh defect correction principle. These correctors appear as not very accurate but sufficiently accurate for our purpose. According to the type of approximation, at least the second one, Defect-Correction is extendable to many models and schemes.

The norm-oriented algorithm is presented as a natural extension of the goal-oriented algorithm which, in our formulation, is itself a natural extension of the Hessian-based algorithm. More precisely, while the Hessian-based algorithm solves only the PDE under study in the mesh-adaptation loop, the goal-oriented algorithm also solves an adjoint system (with linearised operator, transposed). The norm-oriented algorithm solves three systems, a corrector (linearised system with an adhoc RHS), an adjoint (linearised and transposed with the corrector as RHS) and the PDE itself. The three algorithms have in common an anisotropic a priori error analysis and a metric-based mesh parameterisation.

The Hessian-based method produces convergent solution fields but does not take into account the precise equation and discretization. The goal-oriented method takes into account equation and discretization but is too focused on a particular output and does not generally produce convergent solution fields. The norm-oriented method has the advantages of both.

In order to show the improvement obtained with respect to previous methods, we compare, in our experiments, the two field-convergent options, i.e. the Hessian-based method and the norm-oriented method. Our benchmark examines convergence to continuous thanks to the application of a Full-Multigrid (FMG) process. Approximation errors can then be compared as functions of the number of degrees of freedom. Although of larger complexity, the mesh adaptative algorithms are generally more efficient in terms of CPU time than the non-adaptative one. The advantage is renforced by the fact that the same level of error is reached with a much lower number of vertices. Secondly, although the elliptic context is known as favourable to Hessian-based methods, the norm-oriented approach behaves notably better than Hessian-based in terms of approximation error for a given number of vertices. This is observable in particular for singular or stiff contexts. We have presented an example of comparison of computational effort but our algorithm is not optimised. It is, therefore, in several cases, less efficient than Hessian-based adaptative FMG.

The method is rather general and it is now applied to more complex PDE models from CFD (Euler, Navier-Stokes), see [START_REF] Loseille | Anisotropic norm-oriented mesh adaptation for compressible flows[END_REF], involving dominant advection effects for which a Hessian-based approach is much less efficient.

Figure 1 :

 1 Figure 1: Fully 2D Boundary layer test case : sketch of the solution.

Figure 2 :

 2 Figure 2: Fully 2D Boundary layer test case : comparison of error cuts for y = 0.5: plus signs (+) depict the approximation error u -u h and crosses (×) depict the a priori corrector u prio . The corrector is able to correct about 60% of the approximation error.

Figure 3 :

 3 Figure 3: Fully 2D Boundary layer test case : comparison of error cuts for y = 0.5: plus signs (+) depict the approximation error u -u h and crosses (×) depict the Defect-Correction corrector u DC . The corrector is able to correct about 95% of the approximation error.

Figure 4 :

 4 Figure 4: Fully 2D Boundary layer test case: convergence of the error norm |u -u h | L 2 as a function of number of vertices in the mesh for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG, ( * ) norm-oriented adaptative FMG.

Figure 5 :

 5 Figure 5: Fully 2D Boundary layer test case: convergence based on (×) the a priori corrector or on the ( * ) Defect-Correction one, compared with ( ) a virtual adaptation controlled by u -u h .

Figure 6 :

 6 Figure 6: Fully 2D Boundary layer test case: comparison of the predicted error norm proposed by the method, curve with ( ), with the exact error norm (×).

Figure 7 :

 7 Figure 7: Circular layer test case: an adapted mesh and the corresponding numerical solution u h . Palette from -6.4 10 -5 to 0.9996.

Figure 8 :

 8 Figure 8: Thick bubble case: comparison of (+) error u -u h and (×) a priori corrector.

Figure 9 :

 9 Figure 9: Thick bubble case: comparison of (+) error u -u h and (×) Defect-Correction corrector.

Figure 10 :

 10 Figure 10: Thick bubble test case: convergence of the error norm |u -u h | L 2 as a function of number of vertices in the mesh for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG and ( * ) norm-oriented adaptative FMG.

Figure 11 :

 11 Figure 11: Thin bubble case: an adapted mesh and the corresponding numerical solution u h . Palette from 0.000 to 1.0593.

Figure 12 :

 12 Figure 12: Thin bubble case: comparison of (+) error u -u h and (×) a priori corrector.

Figure 13 :

 13 Figure 13: Thin bubble case: comparison of (+) error u -u h and (×) Defect-Correction corrector.

Figure 14 :

 14 Figure 14: Thin bubble test case: convergence of the error norm |u -u h | L 2 as a function of number of vertices in the mesh for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG and ( * ) norm-oriented adaptative FMG.

Figure 15 :

 15 Figure 15: Multiple bubble test case: top: non-adaptative solutions for 441, 1681,6561, 25921 nodes, and, bottom, the norm-oriented mesh adaptative solution for about the same number of nodes (palette retricted to values between the minimum (0.) and the maximum (1.) of the exact solution.

Figure 16 :

 16 Figure 16: Multiple bubble test case: convergence of the error norm |u -u h | L 2 as a function of number of vertices in the mesh for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG and ( * ) norm-oriented adaptative FMG.

Figure 17 :

 17 Figure 17: Poisson problem with discontinuous coefficient: sketch of exact solution definition and a typical computation of it.

Figure 18 :

 18 Figure 18: Poisson problem with discontinuous coefficient: comparison of (+) error u -u h and (×) a priori corrector.

Figure 19 :

 19 Figure 19: Poisson problem with discontinuous coefficient: comparison of (+) error u -u h and (×) Defect-Correction corrector.

Figure 20 :

 20 Figure 20: Poisson problem with discontinuous coefficient: views of final mesh: global view, zoom on right part, zoom of the zoom.

Figure 21 :

 21 Figure 21: Poisson problem with discontinuous coefficient: convergence based on (×) the a priori corrector or on the ( * ) Defect-Correction one, compared with ( ) a virtual adaptation controlled by u -u h .

Figure 22 :

 22 Figure 22: Poisson problem with discontinuous coefficient: convergence of the error norm |u -u h | L 2 as a function of number of vertices in the mesh for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG and ( * ) norm-oriented adaptative FMG.

Figure 23 :

 23 Figure 23: Poisson problem with discontinuous coefficient: efficiency analysis : CPU time for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG and ( * ) norm-oriented adaptative FMG.

Figure 24 :

 24 Figure 24: 1D Boundary Layer: comparison of error cuts for y = 0.5: plus signs (+) depict the approximation error u -u h and crosses (×) depict the a priori corrector u prio . The corrector is able to correct about 55% of the approximation error.

Figure 25 :

 25 Figure 25: 1D Boundary Layer: comparison of error cuts for y = 0.5: plus signs (+) depict the approximation error u -u h and crosses (×) depict the Defect-Correction corrector u DC . The corrector is able to correct about 80% of the approximation error.

Figure 26 :

 26 Figure 26: 1D Boundary Layer: convergence of the error norm |u -u h | L 2 as a function of number of vertices in the mesh for (+) non-adaptative FMG, (×) Hessian-based adaptative FMG, ( * ) norm-oriented adaptative FMG, and (o) norm-oriented adaptative FMG with intersection with RHS Hessian.

Acknowledgements

This work was partly done in the MAIDESC ANR project which is supported by the french ministery of Research under contract ANR-13-MONU-0010. The fellowship of Gautier Brèthes is supported by Lemma and région Provence-Alpes-Côte d'Azur.


Emmanuelle Itam
email: emmanuelle.itam@umontpellier.fr

Stephen Wornom
email: stephen.wornom@inria.fr

Bruno Koobus
email: koobus@math.univ-montp2.fr

Alain Dervieux
email: alain.dervieux@inria.fr

Emmanuelle Itam Name

Hybrid simulation of high-Reynolds number flows relying on a variational multiscale model

We are interested in the study and improvement of the LES component in hybrid RANS-LES formulations. These models are not designed for blunt-body flows with laminar boundary layers, but it is interesting to examine how they behave in that case. A DDES model is compared with a dynamic Variational multi-scale (DVMS) LES model for two subcritical flows past a cylinder. We then propose a hybridation restricting to DVMS in LES regions. The performances of the different options are compared for subcritical flows and for a flow around a tandem cylinder.

Introduction and motivations

This work takes place in a study of numerical simulation methods suited to industrial problems, equipped with turbulence models adapted to the simulation of turbulent flows with massive separations and vortex shedding. Recent contributions concerning this study can be found in [START_REF] Moussaed | Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows[END_REF][START_REF] Moussaed | Simulation of the flow past a circular cylinder in the supercritical regime by blending RANS and variational-multiscale LES models[END_REF]. The numerical approximation is second-order accurate, applicable to unstructured tetrahedral meshes, and involve a very low numerical dissipation.

In order to address high Reynolds number flows, we consider a RANS-LES hybridization similar to the Detached Eddy Simulation (DES) [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF][START_REF] Spalart | A new version of detachededdy simulation, resistant to ambiguous grid densities[END_REF]. While LES is a natural approach for computing subcritical blunt body flows, DES is not designed for this regime but can yet offer reasonable predictions [START_REF] D'alessandro | Detached eddy simulations of the flow over a cylinder at Re = 3900 using OpenFOAM a low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. This can be due to the low impact of the RANS component on the boundary layers at these regimes. In the wake, it remains important to have a LES mode as accurate as possible. In [START_REF] Reddy | A DDES model with a Smagorinsky-type eddy viscosity formulation and log-layer mismatch correction[END_REF], a formulation closer to the standard Smagorinsky model is introduced in a DES formulation by a modification of the energy production term. However, the shear layer between main flow and wake can suffer from the tendancy of a Smagorinsky-like model to reinforce the filtering in these regions. In the present study, we consider a DVMS LES model in which (1) the space-time dependent strength of the filter is controlled by a dynamic process, and (2) the width of the filter is numerically controlled by using a variational multiscale (VMS) formulation.

The present work focuses on the following issues:

• the evaluation of DDES and VMS models for the prediction of subcritical blunt body flows, • the definition of a hybrid model inheriting the accuracy of DVMS,

• the comparison of these models for a subcritical flow and for a well documented supercritical flow, the flow around a tandem cylinder.

Numerical model

The Navier-Stokes system for compressible fluid flow is solved. The spatial discretization is based on a mixed finite-volume/finite-element formulation, with degrees of freedom located at vertices i of the tetrahedrization. The finite volume part is integrated on a dual mesh built in 2D from median (Figure 1), in 3D from median plans. The diffusive fluxes are evaluated by a finite-element method, whereas a finite-volume method is used for the convective fluxes. The numerical approximation of the convective fluxes at the interface of neighboring cells is based on the Roe scheme [START_REF] Roe | Approximate Riemann solvers parameters vectors and difference schemes[END_REF]. In order to obtain second-order accuracy in space, the Monotone Upwind Scheme for Conservation Laws reconstruction method (MUSCL) [START_REF] Leer | Towards the ultimate conservative scheme. iv : a new approach to numerical convection[END_REF] is used, in which the Roe flux is expressed as a function of reconstructed values of the discrete flow variable W h at each side of the interface between two cells. We refer to [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF] for details on the definition of these reconstructed values. We just emphasize that particular attention has been paid to the dissipative properties of the resulting scheme since this is a key point for its successful use in LES. The numerical (spatial) dissipation provided by this scheme is made of sixth-order space derivatives [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF] and is then concentrated on a narrow-band of the highest resolved frequencies. This is expected to limit undesirable damping by numerical dissipation of the large scales. Moreover, this dissipation is O(∆ x 5 ) and a parameter γ S directly controls the amount of introduced viscosity and can be explicitly tuned in order to control the influence of numerical dissipation and, when necessary, reduce it to the minimal amount needed to stabilize the simulation. Time integration uses an implicit second-order backward differencing scheme. In [START_REF] Camarri | A low-diffusion MUSCL scheme for LES on unstructured grids[END_REF], it is shown that when used in combination with moderate time steps, the time dissipation is also quite small.

Turbulence model: VMS-LES

The Variational Multiscale (VMS) model for the large eddy simulation (LES) of turbulent flows has been introduced in [START_REF] Hughes | Large eddy simulation and the variational multiscale method[END_REF] in combination with spectral methods. In [START_REF] Koobus | A variational multiscale method for the large eddy simulation of compressible turbulent flows on unstructured meshes-application to vortex shedding[END_REF], an extension to unstructured finite volumes is defined. That method is adapted in the present work. Let us explain this VMS-LES approach in a simplified context. Assume the mesh is made of two embedded meshes, corresponding to a P 1 -continuous finite-element approximation space V h with the usual basis functions Φ i vanishing on all vertices but vertex i. Let be V 2h its embedded coarse subspace V 2h . Let V h be the complementary space: V h = V 2h ⊕ V h . The space of small scales V h is spanned by only the fine basis functions Φ i related to vertices which are not vertices of V 2h .

We write the compressible Navier-Stokes equations as follows:

∂W ∂t + ∇ • F(W ) = 0
where W = (ρ, ρu, ρE). The VMS-LES discretization writes for W h = ∑WiΦi:

∂W h ∂t , Φ i + (∇ • F(W h ), Φ i ) = -τ LES (W h ), Φ i (1)
.

For a test function related to a vertex of V 2h , the RHS vanishes, which limits the action of the LES term to small scales. In practice, embedding two unstructured meshes V h and V 2h is a constraint that we want to avoid. The coarse level is then built from the agglomeration of vertices/cells as sketched in Figure 2. It remains to Fig. 2 Building the VMS coarse level define the model term τ LES (W h ). This term represents the SGS stress term, acting only on small scales W h , and computed from the small scale component of the flow field by applying either a Smagorinsky or a WALE SGS model, [START_REF] Nicoud | Subgrid-scale stress modelling based on the square of the velocity gradient tensor[END_REF], the constants of these models being evaluated by the Germano-Lilly dynamic procedure [START_REF] Germano | A dynamic subgrid-scale eddy viscosity model[END_REF][START_REF] Lilly | A proposed modification of the Germano subgrid-scale closure method[END_REF][START_REF] Moussaed | Impact of dynamic subgrid-scale modeling in variational multiscale large-eddy simulation of bluff body flows[END_REF]. The main property of the VMS formulation is that the modeling of the unresolved structures is influencing only the small resolved scales, as described in Figure 3, in contrast with, for example the usual Smagorinsky model. This implies two main properties. First, the backscatter transfer of energy to large scales is not damped by the model. Second, the model is naturally unable to produce an artificial viscous layer at a no-slip boundary, or in shear layers.

Turbulence model: RANS-DDES

The RANS component used in the present study is the kε model proposed by Goldberg [START_REF] Goldberg | A wall-distance-free kε model with enhanced near-wall treatment[END_REF]. The turbulent viscosity is limited by the Bradshaw's law in a similar way to Menter's SST model. In what we call DDES/kε/Menter, the above RANS

D DDES k = ρ k 3/2 l DDES
where

l DDES = k 3/2 ε -f d * max(0, k 3/2 ε -C DDES ∆ ) with f d = 1 -tanh((8r d ) 3 ) and r d = ν t + ν max( √ u i, j u i, j , 10 -10 )K 2 d 2 w .
K denotes the von Karman constant (K = 0.41), d w the wall-normal distance, u i, j the x j -derivative of the ith-component of the velocity u, and the model constant C DDES is set to the standard value 0.65 (ν t and ν are the turbulent kinematic viscosity and the kinematic viscosity of the fluid, respectively).

Turbulence model: hybrid RANS/VMS-LES

The LES part of the above RANS-DDES is a kind of Smagorinsky model, but not exactly as remarked in [START_REF] Reddy | A DDES model with a Smagorinsky-type eddy viscosity formulation and log-layer mismatch correction[END_REF]. The central idea of the hybrid VMS model which we propose is to combine the mean flow field obtained with the RANS component with the application of the DVMS model wherever the grid resolution is adequate. First, let us write the semi-discretization of the RANS equations :

∂ W h ∂t , Φ i + (∇ • F( W h ), Φ i ) = -τ RANS ( W h ), Φ i .
A natural hybridation writes:

∂W h ∂t , Φ i + (∇ • F(W h ), Φ i ) = -θ τ RANS ( W h ), Φ i -(1 -θ ) τ LES (W h ), Φ i (5)
where W h denotes now the hybrid variables and θ is defined by:

θ = 1 -f d (1 -tanh(ξ 2 )) with ξ = ∆ l RANS or ξ = µ SGS µ RANS
where f d is defined as in Section 5.

Results

We consider three flows around one or two cylinders. The flow around a cylinder is strongly dependent on the turbulence which exists in the inflow and the turbulence which is created along the wall boundary layer after the stagnation point. In the case where the inflow involves no turbulence, turbulence is created on the wall at a sufficiently high critical Reynolds number, corresponding to the drag crisis. This Reynolds number is between 300, 000 and 500, 000. Under this number, the flow is subcritical and, in short, the boundary layer flow is laminar. Then in principle neither a RANS nor a hybrid model should be applied, but it is interesting to be able to use a unique model for various regimes. A first series of calculations with a Reynolds number of Rey = 3900, is compared with other LES and DES computations [START_REF] Lee | A dynamical subgrid-scale eddy viscosity model with a global model coefficient[END_REF][START_REF] Kravchenko | Numerical studies of flow over a circular cylinder at re=3900[END_REF][START_REF] D'alessandro | Detached eddy simulations of the flow over a cylinder at Re = 3900 using OpenFOAM a low-diffusion MUSCL scheme for LES on unstructured grids[END_REF] and with measurements [START_REF] Parneaudeau | Experimental and numerical studies of the flow over a circular cylinder at Reynolds number 3900[END_REF][START_REF] Norberg | Effects of Reynolds number and low-intensity free-sream turbulence on the flow around a circular cylinder[END_REF]. Three model options are compared. One can notice that DDES behaves very similarly to the no model variant (that is to say without any turbulence model), due to a too low level of turbulence introduced by this model. This behavior is reported in the literature, see e.g. [START_REF] D'alessandro | Detached eddy simulations of the flow over a cylinder at Re = 3900 using OpenFOAM a low-diffusion MUSCL scheme for LES on unstructured grids[END_REF]. In contrast, the VMS-based calculation is closer to experiments, as illustrated by the bulk coefficients shown in Table 1.

Circular cylinder 20000 -subcritical

A second example is the flow around the circular cylinder with a Reynolds number of 20, 000, more than five times larger. 2 Bulk flow parameters at Re=20000. Θ sep is the separation angle, Lv denotes the x-location of the maximum in the turbulent intensity distribution along the wake axis, and Imax holds for the maximum turbulence intensity along the wake axis. The other symbols are the same as in Table 1.

C d C
We use a mesh of 1.8M nodes. Comparisons are done with a rather recent set of measures by Norberg [START_REF] Norberg | Fluctuating lift on a circular cylinder: review and new measurements[END_REF] and Lim-Lee [START_REF] Lim | Flow control of circular cylinders with longitudinal grooved surfaces[END_REF]. We first verify that a computation without any turbulence model produces a fair result with a nearly 10% error on mean drag, an overprediction of rms of lift coefficient of nearly 30%. The application of our URANS gives an even less accurate of bulk coefficients, except pressure at rear stagnation point C P b . Our DDES model, built with the Golberg/Menter RANS model, is not so bad but underpredicts several bulk coefficients (rms of lift coefficient, recirculation length and x-location of the maximum turbulent intensity along the wake axis). In contrast, a LES calculation with the VMS-Wale-dynamic options gives reasonable predictions. As for our hybrid VMS model, a good overall agreement with the experimental data is observed. The accuracy of the prediction is lower than for the LES for the recirculation length and the Strouhal number, but higher (in particular for C Lrms and Cp b ) than for DDES. The calculation of these two subcritical flows tends to show that, for this Reynolds number interval, DDES (1) behaves in a way very close to a no model formulation for Re = 3900, (2) is less predictive than LES models, and in particular the VMS dynamic model which we consider, (3) is less accurate than the proposed Hybrid model for Re = 20000.

Comparisons for subcritical high Reynolds numbers like 140000 introduce a second difficulty for DES/DDES computations since they predict a supercritical (or close to supercritical) flow, see e.g. [START_REF] Travin | Detached-eddy simulations past a circular cylinder[END_REF]. We retain that these results tend to confirm the interest to re-introduce a sophisticated LES model in hybrid formulations.

Tandem cylinder

The last case concerns the flow around two cylinders in tandem, at Re=1.66 × 10 5 . This test case was studied in an AIAA workshop, see [START_REF] Lockard | Summary of the tandem cylinder solutions from the benchmark problems for airframe noise computations-i[END_REF]. Among the conclusions of the workshop, addressing this case with LES was considered as a too difficult task with existing computers. Results of DES-based computations were much closer to measurements. We have computed the case applying our DDES version (with a mesh of 2.7M nodes) and applying the hybrid VMS-LES-Smagorinsky-dynamic combination (with a mesh of 2.59M nodes) . Some comparison of bulk coefficients is given in Table 3. Experimental outputs are taken from [START_REF] Neuhart | Measurements of the flowfield interaction between tandem cylinders[END_REF]. The drag coefficient for the first cylinder shows a small scatter with the various hybrid calculations. The simulation around the second cylinder is more challenging as also illustrated by the pressure and the resolved turbulent kinetic distributions (Figure 4). One can observe, first, that the drag prediction obtained by the DDES model for the second cylinder is in accordance with that of Aybay [START_REF] Lockard | Summary of the tandem cylinder solutions from the benchmark problems for airframe noise computations-i[END_REF], who used a DES/kω/SST model, and second, that this prediction is improved by our hybrid VMS-LES model with a deviation to measurements which is two times smaller.

A short provisional synthesis

Applying a hybrid model on subcritical cylinder flows provokes two kind of prediction losses. First, the turbulent treatment of the boundary layers deviates from the real flow characteristics. It is known that close to critical Reynolds number, this deteriorates the prediction. For lower Reynolds number (3900) this effect is rather small and the model remains a predictive one. Second, the treatment of wake by the hybrid model is of lower accuracy than a good LES model. We propose a way to reintroduce a better LES capturing which tends to improve the results for the two last flows of our study.

Fig. 1

 1 Fig. 1 Dual cell in 2D

Fig. 3

 3 Fig. 3 VMS principle.

Fig. 4

 4 Fig. 4 Tandem cylinder -Left : Mean pressure coefficient distribution on the upstream cylinder (top frame) and on the downstream cylinder (bottom frame) -Right : Resolved turbulent kinetic energy along the centerline between the cylinders (top frame) and behind the second cylinder (bottom frame).

Table 1

 1 Bulk quantities for Re=3900 flow around a cylinder. C d holds for the mean drag coefficient, C p b for the mean pressure coefficient at cylinder basis, C Lrms for the root mean square of lift time fluctuations, Lr is the recirculation length, and St the Strouhal number.

	7.1 Circular cylinder 3900 -subcritical
		Mesh	C d	-C p b C Lrms	L r	S t
	Experiments
	Norberg [18]		[0.94-1.04] [0.83-0.93]	-	-	-
	Parnaudeau [20]		-	-	-	1.51	0.210
	Present simulations
	No model	1.4M	0.87	0.73	0.04	2.11	0.209
	DDES/k -ε/Menter	1.4M	0.88	0.74	0.03	2.07	0.218
	VMS dyn.	1.4M	0.96	0.84	0.12	1.54	0.215
	Other simulations
	Lee (LES) [9]	7.7M	[0.99-1.04] [0.89-0.94]		[1.35-1.37] [0.209-0.212]
	Kravchenko (LES) [8] [0.5M-2.4M] [1.04-1.38] [0.93-1.23]		[1.-1.35] [0.193-0.21]
	SA-IDDES [2]	3.9M	0.98	0.83	0.109	1.67	0.214
	v2 -f DES [2]	3.9M	1.02	0.87	0.14	1.42	0.222

Table 3

 3 Tandem cylinder: Mean drag coefficient (experimental coefficients are computed by integrating experimental pressure). Circled superscripts hold for cylinder 1 and cylinder 2.

		Mesh	C	1	Lrms	St 1
	Experiment
	Neuhart [16] and Lockard [13]
	(BART facility)		0.64	0.31	0.231
	Present simulations
	DDES/k -ε/Menter	2.7M	0.65	0.44 0.156 0.59 0.214
	Hybrid dyn.	2.59M 0.64	0.38 0.077 0.79 0.214
	Other simulations
	DES k -ω SST, Aybay [26]	6.7M	0.64	0.44	0.223
	HRLES, Vatsa[26]	8.7M	0.64	0.45	0.208
	S-A MDDES, Lockard (2014)	16M	0.50	0.45 0.072 0.643 0.229
	S-A DDES, Garbaruk [3]	11M	0.48	0.42 0.078 0.612 0.244

d C 2 d C 1 Lrms C 2

Acknowledgements This work has been supported by French National Research Agency (ANR) through project MAIDESC n o ANR-13-MONU-0010. HPC resources from GENCI-[CINES] (Grant 2010-x2010026386 and 2010-c2009025067) are also gratefully acknowledged.