R Abgrall
email: remi.abgrall@inria.fr

D De Santis

M Ricchiuto

Construction of a high order Residual Distribution scheme for complex viscous flows

Keywords: Navier-Stokes equations, High order schemes, Residual distribution, Computational Fluid Dynamics

In this work we describe the use of the Residual Distribution schemes applied to the discretization of conservation laws. In particular, emphasis is put on the construction of a third order accurate scheme. We first recall the properties of a Residual Distribution scheme and we show how to construct a high order scheme for advection problems. Furthermore, we show how to speed up the convergence of the 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 numerical scheme is used to discretize the advection-diffusion scalar problem and the compressible Navier-Stokes equations.

Introduction

The high order discretization of the Navier-Stokes equation is still a challenging task, especially when complex flow regimes are considered like turbulent flows. The first step to obtain an accurate discretization of these problems is the construction of an accurate and robust solver for the laminar Navier-Stokes problems.

In the last years different high order schemes have been considered to obtain high order (more than two) solutions 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 represent an 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 with the cell polynomial degree and also the limiting strategy is not completely clear in the contest of the DG schemes. 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.

One of the first strategy considered [START_REF] Paillère | Multidimensional upwind residual distribution schemes for the convection-diffusion equation[END_REF] for the discretization of the Navier-Stokes equations was the coupling of the RD scheme for the advection with the standard Galerkin discretization for the diffusion, unfortunately the simply combination of the two schemes results in a scheme which is in general only first order accurate, unless advection or diffusion dominated problems are considered. An improved discretization consists of using a blending function to couple the advective and diffusive scheme, and has been use to obtain a third order discretization of the Navier-Stokes equations on triangular meshes [START_REF] Villedieu | Third order residual distribution schemes for the Navier-Stokes equations[END_REF]. Another approach treats the advection and the diffusion together with the same numerical scheme and has shown to obtain a third order discretization on grids composed by triangles (2D) and tetrahedrons (3D) [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF]. In a more recent approach [START_REF] Nishikawa | A first-order system approach for diffusion equation. II: Unification of advection and diffusion[END_REF] the advection-diffusion problem is rewritten as a first order, hyperbolic, system in which the unknowns are the variables and their derivatives. This approach has been already considered in the framework of the DG scheme [START_REF] Bassi | Discontinuous Galerkin solution of the reynolds-averaged Navier-Stokes and k -ω turbulence model equations[END_REF], where the use of discontinuous solutions allows one to update explicitly the gradients of the solution without considering additional equations. This approach, although promising, needs further investigations and must still be tested on complex test cases.

In this work we extend the RD scheme proposed in [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] for advection problems to the case of advectiondiffusion problems in order to construct a numerical scheme with an uniform order of accuracy for the Navier-Stokes equations. Furthermore the procedure used to imposed the boundary condition is detailed. The paper is organized as follows. In the section 2, we present the basic principle of the RD scheme, and the accuracy properties of the scheme for the scalar advection problem. In the section 3, we extend the RD scheme to the case of the advection-diffusion problem and in the section 4 the ideas introduced for the scalar problem are applied to the case of the system of equations and an efficient iterative procedure for the solution of non linear steady-state problem is presented. In the section 6 are presented the results of the discretization of the scalar advection-diffusion problem in term of the accuracy of the scheme. The proposed scheme is also used to obtain the solution of laminar flows in two and three spatial dimensions. In the last section, we give concluding remarks.

High order RD schemes for hyperbolic conservation laws

In this section we describe the RD scheme for the discretization of hyperbolic conservation laws in order to report the main properties of the scheme and introduce the nomenclature used through the paper. For simplicity a scalar problem is considered first.

The model equation, for a multidimensional problem, reads

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

u| ∂Ω -= g(s), s ∈ ∂Ω -,
where ∂Ω -= {x ∈ ∂Ω | n • ∇ u f < 0}, 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 ∂Ω -. The domain Ω is first discretized with N E non-overlapping elements with characteristic length h E . The set of all the elements is denoted with E h . In the RD scheme the DOF are associated with the points of the mesh and not with the control volumes as in the Finite Volume or DG methods. We denote by {σ l } l=1,N E dof the list of the DOF of the generic element E. When a linear interpolation of the solution is used, the DOF of each element coincide with the vertices. The higher order interpolation of the solution, necessary to construct an high order RD scheme, is obtained by adding extra DOF on each element. In order to keep a local formulation, the extra DOF are added inside the element. Since we use a continuous interpolation, all the DOF on the elements 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 increases more rapidly than in the discontinuous one and both cases become asymptotically similar.

To obtain a discrete counterpart of the Eq. ( 1), the total residual is first computed for each element E of E h as

Φ E = E ∇• f h (u h )dΩ = ∂E f h (u h ) • n d∂Ω, (2)
where u h and f h are approximations of u and f , respectively. Due to the numerical approximation of the solution, the quantity Φ E will be in general not null on each element. In order to handle only nodal values the total residual is distributed to each DOF of the element, the way in which this step is performed characterizes the behavior of the scheme. The quantities distributed to each DOF of the element are indicated for a generic scheme as {Φ E σ } σ=1,N E dof . It is easy to show [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF] that the following conservation constrain must be satisfied

σ∈E Φ E σ = Φ E , ∀ E ∈ E h , (3)
we assume that the residuals Φ E σ depend continuously on the values of {u σ } σ∈E . Imposing the following constrain

E∋σ Φ E σ = 0, ∀ σ ∈ E, ∀ E ∈ E h , (4)
we obtain a relation for each DOF. The numerical solution is obtained resolving the resulting non linear system.

If σ is a DOF belonging to ∂Ω -, the boundary conditions must be taken into account into the Eq. ( 4). We consider a numerical flux f ∂ , which depends on the boundary condition g, the outward normal n Γ to the boundary, and the local state u h . We define the residual on the boundary as Φ Γ σ , with the following conservation constrain

σ∈Γ Φ Γ σ = Γ f ∂ (u h , g; n Γ ) • n Γ -f h (u h ) • n Γ d∂Ω = Φ Γ σ , ∀ Γ ⊂ ∂Ω -, (5)
where Γ is an edge (or a face) on ∂Ω -. The Eq. ( 4) with the contribution of the boundary conditions becomes

