A stabilized DG cut cell method for discretizing the linear transport equation

TL;DR

This paper introduces a new stabilization technique for the discontinuous Galerkin method applied to the linear transport equation on cut cell meshes, enabling explicit time stepping without stability issues.

Contribution

The authors develop and analyze stabilization terms that allow explicit time stepping on cut cell meshes in DG discretizations of the linear transport equation.

Findings

Monotonicity proven for piecewise constant polynomials in 1D

Total variation diminishing stability shown for piecewise linear polynomials

Numerical results confirm theoretical stability and accuracy

Abstract

We present new stabilization terms for solving the linear transport equation on a cut cell mesh using the discontinuous Galerkin (DG) method in two dimensions with piecewise linear polynomials. The goal is to allow for explicit time stepping schemes, despite the presence of cut cells. Using a method of lines approach, we start with a standard upwind DG discretization for the background mesh and add penalty terms that stabilize the solution on small cut cells in a conservative way. Then, one can use explicit time stepping, even on cut cells, with a time step length that is appropriate for the background mesh. In one dimension, we show monotonicity of the proposed scheme for piecewise constant polynomials and total variation diminishing in the means stability for piecewise linear polynomials. We also present numerical results in one and two dimensions that support our theoretical findings.

Tables

Table 1. Table 1: Convergence rates ( L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} norm) for Test 4 for constant ( β C superscript 𝛽 𝐶 \beta^{C} ) and varying ( β V superscript 𝛽 𝑉 \beta^{V} ) velocity field.

Angle $\gamma$:		5∘	10∘	15∘	20∘	25∘	30∘	35∘	40∘	45∘
${\beta }^C$	$L^1$	2.02	2.01	2.01	2.00	2.00	2.01	2.01	2.01	2.02
	$L^{\infty }$	1.88	1.68	1.63	1.60	1.60	1.56	1.54	1.53	1.58
${\beta }^V$	$L^1$	2.02	2.01	2.01	2.01	2.01	2.01	2.01	2.01	2.00
	$L^{\infty }$	1.62	1.60	1.57	1.55	1.51	1.55	1.54	1.53	1.59

Equations

$$u_t+⟨\beta ,\nabla u⟩$$

$$u$$

$${\mathcal{T}}_h:=\{E:=\hat{E}\cap \Omega \left|\hat{E}\in \widehat{{\mathcal{T}}_h}\right\}.$$

$${\Gamma }_{\text{int}}$$

$${\Gamma }_{\text{ext}}$$

$$V_h^k({\mathcal{T}}_h)=\left\{v_h\in L^2(\Omega )|\forall E\in {\mathcal{T}}_h,v_h{|}_E\in P^k(E)\right\},$$

$$\displaystyle\left\llbracket u_{h}\right\rrbracket:={u_{h}}|_{{}_{E_{1}}}n_{E_{1}}+{u_{h}}|_{{}_{E_{2}}}n_{E_{2}},$$

$$\displaystyle\left[u_{h}\right]_{E}:=-\left\llbracket u_{h}\right\rrbracket\cdot n_{E}.$$

$$\left\{\{u_h\right\}\}:=\frac{1}{2}(u_h{|}_{{}_{E_1}}+u_h{|}_{{}_{E_2}}).$$

$$⟨d_t{u}_h(t),w_h⟩+a_h^{\text{upw}}\left(u_h(t),w_h\right)+l_h\left(w_h\right)=0,\forall w_h\in V_h^k({\mathcal{T}}_h),$$

$$\displaystyle\begin{aligned} a_{h}^{\text{upw}}(u_{h},w_{h}):=&-\sum_{E\in\mathcal{T}_{h}}\int_{E}u_{h}\left\langle\beta,\nabla_{h}w_{h}\right\rangle\differential{x}+\sum_{e\in\Gamma_{\text{ext}}}\int_{e}\left\langle\beta,n\right\rangle^{\oplus}u_{h}w_{h}\differential{s}\\ &+\sum_{e\in\Gamma_{\text{int}}}\int_{e}\left(\left\{\mskip-5.0mu\{u_{h}\right\}\mskip-5.0mu\}\left\langle\beta,\left\llbracket w_{h}\right\rrbracket\right\rangle+\frac{1}{2}\left|\left\langle\beta,n_{e}\right\rangle\right|\left\langle\left\llbracket u_{h}\right\rrbracket,\left\llbracket w_{h}\right\rrbracket\right\rangle\right)\differential{s},\end{aligned}$$

$$\begin{array}{rl}{l}_h(w_h):=& -\sum _{e\in {\Gamma }_{\text{ext}}}{\int }_e{⟨\beta ,n⟩}^{\ominus }g{w}_h\text{differential}s,\end{array}$$

$$\begin{array}{rl}{y}^{(1)}& =y^n+\Delta tF(y^n),\\ y^{n+1}& =\frac{1}{2}{y}^n+\frac{1}{2}{y}^{(1)}+\frac{1}{2}\Delta tF(y^{(1)}).\end{array}$$

$${\alpha }_{E,\omega }:=\min \left(\omega \frac{\text{absolutevalue}E}{\Delta t{\int }_{\partial E}{⟨\beta ,n⟩}^{\ominus }\text{differential}s},1\right),\omega \in \mathbb{R},\omega \in (0,1].$$

$$\mathcal{I}=\left\{E\in {\mathcal{T}}_h|E\text{is small cut cell and will be stabilized}\right\}.$$

$${\mathcal{F}}_i(E)$$

$${\mathcal{F}}_o(E)$$

$${\mathcal{N}}_o(E)=\{E^{\prime }\in {\mathcal{T}}_h\left|\partial E\cap \partial E^{\prime }\in {\mathcal{F}}_o(E),E^{\prime }\ne E\right\}.$$

$$\tau (x_0):=\{x\in E\left|\frac{d}{dt}x(t)=\beta (x(t)),x(0)=x_0,t>0\right\}.$$

$$T_E(w)(x):=w({\tau }^{-1}(x)),$$

$$\displaystyle\begin{split}J^{0}_{h}(u_{h},w_{h})=&\sum_{E\in\mathcal{I}}J^{0,E}_{h}(u_{h},w_{h}),\text{\qquad with}\\ J^{0,E}_{h}(u_{h},w_{h})=&-\int_{E}\eta_{E}\,T_{E}(\left[u_{h}\right]_{E})\left\langle\beta,\nabla w_{h}\right\rangle\differential{x}\\ &+\sum_{e\in\mathcal{F}_{o}(E)}\int_{e}\eta_{E}\,T_{E}(\left[u_{h}\right]_{E})\left\langle\beta,\left\llbracket w_{h}\right\rrbracket\right\rangle\differential{s},\end{split}$$

$$\text{and with }{\eta }_E=1-{\alpha }_{E,\frac{1}{2}}.$$

$$\displaystyle\begin{split}J^{1}_{h}(u_{h},w_{h})=&\sum_{E\in\mathcal{I}}J^{1,E}_{h}(u_{h},w_{h}),\text{\qquad with}\\ J^{1,E}_{h}(u_{h},w_{h})=&-\rho\int_{E}\eta_{E}\,\operatorname{len}(T_{E})\,T_{E}(\left[\partial_{\tau}u_{h}\right]_{E})\left\langle\beta,\nabla w_{h}\right\rangle\differential{x}\\ &+\sum_{e\in\mathcal{F}_{o}(E)}\int_{e}\eta_{E}\,\operatorname{len}(T_{E})\,T_{E}(\left[\partial_{\tau}u_{h}\right]_{E})\left\langle\beta,\left\llbracket w_{h}\right\rrbracket\right\rangle\differential{s},\end{split}$$

$$⟨d_t{u}_h(t),w_h⟩+a_h^{\text{upw}}\left(u_h(t),w_h\right)+J_h^0(u_h,w_h)+J_h^1(u_h,w_h)+l_h\left(w_h\right)=0.$$

