# A class of robust numerical schemes to compute front propagation

Nicolas Therme

<figure><img src="image_1.png" /></figure>

To cite this version:

Nicolas Therme. A class of robust numerical schemes to compute front propagation. 2017. hal-01322625v2

HAL Id: hal-01322625

https://hal.science/hal-01322625v2

Preprint submitted on 21 Feb 2017

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A CLASS OF ROBUST NUMERICAL SCHEMES TO COMPUTE FRONT PROPAGATION

## NICOLAS THERME

**ABSTRACT.** In this work a class of finite volume schemes is proposed to numerically solve equations involving propagating fronts. They fall into the class of Hamilton-Jacobi equations. Finite volume schemes based on staggered grids, and initially developed to compute fluid flows, are adapted to the G-equation, using the Hamilton-Jacobi theoretical framework. The designed scheme has a maximum principle property and is consistent an monotonous on Cartesian grids. A convergence property is then obtained for the scheme on Cartesian grids and numerical experiments evidence the convergence of the scheme on more general meshes.

## 1. INTRODUCTION

The work presented here falls into a larger thematic undertaken for several years, which is the development of staggered schemes to simulate all Mach flows. Numerical schemes were proposed for the Navier-Stokes equations [10], and Euler equations [11, 12]. Adaptations of these schemes to more complex models, such as reactive mixture flows, is underway. In this context, equations describing reactive front propagation are involved and need to be discretized using natural extensions of the staggered schemes.

We focus on a particular equation, used in the combustion science to simulate flame front propagation, the so called G-equation, which reads :

$$ (1) \qquad \partial_t(\rho G) + \text{div}(\rho u G) + \rho u_f |\nabla G| = 0, $$

where $\rho$ is the density of the fluid, $G$ stands for the front indicator, $u$ is a convective velocity and $u_f$ is a front propagation speed. The challenging issue is to adapt staggered discretization to the last term $\rho u_f |\nabla G|$ as the convective part of the equation has already been handled previously. When combined with the mass balance equation of the system

$$ \partial_t \rho + \text{div}(\rho u) = 0, $$

the convective part of the equation is a transport operator and we get:

$$ (2) \qquad \partial_t G + u \cdot \nabla G + u_f |\nabla G| = 0, $$

provided that the density never vanish. This is a particular Hamilton-Jacobi equation. The theory of such equations is well known and was vastly developed by J.-L. Lions in [9, 14]. More precisely, consider the following Cauchy problem:

$$ (3) \qquad \begin{cases} \partial_t G + H(\nabla G) = 0, \\ G(0, x) = G_0(x), \end{cases} $$

defined on $[0, T] \times \mathbb{R}^d$, with $H \in C(\mathbb{R}^d)$ and $G_0 \in \text{BUC}(\mathbb{R}^d)$ (BUC($\Omega$) stands for the set of bounded uniformly continuous functions on $\Omega$). There exists exactly one viscosity solution $G \in \text{BUC}([0, T] \times \mathbb{R}^d)$ such that $G(0, x) = G_0(x)$ and $G$ satisfies a particular weak formulation based on the maximum principle (we refer to [14] for

<footnote n="2010">2010 Mathematics Subject Classification. 35F21, 65N08, 65N12.</footnote>

<footnote>Key words and phrases. Finite volumes, Hamilton-Jacobi, Stability, Convergence.</footnote>
more details). Various numerical methods exist to approach such viscosity solutions. A first converging finite difference scheme was developed in [8]. From this point high order extensions to this scheme were given by S.Osher and James A. Sethian in [15], and a simple finite volume scheme was derived in [13], inspired from a unstructured finite difference scheme based on triangular meshes developed by R. Abgrall in [1]. The convergence theory of numerical approximations of Hamilton Jacobi equations, was first proposed for finite difference scheme in [8] and a generalized formulation was given in [3, 21]. Since then, various schemes were presented for Hamilton-Jacobi equations; high-order finite difference schemes in [6, 19, 16] and schemes for unstructured meshes [5, 20, 22, 2]. These methods are difficult to adapt to our problem. Besides, all the existing schemes proposed in the literature are designed to solve very generic Hamilton-Jacobi equations. In this paper, we only deal with a very particular operator, namely, $H(\mathbf{x}) = \mathbf{u} \cdot \mathbf{x} + u_f |\mathbf{x}|$. Consequently, we propose a finite volume discretization of $u_f |\nabla G|$ that is compatible with the staggered discretization of the transport operator $\mathbf{u} \cdot \nabla G$.

For the sake of clarity, we focus on key elements of the discretization and we suppose that $\mathbf{u} = 0$ and $u_f = 1$, so the problem considered here is the unsteady eikonal equation,

$$ (4a) \qquad \partial_t G + |\nabla G| = 0, $$

$$ (4b) \qquad G(0, \mathbf{x}) = G_0(\mathbf{x}), \quad \forall \mathbf{x} \in \mathbb{R}^d. $$

$G_0 \in \text{BUC}(\mathbb{R}^d)$. The choice of such a simplified model is also convenient as its analytical solutions can be computed easily (see appendix A for more details). The scheme proposed to approximate this problem can be defined on unstructured meshes. On Cartesian grids, The scheme is consistent and monotonous and the $L^\infty$ convergence is proved thanks to the theory developed in [3]. Numerical results are given to highlight this convergence results as well as the numerical convergence of the scheme on unstructured discretizations.

The paper is organized as follows. We start by the description of the spatial discretization and the corresponding notations that are used throughout the paper. We present the scheme and its properties in the second part. We finish with some convergence and numerical results.

## 2. SPATIAL DISCRETIZATION

In this section, we focus on the discretization of a multi-dimensional domain (i.e. $d = 2$ or $d = 3$); the extension to the one-dimensional case is straightforward.

Let $\mathcal{M}$ be a mesh of the domain $\Omega$ (which is an open bounded connected subset of $\mathbb{R}^d$ or $\mathbb{R}^d$ itself), supposed to be regular in the usual sense of the finite element literature (e.g. [7]). The cells of the mesh are assumed to be:

- for a general domain $\Omega$, either non-degenerate quadrilaterals ($d = 2$) or hexahedra ($d = 3$) or simplices, both types of cells being possibly combined in a same mesh,

- for a domain whose boundaries are hyperplanes normal to a coordinate axis, rectangles ($d = 2$) or rectangular parallelepipeds ($d = 3$) (the faces of which, of course, are then also necessarily normal to a coordinate axis).