E∋σ Φ E σ + Γ∋σ Φ Γ σ = 0, ∀ σ ∈ E, Γ and ∀ E ∈ E h , ∀ Γ ∈ ∂Ω -.
It can be shown [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 w, such that u h → w when h → 0, then w 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 [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].

Accuracy constraints

On the basis of the work done in [START_REF] Abgrall | High-order fluctuation schemes on triangular meshes[END_REF], are now introduced two important properties that must be considered in the construction of a high-order monotone scheme.

The first property, historically called linearity preserving, establishes the conditions under which the scheme is effectively high order accurate. Assuming u is smooth enough and u h is its P k interpolant, and writing the distributed residual as

Φ E σ = β E σ Φ E and Φ Γ σ = β E σ Φ Γ ,
the numerical scheme is accurate with order h k+1 providing that

Φ E σ = O(h k+d ) and Φ Γ σ = O(h k+d-1 ), (6)
and β E σ is uniformly bounded. The second property is related to the monotonicity of the numerical scheme. If we re-write the total residual on the element as

Φ E σ = σ ′ ∈E c E σσ ′ (u σ -u σ ′ ), (7)
using this formalism, the Eq. ( 4), neglecting the boundary conditions, becomes

E∋σ σ ′ ∈E c E σσ ′ (u σ -u σ ′ ) = 0.
In general, the coefficients c E σσ ′ 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 σ -δ σ E∋σ σ ′ ∈E c E σσ ′ (u σ -u σ ′ ) n , (8)
with δ σ a relaxation parameter. If the scheme satisfies the following positivity conditions

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

min E∋σ min σ ′ ∈E u 0 σ ′ ≤ u n σ ≤ max e∋σ max σ ′ ∈E u 0 σ ′ ,
with u 0 σ ′ the value of the initial solution at the DOF. A scheme that verifies the conditions ( 9) is said monotonicity preserving. These conditions do not imply that the iterative scheme ( 8) 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 unsteady problem[END_REF]. As a consequence, a monotonicity and linearity preserving scheme must be non linear.

A systematic way to construct a non linear scheme which is both linearity and monotonicity preserving has been shown in [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF]. The basic idea consists in distributing the total residual to the DOF of each element by the means of a low-order scheme, monotonicity preserving but not linearity preserving, for which the distribution coefficients

β E σ = Φ E Φ E σ ,
are not bounded. The high-order monotone scheme is constructed by applying a mapping β E σ → βE σ , such that the distributions coefficients βE σ of the high order scheme are bounded and do not violate the conservation and the monotonicity constrains. The high-order distributed residuals are computed as: ΦE σ = βE σ Φ E .

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 scheme1

Φ E σ = Φ E N E DOF + α E (u σ -ū), (10)
with ū = σ∈E uσ N E

DOF

and α E a parameter large enough to guaranty the stability of the scheme. Approximating u with its P k interpolant on the element E, u h = σ∈E u σ ψ σ , the residual reads

Φ E = E ∇ • f (u h ) dΩ = ∂E f (u h ) • n d∂Ω.
If we define now k E σ = E λ • ∇ψ σ dΩ, with λ = ∇ u f , it is easy to see that the Eq. ( 10) can be put in the form of the Eq. [START_REF] Abgrall | High order residual distribution scheme for Navier-Stokes equations[END_REF] with

c E σσ ′ = k E σ -α E N E DOF , and the condition c E σσ ′ ≥ 0 is satisfied if α E ≥ max σ∈E |k E σ |.
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 non-linear technique previously described, and the distributed high-order residual ΦE σ is obtained. 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 [START_REF] Abgrall | Essentially non-oscillatory residual distribution schemes for hyperbolic problems[END_REF] and the solution proposed consists in modifying the scheme as follows

Φ ⋆,E σ = ΦE σ + E (λ • ∇ψ σ ) τ λ • ∇u h dΩ, τ > 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 σ∈E Φ ⋆ σ = Φ E since σ∈E ∇ψ σ = 0. Experimentally, we can see that the non oscillatory properties of the scheme are not spoiled.

Extension to viscous terms

In this section we extend the RD scheme, developed previously for the advection problem, to the case of advection-diffusion problem. The scalar model equation for the advection-diffusion problems reads

∇• f (u) = ∇• (ν∇u), ∀ x ∈ Ω ⊂ R d , d = 2, 3 u = g(s), s ∈ ∂Ω - u = g 0 (s), s ∈ ∂Ω 0 (12)
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 the scheme for the advection equations with the Galerkin approximation of the viscous terms [START_REF] Paillère | Multidimensional upwind residual distribution schemes for the convection-diffusion equation[END_REF]. It has been shown [START_REF] Abgrall | Numerical approximation of parabolic problems by means of residual distribution schemes[END_REF] that a scheme resulting from the coupling of a RD scheme with a Galerkin discretization of the viscous term is still a residual method, but this is true only for P1 elements. Furthermore it is well know [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF] that the scheme obtained as the sum of the RD scheme and the Galerkin scheme is second order accurate in the diffusion and advection limits, but it is only first order accurate when advection and diffusion are equally important. One way to construct a scheme with an uniform order of accuracy has been proposed in [START_REF] Villedieu | Third order residual distribution schemes for the Navier-Stokes equations[END_REF], where the Galerkin scheme for the diffusion is added to the LDA scheme for the advection using a blending function.

A different approach [START_REF] Nishikawa | A first-order system approach for diffusion equation. I: Second-order residual-distribution scheme[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 advection-diffusion equation, Eq. ( 12), is equivalent to a hyperbolic relaxation system at the steady state. The reformulation of the original advectiondiffusion equation as an hyperbolic system, allows one to use all the techniques already known in the field of the hyperbolic system of conservation laws to obtain a solution of the advection-diffusion problem with an uniform order of accuracy on both the solution and its gradient.

We follow an alternative way to discretize the advection-diffusion equation in the framework of the RD scheme. In order to preserve the order of accuracy of the scheme, the total residual must be take into account both the advective and diffusive parts, and a single distribution scheme must be used to distribute the total residual. This idea has been introduced in [START_REF] Caraeni | Compact third-order multidimensional upwind scheme for Navier-Stokes simulations[END_REF], where the LDA scheme is used to distribute the total residual consisting of the advective and the diffusive parts. A third order accurate solution is obtained by considering the Hermite cube interpolation along the edge of the elements (triangles or tetrahedrons). This scheme has been also considered in [START_REF] Nishikawa | On high-order fluctuation-splitting schemes for Navier-Stokes equations[END_REF] where the distribution coefficients of the LDA scheme are modified in a such a way that they make the scheme isotropic in the diffusion limit.

Construction of the RD scheme for the advection-diffusion problem

The first step in the construction of a RD scheme is the computation of the total residual on each element of the mesh. As explained previously, the total residual must contains both the advective and diffusive parts, namely

Φ E = E ∇• f (u) -∇• (ν∇u) dΩ,
thanks to the divergence theorem, the previous equation can be recasted as

Φ E = ∂E f (u) -ν ∇u • n d∂Ω, ( 13
)
where ∇u is the reconstructed gradient of the solution. The use of the reconstructed gradient is necessary because the normal flux to each face of the element must be continuous and in general the components of gradients normal to the element's face are discontinuous across the face shared by two neighboring elements. Different options are available for the gradient reconstruction as, for example, the Green-Gauss formula or the Least-Square procedure. The total residual is distributed to each DOF of the elements according to the procedure described in the sub-section 2.2. The high order distributed residual for the advection-diffusion equation reads

Φ ⋆,E σ = ΦE σ + E λ • ∇ψ σ τ λ • ∇u h -∇ • (ν∇u) dΩ, τ > 0. ( 14
)
In the vanishing viscosity limit, the scheme reduces to the scheme ( 11) developed for the advection problems.

From numerical experiments it has been observed that the previous scheme it not satisfactory because it is characterized by a large numerical error and by a not uniform order of accuracy. A modification of the original scheme has been inspired by the work of Nishikawa [START_REF] Nishikawa | Robust and accurate viscous discretization via upwind scheme I: Basic principle[END_REF], where an accurate viscous discretization is obtained by rewriting the diffusion equation as a hyperbolic first order system which is discretized by a scheme suitable for advection problems. From the discretization of the first order scheme a discrete counterpart of the viscous term can be extracted and applied back to the original advection-diffusion equation. As detailed in the reference [START_REF] Nishikawa | Beyond interface gradient: A general principle for constructing diffusion scheme[END_REF], when the Lax-Wendroff scheme is applied to the diffusion scalar equation the following scheme is obtained for the case of P 1 elements

Φ E σ = Φ E 3 + νδ 4 ∇u E -∇u E • n E σ , (15)
where n E σ is the normal vector to the face opposed to the node σ scaled by the length of the face itself, and we set the parameter δ = 1. The term ∇u E -∇u E represents the difference, on the element, between the gradient of the solution and the reconstructed gradient. In this work the scheme [START_REF] Eisenstat | Choosing the forcing terms in an inexact newton method[END_REF] is written in a more general way as follows

Φ E σ = ΦE σ + E νδ 2 ∇u E -∇u E • ∇ψ E σ dΩ, (16)
where Φ E σ is, as usual, the distributed residual and the integral term vanish when summed over the DOF of the element, thus it does not contribute to the total residual and it acts only as a damping term in the viscous discretization scheme.

The schemes ( 14) and ( 16) can be combined together as follows

Φ ⋆,E σ = ΦE σ + ξ(Re) E λ • ∇ψ σ τ λ • ∇u h -∇ • (ν∇u) dΩ + 1 -ξ(Re) E νδ 2 ∇u E -∇u E • ∇ψ E σ dΩ, (17)
where Re E = λ hE ν is the Reynolds number on the element, with h E characteristic length of the element, and the function ξ(Re) is such that ξ(Re) → 0 in the diffusion limit (Re → 0) and ξ(Re) → 1 in the advection limit (Re → ∞).

Gradient reconstruction

The evaluation of the total residual for the advection-diffusion problem requires a gradient reconstruction in order to have a continuous normal flux at each face. Different techniques have been used, the easiest one is the Green-Gauss formula

∇u |σ = E∋σ ∇u E |σ |E| E∋σ |E| .
Another approach is the Least-Square technique, in which the gradient reconstruction is obtained by solving for the values of the gradients that minimize the following function

N k=1 E 2 ik , ∀i with E 2 ik = -∆u ik + ∂u ∂x i ∆x ik + ∂u ∂y i ∆y ik + ∂ 2 u ∂ 2 x i ∆x 2 ik + ∂ 2 u ∂ 2 y i ∆y 2 ik + ∂ 2 u ∂x∂x i ∆x ik ∆y ik 2 , (18)
with ∆u ik = u(x k , y k ) -u(x i , y k ), ∆x ik = x k -x i , ∆y ik = y k -y i , and where i is the node where the gradient is reconstructed and k = 1 . . . N are the neighboring nodes of i. The resulting linear system is overdetermined and must be solver be the means of a Least-Square procedure. In the case of P 1 elements only the first order derivatives are considered while in the case of P 2 elements also the second order derivatives are considered in the Eq. ( 18)

The last technique considered is the L2 projection in which the reconstructed gradient is obtained solving the following equivalence ∇u = ∇u in a weak sense

Ω ψ i ∇u dΩ = Ω ψ i ∇u dΩ, ∀i.
By expanding the gradient as ∇u ≃ k ψ k ∇u k , ones obtain the following linear system: M { ∇u} = b, with

M ij = Ω ψ i ψ k dΩ, and b i = Ω ψ i ∇u dΩ.

Discretization of system of conservation laws

We consider now the extension of the RD scheme, developed for the scalar advection equation, to the case of the system, namely

∇ • f a (u) = 0, (19)
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 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 , while the the flux function is

f a (u) =       m m ⊗ m ρ + P I (E t + P ) m ρ       ,
where P is the pressure and I ∈ R d×d is the identity matrix. The thermodynamic law of the polytropic ideal gas is used, namely

P = (γ -1) E t - m 2 2ρ .
The advection-diffusion problem for a system of equations with the vector unknown u can be written as follows ∂u ∂t

+ ∇ • f a (u) = ∇ • f v (u, ∇u), (20)
with f v (u, ∇u) the diffusion flux which generally depends on u and ∇u. In the case of Navier-Stokes equations u and f a (u) are the vector of the conservative variables and the advective flux function, respectively, as defined for the Euler equations, while f v (u, ∇u) is the viscous flux function, defines as

f v (u, ∇u) =      0 S S • m ρ -q      ,
where q = -κ∇T is the heat flux, with T the temperature and κ the thermal conductivity, and S is the viscous stress tensor

S = - 2 3 µ∇ • m ρ I + µ ∇ T m ρ + ∇ m ρ ,
with µ the fluid viscosity. It is well know that the viscous flux function f v is homogeneous with respect to the gradient of the conservative variables ∇u

f v (u, ∇u) = K(u)∇u,
with the homogeneity tensor K(u) = ∂f v (u, ∇u) ∂∇u . The discretization of the system of equations can be deduced in a straightforward manner from the corresponding scalar problem, in particular for the advection-diffusion problem the total residual on a generic element E reads

Φ E = e f a (u) -K(u) ∇u • n d∂Ω,
with ∇u the reconstructed gradient of the conservative variables. The total residual is first distributed to all the DOF of the element using the low order Rusanov scheme, from that a non-linear monotone scheme is constructed and the filtering term is added together with the dumping term acting for the viscous part. The schemes shown in the section 3 for the scalar problem applies in the same way to the system, namely

Φ ⋆,E σ = ΦE σ + ξ(Re) E A • ∇φ σ τ A • ∇u h -∇ • K∇u dΩ + 1 -ξ(Re) E 1 2 K ∇u E -∇u E • ∇ψ E σ dΩ.

Boundary conditions

The imposition of the boundary conditions, in the RD framework, has not reached yet a mature level of understanding, for that reason we shown here the way followed to impose the boundary conditions in the case of inviscid and viscous problems.

We propose a flexible way to distribute the residual computed with the corrections fluxes on the boundary elements, that is generalizable to all the spatial dimensions and does not depend on arbitrary parameters. Let us consider the weak form of the conservation law

Ω ψ∇ • f(u h ) dΩ = 0,
an integration by parts gives

∂Ω f(ũ) • n d∂Ω - Ω ∇ψ • f(u h ) dΩ = 0,
where ũ represent the state that take into account the boundary conditions. Replacing the basis function ψ with the standard Lagrangian basis functions ψ i , ∀ i ∈ Ω, and applying again the integration by parts to the second member on the left hand side of the previous equation one obtains

Ei ψ i • f(u h ), dΩ + ∂Ei∩∂Ω ψ i f(ũ) -f(u h ) • n d∂Ω = 0,
where E i is the set of the elements that have the node i in common. The boundary integral on the previous equation can be regarded as the contribution that must be added to the residuals computed on the domain elements without considering the boundary conditions. Note that the correction flux, f(ũ)f(u h ), becomes null as the boundary conditions are correctly enforced into the numerical scheme.

In the case of an a inviscid flow of a solid wall, the boundary conditions require that the normal component of the velocity to the wall is null: v • n = 0, with v the velocity vector. The flux function and the correction flux read

f(ũ wall ) =       0 P n x P n y P n z 0       , f(ũ wall ) -f(u h ) = -v n       ρ ρv x ρv y ρv z E t + P       , respectively, with v n = v • n.
The inflow/outflow boundary conditions can be imposed specifying the state ũin/out . As standard practice the flux function is linearized as

f(ũ in/out ) ≃ A + n (u h ) + A - n (ũ in/out ),
where

A ± n (u) = R(u, n) Λ ± (u, n) L(u, n), with R(u, n), L(u, n
) the matrices of right and left eigenvectors respectively and Λ(u, n) ± the diagonal matrix of the eigenvalues. The operator ± selects the positive and the negative eigenvalues and n is the outward normal versor to the boundary face. The correction flux for the imposition of the inflow/outflow boundary conditions reads

f(ũ in/out ) -f(u h ) = A - n (u h )(ũ in/out -u h ).
The same procedure used to impose the inflow/outflow boundary conditions can be used to specify the far-field boundary conditions, in this case usually the value of pressure at far-field boundary is specified. For a viscous fluid, at the solid surface we apply the adiabatic no-slip boundary condition, which require that the velocity and normal heat flux on the wall are zero: v| wall = 0, q • n| wall = 0. In the RD scheme the adiabatic, no-slip boundary conditions are generally applied by initializing the velocity field such that the velocity components are zero on the wall and enforcing the residual associated to the momentum equation to be also zero on the boundary nodes at each iteration, this corresponds to impose strongly the boundary conditions for the momentum equation. Nothings is explicitly done for the residual associated to the energy equation, this is partially justified by the fact that on the wall

∂Ei∩∂Ω (v • S -q) = - ∂Ei∩∂Ω q,
because the velocity is zero, so if the previous integral is not computed this is equivalent to assume that q is zero on the wall. However this consideration is not true in practice, and imposing the boundary condition in such a way does not guarantee that the heat flux is zero on the wall.

We propose an alternative way to impose the boundary condition for the energy equation in a weak sense, which uses the correction flux technique introduced in the inviscid case. Since the velocity is zero on the wall, the flux function in the direction normal to the wall becomes

f • n w =    0 S • n w -q • n w    ,
with n w the unit vector directed along the normal direction to the wall. In order to have zero heat flux on the wall, a correction flux can be considered. This is exactly the same way how slip boundary conditions are imposed in inviscid simulations. The correction flux on the wall is

f wall =     0 0 0 q • n     .
The components of the heat flux can be easily computed, with no extra effort, from the gradient of the conservative variables (already available from the computation). The heat flux for a perfect gas can be written as

q = -k∇T = -k γ -1 R ∇e,
since for a perfect gas T = γ-1 R e, where γ is the specif heats ratio, R is the gas constant and e is the internal energy: e = E t ρ -1 2 v 2 . Considering, for example the x component of a two dimensional flow, it is easy to see that

∂e ∂x = ∂ ∂x E t ρ - 1 2 v 2 = ∂ ∂x E t ρ - 1 2 ∂ ∂x m 2 x ρ 2 + ∂ ∂x m 2 y ρ 2 = 1 ρ ∂E t ∂x - E t ρ 2 ∂ρ ∂x - m x ρ 2 ∂m x ∂x - m 2 x ρ 3 ∂ρ ∂x + m y ρ 2 ∂m y ∂x - m 2 y ρ 3 ∂ρ ∂x =0,wall = 1 ρ ∂E t ∂x - E t ρ 2 ∂ρ ∂x ,
with m x = ρv x and m y = ρv y , the components of the momentum, which are zero on the wall. The heat flux as function of the conservative variables can be written as

q = - k ρ γ -1 R ∇E t - E t ρ ∇ρ
Now, it is possible to compute the residual on the wall nodes, as

∂Ei∩∂Ω ψ i f wall -f(u h ) • n d∂Ω.
As usual the residual associated to the momentum equation is explicitly imposed to be zero in order to apply the strong boundary conditions for the velocity on the wall.

Curved elements

The main advantage of high order methods over the low order ones is the capability of the high order methods to achieve an higher level of accuracy given the same computational resources. Obviously, since an high order method introduces more degrees of freedom than a lower one, a fair comparison between high and low order methods must be done with the same number of degrees of freedom. In other words the computational grids used by high order methods should be coarser than those used by low order methods. In generating coarse grids for an high order RD scheme, it is critical that curved boundaries are represented with high fidelity otherwise the error produced by a low-fidelity representation of a curved boundary may nullify the benefits of the high order discretization of the governing equations.

In this work the geometry of the high order elements is defined by the following isoparametric transfor-mation from the reference space, E, to the physical space, E

x(ξ) = N E DOF i=1 x i ψ i (ξ) ∀ ξ ∈ E, (21)
where ξ is the independent variable in the reference space and x i are the know coordinates of the nodes defined in the physical space. Note that the nodes of the element edges belonging to the boundary ∂Ω are places on the real geometry of the boundary. A curved boundary is therefore approximated by piecewise polynomial curved segments. The use of Lagrangian function guarantees the geometric continuity between neighboring elements, Figure 1.

The integral that defines the total residual is computed as the sum of the integrals on each face of the element

∂E f(u h ) • n d∂Ω = N E f f =1 ∂Ω f f(u h ) • n d∂Ω,
where N E f is the number of faces of the element E. The integral on each face is approximated with a quadrature formula

∂Ω f f(u h ) • n d∂Ω ≃ N q q=1 ω q f(u q ) • n q J q ,
where N q is the number of quadrature points, ω q are the quadrature weights, n q is the normal versor to the face at the quadrature point q and Jq is the determinant of Jacobian of the transformation from the reference to the physical face, evaluated at the quadrature point q. The quadrature formula must take into account that n, for non-linear elements is not constant on the faces, it is a polynomial one order lower than the transformation polynomials [START_REF] Marques | Comparison of matrix-free acceleration techniques in compressible Navier-Stokes calculations[END_REF].

Solution of the non linear system

The discrete counterpart of the system of equations [START_REF] Leicht | Error estimation and anisotropic mesh refinement for 3d laminar aerodynamic flow simulations[END_REF] or [START_REF] Luo | A fast, matrix-free implicit method for compressible flows on unstructured grids[END_REF] is obtained by assembling for each DOF σ the contributions from all the elements E ∋ σ plus the contribution from the boundary conditions, namely

E∋σ Φ E σ (u h ) + Γ∋σ Φ Γ σ (u h , ũ, n) = 0, ∀σ. (22)
The previous equations represent a system of non linear equations in the form

F (u) = 0, (23)
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 [START_REF] Nishikawa | Beyond interface gradient: A general principle for constructing diffusion scheme[END_REF] is

u k+1 = u k -J(u k ) -1 F (u k ), k = 0, 1, 2, . . . ( 24
)
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 ( 24) 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 ( 22) 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

E∋σ Φ E σ (u h ) + Γ∋σ Φ Γ σ (u h , ũ, n) = 0, -→ |C σ | ∂u h σ ∂t + E∋σ Φ E σ (u h ) + Γ∋σ Φ Γ σ (u h , ũ, n) = 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 ) ( 25
)
u n k+1 = u n k + ∆u n k , k = 0, 1, 2, . . . , (26)
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 ( 24) is retrieved. The evolution of the time step is controlled by the CFL number which is chosen according to the following law [START_REF] Tidiriri | Preconditioning techniques for the newton-krylov solution of compressible flows[END_REF] CFL

n = CFL n-1 F (u n-1 ) ∞ F (u n ) ∞ , (27)
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.

At each time step the linear system (25) must be solved which is recast 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. 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. 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 ) ( 28
)
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. Since in the GMRES algorithm the matrix A is only required in a matrix-vector product and remembering that the matrix contains the Jacobian of the numerical scheme, one can approximate the product of the matrix A with a generic vector w as follows [START_REF] Jackson | Nonlinearly preconditioned krylov subspace methods for discrete newton algorithms[END_REF][START_REF] Brown | Hybrid krylov methods for nonlinear systems of equations[END_REF][START_REF] Knoll | Jacobian-free newton method: a survey of approaches and application[END_REF]

