∎

11institutetext: Jayesh Badwaik22institutetext: Dept. of Mathematics, Univ. of
Würzburg, Würzburg, Germany 22email<EMAIL_ADDRESS>33institutetext:
Praveen Chandrashekar44institutetext: TIFR Center for Applicable Mathematics,
Bangalore, India 44email<EMAIL_ADDRESS>55institutetext:
Christian Klingenberg 66institutetext: Dept. of Mathematics, Univ. of
Würzburg, Würzburg, Germany 66email<EMAIL_ADDRESS>
# Single-step Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for
1-D Euler equations

Jayesh Badwaik  Praveen Chandrashekar  Christian Klingenberg

###### Abstract

We propose an explicit, single step discontinuous Galerkin (DG) method on
moving grids using the arbitrary Lagrangian-Eulerian (ALE) approach for one
dimensional Euler equations. The grid is moved with the local fluid velocity
modified by some smoothing, which is found to considerably reduce the
numerical dissipation introduced by Riemann solvers. The scheme preserves
constant states for any mesh motion and we also study its positivity
preservation property. Local grid refinement and coarsening are performed to
maintain the mesh quality and avoid the appearance of very small or large
cells. Second, higher order methods are developed and several test cases are
provided to demonstrate the accuracy of the proposed scheme.

##  1 Introduction

Finite volume schemes based on exact or approximate Riemann solvers are used
for solving hyperbolic conservation laws like the Euler equations governing
compressible flows. These schemes are able to compute discontinuous solutions
in a stable manner since they have implicit dissipation built into them due to
the upwind nature of the schemes. Higher order schemes are constructed
following a reconstruction approach combined with a high order time
integration scheme. Discontinuous Galerkin methods can be considered as higher
order generalizations of finite volume methods which also make use of Riemann
solver technology but do not need a reconstruction step since they evolve a
polynomial solution inside each cell. While these methods are formally high
order accurate on smooth solutions, they can still introduce too much
numerical dissipation in some situations. Springel springel2010pur  gives the
example of a Kelvin-Helmholtz instability in which adding a large constant
velocity to both states leads to suppression of the instability due to
excessive numerical dissipation. This behaviour is attributed to the fact that
fixed grid methods based on upwind schemes are not Galilean invariant. Upwind
schemes, even when they are formally high order accurate, are found to be too
dissipative when applied to turbulent flows Johnsen:2010:AHM:1670833.1670927
since the numerical viscosity can overwhelm the physical viscosity.

For the linear convection equation $u_{t}+au_{x}=0$, the first order upwind
scheme has the modified partial differential equation

$\frac{\partial u}{\partial t}+a\frac{\partial u}{\partial
x}=\frac{1}{2}|a|h(1-\nu)\frac{\partial^{2}u}{\partial
x^{2}}+O(h^{2}),\qquad\nu=\frac{|a|\Delta t}{h}$

which shows that the numerical dissipation is proportional to $|a|$ which is
the wave speed. In case of Euler equations simulated with a Riemann solver,
e.g., the Roe scheme, the wave speeds are related to the eigenvalues of the
flux Jacobian and the numerical dissipation would be proportional to the
absolute values of the eigenvalues, e.g., $|v-c|,|v|,|v+c|$ where $v$ is the
fluid velocity and $c$ is the sound speed. This type of numerical viscosity is
not Galilean invariant since the fluid velocity depends on the coordinate
frame adopted for the description of the flow. Adding a large translational
velocity to the coordinate frame will increase the numerical viscosity and
reduce the accuracy of the numerical solution. Such high numerical viscosity
can be eliminated or minimized if the grid moves along with the flow as in
Lagrangian methods doi:10.1137/0731002 ; Carre:2009:CLH:1552584.1553086 ;
Maire20092391 . However pure Lagrangian methods encounter the issue of large
grid deformations that occur in highly sheared flows as in the Kelvin-
Helmholtz problem requiring some form of re-meshing. A related approach is to
use arbitrary Lagrangian-Eulerian approach HIRT1974227 ; doneaale  where the
mesh velocity can be chosen to be close to the local fluid velocity but may be
regularized to maintain the mesh quality. Even in the ALE approach, it may be
necessary to perform some local remeshing to prevent the grid quality from
degrading. In springel2010pur , the mesh is regenerated after every time step
based on a Delaunay triangulation, which allows it to maintain good mesh
quality even when the fluid undergoes large shear deformation. However these
methods have been restricted to second order accuracy as they rely on
unstructured finite volume schemes on general polygonal/polyhedral cells,
where achieving higher order accuracy is much more difficult compared to
structured grids.

Traditionally, ALE methods have been used for problems involving moving
boundaries as in wing flutter, store separation and other problems involving
fluid structure interaction LOMTEV1999128 ,Venkatasubban19951743 ,wang2015aa ,
Boscheri2017 , Gaburro2018 . In these applications, the main reason to use ALE
is not to minimize the dissipation in upwind schemes but to account for the
moving boundaries, and hence the grid velocities are chosen based on boundary
motion and with a view to maintain good mesh quality. Another class of methods
solve the PDE on moving meshes where the mesh motion is determined based on a
monitor function which is designed to detect regions of large gradients in the
solution, see ammmbook  and the references therein. These methods achieve
automatic clustering of grid points in regions of large gradients. ALE schemes
have been used to compute multi-material flows as in Luo2004304 , since they
are useful to accurately track the material interface. The mesh velocity was
chosen to be equal to the contact speed but away from the material contact,
the velocity was chosen by linear interpolation and was not close to
Lagrangian. There are other methods for choosing the mesh velocity which have
been studied in Boscheri2013 ; Lohner2004 . Lax-Wendroff type ALE schemes for
compressible flows have been developed in Liu20098872 . Finite volume schemes
based on ADER approach have been developed on unstructured grids
Boscheri_Dumbser_2013 , Boscheri201648 , boscheriIJNMF2014 . The theoretical
analysis of ALE-DG schemes in the framework of Runge-Kutta time stepping for
conservation laws has been done in gero .

In the present work, we consider only the one dimensional problem in order to
set down the fundamental principles with which in an upcoming work, we shall
solve the multidimensional problem. The numerical method developed here will
be usable in the multiple dimensions, but additional work is required in
multiple dimensions to maintain a good mesh quality under fluid flow
deformation. We develop an explicit discontinuous Galerkin scheme that is
conservative on moving meshes and automatically satisfies the geometric
conservation law. The scheme is a single step method which is achieved by
using a predictor computed from a Runge-Kutta scheme that is local to each
cell in the sense that it does not require any data from neighbouring cells
and belongs to the class of schemes called ADER method. Due to the single step
nature of the scheme, the TVD limiter has to be applied only once in each time
step unlike in multi-stage Runge-Kutta schemes where the limiter is applied
after each stage update. This nature of the ADER scheme can reduce its
computational expense especially in multi-dimensional problems and while
performing parallel computations. The mesh velocity is specified at each cell
face as the local velocity with some smoothing. We analyze the positivity of
the first order scheme using Rusanov flux and derive a CFL condition. The
scheme is shown to be exact for steady moving contact waves and the solutions
are invariant to the motion of the coordinate frame. Due to Lagrangian nature,
the Roe scheme does not require any entropy fix. However, we identify the
possibility of spurious contact waves arising in some situations. This is due
to the vanishing of the eigenvalue corresponding to the contact wave. While
the cell averages are well predicted, the higher moments of the solution can
be inaccurate. This behaviour of Lagrangian DG schemes does not seem to have
been reported in the literature. We propose a fix for the eigenvalue in the
spirit of the entropy fix of Harten Harten1983  that prevents the spurious
contact waves from occurring in the solution. The methodology developed here
will be extended to multi-dimensional flows in a future work with a view
towards handling complex sheared flows.

The rest of the paper is organized as follows. Section (2) introduces the
Euler equation model that is used in the rest of the paper. In section (3), we
explain the derivation of the scheme on a moving mesh together with the
quadrature approximations and computation of mesh velocity. The computation of
the predicted solution is detailed in section (4). The TVD type limiter is
presented in section (5) for a non-uniform mesh, section (6) shows the
positivity of the first order scheme and section (7) shows the preservation of
constant states. The grid coarsening and refinement strategy is explained in
section (8) while section (9) presents a series of numerical results.

##  2 Euler equations

The Euler equations model the conservation of mass, momentum and energy, and
can be written as a system of coupled partial differential equations laws of
the form

$\frac{\partial\bm{u}}{\partial t}+\frac{\partial\bm{f}(\bm{u})}{\partial
x}=0$ (1)

where $\bm{u}$ is called the vector of conserved variables and
$\bm{f}(\bm{u})$ are the corresponding fluxes given by

$\bm{u}=\begin{bmatrix}\rho\\\ \rho v\\\
E\end{bmatrix},\qquad\bm{f}(\bm{u})=\begin{bmatrix}\rho v\\\ p+\rho v^{2}\\\
\rho Hv\end{bmatrix}$

In the above expressions, $\rho$ is the density, $v$ is the fluid velocity,
$p$ is the pressure and $E$ is the total energy per unit volume, which for an
ideal gas is given by $E=p/(\gamma-1)+\rho v^{2}/2$, with $\gamma>1$ being the
ratio of specific heats at constant pressure and volume, and $H=(E+p)/\rho$ is
the enthalpy. The Euler equations form a hyperbolic system; the flux Jacobian
$A(\bm{u})=\bm{f}^{\prime}(\bm{u})$ has real eigenvalues and linearly
independent eigenvectors. The eigenvalues are $v-c,\ v,\ v+c$ where
$c=\sqrt{\gamma p/\rho}$ is the speed of sound and the corresponding right
eigenvectors are given by

$r_{1}=\begin{bmatrix}1\\\ v-c\\\ H-vc\end{bmatrix},\qquad
r_{2}=\begin{bmatrix}1\\\ v\\\ \frac{1}{2}v^{2}\end{bmatrix},\qquad
r_{3}=\begin{bmatrix}1\\\ v+c\\\ H+vc\end{bmatrix}$ (2)

The hyperbolic property implies that $A$ can be diagonalized as $A=R\Lambda
R^{-1}$ where $R$ is the matrix formed by the right eigenvectors as the
columns and $\Lambda$ is the diagonal matrix of eigenvalues.

##  3 Discontinuous Galerkin method

###  3.1 Mesh and solution space

Figure 1: Example of a space-time cell in the time interval $(t_{n},t_{n+1})$

Consider a partition of the domain into disjoint cells with the $j$’th cell
being denoted by $C_{j}(t)=[x_{j-\frac{1}{2}}(t),x_{j+\frac{1}{2}}(t)]$. As
the notation shows, the cell boundaries are time dependent which means that
the cell is moving in some specified manner. The time levels are denoted by
$t_{n}$ with the time step $\Delta t_{n}=t_{n+1}-t_{n}$. The boundaries of the
cells move with a constant velocity in the time interval $(t_{n},t_{n+1})$
given by

$w_{j+\frac{1}{2}}(t)=w_{j+\frac{1}{2}}^{n},\qquad t_{n}<t<t_{n+1}$

which defines a cell in space-time as shown in figure (1). The algorithm to
choose the mesh velocity $w_{j+\frac{1}{2}}^{n}$ is explained in a later
section. The location of the cell boundaries is given by

$x_{j+\frac{1}{2}}(t)=x_{j+\frac{1}{2}}^{n}+(t-t_{n})w_{j+\frac{1}{2}}^{n},\qquad
t_{n}\leq t\leq t_{n+1}$

Let $x_{j}(t)$ and $h_{j}(t)$ denote the center of the cell $C_{j}(t)$ and its
length, i.e.,

$x_{j}(t)=\frac{1}{2}(x_{j-\frac{1}{2}}(t)+x_{j+\frac{1}{2}}(t)),\qquad
h_{j}(t)=x_{j+\frac{1}{2}}(t)-x_{j-\frac{1}{2}}(t)$

Let $w(x,t)$ be the continuous linear interpolation of the mesh velocity which
is given by

$w(x,t)=\frac{x_{j+\frac{1}{2}}(t)-x}{h_{j}(t)}w_{j-\frac{1}{2}}^{n}+\frac{x-x_{j-\frac{1}{2}}(t)}{h_{j}(t)}w_{j+\frac{1}{2}}^{n},\qquad
x\in C_{j}(t),\quad t\in(t_{n},t_{n+1})$

