R Abgrall
email: remi.abgrall@inria.fr

D De Santis
email: dante.desantis@inria.fr

Ph.D Candidate

High order residual distribution scheme for Navier-Stokes equations

In this work we describe the use of the Residual Distribution schemes for the discretization of the conservation laws. In particular, emphasis is put on the construction of a third order accurate scheme. We first recall the proprieties of a Residual Distribution scheme and we show how to construct a high order scheme for advection problems, in particular for the system of the Euler equations. Furthermore, we show how to speed up the convergence of implicit scheme to the steady solution by the means of the Jacobian-free technique. We then extend the scheme to the case of advection-diffusion problems. In particular, we propose a new approach in which the residuals of the advection and diffusion terms are distributed together to get high order accuracy. Due to the continuous approximation of the solution the gradients of the variables are reconstructed at the nodes and then interpolated on the elements. The scheme is tested on scalar problems and is used to discretize the Navier-Stokes equations.

I. Introduction

We are interested in the higher order approximation of the multidimensional compressible Navier-Stokes equations. The discretization of the equations is obtained by the means of the Residual Distribution (RD) schemes which have already been shown to be a robust higher order approximation technique for the Euler equations. Residual Distribution schemes represent a very interesting alternative to Discontinuous Galerkin (DG) schemes. [START_REF] Cockburn | Discontinuous Galerkin methods: theory, computation and application[END_REF] While computationally compact and probably more flexible, DG schemes suffer from the serious drawback of a very fast growth of the number of degrees of freedom with the cell polynomial degree. In RD schemes the formulation remains local, as in DG, but the number of degree of freedom growth less quickly. The price to pay is to impose a continuous approximation of the solution, even though some papers report their extension to discontinuous approximation. [START_REF] Abgrall | Development of residual distribution schemes for the discontinuous Galerkin method: the scalar case with linear elements[END_REF][START_REF] Abgrall | A residual distribution method using discontinuous elements for the computation of possibly non smooth flows[END_REF][START_REF] Hubbard | A framework for discontinuous fluctuation distribution[END_REF] Results for RD schemes in the case of order more that two have been presented for the system of the Euler equations in, [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] but the high order discretization of the Navier-Stokes equation is still an open question.

To extend RD scheme to convection-diffusion equations, efforts were made to properly discretize and distribute the diffusion terms. One kind of solution is to use a RD scheme for the convection terms while adding a Galerkin discretization of the diffusion terms. Nishikawa and Roe [START_REF] Nishikawa | On high order fluctuation splitting schemes for Navier-Stokes equations[END_REF] pointed out that the Galerkin discretization for the diffusion terms would lead to a loss of order of accuracy, according to truncation error analysis. They propose to write the governing equations as a first order system to get high order accuracy. Another approach consists in rewriting the RD scheme as a perturbation of the Galerkin scheme weighted by a properly scaled function, so that the scheme gives a correct behavior for the convection or diffusion dominated regimes. The scaling depends on a cell Peclet number.

A natural way to discretise the viscous terms would be to include it in the residual and distribute together with the convective part. This approach has been already used in framework of the DG, that can easily handle the discontinuous nature of the viscous flux across the element boundaries. In this work the idea of distribute the advective and the diffusion term is retrieved and rearranged in the contest of a continuous approximation. Since our formulation requires a continuous approximation, the gradients, which appear in the viscous terms, are reconstructed on each nodes and then interpolated on each elements by the same Lagrangian functions used for the solution. In this way, it is possible to compute the residuals of the advection and diffusion terms and distribute them together.

The paper is organized as follows. In the section II, we present the basic principle of the RD scheme, and the accuracy properties of the scheme for the scalar advection problem. In the section III, the scheme proposed for the scalar equation is extended to the case of the Euler equations. In the section IV we present an implicit iterative solver for the solution of steady non linear problems in combination of the Jacobian-free technique to speed up the convergence rate to steady state. In the section V, we extend the RD scheme to the case of the advection-diffusion problem, and in the section VI we use the proposed approach to discretize scalar problems and to verify the order of accuracy of the scheme. Eventually we use the RD scheme to obtain the solution of the Navier-Stokes equation on the flat plane and over the NACA 0012 airfoil. In the last section, we give concluding remarks.

II. High order RD schemes for hyperbolic conservation laws

Let us consider the multidimensional, steady, hyperbolic scalar equation for conservations laws in the form

∇• f (u) = 0, ∀ x ∈ Ω ⊂ R d , d = 2, 3 (1)
in which u(x) ∈ R and f ∈ R d is the flux function associated to the unknown u. The equation ( 1) must be supplemented by proper inflow boundary conditions

u| ∂Ω -= g(s), s ∈ ∂Ω -,
where Let us consider, now, a tessalation T h of the domain Ω, with h T a characteristic length of the element T of the mesh, the number of the total elements in the domain is N T . In the RD scheme the degrees of freedom (DOF) are associated to the points of the mesh and not to control volumes as in the Finite Volume or DG methods; we denote by {σ l } l=1,N T dof the list of the DOF on each element. In this work the elements considered are triangles and tetrahedrons, in two and three spatial dimensions respectively.

∂Ω -= {x ∈ ∂Ω | n • ∇ u f < 0},
When a linear interpolation of the solution is used, the degrees of freedom of each element coincide with the vertices. To have higher order interpolation of the solution, necessary to construct an high order RD scheme, extra degrees of freedom must me considered on each elements. In order to keep a local formulation, the extra degrees of freedom are added inside the element. When standard Lagrange elements P are used, a quadratic reconstruction of the solution is obtained taking as degrees of freedom the vertices of the element and the edges mid-points. The use of degrees of freedom out of each element 2, 3 destroys the compactness of the computational stencil, resulting in a not computationally efficient scheme. Furthermore the continuity of the standard Lagrange elements requires that all the DOF on the element boundaries are shared by neighboring elements, this results in a number of DOF smaller than the DG scheme. Clearly in the case of the continuous approximation the number of DOF increase less rapidly than in the discontinuous one and both cases become asymptotically similar.

A residual distribution scheme for the Eq. ( 1) reads

