Cosme Louart

Romain Couillet

A RANDOM MATRIX AND CONCENTRATION INEQUALITIES FRAMEWORK FOR NEURAL NETWORKS ANALYSIS

Keywords: Neural networks, random matrix theory, concentration inequalities, extreme learning machines

This article provides a theoretical analysis of the asymptotic performance of a regression or classification task performed by a simple random neural network. This result is obtained by leveraging a new framework at the crossroads between random matrix theory and the concentration of measure theory. This approach is of utmost interest for neural network analysis at large in that it naturally dismisses the difficulty induced by the non-linear activation functions, so long that these are Lipschitz functions. As an application, we provide formulas for the limiting law of the random neural network output and compare them conclusively to those obtained practically on handwritten digits databases.

INTRODUCTION

One of the main popularity features of deep neural networks lies in their (still barely fathomed) performance stability. That is, as the number n and size p of the training data grow large (and so does the network), independently of the random initialization point of the backpropagation learning algorithm, essentially the same performances are ultimately achieved. This characteristic is at the core of a current stream of research, based on tools from statistical physics and random matrix theory [START_REF] Dauphin | Identifying and attacking the saddle point problem in high-dimensional nonconvex optimization[END_REF][START_REF] Choromanska | The loss surfaces of multilayer networks[END_REF][START_REF] Pennington | Geometry of neural network loss surfaces via random matrix theory[END_REF], aiming at theorizing these observations. And, indeed, [START_REF] Choromanska | The loss surfaces of multilayer networks[END_REF] explored a model akin to deep neural networks and concludes that the local minima of the learning cost function become increasingly dense as data and network grow large, having essentially the same associated loss, and that the probability to escape these minima vanishes. However, the statistical physics model of [START_REF] Choromanska | The loss surfaces of multilayer networks[END_REF] is in reality far from satisfying from a neural network perspective as it, for once, breaks all dependence induced by the non-linearities of the activation functions (ReLU non-linearities being replaced by products with independent Bernoulli random variables) and, most importantly, assumes data to be constituted of random independent entries; both conditions ensure that only random uncorrelated scalars propagate through the network, a highly criticizable model for deep nets. Alternative neural network analyses discard the non-linearities altogether, as in [START_REF] Saxe | Exact solutions to the nonlinear dynamics of learning in deep linear neural networks[END_REF]; but in this case, convergence to global minima are studied, which are known not to be achieved by practical deep networks (fortunately so, as this avoids overfitting).

In this article, following our seminal works [START_REF] Louart | Harnessing neural networks: a random matrix approach[END_REF][START_REF] Louart | A random matrix approach to neural networks[END_REF], we propose a different angle of approach to neural network analysis. Rather than modelling a complete deep neural net, we focus here primarily on simple network structures, so far not considering backpropagation learning but accounting for nonlinearities induced when traversing a hidden layer. The main technical driver to this endeavor is the concentration of measure theory, which has the key features of (i) extending many results holding for vectors of independent entries to the wider scope of concentrated random vectors (see definition in subsequent sections) and most importantly of (ii) being a theory "stable to Lipschitz mappings" in that Lipschitz functions f : R p → R q of concentrated random vectors in R p are still concentrated vectors in R q . Feature (ii) notably allows one to accommodate with the non-linear activation functions, so long that these are Lipschitz (e.g., ReLU, sigmoid maps).

In [START_REF] Louart | Harnessing neural networks: a random matrix approach[END_REF], we merely exploited Feature (ii) as a technical means to study the asymptotic (as n, p → ∞) performance of extreme learning machines (ELM) [START_REF] Huang | Extreme learning machine: theory and applications[END_REF] (i.e., single hiddenlayer regression networks with no backpropagation learning), assuming a model encompassing a random connectivity matrix (which induces the concentration of the output vectors) but deterministic data. Under this model, however, while the asymptotic network training performance was readily accessible, the asymptotic generalization performance remained out of technical grasp and only a conjecture under "reasonable" yet unclear assumptions on the deterministic dataset could be proposed. As an answer, the present study strongly suggests that a key property of neural network stability (and likely of many statistical learning methods) lies, not in the (initial) randomness of the inter-layer connections, but rather in a concentration property of the dataset. This property structurally appears when studying the orders of magnitude of the output of an ELM for concentrated versus deterministic data (in the former case the output has a controlled magnitude, while in the latter case the output may diverge as n, p → ∞). Exploiting Features (i) and (ii) together, and therefore working on concentrated input data in the first place, brings us to a more generic analysis framework, where Lipschitz non-linearities need not be explicitly studied as they do not affect the concentration properties of the data.

Under this setting, in the present work, we generalize several results from random matrix theory, by providing notably a deterministic equivalent for the covariance matrix of a kmixture of concentrated vectors along with the asymptotic statistical behavior of a ridge-regression on these vectors. As an immediate aftermath, the asymptotic performance analysis of ELMs, and in passing the spectral characterization of random feature maps, then reduce to mere corollaries by specifying the structural properties of the concentrated vectors.

Notation. In the remainder, C, c > 0 as well as C , c > 0 are constants independent of all other parameters, and C , c > 0 constants dependent only on C, c.

SYSTEM SETTING

Basic notions of concentration of measure

We start our system modelling assumptions with a few definitions and essential notions of the concentration of measure theory [START_REF] Ledoux | The concentration of measure phenomenon[END_REF] that will be used in this article.

Definition 1 (Concentration of measure). The random variable Z ∈ R is said to be concentrated and we denote Z ∈ α( • ) if, for Z an independent copy of Z,

P (|Z -Z | ≥ t) ≤ α(t).
In particular, Z is normally concentrated, denoted

Z ∈ Cα N ( √ c • ), if Z ∈ Ce -c( • ) 2 , and Z is exponentially con- centrated, denoted Z ∈ Cα exp (c • ), if Z ∈ Ce -c • .
Normal and exponential concentrations owe their names to the fact that normal and exponential random variables are respectively normally and exponentially concentrated. The fast decay of their tails allows for additional properties. In particular, both concentrate around their means in the sense that, e.g.,

Z ∈ Cα N ( √ c • ) ⇒ P(|Z -EZ| ≥ t) ≤ C e -c t 2
. Also, these fast concentrations induce moment controls:

Z ∈ Cα N ( √ c • ) ⇒ ∀r > 0, E[|Z -EZ| r ] ≤ C (2r/c) r/2 Z ∈ Cα exp (c • ) ⇒ ∀r > 0, E[|Z -EZ| r ] ≤ C (2r/c) r .
As f (|Z -Z |) ≤ λ|Z -Z | for λ-Lipschitz functions f , we have the following structural property of the theory.

Property 1 (Lipschitz maps). For f : R → R a λ-Lipschitz function and Z ∈ R,

Z ∈ α( • ) ⇒ f (Z) ∈ α( • /λ).
Similarly, linear combinations of concentrated random variables remain concentrated. Products of concentrated random variables are more difficult to handle, but we have the following lemma, of importance in this article.

Lemma 1 (Concentration of squared variables1 ).

If Z ∈ Cα N ( √ c • ), then, with obvious notations, Z 2 ∈ C α exp c 2 • + C α N √ c 4|EZ|
• .

