Geometric Error of Finite Volume Schemes for Conservation Laws on Evolving Surfaces

# Geometric Error of Finite Volume Schemes for Conservation Laws on Evolving Surfaces

Jan Giesselmann and Thomas Müller J. Giesselmann Weierstrass Institute,
Mohrenstr. 39, D-10117 Berlin, Germany
22email: giesselm@wias-berlin.deT. Müller Abteilung für Angewandte Mathematik, Universität Freiburg,
Hermann-Herder-Str. 10, D-79104 Freiburg, Germany
44email: mueller@mathematik.uni-freiburg.de
###### Abstract

This paper studies finite volume schemes for scalar hyperbolic conservation laws on evolving hypersurfaces of . We compare theoretical schemes assuming knowledge of all geometric quantities to (practical) schemes defined on moving polyhedra approximating the surface. For the former schemes error estimates have already been proven, but the implementation of such schemes is not feasible for complex geometries. The latter schemes, in contrast, only require (easily) computable geometric quantities and are thus more useful for actual computations. We prove that the difference between approximate solutions defined by the respective families of schemes is of the order of the mesh width. In particular, the practical scheme converges to the entropy solution with the same rate as the theoretical one. Numerical experiments show that the proven order of convergence is optimal.

###### Keywords:
hyperbolic conservation laws finite volume schemes curved surfaces error bound
###### Msc:
65M08 35L65 58J45

## 1 Introduction

Hyperbolic conservation laws serve as models for a wide variety of applications in continuum dynamics. In many applications the physical domains of these problems are stationary or moving hypersurfaces. Examples of the former are in particular geophysical problems WDHJS92 and magnetohydrodynamics in the tachocline of the sun Gil00; SBG01. Examples of the latter include transport processes on cell surfaces RS05, surfactant flow on interfaces in multiphase flow BPS05 and petrol flow on a time dependent water surface. There are several recent approaches to the numerical computation of such equations. Numerical schemes for the shallow water equations on a rotating sphere can be found in CHL06; Gir06; Ros06. For the simulation of surfactant flow on interfaces we refer to AB09; BS10; JL04. As we are interested in numerical analysis we focus on nonlinear scalar conservation laws as a model for these systems. The intense study of conservation laws posed on fixed Riemannian manifolds started within the last years. There are results on well-posedness BL06; DKM13; LM13 of the differential equations and on the convergence of appropriate finite volume schemes ABL05; Gie09; GW12; LON09. For recent developments on finite volume schemes for parabolic equations we refer to LNR11.

In the previous error analysis for finite volume schemes approximating nonlinear conservation laws on manifolds the schemes were defined on curved elements lying on the curved surface and it was assumed that geometric quantities like lengths, areas and conormals are known exactly. While this is a reasonable assumption for schemes defined on general Riemannian manifolds or even more general structures ALN11; LO08.2 with no ambient space, most engineering applications involve equations on hypersurfaces of and one aims at computing the geometry with the least effort. This is in particular important for moving surfaces where the geometric quantities have to be computed in each time step. Now the question arises to which extent an approximation of the geometry influences the order of convergence of the scheme. We will treat this question in the “embedded“ case where an explicit embedding of the surface under consideration into Euclidean space is known. It is an interesting question for future studies whether our analysis can be extended to errors arising from the discretisation of the geometry of the underlying space in an ”invariant” description, like the very general one analysed in ALN11; LO08.2.

We consider the following initial value problem, posed on a family of closed, smooth hypersurfaces . For a derivation cf. DE07.2; DKM13; Sto90. For some , find with

 ˙u+u∇Γ⋅v+∇Γ⋅f(u,⋅,⋅) =0 in GT, (1) u(⋅,0) =u0 on Γ(0), (2)

where is the velocity of the material points of the surface and are initial data. For every the flux is a smooth vector field tangential to which depends Lipschitz on and smoothly on . Moreover, we impose the following growth condition

 |∇Γ⋅f(¯u,x,t)|≤c+c|¯u|∀¯u∈R,(x,t)∈GT (3)

for some constant . By we denote the material derivative of which is given by

 ˙u(Φt(x),t):=ddtu(Φt(x),t),

