R Abgrall

D De Santis

M Ricchiuto

High order preserving residual distribution schemes for advection-diffusion scalar problems on arbitrary grids

Keywords: Higher order schemes, Residual distribution, Viscous term, Advection-diffusion problems, Gradient recovery Schémas d'ordre élevé, schémas distribuant le résidu, problèmes d'advection diffusion, reconstruction du gradient

This paper deals with the construction of a class of high order accurate Residual Distribution schemes for advection-diffusion problems using conformal meshes. The problems we consider range from pure difusion to pure advection. The approximation of the solution is obtained using standard Lagrangian finite elements and the total residual of the problem is constructed taking into account both the advective and the diffusive terms in order to discretize with the same scheme both parts of the governing equation. To cope with the fact that the normal component of the gradients of the numerical solution is discontinuous across the faces of the elements, the gradient of the numerical solution is recovered at each degree of freedom of the grid and then interpolated with the same shape functions used for the solution. Linear and non-linear schemes are constructed and their accuracy is tested with the discretization of advection-diffusion and anisotropic diffusion problems.

Schémas distribuant le résidu d'ordre élevé pour les problèmes de convection diffusion sur des maillages arbitraires 1 Introduction

In the last years different high order schemes have been developed to obtain an higher order (more than two) discretization of the Navier-Stokes equations. One of the most attractive scheme seems to be the discontinuous Galerkin (DG) scheme [START_REF] Cockburn | Discontinuous Galerkin methods: theory, computation and application[END_REF]. Residual Distribution (RD) schemes [START_REF] Struijs | Fluctuation splitting schemes for the 2d Euler equations[END_REF][START_REF] Abgrall | Toward the ultimate conservative scheme: Following the quest[END_REF][START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF] represent a very interesting alternative to DG schemes. 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 (DOF) with the cell polynomial degree. In RD schemes the formulation remains local, as in DG, but the number of DOFs growths less quickly because the solution is assumed to be continuous. Another difference between RD and DG schemes is that, to date at least, the non oscillatory properties of the RD scheme in the case of discontinuous solutions are probably better understood that for their DG counterpart.

RD schemes have been developed mainly for advection problems due to possibility to construct multidimensional upwind schemes which guarantees a small discretization error compared to the standard Finite Volume schemes, but the discretization of advection-diffusion problems with the RD schemes is still an open problem. One of the main issue concerns the possibility to take into account within the same scheme advective terms, by the means of upwind mechanism, and diffusive phenomena, which on the other hand have an isotropic behavior. To address this problem mixed upwind/central schemes have been developed, in which RD methods for the advection terms are combined with central schemes, usually based on the Galerkin discretization of the diffusion terms. For such type of schemes a proper blending between the RD and the Galerkin schemes must be constructed otherwise the accuracy of the resulting schemes is spoiled [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF]. The approach used in this work is based instead on the construction of a RD method in which the advection and the diffusion are handled within the same scheme. Unfortunately this introduces a new complication because, for polynomial piecewise approximation of the solution, the normal component of the gradient of the numerical solution is discontinuous on the face of two adjacent elements. This would require the introduction of a numerical flux for the viscous term.

Instead of actually considering a numerical flux along the faces of the elements, as happens in the DG or Finite Volume schemes, the approach adopted in this work consists in recovering a unique set of values for the gradient of the numerical solutions at each degree of freedom (DOF) of the grid. Then, these values are interpolated with the same continuous functions used to interpolate the solution. It is evident that gradients have to be recovered with higher order of accuracy to construct a high order scheme: the crucial point is the strategy used to recover the gradients at the DOFs. The problem of the gradient recovery is addressed in the paper together with the construction of accurate and robust, linear and non-linear RD schemes. Note that this issue has already been considered by D. Caraeni [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF], the solution we propose here is different and, in our opinion, easier to implement with a more compact stencil The structure of the paper is as follows. In Section 2, starting from the advection problem, the basic ideas of the RD methods are recalled. The issues related to the discretization of advectiondiffusion problems in the RD framework are discussed in Section 3, while in Section 4 is described in detail the construction of a class of linear and non-linear RD schemes. In Section 5 are discussed and compared different strategies for the gradient recovery. In Section 6 the proposed numerical schemes are extensively tested on linear and non-linear scalar advection-diffusion problems and an anisotropic diffusion problem is also considered. Finally, some concluding remarks are given in the last section.

Basics of the residual distribution method

In this section, the main idea of the RD method for scalar hyperbolic problems is briefly recalled and the fundamental properties of conservation, consistency and accuracy of the numerical scheme are also reported. Furthermore, the notation used through the paper is introduced.

Consider the steady conservation law for the scalar quantity u

∇•f (u) = 0, (1a)
where f (u) ∈ R d is a given flux function of the unknown u(x) ∈ R, x ∈ Ω ⊂ R d , with d the number of the spatial dimensions (here d = 2 or d = 3). The Eq. (1a) must be supplemented with the proper boundary conditions on the inflow portion of the boundary ∂Ω

u| ∂Ω -= g(s), s ∈ ∂Ω -, (1b)
where the function g is known. It represents the boundary condition of the problem on the inflow boundary and ∂Ω -= {x ∈ ∂Ω | a • n < 0}, with n the outward normal vector to the boundary of the domain and a the advection velocity defined by a = df du .

The domain Ω is discretized with N e non-overlapping elements with characteristic length h, the set of all the elements 1 is denoted by E h , the list of the DOFs is denoted by Σ h , the set of all the boundary faces is denoted by F h 2 , and the total number of DOFs is N dof . The solution is approximated on each element by k-th order polynomials which are assumed to be continuous within the elements and accross the faces of the elements. If the standard Lagrangian shape functions are used, the approximated solution u h can be written as

u h (x) = i∈Σ h ψ i (x) u i , x ∈ Ω,
with u i the numerical solution at the generic DOF i, and ψ i the Lagrange basis function at the DOF i.

The approximated solution, in general, will not satisfies the governing equation, which means that the integral of the Eq. (1a) calculated on each elements e will be not null, but will give rise to a residual on each element, namely

Φ e u h = Ωe ∇•f u h dΩ = ∂Ωe f u h • n d∂Ω.
The integral quantity Φ e u h is called the total residual of the element e. Similarly, for any element f on the boundary ∂Ω of Ω, one can define a total residual

Φ f u h = f f u h • n -F u h , g d∂Ω
where F is a numerical flux consistant with (1b).

1 assumed to be polyhedrals 2 a boundary face f is the intersection of an element e and ∂Ω. The mesh is assumed to so that a boundary face f = e ∩ ∂Ω is also a face of e.

Inria

In order to handle only nodal values, the total residual is first distributed, in some way, to the DOFs of the element as follows

Φ e i = β e i u h Φ e u h , ∀i ∈ Σ e h ,
where Σ e h is the list of the DOFs of the element e and β e i are the distribution coefficients, which can be in general function of u h . A similar relation is written on the boundary residuals

Φ f i = β f i u h Φ f u h ∀i ∈ Σ f h ,
where Σ f h is the list of DOFs of the boundary face f . It is easy to see that the following conservation constraints must be satisfied for any element e and any boundary face f [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF]

i∈Σ e h Φ e i = Φ e , ∀e ∈ E h i∈Σ f h Φ f i = Φ f ∀f ∈ F h
in addition to the standard assumptions of the Lax-Wendroff theorem in order to garanty that the limit solution, if it exists, is a weak sokution of (1a). To obtain an equation for each nodal value of the numerical solution the following relations are written for each DOF

e∈E h,i Φ e i u h + f ∈F h,i Φ f i u h = 0, ∀i ∈ Σ h ,
where E h,i (resp. F h,i ) is the set of the elements (resp. faces) which share the DOF i. The previous relations define a set of non-linear equations that must be solved for nodal values of the solution [u i ] i=1,..., N dof . In practice the solution with an RD method is obtained by the means of an iterative method, which in the simplest form reads

u n+1 i -u n i ∆t n i + e∈E h,i Φ e i u h = 0, ∀i ∈ Σ h , (2)
with ∆t n i a scaled pseudo-time step. The change of the nodal values of the solution during the iterative process is driven by the non-zero total residuals on the elements; for n → ∞ the total residual on each element vanishes and the steady state solution is obtained.

Consistency and accuracy

The fundamental properties of consistency and accuracy for RD schemes have been analyzed in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF] and are briefly reported here for sake of completeness.

Assuming that a sequence u h is bounded in L ∞ when h → 0 and if exist w, such that u h → w in L 2 when h → 0, then w is a weak solution of (1). In the proof, the continuity of the interpolant across the faces is assumed, although this constrain may be alleviated and RD schemes with discontinuous elements can be constructed [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].

