# Nonlinear Boundary Conditions for Initial Boundary Value Problems with
Applications in Computational Fluid Dynamics

Jan Nordström<EMAIL_ADDRESS>Department of Mathematics, Applied
Mathematics, Linköping University, SE-581 83 Linköping, Sweden  Department of
Mathematics and Applied Mathematics, University of Johannesburg, P.O. Box 524,
Auckland Park 2006, Johannesburg, South Africa

###### Abstract

We derive new boundary conditions and implementation procedures for nonlinear
initial boundary value problems (IBVPs) with non-zero boundary data that lead
to bounded solutions. The new boundary procedure is applied to nonlinear IBVPs
on skew-symmetric form, including dissipative terms. The complete procedure
has two main ingredients. In the first part (published in [1, 2]), the energy
and entropy rate in terms of a surface integral with boundary terms was
produced for problems with first derivatives. In this second part we
complement it by adding second derivative dissipative terms and bound the
boundary terms.

We develop a new nonlinear boundary procedure which generalise the
characteristic boundary procedure for linear problems. Both strong and weak
imposition of the nonlinear boundary conditions with non-zero boundary data
are considered, and we prove that the solution is bounded. The boundary
procedure is applied to four important IBVPs in computational fluid dynamics:
the incompressible Euler and Navier-Stokes, the shallow water and the
compressible Euler equations. Finally we show that stable discrete
approximations follow by using summation-by-parts operators combined with weak
boundary conditions.

###### keywords:

Nonlinear boundary conditions , Navier-Stokes equations , Euler equations ,
shallow water equations , energy and entropy stability , summation-by-parts

##  1 Introduction

This paper constitutes the second part of the stability theory for nonlinear
IBVPs on skew-symmetric form initiated in [1, 2]. In the first part, energy
rates in terms of surface integrals with boundary terms were derived using
skew-symmetric formulations. Here in the second part we add on second
derivative terms and develop a nonlinear boundary procedure which controls the
boundary terms. In fact, we generalise the well known linear characteristic
boundary procedure [3, 4] and prove that we can bound both strong and weak
boundary conditions in terms of given external data. We stress that this
second part (the new boundary procedure) in itself does not depend on the
first part. It can be used by other preferred formulations of nonlinear IBVPs,
as long as similar settings are obtained. Most of the analysis in this second
part is limited to IBVPs involving first derivatives (essentially hyperbolic
problems). However, the new nonlinear boundary procedure can also be used to
bound dissipative IBVPs involving second derivatives similar to what has been
done in the linear case [5, 6, 7]. We show how to formally extend the analysis
for the first derivative case to the one including second derivatives
(essentially parabolic and incompletely parabolic problems).

We exemplify the procedure on four important IBVPs in computational fluid
dynamics (CFD): the incompressible Euler equations (IEEs) and skew-symmetric
formulations of the shallow water equations (SWEs) and the compressible Euler
equations (CEEs) and the incompressible Navier-Stokes equations (INSEs). It
was shown in [1] that the original form of the velocity-divergence form of the
IEEs equations had the required skew-symmetric form and that it could be
derived for the SWEs. In [2] we showed that also the CEEs could be transformed
to skew-symmetric form with a quadratic mathematical (or generalised) entropy
by using specific variables involving square roots of density and pressure,
also used in [8, 9, 10, 11]. The continuous procedure in [1, 2]), was reused
by discretising the equations in space using summation-by-parts (SBP)
operators [12, 13] which discretely mimic the integration-by-parts (IBP)
procedure used in the continuous analysis. A provably stable nonlinear
approximation was obtained.

The stable boundary procedures that was assumed to exist in [1, 2] is the
target of this paper where we also include IBVPs with second derivatives. With
few exceptions, only nonlinear boundary conditions for solid walls (or
glancing boundaries) have been considered previously as for example in [14,
15, 16, 17, 18, 19]. Solid wall boundary conditions are notoriously simple and
straightforward to apply and implement due to their homogeneous nature, i.e.
no external non-zero data must be considered. In contrast to previous
investigations, we will for the first time (to the best of our knowledge)
treat the general case with non-homogeneous nonlinear boundary conditions and
derive estimates of the solution in terms of given data.

The remaining part of paper is organised as follows: In Section 2 we reiterate
and complement the main content in the previously published first part of
procedure [1, 2]. The new second part: how to formulate and impose general
nonlinear boundary conditions with non-zero data such that a bound is
obtained, is presented in Section 3. In Section 4, we show that nonlinear
inhomogeneous boundary conditions for important IBVPs in CFD can be
implemented in a bounded way using the new theoretical framework. Explicit
examples of boundary conditions and implementation procedures are given for
the IEEs, the SWEs, the CEEs and the INSEs. In Section 5 we return show that
the energy bounded continuous formulation lead to nonlinear stability of an
SBP based scheme, including non-zero boundary data. A summary is provided in
Section 6.

##  2 Part 1: Producing the energy and entropy rates

In this section we repeat the findings in [1, 2] to make the presentation self
contained and to set the stage for the new upcoming boundary condition
analysis. We also complement the governing equations with viscous fluxes to
include e.g. the Navier-Stokes equations. Following [1, 2], we consider the
following IBVP

$PU_{t}+(A_{i}(V)U)_{x_{i}}+B_{i}(V)U_{x_{i}}+C(V)U=\epsilon(D_{i}(U))_{x_{i}},\quad
t\geq 0,\quad\vec{x}=(x_{1},x_{2},..,x_{k})\in\Omega$ (2.1)

augmented with the initial condition $U(\vec{x},0)=F(\vec{x})$ in $\Omega$ and
the non-homogeneous boundary condition

$L(V)U=G(\vec{x},t),\quad t\geq
0,\quad\vec{x}=(x_{1},x_{2},..,x_{k})\in\partial\Omega.$ (2.2)

In (2.2), $L$ is the boundary operator and $G$ the boundary data. In (2.1),
Einsteins summation convention is used and $P$ is a symmetric positive
definite (or semi-definite) time-independent matrix that defines an energy
norm (or semi-norm) $\|U\|^{2}_{P}=\int_{\Omega}U^{T}PUd\Omega$. Note in
particular that $P$ is not a function of the state $U$. The $n\times n$
matrices $A_{i},B_{i},C$ are smooth functions of the $n$ component vector $V$.
The viscous fluxes containing first derivatives are included in $D_{i}$ and
$\epsilon$ (typically the inverse of the Reynolds number) is a parameter
measuring the influence of viscous forces. Furthermore, we assume that $U$ and
$V$ are smooth. The smoothness could be considered either as a result of the
dissipative terms in (2.1) for $\epsilon\neq 0$ or stemming from (2.1) for
short times with smooth initial data when $\epsilon\rightarrow 0$. Note that
(2.1) and (2.2) encapsulates both linear ($V\neq U$) and nonlinear ($V=U$)
problems.

###  2.1 Energy and entropy rates when $\epsilon\rightarrow 0$

We begin with a formal definition and proposition for $\epsilon\rightarrow 0$
that was already published in [1, 2].

######  Definition 2.1.

Firstly, the problem (2.1) with $\epsilon=0$ is energy conserving if
$\|U\|^{2}_{P}=\int_{\Omega}U^{T}PUd\Omega$ only changes due to boundary
effects. Secondly, it is energy bounded if $\|U\|^{2}_{P}\leq\|F\|^{2}_{P}$
for a minimal number of homogeneous $(G=0)$ boundary conditions (2.2).
Thirdly, it is strongly energy bounded if
$\|U\|^{2}_{P}\leq\|F\|^{2}_{P}+\int_{0}^{t}(\oint G^{T}G\ ds)dt$ for a
minimal number of non-homogeneous $(G\neq 0)$ boundary conditions (2.2).

######  Proposition 2.1.

The IBVP (2.1) with $\epsilon=0$ for linear ($V\neq U$) and nonlinear ($V=U$)
is energy conserving if

$B_{i}-A_{i}^{T}=0,\quad C+C^{T}=0,\quad i=1,2,..,k$ (2.3)

holds. It is energy bounded if it is energy conserving and the boundary
conditions (2.2) for $G=0$ lead to

$\oint\limits_{\partial\Omega}U^{T}(n_{i}A_{i})\\\ \ U\\\ \
ds=\oint\limits_{\partial\Omega}\frac{1}{2}U^{T}((n_{i}A_{i})+(n_{i}A_{i})^{T})U\\\
\ ds\geq 0.$ (2.4)

It is strongly energy bounded if it is energy conserving and the boundary
conditions (2.2) for $G\neq 0$ lead to

$\oint\limits_{\partial\Omega}U^{T}(n_{i}A_{i})\\\ \ U\\\ \
ds=\oint\limits_{\partial\Omega}\frac{1}{2}U^{T}((n_{i}A_{i})+(n_{i}A_{i})^{T})U\\\
\ ds\geq-\oint\limits_{\partial\Omega}G^{T}G\\\ \ ds,$ (2.5)

where $G=G(\vec{x},t)$ is independent of the solution $U$.

###### Proof.

The energy method applied to (2.1) with $\epsilon=0$ yields

$\frac{1}{2}\frac{d}{dt}\|U\|^{2}_{P}+\oint\limits_{\partial\Omega}U^{T}(n_{i}A_{i})\\\
\ U\\\ \ ds=\int\limits_{\Omega}(U_{x_{i}}^{T}A_{i}U-U^{T}B_{i}U_{x_{i}})\\\ \
d\Omega-\int\limits_{\Omega}U^{T}CU\\\ \ d\Omega,$ (2.6)

where $(n_{1},..,n_{k})^{T}$ is the outward pointing unit normal. The terms on
the right-hand side of (2.6) are cancelled by (2.3) leading to energy
conservation. If in addition (2.4) or (2.5) holds, an energy bound or a strong
energy bound respectively follows after integration in time. ∎

######  Remark 2.2.

For linear problems, a minimal number of boundary conditions that lead to a
bound is a necessary and sufficient condition for well-posedness [20, 21, 7,
22]. For nonlinear problems, a minimal number of boundary conditions that lead
to a bound is a necessary but not sufficient condition [23].

We aim for a provably stable nonlinear boundary procedure, and focus on smooth
solutions. However, we notice that (2.1) allows for a mathematical (or
generalised) entropy conservation law. The matrix $P$ is not a function of
$U$, such that the mathematical entropy defined below, is a quadratic not a
cubic function of $U$. Note also that the mathematical entropy is not related
to the thermodynamic one [24, 25].

######  Proposition 2.3.

The IBVP (2.1) together with conditions (2.3) leads to the entropy
conservation law

