A Belme
email: belme@dalembert.upmc.fr

F Alauzet

A Dervieux

An a priori anisotropic Goal-Oriented Error Estimate for Viscous Compressible Flow and Application to Mesh Adaptation

Keywords: Viscous compressible flow, goal-oriented mesh adaptation, anisotropic mesh adaptation, adjoint, metric

An a priori anisotropic Goal-Oriented Error Estimate for Viscous Compressible Flow and Application to Mesh Adaptation

A. Belme, F. Alauzet, A.

Introduction

Mesh adaptation progressively plays an increasing role in high-fidelity simulation. Beside the direct quantitative expectations, like an increased accuracy, and an higher e ciency for a given error, an important motivation in the extension of mesh adaptation is the better control of error convergence, for example by a better numerical convergence order, since anisotropic mesh adaptation is known as yielding such progress, see [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF][START_REF] Yano | The importance of mesh adaptation for higher-order discretizations of aerodynamics flows[END_REF].

The case of Navier-Stokes flows takes a particular place since, due to well identified and localized boundary layers, the researcher and the engineer have early managed for concentrating and as far as possible stretching the mesh near the no-slip boundary. However, already for moderate Reynolds number, the automatic mesh adaptation for Navier-Stokes flows, while not a crucial issue, remains a non-trivial task. A major di culty is the choice of an error measure.

