The finite volume method on a Schwarzschild background

TL;DR

This paper develops a finite volume numerical scheme for nonlinear hyperbolic conservation laws on a Schwarzschild black hole background, proving convergence and stability without boundary conditions at the horizon.

Contribution

It introduces a novel finite volume method tailored for curved spacetime, establishing convergence, existence, and stability results for solutions.

Findings

Scheme converges to an entropy weak solution

No boundary conditions needed at the black hole horizon

Provides a new theory for conservation laws on curved backgrounds

Abstract

We introduce a class of nonlinear hyperbolic conservation laws on a Schwarzschild black hole background and derive several properties satisfied by (possibly weak) solutions. Next, we formulate a numerical approximation scheme which is based on the finite volume methodology and takes the curved geometry into account. An interesting feature of our model is that no boundary conditions is required at the black hole horizon boundary. We establish that this scheme converges to an entropy weak solution to the initial value problem and, in turn, our analysis also provides us with a theory of existence and stability for a new class of conservation laws.

Equations

$$v:\mathcal{M}\to [-1,1]$$

$$\nabla_{\alpha}\big{(}X^{\alpha}(v,\cdot)\big{)}=q(v,\cdot)\quad\text{ in }\mathcal{M}.$$

$$\mathcal{M}$$

$$g=-\Big{(}1-{2M\over r}\Big{)}\,dt^{2}+\Big{(}1-{2M\over r}\Big{)}^{-1}dr^{2}+r^{2}\,g_{S^{n-1}}.$$

$$t\in [0,+\infty ),r\in (2M,+\infty ).$$

$$\big{\{}r=2M\big{\}}\subset\overline{\mathcal{M}}$$

$$X=\Big{(}{1\over\sqrt{\det(g)}}{v\over(1-{2M\over r})^{2}},{1\over\sqrt{\det(g)}}{f(v)\over{1-{2M\over r}}},0,\ldots,0\Big{)}$$

$$q(v,x)={2M\over r^{2}\Big{(}1-{2M\over r}\Big{)}^{2}}h(v),$$

$$\partial_{t}\Bigg{(}{v\over(1-{2M\over r})^{2}}\Bigg{)}+\partial_{r}\Bigg{(}{f(v)\over 1-{2M\over r}}\Bigg{)}={2M\over r^{2}\Big{(}1-{2M\over r}\Big{)}^{2}}h(v).$$

$$\Omega:=\big{\{}r>2M\big{\}}\subset\mathbb{R}^{n}.$$

$$\partial_{t}\Bigg{(}{v\over\Big{(}1-{2M\over r}\Big{)}^{2}}\Bigg{)}+\partial_{j}\Bigg{(}{x^{j}\over r\Big{(}1-{2M\over r}\Big{)}}f(v)\Bigg{)}-{(n-1)\over r\Big{(}1-{2M\over r}\Big{)}}f(v)={2M\over r^{2}\Big{(}1-{2M\over r}\Big{)}^{2}}\,h(v).$$

$$\displaystyle\partial_{t}v+\partial_{j}\Big{(}\Big{(}1-{2M\over r}\Big{)}{x^{j}\over r}f(v)\Big{)}=g(v,r),$$

$$\displaystyle g(v,r):=\partial_{j}\Big{(}\Big{(}1-{2M\over r}\Big{)}{x^{j}\over r}\Big{)}f(v)+{2M\over r^{2}}\big{(}f(v)+h(v)\big{)}$$

$$F^{\prime }(v)=f^{\prime }(v)U^{\prime }(v),v\in [-1,1].$$

$$\partial_{t}U(v)+\Big{(}1-{2M\over r}\Big{)}{x^{j}\over r}\partial_{j}F(v)\leq U^{\prime}(v){2M\over r^{2}}\Big{(}f(v)+h(v)\Big{)}.$$

$$v(0,\cdot )=v_0$$

$$\displaystyle\int_{0}^{+\infty}\int_{\Omega}\Big{(}U(v)\partial_{t}\phi+F(v)\partial_{j}\big{(}\Big{(}1-{2M\over r}\Big{)}{x^{j}\over r}\phi\big{)}+U^{\prime}(v){2M\over r^{2}}\Big{(}f(v)+h(v)\Big{)}\phi\Big{)}\,dxdt$$

$$+{\int }_{\Omega }U(v_0)\varphi (0,\cdot )dx\ge 0$$

$$\nabla_{\alpha}\big{(}T^{\alpha\beta}(\rho,u)\big{)}=0,$$