$\Phi_{t}+(\Psi_{i})_{x_{i}}=0,$ (2.7)

where $\Phi=U^{T}PU/2$ is the mathematical (or generalised) entropy and
$\Psi_{i}=U^{T}A_{i}U$ are the entropy fluxes.

###### Proof.

Multiplication of (2.1) from the left with $U^{T}$ in smooth region of the
domain yields

$(U^{T}PU/2)_{t}+(U^{T}A_{i}U)_{x_{i}}=(U_{x_{i}}^{T}A_{i}U-U^{T}B_{i}U_{x_{i}})-U^{T}CU.$
(2.8)

The right-hand side of (2.10) is cancelled by (2.3) leading to the entropy
conservation relation (2.7). ∎

The entropy conservation law (2.7) holds for smooth solutions. In general,
discontinuous solutions can develop for non-linear hyperbolic systems,
regardless of their initial smoothness. In the presence of discontinuities,
the mathematical entropy conservation law (2.7) becomes the entropy inequality

$\Phi_{t}+(\Psi_{i})_{x_{i}}\leq 0.$ (2.9)

If the mathematical entropy is convex (as in this case where $\Phi_{UU}=P$),
then the physically relevant weak solution makes the entropy decrease. The
non-standard compatibility conditions in this case reads

$\partial\Phi/\partial
U=\Phi_{U}=U^{T}P,\quad\Phi_{U}[P^{-1}((A_{i}U)_{x_{i}}+A^{T}_{i}U_{x_{i}}+CU)]=(U^{T}A_{i}U)_{x_{i}}=(\Psi_{i})_{x_{i}}.$
(2.10)

The entropy is identical to the energy, similar to what was found in [26, 23].
We will use energy to denote both quantities, but sometimes remind the reader
by using both notations explicitly.

###  2.2 Energy and entropy rates when $\epsilon>0$

The additional terms that we consider stem from the righthand side of (2.1).
Greens formula leads to

$\epsilon\int\limits_{\Omega}U^{T}(D_{i})_{x_{i}}d\Omega=\epsilon\oint\limits_{\partial\Omega}U^{T}(n_{i}D_{i})ds-\epsilon\int\limits_{\Omega}U_{x_{i}}^{T}D_{i}d\Omega.$
(2.11)

To obtain a bound for $\epsilon>0$, (2.3) is complemented with a condition on
the viscous fluxes which leads to

$B_{i}-A_{i}^{T}=0,\quad C+C^{T}=0,\quad U_{x_{i}}^{T}D_{i}\geq 0,\quad
i=1,2,..,k.$ (2.12)

The viscous boundary terms $U^{T}(n_{i}D_{i})$ can be grouped together with
the inviscid ones on the left hand side in the energy rate. We summarise the
result in the following proposition.

######  Proposition 2.4.

The IBVP (2.1) with $\epsilon>0$ for linear ($V\neq U$) and nonlinear ($V=U$)
solutions is energy bounded if (2.12) holds and the boundary conditions (2.2)
for $G=0$ lead to

$\oint\limits_{\partial\Omega}U^{T}(n_{i}A_{i})U-\epsilon
U^{T}(n_{i}D_{i})ds\geq 0.$ (2.13)

It is strongly energy bounded if it is energy bounded and the boundary
conditions (2.2) for $G\neq 0$ lead to

$\oint\limits_{\partial\Omega}U^{T}(n_{i}A_{i})U-\epsilon
U^{T}(n_{i}D_{i})ds\geq-\oint\limits_{\partial\Omega}G^{T}G\\\ \ ds,$ (2.14)

where $G=G(\vec{x},t)$ is independent of the solution $U$.

###### Proof.

The energy method applied to (2.1) with $\epsilon>0$ together with (2.12 and
(2.14) proves the claim. ∎

##  3 Part 2: The nonlinear boundary conditions and implementation procedure

We aim for a nonlinear boundary procedure that limits the boundary term in
terms of given data. We start with a couple of convenient transformations and
the case $\epsilon\rightarrow 0$. Consider the boundary term

$\oint\limits_{\partial\Omega}U^{T}(n_{i}A_{i})\\\ \ U\\\ \
ds=\frac{1}{2}\oint\limits_{\partial\Omega}U^{T}((n_{i}A_{i})+(n_{i}A_{i})^{T})U\\\
\ ds=\oint\limits_{\partial\Omega}U^{T}\tilde{A}(V)\\\ \ U\\\ \ ds,$ (3.1)

where $\tilde{A}(V)$ is symmetric. Next we rotate the matrix $\tilde{A}$ to
diagonal form as $T^{T}\tilde{A}T=\Lambda=diag(\lambda_{i})$ which gives us
new rotated variables $W=T^{-1}U$ and

$\oint\limits_{\partial\Omega}U^{T}\tilde{A}(V)\\\ \ U\\\ \ ds\\\ \
ds=\oint\limits_{\partial\Omega}W^{T}\Lambda\\\ \ W\\\ \
ds=\oint\limits_{\partial\Omega}(W^{+})^{T}\Lambda^{+}\\\ \
W^{+}+(W^{-})^{T}\Lambda^{-}\\\ \ W^{-}\\\ \
ds=\oint\limits_{\partial\Omega}\lambda_{i}W_{i}^{2}\\\ \ ds,$ (3.2)

where we again use Einsteins summation convention. In (3.2) we use the
indicator matrices $I^{-},I^{+}$ where $I^{-}+I^{+}=I$ to define,
$\Lambda^{+}=I^{+}\Lambda$ and $\Lambda^{-}=I^{-}\Lambda$ which are the
positive and negative parts of $\Lambda$ respectively while $W^{+}=I^{+}W$ and
$W^{-}=I^{-}W$ are the corresponding variables. The new rotated variables
$W=W(U)$ are functions of the solution in both the linear and nonlinear case.
In the nonlinear case, the diagonal matrix $\Lambda(U)$ is solution dependent
and not a priori bounded while in the linear case, $\Lambda(V)$ is bounded by
external data. This difference leads to notable differences in the boundary
condition procedure.

######  Remark 3.1.

For linear problems, the number of boundary conditions is equal to the number
of eigenvalues of $\tilde{A}(V)$ with the wrong (in this case negative) sign
[22]. Sylvester’s Criterion [27], shows that the number of boundary conditions
is equal to the number of $\lambda_{i}(V)$ with the wrong sign if the rotation
matrix $T$ is non-singular. In the nonlinear case where
$\lambda_{i}=\lambda_{i}(U)$ it is more complicated since multiple forms of
the boundary term $W^{T}\Lambda W$ may exist, see Section 4.1 below and [1, 2,
23] for examples. The nonlinear procedure developed here has similarities with
the linear characteristic boundary procedure [3, 4]. Hence, with a slight
abuse of notation we will sometimes refer to $\Lambda(U)$ as "eigenvalues" and
to the rotated variables $W(U)$ as "characteristic" variables, although
strictly speaking they are not.

We will impose the boundary conditions both strongly and weakly. For the weak
imposition we introduce a lifting operator $L_{C}$ that enforce the boundary
conditions in our governing equation (2.1) for $\epsilon\rightarrow 0$ as

$PU_{t}+(A_{i}(V)U)_{x_{i}}+A^{T}_{i}(V)U_{x_{i}}+C(V)U+L_{C}(U)=0,\quad t\geq
0,\quad\vec{x}=(x_{1},x_{2},..,x_{k})\in\Omega.$ (3.3)

The lifting operator for two smooth vector functions $\phi,\psi$ satisfies
$\int\limits\phi^{T}L_{C}(\psi)d\Omega=\oint\limits\phi^{T}\psi ds$ which
enables development of the essential parts of the numerical boundary procedure
in the continuous setting [28, 29].

###  3.1 One significant difference between linear and nonlinear boundary
conditions

Before attacking the nonlinear problem, we digress momentarily to the linear
case to introduce one aspect of our subsequent nonlinear analysis. In the
simplest possible version of (3.5) one can specify $W^{-}=g$ corresponding to
negative $\lambda_{i}(V)$ indicated by $\lambda_{i}^{-}=-|\lambda_{i}(V)|$.
Since $|\lambda_{i}(V)|$ are bounded, we obtain

$\oint\limits_{\partial\Omega}U^{T}\tilde{A}(V)\\\ \ U\\\ \ ds\\\ \
=\oint\limits_{\partial\Omega}W^{T}\Lambda\\\ \ W\\\ \
ds=\oint\limits_{\partial\Omega}(W^{+})^{T}\Lambda^{+}\\\ \
W^{+}+g^{T}\Lambda^{-}\\\ \ g\\\ \
ds\geq-\oint\limits_{\partial\Omega}G^{T}G\\\ \ ds,$ (3.4)

where $G_{i}=\sqrt{|\lambda_{i}^{-}(V)|}g_{i}$. Hence we get a strong energy
bound in terms of external data. However, in the nonlinear case, no estimate
is obtained since $\lambda_{i}^{-}(U)$ is not a priori bounded, see [30, 31]
for IEE examples.

###  3.2 The general form of nonlinear boundary conditions for rotated
variables

The starting point for the derivation of stable general nonlinear (and linear)
boundary conditions (2.2) is the form (3.2) of the boundary term. First we
find the formulation (3.2) with a minimal number of entries in $\Lambda^{-}$
[1, 2, 23] (there might be more than one formulation of the cubic boundary
terms). Next, we specify the characteristic variables $W^{-}$ in terms of
$W^{+}$ and external data. Finally we need a way to combine different ingoing
characteristic variables $W^{-}$ and a scaling possibility. The formulation

$S^{-1}(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+})=G\quad\text{or
equivalently}\quad\sqrt{|\Lambda^{-}|}W^{-}=R\sqrt{\Lambda^{+}}W^{+}+SG$ (3.5)

where $S^{-1}$ is a non-singular matrix combining values of $W^{-}$, $S^{-1}R$
combines values of $W^{+}$ while $G$ is external data will lead to a bound.
The boundary condition (3.5) where we used the notation
$|\Lambda|=diag(|\lambda_{i}|)$ and
$\sqrt{|\Lambda|}=diag(\sqrt{|\lambda_{i}|})$ is general in the sense that it
may involve all components of $W$ combined in an arbitrary way by the matrices
$S$ and $R$. Also,we observe that the data $G$ must contain some nonlinear
interaction since the boundary terms coming from the equations after applying
the energy method are cubic, not quadratic. A weak implementation of (3.5)
using a lifting operator is given by

$L_{C}(U)=L_{C}(2(I^{-}T^{-1})^{T}\Sigma(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+}-SG)),$
(3.6)