T ∋σ Φ T σ = 0, ∀ σ ∈ T, (2)
where the residuals Φ T σ must satisfy the following conservation constraint

4 σ∈T Φ T σ = T λ • ∇u dx = ∂T f h (u h ) • n dℓ = Φ T , ∀ T, (3)
where λ = ∇ u f and f h is an approximation of f . For example, it is possible to construct f h (u h ) as the Lagrange interpolant of f (u) at the degrees of freedom, as well as the flux function evaluated for the Lagrange interpolant of u. We assume that the residuals Φ T σ depend continuously on the values of {u σ } σ∈T . If σ is a DOF belonging to ∂Ω -, the boundary conditions must be taken into account while writing the Eq. ( 3). Indicating with Γ the edge (or the face) of ∂Ω -, we consider a numerical flux F , which depends on the boundary condition g, the inward normal n -to the boundary, and the local state u h . We define the boundary residuals which satisfy the following conservation relation as

σ∈Γ Φ Γ σ = Γ F (u h , g, n -) -f h (u h ) • n dℓ = Φ Γ σ , ∀ Γ ⊂ ∂Ω -. (4)
The Eq. ( 3) with the contribution of the boundary conditions writes

σ∈T Φ T σ + Γ⊂∂Ω -, Γ∋σ Φ Γ σ = 0.
It can be showed [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF] that if the sequence u h is bounded in L ∞ when h → 0 and if exists v, such that u h → v when h → 0, then v is a weak solution of (1). In the proof of this statement one has to assume the continuity of the interpolant across the edges, although this constrains may be alleviated and is possible to define RD schemes on discontinuous elements. 5-7

A. Accuracy constraints

Following the work of, [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF] we report now considerations about the accuracy of the scheme introduced in the previous sub-section. We define the following truncation error, for any smooth function ϕ

E(u h , ϕ h ) = σ∈Ω ϕ(σ)   σ∈T Φ T σ + Γ⊂∂Ω -, Γ∋σ Φ Γ σ   ,
with ϕ h the interpolant of ϕ. The scheme is k-th order accurate if the truncation error is O(h k ) when u h is an interpolant of the exact solution, assumed smooth enough. We have the following result Proposition 1 If the solution u is smooth enough and u h is its P k interpolant, the residuals satisfy

Φ T σ = O(h k+d ) and Φ Γ σ = O(h k+d-1 ), (5)
and if the approximation f h (u h ) is accurate with the order k + 1, then the truncation error satisfies

|E(u h , ϕ h )| ≤ C(ϕ, f , u)h k+1 ,
with C a constant which depends only on ϕ, f , and u.

Under of the hypotheses of the proposition 1, can be showed that

Φ T = O(h k+d ) and Φ Γ = O(h k+d-1 ),
and if exists a constant (in the scalar case) or a matrix (in the system case) β T σ such that

Φ T σ = β T σ Φ T and Φ Γ σ = β T σ Φ Γ ,
then the conditions (5) are satisfied provided that β T σ is uniformly bounded. Such a condition is historically called linearity preserving.

From a numerical point of view the integrals in (3) and ( 4) are computed using quadrature formulas, namely

Φ T = e∈∂K |e| NG p=1 ω p f h (u h ) • n e , ∀ T, and Φ Γ = |Γ| NG p=1 ω p F (u h , g, n -) -f h (u h ) • n e , ∀ Γ ⊂ ∂Ω -,
where e is the edge (or the face) of the element, with |e| the measure of e and n e the outward normal versor to e. We have indicated with N G and ω p respectively the number and the weights of the Gauss points used in the numerical quadrature. The quadrature formulas must be chosen such that the order of accuracy is not spoiled, this give us a criterion for choosing appropriate quadrature formulas since they must be such that

e∈∂K |e| NG p=1 ω p f h (u h ) • n e = ∂T f h (u h ) • n dℓ + O(h k+d ),
and

|Γ| NG p=1 ω p F (u h , g, n -) -f h (u h ) • n e = Γ F (u h , g, n -) -f h (u h ) • n dℓ + O(h k+d-1 ).
The practical approach used in this work is to reconstruct in each element a polynomial flux based on the Lagrange interpolation of the flux values evaluated at the degrees of freedom. The quadrature points coincide with the DOF, and the quadrature weights are easily computed once and for all.

B. Getting high order accuracy and monotonicity preservation

In this sub-section we introduce the conditions that must be satisfied by the numerical scheme in order to have non-oscillatory solutions. We start by re-writing the residuals as

Φ T σ = σ ′ ∈T c T σσ ′ (u σ -u σ ′ ), (6)
using this formalism, the Eq. ( 2) becomes, neglecting the boundary, conditions

T ∋σ σ ′ ∈T c T σσ ′ (u σ -u σ ′ ) = 0.
In general, the coefficients c T σσ ′ depend on the solution, which means that the last expression defines a set of non linear equations that must be solved by an iterative procedure, as for example a Jacobi-like iterations

u n+1 σ = u n σ -δ σ T ∋σ σ ′ ∈T c T σσ ′ (u σ -u σ ′ ) n , (7)
with δ σ a relaxation parameter. If the scheme satisfies the following positivity conditions

T ∋σ σ ′ ∈T c T σσ ′ ≥ 0 ∀σ, σ ′ and 1 -δ σ T ∋σ σ ′ ∈T c T σσ ′ ≥ 0 ∀σ, (8)
then the solution verifies the following discrete maximum principle min

T ∋σ min σ ′ ∈T u 0 σ ′ ≤ u n σ ≤ max T ∋σ max σ ′ ∈T u 0 σ ′ ,
with u 0 σ ′ the value of the initial solution at the DOF points. A more convenient approach is to replace the conditions (8) by a local positivity conditions

c σσ ′ ≥ 0 ∀σ, σ ′ ∈ T and δ σ max T ∋σ C σ C T σ σ ′ ∈T c T σσ ≤ 1 ∀σ, (9)
where the following median dual area has been introduced

C T σ = |T | N T DOF and C σ = T ∋σ C T σ .
A scheme that verifies the conditions ( 9) is said monotonicity preserving. These conditions do not imply that the iterative scheme ( 7) is convergent, but only that the maximum principle is satisfied, i.e., the L ∞ -stability.

