∎

11institutetext: M. Pablo $\cdot$ O. Reula 22institutetext: Instituto de
Física Enrique Gaviola, CONICET, FaMAF - Universidad Nacional de Córdoba,
Ciudad Universitaria, 5000 Córdoba, Argentina
22email<EMAIL_ADDRESS>
22email<EMAIL_ADDRESS>
# New Numerical Interface Scheme for the Kurganov-Tadmor second-order Method

Pablo Montes  Oscar Reula

(Received: date / Accepted: date)

###### Abstract

In this paper, we develop a numerical scheme to handle interfaces across
computational domains in multi-block schemes for the approximation of systems
of conservation laws. We are interested in transmitting shock discontinuities
without lowering the overall precision of the method. We want to accomplish
this without using information from interior points of adjacent grids, that
is, sharing only information from boundary points of those grids. To achieve
this, we choose to work with the second-order Kurganov-Tadmor (KT) method at
interior points, relaxing it to first order at interfaces. This allows us to
keep second-order overall accuracy (in the relevant norm) and at the same time
preserve the TVD property of the original scheme. After developing the method
we performed several standard one and two-dimensional tests. Among them, we
used the one-dimensional advection and Burgers equations to verify the second-
order convergence of the method. We also tested the two-dimensional Euler
equations with an implosion and a Gresho vortexliska2003 . In particular, in
the two-dimensional implosion test we can see that regardless of the
orientation of shocks with respect to the interface, they travel across them
without appreciable deformation both in amplitude and front direction.

###### Keywords:

multigrid method conservative method numerical interface Kurganov-Tadmor
method.

##  1 Introduction

Many physical systems of interest that can be described by partial
differential equations (PDEs) cannot be numerically approximated by using a
single numerical grid. This may be due to several reasons, for example: a) the
topology of the problem does not allow it, as is the case of fields defined on
a sphere, b) the size of the problem calls for block-parallelization, c) the
need for greater resolution in certain parts of the domain. This leads to the
creation of methods that allow for multiple grid patches to cover the
computational domain. One such type of method is that of multi-block domain
decomposition, in which grids are non-overlapping and only connect at boundary
regions called interfaces. On those interfaces, the PDE discretization must be
specially tailored to allow the correct and accurate transmission of numerical
solutions between grids.

The treatment of these interfaces is generally well understood for smooth
solutions. In the finite-difference case satisfying the summation by parts
(SBP) property, simultaneous approximation terms (SAT) schemes can be used in
a very efficient way and with very high accuracy CARPENTER1999341 .
Alternatively, pseudo-spectral methods are also an optimal choice and similar
penalty methods can be applied HESTHAVEN200023 . However, in systems where
shocks or other types of discontinuities are expected, the application of the
above-mentioned schemes is no longer possible. This is mostly due to the
formation of spurious oscillations on the numerical approximation of the
solution near discontinuities.

A very large and important subset of systems that present shocks and which are
ubiquitous in physics are the so-called systems of conservation laws, which
describe the evolution of the density fields of conserved quantities, such as
energy, momentum, or mass; that is, fields whose variation satisfies a flux
law. Their general form is

$\frac{\partial}{\partial t}u+\sum_{j}\frac{\partial}{\partial
x_{j}}F^{j}(u)=0,\quad u|_{\partial\Omega,t=0}=u_{0},$ (1)

where $u$ is a vector of conserved quantities, $F^{j}$ are flux functions for
each coordinate, $\Omega$ is the computational domain, and $\partial\Omega$
its boundary. Although these systems present generic shock discontinuities,
their structure allows for very good and stable approximations of solutions
via so called numerical conservative methods. These are methods that use the
weak formulation of Eq. (1) to numerically enforce the conservation laws,
allowing high orders of convergence in the smooth parts of the solution while
keeping good accuracy near shock discontinuities. They also avoid spurious
oscillations near shocks by using flux limiters that are built to preserve
certain properties of the exact solution. Some of them, for example, preserve
monotonicity, others, guarantee that the total variation of the solution does
not increase LeVeque-1992 . To achieve this, conservative methods determine
the shock region and locally increase the numerical dissipation. In the cases
of interfaces among multi-blocks, the complication resides in the fact that
the information needed to properly determine a shock region and evolve it
accordingly is usually located in an adjacent grid. Nevertheless, if we are
willing to tolerate a first-order scheme in the interface, the complication is
surmountable. Since this region’s size is similar to that of a grid cell, the
overall global order of convergence can be one order higher. There are
previous works that present ways to treat these interfaces. Their common
feature is the use of information from the interior points of adjacent grids;
this requires overlapping between grids, thus needing some interpolation and
more communication among different processes. Some examples of this can be
found in MCCORQUODALE2015181 ; Chesshire1994 ; SHERER2005459 ;
PARTENANDER1994551 . Here we are interested in developing a simple interface
method for conservative schemes which only uses information from the boundary
grid points. This is relevant for implementations in high-performance
computing codes, for which minimal communication between grids, and by
extension computational nodes, is desired. This will also allow the direct
treatment of situations where the traversal grid directions are not aligned
across interfaces, as could be the case of several grid patches covering a
sphere. In this case, it is not straightforward to do a flux reconstruction
using interior points from both sides of an interface.

