[EN]

[EN]

# DoD Stabilization of linear hyperbolic PDEs on general cut-cell meshes

Gunnar Birke111Corresponding author: e-mail<EMAIL_ADDRESS>11 [DE]Applied
Mathematics Münster, Münster University  Christian Engwer 11 Sandra May 22
[SE]Department of Information Technology, Uppsala University  Florian
Streitbürger 33 [DE]Department of Mathematics, Dortmund University

###### Abstract

Standard numerical methods for hyperbolic PDEs require for stability a CFL-
condition which implies that the time step size depends on the size of the
elements of the mesh. On cut-cell meshes, elements can become arbitrarily
small and thus the time step size cannot take the size of small cut-cells into
account but has to be chosen based on the background mesh elements.

A remedy for this is the so called DoD (domain of dependence) stabilization
for which several favorable theoretical and numerical properties have been
shown in one and two space dimensions [4, 9]. Up to now the method is
restricted to stabilization of cut-cells with exactly one inflow and one
outflow face, i.e. triangular cut-cells with a no-flow face (see [4]).

We extend the DoD stabilization to cut-cells with multiple in- and outflow
faces by properly considering the flow distribution inside the cut-cell. We
further prove $L^{2}$-stability for the semi-discrete formulation in space and
present numerical results to validate the proposed extension.

##  1 Introduction

To avoid the mesh generation process of complex geometries, cut-cell methods
are an attractive alternative. The general idea is to start with a simple,
e.g. structured, background mesh and to cut out the desired geometry. This
results in a mesh with unstructured polyhedral cells, called cut-cells. Cut-
cells can have an arbitrary shape and can become arbitrarily small, causing
the small cell problem. To use explicit time stepping schemes for solving
hyperbolic conservation laws, the time step size would need to be chosen based
on the smallest cut-cell in the grid to ensure stability, which is in general
not feasible.

Developing solution approaches to the small cell problem in the context of
discontinuous Galerkin (DG) schemes is a very recent research branch,
including for example the work in [10, 8, 7]. In this contribution we focus on
the domain of dependence (DoD) stabilization, which was introduced in [4] for
the linear transport equation in one and two space dimensions and was extended
to non-linear systems in one space dimension in [9]. It is based on a DG
scheme in space to allow for higher-order approximations and possesses several
desirable theoretical properties. Numerical results show the expected higher-
order behavior in smooth flow and robustness around shocks.

Up to now, the DoD stabilization in two dimensions has only been used to
stabilize small triangular cut-cells for linear advection parallel to a ramp
[4, 11]. In this setup, the small stabilized cut-cells have exactly one inflow
and one outflow face, which was exploited in the design of the stabilization.
When moving to non-linear or coupled linear problems, this does not hold true
anymore and one has to deal with multiple inflow and outflow faces.

In this work, we take the first step in that direction by considering the
linear advection equation on a cut-cell mesh with arbitrary flow directions,
resulting in triangular cut-cells having 2 inflow and 1 outflow neighbor or
reverse. As this causes significant additional complications, we will only
consider piecewise constant polynomials here. We will prove $L^{2}$-stability
for the semi-discrete scheme and present numerical results to validate the new
extension of the stabilization terms.

The outline of the paper is as follows: we will first describe the problem
setup and then introduce the new extended stabilization. Then we will show the
$L^{2}$-stability proof and conclude with numerical results.

##  2 Problem setup

Let $\Omega\subset\mathbb{R}^{2}$ be an open and connected domain. We consider
linear hyperbolic systems of the form

$\displaystyle u_{t}+Au_{x}+Bu_{y}$ $\displaystyle=0$ $\displaystyle\text{in
}\Omega\times(0,T),$ (1a) $\displaystyle\tau u$ $\displaystyle=g$
$\displaystyle\text{on }\partial\Omega\times(0,T),$ (1b) $\displaystyle u$
$\displaystyle=u_{0}$ $\displaystyle\text{on }\Omega\times\\{t=0\\},$ (1c)

where $u(t)\in\mathbb{R}^{m}$, and $A,\,B\in\mathbb{R}^{m\times m}$ constant,
and $\tau$ is an appropriate boundary operator such that we only impose inflow
boundary conditions on incoming waves (and not everywhere on
$\partial\Omega$). We require that for any unit vector
$n=(n_{1},n_{2})^{T}\in\mathbb{S}^{1}$ the matrices $C=(n_{1}A+n_{2}B)$ are
symmetric and simultaneously diagonalizable over the reals, i.e. there is an
orthogonal matrix $O\in\mathbb{R}^{m\times m}$ and diagonal matrices
$\Lambda_{n}\in\text{Diag}(\mathbb{R}^{m\times m})$ such that
$n_{1}A+n_{2}B=O\Lambda_{n}O^{T}\leavevmode\nobreak\ \forall
n\in\mathbb{S}^{1}$.