To analyze the accuracy of the RD the following truncation error is introduced, for any piecewise polynomial function ϕ,

E u h , ϕ = i∈Σ h ϕ(x i ) e∈E h,i Φ e i + f ∈F h,i Φ f i . (3)
If the solution u is smooth enough and the residuals, computed with the numerical solution u h , are such that

Φ e i = O h k+d , ∀e ∈ E h and i ∈ Σ h Φ f i = O h k+d-1 , ∀f ∈ F h and i ∈ Σ h (4)
and if the approximation f u h is accurate with the order k+1, then the truncation error satisfies the following relation

|E u h , ϕ | ≤ C(ϕ, f , u)h k+1 ,
with C a constant which depends only on ϕ, f , and u. It can be shown, under the previous hypothesis, that Φ e = O h k+d and if exists a constant β E i , such that Φ e i = β e i Φ e , then the condition ( 4) is satisfied provided that β e i is uniformly bounded. Such a condition is historically called linearity preserving.

To determine the conditions that must be satisfied by the numerical scheme in order to have non-oscillatory solutions the distributed residual on a generic element e is re-written in the following form

Φ e i = j∈Σ e h j =i c e ij (u i -u j ),
with the coefficients c e ij that in general depend on the solution. The same would hold for the face residuals. By applying the RD scheme [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] with the previous definition of the residual one obtains

u n+1 i -u n i ∆t n i = e∈E h,i j∈Σ e h j =i c e ij u n i -u n j = 0 ∀i ∈ Σ h . (5)
If the scheme satisfies the following positivity conditions

e∈E h,i j∈Σ e h j =i c e ij ≥ 0 and 1 -∆t n i j∈Σ e h j =i c e ij ≥ 0, ∀i ∈ Σ h ,
then the solution verifies the following discrete maximum principle

min e∈E h,i min j∈Σ e h u 0 j ≤ u n i ≤ max e∈E h,i max j∈Σ e h u 0 j , ∀i ∈ Σ h .
A scheme which satisfy the maximum principle is said to be positive. If all the coefficients c ij are independent of the numerical solution a scheme of the form ( 5) is said linear. It is well know from the Godunov's theorem [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for unsteady problem[END_REF] that a linear scheme of the form (5) cannot be simultaneously positive and linearity preserving, which means that a positive and high order scheme must be non-linear.

Extension to the diffusion terms

When in the governing equation (1a) diffusive phenomena are considered together with the advective terms, the following advection-diffusion equation is obtained

∇•f (u) = ∇• ν∇u on Ω ⊂ R d , d = 2, 3 (6)
where ν > 0 is the viscosity, generally function of u. The relative importance of the advection and the diffusion is described by the non-dimensional parameter, Peclet number, Pe = a h/ν.

Inria

In the advection and diffusion limits Pe → ∞, Pe → 0, respectively, while Pe ∼ 1 when advection and diffusion are equally important.

To extend RD methods to advection-diffusion problems, different strategies have been considered to compute and to distribute the residual associated with the diffusion terms. On a first attempt, based on the physical intuition that the diffusion has an isotropic behavior in the space, RD schemes for the advection terms were coupled with the Galerkin discretization of the diffusion terms [START_REF] Van Der Weide | A parallel, implicit, multi-dimensional upwind, residual distribution method for the Navier-Stokes equations on unstructured grids[END_REF][START_REF] Paillere | Multidimensional upwind residual distribution schemes for the Euler and Navier-Stokes equations on unstructured grids[END_REF], but a truncation error analysis revealed that this simple approach results in a first order accurate scheme when advection and diffusion have the same order of magnitude [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF]. A different approach, which was developed for two-dimensional schemes on triangular grids, considered a hybridization of the RD method with a Petrov-Galerkin scheme by the means of a scaling parameter, function of the Peclet number [START_REF] Ricchiuto | On uniformly high-order accurate residual distribution schemes for advection-diffusion[END_REF].

A key aspect that emerges from the work of Nishikawa and Roe [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF] is that a RD scheme with an uniform order of accuracy in all the range of the Peclet numbers should not consider two different distribution schemes for the advection and diffusion terms, but only one distribution process has to be performed for the residual of the whole equation, namely

Φ e = Ωe ∇•f u h -∇• ν∇u h dΩ.
To put the previous expression in term of a boundary integral, one has to cope with the fact that the normal component of the gradient of the numerical solution, ∇u h •n, is in general discontinuous on the faces of the elements and this violates the continuous approximation hypothesis of the numerical scheme. Suppose, now, that an unique value of the gradient is available at each DOF, the gradients can be interpolated with the same shape functions used for the solution and the total residual on the element can be written as follows

Φ e = ∂Ωe f u h -ν ∇u h • n d∂Ω, (7)
where ∇u h is the interpolated gradient of the numerical solution, which is now continuous on the faces of the elements.

Once the total residual is evaluated, it can be distributed to the DOFs of the elements by the distribution coefficients, β e i . This strategy has been adopted in [START_REF] Paillère | Multidimensional upwind residual distribution schemes for the convection-diffusion equation[END_REF] to construct a second order RD scheme for advection-diffusion problems on triangular grids an has been extended to the third order in [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF]. In both works the distribution process is done with purely advective distribution coefficients, which is not appropriate in the diffusion limit. A more general scheme consists in using distribution coefficients which are function of the local Peclet number in order to recover an isotropic scheme in the diffusion limit and an upwind scheme in the advection limit [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF][START_REF] Chou | High order residual distribution conservative finite difference WENO schemes for convection-diffusion steady state problems on non-smooth meshes[END_REF]. Another attempt in that direction is given by [START_REF] Abgrall | Numerical approximation of parabolic problems by means of residual distribution schemes[END_REF], the scheme give satisfactory results except in the region P e ≈ 1, which is typical of a boundary layer. Hence the present contribution can be viewed as an improvement over the previous references.

The key idea of the Eq. ( 7) is the reconstruction of the gradient of the numerical solution at each DOF of the grid and is one of the issue analyzed in this work. Indeed, numerical experiments show that in order to obtain an high order accurate solutions, the gradients must be recovered with the same order of the solution.

An alternative approach has been proposed by Nishikawa for diffusion problems [START_REF] Nishikawa | A first-order system approach for diffusion equation. I: Second-order residualdistribution scheme[END_REF] and advection-diffusion problems [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF], it consists in reinterpreting the advection-diffusion scalar equation as an equivalent hyperbolic first order system, in this way the gradient recovery is no longer necessary, but the price to pay is the increment of the unknowns of the problem due to the fact that a system of equations must be solved instead of a single scalar equation.

Hyperbolic First Order System formulation

The hyperbolic First Order System (FOS) formulation is here recalled for later convenience. The basic idea consists in re-writing the advection-diffusion scalar problem [START_REF] Abgrall | Construction of second-order accurate monotone and stable residual distribution schemes for unsteady problem[END_REF] as an equivalent first order system in witch the second order derivatives of the original problem are replaced by the first order derivatives of auxiliary variables. At the steady state the two formulations will coincide and the value of the auxiliary variables will equal the value of the derivatives of the unknown in the original problem.

Consider the case of a two-dimensional advection-diffusion problem for simplicity, the FOS formulation reads

∂u ∂t + a • ∇u = ν ∂p ∂x + ∂q ∂y ∂p ∂t = 1 T r ∂u ∂x -p ∂q ∂t = 1 T r ∂u ∂y -q (8)
where p and q are the gradient variables and T r is a relaxation time. At the steady state the system ( 8) is equivalent to the original equation ( 6), independently of the parameter T r , and p,q become equivalent to ∂u ∂x , ∂u ∂y respectively. Note that, differently form other schemes which use a first order representation of the advection-diffusion equation, the system (8) is hyperbolic. In vector form one has ∂u ∂t

+ A • ∇u = s, (9) with u
=     u p q     , A x =     a x -ν 0 -1 Tr 0 0 0 0 0     , A y =     a y 0 -ν 0 0 0 -1 Tr 0 0     , s =     0 -p Tr -q Tr     ,
and where, for an arbitrary vector n = (n x , n y ) T , one can write that A•n = A x n x + A y n y , thus

A n = A x n x + A y n y =     a n -νn x -νn y -nx Tr 0 0 - ny Tr 0 0     , (10)
with a n = a • n. The Jacobian matrix [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF] is diagonalizable with real eigenvalues

λ 1 = 1 2 a n -a 2 n + 4ν T r , λ 2 = 1 2 a n + a 2 n + 4ν T r , λ 3 = 0,
and the matrix of the right eigenvectors reads

