Kanokwarun Para

Bubpha Jitsom

Robert Eymard

Surattana Sungnul

Sekson Sirisubtawee

Sutthisak Phongthanapanich

An accuracy comparison of piecewise linear reconstruction techniques for the characteristic finite volume method for two-dimensional convection-diffusion equations

Keywords: Kanokwarun Para, Bubpha Jitsom, Robert Eymard, Surattana Sungnul, Sekson Sirisubtawee, Sutthisak characteristic finite volume method, convection-diffusion equation, piecewise linear reconstruction, two-dimensional triangular grid

published or not. The documents may come   L'archive ouverte pluridisciplinaire

An accuracy comparison of piecewise linear reconstruction techniques for the characteristic finite volume method for two-dimensional convection-diffusion equations 1 Introduction

Research on numerical methods for solving convection-diffusion equations is of great importance since these equations are a basis for describing physical phenomena where particles, energy, or other physical quantities are transferred inside a physical system due to diffusion and convection processes. These equations are also used to describe advection-diffusion processes in environmental science, e.g., pollutant transport in the atmosphere, oceans, rivers, or groundwater [START_REF] Slepchenko | Numerical approach to fast reactions in reaction-diffusion systems: application to buffered calcium waves in bistable models[END_REF][START_REF] Gallo | Finite volume/mixed finite element analysis of pollutant transport and bioremediation in heterogeneous saturated aquifers[END_REF][START_REF] Bertolazzi | A finite volume method for transport of contaminants in porous media[END_REF]. Convection is usually defined as the movement of a fluid due to heat, advection is usually defined as the movement of a substance due to the movement of a fluid transport medium, and diffusion refers to the dispersion of some substance through the physical domain of some other substance, e.g., the dispersion of a pollutant through a river. Generally, the scale of the diffusion is smaller by several orders of magnitude compared to the size of the advective flow field. The accurate simulation of such processes requires numerical methods that are able to compute if sharp layers exist without the occurrence of spurious oscillations.

Many methods have been proposed to overcome the instability and inaccuracy of numerical solutions of convection-diffusion equations. Methods which are usually used to solve convectiondiffusion equations include the finite difference method, the finite element method, and the finite volume method. The principle of the finite difference method is to select some discretization points, then to assign one discrete unknown and write one equation per discretization point [START_REF] Morton | Numerical Solution of Partial Differential Equations[END_REF]. Also, the derivatives of the unknown are replaced by finite differences using a Taylor-series expansion. The finite difference method becomes difficult to use when the coefficients involved in the equation are discontinuous. For the finite element method, the domain of the convectiondiffusion equation is discretized into a set of elements and the method of weighted residuals is usually applied to obtain a set of equations for each element, which are then assembled to form the system of equations. The solution at the element nodes are then obtained by solving the system of equations [START_REF] Nithiarasu | The characteristic-based split (cbs) scheme-a unified approach to fluid dynamics[END_REF][START_REF] Zienkiewicz | A general algorithm for compressible and incompressible flow-part i. the split, characteristic-based scheme[END_REF]. The accuracy and the stability of this method needs to be further improved to achieve an accurate transient convection-diffusion solution [START_REF] Donea | High-order accurate time-stepping schemes for convection-diffusion problems[END_REF]. For the finite volume method, the domain is discretized into a set of regions called control volumes. The divergence theorem is then applied to derive the system of equations. Higher-order temporal and spatial discretizations can also be used to obtain a more accurate solution. In the finite-volume method, discontinuities of coefficients are not a problem if the mesh is chosen so that the discontinuities of the coefficients occur on the boundaries of the control volumes [START_REF] Eymard | Finite volume methods[END_REF]. In recent years, many researchers have used finite volume methods. Some examples of this research are as follows. In 2003, Aboubacar and Webster [START_REF] Aboubacar | Development of an optimal hybrid finite volume/element method for viscoelastic flows[END_REF] proposed a cell-vertex hybrid finite volume-element method that is implemented on triangles and applied to the numerical solution of Oldroyd model fluids in contraction flows. They found that a linear finite volume stress representation with discontinuous stress gradients, and incorporating locally reduced quadrature at the re-entrant corner could enhance the stability properties of the numerical scheme. In 2006, Pudykiewicz [START_REF] Pudykiewicz | Numerical solution of the reaction-advection-diffusion equation on the sphere[END_REF] derived a finite volume algorithm for obtaining the solution of an advectiondiffusion equation on a sphere based on the principle of semi-discretization. The monotonicity of the scheme was achieved with explicit adaptive dissipation. In 2008, Gao and Yuan [START_REF] Gao | The characteristic finite volume element method for the nonlinear convection-dominated diffusion problem[END_REF] presented a characteristic finite volume element method for multiple space, nonlinear, convectiondominated diffusion problems. The calculus of variations, the commutating operator, and the theory of prior estimates was adopted in their proposed scheme. In 2011, Ten Thije Boonkkamp and Anthonissen [START_REF] Ten Thije | The finite volume-complete flux scheme for advection-diffusion-reaction equations[END_REF] proposed a finite volume scheme for advection-diffusion-reaction equations. Their scheme included a new integral representation for the flux of the one-dimensional advection-diffusion-reaction equation which was derived from the solution of a local boundary value problem for the entire equation including the source term. In 2012, Zhang et al. [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF] presented a fourth-order finite volume method for solving the advection-diffusion equation with an AMR-enabled multigrid solver. The fourth-order Runge-Kutta method was used for the time integration. The fourth-order accuracy was used to approximate face-averaged values from cellaveraged values. In 2013, Filimon et al. [START_REF] Filimon | High-order finite volume schemes based on defect corrections[END_REF] proposed an iterated defect correction approach to achieve higher-order accuracy for the approximation of steady-state solutions. The scheme applied a posteriors to estimate the local discretization error of the lower-order finite volume scheme. This estimate was then used to iteratively shift the solution to higher-order accuracy by polynomial reconstruction. In 2013, Chen et al. [START_REF] Chen | A two-grid characteristic finite volume element method for semilinear advection-dominated diffusion equations[END_REF] presented the concept of a two-grid finite volume element method combined with a modified method of characteristics for analyzing semi-linear, time-dependent, advection-dominated diffusion equations in two space dimensions.