In this work, we develop a method which overcomes the aforementioned
challenges and does comply with the desired properties for an extension of the
second-order Kurganov-Tadmor (KT) conservative scheme KT99 . Our method allows
the transmission of shock discontinuities through interfaces while maintaining
the global order of convergence and with minimal communication among grids,
since only the values at the interface boundary from each grid are shared.
Furthermore, it is also built as to preserve the total variation diminishing
(TVD) property of the original semi-discrete KT scheme. Contrary to upwind
methods, KT does not require a characteristic decomposition of the fields,
making it easier to implement.

This paper is organized as follows: In section 2, we give a brief overview of
the KT method and its total variation diminishing (TVD) property. In section
3, we introduce our multi-block version of the KT method and show that the TVD
property is still maintained. It works by falling back to a specific first-
order formulation near the boundary of each grid, and later calculating a cell
average between the values obtained in adjacent grids. Finally, in section 4,
we apply our scheme to the case of the advection and Burgers equations in one
dimension to test convergence and shock propagation. We also corroborate that
the second-order convergence is maintained and shock discontinuities and waves
are completely transmitted through interfaces without loss of neither phase
nor group speed. Finally, we show results for the implementation of this
interface for the two-dimensional Euler equations, solving a Gresho vortex and
an implosion problem from liska2003 .

##  2 Kurganov-Tadmor method

In this section we make a quick overview of the semi-discrete formulation of
the KT method. More information can be found in KT99 .

We start with a simple one-dimensional problem of the form

$\frac{\partial u}{\partial t}(x,t)+\frac{\partial f(u(x,t))}{\partial
x}=0\quad(x,t)\in[0,\infty)\times[0,\infty),$ (2)

where $u$ is a vector of conserved quantities and $f$ are the corresponding
fluxes.

We start by discretizing the spatial dimension in equal intervals of width
$\Delta x$, and denoting $x_{j}\doteq j\Delta x$. We then define the cell
average of $u$ at grid point $j$ as

$\overline{u}_{j}(t)=\frac{1}{\Delta
x}\int_{x_{(j-\frac{1}{2})}}^{x_{(j+\frac{1}{2})}}u(\xi,t)\;d\xi.$ (3)

If we now combine (2) and (3) we obtain an equivalent weak formulation of the
conservation law,

$\frac{\partial\overline{u}_{j}}{\partial
t}(x,t)=f\left(u(x_{j-1/2,},t)\right)-f\left(u(x_{j+1/2,},t)\right).$ (4)

Conservative methods evolve these cell averages of the solution using Eq. (4)
and thus need to approximate the solution at the cell boundaries $x_{j-1/2}$
and $x_{j+1/2}$. The main difference between different schemes lies in the way
that these boundaries values are approximated. From now onwards we denote
$v_{j}(t)$ as the numerical approximation of $\overline{u}(x_{j},t)$.

The semidiscrete KT method takes the conservative form

$\frac{\mathrm{d}}{\mathrm{d}t}v_{j}(t)=-\frac{H_{j+1/2}(t)-H_{j-1/2}(t)}{\Delta
x},$ (5)

$\displaystyle
H_{j+1/2}(t):=\frac{f(v^{+}_{j+1/2}(t))+f(v^{-}_{j+1/2}(t))}{2}-\frac{a_{j+1/2}(t)}{2}\left[v^{+}_{j+1/2}(t)-v^{-}_{j+1/2}(t)\right],$
(6)

$v^{+}_{j+1/2}(t)=v_{j+1}(t)-\frac{\Delta x}{2}(v_{x})_{j+1}(t)\quad
v^{-}_{j+1/2}(t)=v_{j}(t)+\frac{\Delta x}{2}(v_{x})_{j}(t),$ (7)

where $(v_{x})_{j}$ is some approximation of the spatial derivative of $v$,
$a_{j+1/2}(t)$ the maximum propagation speed of the solution between
$v^{-}_{j+1/2}(t)$ and $v^{+}_{j+1/2}(t)$,