In our numerical tests we will choose $\Omega=[0,1]^{2}$ and discretize it by
a structured grid $\widehat{\mathcal{M}_{h}}$. We then introduce an artificial
cut, a straight line going through the square, starting at $(x_{0},0)$ and
having an angle $\gamma$ relative to the $x$-axis. This creates an internal
boundary with two subdomains which we will resolve by a cut-cell mesh
$\mathcal{M}_{h}$. A sketch is contained in Fig. 1. So far we have always
tested with flow parallel to that cut. Here, we consider flow in various
directions, keeping the cut fixed.

${\widehat{\mathcal{M}}}_{h}$$\cap$${\bar{\Omega}}$$(x_{0},0)$$\gamma$$=$${\mathcal{M}_{h}}$$E\in\widehat{\mathcal{M}}_{h}$$E_{1}$$E_{2}$
Figure 1: Construction of the mesh: Out of the structured grid
$\widehat{\mathcal{M}_{h}}$ on the domain $\Omega$ the mesh $\mathcal{M}_{h}$
is constructed by introducing cut-cells $E_{1},E_{2}\subset
E\in\widehat{M_{h}}$ along the cut such that
$\bar{E_{1}}\cup\bar{E_{2}}=\bar{E}$.

We define the sets of internal and external faces as

$\mathcal{F}^{\text{int}}_{h}=\\{F=\partial E_{1}\cap\partial
E_{2}|E_{1},E_{2}\in\mathcal{M}_{h},\;E_{1}\neq
E_{2},\;|F|>0\\},\quad\mathcal{F}^{\text{ext}}_{h}=\\{F=\partial
E\cap\partial\Omega|E\in\mathcal{M}_{h},\;|F|>0\\},$

and the set of faces of an element $E\in\mathcal{M}_{h}$ by
$\mathcal{F}_{h}^{E}=\\{F\in\mathcal{F}^{\text{int}}_{h}\cup\mathcal{F}^{\text{ext}}_{h}|F\subset\partial
E\\}$. We choose a fixed local numbering on each of these sets and denote the
neighbor element of $E$ corresponding to a face $F_{i}\in\mathcal{F}_{h}^{E}$
by $E_{i}$. For the discretization in space we choose the discrete function
space

$\mathcal{V}^{0}_{h}(\mathcal{M}_{h})=\\{v_{h}\in
L^{2}(\Omega)^{m}\>|\>(v_{h})_{i}{|_{E}}\in\mathcal{P}^{0}(E)\>\forall(E,i)\in\mathcal{M}_{h}\times\\{1,..,m\\}\\}.$

For $v_{h}\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h})$ we denote by $v_{h}^{E}$
the value of $v_{h}$ on an element $E\in\mathcal{M}_{h}$.

For interior faces $F\in\mathcal{F}^{\text{int}}_{h}$ we fix an orientation of
the outer unit normal vector $n_{F}$ and denote the inner and outer element of
$F$ by $E_{1}$ and $E_{2}$, respectively. We then define average and jump by

$\textstyle\left\\{\mskip-5.0mu\\{{u_{h}}\right\\}\mskip-5.0mu\\}:=\frac{1}{2}(u_{h}^{{}_{E_{1}}}+u_{h}^{{}_{E_{2}}}),\qquad\left\llbracket{u_{h}}\right\rrbracket:=u_{h}^{{}_{E_{1}}}-u_{h}^{{}_{E_{2}}}.$

For exterior faces $F\in\mathcal{F}^{\text{ext}}_{h}$ we simply choose the
unit outer normal and denote by $u_{h}^{E_{F}}$ the solution on the cell that
lies in the interior of the domain and contains face $F$. We define the flux
matrix on a face $F$ as

$C_{F}=(n_{F})_{1}A+(n_{F})_{2}B=O\Lambda_{F}O^{T},$ (2)

where $(n_{F})_{1,2}$ denote the first and second component of the unit normal
vector $n_{F}$ on face $F$. Based on this, we define matrices which encode the
flux directions as