Figure 2: Example of a discontinuous piecewise polynomial solution

We will approximate the solution of the conservation law by piecewise
polynomials which are allowed to be discontinuous across the cell boundaries
as shown in figure (2). For a given degree $k\geq 0$, the solution in the
$j$’th cell is given by

$\bm{u}_{h}(x,t)=\sum_{m=0}^{k}\bm{u}_{j,m}(t)\varphi_{m}(x,t),\qquad x\in
C_{j}(t)$

where $\\{\bm{u}_{j,m}\in\mathbb{R}^{3},0\leq m\leq k\\}$ are the degrees of
freedom associated with the $j$’th cell. The basis functions $\varphi_{m}$ are
defined in terms of Legendre polynomials

$\varphi_{m}(x,t)=\hat{\varphi}_{m}(\xi)=\sqrt{2m+1}P_{m}(\xi),\qquad\xi=\frac{x-x_{j}(t)}{\frac{1}{2}h_{j}(t)}$

where $P_{m}:[-1,+1]\to\mathbb{R}$ is the Legendre polynomial of degree $m$.
The above definition of the basis functions implies the following
orthogonality property

$\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\varphi_{l}(x,t)\varphi_{m}(x,t)\textrm{d}x=h_{j}(t)\delta_{lm},\qquad
0\leq l,m\leq k$ (3)

We will sometimes also write the solution in the $j$’th cell in terms of the
reference coordinates $\xi$ as

$\bm{u}_{h}(\xi,t)=\sum_{m=0}^{k}\bm{u}_{j,m}(t)\hat{\varphi}_{m}(\xi)$

and we will use the same notation $\bm{u}$ to denote both functions.

###  3.2 Derivation of the scheme

In order to derive the DG scheme on a moving mesh, let us introduce the change
of variable $(x,t)\to(\xi,\tau)$ given by

$\tau=t,\qquad\xi=\frac{x-x_{j}(t)}{\frac{1}{2}h_{j}(t)}$ (4)

For any $0\leq l\leq k$, we now calculate the rate of change of the $l$’th
moment of the solution starting from

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{u}(x,t)\varphi_{l}(x,t)\textrm{d}x$
$\displaystyle=$
$\displaystyle\frac{\textrm{d}}{\textrm{d}\tau}\int_{-1}^{+1}\bm{u}(\xi,\tau)\hat{\varphi}_{l}(\xi)\frac{1}{2}h_{j}(\tau)\textrm{d}\xi$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int_{-1}^{+1}\left[h_{j}(\tau)\frac{\partial\bm{u}}{\partial\tau}+\bm{u}\frac{\textrm{d}h_{j}}{\textrm{d}\tau}\right]\hat{\varphi}(\xi)\textrm{d}\xi$

wherein we used the change of variables given by (4). But we also have the
inverse transform

$t=\tau,\qquad x=x_{j}(\tau)+\frac{\xi}{2}h_{j}(\tau)$

and hence

$\frac{\partial t}{\partial\tau}=1,\qquad\frac{\partial
x}{\partial\tau}=\frac{\textrm{d}x_{j}}{\textrm{d}\tau}+\frac{\xi}{2}\frac{\textrm{d}h_{j}}{\textrm{d}\tau}=\frac{1}{2}(w_{j-\frac{1}{2}}+w_{j+\frac{1}{2}})+\frac{\xi}{2}(w_{j+\frac{1}{2}}-w_{j-\frac{1}{2}})=w(x,t)$

Using the above relations, we can easily show that

$\frac{\partial\bm{u}}{\partial\tau}(\xi,\tau)=\frac{\partial\bm{u}}{\partial
t}(x,t)+w(x,t)\frac{\partial\bm{u}}{\partial x}(x,t)$

Moreover

$\frac{\textrm{d}h_{j}}{\textrm{d}\tau}=w_{j+\frac{1}{2}}-w_{j-\frac{1}{2}}=h_{j}\frac{\partial
w}{\partial x}$

since $w(x,t)$ is linear in $x$ and hence $\partial w/\partial x$ is constant
inside each cell. Hence the $l$’th moment evolves according to

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{u}(x,t)\varphi_{l}(x,t)\textrm{d}x$
$\displaystyle=$
$\displaystyle\int_{-1}^{+1}\left[\frac{\partial\bm{u}}{\partial
t}+w\frac{\partial\bm{u}}{\partial x}+\bm{u}\frac{\partial w}{\partial
x}\right]\hat{\varphi}_{l}(\xi)\frac{1}{2}h_{j}\textrm{d}\xi$ $\displaystyle=$
$\displaystyle\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\left[-\frac{\partial\bm{f}(\bm{u})}{\partial
x}+\frac{\partial}{\partial x}(w\bm{u})\right]\varphi_{l}(x,t)\textrm{d}x$

where we have transformed back to the physical coordinates and made use of the
conservation law (1) to replace the time derivative of the solution with the
flux derivative. Define the ALE flux

$\bm{g}(\bm{u},w)=\bm{f}(\bm{u})-w\bm{u}$ (5)

Performing an integration by parts in the $x$ variable, we obtain

$\displaystyle\frac{\textrm{d}}{\textrm{d}t}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{u}_{h}(x,t)\varphi_{l}(x,t)\textrm{d}x$
$\displaystyle=$
$\displaystyle\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{g}(\bm{u}_{h},w)\frac{\partial}{\partial
x}\varphi_{l}(x,t)\textrm{d}x$
$\displaystyle+\hat{\bm{g}}_{j-\frac{1}{2}}(\bm{u}_{h}(t))\varphi_{l}(x_{j-\frac{1}{2}}^{+},t)-\hat{\bm{g}}_{j+\frac{1}{2}}(\bm{u}_{h}(t))\varphi_{l}(x_{j+\frac{1}{2}}^{-},t)$

where we have introduced the numerical flux

$\hat{\bm{g}}_{j+\frac{1}{2}}(\bm{u}_{h}(t))=\hat{\bm{g}}(\bm{u}_{j+\frac{1}{2}}^{-}(t),\bm{u}_{j+\frac{1}{2}}^{+}(t),w_{j+\frac{1}{2}}(t))$

which provides an approximation to the ALE flux, see Appendix. Integrating
over the time interval $(t_{n},t_{n+1})$ and using (3), we obtain

$\displaystyle h_{j}^{n+1}\bm{u}_{j,l}^{n+1}$ $\displaystyle=$ $\displaystyle
h_{j}^{n}\bm{u}_{j,l}^{n}+\int_{t_{n}}^{t_{n+1}}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{g}(\bm{u}_{h},w)\frac{\partial}{\partial
x}\varphi_{l}(x,t)\textrm{d}x\textrm{d}t$
$\displaystyle+\int_{t_{n}}^{t_{n+1}}[\hat{\bm{g}}_{j-\frac{1}{2}}(t)\varphi_{l}(x_{j-\frac{1}{2}}^{+},t)-\hat{\bm{g}}_{j+\frac{1}{2}}(t)\varphi_{l}(x_{j+\frac{1}{2}}^{-},t)]\textrm{d}t$

The above scheme has an implicit nature since the unknown solution
$\bm{u}_{h}$ appears on the right hand side integrals whereas we only know the
solution at time $t_{n}$. In order to obtain an explicit scheme, we assume
that we have available with us a predicted solution $\bm{U}_{h}$ in the time
interval $(t_{n},t_{n+1})$, which is used in the time integrals to obtain an
explicit scheme. Moreover, the integrals are computed using quadrature in
space and time leading to the fully discrete scheme

$\displaystyle h_{j}^{n+1}\bm{u}_{j,l}^{n+1}$ $\displaystyle=$ $\displaystyle
h_{j}^{n}\bm{u}_{j,l}^{n}$ $\displaystyle+$ $\displaystyle\Delta
t_{n}\sum_{r}\theta_{r}h_{j}(\tau_{r})\sum_{q}\eta_{q}\bm{g}(\bm{U}_{h}(x_{q},\tau_{r}),w(x_{q},\tau_{r}))\frac{\partial}{\partial
x}\varphi_{l}(x_{q},\tau_{r})$ $\displaystyle+$ $\displaystyle\Delta
t_{n}\sum_{r}\theta_{r}[\hat{\bm{g}}_{j-\frac{1}{2}}(\bm{U}_{h}(\tau_{r}))\varphi_{l}(x_{j-\frac{1}{2}}^{+},\tau_{r})-\hat{\bm{g}}_{j+\frac{1}{2}}(\bm{U}_{h}(\tau_{r}))\varphi_{l}(x_{j+\frac{1}{2}}^{-},\tau_{r})]$

where $\theta_{r}$ are weights for time quadrature and $\eta_{q}$ are weights
for spatial quadrature. For the spatial integral, we will use $q=k+1$ point
Gauss quadrature. For the time integral we will use mid-point rule for $k=1$
and two point Gauss quadrature for $k=2,3$. Since the mesh is moving, the
spatial quadrature points $x_{q}$ depend on the quadrature time $\tau_{r}$
though this is not clear from the notation. In practice, the integrals are
computed by mapping the cell to the reference cell, and the basis functions
and its derivatives are also evaluated on the reference cell. The quadrature
points in the reference cell are independent of time due to the linear mesh
evolution.

###  3.3 Mesh velocity

The mesh velocity must be close to the local fluid velocity in order to have a
Lagrangian character to the scheme. Since the solution is discontinuous, there
is no unique fluid velocity at the mesh boundaries. Some researchers,
especially in the context of Lagrangian methods, solve a Riemann problem at
the cell face to determine the face velocity. Since we use an ALE formulation,
we do not require the exact fluid velocity which is anyway not available to
use since we only have a predicted solution. Following the exact trajectory of
the fluid would also lead to curved trajectories for the grid point, which is
an unnecessary complication. In our work, we make two different choices for
the mesh velocities:

  1. 1.

The first choice is to take an average of the two velocities at every face In
the numerical results, we refer to this as ADG.

$\tilde{w}_{j+\frac{1}{2}}^{n}=\frac{1}{2}[v(x_{j+\frac{1}{2}}^{-},t_{n})+v(x_{j+\frac{1}{2}}^{+},t_{n})]$

  2. 2.

The second choice is to solve a linearized Riemann problem at the face at time
$t_{n}$. In the numerical results, we refer to this as RDG. For simplicity of
notation, let the solution to the left of the face $x_{j+\frac{1}{2}}$ be
represented as $u_{j+\frac{1}{2}}^{-}$ and the solution to the right be
represented as $u_{j+\frac{1}{2}}^{+}$, then

$\tilde{w}_{j+\frac{1}{2}}^{n}=\frac{\rho_{j}^{n}c_{j}^{n}v_{j}^{n}+\rho_{j+1}^{n}c_{j+1}^{n}v_{j+1}^{n}}{\rho_{j}^{n}c_{j}^{n}+\rho_{j+1}^{n}c_{j+1}^{n}}+\frac{p_{j}^{n}-p_{j+1}^{n}}{\rho_{j}^{n}c_{j}^{n}+\rho_{j+1}^{n}c_{j+1}^{n}}$

We will also perform some smoothing of the mesh velocity, e.g., the actual
face velocity is computed from

$w_{j+\frac{1}{2}}^{n}=\frac{1}{3}(\tilde{w}_{j-\frac{1}{2}}^{n}+\tilde{w}_{j+\frac{1}{2}}^{n}+\tilde{w}_{j+\frac{3}{2}}^{n})$

Note that our algorithm to choose the mesh velocity is very local and hence
easy and efficient to implement as it does not require the solution of any
global problems. In Springel springel2010pur , the mesh velocity is adjusted
so that the cells remain nearly isotropic which leads to smoothly varying cell
sizes. Such an approach leads to many parameters that need to be selected and
we did not find a good way to make this choice that works well for a range of
problems. Instead, we will make use of mesh refinement and coarsening to
maintain the quality of cells, i.e., to prevent very small or large cells from
occurring in the grid. The use of a DG scheme makes it easy to perform such
local mesh adaptation without loss of accuracy.

###### Remark 1