$${\mathbb{I}}_E(x)=\{\begin{array}{ll}1& \text{for }x\in E,\\ 0& \text{otherwise}.\end{array}$$

$$J_h^{0,E}(u_h,{\mathbb{I}}_{E\cup {\mathcal{N}}_o(E)})+J_h^{1,E}(u_h,{\mathbb{I}}_{E\cup {\mathcal{N}}_o(E)})=0.$$

$$J_h^{0,E}(u_h,{\mathbb{I}}_{E\cup {\mathcal{N}}_o(E)})=J_h^{0,E}(u_h,{\mathbb{I}}_E)+J_h^{0,E}(u_h,{\mathbb{I}}_{{\mathcal{N}}_o(E)})$$

$$=\sum _{e\in {\mathcal{F}}_o(E)}{\int }_e{\eta }_E{T}_E({\left[u_h\right]}_E)⟨\beta ,{\mathbb{I}}_E{n}_E⟩\text{differential}s+\sum _{e\in {\mathcal{F}}_o(E)}{\int }_e{\eta }_E{T}_E({\left[u_h\right]}_E)⟨\beta ,-{\mathbb{I}}_{{\mathcal{N}}_o(E)}{n}_E⟩\text{differential}s$$

$$=0.$$

$$u_t(x,t)+\beta u_x(x,t)=0\text{ in }I\times (0,T),u(0,t)=g(t)\text{ for }t\in (0,T),$$

$$\displaystyle\begin{split}a_{h}^{\text{upw}}(u_{h},w_{h})=&-\sum_{j=1}^{N}\int_{j}\beta u_{h}\partial_{x}w_{h}\differential{x}+\beta u_{h}(x_{N+\frac{1}{2}}^{-})w_{h}(x_{N+\frac{1}{2}}^{-})\\ &+\sum_{j=1}^{N-1}\beta u_{h}(x_{j+\frac{1}{2}}^{-})\left\llbracket w_{h}\right\rrbracket_{j+\frac{1}{2}}.\end{split}$$

$$\displaystyle\left\llbracket w_{h}\right\rrbracket_{j+\frac{1}{2}}=w_{h}(x_{j+\frac{1}{2}}^{-})-w_{h}(x_{j+\frac{1}{2}}^{+}),\quad l_{h}(w_{h}):=-\beta g(0)w_{h}(x_{\frac{1}{2}}^{+}).$$

$$\lambda =\frac{\beta \Delta t}{h}.$$

Full text

1 Introduction

Finite element (FE) and more recently discontinuous Galerkin (DG) schemes have successfully been used on a huge variety of equations and are overall well understood. Therefore, research has advanced from solving model problems on simple geometries to trying to solve real-world problems. As a result, grid generation has become a huge issue. Simulating flow around an airplane, flow in blood vessels, or phase transitions requires to mesh very complicated geometries, which are often given as CAD models or implicitly. The generation of corresponding body-fitted meshes is a very involved process.

In recent years, the usage of embedded boundary or cut cell meshes has become increasingly popular. The details of these approaches vary. In this work, we will focus on the approach of cutting a given geometry out of given background mesh, resulting in so called cut cells along the boundary of the embedded object. These cells can have various shapes and may in particular become arbitrarily small. Special schemes must be developed to guarantee stability on these cells. There already exists a significant amount of literature for stabilizing problems of elliptic and parabolic type on cut cell meshes, see e.g. [11, 24, 18, 13, 2, 6, 1, 38, 29] and the references cited therein. For small cells, stability of higher derivatives is lost and different stabilization techniques have been proposed. A successful approach is the ghost penalty stabilization[12], sometimes referred to as the cutFEM method [13].

For hyperbolic conservation laws on cut cell meshes, different problems arise, compared to solving elliptic and parabolic PDEs. Probably the biggest issue is the small cell problem – that standard explicit schemes are not stable on the arbitrarily small cut cells when the time step is chosen according to the cell size of the background mesh. An additional complication is the fact that there is typically no concept of coercivity that could serve as a guideline for constructing stabilization terms. Furthermore one wants the numerical scheme to satisfy properties such as monotonicity and TVD (total variation diminishing) stability in order to avoid overshoot or oscillations in the presence of discontinuities, which could result in unphysical solutions.

In the context of finite volume schemes, which have traditionally been used for the solution of hyperbolic conservation laws, already several solution approaches exist that solve the small cell problem while dealing with arbitrarily small cut cells; e.g., the flux redistribution method [15, 36, 17], the h-box method [10, 25, 9], and the mixed explicit implicit scheme [33]. For the solution of hyperbolic conservation laws on cut cell meshes using DG schemes, very little work has been done. As one of the first ones, in Bastian et. al [7] use an implicit Euler method to overcome the small cell problem for a linear transport problem, but this approach is bound to first order accuracy. For explicit time stepping schemes some work relies on so called cell merging or cell agglomeration [35, 28]. In this approach, cut cells that are too small are merged with neighboring cells. As this approach puts the complexity back into the mesh generation we do not want to consider it here. Recently, Gürkan and Massing [23] suggested a new scheme for solving the steady advection-reaction problem on a cut cell mesh with potentially arbitrarily small cut cells that uses a ghost penalty approach. However, the authors consider only the steady case and do not discuss properties such as monotonicity or TVD stability. Also, Sticko and Kreiss [40] have worked on using penalty terms to stabilize the solution of the wave equation. In both cases, the penalty term employed has great similarity to the ghost penalty stabilization used for elliptic problems [12].

To the best of our knowledge, the scheme suggested in this work is the first DG scheme for overcoming the small cell problem for time-dependent scalar conservation laws by dealing with the potentially arbitrarily small cut cells while ensuring monotonicity and TVD stability. Therefore, we will focus on the linear advection equation as the standard model problem in the following. Many problems that occur for solving more complex hyperbolic equations already show up for this simple model.

Our approach is based on stabilizing the spatial discretization. One is free in the choice of the time stepping scheme. In particular, our stabilized spatial discretization allows for using standard explicit time stepping schemes everywhere, even on cut cells. We obtain stability on cut cells by adding penalty terms. In that sense the suggested scheme is similar to the ghost penalty approach [12]. However, the details of the terms differ significantly. In particular, the terms are fundamentally different from the ones used in [23, 40]. Conceptually, the terms are designed to reconstruct the proper domain of dependence on small cut cells and their neighbors, similar to the idea of the h-box method, but without an explicit geometry reconstruction. In this work we consider piecewise linear polynomials in one and two space dimensions.

This paper is structured as follows. In section 2.1, we provide background information, such as triangulation and geometry information as well as the standard scheme that we plan to stabilize. Section 3 contains the core of this work, the formulation of the stabilization terms in 2D. In section 4, we focus on the 1D case for a better understanding of the proposed stabilization and provide theoretical results for the case of piecewise constant and piecewise linear polynomials. We discuss implementational aspects of the proposed scheme in section 5. And in section 6, we provide numerical results in both 1D and 2D that support our theoretical findings. We conclude with a short summary and an outlook in section 7.

2 Discretization

2.1 Preliminaries

In this work, we focus on the linear advection equation as the classic model problem for hyperbolic scalar conservation laws, which is given by

[TABLE]

Here, $\Omega\subset\mathbb{R}^{2}$ is a open, connected, polygonal domain, $\partial\Omega$ denotes its boundary, and $\partial\Omega^{\text{in}}:=\{x\in\partial\Omega:\left\langle\beta(x),n(x)\right\rangle<0\}$ its inflow boundary. The velocity field $\beta\in\mathbb{R}^{2}$ is assumed to satisfy $\nabla\cdot\beta=0$ and $\left\langle\cdot,\cdot\right\rangle$ denotes the canonical scalar product in $\mathbb{R}^{2}$.

To construct our discrete approximation, we first consider a larger polygonal domain $\widehat{\Omega}\supset\Omega$, which is easy to mesh, see figure 1.

Definition 2.1** (Cut cell mesh).**

Given a background mesh $\widehat{\mathcal{T}_{h}}$ of $\widehat{\Omega}$, we introduce a cut cell mesh $\mathcal{T}_{h}(\Omega)$. Let $\widehat{\mathcal{T}_{h}}$ be a non-overlapping set of shape-regular elements $\widehat{E}$ such that $\bigcup\limits_{\widehat{E}\in\widehat{\mathcal{T}_{h}}}\overline{\widehat{E}}=\overline{\widehat{\Omega}}$. For simplicity, we will assume in this paper that $\widehat{\mathcal{T}_{h}}$ corresponds to a Cartesian background mesh but this is not necessary. Intersecting $\Omega$ and the background mesh induces the triangulation

[TABLE]

*Note that $\mathcal{T}_{h}$ is a triangulation of $\Omega$ consisting of structured (Cartesian) cells and cut cells. The internal and external skeletons of the partitioning are given by *

[TABLE]

Definition 2.2** (Discrete Function Space).**

Following the usual discontinuous Galerkin formulation we define the discrete space $V_{h}^{k}(\mathcal{T}_{h})\subset L^{2}(\Omega)$ by

[TABLE]

where $P^{k}$ denotes the polynomial space of degree $k$.

In this paper we will only consider the cases $k=0$ and $k=1$, i.e., piecewise constant and piecewise linear polynomials. On the skeleton $\Gamma_{\text{int}}$, functions $u_{h}\in V^{k}_{h}$ are not well-defined. Therefore, we define jump and average as follows.

Definition 2.3** (Jump & Average).**

The jump in normal direction over a face $e_{E_{1},E_{2}}=\partial{E_{1}}\cap\partial{E_{2}}$ between two elements $E_{1}$ and $E_{2}$ is defined as

[TABLE]

where $n_{E_{1}},n_{E_{2}}\in\mathbb{R}^{2}$ denote the unit outward normals of $E_{1}$ and $E_{2}$, respectively. Note that the jump of a scalar quantity is vector-valued, i.e., \left\llbracket u_{h}\right\rrbracket\in\mathbb{R}^{2}. We define a scalar-valued jump on a face $e\in\partial E$ as

[TABLE]

The (scalar-valued) average on a face is defined as

[TABLE]

2.2 Unstabilized DG Formulation

We now introduce the DG scheme as used on structured background cells. Additional stabilization terms will be necessary on small cut cells, as we will discuss in section 3. We will use a method of lines approach and first discretize (1) in space and then in time.

The spatial discretization of (1) corresponds to the standard DG discretization using upwind fluxes [19]: Find $u_{h}\in V^{k}_{h}(\mathcal{T}_{h})$ such that

[TABLE]

with

[TABLE]

where $n\in\mathbb{R}^{2}$ denotes the unit outer normal on $\partial\Omega$ and $n_{e}\in\mathbb{R}^{2}$ denotes a unit normal on a face $e$ (of arbitrary but fixed orientation). The negative and positive components of a quantity $x\in\mathbb{R}$ are defined as $x^{\ominus}:=\frac{\left|x\right|-x}{2}$ and $x^{\oplus}:=\frac{\left|x\right|+x}{2}.$ Note that $x^{\ominus},x^{\oplus}\geq 0$.

The new method does not rely on a particular time stepping scheme. Nevertheless, it is desirable that the scheme is explicit in order to allow for fast operator evaluation and to ease the use of limiters. Furthermore, it should be of the same order of accuracy as the spatial discretization and result in a discretization that is TVD.

While for piecewise constant polynomials in space the explicit Euler scheme suffices, it will diminish the convergence order for $V_{h}^{1}(\mathcal{T}_{h})$. We thus decided to employ the explicit second-order TVD Runge-Kutta (RK) scheme [22] that is given for the ODE $y_{t}=F(y)$ by

[TABLE]

A limiter is necessary to avoid unphysical oscillations close to discontinuities when using piecewise linear polynomials. We have chosen a Barth-Jespersen type limiter [3] that has been extended to the DG setting by exploiting the structure of the local Taylor basis [30]. This limits the gradient in such a way that the local solution of a cell $E$ evaluated at each neighboring centroid does not over/undershoot the maximum/minimum taken over the cell’s $E$ average value and the average values of all of cell’s $E$ face neighbors. The limiter is applied as a postprocessing step to both the intermediate solution $u^{(1)}$ and the solution $u^{n+1}$.

3 Stabilization

There are two problems that need to be dealt with in order to ensure stability on cut cells and to avoid overshoot and oscillations. First, the average mass in each cell must be controlled so that a piecewise constant approximation does not produce oscillations. Second, the gradients must be controlled to avoid oscillations within a single cell and unphysical evaluations at cell faces. Both goals must be reached without violating mass conservation. For this purpose, we introduce two stabilization terms $J^{0}_{h}$ and $J^{1}_{h}$: the first one will control average values and the second one mainly gradients.

One way to think of the stabilization terms is that they ensure that only a certain fraction of the inflow stays in the small cut cell by transporting the remaining part of the cut cell’s inflow directly through the small cut cell into its outflow neighbors.

Definition 3.1** (Capacity).**

We define the capacity of a cut cell $E$ as

[TABLE]

For $\omega=1$, the capacity measures the fraction of the inflow that is allowed to flow into the cut cell $E$ without producing overshoot.

We further denote with $\mathcal{I}$ the set of cut cells on which the stabilization should be applied. The idea is that only small cut cells need stabilization:

[TABLE]

We assume that for every pair of neighboring cut cells, at most one of the two cut cells is an element of $\mathcal{I}$. While it is possible to extend the proposed scheme to several neighboring cut cells, the implementation would be more difficult. We design the stabilization such that the CFL condition for explicit time stepping schemes essentially only depends on the resolution of the background mesh and not on the size of small cut cells from the set $\mathcal{I}$.

Definition 3.2** (Inflow and outflow faces).**

For each cell $E$ we denote the set of inflow faces and outflow faces as

[TABLE]

where $n_{E}$ denotes the unit outer normal of cell $E$ on face $e$.

Definition 3.3** (Outflow neighbors of $E\in\mathcal{I}$).**

The set of outflow neighbors $\mathcal{N}_{o}(E)$ of a cut cell $E\in\mathcal{I}$ is defined by

[TABLE]

Another way to think of the stabilization is that it reconstructs the proper domain of dependence, i.e., it ensures that the outflow neighbors of a small cut cell get information from the inflow neighbors of that cut cell. The stabilization term at a given point on an outflow face thus depends on the flux on the inflow faces, measured upstream along the trajectory.

Definition 3.4** (Trajectory operator).**

For every point $x_{0}\in\mathcal{F}_{i}(E)$ the trajectory $\tau:\mathcal{F}_{i}(E)\to E$ describes the curve traced out by a particle starting at $x_{0}$ and being transported with $\beta$:

[TABLE]

Note that $\tau$ describes an injective mapping. We define $\tau^{-1}$ as the operator that maps back onto $\mathcal{F}_{i}(E)$ and introduce the trajectory operator

[TABLE]

as the evaluation of $w$ on the inflow face, following the trajectory backwards from $x$. An illustration of the trajectory operator is given in figure 2.

Definition 3.5** (Trajectory length).**

*The trajectory length $\operatorname{len}(T_{E})(x)$ measures the length of a trajectory within the cell $E$. It is defined as follows: for a point $x\in E$ we identify the point $x_{0}$ on the inflow face and consider the length of the curve from $x_{0}$ (through $x$) to a point on an outflow face of $E$. *

In section 5.1 we will discuss how to realize $T_{E}$ and $\operatorname{len}(T_{E})$ practically.

The penalty term $J_{h}^{0}$ is now constructed in such a way that any inflow that exceeds the capacity of a small cut cell is moved to the downwind cells and is given by

[TABLE]

The volume contribution can be seen as a transport within the cut cell of the quantity $\eta_{E}=1-\alpha_{E,\frac{1}{2}}$, which should not remain in the cut cell, from its inflow faces to its outflow faces. The second term which is applied on the outflow faces transports this quantity out of the cut cell into its downwind neighbors. We only allow the fraction $\alpha_{E,\frac{1}{2}}$ instead of $\alpha_{E,1}$ to stay within the cut cell as the latter one would result in too restrictive slope limiting, leading to close to zero gradients on the cut cells. The main task of the second stabilization term $J^{1}_{h}$ is to restore control over gradients. Its structure is similar to the one of $J_{h}^{0}$, but it employs the derivative of the discrete solution $u_{h}$ and uses a different scaling and a different relation between skeleton and boundary terms:

[TABLE]

where $\partial_{\tau}u_{h}$ denotes the derivative along the trajectory, which is evaluated $\partial_{\tau}u_{h}=\left\langle\nabla u_{h},\beta\right\rangle/\norm{\beta}$, and $\rho=\frac{1}{2}$. (We introduce the parameter $\rho$ here as we will later examine the stability of the resulting scheme for different values of $\rho$.) The stabilized spatial discretization is then given by: Find $u_{h}\in V^{k}_{h}(\mathcal{T}_{h})$ such that $\forall\>w_{h}\in V_{h}^{k}(\mathcal{T}_{h})$

[TABLE]

Next, we will examine the mass conservation properties of the stabilization terms.

Definition 3.6** (Indicator function).**

The indicator function of a cell or cell patch is defined as

[TABLE]

Lemma 3.1**.**

Let $E\in\mathcal{I}$. Then

[TABLE]

Proof.

We first focus on $J^{0,E}_{h}$. Due to the linearity of $J^{0,E}_{h}$ in $w_{h}$, there holds

[TABLE]

The same argumentation can be used for $J_{h}^{1,E}$. ∎

This result (together with the local mass conservation properties of the standard DG discretization (8)) guarantees local mass conservation in the extended control volume $E\cup\mathcal{N}_{o}(E)$ of a small cut cell $E\in\mathcal{I}$: the same amount of mass that is not allowed to stay in $E$ is redistributed to $E$’s outflow neighbors. This slightly extended setting of local mass conservation is a natural consequence of our stabilization.

4 The 1D case

For a better understanding of the proposed scheme and for the validation of some theoretical properties, we focus on the 1D case in this section. WLOG, we consider the interval $I=(0,1)$ and assume $\beta>0$ to be constant. Our PDE is given by

[TABLE]

with initial data $u(x,0)=u_{0}(x)$. For the analysis in this section we will focus on solving the model problem MP shown in figure 3: we discretize the interval $I$ in $N-1$ cells $E_{j}=[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}]$ of equidistant length $h$. Then we split one cell, the cell $k$, in two cells of lengths $\alpha h$ (referred to as cell $k_{1}$) and $(1-\alpha)h$ (referred to as cell $k_{2}$) with $\alpha\in(0,\frac{1}{2}]$.