$C_{F}^{+}=O\Lambda_{F}^{+}O^{T},\quad
C_{F}^{-}=O\Lambda_{F}^{-}O^{T}\quad\text{with}\quad(\Lambda_{F}^{+})_{i,i}=\max(0,(\Lambda_{F})_{i,i})\quad\text{and}\quad(\Lambda_{F}^{-})_{i,i}=\min(0,(\Lambda_{F})_{i,i}).$

Note that $C_{F}=C_{F}^{+}+C_{F}^{-}.$ We also introduce a generalization of
the absolute value for such flux matrices by $|C_{F}|=C_{F}^{+}-C_{F}^{-}.$

The (unstabilized) upwind semi-discretization in space is then given as: Find
$u_{h}(t)\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h})$ such that

$(\partial_{t}u_{h}(t),v_{h})_{L^{2}(\Omega)}+a^{\text{upw}}_{h}(u_{h}(t),v_{h})+l_{h}(v_{h})=0\quad\forall\>v_{h}\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h})$
(3)

with

$\displaystyle a^{\text{upw}}_{h}(u_{h},v_{h})=$
$\displaystyle\sum_{F\in\mathcal{F}_{h}^{\text{ext}}}\int_{F}\langle
C_{F}^{+}u_{h}^{E_{F}},v_{h}^{E_{F}}\rangle
ds+\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\langle
C_{F}\left\\{\mskip-5.0mu\\{{u_{h}}\right\\}\mskip-5.0mu\\},\left\llbracket{v_{h}}\right\rrbracket\rangle+\frac{1}{2}\langle|C_{F}|\left\llbracket{u_{h}}\right\rrbracket,\left\llbracket{v_{h}}\right\rrbracket\rangle
ds,$ $\displaystyle l_{h}(v_{h})=$
$\displaystyle-\sum_{F\in\mathcal{F}_{h}^{\text{ext}}}\int_{F}\langle
C_{F}^{-}g,v_{h}^{E_{F}}\rangle ds.$

Here, $\langle\cdot,\cdot\rangle$ denotes the scalar product in
$\mathbb{R}^{m}$. We obtain $a^{\text{upw}}_{h}$ and $l_{h}$ by integration by
parts, where the integral over internal edges leads to jump terms (second sum)
and the boundary integral is split into outgoing waves (first sum) and
incoming waves (right hand side). We then discretize in time using the
explicit Euler scheme.

##  3 Stabilization

$F_{1}$$E_{1}$$F_{2}$$E_{2}$$F_{3}$$E_{3}$$F_{4}$$E_{4}$$E_{\text{cut}}$$\beta$
$\langle\beta,n_{F_{2}}\rangle^{+}\int_{F_{2}}\big{(}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}\omega_{1}(u_{h}^{E_{1}}-u_{h}^{E_{\text{cut}}})}\\\
\phantom{\langle\beta,n_{F_{2}}\rangle^{+}\int_{F_{2}}\big{(}}+{\color[rgb]{1,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{1,0,0}\omega_{4}(u_{h}^{E_{4}}-u_{h}^{E_{\text{cut}}})}\big{)}\llbracket
v\rrbracket ds$ Figure 2: Domain of dependence extension illustrated on a
four-sided cut-cell: We introduce a direct mass transport from cells $E_{1}$
and $E_{4}$ into $E_{2}$. The colored regions indicate the coupling between
the faces, the corresponding parts of the stabilization term on $F_{2}$ are
highlighted accordingly. From the graphic one can see that flow orientation
and geometry play a central part in determining the right mass redistribution
and the $\omega_{i}$ will have to be chosen accordingly. The face $F_{3}$
needs to be stabilized as well. Note that there should not be any coupling
between between $F_{2}$ and $F_{3}$ as both are outflow faces.

To deal with small cut-cells, an additional term
$J_{h}^{0}=\sum_{E_{\text{cut}}\in\mathcal{I}}J_{h}^{0,E_{\text{cut}}}$ with
$\mathcal{I}$ being the set of small cut-cells that require stabilization is
added to the space semi-discretization (we will comment on our choice in the
numerical results below). This results in the following scheme: Find
$u_{h}(t)\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h})$ such that

$(\partial_{t}u_{h}(t),v_{h})_{L^{2}(\Omega)}+a^{\text{upw}}_{h}(u_{h},v_{h})+J^{0}_{h}(u_{h},v_{h})+l_{h}(v_{h})=0\quad\forall\>v_{h}\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h}).$
(4)