R n =    -λ 1 T r -λ 2 T r 0 n x n y -n y n y n y n x    .
As usual, the Jacobian matrix can be written as

A n = R n Λ n L n
, where Λ n is the diagonal matrix of eigenvalues and L n = R -1 n is the matrix of the left eigenvectors.

Inria

The parameter T r can be defined as the ratio of a length scale L r to the fastest wave speed of the system, namely

T r = L r |a n | + ν/L r , (11)
while the length scale L r can be determined in order to ameliorate the formulation of the continuous system; for example in [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] L r is chosen such that the magnitude of the biggest and smallest eigenvalues are equal, thus minimizing the stiffness of the system. Since the system of equations ( 9) is hyperbolic, it can be discretized with any scheme already available for hyperbolic problems. If a RD scheme is used, the total residual on a generic element e is defined as follows

Φ e (u h ) =    Φ e u Φ e p Φ e q    = Ωe A•∇u h -s(u h ) dΩ.
The system is written in conservative form by introducing the flux function

f(u) = f x (u), f y (u) T , such that A = ∇ u f(u), f x (u) =     a x u -νp - u T ⋆ r 0     , f y (u) =     a y u -νq 0 - u T ⋆ r     ,
where the parameter T ⋆ r is used instead of the parameter T r , with

T ⋆ r = L r a + ν/L r ,
so that T ⋆ r is constant within the element, this ensure that at the steady state the relations ∂u ∂x = p and ∂u ∂y = q will be satisfied in the integral sense. The total residual can be now written as

Φ e = Ωe ∇•f(u h ) -s(u h ) dΩ = ∂Ωe f(u) • n d∂Ω - Ωe s(u h ) dΩ

Residual distribution discretization of advection-diffusion problems

In the previous sections the distribution process of the total residual is expressed through the use of generic distribution coefficients, in this section is described how actually to perform this step.

In the past years different RD schemes were developed with the objective to construct upwind schemes for linear triangular/tetrahedral elements, however the way how these schemes can be extended to different elements and/or different orders is still not obvious [START_REF] Abgrall | Toward the ultimate conservative scheme: Following the quest[END_REF][START_REF] Abgrall | Residual distribution schemes: Current status and future trends[END_REF]. In this work, the attention is focused on the construction of central schemes which can be formulated on every type of element and which can be easily extended to high order approximations. Linear and non-linear schemes are considered.

Central linear and non-linear RD schemes

Linear scheme

The linear scheme proposed in this work is the extension to the integral formulation of the classical Ni's Lax-Wendroff scheme [START_REF] Ni | A multiple grid scheme for solving the Euler equations[END_REF], namely

Φ e i = Φ e N e dof + Ωe a•∇ψ i τ a • ∇u h -∇• ν∇u h dΩ, (12)
where N e dof is the number of DOFs on the element e, the scaling parameter τ is defined as follows

τ = 1 2 |Ω e | j∈Σ e h max(k j , 0) , with k j = 1 2 ā•n j ,
and where ā represents the arithmetic average of the advection velocity on the element and the vector n j is defined as

n j = Ωe ∇ψ j dΩ.
The scheme ( 12) is linearity preserving but not positive, and due to the integral formulation, it is valid for any type of element and for any order of approximation. The scheme is conservative since i∈Σ e h Φ e i = Φ e , due to the fact that i∈Σ e h ∇ψ i = 0. The scheme is consistent. When the exact solution is injected in the Eq. ( 12) the residual is zero because the total residual vanishes, by definition, and the integral term vanishes as well due to the fact that the term between the brackets is exactly the governing equation.

Non-linear scheme

Non-linear schemes are needed to combine the non-oscillatory behavior of the numerical solution with the high order discretization. The basic idea to construct a non-linear scheme is to start with a first order, positive scheme, and to map its distributed residuals onto a set of positive and non-linear residuals.

To see in practice how to construct of a non-linear scheme, consider the first order accurate and positive Rusanov's scheme (also know as Lax-Friedrichs scheme)3 defined as

Φ e i = Φ e N e dof + 1 N e dof α j∈Σ e h j =i u i -u j ), ∀i ∈ Σ e h ,
with a ≥ max j∈Σ e h k j | + ν 0. Since the Rusanov's scheme is first order accurate, its distribution coefficients, β e i = Φ e i /Φ e , are unbounded. The construction of the non-linear scheme consists in mapping the distribution coefficients of the low order scheme onto non-linear bounded distribution coefficients βe i , this process is generally called limitation. A common choice for the map is the following [START_REF] Struijs | Fluctuation splitting schemes for the 2d Euler equations[END_REF]

βe i = max β e i , 0 j∈Σ e h max β e j , 0
.

The use of a central scheme, like the Rusanov's scheme, in combination with the limiting technique produces undamped spurious modes and a poor iterative convergence to the steady state solution [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF]. The cure to this problem consists in adding a filtering term by means of a streamline dissipation term

Φe i = βe i Φ e + θ e h (u h ) Ωe a•∇ψ i -∇• ν∇ψ i a • ∇u h -∇• ν∇u h dΩ. ( 13
)
The role of the parameter θ e h (u h ) is double. It provides the correct scaling of the streamline filtering and it makes sure that the filtering term is added only in the smooth regions of the solution. The following definition is used here

θ e h (u h ) = ε(u h )      j∈Σ e h |ā•n j | + ν 2|Ω e |      -1
, with ε(u h ) a smoothness sensor which behavior is ε(u h ) ≈ 1 in smooth regions and ε(u h ) ≈ 0 around discontinuities.

Improved discretization of the diffusion terms

Numerical experiments reveal that the schemes ( 12) and ( 13) applied to the discretization of the advection-diffusion problem are unsatisfactory from the point of view of the accuracy and the robustness. In order to obtain a better discretization of the diffusive terms, the advectiondiffusion equation ( 6) is written in the form of a first order system as follows

∇•f (u) -∇•(νq) = 0 q -∇u = 0 (14)
Consider now a numerical scheme for the previous system obtained by writing the weak form of the system plus a streamline stabilization term:

Ωe ψ i ∇•f (u h ) -∇•(νq) q -∇u h dΩ + Ωe A•∇ψ i τ ∇•f (u h ) -∇•(νq) q -∇u h dΩ = 0, (15)
where A = (A x , A y ) with

A x =   a x -ν 0 -1 0 0 0 0 0   and A y =   a y 0 -ν 0 0 0 -1 0 0   , so that A•∇ψ i =        a•∇ψ i -ν ∂ψ i ∂x -ν ∂ψ i ∂y - ∂ψ i ∂x 0 0 - ∂ψ i ∂y 0 0       
.

RR n°8157

The term τ is assumed to be of the following form

τ =   τ a 0 0 0 τ d 0 0 0 τ d   ,
where τ a and τ d are strictly positive coefficients. Supposing, now, that the gradient of the numerical solution has been recovered at each DOF, one can replace the second equation of the system ( 14) with the approximation ∇u u ≃ ∇u h and consider only the first equation, which now reads

Ωe ψ i ∇•f (u h ) -∇•(ν ∇u h ) dΩ + Ωe a•∇ψ i τ c a•∇u h -∇• ν ∇u h dΩ + Ωe ν ∇ψ i • τ d ∇u h -∇u h dΩ = 0. ( 16
)
The first two integrals of the previous equation represent a discretization for the scalar advection-diffusion equation by the means of a central scheme plus a streamline stabilization term, in the same way as shown in the Eq. ( 12). The last integral represents an additional stabilization term, for the diffusive part only, which vanishes in the advective limit and the parameter τ d is dimensionless. It is interesting to note that the additional term penalizes the difference between the discontinuous and the interpolated gradients, on each element.