Consider the application of the proposed ALE-DG scheme to the linear advection
equation $u_{t}+au_{x}=0$. In this case the mesh velocity is equal to the
advection velocity $w_{j+\frac{1}{2}}=a$, i.e., the cells move along the
characteristics. This implies that the ALE flux $g(u,w)=au-wu=0$ and also the
numerical flux $\hat{g}_{j+\frac{1}{2}}=0$. Thus the DG scheme reduces to

$\int_{x_{j-\frac{1}{2}}^{n+1}}^{x_{j+\frac{1}{2}}^{n+1}}u_{h}(x,t_{n+1})\varphi_{l}(x,t_{n+1})\textrm{d}x=\int_{x_{j-\frac{1}{2}}^{n}}^{x_{j+\frac{1}{2}}^{n}}u_{h}(x,t_{n})\varphi_{l}(x,t_{n})\textrm{d}x,\qquad
l=0,1,\ldots,k$

so that the solution at time $t_{n}$ has been advected exactly to the solution
at time $t_{n+1}$. Note that there is no time step restriction involved in
this case and the accuracy of the predicted solution is also not relevant. If
the initial condition has a discontinuity coinciding with a cell face, then
the scheme advects the discontinuity exactly without any diffusion.

##  4 Computing the predictor

The predicted solution is used to approximate the flux integrals over the time
interval $(t_{n},t_{n+1})$ and the method to compute this must be local, i.e.,
it must not require solution from neighbouring cells. Several methods for
computing the predictor have been reviewed in Gassner20114232 . The simplest
approach is to use a Taylor expansion in space and time. Since the cells are
moving, the Taylor expansion has to be performed along the trajectory of the
mesh motion. For a second order scheme, an expansion retaining only linear
terms in $t$ and $x$ is sufficient. Consider a quadrature point
$(x_{q},\tau_{r})$; the Taylor expansion of the solution around the cell
center $x_{j}^{n}$ and time level $t_{n}$ yields

$\displaystyle\bm{u}_{h}(x_{q},\tau_{r})$ $\displaystyle=$
$\displaystyle\bm{u}_{h}(x_{j}^{n},t_{n})+(\tau_{r}-t_{n})\frac{\partial\bm{u}_{h}}{\partial
t}(x_{j}^{n},t_{n})+(x_{q}-x_{j}^{n})\frac{\partial\bm{u}_{h}}{\partial
x}(x_{j}^{n},t_{n})$
$\displaystyle+O(\tau_{r}-t_{n})^{2}+O(x_{q}-x_{j}^{n})^{2}$

and the predicted solution is given by truncating the Taylor expansion at
linear terms, leading to

$\bm{U}(x_{q},\tau_{r})=\bm{u}_{h}(x_{j}^{n},t_{n})+(\tau_{r}-t_{n})\frac{\partial\bm{u}_{h}}{\partial
t}(x_{j}^{n},t_{n})+(x_{q}-x_{j}^{n})\frac{\partial\bm{u}_{h}}{\partial
x}(x_{j}^{n},t_{n})$

Using the conservation law, the time derivative is written as
$\frac{\partial\bm{u}}{\partial t}=-\frac{\partial\bm{f}}{\partial
x}=-A\frac{\partial\bm{u}}{\partial x}$ so that the predictor is given by

$\bm{U}_{h}(x_{q},\tau_{r})=\bm{u}_{h}^{n}(x_{j}^{n})-(\tau_{r}-t_{n})\left[A(\bm{u}_{h}^{n}(x_{j}^{n}))-w_{q}I\right]\frac{\partial\bm{u}_{h}^{n}}{\partial
x}(x_{j}^{n})$ (7)

The above predictor is used for the case of polynomial degree $k=1$. This
procedure can be extended to higher orders by including more terms in the
Taylor expansion but the algebra becomes complicated. Instead we will adopt
the approach of continuous explicit Runge-Kutta (CERK) schemes
Owren:1992:DEC:141072.141096  to approximate the predictor.

Let us choose a set of $(k+1)$ distinct nodes, e.g., Gauss-Legendre or Gauss-
Lobatto nodes, which uniquely define the polynomial of degree $k$. These nodes
are moving with velocity $w(x,t)$, so that the time evolution of the solution
at node $x_{m}$ is governed by

$\displaystyle\frac{\textrm{d}\bm{U}_{m}}{\textrm{d}t}$ $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
t}\bm{U}_{h}(x_{m},t)+w(x_{m},t)\frac{\partial}{\partial
x}\bm{U}_{h}(x_{m},t)$ $\displaystyle=$
$\displaystyle-\frac{\partial}{\partial
x}\bm{f}(\bm{U}_{h}(x_{m},t))+w(x_{m},t)\frac{\partial}{\partial
x}\bm{U}_{h}(x_{m},t)$ $\displaystyle=$
$\displaystyle-[A(\bm{U}_{m}(t))-w_{m}(t)I]\frac{\partial}{\partial
x}\bm{U}_{h}(x_{m},t)=:\bm{K}_{m}(t)$

wherein we have made use of the PDE to write the time derivative in terms of
spatial derivative of the flux. This equation is solved with initial condition

$\bm{U}_{m}(t_{n})=\bm{u}_{h}(x_{m},t_{n})$

Using a Runge-Kutta scheme of sufficient order (see Appendix), we will
approximate the solution at these nodes as

$\bm{U}_{m}(t)=\bm{u}_{h}(x_{m},t_{n})+\sum_{s=1}^{n_{s}}b_{s}((t-t_{n})/\Delta
t_{n})\bm{K}_{m,s},\quad t\in[t_{n},t_{n+1}),\quad m=0,1,\ldots,k$

where $\bm{K}_{m,s}=\bm{K}_{m}(t_{n}+\theta_{s}\Delta t_{n})$,
$\theta_{s}\Delta t_{n}$ is the stage time and $b_{s}$ are certain polynomials
related to the CERK scheme and given in the Appendix. Note that we are
evolving the nodal values but the computation of $\bm{K}_{m,s}$ requires the
modal representation of the solution in order to calculate spatial derivative
of the solution.

Figure 3: Quadrature points for second order scheme Figure 4: Quadrature
points for third order scheme

Once the predictor is computed as above, it must be evaluated at the
quadrature point $(x_{q},\tau_{r})$ as follows. For each time quadrature point
$\tau_{r}\in(t_{n},t_{n+1})$,

  1. 1.

Compute nodal values $\bm{U}_{m}(\tau_{r})$, $m=0,1,\ldots,k$

  2. 2.

For each $r$, convert the nodal values to modal coefficients $\bm{u}_{m,r}$,
$m=0,1,\ldots,k$

  3. 3.

Evaluate predictor
$\bm{U}_{h}(x_{q},\tau_{r})=\sum_{m=0}^{k}\bm{u}_{m,r}\varphi_{m}(x_{q},\tau_{r})$

The conversion from nodal to modal values is accomplished through a
Vandermonde matrix of size $(k+1)\times(k+1)$ which is the same for every cell
and can be inverted once before the iterations start. The predictor is also
computed at the cell boundaries using the above procedure. Figure (3) and (4)
show the quadrature points used in the second, third and fourth order scheme.
For the second order scheme, the values at $\bullet$ and $\Box$ points are
obtained from the predictor based on Taylor expansion as given in equation
(7). For third and fourth order schemes, the nodal values $\times$ are evolved
forward in time by the CERK scheme and evaluated at the $\bullet$ points. The
$\bullet$ point values are converted to modal coefficients using which the
solution at the $\Box$ points is computed.

###### Remark 2

By looping over each cell in the mesh, the predicted solution is computed in
each cell and the cell integral in (3.2) is evaluated. The trace values
$\bm{U}_{j+\frac{1}{2}}^{-}(\tau_{r})$ and
$\bm{U}_{j-\frac{1}{2}}^{+}(\tau_{r})$ at the $\Box$ points needed for
quadrature in time are computed and stored. These are later used in a loop
over the cell faces where the numerical flux is evaluated. Thus the algorithm
is easily parallelizable on multiple core machines and/or using threads.

##  5 Control of Oscillation by Limiting

High order schemes for hyperbolic equations suffer from spurious numerical
oscillations when discontinuities or large gradients are present in the
solution which cannot be accurately resolved on the mesh. In the case of
scalar problems, this is a manifestation of loss of TVD property and hence
limiters are used to satisfy some form of TVD condition. In the case of DG
schemes, the limiter is used as a post processor which is applied on the
solution after the time update has been performed. If the limiter detects that
the solution is oscillatory, then the solution polynomial is reduced to at
most a linear polynomial with a limited slope. In the present scheme, the
limiter is applied after the solution is updated from time $t_{n}$ to time
$t_{n+1}$, i.e., the solution $\bm{u}_{h}^{n+1}$ obtained from (3.2) is post-
processed by the limiter. Since the mesh is inherently non-uniform due to it
being moved with the flow, we modify the standard TVD limiter to account for
this non-uniformity. Also, since we are solving a system of conservation laws,
the limiter is applied on the local characteristic variables which gives
better results than applying it directly on the conserved variables
Cockburn:1989:TRL:69978.69982 .

The solution in cell $j$ can be written as 111We suppress the time variable
for clarity of notation.

$\bm{u}_{h}(x)=\bar{\bm{u}}_{j}+\frac{x-x_{j}}{\frac{1}{2}h_{j}}\bm{s}_{j}+\textrm{higher
order terms}$

where $\bar{\bm{u}}_{j}$ is the cell average value and $\bm{s}_{j}$ is
proportional to the derivative of the solution at the cell center. Let
$R_{j}$, $L_{j}$ denote the matrix of right and left eigenvectors evaluated at
the cell average value $\bar{\bm{u}}_{j}$ which satisfy $L_{j}=R_{j}^{-1}$,
and the right eigenvectors are given in equation (2). The local characteristic
variables are defined by $\bar{\bm{u}}^{*}=L_{j}\bar{\bm{u}}$ and
$\bm{s}^{*}=L_{j}\bm{s}$. We first compute the limited slope of the
characteristic variables from

$\bm{s}_{j}^{**}=h_{j}\
\textrm{minmod}\left(\frac{\bm{s}_{j}^{*}}{h_{j}},\frac{\bar{\bm{u}}^{*}_{j}-\bar{\bm{u}}^{*}_{j-1}}{\frac{1}{2}(h_{j-1}+h_{j})},\frac{\bar{\bm{u}}^{*}_{j+1}-\bar{\bm{u}}^{*}_{j}}{\frac{1}{2}(h_{j}+h_{j+1})}\right)$

where the minmod function is defined by

$\textrm{minmod}\left(a,b,c\right)=\begin{cases}s\min(|a|,|b|,|c|)&\textrm{if
}s=\textrm{sign}(a)=\textrm{sign}(b)=\textrm{sign}(c)\\\
0&\textrm{otherwise}\end{cases}$

If $\bm{s}_{j}^{**}=\bm{s}_{j}^{*}$ then we retain the solution as it is, and
otherwise, the solution is modified to

$\bm{u}_{h}(x)=\bar{\bm{u}}_{j}+\frac{x-x_{j}}{\frac{1}{2}h_{j}}R_{j}\bm{s}_{j}^{**}$

This corresponds to a TVD limiter which is known to lose accuracy at smooth
extrema cockburn1989  since the minmod function returns zero slope at local
extrema. The TVB limiter corresponds to replacing the minmod limiter function
with the following function

$\overline{\textrm{minmod}}\left(a,b,c\right)=\begin{cases}a&\textrm{if
}|a|\leq Mh^{2}\\\
\textrm{minmod}\left(a,b,c\right)&\textrm{otherwise}\end{cases}$ (8)

where the parameter $M$ is an estimate of the second derivative of the
solution at smooth extrema cockburn1989 , and has to be chosen by the user.

##  6 Positivity property

The solutions of Euler equations are well defined only if the density and
pressure are positive quantities. This is not a priori guaranteed by the DG
scheme even when the TVD limiter is applied. In the case of Runge-Kutta DG
schemes, a positivity limiter has been developed in
Zhang:2010:PHO:1864819.1865066  which preserves accuracy in smooth regions.
This scheme is built on a positive first order finite volume scheme. Consider
the first order version of the ALE-DG scheme which is a finite volume scheme
given by