Aw = |C| ∆t n I + J(u n k ) w ≃ |C| ∆t n w + F (u n + ǫw) -F (u n ) ǫ ,
with ǫ a small step size chosen as [START_REF] Pernice | Nitsol: A newton iterative solver for nonlinear systems[END_REF] ǫ

= 1 + u L2 z L2 ǫ rel , ǫ rel = 10 -8 .
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 matrixfree version.

The use a Jacobian-free technique without preconditioning usually leads to a lack of convergence on the iterative process. Among all the type of preconditioning developed for CFD applications, we consider in this work the Jacobi and the LU-SGS [START_REF] Jameson | Lower-upper implicit schemes with multiple grids for the Euler equation[END_REF] preconditioning. They represent a good compromise between robustness and memory requirement, and guarantee a good scalability in parallel algorithms. In the Jacobi preconditioning only the diagonal block of the approximated Jacobian are retained, while in the LU-SGS preconditioning the preconditioning matrix is taken as

P = (L + D)D -1 (D + U ),
where D, is the diagonal block of the approximated Jacobian, L and U are respectively the lower and the upper parts.

The Jacobian-free approach allows quadratic convergence of Newton's 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 iterations number [START_REF] Abgrall | High order residual distribution scheme for Navier-Stokes equations[END_REF].