Remark 4.1** (Notation).**

*Different to 2D, there is a natural order of cells in 1D. In 1D, we will therefore use the notation indicated in figure 3 and refer to ‘cell $j$’ and associate faces with $x_{j\pm 1/2}$ instead of using ‘cell $E$’ or ‘cell $E_{j}$’ or ‘face $e$’. *

Using the notation $x_{j+\frac{1}{2}}^{\pm}=\lim_{\varepsilon\to 0}x_{j+\frac{1}{2}}\pm\varepsilon$, the bilinear form (8) simplifies to

[TABLE]

The definitions of the jump and of $l_{h}$, given by (5) and (9) in 2D, reduce in 1D to

[TABLE]

Further, we define the CFL number

[TABLE]

Note that $\lambda$ is chosen only with respect to the background mesh width $h$ but the volume fraction $\alpha\in\left(0,\frac{1}{2}\right]$ in MP is allowed to become arbitrarily close to 0. For the considered model problem MP, the set of cells that need to be stabilized consists of $\mathcal{I}=\{k_{1}\}$ and therefore $J^{0}_{h}$ and $J^{1}_{h}$ coincide with $J^{0,k_{1}}_{h}$ and $J^{1,k_{1}}_{h}$, respectively. The stabilization terms $J_{h}^{0,k_{1}}$ and $J_{h}^{1,k_{1}}$, given by (14) and (16) in 2D, simplify significantly for the considered setup in 1D:

• •

For cell $k_{1}$, the inflow face is $k-\frac{1}{2}$ and the outflow face is $k_{\text{cut}}$.

• •

For $x=x_{\text{cut}}$, $T_{E}(\left[u_{h}\right]_{E})(x)$ in 2D simply corresponds to \left\llbracket u_{h}\right\rrbracket_{k-\frac{1}{2}} in 1D: the trajectory operator $T_{E}$ is trivial and $\left[u_{h}\right]_{E}$, defined in (6), on an inflow face corresponds to the definition of \left\llbracket u_{h}\right\rrbracket in 1D, given by (22), for $\beta>0$.

• •

The length $\operatorname{len}(T_{E})$ is simply $\alpha h$.

• •

The derivative $\partial_{\tau}u_{h}$ in the definition of $J^{1,k_{1}}_{h}$ simply corresponds to $\partial_{x}u_{h}$.

Summarizing $J_{h}=J_{h}^{0}+J_{h}^{1}$, the stabilization term in 1D for the considered setup is of the form

[TABLE]

with

[TABLE]

Note that $\alpha_{k_{1},\omega}$ is exactly the 1D equivalent of the capacity $\alpha_{E,\omega}$ defined in (11) as using $\lambda=\frac{\beta\Delta t}{h}$ we can reformulate

[TABLE]

The resulting stabilized scheme is then given by: Find $u_{h}\in V_{h}^{k}(\mathcal{T}_{h})$ such that

[TABLE]

Definition 4.1** (MP0 and MP1).**

Both MP0 and MP1 refer to

• •

solving the advection equation (20) using periodic boundary conditions for the model problem MP shown in figure 3,

• •

using the stabilized scheme (26) with CFL $\lambda=\frac{\beta\Delta t}{h}$.

By MP0 we refer to the case of piecewise constant polynomials and by MP1 to the case of piecewise linear polynomials, respectively.

4.1 Behavior of the space discretization without stabilization term $J_{h}$

Before we discuss the properties of our stabilized scheme, we shortly present numerical results that illuminate the different kind of instabilities that occur if one does not use the stabilization term $J_{h}$.

4.1.1 Test 1: 1D single small cut cell

We consider the model problem MP and place the cell $k$ such that $x_{k-1/2}=0.5$. We use discontinuous initial data

[TABLE]

We set $\beta=1$, $\alpha=0.001$, $\lambda=0.5$, and $h=0.05$, and use $V_{h}^{0}(\mathcal{T}_{h})$ as well as periodic boundary conditions.

We compare three scenarios: (1) we do not use stabilization, i.e., we apply (26) without the term $J_{h}$; (2) we apply the stabilization $J_{h}$; (3) instead of the $J_{h}$ suggested in this work, we use a straight-forward adaption of the ghost penalty stabilization [12] to stabilize the problem. The ghost penalty stabilization has been used very successfully to stabilize the solution of elliptic problems on cut cell meshes. A first attempt to transfer the stabilization to the situation considered here would result in using a stabilization term of the form

[TABLE]

(instead of the $J_{h}$ that we suggest). We use $\rho_{1}=\rho_{2}=\eta_{k_{1}}$ (as for $J_{h}$).

