G Saez-Mischlich
email: g.saez-mischlich@isae-supaero.fr

J Sierra-Ausin

G Grondin

J Gressier

On the properties of high-order least-squares finite-volume schemes

Finite-volume

High-order Least-squares k-exact High-order finite-volume schemes based on polynomial least-squares methods are an active research topic for the discretization of hyperbolic equations as they allow to obtain high-order spatial discretization schemes in arbitrary meshes. However, few studies have analyzed their performance in good-quality/near-to-uniform meshes, which are commonly used as a meshing strategy in zones where turbulent effects are important. In this paper, the theoretical numerical properties of commonly used least-squares (LSQ) k-exact high-order finite volume schemes are studied in one-dimensional and in several twodimensional meshes (with some remarks regarding their properties in three-dimensional meshes). These results are compared to those obtained using fully-constrained polynomial reconstructions only compatible with structured meshes. The numerical properties of the schemes are investigated through the von Neumann analysis methodology applied to the one-dimensional and two-dimensional finite volume formulation, including temporal discretization errors. This analysis is also extended to non-uniform and unstructured twodimensional meshes. At last, the schemes are tested with several numerical experiments using the linear advection, the Euler equations and the Navier-Stokes equations. The analysis of both studies yields similar conclusions regarding the numerical errors and stability of the different studied schemes showing that the high-order least-squares finite volume schemes yield stable and robust results across different uniform and non-uniform unstructured meshes. However, their performance is heavily degraded in simulations with low mesh resolution compared to schemes specially catered to structured meshes. On the other hand, the latter schemes lack stability and robustness in general structured meshes and its formulation cannot be straightforwardly extended to unstructured meshes. Moreover, this work shows that the use of weighted-LSQ can drastically improve the results of LSQ schemes in under-resolved simulations.

Introduction