$h_{j}^{n+1}\bar{\bm{u}}_{j}^{n+1}=h_{j}^{n}\bar{\bm{u}}_{j}^{n}-\Delta
t_{n}[\hat{\bm{g}}_{j+\frac{1}{2}}^{n}-\hat{\bm{g}}_{j-\frac{1}{2}}^{n}]$ (9)

The only degree of freedom is the cell average value and the solution is
piecewise constant. We will analyze the positivity of this scheme for the case
of Rusanov flux which is given in A.1. The update equation can be rewritten as

$\displaystyle h_{j}^{n+1}\bar{\bm{u}}_{j}^{n+1}$ $\displaystyle=$
$\displaystyle\left[h_{j}^{n}-\frac{\Delta
t_{n}}{2}(\lambda_{j-\frac{1}{2}}^{n}+w_{j-\frac{1}{2}}^{n}+\lambda_{j+\frac{1}{2}}^{n}-w_{j+\frac{1}{2}}^{n})\right]\bar{\bm{u}}_{j}^{n}$
$\displaystyle+\frac{\Delta
t_{n}}{2}\left[(\lambda_{j-\frac{1}{2}}^{n}-w_{j-\frac{1}{2}}^{n})\bar{\bm{u}}_{j-1}^{n}+\bm{f}_{j-1}^{n}\right]$
$\displaystyle+\frac{\Delta
t_{n}}{2}\left[(\lambda_{j+\frac{1}{2}}^{n}+w_{j+\frac{1}{2}}^{n})\bar{\bm{u}}_{j+1}^{n}-\bm{f}_{j+1}^{n}\right]$
$\displaystyle=$ $\displaystyle a_{j}^{n}\bar{\bm{u}}_{j}^{n}+\frac{\Delta
t_{n}}{2}B_{j}^{n}+\frac{\Delta t_{n}}{2}C_{j}^{n}$

From the definition of the Rusanov flux formula, we can easily see that222We
drop the superscript n in some of these expressions.

$(\lambda_{j-\frac{1}{2}}-w_{j-\frac{1}{2}})+v_{j-1}\geq
c_{j-1}>0,\qquad(\lambda_{j+\frac{1}{2}}+w_{j+\frac{1}{2}})-v_{j+1}\geq
c_{j+1}>0$

Consider the first component of $B_{j}^{n}$

$\displaystyle(\lambda_{j-\frac{1}{2}}-w_{j-\frac{1}{2}})\rho_{j-1}+\rho_{j-1}v_{j-1}$
$\displaystyle\geq$
$\displaystyle(|v_{j-1}-w_{j-\frac{1}{2}}|+c_{j-1}-w_{j-\frac{1}{2}}+v_{j-1})\rho_{j-1}$
$\displaystyle\geq$ $\displaystyle c_{j-1}\rho_{j-1}>0$

Consider the first component of $C_{j}^{n}$

$\displaystyle(\lambda_{j+\frac{1}{2}}+w_{j+\frac{1}{2}})\rho_{j+1}-\rho_{j+1}v_{j+1}$
$\displaystyle\geq$
$\displaystyle(|v_{j}-w_{j+\frac{1}{2}}|+c_{j+1}+w_{j-1}-v_{j+1})\rho_{j+1}$
$\displaystyle\geq$ $\displaystyle c_{j+1}\rho_{j+1}>0$

The pressure corresponding to $B_{j}^{n}$ is

$p_{j-1}\left(-p_{j-1}+\frac{2\rho_{j-1}(v_{j-1}+\lambda_{j-\frac{1}{2}}-w_{j-\frac{1}{2}})^{2}}{\gamma-1}\right)\geq
p_{j-1}\rho_{j-1}c_{j-1}^{2}\frac{\gamma+1}{\gamma-1}\geq 0$

and similarly the pressure corresponding to $C_{j}^{n}$ is non-negative. Hence
if the coefficient in term $a_{j}^{n}$ is positive, then the scheme is
positive. This requires the CFL condition

$\Delta
t_{n}\leq\frac{2h_{j}^{n}}{(\lambda_{j-\frac{1}{2}}^{n}+w_{j-\frac{1}{2}}^{n}+\lambda_{j+\frac{1}{2}}^{n}-w_{j+\frac{1}{2}}^{n})}$

The time step will also be restricted to ensure that the cell size does not
change too much in one time step. If we demand that the cell size does not
change by more than a fraction $\beta\in(0,1)$, then we need to ensure that
the time step satisfies

$\Delta t_{n}\leq\frac{\beta
h_{j}^{n}}{|w_{j+\frac{1}{2}}^{n}-w_{j-\frac{1}{2}}^{n}|}$

Combining the previous two conditions, we obtain the following condition on
the time step

$\Delta t_{n}\leq\Delta
t_{n}^{(1)}:=\min_{j}\left\\{\frac{(1-\tfrac{1}{2}\beta)h_{j}^{n}}{\tfrac{1}{2}(\lambda_{j-\frac{1}{2}}^{n}+\lambda_{j+\frac{1}{2}}^{n})},\frac{\beta
h_{j}^{n}}{|w_{j+\frac{1}{2}}^{n}-w_{j-\frac{1}{2}}^{n}|}\right\\}$ (10)

We can now state the following result on the positivity of the first order
finite volume scheme on moving meshes.

###### Theorem 6.1

The scheme (9) with Rusanov flux is positivity preserving if the time step
condition (10) is satisfied.

###### Remark 3

In this work we have not attempted to prove the positivity of the scheme for
other numerical fluxes. We also do not have a proof of positivity for higher
order version of the scheme. In the computations, we use the positivity
preserving limiter of Zhang:2010:PHO:1864819.1865066  which leads to robust
schemes which preserve the positivity of the cell average value in all the
test cases.

##  7 Preservation of constant states

An important property of schemes on moving meshes is their ability to preserve
constant states for any mesh motion. This is related to the conservation of
cell volumes in relation to the mesh motion. In our scheme, if we start with a
constant state $\bm{u}^{n}_{h}=\bm{c}$, then the predictor is also constant in
the space-time interval, i.e., $\bm{U}_{h}=\bm{c}$. The space-time terms in
(3.2) are polynomials with degree $k+1$ in space and degree one in time and
these are exactly integrated by the chosen quadrature rule. The flux terms at
cell boundaries in (3.2) are of degree one in time and these are also exactly
integrated. Hence the scheme (3.2) can be written as

$\displaystyle h_{j}^{n+1}\bm{u}_{j,l}^{n+1}$ $\displaystyle=$ $\displaystyle
h_{j}^{n}\bm{u}_{j,l}^{n}+\int_{t_{n}}^{t_{n+1}}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{g}(\bm{c},w)\frac{\partial}{\partial
x}\varphi_{l}(x,t)\textrm{d}x\textrm{d}t$
$\displaystyle+\int_{t_{n}}^{t_{n+1}}[\hat{\bm{g}}_{j-\frac{1}{2}}(t)\varphi_{l}(x_{j-\frac{1}{2}}^{+},t)-\hat{\bm{g}}_{j+\frac{1}{2}}(t)\varphi_{l}(x_{j+\frac{1}{2}}^{-},t)]\textrm{d}t$

where $\bm{u}_{j,0}^{n}=\bm{c}$ and $\bm{u}_{j,l}^{n}=0$ for $l=1,2,\ldots,k$.
Due to the constant predictor and by consistency of the numerical flux

$\hat{\bm{g}}_{j+\frac{1}{2}}(t)=\bm{f}(\bm{c})-w_{j+\frac{1}{2}}^{n}\bm{c}$

Moreover, for $l=1,2,\ldots,k$

$\displaystyle\int_{t_{n}}^{t_{n+1}}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\bm{g}(\bm{c},w)\frac{\partial}{\partial
x}\varphi_{l}(x,t)\textrm{d}x\textrm{d}t+\int_{t_{n}}^{t_{n+1}}[\hat{\bm{g}}_{j-\frac{1}{2}}(t)\varphi_{l}(x_{j-\frac{1}{2}}^{+},t)-\hat{\bm{g}}_{j+\frac{1}{2}}(t)\varphi_{l}(x_{j+\frac{1}{2}}^{-},t)]\textrm{d}t$
$\displaystyle=$
$\displaystyle\int_{t_{n}}^{t_{n+1}}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\frac{\partial}{\partial
x}\bm{g}(\bm{c},w)\varphi_{l}(x,t)\textrm{d}x\textrm{d}t+\int_{t_{n}}^{t_{n+1}}[\hat{\bm{g}}_{j-\frac{1}{2}}(t)\varphi_{l}(x_{j-\frac{1}{2}}^{+},t)-\hat{\bm{g}}_{j+\frac{1}{2}}(t)\varphi_{l}(x_{j+\frac{1}{2}}^{-},t)]\textrm{d}t$
$\displaystyle-\int_{t_{n}}^{t_{n+1}}\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\varphi_{l}(x,t)\frac{\partial}{\partial
x}\bm{g}(\bm{c},w)\textrm{d}x\textrm{d}t$ $\displaystyle=$
$\displaystyle-\int_{t_{n}}^{t_{n+1}}\frac{\partial}{\partial
x}\bm{g}(\bm{c},w)\left(\int_{x_{j-\frac{1}{2}}(t)}^{x_{j+\frac{1}{2}}(t)}\varphi_{l}(x,t)\textrm{d}x\right)\textrm{d}t=0$

where we have used the property that $w$ is an affine function of $x$ and
$\varphi_{l}$ are orthogonal. This implies that $\bm{u}_{j,l}^{n+1}=0$ for
$l=1,2,\ldots,k$. For $l=0$, we get

$\displaystyle h_{j}^{n+1}\bm{u}_{j,0}^{n+1}$ $\displaystyle=$ $\displaystyle
h_{j}^{n}\bm{c}+\int_{t_{n}}^{t_{n+1}}[\hat{\bm{g}}_{j-\frac{1}{2}}(t)-\hat{\bm{g}}_{j+\frac{1}{2}}(t)]\textrm{d}t=[h_{j}^{n}+(w_{j-\frac{1}{2}}^{n}-w_{j+\frac{1}{2}}^{n})\Delta
t]\bm{c}$

and since
$h_{j}^{n+1}=h_{j}^{n}+(w_{j-\frac{1}{2}}^{n}-w_{j+\frac{1}{2}}^{n})\Delta t$,
we obtain $\bm{u}_{j,0}=\bm{c}$ which implies that $\bm{u}_{h}^{n+1}=\bm{c}$.

##  8 Grid coarsening and refinement

The size of the cells can change considerably during the time evolution
process due to the near Lagrangian movement of the cell boundaries. Near
shocks, the cells will be compressed to smaller sizes which will reduce the
allowable time step since a CFL condition has to be satisfied. In some
regions, e.g., inside expansion fans, the cell size can increase considerably
which may lead to loss of accuracy. In order to avoid too small or too large
cells from occurring in the grid, we implement cell merging and refinement
into our scheme. If a cell becomes smaller than some specified size $h_{min}$,
then it is merged with one of its neighbouring cells and the solution is
transferred from the two cells to the new cell by performing an $L^{2}$
projection. If a cell size becomes larger than some specified size $h_{max}$,
then this cell is refined into two cells by division and the solution is again
transferred by $L^{2}$ projection. The use of $L^{2}$ projection for solution
transfer ensures the conservation of mass, momentum and energy and preserves
the accuracy in smooth regions. We also ensure that the cell sizes do not
change drastically between neighbouring cells. To keep a track of refinement
of cells, each cell is assigned an initial level equal to $0$. The daughter
cells created during refinement are assigned a level incremented from the
parent cell, while the coarsened cells are assigned a level decremented from
the parent cell.

The algorithm for refinement and coarsening is carried out in three sweeps
over all the active cells. In the first sweep, we mark the cells for
refinement or coarsening based on their size and the level of neighboring
cells. Cells are marked for coarsening if the size is less than a pre-
specified minimum size. They are marked for refinement if either the size of
the cell is larger than the maximum size or if the level of the cell is less
than the level of the neighboring cells. If none of the conditions are
satisfied, the cells are marked for no change. In the second sweep, a cell is
marked for refinement if both the neighboring cells are marked for refinement.
A cell is also marked for refinement if the size of the cell is larger than
twice the size of either of the neighboring cells, and is also larger than
twice the minimum size. The last condition is inserted in order to prevent a
cell being alternately marked for refinement and coarsening in consecutive
adaptation cycles. In the third and final sweep, we again mark cells for
refinement if both the neighboring cells are marked for refinement. Further,
we ensure that a cell marked for refinement does not have a neighboring cell
marked for a coarsening, since this can lead to an inconsistent mesh.