The solution of a nonlinear system on the fine-grid space was reduced to the solution of nonlinear and linear systems on the coarse-grid and the fine-grid, respectively. In 2014, Arachchige and Pettet [START_REF] Arachchige | A finite volume method with linearisation in time for the solution of advection-reaction-diffusion systems[END_REF] introduced a new method for numerically solving the advection-reaction-diffusion systems in one space dimension based on the usual finite volume method with a third order unwinding scheme for discretizing the advection term in space. A new linearization technique for the temporal integration of the nonlinear source terms was proposed to provide a very effective method in terms of computational cost for the simulations. In 2015, Liu et al. [START_REF] Liu | Coupled nonlinear advection-diffusion-reaction system for prevention of groundwater contamination by modified upwind finite volume element method[END_REF] applied a two-grid algorithm to decompose a nonlinear advection-dominated diffusion-reaction equation into a small nonlinear system on a coarse grid and a linear system on a fine grid, as previously proposed by Chen et al. [START_REF] Chen | A two-grid characteristic finite volume element method for semilinear advection-dominated diffusion equations[END_REF]. The modified upwind finite volume element method was then used to approximate the solution of the systems. In 2016, Tambue [START_REF] Tambue | An exponential integrator for finite volume discretization of a reactionadvection-diffusion equation[END_REF] applied a finite volume method with a two-point flux approximation on regular meshes to solve an advection-diffusion-reaction problem with a nonlinear reaction term. Exponential time differencing of order one was then applied for temporal discretization. In 2017, Lin and Zhang [START_REF] Lin | High-order finite-volume solutions of the steady-state advectiondiffusion equation with nonlinear robin boundary conditions[END_REF] proposed a high-order finitevolume scheme for solving a steady-state advection-diffusion equation with nonlinear Robin boundary conditions. A generic algorithm for generating third-order, fourth-order, and even higher-order explicit ghost-filling formulas was used to enforce the nonlinear Robin boundary conditions in multiple dimensions under the framework of finite volume methods. In 2017, Yang and Tine [START_REF] Yang | A hybrid finite volume method for advection equations and its applications in population dynamics[END_REF] proposed a hybrid finite volume scheme where the construction of the numerical flux appeared as a combination of the WENO and Anti Dissipative Method fluxes for transport type equations. This property of the numerical flux is very suitable for deriving the long term asymptotic behavior of the solution of population dynamics models as previously shown in numerical simulations of polymerization/depolymerization type models. In 2018, Xu [START_REF] Xu | A modified finite volume method for convection-diffusion-reaction problems[END_REF] proposed a modified finite volume method to solve convection-diffusion-reaction problems. The scheme was much more accurate and stable than the conventional finite volume method without widening the computational stencil. Some undetermined coefficients were introduced in discretizing the equation and then determined analytically by making the truncation error of the numerical scheme vanish and by expressing the higher order derivatives of the unknown function with lower order ones through the governing equation. Recently, Lan et al. [START_REF] Lan | A new positive finite volume scheme for two-dimensional convection-diffusion equation[END_REF] proposed a new positive finite volume scheme for two-dimensional convection-diffusion equations on deformed meshes. The convective flux was approximated using available information of the diffusive flux to keep the upwind property and obtain second-order accuracy. This method did not need to introduce any slope limiting technique. The outline of this paper is as follows. In section 2, we summarize an explicit characteristicbased finite volume element method developed by Phongthanapanich and others that can provide stabilized numerical solutions for pure advection and advection-dominated diffusion prob-lems [START_REF] Phongthanapanich | A characteristic-based finite volume element method for convection-diffusion-reaction equation[END_REF][START_REF] Gao | The characteristic finite volume element method for the nonlinear convection-dominated diffusion problem[END_REF][START_REF] Phongthanapanich | Combined finite volume and finite element method for convection-diffusion-reaction equation[END_REF][START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF][START_REF] Phongthanapanich | An explicit characteristic finite volume element method for non-divergence free convection-diffusion-reaction equation[END_REF]. The idea of this method is to derive a convection-diffusion equation along a characteristic path. The explicit cell-centered finite volume method is then employed to discretize this characteristic-based equation. The weighted residuals finite element method is then used to approximate cell-face gradient quantities. In Section 3, we analyze the convergence of the method for a one-dimensional example. In section 4, we examine the accuracy and robustness of four piecewise linear reconstruction techniques, e.g., Frink [START_REF] Frink | Fast upwind solver for the euler equations on three-dimensional unstructured meshes[END_REF], Holmes-Connell [START_REF] Holmes | Solution of the 2d navier-stokes equations on unstructured adaptive grids[END_REF], Green-Gauss [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF], and least-squares methods [START_REF] Barth | A 3-d upwind euler solver for on unstructured meshes[END_REF] on structured triangular grids for three test cases. These cases consist of (1) mixing of hot with cold fronts [START_REF] Phongthanapanich | A characteristic-based finite volume element method for convection-diffusion-reaction equation[END_REF][START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF], (2) rotation of a Gaussian pulse [START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF], and (3) triangular wave flow [START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF]. Finally, in section 5, we give conclusions.