It is well know from the Godunov's theorem that a monotonicity preserving scheme with the coefficients c T σσ ′ that do not depend on the solution can not be linearity preserving. [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for steady problems[END_REF] As a consequence, a monotonicity and linearity preserving scheme must be non linear.

There is a systematic way of constructing a non linear scheme which is both linearity and monotonicity preserving. We start from a monotone, first order, scheme

ΦT σ = σ ′ ∈T ĉT σσ ′ (u σ -u σ ′ ).
The coefficients ĉT σσ ′ are all positive and

σ ′ ∈T ΦT σ = Φ T .
If we denote with Φσ the high order residual, it also verifies the relation

σ ′ ∈T ΦT σ = Φ T and ΦT σ = β T σ ∂T f h (u h ) • n dℓ,
with β T σ uniformly bounded and σ∈T β T σ = 1 for the conservation. We introduce, now, the parameter

χ T σ = ΦT σ Φ T with σ∈T χ T σ = 1
, still for the conservation. We write now the high order residual as

ΦT σ = ΦT σ ΦT σ Φσ = σ ′ ∈T ΦT σ ΦT σ ĉT σσ ′ (u σ -u σ ′ ) the modified coefficients čT σσ ′ = ΦT σ ΦT σ ĉT σσ ′ ,
characterize the behavior of the high order scheme and it is easy to see that in order to fulfill the monotonicity preserving condition it is necessary that

ΦT σ ΦT σ ≥ 0, since ĉT σσ ′ ≥ 0.
The previous considerations about the linearity and monotonicity preserving can be rephrased as

σ ′ ∈T β T σ = 1 and σ ′ ∈T χ T σ = 1 conservation χ T σ β T σ ≥ 0, ∀σ ∈ T monotonicity preservation
These relation can be interpreted geometrically 4 and a simple choice, already designed in 9 in the second order case, consists in taking

β T σ = χ T + σ ′ ∈T χ T + σ ′ , χ T + σ = max(χ T σ , 0),
in the previous equation there is no singularity since

σ∈T χ T + σ = σ∈T χ T σ - σ ′ ∈T χ T - σ ≥ 1.

C. Construction of a high order RD scheme

In this sub-section we show how to construct a high order non upwind RD scheme from a first order scheme. The starting point is the Rusanov's scheme a

Φ T σ = Φ T N T DOF + α T (u σ -ū), (10)
with ū = σ∈T uσ N T DOF and α T a parameter large enough to garanty the stability of the scheme. Using the P k interpolant on the element T , u h = σ∈T u σ ψ σ and the residual reads

Φ T = T ∇ • f (u h ) dx = T λ • ∇u h dx = σ∈T u σ T λ • ∇ψ dx.
We define now k T σ = T λ • ∇ψ dx and noting that the Eq. ( 10) can be put in the form of the Eq. ( 6) with

c T σσ ′ = k T σ -α T N T DOF , it is easy to see that the condition c T σσ ′ ≥ 0 is satisfied if α T ≥ max σ∈T |k T σ |.
The scheme is extremely dissipative, but it is very cheap and simple to code and can be easily extended to system case. The high order scheme is constructed from the Eq. ( 10) applying the limitation technique described in the sub-section B.

The use of a central scheme, like the Rusanov's scheme, in combination with the limiting technique may produce a local downwind scheme which results in not accurate and not convergent scheme. The problem has been analyzed 10 and the solution proposed consists in modifying the scheme as follows

Φ ⋆,T σ = ΦT σ + h k T (λ • ∇ψ σ ) τ λ • ∇u h dx, τ > 0. ( 11
)
The last term on the second member of the previous equation is a streamline dissipation term, used in SUPG schemes to suppress the spurious mode of the Galerkin scheme. [START_REF] Brooks | Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations[END_REF] The formal accuracy of the scheme is preserved since the filtering term vanishes when u h is replaced with the exact solution It is worth noting that the conservation property is still preserved because σ∈T Φ ⋆ σ = Φ T since σ∈T ∇ψ σ = 0. Experimentally, we can see that the non oscillatory properties of the scheme are not spoiled.

In the construction of the filtering term a computation of an integral is necessary, this is usually done be the means of a quadrature formula. The use of a consistent quadrature formula for evaluating the integral in (11) may be quite expensive. It has been pointed out [START_REF] Abgrall | A simple construction of very high order non-oscillatory compact schemes on unstructured meshes[END_REF] that a consistent quadrature formula to compute the integral in (11) is not necessary to preserve the order of accuracy of the scheme since in choosing the quadrature formula one should only guaranties that the number of the quadrature points is enough to represent the gradients and the quadrature weights are positive to assure that the filtering term is dissipative.

a Other examples can be considered, such as the rephrasing of standard finite volume schemes in term of RD schemes.

Numerical verification of the accuracy

In order to evaluate the accuracy of the scheme, we perform a convergence study on the advection problem of the form λ • ∇u = 0 on Ω = [0, 1] 2 , with λ = (x, -y) T and on the inlet boundaries the following boundary condition is applied

g(x) = cos 2 (2πx) if x ∈ [0.25, 0.75] 0 else.
The problem corresponds to the rigid rotation of the profile at the inlet. In Table 1 are reported, for a sequence of meshes, the errors between the numerical and the exact solution together with the computed order of accuracy for linear, quadratic and cubic elements. The scheme constructed has proved to be accurate with the right order and has been also used with quadrilateral elements and on hybrid grids. See ref. [START_REF] Abgrall | Construction of very high order residual distribution schemes for steady inviscid flow problems on hybrid unstructured meshes[END_REF] for further details.

h ǫ L2 (P 1 ) ǫ L2 (P 2 ) ǫ L2 (P 3 ) 1/

III. Discretization of the Euler equations

In this section we consider the extension of the Eq. ( 1) to case of a system of the equations, namely

