High order numerical schemes for non-local conservation laws

# High order numerical schemes for one-dimension non-local conservation laws

Christophe Chalons Paola Goatin  and  Luis M. Villada
July 7, 2019
###### Abstract.

This paper focuses on the numerical approximation of the solutions of non-local conservation laws in one space dimension. These equations are motivated by two distinct applications, namely a traffic flow model in which the mean velocity depends on a weighted mean of the downstream traffic density, and a sedimentation model where either the solid phase velocity or the solid-fluid relative velocity depends on the concentration in a neighborhood. In both models, the velocity is a function of a convolution product between the unknown and a kernel function with compact support. It turns out that the solutions of such equations may exhibit oscillations that are very difficult to approximate using classical first-order numerical schemes. We propose to design Discontinuous Galerkin (DG) schemes and Finite Volume WENO (FV-WENO) schemes to obtain high-order approximations. As we will see, the DG schemes give the best numerical results but their CFL condition is very restrictive. On the contrary, FV-WENO schemes can be used with larger time steps. We will see that the evaluation of the convolution terms necessitates the use of quadratic polynomials reconstructions in each cell in order to obtain the high-order accuracy with the FV-WENO approach. Simulations using DG and FV-WENO schemes are presented for both applications.

###### Key words and phrases:
Scalar conservation laws, non-local flux, sedimentation models, traffic flow model, Galerkin discontinuous schemes, finite volume schemes.
###### 2000 Mathematics Subject Classification:
35L65,90B20,76T20,65M08,65M60.
Laboratoire de Mathématiques de Versailles, UMR 8100, Université de Versailles Saint-Quentin-en-Yvelines, UFR des Sciences, bâtiment Fermat, 45 avenue des Etats-Unis, 78035 Versailles cedex, France. E-mail: christophe.chalons@uvsq.fr
Inria Sophia Antipolis - Méditerranée, Université Côte d’Azur, Inria, CNRS, LJAD, 2004, route des Lucioles - BP 93 06902 Sophia Antipolis Cedex, France E-mail: paola.goatin@inria.fr
GIMNAP-Departamento de Matemáticas, Universidad del Bío-Bío, Casilla 5-C, Concepción Chile and CIMA, Universidad de Concepción, Casilla 160-C, Concepción, Chile. E-Mail: lvillada@ubiobio.cl

## 1. Introduction

