# A central-upwind geometry-preserving method for

hyperbolic conservation laws on the sphere

###### Abstract

We introduce a second-order, central-upwind finite volume method for the discretization of nonlinear hyperbolic conservation laws posed on the two-dimensional sphere. The semi-discrete version of the proposed method is based on a technique of local propagation speeds and it is free of any Riemann solver. The main advantages of our scheme are the high resolution of discontinuous solutions, its low numerical dissipation, and its simplicity for the implementation. The proposed scheme does not use any splitting approach, which is applied in some cases to upwind schemes in order to simplify the resolution of Riemann problems. The semi-discrete form of the scheme is strongly linked to the analytical properties of the nonlinear conservation law and to the geometry of the sphere. The curved geometry is treated here in an analytical way so that the semi-discrete form of the proposed scheme is consistent with a geometric compatibility property. Furthermore, the time evolution is carried out by using a total-variation-diminishing Runge-Kutta method. A rich family of (discontinuous) stationary solutions is available for the problem under consideration when the flux is nonlinear and foliated (as identified by the author in an earlier work). We present here a series of numerical examples, obtained by considering non-trivial steady state solutions and this leads us to a good validation of the accuracy and efficiency of the proposed central-upwind finite volume method. Our numerical tests confirm the stability of the proposed scheme and clearly show its ability to capture accurately discontinuous steady state solutions to nonlinear hyperbolic conservation laws posed on the sphere.

A central-upwind geometry-preserving method for

hyperbolic conservation laws on the sphere

Department of Civil Engineering, University of Ottawa, 161 Louis Pasteur, Ottawa,

Ontario, K1N6N5, Canada, Email: abelj016@uottawa.ca

Laboratoire Jacques-Louis Lions & Centre National de la Recherche Scientifique,

Université Pierre et Marie Curie (Paris 6), 4 Place Jussieu, 75258 Paris, France.

Email: contact@philippelefloch.org

Keywords: Hyperbolic conservation law; shock wave; geometry-compatible flux; central-upwind scheme

## 1 Introduction

Solutions to hyperbolic partial differential equations generally develop discontinuities in finite time, even from smooth initial conditions. Various classes of so-called shock-capturing schemes have been proposed. In particular, upwind and central schemes have been used to numerically solve these equations. Generally, it can be stated that the difference between these schemes is that upwind methods use characteristic-related information, while central methods do not. The use of characteristic information in upwind schemes can improve the results but renders these schemes, in some cases, computationally expensive. Central schemes are widely used (see e.g.[28]) after the pioneering work of Nessyahu and Tadmor [26], where a second-order finite volume central method on a staggered grid in space-time was first proposed. This strategy leads to high-resolution and the simplicity of the Riemann-solver free method. As observed in Kurganov and Tadmor [13], this scheme suffers from excessive numerical viscosity when a small time step is considered.

In order to improve the performance of central schemes, some characteristic information can still be used. Kurganov et al. [16] proposed the central-upwind schemes which are based on information obtained from the local speeds of wave propagation. The central-upwind schemes can be considered as a generalization of central schemes originally developed by Kurganov and Tadmor [13, 14], Kurganov and Levy [15], and Kurganov and Petrova [17]. The central-upwind schemes are simple, since they use no Riemann solvers, and they have proven their effectiveness in multiple studies, as shown in [20, 21, 23, 24, 22, 18]. Kurganov and Petrova [19] extended the central-upwind schemes to triangular grids for solving two-dimensional Cartesian systems of conservation laws. Next, Beljadid et al. [5] proposed a two-dimensional well-balanced and positivity preserving cell-vertex central-upwind scheme for the computation of shallow water equations with source terms due to bottom topography.

Several studies have been recently developed for hyperbolic conservation laws posed on curved manifolds. The solutions of conservation laws including the systems on manifolds and on spacetimes were studied in [27, 25] and by LeFloch and co-authors [2, 3, 7, 6] and [11]–[12]. More recently, hyperbolic conservation laws for evolving surface were investigated by Dziuk, Kroöner and Müller [8], Giesselman [10], and Dziuk and Elliott [9]. Earlier on, for such problems, Ben-Artzi and LeFloch [7] and LeFloch and Okutmustur [12] established a general well-posedness theory for conservation laws on manifolds.