where $\Sigma$ is a penalty matrix.

The procedure to arrive at a general stable nonlinear inhomogeneous boundary
condition and implementation consist of the following steps for the unknowns
in $R,S,\Sigma$ in (3.5) and (3.6).

  1. 1.

Derive strong homogeneous ($G=0$) boundary conditions. This lead to conditions
on matrix $R$.

  2. 2.

Derive strong inhomogeneous ($G\neq 0$) boundary conditions. This lead to
conditions on matrix $S$.

  3. 3.

Derive weak homogeneous ($G=0$) boundary conditions. This lead to
specification of matrix $\Sigma$.

  4. 4.

Show that the weak inhomogeneous ($G\neq 0$) case of the boundary conditions
follow from 1-3 above.

The following Lemma (structured as the step-by-step procedure above) is the
main result of this paper.

######  Lemma 3.2.

Consider the boundary term described in (3.1),(3.2) and the boundary
conditions (3.5) implemented strongly or weakly using (3.6).

The boundary term augmented with 1\. strong nonlinear homogeneous boundary
conditions is positive semi-definite if the matrix $R$ is such that

$I-R^{T}R\geq 0.$ (3.7)

The boundary term augmented with 2\. strong nonlinear inhomogeneous boundary
conditions is bounded by the data $G$ if the matrix $R$ satisfies (3.7) with
strict inequality and the matrix $S$ is such that

$I-S^{T}S-(R^{T}S)^{T}\left(\sum^{\infty}_{k=0}(R^{T}R)^{k}\right)(R^{T}S)\geq
0.$ (3.8)

The boundary term augmented with 3\. weak nonlinear homogeneous boundary
conditions is positive semi-definite if the matrix $R$ satifies (3.7) and the
matrix $\Sigma$ is given by

$\Sigma=\sqrt{|\Lambda^{-}|}.$ (3.9)

The boundary term augmented with 4\. weak nonlinear inhomogeneous boundary
conditions is bounded by the data $G$ if the matrix $R$ satisfies (3.7) with
strict inequality, the matrix $S$ satisfies (3.8) and the matrix $\Sigma$ is
given by (3.9).

###### Proof.

We proceed in the step-by-step manner described above.

1\. The homogeneous boundary condition (3.5) implemented strongly (with $G=0$)
lead to

$W^{T}\Lambda
W=(W^{+})^{T}\Lambda^{+}W^{+}-(W^{-})^{T}|\Lambda^{-}|W^{-}=(\sqrt{\Lambda^{+}}W^{+})^{T}(I-R^{T}R)(\sqrt{\Lambda^{+}}W^{+})\geq
0$ (3.10)

if condition (3.7) holds.

2\. The inhomogeneous boundary condition (3.5) implemented strongly (with
$G\neq 0$) lead to

$W^{T}\Lambda
W=(\sqrt{\Lambda^{+}}W^{+})^{T}(\sqrt{\Lambda^{+}}W^{+})-(R(\sqrt{\Lambda^{+}}W^{+})+SG)^{T}(R(\sqrt{\Lambda^{+}}W^{+})+SG).$
(3.11)

Expanding (3.11), adding and subtracting $G^{T}G$ lead to the result

$W^{T}\Lambda W=\begin{bmatrix}\sqrt{\Lambda^{+}}W^{+}\\\
G\end{bmatrix}^{T}\begin{bmatrix}I-R^{T}R&-R^{T}S\\\
-S^{T}R&I-S^{T}S\end{bmatrix}\begin{bmatrix}\sqrt{\Lambda^{+}}W^{+}\\\
G\end{bmatrix}-G^{T}G.$ (3.12)

The matrix in (3.12) can be rotated into block-diagonal form with the upper
left block preserved and the lower right block being
$I-S^{T}S-(R^{T}S)^{T}\left(I-R^{T}R\right)^{-1}(R^{T}S)$. Since (3.7) holds
strictly, (3.8) follows.

3\. The homogeneous boundary condition (3.5) implemented weakly (with $G=0$)
using the lifting operator in (3.6) and applying the energy method lead to the
boundary term

$W^{T}\Lambda
W+2U^{T}(J^{-}T^{-1})^{T}\Sigma(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+}))=W^{T}\Lambda
W+2(W^{-})^{T}\Sigma(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+}).$

Collecting similar terms and making the choice (3.9) for $\Sigma$ transforms
the right hand side to

$(\sqrt{\Lambda^{+}}W^{+})^{T}(\sqrt{\Lambda^{+}}W^{+})+(\sqrt{|\Lambda^{-}|}W^{-})^{T}(\sqrt{|\Lambda^{-}|}W^{-})-2(\sqrt{|\Lambda^{-}|}W^{-})^{T}R(\sqrt{\Lambda^{+}}W^{+}).$
(3.13)

Adding and subtracting
$(R\sqrt{\Lambda^{+}}W^{+})^{T}(R\sqrt{\Lambda^{+}}W^{+})$ and using (3.7)
guarantee positivity of (3.13) since

$(\sqrt{\Lambda^{+}}W^{+})^{T}(I-R^{T}R)(\sqrt{\Lambda^{+}}W^{+})+(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+})^{T}(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+})\geq
0.$ (3.14)

4\. The inhomogeneous boundary condition (3.5) implemented weakly (with $G\neq
0$) using the lifting operator in (3.6) and the choice $\Sigma$ in (3.9) lead
to the boundary terms

$W^{T}\Lambda
W+2(\sqrt{|\Lambda^{-}|}W^{-})^{T}(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+}-SG)).$
(3.15)

By adding and subtracting $G^{T}G$ and rearranging, the boundary terms (3.15)
above can be written as

$\displaystyle(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+}-SG)^{T}(\sqrt{|\Lambda^{-}|}W^{-}-R\sqrt{\Lambda^{+}}W^{+}-SG)+$
$\displaystyle\begin{bmatrix}\sqrt{\Lambda^{+}}W^{+}\\\
G\end{bmatrix}^{T}\begin{bmatrix}I-R^{T}R&-R^{T}S\\\
-S^{T}R&I-S^{T}S\end{bmatrix}\begin{bmatrix}\sqrt{\Lambda^{+}}W^{+}\\\
G\end{bmatrix}-G^{T}G$ (3.16)

The first term is obviously positive semi-definite and the remaining part is
identical to the term in (3.12). ∎

######  Remark 3.3.

Lemma 3.2 can be used to prove that the estimates (2.4) and (2.5) in
Proposition 2.1 holds.

We are now ready to connect the characteristic boundary condition formulation
(3.5) with (2.2) in the original variables. The definitions $W=T^{-1}U$,
$W^{-}=I^{-}W$, $W^{+}=I^{+}W$, $\sqrt{|\Lambda^{-}|}=I^{-}\sqrt{|\Lambda|}$
and $\sqrt{\Lambda^{+}}=I^{+}\sqrt{|\Lambda|}$ transforms (3.5) to

$S^{-1}(I^{-}-RI^{+})\sqrt{|\Lambda|}T^{-1}U=G\quad\text{and hence}\quad
L=S^{-1}(I^{-}-RI^{+})\sqrt{|\Lambda|}T^{-1}.$ (3.17)

This concludes the analysis of the formulation of nonlinear boundary
conditions for first derivative terms.

###  3.3 Extension of the analysis to include second derivative terms

The boundary terms to consider in the case when $\epsilon>0$ are given in
Proposition 2.4. The argument in the surface integral can be conveniently
reformulated. By denoting the symmetric part of $n_{i}A_{i}$ as $\tilde{A}$
and introducing the notation $n_{i}D_{i}/2=\tilde{F}$ for the viscous flux we
rewrite the boundary terms as

$U^{T}(n_{i}A_{i})U-\epsilon U^{T}(n_{i}D_{i})=U^{T}\tilde{A}U-\epsilon
U^{T}\tilde{F}-\epsilon\tilde{F}^{T}U=\begin{bmatrix}U\\\
\epsilon\tilde{F}\end{bmatrix}^{T}\begin{bmatrix}\tilde{A}&-I\\\
-I&0\end{bmatrix}\begin{bmatrix}U\\\ \epsilon\tilde{F}\end{bmatrix}$ (3.18)

where $I$ is the identity matrix. We can now formally diagonalise the boundary
terms in (3.18), apply the boundary condition (3.5) and Lemma 3.2 to obtain
energy bounds.

However, it is instructive the take the analysis one step further. We start by
rotating the boundary terms (3.18) to block-diagonal form (assuming that
$\tilde{A}$ is non-singular)

$\begin{bmatrix}U\\\
\epsilon\tilde{F}\end{bmatrix}^{T}\begin{bmatrix}\tilde{A}&-I\\\
-I&0\end{bmatrix}\begin{bmatrix}U\\\
\epsilon\tilde{F}\end{bmatrix}=\begin{bmatrix}U-\epsilon\tilde{A}^{-1}\tilde{F}\\\
\epsilon\tilde{F}\end{bmatrix}^{T}\begin{bmatrix}\tilde{A}&0\\\
0&-\tilde{A}^{-1}\end{bmatrix}\begin{bmatrix}U-\epsilon\tilde{A}^{-1}\tilde{F}\\\
\epsilon\tilde{F}\end{bmatrix}.$ (3.19)

Next, we insert the previous transformations $\tilde{A}=(T^{-1})^{T}\Lambda
T^{-1}$ and $W=T^{-1}U$ into (3.19) to get

$\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\ \epsilon
T^{T}\tilde{F}\end{bmatrix}^{T}\begin{bmatrix}\Lambda&0\\\
0&-\Lambda^{-1}\end{bmatrix}\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\epsilon T^{T}\tilde{F}\end{bmatrix}.$ (3.20)

Finally we move the factor $\Lambda^{-1}$ from the matrix into the dependent
vectors as shown below

$\tilde{W}^{T}\tilde{\Lambda}\tilde{W}=\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\,\
\quad-\epsilon\Lambda^{-1}T^{T}\tilde{F}\end{bmatrix}^{T}\begin{bmatrix}\Lambda&0\\\
0&-\Lambda\end{bmatrix}\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\,\ \quad-\epsilon\Lambda^{-1}T^{T}\tilde{F}\end{bmatrix},$ (3.21)

which results in the same formal setting as for the first derivative problems.

######  Remark 3.4.

If there are no hidden linear dependencies in (3.21), then we can directly
apply the boundary condition (3.5) and Lemma 3.2 to obtain energy bounds.
However, for incompletely parabolic equations, the vector $\tilde{F}$ is
shorter than $W$ and reduce the number of boundary conditions below the number
indicated by the diagonal entries in $\tilde{\Lambda}$. We will discuss this
phenomenon in detail for the INSEs below.