By $\mathcal{E}$ and $\mathcal{E}(K)$ we denote the set of all $(d-1)$-faces $\sigma$ of the mesh and of the element $K \in \mathcal{M}$ respectively. The set of faces included in the boundary of $\Omega$ is denoted by $\mathcal{E}_{\text{ext}}$ and the set of internal faces (i.e. $\mathcal{E} \setminus \mathcal{E}_{\text{ext}}$) is denoted by $\mathcal{E}_{\text{int}}$; a face $\sigma \in \mathcal{E}_{\text{int}}$ separating the cells $K$ and $L$ is denoted by $\sigma = K|L$. The outward normal vector to a face $\sigma$ of $K$ is denoted by $n_{K,\sigma}$. For $K \in \mathcal{M}$ and $\sigma \in \mathcal{E}$, we denote
by $|K|$ the measure of $K$ and by $|\sigma|$ the $(d-1)$-measure of the face $\sigma$. The mass center of a face is denoted by $\boldsymbol{x}_{\sigma}$.

Finally we denote by $d_{\sigma}$ the measure of $\overline{\boldsymbol{x}_{K}\boldsymbol{x}_{L}}$.

The unknown discrete function $G$ is piecewise constant on the cells $K$. We denote by $H_{\mathcal{M}}$ the space of such piecewise constant functions.

$$G_{\mathcal{M}} \in H_{\mathcal{M}} \iff G_{\mathcal{M}} = \sum_{K \in \mathcal{M}} G_K \chi_K,$$

where $\chi_O$ stands for the characteristic function of the set $O$.

### 3. THE SCHEME

The problem eq. (4) is posed over $\mathbb{R}^d \times (0, T)$, where $(0, T)$ is a finite time interval. Concerning the initial data, we have $G_0 \in \text{BUC}(\mathbb{R}^d)$. According to the known results at the continuous level, the problem has a unique viscosity solution in $\text{BUC}([0, T] \times \mathbb{R}^d)$, that we denote $\tilde{G}$. In order to be able to perform computations, the domain can be reduced to an open bounded subset $\Omega$ of $\mathbb{R}^d$ with zero-flux boundary conditions. We propose three versions of the scheme depending on the regularity of the mesh. The finite volume scheme is written on an alternative form of Equation eq. (4a) :

$$(5) \qquad \partial_t G + \left( \frac{\nabla G}{|\nabla G|} \right) \cdot \nabla G = 0,$$

and makes use of the classical identity:

$$(6) \qquad \boldsymbol{u} \cdot \nabla \phi = \operatorname{div}(\phi \boldsymbol{u}) - \phi \operatorname{div}(\boldsymbol{u}).$$

Let us consider a partition $0 = t_0 < t_1 < \dots < t_N = T$ of the time interval $(0, T)$, which we suppose uniform for the sake of simplicity, and let $\delta t = t_{n+1} - t_n$ for $n = 0, 1, \dots, N-1$ be the (constant) time step. We consider an explicit-in-time scheme, which reads, for $0 \le n \le N-1$ and $K \in \mathcal{M}$:

$$(7) \qquad \partial_t G^n + F_{\mathcal{M}}(G^n) = 0,$$

with,

$$(8) \qquad \partial_t G^n = \sum_{K \in \mathcal{M}} \frac{G_K^{n+1} - G_K^n}{\delta t} \chi_K,$$

and

$$(9) \qquad F_{\mathcal{M}}(G^n) = \operatorname{div} \left( \frac{\nabla \varepsilon G^n}{|\nabla \varepsilon G^n|} G^n \right)_K - G_K^n \operatorname{div} \left( \frac{\nabla \varepsilon G^n}{|\nabla \varepsilon G^n|} \right)_K.$$

The discrete divergence operator is given by:

$$(10) \qquad \text{for } K \in \mathcal{M}, \quad (\operatorname{div} \boldsymbol{u})_K = \frac{1}{|K|} \sum_{\sigma=K|L \in \mathcal{E}(K)} \kappa_{K,\sigma}^{M} |\sigma| \boldsymbol{u}_{\sigma} \cdot \boldsymbol{n}_{K,\sigma},$$

where $\kappa_{K,\sigma}^M$ is a coefficient equal to 1 for unstructured meshes, and equal to $\kappa_{K,\sigma}^M = \frac{2|K|}{|K| + |L|}$ on Cartesian grids. Likewise

$$(11) \qquad \text{for } K \in \mathcal{M}, \quad (\operatorname{div} G \boldsymbol{u})_K = \frac{1}{|K|} \sum_{\sigma=K|L \in \mathcal{E}(K)} \kappa_{K,\sigma}^{M} |\sigma| G_{\sigma} \boldsymbol{u}_{\sigma} \cdot \boldsymbol{n}_{K,\sigma},$$

where $G_{\sigma}$ denotes an interpolation of $G$ on the edge $\sigma$ that is:

$$\text{for } \sigma = K|L \in \mathcal{E}_{\text{int}}, \quad G_{\sigma} = \begin{cases} G_K & \text{if } \boldsymbol{u}_{\sigma} \cdot \boldsymbol{n}_{K,\sigma} \ge 0, \\ G_L & \text{otherwise.} \end{cases}$$
For a face $\sigma \in \mathcal{E}_{\text{ext}}$ one simply take $G_{\sigma} = G_K$ so that

$$ \nabla G \cdot n_{K,\sigma} = \frac{|\sigma|}{|K|} (G_{\sigma} - G_K) = 0. $$

The expression of the discrete spatial operator eq. (9) becomes

$$ (12) \quad F_M(G_M^n) = \sum_{K \in \mathcal{M}} \left[ \sum_{\sigma=K|L \in \mathcal{E}(K)} \kappa_{K,\sigma}^{\mathcal{M}} \frac{|\sigma|}{|K|} \frac{(\nabla_{\mathcal{E}} G^n)_{\sigma}}{|(\nabla_{\mathcal{E}} G^n)_{\sigma}|} \cdot n_{K,\sigma} (G_{\sigma}^n - G_K^n) \right] \chi_K, $$

where $\nabla_{\mathcal{E}}$ refers to a discrete gradient operator defined on every $\sigma \in \mathcal{E}_{\text{int}}$.

**3.1. Unstructured meshes.** For $\sigma = K|L \in \mathcal{E}_{\text{int}}$, we take:

$$ (13) \quad (\nabla_{\mathcal{E}} G)_{\sigma} = \sum_{\sigma \in \partial(K \cup L)} \frac{|\sigma|}{|K \cup L|} \tilde{G}_{\sigma} n_{K \cup L, \sigma}, $$

with $\tilde{G}_{\sigma}$ a second order approximation of $G$ at the barycenter of the face $\sigma$.

**3.2. Cartesian meshes.** When the scheme is based on Cartesian grids, we have for $\sigma = \overrightarrow{K|L}$ (which means the flow goes from $K$ to $L$):

$$ (14) \quad \text{For } \sigma \in \mathcal{E}_{\text{int}}, \quad (\nabla_{\mathcal{E}} G)_{\sigma} = \left[ \frac{G_L - G_K}{d_{\sigma}} n_{K,\sigma} + \nabla_{\parallel\sigma} G \right], $$

where $\nabla_{\parallel\sigma}^C$ is defined by:

$$ (15) \quad (\nabla G)_{\parallel\sigma}^C = \sum_{\substack{i=1, \\ e^{(i)} \cdot n_{K,\sigma}=0}}^{d} \frac{(G_{K_i^+} - G_K)^+}{d_{\sigma_i^+}} - \frac{1}{2} (1 - \operatorname{sgn}(G_{K_i^+} - G_K)^+) \frac{(G_K - G_{K_i^-})^-}{d_{\sigma_i^-}} e^{(i)}, $$

with $\sigma = \overrightarrow{K|L}$. For a cell $K \in \mathcal{M}$, $\sigma_i^+$ and $\sigma_i^-$ stand for the two faces of $K$ normal to $e^{(i)}$. Superscripts - and + refer to the up and down faces of $K$ respectively. We set $\sigma_i^+ = K|K_i^+$ and $\sigma_i^- = K|K_i^-$. We illustrate these notations in the following figure. We recall that $a^+ = \max(a, 0)$ and $a^- = \max(-a, 0)$, for $a \in \mathbb{R}$.

<figure><img src="image_3.png" /><figcaption>FIGURE 1. Notations for the alternative gradient definition on Cartesian grids with $\mathbf{F} = (G_L - G_K)n_{K,\sigma}$.</figcaption></figure>
### 3.3. High order extension.
It is possible to replace the upwind interpolation by a higher order interpolation based on a MUSCL reconstruction. Adopting the same notations as in eq. (11), its important property, based on [18] is stated below. For any $K \in \mathcal{M}$, and for any $\sigma \in \mathcal{E}(K) \cap \mathcal{E}_{\text{int}}$, there exists $\alpha_{K,\sigma} \in [0, 1]$ such that :

$$ (16) \qquad G_{\sigma} - G_{K} = \begin{cases} \beta_{K,\sigma}(G_{K} - G_{M_{\sigma}^{K}}) & \text{if } \frac{\nabla_{\epsilon} G_{\sigma}^{n}}{|\nabla_{\epsilon} G_{\sigma}^{n}|} \cdot n_{K,\sigma} \ge 0, \\ \beta_{K,\sigma}(G_{M_{\sigma}^{K}} - G_{K}) & \text{otherwise.} \end{cases} $$

The procedure is the following:

* We define a tentative value $\tilde{G}_{\sigma}$ based on a high order geometric interpolation.

* The next step is to create a limitation procedure for $\rho_{\sigma}$ and $e_{\sigma}$. Let $\sigma \in \mathcal{E}_{\text{int}}$, $\sigma = K|\overleftarrow{L}$ and $V_K$ a set of neighboring cells to K. We make the two following assumptions :

$$ (H1) \qquad G_{\sigma} - G_{K} \in \left| \left[ 0, \frac{\zeta^{+}}{2} (G_{L} - G_{K}) \right] \right| \\ (H2) \qquad \exists M \in V_{K}, G_{\sigma} - G_{K} \in \left| \left[ 0, \frac{\zeta^{-}}{2} \frac{d_{\sigma}}{d_{K|M}} (G_{K} - G_{M}) \right] \right|, $$

where, for $a, b \in \mathbb{R}$, we denote by $|[a, b]|$ the convex hull of a and b and $K|\overleftarrow{L}$ means that the flow is going from K to L ($\frac{\nabla_{\epsilon} G_{\sigma}^{n}}{|\nabla_{\epsilon} G_{\sigma}^{n}|} \cdot n_{K,\sigma} \ge 0$). The parameters $\zeta^{+}$ and $\zeta^{-}$ lie in $[0, 2]$.

* We compute $G_{\sigma}$ as the nearest point to $\tilde{G}_{\sigma}$ in the limitation interval.

Whenever it is possible (i.e. with a mesh obtained by $Q_1$ mappings from the $(0, 1)^d$ reference element), $V_K$ may be chosen as the opposite cells to $\sigma$ in K. Otherwise $V_K$ is defined as the set of "upstream cells" to K. Note that, for a structured mesh, the first choice allows to recover the usual minmod limiter.

**Remark 3.1 (Cartesian grids).** We impose $\zeta^{+} = \zeta^{-} = 1$ for the Cartesian version of the scheme. This particular choice of parameters is the only one possible if we want to get consistency properties for the discrete spatial operator of the scheme.

## 4. PROPERTIES OF THE SCHEME

We expose in this section the properties of the scheme. Specific paragraph is devoted to its additional properties on Cartesian grids, derived from the convergence theory [3, 21]. This ensures that the given discretization behaves like usual finite difference methods for Hamilton-Jacobi equations.

### 4.1. Stability.
Thanks to the definition of the discrete convective operator, we have the following property:

**Proposition 4.1 (Maximum principle).** Let $G_{\mathcal{M}}^{n} \in H_{\mathcal{M}}$, $n \in [0, N]$, be the solution of the scheme eq. (7). For all $K \in \mathcal{M}$ and $n \in [0, N-1]$, we have:

$$ \min_{L \in \mathcal{M}} G_{L}^{n} \le G_{K}^{n+1} \le \max_{L \in \mathcal{M}} G_{L}^{n}, $$

under the CFL condition:

$$ (18) \qquad \delta t \le \min_{K \in \mathcal{M}} \frac{|K|}{\sum_{\sigma \in \mathcal{E}(K)} |\sigma|} $$
Proof. We have, for $K \in \mathcal{M}$ and $n \in [0, N-1]$:

$$ G_K^{n+1} = \left( 1 - \delta t \sum_{\sigma \in \mathcal{E}(K)} \frac{|\sigma|}{|K|} \left( \frac{\nabla_{\mathcal{E}} G_{\sigma}^n}{|\nabla_{\mathcal{E}} G_{\sigma}^n|} \cdot \mathbf{n}_{K,\sigma} \right)^- \right) G_K^n + \delta t \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{|\sigma|}{|K|} \left( \frac{\nabla_{\mathcal{E}} G_{\sigma}^n}{|\nabla_{\mathcal{E}} G_{\sigma}^n|} \cdot \mathbf{n}_{K,\sigma} \right)^- G_L^n. $$

Consequently, $G_K^{n+1}$ is a convex combination of its neighbors at time $n$ if eq. (18) is verified, which completes the proof. □

**Remark 4.1** (Cartesian grids). The property remains the same with the scheme on Cartesian grids, only the CFL is modified. One must replace $|K|$ by $\frac{|K|+|L|}{2}$ in eq. (18).

**Remark 4.2** (MUSCL interpolation). Concerning the MUSCL interpolation, we use the property eq. (16) and use it in the scheme to get:

$$ G_K^{n+1} = \left( 1 - \delta t \sum_{\sigma \in \mathcal{E}(K)} \frac{|\sigma|}{|K|} \beta_{K,\sigma} \frac{\nabla_{\mathcal{E}} G_{\sigma}^n}{|\nabla_{\mathcal{E}} G_{\sigma}^n|} \cdot \mathbf{n}_{K,\sigma} \right) G_K^n + \delta t \sum_{\sigma \in \mathcal{E}(K)} \frac{|\sigma|}{|K|} \beta_{K,\sigma} \frac{\nabla_{\mathcal{E}} G_{\sigma}^n}{|\nabla_{\mathcal{E}} G_{\sigma}^n|} \cdot \mathbf{n}_{K,\sigma} |G_{M_K^n}|. $$