This paper is concerned with the design of numerical schemes for the one-dimensional Cauchy problem for non-local scalar conservation law of the form

 {ρt+(f(ρ)V(ρ∗ωη))x=0,x∈R,t>0,ρ(x,0)=ρ0(x)x∈R, (1.1)

where the unknown density depends on the space variable and the time variable , is a given velocity function, and is the usual flux function for the corresponding local scalar conservation law. Here, (1.1) is non-local in the sense that the velocity function is evaluated on a “neighborhood” of defined by the convolution of the density and a given kernel function with compact support. In this paper, we are especially interested in two specific forms of (1.1) which naturally arise in traffic flow modelling [9, 17] and sedimentation problems [8]. They are given as follows.

A non-local traffic flow model. In this context, we follow [9, 17] and consider (1.1) as an extension of the classical Lighthill-Whitham-Richards traffic flow model, in which the mean velocity is assumed to be a non-increasing function of the downstream traffic density and where the flux function is given by

 f(ρ)=ρ,V(ρ)=1−ρ,(ρ∗ωη)(x,t)=∫x+ηxωη(y−x)ρ(y,t)dy. (1.2)

Four different nonnegative kernel functions will be considered in the numerical section, namely , , and with support on for a given value of the real number . Notice that the well-posedness of this model together with the design of a first and a second order FV approximation have been considered in [9, 17].

A non-local sedimentation model. Following [8], we consider under idealized assumptions that equation (1.1) represents a one-dimensional model for the sedimentation of small equal-sized spherical solid particles dispersed in a viscous fluid, where the local solids column fraction is a function of depth and time . In this context, the flux function is given by

 f(ρ)=ρ(1−ρ)α,V(ρ)=(1−ρ)n,(ρ∗ωη)(x,t)=∫2η−2ηωη(y)ρ(x+y,t)dy, (1.3)

where and the parameter satisfies or . The function is the so-called hindered settling factor and the convolution term is defined by a symmetric, nonnegative, and piecewise smooth kernel function with support on for a parameter and . More precisely, the authors define in [8] a truncated parabola by

 K(x)=38(1−x24) for |x|<2,K(x)=0 otherwise,

and set

 ωη(x):=η−1K(η−1x). (1.4)

Conservation laws with non-local terms arise in a variety of physical and engineering applications: besides the above cited ones, we mention models for granular flows [4], production networks [21], conveyor belts [18], weakly coupled oscillators [3], laser cutting processes [15], biological applications like structured populations dynamics [23], crowd dynamics [10, 12] or more general problems like gradient constrained equations [6].
While several analytical results on non-local conservation laws can be found in the recent literature (we refer to [7] for scalar equations in one space dimension, [5, 13, 24] for scalar equations in several space dimensions and [2, 14, 16] for multi-dimensional systems of conservation laws), few specific numerical methods have been developed up to now. Finite volume numerical methods have been studied in [2, 17, 22, 25]. At this regard, it is important to notice that the lack of Riemann solvers for non-local equations limits strongly the choice of the scheme. At the best of our knowledge, two main approaches have been proposed in the literature to treat non-local problems: first and second order central schemes like Lax-Friedrichs or Nassyau-Tadmor [2, 7, 8, 17, 22] and Discontinuous Galerkin (DG) methods [19]. In particular, the comparative study presented in [19] on a specific model for material flow in two space dimensions, involving density gradient convolutions, encourages the use of DG schemes for their versatility and lower computational cost, but further investigations are needed in this direction. Besides that, the computational cost induced by the presence of non-local terms, requiring the computation of quadrature formulas at each time step, motivate the development of high order algorithms.

The aim of the present article is to conduct a comparison study on high order schemes for a class of non-local scalar equations in one space dimension, focusing on equations of type (1.1). In Section 2 we review DG and FV-WENO schemes for classical conservation laws. These schemes will be extended to the non-local case in Sections 3 and 4. Section 5 is devoted to numerical tests.

## 2. A review of Discontinuous Galerkin and Finite Volume WENO schemes for local conservation laws

The aim of this section is to introduce some notations and to briefly review the DG and FV-WENO numerical schemes to solve the classical local nonlinear conservation law

 {ρt+g(ρ)x=0,x∈R,t>0,ρ(x,0)=ρ0(x),x∈R. (2.1)

We first consider a partition of . The points will represent the centers of the cells , and the cell sizes will be denoted by . The largest cell size is . Note that, in practice, we will consider a constant space step so that we will have .

### 2.1. The Discontinuous Galerkin approach

In this approach, we look for approximate solutions in the function space of discontinuous polynomials

 Vh:=Vkh={v:v|Ij∈Pk(Ij)},

where denotes the space of polynomials of degree at most in the element . The approximate solutions are sought under the form

 ρh(x,t)=k∑l=0c(l)j(t)v(l)j(x),v(l)j(x)=v(l)(ζj(x)),

where are the degrees of freedom in the element . The subset constitutes a basis of and in this work we will take Legendre polynomials as a local orthogonal basis of , namely

 v(0)(ζj)=1,v(1)(ζj)=ζj,v(2)(ζj)=12(3ζ2j−1),…,ζj:=ζj(x)=x−xjΔx/2,

see for instance [11, 26].
Multiplying (2.1) by and integrating over gives

 ∫Ijρtvhdx−∫Ijg(ρ)vh,xdx+g(ρ(⋅,t))vh(⋅)∣∣xj+12xj−12=0, (2.2)

and the semi-discrete DG formulation thus consists in looking for , such that for all and all ,

 ∫Ijρhtvhdx−∫Ijg(ρh)vh,xdx+^gj+12v−h(xj+12)−^gj−12v+h(xj−12)=0, (2.3)

where is a consistent, monotone and Lipschitz continuous numerical flux function. In particular, we will choose to use the Lax-Friedrichs flux

 ^g(a,b):=g(a)+g(b)2+αa−b2,α=maxu|g′(u)|.

Let us now observe that if is the -th Legendre polynomial, we have , and , . Therefore, replacing by and by in (2.3), the degrees of freedom satisfy the differential equations

 ddtc(l)j(t)+1al(−∫Ijg(ρhj)ddxv(l)(ζj(x))dx+^g(ρh,−j+12,ρh,+j+12)−(−1)l^g(ρh,−j−12,ρh,+j−12))=0, (2.4)

with

 al=∫Ij(v(l)(ζj(x)))2dx=Δx2l+1,l=0,1,…,k.

On the other hand, the initial condition can be obtained using the -projection of , namely

 c(l)j(0)=1al∫Ijρ0(x)v(l)(ζj(x))dx=2l+12∫1−1ρ0(Δx2y+xj)v(l)(y)dy,l=0,…,k.

The integral terms in (2.4) can be computed exactly or using a high-order quadrature technique after a suitable change of variable, namely

 ∫Ijg(ρhj)ddxv(l)(ζj(x))dx=∫1−1g(ρhj(Δx2y+xj,t))(v(l))′(y)dy.

In this work, we will consider a Gauss-Legendre quadrature with nodes for integrals on

 ∫1−1g(y)dy=NG∑e=1weg(ye),

where are the Gauss-Legendre quadrature points such that the quadrature formula is exact for polynomials of degree until [1].

The semi-discrete scheme (2.4) can be written under the usual form

 ddtC(t)=L(C(t)),

where is the spatial discretization operator defined by (2.4). In this work, we will consider a time-discretisation based on the following total variation diminishing (TVD) third-order Runge-Kutta method [28],

 C(1) = Cn+ΔtL(Cn), C(2) = 34Cn+14C(1)+14ΔtL(C(1)), (2.5) Cn+1 = 13Cn+23C(2)+23ΔtL(C(2)).

Other TVD or strong stability preserving SSP time discretization can be also used [20]. The CFL condition is

 ΔtΔxmaxρ|g′(ρ)|≤CCFL=12k+1,

where is the degree of the polynomial, see [11]. The scheme (2.4) and (2.1) will be denoted RKDG.

### 2.2. Generalized slope limiter

It is well-known that RKDG schemes like the one proposed above may oscillate when sharp discontinuities are present in the solution. In order to control these instabilities, a common strategy is to use a limiting procedure. We will consider the so-called generalized slope limiters proposed in [11]. With this in mind and , we first set

 u−j+12:=¯ρj+mm(ρhj(xj+12)−¯ρj,Δ+¯ρj,Δ−¯ρj)

and

 u+j−12:=¯ρj−mm(¯ρj−ρhj(xj−12),Δ+¯ρj,Δ−¯ρj),

where is the average of on , , , and where is given by the minmod function limiter

 mm(a1,a2,a3)={s⋅minj|aj| if s=sign(a1)=sign(a2)=sign(a3)0 otherwise,

or by the TVB modified minmod function

 ¯¯¯¯¯¯¯¯¯mm(a1,a2,a2)={a1 if |a1|≤Mbh2,mm(a1,a2,a3) otherwise, (2.6)

where is a constant. According to [11, 26], this constant is proportional to the second-order derivative of the initial condition at local extrema. Note that if is chosen too small, the scheme is very diffusive, while if is too large, oscillations may appear.
The values and allow to compare the interfacial values of with respect to its local cell averages. Then, the generalized slope limiter technique consists in replacing on each cell with defined by

Of course, this generalized slope limiter procedure has to be performed after each inner step of the Runge-Kutta scheme (2.1).

### 2.3. The Finite Volume WENO approach

In this section, we solve the nonlinear conservation law (2.1) by using a high-order finite volume WENO scheme [27, 28]. Let us denote by the cell average of the exact solution in the cell  :

 ¯ρ(xj,t):=1Δx∫Ijρ(x,t)dx.

The unknowns are here the set of all which represent approximations of the exact cell averages . Integrating (2.1) over we obtain

 ddt¯ρ(xj,t)=−1Δx(g(ρ(xj+12,t))−g(ρ(xj−12,t))),∀j∈Z.

This equation is approximated by the semi-discrete conservative scheme

 ddt¯ρj(t)=−1Δx(^gj+12−^gj−12),∀j∈Z, (2.7)

where the numerical flux is the Lax-Friedrichs flux and and are some left and right high-order WENO reconstructions of obtained from the cell averages . Let us focus on the definition of . In order to obtain a th-order WENO approximation, we first compute reconstructed values

 ^ρ(r)j+12=k−1∑i=0cir¯ρi−r+j,r=0,…,k−1,

that correspond to possible stencils for . The coefficients are chosen in such a way that each of the reconstructed values is th-order accurate [27]. Then, the th-order WENO reconstruction is a convex combination of all these reconstructed values and defined by

 ρlj+12=k−1∑r=0wr^ρ(r)j+12,

where the positive nonlinear weights satisfy and are defined by

 wr=αr∑k−1s=0αs,αr=dr(ϵ+βr)2.

Here are the linear weights which yield the th-order accuracy, are called the “smoothness indicators” of the stencil , which measure the smoothness of the function in the stencil, and is a small parameter used to avoid dividing by zero (typically ). The exact form of the smoothness indicators and other details about WENO reconstructions can be found in [27].
The reconstruction of is obtained in a mirror symmetric fashion with respect to . The semi-discrete scheme (2.7) is then integrated in time using the (TVD) third-order Runge-Kutta scheme (2.1), under the CFL condition

 ΔtΔxmaxρ|g′(ρ)|≤CCFL<1.

.

## 3. Construction of DG schemes for non-local problems

We now focus on the non-local equation (1.1), for which we set . Since , is a weak solution of the non-local problem (1.1), the coefficients can be calculated by solving the following differential equation,

 ddtc(l)j(t)+1al(−∫Ijf(ρh)V(R)ddxv(l)(ζj(x))dx+ˇfj+12−(−1)lˇfj−12)=0, (3.1)

where is a consistent approximation of at interface . Here, we consider again the Lax-Friedrichs numerical flux defined by

 ˇfj+12=12(f(ρh,−j+12)V(Rh,−j+12)+f(ρh,+j+12)V(Rh,+j+12)+α(ρh,−j+12−ρh,+j+12)), (3.2)

with and where and are the left and right approximations of at the interface . Since is defined by a convolution, we naturally set , so that (3.2) can be written as

 ˇfj+12:=ˇf(ρh,−j+12,ρh,+j+12,Rj+12)=12((f(ρh,−j+12)+f(ρh,+j+12))V(Rj+12)+α(ρh,−j+12−ρh,+j+12)). (3.3)

Next, we propose to approximate the integral term in (3.1) using the following high-order Gauss-Legendre quadrature technique,

 ∫Ijf(ρh)V(R)ddxv(l)(ζj(x))dx = ∫1−1f(ρh(Δx2y+xj,t))V(R(Δx2y+xj,t))(v(l))′(y)dy (3.4) = NG∑e=1wef(ρh(^xe,t))V(R(^xe,t))(v(l))′(ye),

where we have set , being the Gauss-Legendre quadrature points ensuring that the quadrature formula is exact for polynomials of order less or equal to .
It is important to notice that the DG formulation for the non-local conservation law (3.1) requires the computation of the extra integral terms in (3.3) and in (3.4) for each quadrature point, which increases the computational cost of the strategy. For , we can compute these terms as follows for both the LWR and sedimentation non-local models considered in this paper.

Non-local LWR model. For the non-local LWR model, we impose the condition for some , so that we have

 Rj+12 = ∫xj+12+ηxj+12ωη(y−xj+12)ρh(y,t)dy=N∑i=1∫Ij+iωη(y−xj+12)ρhj+i(y,t)dy = N∑i=1k∑l=0c(l)j+i(t)∫Ij+iωη(y−xj+12)v(l)(ζj+i(y))dy = N∑i=1k∑l=0c(l)j+i(t)Δx2∫1−1ωη(Δx2y+(i−12)Δx)v(l)(y)dyΓi,l=N∑i=1k∑l=0c(l)j+iΓi,l,

while for each quadrature point we have

 R(^xe,t) = ∫^xe+η^xeωη(y−^xe)ρh(y,t)dy=∫xj+12^xeωη(y−^xe)ρhj(y,t)dyΓa+ N−1∑i=1∫Ij+iωη(y−^xe)ρhj+i(y,t)dyΓib+∫^xe+ηxj+N−12ωη(y−^xe)ρhj+N(y,t)dyΓc.

The three integrals , and can be computed with the same change of variable as before, namely

 Γa = Δx2∫1yeωη(Δx2(y−ye))k∑l=0c(l)j(t)v(l)(y)dy = k∑l=0c(l)j(t)Δx2∫1yeωη(Δx2(y−ye))v(l)(y)dyΓ(e)0,l=k∑l=0c(l)jΓ(e)0,l, Γib = Δx2∫1−1ωη(Δx2(y−ye)+iΔx)k∑l=0c(l)j+i(t)v(l)(y)dy = k∑l=0c(l)j+iΔx2∫1−1ωη(Δx2(y−ye)+iΔx)v(l)(y)dyΓ(e)i,l=k∑l=0c(l)j+iΓ(e)i,l, Γc = Δx2∫ye−1ωη(Δx2(y−ye)+NΔx)k∑l=0c(l)j+Nv(l)(y)dy = k∑l=0c(l)j+N(t)Δx2∫ye−1ωη(Δx2(y−ye)+NΔx)v(l)(y)dyΓ(e)N,l=k∑l=0c(l)j+NΓ(e)N,l.

Finally we can compute as

 R(^xe,t)=N∑i=0k∑l=0c(l)j+i(t)Γ(e)i,l,

in order to evaluate (3.4).

Non-local sedimentation model. Considering now the non-local sedimentation model, we impose for some , so that we have

 Rj+12=∫xj+12+2ηxj+12−2ηωη(y−xj+12)ρh(y,t)dy=N∑i=−N+1k∑l=0c(l)j+i(t)Γi,l,

and for each quadrature point ,

 R(^xe,t)=∫^xe+2η^xe−2ηωη(y−^xe)ρh(y,t)dy=N∑i=−Nk∑l=0c(l)j+i(t)Γ(e)i,l,

with

 Γ(e)−N,l = Δx2∫1yeωη(Δx2(y−ye)−NΔx)v(l)(y)dy, Γ(e)i,l = Δx2∫1−1ωη(Δx2(y−ye)+iΔx)v(l)(y)dy, Γ(e)N,l = Δx2∫ye−1ωη(Δx2(y−ye)+NΔx)v(l)(y)dy.

Remark. In order to compute integral terms in (3.1) as accurately as possible, the integrals and above, and in particular, the coefficients , must be calculated exactly or using a suitable quadrature formula accurate to at least where is the degree of the convolution term . It is important to notice that the coefficients can be precomputed and stored in order to save computational time.
Finally the semi-discrete scheme (3.1) can be discretized in time using the (TVD) third-order Runge-Kutta scheme (2.1), under the CFL condition

 ΔtΔxmaxρ|∂ρ(f(ρ)V(ρ))|≤CCFL=12k+1,

where is the degree of the polynomial.

## 4. Construction of FV schemes for non-local conservation laws

Let us now extend the FV-WENO strategy of Section 2.3 to the non-local case. We first integrate (1.1) over to obtain

 ddt¯ρ(xj,t)=−1Δx(f(ρ(xj+12,t))V(R(xj+12,t))−f(ρ(xj−12,t))V(R(xj−12,t))),∀j∈Z,

so that the semi-discrete discretization can be written as

 ddt¯ρj(t)=−1Δx(ˇfj+12−ˇfj−12),∀j∈Z, (4.1)

where the use of the Lax-Friedrichs numerical flux gives

 ˇfj+12=ˇf(ρlj+12,ρrj+12,Rj+12)=12((f(ρlj+12)+f(ρrj+12))V(Rj+12)+α(ρlj+12−ρrj+12)).

Recall that and are the left and right WENO high-order reconstructions at point .

At this stage, it is crucial to notice that in the present FV framework, the approximate solution is piecewise constant on each cell , so that a naive evaluation of the convolution terms may lead to a loss of high-order accuracy. Let us illustrate this. Considering for instance the non-local LWR model and using that is piecewise constant on each cell naturally leads to

 Rj+12 = ∫xj+12+ηxj+12ωη(y−xj+12)ρ(y,t)dy=N∑i=1∫Ij+iωη(y−xj+12)ρ(y,t)dy = ΔxN∑i=1¯ρj+i∫Ij+