2 The Numerical Scheme

The characteristic unsteady convection-diffusion equation

The governing partial differential equation for the two-dimensional unsteady convection-diffusion problem is

∂φ ∂t + ∇ • (vφ -ε∇φ) = 0 in Ω × T, φ(x, 0) = φ 0 (x) x ∈ Ω, Ω ⊂ R 2 , (1)
where φ is the scalar quantity, v = v(x) is the given convection velocity vector, ε ≥ 0 is the diffusion coefficient, and t ∈ (0, T ] for T < ∞. The boundary conditions are

φ = g D in ∂Ω D , ε ∂φ ∂n = g N in ∂Ω N , (2)
with

∂Ω = ∂Ω D ∪ ∂Ω N with ∂Ω D ∪ ∂Ω N = ∅.
If a moving coordinate is assumed along the path of the characteristic wave, x with a speed of v, the change of variable x to x is expressed by dx = dx -vdt. By applying a characteristic approximation to trace advection in time as explained in Refs. [START_REF] Phongthanapanich | A characteristic-based finite volume element method for convection-diffusion-reaction equation[END_REF][START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF][START_REF] Phongthanapanich | An explicit characteristic finite volume element method for non-divergence free convection-diffusion-reaction equation[END_REF] to Equation (1), we obtain ∂φ ∂t

-∇ • ε∇ φ = 0, (3)
where all terms are to be evaluated at x =x (t). We then discretize the time with time steps ∆t and define φ n = φ(t n ), where t n = n∆t, n = 0, 1, 2, . . .. Then, by carrying out a Taylor-series expansion, the convection term reappears in the equation along with an additional secondorder term. This term acts as a smoothing operator that reduces the oscillations arising from the spatial discretization of the convection term. The fully explicit characteristic unsteady convection-diffusion equation is then given by

φ n+1 -φ n = -∆t (v • ∇φ n ) + ε∇ • ∇φ n + (∆t) 2 2 v • ∇ (v • ∇φ n ) . (4)
Together with the use of the divergence-free assumption (∇ • v = 0), Equation ( 4) can be written finally in the conservation form as

φ n+1 -φ n = -∆t [∇ • (vφ n -ε∇φ n )] + (∆t) 2 2 ∇ • [v (v • ∇φ n )] . (5)

The finite volume method

The computational domain is first discretized into a collection of non-overlapping control volumes Ω i = Ω, i = 1, ..., N, that completely cover the domain such that Ω = ∪ N i=1 Ω i , Ω i = 0 and Ω i ∩ Ω j = 0 if i = j. To obtain the finite volume equation, Equation ( 5) is then integrated over the control volume Ω i to yield

Ω i φ n+1 -φ n dA = Ω i -∆t [∇ • (vφ n -ε∇φ n )] + (∆t) 2 2 ∇ • [v (v • ∇φ n )] dA. (6)
Then the divergence theorem is applied to the spatial terms to yield a fully explicit characteristicbased scheme for solving Equation ( 5) in the form

φ n+1 i = φ n i - ∆t | Ω i | N f j=1 | Γ ij | nij • v n ij φ n ij -ε∇φ n ij - ∆t 2 v n ij (v n i • ∇φ n i ) , (7)
where N f is the number of adjacent cell faces, |Ω i | is the measure of Ω i , Γ ij is the segment of the boundary ∂Ω i between the two adjacent control volumes Ω i and Ω j and |Γ ij | is the measure of Γ ij . The quantities at the time t n are defined by

φ n i = φ i (t n ) = 1 |Ω i | Ω i φ(x, t n )dx, (8)
and

φ n ij = φ ij (t n ). (9)
Finally, the scalar quantities at the cell faces, φ n ij , are approximated by applying a Taylor-series expansion in space such that

φ n ij = φ n i + (x ij -x i ) • ∇φ n i , v n ij • nij ≥ 0 φ n j + (x ij -x j ) • ∇φ n j , otherwise , (10)
where v n ij is the scaled normal velocity at Γ ij , and | Γ c ij | is the cell characteristic length. Moreover, the gradient term, ∇φ n ij , is approximated by the weighted residuals method which is commonly used in the finite element technique [START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF], and the time-step within each control volume i for Equation ( 7) is determined from

∆t = C min i | Ω i | max j=1,...,N f | v n ij | , | Γ c ij | 2 2ε , (11)
where 0 < C ≤ 1.

The upwind schemes require flow states to be specified on the left and the right of a control surface. The first-order accurate spatial discretization gives a constant solution inside each control volume. The piecewise linear reconstruction must be applied to each control volume to achieve second-order accuracy. The main aim of this research is to compare the accuracy and robustness of four piecewise linear reconstruction techniques, e.g., Frink, Holmes-Connell, Green-Gauss, and least-squares methods for a computational domain discretized into a triangular grid.

Frink's linear-reconstruction method [START_REF] Frink | Fast upwind solver for the euler equations on three-dimensional unstructured meshes[END_REF] can explicitly evaluate the gradient for a triangular grid due to the invariant geometric properties of the grid, e.g., a line from a node through the cell-centroid will intersect the midpoint of the opposite face, and the distance from the cellcentroid to the face-midpoint is one-third of the distance from the face-midpoint to the opposite node. The scalar values at the cell-face of two adjacent grids (φ L , φ R ) are calculated by