Let $E_{\text{cut}}$ be a small cut-cell that requires stabilization. The idea
behind $J_{h}^{0,E_{\text{cut}}}$ is the following: When the time step size is
not chosen to respect the small size of $E_{\text{cut}}$, the domain of
dependence of an outflow neighbor $E$ of $E_{\text{cut}}$ will extend beyond
$E_{\text{cut}}$. We therefore extend the numerical DoD of $E$ such that it
receives information directly from the inflow neighbors of $E^{\text{cut}}$.
The amount of mass passed directly between neighbors of $E_{\text{cut}}$ is
chosen such that the update on $E_{\text{cut}}$ becomes stable.

To explain the main concept we consider the scalar linear transport equation
$u_{t}+\nabla\cdot(\beta u)=0$ with constant velocity
$\beta=(\beta_{1},\beta_{2})^{T}\in\mathbb{R}^{2}$, i.e., we use
$A=\beta_{1}$, $B=\beta_{2}$ in (1). Consider Fig. 2 where the domain of
dependence of $E_{2}$ potentially reaches into $E_{1}$ and $E_{4}$. In that
case, for choosing the time step based on the size of background cells, mass
(physically) is moved from $E_{1}$ and from $E_{4}$ to $E_{2}$ in a single
time step and this coupling must be mimicked by the stabilization.

We introduce an extension operator
$\mathcal{L}_{E^{\prime}}^{\text{ext}}(u_{h})(x)=u^{E^{\prime}}_{h}(x)$ for
$E^{\prime}\in\mathcal{M}_{h}$ which acts on functions
$u_{h}\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h})$ and corresponds to evaluating
the (constant) polynomial of cell $E^{\prime}$ outside its original support.
Using this operator, the following stabilization term was introduced in [4]
for triangular cut-cells with single outflow face $F_{2}$ and single inflow
face $F_{1}$:

$J^{0,E_{\text{cut}}}_{h}(u_{h},v_{h})=\eta_{E_{\text{cut}}}\int_{F_{2}}\langle\beta,n_{F_{2}}\rangle(\mathcal{L}_{E_{1}}^{\text{ext}}(u_{h})-u_{h}^{E_{\text{cut}}})\left\llbracket{v_{h}}\right\rrbracket
ds.$ (5)

This term introduces a direct coupling between $E_{1}$ and $E_{2}$. The
parameter $\eta_{E_{\text{cut}}}\in[0,1]$ controls how much mass is
transported via this coupling. In the case of exactly one inflow face $F_{1}$
and one outflow face $F_{2}$ (and one face with no-penetration b.c.), as
considered in [4], this term suffices to ensure stability. Here, however, the
cut-cell is allowed to have two inflow faces as illustrated in Fig. 2.
Furthermore, these inflow faces will in general influence multiple neighbors
of the cut-cell where the degree of influence depends on the geometry and flow
direction. To handle this we propose the following extension of (5)

$J^{0,E_{\text{cut}}}_{h}(u_{h},v_{h})=\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}\langle\beta,n_{F_{j}}\rangle^{+}(\mathcal{L}_{E_{i}}^{\text{ext}}(u_{h})-u_{h}^{E_{\text{cut}}})\left\llbracket{v_{h}}\right\rrbracket
ds.$ (6)

Here the $\omega_{i}\in\mathbb{R}$ provide information about the flow
distribution for incoming flow of the face $F_{i}$. Note that the extended
solutions of all neighbor elements are evaluated on all faces, and
$\omega_{i}=0$ if $F_{i}$ is not an inflow face. We provide a specific formula
of how to choose these weights for triangular cut-cells below in section 3.2.

Going back to the system case we allow $\omega_{i}\in\mathbb{R}^{m\times m}$
and arrive at our final formulation

$J^{0,E_{\text{cut}}}_{h}(u_{h},v_{h})=\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\langle\omega_{i}C^{+}_{F}(\mathcal{L}_{E_{i}}^{\text{ext}}(u_{h})-u_{h}^{E_{\text{cut}}}),\left\llbracket{v_{h}}\right\rrbracket\rangle
ds.$ (7)

Note that in defining $J^{0,E_{\text{cut}}}$ we assume that all normal vectors
$n_{j}$ for $F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}$ correspond to outward
normal vectors with respect to $E_{\text{cut}}$. In order to ensure
consistency and stability the weights $\omega_{i}$ must fulfill