where is a family of diffeomorphisms depending smoothly on , such that is the identity on Obviously this excludes changes of the topology of We will assume that the movement of the surface and also the family is prescribed. A main result of this paper is a bound for the difference between two approximations of . In particular, we will give an estimate for the difference between the flat approximate and the curved approximate solution. By curved approximate solution we refer to a numerical solution given by a finite volume scheme defined on the curved surface, cf. Section 2.2, and by flat approximate solution we refer to a numerical solution given by a finite volume scheme defined on a polyhedron approximating the surface, cf. Section 2.3. We will see that the arising geometry errors can be neglected compared to the error between the curved approximate solution and the exact solution, i.e. both approximate solutions converge to the entropy solution with the same convergence rate. We will present numerical examples showing that the proven convergence rate is optimal under the assumptions for the numerical analysis. However, for most numerical experiments we observe higher orders of convergence.

Our analysis also indicates that the geometry error poses an obstacle to the construction of higher order schemes. To this end we perform numerical experiments underlining in which manner the order of convergence of the higher order scheme is restricted by the approximation of the geometry. This shows that to obtain higher order convergence also the geometry of the manifold has to be approximated more accurately, cf. Dem09 in a finite element context.

The outline of this paper is as follows. In Section 2 we review the definition of finite volume schemes on moving curved surfaces and define finite volume schemes on moving polyhedra approximating the surfaces. The approximation errors for geometric quantities are established in Section 3. Section 4 is devoted to estimating the difference between the curved and the flat approximate solution. Finally, numerical experiments are given in Section LABEL:sec:numerics.

## 2 The Finite Volume Schemes

This section is devoted to the construction of a family of triangulations of the surfaces suitably linked to polyhedral approximations of the surfaces. Afterwards we will recall the definition of a finite volume scheme on which was considered in the hitherto error analysis and define a finite volume scheme on which is an algorithm only relying on easily computable quantities. We would like to point out that our finite volume scheme is applicable to closed, smooth hypersurfaces of arbitrary geometry, which additionally may evolve in time. In the case of simpler geometries of special interest, e.g. considering conservation laws on , one can make use of the additional structures. In BFL09 for instance, studying scalar conservation laws on , the special structure of this setting is exploited by considering a longitude-latitude grid which allows the exact computation of cell areas and edge lengths and the design of a Godunov-type finite volume scheme based on the dimension-wise solution of one-dimensional Riemann problems for a class of analytic flux functions. Additionally, the finite volume scheme in BFL09 is geometry-compatible in the sense that for a divergence-free flux the numerical fluxes are (discretely) divergence-free, as well. Other approaches employing the special knowledge available for the sphere include logically rectangular grids developed in CHL06 and grids in which all edges are geodesic arcs in Gir06.

We mention that our triangulation as well as the definition of the finite volume scheme on is in the same spirit as the one from LNR11 which was developped for the diffusion equation on evolving surfaces.

### 2.1 Triangulation

We start by mentioning that there are neighbourhoods of such that for every there is a unique point such that

 x=a(x,t)+d(x,t)νΓ(t)(a(x,t)), (4)

where denotes the signed distance function to and the unit normal vector to pointing towards the non-compact component of . See DE07 for example.

Let us choose a polyhedral surface which consists of flat triangles such that the vertices of lie on and is the length of the longest edge of In addition we impose that the restriction of is one-to-one. We define as the polyhedral surface that is constructed by moving the vertices of via the diffeomorphism and connecting them with straight lines such that all triangulations share the same grid topology. A triangulation of is automatically given by the decomposition into closed faces. We define the triangulation on as the image of under We will denote the closed curved cells with and the closed curved faces with . A flat quantity corresponding to some curved quantity is denoted by the same letter and a bar, e.g. let be a curved face then In order to reflect the fact that all triangulations share the same grid topology we introduce the following notation. We denote by the family of all (closed) curved triangles relating to the same (closed) triangle on We do the same for Analogously by we denote the family of such families of triangles

For later use we state the following Lemma summarizing geometric properties, whose derivation can be found in DE07.

###### Lemma 1

Let be a polyhedral approximation of as described above then there exists such that for all

1. .

We will use the following notation. By we denote the diameter of each cell, furthermore and are the Hausdorff measures of and the boundary of respectively. When we write we mean to be a face of .

We need to impose the following assumption uniformly on all triangulations There is a constant number such that for each flat cell we have

 αh2≤∣∣¯K(t)∣∣,α∣∣∂¯K(t)∣∣≤h. (5)

Later on, we will see that (5) implies the respective estimate for the curved triangulation, cf. Remark 1. A consequence of (5) is that is a lower bound of the radius of the inner circle of , which implies that the sizes of the angles in are bounded from below. Furthermore we denote by the supremum of the spectral norm of i.e. is a bound on the eigencurvatures. By straightforward continuity and compactness arguments is uniformly bounded in space and time.

