Moment-based multi-resolution HWENO scheme for hyperbolic conservation laws

Jiayin Li111School of Mathematical Sciences, Xiamen University, Xiamen, Fujian
361005, P.R. China. E-mail<EMAIL_ADDRESS>Chi-Wang Shu222Division of
Applied Mathematics, Brown University, Providence, RI 02912. E-mail: Chi-
<EMAIL_ADDRESS>and Jianxian Qiu333School of Mathematical Sciences and
Fujian Provincial Key Laboratory of Mathematical Modeling and High-Performance
Scientific Computing, Xiamen University, Xiamen, Fujian 361005, P.R. China.
E-mail<EMAIL_ADDRESS>
Abstract

In this paper, a high-order moment-based multi-resolution Hermite weighted
essentially non-oscillatory (HWENO) scheme is designed for hyperbolic
conservation laws. The main idea of this scheme is derived from our previous
work [J. Comput. Phys., 446 (2021) 110653], in which the integral averages of
the function and its first order derivative are used to reconstruct both the
function and its first order derivative values at the boundaries. However, in
this paper, only the function values at the Gauss-Lobatto points in the one or
two dimensional case need to be reconstructed by using the information of the
zeroth and first order moments. In addition, an extra modification procedure
is used to modify those first order moments in the troubled-cells, which leads
to an improvement of stability and an enhancement of resolution near
discontinuities. To obtain the same order of accuracy, the size of the stencil
required by this moment-based multi-resolution HWENO scheme is still the same
as the general HWENO scheme and is more compact than the general WENO scheme.
Moreover, the linear weights can also be any positive numbers as long as their
sum equals one and the CFL number can still be 0.6 whether for the one or two
dimensional case. Extensive numerical examples are given to demonstrate the
stability and resolution of such moment-based multi-resolution HWENO scheme.

Key Words: Moment-based scheme; Multi-resolution scheme; HWENO scheme;
Hyperbolic conservation laws; KXRCF troubled-cell indicator; HLLC-flux.
AMS(MOS) subject classification: 65M60, 35L65

##  1 Introduction

In this paper, a high-order moment-based multi-resolution Hermite weighted
essentially non-oscillatory (HWENO) scheme is designed for hyperbolic
conservation laws