We rewrite the described spatial discretization as a system of ODEs of the form $\underline{y}_{t}=L\underline{y}$ with a suitable matrix $L$. In figure 4(a) we present the eigenvalue distribution of $\Delta tL$ where $\Delta t=\lambda h$. Most values lie in the stability region of the explicit Euler scheme. But if we do not stabilize or use the ghost penalty stabilization (28), there is one outlier eigenvalue (corresponding to the small cut cell), leading to stability problems. With our stabilization, all values are in the stability region of explicit Euler.

Next, we consider the usage of an implicit time stepping scheme. The result after one time step for using the implicit Trapezoidal scheme is shown in figure 4(b). Despite the time stepping scheme being stable in an ODE stability sense, we observe a strong overshoot: instead of staying between 0 and 1, the value on the cut cell jumps up to 2 if we do not stabilize. When using the very simple ghost penalty stabilization (28) the overshoot is smaller but still significant. The reason for this behavior is that the implicit Trapezoidal scheme is not unconditionally TVD. In fact, no second- or higher-order scheme has this property [39, 27]. As a result, if the cell fraction $\alpha$ becomes too small, the scheme develops unphysical oscillations and/or overshoot unless one uses a stabilization that prohibits that behavior. Our stabilization has been designed to achieve this.

In the following, we will support these numerical results by mathematical facts: we will prove that our stabilization term $J_{h}$ (a) makes explicit time stepping stable again and (b) guarantees that no unphysical oscillations and/or overshoot occur.

4.2 Piecewise constant polynomials

In the special case of piecewise constant polynomials, the discrete solution $u_{h}$ on cell $j$ at time $t^{n}$, which we will denote by $u_{j}^{n}$, is an approximation to the average of $u$ over cell $j$ and the method is equivalent to a first-order finite volume scheme. For MP0, the stabilization term (24) reduces to (with $\eta_{k_{1}}$ given by (25))

[TABLE]

4.2.1 Monotonicity

One very desirable property of a first-order scheme for hyperbolic conservation laws is to be monotone. This property guarantees for example that an overshoot as seen in figure 4(b) cannot occur.

Definition 4.2** (see [26]).**

A method is called monotone, if

[TABLE]

Since this is challenging to verify we will use the following equivalent definition.

Definition 4.3** (see [41]).**

A method $u^{n+1}_{j}=H(u^{n}_{j-i_{L}},u^{n}_{j-i_{L}+1},...,u^{n}_{j+i_{R}})$ is called monotone, if $\>\forall j$ there holds for every $l$ with $-i_{L}\leq l\leq i_{R}$

[TABLE]

For a linear scheme this implies that all coefficients need to be non-negative. We will verify this property for our stabilized scheme using the concept of M-matrices.

Definition 4.4** (see [37]).**

Let $B\in\mathbb{R}^{n\times n}$ be a Z-Matrix, i.e., for $B=(b_{ij})$ there holds $b_{ij}\leq 0$ for every $i\neq j$. If $b_{ii}>0$ $\forall i$ and there exists a positive diagonal matrix $D$ such that $BD$ is strictly diagonal dominant, we call B an M-matrix.

Lemma 4.1** (see [37]).**

Let $B\in\mathbb{R}^{n\times n}$ be an M-matrix. Then $B^{-1}$ exists and $B^{-1}\geq 0.$

Lemma 4.2**.**

Consider MP0. Discretize the time using the theta scheme. Then, the method is monotone under the CFL condition $\lambda<\frac{1}{2(1-\Theta)}$.

Proof.

The theta scheme for the ODE $y_{t}=F(y)$ is defined as

[TABLE]

with $\Theta\in\left[0,1\right]$ ($\Theta=0:$ explicit Euler, $\Theta=1$: implicit Euler, $\Theta=\frac{1}{2}$: implicit Trapezoidal). Let us first assume $\alpha<2\lambda$, which is the interesting case. Then, $\alpha_{k_{1},\frac{1}{2}}=\frac{\alpha}{2\lambda}$. Rewriting one time step of our stabilized discretization in matrix-vector formulation results in

[TABLE]

Defining $\tau_{1}:=\Delta t\Theta\beta\geq 0$ and $\tau_{2}:=\Delta t\left(1-\Theta\right)\beta\geq 0$, the two matrices are of the form (entries on the main diagonal are boxed)

[TABLE]

and

[TABLE]

According to definition 31, we need to show that $B^{-1}C$ is non-negative. Following lemma 4.1 we will verify that $B$ is an M-matrix and $C$ is non-negative.

Most cells use a standard first-order upwind scheme in space and the theta scheme with standard CFL condition in time and obviously satisfy these conditions. We will focus on the rows corresponding to cells $k_{1}$ and $k_{2}$, which differ from the remaining rows due to the stabilization term $J_{h}$. We start with matrix $B$ and the row corresponding to cell $k_{2}$. The diagonal entry $(1-\alpha)h+\tau_{1}$ is positive while the remaining entries are negative, since $\alpha<2\lambda$. To prove the diagonal dominance we compute:

[TABLE]

Similar calculations for the row corresponding to $k_{1}$ imply that $B$ is an M-matrix.

Next we will show positivity of $C$, considering the CFL condition $\lambda<\frac{1}{2(1-\Theta)}$. For a standard, equidistant cell, the positivity of the entries is guaranteed already by the standard CFL condition and thus also by our slightly stricter condition. Straight-forward calculations further reveal that all entries in rows $k_{1}$ and $k_{2}$ are positive. We note that in particular the positivity of entry $c_{k_{2},k_{2}}$ is guaranteed by the CFL condition and the fact that a lower bound for the size of $k_{2}$ is $\frac{h}{2}$.

Finally, for $\alpha\geq 2\lambda$, there holds $\alpha_{k_{1},\frac{1}{2}}=1$ and therefore the factor $\eta_{k_{1}}=1-\alpha_{k_{1},\frac{1}{2}}$ in the definition of $J_{h}$ evaluates to 0, i.e., we would not stabilize. But as in this case $c_{k_{1},k_{1}}=\alpha h-\tau_{2}\geq\lambda h\Theta\geq 0$, all required properties of $B$ and $C$ would hold true without stabilization. This concludes the proof. ∎

Remark 4.2** (Ghost penalty stabilization).**

Using the stabilization (28) instead of the $J_{h}$ given in (29), this would result in additional $2\times 2$ blocks in the matrices $B$ and $C$ at positions $((k-1):k_{1},(k-1):k_{1})$ and $(k_{1}:k_{2},k_{1}:k_{2})$, respectively. We note that it is not possible to find factors $\rho_{1}$ and $\rho_{2}$ that guarantee that all entries of $C$ are non-negative for $\alpha\to 0$. In contrast, our stabilizing $2\times 2$ block has been shifted to $(k_{1}:k_{2},(k-1):k_{1})$. This shifted stabilization reflects the fact that for the hyperbolic equations there is a designated direction of information transport, whereas for elliptic problems there is not. This shifted stabilization is also a major difference to the stabilization terms used by other authors [23, 40].

4.2.2 $L^{1}$ stability and TVD stability

Next, we will show that MP0 with explicit Euler in time is $L^{1}$ stable as well as TVD stable (to be defined below) under the standard CFL condition, independent of the size of $\alpha$.

On standard cells away from the cut cells $k_{1}$ and $k_{2}$, the scheme given by (26) in combination with explicit Euler in time corresponds to the standard upwind scheme [31]. We will assume for the rest of this subsection that $\alpha<2\lambda$. If this was not the case, we would still have a non-uniform mesh, but would no longer violate the CFL condition on cell $k_{1}$ needed for $L^{1}$ and TVD stability and the results in the following will remain valid. For $\alpha<2\lambda$, we get the following formulae in the neighborhood of the small cut cell $k_{1}$ after a minor reordering and simplification of the terms