∇ • f(u) = 0, (12)
where u(x) ∈ R p , with p the number of unknowns and f = (f 1 , . . . , f d ) is the flux function associated to u, with f i=1,d ∈ R p . In the case of the system of the Euler equations for steady flows, written in the conservative form, the vector unknown is the vector of the conservative variables density, momentum and total energy per unit volume u = (ρ, m, E t ) T , with the momentum m = ρu, were u = (u, v) T is the velocity vector. The components of the flux function are given by

f x (u) =      ρu ρu 2 + P ρuv u(E t + P )      , f y (u) =      ρv ρuv ρv 2 + P v(E t + P )     
with P the thermodynamic pressure, related to the other conservative variables by the equations of state. In this work we use only the model of a polytropic ideal gas, for which

P = (γ -1) E t - 1 2 ρ u 2 .
In the case of the system, the first order Rusanov's scheme reads

Φ T σ = Φ T N T DOF + α T (u h -ū),
where ū is the average state ū = σ∈T uσ NDOF and α T is taken as twice the spectral radius of the Jacobian matrix of the fluxes at the degrees of freedom.

The high order scheme is obtained by the same technique described in the sub-section C, however the ratios Φ T σ /Φ T are not defined in the case of system. In this case, resorting the eigenstructure of the system, it is possible to extend the limiting technique to the case of the system. We define first the local average flow direction as

η = ū ū = (η 1 , η 2 ) T ,
where we have supposed to write the equations in two spatial dimensions for simplicity. The Jacobian matrix in the direction of η reads

A η (ū) = A 1 (ū)η 1 + A 2 (ū)η 2 ,
which is diagonalizable in R with eigenvalues λ i=1,4 . Recalling that a generic vector w ∈ R 4 can be decomposed as

w = 4 i=1 w i r i ,
where r i are the right eigenvectors. The limiting technique can be applied ad follows 1. Decompose the first order residual into the characteristic components:

Φ T σ = 4 i=1 Φ T σ,i r i , ∀σ, with σ∈T Φ T σ,i = Φ T i .
2. Limit the characteristic residual as in the scalar case:

β T σ,i = Φ T σ,i Φ T i + σ∈T Φ T σ,i Φ T i +
3. Construct the limited (high) order residual:

ΦT σ = 4 i=1 β T σ,i Φ T σ,i r i
where r i are the right eigenvectors.

The final form of the high order residual in the case of system reads

Φ ⋆,T σ = ΦT σ + T (A η • ∇φ σ ) τ A η • ∇u h dx, ( 13
)
where τ is a scaling matrix that as been taken as

τ = h -1 T N with N = 2 n∈∂T |A| n -1
,

where n are the inward unit normal vectors to the boundaries ∂T and Ān is the Jacobian matrix evaluated at the average state ū on the direction of n.

IV. An iterative solver for the implicit scheme

The discrete counterpart of the system of equations ( 12) is obtained by assembling for each DOF σ the contributions from all the elements T ∋ σ, namely

T ∋σ Φ T σ (u h ) = 0, ∀σ. (14)
The previous equations represent a system of non linear equations in the form

F (u) = 0, (15)
to be solved by the means of an iterative process. Among all kinds of methods for solving a non linear system of equations, the Newton's method is one of the most popular and has a local quadratic convergence. The general form of the Newton's method for solving (15) is

u k+1 = u k -J(u k ) -1 F (u k ), k = 0, 1, 2, . . . ( 16
)
where u 0 is an initial guess of the solution and J(u k ) = ∂F (u k ) ∂u k , the Jacobian of F , is non singular at each iteration. In practice, the Newton iteration ( 16) is implemented by the following two steps

J(u k )∆u k = -F (u k ) u k+1 = u k + ∆u k
Usually the problem ( 14) is replaced by a pseudo-transient one and the steady solution is the limit, for the time variable that tends to infinity, of the pseudo-unsteady problem

T ∋σ Φ T σ (u h ) = 0 -→ |C σ | ∂u h σ ∂t + T ∋σ Φ T σ (u h ) = 0,
with |C σ | the area of the dual cell associated to the DOF σ. The presence of the time derivative enables a better convergence of the Newton's method, overcoming the harsh start-up phase when the solution is far from an optimal initial guess. Furthermore, the Jacobian associated to the modified problem is better conditioned than the Jacobian of the original problem during the start-up phase. The pseudo-transient scheme can be written in the following form

|C| ∂u ∂t = -F (u).
When the Backward Euler formula is used for the discretization of the time derivative, the fully discrete counterpart of the previous system reads

|C| u n+1 -u n ∆t n = -F (u n+1 ), n = 0, 1, 2, . . .
where n is number of the time steps and u 0 is the initial value of the solution. For each time step n a non linear problem should be solved with the Newton's method, namely

|C| ∆t n I + J(u n k ) ∆u n k = -F (u n k ) u n k+1 = u n k + ∆u n k , k = 0, 1, 2, . . . ,
with I the identity matrix. In practice at each time step only one Newton iteration is performed.

The parameter ∆t n is the discrete time step, with ∆t n → ∞ as n → ∞. Note that for ∆t n → ∞ the iteration of the original Newton's method ( 16) is retrieved. The evolution of the time step is controlled by the CFL number which is chosen according to the following law 15

CFL n = CFL n-1 F (u n-1 ) ∞ F (u n ) ∞ , (17)
starting from a low CFL number. The iterative process is stopped when the residual of the equations becomes small enough respect to the initial residual.

A. Jacobian-free Newton/GMRES method

In the previous sub-section we have seen that at each time step the following linear system must be solved

|C| ∆t n I + J(u n ) ∆u n = -F (u n ),
which is recasted for convenience in the following form: A(u n )∆u n = -F (u n ). The matrix A is non symmetric and has dimension N × N with N = N DOF × p, so the number of the non-zero elements can be very high. Krylov methods can be used to solve this class of linear systems, in particular the GMRES [START_REF] Saad | GMRES: a generalized minimum residual algorithm for solving nonsymmetric linear systems[END_REF] is widely used. This method has the property of minimizing the L 2 -norm of the residual over all vectors in the Krylov subspace. The GMRES method computes a new search vector every iteration. The vector is added to the Krylov subspace to progressively improve the solution. However, more search directions incur higher memory and computational costs. For large problems, this limits the maximum number of iterations that can be used. The restarted version of the algorithm can be used, where the algorithm is restarted from the most recent solution.