$$T^{\alpha \beta }(\rho ,u)=\rho u^{\alpha }{u}^{\beta },$$

$$-1=-\Big{(}1-{2M\over r}\Big{)}(u^{0})^{2}+\Big{(}1-{2M\over r}\Big{)}^{-1}(u^{1})^{2}.$$

$$v:=\frac{1}{(1-2M/r)}\frac{u^1}{u^0},$$

$$(u^0{)}^2=\frac{1}{(1-v^2)(1-2M/r)},(u^1{)}^2=\frac{v^2}{1-v^2}(1-2M/r).$$

$$\displaystyle\partial_{t}\Big{(}{\rho\over{1-v^{2}}}\Big{)}+\Big{(}{1-{2M\over r}}\Big{)}\partial_{r}\Big{(}{\rho v\over{1-v^{2}}}\Big{)}+\rho{{v(2r-2M)}\over{r^{2}(1-v^{2})}}$$

$$\displaystyle\partial_{t}\Big{(}{\rho v\over{1-v^{2}}}\Big{)}+\Big{(}{1-{2M\over r}}\Big{)}\partial_{r}\Big{(}{\rho v^{2}\over{1-v^{2}}}\Big{)}+\rho{{M(1-3v^{2})+2v^{2}r}\over{r^{2}(1-v^{2})}}$$

$$\partial_{t}\Big{(}{v\over(1-{2M\over r})^{2}}\Big{)}+\partial_{r}\Big{(}{{v^{2}/2}\over{1-{2M\over r}}}\Big{)}+{M\over{r^{2}(1-{2M\over r})^{2}}}=0.$$

$$\partial_{t}\Big{(}{v\over(1-2M/r)^{2}}\Big{)}+\partial_{r}\Big{(}{f(v)\over 1-2M/r}\Big{)}={2Mr^{2}\over(1-2M/r)^{2}}h(v).$$

$$f(s)=s^2/2-1/2,h(s)=0.$$

$$f(\pm 1)+h(\pm 1)=0,f^{\prime }(\pm 1)+h^{\prime }(\pm 1)\ne 0.$$

$$\frac{dt}{ds}$$

$$\frac{d{x}^j}{ds}$$

$$u(s)$$

$$u^{\prime }(s)$$

$$\displaystyle=\Bigg{(}{{d-1}\over{r({1-{2M\over r})}}}-{\partial_{j}\Big{(}{x^{j}\over{r({1-{2M\over r})}}}}\Big{)}\Bigg{)}f\big{(}u(s)\big{)}+{2M\over(r-2M)^{2}}\,h\big{(}u(s)\big{)}$$

$$\displaystyle={2M\over(r-2M)^{2}}\Big{(}f\big{(}u(s)\big{)}+h\big{(}u(s)\big{)}\Big{)}.$$

$$\underset{\Omega }{\sup }|v_0|\le 1,$$

Full text

1. 1 Introduction
2. 2 Formulation based on the Schwarzschild coordinates
3. 3 Characteristics and steady states
4. 4 Coordinates covering the black hole interior
5. 5 Finite volume method and convergence analysis

1 Introduction

We design and study a finite volume scheme for a class of nonlinear hyperbolic equations posed on a Schwarzschild black hole background. This paper is the follow-up of earlier investigations by LeFloch and co-authors [2, 10, 13, 14]. As is common in the mathematical theory of hyperbolic balance laws, we consider a (drastically) simplified version of the compressible Euler equations and we describe the fluid evolution by a single scalar unknown function, typically representing the velocity of the fluid. For relativistic problems the velocity is naturally bounded and, after normalization, we seek for solutions

[TABLE]

defined on a “spacetime” $\mathcal{M}$ —explicitly described below in a global coordinate chart222so that the present paper is aimed at a reader interested in the discretization of nonlinear hyperbolic equation with variable coefficients.— and satisfying the following hyperbolic balance law

[TABLE]

Here, $X^{\alpha}=X^{\alpha}(w,\cdot)$ is the so-called flux vector field parametrized by the real variable $w\in[-1,1]$ and defined on $\mathcal{M}$, while $q=q(w,\cdot)$ is a prescribed real-valued function. Structural conditions (even for smooth solutions, as specified later in this text) must be imposed on the vector field in order for the balance law to admit a well-posed initial value formulation.

Our objectives in this paper are as follows:

• •

Choosing $\mathcal{M}$ to be (the outer domain of communication of) a Schwarzschild black hole, we introduce a class of hyperbolic balance laws (1.2) and formulate the associated initial value problem. We then seek for weak solutions $v:\mathcal{M}\to[-1,1]$ possibly containing shock waves which must satisfy a suitable entropy condition (discussed below).