φ L = φ i + Ψ i 3 (φ i1 + φ i2 ) 2 -φ i3 , φ R = φ j + Ψ j 3 (φ j1 + φ j2 ) 2 -φ j3 , (12)
where φ i and φ j are assumed to be the solutions at the cell-centroids, φ ik , φ jk , where k = 1, 2, 3, are the solutions at the nodes of the left and right cells, respectively, and Ψ i and Ψ j are limiter functions [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF][START_REF] Venkatakrishnan | Convergence to steady state solutions of the euler equations on unstructured grids with limiters[END_REF][START_REF] Tamamidis | A new upwind scheme on triangular meshes using the finite volume method[END_REF]. The value at a node may be computed by inverse-distance weighting or by the pseudo-Laplacian method [START_REF] Frink | Fast upwind solver for the euler equations on three-dimensional unstructured meshes[END_REF]. The inverse-distance weighting from a cell-centroid to the cell nodes that preserves the principle of positivity can be expressed as

φ n = Nc j=1 φ j |r j | / Nc j=1 1 |r j | , (13)
where φ j are the surrounding cell-centered values of the node n , |r j | is the distance from the cell-centroid to node n, and N c is the number of the surrounding cells.

Otherwise, the pseudo-Laplacian method proposed by Holmes and Connell [START_REF] Holmes | Solution of the 2d navier-stokes equations on unstructured adaptive grids[END_REF] can be used to determine nodal quantities such that

φ n = Nc j=1 (w j φ j )/ Nc j=1 φ j , (14)
where the weights w j are determined from

w j = 1 + λ x (x j -x n ) + λ y (y j -y n ), (15)
and the Lagrange multipliers, λ x and λ y are defined by

λ x = Nc j=1 (x j -x n )(y j -y n ) Nc j=1 (y j -y n ) - Nc j=1 (y j -y n ) 2 Nc j=1 (x j -x n ) Nc j=1 (x j -x n ) 2 Nc j=1 (y j -y n ) 2 -   Nc j=1 (x j -x n )(y j -y n )   2 , λ y = Nc j=1 (x j -x n )(y j -y n ) Nc j=1 (x j -x n ) - Nc j=1 (x j -x n ) 2 Nc j=1 (y j -y n ) Nc j=1 (x j -x n ) 2 Nc j=1 (y j -y n ) 2 -   Nc j=1 (x j -x n )(y j -y n )   2 . ( 16
)
The Green-Gauss approach [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] approximates the gradient of φ as the surface integral of the product of φ with an outward unit normal vector over a control volume such that

∇φ = 1 Ω ∂Ω φndS. (17)
For the center-centered scheme, the gradient at the center of cell i of the triangular grid [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF] can be obtained from

∇φ i = 1 2Ω 3 k=1 (φ i + φ k )n ik ∆S ik , (18)
where the summation extends over all faces of the cell with volume Ω. The n ik is the outward unit normal vector and ∆S ik is the k-face area of the cell i.

Lastly, the least-squares approach [START_REF] Barth | A 3-d upwind euler solver for on unstructured meshes[END_REF] is based on the use of a Taylor-series expansion for each edge which is incident to the cell-centered i. The linear reconstruction term along the common edge between control volume i and j can be computed as

∇φ i • r ij = φ j -φ i , (19)
where r ij is the vector from cell-centered i to j. By applying Equation ( 19) to three neighboring cells, we obtain the following system of linear equations

  ∆x i1 ∆y i1 ∆x i2 ∆y i2 ∆x i3 ∆y i3   ∂φ i ∂x ∂φ i ∂y =    φ 1 -φ i φ 2 -φ i φ 3 -φ i    . ( 20
)
Finally, it is worth comparing the finite volume method (FVM) and discontinuous Galerkin method (DGM) concerning their efficiency. Moreover, DGM gives error estimates in stronger norms and with less regularity [START_REF] Di Pietro | Mathematical Aspects of Discontinuous Galerkin Methods[END_REF]. They are identical to finite volume schemes in the case of the formal first order (degree k = 0), and their computational costs are the same. However, for higher-order, DGM uses nonconformal ansatz functions whose restrictions to single cells are higher-degree polynomials. DGM with k = 1 may yield similar results as a piecewise linear reconstruction technique, but the effort for computing the fluxes for DGM is higher than for FVM.

Convergence of the characteristic finite volume method for one-dimensional convection

The analysis of the convergence of the scheme defined by ( 7) presents a few difficulties, that we analyse for the case of a one-dimensional problem and a non-uniform mesh. Let Ω = R. We intend to approximate the pure convection equation with a constant convection velocity v > 0,

∂φ ∂t + ∇ • (vφ) = ∂φ ∂t + v ∂φ ∂x = 0. ( 21
)
We remark that, defining φ ini (x) = φ(x, 0) for all x ∈ R, the solution of ( 21) is given by φ(x, t) = φ ini (x -vt).

Although we don't include a diffusion term in the continuous equation, we include in the scheme a diffusion term, which is vanishing with the size of the mesh, and which plays a central role in the convergence proof. We assume that an increasing sequence (x i+ 1

2

) i∈Z ⊂ R is given such that lim i→±∞ x i+ 1 2 = ±∞, we define h i = x i+ 1 2 -x i-1 2
. We assume that there exist h > 0, measuring the size of the mesh, and β ∈ (1, +∞), measuring the regularity of the mesh, such that

βh ≥ h i ≥ 1 β h for all i ∈ Z. ( 22
)
This means that the ratio of the sizes of any pair of control volumes is bounded by below and above. We then define the mesh

(Ω i ) i∈Z by Ω i = (x i-1 2 , x i+ 1 2
). The center of gravity of Ω i is

x i = 1 2 (x i-1 2 + x i+ 1 2
), and there holds |Ω i | = h i . Then, for i ∈ Z, we let