The concentration of measure theory however finds its fullest significance when considering random vectors (rather than scalars) Z ∈ R p . As most random vectors of practical interest do not localize (e.g., large Gaussian vectors tend to spread along a sphere), the notion of concentration of measure for vectors is defined by means of all their Lipschitz "observations".

Definition 2 (Concentration of random vectors). A vector

Z ∈ R p is concentrated, denoted Z ∈ α( • ), if for every 1-Lipschitz map f : R p → R, f (Z) ∈ α( • ).
Normally and exponentially distributed vectors are defined similarly.

A non-trivial result [8, Prop. 1.9.] is that standard normal random vectors Z ∈ R p are indeed normally distributed, in the sense of Definition 2, with parameters C, c independent of p. Precisely, Z ∼ N (0,

I p ) ⇒ Z ∈ 2α N ( • / √ 2).
With these notions at hand, we are in position to present our work setting.

Setting

Keeping in mind that neural networks are mostly used for regression or classification, we consider here a set of inputoutput data pairs (x 1 , y 1 ), . . . , (x n , y n ) with x i ∈ R p and y i ∈ R d (p will be supposed large while d remains small irrespective of p). Our core assumption is that there exists k measures µ 1 , . . . , µ k such that, for each l ≤ n, x l ∼ µ for some ≤ k; besides all x l are independent and

Z ∼ µ ⇒ Z ∈ C α N ( √ c • ).
For Z ∼ µ , we denote M ≡ E[Z] (and M = [ M1 , . . . , Mk ]) and C ≡ E[ZZ T ] (not to be confused with the covariance matrix). To avoid technicalities, we assume that C < max 1≤ ≤k { 1 p tr C } < C for C, C > 0 independent of p. We finally denote n the number of x l 's drawn from µ .

This set of hypotheses is of interest as it notably encompasses the cases where:

1. the x l 's arise from a Gaussian mixture model where µ = N (m , Σ );

2. the x l 's are the output of a random feature map x l = σ(W s l ) for W ∈ R p×q deterministic with W ≤ 1, σ : R → R a 1-Lipschitz map (operating here entrywise), and

s l ∈ R q such that s l ∈ C α N ( √ c • ) (for instance s l ∼ N (m , Σ )).
Item 2 justifies our claim that, under this setting, simple random neural networks analysis reduces to the analysis of concentrated random vectors, disregarding the specificities of the non-linear activation function since

s l ∈ C α N ( √ c • ) ⇒ x l ∈ C α N ( √ c • )
and that only this concentration property will be effectively used in the forthcoming analysis.

MAIN RESULTS

Under the assumptions of Section 2.2, extreme learning machines [START_REF] Huang | Extreme learning machine: theory and applications[END_REF] may be merely seen as a linear ridge-regression with training set (x 1 , y 1 ), . . . , (x n , y n ), where x l = σ(W s l ) for some input observations s 1 , . . . , s n . Since linear ridgeregression involves as a core object the sample covariance matrix

C X ≡ 1 n XX T = 1 n n i=1 x i x T i with X ≡ [x 1 , . . . , x n ] ∈ R p×n
, our first objective is to characterize the eigenspectrum of this matrix for large p and n. This first analysis shall subsequently allow for a full characterization of functionals of the random matrix X, among which the performance of ELM regression and classification.

Sample covariance matrix analysis

To tackle the eigenspectrum analysis of C X , random matrix theory [START_REF] Tao | Topics in random matrix theory[END_REF][START_REF] Couillet | Random Matrix Methods for Wireless Communications[END_REF] provides a quite versatile tool: the resolvent

Q X ≡ (C X + zI p ) -1
of C X , for z > 0. 2 The matrix Q X encapsulates much information about C X ; notably, the so-called Stieltjes transform 1 p tr Q X uniquely characterizes the empirical eigenvalue distribution 1 p p i=1 δ λi(C X ) of C X [13, Th. B.9.], while quadratic forms of the type a T Q X a allow for a characterization of the projections |a T u i (C X )| of deterministic vectors a ∈ R p on the "isolated" eigenvectors u i (C X ) of C X [START_REF] Benaych-Georges | The singular values and vectors of low rank perturbations of large rectangular random matrices[END_REF].

Besides, Q X is a convenient tool for the present article as it naturally transfers concentrations. Indeed, since

Q X ≤ z -1 and Q X X ≤ √ n, it is easily shown that the mapping X → Q X is 2/ √ nz-Lipschitz so that X ∈ α( • ) ⇒ Q X ∈ α( √ nz/2 • ) (with X and Q X respectively seen as vectors in R np and R p 2 ).
Owing to this property, we have the following core result.

Theorem 1. Let Q be defined as

Q ≡ k =1 n n C 1 + δ + zI p -1
2 Usually one defines Q X as Q X = (C X -zIp) -1 for z ∈ C with z > 0 but analyticity arguments along with the nonnegative definiteness of C X justify the equivalence of this definition here.

where (δ 1 , . . . , δ k ) ∈ R k + is the unique solution with nonnegative elements of the system δ = 1 n tr C Q, = 1, . . . , k. Then, for any unit-norm u ∈ R p , the following nonasymptotic inequalities hold

E[Q X ] -Q ≤ C p n max 1, p n √ n P 1 p tr (Q X -Q) ≥ t + C max 1, p n √ n ≤ Ce -cnpt 2 P u T (Q X -Q)u ≥ t + C p n max 1, p n √ n ≤ Ce -cnt 2 .
Sketch of Proof. Given the bound on EQ X -Q , the second and third results follow from Q X -Q = (Q X -EQ X ) + (EQ X -Q) and normal concentration inequalities based on the Lipschitz character of X → Q X . We are then left to prove the first result. Let

X -i = [x 1 , . . . , x i-1 , x i+1 , . . . , x n ], Q -i = Q X-i and Q = Q X . By the identity Qx i (1 + 1 n x T i Q -i x i ) = Q -i x i , we have E[Q] -Q = k =1 n n E[ε 1 + ε 2 ]
where, for any x l ∼ µ , we introduced the matrices

ε 1 ≡ Q -l x l x T l Q(δ -1 n x T l Q -l x l ) (1 + 1 n x T l Q -l x l )(1 + δ ) , ε 2 ≡ Q -l x l x T l Q C Q n(1 + δ ) .
Because of the leading factor 1 n in ε 2 , algebraic manipulations similar to [START_REF] Louart | A random matrix approach to neural networks[END_REF] ensure that this term vanishes faster than ε 1 . As for ε 1 , the main technical part is to control

x T l Q -l x l -nδ = (x T l Q -l x l -tr C Q -l ) + tr C (Q -l -EQ -l ) + (tr C EQ -l -nδ ). Since Q -l ≤ z -1 and x l ∈ C α N ( √ c • ), Lemma 1 ap- plied to x T l Q -l x l = Q 1 2
-l x l 2 ensures that the first righthand side difference is normally-exponentially concentrated. Since X → Q X is Lipschitz, the second difference is normally concentrated. As for the third (deterministic) term, its control follows from pre-established random matrix results (see e.g., [START_REF] Benaych-Georges | Spectral analysis of the gram matrix of mixture models[END_REF]). Since normal and exponential concentrations convert to bounds on moments, E[ε 1 ] can be appropriately bounded, thereby completing the proof.

Ridge-regression and ELM classification