With a slightly different procedure, a similar stabilization term for the diffusive part has been obtain by Nishikawa [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -I: Basic principle[END_REF] for the RD discretization of the diffusion problem.

With the Eq. ( 16) in mind, it is proposed here a modification of the schemes ( 12) and ( 13), previously introduced, in order to include the extra stabilization term for the diffusive part of the equation. In practice, the linear scheme reads

Φ e i = Φ e N e dof + Υ(Pe e ) Ωe a•∇ψ i τ a • ∇u h -∇• ν∇u h dΩ + 1 -Υ(Pe e ) Ωe ν ∇ψ i • ∇u h -∇u h dΩ, (17)
while the non-linear scheme becomes

Φe i = βe i Φ e + Υ(Pe e ) θ e h (u h ) Ωe a•∇ψ i -∇• ν∇ψ i a • ∇u h -∇• ν∇u h dΩ + 1 -Υ(Pe e ) Ωe ν ∇ψ i • ∇u h -∇u h dΩ, (18)
where it has been introduced the local Peclet number, define as Pe e = a h e /ν, with h e the characteristic length size of the element e, the function Υ(Pe e ) is defined such that Υ(Pe e ) → 0 in the diffusive limit and Υ(Pe e ) → 1 in the advective limit. In the numerical simulations the following definition is used

Υ(Pe e ) = max 0, 1 - 1 Pe e .
Note that in the schemes [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] or [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -I: Basic principle[END_REF], the use of the blending function Υ(Pe e ), makes possible to recover, in the case of the pure advection, the same scheme used for the discretization of pure advective problems, while in the case of pure diffusion problems only the stabilization term for the diffusive terms is taken into account.

Inria

Discretization of the hyperbolic first order system

The schemes introduced for the scalar advection-diffusion problem can be easily extended to case of a system of equations meaning that the discretization of the hyperbolic FOS is straightforward, with a simplification: there are no diffusive terms. The construction of non-linear scheme for hyperbolic system of equation has been analyzed in [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF], while the linear scheme for a system of equations reads

Φ e i = Φ e N e dof + Ωe A•∇ψ i Ξ A • ∇u h -s dΩ.
Here the scaling matrix Ξ is defined as

Ξ = 1 2 |Ω e |   j∈Σ e h K + j   -1 , with K + j = 1 2 R nj Λ + nj L nj ,
where the operator Λ + nj selects only the positive eigenvalues and sets to zero the negatives ones. With respect to the original work of Nishikawa, where only strong boundary conditions are considered for the solution and its gradient, here the boundary conditions are imposed in a weak sense as typicaly done for advection problems, like the compressible Euler equations for example. The total residual is first computed without considering the boundary contributions, then a correction residual is added to correctly take into account the boundary conditions. For a node i belonging to the boundary, the residual associated to the boundary conditions can be written as [START_REF] Abgrall | High order residual distribution scheme for rans equations[END_REF]]

Φ e i, ∂ = ∂Ωe∪∂Ω ψ i f(u ∂ ) -f(u h ) • n d∂Ω,
where u ∂ is the state that has to be imposed on the boundary and the term f(u ∂ )f(u h ) is a correction flux which vanishes when the solution on the boundary equals the state u ∂ . In this work the correction flux is taken as follows

f(u ∂ ) -f(u h ) = A - n (u h )(u ∂ -u h ), with A - n (u) = R n (u) Λ - n (u) L n (u).
Note that this correction flux corresponds to the classical Riemann flux used to impose the in/out flow boundary conditions for hyperbolic problems. The weak boundary conditions have been found to be more effective that the strong boundary conditions since the iterative convergence to the steady state solution is much faster.

Implementation details

From a numerical point of view, the total residual is computed by the means of quadrature formulas, for example Eq. ( 7) is approximated as

Φ e i ≃ l∈Γ e   N l quad q=1 f (u q ) -ν q ∇u q • n l q ω l q J l q ,  
where Γ e is the set of the faces of the elements e, N l quad and ω l q are, respectively, the number and the weights of the quadrature points on the l-th face, n l q is the outward normal versor to the face and J l q is the determinant of the Jacobian of the transformation from the reference element to the physical element. In the two-dimensional simulations considered in this work, two and three Gauss points are used on each face of linear and quadratic elements, respectively. Also the stabilization terms are numerically computed by using quadrature formulas; the number of the quadrature points depends on the kind of element considered, and is taken such that the functions are integrated exactly.

The explicit Euler scheme is used to iterate the scalar numerical scheme at the steady state. It uses a local time stepping, the scaled time step is choosen as follows

∆t n i = CFL e∈E h,i max j∈Σ e h |k n j | + ν ,
where CFL < 1 is a safety parameter. The same definition is used also for the hyperbolic FOS, where instead of the parameter k j the maximum value of the eigenvalues is used.

In numerical experiments, it has been observed that the high order discretization of the hyperbolic FOS converges very slowly to the steady state, making the use of a explicit scheme almost impossible. For this reasons, when quadratic elements are used, an implicit Euler scheme is employed in combination with the pseudo-transient continuation strategy for which the local time step is defined as follows

∆t n i = CFL n e∈E h,i max j∈Σ e h |k n j | + ν ,
with the CFL law taken as

CFL n = CFL 0 R n-2 L 2 R n-1 L 2 , with CFL 0 < 1,
where R n-1 L 2 and R n-2 L 2 are the L 2 norms of the residual at the time steps n-1 and n-2, respectively. In the simulation CFL 0 is takes as 0.9 and the maximum value of CFL is limited to 10 6 . The implicit problem is solved by the means of the inexact Newton-Krylov method. The GMRES algorithm with the ILU preconditioner is used to solve the resulting linear systems.

Gradient recovery strategies

As explained in Section 3, in order to construct the total residual for the whole advection-diffusion equation, one has to assume that a continuous value of the gradient of the numerical solution is available on the faces of the elements. The strategy adopted in this work to obtain a continuous approximation of gradient consists in recovering the gradients at every DOF of the grid and then the nodal values of gradients are interpolated with Lagrangian functions on each element.

The key point is the recovery of the gradients at the DOF, for this reason are here recalled some of the most used techniques in the field of the gradient recovery. Attention is focused on the possibility to obtain a high order gradient recovery, e.g., the gradient is recovered with the same order of accuracy of the solution. For simplicity, the description is always limited to the two-dimensional case, but the extension to the three spatial dimensions is straightforward.

Inria

Theory

Area-weighted method

One of the easiest way to recover the gradient at the grid nodes is the area-weight average of the gradients in each element surrounding a node, namely

∇u i = e∈E h,i ∇u h (x i ) |Ω e | e∈E h,i |Ω e | , ∀i ∈ Σ h ,
where |Ω e | is the area of the element e. The previous relation, in the case of linear elements is the so-called Green-Gauss formula.

L 2 -Projection

In the L 2 -Projection, the reconstructed gradients are obtained by solving the following equivalence ∇u h = ∇u h in a weak sense

Ω ψ ∇u h dΩ = Ω ψ∇u h dΩ, ∀ψ ∈ V h . (19)
From a numerical point of view, the weight function ψ is taken in the finite dimensional space of the Lagrangian functions. The gradient is expressed as follows

∇u h ≃ j∈Σ e h ψ j ∇u j ,
where ∇u j is the recovered gradient at the generic DOF j. If the following vectors of unknowns are defined, for the components of the gradient

x ∂ =   ∂u ∂x j   j=1,N dof and y ∂ =   ∂u ∂y j   j=1,N dof ,
than the discrete solution of the problem ( 19) can be obtained by solving the following linear systems

M x ∂ = b x and M y ∂ = b y , with M ij = Ωij ψ i ψ j dΩ, b xi = Ωi ψ i ∂u h ∂x dΩ, and b yi = Ωi ψ i ∂u h ∂y dΩ,
where Ω i is the support of the shape function ψ i and Ω ij = Ω i ∩ Ω j . The gradient ∇u h is computed by resorting the gradient of the shape functions, as standard practice in the Finite Element field. This technique requires the inversion of a global linear system that can be quite expensive for a high number of unknowns. Obviously, since the matrix M depends only on the geometry of the grid, it can be inverted only once and can be used for several calculations on the same grid.

Least-square method

Another approach to recover the gradient at each DOF of the grid is the least-square reconstruction. The technique is unrelated to the mesh topology and it involves only the information associated to the neighboring nodes. Although the stencil is arbitrary, the natural choice involves only the nearest neighboring nodes.

The starting point consists in expanding the solution in a Taylor series around the node i for each node j belonging to the stencil of i, see Fig. 1,

u j = u i + ∂u ∂x i (x j -x i ) + ∂u ∂y i (y j -y i )+ ∂ 2 u ∂x 2 i (x j -x i ) 2 + ∂ 2 u ∂y 2 i (y j -y i ) 2 + ∂ 2 u ∂x∂y i (x j -x i )(y j -y i ) + . . . , (20)
where u i = u(x i ) and u j = u(x j ). The gradient reconstruction is obtained by solving for the values of the gradients that minimize the following function

N j=1 ω 2 ij E 2 ij , ∀ i ∈ Σ h with E 2 ij = -∆u ij + ∂u ∂x i ∆x ij + ∂u ∂y i ∆y ij + ∂ 2 u ∂x 2 i ∆x 2 ij + ∂ 2 u ∂y 2 i ∆y 2 ij + ∂ 2 u ∂x∂x i ∆x ik ∆y ij + . . . 2 ,
where ∆u ij = u j -u i , ∆x ij = x j -x i , ∆y ij = y j -y i , while ω ij is a weight factor. In the case of linear elements, the solution is expanded only up to the first derivatives in the Taylor series and the components of the gradient are obtained by solving the following minimization problems for the first derivatives

∂   N j=1 ω 2 ij E 2 ij   ∂ ∂u ∂x i = 0 and ∂   N j=1 ω 2 ij E 2 ij   ∂ ∂u ∂y i = 0.
By simple algebra, it is easy to see that the previous minimization problems correspond to the solution of following small linear system

        N j=1 ω 2 ij ∆x 2 ij N j=1 ω 2 ij ∆x i ∆y ij N j=1 ω 2 ij ∆x ij ∆y ij N j=1 ω 2 ij ∆y 2 ij                 ∂u ∂x i ∂u ∂y i         =         N j=1 ω 2 ij ∆x ij ∆u ij N j=1 ω 2 ij ∆y ij ∆u ij         .
The weight factor ω ij is generally taken as the inverse of the distance between the nodes i and j.

The extension to the case of quadratic elements is straightforward, it consist in taking also the second derivatives in the Taylor series and the minimization is done respect both first and second derivatives.

Super-convergent patch recovery

In the field of the Finite Element method applied to the mechanical structures, it is know that the stresses (gradients of the displacements) sampled at certain points in a element possess a super-convergent property. This means that the stresses have the same order of accuracy of the displacements [START_REF] Zienkiewicz | Finite Element Method (5th Edition) Volume 1 -The Basis[END_REF]. It can be shown that in the case of a segment element such particular points correspond to the Gauss-Legendre points [START_REF] Herrmann | Interpretation of finite element procedure as stress error minimization procedure[END_REF], obviously by tensor product such points can be defined also in the case of quadrangles and hexahedrons. For triangles or tetrahedrons such property cannot be rigorously shown, at least up to our knowledge, but numerical experiments confirm that the stresses sampled at certain points have high order of accuracy.

Accepting the fact that gradients are sampled with high order accuracy in certain points of the element, it is possible to compute gradients which are high order accurate within all the element. Indeed, if at sampling points the value of gradients is accurate to order k + 1, by using a polynomial of degree k (the same order used to interpolate the solution), it is possible to obtain an approximation which has high order accuracy everywhere within the elements if this polynomial is made to fit the values of the sampled gradients in a least square manner. Such a technique is called super-convergent patch recovery introduced by Zienkiewicz and Zhu (SPR-ZZ) [START_REF] Zienkiewicz | A simple error estimator and adaptive procedure for practical engineerng analysis[END_REF][START_REF] Zienkiewicz | The superconvergent patch recovery and a posteriori error estimates. part 2: Error estimates and adaptivity[END_REF].

Assume that the numerical solution u h of the problem is known at each DOF of the grid to the k+1-th order of accuracy. The aim is to obtain the values of the solution gradient, ∇u h at all the DOFs with same order of accuracy of the solution. The components of the recovered gradient, at the generic DOF i, are written in a polynomial form as follows Assuming that N s sampling points, ξ ℓ , ℓ = 1 . . . N i s , are available for each DOF i, the objective is to minimize the following functions

F x = N i s k=1 ∂u h ∂x (ξ k ) -p T k a x and F y = N i s k=1 ∂u h ∂y (ξ k ) -p T k a y ,
with p k = p(ξ k ). The vectors of the coefficients a x and a y are obtained by solving the following minimization problems ∂F x ∂a x = 0 and ∂F y ∂a y = 0.

We discuss the structure of the sampling points a little bit later in the text for clarity reasons.

It is easy to verify that the minimization problems correspond to the solution of the following linear systems

A T Aa x = A T b h x , and
A T Aa y = A T b h y , (21)
where

b h x =            ∂u h ∂x (ξ 1 ) ∂u h ∂x (ξ 2 )
. . .

∂u h ∂x (ξ N )            , b h y =            ∂u h ∂y (ξ 1 ) ∂u h ∂y (ξ 2 )
. . .

∂u h ∂y (ξ N )            and A =        1 x 1 y 1 . . . y k+1 1 1 x 2 y 2 . . . y k+1 2 . . . . . . . . . . . . . . . 1 x N y N . . . y k+1 N        .
To compute the coefficients a x and a y , a small linear system must be solved for each DOF i of the grid. The dimension of the matrix A i are determined by the number of sampling points N i s and by the degree of the polynomials used to express the recovered gradient, that is

A i ∈ R N i s ×m
, where m is the number of the coefficients in the vector a x or a y . The problems in the Eq. ( 21) admit a unique solution if Rank A i = m, which is always satisfied in the case in which N i s ≥ m. It is worth also noticing that since the matrix A i depends only on the geometry, for a given grid the matrix (A T i A i ) -1 A T i needs to be computed only once.

Generally, the number of elements which share the same node within the domain is such that the condition N s ≥ m is always satisfied, this means that the gradient recovery is compact because it involves only the elements contained within the support of a grid node. For the nodes belonging to the boundary of the grid the condition N s ≥ m might not be satisfied without enlarging the stencil, otherwise the problem is ill conditionned. In this case, to avoid the use of larger stencil for a boundary node -which would break the compactness requirement-it is possible to obtain the value of the recovered gradient with the same polynomial expansion used for nearest domain node. Now we discuss the structure of the sampling set for this gradient reconstruction technique. For each degree of freedom i, the set S i of sampling points is defined by

S i = ∪ e∋i S e i
where the sampling point for the element e depend on the structure of the element. In this paper, we have considered quads and triangles (the extension to 3D is rather obvious), hence the sampling points ξ are as defined on Figre 2. In Fig. 2 are shown examples of patches used to recover the gradient for a domain node in the case of quadrangular and triangular elements.

For a quadrangle the sampling points are defined uniquely. Considering a reference segment defined as x = [-1, 1], the sampling point is the point x = 0 in the case of a linear element, while in the case of a quadratic element the sampling points have coordinates ±1/ √ 3. The sampling points on the reference quadrangle are simply obtained by a tensorial product of the points defined on the reference segment.

For a linear triangle, the sampling point is the point with barycentric coordinated λ = (1/3, 1/3, 1/3), while in the case of a quadratic triangle the sampling points are not unique, different choices are available. In Fig. 3 are shown three examples of sampling points used. In the first option (Fig. 3-(a)) are used three points with barycentric coordinates

λ 1 = 2 3 , 1 3 , 1 3 , λ 2 = 1 3 , 2 3 , 1 3
, and

λ 3 = 1 3 , 1 3 , 2 3 ,
In the second option (Figure 3-(b)) four points are used with barycentric coordinates

λ 1 = 1 3 , 1 3 , 1 3
, λ 2 = (0.6, 0.2, 0.2) , λ 3 = (0.2, 0.6, 0.2) , λ 4 = (0.2, 0.2, 0.6) .

Another option (Fig. 3-(c)) consists in taking as sampling points the three points with barycentric coordinates

λ 1 = 1 2 , 1 2 , 0 , λ 2 = 0, 1 2 , 1 2 , λ 3 = 1 2 , 1 2 , 0 . RR n°8157 (a) (b) (c)

Results and discussion

To study the accuracy of the presented gradient recovery strategies the following function is used

u = -cos(2πη) exp ξ 1 - √ 1 + 16π 2 ν 2 2ν , (22)
with η = a y x -a x y and ξ = a x x + a y . Here a x = 0.5, a y = √ 3/2 and ν = 0.01. The solution, shown in Fig. 4, is infinitely differentiable with continuous gradients. The computations of the recovered gradients are performed on four different kind of grids, shown in Fig. 5, namely unstructured grids of triangles, quadrangles and hybrid elements and highly deformed unstructured meshes of triangles, obtained randomly perturbing a regular grid.

The error of the recovery procedure is computed as the L 2 norm of the difference between the computed gradient, ∇u, and the exact gradient, ∇u ex , for each spatial component, namely

ǫ L 2 = Ω ∇u -∇u ex 2 dΩ Ω ∇u ex 2 dΩ .
Table . 1 shows the L 2 errors and the orders of convergence of different gradient recovery methods on triangular grids. In the case of linear elements, the differences between the recovery methods is small and all the schemes reach almost the second order accuracy. The L2-Projection and the SPR-ZZ methods have the smallest level of error, but the former scheme is much more expensive because it requires the solution of a global linear system. In the case of quadratic elements, the methods have an order of accuracy no more than two except for the SPR-ZZ method which show almost third order accuracy, it is also worth noticing that the errors obtained with this method are one order of magnitude smaller of those obtained with other methods. In Table . 2 are reported the errors obtained with the SPR-ZZ procedure on quadratic triangular elements for the three different sampling strategies shown in Fig. 3. It is evident that the first strategy guaranties the smallest level of error, while the four-points strategy is unsatisfactory. Note that the results reported in the previous paragraph are obtained with the first strategy.

The errors of the recovery methods on unstructured grids of quadrangles and of hybrid elements are reported in Table . 3 and Table . 4, respectively. The behavior of the recovery methods is the same observed in the case of triangular grids. Table . 5 shows the errors computed on a 22) used to test the gradient recovery methods. sequence of highly distorted triangular grids, the performance of the recovery methods is not optimal anymore due to very poor quality of the meshes, nevertheless the errors obtained with the SPR-ZZ methods are always much smaller than those obtained with other procedures.