The maximum principle is still satisfied with the same CFL condition.

## 5. INVARIANCE UNDER TRANSLATION

**Proposition 5.1** (Invariance under Translation with constants). 

$\forall \lambda \in \mathbb{R}$, and $\forall \phi_M \in H_M$,

$$ (19) \qquad F_M(\phi_M + \lambda) = F_M(\phi_M). $$

Proof. Let $\lambda \in \mathbb{R}$ and $\phi_M \in H_M$. Looking at eq. (12), we need to check that $\nabla_{\mathcal{E}}(\phi_M + \lambda) = \nabla_{\mathcal{E}}\phi_M$. We remind that:

$$ \nabla_{\mathcal{E}}(\phi_M + \lambda) = \sum_{\sigma \in \partial K \cup L} \frac{|\sigma|}{|K \cup L|} (\phi_\sigma + \lambda) \mathbf{n}_{K \cup L, \sigma} $$

We have:

$$ \nabla_{\mathcal{E}}(\phi_M + \lambda) = \nabla_{\mathcal{E}}\phi_M + \lambda \sum_{\sigma \in \partial K \cup L} \frac{|\sigma|}{|K \cup L|} \mathbf{n}_{K \cup L, \sigma}. $$

Using the divergence theorem, we get that:

$$ \sum_{\sigma \in \partial K \cup L} \frac{|\sigma|}{|K \cup L|} \mathbf{n}_{K \cup L, \sigma} = \int_{K \cup L} \nabla(1) = 0, $$

which concludes the proof. □

On Cartesian meshes, the result is immediate.

### 5.1. Properties of the Cartesian scheme.
5.1.1. Consistency. We need to define interpolates of test functions on the mesh. Let $\phi \in C_c^\infty(\Omega)$. We set:

$$ (20) \qquad \phi_M = \sum_{K \in \mathcal{M}} \phi_K \chi_K, \qquad \phi_K = \phi(\mathbf{x}_K). $$

We now give the definition of the consistency property.

**Definition 5.1 (Consistency).** Let $F(G)$ be an operator approximated by $F_M(G_M)$. Let $h_M = \max_{K \in \mathcal{M}} \mathrm{diam}(K)$. Let $\mathcal{D}^{(m)} = \{\mathcal{M}^{(m)}, \mathcal{E}^{(m)}, \mathcal{P}^{(m)}\}$ be a sequence of discretizations such that the size $h_M^{(m)}$ tends to zero as $m \to \infty$. The discrete spatial operator $H_M$ is said to be consistent with $H$ if for every $\phi \in C_c^\infty(\Omega)$:

$$ \lim_{m \to \infty} \|F_{M^{(m)}}(\phi_{M^{(m)}}) - F(\phi)\|_{L^\infty(\Omega)} = 0. $$

The next proposition states

**Proposition 5.2.** The spatial operator in the Cartesian case, given by, for $G_M \in H_M$:

$$ (21) \quad F_M(G_M) = \sum_{K \in \mathcal{M}} \left[ \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_{\sigma}} \frac{(G_L - G_K)}{\sqrt{(G_L - G_K)^2 + d_{\sigma}^2 |\nabla_{//\sigma} G_M|^2}} (G_{\sigma} - G_K) \right] \chi_K, $$

is consistent with $|\nabla G|$.

*Proof.* Let $\phi \in C_c^\infty(\Omega)$ and $\phi_M \in H_M$ its interpolation on the mesh. Consider $K \in \mathcal{M}$ and $\mathbf{v}$ a constant vector. Let $\tilde{F}_K(\phi_M, \mathbf{v})$ be:

$$ \tilde{F}_K(\phi_M, \mathbf{v}) = \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_{\sigma}} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})(\phi_{\sigma} - \phi_K). $$

With the upwind interpolation, we get that:

$$ \tilde{F}_K(\phi_M, \mathbf{v}) = - \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_{\sigma}} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^{-}(\phi_L - \phi_K). $$

A simple Taylor expansion leads to:

$$ \tilde{F}_K(\phi_M, \mathbf{v}) = - \sum_{\sigma=K|L \in \mathcal{E}(K)} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^{-} \nabla \phi(\mathbf{x}_K) \cdot \mathbf{n}_{K,\sigma} + \mathcal{O}(h_M), $$

so

$$ \tilde{F}_K(\phi_M, \mathbf{v}) = \nabla \phi(\mathbf{x}_K) \cdot \sum_{\sigma=K|L \in \mathcal{E}(K)} (\mathbf{v} \cdot \mathbf{n}_{L,\sigma})^{+} \mathbf{n}_{L,\sigma} + \mathcal{O}(h_M). $$

Thanks to the Cartesian grid, we have:

$$ \sum_{\sigma=K|L \in \mathcal{E}(K)} (\mathbf{v} \cdot \mathbf{n}_{L,\sigma})^{+} \mathbf{n}_{L,\sigma} = \sum_{i=1}^{d} (\mathbf{v} \cdot \mathbf{e}^{(i)}) \mathbf{e}^{(i)} = \mathbf{v}, $$

so we have:

$$ \tilde{F}_K(\phi_M, \mathbf{v}) = \mathbf{v} \cdot \nabla \phi(\mathbf{x}_K) + \mathcal{O}(h_M). $$
Concerning the MUSCL interpolation, we have:

$$
\begin{aligned}
\tilde{F}_K(\phi_M, \mathbf{v}) ={}& \frac{1}{2} \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^+ \min \left( \phi_K - \phi_{M_K^\sigma} \frac{d_\sigma}{d_{K|M_K^\sigma}}, \phi_L - \phi_K \right) \\
& - \sum_{\sigma=K|L \in \mathcal{E}(K)} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^{-} (\phi_L - \phi_K) \\
& - \frac{1}{2} \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^{-} \min \left( \phi_L - \phi_{M_L^\sigma} \frac{d_\sigma}{d_{L|M_L^\sigma}}, \phi_K - \phi_L \right),
\end{aligned}
$$

where $M_K^\sigma$ refers to the opposite cell to $\sigma$ in $K$. It is easy to see that:

$$ \frac{1}{d_{\sigma}} \min \left( \phi_{K} - \phi_{M_{K}^{\sigma}} \frac{d_{\sigma}}{d_{K|M_{K}^{\sigma}}}, \phi_{L} - \phi_{K} \right) = \nabla \phi(\mathbf{x}_{K}) \cdot \mathbf{n}_{K,\sigma} + \mathcal{O}(h_{M}), $$

and,