[TABLE]

Lemma 4.3**.**

Consider MP0 with explicit Euler in time. Then, for $\lambda<\frac{1}{2}$ there holds

[TABLE]

We note that the required CFL condition is independent of the size of $\alpha$ (but takes into account that the bigger cut cell $k_{2}$ is not stabilized).

Proof.

We define

[TABLE]

On all cells except for cells $k_{1}$ and $k_{2}$, the standard upwind scheme is used (compare also (34a) and (34d)). Plugging in the corresponding formulae and using $\lambda>0$ as well as $1-\lambda>0$ gives

[TABLE]

For cells $k_{1}$ and $k_{2}$, we directly obtain from equations (34b) and (34c) (using $1-\alpha-\lambda\geq 0$ and $2\lambda-\alpha\geq 0$)

[TABLE]

Summing up the estimates for $T_{1},\ldots,T_{4}$ implies the claim. ∎

Definition 4.5** ([16]).**

A DG scheme is called TVDM (total variation diminishing in the means) if

[TABLE]

*holds for all $n\geq 0$. Here, $\overline{u}_{j}^{n}$ denotes the mean of $u$ on cell $j$ at time $t_{n}$. *

For piecewise constant polynomials, the means $\overline{u}^{n}_{j}$ correspond to the unknowns $u^{n}_{j}$ and TVDM coincides with the TVD (total variation diminishing) [31] property. (For later use for piecewise linear polynomials we provide the more general definition here.)

Lemma 4.4**.**

Consider MP0 with explicit Euler in time. Then, the scheme is TVD stable for $\lambda<\frac{1}{2}$.

Proof.

We decompose

[TABLE]

For the standard parts of the scheme, we get

[TABLE]

Using (34a) and (34b), a direct substitution of the formulae results in

[TABLE]

For $T_{3}$ and $T_{4}$, we reorder the terms resulting from the definitions (34b)-(34d) to get

[TABLE]

We emphasize that all factors outside of the absolute values are non-negative due to the made assumptions. Summing up the estimates for $T_{1},\ldots,T_{5}$ implies the claim. ∎

Remark 4.3**.**

On an equidistant mesh, a standard finite volume scheme that is monotone is automatically TVD [26]. As the standard proof for this result does not transfer to our situation, we provide here the proof for TVD stability in addition.

4.3 Piecewise linear polynomials

In this subsection we examine MP1 more closely. In this case, the stabilization term $J_{h}$ is given by (24). For simplicity, we will assume in this subsection that we use a centered, rescaled moment basis, see also section 5.2. In time, we will consider the second-order explicit TVD RK scheme given by (10).

4.3.1 Eigenvalue analysis

We now want to motivate the stability of our stabilized scheme in combination with the second-order explicit RK scheme (10) by means of an eigenvalue analysis using symbolic computations with sympy[34]. We consider the model problem MP for the special case of $N=5$, i.e., we start with four equidistant cells of length $h$ and split the third cell in two cells of length $\alpha h$ and length $(1-\alpha)h$. We use Dirichlet boundary conditions and WLOG $\beta=1$.

We consider the stabilized scheme (26) with the stabilization term $J_{h}$ given by (24). For the construction of the stabilization term $J_{h}$, we have made several design choices, e.g., the choices of $\omega$ and $\eta_{k_{1}}$ and the general structure of the terms. Here, we will focus on the parameter $\rho$. We will show that $\rho=\frac{1}{2}$ is exactly the value that one needs to make explicit time stepping stable again.

We rewrite the variational formulation (26) as a system of ODEs with $M$ being the mass matrix, $A$ being the stiffness matrix incorporating $a_{h}^{\text{upw}}(\cdot,\cdot)$, and $J$ incorporating the corresponding parts of $J_{h}(\cdot,\cdot)$. We then symbolically setup the operator

[TABLE]

which is the stability function used in the ODE stability analysis. As $\alpha\leq\frac{1}{2}$, a lower bound for the size of the unstabilized cut cell $k_{2}$ is $\frac{h}{2}$. The limiting CFL number $\lambda$ for piecewise linear polynomials is then $\lambda=\frac{1}{6}$. We therefore choose $\Delta t=\frac{h}{6}$ to compute the eigenvalues of $R$. We find three pairs of complex eigenvalues:

[TABLE]

For a sequence of decreasing $\alpha$, we visualize the eigenvalues of $R$ and observe that for $\Delta t=\frac{h}{6}$ all eigenvalues stay within the stability region of the explicit second-order TVD RK scheme (10) (figure 5 left), indicating the stability of the explicit time stepping scheme. For fixing $\alpha=0.01$ and testing different choices of $\rho$, we observe (figure 5 right) that two eigenvalues diverge. These observations support our choice of $\rho=\frac{1}{2}$ and are consistent with our numerical experiments, in which we never observed stability issues due to using an explicit time stepping scheme on tiny cut cells.

4.3.2 TVDM stability

We can only expect to obtain a TVDM stability result for piecewise linear polynomials if we apply a limiter. Different to 2D, in 1D we only reconstruct to cell faces instead of to neighboring centroids, which results in an MC like limiter [31]. (In 2D, depending on the mesh, this could result in a limiter that is not linearity-preserving.) By assumption, our solution $u^{*}_{j}$ on cell $j$ resulting from solving (26) with the time stepping scheme (10), where $(\cdot)^{*}$ stands for $(\cdot)^{(1)}$ or $(\cdot)^{n+1}$, uses the centered moment basis

[TABLE]

Then, we enforce

[TABLE]

with

[TABLE]

Lemma 4.5**.**

Consider MP1 with explicit Euler in time. Assume that the moment basis (37) and limiter (38) is used, which has been modified on cell $k-1$ to additionally enforce

[TABLE]

Then, for $\lambda<\frac{1}{4}$ the scheme is TVDM stable.

Proof.

The proof follows that of lemma 4.4. Details are given in appendix A. ∎

Corollary 4.1**.**

Let the assumptions of lemma 4.5 hold true but replace the explicit Euler scheme by the second-order scheme TVD RK scheme (10). Then, for $\lambda<\frac{1}{4}$ the scheme is TVDM stable.

Proof.

The result follows directly from the fact that TVD RK schemes are constructed to be convex combinations of explicit Euler steps [22]. ∎

5 Implementation

The implementation is based on the DUNE [4, 5] framework, a feature rich C++ finite element library. For the cut cell discretization we rely on the dune-udg package[20] and its integration with dune-pdelab[8]. The dune-udg module was originally developed for the unfitted DG method[6]. We modified it to realize the extended stencil of the proposed stabilization terms $J^{0}_{h}$ and $J^{1}_{h}$.

The domain is represented as a discrete level set function that uses vertex values, i.e., as a bi-linear finite element function. Given this representation, cut cells and their corresponding quadrature rules are constructed using the TPMC libary[21]. We note that this setup results in a restriction of the geometry that can be realized: it only allows for a polygonal representation of the domain and for the numerical tests we cannot resolve smooth geometries exactly.

5.1 Implementing the trajectory operator & trajectory length

In 1D, the trajectory operator $T_{E}$ is trivial. In higher space dimensions it becomes more involved to evaluate $T_{E}$. We will now describe a simple approach to compute (approximately) the evaluation of $T_{E}(w)(x)$ for every point $x\in E$.

Let $E\in\mathcal{I}$. From the unsteady problem (1), we derive a localized stationary equation given by

[TABLE]

and natural outflow boundary conditions on $\mathcal{F}_{o}(E)$. The solution $\tilde{w}$ approximates $T_{E}(w)(x)\approx\tilde{w}(x)$. To discretize (40), we employ the bilinear form $a_{h}^{\text{upw}}$ restricted to $\Omega=E$, i.e., to a single element, which we denote by $a_{E}^{\text{upw}}$. We thus search $\tilde{w}_{h}\in P^{m}(E)$ such that