##  9 Numerical results

The numerical tests are performed with polynomials of degree one, two and
three, together with the linear Taylor expansion, two stage CERK and four
stage CERK, respectively, for the computation of the predictor. For the
quadrature in time, we use the mid-point rule, two and three point Gauss-
Legendre quadrature, respectively. The time step is chosen using the CFL
condition,

$\Delta t_{n}=\frac{\textrm{CFL}}{2k+1}\Delta t_{n}^{(1)}$

where $\Delta t_{n}^{(1)}$ is given by equation (10), and the factor $(2k+1)$
comes from linear stability analysis cockburn1989 . In most of the
computations we use $\textrm{CFL}=0.9$ unless stated otherwise. We observe
that the results using average or linearized Riemann velocity are quite
similar. We use the average velocity for most of the results and show the
comparison between the two velocities for some results. The main steps in the
algorithm within one time step $t_{n}\to t_{n+1}$ are as follows.

  1. 1.

Choose mesh velocity $w_{j+\frac{1}{2}}$

  2. 2.

Choose time step $\Delta t_{n}$

  3. 3.

Compute the predictor $\bm{U}_{h}$

  4. 4.

Update solution $\bm{u}_{h}^{n}$ to the next time level $\bm{u}_{h}^{n+1}$

  5. 5.

Apply TVD/TVB limiter on $\bm{u}_{h}^{n+1}$

  6. 6.

Apply positivity limiter on $\bm{u}_{h}^{n+1}$ from
Zhang:2010:PHO:1864819.1865066 .

  7. 7.

Perform grid refinement/coarsening

In all the solution plots given below, symbols denote the cell average value.

###  9.1 Order of accuracy

We study the convergence rate of the schemes by applying them to a problem
with a known smooth solution. The initial condition is taken as

$\rho(x,0)=1+\exp(-10x^{2}),\qquad u(x,0)=1,\qquad p(x,0)=1$

whose exact solution is $\rho(x,t)=\rho(x-t,0)$, $u(x,t)=1$, $p(x,t)=1$. The
initial domain is $[-5,+5]$ and the final time is $t=1$ units. The results are
presented using Rusanov and HLLC numerical fluxes. The $L^{2}$ norm of the
error in density are shown in table (1), (3) for the static mesh and in table
(2), (4) for the moving mesh. In each case, we see that the error behaves as
$O(h^{k+1})$ which is the optimal rate we can expect for smooth solutions. In
table (5), we show that the ALE DG methods preserves its higher order in
presence of a limiter.

The mesh velocity is constant since the fluid velocity is constant. In order
to study the effect of perturbations in mesh velocity, we add a random
perturbation to each mesh velocity, $w_{j+\frac{1}{2}}\leftarrow(1+\alpha
r_{j+\frac{1}{2}})w_{j+\frac{1}{2}}$ where $r_{j+\frac{1}{2}}$ is a uniform
random variable in $[-1,+1]$ and $\alpha=0.05$ and a sample velocity
distribution is shown in figure (5). Note that this randomization is performed
in each time step with different random variables drawn for each face. For
moving mesh, there is no unique cell size, and the convergence rate is
computed based on initial mesh spacing which is inversely proportional to the
number of cells. From table (6) which shows results using HLLC flux, we again
observe that the error reduces at the optional rate of $k+1$ even when the
mesh velocity is not very smooth.

$N$ | $k=1$ | $k=2$ | $k=3$
---|---|---|---
Error | Rate | Error | Rate | Error | Rate
100 | 4.370E-02 | - | 3.498E-03 | - | 3.883E-04 | -
200 | 6.611E-03 | 2.725 | 4.766E-04 | 2.876 | 1.620E-05 | 4.583
400 | 1.332E-03 | 2.518 | 6.415E-05 | 2.885 | 9.376E-07 | 4.347
800 | 3.151E-04 | 2.372 | 8.246E-06 | 2.910 | 5.763E-08 | 4.239
1600 | 7.846E-05 | 2.280 | 1.031E-06 | 2.932 | 3.595E-09 | 4.180
Table 1: Order of accuracy study on static mesh using Rusanov flux $N$ | $k=1$ | $k=2$ | $k=3$
---|---|---|---
Error | Rate | Error | Rate | Error | Rate
100 | 2.331E-02 | - | 3.979E-03 | - | 8.633E-04 | -
200 | 6.139E-03 | 1.9250 | 4.058E-04 | 3.294 | 1.185E-05 | 6.186
400 | 1.406E-03 | 2.0258 | 5.250E-05 | 3.122 | 7.079E-07 | 5.126
800 | 3.375E-04 | 2.0366 | 6.626E-06 | 3.077 | 4.340E-08 | 4.760
1600 | 8.278E-05 | 2.0344 | 8.304E-07 | 3.057 | 2.689E-09 | 4.573
Table 2: Order of accuracy study on moving mesh using Rusanov flux $N$ | $k=1$ | $k=2$ | $k=3$
---|---|---|---
Error | Rate | Error | Rate | Error | Rate
100 | 4.582E-02 |  | 3.952E-03 |  | 3.464E-04 |
200 | 9.611E-03 | 2.253 | 4.048E-04 | 3.287 | 2.058E-05 | 4.073
400 | 2.052E-03 | 2.240 | 4.640E-05 | 3.206 | 1.287E-06 | 4.036
800 | 4.803E-04 | 2.192 | 5.623E-06 | 3.152 | 8.061E-08 | 4.023
1600 | 1.184E-04 | 2.149 | 6.929E-07 | 3.119 | 5.050E-09 | 4.016
Table 3: Order of accuracy study on static mesh using HLLC flux $N$ | $k=1$ | $k=2$ | $k=3$
---|---|---|---
Error | Rate | Error | Rate | Error | Rate
100 | 1.590E-02 |  | 1.626E-03 |  | 1.962E-04 |
200 | 4.042E-03 | 1.977 | 2.072E-04 | 2.972 | 1.269E-05 | 3.950
400 | 1.014E-03 | 1.985 | 2.605E-05 | 2.982 | 7.983E-07 | 3.971
800 | 2.538E-04 | 1.990 | 3.261E-06 | 2.988 | 4.997E-08 | 3.980
1600 | 6.349E-05 | 1.992 | 4.077E-07 | 2.991 | 3.124E-09 | 3.985
Table 4: Order of accuracy study on moving mesh using HLLC flux $N$ | $k=1$ | $k=2$
---|---|---
Error | Rate | Error | Rate
100 | 2.053E-02 | - | 2.277E-03 | -
200 | 4.312E-03 | 2.251 | 3.425E-04 | 2.732
400 | 1.031E-03 | 2.064 | 4.565E-05 | 2.907
800 | 2.550E-04 | 2.015 | 5.812E-06 | 2.973
1600 | 6.356E-05 | 2.004 | 7.315E-07 | 2.990
Table 5: Order of accuracy study on moving mesh using Rusanov flux using Higher Order Limiter Zhong:2013:SWE:2397205.2397446  Figure 5: Example of randomized velocity distribution for smooth test case $N$ | $k=1$ | $k=2$ | $k=3$
---|---|---|---
Error | Rate | Error | Rate | Error | Rate
100 | 1.735E-02 | - | 1.798E-03 | - | 2.351E-04 | -
200 | 4.179E-03 | 2.051 | 2.848E-04 | 2.676 | 1.416E-05 | 4.069
400 | 1.054E-03 | 2.035 | 4.301E-05 | 2.703 | 8.578E-07 | 4.041
800 | 2.615E-04 | 1.943 | 6.012E-06 | 2.838 | 5.476E-08 | 3.958
1600 | 7.279E-05 | 1.852 | 8.000E-07 | 2.909 | 3.505E-09 | 3.966
Table 6: Order of accuracy study on moving mesh using HLLC flux with randomly
perturbed mesh velocity

###  9.2 Smooth Test Case with Non-Constant Velocity

We also test the accuracy of our schemes on a isentropic problem with smooth
solutions. The test case The initial conditions are given by

$\displaystyle\rho(x,0)$ $\displaystyle=1+0.9999995\sin(\pi x)$ $\displaystyle
u(x,0)$ $\displaystyle=0$ $\displaystyle p(x,0)$
$\displaystyle=\rho^{\gamma}(x,0)$ (11)

with $\gamma=3$ and periodic boundary conditions. For this kind of special
isentropic problem, the Euler equations are equivalent to the two Burgers
equations in terms of their two Riemann invariants which can then be used to
derive the analytical solution. The errors are then computed with respect to
the given analytical solution. In contrast to the previous test case, the
velocity and pressure are not constant which makes this a more challenging
test case. We run the simulation with a WENO-type limiter from
Zhang:2010:PHO:1864819.1865066  and positivity limiter enabled. As we can see
from Tables 7, 8, the rate of convergence is maintained for the moving mesh
method with the moving mesh methods exhibiting much lower errors.

$N$ | $k=1$ | $k=2$
---|---|---
Error | Rate | Error | Rate
100 | 8.535E-03 | - | 1.033E-03 | -
200 | 1.958E-03 | 2.124 | 1.221E-04 | 3.08
400 | 4.721E-04 | 2.052 | 1.581E-05 | 2.95
800 | 1.238E-04 | 1.931 | 2.14E-06 | 2.89
1600 | 3.563E-05 | 1.796 | 2.63E-07 | 3.02
Table 7: Order of accuracy study on fixed mesh using Roe flux with Non-Constant Velocity Smooth Test Case $N$ | $k=1$ | $k=2$
---|---|---
Error | Rate | Error | Rate
100 | 4.235E-03 | - | 2.238E-04 | -
200 | 1.058E-03 | 2.001 | 3.255E-05 | 2.87
400 | 2.586E-04 | 2.035 | 4.301E-05 | 3.133
800 | 5.804E-05 | 2.155 | 5.762E-06 | 2.901
1600 | 1.271E-05 | 2.192 | 7.401E-07 | 2.96
Table 8: Order of accuracy study on moving mesh using Roe flux with Non-
Constant Velocity Smooth Test Case

###  9.3 Single contact wave

In this example, we choose a Riemann problem which gives rise to a single
contact wave in the solution that propagates with a constant speed. The
initial condition is given by

$(\rho,v,p)=\begin{cases}(2.0,1.0,1.0)&\textrm{ if }x<0.5\\\
(1.0,1.0,1.0)&\textrm{ if }x>0.5\end{cases}$

and the contact wave moves with a constant speed of 1.0. The solution on
static and moving meshes are shown in figure (6) at time $t=0.5$ using Roe
flux. The moving mesh is able to exactly resolve the contact wave while the
static mesh scheme adds considerable numerical dissipation that smears the
discontinuity over many cells. The accurate resolution of contact waves is a
key advantage of such moving mesh methods, which are capable of giving very
good resolution of the contact discontinuity even on coarse meshes.

|
---|---
(a) | (b)
Figure 6: Single contact wave using Roe flux and 100 cells: (a) static mesh,
(b) moving mesh

###  9.4 Sod problem

The initial condition for the Sod test case is given by sod1978

$(\rho,v,p)=\begin{cases}(1.0,0.0,1.0)&\textrm{ if }x<0.5\\\
(0.125,0.0,0.1)&\textrm{ if }x>0.5\end{cases}$

and the solution is computed upto a final time of $T=0.2$ with the
computational domain being $[0,1]$. Since the fluid velocity is zero at the
boundary, the computational domain does not change with time for the chosen
final time. The exact solution consists of a rarefaction fan, a contact wave
and a shock wave. In figure (7), we show the results obtained using Roe flux
with 100 cells and TVD limiter on static and moving mesh. The contact wave is
considerably well resolved on the moving mesh as compared to the static mesh
due to reduced numerical dissipation on moving meshes.

|
---|---
(a) | (b)
Figure 7: Sod problem using Roe flux, 100 cells and TVD limiter: (a) static mesh (b) moving mesh |
---|---
(a) | (b)
Figure 8: Sod problem using Roe flux, 100 cells and TVD limiter: (a) static
mesh (b) moving mesh