Burgers equation provides a simple, yet challenging equation, which admits discontinuous solutions and it provides a simplified setup for the design and validation of shock-capturing numerical methods. Burgers equation and its generalizations to a curved manifold have been widely used in the physical and mathematical literature. In [4], we have used a class of Burgers-type equations on the sphere and adopted the methodology first proposed by Ben-Artzi, Falcovitz, and LeFloch [6], which uses second–order approximations based on generalized Riemann problems. In [4], a scheme was proposed which uses piecewise linear reconstructions based on solution values at the center of the computational cells and on values of Riemann solutions at the cell interfaces. A second-order approximation based on a generalized Riemann solver was then proposed, together with a total variation diminishing Runge-Kutta method (TVDRK3) with operator splitting for the temporal integration.

The finite volume methods developed in [4] and [6] are strongly linked to the structure of governing equation. The geometric dimensions are considered in an analytical way which leads to discrete forms of schemes that respect exactly the geometric compatibility property. The splitting approach which is used in these schemes simplifies the resolution of the Riemann problem but it increases the computational cost.

In the present study, we propose a new finite volume method which is less expensive in terms of computational cost. This scheme is free of any Riemann solver and does not use any splitting approach, while such a splitting is widely used in upwind schemes when one needs to simplify the resolution of Riemann problems. The proposed paper provides the first study of central-upwind schemes for conservation laws on a curved geometry.

Burgers equation and its generalizations will be used in the present paper in order to develop and validate the new finite volume method. We design in full detail a geometry-compatible central-upwind scheme for scalar nonlinear hyperbolic conservation laws on the sphere. This system has a simple appearance but it generates solutions that have a very rich wave structure (due to the curved geometry) and its solutions provide an effective framework for assessing numerical methods. Our goal is to develop and validate a finite volume method which is free of any Riemann problem and is consistent with the geometric compatibility (or divergence free) condition, at the discrete level. As we prove it, the proposed scheme is efficient and accurate for discontinuous solutions and implies only negligible geometric distorsions on the solutions.

An outline of the paper is as follows. In Section 2, the governing equations related to this study are presented. Section 3 is devoted to the derivation of the semi-discrete version of our scheme. In Section 4, the coordinate system and the non-oscillatory reconstruction are described. In Section 5, we present the geometry-compatible flux vectors and some particular steady state solutions as well as confined solutions, which will be used to validate the performance of the proposed method. In Section 6, we demonstrate the high-resolution of the proposed central-upwind scheme thanks to a series of numerical experiments. Finally, some concluding remarks are provided.

## 2 Governing equations

We consider nonlinear hyperbolic equations posed on the sphere and based on the flux-vector , depending on the function and the space variable . This flux is assumed to satisfy the following geometric compatibility condition: for any arbitrary constant value ,

(1) |

and that the flux takes the form

(2) |

where is the unit normal vector to the sphere and the function is a vector field in , restricted to and defined by

(3) |

Here, is a smooth function depending on the space variable and the state variable . Observe that (for instance by Claim 2.2 in [6]) the conditions and for the flux vector are sufficient to ensure the validity of the geometric compatibility condition .

We are going here to develop and validate a new geometry-preserving central-upwind scheme which approximates solutions to the hyperbolic conservation law

(4) |

where is the divergence of the vector field . Given any data prescribed on the sphere, we consider the following initial condition for the unknown function

(5) |

Equation can be rewritten, using general local coordinates and the index of summation , in the form

(6) |

or

(7) |

where in local coordinates , the derivatives are denoted by , are the components of the flux vector and is the metric.

The conservation law becomes as follows

(8) |

where . This form will be used in the derivation of the semi-discrete form of the proposed scheme. For latitude-longitude grid on the sphere, the divergence operator of the flux vector is

(9) |

where and are the flux components in the latitude () and longitude directions on the sphere respectively.

## 3 Derivation of the proposed method

### 3.1 Discretization of the divergence operator