##  4 Application of the general theory to initial boundary value problems in
CFD

We begin by considering boundary conditions for the first derivative cases
including the IEEs, the SWEs and the CEEs. Next, we extend the analysis to
include second derivatives and exemplify with the INSEs.

###  4.1 The 2D incompressible Euler equations

The incompressible 2D Euler equations in split form are

$PU_{t}+\frac{1}{2}\left[(AU)_{x}+AU_{x}+(BU)_{y}+BU_{y}\right]=0.$ (4.1)

where $U=(u,v,p)^{T}$ and

$P=\begin{bmatrix}1&0&0\\\ 0&1&0\\\ 0&0&0\end{bmatrix},\quad
A=\begin{bmatrix}u&0&1\\\ 0&u&0\\\ 1&0&0\end{bmatrix},\quad
B=\begin{bmatrix}v&0&0\\\ 0&v&1\\\ 0&1&0\end{bmatrix}.$ (4.2)

Since the matrices $A,B$ are symmetric, the formulation (4.1) is in the
required skew-symmetric form (3.3) We obtain an estimate in the semi-norm
$\|U\|^{2}_{P}=\int_{\Omega}U^{T}PUd\Omega$ involving only the velocities.
Note that the pressure $p$ includes a division by the constant density, and
hence has the dimension velocity squared.

By applying the transformations in Section 3.3, the boundary term gets the
form

$U^{T}(n_{1}A+n_{2}B)U=U^{T}_{n}\tilde{A}U_{n}=W^{T}\Lambda
W=(W^{+})^{T}\Lambda^{+}W^{+}+(W^{-})^{T}\Lambda^{-}W^{-}$ (4.3)

where the solution vector is rotated into the normal $(u_{n}=n_{1}u+n_{2}v)$
and tangential ($u_{\tau}=-n_{2}u+n_{1}v)$ direction, and $W=T^{-1}U_{n}$ such
that

$U_{n}=\begin{bmatrix}u_{n}\\\ u_{\tau}\\\
p\end{bmatrix},\tilde{A}\begin{bmatrix}u_{n}&0&1\\\ 0&u_{n}&0\\\
1&0&0\end{bmatrix},T^{-1}=\begin{bmatrix}1&0&1/u_{n}\\\ 0&1&0\\\
0&0&1\end{bmatrix},W=\begin{bmatrix}u_{n}+p/u_{n}\\\ u_{\tau}\\\
p\end{bmatrix},\Lambda=\begin{bmatrix}u_{n}&0&0\\\ 0&u_{n}&0\\\
0&0&-1/u_{n}\end{bmatrix}.$ (4.4)

At inflow where $u_{n}<0$

$W^{-}=\begin{bmatrix}u_{n}+p/u_{n}\\\
u_{\tau}\end{bmatrix},\quad\Lambda^{-}=\begin{bmatrix}u_{n}&0\\\
0&u_{n}\end{bmatrix},\quad W^{+}=p,\quad\Lambda^{+}=-1/u_{n}$ (4.5)

while at outflow with $u_{n}>0$ we get the reversed situation with

$W^{+}=\begin{bmatrix}u_{n}+p/u_{n}\\\
u_{\tau}\end{bmatrix},\quad\Lambda^{+}=\begin{bmatrix}u_{n}&0\\\
0&u_{n}\end{bmatrix},\quad W^{-}=p,\quad\Lambda^{-}=-1/u_{n}.$ (4.6)

By using the definitions in (4.5) and (4.6), it is straightforward to check
for boundedness using Lemma 3.2.

######  Example 4.1.

Recall that we must allow for some nonlinear scaling or interaction since the
boundary terms coming from the equations after applying the energy method are
typically cubic, not quadratic.

Dirichlet like inflow conditions on the velocities using the boundary
condition (3.5) are obtained by

$S^{-1}\left(\begin{bmatrix}\sqrt{|u_{n}|}&0\\\
0&\sqrt{|u_{n}|}\end{bmatrix}\begin{bmatrix}u_{n}+p/u_{n}\\\
u_{\tau}\end{bmatrix}-\begin{bmatrix}R_{1}\\\
R_{2}\end{bmatrix}p/\sqrt{|u_{n}|}\right)=G$ (4.7)

where $R=(1,0)^{T}$ removes the influence of the pressure $p$. This lead to
$I-R^{T}R=0$ and hence condition (3.7) is satisfied, but not strictly. This
makes the choice of $S$ in (4.7) (and (3.8)) irrelevant since only the
homogeneous case is possible. This leads to boundedness, but not strong
boundedness as defined in Proposition 2.1. A weak implementation require
$\Sigma=\sqrt{|\Lambda^{-}|}=diag(\sqrt{|u_{n}|},\sqrt{|u_{n}|})$ as specified
in (3.9).

For an outflow condition on the pressure, the boundary condition (3.5) becomes

$S^{-1}\left(p/\sqrt{|u_{n}|}-(R_{1},R_{2})\begin{bmatrix}\sqrt{|u_{n}|}&0\\\
0&\sqrt{|u_{n}|}\end{bmatrix}\begin{bmatrix}u_{n}+p/u_{n}\\\
u_{\tau}\end{bmatrix}\right)=G.$ (4.8)

The characteristic boundary condition for the pressure $p$ is recovered with
$R=(0,0)$ and hence condition (3.7) holds strictly. This leads to a strongly
energy bounded solution for $S$ such that $I-S^{T}S\geq 0$ as required in
(3.8) for a vanishing $R$. A weak implementation require
$\Sigma=\sqrt{|\Lambda^{-}|}=1/\sqrt{|u_{n}|}$ as specified in (3.9).

###  4.2 The 2D shallow water equations

The 2D SWEs on the skew-symmetric form presented in Proposition 2.1 were
derived in [1]. They are

$U_{t}+(AU)_{x}+A^{T}U_{x}+(BU)_{y}+B^{T}U_{y}+CU=0,$ (4.9)

where $U=(U_{1},U_{2},U_{3})^{T}=(\phi,\sqrt{\phi}u,\sqrt{\phi}v))^{T}$,
$\phi=gh$ is the geopontential [32], $h$ is the water height, $g$ is the
gravitational constant and $(u,v)$ is the fluid velocity in $(x,y)$ direction
respectively. The Coriolis forces are included in the matrix $C$ with the
function $f$ which is typically a function of latitude [33, 34]. Note that
$h>0$ and $\phi>0$ from physical considerations. The matrices in (4.9)
constitute a two-parameter family

$A=\begin{bmatrix}\alpha\frac{U_{2}}{\sqrt{U_{1}}}&(1-3\alpha)\sqrt{U_{1}}&0\\\
2\alpha\sqrt{U_{1}}&\frac{1}{2}\frac{U_{2}}{\sqrt{U_{1}}}&0\\\
0&0&\frac{1}{2}\frac{U_{2}}{\sqrt{U_{1}}}\end{bmatrix},B=\begin{bmatrix}\beta\frac{U_{3}}{\sqrt{U_{1}}}&0&(1-3\beta)\sqrt{U_{1}}\\\
0&\frac{1}{2}\frac{U_{3}}{\sqrt{U_{1}}}&0\\\
2\beta\sqrt{U_{1}}&0&\frac{1}{2}\frac{U_{3}}{\sqrt{U_{1}}}\end{bmatrix},C=\begin{bmatrix}0&0&0\\\
0&0&-f\\\ 0&+f&0\end{bmatrix}$ (4.10)

where the parameters $\alpha,\beta$ are arbitrary. (Symmetric matrices are
e.g. obtained with $\alpha=\beta=1/5$.)

By computing the boundary term, we find

$U^{T}(n_{1}A+n_{2}B)U=U^{T}\begin{bmatrix}\frac{\alpha+\beta}{2}u_{n}&\frac{1-\alpha}{2}n_{1}\sqrt{U_{1}}&\frac{1-\beta}{2}n_{2}\sqrt{U_{1}}\\\
\frac{1-\alpha}{2}n_{1}\sqrt{U_{1}}&\frac{1}{2}u_{n}&0\\\
\frac{1-\beta}{2}n_{2}\sqrt{U_{1}}&0&\frac{1}{2}u_{n}\end{bmatrix}U=U^{T}\begin{bmatrix}u_{n}&0&0\\\
&\frac{1}{2}u_{n}&0\\\ 0&0&\frac{1}{2}u_{n}\end{bmatrix}U$ (4.11)

and the (somewhat mysterious) dependency on the free parameters $\alpha$ and
$\beta$ vanishes for algebraic reasons. The relation (4.11) seemingly indicate
that we need three boundary conditions at inflow ($u_{n}<0$), and zero at
outflow ($u_{n}>0$). However, this is a nonlinear problem and as shown in
[23], it can be rewritten by changing variables and observing that
$u_{n}=(n_{1}U_{2}+n_{2}U_{3})/\sqrt{U_{1}}$. Reformulating (4.11) in new
variables give

$U^{T}(n_{1}A+n_{2}B)U=U^{T}\begin{bmatrix}u_{n}&0&0\\\ &\frac{1}{2}u_{n}&0\\\
0&0&\frac{1}{2}u_{n}\end{bmatrix}U=W^{T}\begin{bmatrix}-\frac{1}{2U_{n}\sqrt{U_{1}}}&0&0\\\
&\frac{1}{2U_{n}\sqrt{U_{1}}}&0\\\
0&0&\frac{1}{2U_{n}\sqrt{U_{1}}}\end{bmatrix}W,$ (4.12)

where
$W^{T}=(W_{1},W_{2},W_{3})=(U_{1}^{2},U_{1}^{2}+U_{n}^{2},U_{n}U_{\tau})$. The
variables $(U_{1},U_{n},U_{\tau})=(\phi,\sqrt{\phi}u_{n},\sqrt{\phi}u_{\tau})$
are directed in the normal ($U_{n}$) and tangential ($U_{\tau}$) direction
respectively. Since we search for a minimal number of boundary conditions, we
consider the formulation (4.11) for outflow, where no boundary conditions are
required. The relation (4.12) for inflow indicate that only two boundary
conditions are needed at inflow when $U_{n}<0$. To be specific, at inflow
where $U_{n}<0$ we find

$W^{-}=\begin{bmatrix}U_{1}^{2}+U_{n}^{2}\\\
U_{n}U_{\tau}\end{bmatrix},\quad\Lambda^{-}=\begin{bmatrix}\frac{1}{2U_{n}\sqrt{U_{1}}}&0\\\
0&\frac{1}{2U_{n}\sqrt{U_{1}}}\end{bmatrix},\quad
W^{+}=U_{1}^{2},\quad\Lambda^{+}=-\frac{1}{2U_{n}\sqrt{U_{1}}}.$ (4.13)