### 2.2 The Finite Volume Scheme on Curved Elements

In this section we will briefly review the notion of finite volume schemes on moving curved surfaces. We consider a sequence of times and set Moreover we assign to each and the term approximating the mean value of on and to each and face a numerical flux function , which should approximate

 \fintIn\finte(t)⟨f(u(x,t),x,t),μK(t),e(t)(x)⟩de(t)dt, (6)

where is the line element, is the unit conormal to pointing outwards from and is the standard Euclidean inner product. Please note that is tangential to Then the finite volume scheme is given by

 u0K:=\fintK(0)u0(x)dΓ(0),un+1K:=|K(tn)||K(tn+1)|unK−|In||K(tn+1)|∑e⊂∂K|e(tn)|fnK,e(unK,unKe),uh(x,t):=unK for t∈[tn,tn+1), x∈K(t), (7)

where denotes the cell sharing face with and is the surface element. As usual in corresponding convergence analysis Gie09; LON09 we assume that the used numerical fluxes are consistent, i.e.

 |e(tn)|fnK,e(u,u)=\fintIn∫e(t)⟨f(u,x,t),μK(t),e(t)(x)⟩de(t)dt∀u∈R, (8)

conservative, i.e.

 fnK,e(u,v)=−fnKe,e(v,u)∀u,v∈R, (9)

monotone, i.e.

 ddufnK,e(u,v)≥0,ddvfnK,e(u,v)≤0∀u,v∈R, (10)

and uniformly Lipschitz continuous. Let denote the Lipschitz constant of the numerical fluxes, then additionally the CFL condition

 tn+1−tn≤α2h8L (11)

has to be imposed to ensure stability, i.e. Lemmas 7 and 10. As an example for numerical flux functions satisfying these conditions we introduce Lax-Friedrichs fluxes

 LFfnK,e(u,v):=\fintIn12|e(tn)|∫e(t)⟨f(u,x,t)+f(v,x,t),μK(t),e(t)(x)⟩de(t)dt+λn(u−v), (12)

where is an artificial viscosity coefficient ensuring the monotonicity of and stabilizing the scheme.

### 2.3 The Finite Volume Scheme on Flat Elements

In this section we define finite volume schemes on which are in the same spirit as (7) but only rely on easily accessible geometrical information. We want to point out that the calculation of areas and lengths is straightforward for flat elements. As well, the approximation of integrals can be achieved using quadrature formulas by mapping cells and edges to a standard triangle and the unit interval, respectively, using affine linear maps. In this fashion we obtain for every time quadrature operators and of order respectively. In addition for any compact interval the term denotes a quadrature operator of order For Lipschitz continuous numerical flux functions we define the finite volume scheme on flat elements according to

 ¯u0¯K:=1|¯K(0)|Q¯K(0)(u0(a(⋅,0))),¯un+1¯K:=|¯K(tn)||¯K(tn+1)|¯un¯K−|In||¯K(tn+1)|∑¯e⊂∂¯K|¯e(tn)|¯fn¯K,¯e(¯un¯K,¯un¯K¯e),¯uh(x,t):=¯un¯K,% for t∈[tn,tn+1), x∈K(t). (13)

Note that by (13) the function is defined on . For the numerical analysis we need to impose the following estimate for the (geometric) error between the numerical fluxes and :

 ∣∣fnK,e(u,v)−¯fn¯K,¯e(u,v)∣∣≤Ch2∀ (u,v)∈K, K∈Th, e⊂∂K, (14)

where is a compact subset of and a constant depending only on and .

As an example for easily computable numerical flux functions for the flat scheme we define Lax-Friedrichs flux functions below. We will see in Lemma 5 that assumption (14) is valid for them. Before we can use the quadrature operators to define the numerical fluxes we need to determine the "discrete" conormals. To each flat triangle we fix a unit normal by imposing

 ⟨¯ν¯K(t),νΓ(t)(y)⟩>0, (15)

where is the barycentre of We will see in Lemma 2 that converges to for To each face and adjacent cell there is a unique unit tangent vector such that is a conormal to pointing outward from Hence this vector product is one candidate for However in general

 ¯ν¯K(t)×¯t¯K(t),¯e(t)≠±(¯ν¯K¯e(t)×¯t¯K¯e(t),¯e(t)) (16)

such that a choice like

 ¯μ¯K(t),¯e(t)=¯ν¯K(t)×¯t¯K(t),¯e(t)

