2 Continuous analysis

The one-dimensional (1D) SWWEs are

(2.1) $$ \begin{align} \frac{\partial h}{\partial t} + \frac{\partial (uh)}{\partial x} = 0, \quad \frac{\partial (uh)}{\partial t} + \frac{\partial (u^2 h + gh^2/2) }{\partial x} = 0, \end{align} $$

where $x\in \mathbb {R}$ is a spatial variable, $t\ge 0$ is time, $h(x,t)>0$ and $u(x,t)$ are the water depth and the depth averaged fluid velocity, respectively, and $g>0$ is the gravitational acceleration.

To make our analysis tractable, we linearise the SWWEs by substituting $h = H + \widetilde {h}$ and $u = U+\widetilde {u}$ into (2.1), where $\widetilde {h}$ and $\widetilde {u}$ denote perturbations of the constant water depth $H>0$ and fluid velocity U, respectively.

After simplifying, the linearised SWWEs are

(2.2) $$ \begin{align} &\frac{\partial \tilde{h}}{\partial t} + U \frac{\partial \tilde{h}}{\partial x} + H \frac{\partial \tilde{u}}{\partial x} = 0, \quad \frac{\partial \tilde{u}}{\partial t} + g \frac{\partial \tilde{h}}{\partial x} + U \frac{\partial \tilde{u}}{\partial x} = 0. \end{align} $$

Introducing the unknown vector field $\mathbf {p} = \begin {bmatrix} {\widetilde {h}}&{\widetilde {u}} \end {bmatrix}^\top $ , the linear equation (2.2) can be rewritten in a more compact form as

(2.3) $$ \begin{align} \frac{\partial \mathbf{{p}}}{\partial t} = \mathcal{D}\textbf{p}, \quad \mathcal{D} = - {\mathcal{M}}\frac{\partial}{\partial x}, \quad \mathbf{p}=\begin{bmatrix}\widetilde{h} &\widetilde{u} \end{bmatrix}^\top, \quad {\mathcal{M}}=\begin{bmatrix}U&H\\g&U\end{bmatrix}. \end{align} $$

To avoid inconsistencies in the units of the matrix entries, we rescale the variables as follows:

(2.4) $$ \begin{align} \mathbf{q} = \begin{bmatrix} \widetilde{\widetilde{h}}\\ \widetilde{\widetilde{u}} \end{bmatrix} = \begin{bmatrix} {\widetilde{h}}/{H}\\ {\widetilde{u}}/{c} \end{bmatrix} = W\mathbf{p}, \quad W = \begin{bmatrix} 1/H & 0\\ 0 & 1/c \end{bmatrix}. \end{align} $$

First, we multiply (2.3) by W, and introduce the dimensionless variable $\mathbf {q}$ and the symmetric matrix $M = W\mathcal {M}W^{-1}$ . After simplifying and dropping the double-tilde, for simplicity, the dimensionless equation can be written as

(2.5) $$ \begin{align} \frac{\partial \textbf{q}}{\partial t} = D\textbf{q}, \quad {D} = - {M}\frac{\partial}{\partial x}, \quad \mathbf{q}=\begin{bmatrix}h &u \end{bmatrix}^\top, \quad {M}=\begin{bmatrix}U&c\\c&U\end{bmatrix}. \end{align} $$

Note that in (2.5), the constant coefficient matrix ${M}$ is symmetric. This will simplify the following analysis.

We will consider (2.5) in a bounded domain, and augment it with initial and boundary conditions. Let our domain be $\Omega =[0,1]$ and $\Gamma =\{0,1\}$ be the boundary points. We consider the IBVP

(2.6a) $$ \begin{align} &\frac{\partial \textbf{q}}{\partial t} = D\textbf{q} , \quad x\in \Omega, \ t\ge 0, \end{align} $$

(2.6b) $$ \begin{align} &\textbf{q}(x,0) = \textbf{f}(x), \quad x\in \Omega, \end{align} $$

(2.6c) $$ \begin{align} &\mathcal{B}\textbf{q} = \textbf{b}(t), \quad x\in \Gamma, \ t\ge 0, \end{align} $$

where $\mathcal {B}$ is a linear boundary operator, $\textbf {b}$ is the boundary data and $\textbf {f} \in L^2(\Omega )$ is the initial condition. One objective of this study is to investigate the choice of the boundary operator $\mathcal {B}$ which ensures that the IBVP (2.6) is well-posed. To simplify the coming analysis, we will consider zero boundary data $\textbf {b}=0$ , but the results can be extended to nontrivial boundary data $\textbf {b} \ne 0$ . Furthermore, numerical experiments performed later in this paper confirm that the analysis is valid for nonzero boundary data.

Let $\textbf {p}$ and $\textbf {q}$ be real-valued vector functions, and define the standard $L^2(\Omega )$ scalar product and the norm

(2.7) $$ \begin{align} (\textbf{p},\textbf{q}) = \int_\Omega \textbf{p}^\top \textbf{q} \, dx, \quad \Vert \textbf{q} \Vert^2 = (\textbf{q},\textbf{q})>0, \end{align} $$

for all nonzero $\textbf {q} \in \mathbb {R}^2$ .

Definition 2.1. The IBVP (2.6) is well-posed if a unique solution $\textbf {q}$ satisfies

$$ \begin{align*} \Vert \textbf{q}(\cdot,t) \Vert \leq \kappa e^{\nu t} \Vert \textbf{f} \Vert, \quad \Vert \textbf{f} \Vert < \infty \end{align*} $$

for some constants $\kappa>0$ and $\nu \in \mathbb {R}$ independent of $\textbf {f}$ .

The well-posedness of the IBVP (2.6) can be related to the boundedness of the differential operator ${D}$ . We introduce the function space

$$ \begin{align*} \mathbb{V} =\{\textbf{p}| \ \textbf{p}(x)\in \mathbb{R}^2, \ \Vert \textbf{p} \Vert < \infty, \ 0\le x\le 1, \ \{\mathcal{B}\textbf{p}=0, \ x\in \Gamma \} \}. \end{align*} $$

The following two definitions are useful.