Numerical experiments

The section presents an extensive evaluation of the numerical schemes proposed. The objective is to show that 1. the high order RD schemes previously proposed can be successfully used in the discretization of the advection-diffusion equation, 2. the high order accuracy is preserved in all the range of the Peclet number. This is contrast with the method proposed in [START_REF] Abgrall | Numerical approximation of parabolic problems by means of residual distribution schemes[END_REF] where the region P e ≈ 1 was problematic.

In all the simulations, the steady state is considered to be reached when the L 2 norm of the initial residual is reduced by ten orders of magnitude. If the residual of the scheme stagnates at a high level, it is marked that the simulation is not converged. The CFL number is taken as 0.9 and 0.6, respectively for the linear and the non-linear schemes. The same kind of grids shown in Fig. 5 are considered in the numerical simulations.

Linear advection-diffusion equation

To verify the order of accuracy of the linear and non-linear schemes, as well as to study the influence of the accuracy of the gradient recovery methods on the accuracy of the numerical solution, the linear advection-diffusion problem with constant viscosity is considered here, namely

a•∇u = ν∇•∇u, on Ω = [0, 1] 2 ,
the exact solution of the problem reads

u = -cos(2πη) exp ξ 1 - √ 1 + 16π 2 ν 2 2ν ,
with η = a y x -a x y and ξ = a x x + a y y. Here a = (0, 1) T and ν = 0.01, which is the most critical case because the advection and the diffusion have similar orders of magnitude and traditional high order RD schemes generally loose an order of accuracy in this regime.