$$ \frac{1}{d_{\sigma}} \min \left( \phi_{L} - \phi_{M_{L}^{\sigma}} \frac{d_{\sigma}}{d_{L|M_{L}^{\sigma}}}, \phi_{K} - \phi_{L} \right) = \nabla \phi(\mathbf{x}_{K}) \cdot \mathbf{n}_{L,\sigma} + \mathcal{O}(h_{M}). $$

Therefore,

$$ \tilde{F}_K(\phi_M, \mathbf{v}) = \frac{1}{2} \sum_{\sigma \in \mathcal{E}(K)} (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^+ \nabla \phi(\mathbf{x}_K) \cdot \mathbf{n}_{K,\sigma} + \\ \frac{1}{2} \sum_{\sigma \in \mathcal{E}(K)} (\mathbf{v} \cdot \mathbf{n}_{L,\sigma})^+ \nabla \phi(\mathbf{x}_K) \cdot \mathbf{n}_{L,\sigma} + \mathcal{O}(h_M), $$

which leads to:

$$
\begin{aligned}
\tilde{F}_K(\phi_M, \mathbf{v}) &= \nabla\phi(\mathbf{x}_K) \cdot \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{2} \left( (\mathbf{v} \cdot \mathbf{n}_{K,\sigma})^+ \mathbf{n}_{K,\sigma} + (\mathbf{v} \cdot \mathbf{n}_{L,\sigma})^+ \mathbf{n}_{L,\sigma} \right) \\
&\quad + \mathcal{O}(h_M) \\
&= \nabla\phi(\mathbf{x}_K) \cdot \mathbf{v} + \mathcal{O}(h_M).
\end{aligned}
$$

Noticing, thanks to the consistency of $\nabla_{\mathcal{E}}$, that:

$$ F_M(\phi_M) = \sum_{K \in \mathcal{M}} \tilde{F}_K \left( \phi_M, \frac{\nabla \phi(\mathbf{x}_K)}{|\nabla \phi(\mathbf{x}_K)|} \right) \mathcal{X}_K + \mathcal{O}(h_M), $$

we can conclude that:

$$ \lim_{m \to \infty} F_M(\phi_M) = |\nabla \phi|, $$

which concludes the proof. $\square$

<figure><img src="image_4.png" /></figure>

5.1.2. *Monotonicity.* Let $(\phi_M, \psi_M) \in H_M$. Let us define the following partial order

$$ (22) \qquad \phi_M \le \psi_M \iff \forall K \in \mathcal{M}, \quad \phi_K \le \psi_K. $$

Then we get the following result with the Cartesian upwind scheme only.

**Proposition 5.3 (Monotonicity of the upwind Cartesian scheme).**

Suppose that the following CFL condition is satisfied

$$ (23) \qquad \delta t \le \frac{1}{\sum_{\sigma \in \mathcal{E}(K)} \frac{1+\frac{1}{2}\sqrt{1+r^2}}{d_{\sigma}}}, \qquad r = \max_{(\sigma, \sigma') \in \mathcal{E}(K)} \frac{d_{\sigma}}{d_{\sigma'}}. $$

Then he have the following result:

$$ \forall (\phi_M, \psi_M) \in H_M, \quad \phi_M \le \psi_M \implies \phi_M + \delta t F_M(\phi_M) \le \psi_M + \delta t F_M(\psi_M). $$
Proof. For the sake of clarity we prove the result in 2D. The extension to all dimension can be done at the cost of heavier notations and CFL conditions. We can equivalently check that SCH is a non decreasing function of each variable. Let $K \in \mathcal{M}$ and $\phi_M \in H_M$. We have:

$$ \phi_M + \delta t F_M(\phi_M)|_K = \phi_K + \delta t \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} f_{K,\sigma}(\phi_M), $$

with,