Numerical results

In this section, the RD scheme shown in the previous sections, is used first to obtain the discretization of the simple scalar advection-diffusion problem for which an exact solution exists so that the numerical accuracy of the scheme can determined. Subsequently the scheme is also applied to the discretization of the Navier-Stokes equations. In particular two test cases are considered: the standard two dimensional laminar flow over a flat plate and then the more challenging three dimensional flow over a delta wing. All the simulations are considered to be converged when the initial residual is reduced by ten orders of magnitude.

Scalar advection-diffusion

To test the accuracy of the discretization of the advection-diffusion equations, the following two-dimensional scalar problem is solved with the RD scheme:

a ∇u = ν∇ • ∇u , in Ω = [0, 1] × [0, 1]
with a = a x , a y T and the viscosity, ν, constant. The exact solution 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 different values of ν are considered. The problem has been discretized by the schemes presented in section 3.1, namely Eq. ( 14) and Eq. ( 17), furthermore a modified scheme is also considered which is obtained imposing α = 0 in [START_REF] Jameson | Lower-upper implicit schemes with multiple grids for the Euler equation[END_REF], this corresponds to use a linear scheme.

We performed an accuracy study on a sequence of unstructured triangular grids and as error estimation the normalized L 2 norm of the difference between the exact and the numerical solution has been used

e L2(Ω) = Ω (u -u e ) 2 dΩ Ω u 2 e dΩ
, where u e is the exact solution. In all the simulations we initialize the solutions with a value zero everywhere except on inflow boundaries where the exact solution is imposed. From the numerical simulations it has been observed that the direct use of the scheme (14) generally does not guarantee that the iterative process reaches the steady state solution and that the scheme is not consistent, for this reason the results for this scheme are not shown.

In the first test case the viscosity coefficient is ν = 10 -6 such that an advection dominated problem is considered. In Figure 2-(a) is reported the L2 norm of the error as function of the parameter h = 1/ √ N dof , with N dof the total number of the DOF. It can be observed that the linear scheme is capable to preserve the right order of accuracy for both the second and the third order scheme. The non-linear scheme is second order accurate with the P 1 elements although it has a bigger error constant than the linear scheme. In the case of the P 2 elements the non-linear scheme loses the third order accuracy with the gird refinement, deteriorating to a first order only accurate scheme.

To investigate the behavior of the scheme when the advection and the diffusion are equally important we consider the advection-diffusion problem with ν = 0.01. The accuracy study is reported in Figure 2-(b). It must be observed that the schemes have at the most second order accuracy with the P 2 elements, and the linear scheme has lower error constant than the non-linear one. The second order schemes has the right order of accuracy. In Figure 3 is reported the accuracy study of the same advection-diffusion problem with ν = 0.01, as before, with the linear scheme and the special case in which the exact gradient of the solution is used instead of the reconstructed one. Since the exact solution of the problem is infinitely differentiable, the exact gradient is continuous. In this case the scheme has the right order of accuracy for both the second and third order scheme. Since the problem of loss of accuracy for the third order scheme seems related to the accuracy of the gradient reconstruction, we have performed an accuracy study of gradient reconstruction, on the same type of problem, for different techniques and two values of viscosity coefficient: ν = 10 -6 and ν = 0.01. The results are shown in Figure 4, it is worth noticing that the accuracy of gradient reconstruction does not depend on the value of the viscosity coefficient and that in the case of P 1 element the L2 projection technique has the smallest error while in the P 2 case the Green-Gauss formula and L2 projection has the same accuracy. The Least-Square technique is the less accurate. It is important to note that in the case P 1 the accuracy of the gradient reconstruction is much closer to the behavior of second order accurate than to the first order accurate technique, while in the case P 2 the accuracy of the gradient reconstruction is never more than two.

Navier-Stokes equations

We consider here the discretization of the laminar Navier-Stokes equations in order to verify the accuracy and the robustness of the proposed approach. The first test case that we consider is a laminar flow over a flat plate with a free stream Mach number, M = 0.3 a Reynolds number based on the free stream condition and the flat plate length Re = 5 000. The length of the plate is L = 1. The range of computational domain extends in the x-direction is [-1, 1] with the leading edge of the flat plate at x = 0. The size of the computational domain in the y-direction is 1, which is 10 times the boundary layer thickness at x = 1.0. At the inlet the inflow boundary condition in imposed, at the top and the exit the static pressure is imposed. Along the plane y = 0, the symmetry boundary condition is imposed for -1 ≤ x ≤ 0 and the wall no-slip boundary condition is imposed for 0 ≤ x ≤ 1. The coarsest mesh used is shown if Figure 5, if consists of 187 DOF. The medium and the fine grids are obtained successively dividing each triangle in four triangles.

We first consider the effect of the boundary conditions on the solution. In Figure 6 are compared the temperature and the axial velocity profiles at x = 1 for the second order schemes using the old boundary condition procedure and the new proposed approach on the medium grid. The results are compared against the exact solution of the thermal boundary layer [START_REF] Schlichting | Boundary-layer theory[END_REF]. Clearly the new approach is able to correctly impose the adiabatic condition, while the old approach completely fail to predict the temperature profile and also the heat flux is not zero on the wall. The velocity profile is instead in good agreement with the exact solution since the no-slip condition is strongly imposed in both the procedures. From the Figure 7-(a) can be also deduced that the new way to handle the adiabatic condition assure a fast convergence to the steady state solution. The third order scheme is much more sensible to the imposition of the boundary conditions as can be seen in Figure 7-(b), only the new way to impose the boundary condition on the wall makes the scheme converge. In Figure 8 is shown the axial velocity profile at x = 1 for the second and third order schemes on three different grids. The numerical results are in good agreement with the analytical solution even on the coarsest grid. If Figure 9 is shown the skin friction coefficient along the flat plate for the second and third order schemes on three different grids. Note that the agreements becomes better refining the grid and increasing the order of the scheme. The last test case we consider is the steady laminar three dimensional flow around a delta wing, the same test case has been considered in [START_REF] Leicht | Error estimation and anisotropic mesh refinement for 3d laminar aerodynamic flow simulations[END_REF]. The delta wing has sharp edges, the geometry of the wing is shown is Figure 10, together with the initial mesh which consist of 21 035 nodes (110 629 elements). The free-stream Mach number is 0.3, the Reynolds number is 4 000 and the angle of attack is 12.5 • . As the flow reaches the wing it rolls up creating a main vortex together with a secondary small vortex over the wing. The resulting vortex system remains over long distance behind the wing. Considering the complex flow patter, the numerical solution of the problem has been combined with a sequence of mesh adaptations. The grids are adapted using Mmg3d library [START_REF] Dobrzynski | Mmg3d: User Guide[END_REF] with the vorticity magnitude as criterion for the adaptation. In Figure 11 are compared the contours of Mach number and of the vorticity magnitude together with the mesh, at the end of the wing, for the initial grid and the adapted grid after one cycle. In Figure 12 are compared the contours of the Mach number and of the vorticity magnitude together with the mesh, at the end of the wing and on the wake, for the adapted grid after one cycle and two cycles of adaptation. The final grid consist of 102 413 nodes (594 595 elements). The grid adaptation is able to refine the mesh in the zone that contains the main vortex, resulting in a better resolution of the flow patter behind the wing. This can be observed in Figure 13, where are shown the slices of the Mach number and of the vorticity magnitude for the original grid and the adapted grids, at different positions along the wing and the wake. It is worth noticing how adapting mesh allows to capture better the vortex structure behind the wing, as it can be also seen in Figure 14 where the Mach iso-surface for the original and final adapted grid is shown.