To accelerate the convergence of the iterative linear solver, preconditioning of the matrix A is used. This consists in solving a modified linear system AP -1 P ∆u = -F , with P a preconditioning matrix. When the right preconditioning is used, one first solves

AP -1 w = -F ,
for w, and the solves ∆u = P -1 w, for ∆u. Only P -1 is required. The right preconditioned GMRES pseudo-code is reported for clarity in the Algorithm 1.

Complete solving of the linear system is unnecessary for the convergence of the scheme, usually inexact Newton's method is used to reduce the computational effort and avoid over-solving of the system. [START_REF] Eisenstat | Choosing the forcing terms in an inexact Newton method[END_REF] The linear system is solved until

F (u n ) + A(u n )∆u n ≤ η n F (u n ) ( 18
)
with a tolerance η n < 1.

The construction of the matrix A requires to compute the Jacobian of F . In order to obtain the quadratic convergence of Newton's method, the linearization of the residuals must be exact. Unfortunately, explicit formation of the Jacobian matrix resulting from the linearization of the high order residuals is extremely expensive, if not impossible. The Jacobian of the low order residual is generally used in the construction of the matrix A, but the quadratic convergence of the Newton's method is lost due to the inconsistency between the right hand side of the linear system, constructed with the high order residuals, and the matrix, constructed with the low order residuals.

Looking at the Algorithm 1, the matrix A is required only in the form of a matrix-vector products as

Av = |C| ∆t n v + J(u n )v,
the matrix-vector products which involve the Jacobian can be approximated by using a finite difference approximation of the directional derivative 18, 19

J(u n )v ≃ F (u n + ǫv) -F (u n ) ǫ , (19)
with ǫ a small step size chosen as 20

ǫ = 1 + u L2 v L2 ǫ rel Algorithm 1 Right preconditioned GMRES. Solve Ax = b 1: Choose x 0 2: r 0 = b -Ax 0 3: β = r 0 2 4: v 1 = r 0 /β 5: Define H m = {h i,j } 1≤i≤m+1,1≤j≤m
6: H m = 0 7: for j = 1 to m do 8:

w 1 = AP -1 v j 9:
for i = 1 to j do 10:

h i,j = (w j , v i )

11:

w j = w j -h ij v i 12:
h j+1,j = w j 2

13:

if h j+1,j = 0 then v j+1 = w j /h j+1,j

18:

end for 19: end for 20: y m = argmin

y βe 1 -H m y 2 21: V m = [v 1 , . . . , v m ] 22: x m = x 0 + P -1 V m y m 23: if βe 1 -H m y 2 < η n b 2 then 24:
J(u n )P -1 v ≃ F (u n + ǫP -1 v) -F (u n ) ǫ .
Since there is no need to compute explicitly the Jacobian, this approach is called Jacobian-free. However, a rough approximation of the Jacobian is always computed at each step, this matrix is used as a preconditioner in the GMRES algorithm. In other works [START_REF] Luo | A Fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF][START_REF] Marques | Comparison of matrix-free acceleration techniques in compressible Navier-Stokes calculations[END_REF] the Jacobian-free technique is implemented in a full matrix-free version, this means that even the preconditioning matrix is computed by the means of the Eq. ( 19).

The Jacobian-free approach allows quadratic convergence of Newtons method because the matrix of the linear system is a complete linearization of the residual vector. The price to pay for using this technique is an increment of the computational effort, because at each time step it is necessary to compute several times the residual F (u) on the whole domain. This is largely compensated by a drastic diminution of the iteration number, as it can be seen in the experiments that are now reported.

B. Verification of the effectiveness of the Jacobian-free technique

We tested the implicit method with the Jacobian-free technique on two dimensional problems and we compared the convergence of the resulting scheme against the implicit scheme that makes use of the Jacobian computed explicitly. The Jacobian is computed considering a first order scheme, not limited and not stabilized. The solution tolerance for the GMRES is set to 0.3 (i.e., η n = 0.3 in the Eq. ( 18)), with maximum 300 iterations. The CFL number is allowed to variate according to the law (17), starting from an initial value of 10 to a maximum of 10 9 .

The first problem we consider is a subsonic flow in channel with a 10% thick bump on the bottom, the inlet Mach number is 0.5. The mesh used consists of 495 nodes (888 triangles). In Fig. 1 are showed the Mach isolines of the solution obtained with the third order scheme. In Fig. 2 the convergence curves are reported in terms of number of iterations and CPU time (in seconds) for the Jacobian-free technique with different type of preconditioners. The residual is computed as the norm L 2 of the continuity equation and is normalized respect to the residual at the first iteration. In all the computations 40 search directions has been used in the GMRES algorithm. It can be seen that the GMRES algorithm without preconditioning does not perform well even for a simple problem, as that considered here. The ILU(0) preconditioner allows the residual to reach the zero machine as well as the Jacobi and the LU-SGS [START_REF] Yoon | Lower-upper implicit schemes with multiple grids for the Euler equations[END_REF] preconditioners, but the last two types of preconditioning allow the implicit scheme to converge faster. In the Jacobian-free technique the preconditioner matrix is computed using the Jacobian of the first order scheme. The same matrix used directly as Jacobian in the implicit scheme does not allow the implicit scheme to converge to the zero machine, whatever the preconditioning strategy is. Different preconditioners has been used for the Jacobian-free technique, the convergence history for the implicit scheme that uses the Jacobian of the first order scheme is also reported.

Iterations

In Fig. 3 is evaluated the influence of the number of the search directions in the GMRES algorithm in combination with the LU-SGS preconditioner. As expected the use of the implicit method with the Jacobian of the first order scheme does not make the scheme converge, whatever the number of search directions is.