[TABLE]

Remark 5.1** (Test spaces).**

Depending on the shape of the streamlines of the vector field $\beta$ and the domain setup, we might choose $m\neq k$, i.e., use a different polynomial degree for this local subproblem. For $\beta|_{E}\in[P^{l}(E)]^{2}$ and $w|_{e}\in P^{k}(e),e\in\mathcal{F}_{i}(E)$, the solution $\tilde{w}$ is also polynomial and can thus be computed exactly. If we assume $\beta$ to depend linearly on $x$ and use $V_{h}^{0}(\mathcal{T}_{h})$ to solve (1), then $\left[u_{h}\right]_{e}\in P^{0}(e)$ and $\tilde{w}_{h}(x)=T_{E}(w_{h})(x)$ is exact for $\tilde{w}_{h}(x)\in P^{1}(E)$.

Remark 5.2** (Fast evaluation).**

Note that (40) is a linear problem. We can interpolate $w|_{\mathcal{F}_{i}}$ into a polynomial basis and also $\tilde{w}$ is represented in a local polynomial basis. We can therefore precompute local (per element) trajectory operators, mapping the coefficients $w_{i}$ onto coefficient $\tilde{w}_{j}$, and store this operator as a small dense matrix.

The trajectory length can be computed in a very similar way. It consists of two steps. First, we solve a local problem, which yields the trajectory length $\operatorname{len}(T_{E})(x)$ for every point $x\in\mathcal{F}_{i}(E)$. In the second step we compute $T_{E}(\operatorname{len}(T_{E}))$, which allows evaluating the local length of the corresponding trajectory for every point in $\bar{E}$.

To compute $\operatorname{len}(T_{E})(x)$ we consider the following local problem:

[TABLE]

and natural outflow boundary conditions on $\mathcal{F}_{i}(E)$. It basically computes the path integral over 1 along the trajectory and thus yields the trajectory length $\operatorname{len}(T_{E})(x)=\tilde{l}(x)$. Again we discretize using the weak form as derived for (1) and then we extend $\operatorname{len}(T_{E})(x)$ from $x\in\mathcal{F}_{i}(E)$ into the cell by computing $\tilde{w}$ in (40).

5.2 Construction of local shape-functions

Although in principal the choice of the local basis should not have impact on the method itself, there are practical aspects to consider. In particular the implementation of the limiter becomes significantly simpler, if one can separate the average value in the cell from the local fluctuations. Thus we choose to use a moment basis. In 1D we have chosen the local basis $\Phi_{k}$ in cell $k$ (in global coordinates) as

[TABLE]

This allows modifying the gradient without changing the average mass in the cell.

In 2D we don’t have an explicit formula, but we again construct the basis via moments in such a way, i.e., we split into a constant function and linear functions with average zero. The linear functions are centered around the center of mass of the cut cell, but in contrast to the 1D case, we don’t do the additional rescaling, so that the gradients are the same as on the background mesh.

6 Numerical results

In this section we support our theoretical findings by numerical experiments in 1D and 2D.

6.1 Numerical results for the 1D case

For the numerical results in 1D, we solve (26) with the stabilization term $J_{h}$ defined in (24) using $V_{h}^{1}(\mathcal{T}_{h})$.

6.1.1 Test 2: 1D, smooth initial data

We consider a variant of MP1 by changing the model problem MP slightly: we split all cells between $x=0.1$ and $x=0.9$ in cut cell pairs of length $\alpha_{k}h$ and $(1-\alpha_{k})h$ and compare 3 different scenarios:

S1: the fraction $\alpha_{k}$ is the same for all cells, i.e., $\alpha_{k}=\alpha$ (‘$\alpha=$’); 2. 2.

S2: the cell fraction $\alpha_{k}$ varies and is computed as $\alpha_{k}=0.1X_{k}+10^{-6}$ with $X_{k}$ being a uniformly distributed random number in $(0,1)$ (‘random $\alpha$’); 3. 3.

S3: we do not split in cut cell pairs (‘equidistant’).

We use smooth initial data $u_{0}(x)=\sin(2\pi x)$, set $\beta=1$, $\lambda=\frac{1}{6}$, and use the time stepping scheme (10).

In figure 6, the error at time $T=1$, measured in the $L^{1}$ and $L^{\infty}$ norm, are shown. We observe second-order convergence for all tested scenarios. In particular, the errors for pairs with $\alpha=0.1$ and $\alpha=10^{-4}$ are almost the same and fairly comparable to the case of not cutting the cells.

Figure 7 (left) shows the solution at time $T=1$ for S2 without and with using the limiter described in lemma 4.5 for a coarse mesh with $h=0.05$. Without limiter, the solution matches very well with the exact solution. With limiter, we observe the expected peak clipping but there is no form of staircasing or other unwanted interaction between the limiter and the gradient stabilization.

6.1.2 Test 3: 1D, discontinuous initial data

We use the setup of Test 2 with the scenario S2 but use discontinuous initial data $u_{0}$ given by (27). In figure 7 (right) we show the solution at time $T=1$ for $h=0.05$ with and without limiter. As expected, the solution does not show overshoot when the limiter is applied.

Additionally, we measure $\text{TV}(\overline{u}^{n})$, defined in (35), and observe that it is non-increasing during the simulation if the limiter is used. We note that this is not the case if we do not apply the additional condition (39) for limiting on the inflow neighbor of a small cut cell of type $k_{1}\in\mathcal{I}$.

6.2 Numerical results for the 2D case

For the numerical experiments in 2D we focus on the geometry of a ramp with angle $\gamma$, see figure 2(a): we use $\widehat{\Omega}=(0,1)^{2}$ and discretize it with $N\times N$ Cartesian cells, resulting in the mesh width $h=\frac{1}{N}$. Afterwards we cut out a ramp at $x=0.2001$ with angle $\gamma$ with a straight intersection, resulting in the computational domain $\Omega$. Our initial data $u_{0}$ are defined with respect to a standard Cartesian coordinate system $(x,y)$ and then transformed to our rotated and shifted coordinate system $(\hat{x},\hat{y})$ by using the transformation

[TABLE]

The cut cells are located along the ramp and have various shapes and sizes. In our experiments, the smallest volume fractions of cut cells $E$, computed as $\frac{\absolutevalue{E}}{h^{2}}$, varied typically between $10^{-5}$ and $10^{-8}$. Generally the set $\mathcal{I}$ should include all cells that are small in direction of $\beta$. In this particular setup, this is equivalent to choosing triangular cut cells, where the cut (i.e., the hypotenuse) is shorter than $\frac{h}{2}$. For ease of implementation we decided to simply use $\mathcal{I}=\{E\in\mathcal{T}_{h}\,|\,\frac{\absolutevalue{E}}{h^{2}}<0.1\}$. Using this definition of $\mathcal{I}$ for $5^{\circ}\leq\gamma\leq 45^{\circ}$, the shortest hypotenuse of a triangle cut cell that is not stabilized corresponds to roughly $0.6h$.

When choosing the time step size $\Delta t$, we need to take into account that the step size needs to be appropriate for $V_{h}^{k}(\mathcal{T}_{h}),k=0,1,$ for both Cartesian cells of size $h^{2}$ and cut cells $E$ that are not stabilized, i.e., for cut cells $E\not\in\mathcal{I}$. Finding a tight bound for the latter category is non-trivial and a project in itself. We have chosen the constraint

[TABLE]

which has worked well in our numerical experiments and which is in good agreement with a recently suggested CFL condition for DG schemes on triangular meshes [14].

For $\beta$, we choose

[TABLE]

This incompressible velocity field transports the mass parallel to the ramp with decreasing speed for increasing distance to the ramp. For comparison, we also use the constant velocity field