The definitions in (4.13) can be used to check any inflow conditions for
boundedness using Lemma 3.2.

######  Example 4.2.

Consider the general form of nonlinear boundary condition in (3.5).

With Dirichlet like inflow conditions on $U_{n},U_{\tau}$ we must remove the
influence of $U_{1}$. Using (3.5) we find

$S^{-1}\left(\begin{bmatrix}\frac{1}{\sqrt{2|U_{n}|\sqrt{U_{1}}}}&0\\\
0&\frac{1}{\sqrt{2|U_{n}|\sqrt{U_{1}}}}\end{bmatrix}\begin{bmatrix}U_{1}^{2}+U_{n}^{2}\\\
U_{n}U_{\tau}\end{bmatrix}-\begin{bmatrix}R_{1}\\\
R_{2}\end{bmatrix}\frac{1}{\sqrt{2|U_{n}|\sqrt{U_{1}}}}U_{1}^{2}\right)=G.$
(4.14)

In a similar way as for the Dirichlet inflow condition for the IEE we find
that this require $R=(1,0)^{T}$ which in turn lead to $I-R^{T}R=0$. Hence
condition (3.7) is satisfied, but not strictly, which makes the choice of $S$
in (4.14) (and (3.8)) irrelevant (similar to the inflow case in Example 4.4).
This leads to boundedness, but not strong boundedness as defined in
Proposition 2.1. A weak implementation require $\Sigma=\sqrt{|\Lambda^{-}|}$
in (4.13).

By instead specifying the characteristic variable $W^{-}$ directly (similar to
the outflow case in Example 4.4) we have $R=(0,0)^{T}$ and (3.7) holds
strictly. This lead to a strongly bounded solution for $S$ such that
$I-S^{T}S\geq 0$ as required in (3.8). A weak implementation require
$\Sigma=\sqrt{|\Lambda^{-}|}$ as given in (4.13), see (3.9).

###  4.3 The 2D compressible Euler equations

The 2D CEEs on the skew-symmetric form presented in Proposition 2.1 were
derived in [2]. They are

$P\Phi_{t}+(A\Phi)_{x}+A^{T}\Phi_{x}+(B\Phi)_{y}+B^{T}\Phi_{y}=0,$ (4.15)

where $\Phi=(\sqrt{\rho},\sqrt{\rho}u,\sqrt{\rho}v,\sqrt{p})^{T}$,
$P=diag(1,(\gamma-1)/2,(\gamma-1)/2,1)$ and

$A=\frac{1}{2}\begin{bmatrix}u&0&0&0\\\ 0&\frac{(\gamma-1)}{2}u&0&0\\\
0&0&\frac{(\gamma-1)}{2}u&0\\\
0&2(\gamma-1)\frac{\phi_{4}}{\phi_{1}}&0&(2-\gamma)u\end{bmatrix},\quad
B=\frac{1}{2}\begin{bmatrix}v&0&0&0\\\ 0&\frac{(\gamma-1)}{2}v&0&0\\\
0&0&\frac{(\gamma-1)}{2}v&0\\\
0&0&2(\gamma-1)\frac{\phi_{4}}{\phi_{1}}&(2-\gamma)v\end{bmatrix}.$ (4.16)

By rotating the Cartesian velocities to normal and tangential velocities at
the boundary, we obtain

$\Phi^{T}(n_{1}A+n_{2}B)\Phi=\Phi^{T}_{r}\begin{bmatrix}u_{n}&0&0&0\\\
0&\frac{(\gamma-1)}{2}u_{n}&0&(\gamma-1)\frac{\phi_{4}}{\phi_{1}}\\\
0&0&\frac{(\gamma-1)}{2}u_{n}&0\\\
0&(\gamma-1)\frac{\phi_{4}}{\phi_{1}}&0&(2-\gamma)u_{n}\end{bmatrix}\Phi_{r},$
(4.17)

where
$\Phi_{r}=(\phi_{1},\phi_{2},\phi_{3},\phi_{4})^{T}=(\sqrt{\rho},\sqrt{\rho}u_{n},\sqrt{\rho}u_{\tau},\sqrt{p})^{T}$.

The boundary term (4.17) can be rotated to diagonal form which yield the
boundary term $W^{T}\Lambda W$ where

$W=\begin{bmatrix}\phi_{1}\\\ \phi_{2}+2\phi_{4}^{2}/\phi_{2}\\\ \phi_{3}\\\
\phi_{4}\end{bmatrix}\quad\Lambda=\begin{bmatrix}u_{n}&0&0&0\\\
0&\frac{(\gamma-1)}{2}u_{n}&0&0\\\ 0&0&\frac{(\gamma-1)}{2}u_{n}&0\\\
0&0&0&(2-\gamma)u_{n}\Psi(M_{n})\end{bmatrix}.$ (4.18)

The last diagonal entry contain the factor $\Psi(M_{n})$ is a function of the
normal Mach number $M_{n}=u_{n}/c$. Explicitly we have

$\Psi(M_{n})=1-\frac{2(\gamma-1)}{\gamma(2-\gamma)}\frac{1}{M_{n}^{2}},$
(4.19)

which switches sign at $M_{n}^{2}=2(\gamma-1)/(\gamma(2-\gamma))$. As shown in
[2], this yields $|M_{n}|=1$ for $\gamma=\sqrt{2}$, while for $\gamma=1.4$ we
get $|M_{n}|\approx 0.95$.

Due to the sign shift in $\Psi$ at $M_{n}^{2}=\gamma(2-\gamma)/(2(\gamma-1))$
we get different cases for inflow when $u_{n}<0$. We find that for inflow when
$u_{n}<0,\Psi<0$, the relation (4.18) leads to

$W^{-}=\begin{bmatrix}\phi_{1}\\\ \phi_{2}+2\phi_{4}^{2}/\phi_{2}\\\
\phi_{3}\end{bmatrix},\quad\Lambda^{-}=\begin{bmatrix}u_{n}&0&0\\\
0&\frac{(\gamma-1)}{2}u_{n}&0\\\ 0&0&\frac{(\gamma-1)}{2}u_{n}\\\
\end{bmatrix},\quad
W^{+}=\phi_{4},\quad\Lambda^{+}=(2-\gamma)u_{n}\Psi(M_{n}).$ (4.20)

For inflow $u_{n}<0,\Psi>0$ we get $W^{-}=W$ and $\Lambda^{-}=\Lambda$ from
relation (4.18), i.e. all eigenvalues are negative.

In the outflow case, the shift in speed can be ignored since an alternate form
of (4.17) different from (4.18) exist. By contracting (4.17) we find that

$\Phi^{T}(n_{x}A+n_{y}B)\Phi=u_{n}(\phi_{1}^{2}+\frac{(\gamma-1)}{2}(\phi_{2}^{2}+\phi_{3}^{2})+\gamma\phi_{4}^{2})=\Phi^{T}_{r}\begin{bmatrix}u_{n}&0&0&0\\\
0&\frac{(\gamma-1)}{2}u_{n}&0&0\\\ 0&0&\frac{(\gamma-1)}{2}u_{n}&0\\\
0&0&0&\gamma u_{n}\end{bmatrix}\Phi_{r},$ (4.21)

which proves that no boundary conditions are necessary in the outflow case.

######  Example 4.3.

Consider the general form of boundary condition in (3.5).

With Dirichlet inflow conditions on $\phi_{1},\phi_{2},\phi_{3}$ for
$\Psi(M_{n})<0$ we find using (4.20)

$S^{-1}\left(\begin{bmatrix}\sqrt{|u_{n}|}&0&0\\\
0&\sqrt{\frac{(\gamma-1)}{2}|u_{n}|}&0\\\
0&0&\sqrt{\frac{(\gamma-1)}{2}|u_{n}|}\\\
\end{bmatrix}\begin{bmatrix}\phi_{1}\\\ \phi_{2}+2\phi_{4}^{2}/\phi_{2}\\\
\phi_{3}\end{bmatrix}-\sqrt{(2-\gamma)u_{n}\Psi(M_{n})}\begin{bmatrix}R_{3}\end{bmatrix}\phi_{4}\right)=G$
(4.22)

which lead to $R=(0,R_{2},0)^{T}$ where
$R_{2}=\sqrt{2(\gamma-1)/((2-\gamma)|\Psi|)}$. Hence condition (3.7) is
violated since $I-R^{T}R=1-R_{2}^{2}=1+(\Psi-1)/|\Psi|)=-1/|\Psi|)<0$ and no
bound can be found.

By instead specifying the characteristic variable $W^{-}$ directly (as for the
outflow case in Example 4.4 and inflow case in Example 4.2) we have
$R=(0,0,0)^{T}$ and (3.7) holds strictly. This lead to a strongly bounded
solution if $S$ is such that $I-S^{T}S\geq 0$ as required in (3.8). A weak
implementation require $\Sigma=\sqrt{|\Lambda^{-}|}$ as given in (4.20) and
(3.9).

###  4.4 The 2D incompressible Navier-Stokes equations

Here we extend the analysis of the IEEs (4.1) to the INSEs by adding on the
viscous fluxes to obtain

$PU_{t}+\frac{1}{2}\left[(AU)_{x}+AU_{x}+(BU)_{y}+BU_{y}\right]=\epsilon((PU_{x})_{x}+(PU_{y})_{y})$
(4.23)

where $U=(u,v,p)^{T}$, $\epsilon$ is the non-dimensional viscosity and the
matrices $P,A,B$ are given in (4.2). The righthand side of (4.23) identifies
the viscous fluxes. The energy method applied to (4.23) yields

$\frac{1}{2}\frac{d}{dt}\|U\|^{2}_{P}+\epsilon(\|U_{x}\|^{2}_{P}+\|U_{y}\|^{2}_{P})+\frac{1}{2}\oint\limits_{\partial\Omega}U^{T}(n_{1}A+n_{2}B)U-2\epsilon
U^{T}(n_{1}PU_{x}+n_{2}PU_{y})\\\ \ ds=0.$ (4.24)

By using the notations of (3.18) in Section (3.3), the argument of the surface
integral in (4.24) becomes

$U^{T}_{n}\tilde{A}U_{n}-\epsilon(U^{T}_{n}\tilde{F}+\tilde{F}^{T}U_{n})\quad\text{where}\quad\tilde{F}=(\tilde{F}_{n},\tilde{F}_{\tau},0)^{T}$
(4.25)