$\displaystyle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}=$
$\displaystyle\,\text{Id}_{m\times m},$ (8)
$\displaystyle\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\omega_{i}C^{+}_{j}ds=$
$\displaystyle-\int_{F_{i}}C^{-}_{i}ds\qquad\forall\>F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}.$
(9)

Additionally we require that $\omega_{i}C_{F_{j}}^{+}$ is always symmetric and
positive semi-definite. Equation (8) can be understood as an assurance that
the overall amount of mass moved over a face by our stabilization is correct.
For the scalar case, we require the $\omega_{i}$ to build a convex
combination. Equation (9) means that a portion of the inflow is exactly
redistributed over all outflow face candidates. This in particular prevents
overshoots on small cut-cells for appropriate choices of
$\eta_{E_{\text{cut}}}$.

###  3.1 $L^{2}$-stability

Equipped with the aforementioned properties of our stabilization we can show
$L^{2}$-stability for the semi-discrete scheme in space. For brevity we
consider homogeneous inflow boundary conditions and assume that exact and
discrete solution vanish on the boundary, so that we can ignore any domain
boundary terms and focus on the situation of cut-cells. For $P^{0}$ functions
this means that they are zero in all boundary cells.

######  3.1

Consider (1) with homogeneous boundary conditions. Assume that the discrete
solution $u_{h}(t)$ vanishes on the boundary $\partial\Omega$ for all
$t\in(0,T)$. Let $u_{h}(t)\in\mathcal{V}^{0}_{h}(\mathcal{M}_{h})$ be the
solution to the semi-discrete problem (4). Then it holds

$||u_{h}(t)||_{L^{2}(\Omega)}\leq||u_{h}(0)||_{L^{2}(\Omega)}\qquad\forall\;t\in(0,T).$

  * Proof.

We choose $v_{h}=u_{h}(t)$ in (4). Any boundary terms vanish as
$\operatorname{trace}(u_{h}(t))=0$ on $\partial\Omega$. This yields

$(\partial_{t}u_{h}(t),u_{h}(t))_{L^{2}(\Omega)}+a^{\text{upw}}_{h}(u_{h}(t),u_{h}(t))+J_{h}^{0}(u_{h}(t),u_{h}(t))=0.$

By the fundamental theorem of calculus

$\int_{0}^{t}(\partial_{\tau}u_{h}(\tau),u_{h}(\tau))_{L^{2}(\Omega)}d\tau=\int_{0}^{t}\frac{d}{d{\tau}}\frac{1}{2}||u_{h}(\tau)||_{L^{2}(\Omega)}^{2}d\tau=\frac{1}{2}||u_{h}(t)||^{2}_{L^{2}(\Omega)}-\frac{1}{2}||u_{h}(0)||^{2}_{L^{2}(\Omega)}.$

To ease notation we will write $u=u_{h}(t)$ in the following. The goal now is
to show that $a^{\text{upw}}_{h}(u,u)+J_{h}^{0}(u,u)\geq 0$. We first consider
$a^{\text{upw}}_{h}(u,u)$, which expands into

$\displaystyle a^{\text{upw}}_{h}(u,u)=$
$\displaystyle\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\langle
C_{F}\\{\\!\\{u\\}\\!\\},\llbracket
u\rrbracket\rangle+\langle\tfrac{1}{2}|C_{F}|\llbracket u\rrbracket,\llbracket
u\rrbracket\rangle ds=\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\langle
C_{F}^{+}u^{E_{1}}+C_{F}^{-}u^{E_{2}},u^{E_{1}}-u^{E_{2}}\rangle ds$
$\displaystyle=$
$\displaystyle\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\langle
C_{F}^{+}u^{E_{1}},u^{E_{1}}\rangle-\langle
C_{F}^{+}u^{E_{1}},u^{E_{2}}\rangle+\langle
C_{F}^{-}u^{E_{2}},u^{E_{1}}\rangle-\langle
C_{F}^{-}u^{E_{2}},u^{E_{2}}\rangle ds.$ Now we add zeros (in form of
$\pm\tfrac{1}{2}\langle C_{F}^{-}u^{E_{1}},u^{E_{1}}\rangle$ and
$\pm\tfrac{1}{2}\langle C_{F}^{+}u^{E_{2}},u^{E_{2}}\rangle$) to get
$\displaystyle=$
$\displaystyle\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\frac{1}{2}\langle(C_{F}^{+}-C_{F}^{-})u^{E_{1}},u^{E_{1}}\rangle+\frac{1}{2}\langle(C_{F}^{+}+C_{F}^{-})u^{E_{1}},u^{E_{1}}\rangle-\langle
C_{F}^{+}u^{E_{1}},u^{E_{2}}\rangle$
$\displaystyle\phantom{\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}}+\langle
C_{F}^{-}u^{E_{2}},u^{E_{1}}\rangle+\frac{1}{2}\langle(C_{F}^{+}-C_{F}^{-})u^{E_{2}},u^{E_{2}}\rangle-\frac{1}{2}\langle(C_{F}^{+}+C_{F}^{-})u^{E_{2}},u^{E_{2}}\rangle
ds$ $\displaystyle=$
$\displaystyle\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\frac{1}{2}\langle|C_{F}|(u^{E_{1}}-u^{E_{2}}),u^{E_{1}}-u^{E_{2}}\rangle+\frac{1}{2}\langle
C_{F}u^{E_{1}},u^{E_{1}}\rangle-\frac{1}{2}\langle
C_{F}u^{E_{2}},u^{E_{2}}\rangle$
$\displaystyle\phantom{\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}}-\frac{1}{2}\langle
C_{F}u^{E_{1}},u^{E_{2}}\rangle+\frac{1}{2}\langle
C_{F}u^{E_{2}},u^{E_{1}}\rangle ds.$