Definition 2.2. The operator ${D}$ is said to be semi-bounded in the function space $\mathbb {V}$ if it satisfies

$$ \begin{align*} (\textbf{q},{D}\textbf{q}) \leq \nu \Vert \textbf{q} \Vert^2, \quad \nu \in \mathbb{R}. \end{align*} $$

It is well known that the maximally semi-boundedness of differential operator ${D}$ is a necessary and sufficient condition for the well-posedness of the IBVP (2.6) [Reference Gustafsson, Kreiss and Oliger6].

Thus, to ensure that the IBVP (2.6) is well-posed, we need: (a) the differential operator ${D}$ to be semi-bounded; and (b) the minimal number of boundary conditions such that ${D}$ is maximally semi-bounded.

To begin with, we will show that the differential operator ${D}$ is semi-bounded in $L_2(\Omega )$ .

Proof. We consider $(\textbf {q},{D}\textbf {q})$ and use integration-by-parts, then

$$ \begin{align*} (\textbf{q},{D}\textbf{q}) = -\int_\Omega \textbf{q}^\top {M}\frac{\partial \textbf{q}}{\partial x} \, dx = -\frac{1}{2}\int_\Omega \frac{\partial }{\partial x} (\textbf{q}^\top {M} \textbf{q}) \, dx = - \frac{1}{2}(\textbf{q}^\top {M} \textbf{q})|_0^1. \end{align*} $$

Thus, if the boundary term $ (\textbf {q}^\top {M} \textbf {q})|_0^1 \ge 0$ , then $(\textbf {q},\textbf {D}\textbf {q}) \le 0$ . In particular, the upper bound $(\textbf {q},\textbf {D}\textbf {q}) = 0$ satisfies Definition 2.2 with $\nu =0$ .

The next step will be to derive boundary operators $\{\mathcal {B}\textbf {p}=0, \ x\in \Gamma \}$ with minimal number of boundary conditions such that the boundary term is never negative, $ (\textbf {q}^\top {M} \textbf {q})|_0^1 \ge 0$ .

Noting that the norm is related to the dimensionless mechanical energy $\tilde {E}$ , that is,

$$ \begin{align*} \frac{1}{2}\Vert \textbf{q} \Vert^2={E} := \int_{\Omega} \frac12 (h^2 + u^2) \,dx> 0 \quad \text{for all}\, \textbf{q} \in \mathbb{R}^2\backslash \{\textbf{0}\}. \end{align*} $$

The mechanical energy in the physical units can be recovered through the scaling $\widetilde {E}=c^2H {E}$ .

We introduce the boundary term

$$ \begin{align*} \text{BT}:= -(\textbf{q}^\top {M} \textbf{q})|_{0}^{1}. \end{align*} $$

By using the eigen-decomposition ${M} = S\Lambda S^{T}$ given by

(2.8) $$ \begin{align} S = \frac{1}{\sqrt{2}}\begin{bmatrix} 1&1\\1&-1 \end{bmatrix}, \quad \Lambda = \begin{bmatrix} \lambda_1&0\\0&\lambda_2 \end{bmatrix}, \quad \lambda_1 = U+c, \quad \lambda_2 = U-c, \end{align} $$

with the linear transformation

(2.9) $$ \begin{align} &\begin{bmatrix} w_1\\w_2 \end{bmatrix} = {S}^{\top} \textbf{q} = \dfrac{1}{\sqrt{2}} \begin{bmatrix} h+u \\ h-u \end{bmatrix}, \end{align} $$

the boundary term can be re-written as

(2.10) $$ \begin{align} \text{BT}=(\lambda_1 w_1^2 + \lambda_2 w_2^2)|_{x=0} - (\lambda_1 w_1^2 + \lambda_2 w_2^2)|_{x =1}. \end{align} $$

The number of boundary conditions will depend on the signs of the eigenvalues $\lambda _1$ , $\lambda _2$ , which in turn depend on the magnitude of the flow U relative to the characteristic wave speed c, and are determined by the Froude number $Fr = |U|/c$ . If $0\le Fr <1$ , then $\lambda _1>0$ and $\lambda _2 <0$ for any U. In this case, we have one boundary condition each in the inflow and outflow. For $Fr> 1$ , the sign of the eigenvalues $\lambda _1$ and $\lambda _2$ will take the sign of U. In this case, we have either two boundary conditions on the left if $U>0$ (inflow at $x=0$ ) or two boundary conditions on the right if $U<0$ (inflow at $x=1$ ). The case where $Fr=1$ , we have $\lambda _1>0$ and $\lambda _2 = 0$ if $U>0$ , and $\lambda _1=0$ and $\lambda _2 < 0$ if $U<0$ . That is, we only have one boundary condition at the inflow at $ x = 0$ if $U>0$ or at $x = 1$ if $U<0$ .

Sub-critical flow. The flow is sub-critical when $\operatorname {Fr} < 1$ , which implies $\lambda _1>0$ and $\lambda _2 <0$ . We need one boundary condition at $x = 0$ and another boundary condition at $x = 1$ . Therefore, for the sub-critical flow regime, we always need an inflow boundary condition and an outflow boundary condition for any U. We formulate the boundary conditions by

(2.11) $$ \begin{align} \{\mathcal{B}\textbf{p}=\textbf{b}, \ x\in \Gamma \} \equiv \{ {w_1 - \gamma_0 w_2 = b_1(t), \ x =0}; \ {w_2 - \gamma_N w_1 = b_2(t), \ x = 1}\}, \end{align} $$

where $\gamma _0, \gamma _N \in \mathbb {R}$ are boundary reflection coefficients. The following lemma constrains the parameters $\gamma _0, \gamma _N$ .

Proof. Let $w_1 = \gamma _0 w_2$ at $ x =0$ and $w_2 = \gamma _N w_1$ at $ x =1$ , and consider