include the normal $(\tilde{F}_{n})$ and tangential $(\tilde{F}_{\tau})$ shear
stresses.

We can now formulate (4.25) in the matrix-vector form presented in Section
3.3, which leads to

$\tilde{W}^{T}\tilde{\Lambda}\tilde{W}=\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\,\
\quad-\epsilon\Lambda^{-1}T^{T}\tilde{F}\end{bmatrix}^{T}\begin{bmatrix}\Lambda&0\\\
0&-\Lambda\end{bmatrix}\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\,\ \quad-\epsilon\Lambda^{-1}T^{T}\tilde{F}\end{bmatrix}.$ (4.26)

To determine the minimal number of boundary conditions needed for a bound, we
consider the terms related to the lower right block
$(\Lambda^{-1}T^{T}\tilde{F})^{T}\Lambda(\Lambda^{-1}T^{T}\tilde{F})$ in
detail. That lower block is equal to

$(T^{T}\tilde{F})^{T}\Lambda^{-1}(T^{T}\tilde{F})=\begin{bmatrix}\tilde{F}_{n}\\\
\tilde{F}_{\tau}\\\
-\tilde{F}_{n}/u_{n}\end{bmatrix}^{T}\begin{bmatrix}-1/u_{n}&0&0\\\
0&-1/u_{n}&0\\\ 0&0&u_{n}\end{bmatrix}\begin{bmatrix}\tilde{F}_{n}\\\
\tilde{F}_{\tau}\\\
-\tilde{F}_{n}/u_{n}\end{bmatrix}=-\tilde{F}_{\tau}^{2}/u_{n}$ (4.27)

which reduces the $6\times 6$ diagonal matrix in relation (4.26) to a $4\times
4$ diagonal matrix. We find that

$\tilde{W}^{T}\tilde{\Lambda}\tilde{W}=\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\,\ \quad-\epsilon\tilde{F}_{\tau}\end{bmatrix}^{T}\begin{bmatrix}\Lambda&0\\\
0&-1/u_{n}\end{bmatrix}\begin{bmatrix}W-\epsilon\Lambda^{-1}T^{T}\tilde{F}\\\
\,\ \ -\epsilon\tilde{F}_{\tau}\end{bmatrix}.$ (4.28)

In full detail (with a slightly different scaling from the one used in Section
4.1), the result (4.28) becomes

$\tilde{W}^{T}\tilde{\Lambda}\tilde{W}=\begin{bmatrix}u_{n}^{2}+P-\epsilon\tilde{F}_{n}\\\
u_{n}u_{\tau}\ \ \ -\epsilon\tilde{F}_{\tau}\\\ \ \ \ \ \ \ P\
-\epsilon\tilde{F}_{n}\\\ \ \ \ \ \ \ \ \ \ \ -\
\epsilon\tilde{F}_{\tau}\end{bmatrix}^{T}\begin{bmatrix}1/u_{n}&0&0&0\\\
0&1/u_{n}&0&0\\\ 0&0&-1/u_{n}&0\\\
0&0&0&-1/u_{n}\end{bmatrix}\begin{bmatrix}u_{n}^{2}+P-\epsilon\tilde{F}_{n}\\\
u_{n}u_{\tau}\ \ \ -\epsilon\tilde{F}_{\tau}\\\ \ \ \ \ \ \ P\
-\epsilon\tilde{F}_{n}\\\ \ \ \ \ \ \ \ \ \ \ -\
\epsilon\tilde{F}_{\tau}\end{bmatrix}$ (4.29)

which is exactly what was obtained in [30, 31] using a slightly different
procedure. An inspection of (4.29) reveals that 2 boundary conditions are
required for a bound, both at inflow and outflow.

At inflow where $u_{n}<0$

$W^{-}=\begin{bmatrix}u_{n}^{2}+P-\epsilon\tilde{F}_{n}\\\ u_{n}u_{\tau}\ \ \
-\epsilon\tilde{F}_{\tau}\end{bmatrix},\quad\Lambda^{-}=\begin{bmatrix}1/u_{n}&0\\\
0&1/u_{n}\end{bmatrix},\quad W^{+}=\begin{bmatrix}P-\epsilon\tilde{F}_{n}\\\ \
\ \ \
-\epsilon\tilde{F}_{\tau}\end{bmatrix},\quad\Lambda^{+}=\begin{bmatrix}-1/u_{n}&0\\\
0&-1/u_{n}\end{bmatrix},$ (4.30)

while at outflow with $u_{n}>0$ we get the reversed situation with

$W^{+}=\begin{bmatrix}u_{n}^{2}+P-\epsilon\tilde{F}_{n}\\\ u_{n}u_{\tau}\ \ \
-\epsilon\tilde{F}_{\tau}\end{bmatrix},\quad\Lambda^{+}=\begin{bmatrix}1/u_{n}&0\\\
0&1/u_{n}\end{bmatrix},\quad W^{-}=\begin{bmatrix}P-\epsilon\tilde{F}_{n}\\\ \
\ \ \
-\epsilon\tilde{F}_{\tau}\end{bmatrix},\quad\Lambda^{-}=\begin{bmatrix}-1/u_{n}&0\\\
0&-1/u_{n}\end{bmatrix}.$ (4.31)

By using the definitions (4.30) and (4.31), it is straightforward to check for
boundedness using Lemma 3.2.

######  Example 4.4.

Consider the general form of nonlinear boundary condition in (3.5).

Dirichlet like inflow conditions on the velocities using the boundary
condition (3.5) and (4.30) become

$S^{-1}\left(\begin{bmatrix}\sqrt{|u_{n}|}&0\\\
0&/\sqrt{|u_{n}|}\end{bmatrix}^{-1}\begin{bmatrix}u_{n}^{2}+P-\epsilon\tilde{F}_{n}\\\
u_{n}u_{\tau}\ \ \
-\epsilon\tilde{F}_{\tau}\end{bmatrix}-\begin{bmatrix}R_{11}&R_{12}\\\
R_{21}&R_{22}\end{bmatrix}\begin{bmatrix}\sqrt{|u_{n}|}&0\\\
0&/\sqrt{|u_{n}|}\end{bmatrix}^{-1}\begin{bmatrix}P-\epsilon\tilde{F}_{n}\\\ \
\ \ \ -\epsilon\tilde{F}_{\tau}\end{bmatrix}\right)=G.$ (4.32)

The matrix $R$ with elements $R_{11}=R_{22}=1$ and $R_{12}=R_{21}=0$ removes
the influence of the pressure $p$ and shear stresses. This lead to
$I-R^{T}R=0$ and hence condition (3.7) is satisfied, but not strictly. This
makes the choice of $S$ in (4.7) (and (3.8)) irrelevant since only the
homogeneous case is possible. This leads to boundedness, but not strong
boundedness as defined in Proposition 2.1. A weak implementation require
$\Sigma=\sqrt{|\Lambda^{-}|}=diag(1/\sqrt{|u_{n}|},1/\sqrt{|u_{n}|})$ as
specified in (3.9).

The outflow conditions on the characteristic variables containing the pressure
and shear stresses, the boundary condition (3.5) and (4.31) leads to the
reversed formulation

$S^{-1}\left(\begin{bmatrix}\sqrt{|u_{n}|}&0\\\
0&/\sqrt{|u_{n}|}\end{bmatrix}^{-1}\begin{bmatrix}P-\epsilon\tilde{F}_{n}\\\ \
\ \ \ -\epsilon\tilde{F}_{\tau}\end{bmatrix}-\begin{bmatrix}R_{11}&R_{12}\\\
R_{21}&R_{22}\end{bmatrix}\begin{bmatrix}\sqrt{|u_{n}|}&0\\\
0&/\sqrt{|u_{n}|}\end{bmatrix}^{-1}\begin{bmatrix}u_{n}^{2}+P-\epsilon\tilde{F}_{n}\\\
u_{n}u_{\tau}\ \ \ -\epsilon\tilde{F}_{\tau}\end{bmatrix}\right)=G.$ (4.33)

Here we can choose the matrix $R$ to be the null matrix and and hence
condition (3.7) holds strictly. This leads to a strongly energy bounded
solution for $S$ such that $I-S^{T}S\geq 0$ as required in (3.8) for a
vanishing $R$. A weak implementation require $\Sigma=\sqrt{|\Lambda^{-}|}$ as
specified in (3.9).

###  4.5 Open questions for nonlinear boundary conditions

We will end this section by discussing some open questions stemming from the
nonlinear analysis above.

####  4.5.1 The number of boundary conditions in nonlinear IBVPs required for
boundedness

The boundary conditions for the SWEs and CEEs are similar in the sense that at
least two different formulations of the boundary terms can be found. The
minimal number of required conditions differ both in the inflow and outflow
cases. One common feature is that that no outflow conditions seem to be
necessary. Another similar feature is that the number of outflow conditions is
independent of the speed of sound for the CEEs and the celerity in the SWE
case. Both these effects differ from what one finds in a linear analysis.

######  Remark 4.5.

By substituting the IEE variables $W=(u_{n}+p/u_{n},u_{\tau},p/u_{n})^{T}$ in
(4.3) with $W=(u_{n},u_{\tau},\sqrt{p})^{T}$ (similar to the ones used in the
CEE and SWE cases) one obtains a similar situation also for the IEEs. The
eigenvalues for the IEEs transform from $\Lambda=diag(u_{n},u_{n},-u_{n})$ to
$\Lambda=diag(u_{n},u_{n},2u_{n})$ which leads to different number of boundary
conditions.

####  4.5.2 The effect of nonlinear boundary conditions on uniqueness and
existence

Roughly speaking, a minimal number of dissipative boundary conditions in the
linear case leads to uniqueness by the fact that it determines the normal
modes of the solution [20, 21]. The minimal number of boundary conditions can
also be obtained using the energy method, see [22]. If uniqueness and
boundedness for a minimal number of boundary conditions are given, existence
can be shown (e.g. using Laplace transforms or difference approximations [35,
36]). For linear IBVPs, the number of boundary conditions is independent of
the solution and only depend on known external data. For nonlinear IBVPs, that
is no longer the case, and the number may change in an unpredictable way as
the solution develops in time. In addition, as we have seen above, it also
varies depending on the particular formulation choosen. This is confusing and
raises a number of questions that we will speculate on below.

Let us consider the SWEs as an example. The two forms of the boundary terms
given in (4.12) were

$U^{T}\begin{bmatrix}u_{n}&0&0\\\ &\frac{1}{2}u_{n}&0\\\
0&0&\frac{1}{2}u_{n}\end{bmatrix}U=W^{T}\begin{bmatrix}-\frac{1}{2U_{n}\sqrt{U_{1}}}&0&0\\\
&\frac{1}{2U_{n}\sqrt{U_{1}}}&0\\\
0&0&\frac{1}{2U_{n}\sqrt{U_{1}}}\end{bmatrix}W,$ (4.34)