would lead to a loss of conservativity of the resulting numerical fluxes. Therefore we choose

 ¯μ¯K(t),¯e(t):=12(¯ν¯K(t)×¯t¯K(t),¯e(t)+¯ν¯K¯e(t)×¯t¯K(t),¯e(t)).

We define numerical Lax-Friedrichs fluxes by

 LF¯fn¯K,¯e(u,v):=1|In|QIn[12|¯e(tn)|Q¯e(⋅)(⟨f(u,⋅,⋅)+f(v,⋅,⋅),¯μ¯K(⋅),¯e(⋅)⟩)]+λ(u−v) (17)

for some sufficiently large and being smoothly extended from to the whole of . Note, that here and in the following the quadrature operator is applied to the time dependence while is applied to the space dependence. In particular, for each time quadrature point , specified by , the space quadrature points lie on the (moving) flat edge .

## 3 Geometrical Estimates

In this section we derive estimates for the approximation errors of the geometric quantities. Throughout this section we suppress the time dependence of all quantities. All the estimates can be derived uniformly in time. To obtain the geometrical estimates, we introduce the following lift operator.

###### Definition 1

Let and a function on then we define a function on as

 ¯gl=¯g∘a|−1Γh.

Similarly we define the inverse of this lift operator by

 g−l=g∘a|Γh

for a function defined on some .

We begin our investigation with the differences between the normal vectors of the flat and curved elements.

###### Lemma 2

There is a constant such that for all flat cells and every we have

 ∥∥ν−lΓ(y)−¯ν¯K∥∥ ≤Ch. (18)

The constant depends on derivatives of in particular on

###### Proof

Without loss of generality we can assume that is a subset of such that is one of its faces and for some . We start by showing that there exists some constant such that

 |(νΓ)i|≤Ch,for i=1,2. (19)

We recall that where is the signed distance function to As the vertices of lie on we know that, if we denote the third vertex by it holds

 d(0,0,0)=0, d(h¯e,0,0)=0, d(x,y,0)=0.

Hence, the directional derivatives of with respect to and need to vanish somewhere in Thus their absolute value is of order on Due to the angle condition (5) an analogous inequality also holds for the directional derivative of with respect to . As the directional derivative of with respect to coincides with respectively, this proves (19). This immediately implies By assumption everywhere and by (15) we have which proves (18). ∎

###### Lemma 3

For the difference between the length of a curved edge and the corresponding flat edge we have

 ∣∣∣|e||¯e|−1∣∣∣≤Ch2, (20)

and for the difference between the area of a curved cell and the corresponding flat cell we have

 ∣∣ ∣∣|K|∣∣¯K∣∣−1∣∣ ∣∣≤Ch2, (21)

where does not depend on but on

Furthermore let be the parametrization of over given by then we have

 |∥c′e(s)∥−1|≤Ch2. (22)
###### Proof

We assume without loss of generality that . For small enough we can parametrize the curved cell according to (4) by a parametrization with

 c(x1,x2)=(x1,x2,0)−d(x1,x2,0)νΓ(c(x1,x2)),

where we suppressed the third coordinate in . The ratio of volume elements of and with respect to the parametrization is given by

 √|g|:=√det(g),

where the matrix is defined by

 g=(gij)1≤i,j≤2:=(⟨∂ic,∂jc⟩)1≤i,j≤2.

For the parametrization of we have

 ∂ic=ei−⟨∇d,ei⟩νΓ∘c−d∂ic(∇νΓ)T∘cfor i=1,2,

where denotes the -th standard unit vector. Due to the bounded curvature of and Lemma 1 we can show that

 ∂ic=ei−((νΓ)iνΓ)∘c+O(h2)for i=1,2. (23)

Applying (18) we see that

 νΓ=±(0,0,1)+O(h) and ⟨ei,νΓ⟩=(νΓ)i=O(h)for i=1,2.

Thus, for the matrix we have

 g=(1+O(h2)O(h2)O(h2)1+O(h2))

which implies for the volume element

 dK=√|g|d¯K=√1+O(h2)d¯K=d¯K+O(h2)d¯K. (24)

Therefore, we arrive at

 ∣∣|K|−∣∣¯K∣∣∣∣=∣∣∣∫¯K√|g|−1d¯K∣∣∣≤C∣∣¯K∣∣h2

for the error of the cell area which proves (21).