On the final grid we performed also a third order simulation, the results of the second and third order scheme are compared in Figure 6.2. Finally, in Table 6.2 are reported the values of the aerodynamic coefficients for the different simulations. The values obtained on the final grid are not far from the reference values, the difference is probably due to the fact that during the adaptation process the surface discretization of the wing has not be modified respect to the original coarse grid, since the library used is able to adapt only on the volume. Nevertheless, it is worth noticing that the value of the drag coefficient is very close to the reference value, because in this case the main contribution to the drag is due to the flow separation over and behind the wing, and this feature is well captured by the grid adaptation.

C L C D
Original Grid (Order 2) 0.32737699 0.144857313 Adapted grid 1 (Order 2) 0.33456865 0.145790923 Adapted grid 2 (Order 2) 0.33923147 0.147307554 Adapted grid 2 (Order 3) 0.33667129 0.161839115 Reference [START_REF] Leicht | Error estimation and anisotropic mesh refinement for 3d laminar aerodynamic flow simulations[END_REF] 0.347 0.165

Conclusion

A high order accurate Residual Distribution scheme for the solution of the compressible Navier-Stokes equations has been presented. The scheme avoid the use of the combined discretization with the RD scheme for the advection and the Galerkin scheme for the diffusion. The contribution of the advection and of the diffusion are handled with the same scheme. It has been also proposed a Jacobian-free technique to accelerate the convergence of the implicit iterative scheme to the steady state solution. The proposed scheme, in combination with an improvement of the adiabatic no-slip boundary conditions, has shown to perform the standard test cases used for laminar compressible flows as well as complex three dimensional aerodynamic flows.

Figure 1 :

 1 Figure 1: Visual representation of the mapping from the reference element the physical element with curved faces. On the boundary the normal vectors to the face at the quadrature points are also shown.

Figure 2 :

 2 Figure 2: Normalized L2 solution error for the advection-diffusion scalar problem with ν = 10 -6 (left) and ν = 0.01 (right).

3 Figure 3 : 1 Figure 4 :

 3314 Figure 3: Normalized L2 error for the advection-diffusion scalar problem with ν = 0.01 in which the exact gradient is used instead of the reconstructed one.

Figure 5 :Figure 6 :

 56 Figure 5: The coarsest mesh used for the flat plate boundary layer case.

Figure 7 :

 7 Figure 7: L2 iterative residual of axial momentum equation for the old and the improved way to impose the adiabatic boundary conditions for the second (left) and third (right) order scheme.

Figure 8 :

 8 Figure 8: Axial velocity profile for the second (left) and third (right) order schemes at x = 1 as function of the dimensionless wall distance (η) on three different grids. The numerical solutions are compared against the theoretical solution.

Figure 9 :

 9 Figure 9: Skin friction coefficient along the flat plate fort he second (left) and third (right) order schemes on three different grids. The numerical solutions are compared against the theoretical solution.

Figure 10 :

 10 Figure 10: Geometry of the delta wing. The following parameter has been used: Λ = 75 • , σ = 60 • and t/c = 0.024. The initial coarse mesh is also shown.

Figure 11 :

 11 Figure 11: Contours of the Mach number (left) and of the vorticity magnitude (right) at the end of the wing. In each figure are shown the solution and the grid before and after the first grid adaptation. Second order scheme.

Figure 12 :Figure 13 :

 1213 Figure 12: Contours of the Mach number and of the vorticity magnitude at the end of the wing (first row) and behind the wing (second row). In each figure are shown the solution and the grid before and after the second grid adaptation. Second order scheme.

Figure 14 :Figure 15 :

 1415 Figure 14: Mach iso-surface; on the left part of the wing solution on the original grid on the right part of the wing solution on the final adapted grid. Second order scheme.

Table 1 :

 1 Values of the lift and drag coefficients for the different simulations of the delta wing test case.

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


Benoit Combemale

Xavier Thirioux

Benoit Baudry

Formally Defining and Iterating Infinite Models

The wide adoption of MDE raises new situations where we need to manipulate very large models or even infinite model streams gathered at runtime. These new uses cases for MDE raise challenges that had been unforeseen by the time standard modeling framework were designed. This paper proposes a formal definition of an infinite model, as well as a formal framework to reason on queries over infinite models. This formal query definition aims at supporting the design and verification of operations that manipulate infinite models. First, we precisely identify the MOF parts which must be refined to support infinite structure. Then, we provide a formal coinductive definition dealing with unbounded and potentially infinite graph-based structure.

Introduction

The growing adoption of Model-Driven Engineering (MDE) at all steps of software development comes with new requirements for MDE theories and tools. In particular, this work focuses on the need to process (i) finite but very large models, and (ii) infinite models. A major challenge to process these categories of models consists in understanding the exact meaning of a query over a model for which the interpretation does not know the size at a given point in time.

To illustrate the need to process finite but very large models, let us consider the complete model representing the entire Eclipse platform (the minimal workbench with OSGi). This model includes about 5 million model elements. Current model processing tools require all the model elements in memory (e.g., Eclipse Modeling Framework (EMF) [START_REF] Steinberg | EMF: Eclipse Modeling Framework (2nd Edition)[END_REF]). With EMF, the model of the Eclipse platform Java code requires 900MB in RAM memory. Steel et al. [START_REF] Steel | Model interoperability in building information modelling[END_REF] provide an even bigger example: when they adopted a MDE approach to analyze civil engineering models, they had to deal with more than 7.3 million computational objects. Programming languages provide a good source of inspiration to deal with these issues. Through the notion of lazy evaluation, programming languages allow (lazy) iterations on potentially infinite data-structures. Even in Java, which does not support laziness natively, skilled programmers tend to manually postpone object instantiations as much as possible, i.e. only when needed, in order to reduce instantaneous memory consumption. Recent work were inspired by this approach to propose lazy model transformations to process very large models [START_REF] Tisi | Lazy execution of model-to-model transformations[END_REF], or NoSQLbased approaches for model persistence [START_REF] Espinazo-Pagán | Morsa: A scalable approach for persisting and accessing large models[END_REF].

Beyond the problem of very large models, stands the issue of processing infinite models. This requirement becomes more and more critical with the growing adoption of the models@runtime paradigm [START_REF] Blair | Models@ run.time[END_REF]. For instance, if we consider a monitoring autonomous system which relies on a model at runtime to abstract information from a complex event processing (the CEP) engine. The CEP will indefinitely provide information about the environment and thus one cannot consider that the model at runtime will be bounded as it could grow indefinitely. Another illustration can be found in the realm of reactive systems, which are modeled by transition systems intended to run forever. In this case, the model at runtime that records the trace of states and events triggered during execution is another form of infinite model. To deal with infinite models, we could leverage mechanisms established in the area of web feeds, such as RSS syndication. In these cases, data following a predefined format is timely and infinitely appended to an initial model (even though the RSS file usually corresponds to a sliding window because older elements are removed). However, as far as we know, no model processing solution adopts the notion of sliding window over an infinite flow of model elements in order to deal with infinite models.

If we look at the current state of MDE theories and tools to deal with large and infinite models, we make two observations. First, metamodeling formalisms and most of the tools are deeply rooted on the assumption that models include a bounded number of model elements and that this bound is known when computing a query over the model. Second, there exist some ad-hoc implementations to deal with these issues, but there is no formal definition of infinite models and no reference formal semantics for a query over an infinite model. The consequence of these two observations is that it is currently impossible to formally verify operations that process very large or infinite models. This is a major challenge for the adoption of MDE in the use cases discussed above.

In this paper, we tackle the two limitations listed above through two major contributions.

a detailed analysis of current metamodeling standards and a precise identification on how they prevent the definition of infinite models. In order to face this, we propose an extension to the MOF formalism to enable the local definition of an infinite part of a model. a coinductive semantics for a query operation over an infinite part of a model. This semantics relies on a formal definition of infinite models, and provides both a reference for various implementations and the foundations for the verification of operations that must process models of unknown size.