To study the Galilean invariance or the dependence of the solution on the
choice of coordinate frame, we add a boost velocity of $V=10$ or $V=100$ to
the coordinate frame, while implies the initial fluid velocity is $v(x,0)=V$
and the other quantities remain as before. Figure (10a) shows that the
accuracy of the static mesh results degrades with increase in velocity of the
coordinate frame, particularly the contact discontinuity is highly smeared.
The results given in figure (10b) clearly show the independence of the results
on the moving mesh with respect to the coordinate frame velocity. The allowed
time step from CFL condition decreases with increase in coordinate frame speed
for the static mesh case, while in case of moving mesh, it remains invariant.
This means that in case of static mesh, we have to perform more time steps to
reach the same final time as shown in table (9), which increases the
computational time. Thus the moving mesh scheme has the additional advantage
of allowing a larger time step compared to the fixed mesh scheme.

Figure 9: Sod problem using Roe flux, 100 cells and TVD limiter. ADG : Average Velocity, RDG : Linearized Riemann Velocity $V$ | $0$ | $10$ | $100$
---|---|---|---
static mesh | 144 | 810 | 6807
moving mesh | 176 | 176 | 176
Table 9: Number of iterations required to reach time $t=0.2$ for Sod test for different boost velocity of the coordinate frame |
---|---
(a) | (b)
Figure 10: Effect of coordinate frame motion on Sod problem using Roe flux,
100 cells and TVD limiter: (a) static mesh (b) moving mesh

Finally, we compute the solutions using quadratic and cubic polynomials and
the results are shown in figure (11). The solutions look similar to the case
of linear polynomials and have the same sharp resolution of discontinuities.

|
---|---
(a) | (b)
Figure 11: Sod problem on moving mesh using Roe flux, 100 cells and TVD
limiter: (a) Degree = 2 (b) Degree = 3

###  9.5 Lax problem

The initial condition is given by

$(\rho,v,p)=\begin{cases}(0.445,0.698,3.528)&\textrm{if }x<0\\\
(0.5,0.0,0.571)&\textrm{if }x>0\end{cases}$

The computational domain is $[-10,+10]$ and we compute the solution up to a
final time of $T=1.3$. This problem has a strong shock and a contact wave that
is difficult to resolve accurately. The zoomed view of density is shown at the
final time in figure (12), and we observe the moving mesh results are more
accurate for the contact wave, which is the first discontinuity in the figure.
The second discontinuity is a shock which is equally well resolved in both
cases. We can observe that the grid is automatically clustered in the region
between the contact and shock wave, but no explicit grid adaptation was used
in this simulation.

|
---|---
(a) | (b)
Figure 12: Lax problem using HLLC flux, 100 cells and TVD limiter: (a) static
mesh (b) moving mesh Figure 13: Lax problem using HLLC flux, 100 cells and TVD
limiter. ADG : Average Velocity, RDG : Linearized Riemann Velocity

###  9.6 Shu-Osher problem

The initial condition is given by Shu1988439

$(\rho,v,p)=\begin{cases}(3.857143,2.629369,10.333333)&\textrm{if }x<-4\\\
(1+0.2\sin(5x),0.0,1.0)&\textrm{if }x>-4\end{cases}$

which involves a smooth sinusoidal density wave which interacts with a shock.
The domain is $[-5,+5]$ and the solution is computed up to a final time of
$T=1.8$. The solutions are shown in figure (14a)-(14b) on static and moving
meshes using 200 cells and TVD limiter. The moving mesh scheme is considerably
more accurate in resolving the sinusoidal wave structure that arises after
interaction with the shock. In figure (14c) we compute the solution on static
mesh with TVB limiter and the parameter $M=100$ in equation (8). In this case
the solutions on static mesh are more accurate compared to the case of TVD
limiter but still not as good as the moving mesh results. The moving mesh
result has more than 200 cells in the interval $[-5,+5]$ at the final time
since cells enter the domain from the left side. Hence in figure (14d), we
show the static mesh results with 300 cells and using TVB limiter. The results
are further improved for the static mesh case but still not as accurate as the
moving mesh case. The choice of parameters in the TVB limiter is very critical
but we do not have a rigorous algorithm to choose a good value for this. Hence
it is still advantageous to use the moving mesh scheme which gives improved
solutions even with TVD limiter.

|
---|---
(a) | (b)
|
(c) | (d)
Figure 14: Shu-Osher problem using Roe flux: (a) static mesh, 200 cells, $M=0$
(b) moving mesh, 200 cells, $M=0$ (c) static mesh, 200 cells, $M=100$ (d)
static mesh, 300 cells, $M=100$ Figure 15: Shu-Osher problem using Roe flux on
moving mesh

The above results show the ALE method is very accurate in terms of the cell
averages. In figure (15), we show a zoomed view of density and pressure, where
we also plot the linear polynomial solution. The slope of the solution is not
accurately predicted with the Roe scheme and there are spurious contact
discontinuities as the pressure and velocity are nearly continuous. This
behaviour is observed with all contact preserving fluxes like Roe, HLLC and
HLL-CPS but not with the Rusanov flux. Due to the almost Lagrangian character
of the scheme, the eigenvalue corresponding to the contact wave,
$\lambda_{2}=v-w$, is nearly zero, which leads to loss of dissipation in the
corresponding characteristic field. If a spurious contact wave is generated
during the violent dynamics, then this wave will be preserved by the scheme
leading to wrong solutions. We modify the Roe scheme by preventing this
eigenvalue from becoming too small or zero, which is similar to the approach
used for the entropy fix. The eigenvalue $|\lambda_{2}|$ used in the
dissipative part of the Roe flux is determined from

$|\lambda_{2}|=\begin{cases}|v-w|&\textrm{if }|v-w|>\delta=\alpha c\\\
\frac{1}{2}(\delta+|v-w|^{2}/\delta)&\textrm{otherwise}\end{cases}$

With this modification and using $\alpha=0.1$, the solution on moving mesh is
shown in figure (16) and we do not observe the spurious contact
discontinuities which arise with the standard Roe flux, while at the same
time, the solution accuracy compares favourably with the previous results that
did not use the eigenvalue fix.

Figure 16: Shu-Osher problem using modified Roe flux on moving mesh

We next compute the solutions using quadratic polynomials. Figure (17) shows
the results obtained with the TVD limiter which shows the dramatically better
accuracy that is achieved on moving mesh compared to static mesh. In figure
(18) we perform the same computation with a WENO limiter taken from
Zhong:2013:SWE:2397205.2397446 . The static mesh results are now improved over
the case of TVD limiter but still not as good as the moving mesh results in
terms of capturing the extrema. In figure (19) we show a zoomed view of the
results on moving mesh with TVD and WENO limiter. We see that the TVD limiter
is also able to capture all the features and is almost comparable to the WENO
limiter.

|
---|---
(a) | (b)
Figure 17: Shu-Osher problem using modified Roe flux, TVD limiter, quadratic polynomials and 150 cells. (a) static mesh, (b) moving mesh |
---|---
(a) | (b)
Figure 18: Shu-Osher problem using modified Roe flux, WENO limiter, quadratic polynomials and 150 cells. (a) static mesh, (b) moving mesh |
---|---
(a) | (b)
Figure 19: Shu-Osher problem using modified Roe flux, moving mesh, quadratic
polynomials and 150 cells. (a) TVD limiter, (b) WENO limiter

###  9.7 Titarev-Toro problem

Titarev-Toro problem is an extension of the Shu-Osher problem Titarev2014  to
test a severely oscillatory wave interacting with a shock wave. It aims to
test the ability of higher-order methods to capture the extremely high
frequency waves. The initial condition is given by

$\displaystyle(\rho,v,p)=\begin{cases}(1.515695,0.523346,1.805),&-5<x\leq-4.5\\\
(1+0.1\sin(20\pi x),0,1),&-4.5<x\leq 5\end{cases}$ (12)

The computation is carried out on a mesh of 1000 cells with the final time
$T=5$ and the density at this final time is shown in Figures (20), (21). The
fixed mesh is not able to resolve the high frequency oscillations due to
dissipation in the fluxes and the TVD limiter, but the ALE scheme gives an
excellent resolution of these high frequency oscillations. Note that the ALE
scheme also uses the same TVD limiter but it is still able to resolve the
solution to a very degree of accuracy. This result again demonstrates the
superior accuracy that can be achieved by using a nearly Lagrangian ALE scheme
in problems involving interaction of shocks and smooth flow structures.

Figure 20: Titarev Problem with HLLC flux, 1000 cells and TVD limiter Figure
21: Titarev Problem with HLLC flux, 1000 cells and TVD limiter. (Zoomed
Version)

###  9.8 123 problem

The initial condition is given by torobook

$(\rho,v,p)=\begin{cases}(1.0,-2.0,0.4)&x<0.5\\\
(1.0,+2.0,0.4)&x>0.5\end{cases}$

The computational domain is $[0,1]$ and the final time is $T=0.15$. The
density using 100 cells is shown in figure (22) with static and moving meshes.
The mesh motion does not significantly improve the solution compared to the
static mesh case since the solution is smooth. On the contrary, the mesh
becomes rather coarse in the expansion region, though the solution is still
well resolved. However, severe expansion may lead to very coarse meshes which
may be undesirable. To prevent very coarse cells, we switch on the mesh
refinement algorithm as described before and use the upper bound on the mesh
size as $h_{\textrm{max}}=0.05$. The resulting solution is shown in figure
(23) where the number of cells has increased to 108 at the time shown. The
central expansion region is now resolved by more uniformly sized cells
compared to the case of no grid refinement.

|
---|---
(a) | (b)
Figure 22: 123 problem using HLLC flux and 100 cells: (a) static mesh, (b) moving mesh |
---|---
(a) | (b)
Figure 23: 123 problem using HLLC flux and grid refinement: (a) static mesh,
(b) moving mesh with mesh adaptation ($h_{max}=0.05$) leading to 108 cells at
final time Figure 24: 123 problem using HLLC flux, 100 cells and TVD limiter.
ADG : Average Velocity, RDG : Linearized Riemann Velocity

###  9.9 Blast problem

The initial condition is given by Woodward1984115

$(\rho,v,p)=\begin{cases}(1.0,0.0,1000.0)&x<0.1\\\ (1.0,0.0,0.01)&0.1<x<0.9\\\
(1.0,0.0,100.0)&x>0.9\end{cases}$

with a domain of $[0,1]$ and the final time is $T=0.038$. A reflective
boundary condition is used at $x=0$ and $x=1$. A mesh of 400 cells is used for
this simulation and in case of moving mesh, we perform grid adaptation with
$h_{min}=0.001$ since some cells become very small during the collision of the
two shocks. The positivity preserving limiter of
Zhang:2010:PHO:1864819.1865066  is applied together with TVD limiter and HLLC
flux. The static mesh results shown in figure (25a) indicate too much
numerical viscosity in the contact wave around $x=0.6$. This wave is more
accurately resolved in the moving mesh scheme as seen in figure (25b) which is
an advantage due to the ALE scheme and is a very good indicator of the scheme
accuracy as this is a very challenging feature to compute accurately.

|
---|---
(a) | (b)
Figure 25: Blast problem using HLLC flux and 400 cells. (a) static mesh, (b)
moving mesh with adaptation ($h_{min}=0.001$) leading to 303 cells at final
time.

We next compute the same problem using quadratic polynomials with all other
parameters being as before. The solutions are shown in figure (26) and
indicate that the Lagrangian moving mesh scheme is more accurate in resolving
the contact discontinuity. The higher polynomial degree does not show any
major improvement in the solution compared to the linear case, which could be
a consequence of the strong shock interactions present in this problem, see
figure (4.11-4.12) in Zhong:2013:SWE:2397205.2397446  and figure (3.7) in
Zhu2016110  in comparison to current results.

|
---|---
(a) | (b)
Figure 26: Blast problem using HLLC flux, quadratic polynomials and 400 cells.
(a) static mesh, (b) moving mesh with adaptation ($h_{min}=0.001$) leading to
293 cells at final time. Figure 27: Blast problem using HLLC flux, 100 cells
and TVD limiter. ADG : Average Velocity, RDG : Linearized Riemann Velocity