• •

Next, we design a finite volume scheme that allows us to numerically approximate these weak solutions and we derive several fundamental properties of interest: maximum principle, entropy inequalities, etc. We establish the strong convergence of this scheme toward a weak solution of the initial value problem.

Our arguments are based on a generalization of DiPerna’s theory of measure-valued solutions [6] and require us to cope with the effects of the curved black hole geometry.

An outline of this paper is as follows. In Section 2 we introduce the class of hyperbolic equations of interest and provide a motivation from pressureless fluid dynamics. In Section 3, we analyze the geometry of the curved characteristics in the black hole geometry and the class of steady state solutions which represent a fluid at rest. In Section 4, we discuss an alternative choice of slicing and which illustrate how the balance gets transform under change of coordinates. In Section 5, we introduce our finite volume scheme and state the convergence theory. The entropy inequalities satisfied by the weak solutions and their discrete version are also derived, and the proof of convergence is completed.

2 Formulation based on the Schwarzschild coordinates

2.1 The choice of coordinates

[TABLE]

Remark 2.1**.**

Passing to the so-called Eddington-Finkelstein coordinates would allow us to eliminate this singularity, but at the expense of adding further complexity in the algebraic expressions. Fortunately, the coordinates in (2.1a) are suitable for our purpose of analyzing the dynamics of a fluid outside the horizon. See Section 4 for a different choice of coordinates.

2.2 The model of interest

Choosing the vector field in the left-hand side of the balance law (1.2) to be

[TABLE]

and the source term to be

[TABLE]

we arrive at the following hyperbolic balance law:

[TABLE]

Here, the functions $f=f(w)$ and $h=h(w)$ are prescribed functions, while the unknown scalar field is $v:\mathbb{R}_{+}\times\Omega\mapsto[-1,1]$, defined for all $t\geq 0$ and $r\geq 2M$, and we work in the exterior of the ball with radius $2M$, that is

[TABLE]

In our model the unknown $v$ need not be spatially symmetric, so it convenient to rewrite (2.2) in Cartesian coordinates, i.e.

[TABLE]

Finally, in order to eliminate the singularity ${1\over 1-{2M/r}}$, we propose an equivalent form, as follows.

Definition 2.2**.**

The equation with unknown $v:\mathbb{R}_{+}\times\Omega\mapsto[-1,1]$

[TABLE]

is referred to as a hyperbolic balance law on a Schwarzschild black hole.

At this juncture, it should be emphasized that further conditions (presented in Section 3) will be required on the flux function $f$ in order for the interval $[-1,1]$ to be an invariant domain.

Definition 2.3**.**

A pair of functions $(U,F):[-1,1]\to\mathbb{R}\times\mathbb{R}$ is called a convex entropy pair for the equation (2.5) if the function $v\in[-1,1]\mapsto U(v)$ is convex and

[TABLE]

We always tacitly assume that an entropy $U$ is normalize to satisfy $U(0)=0$. Then, by definition an entropy solution to the equation (2.5) must satisfy, for all convex entropy pair $(U,F)$,

[TABLE]

We prescribe an initial data $v_{0}$ at the time $t=0$, that is,

[TABLE]

and we formalize our notion of solution as follows.

Definition 2.4**.**

Given a measurable function $v_{0}:\Omega\to[-1,1]$, a measurable function $v:\mathbb{R}_{+}\times\Omega\mapsto[-1,1]$ is called an entropy solution to the Cauchy problem (2.5) and (2.8) if the following inequalities hold

[TABLE]

for all convex entropy pairs $(U,F)$ and all compactly supported test-functions $\phi\geq 0$.

2.3 Derivation from the relativistic Euler system

Our motivation for introducing the above class of balance laws comes from a formal derivation made from the Euler equations for a relativistic compressible fluid, which read

[TABLE]

in which $\nabla$ denotes the Levi-Civita connection associated with the Schwarzschild metric (2.1a). We are interested here in the energy-momentum tensor of a pressureless fluid, given by

[TABLE]

where $\rho:M\mapsto(0,+\infty)$ denotes the density of the fluid and the velocity field $u=(u^{\alpha})$ is normalized to be future-oriented, unit and timelike $u^{\alpha}u_{\alpha}=g_{\beta\beta^{\prime}}u^{\alpha}u^{\beta^{\prime}}=-1$ with $u^{0}>0$ and, therefore,

[TABLE]

By assuming spherical symmetry, we can derive from the above system a single equation satisfied by a suitably normalized component of the velocity field, denoted below by $v\in(-1,1)$.