$a_{j+1/2}(t)=\max\left(\rho\left(\left|\frac{\partial f}{\partial
v}\left(v^{+}_{j+1/2}(t)\right)\right|\right),\rho\left(\left|\frac{\partial
f}{\partial v}\left(v^{-}_{j+1/2}(t)\right)\right|\right)\right),$ (8)

and $\rho(|\frac{\partial f}{\partial v}(v)|)$ represents the spectral radius
of the matrix $|\frac{\partial f}{\partial v}(v)|$.

While this scheme may look complicated at first glance, a good starting point
to understand it is to set $(v_{x})_{j}=0$. This assumption gives us the first
order KT method, which can be expressed in the dissipative form

$\begin{split}\frac{d}{dt}(v_{j})=&-\frac{1}{2\Delta
x}\left(f(v_{j+1})-f(v_{j-1})\right)+\\\ &{\frac{1}{2\Delta
x}\left(a_{j+1/2}(v_{j+1}-v_{j})-a_{j-1/2}(v_{j}-v_{j-1})\right)}.\end{split}$
(9)

While the first term in the right-hand side of Eq. 9 corresponds to the
approximation of the numerical flux, the second one corresponds to a numerical
dissipation, being similar to a numerical approximation of the second
derivative of $u$. Because $(v_{x})_{j}=0$, we can interpret our numerical
approximation as being piecewise constant in the intervals
$x_{j-1/2},x_{j+1/2}$.

In contrast, the second-order KT method does not set $v_{x}$ to zero and is
piecewise linear in those same intervals. The approximation of $v_{x}$ is
chosen as to preserve the total variation diminishing property (TVD). Recall
that the total-variation (TV) in the case of a numerical solution is defined
as

$TV(v(t))=\sum_{j}|v_{j+1}(t)-v_{j}(t)|,$ (10)

while in the case of an exact solution, the TV is defined as

$TV(u(t))=\int_{\Omega}\left|\frac{\partial u}{\partial{x}}(x,t)\right|dx.$
(11)

The diminishing of the TV quantity is called the TVD property. In the case
that the quantity is differentiable, it implies that $\frac{d}{dt}TV(v)\leq
0$. This property holds for exact solutions of scalar one-dimensional
conservation laws and it ensures that no new extrema are created in the
evolution of the numerical solution, thus preventing spurious oscillations
near shock discontinuitiesLeVeque-1992 . Note that so far we are referring to
the semi-discrete formulation of the KT scheme; to preserve the TVD property
for the fully discretized system we also need to require the use of special
temporal integrators tvd-runge-kutta .

The following choice for $v_{x}$ satisfies the TVD property TVDproperty .

$(v_{x})_{j}=\textrm{minmod}\left(\theta\frac{v_{j}-v_{j-1}}{\Delta
x},\frac{v_{j+1}-v_{j-1}}{2\Delta x},\theta\frac{v_{j+1}-v_{j}}{\Delta
x}\right),$ (12)

with the function minmod defined as