###  9.10 Le Blanc shock tube test case

The Le Blanc shock tube test case is an extreme shock tube problem where the
initial discontinuity separates a region of high energy and density from one
of low energy and density. This is a much more severe test than the Sod
problem and hence more challenging for numerical schemes. The computational
domain is $0\leq x\leq 9$ and is filled with an ideal gas with $\gamma=5/3$.
The gas is initially at rest and we perform the simulation up to a time of
$T=6$ units. The initial discontinuity is at $x=3$ and the initial condition
is given by

$\displaystyle(\rho,v,p)=\begin{cases}(1.0,0.0,0.1)&\text{if }x<3\\\
(0.001,0.0,10^{-7})&\text{if }x>3\end{cases}$ (13)

Note that both the density and pressure have a very large jump in the initial
condition. The solution that develops from this initial condition consists of
a rarefaction wave moving to the left and a contact discontinuity and a strong
shock moving to the right. In Figure (28), we show the comparison of the
internal energy profile at final time between a fixed mesh solution and moving
mesh solutions with two different mesh velocities as described before. Most
methods tend to generate a very large spike in the internal energy in the
contact region, e.g., compare with Figure (11) in Shashkov2005 , while the
present ALE method here is able to give a better profile. We plot the pressure
profile in Figure (29) which shows that the ALE scheme is able to better
represent the region around the contact wave as compared to fixed mesh method.

Figure 28: Internal energy for Le Blanc Shock Tube with Rusanov flux, 1400
cells and TVD limiter, ADG : Average Velocity, RDG : Linearized Riemann
Velocity Figure 29: Pressure for Le Blanc Shock Tube with Rusanov flux, 1400
cells and TVD limiter, ADG : Average Velocity, RDG : Linearized Riemann
Velocity

###  9.11 Two Dimensional Isentropic Vortex Test Case

The extension to two dimensions involves two aspects that need to be
addressed. The first issue is how to handle the grid motion and the second is
how to formulate the ALE-DG scheme. The second part is a natural
generalization of the DG scheme we have described for the 1-D case in this
paper, except that we have to construct basis functions on triangles and
perform some numerical quadrature. The first part involving grid movement is
more complicated and we present only very preliminary results in this section
to demonstrate that the idea has merit. We now consider the two dimensional
Euler equations written as

$\displaystyle\frac{\partial u}{\partial t}+\frac{\partial f(u)}{\partial
x}+\frac{\partial g(u)}{\partial y}$ $\displaystyle=0$ (14)

where

$\displaystyle u$ $\displaystyle=\begin{bmatrix}\rho\\\ \rho u\\\ \rho v\\\
E\end{bmatrix}$ $\displaystyle f(u)$ $\displaystyle=\begin{bmatrix}\rho u\\\
p+\rho u^{2}\\\ \rho uv\\\ (E+p)u\end{bmatrix}$ $\displaystyle g(u)$
$\displaystyle=\begin{bmatrix}\rho v\\\ \rho uv\\\ p+\rho v^{2}\\\
(E+p)v\end{bmatrix}$ (15) $\displaystyle
p=(\gamma-1)\left[E-\frac{1}{2}\rho(u^{2}+v^{2})\right]$ (16)

The test case we consider involves an isentropic vortex that is advecting with
constant velocity and is a smooth solution for which error norms can be
calculated. The test is carried out on a square domain
$[-10,10]\times[-10,10]$ with periodic boundary conditions. The initial
conditions is an isentropic vortex

$\displaystyle T$
$\displaystyle=1-\frac{(\gamma-1)\beta^{2}}{8\gamma\pi^{2}}e^{1-r^{2}}$ (17)
$\displaystyle\rho$ $\displaystyle=T^{\frac{1}{\gamma-1}}$ (18) $\displaystyle
u$ $\displaystyle=u_{\infty}-\frac{\beta}{2\pi}ye^{\frac{1-r^{2}}{2}}$ (19)
$\displaystyle v$
$\displaystyle=v_{\infty}-\frac{\beta}{2\pi}ye^{\frac{1-r^{2}}{2}}$ (20)
$\displaystyle p$ $\displaystyle=\rho^{\gamma}$ (21)

with $u_{\infty}=1,v_{\infty}=0,\gamma=1.4,\beta=10$. As the solution evolves
in time, the mesh becomes quite deformed because the vortex is continually
shearing the mesh, which can lead to degenerate meshes, as shown in figure 30.
We avoid the occurence of badly shaped triangles by using a combination of
face swapping and mesh velocity smoothing algorithms. The mesh modification is
a very local procedure and does not require global remeshing which is a costly
process. With these techniques, we are able to maintain a good mesh quality
even after the vortex has rotated 4 times around its center as shown in
figures 31. As the vortex is translating, we plot the solution in a window
centered at the vortex center. We can see that the method maintains its high
order of accuracy from the convergence rates of the error shown in table 11;
using linear basis functions yields second order convergence while quadratic
basis functions lead to third order convergence.

$N$ | $k=1$ | $k=2$
---|---|---
Error | Rate | Error | Rate
50x50 | 2.230e-03 |  | 1.762E-04 |
100x100 | 5.987E-04 | 1.945 | 2.305E-05 | 2.934
200x200 | 1.498E-04 | 1.998 | 2.973E-06 | 2.955
400x400 | 3.786E-05 | 1.984 | 3.762E-07 | 2.982
800x800 | 9.617E-06 | 1.977 | 3.474E-08 | 2.991
Table 10: Isentropic Vortex in 2D: Order of accuracy study on two dimensional static mesh $N$ | $k=1$ | $k=2$
---|---|---
Error | Rate | Error | Rate
50x50 | 2.230e-03 |  | 1.762E-04 |
100x100 | 5.987E-04 | 1.945 | 2.305E-05 | 2.934
200x200 | 1.498E-04 | 1.998 | 2.973E-06 | 2.955
400x400 | 3.786E-05 | 1.984 | 3.762E-07 | 2.982
800x800 | 9.617E-06 | 1.977 | 3.474E-08 | 2.991
Table 11: Isentropic Vortex in 2D: Order of accuracy study on two dimensional moving mesh Figure 30: Isentropic Vortex in 2-D: Skewed Mesh without Remeshing $t=2.660534$ |
---|---
(a) | (b)
|
(c) | (d)
Figure 31: Isentropic vortex in 2-D: Mesh and pressure solution at various
times (a) $t=0$ (b) $t=6$ (c) $t=12$ (d) $t=20$

##  10 Summary and conclusions

We have developed an explicit DG scheme on moving meshes using ALE framework
and space-time expansion of the solutions within each cell. The near
Lagrangian nature of the mesh motion dramatically reduces the numerical
dissipation especially for contact waves. Even moving contact waves can be
exactly computed with a numerical flux that is exact for stationary contact
waves. The scheme is shown to yield superior results even in the presence of
large boost velocity of the coordinate system indicating its Galilean
invariance property. The standard Roe flux does not suffer from entropy
violation when applied in the current nearly Lagrangian framework. However, in
some problems with strong shocks, spurious contact waves can appear and we
propose to fix the dissipation in Roe-type schemes that eliminates this issue.
The method yields accurate solutions even in combination with standard TVD
limiters, where fixed grid methods perform poorly. The mesh motion provides
automatic grid adaptation near shocks but may lead to very coarse cells inside
expansion waves. A grid adaptation strategy is developed to handle the problem
of very small or very large cells. The presence of the DG polynomials makes it
easy to transfer the solution during grid adaptation without loss of accuracy.
The proposed methodology is general enough to be applicable to other systems
of conservation laws modelling fluid flows. The basic idea can be extended to
multi-dimensions but additional considerations are required to maintain good
mesh quality under fluid deformations. The preliminary results shown for the
isentropic vortex is very promising for the 2-D case.

## Acknowledgements

This work was started when Jayesh Badwaik was a project assistant at TIFR-CAM,
Bangalore. The first two authors gratefully acknowledge the financial support
received from the Airbus Foundation Chair on Mathematics of Complex Systems
established in TIFR-CAM, Bangalore, for carrying out this work. Christian
Klingenberg acknowledges the support of the Priority Program 1648: Software
for Exascale Computing by the German Science Foundation.On behalf of all
authors, the corresponding author states that there is no conflict of
interest.

##  Appendix A Numerical flux

The ALE scheme requires a numerical flux
$\hat{\bm{g}}(\bm{u}_{l},\bm{u}_{r},w)$ which is usually based on some
approximate Riemann solver. The numerical flux function is assumed to be
consistent in the sense that

$\hat{\bm{g}}(\bm{u},\bm{u},w)=\bm{g}(\bm{u},w),\qquad\forall\
\bm{u}\in\mathbb{R}^{3},w\in\mathbb{R}$

Since the ALE versions of the numerical fluxes are not so well known, here we
list the formulae used in the present work.

###  A.1 Rusanov flux

The Rusanov flux is a variant of the Lax-Friedrich flux and is given by

$\hat{\bm{g}}(\bm{u}_{l},\bm{u}_{r},w)=\frac{1}{2}[\bm{g}(\bm{u}_{l},w)+\bm{g}(\bm{u}_{r},w)]-\frac{1}{2}\lambda_{lr}(\bm{u}_{r}-\bm{u}_{l})$

where $\lambda_{lr}=\lambda(\bm{u}_{l},\bm{u}_{r},w)$

$\lambda(\bm{u}_{l},\bm{u}_{r},w)=\max\\{|v_{l}-w|+c_{l},|v_{r}-w|+c_{r}\\}$

which is an estimate of the largest wave speed in the Riemann problem. Since
the mesh velocity is close to the fluid velocity, the value of $\lambda$ is
close to the local sound speed. Thus the numerical dissipation is independent
of the velocity scale.

###  A.2 Roe flux

The Roe scheme Roe1981357  is based on a local linearization of the
conservation law and then exactly solving the Riemann problem for the linear
approximation. The flux can be written as

$\hat{\bm{g}}(\bm{u}_{l},\bm{u}_{r},w)=\frac{1}{2}[\bm{g}(\bm{u}_{l},w)+\bm{g}(\bm{u}_{r},w)]-\frac{1}{2}|A_{w}|(\bm{u}_{r}-\bm{u}_{l})$

where the Roe average matrix $A_{w}=A_{w}(\bm{u}_{l},\bm{u}_{r})$ satisfies

$\bm{g}(\bm{u}_{r},w)-\bm{g}(\bm{u}_{l},w)=A_{w}(\bm{u}_{r}-\bm{u}_{l})$

and we define $|A_{w}|=R|\Lambda-wI|R^{-1}$. This matrix is evaluated at the
Roe average state $\bm{u}(\bar{\bm{q}})$,
$\bar{\bm{q}}=\frac{1}{2}(\bm{q}_{l}+\bm{q}_{r})$, where
$\bm{q}=\sqrt{\rho}[1,\ v,\ H]^{\top}$ is the parameter vector introduced by
Roe.

###  A.3 HLLC flux

This is based on a three wave approximate Riemann solver and the particular
ALE version we use can also be found in Luo2004304 . Define the relative
velocity $q=v-w$; then the numerical flux is given by

$\hat{\bm{g}}(\bm{u}_{l},\bm{u}_{r},w)=\begin{cases}\bm{g}(\bm{u}_{l},w)&S_{l}>0\\\
\bm{g}^{*}(\bm{u}_{l}^{*},w)&S_{l}\leq 0<S_{M}\\\
\bm{g}^{*}(\bm{u}_{r}^{*},w)&S_{M}\leq 0\leq S_{r}\\\
\bm{g}(\bm{u}_{r},w)&S_{r}<0\end{cases}$

where the intermediate states are given by

$\bm{u}_{\alpha}^{*}=\frac{1}{S_{\alpha}-S_{M}}\begin{bmatrix}(S_{\alpha}-q_{\alpha})\rho_{\alpha}\\\
(S_{\alpha}-q_{\alpha})(\rho v)_{\alpha}+p^{*}-p_{\alpha}\\\
(S_{\alpha}-q_{\alpha})E_{\alpha}-p_{\alpha}q_{\alpha}+p^{*}S_{M}\end{bmatrix},\qquad\alpha=l,r$