where
$W^{T}=(W_{1},W_{2},W_{3})=(U_{1}^{2},U_{1}^{2}+U_{n}^{2},U_{n}U_{\tau})$ and
$(U_{1},U_{n},U_{\tau})=(\phi,\sqrt{\phi}u_{n},\sqrt{\phi}u_{\tau})$. Based on
the two formulations in (4.34), one may base the boundary procedure on one of
the following four scenarios.

  1. 1.

The left formulation with variable $U$ at both inflow and outflow boundaries.

  2. 2.

The right formulation with variable $W$ at both inflow and outflow boundaries.

  3. 3.

The left formulation with variable $U$ at inflow and the right formulation
with $W$ at outflow boundaries.

  4. 4.

The right formulation with variable $W$ at inflow and the left formulation
with $U$ at outflow boundaries.

Scenario 1 would in a one-dimensional setting lead to three boundary
conditions all applied on the inflow boundary. Scenario 2 would also give
three boundary conditions, but now two would be applied on the inflow boundary
and one on the outflow boundary. Scenario 3 would lead to four boundary
conditions, three on the inflow and one on the outflow boundary. Scenario 4
would only give two boundary conditions, both applied on the inflow boundary.

If the above scenarios were interpreted in the linear sense, both Scenario 1
and 2 would determine the solution uniquely. (One of them would be a better
choice than the other depending on the growth or decay of the solution away
from the boundary [20, 21].) In scenario 3, the solution would be
overspecifed, leading to loss of existence. In scenario 4, the solution would
be underspecified, leading to loss of uniqueness. In summary: Scenario 1 and 2
may lead to acceptable solutions, Scenario 3 give no solution at all, while
scenario 4 yield a bounded solution with limited (or no) accuracy.

However, since these results are nonlinear, the above summary is merely
speculative. We do not know exactly how to interpret them, since the present
nonlinear theory is incomplete. We only know that boundedness is required. It
also seems likely though that scenario 1 and 2 are should be preferred over
scenario 3 and 4. The speculations in this section are of course equally valid
(or not valid) for the CEEs and IEEs.

######  Remark 4.6.

Whether a reformulation that leads to different minimal numbers of boundary
conditions can be done also in the second derivative cases, e.g. in the INSEs,
is not clear. We think not, since the added viscous terms are linear and do
not contribute to the diagonal entries.

##  5 Semi-discrete nonlinear stability

Consider the extended version (3.3) of (2.1) rewritten (using Einsteins
summation convention) for clarity

$PU_{t}+(A_{i}U)_{x_{i}}+A^{T}_{i}U_{x_{i}}+CU+L_{C}=0,\quad t\geq
0,\quad\vec{x}=(x_{1},x_{2},..,x_{k})\in\Omega.$ (5.1)

Equation (5.1) is augmented with the initial condition
$U(\vec{x},0)=F(\vec{x})$ in $\Omega$ and boundary conditions of the form
(3.5) on $\partial\Omega$. Furthermore $A_{i}=A_{i}(U)$, $C=C(U)$ and $P$ are
$n\times n$ matrices while $U$ and $L_{C}$ are $n$ vectors. $L_{C}$ is the
continuous lifting operator of the form (3.6) implementing the boundary
conditions weakly. A straightforward approximation of (5.1) on summation-by-
parts (SBP) form in $M$ nodes is

$(P\otimes I_{M})\vec{U}_{t}+{\bf D_{x_{i}}}{\bf A_{i}}\vec{U}+{\bf
A_{i}^{T}}{\bf D_{x_{i}}}\vec{U}+{\bf
C}\vec{U}+{\vec{L}_{D}}=0,\quad\vec{U}(0)=\vec{F}$ (5.2)

where $\vec{U}=(\vec{U}_{1}^{T},\vec{U}_{2}^{T},...,\vec{U}_{n}^{T})^{T}$
include approximations of $U=(U_{1},U_{2},...,U_{n})^{T}$ in each node. The
discrete lifting operator ${\vec{L}_{D}}(\vec{U})$ implements the boundary
conditions in a similar way to $L_{C}(U)$ and $\vec{F}$ denotes the discrete
initial data with the continuous initial data injected in the nodes. The
matrix elements of ${\bf A_{i}},{\bf C}$ are matrices with node values of the
matrix elements in $A_{i},C$ injected on the diagonals as exemplified below

$A_{i}=\begin{pmatrix}a_{11}&\ldots&a_{1n}\\\ \vdots&\ddots&\vdots\\\
a_{n1}&\ldots&a_{nn}\end{pmatrix},\quad{\bf A_{i}}=\begin{pmatrix}{\bf
a_{11}}&\ldots&{\bf a_{1n}}\\\ \vdots&\ddots&\vdots\\\ {\bf
a_{n1}}&\ldots&{\bf a_{nn}}\end{pmatrix},\quad{\bf
a_{ij}}=diag(a_{ij}(x_{1},y_{1}),\ldots,a_{ij}(x_{M},y_{M})).$ (5.3)

Moreover ${\bf D_{x_{i}}}=I_{n}\otimes D_{x_{i}}$ where $\otimes$ denotes the
Kronecker product, $I_{n}$ is the $n\times n$ identity matrix,
$D_{x_{i}}=P_{\Omega}^{-1}Q_{x_{i}}$ are SBP difference operators,
$P_{\Omega}$ is a positive definite diagonal volume quadrature matrix that
defines a scalar product and norm such that

$(\vec{U},\vec{V})_{\Omega}=\vec{U}^{T}P_{\Omega}\vec{V}\approx\int\limits_{\Omega}U^{T}Vd\Omega,\quad\text{and}\quad(\vec{U},\vec{U})_{\Omega}=\|\vec{U}\|^{2}_{\Omega}=\vec{U}^{T}P_{\Omega}\vec{U}\approx\int\limits_{\Omega}U^{T}Ud\Omega=\|U\|^{2}_{\Omega}.$
(5.4)

Following [37] we introduce the discrete normal
$\mathbf{N}=(N_{1},N_{2},...,N_{k})$ approximating the continuous normal
$\mathbf{n}=(n_{1},n_{2},...,n_{k})$ in the $N$ boundary nodes and a
restriction operator $E$ that extracts the boundary values $E\vec{U}$ from the
total values. We also need a positive definite diagonal boundary quadrature
$P_{\partial\Omega}=diag(ds_{1},ds_{2},...,ds_{N})$ such that
$\oint_{\partial\Omega}U^{T}Uds\approx(E\vec{U})^{T}P_{\partial\Omega}(\vec{U})=(EU)_{i}^{2}ds_{i}$.
With this notation in place (again using Einsteins summation convention), the
SBP constraints for a scalar variable becomes

$Q_{x_{i}}+Q_{x_{i}}^{T}=E^{T}P_{\partial\Omega}N_{i}E,$ (5.5)

which leads to the scalar summation-by-parts formula mimicking integration-by-
parts

$(\vec{U},D_{x_{i}}\vec{V})=\vec{U}^{T}P_{\Omega}(D_{x_{i}}\vec{V})=-(D_{x_{i}}\vec{U},\vec{V})+(E\vec{U})^{T}P_{\partial\Omega}N_{i}(E\vec{V}).$
(5.6)

The scalar SBP relations in (5.5),(5.6), correspond to the SBP formulas for a
vector with $n$ variables as

$(\vec{U},{\bf D_{x_{i}}}\vec{V})=\vec{U}^{T}(I_{n}\otimes P_{\Omega})({\bf
D_{x_{i}}}\vec{V})=-({\bf
D_{x_{i}}}\vec{U},\vec{V})+(E\vec{U})^{T}(I_{n}\otimes
P_{\partial\Omega})N_{i}(E\vec{V}).$ (5.7)

It remains to construct the discrete lifting operator $\vec{L}_{D}$ (often
called the SAT term [12, 13]) such that we can reuse the continuous analysis.
We consider an operator of the form $\vec{L}_{D}=(I_{n}\otimes
P_{\Omega})(DC)\vec{L}_{C}$. The transformation matrix $DC$ first extracts the
boundary nodes from the volume nodes, secondly permute the dependent variables
from being organised as $(E\vec{U}_{1},E\vec{U}_{2},..,E\vec{U}_{n})^{T}$ to
$((E\vec{U})_{1},(E\vec{U})_{2},..,(E\vec{U})_{N})^{T}$ using the permutation
matrix $P_{erm}$ and thirdly numerically integrate the resulting vector
against the continuous lifting operator $\vec{L}_{C}$ (now applied to the
discrete solution). More specifically we have

$\displaystyle\vec{L}_{D}=(I_{n}\otimes P_{\Omega})^{-1})(DC)\vec{L}_{C}$
$\displaystyle\ \ \ DC=(I_{n}\otimes
E^{T})(P_{erm})^{T}(P_{\partial\Omega}\otimes I_{n}),$ (5.8)
$\displaystyle\vec{L}_{C}=((L_{C})_{1},(L_{C})_{2},...,(L_{C})_{N})^{T}$
$\displaystyle(L_{C})_{j}=(2(J^{-}T^{-1})^{T}\Sigma(\sqrt{|\Lambda^{-}|}\vec{W}^{-}-R\sqrt{\Lambda^{+}}\vec{W}^{+}-S\vec{G}))_{j}.$
(5.9)

We can now prove the semi-discrete correspondence to Proposition 2.1.

######  Proposition 5.1.

Consider the nonlinear scheme (5.2) with $\vec{L}_{D}$ defined in (5.8) and
(5.9).

It is nonlinearly stable for $\vec{G}=0$ if the relations (3.7) and (3.9) in
Lemma 3.2 hold and the solution satisfies the estimate

$\|\vec{U}\|_{P\otimes P_{\Omega}}^{2}\leq\|\vec{F}\|_{P\otimes
P_{\Omega}}^{2}.$ (5.10)

It is strongly nonlinearly stable for $\vec{G}\neq 0$ if the relations
(3.7),(3.8) and (3.9) in Lemma 3.2 hold and the solution satisfies the
estimate

$\|\vec{U}\|_{P\otimes P_{\Omega}}^{2}\leq\|\vec{F}\|_{P\otimes
P_{\Omega}}^{2}+2\int_{0}^{t}\sum_{j=1,N}[\vec{G}^{T}\vec{G})]_{j}ds_{j}\ dt.$
(5.11)