To prove (20) and (22) we consider without loss of generality an edge , where denotes the length of . Considering the derivation of the parametrization

 ce(s)=c(s,0)=(s,0,0)−d(s,0,0)νΓ(ce(s)) (25)

of the curved edge and applying the same arguments as we used to prove (21) completes the proof. ∎

###### Remark 1

Let us note that an analogous estimate to (5) for curved elements is an easy consequence of (5), (20), (21) and the fact which is a consequence of Lemma 3.

###### Lemma 4

There is a constant (depending on ) such that for all flat cells all flat edges and every we have

 ∣∣⟨¯μ¯K,¯e,t−l(x)⟩∣∣ ≤Ch2, (26) ∣∣⟨¯μ¯K,¯e,ν−lΓ(x)⟩∣∣ ≤Ch, (27) ∣∣⟨¯μ¯K,¯e,μ−lK,e(x)⟩−1∣∣ ≤Ch2, (28)

where denotes a unit tangent vector to . We want to point out that this estimate is independent of the sign of

###### Proof

It is sufficient to show versions of (26) - (28) where is substituted by Then analogous results for are immediate and therefore estimates (26) - (28) follow because is the mean of the vectors and . Firstly, we address the proof of (26). Let the same assumptions as in the proof of Lemma 2 hold and in addition let be given by . We obviously have

 ¯ν¯K×¯t¯K,¯e=(0,1,0). (29)

Note that the assumptions of the proof of Lemma 3 are satisfied. Hence we can use (25), i.e. the parametrization of given by satisfies

 c′e(s)=(1,0,0)−νΓ(ce(s))(νΓ)1(ce(s))+O(h2), (30)

and, by definition of , it holds Hence, in view of (22) we obtain

 t−l(x)=(1,0,0)−νΓ(ce(s))(νΓ)1(ce(s))+O(h2) (31)

for some . Combining (29) and (31) we find using (19)

 ∣∣⟨¯ν¯K×¯t¯K,¯e,t−l(x)⟩∣∣=|(νΓ)2(ce(s))(νΓ)1(ce(s))|+O(h2)≤Ch2,

which is (26). Concerning (27),

 ∣∣⟨¯ν¯K×¯t¯K,¯e,ν−lΓ(x)⟩∣∣≤∣∣(ν−lΓ)2(x)∣∣≤Ch

holds because of (29) and (19). Thus, it remains to show (28). By definition form an orthonormal basis of and the vector is of unit length. This means that for every in there exist satisfying such that

 ¯ν¯K×¯t¯K,¯e=b1(¯x)t−l(¯x)+b2(¯x)ν−lΓ(¯x)+b3(¯x)μ−lK,e(¯x). (32)

We know from (26) and (27) that for some which implies using Taylor expansion

 b3(¯x)=±√1+O(h2)=±1+O(h2). (33)

Note that it only remains to show that in (33) the “+” holds. As depends continuously on it is sufficient to find one such that To that end we consider the curve for and small where is the parametrization of from Lemma 3. As is leaving through we have

 0<⟨γ′(0),μK,e(γ(0))⟩. (34)

Due to (32), (33) and the fact that is of unit length we already know that

 μK,e≡±(0,1,0)+O(h). (35)

Inserting (23) for and (35) in (34) leads, together with Lemma 2, to the “+” in (33), which completes the proof. ∎

## 4 Estimating the Difference Between Both Schemes

This section is devoted to establishing a bound for the difference between the curved and flat approximate solutions. To start with we show that the Lax-Friedrichs numerical fluxes from Section 2 satisfy assumption (14).

###### Lemma 5

Let be some compact subset of . Then there is a constant depending only on and such that for the Lax-Friedrichs fluxes (12) and (17) with the same diffusion rate the following inequality holds

 ∣∣LFfnK,e(u,v)−LF¯fn¯K,¯e(u,v)∣∣≤Ch2∀ (u,v)∈K, K∈Th, e⊂∂K.
###### Proof

We start by observing that the diffusive terms drop out, such that

 2∣∣LFfnK,e(u,v)−LF¯fn¯K,¯e(u,v)∣∣=|EnK,e(u)+EnK,e(v)| (36)

with

 EnK,e(u):=\fintIn1|e(tn)|∫e(t)⟨f(u,x,t),μK(t),e(t)(x)⟩de(t)dt−1|In|QIn[1|¯e(tn)|Q¯e(⋅)[⟨f(u,⋅,⋅),¯μ¯K(⋅),¯e(⋅)⟩]].

As and appear symmetrically in (36), we focus on the error analysis of only .