As usual, by taking

[TABLE]

we get

[TABLE]

Elementary computations (following [13]) yield us

[TABLE]

Combining these two equations together, we get

[TABLE]

We now compare (2.15) with (2.4). Restricting now attention to radially symmetric solutions, then (2.4) is equivalent to

[TABLE]

Clearly, this latter equation includes (2.15) as a special case, obtained by taking

[TABLE]

Hence, we can regard (2.4) as a generalization to (2.15).

3 Characteristics and steady states

3.1 Maximum principle

The method of characteristics allows us to obtain a first insight about the properties of (sufficiently regular solutions) to our balance law (2.4). It leads to ordinary differential equations along characteristic curves parametrized with respect to some time parameter (denoted by $s$ below). We would like to deduce some properties of solution $v$ by proposing the following assumption on the flux $f$ and the source $h$.

Assumption 1**.**

The flux and source functions are assumed to satisfy

[TABLE]

We motivate our condition by the following analysis along characteristic curves. So, we consider the coupled system

[TABLE]

A straightforward computation shows that

[TABLE]

This equation tells us how the values of a solution evolve along characteristics, and we use it in order to establish a maximum principle. It is convenient to assume a strict inequality in the data.

Proposition 3.1** (Maximum principle).**

Consider the balance law (2.5) under the condition (3.1). Then, given any initial data (2.8) satisfying

[TABLE]

the solution $v=v(t,x)$ satisfies the same bound for all times

[TABLE]

as long as they remain sufficiently regular.

Proof.

Observe first that if $v_{0}=\pm 1$ initially then it remains so for all times. It is sufficient to show that $v\leq 1$, since exactly the same arguments apply to showing $v\geq-1$.

Consider first the case of $H^{\prime}(1)>0$ with $H(s):=f(s)+h(s)$. By continuity, we have $H(s)<0$ for all $s\in(1-\epsilon,1)$ and some $\epsilon$. Hence, if $u\in[1-\epsilon,1)$, (3.3) implies that $u^{\prime}<0$ and, consequently, $u\leq 1$ for all times.

In the case $H^{\prime}(1)<0$, we have $H(s)<0$ for all $s\in(1,1+\epsilon)$ and some $\epsilon$. Recall $\sup|v_{0}|<1$, if $u(s_{2})>1$ for some $s_{2}$, then $u(s)>1$ for $s\in(s_{1},s_{2})$ with $s_{1}:=\sup\{s:u(s)\leq 1,s<s_{2}\}$. We see that $u(s_{1})=1$. However, if we integrate (3.3) in $[s_{1},s_{2}]$, the right-hand side would be negative ,while the left-hand side would be positive. Hence, $u\leq 1$ (actually, if $u$ could reach 1 at some $s$, then it must remain identically $1$ afterwards). ∎

3.2 Geometry of the characteristic curves

[TABLE]

and, more explicitly,

[TABLE]

where $r_{0}=r(s_{0})$ and $u_{0}=u(s_{0})$ are given data at some time $s_{0}$.

Concerning the global behavior of the characteristics in the special case $f(w)=w^{2}/2-1/2$ and $h(w)=0$, which is the Burgers equation posed on the Schwarzschild background, the weak solutions in the $(u,r)$ plane can be expressed in terms of the initial data via a minimisation formulation based on characteristics; see [2].

Here, to proceed with the study of the characteristic curves and for the sake of definitness, we assume some specific signs about the functions $f$ and $h$.

Assumption 2**.**

[TABLE]

and

[TABLE]

A direct consequence from (3.6) and (3.7) is that

[TABLE]

We rewrite (3.5) as

[TABLE]

and, by solving for $u$, the ordinary differential equation (3.4a) for the radius function $r(s)$ can be written as

[TABLE]

Here, $\widehat{F}_{+}^{(-1)}$ and $\widehat{F}_{-}^{(-1)}$ are the inverse functions of $\widehat{F}_{+}$ and $\widehat{F}_{-}$, respectively, and

[TABLE]

Note that $\widehat{F}_{\pm}$ are single-valued functions within the domain of interest.

We follow [2] and introduce the escape velocity (whenever it exists)

[TABLE]

which satisfies the property

[TABLE]

Replacing the radius $r_{0}$ by the escape velocity parameter $u_{0}^{E}$ in (3.12), we obtain

[TABLE]

[TABLE]

The late-time behavior of $u=u(s)$ can be checked to be described as follows:

• •