In (5.10) and (5.11), $\vec{F}$ and $\vec{G}$ are external data from $F$ and
$G$ injected in the nodes.

###### Proof.

The discrete energy method (multiply (5.2) from the left with
$\vec{U}^{T}(I_{n}\otimes P_{\Omega}$) yields

$\vec{U}^{T}(P\otimes P_{\Omega})\vec{U}_{t}+(\vec{U},{\bf D_{x_{i}}}{\bf
A_{i}}\vec{U})+({\bf A_{i}}\vec{U},{\bf
D_{x_{i}}}\vec{U})+(\vec{U},{\vec{L}_{D})}=0,$ (5.12)

where we have used that $(I_{n}\otimes P_{\Omega})$ commutes with ${\bf
A_{i}}$ (since the matrices have diagonal blocks) and that the symmetric part
of $C$ is zero. The SBP constraints (5.7) and the notation
$\vec{U}^{T}(P\otimes P_{\Omega})\vec{U}=\|\vec{U}\|_{P\otimes
P_{\Omega}}^{2}$ simplifies (5.12) to

$\dfrac{1}{2}\dfrac{d}{dt}\|\vec{U}\|_{P\otimes
P_{\Omega}}^{2}+\vec{U}^{T}(I_{n}\otimes E^{T}P_{\partial\Omega}N_{i}E){\bf
A_{i}}\vec{U}+(\vec{U},{\vec{L}_{D}})=0.$ (5.13)

The semi-discrete energy rate in (5.13) mimics the continuous energy rate in
the sense that only boundary terms remain. To make use of the already
performed continuous energy analysis, we expand the boundary terms and exploit
the diagonal form of $P_{\partial\Omega}$ and the relation
$W_{i}=(T^{-1}(E\vec{U}))_{i}$. The result is

$\vec{U}^{T}(I_{n}\otimes E^{T}P_{\partial\Omega}N_{x_{i}}E){\bf
A_{i}}\vec{U}=\sum_{j=1,N}[(E\vec{U})^{T}(N_{i}{\bf
A_{i}})(E\vec{U})]_{j}ds_{j}=\sum_{j=1,N}[\vec{W}^{T}\Lambda\vec{W}]_{j}ds_{j}.$
(5.14)

The discrete boundary terms in (5.14) now have the same form as the continuous
ones in (3.2). By using (5.8) and (5.9) we find that

$(\vec{U},{\vec{L}_{D}})=\sum_{j=1,N}[2(\vec{W}^{-})^{T}\Sigma(\sqrt{|\Lambda^{-}|}\vec{W}^{-}-R\sqrt{\Lambda^{+}}\vec{W}^{+}-S\vec{G})]_{j}ds_{j}.$
(5.15)

The combination of (5.13)-(5.15) leads to the final form of the energy rate

$\dfrac{1}{2}\dfrac{d}{dt}\|\vec{U}\|_{P\otimes
P_{\Omega}}^{2}+\sum_{j=1,N}[\vec{W}^{T}\Lambda\vec{W}+2(\vec{W}^{-})^{T}\Sigma(\sqrt{|\Lambda^{-}|}\vec{W}^{-}-R\sqrt{\Lambda^{+}}\vec{W}^{+}-S\vec{G})]_{j}ds_{j}=0.$
(5.16)

By using (3.7)-(3.9) in Lemma 3.2, the estimates (5.10) and (5.11) follow by
using the same technique that was used for the continuous estimates in the
proof of Lemma 3.2. ∎

##  6 Summary

In this paper we have complemented the stability analysis of nonlinear initial
boundary value problems on skew-symmetric first derivative form initiated in
[1, 2], by adding the analysis of nonlinear boundary conditions. We also
extend the analysis to parabolic and incompletely parabolic problems and
derive these boundary conditions with non-zero data explicitly. The general
procedure is exemplified using: the shallow water equations, the
incompressible Euler equations and Navier-Stokes equations and the
compressible Euler equations. Once the bound of the continuous problem is
obtained, we show how to implement the boundary conditions in a provable
stable way using summation-by-parts formulations and weak boundary procedures.

## Acknowledgment

Jan Nordström was supported by Vetenskapsrådet, Sweden [award no. 2021-05484
VR].

## References

  * [1] J. Nordström, Nonlinear and linearised primal and dual initial boundary value problems: When are they bounded? how are they connected?, Journal of Computational Physics 455 111001 (2022).
  * [2] J. Nordström, A skew-symmetric energy and entropy stable formulation of the compressible Euler equations, Journal of Computational Physics 470 111573 (2022).
  * [3] H.-O. Kreiss, Initial boundary value problems for hyperbolic systems, Commun. Pur. Appl. Math. 23 (3) (1970) 277–298.
  * [4] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Math. Comp. 31 (139) (1977) 629–651.
  * [5] J. Nordström, The use of characteristic boundary conditions for the Navier-Stokes equations, Comput. & Fluids 24 (5) (1995) 609–623.
  * [6] J. Nordström, M. H. Carpenter, Boundary and interface conditions for high-order finite-difference methods applied to the Euler and Navier-Stokes equations, J. Comput. Phys. 148 (2) (1999) 621–645.
  * [7] J. Nordström, M. Svärd, Well-posed boundary conditions for the Navier-Stokes equations, SIAM J. Numer. Anal. 43 (3) (2005) 1231–1255.
  * [8] P. N. Vabishchevich, On the form of the hydrodynamics equations, in: West-East High Speed Flow Field Conference, Moscow, Russia, 19–22 November 2007.
  * [9] W. Rozema, J. Kok, R. Verstappen, A. Veldman, A symmetry-preserving discretisation and regularisation model for compressible flow with application to turbulent channel flow, Journal of Turbulence 15 (6) (2014) 386–410.
  * [10] J. Reiss, J. Sesterhenn, A conservative, skew-symmetric finite difference scheme for the compressible Navier–Stokes equations, Computers and Fluids 101 (2014) 208–219.
  * [11] F. D. Halpern, R. E. Waltz, Anti-symmetric plasma moment equations with conservative discrete counterparts, Physics of Plasmas 25 (6) (2018) 060703.
  * [12] M. Svärd, J. Nordström, Review of summation-by-parts schemes for initial–boundary-value problems, Journal of Computational Physics 268 (2014) 17–38.
  * [13] D. C. D. R. Fernández, J. E. Hicken, D. W. Zingg, Review of summation-by-parts operators with simultaneous approximation terms for the numerical solution of partial differential equations, Computers & Fluids 95 (2014) 171–196.
  * [14] M. Svärd, H. Özcan, Entropy-stable schemes for the Euler equations with far-field and wall boundary conditions, J. Sci. Comput. 58 (1) (2014) 61–89.
  * [15] M. Parsani, M. H. Carpenter, E. J. Nielsen, Entropy stable wall boundary conditions for the three-dimensional compressible Navier–Stokes equations, J. Comput. Phys. 292 (2015) 88–113.
  * [16] M. Svärd, M. H. Carpenter, M. Parsani, Entropy stability and the no-slip wall boundary condition, SIAM Journal on Numerical Analysis 56 (1) (2018) 256–273.
  * [17] M. Svärd, Entropy stable boundary conditions for the Euler equations, Journal of Computational Physics 426 109947 (2021).
  * [18] J. Chan, Y. Lin, T. Warburton, Entropy stable modal discontinuous Galerkin schemes and wall boundary conditions for the compressible Navier-Stokes equations, Journal of Computational Physics 448 110723 (2022).
  * [19] A. Gjesteland, M. Svärd, Entropy stability for the compressible Navier-Stokes equations with strong imposition of the no-slip boundary condition, Journal of Computational Physics 470 111572 (2022).
  * [20] H. Kreiss, Initial boundary value problems for hyperbolic systems, Communications on Pure and Applied Mathematics 23 (3) (1970) 277–298.
  * [21] J. Strikwerda, Initial boundary value problems for incompletely parabolic systems, Commun. Pure Appli. Math. 30 (6) (1977) 797–822.
  * [22] J. Nordström, T. M. Hagstrom, The number of boundary conditions for initial boundary value problems, SIAM Journal on Numerical Analysis 58 (5) (2020) 2818–2828.
  * [23] J. Nordström, A. R. Winters, A linear and nonlinear analysis of the shallow water equations and its impact on boundary conditions, Journal of Computational Physics 463 111254 (2022).
  * [24] R. J. LeVeque, Numerical Methods for Conservation Laws, Birkhäuser, 1992.
  * [25] E. Godlewski, P.-A. Raviart, Numerical Approximation of Hyperbolic Systems of Conservation Laws, Springer, 1996.
  * [26] A. R. Winters, C. Czernik, M. B. Schily, G. J. Gassner, Entropy stable numerical approximations for the isothermal and polytropic Euler equations, BIT 60 (3) (2020) 791–824.
  * [27] R. A. Horn, C. R. Johnson, Matrix Analysis, Cambridge University Press, 2012.
  * [28] D. Arnold, F. Brezzi, B. Cockburn, L. Donatella Marini, Unified analysis of discontinuous Galerkin methods for elliptic problems, SIAM Journal on Numerical Analysis 39 (5) (2001) 1749–1779.
  * [29] J. Nordström, A roadmap to well posed and stable problems in computational physics, J. Sci. Comput. 71 (1) (2017) 365–385.
  * [30] J. Nordström, C. L. Cognata, Energy stable boundary conditions for the nonlinear incompressible Navier–Stokes equations, Math. Comput. 88 (316) (2019) 665–690.
  * [31] J. Nordström, F. Laurén, The spatial operator in the incompressible Navier–Stokes, Oseen and Stokes equations, Computer Methods in Applied Mechanics and Engineering 363 (2020) 112857.
  * [32] J. Oliger, A. Sundström, Theoretical and practical aspects of some initial boundary value problems in fluid dynamics, SIAM J. Appl. Math. 35 (3) (1978) 419–446.
  * [33] C. B. Vreugdenhil, Numerical Methods for Shallow-Water Flow, Vol. 13, Springer, 2013\.
  * [34] G. B. Whitham, Linear and Nonlinear Waves, JWS, 1974.
  * [35] B. Gustafsson, H.-O. Kreiss, J. Oliger, Time dependent problems and difference methods, Vol. 24, JWS, 1995.
  * [36] H.-O. Kreiss, J. Lorenz, Initial-boundary value problems and the Navier-Stokes equations, Vol. 47, SIAM, 1989.
  * [37] T. Lundquist, A. Malan, J. Nordström, A hybrid framework for coupling arbitrary summation-by-parts schemes on general meshes, Journal of Computational Physics 362 (2018) 49 – 68.