Due to the symmetry of $C_{F}$, the terms in the last line cancel each other.
For the last two terms in the second to last line we use the divergence
theorem. Since $u_{h}(t)$ is elementwise constant, on $E\in\mathcal{M}_{h}$ it
holds that

$0=\int_{E}\nabla\cdot(\langle Au^{E},u^{E}\rangle,\langle
Bu^{E},u^{E}\rangle)dx=\sum_{F\in\mathcal{F}^{\text{int}}_{h}\cup\mathcal{F}_{h}^{\text{ext}},F\cap\partial
E\neq\emptyset}\int_{F}\langle C_{F}u^{E},u^{E}\rangle ds,$

and therefore, these terms vanish as well. Finally, due to $|C_{F}|$ being
positive semi-definite, we obtain positivity of $a^{\text{upw}}_{h}(u,u)$:

$a^{\text{upw}}_{h}(u,u)=\sum_{F\in\mathcal{F}_{h}^{\text{int}}}\int_{F}\frac{1}{2}\langle|C_{F}|(u^{E_{1}}-u^{E_{2}}),u^{E_{1}}-u^{E_{2}}\rangle
ds\geq 0.$

We now investigate $J^{0}_{h}$. For a small cut-cell
$E_{\text{cut}}\in\mathcal{I}$ we have due to (8)