x ij = x i± 1 2 for j = i ± 1. The initial condition is approximated by ∀i ∈ Z, φ 0 i = φ ini (x i ), (23)
and Equation ( 7) becomes, in this 1D case,

∀i ∈ Z, ∀n ∈ N, |Ω i | φ n+1 i -φ n i ∆t + j=i±1 |Γ ij |n ij v n ij φ n ij - ∆tv n i 2 ∇φ n i -ε∇φ n ij = 0. ( 24
)
with

|Γ ij | = 1, v n ij = v, v n i = v, nij = ±1 for j = i ± 1, and
∇φ n i := φ n i+1 -φ n i-1 x i+1 -x i-1 , ∇φ n ij := φ n j -φ n i x j -x i .
The vanishing diffusion term is defined by ε = vβh. For j = i + 1, using the fact that the upwinding direction is from i to i + 1, we can write

φ n ij = φ n i + (x ij -x i )∇φ n i = φ n i + h i 2 φ n i+1 -φ n i-1 x i+1 -x i-1 .
We then set

φ n i+ 1 2 := φ n ij - ∆tv 2 ∇φ n i - ε v ∇φ n ij = φ n i + α i (φ n i+1 -φ n i-1 ) + θ i+ 1 2 (φ n i -φ n i+1 ), (25)
with α i = h i -∆tv 2(x i+1 -x i-1 ) and θ i+ 1 2 = βh x i+1 -x i . In the following, we will assume that

v∆t ≤ C cfl h, (26)
for the value C cfl (which only depends on β) defined below in [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF].

Using C cfl ≤ 1 β , x i+1 -x i = 1 2 (h i + h i+1 ) and x i+1 -x i-1 = 1 2 (h i-1 + h i+1 ) + h i , we obtain 0 ≤ α i ≤ 1 2 and 1 ≤ θ i+ 1 2 ≤ β 2 . ( 27
)
From ( 25) and ( 24), we get

∀i ∈ Z, ∀n ∈ N, h i φ n+1 i -φ n i ∆t + v(φ n i+ 1 2 -φ n i-1 2 ) = 0. ( 28
)
Note that (28) makes clear the finite volume structure of Scheme [START_REF] Pudykiewicz | Numerical solution of the reaction-advection-diffusion equation on the sphere[END_REF] and that, due to the negative coefficient in the dependence of φ n i-1 2 with respect to φ n i-2 , no L ∞ bound can be proved on the scheme. We can state the following convergence result. Theorem 3.1 Let us assume that, letting φ be the solution to [START_REF] Phongthanapanich | Combined finite volume and finite element method for convection-diffusion-reaction equation[END_REF] and defining φ ini (x) = φ(x, 0) for all x ∈ R, the function φ ini ∈ C 2 (R) has a compact support. Let T > 0 and β ∈ (1, +∞) be given, and let h ∈ (0, 1) and N ∈ N be given such that ∆t = T /N satisfies (26) where C cfl is defined by [START_REF] Zhang | A fourth-order accurate finite-volume method with structured adaptive mesh refinement for solving the advection-diffusion equation[END_REF]. For a given mesh defined as above, and satisfying [START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF], let φ n i , for all i ∈ Z and n ∈ N, be the solution of Scheme (23), ( 25) and [START_REF] Ten Thije | The finite volume-complete flux scheme for advection-diffusion-reaction equations[END_REF]. Then there exists C > 0, only depending on T , v, φ ini and β such that

e N L 2 ≤ Ch, (29)
where e N L 2 , the approximate L 2 norm of the error at time T , is defined by

e N 2 L 2 = i∈Z h i (φ N i -φ(x i , T )) 2 .

Proof

Let us compute the consistency error, obtained by replacing φ n i in the scheme by φ n i := φ(x i , n∆t). We obtain

∀i ∈ Z, ∀n ∈ N, h i φ n+1 i -φ n i ∆t + v( φ n i+ 1 2 -φ n i-1 2 ) = vh i R n i , (30)
with

φ n i+ 1 2 = φ n i + α i ( φ n i+1 -φ n i-1 ) + θ i+ 1 2 ( φ n i -φ n i+1 ), φ n i-1 2 = φ n i-1 + α i-1 ( φ n i -φ n i-2 ) + θ i-1 2 ( φ n i-1 -φ n i ).
Expressing φ n j = φ ini (x j -vn∆t), for any j = i -2, . . . , i + 1, as a Taylor expansion from x i,0 = x i -vn∆t, we find that

φ n i+ 1 2 = φ ini (x i,0 ) + ( 1 2 (h i -∆tv) -βh)φ ini (x i,0 ) + h i R n i,+ and φ n i-1 2 = φ ini (x i,0 ) + ( 1 2 (-h i -∆tv) -βh)φ ini (x i,0 ) + h i R n i,-, with ∀i ∈ Z, ∀n ∈ N, |R n i,± | ≤ C 0,
x h, for some C 0,x > 0, depending only on a bound of |φ ini |, on v and on β. This leads to

v( φ n i+ 1 2 -φ n i-1 2 ) = h i v(φ ini (x i,0 ) + R n i,+ -R n i,-).
We have

φ n+1 i -φ n i ∆t = -vφ ini (x i,0 ) + vR n,+ i , with ∀i ∈ Z, ∀n ∈ N, |R n,+ i | ≤ C 0,t
∆t, for some C 0,t > 0, depending only on a bound of |φ ini | and on v. Gathering the above results, we get that there exists C 0 > 0, depending only on a bound of |φ ini |, on v and on β, such that