The second test case considered is a transonic flow, with Mach number M ∞ = 0.8, over NACA 0012 airfoil at 1.25 • of incidence. The mesh used for the simulation consists of 6 199 nodes (12 190 triangles), a particular of the mesh is depicted in Fig. 4 on the left side, while in the righ side are reported the isolines of the Mach number for the solution computed with the third order scheme. In Fig. 5 Different preconditioners has been used for the Jacobian-free technique, the convergence history for the implicit scheme that uses the Jacobian of the first order scheme is also reported.

of the implicit scheme with the Jacobian-free technique is compared against that of the implicit scheme constructed with the Jacobian of the first order method. It is clear that the Jacobian-free technique allows the solution to converge at the steady state.

V. Extension to viscous term

In this section we extend the RD scheme, developed previously for purely advection problem, to the case of advection-diffusion problems. The scalar model equation for the advection-diffusion problems reads where ν > 0 is the viscosity, generally function of u and ∇u. The portion of the boundary ∂Ω -is the inflow boundary where the weak boundary conditions are specified and ∂Ω 0 represents the part of the boundary where strong boundary conditions are imposed. Within the RD framework, the discretization of the viscous terms has been traditionally obtained by coupling, for P 1 elements, the scheme for the advection equations with the Galerkin approximation of the viscous terms. [START_REF] Paillere | Multidimensional upwind residual distribution schemes for the convectiondiffusion equation[END_REF] This technique has been explored further recently, [START_REF] Villedieu | Third order residual distribution schemes for the NavierStokes equations[END_REF] in particular most of the work concerns how to obtain an uniform order of accuracy. A different approach [START_REF] Nishikawa | A first-order system approach for diffusion equation. I: Second-order residual-distribution schemes[END_REF][START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] is based on the idea that the (steady) diffusion equation Eq. ( 20) is equivalent to a hyperbolic relaxation system at the steady state. This approach give a new point of view for the discretization of the diffusion terms.

∇• f (u) = ∇• (ν∇u), ∀ x ∈ Ω ⊂ R d , d = 2, 3 u = g(s), s ∈ ∂Ω - u = g 0 (s), s ∈ ∂Ω 0 (20)
In this section we first explain why a scheme resulting from the coupling of a RD scheme with a Galerkin discretization of the viscous term is still a residual method, and then we show that this approach works only for P 1 elements. We propose an alternative way to discretize the advection-diffusion equation in the framework of the residual method.

A. Approximation with P 1 elements

Let us consider the RD scheme for the advection problem in the form

T ∋σ Φ T σ (u h ) = 0, ∀σ, with Φ T σ = β σ i Φ T .
Using the P 1 shape functions, ψ σ , we can write Φ T σ , in a Petrov-Galerkin manner, as

Φ T σ = T ψ σ ∇• f (u h ) + T β T σ - 1 3 ∇• f (u h ) = T ω T σ ∇• f (u h ),
where the differentials have been omitted in the integrals for sake of brevity, as hereafter. Unfortunately ω T σ is not continuous across the edges and cannot be used to approximate the problem.

The same scheme can be written in a different way. Denote b T the bubble function on the element T which is 0 on ∂T and 1 in the gravity center of T , we can write

Φ T σ = β T σ Φ T = T ψ σ ∇• f (u h ) + γ T σ T b T ∇• f (u h ) with γ T σ T b T = β T σ - 1 3 |T |,
it is possible to define a new continuous function

ω σ | T = ψ σ + γ T σ b T .
Denoting by W h = span(ω σ ) and V h = span(ψ σ ), the problem for the advection-diffusion equation becomes to find

u h ∈ V h such that for all ω ∈ W h T ω ∇• f -ν∆u h = 0,
where for simplicity the viscosity ν has be assumed to be constant. The fully discrete form of the previous equation reads

T ∋σ T ω σ ∇• f (u h ) + T ∋σ T ∇ω σ ∇u h = 0.
The term on the left hand side gives back β T σ Φ T , while the term on the right hand side reads

T ∇ω σ ∇u h = T ∇ψ σ ∇u h + γ T σ ∇u h • T ∇b T ,
since ∇u h is constant for linear elements. Furthermore by the Green formula T ∇b T = ∂T b T n = 0, so the variational formulations is: find u h such that for any σ

T ∋σ β T σ Φ T + ν T ∇ψ σ • ∇u = 0,
showing that the scheme corresponds to the RD discretization of the advection term plus the Galerkin approximation for the diffusion. This explain why the method can be see as a RD scheme, although this is not clear from the initial formulation. However, the coercivity of the scheme cannot be proved and it is necessary to introduce a blending parameter between the RD and the Galerkin scheme in order to get and uniform order of accuracy, [START_REF] Villedieu | Third order residual distribution schemes for the NavierStokes equations[END_REF] as it is also standard for the SUPG scheme at low Reynolds.

B. Extension to higher degree elements

In the previous sub-section we showed how the scheme resulting from the blending between the RD and the Galerkin scheme can be see still as a RD scheme in the case of P 1 . We try now to extend the same consideration to the case of higher order elements.

Let us define a function γ T σ ∈ H 1 (T ) such that 1. when we use a P k Lagrange interpolant

T (ψ σ + γ T σ )∇• f (u h ) = β T σ T ∇• f (u h ) rearranged as T γ T σ ∇• f (u h ) = β T σ T ∇• f (u h ) - T ψ σ ∇• f (u h ), (21)
2. enable to construct the H 1 basis functions

γ T σ | ∂T = 0 (22)
3. have no role on the viscous terms

T γ T i ∆u = T ∇• γ T σ ∇u - T ∇γ T σ • ∇u = 0 that means T ∇γ T σ • ∇u = 0, (23)
because of the condition (22).

The conditions ( 21)-( 23) are affine conditions of the type

ℓ p (γ T σ ) = a p ,
where the linear functions ℓ p are defined such that

ℓ 1 P (ω) = T ω∇• f (u h ) and a p = β T i T ∇• f (u h ) - T ψ σ ∇• f (u h )
form the condition (21), and

ℓ 2 P (ω) = T ∇ω∇ψ σ and a p = 0
In general, there is no solution to this problem as one can easily see considering a simple one dimensional case with quadratic elements.

C. Variational method based on the reconstruction of the gradient

Let us consider the variational formulation of the residual distribution scheme for the advection equation