Many researchers have chosen to focalize on the interpolation error committed on the unknown(s) or on several userprescribed sensors depending on the unknowns. Mesh adaptation based on P 1 interpolation error estimate is referred as Hessianbased mesh adaptation. Most works tend to equidistribute the error which consists in minimizing the maximum of the interpolation error. Pioneering works have shown a fertile development of Hessian-based and metric-based methods [START_REF] Belhamadia | Three-dimensional anisotropic mesh adaptation for phase change problems[END_REF][START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF][START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent CFD[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF][START_REF] Gruau | 3D tetrahedral, unstructured and anisotropic mesh generation with adaptation to natural and multidomain metric[END_REF][START_REF] Hecht | Mesh adaptation by metric control for multi-scale phenomena and turbulence[END_REF][START_REF] Li | 3D anisotropic mesh adaptation by mesh modification[END_REF][START_REF] Michal | Comparison of fixed and adaptive unstructured grid results for drag prediction workshop 6[END_REF][START_REF] Pain | Tetrahedral mesh optimisation and adaptivity for steady-state and transient finite element calculations[END_REF][START_REF] Tam | Anisotropic mesh adaptation for 3D flows on structured and unstructured grids[END_REF]. In contrast to these equidistribution methods, the "multiscale" variant relies on the optimization of a L p norm of the interpolation error [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF][START_REF] Loseille | Optimal 3D highly anisotropic mesh adaptation based on the continuous mesh framework[END_REF]. The L p formulation allows the mesh adaptive process to approximate discontinuous solutions with higher-order convergence [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF]. However, these methods are limited to the minimization of some interpolation errors for some solution fields, the "sensors". If for many applications, this simplifying standpoint is an advantage, there are also many applications where Hessian-based mesh adaptation is far from optimal regarding the way the degrees of freedom are distributed in the computational domain. Indeed, Hessian-based methods aim at controlling the interpolation error but this purpose is not often so close to the objective that consists in obtaining the best approximate solution of the PDE. Further, in many engineering applications, one (or several) specific scalar output needs to be accurately evaluated, e.g. lift, drag, or heat flux. Hessian-based mesh adaptation methods are not designed to address this issue.

A generation of a posteriori error estimates has renewed the existing answers to that issue. A posteriori error estimates take the generic form:

H 7 ! E(u h ) ⇡ E exact (u u h ) ,
where u h = u h (H) is the discrete solution obtained using mesh H. In the Dual Weighted Residual (DWR) method (see [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF]), an important element is introduced, the functional error, which carries a mathematical formulation of the purpose of the mesh adaptation e↵ort. Let u, u h , u ⇤ be the state, the discrete state, and the adjoint state such that a(u, ') = ( f, ') ; a(u h , ' h ) = ( f, ' h ) ; a( , u ⇤ ) = (g, ) ,

the DWR method starts with:

E exact (u u h ) = (g, u u h ) = ( f, u ⇤ ⇧ h u ⇤ ) a(u h , u ⇤ ⇧ h u ⇤ ) ,
where ⇧ h is the P 1 interpolation operator on mesh H. We observe that the right-hand side is zero if u h solves the continuous equation, or if u ⇤ is a ne, or both. By decomposing the integrals over mesh elements and applying a Green formula, the DWR method identifies the local error which the mesh update will attempt to make uniform. Another way is to evaluate the local error with a finer mesh or by increasing the order of the numerical scheme. These goal-oriented methods have been developed in a series of papers, see e.g. [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF][START_REF] Giles | Adjoint methods for PDEs: a posteriori error analysis and postprocessing by duality[END_REF][START_REF] Jones | Validation of 3D adjoint based error estimation and mesh adaptation for sonic boom reduction[END_REF][START_REF] Leicht | Error estimation and anisotropic mesh refinement for 3d laminar aerodynamic flow simulations[END_REF][START_REF] Park | Adjoint-based, three-dimensional error prediction and grid adaptation[END_REF][START_REF] Venditti | Anisotropic grid adaptation for functional outputs: application to two-dimensional viscous flows[END_REF][START_REF] Verfürth | A review of A Posteriori Error Estimation and Adaptative Mesh-Refinement techniques[END_REF][START_REF] Wintzer | Adjoint-based adaptive mesh refinement for sonic boom prediction[END_REF][START_REF] Yano | An optimization-based framework for anisotropic simplex mesh adaptation[END_REF].

However, the initial form of the DWR method does not give a good access to anisotropic mesh adaptation. This di culty has been circumvented by combining the DWR method and interpolation error criteria as in [START_REF] Venditti | Anisotropic grid adaptation for functional outputs: application to two-dimensional viscous flows[END_REF]. It can also be solved by introducing local deformation maps as in [START_REF] Formaggia | New anisotropic a priori error estimates[END_REF][START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-di↵usion-reaction and the Stokes problems[END_REF]. It is also addressed in [START_REF] Yano | An optimization-based framework for anisotropic simplex mesh adaptation[END_REF], where it is proposed to perform a "what if" study which reevaluate locally the error after an anisotropic division has been applied to an element. The resulting information is then projected into an anisotropic refinement criterion involving a stretching direction and anisotropy strength. Such approaches are only locally optimal, they do not provide a global optimum. As a consequence, the adaptive process may require a large number of iterations.

Our proposal relies on a priori estimates. An a posteriori estimate depends on the mesh through the discrete solution. In contrast, an a priori estimate is expressed directly as a function of the mesh and of the continuous solution u. Typically, it writes:

H 7 ! E(u, H) ⇡ E exact (u u h ).
A blemish of a priori error estimates is the fact that the main ingredients in right hand-side are high-order derivative (the order is higher than or equal to the order of the equation) of the continuous solution which is unknown. A typical example of mesh adaptation method based on a priori error estimate is the very popular family of Hessian-based methods, also called feature-based mesh adaptation, using estimates such as:

ku u h k L 2  C h 2 kuk H 2 .
Many strategies have been developed to evaluate these second derivatives from the discrete solution u h such as the L 2 projection method [START_REF] Ph | Approximation by finite element functions using local regularization[END_REF][START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF], least-square methods, patch recovery process [START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 1: The recovery technique[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. Part 2: Error estimates and adaptivity[END_REF], hierarchical approaches [START_REF] Bank | A posteriori error estimate based on hierarchical bases[END_REF], ... All these approaches are equivalent to reconstruct a higher-order solution from the piecewise linear discrete solution. Mathematically Hessian reconstruction methods converge at zeroth order, nevertheless in practice super-convergence appears and a higher rate of convergence is obtained. Even if we can find very specific exemple where some methods do not converge, these approaches are very robust and have been intensively validated in the context of mesh adaptation by the numerous studies that have been published on the subject.

In compensation, these techniques have some attracting features. Among them, we remark that for a scalar error (goal oriented formulation), the mesh adaptation problem can be expressed as an inverse problem (u fixed):

Find H opt such that E(u, H opt ) = min H (E(u, H)).

(

) 1
In Problem [START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF], it is mandatory that the "min" is sought in a set of meshes enjoying some compactness in order to have a wellposed inverse problem. At least, the number of degrees of freedom of the meshes of this set has to be bounded. But, a second important feature of a priori estimates is that we can specify in a more detailed way the class of considered meshes, imposing some regularity, which will lead to a more accurate error analysis. Lastly, the right-hand side of the a priori estimate which we propose here is somewhat dual to the DWR one. Indeed, it writes:

E exact (u u h ) = (g, u u h ) = a(⇧ h u u, u ⇤ h ) + (g, u ⇧ h u) , (2)
where u ⇤ h is the discrete adjoint state. We observe that the error is zero for an a ne exact solution: this is the P 1 -exactness. In contrast to the DWR method, this formulation gives the priority to the reduction of the interpolation error for the state variables.

Right hand-side of Relation (2) involves derivatives of ⇧ h u u leading to a H 1 analysis of the interpolation error which is much more complex. However, due to the scalar (i.e., integral) form of the functional, integration by parts combined with a smoothness assumption for the adjoint solves this obstacle. Indeed, using integration by part, derivatives can be transposed onto the adjoint state to obtain interpolation errors in L 1 norm weighted by derivatives of the adjoint state. It is then possible to apply the interpolation error theory to get a majoration of the error. It is also possible to consider the continuous mesh framework [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] to derive the optimal mesh minimizing the considered error for a given number of vertices. This gives an optimal answer to Problem [START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF]. This stresses the central role of P k -exactness of the approximation: if the above interpolation error vanishes, then the approximation error is also zero. And, according to the standard finite element analysis, if the H 1 norm of this interpolation error is small, so is the approximation error (in the same norm).

The first contribution of this paper is to propose an adjoint-based a priori analysis for elliptic models. A preliminary formulation of this analysis was given in [START_REF] Belme | Unsteady aerodynamic and adjoint method[END_REF]. It was tested for elliptic models in [START_REF] Brèthes | A mesh-adaptive metric-based full-multigrid for the Poisson problem[END_REF]. In the present paper, we establish in a more direct manner the error estimation for the elliptic terms.

The second contribution is to extend this analysis to the complete compressible Navier-Stokes system, which involves nonlinear parabolic terms. We can demonstrate by manipulating the non-linear viscous terms that each of them can be written as a combination of an elliptic terms (on which the above mentioned error estimation applies) and higher order error terms (that can be neglected). By addressing viscous terms, this analysis complements the previous inviscid analysis performed for the compressible Euler equations [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF].

The paper is outlined as follows. Section 2 recalls the continuous mesh framework which provides a duality between discrete meshes and entities and Riemannian metric space and metric tensor. The continuous interpolation error local model is also provided. Section 3 proposes a new a priori estimate for the Poisson problem. Then, Section 4 gives the Navier-Stokes system for compressible gas and the considered variational discrete formulation. The goal-oriented error estimate for the Navier-Stokes equations is given in Section 5, and Section 6 formulates the mesh optimization problem leading to the expression of the optimal continuous mesh. Finally, Section 7 states how the optimal discrete meshes is obtained and Section 8 gives few numerical experiments illustrating the optimality of the proposed adaptation process.

Continuous mesh model

We propose to work in the continuous mesh framework introduced 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. In that context, a continuous mesh M of computational domain ⌦ ⇢ R d is identified to a Riemannian metric field (M(x)) x2⌦ . It defines proper di↵erentiable optimization [START_REF] Arsigny | Log-Euclidean metrics for fast and simple calculus on di↵usion tensors[END_REF], in other words calculus of variations can be used on continuous meshes while it cannot be applied on discrete meshes. For d = 3 and any x of ⌦, M(x) is a symmetric 3 ⇥ 3 matrix having ( i (x)) i=1,3 as eigenvalues along the principal directions R(x) = (v i (x)) i=1,3 . Metric tensor M(x) is a continuous element modeling discrete element K (triangle in 2D and tetrahedron in 3D). For i = 1, 3, sizes along principal directions v i (x) are given by h i (x) = 1 2 i (x) and anisotropy quotients r i are defined by: r i = h 3 i (h 1 h 2 h 3 ) 1 . Anisotropic quotient represent the anisotropy of the continuous element. The diagonalisation of M(x) writes:

M(x) = d 2 3 M (x) R(x) 0 B B B B B B B B B B B B @ r 2 3 1 (x) r 2 3 2 (x) r 2 3 3 (x) 1 C C C C C C C C C C C C A t R(x) . (3)
The continuous mesh density d M (x) at node x is equal to:

d M = (h 1 h 2 h 3 ) 1 = ( 1 2 3 ) 1 2 = p det(M)
. By integrating it, we define the continuous mesh complexity C(M):

C(M) = Z ⌦ d M (x) dx = Z ⌦ p det(M(x)) dx ,
which enables the user to control the level of accuracy of the mesh, and thus, to implicitly control the number of vertices of the resulting discrete mesh.

The main idea of metric-based mesh adaptation, initially introduced in [START_REF] George | Creation of internal points in Voronoi's type method. Control and adaptation[END_REF], is to generate a unit mesh in the prescribed Riemannian metric space M, e.g. a mesh of ⌦ ⇢ R d such that each edge e has a unit length and each element K -tetrahedron or triangle -is regular (or equilateral) with respect to (M(x)) x2⌦ : 8e, `M(e) = 1 and 8K,

|K| M = p 3 4 in 2D or |K| M = p 2 12 in 3D ,
where the length of edge e = ab with respect to M is computed using the straight line parameterization (t) = a + t ab, where t 2 [0, 1]:

`M(ab) = Z 1 0 k 0 (t)k M dt = Z 1 0 q ab T M(a + t ab) ab dt ,
and the volume of element K with respect to M is:

|K| M = Z K p det M(x) dx .
The resulting discrete mesh in the canonical Euclidean space will be anisotropic and adapted. We want to emphasize that the set of all the discrete meshes that are unit meshes with respect to a unique continuous mesh M contains an infinite number of discrete meshes. Given a smooth function u defined on ⌦, to each unit mesh H with respect to continuous mesh M corresponds a local discrete interpolation error |u ⇧ h u| on each element. 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], it is shown that all these interpolation errors are well represented by the continuous interpolation error related to M which is locally expressed in term of the Hessian H u of function u as follows:

|u ⇡ M u|(x) = c d trace ⇣ M 1 2 (x) |H u (x)| M 1 2 (x) ⌘ , (4)
where c d equals 1 8 in 2D and 1 20 in 3D, and |H u | is deduced from H u by taking the absolute values of its eigenvalues.

Remark 2.1. The above metric density d M is generally a smooth function representing the discrete vertex density d H of the unit mesh H with respect to M. In contrast, the density d H is a sum of Dirac located at vertices and can be therefore considered as highly oscillating. The metric density is a homogenization of the unit mesh density. Indeed, a continuous mesh convergent sequence can be defined according to the continuous mesh complexity C(M) = N and a reference metric M 1 of complexity equal to one C(M 1 ) = 1:

M N = C(M) M 1 = N M 1 , N ! 1 .
To fix the idea, we can consider an average mesh size h corresponding to complexityN, and reformulate the continuous mesh convergent sequence as:

M h = h d M 1 , h ! 0
, where d is the domain dimension. The continuous and discrete densities are close to each other in the sense that:

8 2 D(⌦), (d M N , ) (d H N , ) ! 0 as N ! 1 ,
where H N is a unit mesh for M N , and D(⌦) is the subset of C 1 ( ⌦) functions with compact support.

Similarly, given a smooth function u defined on ⌦, and its P 1 interpolation ⇧ h u on unit mesh H for M h , the interpolation error u ⇧ h u is also a highly oscillatory function. Now, according to [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], it is possible to homogenize the discrete interpolation error using the continuous interpolation error of Relation (4) expressed in terms of the Hessian H u of u in such a way that, for h ! 0,

|u ⇧ h u| = h 2 e homog interp (u) + h 2 e oscill interp (u) + h 2 o(h) L 2 (⌦) with e homog interp (u) = |u ⇡ M 1 u|
where e homog interp (u) is the homogenized interpolation error, h 2 o(h) L 2 (⌦) is an error term of order strictly higher than two according to L 2 analysis, and h 2 e oscill interp (u) is the oscillatory component of the interpolation error which is also a error term of order strictly higher that two since e oscill interp ! 0 in D 0 (⌦).

A priori finite-element analysis

Standard a priori estimates have been early derived in H 1 (⌦) ("projection property"), and in L 2 (⌦) (Aubin-Nitsche analysis), but only by means of inequalities. Moreover, the leading term of the error is generally not exhibited (only bounds of it are proposed). In this section, we go a little further in the Aubin-Nitsche a priori analysis to be able to consider the continuous mesh framework and to be able to exhibit the optimal adapted mesh such as in [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF][START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. To this end, we consider an important simplification in the analysis by neglecting the boundary error terms. It has already been done and discussed in previous works [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. In short, this is possible because in our type of approximation, close to FEM, the accuracy of the boundary discretization can degrade of one order of accuracy without changing the asymptotic convergence in H 1 (⌦) or L 2 (⌦). This simplification avoids to generate adapted meshes which would be strongly inhomogeneous close to boundaries which may degrades numerical scheme stability and accuracy. In this section, we focus on the Poisson problem which is set on domain ⌦ ⇢ R d :

u = f on ⌦ ; u = 0 on @⌦ . (5)
Its variational form writes:

a(u, v) = Z ⌦ ru.rv dx = ( f, v) , 8 v 2 V , (6)
where V holds for the Sobolev space V = H 1 0 (⌦) = u 2 L 2 (⌦), ru 2 (L 2 (⌦)) d , u | @⌦ = 0 . In order to derive an a priori estimate, we assume that solution u has some extra smoothness:

u 2 V = V \ C 3 ( ⌦).
where C 3 ( ⌦) is the set of functions of class C 3 on ⌦ [ @⌦. Let H be a mesh of ⌦ made of simplices (triangles in 2D and tetrahedrons in 3D): H = S k K k . Let V h be the subspace of V of continuous functions that are P 1 on each element of the mesh:

V h = n ' h 2 V ' h|K is a ne 8K 2 H o .
The discrete variational problem is then defined by:

a(u h , v h ) = ( f, v h ) , 8 v h 2 V h . (7)
Let us introduce the linear interpolation operator ⇧ h from vertices values:

⇧ h : V ! V h ; u 7 ! ⇧ h u such that ⇧ h u | K is a ne 8K 2 H and ⇧ h u(x) = u(x)
, for all vertices x of mesh H .

In a goal-oriented analysis, we are interested in estimating (g, u h u) which can be split into two components:

(g, u h u) = (g, u h ⇧ h u) + (g, ⇧ h u u) , (8)
where we recognize in the second di↵erence ⇧ h u u the interpolation error, and the first di↵erence u h ⇧ h u is referred as the implicit error. Introducing the continuous and the discrete adjoint states u ⇤ and u ⇤ h verifying:

a( , u ⇤ ) = (g, ) , 8 2 V and a( h , u ⇤ h ) = (g, h ) , 8 h 2 V h , we get (g, u h u) = a(u h ⇧ h u, u ⇤ h ) + (g, ⇧ h u u) . (9)
The second term of the right hand-side can be estimated without any di culty using Relation (4) :

|(g, ⇧ h u u)| Z ⌦ |g| |u ⇧ h u| d⌦ Z ⌦ c d |g| trace ⇣ M 1 2 |H u | M 1 2
⌘ d⌦ , while it would have been a lot more complicated using the equality:

(g, ⇧ h u u) = a(⇧ h u u, u ⇤ ) .
Analyzing the first term of the right hand-side of Relation (9) comes to study the following term:

a(u h ⇧ h u, ⇧ h ') ,
where ' is any su ciently smooth function. To this end, we first express the implicit error term as a function of the interpolation error. It is useful to remark that the discrete statement is equivalently written:

a(u h , ⇧ h ') = ( f, ⇧ h ') , 8 ' 2 V. ( 10
)
Using Relation [START_REF] Belme | Unsteady aerodynamic and adjoint method[END_REF] and then Relation (6), we get:

a(u h ⇧ h u, ⇧ h ') = a(u h , ⇧ h ') a(⇧ h u, ⇧ h ') = ( f, ⇧ h ') a(⇧ h u, ⇧ h ') = a(u, ⇧ h ') a(⇧ h u, ⇧ h ') , (11)
which gives:

a(u h ⇧ h u, ⇧ h ') = a(u ⇧ h u, ⇧ h ') , 8 ' 2 V. (12)
Note that u ⇧ h u is not solution of the discrete adjoint system because u is not in V h . Therefore, we propose the following main result:

Lemma 3.1. For any couple of smooth functions (u, '), where u is not necessarily a solution of Problem (5), we have the following bounds:

| Z ⌦ @ @x i (u ⇧ h u) @ @x j ⇧ h ' d⌦| K d Z ⌦ | ⇢ H (')| |u ⇧ h u| d⌦ + BT | Z ⌦ @ @x i u(u ⇧ h u) @ @x j ⇧ h ' d⌦| K d Z ⌦ | ⇢ H (')| |u| |u ⇧ h u| d⌦ + BT
where K d = 3 in two dimensions, K d = 6 in three dimensions, and A B holds for a majoration asymptotically valid, i.e. A  B + o(A) when mesh size tends to zero. Expression | ⇢ H (') | holds for spectral radius of H(') which is the Hessian of ', i.e., the largest (in absolute value) eigenvalue of H('). The boundary terms BT are not used in the sequel.

To prove this result, we analyze the right-hand side term of Equality [START_REF] Bottasso | Anisotropic mesh adaption by metric-driven optimization[END_REF]:

a(u ⇧ h u, ⇧ h ') = Z ⌦ r(u ⇧ h u) • r⇧ h ' d⌦ = X K2H Z K r(u ⇧ h u) • r⇧ h ' dK ,
where the sum ⌃ is taken over any element K of mesh H. As the aim is to extend this analysis to the Navier-Stokes equations, we analyze the following generic integral term (in which index i is not necessary equal to index j contrary to the Poisson problem):

I = X K2H Z K @ @x i (u ⇧ h u) @ @x j ⇧ h ' ! dK = X K2H Z @K (u ⇧ h u) @ @x j ⇧ h ' ! n K x i d Z K (u ⇧ h u) @ 2 @x i @x j ⇧ h ' ! dK ! ,
after integration by parts, and, as ⇧ h ' is linear over K, it reduces to:

I = X K2H Z K @ @x i (u ⇧ h u) @ @x j ⇧ h ' ! dK = X K2H Z @K (r⇧ h ') | K • e j n K x i (u ⇧ h u) d .
where (e j ) j=1..d stands for the canonical basis of R d , and n K x i is the i th component of unit outward normal

n K = (n K x 1 , n K x 2 , n K
x 3 ) T to the element boundary. As stated above, boundary integrals on @⌦ do not contribute to the volume estimate and are thus discarded in our analysis. The above integral can also be written as an integral over the edges in 2D or the faces in 3D of the mesh. Hence, to develop further this analysis, we consider the integral on the edge or the face sharing the two neighboring triangles or tetrahedrons K + and K . Thanks to the continuity of u ⇧ h u, it writes:

I = X @K + \@K Z @K + \@K ⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ • e j n K + x i (u ⇧ h u) d , where n K + x i is the i th component of n K + = (n K + x 1 , n K + x 2 , n K +
x 3 ) T the unit normal to the considered edge/face. Introducing the following notation for the jump of the '-derivative over the edge/face:

⇥ r⇧ h ' ⇤ + = ⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ ,
we observe that the gradient component of ⇧ h ' in direction t K + tangent to the common edge/face is continuous, therefore we have:

⇥ r⇧ h ' ⇤ + • t K + =
0 , from which we deduce the following relations with the unit edge/face normal:

⇥ r⇧ h ' ⇤ + = ⇥ r⇧ h ' ⇤ + n K + and ⇥ r⇧ h ' ⇤ + = ⇥ r⇧ h ' ⇤ + • n K + .
From the above relations, we get:

⇥ r⇧ h ' ⇤ + • e j = ⇥ r⇧ h ' ⇤ + n K + x j = ⇥ r⇧ h ' ⇤ + • n K + n K + x j ,
and, finally, we obtain:

⇥ r⇧ h ' ⇤ + • e j n K + x i = ⇥ r⇧ h ' ⇤ + • n K + n K + x j n K +
x i . Thus, the above integral becomes:

I = X @K + \@K Z @K + \@K ⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ • n K + n K + x j n K + x i (u ⇧ h u) d ,
|I| = Z ⌦ @ @x i (u ⇧ h u) @ @x j ⇧ h ' ! d⌦ X @K + \@K ⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ • n K + Z @K + \@K |u ⇧ h u| d ,
To pursue our analysis, this integral can be transformed into an integral over ⌦ by considering a diamond cell partitioning:

H = S k K k = S i D i
, where each diamond cell D i is associated to an edge e in 2D or a face f in 3D. In two dimensions, a diamond cell D e associated to an edge e is the union of the two sub-triangles that are build by joining the centers of gravity of the two triangles sharing the edge e to the edge e, see Figure 1 (left). In three dimensions, a diamond cell D f attached to a face f is the union of the two sub-tetrahedrons that are build by joining the centers of gravity of the two tetrahedrons sharing the face f to the face f , see Figure 1

|D e | = |K + | + |K | 3 = |e| |h + | + |h | 6 , (13)
and for the volume in 3D:

|D f | = |K + | + |K | 4 = | f | |h + | + |h | 12 . (14)
Moreover, taking the mean of |u ⇧ h u| over the edge e = @K + \ @K or the face f = @K + \ @K is a consistent quadrature for the integration of |u ⇧ h u| over the diamond cell D + , thus:

1

|@K + \ @K | Z @K + \@K |u ⇧ h u| d ⇡ 1 |D + | Z D + |u ⇧ h u| dD + , (15)
And finally, we observe that the term:

⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ • n K + = @ @n K + ⇧ h '| K + @ @n K ⇧ h '| K
appears as a second derivative of ' in the normal direction to the edge/face weighted by the inverse of the average height of the two neighboring simplex. If we introduce | ⇢ H (') | the spectral radius of H(') which is the Hessian of ', the spectral radius being the largest (in absolute value) eigenvalue of H('), we get the following bound:

⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ • n K + |h + | + |h | 2 | ⇢ H (')| . (16)
In two dimensions, using Relations ( 13), [START_REF] Brèthes | A mesh-adaptive metric-based full-multigrid for the Poisson problem[END_REF], and ( 16), we obtain the following upper bounds:

|I| X @K + \@K ⇥ (r⇧ h ') | K + (r⇧ h ') | K ⇤ • n K + Z @K + \@K |u ⇧ h u| d X e2H |h + | + |h | 2 | ⇢ H (')| ! 0 B B B B B @ |e| |e| |h + |+|h | 6 Z D e |u ⇧ h u| dD e 1 C C C C C A ,
and after canceling terms, we finally obtain:

|I| X e2H 3 | ⇢ H (')| Z D e |u ⇧ h u| dD e = 3 Z ⌦ | ⇢ H (')| |u ⇧ h u| d⌦ , in 2D .
Similarly, in three dimensions, using Relations ( 14), [START_REF] Brèthes | A mesh-adaptive metric-based full-multigrid for the Poisson problem[END_REF], and ( 16) and canceling terms, we obtain the following upper bounds:

|I| X f 2H 6 | ⇢ H (')| Z D f |u ⇧ h u| dD f = 6 Z ⌦ | ⇢ H (')| |u ⇧ h u| d⌦ , in 3D .
The same demonstration holds if we replace the term (u ⇧ h u) by the term u (u ⇧ h u). This concludes the proof of the Lemma 3.1. The estimate of Lemma 3.1 has been successfully tested for mesh adaptation in [START_REF] Brèthe | Anisotropic norm-oriented mesh adaptation for a Poisson problem[END_REF] and also compared with another estimate in [START_REF] Brèthe | A tensorial-based mesh adaptation for a Poisson problem[END_REF].

Navier-Stokes Model

Continuous state system

The steady compressible Navier-Stokes system for a perfect gas is set in computational domain ⌦ ⇢ R 3 and can be written under a compact form as:

r • F E (W) + r • F V (W) = 0 on ⌦ ( 17
)
where W = t (⇢, ⇢u, ⇢E) is the conservative flow variables vector and vector F E represents the Euler fluxes:

F E (W) = t (⇢u, ⇢u 1 u + pe 1 , ⇢u 2 u + pe 2 , ⇢u 3 u + pe 3 , ⇢uH) .
We have noted ⇢ the density, u = (u 1 , u 2 , u 3 ) the velocity vector, H = E + p/⇢ is the total enthalpy, E = T + kuk 2 2 the total energy and p = ( 1)⇢T the pressure with = 1.4 the ratio of specific heat capacities, T the temperature, and (e 1 , e 2 , e 3 ) the canonical basis.

We describe in short the viscous fluxes and viscous stress tensor as follows:

F V (W) = [0, , (q u. )] T ; = µ(ru + ru T ) 2 3 µr.uI,
with µ representing the constant viscosity. The heat flux q is given by Fourier's law q = rT , where is the heat conduction (assumed here to be constant). The five unknowns are gathered in functional product space (we use the same notation):

V = h H 1 (⌦) \ C 3 ( ⌦) i ⇥ h H 1 0 (⌦) \ C 3 ( ⌦) i 3 ⇥ h H 1 (⌦) \ C 3 ( ⌦) i , (18)
assuming adiabatic conditions on walls. We formulate the Navier-Stokes model in a compact variational formulation:

Find W 2 V such that 8 2 V, ( (W) , ) = 0 with = E + + V , (19) where
= ⇣ ⇢ , ⇢u 1 , ⇢u 2 , ⇢u 3 , ⇢E ⌘ T .
The Euler term E relies on the usual Euler fluxes F E :

⇣ E (W) , ⌘ = Z ⌦ • r • F E (W) d⌦. ( 20
)
Term holds for boundary fluxes contribution which we denote shortly:

⇣ (W) , ⌘ = Z • F (W) • n d .
Viscous fluxes F V provide seven terms:

( V (W), ) = Z ⌦ • r • F V (W) d⌦ = 7 X k=1 T V k . (21)
The first three terms come from moment equations and depend only on ⇢u = ( ⇢u 1 , ⇢u 2 , ⇢u 3 ) T :

T V 1 = Z ⌦ ⇢u • r • (µru) d⌦ T V 2 = Z ⌦ ⇢u • r • ⇣ µ(ru) T ⌘ d⌦ T V 3 = 2 3 Z ⌦ ⇢u • r • (µ (r • u) I) d⌦ .
The last four terms are derived from the energy equation:

T V 4 = Z ⌦ ⇢E r • ( rT ) d⌦ T V 5 = Z ⌦ ⇢E r • ⇣ µ u • (ru) T ⌘ d⌦ T V 6 = Z ⌦ ⇢E r • (µ u • ru) d⌦ T V 7 = 2 3 Z ⌦ ⇢E r • (µ u • ((r • u) I)) d⌦ .

Variational discrete formulation

For the spatial semi-discrete model, we consider the mixed Finite Element -Finite Volume formulation [START_REF] Cournède | Positivity statements for a Mixed-Element-Volume scheme on fixed and moving grids[END_REF][START_REF] Menier | CFD validation and adaptivity for viscous flow simulations[END_REF]. As in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF], we reformulate it under the form of a finite element variational formulation. We assume that ⌦ is covered by a unit mesh H with respect to a given Riemannian metric space M = (M(x)) x2⌦ composed of simplicial elements, denoted K. Let us introduce the following approximation space:

V h = n h 2 V h|K is a ne 8K 2 H o
, where V is defined by Relation ( 18).

The weak discrete formulation writes:

Find W h 2 V h such that 8 h 2 V h , h (W h ) , h = 0, with h (W h ) , h = Z ⌦ h • r • F E h (W h ) d⌦ Z h • F h (W h ).n d + Z ⌦ h • r • F V h (W h ) d⌦ + Z ⌦ h • D h (W h ) d⌦ . ( 23
)
The fourth term in Formulation ( 23) is an added numerical di↵usion dedicated to numerical stability. In short, the D h term involves the di↵erence between the Galerkin central-di↵erences approximation and a second-order Godunov approximation defined as in [START_REF] Cournède | Positivity statements for a Mixed-Element-Volume scheme on fixed and moving grids[END_REF]. In the present study, we only need to know that for smooth fields, the D h term is a third-order term with respect to the mesh size parameter h and will not be considered in our estimate.

Taking in System (23) the P 1 interpolation of the fluxes F as discretization principle, produces a finite-element scheme which is identical to the central-di↵erenced finite-volume scheme built on the so-called median dual cells [START_REF] Menier | CFD validation and adaptivity for viscous flow simulations[END_REF], thus:

8 W 2 V h [ V , F E h (W) = ⇧ h F E (W) and F h (W) = ⇧ h F (W) . (24)
However, concerning the viscous term, we cannot apply a similar treatment because we have fluxes of order two (second derivatives) and we cannot apply directly the interpolation operator. Let f be the transformation function from primitive variables U = (⇢, u, T ) into conservatives ones W = (⇢, ⇢u, ⇢E), we set:

F V h (W) = F V h ( f (U)) = F V ( f (⇧ h U)) = F V ( f (⇧ h f 1 (W))
). In other words, our discretization consists of P 1 interpolating the primitive variables in the above viscous terms {T V i } i=1..7 for our discretization principle. This completes the definition of the discrete system under study.

Mesh adaptation : discrete problem statement

Let g be a function of V. We assume that the purpose of the numerical problem is to evaluate the output functional:

j = (g, W)
where W is the solution of Problem [START_REF] Castro | Systematic continuous adjoint approach to viscous aerodynamic design on unstructured grids[END_REF].

The problem addressed in this paper is to find the discrete mesh which minimizes the following functional error given a fixed number of vertices N: j = (g, W W h ) where W is the solution of Problem ( 19) and W h is the solution of Problem [START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-di↵usion-reaction and the Stokes problems[END_REF].

The next sections are devoted to this error analysis and the optimal formulation of the mesh adaptation problem.

Error analysis for Navier-Stokes problem

The Navier-Stokes equations are a non-linear system, thus some extra justifications are required to be able to perform an error analysis similar to the linear case (done in Section 3). The following justify why it is possible to apply a similar analysis in the non-linear case.

Let be a smooth test function of V defined in Relation [START_REF] Carabias | Analyse et adaptation de maillage pour des schémas non-oscillatoires d'ordre élevé[END_REF]. Let W be the solution of Problem [START_REF] Castro | Systematic continuous adjoint approach to viscous aerodynamic design on unstructured grids[END_REF] and W h the solution of Problem [START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-di↵usion-reaction and the Stokes problems[END_REF], the continuous and discrete state equation write:

( (W) , ) = 0 and ( h (W h ) , ⇧ h ) = 0 ,
where ⇧ h lies in V h defined by Relation [START_REF] Dompierre | Anisotropic mesh adaptation: towards a solver and user independent CFD[END_REF]. We also introduce the continuous and the discrete adjoint states: W ⇤ and W ⇤ h . The continuous adjoint system related to the objective functional writes:

W ⇤ 2 V , 8 2 V : @ @W (W) , W ⇤ ! = (g, ) . (25)
From functional analysis theory, a well-posed continuous adjoint system can be derived for any functional output as far as the linearized system is well posed. This however does not mean that any output functional leads to properly defined adjoint boundary conditions. Several works in the literature [START_REF] Anderson | Aerodynamic design optimization on unstructured grids with a continuous adjoint formulation[END_REF][START_REF] Arian | Admitting the inadmissible: Adjoint formulation for incomplete cost functionals in aerodynamic optimization[END_REF][START_REF] Bueno-Orovio | Continuous adjoint approach for the Spalart-Allmaras model in aerodynamic optimization[END_REF][START_REF] Castro | Systematic continuous adjoint approach to viscous aerodynamic design on unstructured grids[END_REF] illustrate this problem and propose solutions, usually by adding auxiliary boundary terms to the Lagrangian functional. In [START_REF] Castro | Systematic continuous adjoint approach to viscous aerodynamic design on unstructured grids[END_REF], it is concluded that for the compressible Navier-Stokes system, only functionals which involve the entire stress tensor at obstacle boundary are admissible. We assume here that ( 25) is well-posed and gives a su ciently smooth continuous adjoint state. The discrete adjoint system writes:

W ⇤ h 2 V h , 8 h 2 V h : @ h @W (W h ) h , W ⇤ h ! = (g, h ) .

Linearized error system

In our error estimation problem presented in Section 4.3, the approximation error can be decomposed into an implicit error term and an interpolation error term:

j = (g, W W h ) = (g, W ⇧ h W) + (g, ⇧ h W W h ) .
The interpolation error can be easily estimated, see Relation (4), while the implicit error is solution of a discrete system that we derive in the following. Now, we assume that W h can be made close enough to ⇧ h W when h ! 0 in such a way that we can identify the main term of the left-hand side as a Jacobian times the di↵erence:

@ h @W (W h )(W h ⇧ h W) , ⇧ h ! ⇡ ⇣ h (W h ) h (⇧ h W) , ⇧ h ⌘ as h ! 0.
Then, combining continuous and discrete systems, we can write similarly to Relation [START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] an equality linking implicit and interpolation errors which is valid for all :

( h (W h ) h (⇧ h W) , ⇧ h ) = ( h (W h ) , ⇧ h ) ( h (⇧ h W) , ⇧ h ) = ( (W) , ⇧ h ) ( h (⇧ h W) , ⇧ h ) = ( (W) h (⇧ h W) , ⇧ h )
We are then interested in the following error on functional:

|(g, W h ⇧ h W)| ⇡ ⇣ h (W h ) h (⇧ h W) , W ⇤ h ⌘ = ⇣ (W) h (⇧ h W) , W ⇤ h ⌘ . ( 26
)
The right-hand side of the above relation is composed of (see Relation [START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-di↵usion-reaction and the Stokes problems[END_REF]) the Euler term, the Euler boundary term, the viscous term, and the stabilization term. In the following, we neglect the boundary term (as in the previous analysis) and the stabilization term thanks to smoothness of functions W and W ⇤ . The boundary term analysis is given in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF].

The method proposed here involves some heuristics. In particular, we assume that the interpolate of the adjoint is close to the discrete adjoint:

⇧ h W ⇤ ⇡ W ⇤ h . Therefore: |(g, W h ⇧ h W)| ⇡ ( (W) h (⇧ h W) , ⇧ h W ⇤ ) .

A priori estimate for Navier-Stokes problem

We can now give the main result of this paper which is the following a priori estimate:

Proposition 5.1. Let us assume that W 2 V and 2 V where V is defined by Relation ( 18) and = ⇣ ⇢ , ⇢u 1 , ⇢u 2 , ⇢u 3 , ⇢E ⌘ T . Then, we have the following error bound for h su ciently small:

|( (W) h (⇧ h W), ⇧ h )| E
with: The proof of this Proposition is given in Appendix A. Le last term of this estimation requires additional manipulations (done in the Appendix) to make its implementation easier. This term is expressed in the following by means of the vector !.

E = Z ⌦ | r | • F E (W) ⇧ h F E (W) d⌦ + d + 1 3 ! K d d X i=1 Z ⌦ µ | ⇢ H ( ⇢u i )| |u i ⇧ h u i | d⌦ + 1 3 K d d X i=1 d X j=1 j,i Z ⌦ µ | ⇢ H ( ⇢u i )| |u j ⇧ h u j | d⌦ + K d Z ⌦ | ⇢ H ( ⇢E )| |T ⇧ h T | d⌦ + (d + 1 3 ) K d d X i=1 Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u i ⇧ h u i | d⌦ + 1 3 K d d X i=1 d X j=1 j,i Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u j ⇧ h u j | d⌦ + 5 3 d X i=1 d X j=1 j,i Z ⌦ µ @(⇧ h ⇢E ) @x j (u j ⇧ h u j ) @u i @x i (u i ⇧ h u i ) @u j @x i
This main result can be written in a more convenient way to facilitate its implementation by gathering interpolation error terms on the primitive variables. Here, we give the 2D and the 3D re-writing of it, using notations introduced in Appendix A.

Corollary 5.1. Let us assume that W 2 V and 2 V where V is defined by Relation (18) and =

⇣ ⇢ , ⇢u 1 , ⇢u 2 , ⇢E ⌘ T .
Then, in two dimensions (K d = 3), we have the following error bound for h su ciently small:

|( (W) h (⇧ h W), ⇧ h )| Z ⌦ G F E F E (W) ⇧ h F E (W) d⌦ + Z ⌦ G u 1 |u 1 ⇧ h u 1 | d⌦ + Z ⌦ G u 2 |u 2 ⇧ h u 2 | d⌦ + Z ⌦ G T |T ⇧ h T | d⌦ ,
with the coe cients:

G F E = | r | G u 1 = µ 7 | ⇢ H ( ⇢u 1 )| + | ⇢ H ( ⇢u 2 )| + ⇣ 7|u 1 | + |u 2 | ⌘ | ⇢ H ( ⇢E )| + 5 3 |! u 2 ,z | ! G u 2 = µ | ⇢ H ( ⇢u 1 )| + 7 | ⇢ H ( ⇢u 2 )| + ⇣ |u 1 | + 7|u 2 | ⌘ | ⇢ H ( ⇢E )| + 5 3 |! u 1 ,z | ! G T = 3 | ⇢ H ( ⇢E )|
and the ! vector defined by:

ru i ⇥ r ⇢E = ! u i = ⇣ ! u i ,x , ! u i ,y , ! u i ,z ⌘ T ,
where only the last component is not zero.

Corollary 5.2. Let us assume that W 2 V and 2 V where V is defined by Relation ( 18) and = ⇣ ⇢ , ⇢u 1 , ⇢u 2 , ⇢u 3 , ⇢E ⌘ T . Then, in three dimensions (K d = 6), we have the following error bound for h su ciently small:

|( (W) h (⇧ h W), ⇧ h )| Z ⌦ G F E F E (W) ⇧ h F E (W) d⌦ + Z ⌦ G u 1 |u 1 ⇧ h u 1 | d⌦ + Z ⌦ G u 2 |u 2 ⇧ h u 2 | d⌦ + Z ⌦ G u 3 |u 3 ⇧ h u 3 | d⌦ + Z ⌦ G T |T ⇧ h T | d⌦ ,
with the coe cients:

G F E = | r | G u 1 = µ 20 | ⇢ H ( ⇢u 1 )| + 2 | ⇢ H ( ⇢u 2 )| + 2 | ⇢ H ( ⇢u 3 )| + ⇣ 20|u 1 | + 2|u 2 | + 2|u 3 | ⌘ | ⇢ H ( ⇢E )| + 5 3 |! u 3 ,y ! u 2 ,z | ! G u 2 = µ 2 | ⇢ H ( ⇢u 1 )| + 20 | ⇢ H ( ⇢u 2 )| + 2 | ⇢ H ( ⇢u 3 )| + ⇣ 2|u 1 | + 20|u 2 | + 2|u 3 | ⌘ | ⇢ H ( ⇢E )| + 5 3 |! u 1 ,z ! u 3 ,x | ! G u 3 = µ 2 | ⇢ H ( ⇢u 1 )| + 2 | ⇢ H ( ⇢u 2 )| + 20 | ⇢ H ( ⇢u 3 )| + ⇣ 2|u 1 | + 2|u 2 | + 20|u 3 | ⌘ | ⇢ H ( ⇢E )| + 5 3 |! u 2 ,x ! u 1 ,y | ! G T = 6 | ⇢ H ( ⇢E )|
and the ! vector defined by:

ru i ⇥ r ⇢E = ! u i = ⇣ ! u i ,x , ! u i ,y , ! u i ,z ⌘ T .
For the next section where we seek for the optimal mesh that minimize the above error model, it is useful to note that the error model can be written under the compact form:

|( (W) h (⇧ h W), ⇧ h )| Z ⌦ X k G k (µ, , U, r , | ⇢ H ( )|) S k (W) ⇧ h S k (W) d⌦ , (27)
in other words, the error model is a sum of interpolation errors weighted by algebraic functions.

6. Goal-oriented optimal continuous mesh 6.1. Mesh adaptation : continuous problem statement According to the continuous mesh theory [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] presented in Section 2, we transform the a priori error estimate of Proposition 5.1 into a continuous error estimate. To this end, the following modifications are done:

• mesh H is replaced by continuous mesh M = (M(x)) x2⌦
• number of vertices N is replaced by complexity C(M) = N

• the P 1 projection operator ⇧ h defined on mesh H is replaced by the P 1 continuous mesh operator ⇡ M defined on continuous mesh M = (M(x)) x2⌦

• discrete interpolation error u ⇧ h u defined on mesh H is replaced by continuous interpolation error u ⇡ M u defined on continuous mesh M = (M(x)) x2⌦ .

Then, using Relations ( 26) and ( 27), we can reformulate the discrete mesh adaptation problem introduced in Section 4.3 in the continuous mesh framework: find the continuous mesh M opt which minimizes the following functional error given a fixed complexity N:

j = (g, W W h ) ⇡ (g, W ⇧ h W) + ( (W) h (⇧ h W) , ⇧ h W ⇤ ) Z ⌦ X k |g k | W k ⇡ M W k d⌦ + Z ⌦ X k G k (µ, , U, rW ⇤ , | ⇢ H (W ⇤ )|) S k (W) ⇡ M S k (W) d⌦ , ( 28
)
where W is the solution of Problem [START_REF] Castro | Systematic continuous adjoint approach to viscous aerodynamic design on unstructured grids[END_REF], W h is the solution of Problem [START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-di↵usion-reaction and the Stokes problems[END_REF] and W ⇤ is the solution of Problem [START_REF] Frey | About surface remeshing[END_REF].

As the above relation is a weighted sum of interpolation errors, by introducing the following positive symmetric matrix

H go (x) = X k |g k | |H W k | + X k G k (µ, , U, rW ⇤ , | ⇢ H (W ⇤ )|) |H S k (W) | (29)
where |H W k | and |H S k (W) | are the absolute value of the Hessian of fields W k and S k (W), and using the definition of the continuous interpolation error given by Relation (4), we can state the following error estimate on continuous mesh M = (M(x)) x2⌦ :

j ⇡ E go (M) = c d Z ⌦ trace ⇣ M 1 2 (x) H go (x) M 1 2 (x) ⌘ d⌦ .
It is then possible to set the well-posed global optimization problem of finding the optimal continuous mesh M go minimizing continuous interpolation error E go (M):

Find M go = min M E go (M) under the constraint C(M) = N . (30)

Optimal goal-oriented continuous mesh

We seek for the optimal continuous mesh M go solution of Problem [START_REF] Hecht | Mesh adaptation by metric control for multi-scale phenomena and turbulence[END_REF]. Similarly to [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF], solving the optimality conditions provides the optimal goal-oriented continuous mesh M go = (M go (x)) x2⌦ defined point-wise by:

M go (x) = N 2 d Z ⌦ (det H go (x)) 1 d+2 dx ! 2 d ⇣ det H go (x) ⌘ 1 d+2 H go (x). ( 31
)
where d is the dimension. The corresponding optimal error writes:

E go (M go ) = = 3 N 2 d Z ⌦ (det H go (x)) 1 d+2 dx ! 2+d d , (32)
where the exponent of N illustrates the second-order accuracy of the method.

From theory to practice

The continuous mesh adaptation problem takes the form of the following continuous optimality system:

W 2 V , 8 2 V , ( (M, W), ) = 0 "Navier-Stokes system" W ⇤ 2 V , 8 2 V , @ @W (M, W) , W ⇤ ! = (g, ) "Adjoint system" M = M go "Adapted continuous mesh".
In practice, it is necessary to approximate the above three-field coupled system by the discrete optimality system:

W 2 V h , 8 h 2 V h , h (H, W h ), h = 0 "Discrete Navier-Stokes system" W ⇤ h 2 V h , 8 h 2 V h , @ @W (H, W h ) h , W ⇤ h ! = (g, h ) "Discrete Adjoint system" H = H go "Discrete adapted mesh".
The discrete optimality system is solved using the following fixed-point mesh adaptation algorithm:

Algorithm 1 Viscous Goal-Oriented Mesh Adaptation Loop for Steady Flows Initial mesh and solution (H 1 go , S 1 0 ) and set targeted functional j and complexity N # Adaptive loop to converge the steady-state solution and its associated optimal adapted discrete mesh For i = 1, n adap The following sections detail each step of the above algorithm. The state and adjoint solutions are given in Section 7.1. The computation of the discrete metric field is given in Section 7.2. And, the methodology to generate the adapted discrete mesh is discussed in Section 7.3. The last step is not presented here, we refer to [START_REF] Alauzet | A parallel matrix-free conservative solution interpolation on unstructured tetrahedral meshes[END_REF] for more details.

The flow solver: Wolf

The flow is modeled by the conservative laminar Navier-Stokes equations given by Relation [START_REF] Bueno-Orovio | Continuous adjoint approach for the Spalart-Allmaras model in aerodynamic optimization[END_REF].

The spatial discretization of the fluid equations is based on a vertex-centered finite element/finite volume formulation on unstructured meshes composed of triangles in 2D and tetrahedrons in 3D. It combines a HLLC upwind schemes for computing the convective fluxes and the Galerkin centered method for evaluating the viscous terms. Second order space accuracy is achieved through a piecewise linear interpolation based on the Monotonic Upwind Scheme for Conservation Law (MUSCL) procedure which uses a particular edge-based formulation with upwind elements. A specific low dissipation scheme is considered using combination of centered and upwind gradients. A dedicated slope limiter is employed to damp or eliminate spurious oscillations that may occur in the vicinity of discontinuities.

The temporal discretization considers implicit BDF1 time advancing scheme. To solve the non-linear system, we follow the approach based on Symmetric Gauss-Seidel (SGS) implicit solver. The Newton's method based on the SGS relaxation is very attractive because they use an edge-based data structure which can be e ciently parallelized with p-threads. From our experience, we have made the following -crucial -choices to solve the compressible Navier-Stokes equations.

• the SGS relaxation iterates until the residual of the linear system is reduced by two orders of magnitude

• the Breadth-first search renumbering proves to be the more e↵ective for the convergence of the implicit method and the overall e ciency

• we found very advantageous to fully di↵erentiate the HLLC approximate Riemann solver and the FEM viscous terms

• to achieve high e ciency, automation and robustness in the resolution of the non-linear system of algebraic equations to steady-state, it is mandatory to have a clever strategy to specify the time step.

The flow solver Wolf is thoroughly detailed in [START_REF] Menier | CFD validation and adaptivity for viscous flow simulations[END_REF] with all the associated bibliography.

As regards the adjoint state computation, the matrix of the linear system is simply the implicit matrix transposed and the right hand-side of the system is the chosen functional (for instance, drag, lift, ...) exactly di↵erentiated:

@ h @W h T W ⇤ h = @ j @W h .
In particular, for viscous flows, µ and the stress tensor ⌧ are exactly di↵erentiated. To solve the adjoint system, we use a restarted GMRES preconditioned with LUSGS relaxation. Note that, it is important to solve the adjoint linear system to machine precision.

Computation of the optimal continuous mesh

The expression of the optimal goal-oriented metric field is given by Relations ( 29) and [START_REF] Jones | Validation of 3D adjoint based error estimation and mesh adaptation for sonic boom reduction[END_REF] where algebraic functions G k and S k are given in Corollaries 5.1 and 5.2 where is replaced by the adjoint state W ⇤ . This expression involves the state and the adjoint state, and their derivatives (first and second). In practice, these terms are approximated by the discrete states and a derivative recovery is applied to get gradients and Hessians. In this paper, the recovery method is based on the L 2 -projection formula. Its formulation along with some comparisons to other methods is available in [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF]. We then obtain a discrete metric field defined at vertices of the mesh.

The discrete adjoint state W ⇤ h is taken to represent the adjoint state W ⇤ . The gradient of the adjoint state rW ⇤ is replaced by r R W ⇤ h where r R stands for the operator that recovers numerically the first derivatives of an initial piecewise linear solution field. The Hessian is obtained by applying the recovery operator two times:

H R = r R r R . Then, | ⇢ H (W ⇤ )| is obtained from | ⇢ H (W ⇤
h )| which is evaluated by computing the maximal, in absolute value, eigenvalue of H R (W ⇤ h ). Similarly, the discrete state W h is taken to represent the state W and to compute the Euler fluxes, velocity components and temperature (i.e., each term S k (W) is replaced by S k (W h )). The Hessians of the S k (W h ) are obtained using the recovery operator H R and the absolute values of the Hessians are computed by taking the absolute values of the eigenvalues.

The local adaptive remesher: Feflo.a

The adaptive remesher used in this paper is based a combination of generalized standard operators (insertion, collapse, swap of edges and faces). The generalized operators are based on recasting the standard operators in a cavity framework [START_REF] Loseille | Cavity-based operators for mesh adaptation[END_REF]. Additional modifications on the cavity allow to either favor a modification, that would have been rejected with the standard operator, or to improve the final quality by combining automatically many standard operators at once. In addition, the CPU time is also improved and becomes independent of the current modification. The unit speed is around 20, 000 points inserted or removed per second on Intel i7 architecture at 2.7 GHz. For robustness purpose, both the surface and the volume mesh are adapted simultaneously, and each local modification is checked to verify that a valid mesh is obtained. For the volume, the validity consists in checking that each newly created element has a strictly positive volume. For the surface, the validity is checked by ensuring that the deviation of the geometric approximation with respect to a reference surface mesh remains within a given tolerance [START_REF] Frey | About surface remeshing[END_REF].

The boundary layer region involves highly directional features requiring the generation of highly stretched meshes (anisotropic ratio from thousands to hundreds thousands). To control the overall mesh quality and to avoid the creation of large dihedral angles, quasi-structured mesh generation in boundary layer has been favored (for instance [START_REF] Marcum | Adaptive unstructured grid generation for viscous flow applications[END_REF][START_REF] Pirzadeh | Viscous unstructured three dimensional grids by the advancing-layers method[END_REF]) because the surface of the geometry can be used as a support. Such quasi-structured meshes impacts favorably the stability and the accuracy of the numerical schemes. The generation of such quasi-structured meshes has been generalized in a metric-based mesh adaptation context [START_REF] Loseille | Metric-orthogonal anisotropic mesh generation[END_REF][START_REF] Loseille | Alignment and orthogonality in anisotropic metric-based mesh adaptation[END_REF][START_REF] Marcum | Aligned metric-based anisotropic solution adaptive mesh generation[END_REF] allowing quasi-structured stretched elements regions to be generated in the vicinity of any physical features (shocks, wakes, shear layers, ...). To this end, a procedure to generate metric-aligned and metric-orthogonal anisotropic adapted meshes has been proposed. The metric-aligned approach uses a point placement that aligns the elements with the eigenvectors of the underlying metric field. This simple strategy produces solution-adapted anisotropic meshes that are optimal in size and alignment with the solution gradients upon which the metric field is based. When the goal is to generate locally structured meshes in an anisotropic context and in the presence of complex geometries, a metric-orthogonal approach is used. In that case, the point placement proposes points that align with the eigenvectors frame of the underlying metric field. The adaptive algorithm relies a frontal approach to propose the points which are iteratively inserted using the cavity-based point insertion operator [START_REF] Loseille | Cavity-based operators for mesh adaptation[END_REF]. The local alignment and orthogonality are naturally inherited from the eigenvectors and eigenvalues of an input metric field.

This strategy has been considered in all the numerical examples presented hereafter.

Numerical Experiments

In order to validate the aforementioned error estimate for viscous flows we have performed a series of simulations of compressible flows dominated by viscous e↵ects around a NACA0012 airfoil. The first three test cases studied hereafter are rather well studied cases by the aerodynamic community, providing thus (for some examples) useful comparison support.

The initial mesh employed for these simulations is triangular, unstructured of about 16 000 vertices in a circular domain of radius 100 and the airfoil chord length is 1. We are interested in deriving the best goal-oriented adapted anisotropic mesh to observe some quantity of interest. For the first three problems, we have considered the total (induced plus friction) drag coe cient C tot D integrated over the airfoil profile:

C tot D = C p + C f = 1 1 2 ⇢ 1 ku 1 k 2 L re f Z (p p 1 )(n x cos ↵ + n y sin ↵) ⇣ cos ↵(⌧ xx n x + ⌧ xy n y ) + sin ↵(⌧ xy n x + ⌧ yy n y ) ⌘ .
where n = (n x , n y ) is the outward unit normal vector to the surface , 1 2 ⇢ 1 ku 1 k 2 is the free-stream dynamic pressure, ↵ the angle of attack, and L re f is the chord reference length chosen equal to 1.0089. Note that after each mesh adaptation, vertices located on the NACA0012 geometry are reprojected exactly on the geometry.

It is interesting to see how the di↵erent terms of the error estimate influence the generation of the adapted meshes for the computation of the total drag. To this end, for each test case, we compare the viscous and the inviscid [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF] (i.e., without the viscous terms) error estimates by confronting the corresponding meshes in order to emphasize the impact of the error estimate viscous terms in the mesh adaptation process. Moreover, as the proposed error estimate should be optimal, it must at least do better than uniform refinement and the inviscid error estimate. Thus, we compare the drag convergence using the viscous and the inviscid error estimates to assess the optimality new approach. Remark 8.1. We have now to recall that viscous terms are integrated with the P 1 -Galerkin discretization. In contrast to many structured finite-volume methods no special high-order interpolation is applied near the wall. Instead, the P 1 -Galerkin approximation stencil reduces to solely the wall vertex and its direct neighbors. According to the finite-element theory, convergence is second order in L 2 and first-order in H 1 , that is for derivatives of unknowns. Taking the trace at wall of the derivative may degrade this order, but the integral over the wall may upgrade it (through a Green formula). This short a priori analysis leads to a first-order convergence for the drag. Then, we have to rely on a super-convergence e↵ect if we want a convergence of the drag coe cient better than one. In our test cases, drag levels are in some cases compared with reference values obtained with hundreds-vertices finer meshes. Drag numerical convergence order is never measured with these reference values but from the numerical outputs.

Note that there is only few accuracy order statements for mesh adaptation methods, see [START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF] for second order method and [START_REF] Yano | The importance of mesh adaptation for higher-order discretizations of aerodynamics flows[END_REF] for higher order methods. A numerical convergence order ↵ can be measured according to a power-of-N law, where N is the number of nodes:

error = O(N ↵ d )
where d holds for the dimension of the computational domain . When the exact solution is available, the error is a deviation with respect to it. When the exact solution is not available, either the above relation is computed with a very accurate reference solution, or it must be injected in a formula using several solutions with di↵erent mesh complexity N. This measure can be applied to essentially any mesh-based series of calculations with an increasing N (structured, unstructured) and its use for mesh adaptation dates back to the beginning of h p adaptive methods. The above measure needs to be interpreted, this can be done with respect to some principles like:

-a mesh adaptation method should show a convergence order close to the uniform mesh convergence order, possibly with a much lower number of nodes -the convergence order measured for a mesh adaptation method should be better that uniform mesh convergence order when applied to fields with singularities.

In two dimensions, the convergence order ↵ of our quantity of interest j for each case is measured according to the relation:

↵ i = log(| j i j re f |) log(| j i 1 j re f |) log(h i ) log(h i 1 ) with h i = 1 p N i ,
where N is the number of vertices of the mesh and j re f is the quantity of interest for the reference solution. For each case, the reference solution is the solution obtained by the mesh adaptation process on the finest adapted mesh. As in a mesh adaptation context the mesh size is varying in the domain, it is meaningless to analyze the convergence using the size h of the elements.

It is a lot more appropriate to consider the mesh size N because it represents the number of degrees of freedom. The first test case is the laminar transonic flow around the NACA0012 airfoil. For this case, the free-stream Mach number is 0.8, the angle of attack is 10 and the Reynolds number is 73. This problem has been investigated by several authors for the GAMM Workshop (see [START_REF] Bristeau | Presentation of problems and discussion of results[END_REF]) on the Solution of Compressible Navier-Stokes flows. The Mach number contours of the computed solution are depicted in Figure 2, where a thick boundary layer along the airfoil surface is observed, and locally supersonic flow is attained only in a small pocket outside the edge of the viscous layer on the upper surface.

The considered functional of interest is the total drag coe cient integrated along the airfoil profile. We have considered our new a priori error estimate to compute the associated optimal anisotropic adapted meshes for varying mesh complexity. Multiple views of the resulting adapted mesh composed of about 38 000 vertices obtained after performing 6 fixed-point iterations are shown in Figure 4 (right column). We observe that most of the refinement is located around the airfoil (boundary layer) and outside the edge of the viscous layer on the upper and lower regions of the domain. We also notice strong refinement at the trailing edge of the airfoil.

We also compared adapted meshes obtained using the error estimate with and without the viscous terms. In Figure 4, we compare the goal-oriented adapted mesh obtained by adapting only for the convective flux with the goal-oriented adapted mesh where both fluxes (convective and di↵usive) are considered. Both meshes have a similar number of vertices. In the far field, we observe similar anisotropic refinement, whereas close to the body we notice that the viscous error estimate adds a notable supplement of resolution in the boundary layer close to the wall and imposes a finer resolution in the vicinity of the trailing edge. As expected, the viscous terms of the error estimate govern the mesh density in the boundary layer. These di↵erences strongly impact the drag estimation as confirmed by the convergence curves of the drag presented in Figure 3. We observe that our new adaptive method converges quicker to the reference drag value. For instance, the adapted mesh obtained with the viscous error estimate composed of 18 943 vertices achieves a drag value close to the adapted mesh obtained with the inviscid error estimate composed of 141 135 vertices. In other words, seven times less vertices are needed to achieve the same accuracy.

The accuracy of the computed drag values and the numerical convergence order are plotted in Figure 3 (right) and listed in Table 1. To compute the convergence order, we have chosen as reference drag the drag obtained on the finest adapted mesh (149 465 vertices) obtained with the viscous goal-oriented error estimate: C Re f D = 0.6648091. We notice a higher order of convergence with the viscous goal-oriented error estimates (⇡ 2.25) than with the inviscid goal-oriented error estimate (⇡ 1.75). The observed convergence order of the total drag is a reasonable value for this problem. We present a supersonic laminar viscous flow around the same NACA0012 airfoil starting from the same initial mesh as before. The free-stream Mach number considered here is 2, with an angle of attack of 10 degrees and a laminar Reynolds number of 106. A thick bow shock appears upstream of the airfoil profile as shown in Figure 5 where the Mach field is pictured. The close-up view of the airfoil shows a thick boundary layer and a thick wake for this low Reynolds number flow.

In order to perform the convergence study, several adapted meshes with di↵erent complexity have been obtained for the inviscid (no viscous terms) and the viscous goal-oriented error estimates targeting the total drag as output. Multiple views for the adapted meshes containing about 38 000 vertices are shown in Figure 7 for both approaches. These adapted meshes are obtained after 6 iterations of the mesh adaptation loop. We observe again that the viscous error estimate increases the mesh resolution in the boundary layer close to the airfoil and at the trailing edge, while the inviscid one add more resolution in the bow shock vicinity.

As the targeted output functional is the drag, meshes obtained with the new error estimates are closer to our expectation. This is confirmed by the drag and error convergence curves, see Figure 6 and Table 2. In this example too, an earlier convergence to the reference value is achieved with the new error estimate. For instance, the adapted mesh obtained with the viscous error estimate composed of 19 103 vertices achieves a similar drag value as the adapted mesh obtained with the inviscid error estimate composed of 140 924 vertices. In other words, seven times less vertices are needed to achieve the same accuracy. To analyze the convergence order, we consider as reference drag, the drag obtained on the finest mesh (149 514 vertices) with the viscous goal-oriented mesh adaptation: C Re f D = 0.5549227. Then, we observe a convergence order around 2.4 for the viscous approach and 1.8 for the inviscid one. This confirms the superiority of the viscous goal-oriented error estimates and the large influence of the viscous terms in the error estimates. It also reflects the optimality of the obtained meshes. We present a subsonic laminar viscous flow at a higher Reynolds number around the same NACA0012 airfoil still starting from the same initial mesh. The free-stream Mach number considered here is 0.5, with an angle of attack of 3 degrees and a Reynolds number of 5 000. This test case have been initially proposed by Swanson and Turkel [START_REF] Swanson | A multistage time-stepping scheme for the Navier-Stokes equations[END_REF] without any angle of attack to evaluate compressible flow solver. This particular case was chosen because it has a small amount of trailing edge separation and was a good test case to check the levels of dissipation being produced by the numerical scheme. Then, this test case has been generalized to several angles of attack. These test cases require a robust flow solver to converge the solution to machine zero. Indeed, it has been shown that a non robust numerical resolution leads to unsteady solutions due to a lack of numerical convergence. This problem has been deeply investigated by Swanson et al. [START_REF] Swanson | Steady-state laminar flow solutions for naca 0012 airfoil[END_REF] where the authors pointed out some possible reason for failing to obtain steady solutions. They mentioned the case presented by Taube et al. in [START_REF] Kroll | Adigma -a european initiative on the development of adaptive higherorder variational methods for aerospace applications[END_REF] where the angle of attack ↵ = 2 . The authors included adaptation in their calculation, but on a rather coarse mesh where the adaptive criterion pointed out the wake as the refinement region, resulting thus in an insu cient resolution in the vicinity of the separation point, which, as pointed out by Swanson et al., might led to unsteadiness. However, the studies conducted in [START_REF] Swanson | Steady-state laminar flow solutions for naca 0012 airfoil[END_REF] and on adapted meshes in [START_REF] Venditti | Anisotropic grid adaptation for functional outputs: application to two-dimensional viscous flows[END_REF] showed that this test case is indeed a steady problem.

The Mach contour solution is shown in Figure 8. The flow separates before mid-chord location on the upper surface of the airfoil. From the streamline patterns, we can observe primary and secondary recirculating regions, with the secondary region attached to the trailing edge.

As it has been pointed out in [START_REF] Swanson | Steady-state laminar flow solutions for naca 0012 airfoil[END_REF] the computed aerodynamic values variation should vanish with mesh refinement. We consider the total drag coe cient C tot D integrated over the airfoil profile as our quantity of interest and we perform convergence studies with three di↵erent error estimates: the classical Hessian-based error estimates controlling the L 2 -norm of the interpolation error of the local Mach number field [START_REF] Alauzet | High order sonic boom modeling by adaptive methods[END_REF][START_REF] Loseille | Achievement of global second-order mesh convergence for discontinuous flows with adapted unstructured meshes[END_REF], the inviscid goal-oriented error estimate and the viscous error estimate. For each estimate, we perform multiple adaptations at di↵erent mesh complexity (from 2 000 to 140 000). For each complexity, the final adapted mesh is obtained after 6 adaptive iterations.

The resulting adapted meshes composed of about 38 000 vertices are shown in Figures 10 and11. The top pictures represent the adapted mesh obtained with the Hessian-based error estimate. We observe that the remeshing e↵ort is very poor at the leading edge, at the trailing edge, and around the airfoil profile in the boundary layer, which is what expected from this method based on minimizing the Mach errors values and not directly on improving the total drag. We can also observe that with this method the wake is highly refined (see the second column of Figure 11). The middle pictures show the adapted mesh obtained with the inviscid goal-oriented error estimate. We observe that the leading edge, trailing edge and the boundary layer are more refined than using the Hessian-based error estimate while the wake has been adapted slightly. However, a closer look near the body shows that the boundary layer is not highly refine as this estimate does not take into account viscous e↵ects. The bottom pictures illustrate the adapted mesh obtained with the viscous goal-oriented error estimate. We notice that the leading edge, trailing edge and the boundary layer are even more refined and again less e↵ort has been put in the wake. A closer look near the body points out the increase refinement in the boundary layer region.

The conclusions regarding the comparison between the three adaptive methods are reinforced by Figure 9 and Table 3 where the total drag value convergence with increasing mesh size is given. We can see that the Mach-adapted method shows only a late convergence for meshes above 70 000 vertices. Both goal-oriented methods shows earlier convergence and again the viscous estimates proves to be superior to the other methods.

8.4. NACA0012 shock-boundary layer interaction M = 1.4, ↵ = 0 and Re = 1 000

The next example is a shock-boundary layer interaction problem in a rectangular domain of size [START_REF] Belme | Unsteady aerodynamic and adjoint method[END_REF][START_REF] Ph | Approximation by finite element functions using local regularization[END_REF]⇥ [START_REF] Belme | Unsteady aerodynamic and adjoint method[END_REF][START_REF] Belme | Unsteady aerodynamic and adjoint method[END_REF] containing a NACA0012 airfoil. The left, right and top parts of the domain are considered as inflow and outflow boundary conditions. The bottom part of the domain is the ground. It is considered as a slip surface except a region between x = 4 and x = 20 which is a no-slip surface where develops a boundary layer. The leading edge of the airfoil is located at (0, 0); this represents a distance of 10 chords length from the ground. The airfoil is also considered as a slip surface.

As the supersonic flow conditions are a free-stream Mach number of 1.4, an angle of attack of 0 and a Reynolds number of 1 000; the flow is strongly shock dominated. It is then interesting to investigate the interaction between the bow shock emitted by the airfoil and the boundary layer on the ground. For the observation, we consider the pressure deviation on a small part of no-slip ground = [START_REF] Becker | A feed-back approach to error control in finite element methods: basic analysis and examples[END_REF][START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF] ⇥ 0 (shown in Figure 12):

j(w) = 1 2 Z (p p 1 ) 2 p 1 d .
We run this simulation on a quasi-uniform mesh composed of about 42 000 vertices. The obtained solution is depicted in Figure 12 where the local Mach number is plotted. The supersonic airfoil generates a bow shock that propagates to the ground. This shock interacts with the reflected shock coming from the edge of the boundary layer and with the boundary layer. At first sight, on this quasi-uniform mesh there is no visible impact resulting from the shock-boundary layer interaction. Physical phenomena are spread out on this non-adapted mesh.

We have applied our goal-based adaptive method starting from a quasi-uniform mesh containing about 16 000 vertices. For this initial mesh, no particular boundary layer mesh has been considered at the ground, and we expect the mesh adaptive method to capture it and mesh it adequately to reduce the error committed on our scalar functional j. It is also interesting to observe the behavior of the goal-oriented mesh adaptation as an indirect physical phenomenon (the bow shock) is interacting with the boundary layer just before the observation area. To perform a convergence study, many complexities have been considered. Several pictures of the resulting adaptive mesh using the viscous goal-oriented error estimate composed of about 72 000 vertices are given in Figure 14. The first view (top pictures) illustrates the resulting adaptive mesh for the entire computational domain with the associated solution. We observe, as expected, a concentration of the refinement e↵ort in the lower half of the domain and to the ground meaning that the flow in these regions contributes the most to the errors committed on our observation. Indeed, even if the upper part of the shocks are refined by the algorithm, the mesh density in that region is a lot coarser than in the lower part of the domain as its impact on the observation is weak. We notice that shocks reflection to the ground are also stopped to be refined after their interactions with the tail shock. The last view (bottom pictures) shows a close-up view of the shock-boundary layer interaction region. Thanks to the mesh adaptation, we observe that we are able to capture accurately the separation bubble and the lift o↵ of the boundary layer due to the impact of the oblique shock, and we also capture the reflected shock in front of the bubble and the expansion fan at the foot of the oblique shock. All these waves are visible in the mesh.

We have also considered a series of adapted meshes for the inviscid error estimator and compared the convergence of the scalar functional with increasing number of vertices. The results are shown in Figure 13 and Table 4. We can see that with the exception of the first computations on a rather coarse mesh of 2 000 vertices, for both cases, the value of the scalar functional converges quickly and remains quite close to final value. Both approaches are almost superposed with increasing mesh size. From far, the adapted meshes obtained with the viscous and the inviscid goal-oriented error estimates look similar, indeed the viscous terms have a little impact in the far field. In fact, the discrepancies between the two goal-oriented adaptive methods for this problem are mostly visible in the boundary layer at the ground level. If we zoom in the boundary layer region near the observation area , we notice that the region close to the ground is more refined with the viscous error estimate as shown in Figure 15. The last example is a three-dimensional turbulent shock-boundary layer interaction problem. The geometry consists in a diamond shape of length 0.2 and thickness 0.02 immersed in a box of size [ 0.1, 1.25] ⇥ [ 0.05, 0.05] ⇥ [0, 0.5]. The leading edge is along the y-axis located at x = 0 and z = 0.2. The supersonic flow conditions are a free-stream Mach number of 1.4, an angle of attack of 0 and a Reynolds number of 3.6 ⇥ 10 6 . The boundary conditions of this test case are presented in Figure 16. The diamond geometry is creating shock waves which impact the flat plate located at the bottom of the domain. The flat plate is also impacted by shock waves reflected by the top wall. This configuration leads to a complex pattern of shock waves with many shock-boundary layer interactions.

The viscous goal-oriented error estimate has been considered to control the error on the total drag on the flat plate. The adaptive process starts from a non-adapted mesh composed of about 170, 000 vertices. This mesh contains a coarse boundary layer mesh with a first layer at y+ = 13 (i.e., h = 10 4 ). The solution on this initial mesh is depicted in Figure 16 (right). Six mesh adaptation iterations are considered to converge the couple mesh solution. The resulting adapted mesh containing 1.6 million vertices and the associated local Mach field number are shown in Figure 17. The complex shock waves pattern is accurately captured by the adaptive process in all the regions where it impacts the flat plate (we observe that at the end of the domain, the mesh has been coarsen). The boundary layer is also highly refined, and a close-up view of the shock boundary layer interaction region shows that we capture accurately the separation bubble and the lift o↵ of the boundary layer due to the impact of the oblique shock. We also capture the reflected shock in front of the bubble and the expansion fan at the foot of the oblique shock. All these waves are visible in the mesh.

Conclusion

This study has presented an error analysis for anisotropic goal-oriented mesh adaptation for the compressible Navier-Stokes equations.

The simplified feature-based/Hessian-based approach consisting in minimizing interpolation errors has early proved its efficiency and ease of use. Unfortunately, the application to systems is penalized by the arbitrary choice of adequate sensors (features). Indeed, the mathematical model is not taken into account and the genuine purpose of reducing the approximation error is not satisfactorily met. In order to reduce e ciently approximation errors, it is mandatory to use an error estimate. By construction, an a posteriori estimate does not directly involve information concerning the behavior of local error when the mesh is stretched. For that reason, we proposed to consider an a priori analysis of a goal-oriented error formulation in order to express in a natural way the error model in terms of interpolation errors. As in the Hessian-based method, a main advantage of this approach is that interpolation errors are then transformed to account for the coupled influence of mesh size and stretching. Moreover, an analytic expression of the optimal mesh is directly obtained thanks to the continuous mesh framework.

The key factor for success was to derive the a priori error estimate. We have first focused on an elliptic model. The extension of this analysis to the complete collection of terms in the compressible Navier-Stokes system was proposed. It re-uses consistently the elliptic error analysis and the previous work dedicated to the Euler system [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. The main result has been re-written in a more comprehensive form showing that the proposed error estimate is simply a sum of interpolation errors weighted by derivatives of the adjoint states and velocity components.

The proposed error estimate has been validated on a set of test cases dealing with steady external Navier-Stokes compressible flows with low Reynolds number at subsonic, transonic and supersonic regime. For all the test cases, a nice mesh convergence of the considered functional has been observed. Moreover, the viscous estimate has been compared to the inviscid one and the Hessian-based approach, and for each case has proven its superiority. This emphasizes the optimality of the viscous goal-oriented adapted mesh. The proposed mesh adaptation brings more fidelity on the output evaluation. The proposed approach has also been applied to a turbulent test case in three dimensions.

The proposed analysis has a good potential for the extension of the method to other (possibly more complex) physical models, in particular to turbulence models such as the Spalart-Allmaras on equation model or the Menter k ! SST two equations model. It illustrates the rather general impact of the a priori strategy, which potentially can be applied to numerical approximations enjoying k exactness: see [START_REF] Carabias | Analyse et adaptation de maillage pour des schémas non-oscillatoires d'ordre élevé[END_REF] for a similar study applied to an ENO approximation. Time discretization error is not considered in this study. Solving this question is not so important for the type of calculation that are shown in this paper, but can be of paramount impact in many other cases, since implicit time advancing is mandatory for most Navier-Stokes calculations. Then in that case, the user needs to prescribe the time step. In a future work, the authors plan to consider a space-time error analysis in the context of the proposed method. Then, the global fixed-point unsteady mesh adaptation algorithm will be a candidate for the extension to unsteady compressible Navier-Stokes [START_REF] Alauzet | Multi-scale anisotropic mesh adaptation for time-dependent problems[END_REF][START_REF] Belme | Time accurate anisotropic goal-oriented mesh adaptation for unsteady flows[END_REF].

Even if a test case has been shown, the most important next issue remains the extension to higher-Reynolds flows where meshing issues are still numerous. Indeed, this would require to address the generation of highly-stretched (aspect ratios of tens to hundreds of thousands) adapted meshes. This extension would use the results of this work, but needs also to solve several other important di culties: the metric-controlled generation of highly stretched meshes (in particular for surface meshing), the definition of a more accurate Hessian recovery, and a consistent interpolation of fields in the boundary layer, three subjects which we are currently studying.

Z ⌦ ⇧ h ⇢E r • ⇣ u.µ(ru) T ⌘ d⌦ = d X i=1 d X j=1 Z ⌦ µ ⇧ h ⇢E @ @x j u i @u i @x j ! d⌦ .
Finding an upper bound of T V 5 = T V 5 (⇧ h , W) T V 5 (⇧ h , ⇧ h W) comes to analyze the following error term:

e i j 5 = Z ⌦ µ ⇧ h ⇢E @ @x j u i @u i @x j ! d⌦ Z ⌦ µ ⇧ h ⇢E @ @x j ⇧ h u i @⇧ h u i @x j ! d⌦.
Again discarding boundary terms, we get by integration by parts:

e i j 5 = Z ⌦ µ @(⇧ h ⇢E ) @x j u i @u i @x j ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j ⇧ h u i @⇧ h u i @x j ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j u i @(u i ⇧ h u i ) @x j ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j (⇧ h u i u i ) @⇧ h u i @x j ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j u i @(u i ⇧ h u i ) @x j ! d⌦ Z ⌦ µ @(⇧ h ⇢E ) @x j (u i ⇧ h u i ) @u i @x j ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j (u i ⇧ h u i ) @u i @x j ! d⌦ Z ⌦ µ @(⇧ h ⇢E ) @x j (u i ⇧ h u i ) @⇧ h u i @x j ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j @ u i (u i ⇧ h u i ) @x j d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j (u i ⇧ h u i ) @(u i ⇧ h u i ) @x j ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j @ u i (u i ⇧ h u i ) @x j d⌦ + 1 2 Z ⌦ µ @(⇧ h ⇢E ) @x j @(u i ⇧ h u i ) 2 @x j d⌦ | {z } I i j
.

Term I i j is homogeneous to a squared interpolation error, thus it is an error term of order four which can be neglected. Then, applying Lemma 3.1 to the first term of the right hand-side, we obtain the following estimate:

| e i j 5 | K d Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u i ⇧ h u i | d⌦ .
Finally, we obtain the following upper bound:

| T V 5 | d K d d X i=1 Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u i ⇧ h u i | d⌦ . (A.6)
In contrast with the previous term, the gradient of the velocity vector is not transposed, thus velocity components will be crossed. We have:

T V 6 (⇧ h , W) = Z ⌦ ⇧ h ⇢E r • (u.µru) d⌦ = d X i=1 d X j=1 Z ⌦ µ ⇧ h ⇢E @ @x j u i @u j @x i ! d⌦.
Again, to find an upper bound of T V 6 = T V 6 (⇧ h , W) T V 6 (⇧ h , ⇧ h W) , we analyze the following error term:

e i j 6 = Z ⌦ µ ⇧ h ⇢E @ @x j u i @u j @x i ! d⌦ Z ⌦ µ ⇧ h ⇢E @ @x j ⇧ h u i @⇧ h u j @x i ! d⌦ .
After integrating by parts both integrals and neglecting the boundary terms, we obtain:

e i j 6 = Z ⌦ µ @(⇧ h ⇢E ) @x j u i @u j @x i ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j ⇧ h u i @⇧ h u j @x i ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j u i @(u j ⇧ h u j ) @x i ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j (⇧ h u i u i ) @⇧ h u j @x i ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j u i @(u j ⇧ h u j ) @x i ! d⌦ Z ⌦ µ @(⇧ h ⇢E ) @x j (u j ⇧ h u j ) @u i @x i ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j (u j ⇧ h u j ) @u i @x i ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x j (⇧ h u i u i ) @⇧ h u j @x i ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x j @ u i (u j ⇧ h u j ) @x i d⌦ | {z } I i j 1 + Z ⌦ µ @(⇧ h ⇢E ) @x j (u j ⇧ h u j ) @u i @x i (u i ⇧ h u i ) @⇧ h u j @x i ! d⌦ | {z } I i j 2
Lemma 3.1 can be applied to term I i j 1 leading to the upper bound:

|I i j 1 | K d Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u j ⇧ h u j | d⌦ .
Term I i j 2 term can be further analyzed, and re-written as follows:

I i j 2 =
Z ⌦ µ @(⇧ h ⇢E ) @x j (u j ⇧ h u j ) @u i @x i (u i ⇧ h u i ) @u j @x i (u i ⇧ h u i ) @(⇧ h u j u j ) @x i ! d⌦ .

The last term of the above relation is a product of interpolation error, thus it is an error term of order four which is neglected. For i = j, the first and the second terms cancel each other out, the remaining error term is of fourth order error term, hence I ii 2 is neglected. For i , j, only the crossed velocity components terms contribute to the error estimation. Therefore, we get the following upper bound:

| T V 6 | = d X i=1 d X j=1 I i j 1 + d X i=1 d X j=1 j,i I i j 2 K d 0 B B B B B B @ d X i=1 d X i= j Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u j ⇧ h u j | d⌦ 1 C C C C C C A + d X i=1 d X j=1 j,i Z ⌦ µ @(⇧ h ⇢E ) @x j (u j ⇧ h u j ) @u i @x i (u i ⇧ h u i ) @u j @x i ! d⌦ . (A.7)
In the above above relation, the second term can be reformulated more conveniently as interpolation error on the velocity terms weighted by cross products of the gradients @(⇧ h ⇢E ) @x j with @u i @x k . To this end, we introduce the following vector:

ru i ⇥ r(⇧ h ⇢E ) = ! u i = ⇣ ! u i ,x , ! u i ,y , ! u i ,z ⌘ T .
In two dimensions, we obtain: This conclude the error estimation of the sixth term.

Appendix A.6. Study of viscous term T V

7

The same remark as previously about crossed velocity components still hold for the seventh term because of multiplication by the identity matrix. This term writes:

T V 7 (⇧ h , W) = 2 3 Z ⌦ ⇧ h ⇢E r • (u.µr.uI) d⌦ = 2 3 d X i=1 d X j=1 Z ⌦ µ ⇧ h ⇢E @ @x i u i @u j @x j ! d⌦.
To obtain an upper bound of T V 7 = T V 7 (⇧ h , W) T V 7 (⇧ h , ⇧ h W) , we analyze the following error term:

e i j 7 = Z ⌦ µ ⇧ h ⇢E @ @x i u i @u j @x j ! d⌦ Z ⌦ µ ⇧ h ⇢E @ @x i ⇧ h u i @⇧ h u j @x j ! d⌦ .
After integrating by parts both integrals and neglecting the boundary terms, we obtain:

e i j 7 = Z ⌦ µ @(⇧ h ⇢E ) @x i u i @u j @x j ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x i ⇧ h u i @⇧ h u j @x j ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x i u i @(u j ⇧ h u j ) @x j ! d⌦ + Z ⌦ µ @(⇧ h ⇢E ) @x i (⇧ h u i u i ) @⇧ h u j @x j ! d⌦ = Z ⌦ µ @(⇧ h ⇢E ) @x i @(u i (u j ⇧ h u j )) @x j d⌦ | {z } I i j 1 + Z ⌦ µ @(⇧ h ⇢E ) @x i (u j ⇧ h u j ) @u i @x j (u i ⇧ h u i ) @⇧ h u j @x j ! d⌦ | {z } I i j 2 .
Similarly to term T V 6 , we can directly apply Lemma 3.1 to term I i j 1 , thus the following estimation holds:

|I i j 1 | K d Z ⌦ µ | ⇢ H ( ⇢E )| |u i | |u j ⇧ h u j | d⌦ .
Regarding the integral I i j 2 , the following equality holds:

I i j 2 =
Z ⌦ µ @(⇧ h ⇢E ) @x i (u j ⇧ h u j ) @u i @x j (u i ⇧ h u i ) @u j @x j (u i ⇧ h u i ) @(⇧ h u j u j ) @x j ! d⌦ .

As previously, for i = j, the first and the second terms cancel each other out. The remaining term is a fourth order error term because it is the product of interpolation errors. Therefore, I ii 2 can be neglected. For i , j, the last term is neglected as it is a fourth order error term, then only the crossed velocity components terms contribute to the error estimation. We have the following upper bound: µ @(⇧ h ⇢E ) @x i (u j ⇧ h u j ) @u i @x j (u i ⇧ h u i ) @u j @x j ! d⌦ (A.8)

| T V 7 | =
The second term can be reformulated more conveniently using the following vector:

ru i ⇥ r(⇧ h ⇢E ) = ! u i = ⇣ ! u i ,x , ! u i ,y , ! u i ,z ⌘ T .

Figure 1 :

 1 Figure 1: Super-convergent molecule for a vertical derivative. Diamond decomposition of the mesh over each edge in 2D (left) and each face in 3D (right). Now, we can provide a first upper bound:

  (right). If we denote by | • | the measure (length, area, volume) of geometric entity and h the height of a simplex K, by definition, we have the following relations for the area in 2D:

  ! d⌦ , and K d = 3 in two dimensions, K d = 6 in three dimensions, and A B holds for a majoration asymptotically valid, i.e. A  B + o(A) when mesh size tends to zero. Expression | ⇢ H (') | holds for spectral radius of the Hessian of '.

8. 1 .

 1 Laminar transonic flow around NACA0012: M = 0.8, ↵ = 10 and Re = 73

8. 2 .

 2 Laminar supersonic flow around NACA0012: M = 2, ↵ = 10 and Re = 106

8. 3 .

 3 Laminar subsonic flow around NACA0012: M = 0.5, ↵ = 3 and Re = 5 000

Figure 2 :

 2 Figure 2: NACA0012 at M = 0.8, ↵ = 10 and Re = 73: Local Mach number solution field and iso-contours for the transonic laminar viscous flow obtained with the goal-oriented viscous adapted mesh composed of 37 838 vertices.

Figure 3 :

 3 Figure 3: NACA0012 at M = 0.8, ↵ = 10 and Re = 73: Convergence of the total drag with respect to the number of vertices (left) and error convergence (right) for the transonic laminar viscous flow computed with the inviscid (red curves) and the viscous (blue curves) goal-oriented anisotropic mesh adaptation.

Figure 4 :

 4 Figure 4: NACA0012 at M = 0.8, ↵ = 10 and Re = 73: Comparison of adapted meshes composed of about 38 000 vertices for the transonic laminar viscous flow. Left, goal-oriented adapted mesh without viscous flux contribution to optimal metric (only the Euler fluxes criterion). Right, goal-oriented adapted mesh with viscous flux contribution to optimal metric.

Figure 5 :

 5 Figure 5: NACA0012 at M = 2, ↵ = 10 and Re = 106: Local Mach number solution field and iso-contours for the supersonic laminar viscous flow obtained with the goal-oriented viscous adapted mesh composed of 37 842 vertices.

Figure 6 :

 6 Figure 6: NACA0012 at M = 2, ↵ = 10 and Re = 106: Convergence of the total drag with respect to the number of vertices (left) and error convergence (right) for the transonic laminar viscous flow computed with the inviscid (red curves) and the viscous (blue curves) goal-oriented anisotropic mesh adaptation.

Figure 7 :

 7 Figure 7: NACA0012 at M = 2, ↵ = 10 and Re = 106: Comparison of adapted meshes composed of about 38 000 vertices for the supersonic laminar viscous flow. Left, goal-oriented adapted mesh without viscous flux contribution to optimal metric (only the Euler fluxes criterion). Right, goal-oriented adapted mesh with viscous flux contribution to optimal metric.

Figure 8 :

 8 Figure 8: NACA0012 at M = 0.5, ↵ = 3 and Re = 5 000: Local Mach number solution field and iso-contours obtained on the goal-oriented viscous adapted mesh composed of 37 353 vertices.

Figure 9 :

 9 Figure9: NACA0012 at M = 0.5, ↵ = 3 and Re = 5 000: Convergence of the total drag with respect to the number of vertices for the subsonic laminar viscous flow computed with Mach Hessian-based error estimate (pink curve), the inviscid (red curves) and the viscous (blue curves) goal-oriented error estimates.

Figure 10 :

 10 Figure 10: NACA0012 at M = 0.5, ↵ = 3 and Re = 5 000: Comparison of adapted meshes composed of about 38, 000 vertices for the subsonic laminar viscous flow. Top, Hessian-based error estimates controlling the L 2 -norm of the interpolation error of the local Mach number field. Middle, goal-oriented adapted mesh without viscous flux contribution to optimal metric (only the Euler fluxes criterion). Right, goal-oriented adapted mesh with viscous flux contribution to optimal metric.

Figure 11 :

 11 Figure 11: NACA0012 at M = 0.5, ↵ = 3 and Re = 5 000: Comparison of adapted meshes composed of about 38, 000 vertices for the subsonic laminar viscous flow. Top, Hessian-based error estimates controlling the L 2 -norm of the interpolation error of the local Mach number field. Middle, goal-oriented adapted mesh without viscous flux contribution to optimal metric (only the Euler fluxes criterion). Right, goal-oriented adapted mesh with viscous flux contribution to optimal metric.

8. 5 .

 5 Three dimensional shock-boundary layer interaction M = 1.4, ↵ = 0 and Re = 3.6 ⇥ 10 6

Figure 12 :Figure 13 :

 1213 Figure 12: NACA0012 at M = 1.4, ↵ = 0 and Re = 1 000: Local Mach number solution field and iso-contours for the shock-boundary layer interaction problem obtained on a quasi-uniform mesh of about 42, 000 vertices. Left, visualization of the observation region on the ground.

Figure 14 :

 14 Figure 14: NACA0012 at M = 1.4, ↵ = 0 and Re = 1 000: Adapted meshes (left) composed of about 72 000 vertices and associated Mach number field (right) obtained using the goal-oriented viscous error estimate. From top to bottom, global view of the domain, close-up view of the airfoil and zoom in the boundary layer region where the shock-boundary layer interaction occurs.

Figure 15 :Figure 16 :

 1516 Figure 15: NACA0012 at M = 1.4, ↵ = 0 and Re = 1 000: Comparison of adapted meshes composed of about 72 000 vertices in the boundary layer region. Left, goal-oriented adapted mesh without viscous flux contribution to optimal metric (only the Euler fluxes criterion). Right, goal-oriented adapted mesh with viscous flux contribution to optimal metric.

Figure 17 :

 17 Figure 17: Diamond at M = 1.4, ↵ = 0 and Re = 3.6 ⇥ 10 6 : Adapted meshes (left) composed of about 1.6 million vertices and associated Mach number field (right) obtained using the goal-oriented viscous error estimate. From top to bottom, global view of the domain surface, global view of the domain with a cut plane at y = 0, and close-up view in the boundary layer region where the shock-boundary layer interaction occurs.

|u 1 ⇧|u 1 ⇧

 11 h u 1 | |! u 2 ,z | + |u 2 ⇧ h u 2 | |! u 1 ,z | ⌘ d⌦ ,and in three dimensions, we get:h u 1 | |! u 3 ,y ! u 2 ,z | + |u 2 ⇧ h u 2 | |! u 1 ,z ! u 3 ,x | + |u 3 ⇧ h u 3 | |! u 2 ,x ! u 1 ,y | ⌘ d⌦ .

  H ( ⇢E )| |u i | |u j ⇧ h u j | d⌦

Table 1 :

 1

	Viscous Error Estimate		Inviscid Error Estimate
	Mesh size 2 384 0.6564739 8.3352 ⇥ 10 3 C tot D Error 4 840 0.6610405 3.7686 ⇥ 10 3 9 562 0.6630579 1.7512 ⇥ 10 3 18 943 0.6639960 8.1310 ⇥ 10 4 37 838 0.6644963 3.1280 ⇥ 10 4 75 241 0.6647139 9.5200 ⇥ 10 5 149 465 0.6648091 Reference	Conv. order -2.24 2.25 2.24 2.76 3.46 -	Mesh size 2 330 0.6394387 2.5370 ⇥ 10 2 C tot Error D 4 631 0.6511749 1.3634 ⇥ 10 2 9 139 0.6576796 7.1295 ⇥ 10 3 17 987 0.6608486 3.9605 ⇥ 10 3 35 881 0.6625947 2.2144 ⇥ 10 3 71 038 0.6635034 1.3057 ⇥ 10 3 141 135 0.6640958 7.1330 ⇥ 10 4	Conv. order -1.81 1.91 1.74 1.68 1.55 1.76

NACA0012 at M = 0.8, ↵ = 10 and Re = 73: Total drag coe cient convergence, error convergence and convergence order for the viscous and inviscid error estimates.

Table 2 :

 2

NACA0012 at M = 2, ↵ = 10 and Re = 106: Total drag coe cient convergence, error convergence and convergence order for the viscous and inviscid error estimates.

Table 3 :

 3 NACA0012 at M = 0.5, ↵ = 3 and Re = 5 000: Total drag coe cient convergence for the viscous, the inviscid, and the Hessian-based error estimates.

	Viscous Error Estimate Mesh size C tot D 2 397 0.05829608 4 772 0.05772547 9 473 0.05768550 18 861 0.05770964 37 353 0.05773110 73 910 0.05774073 147 257 0.05774982 293 706 0.05775450	Inviscid Error Estimate Mesh size C tot D 2 412 0.05819562 4 739 0.05756836 9 342 0.05757901 18 479 0.05764147 36 487 0.05767100 71 848 0.05769866 142 969 0.05772853	Mach Hessian L 2 Error Estimate Mesh size C tot D 2 343 0.06264318 4 662 0.05826357 9 303 0.05836721 18 337 0.05751269 35 956 0.05752201 71 395 0.05740691 139 523 0.05749961

Table 4 :

 4

	Viscous Error Estimate	Inviscid Error Estimate
	Mesh size	C tot D	Mesh size	C tot D
	2 319	0.004462248	2 324	0.004258559
	4 601	0.004662791	4 598	0.004644956
	9 200	0.004678465	9 187	0.004652647
	18 107	0.004678565	18 476	0.004674319
	36 090	0.004675225	36 104	0.004672913
	71 887	0.004670773	71 822	0.004671156

NACA0012 at M = 1.4, ↵ = 0 and Re = 1 000: Pressure functional j convergence for the viscous and inviscid goal-oriented error estimates.

Acknowledgements

The authors want to acknowledge Adrien Loseille for providing software Feflo.a. They also acknowledge Loïc Frazza for his constructive comments on the error estimate.

This work has been supported by French National Research Agency (ANR) through COSINUS program (project Ecinads n o ANR-09-COSI-003) and Maidesc (ANR-13-MONU-0010). HPC resources from GENCI-[CINES] are also gratefully acknowledged.

Appendix A. Proof of Proposition 5.1

In this section, we demonstrate the result of Proposition 5.1. Thanks to Relation [START_REF] Frey | Anisotropic mesh adaptation for CFD computations[END_REF], to obtain an error estimate of j, we aim at giving an upper bound of:

where h is composed of (see Relation [START_REF] Formaggia | Anisotropic mesh adaptation in computational fluid dynamics: Application to the advection-di↵usion-reaction and the Stokes problems[END_REF]) the Euler term, the Euler boundary term, the viscous term, and the stabilization term. In the following, we neglect the boundary term and the stabilization term thanks to smoothness of functions W and . The analysis of the boundary term is given in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF].

Appendix A.1. Study of the inviscid term Thanks to Relation [START_REF] Ph | Approximation by finite element functions using local regularization[END_REF], the error for the Euler term (without boundary term) writes:

Here, we give an estimation of the above term with a more concise analysis than the one proposed in [START_REF] Loseille | Fully anisotropic goal-oriented mesh adaptation for 3D steady Euler equations[END_REF]. First, we add and subtract a term:

Assuming smoothness of and F E (W), we deduce that on ⌦, interpolation errors are error terms of order two and interpolation error gradients are error terms of order one. Therefore, the first term of the right-hand side is an error term of order three and can be thus neglected. The second term of the right-hand side can be transformed into:

Again, discarding the boundary integrals, we finally get:

For the viscous term, we consider:

which is the sum of seven terms, see Relation [START_REF] Cournède | Positivity statements for a Mixed-Element-Volume scheme on fixed and moving grids[END_REF]. We first focus on the first three terms of this summation coming from the momentum equation:

with ⇢u = ( ⇢u 1 , ⇢u 2 , ⇢u 3 ) T . We remark that T V 2 and T V 3 expressions can be directly added with an exchange of i and j derivatives (i.e., @ 2 u l @x i @x j = @ 2 u l @x j @x i ). Then, the sum of these terms writes:

To give a bound of Relation (A.3), we analyse:

It is su cient to consider the general expression e i jkl to analyze all the previous terms:

After an integration by part and neglecting the boundary terms, the error term writes:

Now, according to Lemma 3.1 from the elliptic error analysis such a volume contribution is overestimated as:

with K d = 3 in two dimensions and K d = 6 in three dimensions. Finally, going back to our initial summation, the following a priori estimate holds for the first three terms of viscous flux contribution:

If we split the second term of the right hand-side into a sum with i = j and a sum with i , j, the previous bound can be re-written:

We analyse the first term coming from the energy equation:

We focus on the error term:

After an integration by parts and by neglecting the boundary term as mentioned above, we get:

Then, Lemma 3.1 can be applied and we obtain the following estimate:

However, this error analysis cannot be directly applied to the next three remaining terms, T V 5 , T V 6 and T V 7 , because of their non-linearity. A slightly di↵erent analysis is considered in which higher order error terms, that can be neglected, appear. The T V 5 viscous term can be written under the following development:

In two dimensions, we obtain:

and in three dimensions, we get:

In fact, we can get a sharper estimate using the following observations on T V 6 and T V 7 terms before considering the absolute value into the sum: Then, we deduce that:

and, again, we can re-formulate the second right hand-side term in a more convenient way.

This concludes the proof of Proposition 5.1, Corollary 5.1 and Corollary 5.2.