$\displaystyle J_{h}^{0,E_{\text{cut}}}(u,u)=$
$\displaystyle\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\langle(\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}})-C_{F_{j}}^{+}u^{E_{\text{cut}}},u^{E_{\text{cut}}}-u^{E_{j}}\rangle
ds$ $\displaystyle=$
$\displaystyle\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{\text{cut}}}\rangle-\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{j}}\rangle$
$\displaystyle\phantom{\eta_{E_{\text{cut}}}\int_{F_{j}}}-\langle
C_{F_{j}}^{+}u^{E_{\text{cut}}},u^{E_{\text{cut}}}\rangle+\langle
C_{F_{j}}^{+}u^{E_{\text{cut}}},u^{E_{j}}\rangle ds.$ Adding again zeros (in
form of
$\pm\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{i}}\rangle$
and $\pm\langle C_{F_{j}}^{+}u^{E_{j}},u^{E_{j}}\rangle$) and reordering gives
$\displaystyle=$
$\displaystyle-\frac{1}{2}\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\langle
C_{F_{j}}^{+}u^{E_{\text{cut}}},u^{E_{\text{cut}}}\rangle-2\langle
C_{F_{j}}^{+}u^{E_{\text{cut}}},u^{E_{j}}\rangle+\langle
C_{F_{j}}^{+}u^{E_{j}},u^{E_{j}}\rangle$
$\displaystyle\phantom{-\frac{1}{2}\int_{F_{j}}}+\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{i}}\rangle-2\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{\text{cut}}}\rangle+\langle
C_{F_{j}}^{+}u^{E_{\text{cut}}},u^{E_{\text{cut}}}\rangle$
$\displaystyle\phantom{-\frac{1}{2}\int_{F_{j}}}-\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{i}}\rangle+2\langle\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\omega_{i}C_{F_{j}}^{+}u^{E_{i}},u^{E_{j}}\rangle-\langle
C_{F_{j}}^{+}u^{E_{j}},u^{E_{j}}\rangle ds$ $\displaystyle=$
$\displaystyle-\frac{1}{2}\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\langle
C_{F_{j}}^{+}(u^{E_{\text{cut}}}-u^{E_{j}}),u^{E_{\text{cut}}}-u^{E_{j}}\rangle$
$\displaystyle\phantom{-\frac{1}{2}\int_{F_{j}}}+\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\langle\omega_{i}C_{F_{j}}^{+}(u^{E_{i}}-u^{E_{\text{cut}}}),u^{E_{i}}-u^{E_{\text{cut}}}\rangle\qquad\text{($C_{F_{j}}^{+},\omega_{i}C_{F_{j}}^{+}$
symm. and \eqref{eq:flux-normalization})}$
$\displaystyle\phantom{-\frac{1}{2}\int_{F_{j}}}-\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\langle\omega_{i}C_{F_{j}}^{+}(u^{E_{i}}-u^{E_{j}}),u^{E_{i}}-u^{E_{j}}\rangle
ds$ $\displaystyle\overset{\eqref{eq:flux-distribution}}{=}$
$\displaystyle-\frac{1}{2}\eta_{E_{\text{cut}}}\sum_{F_{j}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F_{j}}\langle|C_{F_{j}}|(u^{E_{\text{cut}}}-u^{E_{j}}),u^{E_{\text{cut}}}-u^{E_{j}}\rangle$
$\displaystyle\phantom{-\frac{1}{2}\int_{F_{j}}}-\sum_{F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}}\langle\omega_{i}C_{F_{j}}^{+}(u^{E_{i}}-u^{E_{j}}),u^{E_{i}}-u^{E_{j}}\rangle
ds.$

Since $\eta_{E_{\text{cut}}}\in[0,1]$ the first term inside the sum can be
compensated with terms from $a_{h}^{\text{upw}}$. The second term is always
non-negative since $\omega_{i}C_{F_{j}}^{+}$ is always positive semi-definite.
Note that the second term corresponds to dissipation introduced by an extended
jump. This concludes the proof. ∎

###  3.2 Choice of parameters

To perform actual computations we need to select concrete $\omega_{i}$ in (7)
that fulfill properties (8) and (9). For the situation of linear
simultaneously diagonalizable hyperbolic systems and triangular cut-cells we
suggest
$\omega_{i}=|F_{i}|C_{F_{i}}^{-}(\sum_{F_{k}\in\mathcal{F}_{h}^{E_{\text{cut}}}}|F_{k}|C_{F_{k}}^{-})^{-1}$
for each $F_{i}\in\mathcal{F}_{h}^{E_{\text{cut}}}$. For linear advection,
this would result in $\omega_{i}=0$ for an outflow edge $F_{i}$ and
$\omega_{i}$ corresponding to some sort of weighted proportion of the total
inflow for an inflow edge $F_{i}$. We also need to set $\eta_{E_{\text{cut}}}$
for $E_{\text{cut}}\in\mathcal{I}$. A stable but not necessarily optimal
choice is
$\eta_{E_{\text{cut}}}=||\sum_{F\in\mathcal{F}_{h}^{E_{\text{cut}}}}\int_{F}\Lambda_{F}^{-}ds||_{\infty}$.

##  4 Numerical results

For the numerical tests we select angles
$\gamma,\theta,\rho_{1},\rho_{2}\in[0,2\pi)$ where $\gamma$ is the angle of
the cut, see Fig. 1, and set

$\Lambda_{1}=\begin{pmatrix}\scriptsize\cos(\rho_{1})&0\\\
0&\cos(\rho_{2})\end{pmatrix},\;\Lambda_{2}=\begin{pmatrix}\scriptsize\sin(\rho_{1})&0\\\
0&\sin(\rho_{2})\end{pmatrix},\;O=\begin{pmatrix}\scriptsize\cos(\theta)&-\sin(\theta)\\\
\sin(\theta)&\cos(\theta)\end{pmatrix}.$

Then our system matrices will be given by $A=O\Lambda_{1}O^{T}$ and
$B=O\Lambda_{2}O^{T}$.