∀i ∈ Z, ∀n ∈ N, |R n i | ≤ C 0 h.
Accounting for the fact that the support of φ ini is compact, and assuming that h ≤ 1, we get that there exists C 1 > 0, only depending on a bound of |φ ini |, on the support of φ ini , on v and on β, such that ∀n ∈ N,

i∈Z h i (R n i ) 2 ≤ C 1 h 2 . ( 31
)
We then denote e n i = φ n i -φ n i , and we obtain, by subtracting ( 28) from [START_REF] Venkatakrishnan | Convergence to steady state solutions of the euler equations on unstructured grids with limiters[END_REF],

∀i ∈ Z, ∀n ∈ N, h i e n+1 i -e n i ∆t + v(e n i+ 1 2 -e n i-1 2 ) = vh i R n i , (32)
with

e n i+ 1 2 = e n i + α i (e n i+1 -e n i-1 ) + θ i+ 1 2 (e n i -e n i+1
), e n i-

1 2 = e n i-1 + α i-1 (e n i -e n i-2 ) + θ i-1 2 (e n i-1 -e n i ), and ∀i ∈ Z, e 0 i = 0. ( 33
)
Let us notice that, due to the compact support assumption, for any n ∈ N, the number of i ∈ Z such that e n i = 0 or R n i = 0 is finite. We multiply (32) by ∆te n i , and we sum on i ∈ Z and on n = 0, . . . , N -1. This leads to T 1 + T 2 = T 3 with

T 1 = N -1 n=0 i∈Z h i (e n+1 i -e n i )e n i , T 2 = ∆tv N -1 n=0 i∈Z (e n i+ 1 2 -e n i- 1
2

)e n i ,
and

T 3 = ∆tv N -1 n=0 i∈Z h i e n i R n i .
We have, thanks to a discrete integration by parts, that

T 2 = ∆tv N -1 n=0 i∈Z e n i+ 1
2 (e n i -e n i+1 ).

Since there holds

e n i (e n i -e n i+1 ) = 1 2 (e n i ) 2 - 1 2 (e n i+1 ) 2 + 1 2 (e n i -e n i+1 ) 2 ,
we get

T 2 = ∆tv N -1 n=0 i∈Z ( 1 2 + θ i+ 1 2 -α i )(e n i -e n i+1 ) 2 + α i (e n i -e n i-1 )(e n i -e n i+1
) .

Using the inequality

(e n i -e n i-1 )(e n i -e n i+1 ) ≥ - 1 2 ((e n i -e n i-1 ) 2 + (e n i -e n i+1 ) 2 ),
we obtain, using ( 27),

T 2 ≥ ∆tv N -1 n=0 i∈Z ( 1 2 + θ i+ 1 2 -α i - 1 2 (α i + α i+1 ))(e n i -e n i+1 ) 2 ≥ 1 2 ∆tv N -1 n=0 i∈Z (e n i -e n i+1 ) 2 .
Using the relation

(e n+1 i -e n i )e n i = 1 2 (e n+1 i ) 2 - 1 2 (e n i ) 2 - 1 2 (e n+1 i -e n i ) 2 ,
we get

T 1 = i∈Z h i 2 (e N i ) 2 - N -1 n=0 i∈Z h i 2 (e n+1 i -e n i ) 2 .
We now use [START_REF] Xu | A modified finite volume method for convection-diffusion-reaction problems[END_REF] and [START_REF] Tambue | An exponential integrator for finite volume discretization of a reactionadvection-diffusion equation[END_REF]. We get

h i |e n+1 i -e n i | ≤ v∆t β 2 |e n i -e n i+1 | + ( 3 2 + β 2 )|e n i -e n i-1 | + 1 2 |e n i-1 -e n i-2 | + h i |R n i | .
Applying the Cauchy-Schwarz inequality, we obtain

-h i (e n+1 i -e n i ) 2 ≥ - v 2 ∆t 2 h i (2 + 2β 2 + h i ) × β 2 (e n i -e n i+1 ) 2 + ( 3 2 + β 2 )(e n i -e n i-1 ) 2 + 1 2 (e n i-1 -e n i-2 ) 2 + h i (R n i ) 2 .
This yields, under the condition h < 1,

T 1 ≥ i∈Z h i 2 (e N i ) 2 - βv 2 ∆t 2 2h (2 + 2β 2 + β) (2 + 2β 2 ) N -1 n=0 i∈Z (e n i -e n i+1 ) 2 + N -1 n=0 i∈Z h i (R n i ) 2 .
We now choose C cfl in [START_REF] Tamamidis | A new upwind scheme on triangular meshes using the finite volume method[END_REF] such that the terms in (e n i -e n i+1 ) 2 in T 1 + T 2 are non-negative. Indeed, there holds

v∆t 2 - βv 2 ∆t 2 2h (2 + 2β 2 + β)(2 + 2β 2 ) ≥ 0,
under condition [START_REF] Tamamidis | A new upwind scheme on triangular meshes using the finite volume method[END_REF], provided that

C cfl = 1 β(2 + 2β 2 + β)(2 + 2β 2 ) which is such that C cfl < 1 β . (34)
We thus obtain the existence of C 2 , only depending on β such that

T 1 + T 2 ≥ i∈Z h i 2 (e N i ) 2 -C 2 v∆t N -1 n=0 i∈Z h i (R n i ) 2 .
Gathering the above relations, we obtain

i∈Z h i 2 (e N i ) 2 -C 2 v∆t N -1 n=0 i∈Z h i (R n i ) 2 ≤ v∆t N -1 n=0 i∈Z h i |e n i ||R n i |.
We apply the inequality ab ≤ 1 2 (a 2 + b 2 ) to the right-hand-side of the above inequality, and we multiply by 2, which yields