Linear ridge-regression for the training data pairs (x i , y i ) previously defined consists in determining the vector β ∈ R p×d that minimizes, for z > 0, the cost

E train (β) ≡ 1 n n i=1 y i -β T x i 2 F + z β 2 F
(with • F the Frobenius norm). Letting Y = [y 1 , . . . , y n ] T ∈ R n×d , the solution is explicitly given by β = 1 n Q X XY . For a test data pair (x, y) ∈ R p×d , the regression output of x is then given by

S(x) = x T β = 1 n x T Q X XY.
For k-class classification purposes, one naturally takes y l = e ∈ R k when x l ∼ µ (hence d = k), with [e ] a = δ a the indicator vector of class . This procedure strongly relates to the kernel LS-SVM approach [START_REF] Suykens | Least squares support vector machine classifiers[END_REF][START_REF] Liao | A large dimensional analysis of least squares support vector machines[END_REF].

The associated regression and classification performance measures are the mean-square error E[ S(x) -y 2 ] and the misclassification rate, respectively. In both cases, these relate to the probability distribution of S(x), which is then our object of present interest. Our main result, restricted for readability to the classical random matrix regime on n, p, reads:

Theorem 2. Let x ∼ µ and assume max 1≤i≤n y i < C. Then, as n, p → ∞ with p/n → γ ∈ (0, ∞), V -1 2 S(x) -S L -→ N (0, I d ) where V = V -S ST , S ≡ 1 n M T Q M ∆J T Y V ≡ 1 + δ n 2 Y T J∆ [Φ] ,•,• + D Ψ - n n (D Θ + ΘD ) ∆J T Y where J = [j 1 , . . . , j k ] ∈ R n×k with [j ] i = δ xi∼µ , Φ = (I k -Ψ) -1 Φ ∈ R k×k×k , 3 Ψ = (I k -Ψ) -1 Ψ ∈ R k×k , D = diag(Ψ •, ) ∈ R k×k , D Ψ = diag( n 2 (1+δ•) n 2 • Ψ ,• ) ∈ R k×k , ∆ = diag((1 + δ • ) -1 ), and
Θ = M T Q M Φ = M T i Cl Q Mj 1 + δ l 1≤l,i,j≤k Ψ = n j n 2 tr Q Ci Q Cj (1 + δ i )(1 + δ j ) 1≤i,j≤k .
Proof. Asymptotic means and covariances follow from a concentration inequality-based analysis. The central limit is then obtained from a refined version of [START_REF] Chatterjee | Fluctuations of eigenvalues and second order poincaré inequalities[END_REF] adapted to our present setting. Details are provided in the extended article [START_REF] Louart | Sample covariance of concentrated random vectors[END_REF].

As discussed previously, letting x l = σ(W s l ) for s l ∈ R q with q ∼ p such that s

l ∈ C α N ( √ c • ), σ(•) 1-Lipschitz and W ∈ R p×q with W ≤ 1 (thereby ensuring that x l ∈ C α N ( √ c • )), S(x)
is the output of an ELM. The asymptotic statistics of S(x) in Theorem 2 therefore directly translate in terms of simple neural network performances. 3 Here we understand the product AB for A ∈ R k×k and B ∈ R k×k×k by [AB] abc = d A ad B dbc . Figure 1 provides the simulated S(x) ∈ R 2 versus theoretical average S output of an ELM for 2-class (C 1 , C 2 ) classification with σ(t) = ReLU(t) = max(t, 0) and σ(t) = erf(t) = 2 √ π t 0 e -u 2 du, W random unitary, and s l extracted from the MNIST handwritten digits dataset [START_REF] Lecun | The MNIST database of handwritten digits[END_REF]. Here Y ∈ R n×k is defined by Y l = δ s l ∈C . Note the accurate fit between theory and practice, suggesting that the MNIST data are conveniently modelled as concentrated random vectors.

Figure 2 proceeds similarly with Gaussian s l ∼ N (0, α I p ) with α ∈ {1, 2} according to the class. Being non-linearly separable classes, a straightforward application of Theorem 2 ensures that both classes are non-discriminable for σ(t) such that σ(-t) = -σ(t), which Figure 2 visually confirms; indeed, for such symmetric σ, for all , M = E[σ(W s l )] = 0 where s l ∼ N (0, α I p ), and we thus find that S = 0 indeed not allowing for class discrimination.

10 - 3 10 -1 10 1 10 3 0 1 10 - 3 10 -1 10 1 10 3 Fig. 1 . 10 - 3 10 -1 10 1 10 3 0 1 2 10 - 3 10 -1 10 1 10 3 Fig. 2 .

 33133110312332 Fig. 1. Scores [S(x)] 1 (blue dashed) and [S(x)] 2 (blue solid) of 2-class MNIST ELM as a function of regularization z ∈ [10 -4 , 10 4 ] (digits 3 for C 1 and 8 for C 2 ) for x ∈ C 1 , versus theory (red dashed and solid). Based on n = 2048 samples, p = q = 784, W random unitary.

This result is similar to the Hanson-Wright inequality [9, Th. 6.2.1], but more adequate to our present setting; the proof, provided in an extended version of this article[START_REF] Louart | Sample covariance of concentrated random vectors[END_REF], is also structurally simpler as it relies on more elementary properties.

This work is supported by the ANR RMT4GRAPH Project (ANR-14-CE28-0006).


Loïc Frazza

Adrien Loseille

Frédéric Alauzet

Alain Dervieux

Nonlinear corrector for RANS equations

The scope of this paper is to present a nonlinear error estimation and correction for Navier-Stokes and RANS equations. This correction is obtained by deducing a source term from the evaluation of the residual of the solution interpolated on the h/2 mesh. To avoid the generation of the h/2 mesh (which is prohibitive for realistic applications), the residual at each vertex is computed by local refinement only in the neighborhood of the considered vertex. It successfully improves solution predictions and yields a sharp estimate of the numerical error.

I. Introduction

The numerical simulation of engineering problems involves a discrete solution which is alleged to converge toward the continuous solution of the problem as the mesh-size decreases. As the exact analytical solution of the problem is sought, the difference between the numerical and the analytical solution can be seen as a noise. In an engineering context, this noise can have disastrous consequences on the numerical prediction which can, in turn, lead to actual accidents or misconceptions.

It is thus essential to reduce and estimate this error. This is usually done by leading a grid convergence study, where the error is estimated with the difference between two solutions computed on two successives grids. However, this approach requires an additional finer grid (thus expensive) to appropriately estimate the error. Moreover, it requires the expertise of an engineer. In a mesh adaptation context multiple successive grids are automatically generated following a sensor which estimates the error related to the local mesh size. It does not requires any mesh prescription, but it needs a relatively accurate estimation of the numerical error to be efficient and to stop when the mesh and the solution are converged.