The cut starts at $(x_{0},0)=(0.2001,0)$. As initial conditions we choose
$u_{0}(x)=O\begin{pmatrix}\scriptsize\sin(2\pi(x_{1}\cos(\rho_{1})+x_{2}\sin(\rho_{1})))\\\
\cos(2\pi(x_{1}\cos(\rho_{2})+x_{2}\sin(\rho_{2})))\end{pmatrix}$. Boundary
conditions are given by the exact solution. Let $N$ denote the number of
background cells in either coordinate direction. We compute the time step size
via $\Delta
t=0.4\frac{h}{\max_{n\in\mathbb{S}^{1}}||n_{1}\Lambda_{1}+n_{2}\Lambda_{2}||_{\infty}}$,
where $h=\frac{1}{N}$. The factor of $0.4$ allows that bigger cut-cells do not
need to be stabilized. We then choose
$\mathcal{I}=\\{E\in\mathcal{M}_{h}|\frac{|E|}{h^{2}}<0.4\\}$ as the set of
stabilized cut-cells. The final time of our simulation is $T=0.5$.

Our implementation is based on the DUNE framework (see [2], [1]), in
particular the dune-udg (see [3], [5]) and dune-pdelab modules. The local
subtriangulations for the cut-cells are computed by the TPMC library (see
[6]).

Figure 3: Plots of the error, measured in the $L^{1}$\- (blue, red) and
$L^{\infty}$-norm (orange, transparent purple). Left plot has values
$\rho_{1}=\frac{7\pi}{4}$, $\rho_{2}=\pi$ and $\gamma=40$°, right plot has
$\rho_{1}=\frac{7\pi}{4}$, $\rho_{2}=\frac{3\pi}{2}$ and $\gamma=30$°. In all
cases we have chosen $\theta=\frac{4\pi}{3}$. The green line is for reference.

Fig. 3 shows convergence plots for two particular setups. Note that the flow
directions have been chosen to not be parallel to the ramp angle $\gamma$. We
observe the expected order of convergence for a first order scheme in both the
$L^{1}$\- and the $L^{\infty}$-norm. In addition, all numerically computed
solution values, including those on small cut-cells, stayed within the bounds
of the initial state during the simulation, confirming the added stability of
the DoD stabilization.

##  5 Discussion and Outlook

We have extended the DoD stabilization to cut-cells with multiple
inflow/outflow faces for the case of component-wise and piecewise constant
trial and test functions and linear, simultaneously diagonalizable systems in
two dimensions. We have proven $L^{2}$-stability for the semi-discrete
setting. Numerically we observe full first-order convergence in different
numerical tests and no over/undershoot on cut-cells. In future work, we plan
to extend our method to more general systems, e.g., the acoustics and Euler
equations. An extension of the presented formulation to higher-order
approximations is ongoing research.

This work has been partially funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) as project 439956613 under contract numbers
EN 1042/5-1 and MA 7773/4-1/2.

## References

  * [1] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, R. Kornhuber, M. Ohlberger, and O. Sander, Computing, 82, 121-138 (2008)
  * [2] P. Bastian, M. Blatt, A. Dedner, C. Engwer, R. Klöfkorn, M. Ohlberger, and O. Sander, Computing, 82, 103-119 (2008)
  * [3] P. Bastian and C. Engwer, Int. Jour. for Num. Meth. in Eng., 79, 12, 1557-1576 (2009)
  * [4] C. Engwer, S. May, A. Nüßing, and F. Streitbürger, SIAM J. Sci. Comput. 42, 6, A3677-A3703 (2020).
  * [5] C. Engwer and F. Heimann, Proceedings of the DUNE user meeting, Stuttgart, Germany, Advances in DUNE (Springer Berlin, Heidelberg, 2012) pp. 89-100
  * [6] C. Engwer and A. Nüßing, ACM Trans. on Math. Soft., 44, 2, Art. No. 14 (2018)
  * [7] P. Fu and G. Kreiss, SIAM J. Sci. Comput., 43,4, A2404–A2424, 2021.
  * [8] A. Giuliani, SIAM J. Sci. Comput. 44, 1, A389-A415 (2022)
  * [9] S. May, F. Streitbürger, Appl. Math. Comput. 419, Art. 126854 (2022).
  * [10] S. Schoeder, S. Sticko, G. Kreiss and M. Kronbichler, Int. J. Numer. Meth. Engrg. 121, 13, 2979-3003 (2020)
  * [11] F. Streitbürger, G. Birke, C. Engwer, and S. May, Spectral and High Order Methods for Partial Differential Equations ICOSAHOM 2020+1., 137 (Springer International Publishing, 2023).