The paper is organized as follows. Section 2 illustrates through a concrete example how the current metamodeling formalisms such as MOF prevent the definition and manipulation of infinite models. Based on this observation, Section 3 introduces MOF extensions supported by a formal definition of infinite models, and Section 4 proposes a coinductive semantics supporting the manipulation of such infinite models. Since the proposed formal semantics is independent of any implementation choice, we discuss in Section 5 the various existing and possible implementations of coinductive operators. Finally, we conclude and outline our perspectives in Section 6.

Illustrative Example: How MOF Does Not Support Infinite Models

Model Driven Engineering (MDE) considers software artifacts as abstract typed graphs (i.e., models conforming to precisely defined metamodels). As discussed in the introduction, we have to deal with increasingly large models. In many cases these models may even be considered of unbounded and infinite size (i.e., their size is a priori unknown). Since models are conforming to metamodels, such situations must be considered in the definition of metamodels. These metamodels are themselves implemented using a meta-language, usually compliant with the Meta Object Facility (MOF) [START_REF]Meta Object Facility (MOF) 2.4.1 Core Specification[END_REF], such as Ecore [START_REF] Steinberg | EMF: Eclipse Modeling Framework (2nd Edition)[END_REF]. This section illustrates how a meta-language such as MOF ties MDE practices to a vision of finite and bounded models and thus prevents the definition and manipulation of infinite models. We illustrate these issues with the UML2 State Machine formalism.

UML2 State Machine as an illustrative example

The state machine sketched in the bottom left corner of Figure 1 conforms to the State Machine metamodel displayed in the middle left of Figure 1 (see the metamodel level). This metamodel defines a StateMachine as composed of several State elements, including an initial state, as well as several Event elements which it may react to. States are pairwise linked through Transition elements as source and target states. Each transition is triggered by a set of events (Trigger) and in return sends events (SendEvent) as the e f f ect of its firing.

The execution semantics of such a state machine processes as follows: The first RuntimeEvent to be processed (that is, in our case, an In jectEvent element) is popped from the eventToProcess queue belonging to the state machine. This element represents an EventOccurrence of some event. A RuntimeEvent is either locally raised by the state machine transitions (endogenous) or brought by the environment (exogenous) and several other preceding events, which have been previously processed, constitute the cause for which it occurs. Runtime events keep on being popped and put aside as unhandled, until this set contains enough events, for some outgoing transitions of the current state to be triggered. Then, actually firing the transition pushes some new events at the end of the event queue, changes the current state of the machine and removes the triggering events from the unhandled ones.

Since a model conforms to a precisely defined metamodel, the underlying model of a state machine (graph of objects) follows the constraints expressed in the metamodel, according to the MOF semantics.

Bounded collections of properties

As illustrated in Figure 1, the upper bound of the collection eventToProcess is relative to its cardinality. In MOF, this cardinality is reified in terms of the lower and upper attributes of the Property construct (cf. top of Fig. 1). According to the OMG MOF Specification, the value of the upper attribute is typed by UnlimitedNatural Consequently, the execution semantics of the iterators (e.g., the ones coming from OCL) on the collection eventToProcess is bounded. This means that the iterators assume that all elements are considered as available at any time of the iteration.

Transitive closure

The underlying cyclic directed graph structure of any MOF-based metamodel (cf. top of Fig. 1) raises the issue of evaluating the transitive closure of a cycle. Such an issue can be shown in Figure 1, with the execution semantics of the closure over the next states (obtained from a given state with outgoing → collect(target)). Since OCL 2.3, the standard includes a closure operation [8, §7.7.5]. This is very useful to specify recursive OCL operations. For instance, in Figure 1, the reachable states from a given state can be specified as proposed in the following OCL expression.

context State :: reachableStates () : Set( State ) body : { self }->closure(outgoing->collect( target ));

As stated by the OCL specification for the operator closure, "the collection type of the result collection is the unique form (Set or OrderedSet) of the original source collection.

If the source collection is ordered, the result is in depth first preorder.". Here again, the underlying semantics refers to the type Collection. Consequently, the execution semantics of the closure is a finite processing, which assumes that the whole model is available for evaluation.

Discussion

Iterating both over the collection eventToProcess or the corresponding closure of the reachable states is thus a finite process for which the whole model is required (e.g., a state machine defined at design time to model the behavior of a given class). The collection eventToProcess is bounded by the semantics of the attribute upper of a MOF property, and sets the "width" (i.e., number of outgoing edges from the same node) of the underlying graph of a conforming model. The reachable states process is finite, as defined by the underlying unfolding semantics as considered in the OCL operator closure, and sets the "depth" (i.e., length of a path with unique nodes) of the underlying graph of a conforming model. These are strong limitations imposed by current metamodeling formalisms since, as illustrated in the introduction, such a state machine should also be considered as locally infinite (e.g., state machine continuously updated to monitor at runtime a running and non-terminating program). In the following sections, we introduce slight modifications in MOF, which broaden its scope. In the context of infinite models, these modifications support the definition of a formal semantics for the MOF attribute upper and the unfolding semantics as used in the OCL operator closure.

Defining Infinite Models

This section starts with two proposals to extend the scope of MOF and allow the identification in metamodels of the infinite parts of the conforming models (i.e., parts which need to be manipulated despite their unknown size). These extensions are then used to provide a formal definition of infinite models.

Intuitive presentation

While the semantics described by the OMG in the MOF specification involves a finite interpretation of models, some situations require an infinite interpretation of the same structure. For instance, a state machine can be considered as infinite (cf. right part of Fig. 1) if it abstracts the execution trace of a non-terminating program. In practice, this execution trace can be lazily built at design time while exploring the graph of reachable states, or continuously built at runtime during the system execution (e.g., monitoring).

As seen in the previous section and illustrated in the example depicted in Figure 1, a model can be infinite in two situations, respectively in the width and in the depth of the underlying graph. In the following, we come back into these two situations and we propose MOF extensions with a concrete syntax to locally characterize in a metamodel the infinite parts of the conforming models.

-The collection defined in Figure 1 by the relation eventToProcess on StateMachine may be considered as infinite in case of a non terminating execution (i.e., an infinite sequence of runtime event). We are noting ω the upper bound of the multiplicity of an infinite collection (cf. right part of Fig. 1), compared to the * value which defines an unbounded but finite collection (cf. left part of Fig. 1). We assume that an infinite collection is ordered and then countable. -A reflexive relation may be indefinitely unfolded and then, the computation of the closure may not terminate (currently, common practices consider that the closure computation terminates). This situation may even be mandatory if we consider in this relation a multiplicity with a lower bound greater than zero. We propose to graphically note such reflexive relations with infinite unfoldings as an arrow with two heads (cf. in Fig. 1, right) 4 . We are aware that infinite unfolding may come from more complex cycles in a graph-based metamodel. For example in Figure 1, we may consider the cycle between State and Transition as an infinite unfolding. Common textual and graphical metamodeling notations cannot easily characterize a cycle. Nevertheless, a reflexive relation may be derived from an OCL expression characterizing the cycle. For example, the reflexive reference nextStates on State in Figure 1 is derived as specified by the following OCL expression. This derived reference characterizes the cycle between State and Transition.

context State :: nextStates : Set( State ) derive : self . outgoing->collect( target );

According to this new notation, the classical example proposed in the left of Figure 1 is modified as shown in the right of Figure 1 in order to locally consider infinite structures (in our case to consider a possible infinite execution of a state machine). Note that several collections and cycles voluntarily keep the initial semantics based on a finite part of the model. For example, the collection of states (resp. transitions) which compose a state machine is unbounded but remains finite, and the transitive closure of the reflexive relation cause always terminate. So this syntax allows to clearly define in a metamodel the parts of a model which should be interpreted as infinite.

Formal presentation

We propose a formal definition relying on the previous intuitive presentation of our extended metamodeling facilities.

In the following formal definition, we assume that we have finite sets of metaelements (MetaElements) and relations (Relations), i.e., a finite metamodel. We are also using Elements as the set of possible model elements without any type information.

Definition 1 (Infinite Model). Let M E ⊆ MetaElements be a bounded set of metaelements. Let R ⊆ { me 1 , r, me 2 | me 1 , me 2 ∈ M E , r ∈Relations} be the bounded set of relations among meta-elements such that ∀me 1 ∈ M E , ∀r ∈Relations,

card{me 2 | me 1 , r, me 2 ∈ R} ≤ 1.
We define an infinite model E, L ∈ Model(M E , R) as a multigraph built over an unbounded set E of typed elements and an unbounded set L of labeled links such that:

E ⊆ { e, me | e ∈ Elements, me ∈ M E } L ⊆ e 1 , me 1 , r, e 2 , me 2 e 1 , me 1 , e 2 , me 2 ∈ E, me 1 , r, me 2 ∈ R
We define the auxiliary type Natural ω = N ∪ { * , ω}. Natural ω is an extension of the UnlimitedNatural type provided by the OMG MOF specification. It is used to represent a range of possible numbers of instances. Unbounded finite ranges can be modeled using the * value whereas unbounded infinite ranges can be modeled using the ω value. The type Natural ω also comes equipped with the following order: m < * < ω, for all m ∈ N.