The derivation of the new central-upwind scheme will be described in detail for the three steps: reconstruction, evolution, and projection. We will develop and give a semi-discrete form of the proposed method for general computational grid used to discretize the sphere. We assume the discretization of the sphere , where are the computational cells with area . We denote by the number of cell sides of and by the neighboring computational cells that share with the common sides , respectively. The length of each cell-interface is denoted by . The discrete value of the state variable inside the computational cell at a point is denoted by at step . The longitude and latitude coordinates of the suitable point to use inside each computational cell are presented in Section 4.2. These coordinates should be chosen according to the reconstruction of the state variable over the computational cells used on the sphere. Finally, we use the notations and for the time step and the time at step , respectively. Note that the development of the first order in time is sufficient to have the exact semi-discrete form of the proposed scheme. The resulting ODE can be numerically solved using a higher-order SSP ODE solver as Runge-Kutta of multistep methods. In the numerical experiments, the third-order TVD Runge-Kutta method proposed by Shu and Osher [29] is used.

In this section, we will present a general form of the discretization of the divergence operator for general computational grid on the sphere. The approximation of the flux divergence can be written using the divergence theorem as

(10) | ||||

where is the unit normal vector to the boundary of the computational cell and is the infinitesimal length along .

The scalar potential function is used to obtain the following approximation along each side of the computational cell .

###### Claim 3.1

For a three-dimensional flux given by , where is a smooth function in the neighborhood of the sphere , the total approximate flux through the cell interface is given by

(11) |

where and are the initial and final endpoints of the side using the sense of integration and is the estimate value of the variable along the side .

Namely, the flux vector is written in the form and we can derive the approximation of the integral along each cell side of

(12) | ||||

where is the unit vector tangent to the boundary .

###### Remark 3.2

Using the discrete approximations based on Claim 3.1, if a constant value of the state variable is considered one obtains

(13) | ||||

This confirms that the discrete approximation of the divergence operator respects the divergence free condition which is the geometric requirement that the proposed scheme should satisfy.

### 3.2 Reconstruction method

In the following, we will present the reconstruction of the proposed central-upwind scheme. The semi-discrete form of the proposed scheme for Equation will be derived by using the approximation of the cell averages of the solution. At each time the computed solution is

(14) |

where .

The discrete values of the solution at time are used to construct a conservative piecewise polynomial function with possible discontinuities at the interfaces of the computational cells

(15) |

where is a polynomial in two variables ( and ), and is the characteristic function which is defined using the Kronecker symbol and for any point of spatial coordinate inside the computational cell we consider .

The maximum of the directional local speeds of propagation of the waves at the kth interface inward and outward of the computational cell are denoted by and , respectively. When the solution evolves over a time step , the discontinuities move inward and outward the kth interface of the computational cell with maximum distances and , respectively. These distances of propagation are used at the computational cells to delimit different areas in which the solution still smooth and the areas in which the solution may not be smooth when it evolves from the time level to .

We define the domain as the part inside the cell in which the solution still smooth, see Figure 1. Two other types of domains are defined, the first type includes the rectangular domains , , along each side of of width and length and the second type includes the domains denoted by , , around the cell vertices of computational cells. These domains are decomposed into two sub-domains and , where the sub-domains with the superscripts plus signs “” and minus signs “” are the domains inside and outside of the cell , respectively. For purely geometrical reasons, the areas of the three types of sub-domains are of orders , and .

We consider the projection of the flux vector according to the orthogonal to the kth cell interface

(16) |

where is the unit normal vector to the cell interface and has the components which are used in Equation (8).

The one-sided local speeds of propagation of the waves at the kth cell interface , inward and outward the computational cell , are estimated by

(17) | ||||

where is the value of the state variable at the midpoint of , which is obtained from the non-oscillatory reconstruction for the computational cell and is the value of at the same point using the non-oscillatory reconstruction for the neighboring cell .

### 3.3 Evolution and projection steps

The computed cell averages of the numerical solution at time step over the computational cells are used to obtain the piecewise linear reconstruction which should satisfy the following conservative requirement

(18) |

The average of the function over the domain is denoted by

(19) |

Note that it is possible to derive the fully discrete form of the proposed scheme but it is impractical to use and for simplicity, we will develop the semi-discrete form of the scheme. The ODE for approximating the cell averages of the solutions is derived by tending the time step to zero. This eliminates some terms because of their orders and we keep the more consistent terms

(20) |

Since the areas of domains with are of order we obtain

(21) |

This approximation allows us to deduce that the third term on the right-hand side of Equation is of order and the result for the limit of this term vanishes for the ODE.

The second term in Equation , in which we use the rectangular domains , will be estimated by using the assumption that the spatial derivatives of are bounded independently of . Under this assumption the following Claim gives an estimation of this term with an error of order for each .