Recently, different adjoint based mesh adaptation strategies [START_REF] Fidkowski | Review of Output-Based Error Estimation and Mesh Adaptation in Computational Fluid Dynamics[END_REF][START_REF] Venditti | Adjoint Error Estimation and Grid Adaptation for Functional Outputs: Application to Quasi-One-Dimensional Flow[END_REF][START_REF] Venditti | Grid adaptation for functional outputs: application to two-dimensional inviscid flows[END_REF][START_REF] Formaggia | Anisotropic error estimates for elliptic problems[END_REF][START_REF] Jones | Validation of 3D Adjoint Based Error Estimation and Mesh Adaptation for Sonic Boom Reduction[END_REF][START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] have been developed in order to optimize the mesh with respect to a specific functional output. These techniques uses the adjoint solution to estimate the error committed on the output and minimize it with respect to the mesh size. However, these approaches are based on a linearization of the equations and can thus be inefficient for strongly nonlinear problems like RANS equations. This is why we proposed a nonlinear equivalent of these approaches in order to efficiently compute a correction in the global solution and in the output functional (lift, drag, heat transfer, sound, ...).

II. Adjoint based error estimation and correction

A. Problem setup and Numerical discretization

Navier-Stokes RANS equations

The compressible Navier-Stokes equations for mass, momentum and energy conservation read:

             ∂ρ ∂t + ∇ • (ρu) = 0 , ∂(ρu) ∂t + ∇ • (ρu ⊗ u) + ∇p = ∇ • T , ∂(ρE) ∂t + ∇ • ((ρE + p)u) = ∇ • (T • u) + ∇ • (λ∇T ) , (1)
where ρ denotes the density (kg/m 3 ), u the velocity (m/s), E the total energy per mass (m 2 .s -2 ), p the pressure (N/m 2 ), T the temperature (K), µ the laminar dynamic viscosity (kg/(m.s)) and λ the laminar conductivity. T the laminar stress tensor:

T = (µ + µ t ) (∇ ⊗ u + t ∇ ⊗ u) - 2 3 ∇. u I ,
where (in 3D) u = (u, v, w) and

∇. u I =    u x + v y + w z 0 0 0 u x + v y + w z 0 0 0 u x + v y + w z    ,
where u x = ∂u ∂x , u y = ∂u ∂y , u z = ∂u ∂z (idem for v and w). The variation of nondimensionalized laminar dynamic viscosity and conductivity coefficients µ and λ as a function of the dimensional temperature T is defined by Sutherland's law:

µ = µ ∞ T T ∞ 3 2 T ∞ + Su T + Su and λ = λ ∞ T T ∞ 3 2 T ∞ + Su T + Su ,
where Su = 110 is the Sutherland constant and the index ∞ denotes reference quantities. The relation linking µ and λ is expressed from the Prandtl laminar number: Pr = µc p λ with Pr = 0.72 for (dry) air , where C p is the specific heat at constant pressure.

Turbulence modeling

According to the standard approach to turbulence modeling based upon the Boussinesq hypothesis, the turbulence is modeled with an eddy viscosity µ t , which is added to the laminar (or dynamic) viscosity µ.

The dynamic viscosity is usually taken to be a function of the temperature, whereas µ t is obtained using a turbulence model. Here we choose the Spalart-Allmaras one equation turbulence model [START_REF] Spalart | A one-equation turbulence model for aerodynamic flows[END_REF] given by the following equation:

∂ ν ∂t + u • ∇ν = c b1 [1 -f t2 ] S ν -c w1 f w - c b1 κ 2 f t2 ν d 2 + 1 σ ∇ • ((ν + ν)∇ν) + c b2 ∇ν 2 + f t1 ∆u 2 , ( 2
)
where ν is the turbulent kinematic viscosity and all the constants are defined below. In the standard model the trip term is being left out, i.e., f t1 = 0. Moreover, some implementations also ignore the f t2 term as it is argued that if the trip is not included, then f t2 is not necessary. [START_REF] Eca | Verification of RANS Solvers with Manufactured Solutions[END_REF] In Wolf, this simplified version has been considered and we prefer to write it under the following form, which is more appropriate for its discretization with the finite element/finite volume method. Indeed, Equation ( 2) can be decomposed into the following terms:

∂ρν ∂t + u • ∇ρν convection = c b1 Sρν production -c w1 f w ρ ν d 2 destruction + ρ σ ∇ • ((ν + ν)∇ν) dissipation + c b2 ρ σ ∇ν 2 dif f usion .
Notice that this is not a conservative model. If a conservative form of the Spalart-Allmaras is foreseen, we have to consider the variation proposed by Catris and Aupoix. [START_REF] Catris | Density corrections for turbulence models[END_REF] The turbulent eddy viscosity is computed from:

µ t = ρνf v1 ,
where

f v1 = χ 3 χ 3 + c 3 v1 and χ = ν ν with ν = µ ρ .
Additional definitions are given by the following equations:

f v2 = 1 - χ 1 + χf v1 and S = Ω + ν κ 2 d 2 f v2 where Ω = ∇ × u .
d is the distance to nearest wall which is computed for each vertex at the beginning of the simulation. The set of closure constants for the model is given by

σ = 2 3 , c b1 = 0.1355 , c b2 = 0.622 , κ = 0.41 , c w1 = c b1 κ + 1 + c b2 σ , c w2 = 0.3 , c w3 = 2 , c v1 = 7.1 .
Finally, the function f w is computed as:

f w = g 1 + c 6 w3 g 6 + c 6 w3 1/6
with g = r + c w2 r 6 -r and r = min ν Sκ 2 d 2 , 10 .

Vector form of the RANS system

We write the RANS system in the following (more compact) vector form:

W t + F 1 (W ) x + F 2 (W ) y + F 3 (W ) z = S 1 (W ) x + S 2 (W ) y + S 3 (W ) z + Q(W ) ,
where S i (W ) a = ∂S i (W ) ∂a (i = 1, 2, 3, a = x, y, z) (idem for F ). W is the nondimensionalized conservative variables vector: W = (ρ, ρu, ρv, ρw, ρE, ρν) T .

F (W ) = (F 1 (W ), F 2 (W ), F 3 (W )
) are the convective (Euler) flux functions:

F 1 (W ) = ρu, ρu 2 + p, ρuv, ρuw, u(ρE + p), ρuν T , F 2 (W ) = ρv, ρuv, ρv 2 + p, ρvw, v(ρE + p), ρvν T , (3)
F 3 (W ) = ρw,
S 3 (W ) = 0, T xz , T yz , T zz , uT xz + vT yz + wT zz + λT z , ρ σ (ν + ν)ν z T , (4)
where T ij are the components of laminar stress tensor T = (µ + µ t ) (∇ ⊗ u + t ∇ ⊗ u) -2 3 ∇. u I :

T xx = (µ + µ t ) 2 3 (2u x -v y -w z ) , T xy = (µ + µ t ) (u y + v x ) , T xz = (µ + µ t ) (u z + w x ) , ... Q(W
) are the source terms, i.e. the diffusion, production and destruction terms from the Spalart-Allmaras turbulence model:

Q(W ) = (0, 0, 0, 0, 0, c b2 ρ σ ∇ν 2 + ρc b1 S ν + c w1 f w ρ ν d 2 ) T . (5)
Note that Q = 0 in the case of the laminar Navier-Stokes equations, unless additional source terms are added (to take into account gravity, for instance).

Spatial discretization

The spatial discretization of the fluid equations ( 1) and ( 2) is based on a vertex-centered finite element/finite volume formulation on unstructured meshes. The equations are integrated on each finite volume cell C i (using the Green formula):

|C i | dW i dt + F i = S i + Q i , (6)
where W i is the mean value of the solution W on cell C i , F i , S i and Q i are respectively the numerical convective, viscous and source flux terms:

F i = ∂Ci F (W i ) • n i dγ , S i = ∂Ci S(W i ) • n i dγ , Q i = Ci Q(W i ) dΩ ,
where n i is the outer normal to the finite volume cell surface ∂C i , and F , S and Q are respectively the convective, viscous and source terms flux functions as defined previously in Relations (3), ( 4) and (5). The numerical scheme combines a HLLC approximate Riemann solver 10 to compute the convective fluxes and the Galerkin centered method to evaluate the viscous terms. Second order space accuracy is achieved through a piecewise linear extrapolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) procedure [START_REF] Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF] which uses a particular edge-based formulation with upwind elements.

HLLC approximate Riemann solver. The idea of the HLLC flow solver is to consider locally a simplified Riemann problem with two intermediate states depending on the local left and right states. The simplified solution to the Riemann problem consists of a contact wave with a velocity S M and two acoustic waves, which may be either shocks or expansion fans. The acoustic waves have the smallest and the largest velocities (S i and S j , respectively) of all the waves present in the exact solution. If S i > 0 then the flow is supersonic from left to right and the upwind flux is simply defined from F (W i ) where W i is the state to the left of the discontinuity. Similarly, if S j < 0 then the flow is supersonic from right to left and the flux is defined from F (W j ) where W j is the state to the right of the discontinuity. In the more difficult subsonic case when S i < 0 < S j we have to calculate F (W * i ) or F (W * j ). Consequently, the HLLC flux is given by:

Φ HLLC ij (W i , W j , n ij ) =          F (W i ) • n ij if S i > 0 F (W i ) • n ij if S i ≤ 0 < S M F (W j ) • n ij if S M ≤ 0 ≤ S j F (W j ) • n ij if S j < 0 .
W * i and W * j are evaluated as follows. We denote by η = u • n. Assuming that η * = η * i = η * j = S M , the following evaluations are proposed 10 (the subscripts i and j are omitted for clarity):

W * = 1 S -S M    ρ (S -η) ρu (S -η) + (p * -p)n ρE (S -η) + p * S M -pη    where p * = ρ (S -η)(S M -η) + p .
A key feature of this solver is in the definition of the three waves velocity. For the contact wave we consider:

S M = ρ j η j (S J -η j ) -ρ i η i (S I -η i ) + p i -p j ρ j (S J -η j ) -ρ i (S I -η i ) ,
and the acoustic wave speeds based on the Roe average:

S I = min(η i -c i , η -c) and S J = max(η j + c j , η + c) .
With such waves velocities, the approximate HLLC Riemann solver has the following properties. It automatically (i) satisfies the entropy inequality, (ii) resolves isolated contacts exactly, (iii) resolves isolated shocks exactly, and (iv) preserves positivity.

2nd-order accurate version. The MUSCL type reconstruction method has been designed to increase the order of accuracy of the scheme. [START_REF] Leer | Towards the ultimate conservative difference scheme I. The quest of monotonicity[END_REF] The idea is to use extrapolated values W ij and W ji instead of W i and W j at the interface ∂C ij to evaluate the flux with the approximate Riemann solver. Note that, in the implementation, the primitive variables (ρ, u, p) are extrapolated to guarantee the positivity of the density and the pressure, then the conservative variables are reconstructed from these values. Thus, the gradients of the primitive variables are evaluated. However, in the following, we still denote by W the primitive variables vector. The numerical flux becomes:

Φ ij = Φ HLLC ij (W ij , W ji , n ij ) ,
where W ij and W ji are linearly extrapolated as:

W ij = W i + 1 2 (∇W ) i • --→ P i P j and W ji = W j + 1 2 (∇W ) j • --→ P j P i .
In contrast to the original MUSCL approach, the approximate "slopes" (∇W ) ij and (∇W ) ji are defined for any edge and obtained using a combination of centered, upwind and nodal gradients. The centered gradient, which is related to edge P i P j , is implicitly defined along edge P i P j by the relation:

(∇W ) C ij • --→ P i P j = W j -W i and (∇W ) C ij • --→ P j P i = W i -W j .
Upwind and downwind gradients, which are also related to edge P i P j , are computed according to the definition of upwind and downwind tetrahedra of edge P i P j . These tetrahedra are respectively denoted K ij and K ji . K ij (resp. K ji ) is the unique tetrahedron of the ball of P i (resp. P j ) the opposite face of which is crossed by the line defined by the edge P i P j , see Figure 1. Upwind and downwind gradients are then defined for vertices P i and P j as:

(∇W ) U ij = (∇W )| Kij and (∇W ) D ij = (∇W )| Kji .
where (∇W )| K = P ∈K W P ∇φ P | K is the P 1 -Galerkin gradient on tetrahedron K. Parametrized gradients are built by introducing the β-scheme:

(∇W ) i • --→ P i P j = (1 -β)(∇W ) C ij • --→ P i P j + β (∇W ) U ij • --→ P i P j (∇W ) j • --→ P j P i = (1 -β)(∇W ) C ij • --→ P j P i + β (∇W ) D ij • --→ P j P i ,
where β ∈ [0, 1] is a parameter controlling the amount of upwinding. For instance, the scheme is centered for β = 0 and fully upwind for β = 1.

Fourth-order numerical dissipation: V4-scheme. The most accurate β-scheme is obtained for β = 1/3. Indeed, it can be demonstrated that this scheme is third-order for the two-dimensional linear advection on structured triangular meshes. On unstructured meshes, a second-order scheme with a fourthorder numerical dissipation is obtained. These high-order gradients are given by:

(∇W ) V 4 i • --→ P i P j = 2 3 (∇W ) C ij • --→ P i P j + 1 3 (∇W ) U ij • --→ P i P j (∇W ) V 4 j • --→ P j P i = 2 3 (∇W ) C ij • --→ P j P i + 1 3 (∇W ) D ij • --→ P j P i .
The MUSCL schemes are not monotone and can be a source of spurious oscillations especially in the vicinity of discontinuities. [START_REF] Cournède | Positivity statements for a Mixed-Element-Volume scheme on fixed and moving grids[END_REF] These oscillations can affect the accuracy of the final solution or simply end the computation because (for instance) of negative pressures. A widely used technique for addressing this issue is to guarantee the TVD property in 1D [START_REF] Harten | On upstream differencing and godunov-type schemes for hyperbolic conservation laws[END_REF] or the LED property in 2D/3D of the scheme, which ensures that the extrapolated values W ij and W ji are not invalid. To guarantee the TVD or the LED properties, limiting functions are coupled with the previous high-order gradient evaluations. The gradient is substituted by a limited gradient denoted (∇W ) lim ij . The choice of the limiting function is crucial as it directly affects the convergence of the simulation. In this work, we use the Piperno limiters [START_REF] Piperno | Criteria for the design of limiters yielding efficient high resolution TVD schemes[END_REF][START_REF] Piperno | Schéma TVD d'ordre élevé pour la résolution de l'équation de Burgers[END_REF] which is expressed in a factorized form,

(∇W ) Lim = ∇W C ψ P I ∇W C ∇W V 4 ,
with

ψ P I (R) = ( 1 3 + 2 3 R)        3 1 R 2 -6 1 R + 19 1 R 3 -3 1 R + 18 if R < 1 1 + ( 3 2 1 R + 1)( 1 R -1) 3 if R ≥ 1 , where R = ∇W C ∇W V 4 .