σ v σ T ∋σ = T σ∈T β T σ v h σ ∂T f (u h ) • n, with v h ∈ V h ⊂ H 1 .
We introduce now the space W h of the functions that are piecewise constant on the element T and the mapping

v h ∈ V h → π h β (v h ) ∈ W h , with π h β (v h ) = σ∈T β T σ v h σ .
We can reformulate the RD scheme as finding u h ∈ V h such as for any

v h ∈ V h a(u h , v h ; u h ) = 0, with a(u h , v h ; w h ) = T ∂T π h β (w h )(v h )f (u h ) • n - T ∇π h β (w h )(v h ) • f (u h ),
where we have indicated the possible dependence of β T σ on w h by the mechanism of the limitation. We note that the exact solution u of the problem, if it is smooth enough, also satisfies the residual condition a(u, v h ; w h ) = 0 for any v h , w h ∈ V h because u is smooth enough, we can also write

T T π h β (w h )(v h )∇•f (u) = 0.
The last relation can be used to obtain the formulation for the viscous problem. Assuming u smooth enough, we can write

T T π h β (w h )(v h ) [∇•(f (u) -ν∇u)] = 0,
from which we obtain

T ∂T π h β (w h )(v h )f (u) • n + ∂T π h β (w h )(v h )ν∇u • n = 0.
Using the average operator, we can rewrite the previous relation as

T ∂T π h β (w h )(v h )f (u) • n + ∂T π h β (w h )(v h ){ν∇u • n} = 0,
so the variational formulation of the RD scheme for the advection-diffusion reads

T ∂T π h β (w h )(v h )f (u) • n + ∂T π h β (w h )(v h ){ν∇u • n} + h T T λ∇v h -∇ • (ν∇v h ) λ∇u h -∇ • (ν∇u h ) = 0.
The residual distribution process is the same as for the advection problem, where now the total residual on the element contains both the advection and diffusion term, namely

Φ T = ∂T f (u h ) -{ν∇u} • n
By evaluating the nodal gradient by a Green-Gauss-like formula

{∇u h } σ = T ∋σ |T |∇u h T ∋σ |T |
, it is possible to interpolate the gradient on the elements by the same Lagrangian test function as

∇u h = σ∈T {∇u} σ ψ σ
and the term ν∇u h can be approximated with the right order as

ν∇u h ≃ σ∈T {ν∇u} σ ψ σ .

VI. Numerical results for the advection-diffusion problem

In this section we report some numerical results for the discretization of the advection-diffusion problem. First we consider a scalar problem, for which it is reported a convergence study of the scheme. The second kind of results concerns the discretization of the Navier-Stokes equations.

A. Scalar equation

In order to check the accuracy of the scheme, we consider the following scalar problem ∂u ∂y -ǫ ∂ 2 u ∂x 2 = 0, on the square [0, 1] × [0, 1] with the following boundary conditions u(x, 0) = sin(πx) on y = 0 u(x, 0) = ϕ(x, y) on x = 0, x = 1, with ϕ(x, y) = e -ǫy sin(πx), that is also the exact solution of the scalar equation. We performed calculations on a series of four meshes with 10, 20, 40 and 80 nodes on each side. The results are showed in the table A. Results confirm the third order accuracy of the scheme, however the formal accuracy degrades when the Reynolds becomes smaller.

We consider now the Smith-Hutton test case. The problem consists in the solution of an advectiondiffusion problem on [-1, 1] × [0, 1] with the advection vector defined as

λ = (2y(1 -x 2 ), -2x(1 -y 2 )) T , the inflow profile is u(x, 0) = 1 + tanh(α(2x + 1)), ∆x (Log)
     u(1, y) = 1 -tanh(α) u(x, 1) = 1 -tanh(α) u(x, -1) = 1 -tanh(α)
This is a standard case where a sharp transition between two constant states occurs. In most of the reported calculations, if not all, the parameter α that drives the sharpness of the transition is set to α = 10. Here we take α = 100 because we also want to test the properties of the scheme with respect to the maximum principle. No exact solution in closed form is knowns, except in the case ν = 0 of course.

In Fig. 6 are reported the contour of the solution for two different values of the viscosity coefficient. This test case is useful since it allows to verify the non-oscillatory character of scheme near zones with strong gradients. As can be seen from the Fig. 7, which display the profile on the outflow boundary, there is no oscillation on the solution.

B. Navier-Stokes equations

We consider, now, the two dimensional Navier-Stokes equations written in conservation form

∂u ∂t + ∇•f E -∇•f V = 0, 18 of 24
American Institute of Aeronautics and Astronautics where u and f E are the conservative variables and the advection flux, respectively, as defined in the Euler equations. The components of the viscous flux f V are given by and where µ is the dynamic viscosity, C p is the specific heat at constant pressure, P r is the Prandtl number and T is the temperature. The system of the Navier-Stokes equations can be easily rephrased as follows

f V x = µ      0 τ xx τ xy uτ xx + vτ xy - cp Pr ∂T ∂x      , f V y = µ      0 τ xy τ yy uτ xy + vτ yy - cp Pr ∂T ∂y      , with τ xx = µ 4 3 ∂u ∂x - 2
∂u ∂t + ∇•f E + ∇•(K∇u) = 0
with K∇u = K i,j ∂u ∂x j .

In the first test case we consider the laminar flow on an adiabatic flat plate characterized by a free stream Mach number M = 0.3 and by a Reynolds number based on the free stream conditions and on the plate length Re = 1 000. The length of the plate is set to be L = 1.0. The range of the computational domain in the x-direction is [-1, 1] and in the y-direction is [0, 1.5] such that the size of the domain in the y-direction is almost 10 times the boundary layer thickness at the end of the plate. At the inlet, the free stream condition is imposed. At the top and exit boundaries, the static pressure is imposed. At the bottom boundary, the adiabatic wall boundary condition is imposed on the plate and the symmetry plane boundary condition is imposed on the other portion of the boundary.