Negative initial data. The function $u(s)$ decreases (as follows from (3.3) and the assumption (3.6)). If $u_{0}\in(-1,0]$ with initial data $(s_{0},r_{0})$, then $u(s)$ remains negative and decreasing and $r(s)$ decreases towards $2M$. More precisely, we have

[TABLE]

• •

Positive initial data with $0<u_{0}<u_{0}^{E}$. The positivity of $u_{0}$ initially ensures $dr/ds>0$, that is, the characteristic curve initially moves away from the black hole. However, $u(s)$ keeps decreasing and eventually reaches [math] at some time $s_{0}$. The dynamics then coincides with that for negative initial data. We conclude that

[TABLE]

• •

Positive initial data with $u_{0}\geq u_{0}^{E}$. In this case, the characteristic curve moves away from the black hole for all times and the asymptotic behavior is

[TABLE]

3.3 Steady state solutions

Finally, let us consider solutions that are steady states representing a fluid at rest in the curved black hole geometry. This is a special class of solutions of interest, for instance, in designing (well-balanced) numerical schemes and in finding test cases. In view of the radially symmetric form of our equation (2.16) (but possibly for non-radially symmetric solutions), for a time-independent solution we obtain the ordinary differential equation

[TABLE]

and, once again, we get the same ordinary differential equation as (3.4b)

[TABLE]

Under Assumptions (1) and (2), for any given data $(r_{0},u_{0})$ we can distinguish between two cases:

• •

Negative $u_{0}$. Then $u$ is increasing and

[TABLE]

• •

Positive $u_{0}$. Then $u$ is decreasing and

[TABLE]

4 Coordinates covering the black hole interior

4.1 An alternative choice of time slicing

In this section, we illustrate the fact that coordinates can be chosen in many different manners. While, for Schwarschild spacetime, this leads to significantly more involved algebraic expressions, such alternative coordinates may allow one to cover a larger region of the spacetime. For definiteness, in this section we take $n=3$. Hence, we now introduce a nonlinear hyperbolic equation posed in a larger domain of the Schwarzschild geometry, obtained by “crossing” the horizon and we study the interior of the black hole. We follow [5] and introduce the following metric:

[TABLE]

Here $R_{0}\in(0,M]$ is a parameter that is fixed, and we observe that the above expression is identical to the metric (2.1a) in the limit $R_{0}\to 0$, for which the radial variables $R$ and $r$ would then coincide. This new slicing (4.1a) allows us to go inside of the black hole (when $R_{0}>0$), and we cover the region $\{r:r+R_{0}-2M>0\}$, within which the metric remains of a definite Lorenztian signature.

In fact, we can transform (4.1a) (for a restricted domain of the variables, only) into the metric (2.1a), by setting

[TABLE]

In the following, it will be convenient to rely on the vector fields

[TABLE]

We rewrite (4.1a) in the matrix form

[TABLE]

with inverse

[TABLE]

After a tedious computation, the (non-vanishing) Christoffel symbols $\Gamma^{\mu}_{\alpha\beta}={1\over 2}\widehat{g}^{\mu\nu}(\widehat{\partial}_{\alpha}\widehat{g}_{\beta\nu}+\widehat{\partial}_{\beta}\widehat{g}_{\alpha\nu}-\widehat{\partial}_{\nu}\widehat{g}_{\alpha\beta})$ are found to be

[TABLE]

4.2 Formulation of the balance law

We follow the strategy in the previous section and derive our equation from the pressureless Euler system. For the (normalized) vector

[TABLE]

we find

[TABLE]

Proposition 4.1**.**

From the pressureless Euler equations, the velocity vector

[TABLE]

satisfies the nonlinear hyperbolic equation

[TABLE]

Proof.

In view of the notation (4.5), we get

[TABLE]

which is one representation of $(\widehat{u}^{0},\widehat{u}^{1})$. Plugging these expressions into (2.11), we obtain

[TABLE]

From (2.10), thus $\widehat{\partial}_{\alpha}T^{\alpha\beta}+\Gamma^{\alpha}_{\alpha\gamma}T^{\gamma\beta}+\Gamma^{\beta}_{\alpha\gamma}T^{\alpha\gamma}=0$, we get

[TABLE]

[TABLE]

Now we set

[TABLE]

and our calculation leads us to

[TABLE]

Finally, further cumbersome calculations give usthe final form (4.6). ∎

4.3 Characteristics and maximum principle