i∈Z h i (e N i ) 2 ≤ v∆t N -1 n=0 i∈Z h i (e n i ) 2 + (1 + 2C 2 )v∆t N -1 n=0 i∈Z h i (R n i ) 2 .
Applying (31), we deduce

i∈Z h i (e N i ) 2 ≤ v∆t N -1 n=0 i∈Z h i (e n i ) 2 + (1 + 2C 2 )vT C 1 h 2 .
We remark that all the preceding inequalities hold, replacing N by any value in {1, . . . , N }.

Therefore we can apply Lemma 3.1, and we obtain

i∈Z h i (e N i ) 2 ≤ e vT (1 + 2C 2 )vT C 1 h 2 .
This proves [START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF], and concludes the proof of the theorem.

Theorem 3.1 (Discrete Gronwall lemma) Let N ∈ N , T > 0 and ∆t = T /N . Let b ≥ 0 and v > 0 be given, and let (a n ) n=0,...,N be non-negative real values with a 0 = 0 and

a n ≤ v∆t n-1 k=0 a k + b for n ∈ {0, . . . , N } Then a N ≤ be vT .

Proof

We set A n = n-1 k=0 a k for n = 1, . . . , N and A 0 = 0. We then have

A n -A n-1 ≤ v∆tA n-1 + b,
and therefore

A n ≤ (1 + v∆t) n A 0 + b((1+v∆t) n -1) 1+v∆t-1 = b((1+v∆t) n -1) v∆t ≤ b(e vT -1)
v∆t since A 0 = 0 (recall that (1 + x N ) n ≤ e x for x ≥ 0 and n ≤ N ). This finally leads to a N ≤ v∆tA N + b ≤ b(e vT -1) + b = be vT .

Test Problems

In this section, three examples of pure convection and convection-dominated diffusion problems are examined to evaluate the accuracy and robustness of four piecewise linear reconstruction methods. These examples are (1) the mixing of hot with cold fronts, (2) the rotation of a Gaussian pulse, and (3) triangular wave flow. All examples in this section were tested on structured triangular grids. For all numerical solutions, we use a triangular grid discretization of the domain, and assume that the limiter function Ψ = 1 (equation ( 12)).

Mixing of hot with cold fronts

The first example is a pure-convection problem of the mixing of hot with a cold fronts [START_REF] Phongthanapanich | A characteristic-based finite volume element method for convection-diffusion-reaction equation[END_REF][START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF]. The computational domain is Ω = (-4, -4) × (4, 4) and the velocity field is given by

v = - y r f t 0.385 i + x r f t 0.385 j, (35)
where r = x 2 + y 2 is the distance from the origin and f t = tanh(r) cosh 2 (r) . The initial condition is specified by

φ 0 (x) = -tanh( y 2 ). (36)
The exact solution is

φ(x, t) = -tanh y 2 cos f t 0.385r t - x 2 sin f t 0.385r t . (37)
The 2D and 3D plots of the exact solution at the final time t = 4 are shown in Figures 1 and2, respectively. The corresponding numerical solutions obtained from the three uniform grids S1 (16×16), S2 (32 × 32), and S3 (64 × 64) compared to the exact solution along the line y = 0 are depicted in Figures 345, respectively. The L 2 -norm error of the solutions obtained by using four piecewise linear reconstruction techniques are shown in Table 1. We can see that the Frink and Holmes-Connell methods gave the most accurate solutions on grids S1 and S2, while the least-squares method gave the most accurate solution on grid S3. The Green-Gauss method gave the worst solution on all grid sizes.

Rotation of Gaussian pulse

The second example is the rotation of a Gaussian pulse around the domain Ω = (-0.5, -0.5) × (0.5, 0.5). We started the test using the pure convection problem with the initial condition φ(x, 0) given as [START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF] φ(x, 0) = exp - The 2D and 3D plots of the exact solution at the final time t = π/2 are shown in Figures 6 and7, respectively. The comparison between exact and numerical solutions obtained from the grids S2, S3, and S4 (128×128) along the line y = 0 are illustrated in Figures 8-10, respectively. The comparison of the L 2 -norm errors in the solutions is given in Table 2. We may conclude that the last-squares method gave the most accurate solution and the Green-Gauss method gave the least accurate solution. On grid S4, the L 2 -norm error value of the solution obtained using the least-square methods is about 100% better than the Frink and the Holmes-Connell methods, and 440% better than the Green-Gauss method. This example is repeated again as a convection-dominated diffusion problem by setting the diffusion coefficient of = 10 -4 . The exact solution of this test at the same final time is determined by [START_REF] Wang | Development of cfl-free, explicit schemes for multidimensional advection-reaction equations[END_REF][START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF] as 3. As we found for the pure convection case, the least-squares method gave the most accurate solution for all grids. On grid S4, the L 2 -norm error values of the solutions obtained by the least-squares method is about 123% better than the Frink and the Holmes-Connell methods, and 530% better than the Green-Gauss method. Finally, we note that the Frink and the Holmes-Connell methods

(x -x c ) 2 + (y -y c ) 2 2σ 2 , (38)
φ(x, 0) = 2σ 2 2σ 2 + 4tε exp - (x -x c ) 2 + (y -y c ) 2 2σ 2 + 4tε . (40)

Triangular wave flow

The last example is the so-called triangular wave flow problem [START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF]. This is a pure-convection problem with an initial discontinuous flow profile. The computational domain is a unit square