The mesh used for the simulation is depicted in the Fig. 8, it consist of 779 nodes (1 440) triangles, with 21 points on the plate and 19 points on the y-direction of which 6 points are in the boundary layer. In Fig. 9 are reported the velocity profiles, with comparison to the Blasius solution, computed with the second and third order scheme. In the same figure is reported also the simulation with the scheme at second order on a finer grid which has the same number of DOF as the third order scheme on the coarser mesh. The non-dimensional wall distance is defined as y * = ρ∞U∞ µx . The computed skin friction profiles are reported in Fig. 10. The agreement between the computed and the exact solution is very good for the third order scheme. It is worth noting that the third order scheme performs better that the second order scheme with the same number of DOF.

The second test case we consider is a subsonic, laminar flow over a Naca airfoil at zero incidence. The free stream Mach number is M = 0.5 and the Reynolds number, based on the airfoil cord, is Re = 500. The computational grid is displayed in Fig. 11 and consists of 11 959 points (21 591) triangles. In Fig. 12 are displayed the Mach contours for the second and third order simulation and in Fig. 13 are reported the computed Cp profiles.

VII. Conclusion

We have presented the general framework of the RD scheme for the high order discretization of the conservation laws. First we described the accuracy proprieties of the scheme in the case of the scalar advection problem. Then we showed how to extend the scheme to the case of the Euler equations and we showed how to speed up the convergence of the method by the Jacobian-free technique. We extended the scheme to the case of the advection-diffusion problem, and we verified the accuracy of the scheme on scalar problems and we used the scheme to obtain the solution of the Navier-Stokes equations on the flat plane and the NACA 0012 airfoil.

x

  Vector of the spatial dimensions N T DOF Number of the degree of freedom of the element Φ T Total residual on the cell Φ

  with ∂Ω boundary of Ω and n the outward normal vector to the boundary of the domain. The function g is known and represents the weak Dirichlet boundary conditions of the problem on the boundary ∂Ω -.

  end if with ǫ rel = 10 -10 . The right preconditioned version of Eq. (19) reads

Figure 1 .

 1 Figure1. Mach isolines (∆M = 0.01) for the bump simulation using a third order scheme.

Figure 2 .

 2 Figure 2. Convergence curves as function of the number of iterations (left) and the CPU time (right) for the bump problem. Different preconditioners has been used for the Jacobian-free technique, the convergence history for the implicit scheme that uses the Jacobian of the first order scheme is also reported.

Figure 3 .

 3 Figure 3. Convergence curves as function of the number of iterations (left) and the CPU time (right) for the bump problem. Different preconditioners has been used for the Jacobian-free technique, the convergence history for the implicit scheme that uses the Jacobian of the first order scheme is also reported.

Figure 4 .

 4 Figure 4. Transoni flow over a NACA airfoil at Mach 0.8 and incidence 1.25 • . Particular of the mesh (left) and Mach isolines (right) with ∆M = 0.05.

Figure 5 .

 5 Figure 5. Convergence curves as function of the number of iterations (left) and the CPU time (right) for the simulation around the NACA airfoil.

Figure 6 .Figure 7 .

 67 Figure 6. Smith-Hutton problem: contour of the solution (left) ν = 0.01, (right) ν = 0.0001.

Figure 8 .

 8 Figure 8. Mesh used for the simulation of the flat plate boundary layer.

Figure 9 .Figure 10 .

 910 Figure 9. Non-dimensional axial velocity for the plate plane case.

Figure 11 .Figure 12 .Figure 13 .

 111213 Figure 11. Computational grid for subsonic viscous flow over the NACA0012 airfoil

Table 1 .

 1 Computer errors and order of accuracy for the advection problem with linear, quadratic and cubic elements.

	25	0.50493E-02 0.32612E-04 0.12071E-05
	1/50	0.14684E-02 0.48741E-05 0.90642E-07
	1/75	0.74684E-03 0.13334E-05 0.16245E-07
	1/100 0.41019E-03 0.66019E-06 0.53860E-08
		O = 1.790	O = 2.848	O = 3.920

Table 2 .

 2 Convergence results for the advection-diffusion scalar problem and the boudary condtions are

		L ∞ error (Log)	Order	L 2 error (Log)	Order
	-0.532115666963180 -2.41783092916874		-2.42918055471673
	-0.846872634669396 -3.22731516724477	2.57	-3.15553327436338	2.30
	-1.08957273264021	-4.05793545206366	3.42	-3.87000533630969	2.94
	-1.36540918681519	-4.90199016882381	3.06	-4.60138273684662	2.65
		ǫ = 0
	∆x (Log)	L ∞ error (Log)	Order	L 2 error (Log)	Order
	-0.532115666963180 -2.42235466229356		-2.43644370369152
	-0.846872634669396 -3.24877046688954	2.62	-3.21509140129168	2.47
	-1.08957273264021	-4.09492244395854	3.48	-3.95823335106917	3.06
	-1.36540918681519	-4.99047469215026	3.24	-4.85559507238436	3.25
		ǫ = 0.0001
	∆x (Log)	L ∞ error (Log)	Order	L 2 error (Log)	Order
	-0.532115666963180 -2.45230965825349		-2.52191658082643
	-0.846872634669396 -3.29453851242374	2.67	-3.26021775685192	2.34
	-1.08957273264021	-4.01681756317218	2.97	-3.74468087319104	1.99
	-1.36540918681519	-4.71151297471185	2.51	-4.48933815669847	2.7
		ǫ = 0.001
	∆x (Log)	L ∞ error(Log)	Order	L 2 error (Log)	Order
	-0.532115666963180 -2.12079249189368		-2.07369114240901
	-0.846872634669396 -2.56866661478255	1.42	-2.55012986275973	1.51
	-1.08957273264021	-3.19486137685157	2.58	-3.19339332748624	2.65
	-1.36540918681519	-4.13416777580946	3.40	-3.82076591805969	2.27
		ǫ = 0.01

of 24 American Institute of Aeronautics and Astronautics

of 24 American Institute of Aeronautics and Astronautics

American Institute of Aeronautics and Astronautics

Acknowledgments

R. Abgrall has been financed in part by the EU ERC Advanced grant "ADDECCO" #226616. D de Santis has been financed by the EU STREP IDIHOM #265780.