As we mentioned in the beginning of this section, the new metric $\widehat{g}$ coincides with the Schwarzchild metric $g$ when $r$ is replaced by $R=r+R_{0}$ in (2.1a), hence it is not surprising to have the following result. Namely, if we replace $r$ by $R$ throughout Section 2, then Burgers equation (2.15) is equivalent to (4.6). This is easy to check with

[TABLE]

From the equation (2.15) we have

[TABLE]

which coincides with (4.6). We can now restate our previous results in the new coordinates, and we only discuss in detail the new features.

From our equation (4.6), the characteristic curves with $\widehat{u}(s):=\widehat{v}(t(s),r(s))$ are given by

[TABLE]

Moreover, we have

[TABLE]

Similarly to what we did in Proposition 3.1, we can check the following result.

Proposition 4.2** (Maximum principle).**

Consider the equation (4.6). If the initial data satisfies

[TABLE]

then any smooth solution to (4.6) also satisfies

[TABLE]

Now, in the $(\widehat{u},r)$–plane let us observe that

[TABLE]

which follows from (4.13) and (4.14). From this, we obtain

[TABLE]

where $(\widehat{u}_{0},R_{0})$ is the initial location. Thus, our previous conclusions concerning the characteristic curves are recovered here. This is of course not surprising, since we are treating the same differential equation expressed in different coordinates. This second formulation however may have some numerical advantage when the horizon, instead of being fixed as it is in the present model, is dynamical.

5 Finite volume method and convergence analysis

5.1 Formulation of the finite volume scheme

Having considered the formulation (2.5) (in Section 2) and the formulation (4.6) (in Section 4), we now study the numerical approximation of the general balance law (1.2) in a setting that, in principle, may encompass both formulations. For definiteness, we treat the outer domain of communication so that the space variable varies in a half-line and no boundary condition is required at the boundary. The hyperbolic model of interest reads

[TABLE]

in which $\mathcal{M}$ denotes the outer domain of communication of a Schwarzschild black hole with radius $2M$, as we described earlier. Here $X^{\alpha}=X^{\alpha}(w,\cdot)$ is a smooth vector field on $\mathcal{M}$, depending upon the real variable $w\in[-1,1]$. An hyperbolicity condition and a condition at the boundary will be made explicit below.

We are going to formulate a finite volume scheme for the equation (5.1) and establish its convergence by generalizing the technique of proof in [1]. In contrast with this later work, the spacelike slices in $\mathcal{M}$ are non-compact and the flux vector $X^{\alpha}(\bar{v},x)$ is no longer assumed to be geometry compatible (i.e. $X^{\alpha}(\bar{v},x)$ does not satisfy divergence free condition), and at the (horizon) boundary, no boundary data is needed.

In the class of interest in the present paper, the flux vector field satisfies the following property, which implies that no boundary condition is needed: the spatial components of the vector field $X(\cdot,x)$ vanishes on the boundary, i.e.

[TABLE]

Following [1] we design a finite volume scheme for (5.1) as follows. We introduce a spacetime triangulation $\mathcal{T}_{h}=\bigcup_{K\in\mathcal{T}_{h}}K$ of $\mathcal{M}$ such that the boundary $\partial K$ of each element $K$ is the union of three possible types of faces:

• •

A face $e_{K}^{+}$ is spacelike and we denote its future–oriented outward unit normal by $n_{K,e_{K}^{+}}$.

• •

A face $e_{K}^{-}$ (in the past of $e_{K}^{+}$) is also spacelike, and we denote its past–oriented outward unit normal by $n_{K,e_{K}^{-}}$.

• •

A vertical face denoted by $e^{0}$ is timelike, and whose inward unit normal is denoted by $n_{K,e^{0}}$, and the union of all of such faces is denoted by $\partial^{0}K:=\partial K\setminus\{e_{K}^{+},e_{K}^{-}\}$.

By definition, for every pair of distinct elements $K,K^{\prime}\in\mathcal{T}_{h}$, the intersection $K\cap K^{\prime}$ is either a common face of $K,K^{\prime}$ or a submanifold of co-dimension at least $2$. We use $K_{e}$ to denote the unique neighbor of $K$ sharing the same edge $e$. We denote by $K^{\pm}$ the neighbors of $K$ which share the same edge $e^{\pm}$.

By integrating the equation (5.1) over an arbitrary element $K\in\mathcal{T}_{h}$ and applying the divergence theorem, we obtain

[TABLE]

Here, $n$ denotes the exterior and unit, normal vector field along the boundary face under consideration, while $d\overline{V}$ is the induced measure element on the boundary. Our finite volume scheme is based on the following approximation formulas, in which $e^{0}$, etc. denotes an edge of $K$:

[TABLE]