Finite volume (FV) is a well-established numerical technique for the reconstruction of solutions of hyperbolic partial differential equations (PDEs), e.g. Euler's and Maxwell's equations. Its success relies on two key properties: these schemes fulfill the discrete global conservation property and L ∞ -stable schemes are relatively easy to design, e.g. non-oscillatory methods (ENO [START_REF] Liu | Wls-eno: weighted-least-squares based essentially non-oscillatory schemes for finite volume methods on unstructured meshes[END_REF], WENO [START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF][START_REF] Tsoutsanis | Stencil selection algorithms for weno schemes on unstructured meshes[END_REF], MOOD [START_REF] Farmakis | Weno schemes on unstructured meshes using a relaxed a posteriori mood limiting approach[END_REF]). Finite volume is a particular case of a Bubnov-Galerkin method, i.e. test and solution functions do not belong to the same functional space, where test functions are constant along each element and solution functions are polynomial in each element. Therefore, higher-order reconstructions are possible if one considers polynomial reconstructions of higher degree. In particular, the quality of a finite volume reconstruction depends on two operations: variable reconstruction from neighboring element-averaged and face average normal flux reconstruction. The former is usually realized by interpolation and it can lead to the well-known k-exact schemes [START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF][START_REF] Martin | Implementation and validation of semi-implicit weno schemes using openfoam®[END_REF][START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Dumbser | Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems[END_REF], whereas the latter is commonly evaluated in high-order formulations via quadrature integration [START_REF] Zhang | On the order of accuracy and numerical performance of two classes of finite volume WENO schemes[END_REF][START_REF] Titarev | Finite-volume weno schemes for three-dimensional conservation laws[END_REF][START_REF] Haider | A high-order interpolation for the finite volume method: the coupled least squares reconstruction[END_REF][START_REF] Tsoutsanis | Comparison of structured-and unstructured-grid, compressible and incompressible methods using the vortex pairing problem[END_REF]. There exist other face average normal flux reconstruction methods based on successive error corrections via Taylor expansion series. Interested reader is referred to the studies of [START_REF] Mccorquodale | A high-order finite-volume method for conservation laws on locally refined grids[END_REF][START_REF] Du | A high-order modified finite volume WENO method on 3d Cartesian grids[END_REF][START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible rans-les simulations on fully unstructured grids[END_REF] for more information on these methods.

Connectivity characteristics of the discretization of the physical domain play an important role in the numerical properties of finite volume schemes. Structured numerical grids possess a constant number of face and vertex neighbors to each element and there exists a labeling function relating N-tuples of indexes to element centroids, e.g. in two dimensional domains, there is a labeling (i, j) ∈ Z 2 for each element centroid. Under certain conditions, these properties allow to determine coefficients of the reconstruction scheme by solving a non-singular linear system of equations [START_REF] Arbogast | Finite volume weno schemes for nonlinear parabolic problems with degenerate diffusion on non-uniform meshes[END_REF][START_REF] Fosso Pouangué | Curvilinear finite-volume schemes using high-order compact interpolation[END_REF]. Nonetheless, unstructured grids do not have, in general, a constant number of neighbors for each element and the labeling of each element centroid cannot be decomposed onto a Cartesian product of N-sets of indexes for N-dimensional domains. For such reason, the stencils employed for the reconstruction of the solution are not unique and its selection is crucial in the numerical properties of the scheme [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Jalali | Higher-order finite volume solution reconstruction on highly anisotropic meshes[END_REF][START_REF] Vasconcelos | A very high-order finite volume method based on weighted least squares for elliptic operators on polyhedral unstructured grids[END_REF]. The determination of the variable reconstruction at the faces may be determined by an iterative procedure of successive Green-Gauss [START_REF] Mccorquodale | A high-order finite-volume method for conservation laws on locally refined grids[END_REF][START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible rans-les simulations on fully unstructured grids[END_REF] approximations or by the selection of a stencil and then solving a L 2 -minimization problem using a least-squares procedure [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Titarev | Finite-volume weno schemes for three-dimensional conservation laws[END_REF][START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF].

The present study aims to analyze numerical properties, e.g. order of accuracy, stability, dispersion and diffusion errors, of finite volume schemes on structured and unstructured grids with a quadrilateral discretization applied to hyperbolic problems with smooth solutions. In particular the emphasis is on properties of least-squares schemes with respect to their structured counter-part. Firstly, the general framework of finite volume method is synthesized in Section 2. Then in Section 3 the attention is paid on a non-compact reconstruction method on unstructured grids by a least-squares procedure with a given stencil based on the superposition of topological layers of neighbors. Such a methodology is particularized in Section 3.2 for structured grids. An analytical study of the truncation error, i.e. order of accuracy, is carried out in Section 4.1 for one and two-dimensional structured grids. Additionally, a numerical test based on the k-polynomial reconstruction property is employed to determine the order of accuracy of least-squares schemes on unstructured grids. Spectral properties of schemes such as -L 2 -stability, dispersive and diffusion errors -are studied in Section 5. Finally, a set of numerical examples is reported in Section 6 to validate the analytical findings.

Finite volume applied to conservation laws

Finite volume methods naturally apply to conservation laws expressed in divergence form:

∂u(x, t) ∂t + ∇ • F (u, ∇u) = 0 ∀ x in ⊂ R N d , N d = 1, 2, 3 ( 1 )
where u(x, t) is a conserved scalar or vector quantity, x is the position vector, F is the flux operator and N d is the number of spatial dimensions of the problem. The flux operator is usually decomposed in two components: the convective flux F c (u) which depends only on the conserved variable and not its gradient, and the viscous flux F v (u, ∇u) which depends both on the conserved variable and its gradient. The numerical resolution of Eq. ( 1) is carried out via the finite volume method.

Partitioning of the physical domain

Several types of unstructured meshes are common in finite volume methods. Most polygonal meshes are composed of N d -simplexes (triangles or tetrahedrons) and tensor-product elements (quadrilaterals or hexahedrons). In this work the meshes will be decomposed in N V linear tensor-product elements or cells. As it is common in the literature, independently of the shape of the chosen polygon covering the physical domain , the partitioning of the numerical domain h produces the partition V h , which possesses the usual required properties for finite element method, see [START_REF] Feistauer | Mathematical and Computational Methods for Compressible Flow[END_REF]. Cells V i ∈ V h are closed convex volume elements, satisfying the following conformal property:

if V i , V j ∈ V h , then either V i ∩ V j = ∅, or V i ∩ V j is a common ⎧ ⎨ ⎩ vertex, edge or face of V i and V j , or V i = V j .
(2)

In addition, let us consider the characteristic length h i of each volume element, which proves to be useful in the analysis of numerical properties of the space discretization. This parameter is defined [START_REF] Ndjinga | Spectral stability of finite volume schemes for linear hyperbolic systems[END_REF] as

h i = |L i | ∈L i |S | |V i | , (3)
where L i is the set of indices of the faces S of element V i , with respective cardinality |L i | and measures |S | and |V i |.

The global edge size h is defined as the minimum edge size h = min i h i .

(4) Finally, the centroid x i of each element is defined as

x i = 1 |V i | V i
x dV .

(5)

Finite volume semi-discretization

The finite volume formulation of Eq. ( 1) describes the evolution of the cell average u i of the conserved quantity on cell V i defined as

u i = 1 |V i | V i u dV i , (6)
and is obtained by integration of Eq. (1) on volume element V i , assuming a static mesh with no deformation or movement of volume elements, which yields

du i dt + 1 |V i | ∈L i S F (u) • n i, dS = 0 (7)
where n i, is the unitary outward-facing face normal of face S with respect to cell V i computed as in [START_REF] Naff | Shape functions for velocity interpolation in general hexahedral cells[END_REF]. In this work, only planar faces are considered, which implies that face normals are constant on each face. Let us define the face average operator on a given face S as

= 1 |S | S (x) dS . ( 8
)
Face average normal fluxes of Eq. ( 7) are also computed as

F i, = 1 |S | S F (u) • n i, dS . (9)
For the sake of brevity, face average normal fluxes will also be referred as face average fluxes in this work. Eq. ( 7) is then transformed into the following discrete system by spatial discretization of the fluxes into the semi-discrete system

du δ i dt = - 1 |V i | ∈L i F δ i, |S | , ( 10
)
where the super-index δ symbolizes a numerical approximation. It can then be clearly observed that the numerical errors in the FV formulation are introduced through the discretization of the face average fluxes and the time integrator. A review of the most common approaches for the reconstruction of high-order face average fluxes (higher than second order) is carried out in Section 2.3. In this study, the time scheme used for time integration of Eq. ( 10) is a third order Runge-Kutta (RK) propagator method (RK33) [START_REF] Butcher | Numerical Methods for Ordinary Differential Equations[END_REF].

Finite volume face average flux reconstruction methods

In this work, high-order discretization of the face average fluxes (higher than second order) will be used to discretize the FV equations. In high-order FV methods, the most common techniques for the numerical computation of the flux face integrals in Eq. ( 9) consist of quadrature methods using local values of the conserved variable at the given quadrature points (see [START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF][START_REF] Tsoutsanis | Comparison of structured-and unstructured-grid, compressible and incompressible methods using the vortex pairing problem[END_REF][START_REF] Tsoutsanis | Extended bounds limiter for high-order finite-volume schemes on unstructured meshes[END_REF][START_REF] Tsoutsanis | Stencil selection algorithms for weno schemes on unstructured meshes[END_REF][START_REF] Zhang | On the order of accuracy and numerical performance of two classes of finite volume WENO schemes[END_REF][START_REF] Farmakis | Weno schemes on unstructured meshes using a relaxed a posteriori mood limiting approach[END_REF][START_REF] Titarev | Finite-volume weno schemes for three-dimensional conservation laws[END_REF]) where G is the set of indices of the weights w g associated with the quadrature points x g of face S (see Fig. 1). In these methods, a polynomial representation of the conserved variables in each cell is mandatory to obtain a high-order approximation of the values of the conserved variables at quadrature points of each cell. To abbreviate the formulation, the values of the conserved variables and their gradients at a given quadrature point with index g are referred to u g and ∇u g respectively. Eq. ( 11) can be rewritten in terms of the numerical values of the conserved variables at the quadrature points as

F i, ≈ F δ i, = g∈G w g F u x g , ∇u x g • n i,g , (11)
F δ i, = g∈G w g F δ u δ,- g , ∇u δ,- g , u δ,+ g , ∇u δ,+ g , n g = g∈G w g F δ c,i,g + F δ v,i,g , ( 12
)
where the numerical convective and viscous flux operators F δ c,i,g and F δ v,i,g have been introduced. The symbols ± represent the two possible numerical reconstructions (commonly referred to as left and right states) of the conserved variables at a given quadrature point of a non-boundary face. The left and right reconstructions, u δ,g and u δ,+ g , are respectively computed using the stencils of the owner cell of the face V ,o and V ,n of face S . Meanwhile, the right reconstruction u δ,+ g is computed using the stencil of the neighbor cell of the face V ,n . The face normal at each quadrature point n g of face S with respect to their owner and neighbor cells is defined such that n g • (x nx o ) ≥ 0. The reader may refer to Section 3 for a detailed explanation of the reconstruction of the conserved variables at the quadrature points.

Remark 1. We are aware that there exist other high-order finite-volume flux reconstruction methods for structured and unstructured meshes which were initially introduced for structured grids by [START_REF] Mccorquodale | A high-order finite-volume method for conservation laws on locally refined grids[END_REF]. These methods rely on the computation of the Taylor expansion of the numerical face average flux. A 6 th -order FV WENO formulation was recently proposed in [START_REF] Du | A high-order modified finite volume WENO method on 3d Cartesian grids[END_REF].

An extension of the method for unstructured grids was developed in [START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible rans-les simulations on fully unstructured grids[END_REF].

Processing of face average convective fluxes

To compute the numerical convective flux at each internal quadrature point F δ c,i,g , the following weighted expression is used

F δ c,i,g = f g F c,upwind u δ,- g , u δ,+ g , n g + (1 -f g )F c,central u δ,- g , u δ,+ g , n g sgn(n i,g • n g )) , ( 13
)
where f g ∈ [0, 1] is a weighting factor, F c,upwind is an upwind flux whose formulation is usually especially tailored for a given non-linear problem using Riemann solvers. For the Euler equations (Section 2.5.2) one can find Riemann solvers such as Roe [START_REF] Roe | Approximate Riemann solvers, parameter vectors, and difference schemes[END_REF], HLL or HLLC [START_REF] Toro | The HLL and HLLC Riemann solvers[END_REF]. F c,central is a central flux which is computed in this work as

F c,central u δ,- g , u δ,+ g , n g = 1 2 F c u δ,- g + F c u δ,+ g • n g . ( 14
)
For f = f g = 1 upwind flux reconstruction schemes are obtained while the case f = f g = 0 is representative of central schemes. In this work, both upwind and central schemes will be studied and they will be denoted with the initial letter U or C respectively. The order of accuracy of these face average flux reconstruction methods is determined by two parameters: the order of accuracy of the interpolation method of the conserved variable at the quadrature points (analyzed in our studies in Section 4.2) and that of quadrature method used to evaluate the face average flux (discussed in Section 3.3).

Processing of face average viscous fluxes

In contrast with the convective terms, the numerical viscous fluxes involve reconstruction of gradients at the quadrature points on each side of the face. In this work, the numerical viscous flux F v,i,g at each quadrature point is directly computed as

F δ v,i,g = F v u δ g , ∇u δ g • n i,g , (15)
where u δ g and ∇u δ g are computed with a central scheme as

u δ g = 1 2 u δ,- g + u δ,+ g and ∇u δ g = 1 2 ∇u δ,- g + ∇u δ,+ g . ( 16
)
One could also consider other types of numerical viscous fluxes such as the Interior-Penalty (IP) method or the Rebay-Bassi (RB) approach. The interested reader is referred to [START_REF] Gassner | A contribution to the construction of diffusion fluxes for finite volume and discontinuous Galerkin schemes[END_REF][START_REF] Jalali | Accuracy analysis of unstructured finite volume discretization schemes for diffusive fluxes[END_REF] for more information on numerical viscous fluxes formulations. At last, Section 3.4 discusses the computation of the gradient of the conserved variables at quadrature points, using the numerical schemes of choice of this work.

Finite volume variables reconstruction methods

This section is devoted to the description of the state of the art of the FV numerical schemes which provide high-order reconstructions of the conserved variables at the face quadrature points.

k-exact schemes for structured meshes

k-exact methods for structured meshes have been developed in [START_REF] Zhang | On the order of accuracy and numerical performance of two classes of finite volume WENO schemes[END_REF][START_REF] Titarev | Finite-volume weno schemes for three-dimensional conservation laws[END_REF]. They allow to compute high-order face average flux reconstructions using Eq. [START_REF] Mccorquodale | A high-order finite-volume method for conservation laws on locally refined grids[END_REF]. They are based on the use of tensor-product of one-dimensional polynomials to represent a modal basis of the solution. The extension of this method for non-uniform structured grids was proposed in [START_REF] Arbogast | Finite volume weno schemes for nonlinear parabolic problems with degenerate diffusion on non-uniform meshes[END_REF] and will be analyzed in this work in Section 3.2.

k-exact schemes for arbitrary meshes

There exists a variety of k-exact methods which are applicable to any kind of mesh and mesh elements. These methods are based on the computation of a p-degree polynomial to approximate the conserved variables in each cell. The polynomial approximation enables a high-order interpolation of the conserved variable at quadrature points. These interpolations are used to compute numerical fluxes using quadrature integrals Eq. [START_REF] Mccorquodale | A high-order finite-volume method for conservation laws on locally refined grids[END_REF]. Several types of k-exact reconstruction methods can be listed:

1. The compact methods of [START_REF] Wang | Compact high order finite volume method on unstructured grids i: basic formulations and one-dimensional schemes[END_REF][START_REF] Haider | A high-order interpolation for the finite volume method: the coupled least squares reconstruction[END_REF] in which coefficients of the polynomial are computed iteratively using data from the closest neighbors via successive least-squares problems. 2. The non-compact method of [START_REF] Pont | Multiple-correction hybrid k-exact schemes for high-order compressible rans-les simulations on fully unstructured grids[END_REF], where the coefficients of a polynomial approximation of the conserved variables are approximated through several LSQ problems on each layer of neighbors. 3. The non-conservative moving least-squares (MLS) method proposed by [START_REF] Cueto-Felgueroso | Finite volume solvers and moving least-squares approximations for the compressible Navier-Stokes equations on unstructured grids[END_REF] which can also be used in mesh-free formulations. 4. The non-compact least-squares (LSQ) method proposed by [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF][START_REF] Tsoutsanis | Comparison of structured-and unstructured-grid, compressible and incompressible methods using the vortex pairing problem[END_REF][START_REF] Tsoutsanis | Stencil selection algorithms for weno schemes on unstructured meshes[END_REF] which is one pf the main topic of research of this work. In these methods the conserved variables are spatially estimated by a polynomial of a given degree whose coefficients are obtained from neighbor cells through the solution of a single least-squares problem in a non-compact manner. These methods are commonly utilized in conjunction with ENO or WENO [START_REF] Farmakis | Weno schemes on unstructured meshes using a relaxed a posteriori mood limiting approach[END_REF] methods by reconstructing a plethora of polynomial approximations of the conserved variables at each cell using different stencil topologies. These polynomials can then be combined by measuring their smoothness ensuring the correct treatment of solution discontinuities.

Choice of conservative equations

Three systems of conservation laws are considered in this study: the linear advection equation, the Euler equations and the Navier-Stokes equations.

Linear advection equation

The linear advection equation is a single equation conservation law in which the flux operator F of the conservation law is linear F = cu, where c is the linear advection velocity and it is supposed to be constant in this work. In this case, Eq. (1) simplifies to ∂u ∂t

+ c • ∇u = 0 . ( 17
)
Supposing planar faces, the finite volume formulation of Eq. ( 17) can be written as

du i dt + 1 |V i | ∈L i u c • n i,l |S | = 0 . ( 18
)
In future references to the linear advection equation in two-dimensional context, the linear advection velocity will be expressed as

c = c(cos ζ, sin ζ ) ( 19
)
where ζ is the advection velocity angle.

Euler equations

The Euler equations are written in conservative form as

∂ U ∂t + ∇ • F (U ) = 0 , ( 20
)
where U is the vector of conservative variables given by

U = ⎛ ⎝ ρ ρu ρ E ⎞ ⎠ , ( 21
)
ρ is the fluid density, ρu the momentum, u the velocity and ρ E is the total energy. The flux operator of the Euler equations can be written as

F (U ) = ⎛ ⎝ ρu ρu ⊗ u + P I N d (ρ E + P )u ⎞ ⎠ , ( 22
)
where the symbol ⊗ represents the dyadic operator, P is the pressure and I N d ∈ R N d ×N d is the identity matrix. The nonlinear system is closed with the equation of state

ρ E = P γ -1 + 1 2 ρu • u , ( 23
)
where γ is the ratio of specific heats. The primitive variables of the Euler equations ρ, u, P and the temperature can be computed from the conservative variables using the equation of state Eq. ( 23) and the following additional expressions:

u = ρu ρ , P = (γ -1) ρ E - 1 2 ρu • u , = P ρr , c p = γ r (γ -1) , ( 24
)
where r is the gas constant and c p is the heat capacity at constant pressure. The flux operator of the Euler equations is non-linear and therefore the only possibility of obtaining a formally high-order finite volume numerical solution is by using the high-order face average flux reconstructions such as Eq. ( 12).

Compressible Navier-Stokes equations

The compressible Navier-Stokes equations are a system of non-linear equations similar to the Euler equations. Both sets of equations share the same conserved (Eq. ( 21)) and primitive variables (Eq. ( 24)) and have similar flux operators. However, the Navier-Stokes equations present additional viscous terms which depend on the values and the gradients of the conserved and/or the primitive variables. These viscous terms are denoted as F v (U , ∇U ). The flux operator of the Navier-Stokes equations can be written as

F (U ) = ⎛ ⎝ ρu ρu ⊗ u + P I N d (ρ E + P )u ⎞ ⎠ + ⎛ ⎝ 0 -σ - μc p Pr ∇ -σ • u ⎞ ⎠ , ( 25
)
where μ is the kinematic viscosity, Pr is the Prandtl number and σ is the viscous stress tensor defined as

σ = μ ∇u + ∇u T - 2 3 ∇ • u I N d . ( 26
)
The velocity and temperature gradients can be computed from the gradients of the conservative variables through the chain rule and the relationship between the conservative and primitive variables expressed in Eq. (24).

k-exact finite volume schemes

This section addresses the different non-compact k-exact schemes discussed in this work. Section 3.1 discusses leastsquares k-exact schemes compatible with both structured and unstructured meshes. Section 3.2 discusses fully constrained k-exact schemes only compatible with structured meshes. The formulations are written supposing two-dimensional problems. The extension to one-dimensional and three-dimensional problems is trivial. However, remarks regarding the extension to different number of dimensions will be added when needed.

k-exact least-squares schemes

The present section addresses the least-squares k-exact reconstruction approach in unstructured meshes used in this work. The methodology has been developed by [START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF] and [START_REF] Martin | Implementation and validation of semi-implicit weno schemes using openfoam®[END_REF] whose work is a continuation of [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF] and [START_REF] Dumbser | Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems[END_REF]. Firstly, let us consider a reference coordinate space for each cell V i denoted as ξ = (ξ, η) related to the global coordinates x = (x, y) via a linear mapping:

x y = x 0 y 0 + J • ξ η , ( 27
)
where (x 0 , y 0 ) are the coordinates of an arbitrary reference node of the cell and J is the Jacobian matrix computed as

J = x 1 -x 0 x 2 -x 0 y 1 -y 0 y 2 -y 0 . ( 28
)
In the latter expression, sub-indices 1 and 2 refer to two nodes of cell i which are connected to the chosen reference node. The choice of reference node is arbitrary and therefore, the properties of these schemes are dependent on the vertexcell connectivity of the mesh. It is important to note that any integral evaluation in the global coordinate space can be transformed into the reference space as

1 |V i | V i T (x) dV = 1 |V i | V i T (x) dx d y = 1 |V i | V i T (ξ ) dξ dη , ( 29
)
where

|V i | = |V i ||J -1 |
is the volume measure of cell i in the reference space and |J -1 | is the determinant of the inverse Jacobian matrix. In these methods a polynomial P i (ξ ) is built for each cell i of the domain as

P i (ξ ) = K γ =1 a γ φ γ (ξ ) = u δ i (ξ ) , (30)
where the superscript δ denotes a numerical approximation of a given function over a cell and the functions {φ γ } are a family of polynomial functions of the form

∀ γ ∈ [1, K ] φ γ (ξ ) = γ (ξ ) -C γ , ( 31
) such that { γ (ξ ), γ ∈ [1, K ]} = {(ξ -ξ i ) q (η -η i ) r , q ≥ 0, r ≥ 0, q + r ≤ p} . ( 32
)
In Eqs. [START_REF] Jaśkowiec | High-order cubature rules for tetrahedra[END_REF] and [START_REF] Schlömer | nschloe/quadpy v0.16[END_REF], (ξ i , η i ) is the centroid of cell i in the reference space, p is the degree of the polynomial, K is the cardinality of the bases {φ γ (ξ )} and { γ (ξ )} and depends on the spatial dimension N d as

K = p + N d p , ( 33
)
and {C γ } is a family of constants whose values will be determined further into this section.

In general, the coefficients a γ of Eq. ( 30) are determined by setting integral equations for the root cell i and a number of its closest neighbors, which together with the root cell form the stencil named C i with cardinality M = |C i |. A discussion on the definition of the stencil for these LSQ methods is developed in Section 3.1.4. In the least-squares schemes studied in this section the number of neighbors is chosen such that M > K . The setting M = K is specific to a small family of k-exact schemes for structured meshes as discussed in Section 3.2.

Definition of the LSQ problem

The LSQ problem is generated by equating the average of the reconstructed polynomial P i on each cell V j of the stencil C i with the cell average u j of the solution as in [START_REF] Lawson | Mathematics, Solving Least Squares Problems[END_REF] which can be rearranged as

∀ j ∈ C i \ {i} 1 |V j | V j P i (ξ ) dξ dη = u j ,
∀ j ∈ C i \ {i} K γ =1 a γ 1 |V j | V j φ γ (ξ ) dξ dη = u j . ( 35
)
The least-squares problem is constrained on cell i by the equation 1

|V i | V i P i (ξ ) dξ dη = u i , ( 36
)
ensuring that the cell average is exactly reconstructed by the polynomial P i .

Let us impose that the constant polynomial obtained for q = r = 0 in Eq. ( 32) be assigned to index γ = 1 with a null constant C γ =1 = 0 and that the integral of φ γ on cell i be zero for every other γ :

∀ γ > 1 1 |V i | V i φ γ (ξ ) dξ dη = 0 . ( 37
)
By inserting the previous assumptions in Eq. [START_REF] Petrovskaya | The accuracy of least-squares approximation on highly stretched meshes[END_REF], the following values for the constants C γ are obtained (as in the studies of [START_REF] Dumbser | Arbitrary high order non-oscillatory finite volume schemes on unstructured meshes for linear hyperbolic systems[END_REF][START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF]): [START_REF] Epperson | An Introduction to Numerical Methods and Analysis[END_REF] from which ensues the values of a 1 , as described in [START_REF] Tsoutsanis | Weno schemes on arbitrary mixed-element unstructured meshes in three space dimensions[END_REF][START_REF] Tsoutsanis | Comparison of structured-and unstructured-grid, compressible and incompressible methods using the vortex pairing problem[END_REF]:

∀ γ > 1 C γ = 1 |V i | V i γ (ξ ) dξ dη ,
a 1 = a 1 φ 1 = a 1 1 = u i . ( 39
)
We stress out that this choice is arbitrary and will indeed modify the numerical properties and stability of the schemes, although the analysis of this topic will not be addressed in this work.

From the above expressions and the constraint imposed in Eq. ( 36), the over-constrained problem may be recast in matrix form as AX = b , [START_REF] Manzanero | Insights on aliasing driven instabilities for advection equations with application to Gauss-Lobatto discontinuous Galerkin methods[END_REF] where A is a matrix in R (M-1)×(K -1) with elements

A jγ = 1 |V j | V j φ γ +1 (ξ ) dξ dη = 1 |V j | V j γ +1 (ξ ) dξ dη -C γ +1 . ( 41
)
Vectors X and b belong to R K -1 and R M-1 respectively and are set as

∀ γ ∈ [1, K -1] X γ = a γ +1 , ( 42
) and ∀ j ∈ C i \ {i} b j = u j -u i . ( 43
)
Remark 2. In this work, all terms in matrix A are computed using an exact numerical integration quadrature of the polynomial basis up to machine precision. In one-dimensional meshes, this is performed using Gauss-Legendre quadratures of appropriate order. In two-dimensional and three-dimensional meshes, the volume integrals are evaluated by decomposing the tensor-product elements (quadrilaterals or hexahedrons) in N d -simplex elements (triangles or tetrahedrons) and evaluating the volume integrals on those elements. To do so, this work relies on the quadratures schemes described in [START_REF] Xiao | A numerical algorithm for the construction of efficient quadrature rules in two and higher dimensions[END_REF] for triangles (2-simplexes) and the quadratures scheme of [START_REF] Jaśkowiec | High-order cubature rules for tetrahedra[END_REF] for tetrahedrons (3-simplexes). The quadrature points positions and weights of these schemes were obtained using the QUADPY software [START_REF] Schlömer | nschloe/quadpy v0.16[END_REF]. The decomposition of tensor-product elements in simplex cells does not incur geometric discretization errors if one only considers tensor-product elements with planar faces. Therefore, the volume integrals are evaluated by subdividing the cells into a set of N d -simplexes and applying quadrature rules on those subelements. As the determinant of the Jacobian of the isoparametric transform of simplex elements is constant, the integrals can be exactly evaluated up to machine round-off precision with an appropriate quadrature. It is important to note that the exact evaluation of volume integrals in arbitrary tensor-product elements with Gauss-Legendre quadratures is also possible but one needs to take into account that the determinant of the isoparametric Jacobian transformation (needed to evaluate Gauss-Legendre quadratures) is a 3 rd degree polynomial [START_REF] Abe | On the freestream preservation of high-order conservative flux-reconstruction schemes[END_REF] in the worst case scenario (hexahedron cells). Therefore, the Gauss-Legendre quadrature should be capable of exactly integrating a p + 3 degree polynomial.

Resolution of the least-squares problem

The over-constrained problem described in Eq. ( 40) is solved with the least-squares L 2 -minimization approach, which implies the computation of the pseudoinverse A + of the matrix A. In this work, the pseudoinverse is obtained using the SVD (Singular Value Decomposition) method [START_REF] Lawson | Mathematics, Solving Least Squares Problems[END_REF][START_REF] Martin | Implementation and validation of semi-implicit weno schemes using openfoam®[END_REF], implemented in the EIGEN library [START_REF] Guennebaud | Eigen v[END_REF]. This process is only performed during the simulation's initialization process, because A + uniquely depends on the geometrical properties of the mesh. The solution vector can be then obtained as

X = A + b . ( 44
)
The choice of stencil and the polynomial order will determine the properties of the least-squares spatial interpolation schemes. This issue will be discussed in Section 3.1.4.

Distance weighting in the LSQ problem

Several studies [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Petrovskaya | The accuracy of least-squares approximation on highly stretched meshes[END_REF][START_REF] Jalali | Higher-order finite volume solution reconstruction on highly anisotropic meshes[END_REF] have pointed out the influence of the addition of distance weighting methods to the leastsquares problem on the numerical properties of the k-exact schemes. The distance-weighted least-squares problem is obtained by rewriting the previous least-squares problem as A w X = b w , [START_REF] Alhawwary | Fourier analysis and evaluation of dg, fd and compact difference methods for conservation laws[END_REF] with

A w jγ = 1 ||x j -x i || α A jγ , ( 46
)
and

b w j = 1 ||x j -x i || α b j , ( 47
)
where, α is the weighting exponent, and x i and x j are the vector coordinates of the centers of cells V i and V j . The solution vector can then be computed as

X = A w,+ b w , ( 48
)
where A w,+ is the pseudoinverse matrix of A w . Schemes obtained with distance weighting least-squares method will be designated as LSQpW. In this work the weighting exponent is fixed as α = 5

2 , chosen as a compromise between the scheme accuracy and the conditioning of matrix A. Constraint equation [START_REF] Petrovskaya | The accuracy of least-squares approximation on highly stretched meshes[END_REF] and constants C γ from Eq. ( 38) remain unchanged with respect to the non-weighted least-squares problem.

Definition of the k-exact LSQ schemes

The accuracy of polynomial P i used for the reconstruction of the conservative variable greatly depends on the choice of stencil C in general. The choice of stencil was proven crucial for the accuracy of k-exact schemes in unstructured meshes [START_REF] Tsoutsanis | Stencil selection algorithms for weno schemes on unstructured meshes[END_REF]. Nonetheless, few studies have analyzed the effects of the stencil shape and topology on the accuracy of these schemes for smooth meshes using quadrilateral or hexahedral elements. Most studies use stencils by gathering a given number of topological neighbors with connected faces or vertices with each other. The number of neighbor elements in the stencil is commonly computed as a function of the polynomial cardinal functions using the β coefficient as

β = M -1 K -1 . ( 49
)
Please note that M -1 is the number of neighbors of a given root cell used to reconstruct a least-squares polynomial and K -1 is the size of the polynomial basis for the least-squares procedure given that the constant polynomial was specifically constrained as stated in Eq. [START_REF] Burden | Numerical Analysis[END_REF].

In practice, the value β = 2 is usually adopted. However, this condition generates ill-formed stencils in smooth meshes with quadrilateral and hexahedral meshes, because it is not possible to ensure topologically symmetric stencils for every cell if the value of β is fixed to a value which is not representative of the number of topological neighbors of a given cell.

If the connectivity of the stencil of neighboring cells is dissimilar, the truncation error will not behave regularly, therefore some error canceling properties of the truncation error will be lost even in smoothly varying grids.

In this study, it has been decided to associate a number of neighbors M and a parameter β to a given topological shape of the neighbor cells of the root cell which share a given face with it (face neighbor cells) similarly to the methods described in [START_REF] Ollivier-Gooch | A high-order-accurate unstructured mesh finite-volume scheme for the advection-diffusion equation[END_REF][START_REF] Jalali | Higher-order finite volume solution reconstruction on highly anisotropic meshes[END_REF][START_REF] Vasconcelos | A very high-order finite volume method based on weighted least squares for elliptic operators on polyhedral unstructured grids[END_REF]. This ensures geometric symmetry of the stencils in structured meshes and a smooth variation of the truncation error, provided that the mesh is also locally smooth. The choice to generate the reconstruction stencil using cell face neighbors is motivated by the will to ensure an almost constant value of β.

Geometric symmetry of the stencil is ensured by associating a number of neighbors equal to the actual number of neighbors of a face. Such a property provides a smooth variation of the truncation error. Therefore, it is possible to generate what we call face neighbor layers with the following algorithm. First, for each root cell control volume V i , all face neighbors are added to the stencil C i generating the first face neighbor layer. If further neighbor layers are needed, the previous process is repeated for each cell added in the previous iteration (ensuring unique indices in C i ) until the stencil contains a number of face neighbor layers equal to the polynomial order p of the scheme.

Regarding the LSQ schemes that will be analyzed in this work, we will restrict the study to those with polynomial degree such that p = 1, 2, 3 and 4. The resulting schemes are denoted as LSQp, with p = 1, 2, 3 and 4. Moreover, the effect of the weighting method in the least-squares problem resolution will be studied for p = 3. The subsequent weighted scheme will be referred as LSQ3W.

Least-squares stencils and face neighbor layers for schemes LSQp with p = 1, 2 and 3 are represented for twodimensional problems in Fig. 2. These figures allow to identify the stencil's geometry of some of these schemes in twodimensional meshes.

To better understand the value of the parameter β with the different combinations of polynomial order p and number of dimensions N d , Table 1 shows the number of neighbors of a given stencil for the different LSQ schemes and Table 2 represents the value of the β parameter as a function of the dimensions of the problem and the different LSQ schemes. In these tables, it is possible to observe that the parameter β increases with the number of dimensions of the problem for a given polynomial order. It can even reach values higher than β > 3 in some three-dimensional cases. With such high values of β, it could be interesting to reduce by one the neighbor layer index used to build the k-exact scheme. However, to remain consistent for all spatial dimensions, we decided not to change the index of the neighbor layer used to build the schemes, independently of the number of dimensions. It is also worth noting that reducing by one the index of the neighbor layer yields always β < 2 in the cases presented in Table 2.

k-exact schemes for structured meshes of tensor-product elements

Recent developments made by [START_REF] Arbogast | Finite volume weno schemes for nonlinear parabolic problems with degenerate diffusion on non-uniform meshes[END_REF] allow to obtain multidimensional k-exact schemes using a fully constrained problem for structured non-uniform meshes made up of tensor elements (quadrilaterals and hexahedrons). To develop these methods, a structured mesh with principal directions xy will be supposed. The extension to one and three-dimensional problems is straightforward. These schemes use the following linear mapping

x y = x i y i + J • ξ η , ( 50
)
where (x i , y i ) are the coordinates of the root cell centroid and

J = h i 0 0 h i . ( 51
)
The polynomial basis used by the schemes belongs to the polynomial basis of tensor product one dimensional polynomials of degree p

{ (ξ )} = ⎧ ⎨ ⎩ N d i=1 (1, ξ 1 i , . . . , ξ p i ) ⎫ ⎬ ⎭ . ( 52
)
Please note that the 0 th degree of freedom is included polynomial basis, opposed to the choice made in k-exact leastsquares schemes where this polynomial basis is supposed equal to the cell average (see Section 3.1). Moreover, the maximum monomial degree is N d p. The cardinal of the polynomial basis is redefined for tensor product polynomials as

K = (p + 1) N d . ( 53
)
Regarding the neighbors used to obtain the polynomial coefficients, as the 0 th degree polynomial basis is unknown, the root cell i will also be added to the central stencil C i . These k-exact schemes build their stencils using cell vertex neighbors instead of cell face neighbors. For the sake of brevity, the discussion of the generation of the vertex neighbors based stencils is avoided as it can be straightforwardly understood from the generation face neighbors based stencils discussed in Section 3.1.4. The number of neighbors M = |C i | for a given cell vertex neighbor layer index N v in structured meshes is

M = (2N v + 1) N d . ( 54
)
In the latter expression it is possible to observe that the condition K = M (which is needed to ensure a fully-constraint problem) imposes the following condition for even values of p

N v = p 2 . ( 55
)
Therefore, a given tensor product polynomial of even degree p can be computed without using least-squares methods in a structured mesh by using the correct number of vertex neighbor layers. Odd degree polynomials can also be considered, although one needs to consider a different family of polynomials and face-based stencils as explained in [START_REF] Arbogast | Finite volume weno schemes for nonlinear parabolic problems with degenerate diffusion on non-uniform meshes[END_REF]. Therefore, we particularize our study for even degree polynomials. In these cases, the least-squares problem defined in Section 3.1 can be rewritten as a linear system as M = K . For each neighbor cell in the C i the following constraint equation is used to build the matrix A ∈ R K ×K as

∀ j ∈ C i K γ =1 a γ |V j | V j φ γ (ξ ) dξ dη = u j . ( 56
)
The expression of the terms of matrix A is

A nm = 1 |V n | V n φ m (ξ ) dξ dη . ( 57
)
Remark 3. In this work, all the terms from matrix A are computed using a numerical integration quadrature which integrates exactly all the polynomial basis up to machine precision. Compared to the LSQ schemes, these methods require higher order numerical quadrature if one wants to integrate exactly all the polynomial basis. For example, the maximum polynomial degree is pN d , which has a value of 12 for a 4 th degree polynomial.

Meanwhile, the right-hand-side vector b ∈ R M is equal to the cell averages of all the cells present in C i . The polynomial coefficients X ∈ R K can then be obtained by building the inverse matrix A -1 and using the previous definition of the right-hand-side vector: In uniform structured meshes these schemes are equivalent to the finite-volume k-exact schemes for uniform structured meshes studied in [START_REF] Zhang | On the order of accuracy and numerical performance of two classes of finite volume WENO schemes[END_REF]. When used in the discretization of linear equations, the schemes presented in this section are equivalent to those described in [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF].

X = A -1 b . (58)
It is important to note that these schemes have only been developed for uniform and non-uniform structured meshes and therefore, they will only be used for comparison purposes in those cases. Remark 4. The extension of the S schemes to quadrilateral or hexahedral arbitrarily-oriented structured meshes or unstructured grids is not straightforward due to the difficulty of choosing a suited reference frame to compute the polynomial. In our studies, we have observed that the numerical properties such as dissipation, dispersion and stability of the S schemes are heavily dependent on the reference frame chosen to generate the transformation to the reference space. Moreover, quadrilateral or hexahedral unstructured grids might have different numbers of vertex neighbors for a given layer index than structured meshes. This implies that one needs to reformulate the polynomial coefficients computation problem in these cases. Several choices can be made to deal with those issues such as recasting the fully constrained problem into a least-squares problem by modifying the polynomial basis, adding more cell face/vertex neighbors, etc. These choices will be analyzed in future studies.

Definition of the S schemes

In this study, the k-exact interpolation schemes described in Section 3.2 will be denoted with the S symbol. We particularize our study for S schemes with tensor-product polynomial degrees p = 2 and 4. The resulting schemes are denoted as S2 and S4 respectively. Fig. 3 represents the stencil and neighbor cells used to build the interpolation operators of schemes S2 and S4. Table 3 shows the number of neighbors M used by S2 and S4 as a function of the number of spatial dimensions of the problem. This table can be compared with the number of neighbors of LSQ schemes showed in Table 1, to assess the equivalence in terms of number of neighbors (which are directly related to the computational cost) between LSQ and S schemes. The S2 scheme uses same number of neighbors than the LSQ1 scheme in one-dimensional problems, although for twodimensional and three dimensional problems its number of neighbors is much closer to that of the LSQ2 scheme. Regarding the S4 scheme, it presents the same number of neighbors than the LSQ2 scheme in one dimensional problems. Meanwhile, in two-dimensional and three three-dimensional problems, the S4 scheme uses a very similar number of neighbors than the LSQ3 and LSQ4 schemes respectively.

Quadrature points and conserved variables face averages

This section is devoted to the analysis of the number of quadrature points needed to compute the face average fluxes in the FV formulation using tensor-product elements.

Gauss-Legendre quadratures [START_REF] Epperson | An Introduction to Numerical Methods and Analysis[END_REF][START_REF] Burden | Numerical Analysis[END_REF] allow the integration in a given face with index of a tensor-product element of a tensor-product polynomial of degree p by using a quadrature such that

= 1 p = 2 p = 3 p = 4 1 D 1 1 1 1 2 D 1 2 2 3 3 D 1 4 4 9 |G | = ceil p + 1 2 N d -1 (59)
where |G | is the number of quadrature points of face . The use of Gauss-Legendre quadratures allows to numerically compute the polynomial basis integrals along the faces as

1 |S | S φ γ (ξ ) dξ dη = g∈G w g φ γ ξ g ( 60
)
where ξ g denotes the coordinates of quadrature points in the reference space and w g are the Gauss-Legendre quadrature weights which also include the isoparametric Jacobian transform determinant of the given face. It is important to note the polynomial basis group of LSQ schemes is contained in the tensor-product polynomial basis of S schemes. Eq. ( 60) holds true provided that the determinant of the isoparametric transformation of the given tensor-product element face is constant in all of its points. If this property is not fulfilled (which might be the case in hexahedral elements), numerical errors will be introduced in the evaluation of Eq. ( 60) through Gauss-Legendre quadratures. In these cases, Eq. ( 60) has to be rewritten as an approximation when evaluated in an arbitrary face of a hexahedron. These numerical errors could be reduced by increasing the number of quadrature points, knowing that, in the worst case scenario, the metric terms of an arbitrary face are a 2 nd degree polynomial [START_REF] Naff | Shape functions for velocity interpolation in general hexahedral cells[END_REF] for hexahedrons.

When dealing with non-linear fluxes, the number of quadrature points is usually obtained from Eq. ( 59) (see [START_REF] Tsoutsanis | Extended bounds limiter for high-order finite-volume schemes on unstructured meshes[END_REF]) and it is shown in Table 4. It is important to note that, with non-linear fluxes, this choice of quadrature points will produce under-integration as in non-linear cases the flux polynomial representation is a product of the conserved variables. This can cause additional numerical errors and/or aliasing effects as described in [START_REF] Manzanero | Insights on aliasing driven instabilities for advection equations with application to Gauss-Lobatto discontinuous Galerkin methods[END_REF].

Gradient computation

In order to compute the high-order viscous face average fluxes as described in Section 2.3 and Eq. ( 15), one needs to compute the gradients of the polynomial representative of the conserved variables at each quadrature point. To do so, the gradients of the polynomial P i of each conserved variables at each cell i are computed in the reference space (ξ, η) yielding ∇ δ,± ξ u. The polynomial gradient in the reference space is then extrapolated to the quadrature points of each face of the cell by transforming the quadrature points coordinates into the reference space using the Jacobian matrix J i associated to each cell. Next, this gradient in the reference space is transformed to the physical reference frame following the expression from [

2] ∇ δ,± u g = J -1 T ∇ δ,± ξ u g . ( 61
)
At last, the gradient at each quadrature point is used to build the face average viscous flux term using Eq. ( 15) and Eq. ( 16).

Order of accuracy

This section is devoted to the analysis of the order of accuracy of LSQ and S schemes via the analysis of the truncation error and the k-exact properties of the reconstruction polynomials.

Linear advection spatial truncation error analysis in uniform meshes

The spatial truncation error is defined as the difference between the discretized and analytical equation. Truncation error is computed in the linear advection equation [START_REF] Vasconcelos | A very high-order finite volume method based on weighted least squares for elliptic operators on polyhedral unstructured grids[END_REF] neglecting time integration errors and supposing planar faces with constant metrics as

E i = 1 |V i | ∈L i u -u δ c • n i, |S | , ( 62
)
where the numerical face average u δ is computed as

∂x 3 i + O h 3 2 nd CLSQ1 -0.083ch 2 ∂ 3 u ∂x 3 i + O h 4 2 nd ULSQ2 0.283ch 3 ∂ 4 u ∂x 4 i + O h 4 3 rd CLSQ2 -0.114ch 4 ∂ 5 u ∂x 5 i + O h 6 4 th ULSQ3 0.058ch 4 ∂ 5 u ∂x 5 i + O h 5 4 th CLSQ3 0.058ch 4 ∂ 5 u ∂x 5 i + O h 6 4 th ULSQ3W 0.017ch 4 ∂ 5 u ∂x 5 i + O h 5 4 th CLSQ3W 0.017ch 4 ∂ 5 u ∂x 5 i + O h 6 4 th ULSQ4 0.331ch 5 ∂ 6 u ∂x 6 i + O h 6 5 th CLSQ4 0.139ch 6 ∂ 7 u ∂x 7 i + O h 8 6 th u δ = f 1 + sign(c • n i, ) 2 u δ,-+ 1 -sign(c • n i, ) 2 u δ,+ + (1 -f ) 1 2 u δ,-+ u δ,+ . ( 63
)
The first term of the expression is representative of an upwind numerical flux while the second term is representative of a central numerical flux. We recall that upwind schemes (denoted with a starting U symbol) and central schemes (denoted with a starting C symbol) use f = 1 and f = 0 respectively in all quadrature points. Eq. ( 62) can be further simplified in uniform meshes in which the characteristic edge size h i is constant h i = h as

E i = 1 h ∈L i u -u δ c • n i, . (64)
The latter expression allows to obtain the order of accuracy and the leading term of the truncation error in uniform meshes by using face and cell average Taylor expansion series of the conserved variable u. In addition, numerical values u δ can be computed by integrating analytically the polynomials resulting from the LSQ reconstruction.

Table 5 shows the leading truncation error term and order of accuracy of different k-exact least-squares schemes in the one-dimensional linear advection problem using uniform structured meshes. In these results it is possible to observe that all ULSQp schemes show at least p + 1 order of accuracy, meaning that they respect the k-exact reconstruction condition. Moreover, the CLSQp schemes also show at least p + 1 order of accuracy. However, it is important to note that the central variants of the LSQ schemes yield p + 2 order of accuracy if p is even. Furthermore, central schemes always present truncation error terms which have even exponents with h. These assessments are commonly observed in the finite-difference and finite-volume literature. At last, compared to the default LSQ3 scheme the weighting distance LSQ3W scheme allows to reduce the amplitude of the leading truncation error term by a factor close to three.

Regarding the comparison of one-dimensional S and LSQ schemes, the S2 scheme has the same stencil size as the LSQ1 scheme while the S4 scheme has the same stencil size as the LSQ2 scheme. Table 6 shows the leading truncation error terms and order of accuracy of these two schemes. In these results it is possible to observe that the order of accuracy of the S2 scheme variants is similar to that of the LSQ2 scheme variants. Nonetheless, the latter has a bigger stencil and its leading truncation error term is bigger. Moreover, the S4 scheme has the same order of accuracy than the LSQ4 schemes although the LSQ4 scheme has a larger stencil and the leading truncation error term is bigger. At last, the central variants of these schemes have similar properties to the LSQ schemes built with even polynomial order, meaning that their truncation error only shows terms with even exponent with h and that the central variants order of accuracy is one order higher accurate than their upwind counterparts.

The truncation error analysis of the two-dimensional linear advection equation yields similar results as the ones shown in Table 5 and Table 6 regarding the order of accuracy of the LSQ and S schemes. Nonetheless, LSQ schemes present additional cross-derivatives xy terms present in its expression of the truncation error. To give an example of this issue the truncation error of the LSQ2 scheme in the two-dimensional linear advection equation is

E ULSQ2 (h) = c x h 3 19 84 ∂ 4 u ∂x 4 i, j + 1 14 ∂ 4 u ∂x 2 ∂ y 2 i, j + c y h 3 19 84 ∂ 4 u ∂ y 4 i, j + 1 14 ∂ 4 u ∂x 2 ∂ y 2 i, j + O h 4 , ( 65
)
where the presence of cross-derivative xy terms can be clearly observed. This fact is exacerbated as the order of accuracy of the LSQ scheme increases. This allows to conclude that the LSQ schemes in two-dimensional uniform meshes can present increased numerical errors compared to the one-dimensional LSQ schemes counterparts in the two-dimensional linear advection equation.

It is important to note that the S schemes do not present this cross-derivative terms in the expression of their truncation error. For example, the truncation error of the two-dimensional S2 scheme reads

E US2 (h) = - c x h 3 12 ∂ 4 u ∂x 4 i, j - c y h 3 12 ∂ 4 u ∂ y 4 i, j + O h 4 . ( 66
)
It can be observed in the previous expression that S schemes recover the one-dimensional truncation error expression if the solution is uniform along a given principal direction of the structured mesh.

For the sake of brevity, we have decided to not include the leading truncation error terms of the rest of the twodimensional LSQ schemes, although their order of accuracy is the same as that of their one-dimensional counterparts.

Two dimensional order of accuracy in non-uniform and unstructured meshes

The order of accuracy of least-squares schemes in unstructured meshes is assessed numerically. The p-th order reconstruction of a given scheme is evaluated by the definition of the order of accuracy or k-exact property. Such an approach is inspired on the work of [START_REF] Petrovskaya | The accuracy of least-squares approximation on highly stretched meshes[END_REF][START_REF] Tsoutsanis | Stencil selection algorithms for weno schemes on unstructured meshes[END_REF]. Previous works studied the k-exact property of the least-squares polynomial and polynomial gradient reconstruction particularized at the cell centers. Nonetheless, the p-th order reconstruction property is extended to analyze the polynomial reconstruction property at quadrature points and face averages reconstructions. A test-polynomial T is defined as

T (x) = (ax + by) p , ( 67
)
where a, b ∈ R and p is the polynomial order of the test error function. In particular, the use of different coefficients a and b is justified to avoid the possible mesh principal direction x = y. The order of accuracy of a polynomial k-exact scheme can be computed from the greatest value of p which ensures that the reconstructed polynomial (built from T ) at each cell i follows

P i (ξ ) = T (x (ξ )) .
(68) A k-exact scheme which uses a polynomial of degree p should be capable of reconstructing this test error function with the same or less order p ≤ p at each point of the cell. In this case the scheme is said to be p + 1 order accurate. As explained beforehand, the focus is on the k-exact properties of schemes with regards to the interpolation to the quadrature points and the face average reconstructions. To study these properties, two different relative errors will be defined. The first is defined at each face l of the domain and it is computed as

E l,u = max P --T T , P + -T T , ( 69
)
if the studied scheme is upwind and

E l,c = 1 2 P -+ P + -T T , ( 70
)
for central schemes. In the latter expressions, P -and P + refer to the face average of the polynomial build using the central stencil of the face owner and face neighbor cells respectively. These error estimators are related to the numerical error of the computation of the face average conserved variables.

The second relative error used to assess the k-exact properties of the spatial schemes is defined at each quadrature point for upwind schemes as

E g,u = max P -ξ g -T x ξ g T x ξ g , P + ξ g -T x ξ g T x ξ g , ( 71
)
and for central schemes as

E g,c = 1 2 P -ξ g + P + ξ g -T x ξ g T x ξ g . ( 72
)
This second error estimator is related to the error in the computation of the fluxes at the quadrature points.

Remark 5. It is worth noting that special care should be taken when defining the polynomial test error function to ensure that the minimum absolute value of the polynomial in the mesh is bigger than the machine rounding-off errors by an important margin to avoid divisions by zero in the computation of the numerical error. This is done by displacing the xy reference frame of the test error function.

Remark 6. Although the test error function is not periodic in the meshes that will be used in this study, the error estimators computation can still be performed with meshes that present periodic boundary faces by ensuring that the periodic ghost cells are initialized with the cell average value of the test error function computed with the geometry of each periodic ghost cell, avoiding the face and quadrature points data exchange due to the mesh periodicity condition. The only drawback of this approach is that the polynomial of the neighbor cells of the faces located in the periodic boundaries cannot be used to compute error estimators. This is due to the fact that the assessment of this polynomial requires periodic data exchange, and as the test error function is not periodic the polynomial reconstruction shows wrong values.

Remark 7. The previous two-dimensional test can be easily expanded to three dimensional analysis by adding the dependence on the z coordinate to the error function. Moreover, it can also be used in one-dimensional cases removing the dependence on the y coordinate.

Three two-dimensional meshes (depicted in Fig. 4) defined in the domain x ∈ [0, L] 2 with periodic boundary condition on all boundaries, are used to compute the previously described error estimators.

1. A uniform structured mesh of edge size h 2. A uniform structured mesh of initial edge size h with a random mesh-nodes perturbation following the expression:

x = x + hK rand(-1, 1) ∀ x ∈ (0, L) 2 , ( 73
)
where K = 0.25 is a constant and the rand function is a random floating point number generator in the domain [-1, 1] ∈ R. Note that the points which lay on the boundaries of the domain are not modified to respect the periodicity conditions. 3. A fully quad-unstructured mesh with a characteristic edge size of h.

Table 7 reports error estimators in the uniform mesh configuration (Mesh 1) for all the schemes studied in this work and as a function of the test error function order p. The results show that, as presented in Section 4.1, all schemes are k-exact regarding the interpolation to the quadrature points and the face average reconstruction as they are capable of approximating a test error function with p = p up to machine precision rounding-off errors. Moreover, some central scheme variants are k+1-exact in the sense of the face average reconstruction. This property is characteristic of schemes with even polynomial degree p and was also observed in the truncation error analysis (Section 4.1). It is important to note that the previous k+1-exactness is not observed when analyzing the error estimator at the quadrature points. Moreover, every central scheme shows lower error values at the quadrature points than their upwind counterpart. As it will be shown afterwards, this k+1-exactness of some central schemes is lost in non-uniform and unstructured meshes due to the loss of truncation error canceling. Table 8 presents the values of the different error estimators in Mesh 2 configuration for every scheme studied in this work and as a function of the test error function order p. As for uniform structured meshes, all upwind LSQ schemes are k-exact schemes regarding the interpolation to the quadrature points and the face average reconstruction, as they are capable to approximate a test error function with p = p up to machine precision rounding-off errors in this non-uniform mesh configuration. Moreover, it is possible to observe that the k+1-exactness property of certain central schemes regarding the face average reconstruction is lost in this non-uniform mesh configuration. Nonetheless, it must be noted that the face average error estimator for test error functions such that p = p + 1 is much lower for the central schemes with even polynomial order p with respect to their upwind counterpart. This difference is less pronounced for least-squares schemes built with odd degree polynomials. Table 9 shows the different error estimators in the Mesh 3 configuration. The results are very similar compared to those showed in the previous analysis of Table 8 therefore, for the sake of brevity, we avoid its discussion.

Spectral properties of finite volume schemes

In the following, spectral properties of the numerical reconstruction of the conservation equation ( 10) will be addressed.

Let us assume a linear flux operator F (u) = cu, where c is the convective speed. In particular, Eq. ( 10) discretized in a mesh of N V cells, may be recast to a dynamical system [START_REF] Haider | A high-order interpolation for the finite volume method: the coupled least squares reconstruction[END_REF] of

du δ - dt = J u δ - , (74)
where J ∈ R N V ×N V is a matrix which depends on the upwinding parameter f of each numerical scheme and it is constructed using the interpolation polynomials. Therefore, it uniquely depends on geometric properties of the mesh discretization. In addition, u δ

-∈ C N V = u δ 1 , u δ 2 , . . . , u δ i , . . . u δ N V
is a vector containing cell average values of the solution at each cell. Von Neumann analysis takes advantage of the particular structural properties of the matrix J , i.e. circulant, orthogonal, etc. For instance in the Toeplitz case, its symbol provides access to the spectrum and the pseudo-spectrum (see [START_REF] Reichel | Eigenvalues and pseudo-eigenvalues of Toeplitz matrices[END_REF]). Nonetheless, the general case is slightly more complex. Toeplitz property is lost whenever the underlying numerical triangulation is not uniform, in such a case entries of the matrix J cease to be in a particular ordering. Integration of the dynamical system Eq. ( 74) can be analytically performed with a semi-group propagator, i.e. the exponential matrix, so one ends up with a discrete map as follows

u δ -(t n+1 ) = R(τ , J ) u δ -(t n ) , (75)
where R(τ , J ) = R δ = e τ J is the exponential of the matrix J , t n is the solution at a given discrete time and t n+1t n is the time step τ . Similarly, one could choose an approximation of the actual semi-group, i.e. a Runge-Kutta multistage integrator R δ with a given number of stages and order. Numerical examples are carried out with a RK33 temporal discretization method. This method allows to compute R δ as

R δ = I + τ J + τ 2 2 J 2 + τ 3 6 J 3 . ( 76
)
Furthermore, let us remind that the semi-group of the linear advection equation with periodic boundary conditions and initial condition: u(x, t = 0) = û e Iκ x cos θ +y sin θ , (77) where κ is the spatial wavenumber, is analytical with the following expression: R = e -Iωτ I , such that ω = κ c x cos θ + c y sin θ is the analytical temporal wavenumber. In this work we will only consider advection velocity vectors such that ζ = π/4. Remark 8. If the exponential operator is used to build the matrix operator R, the results will only take into account spatial discretizations defects. On the contrary, the consideration of Runge-Kutta multistage integrators to build R δ imply that the numerical solution will include both spatial and temporal defects.

Von Neumann analysis in structured uniform meshes

Von Neumann analysis on structured uniform computational domains with periodic boundary conditions takes advantage of the circulant property of the linear operator J , i.e. the action matrix-vector product is equivalent to the convolution of two vectors, thus the Fourier transform is a pointwise multiplication. These properties can only be obtained if the mesh is uniform, the stencils of each cell have the same topological shape and the boundary conditions are periodic. In addition, it is possible to show the equivalency between an analysis of Eq. ( 75) in the finite-volume formulation and in the finitedifference sense. Thus, Eq. ( 75) can be recast into:

u -(t n+1 ) = R δ u -(t n ) , ( 78
)
where u -∈ C N V represents the point-value solution at the cell centroids. As matrix J is circulant, so is R δ . Von Neumann analysis aims to determine the deviation of the numerical solution to Eq. ( 78) and the solution to the advection equation.

Let us consider the solution to the advection equation u(x, t) = û e Iκ(x cos θ +y sin θ )-Iωt , (

which is injected to Eq. ( 78), to obtain the value of the numerical temporal wavenumber of each cell

ω δ i = ω δ = I τ ln ⎛ ⎝ j R δ i, j e Iκ (x j -x i ) cos θ +(y j -y i ) sin θ ⎞ ⎠ = I ln Q τ , ( 80
)
where ω δ is the numerical temporal wavenumber and Q is an amplification factor. It is emphasized that this numerical temporal wavenumber has the same value for each cell due to the properties of the R δ matrix. This uniqueness will be lost in unstructured or non-uniform meshes as it will be shown afterwards.

To assess the dissipation and dispersion errors, a non-dimensional temporal wavenumber error function is defined as

D(κh, CFL, θ) = 1 ω |Re ω δ -ω | + I|Im ω δ -ω | , (81)
where CFL = cτ h . The comparison of the numerical temporal wavenumber and its analytical value is commonly used to assess the dissipation (Im(D)) and dispersion (Re(D)) errors of the numerical schemes [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF][START_REF] Vincent | Insights from von Neumann analysis of high-order flux reconstruction schemes[END_REF][START_REF] Vermeire | On the properties of energy stable flux reconstruction schemes for implicit large Eddy simulation[END_REF][START_REF] Alhawwary | Fourier analysis and evaluation of dg, fd and compact difference methods for conservation laws[END_REF][START_REF] Trojak | Effect of mesh quality on flux reconstruction in multi-dimensions[END_REF]. The one-dimensional case can be particularized by imposing θ = ζ = 0 and by using a one-dimensional mesh.

If temporal discretization errors are neglected (by assuming exponential temporal integration), it is possible to obtain an analytical expression of the non-dimensional error function. For example for the one-dimensional ULSQ1 scheme this error function can be written as

D(κh) = - (κh) 2 12
+ I (κh) 3 8

+ O (κh) 4 . ( 82
)
In the latter formula, it is possible to differentiate the dissipation (imaginary part) and dispersion errors (real part) of the given spatial discretization scheme. Moreover, the leading error term has an exponent which is of the same order of the truncation error term in Table 5. Therefore, the truncation error terms with odd exponent produce dissipation and the ones with even exponent produce dispersion. It is also important to note that central schemes do not present dissipation errors when using analytical integration methods due to the fact that their matrix J is anti-symmetric. This implies that the eigenvalues of their matrix J are purely imaginary. Hence, dissipation of central schemes can uniquely appear due to the temporal numerical integration method. This can also be observed by analyzing the truncation errors of central schemes showed in Table 5.

One-dimensional von Neumann analysis

The one-dimensional von Neumann analysis can be performed by supposing a plane wave and a uniform periodic mesh along the x direction with constant cell size h. This analysis will be performed using the cells per wavelength parameter λ/h = 2π κ (similar to the points per wavelength parameter used in [START_REF] Vanharen | Revisiting the spectral analysis for high-order spectral discontinuous methods[END_REF]) instead of the more broadly used non-dimensional spatial wavenumber parameter κh.

Dissipation & dispersion

Dissipation and dispersion of finite volume schemes are analyzed by considering the real and imaginary parts of D. The study will be performed with three different CFL values. The lowest CFL = 0.1 is chosen to ensure that the spatial discretization errors are dominant; the highest CFL = 1.0 is chosen to evaluate the numerical errors of simulations dominated by temporal discretization errors. At last, the case with CFL = 0.5 will also be presented to assess the interactions between the spatial and temporal discretization errors.

Fig. 5 and Fig. 6 display the dispersion error of the upwind and central schemes respectively. The figures show that for all CFL values the upwind and central scheme variants present a similar trend. Let us start by considering dispersive errors where n is an exponent ranging from 2 to 6 and are used to assess the leading term order of the dispersion errors. where n is an exponent ranging from 2 to 6 and are used to assess the leading term order of the dispersion errors.

at CFL = 0.1. The dispersion error evolves as a power law with the number of cells per wavelength which is representative of even exponents of the truncation error terms analyzed in Section 5.1. Nonetheless, it is worth noting that such a trend appears at values of λ/h > 10, whereas at lower values most least-squares schemes do not exhibit this behavior. Let us stress out that even if the dispersion errors follow the correctly predicted trend for high values of λ/h, error values of most LSQ schemes are higher than those of the corresponding structured (S) schemes. An exception to this behavior is observed when analyzing the LSQ3W scheme which shows lower dispersion errors than the S2 scheme with the same order of accuracy.

For high CFL values most upwind and central schemes show similar dispersion errors and a 4 th order trend with h. This can be explained due to the temporal discretization errors, which in our case, due to constant CFL condition, behave as O τ 3 ∼ O h 3 . As dispersion errors are caused by terms with even exponent truncation errors with h, the dispersion error will behave as the 4 th order truncation error of the temporal discretization method O τ 4 ∼ O h 4 .

It is also interesting to comment the special behavior of the LSQ1 and S2 schemes as the CFL increases. Scheme LSQ1, does not display a 4 th order trend with h because this scheme has a 2 nd order in space which is lower than the order in time, therefore its dispersive errors are not influenced by temporal discretization. On the other hand, the S2 scheme shows a decrease of its dispersion errors as the CFL increases. This fact is exacerbated in its central variant scheme CS2, in which we observe a very high-order trend of the dispersion errors with the cell size parameter. The reason for such superconvergence behavior for CFL = 1 was not identified. This superconvergence has not been observed when combining the S2 scheme with other RK time integrators at CFL = 1. Therefore, further studies are needed to better understand this observation. Fig. 7 displays the evolution of the dissipation error with respect to λ/h for the upwind schemes. Dissipation error of each upwind scheme shows a similar trend to dispersion error. For low CFL values, dissipation is dominated by spatial where n is an exponent ranging from 2 to 6 and are used to assess the leading term order of the dissipation errors.

discretization errors, whereas at large CFL the opposite holds. At low CFL, a power-law trend of the dissipation error is observed for high values of the number of cells per wavelength. The exponent of the power law corresponds to the odd exponent of the truncation error of each scheme. For low values of λ/h the least-squares schemes underperform compared to the US2 and US4 schemes. In particular, we remark that the ULSQ3W over-performs every other least-squares scheme and the S2 scheme for each λ/h, which shows that the addition of distance weighting to the least-squares schemes may improve the dissipation properties of these schemes. As the CFL number increases temporal discretization errors dominate for high resolved waves. Under those circumstances, every scheme displays a O h 3 error trend, similar to what it has been observed with the dissipation of central schemes.

The dissipation errors of central schemes are not shown due to the fact they are just given by the 3 rd numerical error of the RK33 time integrator, due to the choice of constant CFL. Such behavior is equally observed for high values of the CFL number and high values of the cells per wavelength for upwind schemes, i.e. in these cases the temporal discretization error dominates.

Temporal stability

This section is devoted to the study of the L 2 -stability of the linear advection equation, thus maximum CFL values will be reported for every spatial numerical scheme presented in this study. Maximum CFL is assessed from Eq. (80), by taking into account an amplification factor lower than one

Q ≤ 1 , ( 83
)
which is equivalent to stating that the absolute part of every eigenvalue of operator R δ is lower than one.

The eigenvalues of R δ vary with the CFL and λ/h. Therefore, for each value of κh there is a given CFL which ensures linear stability. The lowest value of these CFL conditions, denoted as CFL MAX and is the maximum CFL value and ensuring linear stability for a given scheme.

Table 10b reports the maximum CFL number ensuring linear stability for the S2 and S4 schemes in a uniform mesh using the RK33 time integration method. Let us remark the validity of the approach which has been validated via a crosscomparison of the maximum CFL of US4 scheme, previously reported by [START_REF] Liu | Wls-eno: weighted-least-squares based essentially non-oscillatory schemes for finite volume methods on unstructured meshes[END_REF]. It can be observed that S2 schemes present higher CFL MAX than S4 schemes. Additionally, structured central schemes display a lower CFL MAX than their upwind counterparts. Such a behavior is expected as the central schemes have much lower dissipation values than the upwind schemes and therefore, the linear stability can only be ensured for lower values of the CFL number.

Table 10a reports the maximum CFL number ensuring linear stability for the different one dimensional least-squares schemes in a uniform mesh using the RK33 time integration method. The linear stability CFL criterion of these schemes behaves differently to that of the S2 and S4 schemes and does not follow a coherent trend with respect to their order of accuracy and stencil size. Moreover, one may also observe that for most LSQ schemes, the central variant displays a higher value of the CFL MAX than the upwind variant.

Two-dimensional analysis

The two-dimensional von Neumann analysis is carried out with the same tools developed for the one-dimensional case.

However, the dissipation and dispersion maps depend on the wave angle θ and the advection velocity vector. For these studies, the convection velocity angle ζ is set as ζ = 45 • . This convection velocity angle ensures that the dissipation and dispersion spectrums are computed in the preferential direction of propagation for quadrilateral structured uniform meshes (see [START_REF] Voth | Generalized Fourier analyses of the advection-diffusion equation-part II: two-dimensional domains[END_REF][START_REF] Vichnevetsky | Fourier Analysis of Numerical Approximations of Hyperbolic Equations[END_REF]). Moreover, as dissipation and dispersion depend on the wave angle θ , these parameters need to be analyzed for each wave angle. The temporal discretization time step is fixed such that cτ /h = 0.05, as we have observed that the influence of the temporal discretization methods in the von Neumann analysis is similar to the one-dimensional analysis. We will consider that the computation of the CFL MAX condition of the schemes in two-dimensional problems is beyond the scope of this work, although the methods developed in the previous sections are still valid to compute CFL MAX . The interested reader is referred to [START_REF] Trojak | Effect of mesh quality on flux reconstruction in multi-dimensions[END_REF] for further comments regarding the extension of the CFL MAX condition to two-dimensional problems. Fig. 8a show the contours of the dissipation for the upwind US2, US4, ULSQ2 and ULSQ4 schemes as a function of the number of cells per wavelength and the wave angle. These contours show clearly that dissipation is always minimum at θ = 45 • = ζ for all schemes. This behavior was also described in [START_REF] Lele | Compact finite difference schemes with spectral-like resolution[END_REF]. Moreover, the dissipation values show similar trends as in the one-dimensional studies. The ULSQ2 is more dissipative than the US2 scheme, although they have the same order of accuracy. The previous statement can also be applied when comparing the ULSQ4 and US4 schemes. At last, we observe that the ULSQ4 scheme shows very high values of the dissipation error compared to the other schemes for under-resolved waves in this two-dimensional analysis. This behavior was also noted in the one-dimensional dissipation analysis. Similar conclusions can be obtained regarding the dispersion errors as shown in Fig. 8b.

For comparison purposes, Fig. 9 shows the dissipation and dispersion errors of all two-dimensional upwind schemes at θ = 45 • . The results obtained are very close to those obtained in the one-dimensional study at low CFL number (Fig. 7 and Fig. 5), therefore we avoid its discussion. It is also observed that the two-dimensional central schemes behave similarly to their one-dimensional counterparts. This implies that the dissipation of these schemes is dominated by the temporal discretization errors and that their dispersion is equivalent to the dispersion of the upwind schemes for low CFL numbers. We avoid adding more two-dimensional comparisons for the sake of brevity. This is justified due to the fact that the one-dimensional conclusions regarding the numerical dissipation and dispersion can also be applied to the two-dimensional analysis for all the schemes.

Diffusion & dispersion errors in two-dimensional unstructured grids

The von Neumann analysis of Section 5.1 can be further expanded to study dissipation and dispersion errors in unstructured meshes. It is important to note that in non-uniform meshes, the J matrix is not Toeplitz nor circulant, therefore the properties that allowed to simplify the numerical error analysis in the von Neumann method cannot be used. This is due to the fact that circulant matrices' eigenvectors are Fourier modes and therefore, the dissipation and dispersion properties of a wave described by a Fourier mode are just given by a single eigenvalue of J . To extend the von Neumann analysis for unstructured meshes, we first define the cell average of a given two-dimensional wave function

u i (t) = E i ûi (t) , ( 84
)
where ûi (t = 0) = cte and E ∈ C N V ×N V is a diagonal matrix whose entries E i are defined as

E i = 1 |V i | V i e Iκ(x cos θ +y sin θ) dV . ( 85
)
The linear system in Eq. ( 75) is then expressed as

û -(t n+1 ) = E -1 R δ E û -(t n ) , (86)
where û -∈ C N V ×N V is a vector containing the values of ûi at each cell. If the initial vector û -(t n ) is an eigenvector of E -1 R δ E, then the solution of the previous equation is given by û

-(t n+1 ) = û -(t n ),
where is an eigenvalue. This is typically the case for the von Neumann analysis in uniform meshes. However, in non-uniform meshes the previous statements do not hold true and therefore, the vector û -(t n+1 ) will not only depend on a single eigenvalue, but on a combination of them. This also implies that û -(t n+1 ) won't be uniform even if the initial û -(t n ) is uniform. Therefore discrete system response described in Eq. (86) will present a transient regime. As a result, the dissipation and dispersion analysis in non-uniform meshes has to be performed using global estimators under asymptotic behavior assumptions. To do so, one can build a numerical estimator of the imaginary (dissipation) and real (dispersion) parts of the temporal wavenumber as

Im(ω δ n ) ≈ 1 τ ln û -(t n+1 ) û -(t n ) , (87) and Re
(ω δ n ) ≈ - 1 τ arg û -(t n+1 ), û -(t n ) . (88)
A posteriori analyses of the dissipation and dispersion estimators have confirmed that the values of the dissipation and dispersion estimators at n = 0 yield accurate predictions of the numerical errors observed in the simulation of the linear advection in two-dimensional non-uniform and/or unstructured meshes (see Section 6.3). This fact was also recently discussed in [START_REF] Fernandez | Non-modal analysis of spectral element methods: towards accurate and robust large-Eddy simulations[END_REF], where it was stated that the short-term numerical errors of Spectral Element Methods could be predicted by analyzing the dissipation estimator (Eq. ( 87)) at n = 0. However, it is worth noting that the dissipation and dispersion properties of badly-resolved wavelengths when increasing the iteration number n. This issue appears due to the initial high dissipation values of under-resolved waves and machine rounding-off errors. The previous observations give enough evidence that the analysis of the global estimator of the numerical temporal wave number can be made for the first iteration n = 0 at least for numerical schemes and meshes studied in this work. As explained earlier, these meshes present close-touniform cell sizes with moderate cell aspect ratio changes. It has been observed, a posteriori, that meshes with substantial cell aspect ratio variation exhibit a considerably different eigenvalue spectrum from that obtained with the meshes used in this work. Such eigenvalues distribution has an important impact the linear dissipation and dispersion transient behavior, possibly invalidating the conclusions of this work. At last, we stress out that these observations only apply to the mesh of this study. The conclusions obtained from this section, may not be the extrapolated to meshes with abruptly varying aspect ratios or bigger skewness factor or with non-linear spatial discretization schemes.

Under the previous assumptions, equations (87) and (88) (particularized for n = 0) can be introduced in the temporal wavenumber error function D (Eq. ( 81)) to compute the dissipation and dispersion errors of the global estimator of numerical temporal wavenumber.

It should be noted that the previous method is analogous to the von Neumann analysis in uniform meshes, as in those meshes there are no transient effects due to the properties of circulant matrices.

The previously discussed method can be applied to study dissipation and dispersion properties of unstructured meshes. It is important to note that this method cannot be applied to study arbitrary waves, as the waves need to respect the following periodicity constraints of the domain

κ cos θ = 2πn x L and κ sin θ = 2πn y L , ( 89
)
where n x and n y are positive integer numbers. In uniform meshes, this particularization was not needed as one can assume that the meshes have arbitrary domain and cell sizes which respect the periodicity condition. This leads to the conclusion that high-order least-squares schemes are robust regarding their accuracy and numerical properties when meshes lose their global uniformity. This is not generally the case for other linear reconstruction schemes, which might lose accuracy and spectral properties with mesh deformation or perturbation. Moreover, this reveals that most conclusions obtained in the one-dimensional analysis remain valid. This allows to conclude that the one-dimensional analysis of least-squares schemes is representative of their behavior in the non-uniform two-dimensional meshes described in this work. Such meshes present a certain degree of cell size uniformity which, as described previously, does not substantially modify the eigenvalues spectrum with respect to that of uniform meshes. We omit central schemes from the comparison due to the fact that their J operator has eigenvalues with positive real part in the aforementioned configurations and therefore those schemes are inherently unstable.

Numerical experiments

In this section, numerical simulations will be performed to validate the previous analytical analysis in linear and nonlinear cases.

Error computation

To assess the numerical errors of the simulations the L ∞ norm based on the cell averages of the numerical will be used. In our study, this norm is defined as

L ∞ (u, t) = N V max i=1 |u δ i (t) -u i (t)| . ( 90
)
The computation of this error norm with different number of cells in the domain will allow to perform an a posteriori error norm convergence order study. The order of the error norm convergence should be consistent with the order of the truncation errors studied in Section 4.1.

Linear advection in one-dimensional mesh

In this section, the one-dimensional linear advection equation (Section 2.5.1) is solved using the LSQ and S based finitevolume schemes described in Section 3.1 and Section 3.2. In this analysis, the initial condition is defined by a smooth sine wave [START_REF] Alhawwary | Fourier analysis and evaluation of dg, fd and compact difference methods for conservation laws[END_REF]:

u(x, t = 0) = sin 2π x λ , x ∈ [0, L] . ( 91
)
In the previous initial condition, we fix L λ = 4, L = 1 and c = 1. Moreover, we define the non-dimensional parameter cT L = 1 where T is the wave temporal period in the domain. The domain is uniformly discretized with cell size h and periodic boundary conditions at x = 0 and x = L are applied. To assess the numerical properties of each numerical scheme, we run simulations at fixed CFL = cτ /h = 0.1, 0.5 and 1.0 while decreasing the cell size h (which implies implicitly a reduction of the time step τ ). The numerical errors will be analyzed at t = T using the L ∞ norm (Eq. ( 90)).

Fig. 12 shows the temporal evolution of the L ∞ norm of the upwind schemes as a function of the non-dimensional time t/T , different values of λ/h and the three CFL values mentioned. These results are representative of the combination of the dissipation and dispersion of the spatial and temporal discretization schemes. In this figure, it is possible to observe that error norm of all the spatial discretization schemes shows their predicted theoretical order of accuracy trend (as discussed in Table 5 and Table 6 for CFL = 0.1). However, as the CFL value increases, the higher order schemes show much higher errors due to the temporal discretization schemes, degrading their order of accuracy to a maximum of three as predicted in Section 5.1. To assess the difference spatial discretization schemes, we focus on the results with the lower CFL value. Among LSQ schemes, the ULSQ3W scheme outperforms all others in the studied range of λ/h, although the ULSQ4 (which has a bigger stencil) should yield better results than the ULSQ3W for λ/h > 80 due to its order of accuracy which is one order higher. The error values of the ULSQ3W are three times lower than that of its counterpart ULSQ3 schemes for well-resolved waves. This is consistent with the reduction of the magnitude of the truncation error shown in Table 5. It is worth noting that the ULSQ4 scheme behavior with badly resolved waves is closer to that of the lower order ULSQ1 and ULSQ2 schemes. This might indicate that ULSQ schemes with large stencils might need additional treatment (such as the addition of distance weighting methods described in Section 3.1.3) to better handle low resolution waves.

Regarding the S schemes, they present much lower error norm values than LSQ schemes with the same stencil or with the same order of accuracy. In particular, the US2 scheme, which is comparable to the ULSQ1 and ULSQ2 schemes in terms of stencil size and order of accuracy respectively, shows one more order of accuracy and five times less error than the ULSQ1 and ULSQ2 schemes respectively. Moreover, the US4 scheme, which can be compared to the ULSQ2 and ULSQ4 schemes in terms or order of accuracy, clearly outperforms both by large orders of magnitude. It is also worth noting that the S schemes outmatch the LSQ schemes in the prediction of the propagation of low resolution waves.

Similar conclusions can also be made regarding the errors obtained when using central schemes as shown on Fig. 13. However, due to the truncation error cancellation properties of the central schemes built with even polynomial order, it is possible to note that the CLSQ2, CLSQ4, CS2 and CS4 schemes clearly outperform their upwind variants. Nonetheless, the CLSQ1, CLSQ3 and CLSQ3W show similar results to their upwind variants.

Linear advection in two-dimensional mesh

To test the numerical properties of the finite volume schemes in 2D linear cases the two-dimensional linear advection equation is solved with the following initial wave condition used in [START_REF] Mengaldo | On the connections between discontinuous Galerkin and flux reconstruction schemes: extension to curvilinear meshes[END_REF] u(x, y, t = 0) = sin 2π

x λ cos 2π y λ , x ∈ [0, L] 2 , ( 92
)
where λ = L/4 is a given characteristic wavelength, L = 1 and c x = c y = 1. We define the non-dimensional parameter c x T L = 1 where T is the wave temporal period. To study this case, the three difference meshes discussed in Section 4.2 with characteristic cell size h will be used. Periodic boundary conditions are applied at each edge of the domain. To assess the numerical properties of each numerical scheme we run simulations at fixed c x τ /h = 0.05 while decreasing the cell size parameter h, so as to ensure that spatial discretization errors are dominant with respect to temporal discretization errors. Fig. 15a and Fig. 14b show the error norm for Mesh 1 configuration (uniform mesh) and the upwind and central schemes respectively. The results are very close to those shown in the one-dimensional analysis of Section 6.2, as was also predicted with the von Neumann analysis described in Section 5.1.2. Hence, the analysis of the aforementioned results is not included for the sake of brevity.

All the previous conclusions can also be applied to the error norms obtained using Mesh 2 and Mesh 3 configurations with the upwind schemes shown in Fig. 15a and Fig. 15b respectively. The results of central schemes are not presented as they are unstable in these mesh configurations. In these figures, it is possible to observe that LSQ schemes show similar error trends in non-uniform and unstructured meshes, maintaining the same order of accuracy as in structured meshes. This was also predicted in Section 5.2. S schemes results for Mesh 3 configuration are not presented due to the fact that they are only compatible with structured meshes.

To assess the behavior of the S schemes with meshes whose cell stencils principal directions present certain rotation with respect to those of the mesh, a new mesh configuration called Mesh 4 will be studied. This mesh is produced with the following transformation

x = ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ x + K θ L 0 sin 4π x L + Kh rand(-1, 1) x ∈ (0, L) 2 x + K θ L 0 sin 4π x L x / ∈ (0, L) 2 , ( 93
)
where K = K θ = 1 10 for this mesh configuration. Fig. 16a shows this mesh configuration for h = L/20. Fig. 16b shows the error norm for the Mesh 4 configuration and the upwind central schemes respectively. It can be observed that the US2 scheme errors do not present a smooth behavior as the mesh is refined. In our observations, the increase of the error values is related to mesh regions with maximum curvature. This implies that the S schemes stability is sensitive to mesh curvature. Our own studies have shown that the A matrix of the S2 scheme becomes singular in mesh regions with a relative orientation to the xy principal axis close to 45 • . Moreover, in this mesh configuration, the US4 schemes (results not shown in the figure) are only stable for the coarsest mesh, showing that the sensitivity of S schemes to mesh curvature increases with the polynomial degree of these schemes. At last, in this mesh configuration, the LSQ schemes present similar behavior to that reported in the previous analysis of Mesh 1, Mesh 2 and Mesh 3 configurations.

The performance of the S schemes could be improved in the Mesh 4 configuration by using a local reference frame for each stencil. However, due to the random perturbation added to each node of the mesh, the computation of a local reference frame is not straightforward. In our studies, we have tried several different reference frames described in [START_REF] Fosso Pouangué | Curvilinear finite-volume schemes using high-order compact interpolation[END_REF] and the linear mapping proposed in Eq. [START_REF] Wang | Compact high order finite volume method on unstructured grids i: basic formulations and one-dimensional schemes[END_REF]. These coordinate transformations tended to increase numerical dissipation and dispersion and could not always yield stable solutions in the proposed structured meshes. Moreover, several non-linear mappings were also tested. With these transformations the schemes presented better stability properties although their accuracy was degraded due to undetermined reasons. Further studies are required to generate local reference frames for S schemes to ensure accuracy and stability in arbitrary structured meshes.

Isentropic Euler-vortex in two-dimensional mesh

The two-dimensional static isentropic Euler vortex problem is commonly used to test the order of accuracy of numerical methods for conservation laws with non-linear fluxes [START_REF] Spiegel | A survey of the isentropic Euler vortex problem using high-order methods[END_REF], in particular for the Euler equations. The analytical solution of this problem in primitive variables is defined as follows:

ρ(x, y, t) = 1 - γ -1 2γ 2 1 γ -1 , ( 94
) u(x, y, t) = R -y x , ( 95
) P (x, y, t) = 1 - γ -1 2γ 2 γ γ -1 , ( 96
)
where R is the radius of the vortex, and is defined as

= 0 e -1 2σ 2 x 2 +y 2 R 2 . ( 97
)
In the previous expression 0 and σ are constants which are set to 1. We define the characteristic time as t c = R 2 0 . The initial conditions of the conservative variables of the Euler equations are computed from the previously described primitive variables. The computational domain is set to x ∈ [-L, L] 2 . Periodic boundary conditions are applied at each edge of the domain. It is worth noting that the choice of periodic boundary conditions can have an important influence on the order of accuracy analysis due to the non-periodic nature of the analytical and initial solutions. To avoid spurious numerical errors due to periodic boundary conditions, [START_REF] Spiegel | A survey of the isentropic Euler vortex problem using high-order methods[END_REF] stated the length of the domain must be extended enough to ensure that the amplitude of the analytical solution the edges of the domain is below or close to machine precision which is reached with L/R = 10 for double precision operations.

To reduce temporal discretization errors, the time step of the simulations is chosen as

τ t c = 1 25 h R . ( 98
)
To analyze the accuracy of the different schemes, the L ∞ error norm of the density as a function of the number of cells per radius of the vortex R/h is considered at physical time T = t c for the different schemes and Mesh 1, 2 and 3 configurations. Fig. 17a, 18a shows the results of the upwind schemes in Mesh 1 and Mesh 2 configuration respectively. These results are similar to those obtained in the one and two-dimensional linear advection analysis. The ULSQ3W outperforms all LSQ schemes for the studied range of R/h and the S schemes show lower errors than those of their LSQ counterparts. All schemes show the predicted order of accuracy. Similar conclusions can be obtained in Mesh 3 configuration as shown in Fig. 18b, although in this mesh configuration the S schemes cannot be used due to the unstructured properties of the mesh.

Fig. 17b shows the error values of the central schemes in the static Euler Vortex test case. As for upwind schemes, there results are similar to those showed in the one and two-dimensional linear advection analysis. CLSQ3W scheme shows less error for under-resolved simulations while CLSQ4 presents the lowest error values of LSQ schemes for well-resolved meshes. S schemes outperform all their LSQ counterparts. All central schemes built with an even degree polynomial show k+1exactness. An important conclusion that can be obtained from these results is the fact that the analytical findings regarding the order of accuracy of finite-volume schemes in the linear-advection equation for uniform meshes can also be applied for this non-linear problem when using the analyzed schemes (in particular the high-order face average fluxes reconstruction). For even degree polynomial reconstruction, this conclusion cannot be straightforwardly drawn from our analytical findings, as for non-linear problems it is not generally possible to obtain an analytical expression of the truncation error.

Navier-Stokes decaying compressible Homogeneous Isotropic Turbulence (HIT)

To assess the performance of the studied numerical schemes in non-linear problems in which there exists an interaction between different wavenumbers or Fourier modes, the analysis of a test case of decaying Compressible Homogeneous Turbulence in the Navier-Stokes equations is proposed [START_REF] Lee | Eddy shocklets in decaying compressible turbulence[END_REF][START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF][START_REF] Kotov | Numerical dissipation control in high order shock-capturing schemes for LES of low speed flows[END_REF][START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF][START_REF] Fernández-Fidalgo | A reduced-dissipation WENO scheme with automatic dissipation adjustment[END_REF].

The simulation is performed in a three-dimensional box such that x ∈ [0, 2π ] 3 with periodic boundary conditions on all the boundaries of the domain containing 64 3 cells. The mesh used to discretize the domain is structured and uniform. The density and pressure fields are initialized with constant values ρ 0 and P 0 while the initial velocity field, denoted as u 0 , corresponds to a random solenoidal vector field whose energy spectrum is given by

E(κ) = 16u 2 RMS,0 2 π κ 4 κ 5 0 e -2 κ 2 κ 2 0 , ( 99
)
where κ 0 = 4 and u RMS,0 is the root mean squared value of the initial velocity field. The value of the latter variable is related to the initial energy spectrum as

3 2 u 2 RMS,0 = u 0 • u 0 2 = ∞ κ=0 E(κ) dκ , (100)
where symbolizes an ensemble average over all the domain

= 1 | | d . ( 101
)
The root mean squared value of the initial velocity field can be computed from the following simulation parameters from [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF] and [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF] DNS data digitized from [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF].

Re λ 0 = ρ 0 u RMS,0 λ 0 μ = 100 , t,0 = √ u 0 • u 0 γ P 0 ρ 0 = 0.6 , ( 102
)
and λ 0 refers to the initial length of the Taylor's scales in this simulation. Parameters Re λ 0 and M t,0 are the initial Reynolds number based on the Taylor scales and the initial turbulent Mach number respectively. The viscosity μ is considered as a constant following [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF]. Please note that the choice of energy spectrum implies λ 0 = 2 κ 0 . The initial velocity field is generated following the procedure described in [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF] using a tool provided by the suite PeleC [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF]. As the analysis of this section is not focused on the asymptotic behavior of the numerical error, no distinction between cell averages and point values will be made in the initialization process.

The temporal evolution of the initial prescribed velocity spectrum generates an energy cascade which produces energy exchanges between larger and the small dissipative structures. This implies that the numerical discetization errors of underresolved waves have an impact on the global numerical errors of the simulations. This test case differs drastically from the non-linear isentropic Euler Vortex as in the latter there exists an analytical solution which shows that there is no interaction between different wavenumbers of the Fourier modes of the solution. Nonetheless, there is no guarantee that the numerical solution of the Euler Vortex doesn't present interaction between different wavenumbers, although the error analysis performed shows that this interaction doesn't have an important influence in the solution as the mesh resolution increases, at least for the configuration tested in this work. Therefore, this study should allow to determine the influence of the numerical errors due to under-resolved waves in the prediction of turbulence flows.

At last, as a result of the chosen turbulent Mach number, the simulations performed using central schemes tend to be unstable due to the appearance of eddy shocklets [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF][START_REF] Fernández-Fidalgo | A reduced-dissipation WENO scheme with automatic dissipation adjustment[END_REF]. Therefore, to stabilize the simulations using central schemes, the shock sensor described [START_REF] Larsson | Effect of numerical dissipation on the predicted spectra for compressible turbulence[END_REF] is used to modify the value of f g at each quadrature point g such that

f g = ⎧ ⎨ ⎩ 1 if -∇•u g ∇•u g + ∇×u g ≥ 0.65 0 if -∇•u g ∇•u g + ∇×u g < 0.65 ⎫ ⎬ ⎭ . ( 103
)
It is worth mentioning, that central schemes could be stabilized using other limiters, shock sensors, schemes, etc. Interested reader is referred to the works of [START_REF] Johnsen | Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves[END_REF][START_REF] Fernández-Fidalgo | A reduced-dissipation WENO scheme with automatic dissipation adjustment[END_REF] for more information on this topic. Fig. 19a shows the temporal evolution of the enstrophy ensemble average with upwind schemes computed as

ϒ = ∇ × u 2 . ( 104
)
We note that the post-processing of the results of each scheme is performed using the velocity gradient of the studied scheme. Due to the different precisions of schemes regarding the computation of this gradient, the initial values of variables which depend on the velocity gradient (such as the enstrophy) differ between schemes. The non-dimensionalization factor could be taken as the initial numerical value of the variable considered. However, this option was not chosen in this work because it would mask the disparity in accuracy of the gradient reconstruction between the different schemes. For comparison purposes, reference DNS data from [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF] obtained with very-high-order finite-different method code SMC [START_REF] Emmett | High-order algorithms for compressible reacting flow with complex chemistry[END_REF] and 512 3 points is added to the comparison. The characteristic time t c = λ 0

The results show that all schemes under-predict the production of enstrophy by a large margin in the studied mesh. If the schemes are sorted according to their error with respect to the reference DNS data, it is possible to observe that their ordering closely resembles the one obtained from the errors shown in figure 17a for under-resolved simulations. It can be observed that the results from ULSQ4 are far less accurate than those of the US4 scheme although both have the same order of accuracy and similar stencil sizes. These conclusions are also obtained when comparing ULSQ2 and USLQ3 with US2 schemes. Moreover, the use of least-squares weighting in the ULSQ3W scheme significantly improves the accuracy compared to its counterpart ULSQ3 scheme, yielding the closest results to the US4 scheme while using a substantially smaller stencil. These results seem to indicate that this simulation is drastically under-resolved although it is representative of real ILES turbulent flow test cases in which the resolution of all scales is not computationally possible. One conclusion to be obtained from these results is that the asymptotic behavior of the studied schemes cannot be directly correlated with the expected accuracy of the schemes in under-resolved turbulent flow cases.

Fig. 19b shows the temporal evolution of the enstrophy ensemble average with central schemes using the shock sensor described in Eq. ( 103). All central LSQ schemes except the CLSQ3W scheme showed unstable results, either by divergence of the simulation or by over-predicting the total kinetic energy of the flow (this variable is not shown in this section for the sake of brevity). To ensure that these unstable results are not related to the use of a shock sensor, the same simulations were performed without shock sensor and using a lower Mach number M = 0.2 and the unstable CLSQ schemes (this results from these simulations are not shown for the sake of brevity). In the latter simulations there was no divergence of the solution, although all the LSQ schemes showed an important over-prediction of the temporal evolution of the total kinetic energy. This implies that the behavior of the default CLSQ schemes tends to be unstable in non-linear under-resolved turbulent flow cases. The addition of LSQ weighting distance methods as in the CLSQ3W scheme drastically improves the numerical stability compared to the default CLSQ3 scheme. At last, it is important to note the improvement in the prediction of the enstrophy when using central schemes compared to the results of upwind schemes. Nonetheless, these improvements are partially limited due to the use of shock sensors. One should expect even more differences between central and upwind schemes for lower Mach number cases.

Conclusions

The numerical properties of the finite volume least-squares spatial discretization schemes have been studied and compared to finite-volume schemes especially catered for structured meshes. The k-exactness properties of all schemes in uniform meshes have been investigated through one-dimensional and two-dimensional analysis of the truncation error. A novel method has been developed to assess the k-exact properties of the schemes in all kind of meshes when dealing with the discretization of conservative laws through the finite volume method. This method shows that the order of accuracy of some of these schemes can be improved by using central reconstructions, although the use of central schemes in nonuniform meshes yields unstable simulations. Moreover, an extensive analysis of their dissipation and dispersion properties in uniform and unstructured meshes has been performed, extending the von Neumann's formulation to certain non-uniform and unstructured meshes. This analysis showed that, in uniform one-dimensional and two-dimensional meshes, the leastsquares schemes have worse dissipation and dispersion properties than other spatial discretization schemes developed for structured uniform meshes, specially in the under-resolved part of the spectrum. Moreover, the dissipation and dispersion study in two-dimensional non-uniform meshes proved that the dissipation and dispersion properties of the least-squares schemes are not heavily influenced by the non-uniformity of the studied meshes of this work. Furthermore, we have theoretically and numerically demonstrated that the addition of weighting distance to the least-squares schemes can improve the dissipation and dispersion properties of these schemes. To validate the analytical findings, several one-dimensional and two-dimensional simulations have been performed with linear and non-linear equations and with uniform and non-uniform meshes proving that the analytical findings can also be observed in the simulations. Special emphasis is placed on the behavior of these LSQ schemes in under-resolved simulations of turbulent flows. In these cases, it has been observed that the default LSQ schemes drastically under-perform compared to schemes for structured meshes, although the addition of weighting distance methods to the LSQ formulation notably improves the results. Moreover, the analysis showed that the numerical performance of the schemes follows closely that obtained in less complex non-linear simulations. This demonstrates that the asymptotic analysis of the numerical schemes is not always representative of the behavior of the numerical schemes in coarse turbulent test cases. To conclude, it has been shown that least-squares finite volume schemes allow to easily obtain finite volume k-exact schemes for arbitrary quadrangular meshes in linear and non-linear cases. However, this comes at a cost, the dissipation and dispersion properties of these schemes are not optimized for good-quality meshes and/or under-resolved simulations, which is critical for ILES simulations in which the important flow regions are usually discretized using near-to-uniform meshes and in which not all scales can be correctly captured due to the prohibitively associated computational cost. The addition of weighting distance methods to the LSQ formulation has been shown to greatly improve the numerical performance of the LSQ schemes. Further studies are been carried out to improve the numerical properties of least-squares schemes in near-to uniform meshes. 103). Results obtained with 256 3 cells. DNS data digitized from [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF]. WENO5 (5 th order) data was obtained using 256 3 and it was digitized from [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF].

Fig. 1 .

 1 Fig. 1. Sketch of a two-dimensional mesh with quadrilateral partitioning and the different sets used to identify faces of a given cell and quadrature points of a given face.

Fig. 2 .

 2 Fig. 2. Two-dimensional face neighbors stencils of the LSQ1-2-3 schemes.

Fig. 3 .

 3 Fig. 3. Two-dimensional vertex neighbors stencils of the S2 and S4 schemes.

Fig. 4 .

 4 Fig. 4. Grids used to test the k-exact properties of the LSQ schemes.

Fig. 5 .

 5 Fig. 5. Dispersion errors of 1D upwind schemes as a function of λ/h and for CFL ∈ [0.1, 0.5, 1.0]. Dashed lines are related to the function Im(D) ∼ (λ/h) n

Fig. 6 .

 6 Fig. 6. Dispersion errors of 1D central schemes as a function of λ/h and for CFL ∈ [0.1, 0.5, 1.0]. Dashed lines are related to the function (D) ∼ (λ/h) n

Fig. 7 .

 7 Fig. 7. Dissipation errors of 1D upwind schemes as a function of λ/h and for CFL ∈ [0.1, 0.5, 1.0]. Dashed lines are related to the function Im(D) ∼ (λ/h) n

Fig. 8 .

 8 Fig. 8. Diffusion (a) and Dispersion (b) errors of the two-dimensional US2, US4, ULSQ2 and ULSQ4 schemes at ζ = 45 • . The radial coordinate of the figures represents the cells per wavelength λ/h parameter, while the azimuthal coordinate is the wave angle. Results are obtained for CFL = 0.05. (For interpretation of the colors in the figure(s), the reader is referred to the web version of this article.)

Fig. 9 .

 9 Fig. 9. Dissipation and dispersion errors of 2D upwind schemes of waves with θ = 45 • as a function of λ/h and for cτ /h = 0.05. Dashed lines are related to the function Im(D) ∼ (λ/h) n where n is an exponent ranging from 2 to 6 and are used to assess the leading term order of the dissipation errors.

Fig. 10 .

 10 Fig. 10. Dissipation and dispersion errors of 2D upwind schemes of waves with θ = 45 • in Mesh 2 configuration as a function of λ/h and for cτ /h = 0.05. Dashed lines are related to the function Im(D) ∼ (λ/h) n where n is an exponent ranging from 2 to 6 and are used to assess the leading term order of the dissipation and dispersion errors.

Fig. 10 and

 10 Fig.11show the dissipation and dispersion errors of the upwind least-squares schemes for the Mesh 2 and Mesh 3 configurations at n = 0 and cτ /h = 0.05. Both figures show similar results. The dissipation and dispersion properties of the least-squares schemes in the non-uniform meshes of this study are very similar to those found in uniform meshes.

Fig. 11 .

 11 Fig. 11. Dissipation and dispersion errors of 2D upwind schemes of waves with θ = 45 • in Mesh 3 configuration as a function of λ/h and for cτ /h = 0.05. Dashed lines are related to the function Im(D) ∼ (λ/h) n where n is an exponent ranging from 2 to 6 and are used to assess the leading term order of the dissipation and dispersion errors.

Fig. 12 .

 12 Fig. 12. L ∞ error norm of the solution in the 1D test case at t = T using upwind schemes for three different values of the CFL and as a function of λ/h.

Fig. 13 .

 13 Fig. 13. L ∞ error norm of the solution in test case at t = T using central schemes for three different values of the CFL and as a function of λ/h.

Fig. 14 .

 14 Fig. 14. L ∞ error norm of the solution in the two-dimensional linear advection test case with uniform meshes at t = T using upwind (a) and central (b) schemes with cτ /h = 0.05 and as a function of λ/h.

Fig. 15 .

 15 Fig. 15. L ∞ error norm of the solution in the two-dimensional linear advection test case with Mesh 2 (a) and 3 (b) configurations at t = T using upwind schemes, cτ /h = 0.05 and as a function of λ/h.

Fig. 16 .

 16 Fig. 16. L ∞ error norm of the solution in the two-dimensional linear advection test case with Mesh 4 configuration (a) at t = T using upwind schemes, cτ /h = 0.05 and as a function of λ/h.

Fig. 17 .

 17 Fig. 17. L ∞ error norm of the density in the Euler Vortex test case with uniform meshes (Mesh 1) at T = t c for upwind (a) and central (b) schemes and as a function of R/h.

Fig. 18 .

 18 Fig. 18. L ∞ error norm of the density in the Euler Vortex test case with Mesh 2 (a) and Mesh 3 (b) configuration at T = t c for different upwind schemes and as a function of R/h. Note that the S schemes are missing as they are not compatible with unstructured meshes.

Fig. 19 .

 19 Fig. 19. Non-dimensional enstrophy temporal evolution in the decaying compressible HIT test case with M t,0 = 0.6 for upwind (a) and central (b) schemes.

Fig. A. 20 .

 20 Fig. A.20.Non-dimensional enstrophy temporal evolution in the decaying compressible HIT test case with M t,0 = 0.6 using S central schemes together with the shock sensor described in Eq. (103). Results obtained with 256 3 cells. DNS data digitized from[START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF]. WENO5 (5 th order) data was obtained using 2563 and it was digitized from[START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF].

Table 1

 1 Number of neighbors M of the different LSQ schemes in uniform and structured meshes as a function of the number of dimensions of the problem.

		LSQ1	LSQ2	LSQ3	LSQ4
	1 D	3	5	7	9
	2D	5	13	25	41
	3D	7	25	63	129

Table 2 β

 2

		LSQ1	LSQ2	LSQ3	LSQ4
	1D	2	2	2	2
	2D	2	12/5	8/3	20/7
	3D	2	8/3	62/19	128/34

parameter of the different LSQ schemes in uniform and structured meshes as a function of the number dimensions of the problem.

Table 3

 3 Number of neighbors M of the different Q schemes in uniform and structured meshes as a function of the number of dimensions of the problem.

		S2	S4
	1D	3	5
	2D	9	27
	3D	25	125

Table 4

 4 Number of quadrature points used with LSQ schemes as a function of the dimensions of the problem N d and the polynomial order p.

	p

Table 5

 5 Leading truncation error term and order of accuracy of the different one-dimensional least-squares schemes.

		E	Order of accuracy
	ULSQ1	-0.083ch 2 ∂ 3 u

Table 6

 6 Leading truncation error term and order of accuracy of different schemes described in Section 3.2.

		E			Order of accuracy
	US2	-ch 3 12	∂ 4 u ∂x 4 i	+ O h 4	3 rd
	CS2	-ch 4 30	∂ 5 u ∂x 5 i	+ O h 6	4 th
	US4	-ch 5 60	∂ 6 u ∂x 6 i	+ O h 6	5 th
	CS4	-ch 6 140	∂ 7 u ∂x 7 i	+ O h 8	6 th

Table 7

 7 Maximum face average and quadrature points error in the uniform mesh case (Mesh 1) with h = L/20. Symbol -denotes errors considered as machine round-off errors (< 10 -14 ).

		Maximum Face Average Error	Maximum Quadrature Points Error
		p = p	p = p + 1	p = p + 2	p = p	p = p + 1	p = p + 2
	ULSQ1	-	5.572E-5	1.717E-4	-	5.205E-5	1.607E-4
	CLSQ1	-	5.373E-5	1.612E-4	-	5.373E-5	1.655E-4
	ULSQ2	-	8.640E-6	3.480E-5	-	9.164E-6	3.696E-5
	CLSQ2	-	-	2.276E-7	-	8.183E-6	3.296E-5
	ULSQ3	-	1.379E-7	7.198E-7	-	2.017E-7	1.046E-6
	CLSQ3	-	1.282E-7	6.411E-7	-	1.282E-7	6.690E-7
	ULSQ3W	-	5.042E-8	2.613E-7	-	9.158E-8	4.697E-7
	CLSQ3W	-	4.689E-8	2.345E-7	-	4.689E-8	2.429E-7
	ULSQ4	-	7.070E-8	4.274E-7	-	8.034E-8	4.866E-7
	CLSQ4	-	-	2.884E-9	-	6.458E-8	3.904E-7
	US2	-	1.968E-7	7.894E-7	-	2.036E-7	8.166E-7
	CS2	-	-	2.241E-9	-	1.925E-7	7.722E-7
	US4	-	4.227E-11	2.544E-10	-	4.272E-11	2.571E-10
	CS4	-	-	7.606E-13	-	4.075E-11	2.453E-10

Table 8

 8 Maximum face average and quadrature points error in the Mesh 2 configuration with

		Maximum Face Average Error	Maximum Quadrature Points Error
		p = p	p = p + 1	p = p	p = p + 1
	ULSQ1	-	1.502E-4	-	1.461E-4
	CLSQ1	-	8.090E-5	-	1.273E-4
	ULSQ2	-	9.648E-6	-	1.019E-5
	CLSQ2	-	2.381E-6	-	1.004E-5
	ULSQ3	-	3.772E-7	-	4.517E-7
	CLSQ3	-	2.051E-7	-	4.047E-7
	ULSQ3W	-	1.116E-7	-	1.615E-7
	CLSQ3W	-	8.781E-8	-	1.026E-7
	ULSQ4	-	8.643E-8	-	9.424E-8
	CLSQ4	-	1.650E-8	-	6.978E-8
	US2	-	6.574E-8	-	6.653E-8
	CS2	-	1.510E-8	-	6.183E-8
	US4	-	5.151E-12	-	5.251E-12
	CS4	-	1.127E-12	-	5.151E-12

h ∼ L/20. Symbol -denotes errors considered as machine round-off errors (< 10 -14 ).

Table 9

 9 Maximum face average and quadrature points error in the Mesh 3 configuration with h ∼ L/20. Symbol -denotes errors considered as machine round-off errors (< 10 -14 ).

		Maximum Face Average Error	Maximum Quadrature Points Error
		p = p	p = p + 1	p = p	p = p + 1
	ULSQ1	-	1.584E-4	-	1.471E-4
	CLSQ1	-	7.982E-5	-	1.584E-4
	ULSQ2	-	1.139E-5	-	1.180E-5
	CLSQ2	-	1.684E-6	-	9.351E-6
	ULSQ3	-	4.684E-7	-	6.311E-7
	CLSQ3	-	1.936E-7	-	4.684E-7
	ULSQ3W	-	9.529E-8	-	1.536E-7
	CLSQ3W	-	7.387E-8	-	7.984E-8
	ULSQ4	-	6.802E-8	-	8.250E-8
	CLSQ4	-	1.662E-8	-	6.450E-8

Table 10

 10 Maximum CFL value to ensure temporal stability in the one-dimensional linearadvection as a function of the different schemes studied in this work. Results obtained with the RK33 time integration scheme.

	Scheme	CFL MAX
	ULSQ1	1.176
	CLSQ1	1.101
	ULSQ2	1.353	Scheme	CFL MAX
	CLSQ2	1.666	US2	1.626
	ULSQ3	1.312	CS2	1.262
	CLSQ3	1.342	US4	1.435
	ULSQ3W	1.434	CS4	1.092
	CLSQ3W	1.111	(b) One-dimensional S schemes
	ULSQ4	1.292
	CLSQ4	1.718
	(a) One-dimensional LSQ schemes

u RMS,0 is used to normalize the temporal evolution.

Acknowledgements

This work was performed using HPC resources from GENCI-IDRIS and GENCI-CINES on Jean Zay, Occigen (Grant A0082A07178) and CALMIP on Olympe (Grant 2020-p1425).

CRediT authorship contribution statement G. Saez-Mischlich: Conceptualization, Methodology, Software, Validation, Writing -original draft. J. Sierra-Ausin: Conceptualization, Methodology, Writing -original draft. G. Grondin: Supervision, Writing -review & editing. J. Gressier: Supervision, Writing -review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Appendix A. HIT initial condition validation

To ensure that the initial condition to simulate the HIT test case yields results which can be compared with DNS reference data of [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF], Fig. A.20 shows the ensemble enstrophy as a function of the time for the central S schemes in a 256 3 mesh. The results are compared with DNS and FV-WENO5 (5 th order scheme) data that obtained from [START_REF] Motheau | Investigation of finite-volume methods to capture shocks and turbulence spectra in compressible flows[END_REF]. FV-WENO5 was obtained using the same mesh discretization (256 3 cells). This FV-WENO5 scheme is a finite-volume scheme which uses the S4 scheme formulation described in this work (Section 3.2) in combination with a WENO method to better treat solution discontinuities. The results show almost good agreement with DNS data and that the CS4 scheme with the Ducros shock sensor yields more accurate results than those of the FV-WENO5 scheme. This figure also serves as a validation for both the initial condition of the HIT test case and the implementation of the high-order viscous fluxes discretization described in Section 3.4 and Section 2.3.