B. Adjoint based correction

Linear correction:

Despite some differences, the usual adjoint based approaches relies on the following implicit correction: [START_REF] Fidkowski | Review of Output-Based Error Estimation and Mesh Adaptation in Computational Fluid Dynamics[END_REF][START_REF] Venditti | Adjoint Error Estimation and Grid Adaptation for Functional Outputs: Application to Quasi-One-Dimensional Flow[END_REF][START_REF] Venditti | Grid adaptation for functional outputs: application to two-dimensional inviscid flows[END_REF][START_REF] Formaggia | Anisotropic error estimates for elliptic problems[END_REF][START_REF] Jones | Validation of 3D Adjoint Based Error Estimation and Mesh Adaptation for Sonic Boom Reduction[END_REF][START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] given an analytical problem defined as Ψ(W ) = 0, a numerical discretization of the problem is determined (using here finite volume approach)

Ψ h (W h ) = 0.
From the numerical solution an output functional j h (W h ) is computed. Assuming the numerical solution W h is close to the analytical solution W , the functional j can be linearize as

j h (Π h W ) ≈ j h (W h ) + J • (Π h W -W h ),
where Π h W stands for the projection of the exact solution on the mesh and J = ∂j ∂W is the gradient of j with respect to the solution. The numerical equation can also be linearized as in the Newton's method as

Ψ h (Π h W ) ≈ Ψ h (W h ) + A • (Π h W -W h ), where A = ∂Ψ h
∂W is the jacobian of Ψ h . Put together we obtain

j h (Π h W ) ≈ j h (W h ) + W * • Ψ h (Π h W ),
where W * = J • A -1 = A -T • J is the adjoint of the problem. This formulation is convenient for mesh adaptation as W * can be precomputed separately and Ψ(Π h W ) can be estimated with local interpolation errors. However, we can see that intrinsically this is equivalent to compute a corrected solution as

W = W h + A -1 Ψ h (Π h W )
and then compute the functional j

j h ( W ) ≈ j h (W h ) + J • ( W -W h ),
so that the correction on itself is linear. Similarly, error transport equations are based on the linearization of the given problem: [START_REF] Derlaga | Application of Exact Error Transport Equations and Adjoint Error Estimation to AIAA Workshops[END_REF][START_REF] Hay | Error estimation using the error transport equation for finite-volume methods and arbitrary meshes[END_REF][START_REF] Layton | A defect-correction method for the incompressible Navier-Stokes equations[END_REF][START_REF] Pierce | Adjoint and defect error bounding and correction for functional estimates[END_REF] a linear correction problem is deduced from the initial problem and thus provide a linear correction of the solution. The linearization of the problem is thus a determining factor of the quality of the corrected solution. [START_REF] Yan | Accuracy of Discretization Error Estimation by the Error Transport Equation on Unstructured Meshes -Nonlinear Systems of Equations[END_REF]

Nonlinear correction:

This is why we propose a nonlinear approach for finite volume solvers based on the fact that

Ψ h (Π h W ) = ε h = 0 and Ψ(W h ) = ε = 0.
In the finite volume approach, the residual is defined as the integral of the continuous residual of a reconstructed solution over a cell or equivalently for inviscid flow as the integral of the fluxes through the boundary of the cell. In that sense, for a given mesh an infinite family of continuous fields W verifies Ψ h ( W ) = 0, the reconstructed continuous field based on the discrete solution is one of them. However, if we consider finer meshes, in particular the n th subdivisions of the initial mesh, the only field verifying Ψ h/n ( W ) = 0 for any submesh is the solution of the continuous problem: W = W . Thus, the residual of the discrete solution interpolated on a finer mesh Ψ h/2 (I h/2 W h ) = ε h/2 gives us a usefull indication of the quality of the solution.

We propose the following approach, inspired by multigrid methods, to compute a corrected solution W h from the residual Ψ h/2 (I h/2 W h ). The mesh is initially divided by two (or a power of two) and the solution W h is linearly interpolated on it. The residual is then computed on the finer mesh and accumulated back on the initial mesh.

The linear accumulation process can be seen as an inverse operation of the linear interpolation: The interpolated value of a vertex of a finer mesh is a linear combination of the values of some vertices of the coarse mesh. In the accumulation process, the residual of this vertex of the fine mesh will be added to the contribution of the vertices of the coarse mesh weighted by the same coefficients. This forms a source term

S h = I h/2→h Ψ h/2 I h→h/2 (W h ) ,
which is then introduced in the equation to compute the correction as

Ψ h ( W h ) = S h .
The corrected solution W h is computed performing a few iterations with the same solver, with the source term S h , which nonlinearly propagates the errors. This approach automatically encompasses all features of the flow solver (order, flux reconstruction, limiters, ... ) both in the flux computation and the resolution. All of the operations involved here are local, so that we don't actually need to generate the actual h/2 mesh, which would be far too expensive. Instead, we extract for each point its vicinity mesh and divide it by two to compute and accumulate the source terms fluxes. This approach yields no memory over cost and can be done efficiently in parallel. Finer subdivision (h/4, h/8, ...) of the initial mesh can also be used to improve predictions. The overall procedure is thus:

• Compute the solution on the initial mesh: W h

• Linearly interpolate the solution on the h/2 mesh:

Π h/2 W h • Compute the h/2 residual: Ψ h/2 (I h/2 W h ) = ε h/2
• Accumulate the h/2 residual on the h mesh:

S h = A h/2 h (ε h/2 )
• Perform a few iterations of the flow solver with S h as source term: Ψ h ( W h ) = S h .

Numerical Implementation: Dividing the element size by two multiplies the number of elements by 4 in 2D and 8 in 3D, it is thus particularly expensive to fully generate the submeshes for large practical cases. As the residual computation and accumulation is a local operation, it can be done at a local level, by generating virtually the submeshes of the ball around each vertex. That way no additional memory cost is implied and the process can straightforwardly done in parallel. Similarly, multiple computations of the same element can be avoided with a careful ordering detailed here after. Figure 2 shows an example of 2D mesh and it's subdivision by two. Starting with the solution on the coarse mesh (circle dots and dark lines), the solution is interpolated on the divided mesh (grey lines and square dots). The residual is computed on the fine mesh and accumulated on the coarse mesh (weighted by the barycentric coordinates).

A careful analysis of the numerical discretization presented in Section II.A shows that the computation of the residual of a given vertex requires the values of the vertices of its third order ball (neighbours of the neighbours of its neighbours). Similarly, to compute the accumulated source term of a given vertex, we only need the vertices of the fine mesh encompassed by the first oder ball of the given vertex in the coarse mesh. This is why, we can generate separately the subdivided mesh of the second order ball for each vertex of the coarse mesh to compute its source term.

However, it can be reduced to the subdivision of the first order ball if we consider carefully the different terms involved in the computation of the residual. In our vertex centered finite volume solver Wolf , the advection fluxes are computed and summed for each edge, while the viscous terms are computed for each triangle. We can thus compute them separately in any order as explained below.