and

$\bm{g}^{*}(\bm{u},w)=S_{M}\bm{u}+\begin{bmatrix}0\\\ p^{*}\\\
(S_{M}+w)p^{*}\end{bmatrix}$

where

$p^{*}=\rho_{l}(q_{l}-S_{l})(q_{l}-S_{M})+p_{l}=\rho_{r}(q_{r}-S_{r})(q_{r}-S_{M})+p_{r}$

which gives $S_{M}$ as

$S_{M}=\frac{\rho_{r}q_{r}(S_{r}-q_{r})-\rho_{l}q_{l}(S_{l}-q_{l})+p_{l}-p_{r}}{\rho_{r}(S_{r}-q_{r})-\rho_{l}(S_{l}-q_{l})}$

The signal velocities are defined as

$S_{l}=\min\\{q_{l}-c_{l},\hat{v}-w-\hat{c}\\},\qquad
S_{r}=\max\\{q_{r}+c_{r},\hat{v}-w+\hat{c}\\}$

where $\hat{v}$, $\hat{c}$ are Roe’s average velocity and speed of sound.

##  Appendix B Continuous expansion Runge-Kutta (CERK) schemes

We use a Runge-Kutta scheme to compute the predicted solution used to compute
all the integrals in the DG scheme. In this section, we list down the CERK
scheme for the following ODE

$\frac{\textrm{d}u}{\textrm{d}t}=f(u,t)$

Given the solution $u^{n}$ at time $t_{n}$, the CERK scheme gives a polynomial
solution in the time interval $[t_{n},t_{n+1})$ of the form

$u(t_{n}+\theta
h)=u^{n}+h\sum_{s=1}^{n_{s}}b_{s}(\theta)k_{s},\qquad\theta\in[0,1]$

where $n_{s}$ is the number of stages and $h$ denotes the time step.

###  B.1 Second order (CERK2)

The number of stages is $n_{s}=2$ and

$b_{1}(\theta)=\theta-\theta^{2}/2,\qquad b_{2}(\theta)=\theta^{2}/2$

and

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle f(u^{n},t_{n})$
$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle f(u^{n}+hk_{1},t_{n}+h)$

###  B.2 Third order (CERK3)

The number of stages is $n_{s}=4$ and

$\displaystyle k_{1}$ $\displaystyle=$ $\displaystyle f(u^{n},t_{n})$
$\displaystyle k_{2}$ $\displaystyle=$ $\displaystyle
f(u^{n}+(12/23)hk_{1},t_{n}+12h/23)$ $\displaystyle k_{3}$ $\displaystyle=$
$\displaystyle f(u^{n}+h((-68/375)k_{1}+(368/375)k_{2}),t_{n}+4h/5)$
$\displaystyle k_{4}$ $\displaystyle=$ $\displaystyle
f(u^{n}+h((31/144)k_{1}+(529/1152)k_{2}+(125/384)k_{3}),t_{n}+h)$

and

$\displaystyle b_{1}(\theta)$ $\displaystyle=$
$\displaystyle(41/72)\theta^{3}-(65/48)\theta^{2}+\theta$ $\displaystyle
b_{2}(\theta)$ $\displaystyle=$
$\displaystyle-(529/576)\theta^{3}+(529/384)\theta^{2}$ $\displaystyle
b_{3}(\theta)$ $\displaystyle=$
$\displaystyle-(125/192)\theta^{3}+(125/128)\theta^{2}$ $\displaystyle
b_{4}(\theta)$ $\displaystyle=$ $\displaystyle\theta^{3}-\theta^{2}$

## References

  * (1) Walter Boscheri.High Order Accurate Direct Arbitrary-Lagrangian-Eulerian ADER-MOOD Finite Volume Schemes for Non-Conservative Hyperbolic Systems with Stiff Source Terms. Communications in Computational Physics, vol. 21, issue 01, pp. 271-312.
  * (2) Olivier, G., and Alauzet, F. (2011). A New Changing-Topology ALE Scheme for Moving Mesh Unsteady Simulations. In 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition (p. 474).
  * (3) Elena Gaburro, Manuel J Castro, Michael Dumbser; Well-balanced Arbitrary-Lagrangian-Eulerian finite volume schemes on moving nonconforming meshes for the Euler equations of gas dynamics with gravity, Monthly Notices of the Royal Astronomical Society, Volume 477, Issue 2, 21 June 2018, Pages 2251–2275, https://doi.org/10.1093/mnras/sty542
  * (4) Boscheri W., Dumbser M., and Balsara D. S. (2014), High‐order ADER‐WENO ALE schemes on unstructured triangular meshes—application of several node solvers to hydrodynamics and magnetohydrodynamics, Int. J. Numer. Meth. Fluids, 76; pages 737–778, doi: 10.1002/fld.3947
  * (5) Hong Luo, Joseph D. Baum, and Rainald Löhner. 2004. On the computation of multi-material flows using ALE formulation. J. Comput. Phys. 194, 1 (February 2004), 304-328. DOI=http://dx.doi.org/10.1016/j.jcp.2003.09.026
  * (6) V. Springel, E pur si muove: Galilean-invariant cosmological hydrodynamical simulations on a moving mesh, Monthly Notices of the Royal Astronomical Society 401 (2) (2010) 791–851.
  * (7) E. Johnsen, J. Larsson, A. V. Bhagatwala, W. H. Cabot, P. Moin, B. J. Olson, P. S. Rawat, S. K. Shankar, B. Sjögreen, H. C. Yee, X. Zhong, S. K. Lele, Assessment of high-resolution methods for numerical simulations of compressible turbulence with shock waves, J. Comput. Phys. 229 (4) (2010) 1213–1237.
  * (8) C. D. Munz, On Godunov-type schemes for Lagrangian gas dynamics, SIAM Journal on Numerical Analysis 31 (1) (1994) 17–42.
  * (9) G. Carré, S. Del Pino, B. Després, E. Labourasse, A cell-centered Lagrangian hydrodynamics scheme on general unstructured meshes in arbitrary dimension, J. Comput. Phys. 228 (14) (2009) 5160–5183.
  * (10) P.-H. Maire, A high-order cell-centered Lagrangian scheme for two-dimensional compressible fluid flows on unstructured meshes, Journal of Computational Physics 228 (7) (2009) 2391 – 2425.
  * (11) C. Hirt, A. Amsden, J. Cook, An arbitrary Lagrangian-Eulerian computing method for all flow speeds, Journal of Computational Physics 14 (3) (1974) 227 – 253.
  * (12) J. Donea, A. Huerta, J.-P. Ponthot, A. Rodríguez-Ferran, Arbitrary Lagrangian–Eulerian Methods, John Wiley & Sons Ltd, 2004.
  * (13) I. Lomtev, R. Kirby, G. Karniadakis, A discontinuous galerkin ale method for compressible viscous flows in moving domains, Journal of Computational Physics 155 (1) (1999) 128 – 159.
  * (14) C. S. Venkatasubban, A new finite element formulation for {ALE} (arbitrary lagrangian eulerian) compressible fluid mechanics, International Journal of Engineering Science 33 (12) (1995) 1743 – 1762.
  * (15) L. Wang, P.-O. Persson, High-order Discontinuous Galerkin Simulations on Moving Domains using ALE Formulations and Local Remeshing and Projections, American Institute of Aeronautics and Astronautics, 2015.
  * (16) H. Weizhang, R. D. Russell, Adaptive Moving Mesh Methods, Applied Mathematical Sciences, Springer, 2011.
  * (17) H. Luo, J. D. Baum, R. Löhner, On the computation of multi-material flows using ALE formulation, Journal of Computational Physics 194 (1) (2004) 304 – 328.
  * (18) W. Liu, J. Cheng, C.-W. Shu, High order conservative Lagrangian schemes with Lax–Wendroff type time discretization for the compressible Euler equations, Journal of Computational Physics 228 (23) (2009) 8872 – 8891.
  * (19) W. Boscheri, M. Dumbser, Arbitrary-lagrangian-eulerian one-step WENO finite volume schemes on unstructured triangular meshes, Communications in Computational Physics 14 (5) (2013) 1174–1206.
  * (20) W. Boscheri, M. Dumbser, High order accurate direct arbitrary-Lagrangian-Eulerian ADER-WENO finite volume schemes on moving curvilinear unstructured meshes, Computers & Fluids 136 (2016) 48 – 66.
  * (21) W. Boscheri, M. Dumbser, D. S. Balsara, High-order ADER-WENO ale schemes on unstructured triangular meshes—application of several node solvers to hydrodynamics and magnetohydrodynamics, International Journal for Numerical Methods in Fluids 76 (2014) 737–778.
  * (22) C. Klingenberg, G. Schnücke, Y. Xia, Arbitrary lagrangian-eulerian discontinuous galerkin method for conservation laws: analysis and application in one dimension, Mathematics of Computation.
  * (23) G. Gassner, M. Dumbser, F. Hindenlang, C.-D. Munz, Explicit one-step time discretizations for discontinuous galerkin and finite volume schemes based on local predictors, Journal of Computational Physics 230 (11) (2011) 4232 – 4247\.
  * (24) B. Owren, M. Zennaro, Derivation of efficient, continuous, explicit runge-kutta methods, SIAM J. Sci. Stat. Comput. 13 (6) (1992) 1488–1501.
  * (25) B. Cockburn, S.-Y. Lin, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws III: One-dimensional systems, J. Comput. Phys. 84 (1989) 90–113.
  * (26) B. Cockburn, C.-W. Shu, TVB Runge-Kutta local projection discontinuous Galerkin finite element method for conservation laws II: General framework, Mathematics of Computation 52 (186) (1989) pp. 411–435.
  * (27) X. Zhang, C.-W. Shu, On positivity-preserving high order discontinuous galerkin schemes for compressible Euler equations on rectangular meshes, J. Comput. Phys. 229 (23) (2010) 8918–8934.
  * (28) G. Sod, A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws, Journal of Computational Physics 27 (1978) 1–31.
  * (29) C.-W. Shu, S. Osher, Efficient implementation of essentially non-oscillatory shock-capturing schemes, Journal of Computational Physics 77 (2) (1988) 439 – 471.
  * (30) X. Zhong, C.-W. Shu, A simple weighted essentially nonoscillatory limiter for runge-kutta discontinuous galerkin methods, J. Comput. Phys. 232 (1) (2013) 397–415.
  * (31) E. F. Toro, Riemann Solvers and Numerical Methods for Fluid Dynamics: A Practical Introduction, Springer, 1999.
  * (32) P. Woodward, P. Colella, The numerical simulation of two-dimensional fluid flow with strong shocks, Journal of Computational Physics 54 (1) (1984) 115 – 173\.
  * (33) J. Zhu, J. Qiu, A new fifth order finite difference {WENO} scheme for solving hyperbolic conservation laws, Journal of Computational Physics 318 (2016) 110 – 121.
  * (34) P. L. Roe, Approximate Riemann solvers, parameter vectors, and difference schemes, Journal of Computational Physics 43 (2) (1981) 357 – 372.
  * (35) A. Harten and J. M. Hyman, Self adjusting grid methods for one-dimensional hyperbolic conservation laws, J. Comput. Phys., 50 (1983), 253– 269.
  * (36) V. A. Titarev and E. F. Toro, Finite volume WENO schemes for three-dimensional conservation laws, J. Comput. Phys. 201 (2014) 238-260.
  * (37) Loubere, R., and Shashkov, M. J. (2005). A subcell remapping method on staggered polygonal grids for arbitrary-Lagrangian-Eulerian methods. Journal of Computational Physics, 209(1), 105-138.
  * (38) He, P., and Tang, H. (2012). An Adaptive Moving Mesh Method for Two-Dimensional Relativistic Hydrodynamics. Communications in Computational Physics, 11(1), 114-146.
  * (39) Frederic Alauzet. A changing-topology moving mesh technique for large displacements. Engineering with Computers, Springer Verlag, 2014, 30 (2), pp.175-200.
  * (40) Olivier, G., and Alauzet, F. (2011). A New Changing-Topology ALE Scheme for Moving Mesh Unsteady Simulations. In 49th AIAA Aerospace Sciences Meeting including the New Horizons Forum and Aerospace Exposition (p. 474).