We aim at characterizing the presence of infinity both at the level of collections (i.e. the width of the underlying graph) and reflexive relations unfolding (i.e. the depth of the underlying graph).

First, regarding the width of the graph, we define the upper property which aims at distinguishing finite collections from infinite ones. Either for attributes or references (i.e., relation), a maximum number of instances of a target concept can be defined using the upper attribute, which value n is reflected in the following definition. Whether the upper bound of a relation is finite or infinite impacts the semantics (and implementation as well) of the model elements fetching operator get (cf. Section 4.2).

Definition 2 (Upper). The upper property characterizes an upper bound n of a multiplicity of a given relation, this bound been taken from Natural ω .

upper( me 1 , r, me 2 ∈ R, n ∈ Natural ω ) E, L → ∀ e, me ∈ E, me = me 1 ⇒ card({m 2 ∈ E | e, me 1 , r, m 2 ∈ L}) ≤ n
where the card function returns either m ∈ N or ω.

Second, regarding the depth of the graph, we introduce the property ru_unstable applying on primitive as well as derived relations. It requires the definitions of (maximal) model paths as the results of relation unfolding. The property ru_unstable states that a given relation only gives rise to finite unfoldings, whatever the maximal model path considered. We need to focus on maximal paths as we express properties about possibly infinite unfoldings. Whether a given relation of a model satisfies the ru_unstable property or not impacts the semantics and implementation of the OCL closure and other iteration operators applied to this relation (cf. Section 4.3).

Definition 5 (Unstable Reflexive Unfolding).

Considering model paths as creation traces (cf. definition 3), a relation has only finite unfoldings if and only if it is unstable in any maximal model path π. This condition is rephrased as the following Linear Temporal Logic (LTL) property: ♦¬r or equivalently ¬♦ r. Relying on our model definition, it amounts to directly defining the following property: ru_unstable( me, r, me ∈ R) E, L → ∀π ∈ maximal paths( E, L )

I π = N ⇒ ∃i ∈ I π , me i = me ∨ r i = r
where I π is the set of indexes of π and satisfies definition 3.

A Coinductive Semantics to Iterate Infinite Models

We discuss in this section the ways to manage model elements in the context of an infinite model (Subsection 4.1) and we propose a formal definition of the operators needed to querying such models. First we formalize the common operators for getting model elements (Subsection 4.2), and then we rely on them to formalize an alternative of the main iterators inspired from the OCL language (Subsection 4.3). We finally put into practice the proposed operators, among others for the manipulation of infinite state machines as defined in the right of Figure 1 (Subsection 4.4).

Reasoning on model elements: from finite to infinite model

Standard inductive semantics aims at defining finite data-structures as well as reasoning and programming with them. Therefore, an inductive structure comes naturally equipped first with an induction principle, allowing proofs by induction on the elements of this data-structure ; and second with a generic reduce programming primitive, allowing terminating recursive traversal of these data-structures. Induction principles are commonplace for reasoning about terminating algorithms and finite data-structures. It comes in many flavours such as induction on natural numbers, lists, binary trees, etc. As for the reduce and alike operators, they are also pervasive in programming paradigms, mostly in functional languages, but also in Java (e.g. the Iterator interface) and OCL (the iterate construction that operates on a finite collection).

Dually, coinductive semantics aims at defining potentially infinite data-structure. Therefore, a coinductive data-structure comes equipped with a coinduction principle and also a produce operator. The coinduction principle, among various usages, is at work when typing compilation units in languages supporting separate compilation. In Java for instance, you can type-check a bunch of classes, even if they are totally abstract and don't contain any piece of code. In this respect, checking type safety of two mutually dependent classes A and B works as if you were producing an infinite proof under the form: A is type-safe if B is type-safe if A is type safe, etc. There, type-safety is only proved to be a stable relation, with no base case at all. On the contrary, type-checking a concrete method amounts to reason inductively on the code structure, with assumed well-typed parameters and local object creations as base cases.

The produce operator (dual of the reduce operator) aims at producing potentially infinite data (as it cannot obviously perform a terminating recursive full traversal of an infinite structure), through the repeated execution of a piece of code that generates new values each time, appended to the growing structure. The resulting structure is the limit of this maybe infinite creation process, much akin to the fractal structures resulting from infinite iterations of a subdivision process. Often enough, as it is the case for instance in stream processing languages, the output structure is produced by a piece of code that consumes/explores in turn another corecursive input structure, one element at a time.

These coinductive concepts are mandatory in order to define the formal semantics of our MOF extensions (as used in the right part of Figure 1) independently of a particular implementation. The infinite state sequences of a state machine (cf. right lower part) should follow the semantics of the metamodel, just as plain finite models follow the inductive semantics of the metamodel. These state sequences may be defined as being produced from infinite sequences of transition-enabling events, following the execution semantics of the state machine. A standard inductive viewpoint on these sequences would be for instance to define an allInstances() OCL constraint that checks that each n + 1 th state is the result of executing a transition from a n th state. As each OCL operator is supposed to work on a collection of states taken as a whole, evaluating the constraint on an infinite model would yield a non-terminating behavior and no outcome at all. Moreover, how such a collection may be produced still remains an open question in this case.

As can be seen, coinductive semantics, which amounts to producing infinite proofs and data, may be found in various areas, even though not always presented as such. Defining and formalizing a (coinductive) structure of models will help at elaborating important and practical tools for infinite model manipulation. A coinductive semantics may be implemented in several ways: in a programming language with lazy constructs that will evaluate only the finite browsed part of the model, or with data stream primitives randomly producing new model elements consumed afterwards by the execution semantics when possible, or else with a prefetch semantics that estimates how many model elements should be evaluated in advance, even if not needed. Each such imple-mentation is interesting in its own right as it corresponds to a well-known class of applications. In the remainder of this section, we propose a formalization of a coinductive model semantics and the implementation standpoint is discussed in the next section.

Getting model elements

On a model m = E, L , we first assume the operator getRoots() which corresponds to a minimal set of model elements from which any other model element can be accessed (i.e., a covering set). The accessibility predicate is defined as the existence of a finite model path from a model element to another. A set of model elements is defined as minimal when it is a covering set such that no proper subset is also covering.

accessibility(e, e ∈ E) ∃{ e i , r i , e i+1 } i∈I ∈ model path( E, L ), e 0 = e ∧ ∃i ∈ I, e i = e covering(S ⊆ E)

∀e ∈ E, ∃e ∈ S, accessibility(e, e ) minimality(S ⊆ E) covering(S) ∧ ∀e, e ∈ S, ¬accessibility(e, e )

We can now specify the getRoots() operator which corresponds to the entry point facility on the model:

m.getRoots() ∈ {S ⊆ E | minimality(S)}
We assume that the getRoots() operator return a finite set of roots. We note also that an alternative model definition based upon rooted multigraph may be adopted and would yield a non under-specified definition. In this case, the set of roots is uniquely defined.

Relying on the entry point previously defined by the getRoots() operator, we mainly consider the get() operator for getting model elements from a model E, L . Usually, this operator allows to access to a property r from a model element e (written e.r in OCL).

e.get(r ∈ Relations) {e ∈ E | e, r, e ∈ L} where e ∈ E. In our formalization, we assume that the return value of the get() operator is always a collection of model elements, tagged as either finite or infinite in our MOF extension and processed accordingly with the appropriate iterator. We are aware that in most model management APIs (e.g., EMF) or in the OCL query language, the return value may exhibit different types according to the known multiplicity of the relation.

Iterating model elements

We propose in this section a formalization of an alternative version -called coiterateof the main generic OCL iterator iterate [8, § 7.6.5], from which all other iterators may be defined. Both iterators allow to browse a collection returned by the get() operator and to unfold a reflexive relation, for instance in order to compute its closure. A particularity of OCL finite collections is that they are implicitly ordered. Indeed, the n th element may easily be retrieved with the help of the iterate operator. It may merely appear as a design clumsiness in the API for collections in the finite case. Yet, the ordering of elements is mandatory in the infinite case, as elements will be processed sequentially, one by one, and the user may observe intermediate results of this infinite computation. So we assume that collections, whether finite or infinite, are ordered.

The iterate operator processes the successive values of an (ordered) finite collection in order to compute the final value of an accumulator of any type. Dually to this semantics, the coiterate operator starts from an initial value which processing produces a potentially infinite collection of new values. This approach, which applies the coinduction principles, allows to produce a new collection either from an already existing infinite collection or more generally as the sequence of values built from successive assignments of a variable of any type.

In order not to depart too much from the iterate operator, from a syntactic viewpoint, we define coiterate as follows, first recalling the definition of iterate. We require that coll:Collection(A) possesses the three basic operations provided on Collection by OCL: isEmpty() which tests a collection for emptiness, first() which returns the first element of a collection, and append(elem:A) which appends a new element elem at the end of the collection. For the sake of readability, we also define the operation tail() which returns the collection without the first element 5 . Note that, as e 1 , e 2 and e 3 are expressions and not values in the following definitions, we must use substitutions 6 in order to define recursively (resp. corecursively) these iterators. The coiterate operator starts from a finite collection coll, to which it will append a potentially infinite sequence of elements. For that purpose, it considers a variable acc, initialized with value e 1 . From the current value of acc, if not equal to null, a new element elem with value e 3 built from acc must be appended to the current resulting collection, and the next value of acc is given by e 2 . This process is repeated as long as acc is not null. So the resulting collection is finite if and only if acc finally becomes null, otherwise the collection is infinite. Moreover, both iterators are able to handle several accumulating variables of any name at once, provided the variable elem is defined.