###### Claim 3.3

Consider the reconstruction given by , its evolution over the global domain, and the definitions given in Section 3.2 for the domains and . If we assume that the spatial derivatives of are bounded independently of , then

(22) |

The proof is as follows. It is obvious that for the case or equation is valid. We assume that and we consider

We have

(23) |

where is the length of the domain and is a variable according to the orthogonal outward axis to the kth cell interface, see Figures 1 and 3.

One obtains after the change of variable in the first integral of the last equality in

Using the mean value theorem to the function we obtain

where .

We denote by the upper bound value of the spatial derivative of the function over the domain . Therefore we obtain

Since , and both the areas and are of order we obtain . This completes the proof.

Using Equation in Claim 3.3 one obtains

(24) |

Therefore Equation can be written as

(25) |

where

(26) |

In order to derive the semi-discrete form of the proposed scheme from Equation , one needs to compute the average values and . To compute , Equation is integrated over the space-time control volume . After integration by parts and applying the divergence theorem to transform the surface integral of the divergence operator to the boundary integral and using the approximation (11) of the flux through the cell interfaces, the following equations are obtained

(27) |

and

(28) |

where , , are the four edges of the rectangular domain , and are the initial and final endpoints of the cell interface , and as mentioned before and are the piecewise polynomial reconstructions in the computational cells and respectively at time .

The term on the right-hand side of Equation of order corresponds to the global result of the integration along the two edges of the domain having the length , and the rest of the integration due to the difference between the length of the domain and the length of the cell interface .

In order to compute the spatial integrals in Equation , the Gaussian quadrature can be applied. In our case, the midpoint rule is used for simplicity

(29) |

(30) | ||||

Therefore

(31) | ||||

Now, the average value will be computed. Equation is integrated over the space-time control volume and after integration by parts, using the divergence theorem to transform the surface integral to boundary integral and Equation (11) one obtains

(32) |

Using the previous equality we obtain

(33) |

which leads to

(34) |

Equations and are used together to obtain the following semi-discrete form

(35) |

This equation can be rewritten in the following form

(36) |

which can be rewritten as follows

(37) |

where and are given by

(38) |

The function is defined in the form (38) in order to be consistent with the total approximate flux through the cell interface
as presented by Equation (11) in Claim 3.1.

Remark 2. If the value of in Equation (37) is zero or very close to zero (smaller than in our numerical experiments), we avoid division by zero or by a very small number using the following approximations

(39) |

These approximations are obtained using similar extreme distances of the propagation of the waves at the cell interface inward and outward the computational cell to define the domains , and .

### 3.4 The geometry-compatible condition

In the semi-discrete form (37) and (38) of the proposed scheme, if we consider a constant value of the function , the second term in the right-hand side of Equation (37) vanishes. For this constant function we obtain for each interface cell

(40) |

and

(41) |

The first term in the right-hand side of Equation (37) becomes

(42) |

Since we have

(43) |

we conclude that the first term on the right-hand side of Equation (37) will be canceled which confirms that the proposed scheme respects the geometry-compatibility condition.

###### Remark 3.4

In the proposed central-upwind finite volume method, the midpoint rule was used to compute the spatial integrals. In order to improve the accuracy of the proposed scheme, the Gaussian quadrature can be used instead of the midpoint rule. The Gaussian quadrature will not have any impact on the geometry-compatibility condition of the proposed scheme.

## 4 Formulation using the latitude-longitude grid on the sphere

### 4.1 Computational grid on the sphere

The geometry-compatible scheme was developed in the previous section for scalar nonlinear hyperbolic conservation laws using general grid on the sphere. However, in order to prevent oscillations an appropriate piecewise linear reconstruction should be proposed according to the computational grid used in the proposed method. In the following, we will describe the computational grid and the non-oscillatory piecewise linear reconstruction used in our numerical experiments.

The position of each point on the sphere can be represented by its longitude and its latitude . The grid considered in our numerical examples is shown in Figure 2. The coordinates are singular at the south and north poles, corresponding to and , respectively. The Cartesian coordinates are denoted by for a standard orthonormal basis vectors , and .

The unit tangent vectors in the directions of longitude and latitude at each point on the sphere of coordinates are given as follows

(44) | ||||