Florian Haider
email: florian.haider@onera.fr.

Bernard Courbet
email: bernard.courbet@onera.fr

Jean-Pierre Croisille
email: jean-pierre.croisille@univ-lorraine.fr

A High-Order interpolation for the Finite Volume Method: the Coupled Least Squares Reconstruction

Keywords:

We consider the Finite Volume Method on irregular grids for conservation laws with a particular polynomial reconstruction. This reconstruction is based on a least squares approach to compute consistent approximations of the first, second and third order derivatives. The resulting reconstruction is cubic. It is called the Coupled Least Squares reconstruction. It is obtained by a three stages iteration. At each iteration, only data located in a compact stencil, and not beyond, are accessed. A linear stability analysis is given in the case of regular and irregular one-dimensional grids. Numerical results for various problems, including shock tubes, vortex and accoustic wave propagation, support the interest of this approach. The reconstruction algorithm is presented in detail in the one dimensional case. An outline of the multidimensional case in also given.

.

Introduction

Finite volume schemes with high order reconstruction

In [START_REF] Haider | Efficient implementation of high order reconstruction in finite volume methods[END_REF][START_REF] Haider | Parallel implementation of k-exact finite volume reconstruction on unstructured grids[END_REF], a high order reconstruction approach for finite volume approximations was introduced for the purpose of gas dynamics simulations on general polyhedral grids. This approach, which is called the Coupled Least Squares reconstruction, belongs to the k-exact methods [START_REF] Barth | Higher order solution of the euler equation on unstructured grids using quadratic reconstruction[END_REF]. The purpose of this paper is to give a detailed account of a particular algorithm of the reconstruction. The presentation is mainly restricted to the one dimensional context to emphasize the logic of the numerical procedure.

Consider the conservation law

∂ t u + ∂ x f (u) = 0. (1.1)
We start from the discrete in space, continuous in time, integral version of (1.1):

dv α (t) dt = - 1 |T α | f α+1/2 (t) -f α-1/2 (t) . (1.2)
The grid is made of cells T α and is depicted on Fig 1 .1. In (1.2) t → v α (t) approximates t → ūα (t) defined by

ūα (t) = 1 |T α | Tα u(x, t)dx, (1.3)
T α-1

T α T α+1 x α-1 x α+1 x α h α-1,α h α,α+1
Figure 1.1: The cell T α and the two neighboor cells T α±1 . The barycenters of the cells T α±1 and T α are the points x α±1 and x α , repectively. In addition, h α-1,α = x α -x α-1 and h α,α+1 = x α+1 -x α .

where u(x, t) is the solution of (1.1). Throughout the paper, the numerical flux t → f α+1/2 (t) is of the form:

f α+1/2 (t) = F w α [V (t)](x - α+1/2 ), w α+1 [V (t)](x + α+1/2 ) . (1.4)
where:

• The function (u L , u R ) → F (u L , u R ) denotes some numerical flux function, typically the HLLC flux in the case of the Euler equations [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF].

• The reconstruction x → w α [V (t)](x) is an interpolant based on the vector V (t) = [v 1 (t), v 2 (t), . . . , v N (t)] T .

The reconstruction expounded in this paper has the form of a Taylor expansion of u(x) around the barycenter x α of the cell T α , i.e. (we take out the time dependency for clarity):

w α [ū](x) u(x α ) + u (x α )(x -x α ) + 1 2 u (x α )(x -x α ) 2 + 1 6 u (x α )(x -x α ) 3 . (1.5)
where ū = [ū α ] is the vector of averages of u(x) over the cells T α . The relation (1.5) is a 3-exact reconstruction. The main task is to define a consistent approximation to u (x α ), u (x α ) and u (x α ), where u(x) is supposed to be known by ū only.

The main observation is the following: to preserve fourth order accuracy in (1.5), it is sufficient to calculate approximations of u (x α ), u (x α ) and u (x α ) with order 3, 2 and 1, respectively. It is the purpose of this paper to show how this can be obtained in several steps of an iterative procedure along which the accuracy is progessively enhanced. Moreover, at each step, there is an access only to data located in a small neighborhood of any given cell, and therefore it is called a compact k-exact reconstruction.

This paper is the companion paper of [START_REF] Haider | Parallel implementation of k-exact finite volume reconstruction on unstructured grids[END_REF] where the multidimensional aspect on general polyhedral grids of our algorithm is addressed. This study has been initiated in [START_REF] Haider | Discrétisation en maillage non structuré général -Applications à la LES[END_REF], in the context of the package CEDRE, a parallel multisolver CFD code for aerothermochemistry 1 .

Our point of view is thus the one of a generic computational procedure for any finite volume code and for general compressible flows. For an overview on the k-exact approach for high order computations in CFD, compared to other approaches, we refer to the review paper [START_REF] Vincent | Facilitating the adoption of unstrutured hig-order methods amongst a wider community of fluid dynamicists[END_REF].

High order finite volume schemes

We consider the finite volume method (1.2) for a conservation law (1.1) on a general irregular grid made of cells T α , which support approximations v α (t) of the averages ūα (t). As explained in Section 1.1, we want to define the reconstructed polynomial w α [V (t)](x) which is used in the numerical flux f α+1/2 (t) in (1.4). This polynomial must have the following properties (we drop the t dependency).

First w α [V ](x) depends linearly on the data v β in the cells T β belonging to a small neighborhood of T α . In particular, it does not use gradient data as for example in the Discontinuous Galerkin approximation, but only finite volume data. Furthermore, no flow based local parameters are used in a nonlinear way to define w α [V ](x). Consequently, our reconstruction is decoupled from all nonlinear treatments such as slope limiting or monotonicity preserving procedures. Thus the reconstruction step is decoupled from the limitation procedure. The reconstruction step can be applied in principle to any finite volume code including complex fluid modelling such as turbulent, multiphasic, reactive, spherical flows, etc. With this approach, the nonlinear treatments are performed in a second step, after the reconstruction step. Having in view general CFD codes, the design of the reconstruction must of course be multidimensional. Moreover, it must allow an easy parallel implementation. This is obtained by mean of the procedure in several stages presented in Section 2 and Section 3.

Related works

Reconstruction beyond linear in the finite volume method on irregular grids has been a longstanding issue of great importance, both practical and theoretical in CFD. It takes place in the challenge to go beyond second order with the finite volume method. The problem was introduced in [START_REF] Harten | Multidimensional ENO schemes for general geometries[END_REF][START_REF] Abgrall | On the use of Mühlbach expansions in the recovery step of ENO methods[END_REF] where various polynomial interpolations on triangles were studied. Quadratic reconstruction is undertaken in [START_REF] Barth | Higher order solution of the euler equation on unstructured grids using quadratic reconstruction[END_REF][START_REF] Delanaye | Quadratic-reconstruction finite volume scheme for compressible flows on unstructured adaptive grids[END_REF]. The series of studies [START_REF] Ollivier-Gooch | A high-order accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Michalak | Unstructured high-order accurate finite-volume solutions of the Navier-Stokes equations[END_REF] explores the finite volume method with quadratic reconstruction for compressible Fluid Dynamics. In many works, the reconstruction is analyzed in conjunction with the limiter and with the positivity preserving property [START_REF] Barth | Finite Volume Methods: Foundation and Analysis[END_REF][START_REF] Barth | Numerical methods for conservation laws on structured and unstructured meshes[END_REF], which are considered as constraints on the reconstruction. This is for example the case in the well known ENO/WENO procedure. In this approach, local high order reconstruction functions are based on stencils which vary according to local sensors related the solution, e.g. the total variation. We refer to [START_REF] Harten | ENO Schemes with subcell resolution[END_REF][START_REF] Jiang | Efficient implementation of weighted ENO schemes[END_REF][START_REF] Shu | High order ENO and WENO schemes[END_REF][START_REF] Sonar | On the construction of essentially non-oscillatory finite volume approximations to hyperbolic conservation laws on general triangulations: polynomial recovery, accuracy and stencil selection[END_REF]. Another option [START_REF] Dumbser | Arbitrary high order P N P M schemes on unstructured meshes for the compressible Navier-Stokes equations[END_REF] is where the reconstruction in finite volumes is considered in the general framework of the Discontinuous Galerkin approach. We also refer to the studies [START_REF] Khosla | Using fourth order spatial integration on unstructured meshes to reduce LES run time[END_REF][START_REF] Winkler | The effect of unstructured grid topology and resolution on simulations of decaying turbulence[END_REF], with interest into application oriented simulations, where reconstructions beyond linear are presented. Finally we mention the recent study [START_REF] Pont | Multiple-correction hybrid kexact schemes for high-order compressible RANS-LES simulations on fully unstructured grids[END_REF] where a variant of compact k-exact reconstruction is considered. In this study the authors use an upwinded scheme only up to third order and switch to a fourth order centered scheme in vorticity dominated regions. This is in contrast to our scheme which is fully upwinded up to fourth order independently of the flow regime.

Outline

The outline of the paper is as follows. In Section 2 we present the principle of our reconstruction algorithm. In Section 2.1, the general d-dimensional framework is set up for the reconstructed cubic polynomial w α [ū](x). From Section 2.2 on, we switch to the one dimensional case to present how the reconstruction uses at each step only data in the direct neighborhood of T α . In Section 3 and Section 4 the detail of the algorithm to calculate the cubic reconstruction is summarized (algorithm 1). This is based on fully centered approximate derivatives σα , θα and ψα in the cubic reconstruction (1.5) giving rise to a linearly stable semidiscrete system (1.2), even on irregular grids. Our analysis is based on the numerical computation of discrete spectra. Finally we show in Section 5 numerical results in one and multidimension illustrating our approach.

Coupled Least Squares Reconstruction

The d-dimensional case

The Coupled Least Squares method is a particular method to implement a k-exact reconstruction. This method has been already been introduced in [START_REF] Haider | Parallel implementation of k-exact finite volume reconstruction on unstructured grids[END_REF]. Here we specifically focus on the piecewise cubic reconstruction (i.e k = 3). Let Ω ⊂ R d , d ≥ 1 be a domain with periodic boundary conditions. We consider a grid, a priori irregular, consisting of N general polyhedral cells T α as on Fig. 2.1. The tensor notation used in this section is summarized in Appendix 8. The reconstruction algorithm proceeds in several steps. At each step, the only data that are used are located in a small neighborhood W α of the cell T α . The set W α is a parameter of the method. A typical choice of W α is the first neighborhood 2 defined by

x α T T α β h αβ x β
V (1) α = {β ∈ [1 : N ], | T β is adjacent to T α }. (2.2)
Let x ∈ Ω → u(x) be a given regular function. The value of u and of the first, second and third order derivatives at x α are denoted by 1) u α = D (1) x u(x α ), D (2) u α = D (2) x u(x α ), D (3) u α = D (3) x u(x α ).