Putting the coiterate iterator into practice

We illustrate how the coiterate iterator may be used, through the following two basic examples. In the first one, the infinite collection of the natural numbers is built and in the second one, the infinite collection of the squares of even numbers is built from the first collection. Similarly, the operator coiterate may be used to browse the infinite collection of events to process (cf. eventToProcess in the right part of Fig. 1) and then to specify the simulation of a state machine. For instance, the following listing define the body of an operation simulate() on StateMachine which creates the trace (i.e., a sequence of states) according to the events to process 7 . Our definition of the iterate (resp. coiterate) iterator was shown to browse a collection (resp. create a collection). However, it seems to be an elegant way to use these operators in order to also (co)iterate over unfoldings of a reflexive relation. In the following listing, the coclosure operator of a reflexive relation is defined in terms of coiterate and corresponds to a breadth-first traversal of the underlying graph. From this generic coclosure operator we also specify an operation on State computing the potentially infinite set of reachable states from a given state: We have proposed in the previous section a formal specification to define and evaluate infinite models relying on coinduction principles. Stemming from this first milestone, our immediate next goal will be to explore pragmatical solutions to implement such principles. In all cases, and by definition, it is not possible to store the entire infinite model in memory. When a causal dependency exists between the model producer and the model consumer, the model can be lazily interpreted (from the producer) in order to build, on demand (by the consumer), only the necessary model elements. In our example, the state machine may be lazily executed during a step-by-step simulation. Such lazily interpreted constructs have been recently proposed at the metamodel level in ATL [START_REF] Tisi | Lazy execution of model-to-model transformations[END_REF]. In this case, infinite model processing is close to existing lazy interpretation (also called call-by-need) of a data structure [START_REF] Henderson | A Lazy Evaluator[END_REF] and relies on the following properties :

the model elements in the infinite model are only built when it is necessary, once an infinite model element has been built, it is never built twice. However, in some situations it may be necessary to relax these properties. For example, the last property implies memoization8 and may be relaxed in some cases because it may render certain things impossible (e.g., because more elements may have to be remembered than system memory permits) or very inefficient (e.g., because of synchronization issues in concurrent systems). A non-memoized lazy evaluation also called call-by-name [START_REF] Douence | A systematic study of functional language implementations[END_REF] may then be of interest. For example, such strategies have been explored in [START_REF] Espinazo-Pagán | Morsa: A scalable approach for persisting and accessing large models[END_REF] as part of an approach for model persistence based on NoSQL databases.

Even if lazy interpretation seems promising, specifying the level of laziness is also important. In the context of the MDE, different solutions are possible. For instance, does model element creation entails creation of its references atomically? Or are these references also created on demand?

When the model consumption is disconnected from the model production, the order in which model elements arrive is out of control. Such an infinite model relies on the following properties (inherited from data stream [START_REF] Babcock | Models and issues in data stream systems[END_REF]):

model elements in the infinite model arrive online, once an element from an infinite model has been processed it is discarded or archived. Different semantics may be considered depending on several options. For instance, model elements delivery must take into account the following options from [START_REF] Franklin | A framework for scalable dissemination-based systems[END_REF]: Pull vs. Push, Aperiodic vs Periodic and Unicast Vs. 1-to-N. Then, mechanisms for model elements delivery may be inspired from known protocols in data stream, such as Request/Response (Aperiodic Pull), polling (Periodic Pull), Publish / Subscribe (Aperiodic Push) or Broadcast Disks (Periodic Push). Moreover, an info gatherer may be used to change the appearance of the real model elements delivery at the query language level. In this case, we have also to take into account the conformity points when model interpretation is coherent, for instance using an operation-based model representation [START_REF] Blanc | Detecting model inconsistency through operation-based model construction[END_REF].

The contributions of this paper are motivated by the need to define and manipulate infinite models (i.e. models whose comprehensive set of model elements is too large to be loaded or even not available). After pinpointing how current metamodeling formalisms prevent such situations, we first propose MOF extensions to locally characterize in a metamodel the infinite parts of the conforming models. Then we introduce a formal alternative semantics of the OCL operator iterate, called coiterate, providing an implementation-independent semantics for manipulating infinite models, and supporting specification and verification of complex situations involving infinite models. The coiterate operator can be used both to browse an infinite collection, and to compute the infinite closure of a transitive relation (or more generally a cycle in a model).

Such a coiterate operator can support the formal verification of operations manipulating infinite models. This verification activity can be partially automated by proof assistant which supports coinductive semantics. Among others, Coq is a valuable candidate in that respect [START_REF] Bertot | Coinduction in coq[END_REF]. This paper also discusses various possible implementations of the coiterate operator, whether for reasons of partial availability or partial loading of models.

More generally, iterate and coiterate iterators may be used for model transformation. Indeed, the accumulator may be a model (i.e., modeling in the large) which is finitely extended and returned (iterate) or indefinitely browsed in order to produce a new model (coiterate). If we abstract away the pieces of code used as arguments of these iterators, it turns out that iterate needs a function of type 9 1 + A × B → B whereas coiterate needs a function of type B → 1 + A × B. In order to apply these iterators to model production or browsing, the types A and B may be generalized to dependent types that denote arbitrary predicates over a (mega-)model structure, aiming at identifying patterns of interest and providing entry points in these patterns. In this context, we would need functions with respective types ∀e ∈ E.1 + A(e) × B(e) → B(e) and ∀e ∈ E.B(e) → 1 + A(e) × B(e), where e are elements of the (mega-)model. Here, A(e) and B(e) represent query or creation patterns, the functions induce transformation rules and the iterators represent the execution of a global transformation engine. Finally, specifying models and model transformations with type predicates allows to talk about their respective properties within a single language of types. Turning questions about models into typing problems also brings a rich amount of results in the scope, about type checking, type inference and subtyping issues for models and model transformations at once, as investigated in [START_REF] Steel | On model typing[END_REF].

To conclude with it, some substantial amount of work will be necessary to draw all the consequences of this promising approach for model transformation and to thoroughly compare and cross-fertilize it with existing solutions.

Fig. 1 .

 1 Fig. 1. Finite Model vs. Infinite Model

Definition 3 (

 3 Model Path). Let E, L ∈ Model(M E , R) be an infinite model. A model path is a relation path through model elements, i.e. a sequence of triples { e i , me i , r i , e i+1 , me i+1 } i∈I ∈ L, where either I = [0, sup[ is finite (sup ∈ N) or I = N is infinite. Relying on the previous model path definition, we also assimilate a model path to the (behavioral) trace of the model element creations. Definition 4 (Maximal Path). Let E, L ∈ Model(M E , R) be an infinite model. A maximal path is a model path, such that if it is finite, then the final element of the sequence has no relation to any element of the model.

  coll->iterate(elem: A; acc : B = e 1 | acc = e 2 ) i f coll->isEmpty() then e 1 else let e 1 = e 2 [coll->first() | elem][e 1 | acc] in coll->tail()->iterate(elem: A; acc : B = e 1 | acc = e 2 ) coll->coiterate(acc : B = e 1 | acc = e 2 ; elem: A = e 3 ) i f e 1 = null then coll else let e 1 = e 2 [e 1 | acc] in coll->append(e 3 [e 1 | acc])->coiterate(acc : B = e 1 | acc = e 2 ; elem: A = e 3 )

1 .

 1 Naturals Set{}->coiterate (acc: Integer = 0 | acc = acc+1; elem = acc) 2. Squares Set{}->coiterate (acc: Collection (Integer) = Naturals | acc = acc->tail()-> tail (); elem = acc->first () * acc->first ())

context

  StateMachine :: simulate () : Sequence(State) body : Set{}->coiterate ( acc : Sequence(EventOccurence) = eventToProcess ; current : State = self . initial | current = self . step ( current , acc->first ()) ; acc = acc->tail () ; elem = current )

  context T:: coclosure ( relation ) : Sequence(T) body : Set{}->coiterate ( acc : Sequence(State) = Sequence{self} | acc = acc->tail()->union(acc->first (). relation ->asSequence()) elem = acc->first () ; ) context State :: reachableStates () : Sequence(State) derive : { self }->coclosure( nextStates );

Finite unfoldings where an explicit bound is known could be interesting. We could then add an annotation on relations, belonging to the same type as the upper bound of the multiplicity of collections. This extension is not taken into account in the scope of this paper.

The operation tail() is defined in OCL by: coll->excluding(coll->first()).

e[u | u ] is the expression e where any occurrence of sub-term u has been replaced by u.

We consider the operation step() which returns the current state according to the event given as a parameter.

With memoization, computation of a model element is not repeated: the element is kept in a cache after its first usage.

In abstract set algebra settings, "1" denotes the set with only one element, "+" denotes the disjoint sum and "×" the cartesian product.

This work is partially supported by the EU FP7-ICT-2009.1.4 Project N • 256980, NESSoS: Network of Excellence on Engineering Secure Future Internet Software Services and Systems.