• •

Discretization of the main variable:

[TABLE]

• •

Discretization of the flux:

[TABLE]

• •

Discretization of the source term:

[TABLE]

Here, the numerical flux $f_{K,e}:\mathbb{R}^{2}\to\mathbb{R}$ is chosen to satisfy the properties of consistency, conservation and monotonicity:

[TABLE]

• •

Consistency property:

[TABLE]

• •

Conservation property:

[TABLE]

• •

Monotonicity property:

[TABLE]

Also we have written

[TABLE]

and

[TABLE]

Finally, the finite volume approximations are defined by

[TABLE]

5.2 Convergence and existence theory

Based on the geometric formulation of a finite volume method above, we can now proceed with the analysis of our model problem (2.5). We integrate (2.5) over an element $K$ and by applying the divergence theorem

[TABLE]

where $n$ denotes the outward unit normal vector. We apply the following approximations:

[TABLE]

• •

Discretization of the main variable:

[TABLE]

• •

Discretization of the flux:

[TABLE]

where, with $x_{e}$ being the center of $e$ and $r_{e}=|x_{e}|$,

[TABLE]

• •

Discretization of the source term:

[TABLE]

where, $r_{K}$ being the radial variable evaluated at the center of $K$,

[TABLE]

In the above, we denoted by $e$ some edge of $K$, and the numerical flux $f_{K,e}:~{}\mathbb{R}^{2}\to\mathbb{R}$ is chosen to satisfy the properties of consistency, conservation and monotonicity, that is, in our case

[TABLE]

• •

Consistency property:

[TABLE]

• •

Conservation property:

[TABLE]

• •

Monotonicity property:

[TABLE]

The finite volume approximations are then given by the explicit scheme

[TABLE]

For the sake of stability, we impose the CFL stability condition

[TABLE]

as well as the source stability condition

[TABLE]

Now we are ready to state our convergence result.

Theorem 5.1**.**

Consider the Cauchy problem for the balance law (2.5) posed on the domain $\Omega$ under the assumption (1) and (2). Impose the initial condition (2.8) with $v:[0,+\infty)\times\Omega\to[-1,1]$ in $L^{1}(\Omega)$. Let $\mathcal{T}_{h}$ be a triangulation and $\tau=\tau(h)$ be the time increment, satisfying

[TABLE]

Let $f_{K,e}$ be a family of numerical flux satisfying the consistency, conservation and monotonicity conditions in (5.11a)–(5.11c) and satisfies the CFL condition (5.13) and the stability condition (5.14). Then the discrete scheme (5.12) uniquely defines the family of approximate solution $v^{n}_{K}$. By defining a piecewise constant function $v^{h}:\mathbb{R}_{+}\times\Omega\to\mathbb{R}$ by

[TABLE]

then the sequence $v_{h}:[0,+\infty)\times\Omega\to[-1,1]$ is uniformly bounded in L_{\text{loc}}^{\infty}\big{(}[0,+\infty),L^{1}(\Omega)\big{)} and converges almost everywhere to an entropy solution $v:[0,+\infty)\times\Omega\to[-1,1]$ (in the sense of Definition 2.4) v\in L_{\text{loc}}^{\infty}\big{(}[0,+\infty),L^{1}(\Omega)\big{)}.

The above theorem implies the existence and stability of weak solutions for our model.

Corollary 5.2**.**

Consider the Cauchy problem for the balance law (2.5) posed on the domain $\Omega$ under the assumption (1) and (2). Impose the initial condition (2.8) with $v_{0}\in L^{1}(\Omega)$ and $\|v_{0}\|_{L^{\infty}(\Omega)}\leq 1$. Then there exists an entropy solution $v:[0,+\infty)\times\Omega\to[-1,1]$ in L_{\text{loc}}^{\infty}\big{(}[0,+\infty),L^{1}(\Omega)\big{)} to this problem.

5.3 Discrete entropy inequalities

Entropy inequalities play a key role in the proof of Theorem 5.1.

Proposition 5.3**.**

Under the assumptions in Theorem 5.1, the finite volume approximations satisfy the discrete maximum principle:

[TABLE]

Proof.

We first assume $\max_{\mathcal{T}_{h}}|v^{n}_{K}|\leq 1$ for all elements $K\in\mathcal{T}_{h}$. We rewrite (5.12) as

[TABLE]

which with an obvious notation we express in the form

[TABLE]

We have observed here that $A^{n}_{K}+\sum_{e\in\partial K}A^{n}_{K,e}=1$.

The monotonicity of $f_{K,e}$ implies