$\displaystyle\left\\{\begin{aligned} &u_{t}+\nabla\cdot f(u)=0,\\\
&u(x_{1},...,x_{d},0)=u_{0}(x_{1},...,x_{d}).\end{aligned}\right.$ (1.1)

We concentrate our attention on the one and two dimensional cases ($d=1$ or
2), and in these cases we denote $x_{1}$ as $x$ and $x_{2}$ as $y$.

Conservation laws arise from the physical law that the conservative quantity
in any control body can change only due to the flux passing through its
boundaries, which naturally hold for many fundamental physical quantities,
such as the mass, momentum, energy and so on. Such conservation laws are
widely used in a broad spectrum of disciplines where wave motion or advective
transport is important: gas dynamics, acoustics, elastodynamics, optics,
geophysics, and biomechanics, to name but a few.

The differential equation (1.1) can be derived from the integral equation by
simple manipulations provided that the conservative quantity and its
corresponding flux are sufficiently smooth. This proviso is important because
in practice many interesting solutions are not smooth, but contain
discontinuities such as shock waves. A fundamental feature of nonlinear
conservation laws is that discontinuities can easily develop spontaneously
even from smooth initial data, and must be dealt with carefully both
mathematically and computationally. At a discontinuity in the conservative
quantity, the differential equation does not hold in the classical sense and
it is important to remember that the integral form of the conservation laws
does continue to hold which is more fundamental. This is also why we choose
conservative schemes, such as the finite volume method considered in this
paper, which is based on the integral form of the conservation laws.

Since conservation laws have a very wide range of applications and it is
almost impossible in general to get their exact solutions, many scholars have
explored and proposed a series of numerical methods and are still trying to
improve the performance of these algorithms. In 1994, Liu et al. proposed the
first finite volume WENO scheme in [17], and then, in 1996, Jiang and Shu
improved this WENO scheme to fifth order and to conservative finite difference
formulation (which is more efficient in multi-dimensions), and gave a general
definition of the smoothness indicators and nonlinear weights in [12]. The
methodology of such WENO schemes is to use a nonlinear convex combination of
all the candidate stencils to improve the order of accuracy in smooth regions
without destroying the non-oscillatory behavior near discontinuities. This is
also the difference of such WENO schemes from the ENO schemes in [10, 22, 23],
which only choose the locally smoothest stencil automatically among all the
central and biased spatial stencils. Thereafter, different kinds of WENO
schemes have been developed in, e.g. [3, 18, 6, 8, 28, 29, 2, 1]. Although
these WENO schemes work well for most of the problems we encountered, there is
still room for improvement. For example, if we want to obtain a higher order
scheme, we must further expand the stencil. This will make our scheme not very
compact and will also bring trouble to the processing of the boundary
conditions. In order to overcome this drawback, Qiu and Shu proposed the first
HWENO scheme for one-dimensional problems in [19] and then, in 2005, they
extended this HWENO scheme to two-dimensional problems in [20], where two
different stencils were used to reconstruct the function and its first order
derivative values, respectively. The main difference of such HWENO scheme from
the WENO scheme is that both the function and its first order derivative
values are evolved in time and used in the reconstruction process, not like
the WENO scheme in which only the function values are evolved and used. This
allows the HWENO scheme to obtain the same order of accuracy as the WENO
scheme with relatively narrower stencils. But there occurs a new issue, that
is this HWENO scheme is not stable enough when simulating certain severe
problems with strong discontinuities, including the double Mach and forward
step problems. This difficulty is largely due to the fact that the first order
derivative values may become very large near these discontinuities. Thus, the
stability issue may arise, if these large values are used straightforwardly
without any modification. Driven by the goal of solving this issue, many
effective methods based on the idea of the original HWENO scheme have emerged.
For example, the scheme with a new procedure to reconstruct the first order
derivative values by Zhu and Qiu in [27] in 2008, the scheme with an
additional positivity-preserving limiter by Liu and Qiu in [15, 16] in 2015
and 2016, the scheme with a troubled-cell indicator to modify the first order
moments near the discontinuities before the reconstruction algorithm by Zhao
et al. in [25, 26] in 2020, the scheme with a hierarchy of nested central
spatial stencils by Li et al. in [14] in 2021 and so on, have been developed.

In 2018 and 2019, Zhu and Shu proposed a new type of high-order finite
difference and finite volume multi-resolution WENO schemes in [28, 29], which
only use the information on a hierarchy of nested central spatial stencils and
can not only obtain high order accuracy in smooth regions but also allow the
accuracy near discontinuities to degrade gradually. Following the idea of
these multi-resolution WENO schemes, we designed a new type of high-order
finite volume and finite difference multi-resolution HWENO schemes in [14] in
2021, for which only the function values need to be reconstructed by the HWENO
schemes, and the first order derivative values are obtained directly from the
polynomial with the highest degree in the hierarchy. This can improve the
resolution of the scheme but does not have much effect on its stability, since
the first two layers in the hierarchy do not contain the information of the
first order derivative. In this paper, the function and its first order moment
values are used in our reconstruction algorithm, and only the function values
are needed to be reconstructed according to our control equations. Thus, the
coupling between the function and its first order moment is stronger.

There are two issues to be addressed. The first one is that, although the
bigger the linear weights are for higher degree polynomials, the steeper the
shock transitions become near the discontinuities, the gap between these
linear weights cannot be too large, otherwise the corresponding nonlinear
weights will still be too close to the linear weights near the
discontinuities. This will cause problems, since the higher degree
polynomials, which require the information of the first order derivative or
moment, account for too much in the final reconstruction polynomial, but such
first order derivative or moment values may be very large near the
discontinuities. The other one is that with the suitable choice of the linear
weights, the order of the final reconstruction polynomial will degrade
gradually near the discontinuities until it drops to first order, which will
smooth out these discontinuities to a certain extent. Guided by the idea of
Zhao et al. in [25, 26], we first perform the reconstruction algorithm, and
then modify the first order moments of the troubled-cells and repeat the
reconstruction algorithm for these troubled-cells to update the corresponding
Gauss-Lobatto point values. After such a modification procedure, the
proportion of the last two layers will become much higher than that of the
first two layers in the reconstruction process, thus the resolution near the
discontinuities can be increased significantly. In the meantime, this
modification procedure can also improve the stability of the scheme by
reducing the magnitude of the first order moments near the discontinuities. In
order to better adapt to our high order scheme, we choose the HLLC-flux
(Harten-Lax-van Leer-contact flux) to be our numerical flux, which is an
approximate Riemann solver by assuming that there are four states in the
transition from the left to the right states, thus staying closer to the real
physical situation.

The organization of this paper is as follows: In Section 2, we describe the
reconstruction procedure of the moment-based multi-resolution HWENO scheme for
hyperbolic conservation laws in the one and two dimensions in detail. In
Section 3, we propose a number of numerical examples to illustrate the
accuracy and resolution of our HWENO scheme. Concluding remarks are given in
Section 4.

##  2 Moment-based multi-resolution HWENO scheme

In this section, we describe the reconstruction procedure of the moment-based
multi-resolution HWENO scheme for the one and two dimensional hyperbolic
conservation laws, which has sixth order of accuracy in smooth regions and
high resolution near discontinuities. Here sixth order is simply taken as an
example, arbitrarily high order HWENO schemes can be designed following the
same lines.

###  2.1 One dimensional case

In this subsection, we first consider the following relatively simple one-
dimensional hyperbolic conservation laws

$\displaystyle\left\\{\begin{aligned} &u_{t}+f(u)_{x}=0,\\\
&u(x,0)=u_{0}(x).\end{aligned}\right.$ (2.1)

For simplicity, the computational domain is divided by a uniform mesh
$I_{i}=[x_{i-1/2},x_{i+1/2}]$ with the uniform mesh size $\triangle
x=x_{i+1/2}-x_{i-1/2}$, and the corresponding cell center is denoted by
$x_{i}=\frac{1}{2}(x_{i-1/2}+x_{i+1/2})$.

Firstly, we multiply the governing equation (2.1) by $\frac{1}{\triangle x}$
and $\frac{1}{\triangle x}\frac{x-x_{i}}{\triangle x}$, respectively. Then, we
integrate the resulting equations over the target cell $I_{i}$ and perform the
integration by parts, obtaining the following equations

$\displaystyle\left\\{\begin{aligned} &\frac{1}{\triangle
x}\int_{I_{i}}u_{t}dx=-\frac{1}{\triangle
x}\Big{[}f\big{(}u(x_{i+1/2},t)\big{)}-f\big{(}u(x_{i-1/2},t)\big{)}\Big{]},\\\
&\frac{1}{\triangle x}\int_{I_{i}}u_{t}\frac{x-x_{i}}{\triangle
x}dx=-\frac{1}{2\triangle
x}\Big{[}f\big{(}u(x_{i+1/2},t)\big{)}+f\big{(}u(x_{i-1/2},t)\big{)}\Big{]}+\frac{1}{(\triangle
x)^{2}}\int_{I_{i}}f(u)dx.\end{aligned}\right.$ (2.2)

Next, we define the zeroth and first order moments as follows

$\displaystyle\left\\{\begin{aligned} &\overline{u}_{i}(t)=\frac{1}{\triangle
x}\int_{I_{i}}u(x,t)dx,\\\ &\overline{v}_{i}(t)=\frac{1}{\triangle
x}\int_{I_{i}}u(x,t)\frac{x-x_{i}}{\triangle x}dx,\end{aligned}\right.$ (2.3)

and then, we swap the spatial integration and time derivation and approximate
the flux by a numerical flux, obtaining the following semi-discrete
conservative scheme

$\displaystyle\left\\{\begin{aligned}
&\frac{d\overline{u}_{i}(t)}{dt}=-\frac{1}{\triangle
x}(\hat{f}_{i+1/2}-\hat{f}_{i-1/2}),\\\
&\frac{d\overline{v}_{i}(t)}{dt}=-\frac{1}{2\triangle
x}(\hat{f}_{i+1/2}+\hat{f}_{i-1/2})+\frac{1}{\triangle
x}F_{i}(u).\end{aligned}\right.$ (2.4)

Here, the numerical flux is taken to be the HLLC flux
$\hat{f}_{i+1/2}=\hat{f}^{HLLC}(u_{i+1/2}^{-},u_{i+1/2}^{+})$

$\displaystyle\hat{f}^{HLLC}(u_{L},u_{R})=\left\\{\begin{aligned}
&f_{L},\,\,\,0\leq s_{L},\\\ &f_{L}^{*},\,\,\,s_{L}\leq 0\leq s^{*},\\\
&f_{R}^{*},\,\,\,s^{*}\leq 0\leq s_{R},\\\ &f_{R},\,\,\,s_{R}\leq
0,\end{aligned}\right.$ (2.5)

where $f_{L/R}=f(u_{L/R})$,
$f_{L/R}^{*}=f_{L/R}+s_{L/R}*(u_{L/R}^{*}-u_{L/R})$;
$s_{L}=\mu_{L}-c_{L}*coef_{L}$, $s_{R}=\mu_{R}+c_{R}*coef_{R}$,

$\displaystyle coef_{L/R}=\left\\{\begin{aligned} &1,\,\,\,p^{*}\leq
p_{L/R},\\\
&\sqrt{1+\frac{(\gamma+1)*(p^{*}/p_{L/R}-1)}{2\gamma}},\,\,\,otherwise,\end{aligned}\right.$
(2.6) $\displaystyle
u_{L/R}^{*}=\rho_{L/R}*\frac{s_{L/R}-\mu_{L/R}}{s_{L/R}-s^{*}}*\left(\begin{array}[]{c}1\\\
s^{*}\\\
\frac{E_{L/R}}{\rho_{L/R}}+(s^{*}-\mu_{L/R})*\left(s^{*}+\frac{p_{L/R}}{\rho_{L/R}*(s_{L/R}-\mu_{L/R})}\right)\end{array}\right),$
(2.10)

where $p^{*}=\frac{p_{L}+p_{R}+(\mu_{L}-\mu_{R})*(c^{*}*\rho^{*})}{2}$,
$s^{*}=\frac{\mu_{L}+\mu_{R}+(p_{L}-p_{R})/(c^{*}*\rho^{*})}{2}$;
$\rho^{*}=\frac{\rho_{L}+\rho_{R}}{2}$, $c^{*}=\frac{c_{L}+c_{R}}{2}$. The
integral term $F_{i}(u)$ is approximated by a four-point Gauss-Lobatto
integration

$F_{i}(u)=\frac{1}{\triangle
x}\int_{I_{i}}f(u)dx\approx\sum_{l=1}^{4}\omega_{l}f\big{(}u(x_{l}^{GL},t)\big{)},$
(2.11)

where the weights are $\omega_{1}=\omega_{4}=\frac{1}{12}$,
$\omega_{2}=\omega_{3}=\frac{5}{12}$ and the quadrature points on the target
cell $I_{i}$ are

$\displaystyle
x_{1}^{GL}=x_{i-1/2},x_{2}^{GL}=x_{i-\sqrt{5}/10},x_{3}^{GL}=x_{i+\sqrt{5}/10},x_{4}^{GL}=x_{i+1/2},$

where $x_{i+\xi}=x_{i}+\xi\triangle x$.

Finally, we present the spatial reconstruction procedure of the point values
$\\{u_{i-1/2}^{+}$, $u_{i-\sqrt{5}/10}$, $u_{i+\sqrt{5}/10}$,
$u_{i+1/2}^{-}\\}$ from the given moment values
$\\{\overline{u},\overline{v}|_{i-1};\overline{u},\overline{v}|_{i};\overline{u},\overline{v}|_{i+1}\\}$
in detail as follows:

The 1D Reconstruction Algorithm:

Step 1. Reconstruct the Gauss-Lobatto point values of $u$.

Step 1.1. Reconstruct a series of polynomials of different degrees.

First, we select a series of central spatial stencils and reconstruct a zeroth
degree polynomial $q_{1}(x)$, a quadratic polynomial $q_{2}(x)$, a cubic
polynomial $q_{3}(x)$ and a quintic polynomial $q_{4}(x)$, respectively, such
that

$\displaystyle\left.\begin{aligned} \frac{1}{\triangle
x}\int_{I_{k}}q_{1}(x)dx&=\overline{u}_{k},k=i;\\\ \frac{1}{\triangle
x}\int_{I_{k}}q_{2}(x)dx&=\overline{u}_{k},k=i-1,i,i+1;\\\ \frac{1}{\triangle
x}\int_{I_{k}}q_{3}(x)dx&=\overline{u}_{k},k=i-1,i,i+1;\quad\frac{1}{\triangle
x}\int_{I_{k_{x}}}q_{3}(x)\frac{x-x_{k_{x}}}{\triangle
x}dx=\overline{v}_{k_{x}},{k_{x}}=i;\\\ \frac{1}{\triangle
x}\int_{I_{k}}q_{4}(x)dx&=\overline{u}_{k},k=i-1,i,i+1;\quad\frac{1}{\triangle
x}\int_{I_{k_{x}}}q_{4}(x)\frac{x-x_{k_{x}}}{\triangle
x}dx=\overline{v}_{k_{x}},{k_{x}}=i-1,i,i+1.\end{aligned}\right.$ (2.12)

Then, we obtain the following equivalent expressions for the above
reconstructed polynomials

$\displaystyle p_{l_{2}}(x)=\left\\{\begin{aligned}
&q_{1}(x),\,\,\,l_{2}=1,\\\
&\frac{1}{\gamma_{l_{2},l_{2}}}q_{l_{2}}(x)-\sum_{l=1}^{l_{2}-1}\frac{\gamma_{l,l_{2}}}{\gamma_{l_{2},l_{2}}}p_{l}(x),\,\,\,l_{2}=2,3,4,\end{aligned}\right.$
(2.13)

with $\sum_{l=1}^{l_{2}}\gamma_{l,l_{2}}=1,\gamma_{l_{2},l_{2}}\neq
0,l_{2}=2,3,4$, where these $\gamma_{l_{1},l_{2}}$ for
$l_{1}=1,...,l_{2};l_{2}=2,3,4$ are the linear weights and are defined as

$\gamma_{l_{1},l_{2}}=\frac{\overline{\gamma}_{l_{1},l_{2}}}{\sum_{l=1}^{l_{2}}\overline{\gamma}_{l,l_{2}}};\quad\overline{\gamma}_{l_{1},l_{2}}=10^{l_{1}-1};\quad
l_{1}=1,...,l_{2};\,\,l_{2}=2,3,4.$ (2.14)

Putting these linear weights into (2.13), we obtain the following relations

$\displaystyle\left.\begin{aligned} p_{1}(x)&=q_{1}(x),\\\
p_{2}(x)&=\frac{11}{10}q_{2}(x)-\frac{1}{10}q_{1}(x),\\\
p_{3}(x)&=\frac{111}{100}q_{3}(x)-\frac{11}{100}q_{2}(x),\\\
p_{4}(x)&=\frac{1111}{1000}q_{4}(x)-\frac{111}{1000}q_{3}(x).\end{aligned}\right.$
(2.15)

Step 1.2. Compute the corresponding nonlinear weights of the above
polynomials.

First, we compute the smoothness indicator $\beta_{l_{2}}$ of function
$p_{l_{2}}(x)$ in the interval $I_{i}$:

$\beta_{l_{2}}=\sum_{\alpha=1}^{\kappa}\int_{I_{i}}{\triangle
x}^{2\alpha-1}\left(\frac{d^{\alpha}p_{l_{2}}(x)}{dx^{\alpha}}\right)^{2}dx,\quad
l_{2}=2,3,4,$ (2.16)

where $\kappa=2,3,5$ for $l_{2}=2,3,4$. It is important to note that the
definition of $\beta_{1}$ is different, where a new polynomial
$q_{1}^{new}(x)$ is required and is defined as follows:
(1) We reconstruct two polynomials $q_{1L}(x)$ and $q_{1R}(x)$, such that

$\displaystyle\left.\begin{aligned} \frac{1}{\triangle
x}\int_{I_{k}}q_{1L}(x)dx&=\overline{u}_{k},\,\,\,k=i-1,i;\\\
\frac{1}{\triangle
x}\int_{I_{k}}q_{1R}(x)dx&=\overline{u}_{k},\,\,\,k=i,i+1,\end{aligned}\right.$
(2.17)

and then, we obtain their associated smoothness indicators

$\beta_{1L}=(\overline{u}_{i}-\overline{u}_{i-1})^{2},\quad\beta_{1R}=(\overline{u}_{i+1}-\overline{u}_{i})^{2},$
(2.18)

and the absolute difference $\tau_{1}$ between $\beta_{1L}$ and $\beta_{1R}$

$\tau_{1}=|\beta_{1R}-\beta_{1L}|^{2},$ (2.19)

where the selection of the power is to be consistent with the definition
(2.25) of $\tau_{4}$ later.
(2) We give these two polynomials $q_{1L}(x)$ and $q_{1R}(x)$ the same linear
weights $\gamma_{1L}=\gamma_{1R}=\frac{1}{2}$ and calculate the corresponding
nonlinear weights as

$\omega_{1L}=\frac{\overline{\omega}_{1L}}{\overline{\omega}_{1L}+\overline{\omega}_{1R}},\quad\omega_{1R}=\frac{\overline{\omega}_{1R}}{\overline{\omega}_{1L}+\overline{\omega}_{1R}},$
(2.20)
$\overline{\omega}_{1L}=\gamma_{1L}\left(1+\frac{\tau_{1}}{\beta_{1L}+\varepsilon}\right),\quad\overline{\omega}_{1R}=\gamma_{1R}\left(1+\frac{\tau_{1}}{\beta_{1R}+\varepsilon}\right),$
(2.21)

where $\varepsilon=10^{-10}$ is applied to avoid the denominator of (2.21) to
be zero.
(3) We obtain a new polynomial

$q_{1}^{new}(x)=\omega_{1L}q_{1L}(x)+\omega_{1R}q_{1R}(x),$ (2.22)

and set $\beta_{1}$ to be

$\beta_{1}=\int_{I_{i}}{\triangle
x}\left(\frac{dq_{1}^{new}(x)}{dx}\right)^{2}dx=\Big{(}\omega_{1L}\big{(}\overline{u}_{i}-\overline{u}_{i-1}\big{)}+\omega_{1R}\big{(}\overline{u}_{i+1}-\overline{u}_{i}\big{)}\Big{)}^{2}.$
(2.23)

Then, we still adopt the idea of WENO-Z as shown in [6] to define the
corresponding nonlinear weights

$\omega_{l_{1},4}=\frac{\overline{\omega}_{l_{1},4}}{\sum_{l=1}^{4}\overline{\omega}_{l,4}},\quad\overline{\omega}_{l_{1},4}=\gamma_{l_{1},4}\left(1+\frac{\tau_{4}}{\beta_{l_{1}}+\varepsilon}\right),\quad
l_{1}=1,...,4,$ (2.24)

where $\varepsilon$ is also taken to be $10^{-10}$ and the quantity $\tau_{4}$
is defined to be the absolute difference among above smoothness indicators

$\tau_{4}=\left(\frac{\sum_{l=1}^{3}|\beta_{4}-\beta_{l}|}{3}\right)^{2}.$
(2.25)

Step 1.3. Obtain an approximation polynomial $u_{i}(x)$ of $u(x)$.

The new reconstruction polynomial $u_{i}(x)$ of $u(x)$ is defined as

$u_{i}(x)=\sum_{l=1}^{4}\omega_{l,4}p_{l}(x),$ (2.26)

and the Gauss-Lobatto point values that we need are taken to be

$u_{i-1/2}^{+}=u_{i}(x_{i-1/2}),u_{i\mp\sqrt{5}/10}=u_{i}(x_{i\mp\sqrt{5}/10}),u_{i+1/2}^{-}=u_{i}(x_{i+1/2}).$
(2.27)

Step 2. Update the Gauss-Lobatto point values of $u$ in the troubled-cells.

Step 2.1. Identify the troubled-cells.

The so-called troubled-cells are those cells that may contain discontinuities.
In 2005, Qiu and Shu systematically investigated and compared a few troubled-
cell indicators for the Runge-Kutta discontinuous Galerkin method in [20].
Here, we choose the KXRCF troubled-cell indicator proposed by Krivodonova et
al. in [13] to identify the troubled-cells, and its judgment criterion is that
the target cell $I_{i}$ is identified to be a troubled-cell, if

$\aleph_{i}=\frac{\big{|}\int_{\partial
I_{i}^{-}}\big{(}u_{i}(x)-u_{n_{i}}(x)\big{)}ds\big{|}}{h_{i}^{\frac{l+1}{3}}\big{|}\partial
I_{i}^{-}\big{|}||u_{i}(x)||}>1,$ (2.28)

where $\partial I_{i}^{-}$ is the inflow boundary
($\overrightarrow{v}\cdot\overrightarrow{n}<0$, $\overrightarrow{v}$ is the
velocity of the flow and $\overrightarrow{n}$ is the outer normal vector to
$\partial I_{i}$), $I_{n_{i}}$ is the neighbor of $I_{i}$ on the side of
$\partial I_{i}^{-}$, $h_{i}$ is the length of the cell $I_{i}$, the parameter
$l$ (i.e. the degree of $u_{i}(x)$) is taken to be 5, $u_{i}(x)$ is the
approximation polynomial of $u(x)$ obtained in Step 1 above and the norm is
taken to be the $L^{\infty}$ norm.

Step 2.2. Modify the first order moments in the troubled-cells.

If the target cell $I_{i}$ is identified to be a troubled-cell, we would like
to modify the first order moment $\overline{v}_{i}$ in it. First, we
reconstruct a quartic polynomial, which satisfies

$\frac{1}{\triangle x}\int_{I_{k}}p_{0}(x)dx=\overline{u}_{k},\quad
k=i-1,i,i+1;\quad\frac{1}{\triangle
x}\int_{I_{k_{x}}}p_{0}(x)\frac{x-x_{k_{x}}}{\triangle
x}dx=\overline{v}_{k_{x}},\quad k_{x}=i-1,i+1,$ (2.29)

and then, we modify the first order moment $\overline{v}_{i}$ as

$\overline{v}_{i}=\frac{1}{\triangle
x}\int_{I_{i}}p_{0}(x)\frac{x-x_{i}}{\triangle
x}dx=-\frac{5}{76}\overline{u}_{i-1}+\frac{5}{76}\overline{u}_{i+1}-\frac{11}{38}\overline{v}_{i-1}-\frac{11}{38}\overline{v}_{i+1}.$
(2.30)

Step 2.3. Update the Gauss-Lobatto point values of $u$ for these troubled-
cells.

After modifying the first order moments in the troubled-cells, we repeat the
reconstruction process Step 1 for these troubled-cells to update corresponding
Gauss-Lobatto point values of $u$.

Step 3. Discretize the semi-discrete scheme in time.

After all these Gauss-Lobatto point values are obtained, we substitute them
into the formula of the numerical flux. Then, we discretize (2.4) by a third-
order TVD Runge-Kutta method in time

$\displaystyle\left\\{\begin{aligned} &U^{(1)}=U^{n}+\Delta tL(U^{n}),\\\
&U^{(2)}=\frac{3}{4}U^{n}+\frac{1}{4}U^{(1)}+\frac{1}{4}\Delta tL(U^{(1)}),\\\
&U^{n+1}=\frac{1}{3}U^{n}+\frac{2}{3}U^{(2)}+\frac{2}{3}\Delta
tL(U^{(2)}),\end{aligned}\right.$ (2.31)

to obtain a fully discrete scheme.

Remark 1. In Step 1 above, through a series of Taylor expansion analyses, we
can verify that $\beta_{l}=u^{{}^{\prime}2}\triangle x^{2}+O(\triangle
x^{4}),l=1,2,3,4$, thus $\beta_{4}-\beta_{l}=O(\triangle x^{4})$ for $l<4$ and
$\tau_{4}=O(\triangle x^{8})$, then

$\begin{split}u(x)-u_{i}(x)&=u(x)-\sum_{l=1}^{4}\omega_{l,4}p_{l}(x)\\\
&=\left[u(x)-\sum_{l=1}^{4}\gamma_{l,4}p_{l}(x)\right]+\sum_{l=1}^{4}(\omega_{l,4}-\gamma_{l,4})\Big{(}u\big{(}x\big{)}-p_{l}\big{(}x\big{)}\Big{)}\\\
&=O(\triangle x^{6})+O(\triangle x^{6})*O(\triangle x)\\\ &=O(\triangle
x^{6}).\end{split}$ (2.32)

According to above Taylor expansion analyses, we can prove that our
reconstruction algorithm can obtain sixth order of accuracy in smooth regions.

When there is a discontinuity in the target cell $I_{i}$, its first order
moment will become $\overline{v}_{i}=O(1)$, thus
$\beta_{1}=\beta_{2}=O(\triangle x^{2}),\beta_{3}=\beta_{4}=O(1)$ and
$\tau_{4}=O(1)$, then

$\displaystyle\left.\begin{aligned}
&\overline{\omega}_{1,4}=\gamma_{1,4}\big{(}1+O(\triangle
x^{-2})\big{)},\quad\overline{\omega}_{2,4}=\gamma_{2,4}\big{(}1+O(\triangle
x^{-2})\big{)};\\\
&\overline{\omega}_{3,4}=\gamma_{3,4}\big{(}1+O(1)\big{)},\quad\overline{\omega}_{4,4}=\gamma_{4,4}\big{(}1+O(1)\big{)}.\end{aligned}\right.$
(2.33)

That is to say, the proportion of the first two layers will be magnified
slightly, which will cause the discontinuity to be smoothed out to a certain
extent. But if we modify the first order moment of the troubled-cell $I_{i}$,
its first order moment will become $\overline{v}_{i}=O(\triangle x)$, thus
$\beta_{l}=u^{{}^{\prime}2}\triangle x^{2}+O(\triangle x^{4}),l=1,2,3,4$, and
$\tau_{4}=O(\triangle x^{8})$, then
$\overline{\omega}_{l,4}=\gamma_{l,4}\big{(}1+O(\triangle
x^{6})\big{)},l=1,2,3,4$, i.e. the nonlinear weight of each layer is close to
their corresponding linear weight, thus the resolution near the
discontinuities can be increased significantly.

As in our previous paper, the choice of these linear weights is also not
unique. Even though different choices of the linear weights do not affect the
order of accuracy in smooth regions, they do affect the resolution near
discontinuities. That is to say, the bigger the linear weights are for higher-
degree polynomials, the steeper the shock transitions become, but also the
more unstable the scheme becomes. This modification procedure can also improve
the stability of the scheme by reducing the magnitude of the first order
moments near the discontinuities.

Remark 2. In Step 2 above, we can also choose other indicators to identify
troubled-cells, such as the minmod-based total variation bounded (TVB) limiter
in [7], moment limiter of Biswas, Devine and Flaherty in [4], a modification
of the moment limiter by Burbeau, Sagaut and Bruneau in [5], the monotonicity-
preserving(MP) limiter in [24], a modification of the MP limiter in [21], a
troubled-cell indicator based on Harten’s subcell resolution idea in [9] and
so on. But as shown in Section 3, the KXRCF troubled-cell indicator works
pretty well for our scheme in the one-dimensional case. What needs a special
attention is that, for the one dimensional scalar equation, the solution $u$
is defined as our indicator variable, and then the corresponding
$\overrightarrow{v}=f^{{}^{\prime}}(u)$; for the one dimensional Euler system,
only the density $\rho$ is set to be our indicator variable, and then the
corresponding $\overrightarrow{v}=\mu$ is the velocity of the fluid. In short,
for the one dimensional case, the line integral average in the formula (2.28)
is actually the boundary point value, that is

$\frac{1}{\big{|}\partial I_{i}^{-}\big{|}}\left|\int_{\partial
I_{i}^{-}}\big{(}u_{i}(x)-u_{n_{i}}(x)\big{)}ds\right|=\left|(u_{i-\frac{1}{2}}^{+}-u_{i-\frac{1}{2}}^{-})*sf(\overrightarrow{v}_{i-\frac{1}{2}})+(u_{i+\frac{1}{2}}^{-}-u_{i+\frac{1}{2}}^{+})*sf(-\overrightarrow{v}_{i+\frac{1}{2}})\right|,$
(2.34)

and the norm $||u_{i}(x)||$ is taken to be the maximum norm of all the Gauss-
Lobatto point values in the cell $I_{i}$ (i.e.
$||u_{i}(x)||\approx\max\\{|u_{i-1/2}^{+}|,|u_{i-\sqrt{5}/10}|,|u_{i+\sqrt{5}/10}|,|u_{i+1/2}^{-}|\\}$).
Note that all the values used in the troubled-cell indicator are already
obtained in the reconstruction process Step 1, thus there is no need to
reconstruct an extra polynomial as in [25, 26].

According to [13],

$\displaystyle\int_{\partial I_{i}^{-}}\big{(}u_{i}-u_{n_{i}}\big{)}ds$
$\displaystyle=\int_{\partial I_{i}^{-}}\big{(}u_{i}-u\big{)}ds+\int_{\partial
I_{n_{i}}^{+}}\big{(}u-u_{n_{i}}\big{)}ds$ (2.35)
$\displaystyle=\left\\{\begin{aligned}
&O(h^{p+2})+O(h^{2p+2})=O(h^{p+2}),u|_{\partial I_{i}}\ is\ smooth,\\\
&O(h)+O(h)=O(h),u|_{\partial I_{i}}\ is\ discontinuous.\end{aligned}\right.$

Here, we let the troubled-cell indicator converge to the smooth case twice as
fast as the discontinuous case, then as either $h\rightarrow 0$ or
$p\rightarrow\infty$

$\displaystyle\aleph_{i}=\left\\{\begin{aligned}
&O(h^{\frac{2(l+1)}{3}})\stackrel{{\scriptstyle}}{{\longrightarrow}}0,u|_{\partial
I_{i}}\ is\ smooth,\\\
&\frac{1}{O(h^{\frac{l+1}{3}})}\stackrel{{\scriptstyle}}{{\longrightarrow}}\infty,u|_{\partial
I_{i}}\ is\ discontinuous,\end{aligned}\right.$ (2.36)

thus the troubled-cell indicator can be defined as

$\displaystyle\left\\{\begin{aligned} &\aleph_{i}<1,\,\,\,I_{i}\mbox{ is not a
troubled-cell},\\\ &\aleph_{i}>1,\,\,\,I_{i}\mbox{ is a troubled-
cell}.\end{aligned}\right.$ (2.37)

###  2.2 Two dimensional case

In this subsection, we then consider the following more complicated two-
dimensional hyperbolic conservation laws

$\displaystyle\left\\{\begin{aligned} &u_{t}+f(u)_{x}+g(u)_{y}=0,\\\
&u(x,y,0)=u_{0}(x,y).\end{aligned}\right.$ (2.38)

Still for the sake of simplicity, the computational domain is divided by a
uniform mesh $I_{i,j}=$ $[x_{i-1/2},x_{i+1/2}]\times[y_{j-1/2},y_{j+1/2}]$
with the uniform mesh sizes $\triangle x=x_{i+1/2}-x_{i-1/2}$ in the
$x$-direction and $\triangle y=y_{j+1/2}-y_{j-1/2}$ in the $y$-direction, and
the corresponding cell center is denoted by
$(x_{i},y_{j})=\Big{(}\frac{1}{2}\big{(}x_{i-1/2}+x_{i+1/2}\big{)},\frac{1}{2}\big{(}y_{j-1/2}+y_{j+1/2}\big{)}\Big{)}$.

Firstly, we multiply the governing equation (2.38) by $\frac{1}{\triangle
x\triangle y}$, $\frac{1}{\triangle x\triangle y}\frac{x-x_{i}}{\triangle x}$
and $\frac{1}{\triangle x\triangle y}\frac{y-y_{j}}{\triangle y}$,
respectively. Then, we integrate the resulting equations over the target cell
$I_{i,j}$ and perform the integration by parts, obtaining the following
equations

$\displaystyle\left\\{\begin{aligned} &\frac{1}{\triangle x\triangle
y}\int\limits_{I_{i,j}}u_{t}dxdy=-\frac{1}{\triangle x\triangle
y}\int\limits_{y_{j-1/2}}^{y_{j+1/2}}\Big{[}f\big{(}u(x_{i+1/2},y,t)\big{)}-f\big{(}u(x_{i-1/2},y,t)\big{)}\Big{]}dy\\\
&-\frac{1}{\triangle x\triangle
y}\int\limits_{x_{i-1/2}}^{x_{i+1/2}}\Big{[}g\big{(}u(x,y_{j+1/2},t)\big{)}-g\big{(}u(x,y_{j-1/2},t)\big{)}\Big{]}dx,\\\
&\frac{1}{\triangle x\triangle
y}\int\limits_{I_{i,j}}u_{t}\frac{x-x_{i}}{\triangle
x}dxdy=-\frac{1}{2\triangle x\triangle
y}\int\limits_{y_{j-1/2}}^{y_{j+1/2}}\Big{[}f\big{(}u(x_{i+1/2},y,t)\big{)}+f\big{(}u(x_{i-1/2},y,t)\big{)}\Big{]}dy\\\
&+\frac{1}{(\triangle x)^{2}\triangle
y}\int\limits_{I_{i,j}}f(u)dxdy-\frac{1}{\triangle x\triangle
y}\int\limits_{x_{i-1/2}}^{x_{i+1/2}}\Big{[}g\big{(}u(x,y_{j+1/2},t)\big{)}-g\big{(}u(x,y_{j-1/2},t)\big{)}\Big{]}\frac{x-x_{i}}{\triangle
x}dx,\\\ &\frac{1}{\triangle x\triangle
y}\int\limits_{I_{i,j}}u_{t}\frac{y-y_{j}}{\triangle
y}dxdy=-\frac{1}{\triangle x\triangle
y}\int\limits_{y_{j-1/2}}^{y_{j+1/2}}\Big{[}f\big{(}u(x_{i+1/2},y,t)\big{)}-f\big{(}u(x_{i-1/2},y,t)\big{)}\Big{]}\frac{y-y_{j}}{\triangle
y}dy\\\ &-\frac{1}{2\triangle x\triangle
y}\int\limits_{x_{i-1/2}}^{x_{i+1/2}}\Big{[}g\big{(}u(x,y_{j+1/2},t)\big{)}+g\big{(}u(x,y_{j-1/2},t)\big{)}\Big{]}dx+\frac{1}{\triangle
x(\triangle y)^{2}}\int\limits_{I_{i,j}}g(u)dxdy.\end{aligned}\right.$ (2.39)

Next, we define the zeroth and first order moments in the $x$ and $y$
directions as follows

$\displaystyle\left\\{\begin{aligned}
&\tilde{\overline{u}}_{i,j}(t)=\frac{1}{\triangle x\triangle
y}\int_{I_{i,j}}u(x,y,t)dxdy,\\\
&\tilde{\overline{v}}_{i,j}(t)=\frac{1}{\triangle x\triangle
y}\int_{I_{i,j}}u(x,y,t)\frac{x-x_{i}}{\triangle x}dxdy,\\\
&\tilde{\overline{w}}_{i,j}(t)=\frac{1}{\triangle x\triangle
y}\int_{I_{i,j}}u(x,y,t)\frac{y-y_{j}}{\triangle y}dxdy,\end{aligned}\right.$
(2.40)

and then, we swap the spatial integration and time derivation and approximate
the flux by a numerical flux, obtaining the following semi-discrete
conservative scheme

$\displaystyle\left\\{\begin{aligned}
&\frac{d\tilde{\overline{u}}_{i,j}(t)}{dt}=-\frac{1}{\triangle x\triangle
y}\int_{y_{j-1/2}}^{y_{j+1/2}}\Big{[}\hat{f}\big{(}u(x_{i+1/2},y)\big{)}-\hat{f}\big{(}u(x_{i-1/2},y)\big{)}\Big{]}dy\\\
&-\frac{1}{\triangle x\triangle
y}\int_{x_{i-1/2}}^{x_{i+1/2}}\Big{[}\hat{g}\big{(}u(x,y_{j+1/2})\big{)}-\hat{g}\big{(}u(x,y_{j-1/2})\big{)}\Big{]}dx,\\\
&\frac{d\tilde{\overline{v}}_{i,j}(t)}{dt}=-\frac{1}{2\triangle x\triangle
y}\int_{y_{j-1/2}}^{y_{j+1/2}}\Big{[}\hat{f}\big{(}u(x_{i+1/2},y)\big{)}+\hat{f}\big{(}u(x_{i-1/2},y)\big{)}\Big{]}dy\\\
&+\frac{1}{\triangle x}F_{i,j}(u)-\frac{1}{\triangle x\triangle
y}\int_{x_{i-1/2}}^{x_{i+1/2}}\Big{[}\hat{g}\big{(}u(x,y_{j+1/2})\big{)}-\hat{g}\big{(}u(x,y_{j-1/2})\big{)}\Big{]}\frac{x-x_{i}}{\triangle
x}dx,\\\ &\frac{d\tilde{\overline{w}}_{i,j}(t)}{dt}=-\frac{1}{\triangle
x\triangle
y}\int_{y_{j-1/2}}^{y_{j+1/2}}\Big{[}\hat{f}\big{(}u(x_{i+1/2},y)\big{)}-\hat{f}\big{(}u(x_{i-1/2},y)\big{)}\Big{]}\frac{y-y_{j}}{\triangle
y}dy\\\ &-\frac{1}{2\triangle x\triangle
y}\int_{x_{i-1/2}}^{x_{i+1/2}}\Big{[}\hat{g}\big{(}u(x,y_{j+1/2})\big{)}+\hat{g}\big{(}u(x,y_{j-1/2})\big{)}\Big{]}dx+\frac{1}{\triangle
y}G_{i,j}(u).\end{aligned}\right.$ (2.41)

Here the numerical fluxes are still taken to be the HLLC fluxes, and the
integral terms are also approximated by a four-point Gauss-Lobatto
integration, for instance

$F_{i,j}(u)=\frac{1}{\triangle x\triangle
y}\int_{I_{i,j}}f(u)dxdy\approx\sum_{k=1}^{4}\sum_{l=1}^{4}\omega_{k}\omega_{l}f\big{(}u(x_{k}^{GL},y_{l}^{GL})\big{)},$
(2.42) $\frac{1}{\triangle
y}\int_{y_{j-1/2}}^{y_{j+1/2}}\hat{f}\big{(}u(x_{i+1/2},y)\big{)}dy\approx\sum_{l=1}^{4}\omega_{l}\hat{f}\big{(}u(x_{i+1/2},y_{l}^{GL})\big{)},$
(2.43) $\frac{1}{\triangle
x}\int_{x_{i-1/2}}^{x_{i+1/2}}\frac{x-x_{i}}{\triangle
x}\hat{g}\big{(}u(x,y_{j+1/2})\big{)}dx\approx\sum_{k=1}^{4}\omega_{k}\frac{x_{k}^{GL}-x_{i}}{\triangle
x}\hat{g}\big{(}u(x_{k}^{GL},y_{j+1/2})\big{)},$ (2.44)

where the weights are $\omega_{1}=\omega_{4}=\frac{1}{12}$,
$\omega_{2}=\omega_{3}=\frac{5}{12}$ and the quadrature points on the target
cell $I_{i,j}$ are

$\displaystyle
x_{1}^{GL}=x_{i-1/2},x_{2}^{GL}=x_{i-\sqrt{5}/10},x_{3}^{GL}=x_{i+\sqrt{5}/10},x_{4}^{GL}=x_{i+1/2};$
$\displaystyle
y_{1}^{GL}=y_{j-1/2},y_{2}^{GL}=y_{j-\sqrt{5}/10},y_{3}^{GL}=y_{j+\sqrt{5}/10},y_{4}^{GL}=y_{j+1/2},$

where $x_{i+\xi}=x_{i}+\xi\triangle x$ and $y_{j+\eta}=y_{j}+\eta\triangle y$.

Finally, we present the spatial reconstruction procedure of the point values
$\\{u_{i\mp 1/2,j+\eta_{l}}^{\pm}|l=1,2,3,4\\}$, $\\{u_{i+\xi_{k},j\mp
1/2}^{\pm}|k=1,2,3,4\\}$ and $\\{u_{i+\xi_{k},j+\eta_{l}}|k=2,3;l=2,3\\}$
($\xi_{1}=\eta_{1}=-1/2$, $\xi_{2}=\eta_{2}=-\sqrt{5}/10$,
$\xi_{3}=\eta_{3}=\sqrt{5}/10$ and $\xi_{4}=\eta_{4}=1/2$) from the given
cell-average values
$\\{\tilde{\overline{u}},\tilde{\overline{v}},\tilde{\overline{w}}|_{k_{x},k_{y}};k_{x}=i-1,i,i+1;k_{y}=j-1,j,j+1\\}$
in detail as follows:

The 2D Reconstruction Algorithm:

Step 1. Reconstruct the Gauss-Lobatto point values of $u$.

7 | 8 | 9 | $j+1$
---|---|---|---
4 | 5 | 6 | $j$
1 | 2 | 3 | $j-1$
$i-1$ | $i$ | $i+1$ |

The big stencil and its new labels.

Step 1.1. Reconstruct a series of polynomials of different degrees.

Step 1.1.1. Reconstruct a zeroth degree polynomial $q_{1}(x,y)$ such that

$\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{1}(x,y)dxdy=\tilde{\overline{u}}_{k},\quad k=5.$ (2.45)

Step 1.1.2. Reconstruct a quadratic polynomial $q_{2}(x,y)$ such that

$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{2}(x,y)dxdy=\tilde{\overline{u}}_{k},\quad k=1,...,9.$ (2.46)

Step 1.1.3. Reconstruct a cubic polynomial $q_{3}(x,y)$ such that

$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{3}(x,y)dxdy=\tilde{\overline{u}}_{k},\quad k=1,...,9;$
$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k_{x}}}q_{3}(x,y)\frac{x-x_{k_{x}}}{\triangle
x}dxdy=\tilde{\overline{v}}_{k_{x}},\quad{k_{x}}=5;$ (2.47)
$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k_{y}}}q_{3}(x,y)\frac{y-y_{k_{y}}}{\triangle
y}dxdy=\tilde{\overline{w}}_{k_{y}},\quad{k_{y}}=5.$

Step 1.1.4. Reconstruct a quintic polynomial $q_{4}(x,y)$ such that

$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{4}(x,y)dxdy=\tilde{\overline{u}}_{k},\quad k=1,...,9;$
$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k_{x}}}q_{4}(x,y)\frac{x-x_{k_{x}}}{\triangle
x}dxdy=\tilde{\overline{v}}_{k_{x}},\quad{k_{x}}=1,3,4,5,6,7,9;$ (2.48)
$\displaystyle\frac{1}{\triangle x\triangle
y}\int_{I_{k_{y}}}q_{4}(x,y)\frac{y-y_{k_{y}}}{\triangle
y}dxdy=\tilde{\overline{w}}_{k_{y}},\quad{k_{y}}=1,2,3,5,7,8,9.$

What needs special attention is that the number of equations is greater than
the number of unknowns when we reconstruct the quadratic polynomial
$q_{2}(x,y)$, the cubic polynomial $q_{3}(x,y)$ and the quintic polynomial
$q_{4}(x,y)$. To solve this problem, we require these polynomials must have
the same cell average as $u$ on the target cell $I_{i,j}$ (to maintain
conservation) and match the other conditions in a least square sense as
described in [11].

Step 1.1.5. Then, we further manipulate the above reconstructed polynomials to
obtain the following equivalent expressions

$\displaystyle p_{l_{2}}(x,y)=\left\\{\begin{aligned}
&q_{1}(x,y),\,\,\,l_{2}=1,\\\
&\frac{1}{\gamma_{l_{2},l_{2}}}q_{l_{2}}(x,y)-\sum_{l=1}^{l_{2}-1}\frac{\gamma_{l,l_{2}}}{\gamma_{l_{2},l_{2}}}p_{l}(x,y),\,\,\,l_{2}=2,3,4,\end{aligned}\right.$
(2.49)

with $\sum_{l=1}^{l_{2}}\gamma_{l,l_{2}}=1,\gamma_{l_{2},l_{2}}\neq
0,l_{2}=2,3,4$, where these $\gamma_{l_{1},l_{2}}$ for
$l_{1}=1,...,l_{2};l_{2}=2,3,4$ are still the linear weights and are defined
as (2.14). Likewise putting these linear weights into (2.49), we obtain the
following relations

$\displaystyle\left.\begin{aligned} p_{1}(x,y)&=q_{1}(x,y),\\\
p_{2}(x,y)&=\frac{11}{10}q_{2}(x,y)-\frac{1}{10}q_{1}(x,y),\\\
p_{3}(x,y)&=\frac{111}{100}q_{3}(x,y)-\frac{11}{100}q_{2}(x,y),\\\
p_{4}(x,y)&=\frac{1111}{1000}q_{4}(x,y)-\frac{111}{1000}q_{3}(x,y).\end{aligned}\right.$
(2.50)

Step 1.2. Compute the corresponding nonlinear weights of the above
polynomials.

First, we compute the smoothness indicator $\beta_{l_{2}}$ of the function
$p_{l_{2}}(x,y)$ in the interval $I_{i,j}$:

$\beta_{l_{2}}=\sum_{|\alpha|=1}^{\kappa}\int_{I_{i,j}}|I_{i,j}|^{|\alpha|-1}\left(\frac{\partial^{|\alpha|}}{\partial
x^{{\alpha}_{x}}\partial y^{{\alpha}_{y}}}p_{l_{2}}(x,y)\right)^{2}dxdy,\quad
l_{2}=2,3,4,$ (2.51)

where $\alpha=(\alpha_{x},\alpha_{y})$, $|\alpha|=\alpha_{x}+\alpha_{y}$ and
$\kappa=2,3,5$ for $l_{2}=2,3,4$. Note that the definition of $\beta_{1}$ is
an exception, where a new polynomial $q_{1}^{new}(x,y)$ is required and is
defined as follows:
(1) We reconstruct four polynomials $q_{1k}(x,y)$ for $k=1,2,3,4$, such that

$\displaystyle\left.\begin{aligned} \frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{11}(x,y)dxdy&=\tilde{\overline{u}}_{k},\quad k=4,5,8;\\\
\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{12}(x,y)dxdy&=\tilde{\overline{u}}_{k},\quad k=5,6,8;\\\
\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{13}(x,y)dxdy&=\tilde{\overline{u}}_{k},\quad k=2,5,6;\\\
\frac{1}{\triangle x\triangle
y}\int_{I_{k}}q_{14}(x,y)dxdy&=\tilde{\overline{u}}_{k},\quad
k=2,4,5,\end{aligned}\right.$ (2.52)

and then, we obtain their associated smoothness indicators

$\displaystyle\left.\begin{aligned}
\beta_{11}&=(\tilde{\overline{u}}_{i,j}-\tilde{\overline{u}}_{i-1,j})^{2}+(\tilde{\overline{u}}_{i,j+1}-\tilde{\overline{u}}_{i,j})^{2},\\\
\beta_{12}&=(\tilde{\overline{u}}_{i+1,j}-\tilde{\overline{u}}_{i,j})^{2}+(\tilde{\overline{u}}_{i,j+1}-\tilde{\overline{u}}_{i,j})^{2},\\\
\beta_{13}&=(\tilde{\overline{u}}_{i+1,j}-\tilde{\overline{u}}_{i,j})^{2}+(\tilde{\overline{u}}_{i,j}-\tilde{\overline{u}}_{i,j-1})^{2},\\\
\beta_{14}&=(\tilde{\overline{u}}_{i,j}-\tilde{\overline{u}}_{i-1,j})^{2}+(\tilde{\overline{u}}_{i,j}-\tilde{\overline{u}}_{i,j-1})^{2},\end{aligned}\right.$
(2.53)

and the absolute difference $\tau_{1}$ among these smoothness indicators

$\tau_{1}=\left(\frac{\sum_{k\neq l}|\beta_{1k}-\beta_{1l}|}{6}\right)^{2},$
(2.54)

where the selection of the power is to be consistent with the definition
(2.19) of $\tau_{1}$ in the one dimension.
(2) We give these four polynomials $q_{1k}(x,y)$ the same linear weights
$\gamma_{1k}=\frac{1}{4}$ for $k=1,2,3,4$ and calculate the corresponding
nonlinear weights as

$\omega_{1k}=\frac{\overline{\omega}_{1k}}{\sum_{l=1}^{4}\overline{\omega}_{1l}},$
(2.55)
$\overline{\omega}_{1k}=\gamma_{1k}\left(1+\frac{\tau_{1}}{\beta_{1k}+\varepsilon}\right),\quad
k=1,2,3,4,$ (2.56)

where $\varepsilon$ is still taken to be $10^{-10}$ as in the one dimensional
case.
(3) We obtain a new polynomial

$q_{1}^{new}(x,y)=\sum_{l=1}^{4}\omega_{1l}q_{1l}(x,y),$ (2.57)

and set $\beta_{1}$ to be

$\beta_{1}=\sum_{|\alpha|=1}|I_{i,j}|^{|\alpha|}\left(\frac{\partial^{|\alpha|}}{\partial
x^{{\alpha}_{x}}\partial y^{{\alpha}_{y}}}q_{1}^{new}(x,y)\right)^{2},$ (2.58)

where $\alpha=(\alpha_{x},\alpha_{y})$, $|\alpha|=\alpha_{x}+\alpha_{y}$.

Then, we define the corresponding nonlinear weights as (2.24).

Step 1.3. Obtain an approximation polynomial $u_{i,j}(x,y)$ of $u(x,y)$.

The new reconstruction polynomial $u_{i,j}(x,y)$ of $u(x,y)$ is defined as

$u_{i,j}(x,y)=\sum_{l=1}^{4}\omega_{l,4}p_{l}(x,y),$ (2.59)

and the Gauss-Lobatto point values that we need are taken to be

$\displaystyle\left.\begin{aligned} u_{i\mp
1/2,j+\eta_{l}}^{\pm}&=u_{i,j}(x_{i\mp
1/2},y_{j+\eta_{l}}),\,\,\,l=1,2,3,4;\\\ u_{i+\xi_{k},j\mp
1/2}^{\pm}&=u_{i,j}(x_{i+\xi_{k}},y_{j\mp 1/2}),\,\,\,k=1,2,3,4;\\\
u_{i+\xi_{k},j+\eta_{l}}&=u_{i,j}(x_{i+\xi_{k}},y_{j+\eta_{l}}),\,\,\,k=2,3;\,\,l=2,3,\end{aligned}\right.$
(2.60)

where $\xi_{1}=\eta_{1}=-1/2$, $\xi_{2}=\eta_{2}=-\sqrt{5}/10$,
$\xi_{3}=\eta_{3}=\sqrt{5}/10$ and $\xi_{4}=\eta_{4}=1/2$.

Step 2. Update the Gauss-Lobatto point values of $u$ in the troubled-cells.

Here, we still choose the KXRCF troubled-cell indicator to identify the
troubled-cells as in the one dimensional case, that is the target cell
$I_{i,j}$ is identified to be a troubled-cell, if

$\aleph_{i,j}=\frac{\big{|}\int_{\partial
I_{i,j}^{-}}\big{(}u_{i,j}(x,y)-u_{n_{i,j}}(x,y)\big{)}ds\big{|}}{h_{i,j}^{\frac{l+1}{3}}\big{|}\partial
I_{i,j}^{-}\big{|}||u_{i,j}(x,y)||}>1,$ (2.61)

where $\partial I_{i,j}^{-}$ is the inflow boundary
($\overrightarrow{v}\cdot\overrightarrow{n}<0$, $\overrightarrow{v}$ is the
velocity of the flow and $\overrightarrow{n}$ is the outer normal vector to
$\partial I_{i,j}$), $I_{n_{i,j}}$ is the neighbor of $I_{i,j}$ on the side of
$\partial I_{i,j}^{-}$, $h_{i,j}$ is the length of the cell $I_{i,j}$, the
parameter $l$ (i.e. the degree of $u_{i,j}(x,y)$) is also taken to be 5,
$u_{i,j}(x,y)$ is the approximation polynomial of $u(x,y)$ obtained in Step 1
and the norm is still taken to be the $L^{\infty}$ norm. If the target cell
$I_{i,j}$ is identified to be a troubled-cell, we would like to modify the
first order moment $\tilde{\overline{v}}_{i,j}$ in the $x$ direction by using
the information of
$\\{\tilde{\overline{u}}_{i-1,j},\tilde{\overline{u}}_{i,j},\tilde{\overline{u}}_{i+1,j},\tilde{\overline{v}}_{i-1,j},\tilde{\overline{v}}_{i+1,j}\\}$
and modify the first order moment $\tilde{\overline{w}}_{i,j}$ in the $y$
direction by using the information of
$\\{\tilde{\overline{u}}_{i,j-1},\tilde{\overline{u}}_{i,j},\tilde{\overline{u}}_{i,j+1},\tilde{\overline{w}}_{i,j-1},\tilde{\overline{w}}_{i,j+1}\\}$
in a dimension-by-dimension manner. After modifying the first order moments in
the troubled-cells, we repeat the reconstruction process Step 1 for these
troubled-cells to update the corresponding Gauss-Lobatto point values of $u$.

Step 3. Discretize the semi-discrete scheme in time.

After all these Gauss-Lobatto point values are obtained, we put them into the
formula of the numerical flux. Then, we discretize (2.41) by the third-order
TVD Runge-Kutta method (2.31) in time to complete the entire discretization
process.

Remark 3. In Step 2 above, we still choose the KXRCF troubled-cell indicator
to identify troubled-cells. As shown in Section 3, the KXRCF troubled-cell
indicator works pretty well for our scheme in the two-dimensional case as
well. What needs a special attention is that, for the two dimensional scalar
equation, the solution $u$ is defined as our indicator variable, and then the
corresponding $\overrightarrow{v}=f^{{}^{\prime}}(u)$ in the $x$-direction and
$\overrightarrow{v}=g^{{}^{\prime}}(u)$ in the $y$-direction; for the two
dimensional Euler system, only the density $\rho$ is set to be our indicator
variable, and then the corresponding $\overrightarrow{v}=\mu$ is the velocity
in the $x$-direction of the fluid and $\overrightarrow{v}=\nu$ is the velocity
in the $y$-direction of the fluid. In short, for the two dimensional case, the
line integral average in the formula (2.61) is approximated by a four-point
Gauss-Lobatto integration, that is

$\displaystyle\frac{1}{\big{|}\partial I_{i,j}^{-}\big{|}}\left|\int_{\partial
I_{i,j}^{-}}\big{(}u_{i,j}(x,y)-u_{n_{i,j}}(x,y)\big{)}ds\right|$
$\displaystyle=\Bigg{|}\sum_{l=1}^{4}\triangle
y\omega_{l}(u_{i-\frac{1}{2},j+\eta_{l}}^{+}-u_{i-\frac{1}{2},j+\eta_{l}}^{-})*sf(\overrightarrow{v}_{i-\frac{1}{2},j})$
$\displaystyle+\sum_{l=1}^{4}\triangle
y\omega_{l}(u_{i+\frac{1}{2},j+\eta_{l}}^{-}-u_{i+\frac{1}{2},j+\eta_{l}}^{+})*sf(-\overrightarrow{v}_{i+\frac{1}{2},j})$
$\displaystyle+\sum_{k=1}^{4}\triangle
x\omega_{k}(u_{i+\xi_{k},j-\frac{1}{2}}^{+}-u_{i+\xi_{k},j-\frac{1}{2}}^{-})*sf(\overrightarrow{v}_{i,j-\frac{1}{2}})$
$\displaystyle+\sum_{k=1}^{4}\triangle
x\omega_{k}(u_{i+\xi_{k},j+\frac{1}{2}}^{-}-u_{i+\xi_{k},j+\frac{1}{2}}^{+})*sf(-\overrightarrow{v}_{i,j+\frac{1}{2}})\Bigg{|}\frac{1}{\big{|}\partial
I_{i,j}^{-}\big{|}},$ (2.62)

and the norm $||u_{i,j}(x,y)||$ is taken to be the maximum norm of all the
Gauss-Lobatto point values in the cell $I_{i,j}$(i.e.
$||u_{i,j}(x,y)||\approx\max$$\\{|u_{i\mp 1/2,j+\eta_{l}}^{\pm}|:\,l=1,2,3,4$;
$|u_{i+\xi_{k},j\mp 1/2}^{\pm}|:\,k=1,2,3,4$;
$|u_{i+\xi_{k},j+\eta_{l}}|:\,k=2,3;l=2,3\\}$). Note that all the values used
in the troubled-cell indicator are also already obtained in the reconstruction
process Step 1, thus there is no need to reconstruct an extra polynomial as in
[25, 26].

##  3 Numerical tests

In this section, a number of typical numerical examples are given to
demonstrate the stability and resolution of our moment-based multi-resolution
HWENO scheme which is termed as “HWENO6-M5-I/NI” where “HWENO6” means the
reconstruction process is of the sixth order of accuracy, “M5” means the first
order moments of the troubled-cells are modified by a quartic polynomial, and
“I” means only the first order moments of the troubled-cells are modified
while “NI” means the first order moments of all the cells are modified without
judgment. Before we start to show the results of the examples we have
calculated, let us first explain some of the parameters in particular: the
first one is that we set the CFL number as 0.6 for both the one and two
dimensional cases, but note that for the accuracy tests a suitably reduced
time step is used in order to ensure the dominance of the spatial error; the
second one is that we take the linear weights as $\overline{\gamma}_{1,4}=1$,
$\overline{\gamma}_{2,4}=10$, $\overline{\gamma}_{3,4}=100$ and
$\overline{\gamma}_{4,4}=1000$ both in the one and two dimensions in this
paper.

Example 3.1. One-dimensional scalar Burgers’ equation:

$\mu_{t}+\left(\frac{\mu^{2}}{2}\right)_{x}=0,\ \ 0<x<2,$ (3.1)

with the initial condition $\mu(x,0)=0.5+\sin(\pi x)$ and periodic boundary
condition. As we know, when the time reaches $T=0.5/\pi$, its solution is
still smooth, and now the corresponding errors and orders obtained by the
HWENO6-M5-I and HWENO6-M5-NI schemes are listed in the Table 3.1; but when the
time reaches $T=1.5/\pi$, there will appear discontinuities in its solution,
and now the corresponding numerical solutions obtained by the HWENO6 and
HWENO6-M5-I schemes are drawn on the Fig 3.1 in comparison with the reference
exact solution.

Table 3.1:  $\mu_{t}+\left(\frac{\mu^{2}}{2}\right)_{x}=0$. The initial condition $\mu(x,0)=0.5+\sin(\pi x)$. Periodic boundary condition. $T=0.5/\pi$. HWENO6-M5-I and HWENO6-M5-NI schemes. $L^{1}$ and $L^{\infty}$. | HWENO6-M5-I scheme | HWENO6-M5-NI scheme
---|---|---
grid points |  $L^{1}$ error | order |  $L^{\infty}$ error | order |  $L^{1}$ error | order |  $L^{\infty}$ error | order
10 | 5.24E-03 |  | 2.26E-02 |  | 6.48E-03 |  | 3.01E-02 |
20 | 5.67E-04 | 3.21 | 4.61E-03 | 2.29 | 3.11E-04 | 4.38 | 9.84E-04 | 4.93
40 | 1.23E-06 | 8.85 | 1.37E-05 | 8.40 | 6.97E-06 | 5.48 | 2.55E-05 | 5.27
80 | 2.90E-09 | 8.72 | 3.27E-08 | 8.71 | 5.96E-08 | 6.87 | 9.36E-07 | 4.77
160 | 2.46E-11 | 6.88 | 3.52E-10 | 6.54 | 1.88E-09 | 4.98 | 2.95E-08 | 4.99
320 | 3.95E-13 | 5.96 | 3.69E-12 | 6.57 | 6.02E-11 | 4.97 | 9.53E-10 | 4.95

Figure 3.1:  1D-Burgers’ equation. $T=1.5/\pi$. Left: density; right: density
zoomed in. Solid line: the exact solution; squares: the result of the HWENO6
scheme; triangles: the result of the HWENO6-M5-I scheme. Number of cells: 200.

Example 3.2. Two-dimensional scalar Burgers’ equation:

$\mu_{t}+\left(\frac{\mu^{2}}{2}\right)_{x}+\left(\frac{\mu^{2}}{2}\right)_{y}=0,\
\ 0<x,y<4,$ (3.2)

with the initial condition $\mu(x,y,0)=0.5+\sin\big{(}\pi(x+y)/2\big{)}$ and
periodic boundary condition. As is known to us, when the time reaches
$T=0.5/\pi$, its solution is still smooth, and now the corresponding errors
and orders obtained by the HWENO6-M5-I and HWENO6-M5-NI schemes are listed in
the Table 3.2; but when the time reaches $T=1.5/\pi$, there will appear
discontinuities in its solution, and now the corresponding numerical solutions
at $x=y$ obtained by the HWENO6 and HWENO6-M5-I schemes are drawn on the Fig
3.2 in comparison with the reference exact solution.

Table 3.2:  $\mu_{t}+\left(\frac{\mu^{2}}{2}\right)_{x}+\left(\frac{\mu^{2}}{2}\right)_{y}=0$. The initial condition $\mu(x,y,0)=0.5+\sin\big{(}\pi(x+y)/2\big{)}$. Periodic boundary condition. $T=0.5/\pi$. HWENO6-M5-I and HWENO6-M5-NI schemes. $L^{1}$ and $L^{\infty}$. | HWENO6-M5-I scheme | HWENO6-M5-NI scheme
---|---|---
grid points |  $L^{1}$ error | order |  $L^{\infty}$ error | order |  $L^{1}$ error | order |  $L^{\infty}$ error | order
10$\times$ 10 | 1.03E-02 |  | 2.96E-02 |  | 1.08E-02 |  | 4.31E-02 |
20$\times$ 20 | 7.78E-05 | 7.05 | 2.65E-04 | 6.80 | 1.42E-04 | 6.25 | 8.03E-04 | 5.74
40$\times$ 40 | 2.41E-06 | 5.01 | 1.56E-05 | 4.09 | 6.62E-06 | 4.42 | 3.84E-05 | 4.39
80$\times$ 80 | 3.86E-08 | 5.97 | 5.42E-07 | 4.84 | 1.03E-07 | 6.00 | 8.33E-07 | 5.53
160$\times$160 | 5.62E-10 | 6.10 | 9.83E-09 | 5.78 | 2.87E-09 | 5.17 | 2.30E-08 | 5.18
320$\times$320 | 6.83E-12 | 6.36 | 1.24E-10 | 6.31 | 9.77E-11 | 4.88 | 7.42E-10 | 4.95

Figure 3.2:  2D-Burgers’ equation. $T=1.5/\pi$. Left: density; right: density
zoomed in at $x=y$. Solid line: the exact solution; squares: the result of the
HWENO6 scheme; triangles: the result of the HWENO6-M5-I scheme. Number of
cells: 200$\times$200\.

Example 3.3. One-dimensional Euler equations:

$\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\\ \rho\mu\\\
E\end{array}\right)+\frac{\partial}{\partial
x}\left(\begin{array}[]{c}\rho\mu\\\ \rho\mu^{2}+p\\\
\mu(E+p)\end{array}\right)=0,\ \ 0<x<2\pi,$ (3.3)

where $\rho$ is the density, $\mu$ is the velocity, $E$ is the total energy
and $p$ is the pressure. The initial conditions are

$\begin{array}[]{l}\rho(x,0)=\frac{1+0.2\sin(x)}{2\sqrt{3}},\,\,\mu(x,0)=\sqrt{\gamma}\rho(x,0),\,\,p(x,0)=\rho(x,0)^{\gamma},\end{array}$
(3.4)

and the boundary conditions are periodic. The exact solution of above Euler
equations is given in [14]. When $T=3$, its solution is still smooth, and now
the corresponding errors and orders obtained by the HWENO6-M5-I and
HWENO6-M5-NI schemes are listed in the Table 3.3.

Table 3.3:  1D-Euler equations: The initial condition $\rho(x,0)=\frac{1+0.2\sin(x)}{2\sqrt{3}}$, $\mu(x,0)=\sqrt{\gamma}\rho(x,0)$ and $p(x,0)=\rho(x,0)^{\gamma}$. Periodic boundary condition. $T=3$. HWENO6-M5-I and HWENO6-M5-NI schemes. $L^{1}$ and $L^{\infty}$. | HWENO6-M5-I scheme | HWENO6-M5-NI scheme
---|---|---
grid points |  $L^{1}$ error | order |  $L^{\infty}$ error | order |  $L^{1}$ error | order |  $L^{\infty}$ error | order
10 | 8.54E-04 |  | 2.95E-03 |  | 8.76E-04 |  | 3.67E-03 |
20 | 1.87E-05 | 5.52 | 1.41E-04 | 4.39 | 1.09E-04 | 3.01 | 5.35E-04 | 2.78
40 | 4.99E-07 | 5.23 | 9.01E-06 | 3.97 | 8.44E-06 | 3.69 | 1.05E-04 | 2.35
80 | 7.57E-09 | 6.04 | 1.95E-07 | 5.53 | 2.49E-07 | 5.08 | 5.37E-06 | 4.28
160 | 1.10E-10 | 6.10 | 2.60E-09 | 6.23 | 7.12E-09 | 5.13 | 1.66E-07 | 5.02
320 | 1.69E-12 | 6.03 | 3.65E-11 | 6.15 | 2.15E-10 | 5.05 | 5.03E-09 | 5.05

Example 3.4. Two-dimensional Euler equations:

$\frac{\partial}{\partial t}\left(\begin{array}[]{c}\rho\\\ \rho\mu\\\
\rho\nu\\\ E\end{array}\right)+\frac{\partial}{\partial
x}\left(\begin{array}[]{c}\rho\mu\\\ \rho\mu^{2}+p\\\ \rho\mu\nu\\\
\mu(E+p)\end{array}\right)+\frac{\partial}{\partial
y}\left(\begin{array}[]{c}\rho\nu\\\ \rho\mu\nu\\\ \rho\nu^{2}+p\\\
\nu(E+p)\end{array}\right)=0,\ \ 0<x,y<4\pi,$ (3.5)

where $\rho$ is the density, $\mu$ is the velocity in the $x$-direction, $\nu$
is the velocity in the $y$-direction, $E$ is the total energy and $p$ is the
pressure. The initial conditions are

$\begin{array}[]{l}\rho(x,y,0)=\frac{1+0.2\sin(\frac{x+y}{2})}{\sqrt{6}},\,\,\mu(x,y,0)=\nu(x,y,0)=\sqrt{\frac{\gamma}{2}}\rho(x,y,0),\,\,p(x,y,0)=\rho(x,y,0)^{\gamma},\end{array}$
(3.6)

and the boundary conditions are periodic in both directions. The exact
solution of above Euler equations is given in [14]. When $T=3$, its solution
is still smooth, and now the corresponding errors and orders obtained by the
HWENO6-M5-I and HWENO6-M5-NI schemes are listed in the Table 3.4.

Table 3.4:  2D-Euler equations: The initial condition $\rho(x,y,0)=\frac{1+0.2\sin(\frac{x+y}{2})}{\sqrt{6}}$, $\mu(x,y,0)=\nu(x,y,0)=\sqrt{\frac{\gamma}{2}}\rho(x,y,0)$ and $p(x,y,0)=\rho(x,y,0)^{\gamma}$. Periodic boundary condition. $T=3$. HWENO6-M5-I and HWENO6-M5-NI schemes. $L^{1}$ and $L^{\infty}$. | HWENO6-M5-I scheme | HWENO6-M5-NI scheme
---|---|---
grid points |  $L^{1}$ error | order |  $L^{\infty}$ error | order |  $L^{1}$ error | order |  $L^{\infty}$ error | order
10$\times$ 10 | 6.56E-03 |  | 1.74E-02 |  | 6.70E-03 |  | 2.20E-02 |
20$\times$ 20 | 2.54E-04 | 4.69 | 1.78E-03 | 3.29 | 4.20E-04 | 4.00 | 2.73E-03 | 3.01
40$\times$ 40 | 1.13E-05 | 4.50 | 1.13E-04 | 3.97 | 2.47E-05 | 4.09 | 2.23E-04 | 3.61
80$\times$ 80 | 2.70E-07 | 5.38 | 6.06E-06 | 4.22 | 7.57E-07 | 5.03 | 1.47E-05 | 3.93
160$\times$160 | 4.18E-09 | 6.01 | 1.14E-07 | 5.73 | 2.14E-08 | 5.14 | 4.76E-07 | 4.95
320$\times$320 | 4.67E-11 | 6.49 | 1.37E-09 | 6.39 | 6.24E-10 | 5.10 | 1.30E-08 | 5.20

Comment: According to the results listed in above four tables Table 3.1, Table
3.2, Table 3.3 and Table 3.4, we can see that the HWENO6-M5-I scheme can reach
sixth order of accuracy, this is because the reconstruction process is of the
sixth order and the modification procedure has not been enacted for these
smooth cases. If we modify the first order moments of all the cells without
judgment, we can see that the HWENO6-M5-NI scheme can reach fifth order of
accuracy as expected, this is because the first order moments of all the cells
are modified by a quartic polynomial, which leads to the decrease in the order
of accuracy.

From Fig 3.1 and Fig 3.2, we can observe that both the HWENO6 and HWENO6-M5-I
schemes work well in comparison with the exact solutions and there is not much
difference between the results of these two schemes.

Example 3.5. The Lax problem: one-dimensional Euler equations (3.3) with the
Riemann initial condition:

$(\rho,\mu,p)^{T}=\left\\{\begin{array}[]{c}(0.445,0.698,3.528)^{T},\ \
-0.5<x<0,\\\ (0.5,0,0.571)^{T},\ \ 0<x<0.5.\end{array}\right.$ (3.7)

The computed result of the density $\rho$ and its zoom-in picture obtained by
the HWENO6 and HWENO6-M5-I schemes at the final time $T=0.16$ in comparison
with the reference exact solution, as well as the corresponding locations of
the troubled-cells for the HWENO6-M5-I scheme over time are plotted in the Fig
3.3.

Figure 3.3:  The Lax problem. $T=0.16$. Top: density, density zoomed in;
bottom: the locations of the troubled-cells. Solid line: the exact solution;
squares: the result of the HWENO6 scheme; triangles: the result of the
HWENO6-M5-I scheme. Number of cells: 200.

Example 3.6. The shock density wave interaction problem: one-dimensional Euler
equations (3.3) with a moving Mach=3 shock interaction containing sine waves
in the density:

$(\rho,\mu,p)^{T}=\left\\{\begin{array}[]{c}(3.857143,2.629369,10.333333)^{T},\
\ -5<x<-4,\\\ \big{(}1+0.2\sin(5x),0,1\big{)}^{T},\ \
-4<x<5.\end{array}\right.$ (3.8)

The computational result of the density $\rho$ and its zoom-in picture
obtained by the HWENO6 and HWENO6-M5-I schemes at the final time $T=1.8$ in
comparison with the reference “exact” solution (which is a numerically
converged solution computed by the fifth-order finite difference WENO scheme
[12] with 8000 grid points), as well as the corresponding locations of the
troubled-cells for the HWENO6-M5-I scheme over time are plotted in the Fig
3.4.

Figure 3.4:  The shock density wave interaction problem. $T=1.8$. Top:
density, density zoomed in; bottom: the locations of the troubled-cells. Solid
line: the exact solution; squares: the result of the HWENO6 scheme; triangles:
the result of the HWENO6-M5-I scheme. Number of cells: 400.

Example 3.7. The blast wave problem: one-dimensional Euler equations (3.3)
with the initial condition:

$(\rho,\mu,p)^{T}=\left\\{\begin{array}[]{c}(1,0,10^{3})^{T},\ \ 0<x<0.1,\\\
(1,0,10^{-2})^{T},\ \ 0.1<x<0.9,\\\ (1,0,10^{2})^{T},\ \
0.9<x<1.\end{array}\right.$ (3.9)

The computational result of the density $\rho$ and its zoom-in picture
obtained by the HWENO6 and HWENO6-M5-I schemes at the final time $T=0.038$ in
comparison with the reference “exact” solution(which is a numerically
converged solution computed by the fifth-order finite difference WENO scheme
[12] with 16000 grid points), as well as the corresponding locations of the
troubled-cells for the HWENO6-M5-I scheme over time are plotted in the Fig
3.5.

Figure 3.5:  The blast wave problem. $T=0.038$. Top: density, density zoomed
in; bottom: the locations of the troubled-cells. Solid line: the exact
solution; squares: the result of the HWENO6 scheme; triangles: the result of
the HWENO6-M5-I scheme. Number of cells: 800.

Example 3.8. The Sedov blast wave problem: one-dimensional Euler equations
(3.3) with the initial condition:

$(\rho,\mu,E)^{T}=\left\\{\begin{array}[]{c}(1,0,10^{-12})^{T},\ \
x\in[-2,2]\setminus the\ \ center\ \ cell,\\\ (1,0,\frac{3200000}{\triangle
x})^{T},\ \ x\in the\ \ center\ \ cell.\end{array}\right.$ (3.10)

The computed results of the density $\rho$, the velocity $\mu$ and the
pressure $p$ obtained by the HWENO6 and HWENO6-M5-I schemes at the final time
$T=0.001$ in comparison with the reference exact solution, as well as the
corresponding locations of the troubled-cells for the HWENO6-M5-I scheme over
time are plotted in the Fig 3.6.

Figure 3.6:  The Sedov blast wave problem. $T=0.001$. Top: density, velocity;
bottom: pressure, the locations of the troubled-cells. Solid line: the exact
solution; squares: the result of the HWENO6 scheme; triangles: the result of
the HWENO6-M5-I scheme. Number of cells: 400.

Example 3.9. Double Mach reflection problem: two-dimensional Euler equations
(3.5) in a computational domain $[0,4]$$\times$$[0,1]$ with a reflection wall
lying at the bottom, starting from the position $(x,y)=(\frac{1}{6},0)$,
making a $60^{\circ}$ angle with the $x$-axis. For the bottom of the domain,
at the reflection wall the reflection boundary condition is applied, and at
the rest of the bottom the exact post-shock condition is imposed. For the top
of the domain, the corresponding boundary condition is exactly the motion of a
Mach 10 shock with $\gamma=1.4$. The contours of the computed density $\rho$
and its blow-up region around the double Mach stem obtained by the HWENO6-M5-I
scheme (the reason we do not present the result of HWENO6 scheme is that this
scheme for the double Mach reflection problem will blow up without
modification), as well as the corresponding locations of the troubled-cells
for the HWENO6-M5-I scheme at the final time $T=0.2$ are plotted in the Fig
3.7.

Figure 3.7:  Double Mach reflection problem. $T=0.2$. 30 equally spaced
density contours from 1.5 to 22.7. From top to bottom: density contours, zoom-
in density contours around the Mach stem, the locations of the troubled-cells
at the final time of the HWENO6-M5-I scheme. Number of cells: 1600$\times$400
in the region of $[0,4]$$\times$$[0,1]$.

Example 3.10. Forward step problem: two-dimensional Euler equations (3.5) in a
one length unit wide and three length units long wind tunnel with a 0.2 length
units high step located 0.6 length units from the left side of the tunnel. At
the beginning, we initialize this problem by a right-going Mach 3 flow. Along
the wall of the tunnel the reflection boundary condition is applied, and at
the entrance the inflow boundary condition is imposed, while at the exit the
outflow boundary condition is imposed. The contours of the computed density
$\rho$ obtained by the HWENO6 and HWENO6-M5-I schemes, as well as the
corresponding locations of the troubled-cells for the HWENO6-M5-I scheme at
the final time $T=4$ are plotted in the Fig 3.8.

Figure 3.8:  Forward step problem. $T=4$. 30 equally spaced density contours
from 0.32 to 6.15. From top to bottom: density contours of the HWENO6 scheme,
density contours of the HWENO6-M5-I scheme, the locations of the troubled-
cells at the final time of the HWENO6-M5-I scheme. Number of cells:
600$\times$200 in the region of $[0,3]$$\times$$[0,1]$.

Comment: From above six examples, we can see that the results of the
HWENO6-M5-I scheme have much better resolutions and sharper shock transitions
than those of the HWENO6 scheme without the modification procedure. This might
be due to the fact that after modifying the first order moments of the
troubled-cells, the proportion of the last two layers is much higher than that
of the first two layers in the reconstruction process. Also, this modification
procedure can increase the stability of our scheme according to the double
Mach reflection problem.

##  4 Concluding remarks

In this paper, we have designed a high-order moment-based multi-resolution
HWENO scheme for hyperbolic conservation laws in the one and two dimensional
cases on structured meshes. In comparison with our previous work in [14], the
new feature of this HWENO scheme is that the zeroth and first order moments
rather than the first order derivative are used in the spacial reconstruction
algorithm, and only the function values of the Gauss-Lobatto points in one or
two dimensional case are needed to be reconstructed. Also, after the
reconstruction algorithm, an extra modification procedure is used to modify
those first order moments of the troubled-cells and the corresponding Gauss-
Lobatto point values of these troubled-cells need to be updated by repeating
the reconstruction algorithm, to enhance both resolution and stability. At the
same time, the linear weights can also be any positive numbers as long as
their sum equals one and the CFL number can still be 0.6 for both the one and
two dimensional cases. This HWENO scheme is achieved by reconstructing the
Gauss-Lobatto point values, modifying the first order moments of those cells
which are identified to be troubled-cells by the KXRCF troubled-cell indicator
and repeating the reconstruction algorithm to update the corresponding Gauss-
Lobatto point values of these troubled-cells. In comparison with the multi-
resolution WENO scheme, our major advantages are still the compactness of the
stencils and smaller errors under the same meshes when the same order is
fixed. The framework of this moment-based multi-resolution HWENO scheme would
be particularly efficient and simple on unstructured meshes, the study of
which is our ongoing work.

Acknowledgments: The research of C.-W. Shu is partly supported by AFOSR grant
FA9550-20-1-0055 and NSF grant DMS-2010107. The research of J. Li and J. Qiu
is partly supported by NSFC grant 12071392.

## References

  * [1] D.S. Balsara, S. Garain, V. Florinski and W. Boscheri, An efficient class of WENO schemes with adaptive order for unstructured meshes, J. Comput. Phys., 404 (2020), 109062.
  * [2] D.S. Balsara, S. Garain and C.-W. Shu, An efficient class of WENO schemes with adaptive order, J. Comput. Phys., 326 (2016), 780-804.
  * [3] D.S. Balsara, C.-W. Shu, Monotonicity preserving weighted essentially non-oscillatory schemes with increasingly high order of accuracy, J. Comput. Phys., 160 (2000), 405-452.
  * [4] R. Biswas, K.D. Devine, and J. Flaherty, Parallel, adaptive finite element methods for conservation laws, Appl. Numer. Math., 14 (1994), pp. 255-283.
  * [5] A. Burbeau, P. Sagaut, and C.H. Bruneau, A problem-independent limiter for high-order Runge-Kutta discontinuous Galerkin methods, J. Comput. Phys., 169 (2001), pp. 111-150.
  * [6] M. Castro, B. Costa, W.S. Don, High order weighted essentially non-oscillatory WENO-Z schemes for hyperbolic conservation laws, J. Comput. Phys., 230 (2011), 1766-1792.
  * [7] B. Cockburn and C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Math. Comp., 52 (1989), pp. 411-435.
  * [8] M. Dumbser, W. Boscheri, M. Semplice and G. Russo, Central WENO schemes for hyperbolic conservation laws on fixed and moving unstructured meshes, SIAM J. Sci. Comput., 39 (2017), A2564-A2591.
  * [9] A. Harten, ENO schemes with subcell resolution, J. Comput. Phys., 83 (1989), pp. 148-184.
  * [10] A. Harten, B. Engquist, S. Osher, S. Chakravarthy, Uniformly high order accurate essentially non-oscillatory schemes, III, J. Comput. Phys., 71 (1987), 231-303.
  * [11] C. Hu, C.-W. Shu, Weighted essentially non-oscillatory schemes on triangular meshes, J. Comput. Phys., 150 (1999), 97-127.
  * [12] G.-S. Jiang, C.-W. Shu, Efficient implementation of weighted ENO schemes, J. Comput. Phys., 126 (1996), 202-228.
  * [13] L. Krivodonova, J. Xin, J.-F. Remacle, N. Chevaugeon, J.E. Flaherty, Shock detection and limiting with discontinuous Galerkin methods for hyperbolic conservation laws, Appl. Numer. Math. 48 (2004) 323-338.
  * [14] J. Li, J. Qiu, C.-W. Shu, Multi-resolution HWENO schemes for hyperbolic conservation laws, J. Comput. Phys., 446 (2021) 110653.
  * [15] H. Liu, J. Qiu, Finite difference Hermite WENO schemes for conservation laws, J. Sci. Comput., 63 (2015), 548-572.
  * [16] H. Liu, J. Qiu, Finite difference Hermite WENO schemes for conservation laws, II: an alternative approach, J. Sci. Comput., 66 (2016), 598-624.
  * [17] X.D. Liu, S. Osher, T. Chan, Weighted essentially non-oscillatory schemes, J. Comput. Phys., 115 (1994), 200-212.
  * [18] S. Pirozzoli, Conservative hybrid compact-WENO schemes for shock-turbulence interaction, J. Comput. Phys., 178 (2002), 81-117.
  * [19] J. Qiu, C.-W. Shu, Hermite WENO schemes and their application as limiters for Runge-Kutta discontinuous Galerkin method: one-dimensional case, J. Comput. Phys., 193 (2004), 115-135.
  * [20] J. Qiu, C.-W. Shu, A comparison of troubled-cell indicators for Runge-Kutta discontinuous Galerkin methods using weighted essentially nonoscillatory limiters, SIAM J. Sci. Comput. 27 (2005) 995-1013.
  * [21] W.J. Rider and L.G. Margolin, Simple modifications of monotonicity-preserving limiters, J. Comput. Phys., 174 (2001), pp. 473-488.
  * [22] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, J. Comput. Phys., 77 (1988), 439-471.
  * [23] C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock capturing schemes, II, J. Comput. Phys., 83 (1989), 32-78.
  * [24] A. Suresh and H.T. Huynh, Accurate monotonicity-preserving schemes with Runge-Kutta time stepping, J. Comput. Phys., 136 (1997), pp. 83-99.
  * [25] Z. Zhao, Y. Chen, J. Qiu, A hybrid Hermite WENO scheme for hyperbolic conservation laws, J. Comput. Phys. 405 (2020) 109175.
  * [26] Z. Zhao, J. Qiu, A Hermite WENO scheme with artificial linear weights for hyperbolic conservation laws, J. Comput. Phys., 417 (2020), 109583.
  * [27] J. Zhu, J. Qiu, A class of fourth order finite volume Hermite weighted essentially non-oscillatory schemes, Sci. China Ser. A, Math., 51 (2008), 1549-1560.
  * [28] J. Zhu, C.-W. Shu, A new type of multi-resolution WENO schemes with increasingly higher order of accuracy, J. Comput. Phys., 375 (2018), 659-683.
  * [29] J. Zhu, C.-W. Shu, A new type of multi-resolution WENO schemes with increasingly higher order of accuracy on triangular meshes, J. Comput. Phys., 392 (2019), 19-33.


# Exceptional surgeries in $3$–manifolds

Kenneth L. Baker  Department of Mathematics
University of Miami
Coral Gables, FL 33146
USA<EMAIL_ADDRESS>and  Neil R. Hoffman  Department of Mathematics
Oklahoma State University
Stillwater, OK 74078
USA<EMAIL_ADDRESS>
###### Abstract.

Myers shows that every compact, connected, orientable $3$–manifold with no
$2$–sphere boundary components contains a hyperbolic knot. We use work of
Ikeda with an observation of Adams-Reid to show that every $3$–manifold
subject to the above conditions contains a hyperbolic knot which admits a non-
trivial non-hyperbolic surgery, a toroidal surgery in particular. We conclude
with a question and a conjecture about reducible surgeries.

Myers shows that there are hyperbolic knots in every compact, connected,
orientable $3$–manifold whose boundary contains no $2$–spheres [Mye93]. Might
there be such a $3$–manifold for which every hyperbolic knot has no non-
trivial exceptional surgeries?

One approach to showing the answer is Yes would be to prove that there exists
a $3$–manifold in which every hyperbolic knot has cusp volume larger than $18$
so that the $6$–Theorem [Ago00, Lac00] would obstruct any exceptional surgery.
However, [ACF+06, Corollary 5.2] implies that every closed, connected,
orientable $3$–manifold contains infinitely many hyperbolic knots with cusp
volume at most $9$. So this approach will not work. Furthermore, the knots
constructed in [ACF+06, Corollary 5.2] do not necessarily have any exceptional
surgery, so that work does not address our question.

In this short note we demonstrate the answer to the question is actually No by
constructing hyperbolic knots with a non-trivial toroidal surgery in any
$3$–manifold.

######  Theorem 1.

Let $M$ be a compact, connected, orientable $3$–manifold such that $\partial
M$ contains no $2$–spheres. There exist infinitely many hyperbolic knots in
$M$ that admit a toroidal surgery.

###### Proof.

Let $M$ be a compact, connected orientable $3$–manifold whose boundary
contains no $2$–spheres. In [Ike12], Ikeda shows that $M$ contains an infinite
family of embedded genus $2$ handlebodies in $M$, each with hyperbolic and
anannular complement of its interior where its genus $2$ boundary is totally
geodesic. Let $H$ be any one of these handlebodies.

In Lemma 2 we find a knot $K$ in $H$ that bounds an embedded once-punctured
Klein bottle $\Sigma$ such that $H-K$ is a one-cusped anannular hyperbolic
manifold in which $\partial H$ is totally geodesic. Therefore $M-K$ decomposes
along $\partial H$ into two anannular hyperbolic manifolds. Thus, following an
observation of Adams and Reid [AR93, Observation 2.1], $M-K$ is a hyperbolic
manifold containing a quasi-Fuchsian surface isotopic to $\partial H$, and $K$
is a hyperbolic knot in $M$. (Note that while $\partial H$ is totally geodesic
in both $M-int(H)$ and $H-K$, its hyperbolic structure may not be the same in
these two manifolds. Hence $\partial H$ is not necessarily totally geodesic in
$M-K$.)

Since $K$ bounds the once-punctured Klein bottle $\Sigma$, surgery on $K$
along the slope $\sigma$ of $\partial\Sigma$ produces a manifold
$M_{K}(\sigma)$ containing an embedded Klein bottle $\widehat{\Sigma}$. The
manifold $M_{K}(\sigma)$ will be toroidal unless the torus
$\partial\mathcal{N}(\widehat{\Sigma})$ compresses. However, Lemma 2 shows
that $K$ may be further chosen in $M$ so that $H_{K}(\sigma)-\widehat{\Sigma}$
is also a one-cusped anannular hyperbolic manifold in which $\partial H$ is
totally geodesic boundary. Therefore, as $M_{K}(\sigma)-\widehat{\Sigma}$
decomposes along $\partial H$ into the hyperbolic manifolds $M-int(H)$ and
$H_{K}(\sigma)-\widehat{\Sigma}$, it follows that
$\partial\mathcal{N}(\widehat{\Sigma})$ must be incompressible in
$M_{K}(\sigma)$. ∎

######  Lemma 2.

There is a knot $K$ in a genus $2$ handlebody $H$ that bounds a once-punctured
Klein bottle $\Sigma$ so that $H-K$ is a one-cusped anannular hyperbolic
manifold in which $\partial H$ is totally geodesic. Hence surgery on $K$ along
the slope $\sigma$ of $\partial\Sigma$ produces a manifold $H_{K}(\sigma)$
containing an embedded Klein bottle $\widehat{\Sigma}$. Furthermore, $K$ may
be chosen so that $H_{K}(\sigma)-\widehat{\Sigma}$ is also a one-cusped
anannular hyperbolic manifold in which $\partial H$ is totally geodesic.

###### Proof.

Figure 2(a) shows a surgery description of a trivial $3$–strand tangle in the
ball, along with an arc $k$ that has its endpoints on the tangle strands.
Figure 2(b) shows the result of an isotopy in which the tangle is more
obviously trivial at the expense of elongating the arc $k$. The double
branched cover of this trivial $3$–strand tangle is a handlebody $H$ in which
the arc $k$ lifts to a knot $K$. Figure 2(c),(d), and (e) illustrate the
construction of the knot $K$ in the handlebody $H$. In (c), two caps with red
dual arcs are attached to the $3$–strand tangle to form a trivial $1$–strand
tangle in the ball. After straightening the strand in (d), the double branched
cover is taken in (e). The two caps each lift to $2$–handles attached to $H$.
The two red arcs in (e) are the co-cores of these two $2$–handles, so $H$ is
obtained by drilling them out. For lifting the surgery description, note that
a curve linking the branch locus once with surgery coefficient $1/2a$ lifts to
a single curve with surgery coefficient $1/a$. Thus for each pair of integers
$n,m$, we obtain a knot $K$ in a genus $2$ handlebody $H$.

Figure 2(f) shows a surgery description of the double of $(H,K)$ across
$\partial H$, the link $K\cup\overline{K}$ in $H\cup\overline{H}=S^{1}\times
S^{2}\\#S^{1}\times S^{2}$, obtained by mirroring Figure 2(e) and performing
$0$ surgery on the components formed from the co-cores of the $2$–handles and
their mirrors. A verified computation in SnapPy [CDGW] confirms that the
complement of the link $K\cup\overline{K}$ in $S^{1}\times S^{2}\\#S^{1}\times
S^{2}$ is hyperbolic for choices of $n,m\in\mathbb{Z}$ with $|n|$ and $|m|$
suitably large. (More specifically, a verified computation in SnapPy shows
that after doing the two $0$-surgeries on the two red components in Figure
2(f), the resulting $6$–component link in $S^{1}\times S^{2}\\#S^{1}\times
S^{2}$ has a hyperbolic complement. Then there is a constant $N$ such that the
$2$–component link complement resulting from the surgeries on the green and
purple components will be hyperbolic if both $|n|>N$ and $|m|>N$; see [Koj88,
Lemma 5] or [BHL19, Theorem 3.1].) Since the double has the reflective
symmetry in which $\partial H$ is the fixed set, it must be a totally geodesic
surface. Hence $H-K$ must be a one-cusped anannular hyperbolic manifold in
which $\partial H$ is totally geodesic.

In Figure 2(e) one observes that $K$ bounds a once-punctured Klein bottle
$\Sigma$ in $H$ that is disjoint from the two curves of the surgery
description. As such, Dehn surgery on $K$ in $H$ along the boundary slope
$\sigma=\partial\Sigma$ produces the manifold $H_{K}(\sigma)$ which contains
the Klein bottle $\widehat{\Sigma}$ obtained by capping off $\Sigma$ with a
meridional disk of the surgery.

All that remains is to show that $\widehat{\Sigma}$ is essential in the
filling. First, we may understand the complement of $\widehat{\Sigma}$ through
tangles. As apparent in Figure 2(e), the surface $\Sigma$ may be taken to be
invariant under the involution of $H$ from the branched covering so that the
fixed set intersects $\Sigma$ in two points and an arc. Then $\Sigma$ descends
to a disk $D$ containing the arc $k$ in its boundary and meeting the branch
locus in the remainder of its boundary and two points in its interior. This
disk $D$ may be tracked from its initial quotient of $\Sigma$ in Figure 2(d)
back to Figure 2(a). Now Figure 1(a) shows the exterior of the arc $k$ while
Figure 1(b) shows the rational tangle filling associated to $\sigma$–framed
surgery on $K$. In particular, the disk $D-k$ is completed to a disk
$\widehat{D}$ containing the closed component of the branch locus as its
boundary and meeting the strands of the branch locus in two interior points.
Indeed, the double branched cover of the tangle Figure 1(b) is the manifold
$H_{K}(\sigma)$ in which $\widehat{D}$ lifts to $\widehat{\Sigma}$. Finally,
Figure 1(c) shows the tangle that is the complement of a small regular
neighborhood of $\widehat{D}$.

Figure 3(a) shows a rational tangle filling of Figure 1(c) with the arc
$k^{\prime}$ that is the core of the rational tangle. This $3$–strand tangle
is a trivial tangle as made more apparent in Figures 3(b), (c), and (d) which
isotop the tangle while elongating arc $k^{\prime}$. As before, (c) shows the
attachment of two caps with dual arcs and (d) straightens the resulting
$1$–strand tangle. Figure 3(e) shows the double branched cover which
illustrates the lift of the arc $k^{\prime}$ as the knot $K^{\prime}$ in
another genus $2$ handlebody $H^{\prime}$. Again, the two caps each lift to
$2$–handles attached to $H^{\prime}$, the two red arcs in (e) are the co-cores
of these two $2$–handles, and so $H^{\prime}$ is obtained by drilling them
out. Note that the knot $K^{\prime}$ in $H^{\prime}$ depends on the previously
chosen pair of integers $n,m$ of the surgery description.

It now follows that, by construction, $H_{K}(\sigma)-\widehat{\Sigma}$ is
homeomorphic to $H^{\prime}-K^{\prime}$. We show that $H^{\prime}-K^{\prime}$
is a one-cusped anannular hyperbolic manifold in which $\partial H^{\prime}$
is totally geodesic just as we did for $H-K$. Figure 3(f) shows a surgery
description of the double of $(H^{\prime},K^{\prime})$ across $\partial
H^{\prime}$, the link $K^{\prime}\cup\overline{K^{\prime}}$ in
$H^{\prime}\cup\overline{H^{\prime}}=S^{1}\times S^{2}\\#S^{1}\times S^{2}$,
obtained by mirroring Figure 3(e) and performing $0$ surgery on the components
formed from the co-cores of the $2$–handles and their mirrors. A verified
computation in SnapPy [CDGW] confirms that the complement of the link
$K^{\prime}\cup\overline{K^{\prime}}$ in $S^{1}\times S^{2}\\#S^{1}\times
S^{2}$ is hyperbolic if both $|n|$ and $|m|$ are suitably large. Since the
double has the reflective symmetry in which $\partial H^{\prime}$ is the fixed
set, it must be a totally geodesic surface. Hence $H^{\prime}-K^{\prime}$ is a
one-cusped anannular hyperbolic manifold in which $\partial H^{\prime}$ is
totally geodesic.

Since $H_{K}(\sigma)-\widehat{\Sigma}\cong H^{\prime}-K^{\prime}$, we obtain
the desired results whenever $|n|$ and $|m|$ are large enough to be suitably
large in both situations. ∎

######  Remark 3.

To give a concrete example, taking $n=m=1$ is sufficient for the knots
$K\subset H$ and $K^{\prime}\subset H^{\prime}$ to be hyperbolic. Certainly,
one could verify the hyperbolicity of these knots by hand in the spirit of
what was done in [AR93], but the argument would take longer. Hence we content
ourselves with verified computations in SnapPy [CDGW].

Figure 1. (b) A surgery description of a $3$–strand tangle in the ball with an
unknot component that bounds a disk intersected twice by the strands. (a) The
complement of a rational tangle in this tangle. (c) The complement of a small
neighborhood of the disk bounded by the unknot component. Figure 2. (a) A
rational tangle filling of Figure 1(a) with its core arc $k$. (b) & (c) An
isotopy showing the filled tangle is a rational $3$–strand tangle. The arc $k$
is carried along. (c) Attached to the rational $3$–strand tangle are two caps
with their dual arcs to form a $1$–strand tangle in the ball. (d) The
$1$–strand tangle is straightened. (e) The double branched cover is formed.
Drilling out the red arcs leaves a genus $2$ handlebody $H$ containing the
knot $K$ that covers $k$. (f) A surgery description of the double of $(H,K)$
is formed.  Figure 3. (a) A rational tangle filling of Figure 1(c) with its
core arc $k^{\prime}$. (b) & (c) An isotopy showing the filled tangle is a
rational $3$–strand tangle. The arc $k^{\prime}$ is carried along. (c)
Attached to the rational $3$–strand tangle are two caps with their dual arcs
to form a $1$–strand tangle in the ball. (d) The $1$–strand tangle is
straightened. (e) The double branched cover is formed. Drilling out the red
arcs leaves a genus $2$ handlebody $H^{\prime}$ containing the knot
$K^{\prime}$ that covers $k^{\prime}$. (f) A surgery description of the double
of $(H^{\prime},K^{\prime})$ is formed.

What can be said about other kinds of exceptional surgeries? Considerations of
Betti numbers show that many closed, compact, orientable $3$–manifolds cannot
contain a knot with a Dehn surgery to a lens space or a small Seifert fibered
space. In light of the Cabling Conjecture [GAnS86] whose proof would imply
that no hyperbolic knot in $S^{3}$ has a reducible surgery, it is reasonable
to expect that there are $3$–manifolds in which no hyperbolic knot admits a
reducible surgery. However, we are presently unaware of any $3$–manifold known
to not have a hyperbolic knot with a non-trivial reducible surgery.

######  Question 4.

Which compact, connected, orientable $3$–manifolds do not contain a hyperbolic
knot with a non-trivial reducible surgery?

While non-trivial reducible surgeries on hyperbolic knots in reducible
manifolds do exist, see e.g. [HM03], we suspect that manifolds whose prime
decompositions have at least $3$ summands are candidates.

######  Conjecture 5.

A closed orientable $3$–manifold with at least $3$ summands does not contain a
hyperbolic knot with a non-trivial reducible surgery.

Towards the conjecture, suppose $K$ is a hyperbolic knot in a closed
orientable $3$–manifold $M$ with at least $3$ summands. One may hope that each
planar meridional surfaces in the knot complement $M-K$ arising from $K$
intersecting multiple reducing spheres would contribute a certain amount to
the length of the shortest longitude of $K$. From this, at least if $M$ had
sufficiently many summands, one would be able to use the 6-Theorem to obstruct
a non-trivial reducible surgery. However this would also obstruct a toroidal
surgery contrary to Theorem 1. Indeed, it would also contradict [ACF+06,
Corollary 5.2] which shows that the topology of $M$ cannot force all
longitudes of hyperbolic knots in $M$ to be long.

On the other hand, combinatorial structures in knot complements can induce
obstructions. For instance, hyperbolic alternating knots in $S^{3}$ that have
at least $9$ twist regions (in twist-reduced diagrams) provide an obstruction
the existence of non-trivial exceptional fillings; see [Lac00, Theorem 5.1].

######  Acknowledgments 6.

KB thanks Jacob Caudell for conversations related to [Cau21, Conjecture 5]
that prompted this note. This work was partially supported by Simons
Foundation grant #209184 to Kenneth Baker and by Simons Foundation grant
#524123 to Neil Hoffman.

## References

  * [ACF+06] C. Adams, A. Colestock, J. Fowler, W. Gillam, and E. Katerman.  Cusp size bounds from singular surfaces in hyperbolic 3-manifolds.  Trans. Amer. Math. Soc., 358(2):727–741, 2006.
  * [Ago00] Ian Agol.  Bounds on exceptional Dehn filling.  Geom. Topol., 4(1):431–449, 2000.
  * [AR93] Colin C. Adams and Alan W. Reid.  Quasi-Fuchsian surfaces in hyperbolic knot complements.  Journal of the Australian Mathematical Society. Series A. Pure Mathematics and Statistics, 55(1):116–131, 1993.
  * [BHL19] Kenneth L Baker, Neil R Hoffman, and Joan E Licata.  Jointly primitive knots and surgeries between lens spaces.  arXiv preprint arXiv:1904.03268, 2019.  To appear in Communications in Analysis and Geometry.
  * [Cau21] Jacob Caudell.  Three lens space summands from the Poincaré homology sphere, 2021.  arXiv:2101.01256.
  * [CDGW] Marc Culler, Nathan M. Dunfield, Matthias Goerner, and Jeffrey R. Weeks.  SnapPy, a computer program for studying the geometry and topology of $3$-manifolds.  Available at http://snappy.computop.org (17/01/2021).
  * [GAnS86] Francisco González-Acuña and Hamish Short.  Knot surgery and primeness.  Math. Proc. Cambridge Philos. Soc., 99(1):89–102, 1986.
  * [HM03] James A. Hoffman and Daniel Matignon.  Examples of bireducible Dehn fillings.  Pacific J. Math., 209(1):67–83, 2003.
  * [Ike12] Toru Ikeda.  Hyperbolic spatial graphs in 3-manifolds.  Topology and its Applications, 159(1):279–282, 2012.
  * [Koj88] Sadayoshi Kojima.  Isometry transformations of hyperbolic 3-manifolds.  Topology and its Applications, 29(3):297–307, 1988.
  * [Lac00] Marc Lackenby.  Word hyperbolic Dehn surgery.  Inventiones mathematicae, 140(2):243–282, 2000.
  * [Mye93] Robert Myers.  Excellent $1$-manifolds in compact $3$-manifolds.  Topology Appl., 49(2):115–127, 1993.