$$ f_{K,\sigma}(\phi_M) = \frac{(\phi_L - \phi_K)^-}{\sqrt{(\phi_L - \phi_K)^2 + d_\sigma^2 |\nabla_{//\sigma} \phi_M|^2}} (\phi_L - \phi_K) $$

The monotonicity of $f_{K,\sigma}$ in $\phi_L$ is equivalent to the monotonicity of the function:

$$ f: x \mapsto \frac{x^{-}x}{|x|} = -x^{-}, \quad \forall x \in \mathbb{R} $$

because $\nabla_{//\sigma}\phi_M$ does not depend on $\phi_L$ in the Cartesian case (see eq. (15)). We can conclude that $f_{K,\sigma}$ is a non decreasing function of $\phi_L$. Concerning the monotonicity in $\phi_{K-}$ and $\phi_{K+}$ it is equivalent to the variations of:

$$ f: x \mapsto -\frac{1}{x^{+}}, $$

which is a non decreasing function. We can conclude that $SCH(\phi_M)|_K$ is an increasing function of each $(\phi_M)_{M \neq K}$. Concerning $\phi_K$, we have:

$$ SCH(\phi_M)|_K = g(\phi_K) = \phi_K - \delta t \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} \frac{(\phi_K - \phi_L)^+ (\phi_K - \phi_L)}{\sqrt{(\phi_L - \phi_K)^2 + d_\sigma^2 |\nabla_{//\sigma} \phi_M|^2}} $$

The analysis of this function can be split into three cases. If, $\forall \sigma \in \mathcal{E}(K)$, $\phi_K \le \phi_L$, then $g(\phi) = \phi_K$ which is non decreasing. The second case is when, $\forall \sigma \in \mathcal{E}(K)$, $\phi_K \ge \max(\phi_{K+}, \phi_{K-}, \phi_L)$. We have:

$$ g(\phi_K) = \phi_K - \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{\delta t}{d_\sigma} (\phi_K - \phi_L). $$

which is non decreasing if,

$$ \delta t \le \frac{1}{\sum_{\sigma \in \mathcal{E}(K)} d_\sigma^{-1}}. $$

Finally, suppose that $\forall \sigma \in \mathcal{E}(K)$, $\phi_L \le \phi_K \le \phi_{K+}$ (or $\phi_{K-}$), we have, denoting by $r_\sigma = \frac{d_\sigma}{d_{\sigma+}}$:

$$ g(\phi_K) = \phi_K - \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} \frac{\phi_K - \phi_L}{\sqrt{(\phi_K - \phi_L)^2 + r_\sigma^2 (\phi_K - \phi_{K+})^2}} (\phi_K - \phi_L). $$

Let us derive this function:

$$ g'(\phi_K) = 1 - \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} \frac{\phi_K - \phi_L}{\sqrt{(\phi_K - \phi_L)^2 + r_\sigma^2 (\phi_K - \phi_{K+})^2}} - \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} \frac{r_\sigma^2 (\phi_{K+} - \phi_K)(\phi_K - \phi_L)(\phi_{K+} - \phi_K)}{((\phi_K - \phi_L)^2 + r_\sigma^2 (\phi_K - \phi_{K+})^2)^{3/2}} $$

One can notice directly that:

$$ \sum_{\sigma=K|L \in \mathcal{E}(K)} \frac{1}{d_\sigma} \frac{\phi_K - \phi_L}{\sqrt{(\phi_K - \phi_L)^2 + r_\sigma^2 (\phi_K - \phi_{K+})^2}} \le 1. $$
In order to upper-bound the second sum, we analyze the function

$$h: x \mapsto \frac{r^2 x(a-x)a}{(x^2 + d^2(a-x)^2)^{3/2}},$$

where $a, r$ are strictly positive constants. We split the function in two parts $h(x) = h_1(x)h_2(x)$ with:

$$h_1(x) = \frac{r^2 x(a-x)}{x^2 + r^2(a-x)^2},$$

$$h_2(x) = \frac{a}{\sqrt{x^2 + r^2(a-x)^2}}.$$

Concerning $h_1$ we can equivalently consider the function defined on $\mathbb{R}^+$ by:

$$y \mapsto \frac{r^2}{y + \frac{r^2}{y}} = \frac{r^2 y}{y^2 + r^2}.$$

A quick study of the function shows that,

$$\max_{y \in \mathbb{R}^+} \frac{r^2 y}{y^2 + r^2} = \frac{r}{2} = \max_{x \in [0, a]} h_1(x).$$

The same work is performed with $h_2$ and leads to:

$$\max_{x \in [0, a]} h_2(x) = \frac{\sqrt{1 + r^2}}{r}$$

Gathering the results, we get that:

$$\forall x \in [0, a], \quad h(x) \le \frac{1}{2} \sqrt{1 + r^2}$$

As a result, writing out $r = \max_{(\sigma, \sigma') \in \mathcal{E}(K)} \frac{d_\sigma}{d_{\sigma'}}$, we get that $g'(\phi_K) \ge 0$ provided that eq. (23) is satisfied. This CFL condition ensures that $\phi_M + \delta t F_M(\phi_M)|_K$ is a non-decreasing function of $\phi_K$, which concludes the proof. □

**Remark 5.1.** All the results proved here can be extended with a transport velocity $u \ne 0$ and a front propagation speed $u_f \ne 1$. Only the CFL conditions are modified, the sketch of the proofs is the same. However the monotonicity results cannot be extended to the MUSCL interpolation, and more generally to the non Cartesian case.

## 6. A CONVERGENCE RESULT IN THE CARTESIAN CASE

The previous section ensures that the upwind scheme satisfies the basic properties to seek a convergence result on Cartesian meshes. We first recall the theorem given in [3], adapted to our notations.

**Theorem 6.1.** Let $\mathcal{D}^{(m)} = \{M^{(m)}, \mathcal{E}^{(m)}, \mathcal{P}^{(m)}, \delta t^{(m)}\}$ be a sequence of discretizations such that the space and time steps tend to zero as $m \to \infty$. Let $\tilde{G}$ be the viscosity solution of eq. (4). Consider the following explicit scheme, for $n \in [0, N-1]$:

$$\delta t G_m^n + F_M(G_m^n) = 0,$$

and the complete solution defined by $G_m^{(T)} = \sum_{n=0}^{N-1} G_m^{n+1} \chi_{[t^n, t^{n+1}]}$. We suppose that:

* The spatial operator $F_M$ is strongly consistent with the continuous operator $G \mapsto |\nabla G|$.

* The scheme is invariant under translations: $F_M(G_M + v) = F_M(G_M)$.

* The scheme is monotone.
Then,

$$G_{\mathcal{M}^{(m)}} \to \bar{G} \text{ uniformly as } m \to \infty.$$

Since we have shown the required properties in Theorem 6.1, we can thus conclude to the convergence of the scheme, which we state in the following corollary.

**Corollary 6.2.** Let $\mathcal{D}^{(m)} = \{\mathcal{M}^{(m)}, \mathcal{E}^{(m)}, \mathcal{P}^{(m)}, \delta t^{(m)}\}$ be a sequence of discretizations such that the space and time steps tend to zero as $m \to \infty$. Now suppose there exists $r > 0$, such that $\forall m \in \mathbb{N}, \forall (\sigma, \sigma') \in \mathcal{E}^{(m)}$,

$$\frac{d_{\sigma}}{d_{\sigma'}} \leq r.$$

Suppose that, for any $m \in \mathbb{N}$,

$$\delta t^{(m)} \leq \max_{K \in \mathcal{M}^{(m)}} \frac{1}{\sum_{\sigma' \in \mathcal{E}(K)} \frac{1+\frac{1}{2}\sqrt{1+r^2}}{d_{\sigma}}}.$$

Then the solution of the upwind Cartesian scheme eq. (7)-eq. (21) $G_m^{(T)}$ converges uniformly towards $\bar{G}$.

## 7. NUMERICAL RESULTS

### 7.1. One dimension.
The domain is $\Omega = (0, 1)$. We use zero-flux boundary conditions in $x = 0$ and $x = 1$. We suppose that the time and space steps are constant for simplicity. Consider the following initial data:

$$(24) \qquad G_0(x) = |\sin(4\pi x)|$$

We give the solution at $T = 0.05s$, with an upwind interpolation for the spatial operator, and a fixed CFL equal to 1/10.

<figure><img src="image_5.png" /><figcaption>FIGURE 2. Solution of the G-equation with the upwind scheme at $T = 0.05s$.</figcaption></figure>

It is possible to determine the unique viscosity solution of the eikonal equation for a given bounded uniformly continuous initial data. The expression of the solution is given by eq. (26) and its proof can be found in the appendix A. Consequently we can highlight numerically the theoretical result about the convergence of the
solution of our scheme towards the viscosity solution. The figure below gives the error in $L^1$ norm according to the space step, for a fixed CFL equal to $\frac{1}{10}$.

<figure><img src="image_6.png" /><figcaption>FIGURE 3. L1 norm error at T=0.05s and CFL= $\frac{1}{10}$ – upwind interpolation.</figcaption></figure>

We can also see the behavior of the scheme if we use discontinuous initial data. We consider the following:

$$G_0(x) = \begin{cases} 0, & \text{if } x \le 0.5 \\ 1, & \text{otherwise.} \end{cases}$$

The result at time $T = 0.2s$ is given below, for the upwind scheme and the MUSCL scheme.
<figure><img src="image_7.png" /><figcaption>FIGURE 4. solution at T = 0.05s and CFL= 1/10 with h = 10⁻³.</figcaption></figure>

The MUSCL scheme brings less numerical diffusion, as expected. Normally one can not define a viscosity solution for discontinuous initial data. However one expects the solution to be the same as the general viscosity solution given for BUC initial data (see eq. (26) in the Appendix).

## 7.2. Two Dimensions.

7.2.1. *Unstructured grid.* The computational domain is $\Omega = [-\frac{1}{2}, \frac{1}{2}]^2$. The mesh consists in convex quadrilaterals. We give an example of the discretization below. These grids are built from a regular Cartesian grid for which a random displacement

<figure><img src="image_8.png" /><figcaption>FIGURE 5. Example of a 10 × 10 unstructured grid</figcaption></figure>

of length $\epsilon h$ is applied to each node where $h$ is the space step. We consider zero-flux boundary conditions. The initial data are given in the polar coordinates $(r, \theta)$:

$$G_0(r, \theta) = r \left( 1 + \frac{1}{2} \cos(4\theta) \right).$$
Results obtained at different times are given below. The scheme used is the upwind version for unstructured meshes, with a space step $h = \frac{1}{200}$, and a constant CFL equal to $\frac{1}{10}$.

<figure><img src="image_9.png" /><figcaption>FIGURE 6. G at different times with the upwind scheme on an unstructured mesh - h = $\frac{1}{200}$ - CFL = $\frac{1}{10}$</figcaption></figure>

Another possible test case is the following one:

$$ (25) \qquad G_0(r, \theta) = |\sin(4\pi r)|. $$

Results obtained with different meshes are displayed just below. The scheme used is the upwind version for unstructured meshes, with a space step $h = \frac{1}{400}$, a constant CFL equal to $\frac{1}{10}$ and a final time equal to $T = 0.04s$.
<figure><img src="image_10.png" /><figcaption>FIGURE 7. G on different meshes -T = 0.04s – h = $\frac{1}{400}$ – CFL = $\frac{1}{10}$</figcaption></figure>

Finally we plot some convergence results. Let $G_{\text{visc}}$ be the viscosity solution associated to the initial data eq. (25). We take $G_{\text{visc}}(\cdot, T = 0.01s)$ as a new initial data. The final time is set to $T = 0.04s$. The results are given below, with a constant CFL equal to $\frac{1}{10}$, using three different meshes: an unstructured mesh with a deformation ratio equal to $\epsilon = 0.1$, a triangular mesh which consists of a
square grid where each square is cut in half following the same diagonal, and a Rhomboidal mesh composed of parallelograms with a large angle equal to $2\pi/3$.

<figure><img src="image_11.png" /><figcaption>FIGURE 8. $L^1$ norm error at T=0.05s and CFL= $1/10$ - upwind interpolation.</figcaption></figure>

7.2.2. Cartesian grids. We use the same test to compare the convergence of the MUSCL scheme, the upwind scheme, and an upwind finite difference scheme described in [8] designed for the Hamilton-Jacobi equations. In order to properly observe a difference in the convergence rate we use a Runge-Kutta time discretization of order two.

To conclude, we introduce a test case with a convective velocity $u$ different from zero. Let the computational domain be $\Omega = (0, 1)^2$. Zero-flux boundary conditions are prescribed on the boundary. We consider the following initial data

$$G_0(x) = \begin{cases} 0, & \text{if } \|x - (0.25, 0.8)\| \le 0.15 \\ 1, & \text{otherwise.} \end{cases}$$

The front propagation velocity is equal to $u_f = 0.8$ and the convective velocity corresponds to a vortex centered around $(0, 0)$ with a constant angular speed equal to $2\pi$, namely

$$u = 2\pi r e_{\theta},$$

in polar coordinates.

The upwind scheme is used on a 400 × 400 Cartesian grid with a CFL equal to $1/20$. Results are plotted below.

Numerical simulations performed in this section are in good agreement with the properties verified by the scheme. The discretization proposed in this paper has been implemented on a Computational Fluid Dynamics software called P²REMICS [17]. It will be used on a complex model with reactive flows to simulate the flame front propagation in the explosion phenomenon, for nuclear safety issues. Work on this model and its related numerical simulations is underway.
<figure><img src="image_12.png" /><figcaption>FIGURE 9. L1 norm error at $T = 0.05s$ and CFL= $\frac{1}{10}$</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>FIGURE 10. $G$ at $T = 0s$ (left) – $T = 0.1s$ (right) with a CFL= $\frac{1}{20}$</figcaption></figure>
APPENDIX A. VISCOSITY SOLUTIONS OF THE EIKONAL EQUATION

It is possible to compute the viscosity solution of 4 for every $G_O \in BUC(\mathbb{R}^d)$.
It is then defined on $\mathbb{R}^d \times (0, +\infty)$ by:

$$
(26) \qquad G(\mathbf{x}, t) = \inf_{|\mathbf{x}-\mathbf{y}|\le t} G_0(\mathbf{y}).
$$

The proof of this result can be found in [4], and it is based on the following lemma

**Lemma A.1.** Let us set

$$
S(t)G(\mathbf{x}) = \inf_{|\mathbf{x}-\mathbf{y}|\le t} G(\mathbf{y}).
$$

Then S is a monotonous semigroup on C($\mathbb{R}^d$).

*Proof.* The proof is rather simple as

$$
S(t) \circ S(s) G(\mathbf{x}) = \inf_{|\mathbf{x}-\mathbf{y}|\le t} \left( \inf_{|\mathbf{z}-\mathbf{y}|\le s} G(\mathbf{z}) \right).
$$

This computation is equivalent to seek the infimum in the set

$$
\{ \mathbf{z} \text{ such that } \exists \mathbf{y} \text{ such that } |\mathbf{x} - \mathbf{y}| \le t \text{ and } |\mathbf{z} - \mathbf{y}| \le s \}.
$$

Now, this set is equal to the set

$$
\{z\} \text{ such that } |\mathbf{x} - \mathbf{z}| \leq t + s,
$$

so the infimum are equal and $S(t+s) = S(t) \circ S(s)$. Now consider $G_1$ and $G_2$ two functions of $C(\mathbb{R}^d)$ such that $G_1 \le G_2$ and let $t > 0$. Thanks to the continuity of $G_2$, $\exists \mathbf{y}_{\mathbf{x},t} \in B(\mathbf{x}, t)$ such that $S(t)G_2(\mathbf{x}) = G_2(\mathbf{y}_{\mathbf{x},t})$. Consequently $G_2(\mathbf{y}_{\mathbf{x},t}) \ge G_1(\mathbf{y}_{\mathbf{x},t}) \ge S(t)G_1(\mathbf{x})$, which concludes the proof. $\square$

Now let $\phi \in C^1(\mathbb{R}^d \times (0, +\infty))$ and suppose that $(\mathbf{x}, t)$ is a local maximum of $G - \phi$. Thanks to the semigroup property of $S$ we get that:

$$
G(\mathbf{x}, t) = S(t)G_0(\mathbf{x}) = S(h)S(t-h)G_0(\mathbf{x}) = S(h)G(\mathbf{x}, t-h).
$$

Therefore, for all $0 < h < t$, we have

$$
(27) \qquad G(\mathbf{x}, t) = \inf_{|\mathbf{x}-\mathbf{y}|\le h} G(\mathbf{y}, t-h).
$$

$(\mathbf{x}, t)$ being a local maximum of $G - \phi$, we have, if $h$ is sufficiently small, and
$|\mathbf{x} - \mathbf{y}| \le h$:

$$
G(\mathbf{y}, t-h) - \phi(\mathbf{y}, t-h) \leq G(\mathbf{x}, t) - \phi(\mathbf{x}, t),
$$

which is equivalent to

$$
G(\mathbf{y}, t-h) \leq G(\mathbf{x}, t) - \phi(\mathbf{x}, t) + \phi(\mathbf{y}, t-h).
$$

Injecting this in 27 leads to

$$
\phi(\mathbf{x}, t) \leq \inf_{|\mathbf{x}-\mathbf{y}|\leq h} \phi(\mathbf{y}, t-h).
$$

A first order Taylor expansion at the point $(\mathbf{x}, t)$ leads to

$$
0 \le \inf_{|\mathbf{x}-\mathbf{y}|\le h} \left[ -\partial_t \phi(\mathbf{x}, t) + \nabla \phi(\mathbf{x}, t) \cdot \frac{\mathbf{y}-\mathbf{x}}{h} + o(1) \right].
$$

Using that fact that $-\inf(-) = \sup()$, we have

$$
\partial_t \phi(\mathbf{x}, t) + \sup_{|\mathbf{x}-\mathbf{y}|\le h} \nabla \phi(\mathbf{x}, t) \cdot \frac{\mathbf{x}-\mathbf{y}}{h} + o(1) \le 0.
$$

Thanks to the Cauchy-Schwarz inequality:

$$
|\nabla \phi(\mathbf{x}, t) \cdot \frac{\mathbf{x} - \mathbf{y}}{h}| \leq |\nabla \phi(\mathbf{x}, t)|.
$$
By taking $y = x - \frac{\nabla\phi(x, t)}{|\nabla\phi(x, t)|}h$ we see that the previous upper-bound is reached.

Therefore,

$$\partial_t\phi(x, t) + |\nabla\phi(x, t)| + o(1) \le 0,$$

and passing to the limit when $h \to 0$ leads to the desired result.
REFERENCES

[1] R. Abgrall. Numerical discretization of the first-order Hamilton-Jacobi equation on triangular meshes. *Comm. Pure Appl. Math.*, 49(12):1339–1373, 1996.

[2] Steeve Augoula and Rémi Abgrall. High order numerical discretization for Hamilton-Jacobi equations on triangular meshes. *J. Sci. Comput.*, 15(2):197–229, 2000.

[3] G. Barles and P. E. Souganidis. Convergence of approximation schemes for fully nonlinear second order equations. *Asymptotic Anal.*, 4(3):271–283, 1991.

[4] Guy Barles. *Solutions de viscosité des équations de Hamilton-Jacobi*, volume 17 of *Mathématiques & Applications (Berlin) [Mathematics & Applications]*. Springer-Verlag, Paris, 1994.

[5] Timothy J. Barth and James A. Sethian. Numerical schemes for the Hamilton-Jacobi and level set equations on triangulated domains. *J. Comput. Phys.*, 145(1):1–40, 1998.

[6] Steve Bryson and Doron Levy. High-order central WENO schemes for multidimensional Hamilton-Jacobi equations. *SIAM J. Numer. Anal.*, 41(4):1339–1369 (electronic), 2003.

[7] P. G. Ciarlet. Basic error estimates for elliptic problems. In P. Ciarlet and J.L. Lions, editors, *Handbook of Numerical Analysis, Volume II*, pages 17–351. North Holland, 1991.

[8] M. G. Crandall and P.-L. Lions. Two approximations of solutions of Hamilton-Jacobi equations. *Math. Comp.*, 43(167):1–19, 1984.

[9] Michael G. Crandall and Pierre-Louis Lions. Viscosity solutions of Hamilton-Jacobi equations. *Trans. Amer. Math. Soc.*, 277(1):1–42, 1983.

[10] D. Grapsas, R. Herbin, W. Kheriji, and J.-C. Latché. An unconditionally stable staggered pressure correction scheme for the compressible Navier-Stokes equations. *under revision*, 2015.

[11] R. Herbin, W. Kheriji, and J.-C. Latché. On some implicit and semi-implicit staggered schemes for the shallow water and Euler equations. *ESAIM Math. Model. Numer. Anal.*, 48(6):1807–1857, 2014.

[12] R. Herbin, J.-C. Latché, and T. T. Nguyen. Explicit staggered schemes for the compressible Euler equations. In *Applied mathematics in Savoie—AMIS 2012: Multiphase flow in industrial and environmental engineering*, volume 40 of *ESAIM Proc.*, pages 83–102. EDP Sci., Les Ulis, 2013.

[13] G. Kossioris, Ch. Makridakis, and P. E. Souganidis. Finite volume schemes for Hamilton-Jacobi equations. *Numer. Math.*, 83(3):427–442, 1999.

[14] Pierre-Louis Lions. *Generalized solutions of Hamilton-Jacobi equations*, volume 69 of *Research Notes in Mathematics*. Pitman (Advanced Publishing Program), Boston, Mass.-London, 1982.

[15] Stanley Osher and James A. Sethian. Fronts propagating with curvature-dependent speed: algorithms based on Hamilton-Jacobi formulations. *J. Comput. Phys.*, 79(1):12–49, 1988.

[16] Stanley Osher and Chi-Wang Shu. High-order essentially nonoscillatory schemes for Hamilton-Jacobi equations. *SIAM J. Numer. Anal.*, 28(4):907–922, 1991.

[17] P²REMICS. Computational fluid dynamics software for the simulation of dispersion and explosion. https://gforge.irsn.fr/gf/project/p2remics.

[18] L. Piar, F. Babik, R. Herbin, and J.-C. Latché. A formally second-order cell centred scheme for convection-diffusion equations on general grids. *Internat. J. Numer. Methods Fluids*, 71(7):873–890, 2013.

[19] Susana Serna and Jianliang Qian. Fifth-order weighted power-ENO schemes for Hamilton-Jacobi equations. *J. Sci. Comput.*, 29(1):57–81, 2006.

[20] J. A. Sethian and A. Vladimirsky. Fast methods for the eikonal and related Hamilton-Jacobi equations on unstructured meshes. *Proc. Natl. Acad. Sci. USA*, 97(11):5699–5703, 2000.

[21] Panagiotis E. Souganidis. Approximation schemes for viscosity solutions of Hamilton-Jacobi equations. *J. Differential Equations*, 59(1):1–43, 1985.

[22] Jue Yan and Stanley Osher. A local discontinuous Galerkin method for directly solving Hamilton-Jacobi equations. *J. Comput. Phys.*, 230(1):232–244, 2011.

UNIVERSITÉ DE NANTES, LABORATOIRE DE MATHÉMATIQUES JEAN LERAY, CNRS UMR 6629,
2 RUE DE LA HOUSSINIÈRE, BP 92208, 44322 NANTES, FRANCE
*E-mail address: nicolas.therme@univ-nantes.fr*