[TABLE]

while the CFL condition (5.13) gives us

[TABLE]

Therefore, we have $A^{n}_{K}\geq{1\over 2}$.

On the other hand, since $f(1)+h(1)=0$ and $v^{n}_{K_{e}}\leq 1$ we have

[TABLE]

where we used the source stability condition (5.14). Similarly, we find

[TABLE]

Now, we state a convex decomposition of $v^{n+1}_{K}$, which plays an important role in deriving the discrete entropy inequalities given below. For each $K$ and $e$, we define

[TABLE]

The following lemma provides a standard result concerning the existence of discrete entropy flux terms and an entropy inequality relating $\widetilde{v}^{n+1}_{K,e}$ and $v^{n+1}_{K}$. We omit the proof and refer to [1] and the references therein.

Lemma 5.4** (Discrete entropy inequalities).**

Let $(U,F)$ be a convex entropy pair. Then there exists a family of discrete entropy flux functions $F_{K,e}:\mathbb{R}^{2}\to\mathbb{R}$ satisfying the following conditions:

[TABLE]

• •

Consistency with the entropy flux $F$:

[TABLE]

• •

Conservation property:

[TABLE]

• •

Discrete entropy inequality:

[TABLE]

Equivalently, (5.20c) can be written in terms of $v^{n+1}_{K,e}$ and $v^{n}_{K}$ as

[TABLE]

with $R^{n+1}_{K,e}:=U(v^{n+1}_{K,e})-U(\widetilde{v}^{n+1}_{K,e})$. The entropy dissipation estimate below will serve to establish the convergence result.

Proposition 5.5** (Discrete entropy balance law).**

Let $U:\mathbb{R}\to\mathbb{R}$ be a strictly convex function and set $\alpha:=\inf_{v\in[-1,1]}U^{\prime\prime}(v)$. Then for all $n$ one has

[TABLE]

Proof.

By Lemma 3.5 in [4] and the convex decomposition identity (5.19c), we have

[TABLE]

Next, we multiply by $|e||K|/p_{K}$ in (5.21), and sum up over all $K$ and $e$,

[TABLE]

This leads us to (5.22). ∎

For the proofs of the following lemmas, see [1] and the references therein for details. The local entropy inequalities read as follows.

Lemma 5.6**.**

One has

[TABLE]

The global entropy inequalities read as follows.

Lemma 5.7**.**

Let $(U,F)$ be a convex entropy pair and let $\phi=\phi(t,x)\in C_{c}([0,T)\times\Omega)$ be a test function. For each element $K$ and each face $e\in\partial K$, set

[TABLE]

Then one has

[TABLE]

5.4 Measure-valued solutions and strong convergence

We are now in a position to complete our proof of Theorem 5.1. Based on the entropy inequalities we have established, we are able to pass the limit in the inequality (5.24) as $h\to 0$. Then we associate with a subsequence of $v_{h}$ (which is uniformly bounded in $[0,T)\times\Omega\to\mathbb{R}$ for fixed $T$) a Young measure $\nu:[0,T)\times\Omega\to\text{Prob}(\mathbb{R})$, which is a family of probability measures in $\mathbb{R}$ parametrized by $(t,x)\in[0,T)\times\Omega$. We then show that the Young measure, describing all the weak-star limits of $v_{h}$, is an entropy measure-valued solution in the sense of DiPerna. The strong convergence result follows from the DiPerna’s uniqueness theorem, see [6].

The Young measure allows us to write, for every continuous function $a:\mathbb{R}\to\mathbb{R}$,

[TABLE]

in the $L^{\infty}$ weak-star topology. As presented in [1], it suffices to show that $\nu$ is an entropy measure-valued solution to our balance law, in order to imply that $\nu_{t,x}$ reduce to a Dirac mass $\delta_{v(t,x)}$ if this is true at the initial time $t=0$. The convergence in (5.25) then holds in a strong sense and $v_{h}$ converges to the entropy solution $v$ to the Cauchy problem.

Lemma 5.8**.**

Let $\nu:[0,T)\times\Omega\to\text{Prob}(\mathbb{R})$ be the Young measure associated with the sequence $v_{h}$. Then for every convex entropy pair $(U,F)$ one has

[TABLE]

for all non-negative test functions $\phi:[0,T)\times\Omega\to\mathbb{R}_{+}$.

For all convex entropy pairs, we thus have

[TABLE]

and the proof of of Theorem 5.1 is completed.

Acknowledgments

The authors were supported by the Innovative Training Network (ITN) entitled ModCompShock under the grant 642768.