$\textrm{minmod}(x_{1},x_{2},x_{3},...)=\left\\{\begin{matrix}\max(x_{j})&\textrm{if
}x_{j}<0\;\forall j,\\\ \min(x_{j})&\textrm{if }x_{j}>0\;\forall j,\\\
0&\textrm{otherwise,}\end{matrix}\right.$ (13)

and $\theta$ a parameter between $1$ and $2$ KT99 .

For systems of equations, the TVD property is not valid. For instance,
consider the Euler equations and initial data with zero momentum density in a
compact region. As the system starts to evolve, the fluid starts to move in
different directions, and so the total variation of the momentum density
increases. Nevertheless, enforcing this property for the scalar case singles
out a particular one-parameter family of choices of derivatives. That choice
is then enforced even for systems of equations and even in several dimensions,
where, in practice, shows very nice stability properties and stops the
appearance of local oscillations near shocks. These methods converge to
second-order in the $Lip^{\prime}$ semi-norm KT99 , also known as the
Wasserstein distance, which is especially suited for analyzing the convergence
of conservative methods. The convenience of this norm for weak solutions is
very well explained by U. S. Fjordholm and S. Solem in sec-ord-conv .
Therefore, this is the appropriate norm to test the convergence of our
interface method.

##  3 Kurganov-Tadmor interface generalization

In this section, we explain our interface method. For simplicity, we start
with the one-dimensional case. First, assume that we choose to place our
interface at a point $x_{I}$ of a domain $\Omega$. We label our subdomains to
the left and right of $x_{I}$ with the letters $L$ and $R$ respectively.

We discretize $L$ using $N$ points and $R$ using $M$ points, with
$x_{N-1}^{L}=x_{I}=x_{0}^{R}$, as seen in figure 1; that is, in both grids we
store information of our estimation of $\overline{u}(x_{I},t)$.

Figure 1: We split our domains in two regions $L$ and $R$ around an interface
$x_{I}$, and discretize both regions so that they share a grid point at the
interface. We mark our intermediate points for better understanding of the
scheme. Even though in the figure $\Delta x^{L}=\Delta x^{R}$, this is not a
necessary condition.

In other words, we are covering our one-dimensional domain with two grid
patches that do not overlap each other and share a single interface point
$x_{I}$. For simplicity, given any variable $w$, we will use the notation
$w^{L}_{I-1/2}\doteq w^{L}_{N-3/2}$ and $w^{R}_{I+1/2}\doteq w^{R}_{1/2}$. We
only focus on how to approximate fields in the interface, as the rest of the
grid is solved using the standard second-order Kurganov-Tadmor method
previously explained.

We recall that KT updates the estimated point averages of the solution at each
grid point. We want to estimate the cell average of the solution in the
interval $(x_{I-1/2}^{L},x_{I+1/2}^{R})$. The idea of our interface method is
to use the points of the subgrids $L$ and $R$ to estimate the left and right
averages $v^{L}_{I}(t)$ and $v^{R}_{I}(t)$,

$v^{L}_{I}(t)\approx\frac{1}{\Delta
x^{L}}\int_{x^{L}_{I-1/2}}^{x_{I}}u(\xi,t)\;d\xi,\quad
v^{R}_{I}(t)\approx\frac{1}{\Delta
x^{R}}\int_{x_{I}}^{x^{R}_{I+1/2}}u(\xi,t)\;d\xi.$ (14)

and latter communicate both grids by calculating the full average $u_{I}$,

$v_{I}(t)\doteq\frac{\Delta x^{L}v^{L}_{I}(t)+\Delta x^{R}v^{R}_{I}(t)}{\Delta
x^{L}+\Delta x^{R}}\approx\frac{1}{\Delta x^{L}+\Delta
x^{R}}\int_{x^{L}_{I-1/2}}^{x^{R}_{I+1/2}}u(\xi,t)\;d\xi.$ (15)

Our method then proceeds as follows. At points near the interface, we lower
the approximation to first order. Thus we don’t need to compute the derivative
at the interface, as seen in Figure 2. This way we avoid using information
from the adjacent grid, as normally we would need the values $v_{I-1/2}^{L}$
and $v_{I+1/2}^{R}$ to calculate this derivative. An order reduction is in any
case expected for touching grids methods, as happens with SAT methods
CARPENTER1999341 . In SAT methods this is due to the fact that SBP operators
are lower-order near the interfaces.

Figure 2: The interface is marked with a dashed line. For each point, we also
plot a line of slope $v_{x}$. The left and right extremes of those lines
correspond to the values $v^{+}_{j-1/2}$ and $v^{-}_{j+1/2}$, respectively. We
mark in blue and red the averaging areas of $v^{L}_{I}$ and $v^{R}_{I}$.

At this point, we would still need the values $v_{I-1}^{n}$ and $v_{I+1}^{n}$
to calculate the fluxes to both sides of $x_{I}$. To avoid using the value at
the other side of the interface, we instead choose to calculate partial cell
averages on a smaller cell, for instance, for the left to the interface we
consider the cell between $x_{I-1/2}=x_{N-3/2}$ and $x_{I}$, of width $\Delta
x^{L}/2$. Then, to calculate

$(v_{I}^{L})^{n+1}\approx\int_{x^{L}_{I-1/2}}^{x^{I}}u(x,t^{(n+1)})dx,$ (16)

we use the reconstructed fluxes at $x^{L}_{I-1/2}$ and $x_{I}$. The algorithm
for the left grid becomes

$\dfrac{d}{dt}v^{L}_{I}=-\dfrac{G_{I}^{L+}-G_{I}^{L-}}{\Delta x^{L}},$ (17)

with

$\begin{matrix}G_{I}^{L+}&=&f(v^{L+}_{I-1/2})+a^{L}_{I-1/2}v^{L+}_{I-1/2}\\\
&=&f(v^{L+}_{I-1/2})+a^{L}_{I-1/2}v_{I},\quad\;\;\;\;\;\end{matrix}$ (18)

and

$G_{I}^{L-}=f(v^{L-}_{I-1/2})+a^{L}_{I-1/2}v^{L-}_{I-1/2},$ (19)

where the values $v^{L\pm}_{I-1/2}$ and $a^{L}_{I-1/2}$ are calculated using
the definitions (7) and (8). Because $(v_{x})^{L}_{I}=0$ we also have
$v^{L+}_{I-1/2}=v_{I}$.

Similarly, the algorithm for the right grid becomes

$\dfrac{d}{dt}v^{R}_{I}=-\dfrac{G_{I}^{R+}-G_{I}^{R-}}{\Delta x^{R}},$ (20)

with

$G_{I}^{R+}=f(v^{R+}_{I+1/2})-a^{R}_{I+1/2}v^{R+}_{I+1/2}$ (21)

and

$\begin{matrix}G_{I}^{R-}&=&f(v^{R-}_{I+1/2})-a^{R}_{I+1/2}v^{R-}_{I+1/2},\\\
&=&f(v^{R-}_{I+1/2})-a^{R}_{I+1/2}v_{I}.\quad\;\;\;\;\;\end{matrix}$ (22)

Each of these calculations is done with information present on their
respective side of the interface. Thus, they can be performed without any
communication from the other grid. After obtaining both temporal derivatives,
we integrate in time using a TVD temporal scheme to evolve $v^{L}_{I}$ and
$v^{R}_{I}$. Finally, we communicate both grids by calculating a full grid
average between $x_{I-1/2}$ and $x_{I+1/2}$ by averaging $v^{L}_{I}$ and
$v^{R}_{I}$:

$\begin{array}[]{r c l}v_{I}&=&\dfrac{\Delta x^{L}v^{L}_{I}+\Delta
x^{R}v^{R}_{I}}{\Delta x^{L}+\Delta x^{R}}\\\
v^{L}_{I}&\leftarrow&v^{\phantom{R}}_{I}\\\
v^{R}_{I}&\leftarrow&v^{\phantom{L}}_{I}.\end{array}$ (23)

We notice that the propagation information of the waves is only passed at this
very last instance, namely when we average the new values on each side using
the new values from the other. This way we essentially recover the original
first-order Kurganov-Tadmor method.

We can check that the TVD property still holds. From KT99  and TVDproperty  we
know that this property is guaranteed if

$\left|\frac{\Delta x(v_{x})_{j+1/2}}{\Delta v_{j\pm 1/2}}\right|\leq 2,$ (24)

with $\Delta v_{j+1/2}=v_{j+1}-v_{j}$. This is clearly the case at the
interface point, since $(v_{x})_{I}$ = 0. Although we could, in theory,
approximate $(v_{x})^{L}_{I}=\frac{v_{I}-v^{L}_{N-2}}{\Delta x^{L}}$ and
$(v_{x})^{R}_{I}=\frac{v^{R}_{1}-v_{I}}{\Delta x^{R}}$, doing so would not
guarantee that (24) holds.

As a remark, we note that $\dfrac{dv^{L}_{I-1}}{dt}$ and
$\dfrac{dv^{R}_{I+1}}{dt}$ are also slightly modified by our method. For
example, to calculate $\dfrac{dv^{L}_{I-1}}{dt}$, we need to use equation (7),
which in turns uses the value $(v_{x})^{\phantom{R}}_{I}$ which is set to zero
at the interface. This is only a slight modification to the standard second-
order KT method, thus not altering the TVD condition by the same arguments
used before.

As is the case with the full KT scheme, the extension of this method to
equation systems and several dimensions is done by implementing it to each
vector component and each dimension. We simply need to set the gradient in the
normal direction to zero for the interface values, evolve the truncated KT
scheme for each side, and then calculate a final cell average using the
volumes of the cells in each grid.

##  4 Numerical tests

We present convergence results for one-dimensional scalar equations, as well
as results for the two-dimensional Euler equation. In all cases we use the
third order TVD Runge-Kutta scheme introduced in equation $(3.3)$ in tvd-
runge-kutta . We set the parameter $\theta=2.0$ in the minmod definition (12)
in the scalar cases, and set $\theta=1.2$ for the two dimensional Euler
equation (see discussion of optimal $\theta$ values in KT99 ). The CFL
condition is set so that $\Delta t=0.1\Delta x$ for all cases.

###  4.1 Advection

We consider the one-dimensional advection equation with speed 1,

$\frac{\partial}{\partial t}u(x,t)+\frac{\partial}{\partial x}u(x,t)=0,$ (25)

in the domain $x\in[-2,2)$, with periodic boundary conditions. The initial
data is given by a $C^{7}$ peak function modeled with an eighth order
polynomial,

$u_{0}(x)=\left\\{\begin{matrix}(x-1)^{4}(x+1)^{4}&\textrm{if }|x|<1\\\
0&\textrm{if }|x|\geq 1.\end{matrix}\right.$ (26)

The exact solution to this problem is trivial and corresponds to
$u(x,t)=u_{0}(x-t)$. Since the solution is periodic, we also have
$u(x,4n)=u_{0}(x)\;\forall n\in\mathcal{Z}$.

We evolve (25)-(26) until $T=20$, both with and without a numerical interface
placed at $x_{I}=0.5$. For simplicity in the convergence analysis, we assume
$\Delta x^{L}=\Delta x^{R}$ in the interface simulations. The evolution
results are shown in Figure 3(a), and convergence results are shown in 3(b)
for the seminorm $Lip^{\prime}$. We fit the convergence data with a power
equation

$||u-u_{exact}||_{Lip^{\prime}}=b\Delta x^{p}.$ (27)

(a) Evolution of (26) for the advection equation with different number of
points $N$ and with (w i) and without an interface in $x=0$.

(b) $Lip^{\prime}$ convergence of (26) for the advection equation.

Figure 3: Results for the one-dimensional advection problem at T = 20,
corresponding to 4 periods. The CFL is set to $0.1$ in all cases.

We notice that even though there is slightly greater numerical dissipation in
the presence of interfaces, as seen in the deformation of the initial profile
in Figure 3(a), the second-order global convergence of the method is
maintained. The speed of the wave is also correctly preserved, as evidenced by
the approximations still being centered after 5 periods.

Figure 4: Evolution of the total variation for the advection equation with (w
i) and without an interface. The total variation is constant until the shock
is developed.

Figure 4 shows the total variation evolution of the solution of the advection
equation until $T=40$ for several grid resolutions. We notice that the total
variation of the numerical solution diminishes, even though the total
variation of the exact solution remains constant. We also notice that the
total variation diminishes faster in the interface simulations, particularly
so whenever the peak of the function passes through the interface. This is due
to the fact that the total variation only depends on the extrema of the
function, and is therefore in this case only sensitive to the dissipation of
the peak of the numerical solution. It is well known that the KT schemes are
very dissipative at extrema, thus when a maximum or minimum passes through the
interface, where the scheme is only first order, a substantial diminishing of
the TV is expected. We shall see in the Burgers equation example that in the
case of a shock the situation is much better.

###  4.2 Burgers equation

To study the behavior of shocks passing through an interface we consider now
the one-dimensional Burgers equation

$u_{t}+\left(\frac{u^{2}}{2}\right)_{x}=0,$ (28)

with smooth periodic initial data in the interval $[0,2\pi)$,

$u(x,0)=0.5+\sin(x).$ (29)

The solution of (29) is well known and develops a shock at time $T=1$ with
speed $0.5$. Results of the second-order convergence for the KT scheme for
this problem in the case without interfaces can be found in KT99 .

We first evolve the wave until $T=2.0$ so that the shock has time to develop
and travel across the interface, which is placed at $x=1.25\pi$. The
convergence results are shown in Figure 6. We compare our results with a high
resolution approximation using $N=20480$ points without an interface. We
notice that the shock discontinuity successfully travels across the interface
and the second-order global convergence is maintained .

Figure 5: Evolution of (29) with the Burgers equation at different times for
$N=40$ and $CFL=0.1$ with and without an interface. The interface location is
marked by the dashed black line. Figure 6: Convergence of (29) with the
Burgers equation at time $T=2.0$ with CFL = $0.1$. $u_{hres}$ corresponds to
the approximation with $N=20480$.

We also continue the evolution up to $T=114.0$, to allow the shock to travel
through the interface several times. Given that the speed of the shock is 0.5,
this means that it passes through the interface around $9$ times. Second-order
convergence is maintained even after these 9 crosses, as can be seen in
Figures 7(a) and 7(b).

(a) Evolution of (29) with the Burgers equation with different numbers of
points $N$ and with (dashed line) and without (continuous line) interface.
This is after 10 passes of the shock across the interface.

(b) $Lip^{\prime}$ convergence of (29) for the Burgers equation. $u_{hres}$
corresponds to the approximation with $N=20480$.

Figure 7: Results for the one-dimensional Burgers problem at $T=114$ and with
$CFL=0.1$. Figure 8: Logarithmic plot of the evolution of the total variation
for the Burgers equation with (w i) and without an interface. The total
variation is constant until the shock develops.

Figure 8 shows the total variation evolution of the solution of the Burgers
equation until $T=100$ for several grid resolutions. As expected, the total
variation starts to diminish rapidly as soon as the shock develops. In
contrast with the advection equation, the effect of the interface on the total
variation is almost unnoticeable. Notice that for this case the TV plot is
logarithmic.

###  4.3 Two-dimensional Euler equation

We evolve the two-dimensional Euler equations of gas dynamics LeVeque-1992 ,

$\frac{\partial}{\partial t}\begin{pmatrix}\rho\\\ \rho u\\\ \rho w\\\
E\end{pmatrix}+\frac{\partial}{\partial x}\begin{pmatrix}\rho u\\\ \rho
u^{2}+p\\\ \rho uw\\\ u(E+p)\end{pmatrix}+\frac{\partial}{\partial
y}\begin{pmatrix}\rho u\\\ \rho uw\\\ \rho w^{2}+p\\\ u(E+p)\end{pmatrix}=0,$
(30)

where $\rho$ is the density of the gas, $(u,w)$ the 2D fluid velocity, $E$ the
energy density of the system and $p$ the pressure, which obeys state equation

$E=\frac{1}{2}\rho(u^{2}+w^{2})+\frac{p}{\gamma-1},$ (31)

with $\gamma$ the ratio of specific heats. We use $\gamma=1.4$, which
corresponds to an ideal gas.

####  4.3.1 Implosion

We consider the implosion problem from liska2003 . It consists of a box
$(x,y)\in[-0.3,0.3]\times[-0.3,0.3]$ with reflective boundary conditions, with
a smaller box inside rotated by $\pi/4$ with corners in points $[\pm 0.15,0]$
and $[0,\pm 0.15]$. The fluid velocity is set to zero everywhere, and initial
pressure and density are set at $\rho_{i}=0.125$ and $p_{i}=0.14$ in the inner
box, and $\rho_{o}=1,p_{o}=0$ outside of it, as seen in figure 9.

Figure 9: Initial data for the implosion problem.

We take advantage of the symmetry of the problem and only evolve the upper
right quadrant of the full domain, that is, $(x,y)\in(0,0.3)\times(0,0.3)$.

To impose the boundary conditions we add two extra rows of points outside of
each boundary and set their field values to mirror the ones inside, except for
the normal component of the momentum which is opposite to the mirror one
inside the domain.

Figure 10: Implosion problem on a $400\times 400$ grid with $CFL=0.1$ at
$T=2.5$ with and without interfaces. Density contour lines and velocity arrows
are overlaid over the pressure color maps. Figure 11: Difference between the
approximations with and without interfaces for the implosion problem on a
$400\times 400$ grid with $CFL=0.1$ at $T=2.5$.

We evolve this systems on $400\times 400$ grids until $T=2.5$ for three
different cases: one with no interface, one with a vertical interface at
$x=0.15$, and the last one with two interfaces; one vertical at $x=0.15$ and
another horizontal at $y=0.15$. Results of the implementation of the method
can be seen in figure 10. All three plots are very similar, which shows that
the interface method works as intended. To help distinguish between them, we
subtract the pressure approximations with interfaces from the pressure
approximation without any interface and show them in Figure 11.

We find that shocks are correctly propagated through the interfaces even when
there is an angle between the propagation direction and the interface normal,
and there is no indication of bouncing effects. Most of the differences
between approximations are located around shocks. This is in full accordance
with the differences seen in the one-dimensional case, although the $\theta$
value used here is different than that used in the scalar one-dimensional
cases.

####  4.3.2 Gresho vortex

As a final test we evolve the Gresho vortex liska2003  in a box
$(x,y)\in[-1,1)\times[-1,1)$ with periodic boundary conditions:

$u_{\phi}(r)=\left\\{\begin{matrix}[l]5r\\\ 2-5r\\\
0\end{matrix}\right.,\;p(r)=\left\\{\begin{matrix}[l]5+\frac{25}{2}r^{2}&0\leq
r<0.2\\\ 9+\frac{25}{2}r^{2}-20r+4\ln(5r)&0.2\leq r<0.4\\\ 3+4\ln(2)&0.4\leq
r\end{matrix}\right.$ (32)

Figure 12: Gresho problem on a $400\times 400$ grid with and without
interfaces with $CFL=0.1$. The upper row corresponds to $T=3.0$ and the bottom
row to $T=10.0$. The interfaces are placed slightly off center, so as to break
the symmetry of the problem. Density contours are overlaid over the vorticity
color maps.

The exact solution of this problem is independent of time, as the pressure
gradient is constructed in such a way as to exactly cancel the centrifugal
force given by the velocity profile. We solve this problem in a $400\times
400$ grid with and without interfaces in the same way we did in the implosion
problem but placing the interfaces slightly off-center so as to break the
spatial symmetry of the numerical implementation. This symmetry breaking is
realized by the interface points being treated differently than other points
in the grid. This way we slightly excite unstable perturbations. Otherwise, if
the interfaces are at the center of the vortex, the numerical simulation lasts
a much longer time without unleashing the instabilities.

Results seen in Figure 12 show that the vortex is maintained for small times,
but becomes unstable after a long evolution. This indicates that the usual bar
instabilities of this problem are triggered much earlier due to the interfaces
of this method. To check that we perturbed the initial data with a small
azimuthal function and got similar effects. Figure 13 shows the results of
that perturbation, corroborating that the instability is both physical and
generic.

Figure 13: Gresho problem on a $400\times 400$ grid with initial velocity
perturbed by multiplying it by $(1.0+0.01\cos(\phi))$.

##  5 Conclusions:

We introduced a new simple interface method for dealing with shock propagation
in a multi-block setting. This method preserves the TVD condition while being
second-order convergent. It has minimal communication among multiblocks,
allowing for simple and efficient code parallelization. This method is based
on the Kurganov-Tadmor semidiscrete scheme introduced in KT99 , and works by
lowering the local order of convergence in the interface region to first
order. As expected, the global convergence remains of second-order in the
$Lip^{\prime}$ norm, as was shown to be the case in all one-dimensional
numerical experiments. We also applied our method to the 2D Euler equations
for two different situations: an implosion problem and a Gresho vortex. We
found that the scheme allows the transmission of shock discontinuities without
loss of phase or change in speed regardless of the orientation of the shock
with respect to the interface normal. Finally, we also found that the bar
instabilities of the Gresho vortex problem are triggered much earlier when
interfaces not preserving the quadrant symmetry are present.

## Acknowledgments

The authors would like to thank funding from the Consejo Nacional de
Investigaciones Científicas y Técnicas (CONICET), Fondo para la Investigación
Científica y Tecnológica (FONCYT), and Secretaría de Ciencia y Técnica,
Universidad Nacional de Córdoba and discussions with Manuel Tiglio.

## References

  * (1) R. Liska and B. Wendroff, “Comparison of several difference schemes for the euler equations in 1d and 2d,” in Hyperbolic Problems: Theory, Numerics, Applications (T. Y. Hou and E. Tadmor, eds.), (Berlin, Heidelberg), pp. 831–840, Springer Berlin Heidelberg, 2003.
  * (2) M. H. Carpenter, J. Nordström, and D. Gottlieb, “A stable and conservative interface treatment of arbitrary spatial accuracy,” Journal of Computational Physics, vol. 148, no. 2, pp. 341–365, 1999.
  * (3) J. Hesthaven, “Spectral penalty methods,” Applied Numerical Mathematics, vol. 33, no. 1, pp. 23–41, 2000.
  * (4) R. J. LeVeque, Numerical Methods for Conservation Laws.  Birkhäuser Basel, 1992.
  * (5) P. McCorquodale, M. Dorr, J. Hittinger, and P. Colella, “High-order finite-volume methods for hyperbolic conservation laws on mapped multiblock grids,” Journal of Computational Physics, vol. 288, pp. 181–195, 2015\.
  * (6) G. Chesshire and W. D. Henshaw, “A scheme for conservative interpolation on overlapping grids,” SIAM Journal on Scientific Computing, vol. 15, pp. 819–845, July 1994.
  * (7) S. E. Sherer and J. N. Scott, “High-order compact finite-difference methods on general overset grids,” Journal of Computational Physics, vol. 210, no. 2, pp. 459–496, 2005.
  * (8) E. Pärt-Enander and B. Sjögreen, “Conservative and non-conservative interpolation between overlapping grids for finite volume solutions of hyperbolic problems,” Computers & Fluids, vol. 23, no. 3, pp. 551–574, 1994.
  * (9) A. Kurganov and E. Tadmor, “New high-resolution central schemes for nonlinear conservation laws and convection–diffusion equations,” Journal of Computational Physics, vol. 160, pp. 241–282, 05 2000.
  * (10) S. Gottlieb and C.-W. Shu, “Total variation diminishing runge kutta methods,” Mathematics of Computation, vol. 67, no. 221, pp. 73–85, 1998.
  * (11) S. Osher, “Convergence of generalized muscl schemes,” SIAM Journal on Numerical Analysis, vol. 22, no. 5, pp. 947–961, 1985.
  * (12) S. S. U. S. Fjordholm, “Second-order convergence of monotone schemesfor conservation laws,” SIAM J. Numerical Analysis 54(3) (2016) 1920-1945., 2016.