We can see in Figure 2 that the residual of the fine node C contributes through the accumulation process to the source term of coarse node A and B by a factor 1/2. The residual of node D contributes to the source term of node A but not of node B (barycentric coordinate 0), however, its solution is required to compute the residual of node C and thus contributes indirectly to the source term of B. Finally, node E directly contributes to the source term of A, but it belongs to an upwind element of the edge DC, so that it contributes to the residual of C and the source term of B.

If we want to compute the source term of node B, the residual of node D does not need to be computed, and the residual of node C which contributes to both source terms of A and B can be computed only once. To do so, we can notice that the advection term in the residual of C is the sum of the fluxes across each edge connecting C to one of its neighbours, i.e. here the edges CA, CD, CB, ... As the residual of C is added to the source terms of nodes A and B in the accumulation, we can compute and add separately the contributions of edges CA and BC.

To this end, we generate the local sub-mesh as shown in Figure 3. The grey part of the sub-mesh is not needed and not generated at all, the black part of the sub-mesh is generated for the interpolation of the nodal values and the fluxes are computed only on the red edges and triangles encompassed. That way, only a half of the residual of node C is computed and accumulated on both nodes A and B and the other half of the residual is computed in the vicinity of A. Additionally, the central sub-triangle contribution is computed in the vicinity of the node with the smallest ID in the coarse triangle. Each contribution of each fine edge and triangle are thus computed only once which reduces computational costs. It also ensures that each edge has an upwind and downwind triangles for second order extrapolations in the HLLC flux solver (see Section II.A.4).

III. Numerical Results

In the sequel, given the analytical solution of a problem W and the numerical solution W h , we split the approximation error W -W h into the interpolation error and the implicit error as

W -W h = W -Π h W Interpolation Error + Π h W -W H Implicit Error .
The interpolation error is entirely related to the discretization of the solution and can only be improved by a refinement of the mesh. Here, we are interested in the implicit error as we want to improve the discrete value of our solution. We can thus compute analytically the exact implicit error: the estimated implicit error:

E Imp Exa = 1 0 Π h W -W h 2 dΩ,
E Imp Est = 1 0 W h -W h 2 dΩ,
as they can be expressed element-wise as polynomial expressions. Similarly, we define a function j computed from u, and we define the exact functional error:

E F un Exa = |j(W h ) -j(Π h W )|, the estimated functional error: E F un Est = |j(W h ) -j( W h )|
, and the corrected functional error:

E F un Cor = |j( W h ) -j(Π h W )|.
Note that the interpolation error is not taken into account here as, we compare the numerical functional j(W h ) to j(Π h W ) which is an approximation of the exact j(W ). However this interpolation error can be taken into account in another means. In the sequel, we will use a finer solution as estimation of the exact solution.

A. 2D Laminar NACA0012

We first validate the present corrector in an adaptive context on a NACA0012 geometry, at a Reynolds number of 500 with a purely laminar viscous flow. The Mach number is fixed to M = 0.1 and the angle of attack to α = 3 • . Adapted meshes are generated iteratively with feflo.a 21, 22 using a goal oriented error estimate on the drag coefficient. [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] An example of the resulting meshes and solutions is shown in Figures 4 and5.

Figure 8 shows the convergence of the drag computed for both the initial (red) and corrected (blue) solution. We can see that corrector moderately improves the solution and the computation of the drag coefficient. The corrected solution on a 5 × 10 3 nodes grid provides the same result as the solution on a 8 × 10 3 nodes grid. The moderate improvement is due to, the mesh adaptation context that introduces a side effect: an adapted mesh with four times the number of vertices of a given adapted mesh is not its subdivision by two. Indeed, it has been further improved by the adaptation and provides a better solution than a subdivision. We will try to quantify this effect with the next test case.

Still, the primal and corrected solutions converge to the same value, so that the difference between both of them gives and accurate estimation of the remaining numerical error. This second property is particularly important in both adaptative and fixed grids context as it informs on how confident we can be in a given solution. It also indicates if the grid convergence study should be continued or not.

B. 2D Turbulent NACA0012

The corrector is then validated on a turbulent case with the same geometry, at a Reynolds number of 5 × 10 5 with the Sapalart-Allmaras turbulence model. The Mach number is fixed to M = 0.7 and the angle of attack to α = 1 • . Adapted meshes are generated iteratively with feflo.a using a hessian-based error estimate based on the Mach number field. An example of the resulting meshes and solutions is shown in Figures 6 and7. Additionally, we generate for each mesh its h/2 subdivision and compute the solution on it for comparison.

Figures 9 and10 show the convergence of the drag computed with a first and second order scheme for both the initial (red) and corrected (black) solution. The blue curve represents for each mesh the drag computed on its h/2 mesh subdivision which actual complexity is thus four times larger. Both curves exhibits the same trend as the previous laminar example, but it is now possible to quantify how the mesh adaptation improves the solution on a mesh with four times more vertices compared to the actual h/2-mesh subdivison. This dims the apparent efficiency of the non-linear corrector. In particular we can see that the corrector is very efficient with a first order scheme as it yields the same solution as the h/2-mesh. It is slightly less efficient with a second order scheme but still satisfying.

C. 3D inviscid supersonic case: Sonic Boom Prediction Workshop

The last example is a supersonic case on a realistic aircraft configuration. The considered geometry was proposed for the 2 nd AIAA Sonic-Boom Prediction Workshop which aims at minimizing the noise produced by supersonic aircrafts on the ground. Different geometries, meshes and configurations were proposed, we focus here on the C25D flow through Euler case, see Figures 11 and12. In order to properly lead the shape optimisation process, the noise signature under the aircraft has to be accurate, this is why a grid convergence Figure 10. Convergence of the drag with a second order scheme on NACA0012 for the initial solution (red), the corrected solution (black) and the solution on the mesh divided by two (blue) at Reynolds number of 5 × 10 5 .

study is required to determine the numerical error. To this end, mesh adaptation has proven to be very efficient [START_REF] Alauzet | High Order Sonic Boom Modeling by Adaptive Methods[END_REF] and we expect the corrector to provide a second pertinent indication on the convergence of the solution.

Figures 13, 15 and 17 show the pressure signatures (pressure difference with the free stream pressure normalized) under the plane computed on tailored grids provided for the workshop. Error bars represent the estimated error of each nodal value computed with the nonlinear corrector. We can see that the corrector indicates that the error level remains relatively high, even on the 13 million grid, which has been confirmed by the workshop: the grid convergence is not achieved. This is visible on the leading shock which is still relatively smooth and the shocks in the aft part of the signature (65 < X < 80) which are still forming. Note that the actual subdivision by two of this 13 millions nodes mesh would have yielded a 104 millions nodes mesh. Figures 14, 16 and 18 shows the same computations using goal-oriented mesh adaptation which have been generated to optimally capture the pressure signature. [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] We can se that the error bars are much smaller and tend to reduce during the adaptation process. This is especially visible on the mid part of the signature (50 < X < 60). It is even more enlightening to compare the error levels predicted by the corrector on the tailored grid composed of 13 million nodes in the mid part of the signature to the actual difference between the solutions on the tailored grid and adapted grid composed of 5.9 million nodes. The error predicted by the corrector on tailored grid matches the difference between the solution on the tailored grid (which is not converged) and the solution on the adapted grid (which is converged in this region).

D. 3D ONERA M6 Wing turbulent Case