u α = u(x α ), D (
(2.

3)

The derivative D (l) u α is a symmetric tensor of order l. We also note the mean value

ūα = 1 |T α | Tα u(x)dx. (2.4)
Consider the particular case where u(x) = p(x) is a cubic polynomial. The values pβ -pα and the three derivatives D (l) p α , l = 1, 2, 3, satisfy for all β ∈ W α the following Taylor expansions (see Appendix 8 for the notation),

           pβ -pα = D (1) p α • (z (1) αβ -x (1)
α ) + 1 2 D (2) p α • (z (2) αβ -x (2) α ) + 1 6 D (3) p α • (z (3) αβ -x (3)
α ),

D (1) p β = D (1) p α + D (2) p α • h (1) αβ + 1 2 D (3) p α • h (2) αβ , D (2) p β = D (2) p α + D (3) p α • h (1)
αβ .

(2.5)

2
The first neighborhood of the cell Tα is also called the direct or the Von Neumann neighborhood.

For a general function u(x) these Talor expansions become

           ūβ -ūα = D (1) u α • (z (1) αβ -x (1) α ) + 1 2 D (2) u α • (z (2) αβ -x (2) α ) + 1 6 D (3) u α • (z (3) αβ -x (3) α ) + O(h 4 ), D (1) u β = D (1) u α + D (2) u α • h (1) αβ + 1 2 D (3) u α • h (2) αβ + O(h 3 ), D (2) u β = D (2) u α + D (3) u α • h (1) αβ + O(h 2 ).
(2.6) Therefore, due to the O(h p ) terms in (2.6), the quantities ūβ -ūα , D (1) u α , D (2) u α and D (3) u α are approximate solutions of the system (2.5). Consider now the problem to approximate the three derivatives D (1) u α , D (2) u α and D (3) u α by consistent values. We look for symmetric tensors called σα , θα and ψα

satisfying          σα = D (1) u α + O(h 3 ), θα = D (2) u α + O(h 2 ), ψα = D (3) u α + O(h).
(2.7)

An essential observation is that for all symmetric tensors σα , θα and ψα of order 1, 2 and 3 respectively, satisfying (2.7), the same system than (2.6) holds. For all β ∈ W α we have:

           ūβ -ūα = σα • (z (1)
αβ -x

(1)

α ) + 1 2 θα • (z (2) αβ -x (2) α ) + 1 6 ψα • (z (3) αβ -x (3) α ) + O(h 4 ), σβ = σα + θα • h (1) αβ + 1 2 ψα • h (2) αβ + O(h 3 ), θβ = θα + ψα • h (1) αβ + O(h 2 ). (2.8)
The Coupled Least Squares method consists in approximately solving (2.8), when taking out the O(h p ) terms. Otherwise stated, we solve for ū

= [ū 1 , ū2 , . . . ūN ] → σα [ū], θα [ū], ψα [ū] solution of the system            ūβ -ūα = σα [ū] • (z (1)
αβ -x

(1)

α ) + 1 2 θα [ū] • (z (2)
αβ -x

(2)

α ) + 1 6 ψα [ū] • (z (3) αβ -x (3)
α ), (a),

σβ [ū] = σα [ū] + θα [ū] • h (1) αβ + 1 2 ψα [ū] • h (2)
αβ , (b),

θβ [ū] = θα [ū] + ψα [ū] • h (1)
αβ , (c).

(2.9) Solving (2.9) raises two questions. First the linear system (2.9) does not have a solution in general. This system is rectangular with m α 1 + d + d(d -1)/2 equations and d + d(d + 1)/2 + d(d + 1)(d + 2)/6 unknowns. For example in two dimensions and for a grid made of triangles, there are 18 equations and 9 unknowns. Therefore (2.9) is overdetermined. This leads to seek a solution in the least squares sense. Second, equations (2.9 a ), (2.9 b ) and (2.9 c ) couple unknowns in differents cells. To overcome these problems, the following strategy is used. The system (2.9) is rewritten as:

           h αβ . σα [ū] = ũβ -ūα -1 2 θα [ū] • (z (2) αβ -x (2) α ) -1 6 ψα [ū] • (z (3) αβ -x (3)
α ), (a)

h αβ . θα [ū] = σβ [ū] -σα [ū] -1 2 ψα [ū] • h (2)
αβ , (b)

h αβ . ψα [ū] = θβ [ū] -θα [ū] (c).
(2.10)

Let us now consider in each equation of (2.10) the unknowns in the right-hand-side as parameters. A least squares resolution is applied to each equation in (2.10). Consider first (2.10) a . The right-hand side is denoted by

δu (0) αβ = ũβ -ūα - 1 2 θα [ū] • (z (2) αβ -x (2) α ) - 1 6 ψα [ū] • (z (3) αβ -x (3) α ), (2.11)
Therefore (2.10) a is rewritten as

h αβ . σα [ū] = δu (0) αβ , ∀β ∈ W α . ( 2
min σα∈R d β∈Wα | δu 0 αβ -σ α • h (1) α | 2 . (2.13)
On the other hand, the standard least squares slope in cell T α is solution of

min σα∈R d β∈Wα |ū β -ūα -σ α • h (1) α | 2 . (2.14) Denote C α and H T α the d × m α matrices (m α = #W α ) C α = [c αβ1 , c αβ2 , . . . c αβm α ], H T α = [h αβ1 , h αβ2 , . . . h αβm α ]. (2.15)
Assuming rank(H α ) = d, the unique solution of (2.14) is given by the least squares operator A α defined by

A α : [ū] → σLS α [ū] = β∈Wα c αβ (ū β -ūα ), (2.16)
where

c αβ = (H T α H α ) -1 h (1)
αβ (= O(1/h)).

(2.17)

As a result, σα [ū] solution of (2.13) is expressed as:

σα [ū] = σLS α [ū] -ā(3) α . θα [ū] -b(4) α . ψα [ū]. (2.18)
where the two tensors ā(3) α and b(4) α depend only on the grid and are given by

ā(3) α = z (2) αβ -x (2) α ⊗ c αβ , b(4) α = z (3) αβ -x (3) α ⊗ c αβ . ( 2
αβ = σβ [ū] -σα [ū] - 1 2 ψα [u] • h (2) αβ , β ∈ W α . (2.20)
In (2.10) b , the symmetric tensor θα [ū] is the unknown and δu

(1)
αβ is at the moment considered to be known. There are d(d + 1)/2 unknowns and dm α equations. The linear system is as before overdetermined. It can be expressed as

H (2) α θα [ū] = δu (1) αβ β , (2.21)
where the size of the matrix H

(2) where δu

α is m α d × d(d +
(1)

α = [ δu (1)
α1 , . . . , δu

αm α ] T and where H

(2) α + denotes the standard Moore-Penrose inverse of the matrix H

α . It is defined by [START_REF] Golub | Matrix Computations[END_REF][START_REF] Deuflhard | Numerical Analysis in Modern Scientific Computing[END_REF]:

H (2) α + = H (2) α T H (2) α -1 H (2) α T
.

(2.23)

Therefore the second relation in the space of symmetric tensors of order 2 (i.e. in R d(d+1)/2 ) is

θα [ū] = H (2) α + σβ [ū] -σα [ū] - 1 2 ψα [ū] • h (2)
αβ β .

(2.24)

Proceeding similarly for the third relation (2.10) c provides a relation of the form

ψα [ū] = H (3) α + θβ [ū] -θα [ū] β . (2.25)
The relation (2.25) is an identity in the space of symmetric tensor of order 3 and H

α again is a suitable least squares matrix. The final linear system with unknowns σα [ū], θα [ū] and ψα [ū] is of the form

           σα [ū] + ā(3) α . θα [ū] + b(4) α . ψα [ū] = σLS α [ū], (a) θα [ū] + a (3)
α . ψα [ū] = θLS α [ū], (b) ψα [ū] = ψLS α [ū] (c). (2.26)
In (2.26) we have denoted

   θLS α [ū] = H (2) α + ( σβ [ū] -σα [ū]) β , ψLS α [ū] = H (3) α + ( θβ [ū] -θα [ū]) β .
(2.27)

and where a

α is a third tensor depending only on the grid.

Definition 2.1. The system (2.26) is called the Coupled Least Squares system for the cubic reconstruction.

Two methods to solve (2.26) can be considered. First, the Jacobi (or Gauss Seidel) iteration can be used. This approach is presented in [START_REF] Haider | Parallel implementation of k-exact finite volume reconstruction on unstructured grids[END_REF]. Due to the triangular form of (2.26), a second method consists in using a forward/backward substitution method. The rest of the paper is devoted to the presentation of this second method in the particular case of the one dimensional setting. The interest of the one dimensional case is that it permits to have a better view on the algebra of the algorithm. In particular, it permits to make fully explicit the two steps of the method: first establishing the system (2.26) and second solving it. Note in addition that, even in one dimension, the reconstruction procedure on a nonequispaced grid is by no way straightforward. It clearly deserves an independent study in itself. Having solved (2.26) in one way or another, the piecewise cubic reconstruction is given by

w α [ū](x) = ūα + σα [ū]• (x-x α )-x (1) α + 1 2 θα [ū]• (x-x α ) ⊗2 -x (2) α + 1 6 ψα [ū]• (x-x α ) ⊗3 -x (3) α . (2.28)
Remark 2.2. A full mathematical study in the d-dimensional case, d ≥ 2 is not completed yet. An important (and difficult) issue is to determine the geometric condition on the grid under which the matrices H

(2) α + and H

(3) α + do exist. In the absence of such a condition so far, we rely on numerical experiments.

The one-dimensional case

From now on, we consider the one dimensional case. First we give the detail of the agebra in Section 2.3 in the one dimensional case (d = 1). This permits to emphasize the logic of the reconstruction procedure. In this case, the calculations presented in Section 2.1 are much simpler since all tensors reduce to scalars. The neighborhood

W α is W α = V (1) α with 3 V (1) α = {T α-1 , T α , T α+1 }. (2.

29)

The vector h αβ is defined by (see Fig. 1.1)

h αβ = x β -x α , β = α -1, α, α + 1. (2.30)
In particular, h αα = 0. The hypothesis (8.12) on the irregular grid becomes

ch ≤ h α ≤ Ch, ch ≤ h αβ ≤ Ch, (2.31)
where c and C are positive constants. Let u(x) be a given regular function. The values u α , σ α , θ α and ψ α are (see (2.3)):

u α = u(x α ), σ α = u (x α ), θ α = u (x α ), ψ α = u (x α ). (2.32)
The Taylor expansion of u(x) at x α is:

u(x) = u α + σ α (x -x α ) + 1 2 θ α (x -x α ) 2 + 1 6 ψ α (x -x α ) 3 + O(h 4 ). (2.33)
The average ūα of u(x) over T α satisfies:

ūα = u α + |T α | 2 24 θ α + O(h 4 ), (2.34)
Inserting ūα in (2.33) in replacement of u α leads to expanding u(x) in the form:

u(x) = ūα + σ α (x -x α ) + θ α 1 2 (x -x α ) 2 - 1 24 |T α | 2 + 1 6 ψ α (x -x α ) 3 + O(h 4 ). (2.35)
This is why we look for an approximant w α [ū](x) of u(x) of the form:

w α [ū](x) = ūα + σα (x -x α ) + θα 1 2 (x -x α ) 2 - 1 24 |T α | 2 + 1 6 ψα (x -x α ) 3 . (2.36)
where σα , θα and ψα are suitable approximations to σ α , θ α and ψ α , respectively. Note that (2.36) is the one dimensional form of the multidimensional reconstruction (2.28). Defining the reconstruction w α [ū](x) is therefore equivalent to define σα , θα and ψα and this is the topic of the rest of the paper.

In order to define a fourth order VF scheme, the reconstruction (2.36) is finally inserted in the flux function as indicated in (1.4).

Coupled Least Squares at the continuous level

In this section, we show how to calculate σα , θα and ψα appearing in (2.36). As already mentionned, data are located in the compact stencil W α = V (2.37)

with coefficients c αβ to be defined. The consistency relation is

β∈Wα c αβ h αβ = 1. (2.38)
The relation (2.38) express the property that the slope σ 0 of a linear function u(x) = u 0 + σ 0 x is left invariant by the slope operator (2.37). Our privilegiate choice for c αβ in (2.38) is the least squares slope defined by

σLS α = β∈Wα c LS αβ (ū β -ūα ), (2.39)
where

c LS αβ is (see (2.17)), c LS αβ = h αβ / β∈Wα h 2 αβ .
(2.40)

The least squares slope is consistent since it satisfies (2.38). Henceforth we write c αβ instead of c LS αβ . Note that

c αβ = O(1/h). (2.41)
In the sequel, we also need the least squares values σ LS α , θ LS α and ψ LS α defined by the relations:

               σ LS α = β∈Wα c αβ (u β -u α ), (a) θ LS α = β∈Wα c αβ (σ β -σ α ), (b) ψ LS α = β∈Wα c αβ (θ β -θ α ). (c) (2.42) Lemma 2.3. The approximate derivatives σ LS α , θ LS α , ψ LS α in (2.42) satisfy the consistency relations          σ α + a α θ α + b α ψ α = σ LS α + O(h 3 ), (a) θ α + a α ψ α = θ LS α + O(h 2 ), (b)
ψ α = ψ LS α + O(h), (c) (2.43) with          σ α = u (x α ), θ α = u (x α ), ψ α = u (x α ).
(2.44)

Using the l-momentum H (l)
α in T α defined by:

H (l) α = β∈Wα h l αβ , (2.45)
the coefficients a α , b α in (2.43) are:

a α = 1 2 H (3) α H (2) α , b α = 1 6 H (4) α H (2) α . (2.46)
Furthermore, replacing in (2.43) a the pointwise least squares value σ LS α by σLS α , gives

σ α + āα θ α + bα ψ α = σLS α + O(h 3 ), (2.47)
where the modified coefficients āα and bα are

āα = a α + 1 24 β∈Wα c αβ |T β | 2 -|T α | 2 and bα = b α + 1 24 β∈Wα c αβ |T β | 2 h αβ . (2.48) Proof. For all cells T β ∈ W α , the Taylor expansion of u(x β ) -u(x α ) is u β -u α = σ α h αβ + 1 2 θ α h 2 αβ + 1 6 ψ α h 3 αβ + O(h 4 ). (2.49) Therefore the slope σ LS α in (2.42) a satisfies σ LS α = β∈Wα c αβ σ α h αβ + 1 2 θ α h 2 αβ + 1 6 ψ α h 3 αβ + O(h 3 ). ( 2
| 2 = O(h 2 )), σ LS α = β∈Wα c αβ (ū β -ūα ) - 1 24 β∈Wα c αβ (|T β | 2 θ β -|T α |θ 2 α ) + O(h 4 ) = σLS α - 1 24 β∈Wα c αβ |T β | 2 (θ α + h αβ ψ α ) -|T α | 2 θ α + O(h 4 ) = σLS α - θ α 24 β∈Wα c αβ |T β | 2 -|T α | 2 - ψ α 24 β∈Wα c αβ |T β | 2 h αβ + O(h 4 ).
Replacing in the preceding relation σ LS α by (see (2.43) a )

σ LS α = σ α + a α θ α + b α ψ α + O(h 3 ), (2.51)
we deduce that

σ α + θ α a α + 1 24 β∈Wα c αβ |T β | 2 -|T α | 2 + ψ α b α + 1 24 β∈Wα c αβ |T β | 2 h αβ = σLS α + O(h 3 ).
This proves (2.47), with the coefficients āα , bα given in (2.48).

A Cubic Reconstruction Procedure

In this section the Coupled Least Squares resolution procedure is detailed. We show how to use Lemma 2.3 to define our cubic reconstruction. Let us start from the relations (2.47) and (2.43) bc above relating σ α , θ α , ψ α , σLS α , θ LS α and ψ LS α . These relations are, (compare (2.26)):

         σ α + āα θ α + bα ψ α = σLS α + O(h 3 ), (a) θ α + a α ψ α = θ LS α + O(h 2 ), (b)
ψ α = ψ LS α + O(h). (c) (3.1)
Recall that σ α , θ α and ψ α denote the exact derivatives u (x α ), u (x α ) and u (x α ), respectively. However, consider for a moment the system (3.1) as a 3×3 linear system with σ α , θ α and ψ α as "unknowns".

A forward-backward "resolution" provides Taylor expansions of σ α , θ α and ψ α in terms of the value σLS β . This is expressed in the following

Proposition 3.1. The exact derivatives u (x α ) = σ α , u (x α ) = θ α , u (x α ) = ψ α satisfy the relations                ψ α = β∈Wα c * αβ γ∈V β cβγ σLS γ -σLS β - δ∈Wα cαδ (σ LS δ -σLS α ) + O(h), (a)
θ α = -ã α ψ α + β∈Vα cαβ σLS β -σLS α + O(h 2 ), (b) σ α = -ā α θ α -bα ψ α + σLS α + O(h 3 ), (c) (3.2)
where

• The coefficients āα and bα are given in (2.48),

• The coefficients ãα and cαβ are defined by:

ãα = a α + β∈Wα c αβ ( bβ -bα + āβ h αβ ) 1 + β∈Wα c αβ (ā β -āα ) , (3.3
)

and cαβ = c αβ 1 + β∈Wα c αβ (ā β -āα ) . (3.4)
• The coefficient c * αβ is given by

c * αβ = c αβ 1 + δ∈Wα c αδ (ã δ -ãα ) . (3.5) Proof. First note that Equation (3.2) (c) is identical to Equation (3.1) (a)
and therefore there is nothing to prove. Second, consider the proof of (3.2)

(b) . The relation (3.1) (b) is θ α + a α ψ α = θ LS α + O(h 2 ) (3.6)
Replacing in the right-hand side θ LS α by its value in (2.42) (b) gives

θ α + a α ψ α = β∈Wα c αβ (σ β -σ α ) + O(h 2 ). (3.7)
Consider now the equation (2.47) in all cells T β

σ β + āβ θ β + bβ ψ β = σLS β + O(h 3 ). (3.8)
For all cells T β close to T α , a Taylor expansion gives

θ β = θ α + h αβ ψ α + O(h 2 ), (a) ψ β = ψ α + O(h). (b) (3.9)
Replacing in (3.8) the values of θ β and ψ β given in (3.9) yields (note that āα = O(h) and bα = O(h 2 ))

σ β + āβ (θ α + h αβ ψ α ) + bβ ψ α = σLS β + O(h 3 ), (3.10)
or equivalently

σ β = σLS β -āβ θ α -(ā β h αβ + bβ )ψ α + O(h 3 ). (3.11)
Inserting finally (3.11) in (3.8) gives

θ α + a α ψ α = β∈Wα c αβ σLS α -āβ θ α -(ā β h αβ + bβ )ψ α -σ α + O(h 2 ) = β∈Wα c αβ σLS β -āβ θ α -(ā β h αβ + bβ )ψ α -σLS α + āα θ α + bα ψ α + O(h 2 ).
Collecting the terms in θ α and ψ α in the left-hand-side gives

1 + c αβ (ā β -āα ) θ α + a α + β c αβ ( bβ -bα + āβ h αβ ) ψ α = β∈Wα c αβ σLS β -σLS α + O(h 2 ).
Dividing each side by the coefficient

1 + β∈Wα c αβ (ā β -āα ), (3.12)
this is rewritten as

θ α + ãα ψ α = β∈Wα cαβ σLS β -σLS α + O(h 2 ), (3.13)
where ãα , cαβ are given in (3.4). Therefore (3.2) (b) is proved and we prove (3.2) (c) in a similar way. Using (3.13 ) at x β and the relation

ψ β = ψ α + O(h) gives the following expression for θ β θ β = -ã β ψ α + γ∈W β cβγ σLS γ -σLS β + O(h 2 ). (3.14)
The identity (2.42) (c) is expressed using (3.1), (3.14) and (3.9) (b) as

ψ α = β∈Wα c αβ (θ β -θ α ) = β∈Wα c αβ -ãβ ψ β + γ∈W β cβγ σLS γ -σLS β -θ α + O(h) = β∈Wα c αβ ãα ψ α -ãβ ψ β + γ∈W β cβγ σLS γ -σLS β - δ∈Wα cαδ (σ LS δ -σLS α ) + O(h) = ψ α - β∈Wα c αβ (ã β -ãα ) + β∈Wα c αβ γ∈W β cβγ σLS γ -σ LS β - δ∈Wα cαδ (σ LS δ -σLS α ) + O(h).
Collecting in the left hand side all the terms in ψ α gives:

1 + δ∈Wα c αδ (ã δ -ãα ) ψ α = β∈Wα c αβ γ∈W β cβγ σLS γ -σLS β - δ∈Wα cαδ (σ LS δ -σLS α ) + O(h).
Dividing by the coefficient 1 + δ c αδ (ã δ -ãα ) yields (3.2) (a) where c * αβ is given in (3.5) and the proof is complete. Finally we define the quantities ψα , θα and σα by droping in the relations (3.2) the O(h p ) terms. This suggests the following Corollary 3.2. Let us define the values ψα , θα and σα by

             ψα = β∈Wα c * αβ γ∈W β cβγ σLS γ -σLS β - δ∈Vα cαδ (σ LS δ -σLS α ) , (a) θα = -ã α ψα + β∈Wα cαβ (σ LS β -σLS α ), (b) σα = -ā α θα -bα ψα + σLS α . (c) (3.15)
These three values are approximations to ψ α , θ α and σ α respectively with the following accuracy: Our cubic reconstruction has finally the following form:

         ψα = ψ α + O(h), (a)
Proposition 3.3. Let ū = [ū 1 , ū2 , .
. . , ūN ] be the vector of the averages of a given function u(x) over the cells T α of the grid. Then the cubic reconstruction polynomial defined by (2.36) is a fourth order accurate reconstruction of u(x) near T α .

Proof. This is a simple consequence of (2.35).

Remark 3.4. The difficulty mentioned in Remark 2.2 does not exist here, since the system (2.43) (droping the O(h p ) terms) consists of 3 equations with 3 unknowns σα , θα and ψα . Note however that the logic of the method is the same in the one-dimensional and the multidimensional cases.

Remark 3.5. The approximations σα , θα and ψα depend by construction on data located in cells

T β with β ∈ {α -3, α -2, α -1, α, α + 1, α + 2, α + 3}.
This dependence stencil of the reconstruction is not the accessed stencil at each step of Algorithm 1, which is the compact stencil

W α = {T α-1 , T α , T α+1 }.
Remark 3.6. The preceding analysis does not apply at boundary cells since the neighborood set V α for such cells is not complete. The question of defining a fourth order reconstruction at boundary cells is difficult. It may require one-sided approximations. We do not elaborate further on this point in this paper. In practice, second-order reconstruction is used in boundary cells.

A Fourth Order Finite Volume method

Semi-discrete in space approximation

Let us come back to the conservation law (1.1)

∂ t u + ∂ x f (u) = 0. (4.1)
The semi-discrete scheme in (1.2) is

dv α (t) dt = - 1 |T α | f α+1/2 (t) -f α-1/2 (t) , 1 ≤ α ≤ N, (4.2)
where

• The vector V (t) = [v 1 (t), v 2 (t), • • • , v N (t)] T ∈ R N has components v α (t) defined by: v α (t) ūα (t) = 1 |T α | Tα u(x, t)dx, 1 ≤ α ≤ N, (4.3)
and where u(x, t) is the solution of (4.1).

• The reconstruction operator

V = [v 1 , v 2 , . . . , v N ] → w[V ] α (x)
is given by (2.36).

• The numerical flux function is denoted by

(u L , u R ) → F (u L , u R ).
Consider the particular case of the linear convection equation

∂ t u + c∂ x u = 0, c > 0, x ∈ [0, 1], (4.4)
with periodic boundary conditions. In this case, the function t → ūα (t) satisfies the equation:

dū α (t) dt = - c |T α | u(x α,α+1 , t) -u(x α-1,α , t) , (4.5)
In this case, the semi-discrete scheme (4.2) is:

dv α dt = - c |T α | w α [V (t)](x α,α+1 ) -w α-1 [V (t)](x α-1,α ) . (4.6)
where the polynomial w α [V ](x) in (2.36) is calculated by the Algorithm 1 in Section 3. In matrix form, the scheme (4.6) is expressed as the linear dynamical system:

d dt V (t) = JV (t), J ∈ M N (R), (4.7)
with

V (t) = [v 1 (t), v 2 (t) • • • , v N (t)] T ∈ R N .
The coefficients J αβ of J are such that:

N β=1 J αβ V β = - c |T α | w α [V ](x α,α+1 ) -w α-1 [V ](x α-1,α ) . (4.8)
Note that a proof of the fourth order accuracy in space of the scheme (4.6) is by no way easy in the case of a general irregular grid. In the case of a regular grid, we have the following Proposition 4.1. On a regular grid with stepsize h, the scheme (4.6) is fourth order accurate with respect to the equation

dū α (t) dt = - c h u(x α,α+1 , t) -u(x α-1,α , t) . (4.9)
Proof. We refer to the Appendix for a proof. Further comments in the particular case of an irregular grid are also given.

Statistical analysis for linear stability

An important question concerning the reconstruction (2.36) is whether it provides stability of the linear system (4.7). This is related to the location of the spectrum of the matrix J in the complex plane. Due to the irregularity of the grid, this spectrum is not analytically known in general, and therefore we must rely on some numerical evaluation. We refer to [START_REF] Haider | Stability analysis of the cell centered finite-volume MUSCL method on unstructured grids[END_REF] for the same stability problem in three dimensions and in the case of a piecewise linear reconstruction, Consider a periodic grid with cells of length |T α | = κ α h. The coefficient κ α stands for a measure of the irregularity of the grid. The κ α satisfy α κ α = 1. The spectral abscissa of J is defined by:

Λ = max 1≤k≤N Re(λ k ).
(4.10) The dynamical system (4.7) is stable when Λ ≤ 0. We adopt here an elementary statistical point of view. The spectral absissa Λ is evaluated using a sample of irregular grids. Each grid of the sample has N cells. The irregularity factor κ α , 1 ≤ α ≤ N , is randomly selected so that

κ α ∼ 1 + Cr α > 0, (4.11)
where

C ∈ [0, 1) is a constant and r = [r α ] ∈ [-1, 1]
N is a random array of size N . The magnitude of the constant C determines the irregularity of the grid. Testing a large sample of such grids provides a statistical evaluation of the value of the spectral abscissa Λ. Doing so, no instability was observed, to computer accuracy. On Fig. 4.1 (left panel), we display the spectral abscissa as a function of the index k of 5000 randomly selected grids of size N = 64. This was also performed (not shown) using a sample of 1000 grids of size N = 32. In both cases, the irregularity constant in (4.11) was choosen as C = 0.99, which corresponds to highly irregular grids. As can be observed, Λ 0 to computer accuracy. On Fig. 4.1 (right panel) the distribution of Λ is represented, which appears to be Gaussian. Furthermore a T-test at the confidence level α = 0.05 with null hypothesis on the mean value Λ = 0 (vs Λ = 0) provides a non-rejection and a confidence interval of [Λ min , Λ max ] = [-0.0085(-14), 0.14 (-14)] (for Λ). This clearly indicates that the eigenvalues of the matrix J are located in the left complex half-plane and therefore we can infer that that the linear dynamical system (4.6) is stable. Note that due to the irregularity of the grid, a theoretical localization of Λ seems a difficult problem of spectral analysis.

Dissipation and dispersion analysis

In the particular case of a regular grid with step-size h, the approximation (4.6) can be analytically derived. The coefficients

H (2) α , H (3)
α and H (4) α in (2.45) are

H (2) α = 2h 2 , H (3) α = 0, H (4) α = 2h 4 .
The coefficients a α , ãα , āα , b α , bα in (2.46) are

a α = ãα = āα = 0, b α = h 2 6 , bα = 5 24 h 2 .
The coefficients c α,α+1 , c α,α-1 , c * α,α+1 , c * α,α-1 in (2.42) and (3.5) are

c α,α+1 = 1 2h , c α,α-1 = - 1 2h , c * α,α+1 = 1 2h , c * α,α-1 = - 1 2h .
The least squares slope is

σLS α = ūα+1 -ūα-1 2h . (4.
12)

The values of ψα , θα , ψα in (3.15) are

               ψα = 1 (2h) 2 σLS α+2 + σLS α-2 -2σ LS α , θα = 1 2h σLS α+1 -σLS α-1 , σα = 53 48 σLS α -5 96 σLS α+2 -5 96 σLS α-2 .
(4.13)

In (4.6), w α [V ](x α,α+1
) is given by:

w α [V ](x α,α+1 ) = w α [V ](x α + h/2) (4.14) = ūα + 1 2 hσ α + 1 12 h 2 θα + 1 48 h 3 ψα .
The scheme (4.6) is therefore

dv α dt = - c h v α (t) + h 2 σα (t) + h 2 12 θα (t) + h 3 48 ψα (t) (4.15)
- v α-1 (t) + h 2 σα-1 (t) + h 2 12 θα-1 (t) + h 3 48 ψα-1 (t) .
We now give a dissipation and dispersion analysis of (4.15). Such an analysis is an important indication to interpret the numerical behaviour of the scheme. For each given wavenumber k, -N/1+1 ≤ k ≤ N/2, we consider the periodic initial function u k 0 (x):

u k 0 (x) = exp(2ikπx). (4.16)
The discrete initial data at x α = αh is

u k 0,α = exp(iαφ k ) where φ k = 2πkh ∈ (-π, π]. (4.17)
The solution of the system (4.7) is:

u k α (t) = exp(2ikπ(x α -ĉk t)), (4.18)
where the numerical velocity ĉk , which depends on k, is decomposed into real and imaginary parts as

ĉk = ĉk,R + iĉ k,I . (4.19)
Definition 4.2. The dissipation and the dispersion functions are defined by

   φ k ∈ [0, π[ → D 1 (φ k ) = ĉk,I c , (a) φ k ∈ [0, π[ → D 2 (φ k ) = ĉk,R c (b).
(4.20)

• The dissipation function D 1 (φ k ) is positive when the sinusoidal gridfunction function z k α = exp(iαφ k ) is amplified for the wave number k. This corresponds to an unstable scheme. In the contrary, if D 1 (φ k ) ≤ 0 for all k, then the scheme is stable.

• The dispersive function D 2 (φ k ) ≥ 1 for a propagation at a velocity ĉk,R > c and D 2 (θ) ≤ 1 for a propagation at a velocity ĉk,R < c. In both cases, there is a dispersive behaviour of the scheme.

We proceed by expanding the normalized eigenvalue s k = cλ k /h in powers of h when h → 0. This yields (the value iξ corresponds to the operator ∂ x ): Otherwise stated, the modified equation [START_REF] Yu | The method of differential approximations[END_REF] of the scheme (4.15) is

s k = -c
∂ t u + c∂ x u = -c 17 384 h 4 ∂ (5)
x u(x, t) -

229 3840 h 5 ∂ (6) x u(x, t) + O(h 6 ). (4.22)
The relation (4.22) indicates fourth order accuracy with a leading dispersive error of order 4. Furthermore the first dissipative term is of order 5 with a real part of λ k behaving as: Remark 4.3. Not all cubic reconstruction lead to a stable dynamical system (4.9). Consider for example the consistent approximations σα , θα and ψα to u (x α ), u (x α ) and u (x α ), defined by: An algebraic calculation, not shown here, demonstrates that the spectrum is located on the right of the imaginary axis and therefore that the scheme (4.7) with the reconstruction w α [ū] in (2.36) with σα , θα and ψα given in (4.24) is instable.

Re(λ k ) ∼ - 229 3840 φ 6 k ≤ 0. ( 4
                 ũα + h 2 24 θα =

One dimensional numerical results

In this section, we display numerical results obtained with the scheme (4.2). First we give some results on the linear convection equation. In this case, (4.2) becomes (4.7). Then we show results for nonlinear gasdynamics test cases. In all the results the RK4 time stepping scheme is used. If V n V (t n ) is supposed given, then V n+1 is computed by (in the case of (4.7)):

                 k 0 = JV n , k 1 = J(V n + 1 2 ∆t k 0 ), k 2 = J(V n + 1 2 ∆t k 1 ), k 3 = J(V n + ∆t k 2 ), V n+1 = V n + ∆t 1 6 k 0 + 1 3 k 1 + 1 3 k 2 + 1 6 k 3 .
(5.1)

The scheme (5.1) is separately fourth order in space and time, and therefore it is globally fourth order (assuming a CFL condition for the time step.).

Accuracy and convergence analysis for the linear advection equation

In this section, we present an accuracy and convergence analysis for the linear convection equation

∂ t u + c∂ x u = 0, c > 0, x ∈ [0, 1], (5.2)
approximated by the scheme (4.7)-(5.1). Fourth order accuracy is expected. The grid refinement procedure is as follows. We start from a mildly irregular grid of

Ω = [0, 1], with N cells T α , |T α | = κ α h, α κ α = 1 defined by the law κ α = K 1 + 1 2 cos 2 (8x α π) . (5.3)
The constant K is adjusted so that

N α=1 κ α = 1. (5.4)
After performing a first computation using this first grid, the grid is refined by randomly splitting each cell into two subcells with shapes corresponding to the two values

κ 2α-1 = K 2 κ α (1 -Cr α ), κ 2α = K 2 κ α (1 + Cr α ), (5.5)
where C is a constant monitoring the randomness of the refinement. The grid is refined five times according to (5.5). The result is a highly irregular grid. Fig. 5.1 reports the convergence of the scheme in maximum norm after 1, 5 and 10 periods of propagation and using two different choices of the constant C. We observe very satisfactory error levels and convergence rates matching the expectation of fourth order accuracy. Note that there is little difference between the two refinement patterns.

Shock tube problems

In this section we display several results obtained for shock tube problems of gas dynamics. The one dimensional Euler system has the form (1.1)

∂ t q + ∂ x f = 0, (5.6) with      q = [ρ, ρu, ρ(E + 1 2 u 2 )] T , f = ρu, ρu 2 + p, ρ(E + 1 2 u 2 ) + p u T .
(5.7) The semi discrete fourth order finite volume scheme for t →

V (t) = [v 1 (t), v 2 (t), . . . , v N (t)] T is dv α (t) dt = - 1 |T α | f α+1/2 (t) -f α-1/2 (t) . (5.9)
where v α (t) q α (t), the average of q(x, t) over T α . The numerical flux function is chosen as the HLLC numerical flux [START_REF] Toro | Riemann solvers and numerical methods for fluid dynamics[END_REF]. But other numerical flux can be selected without significant influence on the result. The cubic reconstruction is performed on the physical variables using:

w[V ] α (x) = v α + σα (x -x α ) + θα 1 2 (x -x α ) 2 - 1 24 |T α | 2 + 1 6 ψα (x -x α ) 3 .
(5.10)

where the vector values σα , θα and ψα , are obtained by the Algorithm 1 in Sec. 3. Again, the timestepping scheme is the RK4 scheme. We show the results obtained for two shock tube problems obtained with order the order two and four versions of the scheme, respectively. In both cases, we use the simple limiter (5.11) acting on the the physical variables. It is designed as follows. In each cell T α , the reconstructed poynomial w α [V ](x, t 0 ) is limited so that

w α [V ](x α±1/2 ) ∈ [ min β∈Vα v β , max β∈Vα v β ]. (5.11)
This is obtained by applying the same scalar coefficient λ α to the three values σα , θα and ψα in (2.36).

In other words, the polynomial w α [V ](x) in (5.10) is replaced by

wα (x) = ūα + λ α w α [V ](x) -ūα (5.12)
with λ α > 0 such that (5.11) holds in the cell T α . This limitation is part of the semi-discrete approximation. Therefore it is applied at each substep of the time scheme (5.1). Since this paper is focused on the reconstruction in Algorithm 1, we do not elaborate further on the limiter design and we limit ourselves to the limiter (5.11). Note that this limiter does not a priori prevent new extrema to appear. It is used for the sole purpose to limit nonlinear oscillations. Note also that in such shock tube problems, the physical time of simulation is too short to observe significant differences between the second order and the fourth order schemes. Future work will be devoted to a sharper comparison with a larger physical time of observation. Here a regular grid is used and our purpose is only to show that the fourth order scheme (4.2) is sufficiently robust to handle these gas dynamics test cases. The CFL = 0.3 is selected in order for the time approximation to not significantly modify the accuracy in space.

Sod test case

In all our tests, the dimensions of the tube are:

x 0 = 0.5 m location of the initial discontinuity, x l = 5.0 m lenght of the half of the tube.

(5.13)

For the Sod test, we use the initial left and right states in the unit system of (5. The observation time is t = 0.003s. The results are plotted against the exact solution. The contact discontinuity is slightly better with the scheme of order 4 (right column of Fig. 5.2) than with the one of order 2 (left column of Fig. 5.2). As already mentioned, due to the short observation time, there is no significant difference between the second order and the fourth order versions of the scheme. Observe that the limiter (5.11) (applied to the physical quantities q = [p, T, u] T ) is sufficient to ensure a good behaviour of the scheme near discontinuities.

Lax test case

The Lax test case initial conditions are given using the physical system of units (5. This corresponds to left and right temperature values of T l = 2082K and T r = 300K, respectively. The observation time is t = 0.00575s. Again there is no significant difference between the versions of order 2 (Fig. 5.3) and 4 (right column of Fig. 5.3). We only observe a better resolution of the contact discontinuity with the fourth order scheme. In this case, the limiter (5.11) is not sufficient to limit the overshoot near the contact discontinuity. We oberve this behaviour in the case of a temperature transition of the form T before < T after . Note that no overshoot was observed in the opposite case (T after < T before ). Observe nonetheless that the fourth order scheme is stable.

Multidimensional numerical results

The two test cases displayed hereafter are related to propagation problems in gas dynamics. Both aim to illustrate the efficiency of the cubic reconstruction procedure on an irregular grid when compared to lower order reconstructions.

Vortex transport by a uniform flow

The first test case consists of an isothermal Gaussian vortex moving accross a periodic box with constant velocity. This vortex is an analytical solution of the two dimensional Euler equations. The dimensional characteristics of the vortex are as follows:

                                            
gas constant: R = 8.314 J mol -1 K -1 , molar mass: M g = 28.013410 3 kg mol -1 , σ = 0.001, a = 1 e -r 2 2σ 2 , advection velocity: v x = 50 ms -1 , v y = 20 ms -1 , t start = 0, t end = 0.004, 10 periods along the x direction, 4 periods along the y direction, ∆t = 5.10 -8 s, CFL on finest grid 0.2. (

A piecewise cubic reconstruction of the form (2.28) is used combined with the Roe numerical flux. There is no limitation needed for this case. At final time Fig. 6.2 reports the shape of the vortex when using the irregular grid in the right panel of Fig. 6.1. The vortex spans over approximately 15 cells. As can be observed, after several periods of propagation, there is a drastic difference between the piecewise affine reconstruction and the piecewise cubic reconstruction. For comparison the result using a 200 × 200 Cartesian grid is also reported. Comparing Fig. 6.2 and 6.3, it turns out that the result obtained with the fourth order scheme on the irregular grid is more accurate than the result obtained using the Cartesian grid with the second order scheme.

Acoustic Wave in a 3D Channel

This test case concerns a propagation problem typically occuring in aeroacoustics. Aeroacoustics computations are usually performed by associating two kinds of codes:

• an unsteady CFD code for LES (Large Eddy Simulation) . This first code calculates acoustic sources generated by the turbulent flow and performs the numerical propagation in the acoustic near field region,

• an acoustic code to reconstruct the far field from unsteady CFD results recorded at CFD-acoustic interfaces (e.g. [START_REF] Rahier | Investigation of integral surface formulations for acoustic post-processing of unsteady aerodynamic jet simulations[END_REF]).

It is thus very important to assert the efficiency of the CFD code for accurately propagating acoustic waves in the near field up to the CFD/acoustic interface without excessive damping or phase errors.

Complex geometries are to be considered frequently and unstructured meshes are most of the time mandatory. In general, when 2 nd order CFD schemes are used, fine meshes are required in the near field region to compensate the limited frequency bandwidth of these schemes. In [START_REF] De La Puente | On LAGOON nose landing gear CFD/CAA computation over unstructured mesh using a ZDES approach[END_REF][START_REF] Lorteau | Recent progress in les computation for aeroacoustics of turbulent hot jet. Comparison to experiments and near field analysis[END_REF] numerical results are shown with a 2 nd order scheme for aeroacoustics problems. It is believed that higher order accurate schemes will demonstrate higher bandwidth thus allowing less refined meshes for a fixed precision. A popular measure of the bandwidth of a scheme is simply expressed as the minimum number of grid points (or cells) per wavelength (denoted by ppw ). This is required to properly propagate a given wave pattern over a significant length [START_REF] Bogey | A family of low dispersive and low dissipative explicit schemes for flow noise and noise computations[END_REF]. When using a 2 nd order schemes, a ppw of the order of 20-25 is often mentioned while with recent high order schemes, the ppw parameter is likely of the order 4-6.

Test case outline

Our test case consists of an acoustic wave propagating in a uniform flow. The geometry consists of a square duct with axis Ox and uniform cross section of 0.2 × 0.2m. The length is L=10 m. The duct is filled with a uniform flow of air at V x =1 m•s -1 at ambient conditions: 101 325 Pa and 300 K. Under these conditions the sound velocity is c = 347 m•s -1 . Boundary conditions are set to:

• left edge (x = 0): subsonic inlet conditions set at V x =1 m•s -1 , T =300 K.

• right edge (x = L): subsonic outlet conditions set at P=101 325 Pa.

• lateral surfaces : slip boundary conditions.

We have used a uniform mesh of isotropic tetrahedra of size a =20 mm. At the right edge, four layers of prisms of uniform size of 16 mm were added. This part of the grid intends to allow for a planar acoustic wave generation as described below. Overall, the grid consists of 288876 cells and 600053 faces, (see fig.

6.4).

The simulation is initiated by a sinusoidal perturbation at the right end of the duct: a time harmonic, isothermal (300 K), distributed mass source enters the first layer of prismatic cells on the right of the computational domain. It is expressed as:

ṁ = A sin (2πf t) . (6.2)
In (6.2 ) the amplitude A is set so as to generate a low amplitude acoustic wave (smaller than 20 Pa of physical amplitude) such that non linear effects can be ignored. This source generates a plane acoustic wave propagating upstream, from the right to the left with a phase velocity

V ϕ = c -V x = 346 m•s -1 .
The corresponding wavelength is thus: The physical elapsed time of each run is 20 ms. This corresponds to a propagation length of 6.9 m in the counter flow direction (x < 0). The inlet boundary is therefore not impacted, thus alleviating any acoustic reflection issues at the left end boundary. A total of 7578 time steps are performed at CFL = 0.4. By using a series of values for the source frequency f , the resolution of the acoustic wave for a given grid, as represented by the parameter ppw (number of grid points per wavelength) is estimated by ppw = λ a . (6.4)

λ = V ϕ f . ( 6
Runs were performed for a series of 4 values of the frequency f, which corresponds to a ppw number varying from 6 (coarse mesh) to 59 (extra fine mesh). CFD schemes are used in aeroacoustic, a value of 20-25 ppw is typically used. This corresponds to a specific grid step size, according to the expected frequency range of the acoustic sources. This value was selected for the runs of the R2 series. Then for each run conditions, computations were performed for schemes of various orders, as described in Table 2, and the wave quality was observed and compared at different numerical transducers located along the duct axis. 1 and2) 6.2.2. Amplitude and phase error analysis The analysis of the numerical results was performed using three numerical transducers located on the duct axis at 7 m, 8 m and 9 m, respectively, see Fig. 6.5 top. For example we show on Fig. 6.5 (bottom) the pressure time history for a specific run.

The quality of the schemes was estimated from the wave amplitude damping along the propagation direction. The phase error was measured with respect to an ideal linear acoustic wave propagation solution of the form P (x, t i ) = Re P 0 e i(ωt-kr+ϕ0) (6.5)

with ω = 2πf ∈ R and k = k r + ik i ∈ C. The calculated time series is assumed to have the form

F (x, t i ) = A 0 (x) sin (2πf (x) t i + ϕ (x)) e -α(x)ti . (6.6)
where the four parameter functions A 0 (x), f (x), ϕ(x) and α(x) have to be estimated. Identifying these four functions is obtained by a fitting approximation technique using exponential (in time) sums. This kind of fitting is thoroughly presented e.g. in [START_REF] Wiscombe | Exponential-sum fitting of radiative transmission functions[END_REF]. Here we assume that our numerical data can be represented by the series

F (x, t i ) = N n=1
a n e iσnti with (a n , σ n ) ∈ Z 2 (6.7)

The function F represents a real monochromatic wave, so that we limit the sum to two terms (N = 2 in (6.7)). By identification the functions A 0 (x), f (x), ϕ(x) and α(x) are constants and they are expressed in terms the values a 1 , a 2, σ 1 and σ 2 as

A 0 = √ 4a 1 a 2 , 2πf = σ 1 -σ 2 2 , ϕ = 1 2i log -a 1 a 2 , α = σ 1 + σ 2 2i . (6.8)
An iterative least squares procedure is then used to identify these four constants. Turning back to the damping and dispersion evaluation, we use the now estimated functions A 0 (x) and ϕ(x) as a basis of our measurement 4 . The numerical wave damping for one single wavelength is estimated by

A 0 (x 2 ) A 0 (x 1 ) = e k i (x2-x1) = e λk i ( x 2 -x 1 λ
) ,

e λk i = e λ x 2 -x 1 log A 0 (x 2 )
A 0 (x 1 ) , (6.9) where e λk i represents the wave damping per wavelength traveled. Similarly, the numerical dispersion can be represented by the phase difference between two transducer locations as

∆ϕ 21 = ϕ (x 2 ) -ϕ (x 1 ) , (∆ϕ 21 ) th = -k r (x 2 -x 1 ) , (6.10)
where k r = 2πf

Vϕ . The relative error on the phase velocity error can be finally expressed as

∂V ϕ V ϕ = ∂ (∆ϕ 21 ) 2π (x 2 -x 1 ) λ . (6.11)

Numerical Results

Results using reconstruction of order 2, 3 and 4 (denoted by O2, O3, and O4) are reported. This corresponds to reconstruction polynomials of order 1, 2 and 3 respectively. The comparison analysis of Sec 6.2.2 was carried out with the results obtained from different runs summarized in Table 2.

Run R1 R2 R3 R4 O2 -RK2 x x x x O3 -RK3 x x x O4 -RK4 x x x O4 -RK4, CFL = 0.8
x Table 2: Test Case 6.2: acoustic wave in a channel -Computation matrix Tables 3 and4, and Figs. 6.6 and 6.7 show that the scheme of order 4 provides a drastic improvement in the ability of the code to propagate acoustic waves. Even on a coarse grid (up to 6 points per wavelength) a very small damping and phase error was observed. Regarding a detailed comparison, Fig. 6.6 shows a similar damping between 2 nd order and the 3 rd order scheme whereas Fig. 6.7 shows that the 3 rd order and the 4 th order schemes have a similar phase error. This is easily interpreted in terms of the modified equation (see (4.22) and [START_REF] Yu | The method of differential approximations[END_REF]) associated with the linear transport equation. A 2 nd and a 3 rd order scheme have a first dissipative error h 3 ∂ (4)

x u, whereas a 3 rd and 4 th order schemes have a Run (ppw) O2-RK2 O3-RK3 O4-RK4 O4-RK4(CFL = 0.8) R1 ( 59 -0.0695% -0.0054% 0.0000% R3 [START_REF] Haider | Discrétisation en maillage non structuré général -Applications à la LES[END_REF] -0.3046% -0.0183% -0.0256% -0.0287% R4 ( 6) -1.0516% -0.0466% -0.0744%

x u. Overall the 4 th order scheme detailed in Sections 3 and 4 outperforms the two other schemes.

Table 5 reports numerical values regarding computational cost. The computational complexity of the 3 rd and 4 th schemes was numerically measured against the 2 nd order scheme. Two values are reported. First the intrinsic CPU ratio, which is a numerical evaluation of the cost of one single time iteration. Second the CPU ratio for the full simulation. The latter also depends on the stability properties of each scheme. As an example, the RK4 (explicit) time stepping scheme allows a stability limit of CFL = 0.8 with the 4 th order scheme instead of CFL = 0.4 with the second order scheme. The overall CPU ratio is therefore reduced to a factor, as reported in Table 5. The change from CFL = 0.4 to CFL = 0.8 did not bear any loss in terms of scheme accuracy, as shown in Tables 3 and4 and in Fig. 6.6 and 6.7.

Conclusion

This paper describes the Coupled Least Squares reconstruction procedure to be used in finite volume methods. This piecewide cubic reconstruction belongs to the k-exact reconstruction of [START_REF] Barth | Higher order solution of the euler equation on unstructured grids using quadratic reconstruction[END_REF]. We have investigated the design and implementation of this reconstruction on irregular grids. We have obtained a high order finite volume scheme which can be extended to complex physical models. The resulting scheme exhibits good stability properties.

As explained in the introduction, our main interest is to keep the paradigm of one unknown per cell, avoiding to introduce additional degrees of freedom in the cells. Furthermore, our reconstruction procedure is fully centered. In particular no solution dependent reconstruction stencils are used. Since the algorithm needs data located in a compact neighborhood, the reconstruction can be efficiently implemented regarding parallellism. This is clearly of crucial importance for practical applications. This property will be further analyzed elsewhere in physically realistic situations. The tensor formalism permits to handle the algebra for polynomial reconstruction in multidimension. Refer on this topic to the recent monograph [START_REF] Hackbusch | Tensor Spaces and Numerical Tensor Calculus[END_REF]. Here we review the definitions needed in Section 2.1. We follow [START_REF] Haider | Discrétisation en maillage non structuré général -Applications à la LES[END_REF][START_REF] Haider | Efficient implementation of high order reconstruction in finite volume methods[END_REF]. A tensor a of order k is a linear application

(x 1 , . . . , x k ) ∈ (R d ) k → a (k) (x 1 , . . . , x k ). (8.1)
The components of a

(k) are a i1...i k = a (k) (e i1 , . . . e i k ). (8.2)
The tensor a (k) is symmetric if for all permutation σ of k elements,

a σ(i1)...σ(i k ) = a i1...i k . (8.3)
In Section 2.1 is used the fact that the derivative D (k) u(x) is a symmetric tensor of order k with components

D (k) u i1...i k = ∂u(x) ∂x i1 . . . ∂x i k , 1 ≤ i 1 . . . i k ≤ d (8.4)
Finally for 0 ≤ m ≤ k, the product of the tensors a (k) and b (m) is the tensor of order k -m denoted by a (k) . b (m) with components

(a (k) .b (m) ) i1...i k-m = d j1=1 • • • d jm=1 a i1...i k-m j1...jm b j1...jm . (8.5)
In the particular case m = k, the product of a k and b k is called the contraction and is the scalar denoted by

a (k) • b (k) = d j1=1 • • • d j k =1 a j1...j k b j1...j k . (8.6)
Finally the tensor product of a (k) and b (m) is the tensor c

(k+m) = a (k) ⊗ b (m) with components c i1,...i k+m = a i1...i k b i k+1 ...i k+m . (8.7)
In particular, for x ∈ R d a given vector, the tensor x ⊗k is defined by

x ⊗k i1...i k = x i1 . . . x i k . (8.8)
In our context, specific tensor notation is as follows. Cell indices are α and β ∈ V α .

                                       x α = (x α,1 , . . . , x α,d ) T ∈ R d barycenter of T α , h αβ = x β -x α , with norm h αβ = |h αβ |, h (k)
αβ = h αβ ⊗ • • • ⊗ h αβ k× , z (k) αβ = 1 |T β | T β (x -x α ) ⊗k dx, x (k) α = z (k) αα , |T α | = Tα dx = x (0) α . (8.9)
In particular we have z

(1) α -x (1) α = h (1) αβ = h αβ . (8.10)
The components of the tensor z

(k) αβ are z (k) αβ,i1i2...i k defined for 1 ≤ i 1 . . . i k ≤ d by z (k) αβ,i1i2...i k = T β (x i1 -x α,i1 )(x i2 -x α,i2 ) . . . (x i k -x α,i k )dx. (8.11)
Moreover the consistency of the reconstruction is estimated using the scale h which is such that

ch d ≤ |T α | ≤ Ch d , ch ≤ h αβ ≤ Ch, 1 ≤ α, β ≤ N, (8.12)
for some constants c and C.

Appendix: Spatial accuracy analysis

Here we give some comments on the accuracy analysis of the approximation in space developed in this paper. Let u(x, t) be solution of the advection equation ∂ t u = ∂ x u. Averaging over the cell T α gives (see Sec. [START_REF] Barth | Higher order solution of the euler equation on unstructured grids using quadratic reconstruction[END_REF] Then H α (u) becomes:

H α (u) = 1 h α L - α+1/2 (u) -L - α-1/2 (u) . (9.5)
In the particular case of a regular grid (h α = h for all α), the operator L - α+1/2 is given by (see (4.15)): H α (u) = -17 384 h 4 ∂ (5) x u(x α ) + 229 3840 h 5 ∂ (6) x u(x α ) + O(h 6 ). (9.8) This proves Prop. 4.1, i.e. the fourth order accuracy of the scheme. When h α is non constant the evaluation of the order p > such that H α (u) = O(h p ) is more delicate. In particular, it cannot be performed by a Fourier symbol calculation. We show next how the Peano Kernel Theorem ( [START_REF] Powell | Approximation theory and methods[END_REF], chap. 22, p.270) can provide a set of sufficient conditions on the irregularity of the grid in order for (9.2) to keep a high order accuracy. For simplicity of the presentation, the analysis is restricted to the case of the piecewise linear reconstruction defined by

L - α+1/
w α [V ](x) = v α + σ α [V ](x -x α ) (9.9)
where σ α is the least squares slope [START_REF] Haider | Discrétisation en maillage non structuré général -Applications à la LES[END_REF][START_REF] Haider | Stability analysis of the cell centered finite-volume MUSCL method on unstructured grids[END_REF] given by:

σ α = h α,α+1 h 2 α,α-1 + h 2 α,α+1 (v α+1 -v α ) + h α,α-1 h 2 α,α-1 + h 2 α,α+1
(v α-1 -v α ).

(9.10)

The same analysis can be extended to the cubic reconstruction case. Instead of (9.6) the operator L - α+1/2 (u) is in the case of a linear reconstruction:

L - α+1/2 (u) = u(x α+1/2 ) -ūα + 1 2 h α σα . (9.11)
Since L - α+1/2 (p) = 0 when p(x) is any polynomial of degree 1, the Peano Kernel Theorem implies that L - α+1/2 can be expressed as

L - α+1/2 (u) = x α+1/2
x α-1/2 K - α+1/2 (θ)u (θ)dθ.

(9.12)

The kernel

K - α+1/2 (θ) is K - α+1/2 (θ) = 1 2 L - α+1/2,x {(x -θ) + }, (9.13)
where the operator L - α+1/2 is applied in the x-variable. A simple calculation shows that K - α+1/2 (θ) is expressed as the second order polynomial

K - α+1/2 (θ) = a α (x α+1/2 -θ) 2 + b α (x α+1/2 -θ) + c α , (9.14)
with coefficients (we denote h α+1/2 = (h α + h α+1 )/2): Consider now that α is a continuous parameter such that α → x(α) describes the variation of the grid size. The grid size h α = h(α) is defined by h(α) = x (α) with the additional assumption that h(α) = O(h). (9.16)

                 a α = - 1 4h α + 1 8
h
The relations (9.15) become:

                 a(α) = - 1 4h(α) + 1 8
h(α + 1/2) -h(α -1/2) h 2 (α + 1/2) + h 2 (α -1/2) , b(α) = 1 2 - 1 4
h(α)h(α + 1/2) h 2 (α + 1/2) + h 2 (α -1/2) , c(α) = -1 8

h(α)h(α + 1/2)h(α + 1) h 2 (α + 1/2) + h 2 (α -1/2) .

(9.17)

Using a change of variable, H α (u) in (9.5) is expressed as:

H α (u) = 1 h(α) 1 0 a(α)h(α) 3 τ 2 + b(α)h(α) 2 τ + c(α)h(α) u x(α + 1/2) -τ h(α) -a(α -1)h(α -1) 3 τ 2 + b(α)h(α -1) 2 τ + c(α)h(α -1) u x(α -1/2) -τ h(α -1) dτ. Therefore ∂ σ A = 0 and A = O(h 2 ). Due to the fact that h (σ) = 0, Prop 9.1 proves the second order accuracy of the scheme (9.2) and therefore we recover (9.8.)

Remark 9.3. Proposition 9.1 only claims that the accuracy is preserved when the grid is a slight perturbation of a regular grid in the sense of (9.23). The magnitude of the perturbation is measured by the variation of the function A(σ, τ ). Therefore the problem of the accuracy analysis of the scheme (9.2) on grids with cells presenting large shape ratio remains an open question.

Figure 2 . 1 :

 21 Figure 2.1: Two cells of a general polyhedral grid in the two dimensional case, (d = 2

( 1 )

 1 α only. Let us start by approximating the derivative σ α by a linear relation of the form σ α β∈Wα c αβ (ū β -ūα ).

  θα = θ α + O(h 2 ), (b) σα = σ α + O(h 3 ). (c) (3.16) Proof. The consistency relations (3.16) result directly from comparing (3.15) and (3.2) and from the fact that a α = O(h), ãα = O(h), āα = O(h) and bα = O(h 2 ).

4 .

 4 Using (3.15) (b) , compute the second order derivative θα for all α. 5. Using (3.15) (c) , compute the first order derivative σα for all α.

Figure 4 . 1 :

 41 Figure 4.1: Left panel: spectral abscissa Λ (4.10) for the scheme (4.6) for a sample of 5000 random irregular grids of size 64. The scale along y is 10 -13 . Right panel: distribution of the spectral absissa Λ for a sample of 5000 random grids of size 64: a Gaussian distribution around the value 0 is clearly observed. This indicates a stable scheme.

. 23 )

 23 This suggests a stable behaviour. Actually the normalized full spectrum (see bottom panel of Fig. 4.2) confirms the stability of the scheme (4.15). The dissipation curve φ k ∈ [0, π) → D 1 (φ k ) and the dispersion curve φ k ∈ [0, π) → D 2 (φ k ) are shown on Fig. 4.2 (top panels).

Figure 4 . 2 :

 42 Figure 4.2: Top left panel: dissipation of the scheme (4.15). Top right panel: dispersion of the scheme (4.15). Bottom panel: spectrum of the matrix J in (4.7). The spectrum is located in the left half-plane. This indicates stability. The order of accuracy is reflected by the tangency order at the origin along the vertical axis.

( a )Figure 5 . 1 :

 a51 Figure5.1: Fourth order finite volume scheme for the linear convection equation on a set of consistently refined irregular grids. The shape aspect of the cells is given by (5.3). Left panel: each cell is refined in two subcells. Right panel: random effect in the grid refinement, given in (5.5). In both case, an almost fourth order is observed after 1, 5 and 10 periods.

  3.5 kg . m -3 u l = 0.0 m . s -1 p l = 10 6 Pa    ρ r = 0.4375 kg . m -3 u r = 0.0 m . s -1 p r = 10 5 Pa (5.14) This corresponds to left and right temperature values of T l = 1000K and T r = 800K, respectively.

ρ l = 1 .

 1 038569 kg . m -3 u l = 191.511 m . s -1 p l = 6.17800 10 5 Pa    ρ r = 1.16667 kg . m -3 u r = 0.0 m . s -1 p r = 10 5 Pa (5.15)

Figure 5 . 2 :

 52 Figure 5.2: Sod shock tube . Finite volume scheme with a regular grid of 200 cells. Left column: standard second order reconstruction. Right column: finite volume scheme with fourth order reconstruction (5.10). In both cases, the slope limiter is (5.11) and CFL = 0.3. The exact solution is represented with the continuous line.

Figure 5 . 3 :

 53 Figure 5.3: Lax shock tube . Finite volume scheme with a regular grid of 200 cells. Left column: standard second order reconstruction. Right column: finite volume scheme with fourth order reconstruction (5.10). In both cases, the slope limiter is (5.11). The CFL is CFL = 0.3. The exact solution is represented with the continuous line.

Figure 6 . 1 :

 61 Figure 6.1: Convection of the isothermal vortex (6.1) -Left panel: initial condition. Right panel: irregular grid.

  constant temperature: T = T 0 = 300 K, Gaussian pressure distribution: p = p 0 (1 -ae -r 2 2σ 2 ) Pa, p 0 = 10 5 Pa , (gives p min 99985.96 Pa), vortex angular velocity: v θ =

Figure 6 . 2 :

 62 Figure 6.2: Convection of the isothermal vortex (6.1) with the unstructured grid in Fig. 6.1 at final time -Left panel: piecewise linear reconstruction (second order FV scheme) -Right panel: piecewise cubic reconstruction , (fourth order FV scheme).

Figure 6 . 3 :

 63 Figure 6.3: Convection of the isothermal vortex (6.1) with a Cartesian fine grid at final time -Left panel: piecewise linear reconstruction (second order FV scheme) -Right panel: piecewise cubic reconstruction, (fourth order FV scheme).

. 3 )Figure 6 . 4 :

 364 Figure 6.4: Test Case 1: Acoustic wave in a channel -Duct overview and grid details at the right end

Figure 6 . 5 :

 65 Figure 6.5: Test case 1: Acoustic wave in a channel -Top: illustration of wave propagation -Bottom: pressure time histories at numerical transducers location. Numerical case R3 and scheme O4-RK4 (see Tables1 and 2)

Figure 6 . 6 :Figure 6 . 7 :

 6667 Figure 6.6: Test Case 6.2: acoustic wave in a channel -Wave damping per wavelength

( 9 . 18 )Proposition 9 . 1 .( 9 . 23 ) 9 . 2 .

 9189192392 Define (σ, τ ) ∈ [α -1, α] × [0, 1] → A(σ, τ ) by A(σ, τ ) = a(σ)h(σ) 3 τ 2 + b(σ)h(σ) 2 τ + c(σ)h(σ),(9.19)where, for σ ∈ [α -1, α], the values a(σ), b(σ), c(σ) are given in (9.17) when replacing h α by h(σ) in the right-hand side. Then H α (u) is rewritten asH α (u) = 1 h(α) [α-1,α]×[0,1] ∂ ∂σ A(σ, τ )u x(σ + 1/2) -τ h(σ) dσdτ. (9.20)Therefore by the mean value theorem for integrals, |H α (u)| is estimated by|H α (u)| ≤ 1 h(α) [α-1,α]×[0,1] ∂ ∂σ [A(σ, τ )u x(σ + 1/2) -τ h(σ) ] dτ dσ.(9.21)Finally we have∂ ∂σ A(σ, τ )u (x(σ + 1/2) -τ h(σ)) = ∂A(σ, τ ) ∂σ u x(σ + 1/2) -τ h(σ)(9.22)+ A(σ, τ ) h(σ + 1/2) -τ h (σ) u x(σ + 1/2) -τ h(σ) .This estimate leads to the following result: The scheme (9.2) is second order accurate withy respect to (9.1) under the sufficient condition on the variation of the grid α → x(α), given by the relations A(σ, τ)| ≤ Ch 3 , max σ,τ |A(σ, τ )| |h(σ)| + |h (σ)| ≤ Ch 3 .Remark In the case of a regular grid of size h, the function h(α) is the constant h(α) = h. We haveA(σ, τ ) = h 2

  The outcome of the preceding analysis is the following algorithm to calculate ψα , θα and σα . Suppose given a function u(x) with averages in the cells T α given by ū = [ū 1 , ū2 , . . . , ūN ]. The relations(3.15) translate in the following computational procedure:

	Algorithm 1.	1. Using (2.39), compute σLS α in terms of the averages ūα for all α.
	2. Compute	cαβ (σ LS β -σLS
		β∈Wα

α ) with cαβ given in

(3.4)

, for all α.

3. Using (3.15) (a) , compute the approximate third order derivative ψα for all α.

  ūα ,

	1 6 σα-1 + 2 3 σα + 1 6 σα+1 = δ x ũα ,	(4.24)
	1 12 θα-1 + 10 12 θα + 1 12 θα+1 = 2δ 2 x ũα -δ x σα,
	1 12 ψα-1 + 10 12 ψα + 1 12 ψα+1 = δ 2 x σα .

Table 1 :

 1 Table 1 summarizes the results. When 2 nd order Test Case 1: Acoustic wave in a channel -Values of the wave frequency

	Run	R1	R2	R3	R4
	f (Hz)	300	700	1500	3000
	λ (m)	1.1533 0.4943 0.2307 0.1153
	ppw = λ/a	59	25	12	6

Table 3 :

 3 Test Case 6.2: acoustic wave in a channel -Summary of results for wave damping per wavelength

	)	0.9995
	R2 (25)	0.9961	0.9965	0.9990
	R3 (12)	0.9843	0.9866	0.9942	0.9939
	R4 (6)	0.9188	0.9276	0.9754
	Run (ppw) O2-RK2	O3-RK3	O4-RK4 O4-RK4(CFL = 0.8)
	R1 (59)	-0.0119%
	R2 (25)

Table 4 :

 4 Test Case 6.2: acoustic wave in a channel -Summary of results for the relative phase velocity error(6.11)

	a	O2-RK2 O3-RK3 O4-RK4 O4-RK4(CFL = 0.8)
	Intrinsic CPU ratio	1	4.7	12.0	12.0
	Run CPU ratio	1	4.7	12.0	6.0

Table 5 :

 5 Test Case 6.2: acoustic wave in a channel -Relative CPU cost first dispersive error of the form h 4 ∂

  , t) -u(x α-1/2 , t) , with h α = |T α |. (9.1) The approximation t → v α (t) of ūα (t) is solution of: (t) = [v 1 (t), . . . , v N (t)] and w α [V ](x)is the cubic polynomial (2.36). Let u(x) be any regular function and U = [ū 1 , . . . , ūN ]. The accuracy analysis consists in evaluating p > 0 such that H α (u) = O(h p ) where

	.1):
		dū α (t) dt u(x α+1/2 dv α (t) = 1 h α dt = 1 h α w α [V (t)](x -α+1/2 ) -w α [V (t)](x -α-1/2 ) .	(9.2)
	where V H α (u) =	1 h α	u(x α+1/2 ) -u(x α-1/2 ) -	1 h α	w α [U ](x -α+1/2 ) -w α [U ](x -α-1/2 ) .	(9.3)
	Define the operator L -α+1/2 (u) by
			L -α+1/2 (u) = u(x α+1/2 ) -w α [U ](x -α+1/2 ).	(9.4)

  2 (u) = u(x α+1/2 ) -ūα +The Taylor expansion of H α (u) at x α is (see (4.21)):

		1 2	hσ α +	1 12	h 2	θα +	1 48	h 3	ψα ,	(9.6)
	where the notation is (see (4.12-4.13)):
	                        	σLS α = ψα = (2h) 2 ūα+1 -ūα-1 2h σLS α+2 + σLS , α-2 -2σ LS α 1 σα = 53 48 σLS α -5 96 σLS α+2 -5 σLS α-2 , 96 θα = 2h σLS α+1 -σLS α-1 , 1	.			(9.7)

CEDRE is being developed at Onera, France, http://cedre.onera.fr

We could as well consider a priori the five points dependence stencil Wα = {T α-2 , T α-1 , Tα, T α+1 , T α+2 }

Note that the two other functions f (x) and α(x) are used only to check the accuracy of the exponential fitting.

Acknowledgement: The authors warmly thank F. Vuillot for his help and encouragement which was decisive for the final form of this paper. They also thank the referees for their sharp reading. Their suggestions greatly helped to enhance the presentation of this study.

LaTeX Source Files Click here to download LaTeX Source Files: paper-final-Sep-06-2018.tar.gz


Samy Skander Bahoura
email: samybahoura@yahoo.fr

A compactness result for solutions to an equation with boundary singularity

Keywords: blow-up, boundary, logarithmic singularity, a priori estimate, analytic domain, Lipschitz condition. MSC: 35J60, 35J05, 35B44, 35B45

We give a blow-up behavior for solutions to a problem with boundary singularity and Dirichlet condition. An application, we derive a compactness result for solutions to this Problem with singularity and Lipschitz condition.

Introduction and Main Results

We set ∆ = ∂ 11 + ∂ 22 on an analytic domain Ω ⊂ R 2 .

We consider the following equation:

(P )    -∆u = -log |x| 2d V e u in Ω ⊂ R 2 , u = 0 in ∂Ω.
Here:

d = diam(Ω), 0 ∈ ∂Ω,
and,

0 ≤ V ≤ b < +∞, -log |x| 2d e u ∈ L 1 (Ω), u ∈ W 1,1 0 (Ω).
The previous equation was studied by many authors, with or without the boundary condition, also for Riemannian surfaces, see [START_REF] Aubin | Some Nonlinear Problems in Riemannian Geometry[END_REF][START_REF] Bandle | Isoperimetric Inequalities and Applications[END_REF][START_REF] Bartolucci | sup+Cinf" inequality for Liouville-type equations with singular potentials[END_REF][START_REF] Bartolucci | A 'sup+Cinf' inequality for the equation -∆u = V e u /|x| 2α[END_REF][START_REF] Brezis | A sup+inf inequality for some nonlinear elliptic equations involving exponential nonlinearities[END_REF][START_REF] Brezis | Nonlinear elliptic equations with measures revisited[END_REF][START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF][START_REF] Brezis | Semi-linear second-order elliptic equations in L1[END_REF][START_REF] Chen | A priori estimates for solutions to nonlinear elliptic equations[END_REF][START_REF] Chen | A sharp sup+inf inequality for a nonlinear elliptic equation in R 2[END_REF][START_REF] Hofmann | Geometric and transformational properties of Lipschitz domains, Semmes-Kenig-Toro domains, and other classes of finite perimeter domains[END_REF][START_REF] Yy | Blow-up analysis for solutions of -∆u = V e u in dimension two[END_REF][START_REF] Yy | Harnack Type Inequality: the method of moving planes[END_REF][START_REF] Ma | Convergence for a Liouville equation[END_REF][START_REF] Shafrir | A sup+inf inequality for the equation -∆u = V e u[END_REF], we can find some existence and compactness results.

Among other results, we can see in [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF] the following important Theorem, Theorem.(Brezis-Merle [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF]).If (u i ) i and (V i ) i are two sequences of functions relatively to the problem (P ) with, 0 < a ≤ V i ≤ b < +∞, then, for all compact set K of Ω, sup K u i ≤ c = c(a, b, K, Ω).

If we assume V with more regularity, we can have another type of estimates, a sup + inf type inequalities. It was proved by Shafrir see [START_REF] Shafrir | A sup+inf inequality for the equation -∆u = V e u[END_REF], that, if (u i ) i , (V i ) i are two sequences of functions solutions of the previous equation without assumption on the boundary and, 0 < a ≤ V i ≤ b < +∞, then we have the following interior estimate:

C a b sup K u i + inf Ω u i ≤ c = c(a, b, K, Ω).
Now, if we suppose (V i ) i uniformly Lipschitzian with A the Lipschitz constant, then, C(a/b) = 1 and c = c(a, b, A, K, Ω), see [START_REF] Brezis | A sup+inf inequality for some nonlinear elliptic equations involving exponential nonlinearities[END_REF].

Here we give the behavior of the blow-up points on the boundary and a proof of a compactness result of solutions of a Brezis-Merle type Problem with Lipschitz condition.

Here, we write an extenstion of Brezis-Merle Problem (see [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF]) is:

Problem. Suppose that V i → V in C 0 ( Ω), with, 0 ≤ V i . Also, we consider a sequence of solutions (u i ) of (P ) relatively to (V i ) such that, Ω -log |x| 2d e ui dx ≤ C,
is it possible to have:

||u i || L ∞ ≤ C?
Here, we give a caracterization of the behavior of the blow-up points on the boundary and also a proof of the compactness theorem when the prescribed curvature are uniformly Lipschitzian. For the behavior of the blow-up points on the boundary, the following condition is enough,

0 ≤ V i ≤ b, The condition V i → V in C 0 ( Ω) is not necessary.
But for the proof for the Brezis-Merle type problem we assume that:

||∇V i || L ∞ ≤ A.
We have the following caracterization of the behavior of the blow-up points on the boundary. Theorem 1.1 Assume that max Ω u i → +∞, Where (u i ) are solutions of the probleme (P ) with:

0 ≤ V i ≤ b, and Ω -log |x| 2d e ui dx ≤ C, ∀ i,
then; after passing to a subsequence, there is a finction u, there is a number N ∈ N and there are N points x 1 , x 2 , . . . , x N ∈ ∂Ω, such that,

∂ ν u i → ∂ ν u+ N j=1 α j δ xj , α 1 ≥ 4π, α j ≥ 4π weakly in the sens of measure L 1 (∂Ω).
and,

u i → u in C 1 loc ( Ω -{x 1 , . . . , x N }).
In the following theorem, we have a proof for the global a priori estimate which concern the problem (P ).

Theorem 1.2 Assume that (u i ) are solutions of (P ) relatively to (V i ) with the following conditions:

d = diam(Ω), 0 ∈ ∂Ω, 0 ≤ V i ≤ b, ||∇V i || L ∞ ≤ A, and Ω -log |x| 2d e ui ≤ C,
We have,

||u i || L ∞ ≤ c(b, A, C, Ω),

Proof of the theorems

Proof of theorem 1.1:

We have,

u i ∈ W 1,1 0 (Ω)
. By [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF], e ui ∈ L k , ∀k > 1 and by the elliptic estimates:

u i ∈ W 2,k (Ω) ∩ C 1,ǫ ( Ω), ǫ > 0.
We have,

∂Ω ∂ ν u i dσ ≤ C,
Without loss of generality, we can assume that ∂ ν u i ≥ 0. Thus, (using the weak convergence in the space of Radon measures), we have the existence of a positive Radon measure µ such that,

∂Ω ∂ ν u i φdσ → µ(φ), ∀ φ ∈ C 0 (∂Ω).
We take an x 0 ∈ ∂Ω such that, µ(x 0 ) < 4π. Without loss of generality, we can assume that the following curve, B(x 0 , ǫ) ∩ ∂Ω := I ǫ is an interval.(In this case, it is more simple to construct the following test function η ǫ ). We choose a function η ǫ such that,

           η ǫ ≡ 1, on I ǫ , 0 < ǫ < δ/2, η ǫ ≡ 0, outside I 2ǫ , 0 ≤ η ǫ ≤ 1, ||∇η ǫ || L ∞ (I2ǫ) ≤ C 0 (Ω, x 0 ) ǫ .
We take a ηǫ such that,

-∆η ǫ = 0 in Ω ηǫ = η ǫ on ∂Ω.
Remark: We use the following steps in the construction of ηǫ :

We take a cutoff function η 0 in B(0, 2) or B(x 0 , 2):

We set η ǫ (x) = η 0 (|x -x 0 |/ǫ) in the case of the unit disk it is sufficient. We use a chart [f, B 1 (0)] with f (0) = x 0 .
We can take:

µ ǫ (x) = η 0 (x/ǫ) and η ǫ (y) = µ ǫ (f -1 (y)), we extend it by 0 outside f (B 1 (0)). We have f (B 1 (0)) = D 1 (x 0 ), f (B ǫ (0)) = D ǫ (x 0 ) and f (B + ǫ ) = D + ǫ (x 0 ) with f and f -1 smooth diffeomorphism.            η ǫ ≡ 1, on a the connected set J ǫ = f (I ǫ ), 0 < ǫ < δ/2, η ǫ ≡ 0, outside J ′ ǫ = f (I 2ǫ ), 0 ≤ η ǫ ≤ 1, ||∇η ǫ || L ∞ (J ′ ǫ ) ≤ C 0 (Ω, x 0 ) ǫ .
And,

H 1 (J ′ ǫ ) ≤ C 1 H 1 (I 2ǫ ) = C 1 4ǫ, since f is Lipschitz.
Here H 1 is the Hausdorff measure.

We solve the Dirichlet Problem:

∆η ǫ = ∆η ǫ in Ω ⊂ R 2 , ηǫ = 0 in ∂Ω.
and finaly we set ηǫ = -η ǫ + η ǫ . Also, by the maximum principle and the elliptic estimates we have :

||∇η ǫ || L ∞ ≤ C(||η ǫ || L ∞ + ||∇η ǫ || L ∞ + ||∆η ǫ || L ∞ ) ≤ C 1 ǫ 2 , with C 1 depends on Ω.
We use the following estimate, see [START_REF] Brezis | Semi-linear second-order elliptic equations in L1[END_REF],

||∇u i || L q ≤ C q , ∀ i and 1 < q < 2.
We deduce from the last estimate that, (u i ) converge weakly in W 1,q 0 (Ω), almost everywhere to a function u ≥ 0 and Ω -log |x| 2d e u < +∞ (by Fatou lemma). Also, V i weakly converge to a nonnegative function V in L ∞ . The function u is in W 1,q 0 (Ω) solution of :

   -∆u = -log |x| 2d V e u ∈ L 1 (Ω) in Ω u = 0 on ∂Ω,
As in the corollary 1 of Brezis-Merle result, see [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF], we have e ku ∈ L 1 (Ω), k > 1. By the elliptic estimates, we have u ∈ C 1 ( Ω).

We can write,

-∆((u i -u)η ǫ ) = -log |x| 2d (V i e ui -V e u )η ǫ + 2(∇(u i -u) • ∇η ǫ ). (1)
We use the interior esimate of Brezis-Merle, see [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF],

Step 1: Estimate of the integral of the first term of the right hand side of (1).

We use the Green formula between ηǫ and u, we obtain,

Ω -log |x| 2d V e u ηǫ dx = ∂Ω ∂ ν uη ǫ ≤ 4ǫ||∂ ν u|| L ∞ = Cǫ (2)
We have,

   -∆u i = -log |x| 2d V i e ui in Ω u i = 0 on ∂Ω,
We use the Green formula between u i and ηǫ to have:

Ω -log |x| 2d V i e ui ηǫ dx = ∂Ω ∂ ν u i η ǫ dσ → µ(η ǫ ) ≤ µ(I 2ǫ ) ≤ 4π -ǫ 0 , ǫ 0 > 0 (3)
From (2) and (3) we have for all ǫ > 0 there is

i 0 = i 0 (ǫ) such that, for i ≥ i 0 , Ω -log |x| 2d |(V i e ui -V e u )η ǫ |dx ≤ 4π -ǫ 0 + Cǫ (4)
Step 2: Estimate of integral of the second term of the right hand side of (1).

Let Σ ǫ = {x ∈ Ω, d(x, ∂Ω) = ǫ 3 } and Ω ǫ 3 = {x ∈ Ω, d(x, ∂Ω) ≥ ǫ 3 }, ǫ > 0. Then, for ǫ small enough, Σ ǫ is hypersurface. The measure of Ω -Ω ǫ 3 is k 2 ǫ 3 ≤ µ L (Ω -Ω ǫ 3 ) ≤ k 1 ǫ 3 .
Remark: for the unit ball B(0, 1), our new manifold is B(0, 1 -ǫ 3 ).

(Proof of this fact; let's consider d(x, ∂Ω) = d(x, z 0 ), z 0 ∈ ∂Ω, this imply that (d(x, z 0 )) 2 ≤ (d(x, z)) 2 for all z ∈ ∂Ω which it is equivalent to (z -z 0 ) • (2x -z -z 0 )
≤ 0 for all z ∈ ∂Ω, let's consider a chart around z 0 and γ(t) a curve in ∂Ω, we have;

(γ(t) -γ(t 0 ) • (2x -γ(t) -γ(t 0 )) ≤ 0 if we divide by (t -t 0 ) (with the sign and tend t to t 0 ), we have γ ′ (t 0 ) • (x -γ(t 0 )) = 0, this imply that x = z 0 -sν 0 where ν 0 is the outward normal of ∂Ω at z 0 ))

With this fact, we can say that S = {x, d(x, ∂Ω)

≤ ǫ} = {x = z 0 -sν z0 , z 0 ∈ ∂Ω, -ǫ ≤ s ≤ ǫ}. It is sufficient to work on ∂Ω. Let's consider a charts (z, D = B(z, 4ǫ z ), γ z ) with z ∈ ∂Ω such that ∪ z B(z, ǫ z ) is cover of ∂Ω . One can extract a finite cover (B(z k , ǫ k )), k = 1, .
.., m, by the area formula the measure of S ∩ B(z k , ǫ k ) is less than a kǫ (a ǫ-rectangle). For the reverse inequality, it is sufficient to consider one chart around one point of the boundary).

We write,

Ω |(∇(u i -u) • ∇η ǫ )|dx = Ω ǫ 3 |(∇(u i -u) • ∇η ǫ )|dx+ + Ω-Ω ǫ 3 (∇(u i -u) • ∇η ǫ )|dx. ( 5
)
Step 2.1:

Estimate of Ω-Ω ǫ 2 |(∇(u i -u) • ∇η ǫ )|dx.
First, we know from the elliptic estimates that

||∇η ǫ || L ∞ ≤ C 1 /ǫ 2 , C 1 de- pends on Ω
We know that (|∇u i |) i is bounded in L q , 1 < q < 2, we can extract from this sequence a subsequence which converge weakly to h ∈ L q . But, we know that we have locally the uniform convergence to |∇u| (by Brezis-Merle theorem), then, h = |∇u| a.e. Let q ′ be the conjugate of q.

We have, ∀f ∈ L q ′ (Ω)

Ω |∇u i |f dx → Ω |∇u|f dx
If we take f = 1 Ω-Ω ǫ 3 , we have:

for ǫ > 0 ∃ i 1 = i 1 (ǫ) ∈ N, i ≥ i 1 , Ω-Ω ǫ 3 |∇u i | ≤ Ω-Ω ǫ 2 |∇u| + ǫ 2 . Then, for i ≥ i 1 (ǫ), Ω-Ω ǫ 3 |∇u i | ≤ mes(Ω -Ω ǫ 3 )||∇u|| L ∞ + ǫ 2 = Cǫ.
Thus, we obtain,

Ω-Ω ǫ 3 |(∇(u i -u) • ∇η ǫ )|dx ≤ ǫC 1 (2k 1 ||∇u|| L ∞ + 1) (6)
The constant C 1 does not depend on ǫ but on Ω.

Step 2.2:

Estimate of Ω ǫ 2 |(∇(u i -u) • ∇η ǫ )|dx.
We know that, Ω ǫ ⊂⊂ Ω, and ( because of Brezis-Merle's interior estimates)

u i → u in C 1 (Ω ǫ 3 ). We have, ||∇(u i -u)|| L ∞ (Ω ǫ 3 ) ≤ ǫ 2 , for i ≥ i 3 = i 3 (ǫ).
We write,

Ωǫ 3 |(∇(u i -u) • ∇η ǫ )|dx ≤ ||∇(u i -u)|| L ∞ (Ω ǫ 3 ) ||∇η ǫ || L ∞ ≤ C 1 ǫ for i ≥ i 3 , For ǫ > 0, we have for i ∈ N, i ≥ max{i 1 , i 2 , i 3 }, Ω |(∇(u i -u) • ∇η ǫ )|dx ≤ ǫC 1 (2k 1 ||∇u|| L ∞ + 2) (7)
From ( 4) and (7), we have, for ǫ > 0, there is i

3 = i 3 (ǫ) ∈ N, i 3 = max{i 0 , i 1 , i 2 } such that, Ω |∆[(u i -u)η ǫ ]|dx ≤ 4π -ǫ 0 + ǫ2C 1 (2k 1 ||∇u|| L ∞ + 2 + C) (8)
We choose ǫ > 0 small enough to have a good estimate of (1).

Indeed, we have:

-∆[(u i -u)η ǫ ] = g i,ǫ in Ω, (u i -u)η ǫ = 0 on ∂Ω. with ||g i,ǫ || L 1 (Ω) ≤ 4π -ǫ 0 .
We can use Theorem 1 of [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF] to conclude that there is q ≥ q > 1 such that:

Vǫ(x0) e q|ui-u| dx ≤ Ω e q|ui-u|ηǫ dx ≤ C(ǫ, Ω).
where, V ǫ (x 0 ) is a neighberhood of x 0 in Ω.

Thus, for each x 0 ∈ ∂Ω -{x 1 , . . . , xm } there is ǫ x0 > 0, q x0 > 1 such that:

B(x0,ǫx 0 ) e qx 0 ui dx ≤ C, ∀ i. Now, we consider a cutoff function η ∈ C ∞ (R 2 ) such that: η ≡ 1 on B(x 0 , ǫ x0 /2) and η ≡ 0 on R 2 -B(x 0 , 2ǫ x0 /3).
We write,

-∆(u i η) = -log |x| 2d V i e ui η -2 < ∇u i |∇η > -u i ∆η.
By the elliptic estimates, (u i η) i is uniformly bounded in W 2,q1 (Ω) and also, in C 1 ( Ω).

Finaly, we have, for some ǫ > 0 small enough,

||u i || C 1,θ [B(x0,ǫ)] ≤ c 3 ∀ i. ( 9
)
We have proved that, there is a finite number of points x1 , . . . , xm such that the squence (u i ) i is locally uniformly bounded in Ω -{x 1 , . . . , xm }.

And, finaly, we have:

∂ ν u i → ∂ ν u + N j=1 α j δ xj , α j ≥ 4π weakly in the sens of measure L 1 (∂Ω). ( 10
)
Proof of theorem 1.2:

Without loss of generality, we can assume that 0 is a blow-up point. Since the boundary is an analytic curve γ(t), there is a neighborhood of 0 = x 1 such that the curve γ can be extend to a holomorphic map such that γ ′ (0) = 0 (series) and by the inverse mapping one can assume that this map is univalent around 0. In the case when the boundary is a simple Jordan curve the domain is simply connected. In the case that the domains has a finite number of holes it is conformally equivalent to a disk with a finite number of disks removed. Here we consider a general domain. Without loss of generality one can assume that γ(B + 1 ) ⊂ Ω and also γ(B - 1 ) ⊂ ( Ω) c and γ(-1, 1) ⊂ ∂Ω and γ is univalent. This means that (B 1 , γ) is a local chart around 0 for Ω and γ univalent. (This fact holds if we assume that we have an analytic domain, in the sense of Hofmann see [START_REF] Hofmann | Geometric and transformational properties of Lipschitz domains, Semmes-Kenig-Toro domains, and other classes of finite perimeter domains[END_REF], (below a graph of an analytic function), we have necessary the condition ∂ Ω = ∂Ω and the graph is analytic, in this case γ(t) = (t, φ(t)) with φ real analytic and an example of this fact is the unit disk around the point (0, 1) for example).

By this conformal transformation, we can assume that Ω = B + 1 , the half ball, and ∂ + B + 1 is the exterior part, a part which not contain 0 and on which u i converge in the C 1 norm to u. Let us consider B + ǫ , the half ball with radius ǫ > 0. Also, one can consider a C 1 domain (a rectangle between two half disks) and by charts its image is a C 1 domain).

We know that:

u i ∈ W 2,k (Ω) ∩ C 1,ǫ ( Ω).
Thus we can use integrations by parts (Gauss-Green-Riemann-Stokes formula). The second Pohozaev identity applied around the blow-up 0 = x 1 gives :

B + ǫ -∆u i (x • ∇u i )dx = - ∂ + B + ǫ g(∇u i )dσ, (11) with, g(∇u
i ) = (ν • ∇u i )(x • ∇u i ) -x • ν |∇u i | 2 2 .
After integration by parts, we obtain:

B + ǫ -log |x| 2d V i e ui (1+o (1)
)dx+

B + ǫ (x•∇V i )(-log |x| 2d )V i e ui dx+ ∂B + ǫ -log |x| 2d V i e ui (x•ν) = ∂ + B + ǫ g(∇u i )dσ + o(1), (12)
B + ǫ -log |x| 2d V e u (1+o (1)
)dx+

B + ǫ (x•∇V )(-log |x| 2d )V e u dx+ ∂B + ǫ -log |x| 2d V e u (x•ν) = = ∂ + B + ǫ g(∇u)dσ + o(1), (13)
Thus,

B + ǫ -log |x| 2d V i e ui dx - B + ǫ -log |x| 2d V e u dx+ + B + ǫ (x • ∇V i )(-log |x| 2d )V i e ui dx - B + ǫ (x • ∇V )(-log |x| 2d )V e u dx = = ∂ + B + ǫ g(∇u i ) -g(∇u)dσ + o(1) = o(1),
First, we tend i to infinity after ǫ to 0, we obtain:

lim ǫ→0 lim i→+∞ B + ǫ -log |x| 2d V i e ui dx = 0, (14)
But,

γ(B + ǫ ) -log |x| 2d V i e ui dx = ∂γ(B + ǫ ) ∂ ν u i + o(ǫ) + o(1) → α 1 > 0.
A contradiction.

Here we used a theorem of Hofmann see [START_REF] Hofmann | Geometric and transformational properties of Lipschitz domains, Semmes-Kenig-Toro domains, and other classes of finite perimeter domains[END_REF], which gives the fact that γ(B + ǫ ) is a Lipschitz domain. Also, we can see that γ((-ǫ, ǫ)) and γ(∂ + B + ǫ ) are submanifolds.

We start with a Lipschitz domain B + ǫ because it is convex and by the univalent and conformal map γ the image of this domain γ(B + ǫ ) is a Lipschitz domain and thus we can apply the integration by part and here we know the explicit formula of the unit outward normal it is the usual unit outward normal (normal to the tangent space of the boundary which we know explicitly because we have two submanifolds).

In the case of the disk D = Ω, it is sufficient to consider B(0, ǫ) ∩ D which is a Lipschitz domain because it is convex (and not necessarily γ(B + ǫ )). There is a version of the integration by part which is the Green-Riemann formula in dimension 2 on a domain Ω. This formula holds if we assume that there is a finite number of points y 1 , ..., y m such that ∂Ω -(y 1 , ..., y m ) is a C 1 manifold and for C 1 tests functions.

Remarks about the conformal map : 1-It sufficient to prove that γ 1 ((-ǫ, ǫ)) = ∂Ω ∩ γ1 (B ǫ ) = ∂Ω ∩ γ1 (B ǫ ) ∩ {|abscissa| < ǫ}, for ǫ > 0 small enough. Where γ1 is the holomorphic extension of γ 1 (t) = t + iφ(t). For this, we argue by contradiction, we have for z ǫ ∈ B ǫ , γ1 (z ǫ ) = (t ǫ , φ(t ǫ )) for |t ǫ | ≥ ǫ. Because γ1 is injective on B 1 and γ1 = γ 1 = t + iφ(t) on the real axis, we have necessirely

|t ǫ | ≥ 1. But, by continuity |γ 1 (z ǫ )| → 0 because z ǫ → 0. And, we use the fact that |γ 1 (z ǫ )| = |(t ǫ , φ(t ǫ ))| ≥ |t ǫ | ≥ 1,
to have a contradiction.) (This means that for a small radius when the graph go out from the ball, it never retruns to the ball). (This fact imply that, when we have a curve which cut ∂Ω in γ1 (B ǫ ) then the point have an abscissa such that |abscissa| < ǫ. This fact (by a contradiction with the fact ∂ Ω = ∂Ω and consider paths), imply that the image of the upper part of the ball is one side of the curve and the image of the lower part is in the other side of the curve.

2-Also, we can consider directly the coordinate T and change the function u(z) → u(T ) by z = T /a 1 + x 0 .

Set:

ψ(λ 1 , λ 2 ) → M ∈ Ω such that, --→ x 0 M = λ 1 i ′ 1 + λ 2 j ′ 1 where (i ′ 1 , j ′ 1 ) is a basis such that i ′ 1 = e -iθ i 1 , j ′ 1 = e -iθ j 1 .
And, φ(x 1 , x 2 ) → M such that, --→ OM = x 1 i 1 + x 2 i 2 the canonical basis (i 1 , j 1 ). Then, we have two charts φ and ψ and the complex affix T M = λ 1 + iλ 2 and z M = x 1 + ix 2 are such that (transition map): We write:

T M /e iθ + x 0 = z M = φ -1 oψ(λ 1 , λ 2 ),
∆u(M ) = ∆ λ (uoψ(λ 1 , λ 2 ))
And then we apply the conformal map γ1 which send the affix T M , M in a neighborhood of x 0 ∈ ∂Ω to B ǫ with the fact that send T M , M ∈ ∂Ω to the real axis (-ǫ, ǫ) and the other parts of Ω and Ωc .

3-We can remark that a definition of C k , k ≥ 1 domain, is equivalent to a definition of a submanifold with the condition ∂ Ω = ∂Ω or Ω = Ω.

Remark 2: about a variational formulation. we consider a solutions in the sense of distribution. By the same argument (in the proof of the maximum principle obtained by Kato's inequality W 1,1 0 is sufficient), see [START_REF] Brezis | Nonlinear elliptic equations with measures revisited[END_REF], we prove that the solutions are in the sense C 2 0 of Agmon, see [START_REF] Aubin | Some Nonlinear Problems in Riemannian Geometry[END_REF]. Also, we have corollary 1 of [START_REF] Brezis | Uniform estimates and Blow-up behavior for solutions of -∆u = V (x)e u in two dimension[END_REF]. We use Agmon's regularity theorem. We return to the usual variational formulation in W 1,2 0 and thus we have the estimate of [START_REF] Brezis | Semi-linear second-order elliptic equations in L1[END_REF] or by Stampacchia duality theorem in W 1,2 0 .

  We have:∂ λ1 = cos θ∂ x1 + sin θ∂ x2 , ∂ λ2 = -sin θ∂ x1 + cos θ∂ x2 ,Thus, the metric in the chart ψ or coordinates (λ 1 , λ 2 ) is : g λ ij = δ ij and the Laplacian in the two charts, ψ and φ are the usual Laplacian ∂ λ1λ1 + ∂ λ2λ2 .