The velocity field is given by v = 0.05i. This example is used to evaluate the robustness of the piecewise linear reconstruction techniques against the scalar discontinuity on four grids (S1,S2,S3,S4). For the coarse grid S1 (Figure 16), all methods gave nearly identical solutions without oscillation. For grid size S2, a slight oscillation occurs as shown in Figure 17. For S2, the Frink and the Holmes-Connell methods gave identical solutions along the line y = 0.5 with relatively small oscillations, while the solutions for the Green-Gauss and the least-squares methods gave larger oscillations. For the finer grids S3 and S4, the oscillations were apreciably larger as shown in Figures 18 and19. Finally, the results in Figures 18 and19 show that the Frink and the Holmes-Connell methods gave oscillating solutions downstream of the discontinuity while the Green-Gauss and the least-squares gave an oscillating solution upstream of the discontinuity.

Figure 1 :

 1 Figure 1: Exact solution at the final time of problem 4.1: 2D plot.

Figure 2 :

 2 Figure 2: Exact solution at the final time of problem 4.1: 3D plot.

Figure 3 :Figure 4 :

 34 Figure 3: Comparison of the solutions of problem 4.1 on grid S1

Figure 5 :

 5 Figure 5: Comparison of the solutions of problem 4.1 on grid S3

Figure 6 :

 6 Figure 6: 2D plot of the exact solution at the final time of problem 4.2: Pure convection.

Figure 7 :

 7 Figure 7: 3D plot of the exact solution at the final time of problem 4.2: Pure convection.

Figure 8 :Figure 9 :

 89 Figure 8: Comparison of the solutions of problem 4.2 on grid S2: Pure convection.

Figure 10 :

 10 Figure 10: Comparison of the solutions of problem 4.2 on grid S4: Pure convection.

Figure 11 :

 11 Figure 11: 2D plot of the exact solution at the final time of problem 4.2: Convection-Diffusion.

Figure 12 :Figure 13 :

 1213 Figure 12: 3D plot of the exact solution at the final time of problem 4.2: Convection-Diffusion.

Figure 16 :Figure 17 :

 1617 Figure 16: Comparison of the solutions of problem 4.3 on grid S1.

Figure 18 :

 18 Figure 18: Comparison of the solutions of problem 4.3 on grid S3.

Figure 19 :

 19 Figure 19: Comparison of the solutions of problem 4.3 on grid S4.

Table 1 :

 1 Comparison of L 2 -norm errors for problem 4.1.

	||e|| L 2	S1	S2	S3
	Frink	0.004848248 0.001285789 0.0002236133
	Holmes-Connell 0.004848264 0.001285792 0.0002236136
	Green-Gauss	0.005117222 0.001370439 0.0002785007
	Least-Squares	0.004879665 0.001295393 0.0002080393

Table 2 :

 2 Comparison of the L 2 -norm errors for problem 4.2: Pure convection.

	||e|| L 2	S2	S3	S4
	Frink	0.001100730	0.0002354509 0.00003324699
	Holmes-Connell 0.001100731	0.0002354509 0.00003324699
	Green-Gauss	0.001080318	0.0002938341 0.00007430981
	Least-Squares	0.0008085976 0.0001357213 0.00001689011

Table 3 :

 3 Comparison of L 2 -norm errors for problem 4.2: Convection-Diffusion.

	||e|| L 2	S2	S3	S4
	Frink	0.0009283025 0.0001863148 0.00002602896
	Holmes-Connell 0.0009283033 0.0001863148 0.00002602896
	Green-Gauss	0.0009172858 0.0002428729 0.00006169527
	Least-Squares	0.0006557170 0.0001008259 0.00001165540
	t = 0 are prescribed by
	φ(0, y) =	 2 (y -0.25) , 0.25 ≤ y ≤ 0.50   2 (0.75 -y) , 0.50 ≤ y ≤ 0.75	.
		  0,	otherwise

Acknowledgments

The first author is pleased to acknowledge the Science Achievement Scholarship of Thailand (SAST) for supporting this research work.

of Ω = (0, 1) × (0, 1). For the initial condition, φ 0 (x) is set to zero. The boundary conditions at

Conclusions

This paper presents an explicit characteristic finite volume method for solving the convectiondiffusion equation in a two-dimensional domain. The method combines the characteristic-based scheme developed for the finite element method with the finite volume method for deriving the discretized equations from the convection-diffusion equation. The convergence of numerical solutions of the finite volume equations has been studied for four piecewise linear reconstruction techniques, namely, Frink [START_REF] Frink | Fast upwind solver for the euler equations on three-dimensional unstructured meshes[END_REF], Holmes-Connell [START_REF] Holmes | Solution of the 2d navier-stokes equations on unstructured adaptive grids[END_REF], Green-Gauss [START_REF] Barth | The design and application of upwind schemes on unstructured meshes[END_REF], and least-squares [START_REF] Barth | A 3-d upwind euler solver for on unstructured meshes[END_REF]. An analysis of the convergence of the scheme for a one-dimensional problem has been given in Section 3. The accuracy and robustness of the four linear construction techniques has been tested on three examples, namely, mixing of hot with cold fronts [START_REF] Phongthanapanich | A characteristic-based finite volume element method for convection-diffusion-reaction equation[END_REF][START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF], rotation of a Gaussian pulse [START_REF] Phongthanapanich | Explicit characteristic finite volume method for convection-diffusion equation on rectangular grids[END_REF] and triangular wave flow [START_REF] Theeraek | Solving convection-diffusionreaction equation by adaptive finite volume element method[END_REF], for the special case of a triangular grid discretization of the domain. It was found that the least-squares method provided the most accurate solution for


Kurt Appert

T Hellsten

Hinrich Lütjens

Olivier Sauter

Jan Vaclavik

Laurent Villard

The physics of wave-coupling in low-frequency plasma heating

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.