Inria

Linear elements Quadratic elements On the left, right and bottom boundaries of the domain the exact solution is imposed as Dirichlet boundary condition while on the top boundary nothing is done. The solution is initialized with a zero value everywhere in the domain except on the inflow boundaries where the exact solution is imposed.

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u
In Table . 6 are reported the L 2 errors of the numerical solution obtained on a sequence of triangular grids with the linear scheme [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] and the non-linear one [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme -I: Basic principle[END_REF], for different gradient recovery strategies; linear and quadratic elements are considered. For sake of completeness, in Fig. 6 and in Fig. 7, the errors of the solution are shown together with the errors of the gradients of the numerical solution, for the linear and quadratic approximation, respectively. In the case of linear elements, the accuracy of the schemes with different gradient recovery methods is almost identical, for both linear and non-linear schemes. This is in accordance with the accuracy results observed for the different gradient recovery techniques and it underlines also the fact that high cost of the L2-projection method is not justified, since less expensive methods produce results with the same level of accuracy. This is also in agreement with the observation made by the other researchers.

The situation is very different in the case of quadratic elements, the weighted area and the

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2
Sampling points Fig. 3-( L2-Projection gradient recovery methods produces a sub-optimal scheme with second order only accurate solutions for both linear and non-linear schemes. The use of the SPR-ZZ method allows the construction of an optimal third order accurate scheme and it is worth noticing that also the x-component of the gradient of the numerical solutions is third order accurate, meaning that solution and gradients are computed with the same order of accuracy. On the y-component of the gradient this optimal behavior is lost, and this due to the combined effects of the gradient recovery with the solution error on the outflow boundary, where no boundary condition is imposed. An optimal accuracy on both the components of the gradient has been observed in numerical simulations of the linear advection-diffusion problem with Dirichlet boundary conditions imposed on all the boundaries of the domain.

It is worth noticing that the combination of the non-linear scheme with the least square gradient recovery technique produce an almost optimal scheme although the least square recovery does not allow a high order gradient recovery, by itself.

In Fig. 8 are reported the errors of the solution and of gradient components obtained discretizing the linear advection-diffusion problem on a sequence of unstructured grid of quadrangles with the linear and the non linear schemes. For simplicity, only the weighted area and the SPR-ZZ recovery strategies are used. As previously observed, with linear elements there is no significant difference in the level of accuracy between different gradient recovery techniques, however with quadratic elements only the SPR-ZZ gradient recovery guarantees third order accurate solutions. Of course, the same considerations done for triangular and quadrangular grids are still valid with mesh with hybrid elements, as it is evident from Fig. 9.

In Fig. 10 are reported the errors obtained on a sequence of highly distorted triangular grids. The behavior of the schemes is similar to that observed with more regular meshes, in particular it is important to note that the poor quality of the grids has only a limited influence on the

Inria

Linear elements

Quadratic elements

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u
N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2
Area-eight method Table 4: Accuracy study of different gradient recovery methods with linear and quadratic elements on unstructured grid of hybrid elements.

Inria

Linear elements Quadratic elements Table 5: Accuracy study of different gradient recovery methods with linear and quadratic elements on grids of randomly distrorted triangles.

N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u ∂y O L 2 N dof ǫ L 2 ∂u ∂x O L 2 ǫ L 2 ∂u
Linear elements Quadratic elements

N dof ǫ L 2 (u h ) O L 2 ǫ L 2 (u h ) O L 2 N dof ǫ L 2 (u h ) O L 2 ǫ L 2 (u h ) O L 2
Weighted area (2.21) L2-Projection (

Weighted area (1.78) L2-Projection (1.84) Least square (1.78) SPR-ZZ (

Weighted area (1.16) L2-Projection (1.39) Least square (1.12) SPR-ZZ (1.01)

Figure 6: L 2 error in the solution of the linear advection-diffusion problem on triangular girds with linear elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

accuracy of the scheme. In order to highlight the effectiveness of higher order schemes respect to the second order schemes, in Fig. 11 the discretization errors of the solution is reported versus the number of DOF and versus the CPU time, for brevity only results for the linear scheme on triangular grids are shown. One can see that to get a fixed level of error, let's say 10 -5 , an actual third order scheme requires about 12 000 DOFs and 25 minutes to perform the computation. A second order scheme, on the other hand, requires about 31 000 DOFs and 5 hours to get the same level of error.

The effect of the stabilization term obtained for the viscous part only is now investigate, log

Weighted area (2.02) L2-Projection (2.03) Least square (2.65) SPR-ZZ (

Weighted area (2.03) L2-Projection (2.03) Least square (2.00) SPR-ZZ (2.96)

log(h) log(L 2 error) -2.2 -2 -1.8 -1.6 -1.4 - 1.2 -3.5 -3 -2.5 -2 -1.5
Weighted area (2.19) L2-Projection (2.07) Least square (2.79) SPR-ZZ (

Weighted area (1.01) L2-Projection (1.46) Least square (2.14) SPR-ZZ (

Weighted area (1.95) L2-Projection (1.45) Least square (2.05) SPR-ZZ (2.80)

Figure 7: L 2 error in the solution of the linear advection-diffusion problem on triangular girds with quadratic elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and

h = 1/ √ N dof .
this means that the linear scheme ( 17) is compared against the scheme ( 12) and the non-linear scheme ( 18) is compared against the scheme (13). The comparison is done in term of solution accuracy and number of iterations necessary to reach the steady state, results are reported in Table . 7 for the linear and non-linear schemes, with the SPR-ZZ recovery strategy. It can be observed that in the case of linear elements there is no appreciable difference in term of error between the schemes with and without the stabilization term for the viscous part, however the use of the extra stabilization term makes the linear scheme converge much faster to the steady state. The effect becomes even more important in the case of non-linear scheme, where

Inria log(h) log(L 2 error) -2 -1.5 -1 -5 -4 -3 -2 -1
Weighted area (2.04) SPR-ZZ (2.03) Weighted area (2.10) SPR-ZZ ( it is observed that the absence of the extra dumping term prevents the scheme to converge in most cases. For the quadratic elements, the presence of the extra dumping term has two effects, it improves the convergence of the numerical methods and it introduces a crucial improvement in the level of accuracy, for both linear and non-linear schemes.

For sake of completeness, the linear advection-diffusion problem is also solved with a very small viscosity coefficient, ν = 10 -6 , in order to verify that the numerical schemes are able to preserve the theoretical accuracy in the advection limit. The errors, obtained on a sequence of triangular grids, are reported in Fig. 12 for the linear and the non-linear schemes. As expected, Weighted area ( In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Discretization with the hyperbolic FOS

The linear advection-diffusion problem is now discretized by the means of the hyperbolic FOS scheme described in section 4.3. The objective is to compare the accuracy and the performance of this formulation with the standard scalar discretization. The linear scheme is used to discretize the hyperbolic FOS and the scalar equation, in the latter case the SPR-ZZ gradient recovery strategy is chosen. A sequence of triangular grids is considered, with linear and quadratic elements, and the viscosity coefficient ν is takes as 0.01. In Fig. 13 With linear elements, there is only a small difference in the errors discretization between the scalar and the FOS formulation, however considering the CPU time versus the levels of error it is evident that the scalar scheme is much more effective than the FOS formulation. Note that the CPU time is measured in seconds and is reported in logarithmic scale; for example on the finest mesh the scalar scheme requires about ten minutes to reach the steady state, while the FOS scheme requires about two hours. The slope of the curve CPU time-error is about -2/3, for both schemes. In [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] was observed a slope -2/3 for the FOS scheme while a slope -1/2 was observed for a scalar scheme based on the LDA RD plus a Galerkin scheme.

With quadratic elements, the accuracy on the solution obtained using the FOS scheme is slightly better than that obtained with the scalar scheme but the situation is completely different if one looks at the accuracy of the gradient. As already pointed out in [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF], the discretization of the hyperbolic FOS with RD schemes does not allow to recover the gradients with the same accuracy of the solution, unless the mesh is regular. It can be noticed how the use of the SPR-ZZ strategy allows to obtain a third order accurate gradient while the FOS scheme gives only second order accuracy. The last remark concerns the CPU time. An implicit Euler method has been used for the FOS scheme and an explicit Euler method has been used for the scalar scheme; the scalar scheme is still much more effective than the FOS scheme. Note that in the advection limit the smallest eigenvalue of the hyperbolic FOS vanishes, this means that two of the three eigenvalues are zero and the problem becomes ill conditioned. It has been observed that the FOS scheme is not able to converge for the linear advection-diffusion problem with the viscous coefficient ν taken as 10 -6 .

Inria

Linear elements

Quadratic elements

N dof ǫ L 2 (u h ) Ite ǫ L 2 (u h ) Ite N dof ǫ L 2 (u h ) Ite ǫ L 2 (u h )

Viscous Burger equation

The viscous Burger equation is now considered in order to test the accuracy of the numerical schemes with a non-linear problem. The governing equation reads

∂ ∂x u 2 2 + ∂u ∂y = ν ∂ 2 u ∂x 2 , on Ω = [0, 1] 2 ,
the problem admits the following exact solution u = 2νπ exp(-νyπ 2 ) sin(πx) a + exp(-νyπ 2 ) cos(πx) , with a > 1.

Note that the exact solution of the steady two-dimensional problem is obtained from the unsteady one-dimensional problem, in which the time coordinated is substituted by the y coordinates. In the simulations the parameter a is taken as 1.5 and the viscosity ν coefficient is taken as 0.05. On the bottom, left and right boundaries the exact solution is imposed as Dirichlet boundary condition. The solution is initialized with a zero values everywhere, except on the inflows boundaries where the exact solution is imposed. A sequence of unstructured triangular grid is considered, the weighted area and the SPR-ZZ gradient recovery methods are used. The errors of the solution are reported in Table. 8, while in Fig. 14 the errors of the solution are shown together with the error of the gradients components.

With linear elements, the level of accuracy of the schemes is almost identical, while with quadratic elements the situation is very different. The use of weighted area gradient recovery has a disastrous effect on the accuracy of the solution, indeed the theoretical third order scheme has the same level of accuracy of the second order scheme. On the other had, the use of the SPR-ZZ recovery method allows to construct an optimal third order scheme and also the accuracy of the gradients is improved. , on triangular girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower).

In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Anisotropic diffusion problem

As last test case, a two-dimensional diffusion problem is addressed, the viscosity is not considered to be a scalar anymore but a tensor quantity and is taken to be anisotropic. The aim of this test case is to study the accuracy of the proposed RD schemes with a pure diffusion problem, and the anisotropy of the viscous tensor is introduced to test also the robustness of the numerical scheme. Figure 13: L 2 error in the solution of the linear advection problem with ν = 0.01 on triangular girds, with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) CPU time in seconds versus the error (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

The diffusion problem is formulated as follows

-∇• K∇u = 0, on Ω = [0, 1] 2 , with K = 1 0 0 δ ,
the problem has the following exact solution u = sin(2πx) e -2πy √ 1/δ , Linear elements Quadratic elements and in the numerical simulations δ is takes as 10 3 . A sequence on unstructured triangular grids is considered, the scalar schemes with the weighted area and SPR-ZZ gradient recovery strategies are considered. In Fig. 15 are reported the errors of the solution and of the gradients for linear and quadratic elements. As usual, the second order schemes have the same level of accuracy independently of the gradient recovery method used, but with quadratic elements only the use of the SPR-ZZ method allows to get a third order accuracy on the solution as well as on the gradients. It is interesting to note that the accuracy of the non-linear scheme is severely spoiled by the use of the simple weighted area method with quadratic elements.

N dof ǫ L 2 (u h ) O L 2 ǫ L 2 (u h ) O L 2 N dof ǫ L 2 (u h ) O L 2 ǫ L 2 (u h ) O L 2
In the end, the anisotropic viscous problem is solved on a uniform, structured mesh of quadrangles, results are shown in Fig. 16 and indicated that the theoretical accuracy of the schemes is achieved independently from the gradient recovery technique used, although the accuracy of the gradient is always better with the SPR-ZZ method. This remark is important because for advection-diffusion problems, even the use of uniform, structured grids does not preserve the formal accuracy of the scheme unless the gradients are recovered with high order accuracy.

Conclusion

An high order accurate and robust Residual Distribution scheme for the solution of advectiondiffusion equations has been presented. The method relies on the computation of a total residual for the whole equation without construing two different type of schemes for the advection and diffusion parts. A fundamental aspect of the construction to get an high order approximation of the solution is the high order recovery of the gradient of the numerical solution. Different recovery techniques has been analyzed, in particular the super-convergent patch recovery of Zienkiewicz and Zhu has shown to posses a such level of the flexibility and accuracy to guaranty the construction of third order accurate RD schemes for general unstructured grids. The accuracy of the numerical schemes has been verified with linear, non-linear advection-diffusion problems Figure 16: L 2 error in the solution of the anisotropic diffusion problem on uniform, structured grid of quadrangles, with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower).

In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Journal of Computational Physics, 214:773-808, 2006.

Figure 1 :

 1 Figure 1: Illustration of the stencil for least square gradient recovery at the node i. The symbol (•) indicates the node around which the Taylor series expansion is done while the symbols (•) indicate the node used to construct the least square problem. On the left stencil with linear elements, on the right stencil with quadratic elements for a nodes on the face.

  ∂u h ∂x i = p T a x and ∂u h ∂y i = p T a y , with p T (x) = (1, x, y, x 2 , . . . , x k+1 , x k y, . . . , y k+1 ), a x = (a x1 , a x2 , . . . , a xm ) and a y = (a y1 , a y2 , . . . , a ym ).

Figure 2 :

 2 Figure 2: Interior super-convergent patches for quadrilateral and triangular elements: top linear elements, bottom quadratic elements. The symbols (•) indicate the patch assembly points, the symbols (•) indicate the points where the gradient is recovered and the symbols (△) indicate the super-convergent sampling points.

Figure 3 :

 3 Figure 3: Three different examples of sampling points for quadratic triangles.

Figure 4 :

 4 Figure 4: Graphical representation of the Eq. (22) used to test the gradient recovery methods.

Figure 5 :

 5 Figure 5: Example of different kinds of grid used to test the accuracy of the gradient recovery procedures.

Figure 8 :

 8 Figure8: L 2 error in the solution of the linear advection-diffusion problem on quadrangular girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Figure 10 :

 10 Figure10: L 2 error in the solution of the linear advection-diffusion problem on highly distorted triangular girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

  are shown the errors on the log(N dof )

Figure 11 :

 11 Figure 11: L 2 error in the solution of the linear advection-diffusion problem, on triangular grids, versus the number of DOFs and the CPU time in seconds.

Figure 12 :

 12 Figure12: L 2 error in the solution of the linear advection-diffusion problem, with ν = 10 -6 , on triangular girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Inria

Figure 14 :

 14 Figure14: L 2 error in the solution of the viscous Burger problem on triangular girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Figure 15 :

 15 Figure15: L 2 error in the solution of the anisotropic diffusion problem on triangular girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .

Table 1 :

 1 Accuracy study of different gradient recovery methods with linear and quadratic elements on an unstructured grids of triangles.

	∂y	O L 2

Table 2 :

 2 Accuracy study of SPR-ZZ recovery methods with quadratic triangular elements for different sampling strategies.

	a)
	445 5.4353E -03	-	5.3235E -03	-
	1705 7.1348E -04 3.02 6.9027E -04 3.04
	6673 9.5508E -05 2.94 1.0076E -04 2.82
	26401 1.3631E -05 2.83 1.6522E -05 2.62
	105025 2.2311E -06 2.62 3.0462E -06 2.44
	Sampling points Fig. 3-(b)
	445 5.4270E -03	-	5.9338E -03	-
	1705 7.4236E -04 2.96 7.6624E -04 3.04
	6673 1.1574E -04 2.72 1.2145E -04 2.69
	26401 2.1052E -05 2.47 2.2411E -05 2.45
	105025 4.4137E -06 2.26 4.6939E -06 2.26
	Sampling points Fig. 3-(c)
	445 7.4535E -03	-	7.8030E -03	-
	1705 1.5611E -03 2.32 1.5747E -03 2.38
	6673 3.7352E -04 2.09 3.7461E -04 2.10
	26401 9.2406E -05 2.03 9.2623E -05 2.03
	105025 2.3049E -05 2.01 2.3110E -05 2.01

Table 3 :

 3 Accuracy study of different gradient recovery methods with linear and quadratic elements on unstructured grids of quadrangles.

	Linear elements	Quadratic elements

Table 6 :

 6 L 2 errors and orders of accuracy in the solution of the linear advection-diffusion problem on triangular girds with the linear and non-linear schemes, for different gradient recovery strategies, and with linear and quadratic elements.

	linear scheme		non-linear scheme		linear scheme		non-linear scheme
	Area-weight method
	43 1.2149E -01	-	3.0392E -01	-	149 1.1002E -02	-	2.7329E -02	-
	121 3.0269E -02 2.68	6.9542E -02	2.85	445 3.0503E -03 2.34	7.0412E -03	2.47
	445 7.4341E -03 2.15	1.7157E -02	2.14	1705 8.3842E -04 1.92	1.6639E -03	2.14
	1705 2.3405E -03 1.72	4.2609E -03	2.07	6673 2.2112E -04 1.95	4.2057E -04	2.01
	6673 5.8700E -04 2.02	1.0452E -03	2.05	26401 5.4878E -05 2.02	1.0748E -04	1.98
	L2-Projection method
	43 1.3837E -01	-	3.0722E -01	-	149 1.1434E -02	-	2.7277E -02	-
	121 3.6110E -02 2.59	7.0099E -02	2.85	445 3.0027E -03 2.44	1.1651E -02	1.55
	445 8.6721E -03 2.19	1.7306E -02	2.14	1705 8.4470E -04 1.88	2.4454E -03	2.32
	1705 2.3478E -03 1.94	4.2515E -03	2.09	6673 2.2191E -04 1.95	5.7123E -04	2.13
	6673 5.9575E -04 2.01	1.0418E -03	2.06	26401 5.4986E -05 2.02	not converged
	Least square method
	43 1.2013E -01	-	3.0309E -01	-	149 2.4329E -02	-	2.0318E -02	-
	121 2.8616E -02 2.77	6.9005E -02	2.86	445 5.0608E -03 2.87	3.7344E -03	3.09
	445 6.8597E -03 2.19	1.7308E -02	2.12	1705 8.6129E -04 2.63	4.6880E -04	3.08
	1705 2.3311E -03 1.60	4.2958E -03	2.07	6673 1.1362E -04 2.96	7.6648E -05	2.65
	6673 5.8774E -04 2.01	1.0585E -03	2.05	26401 2.7760E -05 2.04	1.5193E -05	2.35
	SPR-ZZ method
	43 1.2639E -01	-	3.0257E -01	-	149 1.1247E -02	-	1.3349E -02	-
	121 3.3130E -02 2.58	7.2141E -02	2.77	445 1.8777E -03 3.27	1.8975E -03	3.56
	445 8.2461E -03 2.13	1.8068E -02	2.12	1705 1.9648E -04 3.36	2.2616E -04	3.16
	1705 2.3337E -03 1.87	4.4622E -03	2.08	6673 2.3797E -05 3.09	2.9410E -05	2.99
	6673 5.8201E -04 2.03	1.0974E -03	2.05	26401 3.5754E -06 2.75	4.6791E -06	2.67
							Inria

  Figure9: L 2 error in the solution of the linear advection-diffusion problem on hybrid girds with linear (dashed lines) and quadratic (solid lines) elements. Error of the solution (first column), error of the x-component of the gradient (second column) error of the y-component of the gradient (third column). Linear scheme (upper), non-linear scheme (lower). In the legends are reported also the mean slopes of curves and h = 1/ √ N dof .the theoretical accuracy of the schemes is reached independently from the gradient recovery method used, because the diffusive effects are negligible in this regime, nevertheless the level of accuracy of the gradients obtained with the SPR-ZZ technique is better. Note that the error of the gradients are not normalized in this case, due to the very small value of y-component of the gradient.

				Weighted area (1.99)			0		Weighted area (1.33)	1.55)	0	Weighted area (1.08)
				SPR-ZZ (2.01)						SPR-ZZ (1.85)	SPR-ZZ (2.02)			SPR-ZZ (1.01)
		-1		Weighted area (2.11)					Weighted area (2.06) Weighted area (1.99)		Weighted area (1.24)
				SPR-ZZ (2.81)						SPR-ZZ (2.63)	SPR-ZZ (2.80)			SPR-ZZ (1.52)
								-1				-0.5
		-2
													-1
	error)						error)	-2				error)
	log(L 2	-3					log(L 2					log(L 2 -1.5
								-3
		-4											-2
								-4				-2.5
		-5
		-2		-1.5	-1			-2		-1.5	-1		-2		-1.5	-1
				log(h)						log(h)					log(h)
		-1		Weighted area (2.18) SPR-ZZ (2.20) Weighted area (2.07) Weighted area (1.84) SPR-ZZ (2.04) Weighted area (2.25) -1 SPR-ZZ (3.11) SPR-ZZ (2.77)		0		0 Weighted area (1.38) Weighted area (1.66) SPR-ZZ (1.96) Weighted area (2.04) SPR-ZZ (1.80) Weighted area (2.03) SPR-ZZ (2.72) SPR-ZZ (2.43)	0	0 Weighted area (1.06) Weighted area (1.18) SPR-ZZ (0.82) Weighted area (1.96) SPR-ZZ (1.08) Weighted area (1.41) SPR-ZZ (1.96) SPR-ZZ (1.78)
								-1		-1		-0.5	-0.5
				-2
		-2											-1	-1
	error) log(L 2	-3	log(L 2 error)	-3			error) log(L 2	-2	error) log(L 2	-2		error) log(L 2 -1.5	log(L 2 error) -1.5
										-3				-2
				-4				-3					-2
		-4												-2.5
										-4		-2.5
				-5				-4						-3
		-5											-3
		-2	-2 -1.5	-1	-1.5	-1	-2	-2 -1.5	-1.5 -1	-1	-2	-1.5	-2	-1	-1.5	-1
				log(h)	log(h)				log(h)	log(h)				log(h)	log(h)

Table 7 :

 7 L 2 errors and orders of accuracy in the solution of the linear advection-diffusion problem on triangular girds with the linear and non-linear schemes, for different gradient recovery strategies, and with linear and quadratic elements.

	Ite

Table 8 :

 8 L 2 errors and orders of accuracy in the solution of the viscous Burger problem on triangular girds with the linear and non-linear schemes, for different gradient recovery strategies, and with linear and quadratic elements.

RR n°8157

Other low order, non oscillatory schemes can considered, like for example a Finite Volume scheme written as a RD scheme, see e.g.[START_REF] Abgrall | Toward the ultimate conservative scheme: Following the quest[END_REF].Inria

1.0519E -01 -1.0224E -01 -537 1.1979E -02 -1.0599E -02 -

Inria

351, Cours de la Libération Bâtiment A 29

Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399