$$ \begin{align*} (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_0 - (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_1 = w_2^2 ( \lambda_1\gamma_0^2 + \lambda_2 ) \vert_0 - w_1^2 ( \lambda_1 + \lambda_2 \gamma_N^2) \vert_1. \end{align*} $$

Thus, if $0\leq \gamma _0^2 \leq -\lambda _2/\lambda _1$ and $0\leq \gamma _N^2 \le -\lambda _1/\lambda _2$ , then $ \lambda _1\gamma _0^2 {\kern-1pt}+{\kern-1pt} \lambda _2 {\kern-1pt}\le{\kern-1pt} 0 \ \text {and} \ \lambda _1 {\kern-1pt}+{\kern-1pt} \lambda _2 \gamma _N^2 {\kern-1pt}\ge{\kern-1pt} 0,$ and then

$$ \begin{align*} \text{BT} = w_2^2 ( \lambda_1\gamma_0^2 + \lambda_2 ) \vert_0 - w_1^2 ( \lambda_1 + \lambda_2 \gamma_N^2) \vert_1\le 0.\\[-37pt] \end{align*} $$

Super-critical flow. When $\operatorname {Fr}> 1$ , the flow is super-critical, then $\lambda _1$ and $\lambda _2$ both take the sign of the average flow velocity U. That is, if $U>0$ , then $\lambda _1> 0$ , $\lambda _2> 0$ , and if $U<0$ , then $\lambda _1 < 0$ , $\lambda _2 < 0$ . Thus, when $U>0$ and $\operatorname {Fr}> 1$ , we need two boundary conditions at $x =0$ , and no boundary condition at $x =1$ . Similarly, when $U<0$ and $\operatorname {Fr}> 1$ , we need two boundary conditions at $x =1$ and no boundary conditions at $x =0$ . Therefore, for super-critical flows, there are no outflow boundary conditions for any U. We formulate the boundary conditions by

(2.12a) $$ \begin{align} \{\mathcal{B}\textbf{q}=\textbf{b}, \ x\in \Gamma \} \equiv \{ {w_1 = b_1(t), \ w_2 =b_2(t), \ x =0} \ \text{if } U>0 \text{ and } \operatorname{Fr} > 1 \}, \end{align} $$

(2.12b) $$ \begin{align} \{\mathcal{B}\textbf{q}=\textbf{b}, \ x\in \Gamma \} \equiv \{ {w_1 = b_1(t), \ w_2 =b_2(t), \ x =1} \ \text{if } U <0 \text{ and } \operatorname{Fr}> 1\}. \end{align} $$

Proof. Let $U> 0$ with $\lambda _1> 0$ , $\lambda _2> 0$ if $w_1 = 0, \ w_2 =0$ , at $x =0$ , then

$$ \begin{align*} \text{BT} = (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_0 - (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_1 = - \tfrac{1}{2}(\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_1 \le 0. \end{align*} $$

If $U < 0$ with $\lambda _1 < 0$ , $\lambda _2 < 0$ and $w_1 = 0, \ w_2 =0$ , at $x =1$ , then

$$ \begin{align*} \text{BT} = (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_0 - (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_1 = \tfrac{1}{2} (\lambda_1 w^2 + \lambda_2 w_2^2)\vert_0 \le 0.\\[-37pt] \end{align*} $$

Critical flow. The flow is critical when $\operatorname {Fr} = 1$ . Note that this case is degenerate, since there is only one nonzero eigenvalue, that is, $U>0$ implies $\lambda _1> 0$ , $\lambda _2 =0$ and $U <0$ implies $\lambda _1 = 0$ , $\lambda _2 < 0$ . However, it can also be treated by prescribing only one boundary condition for the system. The location of the boundary condition will be determined by the sign of U, similar to the super-critical flow regime. We prescribe the boundary conditions by

(2.13a) $$ \begin{align} \{\mathcal{B}\textbf{q}=\textbf{b}, \ x\in \Gamma \} \equiv \{ {w_1 = b_1(t), \ x =0} \ \text{if } U>0 \text{ and } \operatorname{Fr} = 1\}, \end{align} $$

(2.13b) $$ \begin{align} \{\mathcal{B}\textbf{q}=0, \ x\in \Gamma \} \equiv \{ w_2 =b_2(t), \ x =1 \ \text{if } U <0 \text{ and } \operatorname{Fr} = 1\}. \end{align} $$

Proof. Let $U> 0$ with $\lambda _1> 0$ , $\lambda _2 = 0$ if $w_1 = 0$ , at $x =0$ ,

$$ \begin{align*} \text{BT} = (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_0 - (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_1 = - \tfrac{1}{2}\lambda_1 w_1^2 \vert_1 \le 0. \end{align*} $$

If $U < 0$ with $\lambda _1 = 0$ , $\lambda _2 < 0$ and $w_2 =0$ , at $x =1$ , then also

$$ \begin{align*} \text{BT} = (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_0 - (\lambda_1 w_1^2 + \lambda_2 w_2^2)\vert_1 = \tfrac{1}{2} \lambda_2 w_2^2\vert_0 \le 0. \end{align*} $$

The proof is complete.

We will conclude this section with the theorem that proves the well-posedness of the IBVP (2.6).

Theorem 2.8. Consider the IBVP (2.6) where the boundary operator $\mathcal {B}\textbf {q}=0$ is defined by (2.11) with $\gamma _0^2 \leq -\lambda _2/\lambda _1$ and $\gamma _N^2 \le -\lambda _1/\lambda _2$ for sub-critical flows, $\operatorname {Fr} < 1$ ; by (2.12) for the super-critical flow regime, $\operatorname {Fr}> 1$ ; and by (2.13) for critical flows, $\operatorname {Fr} = 1$ ; we have the energy estimate

(2.14) $$ \begin{align} \frac{1}{2}\frac{d}{dt}\Vert \textbf{q} \Vert_{W}^2 = \mathrm{BT}\le 0. \end{align} $$

Proof. We use the energy method, that is, from the left, we multiply (2.6a) with $\textbf {q}^\top W$ and integrate over the domain. As above, integration-by-parts gives

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \textbf{q} \Vert_{W}^2 = \bigg(\textbf{q},\frac{\partial \textbf{q}}{\partial t}\bigg)_{W} = (\textbf{q},\textbf{D}\textbf{q})_{W} = \mathrm{BT}. \end{align*} $$

Using Lemmas 2.5–2.7 for each flow regime gives $\mathrm {BT}\le 0$ . Then the proof is complete.

This energy estimate (2.14) is what a stable numerical method should emulate.

3 Numerical scheme

We derive a stable finite-volume method for the IBVP (2.6) encapsulated in the SBP framework. Numerical stability is proved by deriving discrete energy estimates analogous to Theorem 2.8.

3.1 The finite-volume method

To begin, the domain, $\Omega = [0, 1]$ , is subdivided into $N+1$ computational nodes having $ x_{i} = x_{i-1} + \Delta {x}_i$ for $i = 1, 2, \ldots N$ , with $x_0 = 0$ , $\Delta {x}_i> 0$ and $\sum _{i=1}^{N}\Delta {x}_i = 1$ . We consider the control cell $I_i = [x_{i-1/2},x_{i+1/2}]$ for each interior node $1\le i \le N-1$ , and for the boundary nodes $\{x_0, x_N\}$ , the control cells are $I_{0} = [x_0, x_{1/2}]$ and $I_{N} = [x_{N-1/2}, x_{N}]$ , see Figure 1. Note that $|I_i| =\Delta {x}_{i}/2 + \Delta {x}_{i+1}/2$ for the interior nodes $1\le i\le N-1$ , and for the boundary nodes $i \in \{0, N\}$ , we have $|I_0| = \Delta {x}_{1}/2$ and $|I_N| = \Delta {x}_{N}/2$ . The control cells $I_i$ are connected and do not overlap, and $\sum _{i=0}^N |I_i| = \sum _{i=1}^{N}\Delta {x}_i = 1.$

Figure 1 Finite-volume nodes $x_i$ and control cells $I_i$ .

Consider the integral form of (2.2) over the control cells $I_i$ :

$$ \begin{align*} &\dfrac{d}{dt}\int_{I_0} \textbf{p}(x,t) \, dx + {\mathcal{M}} \textbf{p}(x_{1/2},t) - {\mathcal{M}}\textbf{p}(x_{0},t) = 0, \\ &\dfrac{d}{dt}\int_{I_i} \textbf{p}(x,t) \, dx + {\mathcal{M}} \textbf{p}(x_{i+1/2},t) - {\mathcal{M}}\textbf{p}(x_{i-1/2},t) = 0, \quad 1\le i\le N-1,\\ & \dfrac{d}{dt}\int_{I_N} \textbf{p}(x,t) \, dx + {\mathcal{M}} \textbf{p}(x_{N},t) - {\mathcal{M}}\textbf{p}(x_{N-1/2},t) = 0 \end{align*} $$

with $\mathcal {M} = \begin {bmatrix} U&H\\g&U \end {bmatrix}$ and $\textbf {p} = \begin {bmatrix} \tilde {h}&\tilde {u} \end {bmatrix}^\top $ .

We introduce the cell-average

$$ \begin{align*} \bar{\textbf{p}}_i = \frac{1}{|I_i|}\int_{I_i} \textbf{p}(x,t)\,dx, \end{align*} $$

and approximate the PDE flux $\mathcal {M}\textbf {p}$ with the local Lax–Friedrich flux

(3.1) $$ \begin{align} \mathcal{M}\textbf{p}(x_{i+1/2},t) \approx \frac{{\mathcal{M}}\bar{\textbf{p}}_{i+1} + {\mathcal{M}}\bar{\textbf{p}}_i}{2} - \frac{\alpha}2(\bar{\textbf{p}}_{i+1} - \bar{\textbf{p}}_i), \quad \alpha \ge 0, \end{align} $$

and

$$ \begin{align*}{\mathcal{M}}\textbf{p}(x_{0},t) \approx \mathcal{M} \bar{\textbf{p}}_0, \quad {\mathcal{M}}\textbf{p}(x_{N},t) \approx \mathcal{M} \bar{\textbf{p}}_N.\end{align*} $$

The evolution of the cell-average is governed by the semi-discrete system

(3.2a) $$ \begin{align} &|I_0| \frac{d\bar{\textbf{p}_0}}{dt} + {\mathcal{M}}\frac{\bar{\textbf{p}}_{1} - \bar{\textbf{p}}_{0}}{2} - \frac{\alpha}{2} (\bar{\textbf{p}}_{1} - \bar{\textbf{p}}_{0}) = 0, \end{align} $$

(3.2b) $$ \begin{align} &|I_i| \frac{d\bar{\textbf{p}_i}}{dt} + {\mathcal{M}}\frac{\bar{\textbf{p}}_{i+1} - \bar{\textbf{p}}_{i-1}}{2} - \frac{\alpha}{2} (\bar{\textbf{p}}_{i+1} - 2\bar{\textbf{p}}_{i} + \bar{\textbf{p}}_{i-1}) = 0, \ 1\le i\le N-1, \end{align} $$

(3.2c) $$ \begin{align} & |I_N| \frac{d\bar{\textbf{p}_N}}{dt} + {\mathcal{M}}\frac{\bar{\textbf{p}}_{N} -\bar{\textbf{p}}_{N-1}}{2} - \frac{\alpha}{2} (\bar{\textbf{p}}_{N-1} - \bar{\textbf{p}}_{N}) = 0. \end{align} $$

Introducing the discrete solution vector $\bar {\textbf {p}} = [\bar {\textbf {p}}_0, \bar {\textbf {p}}_1, \ldots , \bar {\textbf {p}}_N]^{\top }$ and rewriting (3.2) in a more compact form,

(3.3) $$ \begin{align} (I\otimes{P})\frac{d\bar{\textbf{p}}}{dt} + ({\mathcal{M}}\otimes {Q}) \bar{\textbf{p}} - \frac{\alpha}{2} (I\otimes{A})\bar{\textbf{p}}=0, \end{align} $$

where $\otimes $ denotes the Kronecker product and

$$ \begin{align*} {Q} = \begin{pmatrix} - \frac12 & \frac12 & 0 & \cdots & 0 & 0 & 0 \\ -\frac{1}{2} & 0 & \frac{1}{2} & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & -\frac{1}{2} & 0 & \frac{1}{2} \\ 0 & 0 & 0 & \cdots & 0 & -\frac12 & \frac12 \end{pmatrix}, \quad {A} = \begin{pmatrix} {-1} &{1} & 0 & \cdots & 0 & 0 & 0 \\ 1 & -2 & 1 & \cdots & 0 & 0 & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\ 0 & 0 & 0 & \cdots & 1 & -2 & 1 \\ 0 & 0 & 0 & \cdots & 0 & 1 & -1 \end{pmatrix}, \end{align*} $$

and ${P} = \operatorname {diag} ([|I_0|,|I_1|,\ldots ,|I_N|])$ . The matrix Q is related to the spatial derivative operator, A is a numerical dissipation operator and $\alpha \ge 0$ controls the amount of numerical dissipation applied. Note that A is symmetric and negative semi-definite, that is, $A = A^\top $ and $ \mathbf {p}^\top A \mathbf {p} \leq 0 $ for all $\mathbf {p} \in \mathbb {R}^{N+1}$ . The important stability property of the semi-discrete approximation (3.3) is that the associated discrete derivative operator satisfies the SBP property. To see this, we rewrite (3.3) as

$$ \begin{align*} \frac{d\bar{\textbf{p}}}{dt} + ({\mathcal{M}}\otimes {D}_x) \bar{\textbf{p}} - \frac{\alpha}{2} (I\otimes{P}^{-1}{A})\bar{\textbf{p}}=0, \end{align*} $$

where I is the $2\times 2$ identity matrix and

(3.4) $$ \begin{align} {D}_x = {P}^{-1}{Q}, \quad {Q}+{Q}^\top = \operatorname{diag}([-1,0,\dots,0,1]). \end{align} $$

Equation (3.4) is the so-called SBP property [Reference Gustafsson5, Reference Kreiss, Scherer and Boor7] for the first derivative $d/dx$ , which can be useful in proving numerical stability of the discrete approximation (3.3). Note that we have not enforced any boundary condition yet, the boundary condition (2.6c) will be implemented weakly using penalties.

3.2 Numerical boundary treatment and stability

In this section, we will implement the boundary conditions and prove numerical stability. The boundary conditions are implemented using the SAT method, similar terms used as in [Reference Carpenter, Gottlieb and Abarbanel1]; by appending the boundary operators (2.11)–(2.13) to the right-hand side of (3.3) with penalty weights,

(3.5) $$ \begin{align} (I\otimes{P})\frac{d\bar{\textbf{p}}}{dt} + ({\mathcal{M}}\otimes {Q}) \bar{\textbf{p}} - \frac{\alpha}{2} (I\otimes{A})\bar{\textbf{p}}=\mathrm{SAT}. \end{align} $$

With $\mathbf {e}_0 = [1, 0, \ldots 0]^T$ and $\mathbf {e}_N = [0, 0, \ldots 1]^T$ , the SAT for sub-critical flow is

(3.6) $$ \begin{align} \mathrm{SAT}= -\frac{1}{2}(W^{-1}S\otimes\textbf{I}) \begin{bmatrix} {\tau_{01}} \mathbf{e}_0(\bar{w}_1 - \gamma_0 \bar{w}_2 -b_1(t)) + { \tau_{N1}} \mathbf{e}_N (\bar{w}_2 - \gamma_N \bar{w}_1 - b_2(t))\\ {\tau_{02}} \mathbf{e}_0 (\bar{w}_1 - \gamma_0 \bar{w}_2 - b_1(t)) + {\tau_{N2}} \mathbf{e}_N(\bar{w}_2 - \gamma_N \bar{w}_1 - b_2(t)) \end{bmatrix}, \end{align} $$

and for critical/super-critical flow regimes,

(3.7a) $$ \begin{align} \mathrm{SAT}&= -\frac{1}{2}(W^{-1}S\otimes\textbf{I}) \begin{bmatrix} {\tau_{01}} \mathbf{e}_0(\bar{w}_1-b_1(t)) \\ {\tau_{02}} \mathbf{e}_0 (\bar{w}_2-b_2(t)) \end{bmatrix}, \quad U> 0, \end{align} $$

(3.7b) $$ \begin{align} \mathrm{SAT}&= -\frac{1}{2}(W^{-1}S\otimes\textbf{I}) \begin{bmatrix} {\tau_{N1}} \mathbf{e}_N(\bar{w}_1-b_1(t)) \\ {\tau_{N2}} \mathbf{e}_N (\bar{w}_2-b_2(t)) \end{bmatrix}, \quad U < 0. \end{align} $$

Here, S is the orthonormal eigenvector matrix given in (2.8) and W is the diagonal weight matrix given in (2.4). The coefficients $\tau _{01}$ , $\tau _{02}$ , $\tau _{N1}$ , $\tau _{N2}$ are real penalty parameters to be determined by requiring stability. Note that (3.5) is a consistent semi-discrete approximation of the IBVP (2.6) for all nontrivial choices of the penalty parameters. The semi-discrete approximation (3.5), given that the discrete derivative operator satisfies the SBP property (3.4), is often referred to as the SBP-SAT scheme [Reference Gassner3, Reference Lundquist and Nordström8]. We introduce the discrete weighted $L_2$ -norm

$$ \begin{align*} \Vert \bar{\textbf{q}} \Vert^2:= \bar{\textbf{q}}^T(\mathbf{I}\otimes P)\bar{\textbf{q}} \ge 0 \end{align*} $$

for some weighted matrix $\mathcal {W}$ . The semi-discrete approximation (3.5) is stable if the discrete energy norm $\Vert \bar {\textbf {q}} \Vert ^2$ is nonincreasing with time. We will now prove the stability of the semi-discrete approximation (3.5) for sub-critical flows.

Theorem 3.1. Consider the semi-discrete finite-volume approximation (3.5) with the SAT (3.6) and $\textbf {b} =0$ for sub-critical flow regimes, where $\lambda _1>0$ , $\lambda _2 < 0$ and $\gamma _0^2 \le -\lambda _2/\lambda _1$ , $\gamma _N^2 \le -\lambda _1/\lambda _2$ . If the penalty parameters are chosen such that

$$ \begin{align*}\tau_{01} = \lambda_1,\quad \tau_{02} = \gamma_0\lambda_1; \quad \tau_{N2} = -\lambda_2, \quad \tau_{N1} = -\gamma_N\lambda_2,\end{align*} $$

then

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2\le 0 \quad \text{for all}\ t\ge 0. \end{align*} $$

Proof. We use the energy method, that is, from the left, we multiply (3.5) with $(W\otimes \mathbf {I})$ , and using identity $\bar {\mathbf {q}} = W\bar {\mathbf {p}}$ and $M=W\mathcal {M}W^{-1}$ ,

$$ \begin{align*} (I\otimes{P})\frac{d\bar{\textbf{q}}}{dt} + ({{M}}\otimes {Q}) \bar{\textbf{q}} - \frac{\alpha}{2} (I\otimes{A})\bar{\textbf{q}}=(W\otimes\textbf{I})\text{SAT}.\end{align*} $$

We multiply this equation with $\mathbf {q}^\top $ , add its transpose and then simplify further, then:

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2 + \frac{1}{2} \bar{\textbf{q}}^T({M}\otimes ({Q}+{Q}^T)) \bar{\textbf{q}} - \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}=\bar{\textbf{q}}^T(W\otimes \mathbf{I})\mathrm{SAT}. \end{align*} $$

Using the SBP property (3.4) and the eigen-decomposition of ${M}$ ,

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2 - \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}= \frac{1}{2}\mathrm{BT}_{\text{num}}, \end{align*} $$

where

$$ \begin{align*} \mathrm{BT}_{\text{num}} &= (\lambda_1 \bar{w}_1^2 + \lambda_2 \bar{w}_2^2 -(\tau_{01} \bar{w}_1(\bar{w}_1 - \gamma_0 \bar{w}_2) + \tau_{02} \bar{w}_2(\bar{w}_1 - \gamma_0 \bar{w}_2)))|_{i=0} \\ &\quad - (\lambda_1 \bar{w}_1^2 + \lambda_2 \bar{w}_2^2 +(\tau_{N1} \bar{w}_1(\bar{w}_2 - \gamma_N \bar{w}_1 ) + \tau_{N2} \bar{w}_2(\bar{w}_2 - \gamma_N \bar{w}_1 )))|_{i=N}. \end{align*} $$

Thus, if $\tau _{01} = \lambda _1,\, \tau _{02} = \gamma _0\lambda _1; \, \tau _{N2} = -\lambda _2, \, \tau _{N1} = -\gamma _N\lambda _2,$ then

$$ \begin{align*} \mathrm{BT}_{\text{num}} &= ( \lambda_2 + \lambda_1\gamma_0^2)\bar{w}_2^2|_{i=0} - ( \lambda_1 + \lambda_2\gamma_N^2)\bar{w}_1^2|_{i=N}. \end{align*} $$

Since $\lambda _1> 0$ , $\lambda _2 < 0$ and

$$ \begin{align*} (\lambda_2 + \lambda_1\gamma_0^2)\le 0 \!\iff\! \gamma_0^2 \le -\lambda_2/\lambda_1; \quad (\lambda_1 + \lambda_2\gamma_N^2) \ge 0 \!\iff\! \gamma_N^2 \le -\lambda_1/\lambda_2, \end{align*} $$

then it must be true that $\mathrm {BT}_{\text {num}} \le 0$ . Since A is negative semi-definite, then for $\alpha \ge 0$ , we have ${\alpha }/{2} \bar {\textbf {q}}^T(\mathbf {I}\otimes {A})\bar {\textbf {q}}\leq 0$ and

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2 = \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}+ \mathrm{BT}_{\text{num}} \le 0. \end{align*} $$

This completes the proof.

The next theorem will prove the stability of the semi-discrete approximation (3.5) for super-critical flows.

Theorem 3.2. Consider the semi-discrete finite-volume approximation (3.5) with the SAT (3.7) and $\textbf {b}=0$ for super-critical flows. If the penalty parameters are chosen such that $\tau _{01} \ge \lambda _1,\, \tau _{02} \ge \lambda _2; \, \tau _{N1} \ge -\lambda _1, \, \tau _{N2} \ge -\lambda _2, $ then

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2\le 0 \quad \text{for all}\ t\ge 0. \end{align*} $$

Proof. As above, the energy method with the SBP property (3.4) and the eigen-decomposition of ${M}$ yields

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2- \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}= \mathrm{BT}_{\text{num}}, \end{align*} $$

where

$$ \begin{align*} \mathrm{BT}_{\text{num}} &= ((\lambda_1-\tau_{01}) \bar{w}_1^2 + (\lambda_2-\tau_{02}) \bar{w}_2^2)|_{i=0} - (\lambda_1 \bar{w}_1^2 + \lambda_2 \bar{w}_2^2)|_{i=N}, \quad U> 0,\\ \mathrm{BT}_{\text{num}} &= (\lambda_1 \bar{w}_1^2 + \lambda_2 \bar{w}_2^2)|_{i=0} -((\lambda_1+\tau_{N1}) \bar{w}_1^2 + (\lambda_2+\tau_{N2}) \bar{w}_2^2)|_{i=N}, \quad U< 0. \end{align*} $$

Therefore, if $ \tau _{01} \ge \lambda _1, \tau _{02} \ge \lambda _2; \tau _{N1} \ge -\lambda _1, \tau _{N2} \ge -\lambda _2,$ then we have $\mathrm {BT}_{\text {num}} \le 0$ . Noting that $\alpha \ge 0$ and, as previous, we have ${\alpha }/{2} \bar {\textbf {q}}^T(\mathbf {I}\otimes {A})\bar {\textbf {q}}\leq 0$ , which gives us

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2 = \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}+ \mathrm{BT}_{\text{num}} \le 0, \end{align*} $$

so the proof is complete.

Finally, we will prove the stability of the semi-discrete approximation (3.5) for critical flows.

Theorem 3.3. Consider the semi-discrete finite-volume approximation (3.5) with the SAT (3.7) and $\textbf {b}=0$ for critical flows. If the penalty parameters are chosen such that $\tau _{01} \ge \lambda _1,\, \tau _{02} =0; \, \tau _{N1} = 0, \, \tau _{N2} \ge -\lambda _2 $ , then

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2\le 0 \quad \text{for all}\ t\ge 0. \end{align*} $$

Proof. The zero penalties ensure consistency of the SAT, that is, $\tau _{02} =0$ and $\tau _{N1} = 0$ give

$$ \begin{align*} \mathrm{SAT}&= -\frac{1}{2}(W^{-1}S\otimes\textbf{I}) \begin{bmatrix} {\tau_{01}} H\mathbf{e}_0\bar{w}_1 \\ 0 \end{bmatrix}, \quad U> 0,\\ \mathrm{SAT}&= -\frac{1}{2}(W^{-1}S\otimes\textbf{I}) \begin{bmatrix} 0 \\ {\tau_{N2}}g \mathbf{e}_N \bar{w}_2 \end{bmatrix}, \quad U < 0. \end{align*} $$

Again the energy method with the SBP property (3.4) and the eigen-decomposition of ${M}$ yield

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2- \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}= \mathrm{BT}_{\text{num}}, \end{align*} $$

where

$$ \begin{align*} &\mathrm{BT}_{\text{num}} = (\lambda_1-\tau_{01}) \bar{w}_1^2|_{i=0} - \lambda_1 \bar{w}_1^2 |_{i=N}, \quad U> 0, \ \lambda_1 >0, \ \lambda_2 =0, \\ &\mathrm{BT}_{\text{num}} = \lambda_2 \bar{w}_2^2|_{i=0} -(\lambda_2+\tau_{N2}) \bar{w}_2^2|_{i=N}, \quad U< 0, \ \lambda_1 =0, \ \lambda_2 <0. \end{align*} $$

Therefore, if $\tau _{01} \ge \lambda _1$ and $ \tau _{N2} \ge -\lambda _2, $ then we have $\mathrm {BT}_{\text {num}} \le 0$ . Using the fact that $\alpha \ge 0$ and ${\alpha }/{2} \bar {\textbf {q}}^T(\mathbf {I}\otimes {A})\bar {\textbf {q}}\leq 0$ again gives the result that we wanted:

$$ \begin{align*} \frac{1}{2}\frac{d}{dt}\Vert \bar{\textbf{q}} \Vert^2 = \frac{\alpha}{2} \bar{\textbf{q}}^T(\mathbf{I}\otimes{A})\bar{\textbf{q}}+ \mathrm{BT}_{\text{num}} \le 0. \end{align*} $$

This completes the proof.

4 Numerical experiments

In this section, we perform numerical experiments to verify the analysis undertaken in the previous sections. Similar to the theoretical analysis, the numerical experiments cover the three flow regimes, namely the sub-critical, critical and super-critical flow regimes. We used $H=1$ m, $g=9.8$ m/s $^2$ , and $U \in \{\tfrac 12 \sqrt {gH}, \sqrt {gH}, 2 \sqrt {gH}\}$ , which correspond to the three different flow regimes. The interval of interest is $[0,L]$ with $L> 0$ . Note that $U>0$ so that $x=0$ is the inflow boundary and $x = L$ is the outflow boundary. The locations and the number of boundary conditions required are given in Table 1, and the explicit forms of the boundary conditions considered here are given in Table 2, where $g_1(t)$ and $g_2(t)$ are the boundary data.

Table 1 The number and location of the boundary condition in all regime. The boundary at $x=0$ ( $x=1$ ) is inflow (outflow) boundary if $U>0$ and outflow (inflow) boundary if $U<0$ .

Table 2 Transmissive boundary conditions in all regimes with $U>0$ .

The semi-discrete system (3.5) is integrated in time using the classical fourth-order accurate explicit Runge–Kutta method with the time step

$$ \begin{align*} \Delta t = \mathrm{Cr} \, \frac{\Delta x}{|U|+\sqrt{gH}}, \quad \mathrm{Cr} = 0.25, \end{align*} $$

where $\mathrm {Cr}$ is the Courant–Friedrichs–Lewy number, $\Delta x = L/N$ is the uniform cell width and N is the number of finite-volume cells. We will consider a centred numerical flux with $\alpha =0$ and the local Lax–Friedrich’s numerical flux (3.1) with $\alpha>0$ , and verify numerical accuracy. Note that the semi-discrete approximation is energy stable for all $\alpha \ge 0$ .

Nonhomogeneous boundary data. We consider zero initial conditions, that is, $u(x,0)=0$ and $h(x,0)=0$ , and send a wave into the domain through the inflow boundary at $x=0$ . We will consider specifically $g_1(t)\ne 0$ and $g_2(t)=0$ for the boundary conditions given in Table 2, so that the corresponding IBVP has the exact solution

$$ \begin{align*}\widetilde{h}(x, t) = g_1\bigg(t-\frac{x}{U+\sqrt{gH}}\bigg), \quad \widetilde{u}(x, t) = \frac{1}{\sqrt{H/g}}g_1\bigg(t-\frac{x}{U+\sqrt{gH}}\bigg). \end{align*} $$

We will consider a smooth boundary data given by

$$ \begin{align*} g_1(t)= \begin{cases}(\sin(\pi t))^ 4 & \text { if } 0\leq t \leq 1, \\ 0 & \text { otherwise},\end{cases} \quad g_2(t) = 0 \quad \text{for all} \ t\ge 0, \end{align*} $$

and non-smooth boundary data given by

$$ \begin{align*} g_1(t)= \begin{cases}1 & \text { if } 0< t \leq 1, \\ 0 & \text { otherwise},\end{cases} \quad g_2(t) = 0 \quad \text{for all}\ t\ge 0. \end{align*} $$

The boundary data for the boundary conditions in Table 2 can be rewritten as $b_1(t)$ and/or $b_2(t)$ , and in terms of $w_1=b_1(t)=\sqrt {2}/Hg_1(t)$ and $w_2=b_2(t)=\sqrt {2}/Hg_2(t)$ for the given boundary. In the sub-critical case, by using the linear transformation (2.9), the boundary condition can be rewritten in the form (2.11) with $\gamma _0 = 0$ and $\gamma _N = 0$ .

Using the fact that $\lambda _1$ and $\lambda _2$ have different signs, we have $\gamma _0^2 \leq -\lambda _2/\lambda _1$ and $\gamma _N^2 \le -\lambda _1/\lambda _2$ for all $\operatorname {Fr} < 1$ , that is, the condition of Lemma 2.5 is satisfied.

For the critical flow regime, we have $\operatorname {Fr} = 1$ and $\lambda _2=0$ . Only one boundary condition is imposed at the inflow, $\bar {w_1} = b_1(t)=(\sqrt {2}/H)g_1(t)$ . This condition can be rewritten to match the condition in Theorem 3.3 by using the linear transformation (2.9) and the fact that $U^2 = gH$ .

For the super-critical flow regime, two boundary conditions are imposed at the inflow boundary. That is, $\bar {w}_1 = b_1(t)=(\sqrt {2}/H)g_1(t)$ , $\bar {w}_2=b_2(t)=(\sqrt {2}/H)g_2(t)$ as the boundary condition at $x = 0$ . These boundary conditions are equivalent to (2.12a).

The boundary data will generate a pulse from the left boundary at $x =0$ , which will propagate through the domain and leave the domain through $x = L$ .

Figure 2 shows the snapshot of the sub-critical flow solutions at $t = 3.02$ s for both smooth and nonsmooth boundary data, with $\alpha = 0$ and $\alpha = 0.15 \times (U+\sqrt {gH})>0$ . In the plots, we have scaled the horizontal axis by the wave speed $(U+\sqrt {gH})$ so that the solution is spatially invariant for all flow regimes. Note that for the smooth pulse, the numerical solution matches the exact solution excellently well for $\alpha = 0$ and $\alpha = 0.15 \times (U+\sqrt {gH})>0$ , although with $\alpha = 0.15 \times (U+\sqrt {gH})>0$ , the peak of the numerical is slightly dissipated. For the nonsmooth pulse, when $\alpha =0$ , the propagation speed of the pulse is well approximated by the numerical solution. However, there are numerical oscillations generated by the propagating discontinuities. When ${\alpha = 0.15 \times (U+\sqrt {gH})>0}$ , the numerical solution is nonoscillatory, but the discontinuous edges of the solutions are smeared.

Figure 2 The snapshots of the numerical and exact solutions with $\Delta x = L\times 2^{-11}$ m at time $t=3.02$ s for a sub-critical flow regime with smooth and nonsmooth boundary data. For the smooth boundary data, the numerical solution matches the exact solution well for $\alpha = 0$ and $\alpha = 0.15 \times (U+\sqrt {gH})>0$ . Note, however, with $\alpha = 0.15 \times (U+\sqrt {gH})>0$ , the peak of the numerical solution is slightly dissipated. For the nonsmooth boundary data, when $\alpha =0$ , the propagation speed of the pulse is well approximated by the numerical solution. However, there are numerical oscillations generated by the propagating discontinuities. When $\alpha = 0.15 \times (U+\sqrt {gH})>0$ , the numerical solution is nonoscillatory, but propagating discontinuities are smoothed out.

The evolution of the numerical solutions and the exact solutions, at all flow regimes, are shown in Figure 3 for the smooth pulse and in Figure 4 for the nonsmooth pulse. The pulses enter the domain through the inflow boundary at $x=0$ and leave the domain through the outflow at $x=L= (U+\sqrt {gH})\times 5$ . Note that because of the re-scaling of the x-axis to $x/(U+\sqrt {gH})$ , the solutions are invariant for all three flow regimes.

Figure 3 The evolution of the numerical solutions and the exact solutions at all the three flow regimes with smooth boundary data, $\Delta x = L\times 2^{-11}$ m and $\alpha =0$ . The solutions enter the domain through the inflow boundary at $x=0$ and leave the domain through the outflow at $x=L= (U+\sqrt {gH})\times 5$ . Note that because of the re-scaling of the x-axis to $x/(U+\sqrt {gH})$ , the solutions are invariant for all three flow regimes.

Figure 4 The evolution of the numerical solutions and the exact solutions at all three flow regimes with nonsmooth boundary data, $\Delta x = L\times 2^{-11}$ m and $\alpha =0.15\times (U+\sqrt {gH})>0$ . The discontinuous solutions enter the domain through the inflow boundary at $x=0$ and leave the domain through the outflow at $x=L= (U+\sqrt {gH})\times 5$ . Note that because of the re-scaling of the x-axis to $x/(U+\sqrt {gH})$ , the solutions are invariant for all three flow regimes.

Convergence test. Here, we verify the convergence properties of the numerical method. We will use the method of the manufactured solution [Reference Roache11]. That is, we force the system to have the exact smooth solution

(4.1) $$ \begin{align} h(x,t) = \cos(2\pi t) \sin(6 \pi x ), \quad u(x,t) = \sin(2\pi t) \cos(4 \pi x ). \end{align} $$

The initial conditions $h(x,0)$ , $u(x,0)$ and the boundary data $g_1(t)$ and $g_2(t)$ are chosen to match the analytical solution (4.1). We compute the numerical solution on a sequence of increasing number of finite-volume cells, $N=64,128,256,512,1024,2048$ . The $L_2$ -error and convergence rates of the error are shown in Figure 5 and also presented in Table 3. We have performed numerical experiments with no dissipation $\alpha = 0$ and with numerical dissipation set at $\alpha = 0.05$ . From Table 3, we see that the method is second-order accurate $O(\Delta {x}^2)$ when $\alpha =0$ , and first-order accurate $O(\Delta {x})$ when $\alpha>0$ . These are in agreement with the theory.

Figure 5 The error and convergence of the error at final time $t=0.1$ using the manufactured solution for all flow regimes.

Table 3 The error and convergence of the error at final time $t=0.1$ using the manufactured solution for all flow regimes.