We show here the improvements this method can provide on fixed grids in 3D. We choose to use here the ONERA M6 Wing configuration at Mach number M = 0.1, angle of attack α = 3.06 • and Reynolds number of Re = 11 × 10 6 . A series of four grids were generated using AFLR, [START_REF] Marcum | Adaptive unstructured grid generation for viscous flow applications[END_REF][START_REF] Marcum | Unstructured Grid Generation Using Automatic Point Insertion and Local Reconnection[END_REF] with a dimensionless wall spacing based on the Reynolds number of y + = 3, 2, 1, 0.5. The coarsest mesh is shown in Figure 20.

We can see in Figure 19 that here, the corrector significantly improves the drag prediction, it yields on the 5 × 10 5 nodes mesh the same drag as the solution without correction on the 3 × 10 6 nodes mesh. The corrector is even more efficient when the error is not optimally reduced between two grids. We also have a good estimation of the quality of our last solution, without the need of any finer mesh: despite the fine resolution of the boundary layer, we can be pretty confident on the fact that the solution is not converged yet. This is mainly due to the fact that the wake of the wing is poorly discretized, as can be seen in Figure 20. This pollutes the computation of the drag.

IV. Conclusion

We have seen that the nonlinear correction for inviscid, laminar and RANS flow solutions already shows promising results. It improves the overall solution and output functionals, and thus provides a reliable estimation of the numerical error. As it does not requires to effectively generate a finer mesh, it is relatively cheap to use and can provide a pertinent error estimation even on a single grid. It is particularly efficient on fixed grids and can be used to asses the quality of the mesh refinement that has been done. In a mesh adaptation context its apparent efficiency is reduced as the mesh adaptation process efficiently reduces the error. Though, it is still a useful tool to asses the convergence of the overall process. However, the transonic NACA0012 case showed that the correction is less efficient with a second order scheme than with a first order scheme, suggesting that the MUSCL extrapolation requires an additional particular treatment.

Figure 1 .

 1 Figure 1. Downwind K ij and upwind K ji tetrahedra associated with edge P i P j .

Figure 2 .Figure 3 .

 23 Figure 2. Mesh Subdivision

Figure 4 .

 4 Figure 4. Adapted mesh generated using goaloriented error estimate on the drag for NACA0012 at Re = 500, M = 0.1 and α = 3 • : 8009 vertices.

Figure 5 .

 5 Figure 5. Velocity contours around a NACA0012 at Re = 500, M = 0.1 and α = 3 • .

Figure 6 .

 6 Figure 6. Adapted mesh generated using hessianbased error estimate on the Mach number field for NACA0012 at Re = 5 × 10 5 M = 0.7 and α = 1 • : 5394 vertices .

Figure 7 .

 7 Figure 7. Velocity contours around a NACA0012 at Re = 5 × 10 5 , M = 0.7 and α = 1 • .

Figure 8 .

 8 Figure 8. Convergence of the drag on NACA0012 for initial solution (red) and corrected solution (blue) at Reynolds number of 500.

Figure 9 .

 9 Figure 9. Convergence of the drag with a first order scheme on NACA0012 for the initial solution (red), the corrected solution (black) and the solution on the mesh divided by two (blue) at Reynolds number of 5 × 10 5 .

Figure 11 .

 11 Figure 11. C25D Geometry of the 2 nd Sonic-Boom Prediction Workshop.

Figure 12 .

 12 Figure 12. Cut in the volume of the goal-oriented adapted mesh and solution on the C25D geometry.

Figure 13 .

 13 Figure 13. Pressure difference under the plane for tailored grid of 3.4 million vertices (red line) with estimated error provided by nonlinear corrector (black error bars).

Figure 14 .

 14 Figure 14. Pressure difference under the plane for goal-oriented adapted grid of 1.5 million vertices (red line) with estimated error provided by nonlinear corrector (black error bars).

Figure 15 .

 15 Figure 15. Pressure difference under the plane for tailored grid of 6.3 million vertices (red line) with estimated error provided by nonlinear corrector (black error bars).

Figure 16 .

 16 Figure 16. Pressure difference under the plane for goal-oriented adapted grid of 3.1 million vertices (red line) with estimated error provided by nonlinear corrector (black error bars).

Figure 17 .

 17 Figure 17. Pressure difference under the plane for tailored grid of 13 million vertices (red line) with estimated error provided by nonlinear corrector (black error bars).

Figure 18 .

 18 Figure 18. Pressure difference under the plane for goal-oriented adapted grid of 5.9 million vertices (red line) with estimated error provided by nonlinear corrector (black error bars).

Figure 19 .

 19 Figure 19. Convergence of the drag on ONERA M6 Wing for initial solution (red) and corrected solution (green) at Reynolds number of 11 × 10 6 and angle of attack α = 3.06 • .

Figure 20 .

 20 Figure 20. Coarsest M6 Wing mesh generated with a fixed structured boundary layer and frontal approach with AFLR: 110009 vertices.

  T xx , T xy , T xz , uT xx + vT xy + wT xz + λT x , T xy , T yy , T yz , uT xy + vT yy + wT yz + λT y ,

		ρuw, ρvw, ρw 2 + p, w(ρE + p), ρwν	T .
	S(W ) = (S 1 (W ), S 2 (W ), S 3 (W )) are the laminar viscous fluxes:
	S 1 (W ) =	0, ρ σ	(ν + ν)ν x	T	,
	S 2 (W ) =	0, ρ σ	(ν + ν)ν y	T	,


Angélique Villière

Catherine Fillonneau

Carole Prost

published or not. The documents may come

InnOscent chromatographic system : a new potent device for wine aroma analysis

InnOscent chromatographic system : a new potent device for wine aroma analysis Angélique VILLIÈRE Catherine FILLONNEAU Carole PROST A new chromatographic device was developed to go further in aroma analysis. Its innovative configuration not only permits to realize classical analysis by gas chromatography coupled with multiple detectors such as mass spectrometry, flame ionization detector and single or dual olfactometric port but also incorporates preparative fraction collection possibilities to recover bespoke aromatic samples. Accordingly, it becomes possible to evaluate the contribution of any compound to the overall aroma of a product. These experiments enable to point out 33 odorant zones among which a majority were associated to identified compounds. Despite their modest olfactometric detection, the contribution of IPMP and IBMP on the typical pyrazic-green odor of wines made of cabernet franc vineyard was evidenced.

Flexible multi detection/recovery system

Collection of designed extracts

InnOscent device in an innovative selective recovering configuration

Capillary column

These results exemplify the suitability of the approach to improve knowledge of mixture of odorants in a wine aroma context. The extracts are perceived as different (triangular test p<0,05)

The extract without pyrazines is perceived as less "green, pyrazic" and "chemical" than the total extract (p<0,05) Bespoke extracts were recovered to be sensory evaluated (scale 0-10).

Oniris -College of Veterinary Medicine, Food Science and Engineering Nantes, France Flavor team -GEPEA -UMR CNRS 6144

angelique.villiere@oniris-nantes.fr

  potentialities of the InnOscent system is presented through the aromatic analysis of a wine made of Cabernet franc vineyard after extraction of volatile compounds by a Car/PDMS SPME fiber. A particular focus was given to the contribution of pyrazines to the wine aroma.

  Vino Analytica Scientia Symposium -Trento (Italy) -14-17 July 2015