[TABLE]

6.2.1 Test 4: 2D Smooth initial data

We use the velocity field $\beta^{V}$, choose smooth initial data

[TABLE]

and compute the solution at time T = 0.5 using $V_{h}^{1}(\mathcal{T}_{h})$.

Figure 8 shows the error at time $T$ in the $L^{1}$ and $L^{\infty}$ norm. We observe second-order convergence in the $L^{1}$ norm. Due to the irregularity of the cut cells, the $L^{\infty}$ error is not smooth, indicated by the zig-zag behavior of the plotted results. We therefore use a least squares fit to compute the convergence rates. The rates vary, depending on the angle, and lie between 1.52 and 1.63. This slightly reduced order of convergence will need to be examined in more detail in the future. In general it is challenging to achieve full second-order accuracy on cut cells as the sizes of neighboring cells differ by several orders of magnitude and therefore errors of neighboring cells do not cancel the same way as on a structured mesh.

For comparison we show in table 1 the convergence rates for additional angles $\gamma$ as well as for using the constant velocity field $\beta^{C}$. The results for the constant velocity $\beta^{C}$ are very similar to the ones for $\beta^{V}$, except for $5^{\circ}$ degree. We note that it is reasonable that the results for a $5^{\circ}$ degree ramp are better as most cut cells have full length in flow direction, i.e., such a ramp contains significantly fewer ‘problematic’ cut cells than a ramp with a higher angle and the sizes of neighboring cut cells do not differ as strongly.

Finally, in figure 9 (left column), we show the solution for $\beta^{V}$ and $\gamma=30^{\circ}$ on a coarse mesh of $N=30$. The 1D profile along the cut boundary shows the expected sine curve. For the contour plot, we observe that the contour lines are straight lines all the way to the boundary. We note that for the initial data $u_{0}$ the contour lines are perpendicular to the ramp. Due to $\beta^{V}$ having decreasing speed for increasing distance to the ramp, the lines have been rotated during the simulation.

6.2.2 Test 5: 2D Discontinuous initial data

We use the velocity field $\beta^{V}$ and compute until time $T=0.4$ for discontinuous initial data

[TABLE]

Figure 9 (columns 2-4) shows the discrete solution for $N=30$ for $\gamma=30^{\circ}$ for both $V_{h}^{0}(\mathcal{T}_{h})$ and $V_{h}^{1}(\mathcal{T}_{h})$, the latter without and with limiter. For both $V_{h}^{0}(\mathcal{T}_{h})$ and $V_{h}^{1}(\mathcal{T}_{h})$ with limiter, the discrete solution stays between 0 and 1 and we do not observe overshoot. This is also indicated by the 1D line plots along the cut boundary. Further, we verified numerically for $V_{h}^{0}(\mathcal{T}_{h})$ that all matrix indices in the matrix $B^{-1}C$, compare (33) for the 1D version, are non-negative.

Examining the contour lines more closely, we find that for $V_{h}^{0}(\mathcal{T}_{h})$ and $V_{h}^{1}(\mathcal{T}_{h})$ without limiter the contour lines are straight lines all the way to the boundary, indicating that our stabilization does not add extra diffusion on the cut cells. For the case of $V_{h}^{1}(\mathcal{T}_{h})$ with limiter, we observe slightly more diffusion along the boundary. We attribute this to the limiter: on a Cartesian cell, we only limit in $x$- and $y$-direction. On cut cells however we also limit roughly in advection direction when reconstructing to neighboring centroids. We expect this behavior to improve if a more accurate limiter, e.g., the LP limiter [32] is used, but this is not the focus of this work.

7 Conclusions and Outlook

We have presented a new stabilization for DG schemes for linear scalar conservation laws on cut cell meshes that solves the small cell problem and makes explicit time stepping stable again. Our stabilization is designed to only let a certain portion of the inflow of a small cut cell stay in that small cut cell and to transport the remaining portion directly into the cut cell’s outflow neighbors. As a by-product, we reconstruct the proper domain of dependence of the small cut cell’s outflow neighbors. In that sense our stabilization relies on similar ideas as the $h$-box method [25, 9], but without an explicit geometry reconstruction.

The approach for realizing these ideas in a DG setting was inspired by the ghost penalty method [12] but significant changes were necessary to adjust the terms that were developed for elliptic problems to the setting of hyperbolic equations. In this contribution, we have focused on using two standard explicit time stepping schemes (the explicit Euler scheme and the standard second-order TVD RK scheme) and the standard Barth-Jespersen limiter but there are no explicit obstacles for working with other choices of time stepping schemes and limiters.

Our stabilization ensures conservation (in a slightly extended meaning). In one dimension, we have shown that the stabilized scheme is monotone for piecewise constant polynomials and total variation diminishing in the means for piecewise linear polynomials in combination with explicit time stepping schemes. In numerical tests we observed optimal convergence rates in the $L^{1}$ and $L^{\infty}$ norm for 1D. In 2D, the numerical results in the $L^{1}$ norm showed again full second-order convergence but the convergence rate in the $L^{\infty}$ norm was slightly reduced. This will be examined in more detail in future work. We also plan to extend the stabilization to higher-order polynomial degrees and to non-linear conservation laws.

Appendix A Proof of TVDM property

Proof of lemma 4.5.

The structure follows that of the proof of lemma 4.4. We will show that

[TABLE]

We exploit the fact that a moment basis is used. Testing with $w_{h}=\mathbb{I}_{j}$ yields the update for $\bar{u}_{j}$. For the update of the gradient, it is sufficient to consider the information that the postprocessing by applying the limiter provides. Again it is sufficient to consider the case $\alpha<2\lambda$.

In the following $u_{j}(x)$ denotes the full (linear) solution of cell $j$ evaluated at $x$, which can be written as $u_{j}(x)=\bar{u}_{j}+(x-x_{j})\nabla u_{j}$. We first provide the update formulae in the neighborhood of the cut cell $k_{1}$ induced by the stabilization (26):

[TABLE]

We decompose the sum of the TV in the means at time $t^{n+1}$ in terms $T_{1}$ to $T_{5}$:

[TABLE]

In the following, we will estimate each of the terms $T_{1}$ to $T_{5}$ separately. We will use that the MC limiter guarantees

[TABLE]

which implies

[TABLE]

Using this property and $\lambda<\frac{1}{2}$ (to guarantee the non-negativity of the pre-factors), there holds for two equidistant cells away from cells $k_{1}$ and $k_{2}$

[TABLE]

This yields bounds for $T_{1}$ and $T_{5}$:

[TABLE]

For $T_{2}$ we can show using (49a) and (49b)

[TABLE]

For this we need to verify that the two pre-factors are non-negative and that it is therefore allowed to pull them out of the absolute value. The pre-factor of the second term is non-negative due to (51). For the first pre-factor we obtain using (39)

[TABLE]

This implies with (50), $\lambda<\frac{1}{4}$, and (39)

[TABLE]

For $T_{3}$, we get from (49b) and (49c)

[TABLE]

Again we check the pre-factors. The first pre-factor is obviously non-negative because of (51) and $\frac{\lambda}{1-\alpha}\leq\frac{1}{2}$. Further, $\frac{1}{2}+\frac{\alpha-2\lambda}{2(1-\alpha)}\geq 0$ for $\lambda\leq\frac{1}{2}$. Therefore, using (50), the second pre-factor is non-negative as well. Finally, for $T_{4}$ we estimate

[TABLE]

For the last term there holds $\frac{\lambda}{1-\alpha}-\frac{\alpha}{2(1-\alpha)}\geq 0$ due to $2\lambda\geq\alpha.$ Together with other previously shown estimates, this implies that all three pre-factors are non-negative. Finally, summing up the estimates for $T_{1},\ldots,T_{5}$ implies the claim. ∎