Arbitrary Lagrangian-Eulerian discontinuous Galerkin method for conservation laws on moving simplex meshes

TL;DR

This paper extends the ALE-DG method for conservation laws to multiple dimensions on simplex meshes, ensuring geometric conservation, stability, and high-order accuracy, validated through computational experiments.

Contribution

The paper develops a multi-dimensional ALE-DG method on simplex meshes with proven stability, error estimates, and maximum principle, incorporating a limiter that preserves high-order accuracy.

Findings

The method satisfies the geometric conservation law.

L2-stability of the semi-discrete method is proven.

Numerical experiments demonstrate stability and accuracy.

Abstract

In Klingenberg, Schn\"ucke and Xia (Math. Comp. 86 (2017), 1203-1232) an arbitrary Lagrangian-Eulerian discontinuous Galerkin (ALE-DG) method to solve conservation laws has been developed and analyzed. In this paper, the ALE-DG method will be extended to several dimensions. The method will be designed for simplex meshes. This will ensure that the method satisfies the geometric conservation law, if the accuracy of the time integrator is not less than the value of the spatial dimension. For the semi-discrete method the L2-stability will be proven. Furthermore, an error estimate which provides the suboptimal (k+1/2) convergence with respect to the L-infinity-norm will be presented, when an arbitrary monotone flux is used and for each cell the approximating functions are given by polynomials of degree $k$. The two dimensional fully-discrete explicit method will be combined with the bound preserving limiter developed by Zhang, Xia and Shu in (J. Sci. Comput. 50 (2012), 29-62). This limiter does not affect the high order accuracy of a numerical method. Then, for the ALE-DG method revised by the limiter the validity of a discrete maximum principle will be proven. The numerical stability, robustness and accuracy of the method will be shown by a variety of two dimensional computational experiments on moving triangular meshes.

Tables

Table 1. Table 5.3: L 2 superscript L 2 \mathrm{L}^{2} -errors and the rates of convergence for the linear advection equation ( 5.2 ) at final time t = 1 𝑡 1 t=1 on static (right) and moving (left) triangular meshes with cell size h 0 subscript ℎ 0 h_{0} . The bound preserving limiter is not applied.

			$u-u_h^M$				$u-u_h^S$
	$h_0$	$L^2$ norm	order	$\min \left(1.5-u_h^M\right)$	$\min \left(u_h^M-0.5\right)$	$L^2$ norm	order	$\min \left(1.5-u_h^S\right)$	$\min \left(u_h^S-0.5\right)$
$P^1$	1/2	1.30E-01	–	2.22E-02	6.11E-02	1.11E-01	–	2.39E-02	2.49E-02
	1/4	3.09E-02	2.07	-7.38E-02	-8.71E-03	2.09E-02	2.41	-3.50E-02	-1.89E-02
	1/8	6.77E-03	2.19	0.00E+00	9.60E-05	4.43E-03	2.24	-1.29E-02	-1.11E-02
	1/16	1.59E-03	2.09	-7.27E-03	-4.64E-04	1.04E-03	2.09	-3.47E-03	-3.41E-03
	1/32	3.88E-04	2.03	-1.85E-03	-1.10E-04	2.54E-04	2.04	-9.40E-04	-9.27E-04
$P^2$	1/2	2.30E-02	–	-1.01E-01	-2.77E-02	1.79E-02	–	-7.32E-02	-4.33E-02
	1/4	4.88E-03	2.24	-3.11E-02	-4.04E-03	3.00E-03	2.57	-9.71E-03	-5.66E-03
	1/8	7.64E-04	2.68	-4.24E-03	-1.10E-03	4.09E-04	2.88	-8.10E-04	-1.06E-03
	1/16	1.03E-04	2.88	-4.30E-04	-1.14E-04	5.12E-05	3.00	-1.30E-04	-1.87E-04
	1/32	1.31E-05	2.98	-6.00E-05	-2.80E-05	6.28E-06	3.03	-2.00E-05	-2.50E-05
$P^3$	1/2	4.05E-03	–	-2.60E-03	-9.28E-04	2.01E-03	–	-5.50E-04	-2.20E-03
	1/4	3.12E-04	3.70	-1.00E-05	-1.39E-04	1.30E-04	3.95	-6.00E-05	-1.29E-04
	1/8	1.93E-05	4.02	1.00E-05	-6.00E-06	7.85E-06	4.05	0.00E+00	-9.00E-06
	1/16	1.22E-06	3.98	0.00E+00	-1.00E-06	4.79E-07	4.03	0.00E+00	0.00E+00
	1/32	7.71E-08	3.98	0.00E+00	0.00E+00	2.96E-08	4.02	0.00E+00	0.00E+00

Table 2. Table 5.4: L 2 superscript L 2 \mathrm{L}^{2} -errors and the rates of convergence for the moving mesh ALE-DG solution u ~ h M superscript subscript ~ 𝑢 ℎ 𝑀 \widetilde{u}_{h}^{M} with the bound preserving limiter at final time t = 1 𝑡 1 t=1 for the linear advection equation ( 5.2 ) on moving triangular meshes with cell size h 0 subscript ℎ 0 h_{0} .

	$h_0$	$\Vert u-{\tilde{u}}_h^M\Vert$	order	$\min \left(1.5-{\tilde{u}}_h^M\right)$	$\min \left({\tilde{u}}_h^M-0.5\right)$
$P^1$	1/2	1.36E-01	–	4.17E-02	7.82E-02
	1/4	3.31E-02	2.04	0.00E+00	1.99E-03
	1/8	7.94E-03	2.06	0.00E+00	5.05E-03
	1/16	1.84E-03	2.11	0.00E+00	1.09E-03
	1/32	4.41E-04	2.06	0.00E+00	2.72E-04
$P^2$	1/2	6.26E-02	–	0.00E+00	2.20E-02
	1/4	1.07E-02	2.54	0.00E+00	8.60E-04
	1/8	1.18E-03	3.19	0.00E+00	3.64E-05
	1/16	1.23E-04	3.26	0.00E+00	4.52E-06
	1/32	1.46E-05	3.08	0.00E+00	2.13E-08
$P^3$	1/2	5.96E-03	–	0.00E+00	1.63E-03
	1/4	4.69E-04	3.67	2.43E-04	1.98E-04
	1/8	3.02E-05	3.96	4.66E-05	9.56E-06
	1/16	1.76E-06	4.10	1.47E-06	1.62E-06
	1/32	1.01E-07	4.13	0.00E+00	0.00E+00

Table 3. Table 5.5: L 2 superscript L 2 \mathrm{L}^{2} -errors and rates of convergence for Burgers’ equation ( 5.3 ) at final time t = 0.1 𝑡 0.1 t=0.1 on static (left) and moving (center) triangular meshes with cell size h 0 subscript ℎ 0 h_{0} . On the right the L 2 superscript L 2 \mathrm{L}^{2} -errors, rates of convergence and bounds for the ALE-DG solution revised by the bound preserving limiter.

	$h_0$	$\Vert u-u_h^S\Vert$	order	$\Vert u-u_h^M\Vert$	order	$\Vert u-{\tilde{u}}_h^M\Vert$	order	$\min \left(1.5-{\tilde{u}}_h^M\right)$	$\min \left({\tilde{u}}_h^M-0.5\right)$
$P^1$	1/2	6.15E-02	–	6.21E-02	–	6.18E-02	–	0.00E+00	0.00E+00
	1/4	1.78E-02	1.79	1.65E-02	1.91	1.58E-02	1.97	0.00E+00	0.00E+00
	1/8	4.18E-03	2.09	3.89E-03	2.09	3.87E-03	2.03	0.00E+00	0.00E+00
	1/16	1.02E-03	2.04	9.44E-04	2.04	9.82E-04	1.98	0.00E+00	0.00E+00
	1/32	2.49E-04	2.03	2.31E-04	2.03	2.40E-04	2.03	0.00E+00	0.00E+00
$P^2$	1/2	2.54E-02	–	2.54E-02	–	4.71E-02	–	0.00E+00	0.00E+00
	1/4	4.16E-03	2.61	4.10E-03	2.63	1.23E-02	1.93	0.00E+00	0.00E+00
	1/8	7.02E-04	2.57	6.72E-04	2.61	8.18E-04	3.91	0.00E+00	0.00E+00
	1/16	1.14E-04	2.62	1.08E-04	2.64	1.10E-04	2.90	0.00E+00	0.00E+00
	1/32	1.66E-05	2.78	1.59E-05	2.77	1.59E-05	2.78	0.00E+00	0.00E+00
$P^3$	1/2	7.70E-03	–	7.70E-03	–	1.22E-02	–	0.00E+00	0.00E+00
	1/4	8.82E-04	3.12	9.17E-04	3.07	1.07E-03	3.51	0.00E+00	0.00E+00
	1/8	6.44E-05	3.78	6.15E-05	3.90	6.35E-05	4.08	0.00E+00	0.00E+00
	1/16	4.18E-06	3.95	3.93E-06	3.97	4.02E-06	3.98	6.08E-07	9.90E-09
	1/32	2.72E-07	3.94	2.55E-07	3.95	2.59E-07	3.96	6.36E-08	1.48E-08

Table 4. Table 5.6: L 2 superscript L 2 \mathrm{L}^{2} -errors and rates of convergence at final time t = 1 𝑡 1 t=1 for the Euler plain wave problem on static (right) and moving (left) triangular meshes with cell size h 0 subscript ℎ 0 h_{0} .

		$\rho -{\rho }_h^M$		$\rho -{\rho }_h^S$
	$h_0$	$L^2$-norm	order	$L^2$-norm	order
$P^1$	1/2	1.35E-01	–	1.15E-01	–
	1/4	3.04E-02	2.15	1.81E-02	2.67
	1/8	6.06E-03	2.32	3.49E-03	2.37
	1/16	1.40E-03	2.11	8.02E-04	2.12
	1/32	3.41E-04	2.04	1.93E-04	2.05
$P^2$	1/2	2.64E-02	–	2.23E-02	–
	1/4	6.35E-03	2.06	4.58E-03	2.28
	1/8	1.08E-03	2.56	6.74E-04	2.76
	1/16	1.55E-04	2.79	8.62E-05	2.97
	1/32	2.04E-05	2.93	1.04E-05	3.05
$P^3$	1/2	4.75E-03	–	2.37E-03	–
	1/4	3.44E-04	3.79	1.37E-04	4.11
	1/8	2.02E-05	4.09	8.05E-06	4.09
	1/16	1.34E-06	3.92	4.92E-07	4.03
	1/32	8.78E-08	3.93	3.05E-08	4.01

Table 5. Table 5.7: L 2 superscript L 2 \mathrm{L}^{2} -errors and the rates of convergence for the density ρ 𝜌 \rho and the pressure p 𝑝 p at final time t = 10 2 + 5 2 𝑡 superscript 10 2 superscript 5 2 t=\sqrt{10^{2}+5^{2}} for the Euler vortex problem on static (right) and moving (left) triangular meshes with cell size h 0 subscript ℎ 0 h_{0} .

	$h_0$	$\Vert \rho -{\rho }_h^M\Vert$	order	$\Vert p-p_h^M\Vert$	order	$\Vert \rho -{\rho }_h^S\Vert$	order	$\Vert p-p_h^S\Vert$	order
$P^1$	l/2	1.35E-03	–	1.90E-03	–	1.29E-03	–	1.81E-03	–
	l/4	2.84E-04	2.25	3.98E-04	2.25	2.64E-04	2.29	3.72E-04	2.28
	l/8	5.60E-05	2.34	7.83E-05	2.35	5.19E-05	2.35	7.28E-05	2.35
	l/16	1.22E-05	2.19	1.71E-05	2.19	1.16E-05	2.16	1.63E-05	2.16
	l/32	2.81E-06	2.12	3.94E-06	2.12	2.71E-06	2.10	3.80E-06	2.10
$P^2$	l/2	4.75E-04	–	6.59E-04	–	4.32E-04	–	6.00E-04	–
	l/4	6.49E-05	2.87	9.11E-05	2.85	6.34E-05	2.77	8.89E-05	2.75
	l/8	1.16E-05	2.48	1.63E-05	2.48	9.43E-06	2.75	1.32E-05	2.75
	l/16	1.85E-06	2.65	2.59E-06	2.65	1.27E-06	2.90	1.78E-06	2.90
	l/32	2.92E-07	2.66	4.09E-07	2.66	1.77E-07	2.84	2.48E-07	2.84
$P^3$	l/2	1.26E-04	–	1.75E-04	–	1.16E-04	–	1.61E-04	–
	l/4	6.72E-06	4.23	9.36E-06	4.23	5.50E-06	4.40	7.68E-06	4.39
	l/8	3.24E-07	4.37	4.52E-07	4.37	2.20E-07	4.64	3.07E-07	4.64
	l/16	1.54E-08	4.40	2.14E-08	4.40	1.09E-08	4.34	1.51E-08	4.35
	l/32	8.28E-10	4.22	1.15E-09	4.22	5.63E-10	4.27	7.84E-10	4.27

Table 6. Table 5.8: L 2 superscript L 2 \mathrm{L}^{2} -errors for the advection equation and Burgers’ equation at time t = 1 𝑡 1 t=1 with constant initial condition u 0 = 1 subscript 𝑢 0 1 u_{0}=1 on moving triangular meshes with the grid point distribution ( 5.1 ) and cell size h 0 subscript ℎ 0 h_{0} .

	Advection equation	$u-u_h^M$		Burgers’ equation	$u-u_h^M$
$h_0$	$P^1$	$P^2$	$P^3$	$P^1$	$P^2$	$P^3$
1/2	5.71E-16	3.72E-15	8.65E-15	3.03E-16	2.57E-15	7.35E-15
1/4	7.89E-16	7.42E-15	1.99E-14	5.20E-16	5.93E-15	1.56E-14
1/8	2.27E-15	1.24E-14	3.86E-14	1.13E-15	8.86E-15	2.89E-14
1/16	4.21E-15	2.47E-14	7.88E-14	2.44E-15	1.75E-14	5.86E-14
1/32	9.11E-15	5.39E-14	1.67E-13	5.06E-15	3.56E-14	1.19E-13

Table 7. Table 5.9: L ∞ superscript L \mathrm{L}^{\infty} -errors and L 2 superscript L 2 \mathrm{L}^{2} -errors for the advection equation at t = 1.0 𝑡 1.0 t=1.0 with constant initial condition u 0 = 1 subscript 𝑢 0 1 u_{0}=1 on moving triangular meshes with the grid point distribution in Fig. 5.4 .

	Forward Euler		TVD-RK2		TVD-RK3
$N$	$L^{\infty }$	$L^2$	$L^{\infty }$	$L^2$	$L^{\infty }$	$L^2$
32	1.40E-02	8.79E-03	3.55E-15	1.74E-15	1.47E-14	8.21E-15
128	6.21E-03	3.09E-03	1.78E-15	7.29E-16	1.35E-14	9.11E-15
512	2.72E-03	1.30E-03	2.33E-15	4.65E-17	3.57E-12	1.26E-12
2048	1.34E-03	6.22E-04	1.03E-11	2.21E-12	2.39E-11	3.39E-12
8192	6.67E-04	3.07E-04	3.06E-11	1.00E-11	8.45E-11	1.77E-11

Equations

$${\partial }_{t}u+\nabla \cdot \mathbf{f}\left(u\right)=0,\text{in }\Omega \times \left(0,T\right),$$

$$M:=\underset{\mathbf{x}\in \Omega }{\max }{u}_0\left(\mathbf{x}\right)\text{and}m:=\underset{\mathbf{x}\in \Omega }{\min }{u}_0\left(\mathbf{x}\right).$$

$$\overline{\Omega }=\bigcup \left\{\overline{K(t_n)}\mid K(t_n)\in {\mathcal{T}}_{(t_n)}\right\}.$$

$$K_{\text{ref}}:=\left\{\mathbf{\xi }={\left({\xi }_1,...,{\xi }_d\right)}^T\in {\mathbb{R}}^d:{\xi }_{\nu }\ge 0,\forall \nu ,\text{and}\sum _{\nu =1}^d{\xi }_{\nu }\le 1\right\}.$$

$${\mathbf{v}}_{\ell }\left(t\right):={\mathbf{v}}_{\ell }^n+{\mathbf{\omega }}_{K^n,\ell }\left(t-t_n\right),{\mathbf{\omega }}_{K^n,\ell }:=\frac{1}{△t}\left({\mathbf{v}}_{\ell }^{n+1}-{\mathbf{v}}_{\ell }^n\right).$$

$$\begin{array}{rl}K\left(t\right):=& \text{int}\left(\text{conv}\left\{{\mathbf{v}}_1\left(t\right),\dots ,{\mathbf{v}}_{d+1}\left(t\right)\right\}\right),\partial K\left(t\right)=\bigcup _{\nu =1}^{d+1}{F}_{K\left(t\right)}^{\nu },\\ & F_{K\left(t\right)}^{\nu }:=\text{conv}\left(\left\{{\mathbf{v}}_1\left(t\right),\dots ,{\mathbf{v}}_{d+1}\left(t\right)\right\}\setminus \left\{{\mathbf{v}}_{\nu }\left(t\right)\right\}\right),\end{array}$$

$$h:=\underset{t\in \left[0,T\right]}{\max }\underset{K\left(t\right)\in {\mathcal{T}}_{\left(t\right)}}{\max }{h}_{K\left(t\right)}.$$

$$h_{K\left(t\right)}\le \kappa {\rho }_{K\left(t\right)}\text{and}h\le \tau h_{K\left(t\right)},\forall K\left(t\right)\in {\mathcal{T}}_{\left(t\right)}.$$

$${\mathbf{\chi }}_{K\left(t\right)}:K_{\text{ref}}\to \overline{K\left(t\right)},\mathbf{\xi }\mapsto {\mathbf{\chi }}_{K\left(t\right)}\left(\mathbf{\xi },t\right):={\mathbf{A}}_{K\left(t\right)}\mathbf{\xi }+{\mathbf{v}}_1\left(t\right),$$

$${\mathbf{A}}_{K\left(t\right)}:=\left({\mathbf{v}}_2\left(t\right)-{\mathbf{v}}_1\left(t\right),\dots ,{\mathbf{v}}_{d+1}\left(t\right)-{\mathbf{v}}_1\left(t\right)\right).$$

$$J_{K\left(t\right)}=\text{det}\left({\mathbf{A}}_{K\left(t\right)}\right)=d!\left|K\left(t\right)\right|,$$

$$\nabla \cdot \mathbf{f}={\nabla }_{\mathbf{\xi }}\cdot \left[\left({\mathbf{A}}_{K\left(t\right)}^{-1}\right){\mathbf{f}}^{\ast }\right],\nabla u={\mathbf{A}}_{K\left(t\right)}^{-T}{\nabla }_{\mathbf{\xi }}{u}^{\ast },{\mathbf{n}}_{K(t)}=J_{K\left(t\right)}{\mathbf{A}}_{K\left(t\right)}^{-T}{\mathbf{n}}_{K_{\text{ref}}},$$

$$\bm{\omega}_{K\left(t\right)}\left({\bf x},t\right):=\frac{d}{dt}\Big{(}\bm{\chi}_{K\left(t\right)}\left(\bm{\xi},t\right)\Big{)}.$$

$$\partial_{\xi_{i}}\omega_{j}\left(\bm{\chi}_{K\left(t\right)}\left(\bm{\xi},t\right)\Big{)},t\right)=\left[\frac{d}{dt}\left(\mathrm{{\bf A}}_{K\left(t\right)}\right)\right]_{ji},\qquad i,j=1,\dots,d,$$

$$\begin{array}{rl}\frac{d}{dt}\left(J_{K\left(t\right)}\right)=& \text{tr}\left[\text{adj}\left({\mathbf{A}}_{K\left(t\right)}\right)\frac{d}{dt}\left({\mathbf{A}}_{K\left(t\right)}\right)\right]\\ =& \sum _{i=1}^d\sum _{j=1}^d{\left(-1\right)}^{i+j}{M}_{K\left(t\right)}^{ij}{\left[\frac{d}{dt}\left({\mathbf{A}}_{K\left(t\right)}\right)\right]}_{ji}\\ =& \left(\sum _{i=1}^d\sum _{j=1}^d{\left[{\mathbf{A}}_{K\left(t\right)}^{-1}\right]}_{ij}\left({\partial }_{{\xi }_i}{\omega }_j\right)\right)J_{K\left(t\right)}\\ =& \left(\sum _{i=1}^d{\left[{\mathbf{A}}_{K\left(t\right)}^{-1}\left({\partial }_{{\xi }_i}\mathbf{\omega }\right)\right]}_i\right)J_{K\left(t\right)},\end{array}$$

$$\frac{d}{dt}\left(J_{K\left(t\right)}\right)=\left(\sum _{i=1}^d{\left[{\mathbf{A}}_{K\left(t\right)}^{-1}\left({\partial }_{{\xi }_i}\mathbf{\omega }\right)\right]}_i\right)J_{K\left(t\right)}=\left({\nabla }_{\mathbf{\xi }}\cdot \left[{\mathbf{A}}_{K\left(t\right)}^{-1}\mathbf{\omega }\right]\right)J_{K\left(t\right)}.$$

$$\left(\nabla \cdot \mathbf{\omega }\right)J_{K\left(t\right)}=\left({\nabla }_{\mathbf{\xi }}\cdot \left[{\mathbf{A}}_{K\left(t\right)}^{-1}\mathbf{\omega }\right]\right)J_{K\left(t\right)}\in P^{d-1}\left(\left[t_{n-1},t_n\right]\right).$$

$${\mathcal{V}}_h(t):=\left\{v\in {\mathrm{L}}^2\left(\Omega \right):v\circ {\mathbf{\chi }}_{K\left(t\right)}\in P^k\left(K_{\text{ref}}\right),\forall K\left(t\right)\in {\mathcal{T}}_{\left(t\right)}\right\},$$

$$v^{{\text{int}}_{K\left(t\right)}}\left(\mathbf{x}\right):=\underset{\epsilon \to 0^{+}}{\lim }v\left(\mathbf{x}-\epsilon {\mathbf{n}}_{K\left(t\right)}^{\nu }\right),v^{{\text{ext}}_{K\left(t\right)}}\left(\mathbf{x}\right):=\underset{\epsilon \to 0^{+}}{\lim }v\left(\mathbf{x}+\epsilon {\mathbf{n}}_{K\left(t\right)}^{\nu }\right),\forall \mathbf{x}\in F_{K\left(t\right)}^{\nu },$$

$$\{\{v\}\}:=\frac{1}{2}\left(v^{{\text{int}}_{K\left(t\right)}}+v^{{\text{ext}}_{K\left(t\right)}}\right),[[v]]:=v^{{\text{ext}}_{K\left(t\right)}}-v^{{\text{int}}_{K\left(t\right)}}.$$

$$u_h\left(\mathbf{x},t\right)=\sum _{j=1}^r{u}_j^{K\left(t\right)}\left(t\right){\varphi }_j^{K\left(t\right)}\left(\mathbf{x},t\right),\text{for all }t\in \left[t_n,t_{n+1}\right)\text{ and }\mathbf{x}\in K\left(t\right),$$

$${\left(J_{K\left(t\right)}\left({\partial }_t{u}_h\right),v^{\ast }\right)}_{K_{\text{ref}}}+{\left(J_{K\left(t\right)}\left({\nabla }_x\cdot \mathbf{f}\left(u_h\right)\right),v^{\ast }\right)}_{K_{\text{ref}}}=0,$$

$$\frac{d}{dt}{u}_h^{\ast }={\partial }_t{u}_h+\mathbf{\omega }\cdot \nabla u_h={\partial }_t{u}_h+\mathbf{\omega }\cdot {\mathbf{A}}_{K\left(t\right)}^{-T}{\nabla }_{\mathbf{\xi }}{u}_h^{\ast },$$

$$J_{K\left(t\right)}\left(\partial_{t}u_{h}\right)=\frac{d}{dt}\left(J_{K(t)}u_{h}^{*}\right)-J_{K(t)}\nabla_{\bm{\xi}}\cdot\left[{\bf A}_{K\left(t\right)}^{-1}\Big{(}\bm{\omega}u_{h}^{*}\Big{)}\right].$$

$${\left(\frac{d}{dt}\left(J_{K\left(t\right)}{u}_h^{\ast }\right),v^{\ast }\right)}_{K_{\text{ref}}}+{\left(J_{K\left(t\right)}\left({\nabla }_{\mathbf{\xi }}\cdot \tilde{\mathbf{g}}\left(\mathbf{\omega },u_h^{\ast }\right)\right),v^{\ast }\right)}_{K_{\text{ref}}}=0,$$

$$\mathbf{g}\left(\mathbf{\omega },u\right):=\mathbf{f}\left(u\right)-\mathbf{\omega }\left(\mathbf{x},t\right)u.$$

$$\begin{array}{rl}{\left(\frac{d}{dt}\left(J_{K\left(t\right)}{u}_h^{\ast }\right),v^{\ast }\right)}_{K_{\text{ref}}}=& {\left(J_{K\left(t\right)}\tilde{\mathbf{g}}\left(\mathbf{\omega },u_h^{\ast }\right),{\nabla }_{\mathbf{\xi }}{v}^{\ast }\right)}_{K_{\text{ref}}}\\ & -{⟨\hat{g}\left(\mathbf{\omega },u_h^{\ast ,{\text{int}}_{K_{\text{ref}}}},u_h^{\ast ,{\text{ext}}_{K_{\text{ref}}}},J_{K(t)}\tilde{\mathbf{n}}\left(t\right)\right),v^{\ast ,{\text{int}}_{K_{\text{ref}}}}⟩}_{\partial K_{\text{ref}}},\end{array}$$

$${\left(\frac{d}{dt}\left(J_{K\left(t\right)}{u}_h^{\ast }\right),v^{\ast }\right)}_{K_{\text{ref}}}=\frac{d}{dt}{\left(J_{K\left(t\right)}{u}_h^{\ast },v^{\ast }\right)}_{K_{\text{ref}}}=\frac{d}{dt}{\left(u_h,v\right)}_{K\left(t\right)}.$$

$$\begin{array}{rl}\frac{d}{dt}{\left(u_h,v\right)}_{K\left(t\right)}=& {\left(\mathbf{g}\left(\mathbf{\omega },u_h\right),\nabla v\right)}_{K\left(t\right)}\\ & -{⟨\hat{g}\left(\mathbf{\omega },u_h^{{\text{int}}_{K(t)}},u_h^{{\text{ext}}_{K(t)}},{\mathbf{n}}_{K(t)}\right),v^{{\text{int}}_{K(t)}}⟩}_{\partial K(t)}.\end{array}$$

$$\hat{g}\left(\mathbf{\omega },u_h^{{\text{int}}_{K(t)}},u_h^{{\text{ext}}_{K(t)}},{\mathbf{n}}_{K(t)}\right)=-\hat{g}\left(\mathbf{\omega },u_h^{{\text{ext}}_{K(t)}},u_h^{{\text{int}}_{K(t)}},-{\mathbf{n}}_{K(t)}\right).$$

Full text

1 Introduction

The present paper investigates the development and analysis of an arbitrary Lagrangian-Eulerian discontinuous Galerkin (ALE-DG) method for scalar conservation laws in several space dimensions

[TABLE]

with initial condition $u_{0}({\bf x})$ and suitable boundary conditions. The domain $\Omega\subseteq\mathbb{R}^{d}$ is an open convex polyhedron and the flux $\bm{f}\left(u\right):=\left(f_{1}\left(u\right),...,f_{d}\left(u\right)\right)^{T}$ is a suitable vector field. In general the problem (1.1) has no classical solutions. Discontinuities like shock waves could appear in the solution, regardless of the smoothness of the initial data. Hence, the problem needs to be investigated with a class of generalized solutions. The existence of an unique physical relevant solution for the problem was proven by Kružkov in [18]. This solution is called entropy solution. In particular, Kružkov proved that the unique entropy solution satisfies the maximum principle. This means the entropy solution is bounded by the interval $\left[m,M\right]$, where

[TABLE]

The discontinuous Galerkin (DG) method, introduced by Reed and Hill [30] in the context of a neutron transport equation, is a finite element method with discontinuous basis functions. The choice of discontinuous basis functions gives the method a local structure (elements only communicate with immediate neighbors) and the property to handle complex mesh geometries. These features of the DG method are attractive for parallel and high performance computing. Therefore, in particular, the explicit Runge-Kutta DG (RK-DG) method for convection-dominated problems developed and analyzed by Cockburn, Shu and several co-authors in a series of publications (cf. the review article [6] for a summation of their pioneering works) became very popular in the last decades.

The RK-DG method of Cockburn and Shu was developed for a static computational mesh, but in engineering applications like aeroelastic computations of wings (cf. for instance Robinson et al. [31]) numerical methods with a deformable moving mesh are desirable. Nevertheless, a deformable computational domain can lead to strong distortions in the mesh geometry which can be the source of numerical artifacts and instabilities. In the arbitrary Lagrangian-Eulerian (ALE) approach the mesh can move with the fluid like in the Lagrangian specification or the mesh can be static as in the Eulerian specification. This flexibility has a stabilizing effect on an ALE method, since it is possible to switch to the Eulerian specification whenever distortions appear in the mesh geometry. The ALE kinematics were rigorously described by Donea et al. [7]. Moreover, in the literature there are different strategies to combine the ALE approach with the RK-DG method. Among others Lomtev, Kirby, Karniadakis [23], Nguyen [27], Persson et al. [29, 34], Kopriva et al. [17, 26] and Boscheri, Dumbser [2] developed and analyzed ALE-DG methods for convection-dominated problems on a moving domain.

In this paper, an ALE-DG method for solving the problem (1.1) on moving simplex meshes is introduced. This method is an extension of the ALE-DG method developed by Klingenberg et al. [14, 15]. In order to describe the ALE kinematics, we assume that the distribution of the grid points is explicitly given for an upcoming time level by a suitable moving grid methodology. On the basis of this assumption, we can define local affine linear ALE mappings which connect the time-dependent simplex cells with a time-independent reference simplex cell. This simple construction of the ALE mappings ensures that our ALE-DG method has a local structure like the RK-DG method and the discrete geometric conservation law (D-GCL) is satisfied, when a suitable high order accurate Runge-Kutta (RK) method is used. The geometric conservation law (GCL) describes the time evolution of the metric terms in a grid deformation method and has an important influence on the stability and accuracy of a method. The significance of the GCL was first analyzed by Lombard and Thomas [22], thenceforth the GCL was investigated in the context of moving mesh finite volume and finite element methods by Farhat et al. [9, 12, 21], Mavriplis, Yang [25] and Étienne, Garon, Pelletier [8].

Besides the D-GCL, a discrete maximum principle is discussed for our ALE-DG method. In general, even on a static mesh, it is not easy to design a high order method which satisfies a discrete maximum principle for the problem (1.1) without affecting the high order accuracy of the method. Two approaches are commonly used in the literature. The first approach is the flux correction approach. Based on this approach Xu [36] developed a technique to ensure that a high order method satisfies the maximum principle and maintains the high order accuracy. Moreover, an algebraic flux correction approach for finite element methods was introduced by Kuzmin in [20] and [19, Chapter 4]. Another approach was developed by Zhang and Shu [41]. This approach based on a bound preserving limiter which does not affect the high order accuracy of a high order method. In particular, the bound preserving limiter was developed for rectangular meshes by Zhang, Shu [41] and for triangular meshes by Zhang, Xia, Shu [42]. This approach allows the development of high order accurate maximum principle satisfying schemes by a simple investigation of the forward Euler step, since the common convexity argument (cf. Gottlieb and Shu [10]) can be used to extend the result for the forward Euler step to the high order total-variation-diminishing RK (TVD-RK) methods. However, Farhat, Geuzaine and Grandmont [9] proved that for ALE finite volume methods a discrete maximum principle is satisfied, if and only if the D-GCL is satisfied. Unfortunately, for our ALE-DG method the D-GCL is only fulfillment, if the accuracy of the RK method corresponds with the spatial dimension. Hence, we cannot expect that our forward Euler ALE-DG method satisfies a discrete maximum principle, when the bound preserving limiter is applied. Nevertheless, it turns out that the GCL is an ODE in our ALE-DG method. This ODE can be solved exactly by a RK method with an order not less than the value of the spatial dimension. In two dimensions, we use the second and the third order TVD-RK methods developed by Shu in [32] to solve the GCL and the actual ALE-DG method. Then the RK stage solutions for the GCL are used to update the metric terms in the RK stages of the actual ALE-DG method. This time integration strategy allows to develop second and third order accurate fully-discrete ALE-DG methods. We prove that these methods satisfy a discrete maximum principle when the bound preserving limiter is applied. Furthermore, we present numerical experiments which support the expectation that the ALE-DG method also satisfies a discrete maximum principle when the five stage fourth order TVD-RK method developed by Spiteri and Ruuth in [33] and the bound preserving limiter are used.

In addition, we present an a priori error estimate for our ALE-DG method. A priori error analysis to smooth solutions of the second and third order RK-DG method on a static mesh for scalar and symmetrizable systems of conservation laws were mainly done by Zhang, Shu et al. More precisely, Zhang and Shu proved in [37] and [38] that under a slightly more restrictive Courant-Friedrichs-Lewy (CFL) constraint than the commonly used constraint the a priori error of the second order RK-DG method behaves as $\mathcal{O}\left(\triangle t^{2}+h^{k+\frac{1}{2}}\right)$ in the $\mathrm{L}^{2}$-norm, when a local polynomial basis of degree $k\geq 1$ and an arbitrary monotone flux are applied. In this context the quantity $\triangle t$ denotes the time step and $h$ denotes the maximum cell length. Likewise, Zhang and Shu proved in [39] that under the usual CFL constraint the a priori error of the third order RK-DG method for scalar conservation laws behaves as $\mathcal{O}\left(\triangle t^{3}+h^{k+\frac{1}{2}}\right)$ in the $\mathrm{L}^{2}$-norm. This result was extended to symmetrizable systems of conservation laws by Luo, Shu, and Zhang [24]. Furthermore, a priori error estimates for the third order RK-DG method in the context of linear scalar conservation Laws with discontinuous initial data were proven in [40]. Error estimates for fully discrete ALE-DG methods to solve linear conservation laws were proven by Zhou, Xia and Shu in [43]. In this work, we merely prove that for smooth solutions of the problem (1.1) the a priori error of the semi-discrete ALE-DG method behaves as $\mathcal{O}\left(h^{k+\frac{1}{2}}\right)$, when polynomials of degree $k\geq\max\left\{1,\frac{d}{2}\right\}$ are used on the reference cell and an arbitrary monotone flux is applied.

The rest of the paper is organized as follows. In Section 2, we introduce the local affine linear ALE mappings, a time-dependent test function space and our semi-discrete ALE-DG method. In Section 3, we present some theoretical results for the semi-discrete ALE-DG method. In particular, the $\mathrm{L}^{2}$-stability is proven. Afterward, in Section 4, the fully-discrete ALE-DG method is investigated. We prove that the D-GCL is satisfied under certain conditions which are related to the spatial dimension. Furthermore, in two dimensions second and third order accurate fully-discrete ALE-DG methods on moving triangular meshes are presented. We prove that these methods satisfy a discrete maximum principle, when the bound preserving limiter for triangular meshes developed by Zhang, Xia, Shu in [42] is applied. In Section 5, we validate the theoretical results by some computational examples and show that the ALE-DG method is numerically stable and high order accurate. Finally, we give some concluding remarks in Section 6.

Constants and notation

In the present paper, vectors, vector valued functions and matrices are denoted by bold letters. Scalar quantities are denoted by regular letters. The set $K\left(t\right)$ denotes a time-dependent open simplex cell in a $d$ dimensional domain with the edges $F_{K\left(t\right)}^{\nu}$, $\nu=1,\dots,d+1$. Volume integrals with respect to the open set $K\left(t\right)$ and surface integrals with respect to the edges $F_{K\left(t\right)}^{\nu}$, $\nu=1,\dots,d+1$, are denoted by the bracket notation. Hence, for all $v,w\in\mathrm{L}^{2}\left(K\left(t\right)\right)\cup\mathrm{L}^{2}\left(\partial K\left(t\right)\right)$ and $\nu=1,\dots,d+1$ the notations $\left(v,w\right)_{K\left(t\right)}:=\int_{K\left(t\right)}vw\,d\textbf{x}$, $\left\langle v,w\right\rangle_{F_{K\left(t\right)}^{\nu}}:=\int_{F_{K\left(t\right)}^{\nu}}vw\,d\bm{\Gamma}$ and $\left\langle v,w\right\rangle_{\partial K\left(t\right)}:=\sum_{\nu=1}^{d+1}\left\langle v,w\right\rangle_{F_{K\left(t\right)}^{\nu}}$ are applied. Furthermore, to avoid confusion with different constants, we denote by $C$ a positive constant, which is independent of the mesh size and the numerical solutions for the conservation law (1.1), but it may depend on the solution of the PDE and may have a different value in each occurrence.

2 The ALE-DG discretization

In this section, we present the semi-discrete ALE-DG discretization of the problem (1.1). At first, the ALE framework to derive the ALE-DG method in several dimensions is briefly listed. Afterward, the ALE framework is used to derive the semi-discrete ALE-DG method for solving the problem (1.1).

2.1 The ALE-DG setting

In this section, we present the time dependent cells and introduce some identities for the metric quantities to transform derivatives on a reference cell.

2.1.1 The time-dependent simplex mesh

We assume that there exists a regular mesh $\mathcal{T}_{(t_{n})}$ of simplices at any time level $t_{n}$, $n=0,\dots,\mathcal{N}$, which covers exactly the convex polyhedron domain $\Omega$ such that

[TABLE]

The mesh topology of $\mathcal{T}_{(t_{n})}$ and $\mathcal{T}_{(t_{n+1})}$ is assumed to be the same. This means:

• a)

$\mathcal{T}_{(t_{n})}$ and $\mathcal{T}_{(t_{n+1})}$ are simplex meshes of the domain $\Omega$.

• b)

$\mathcal{T}_{(t_{n})}$ and $\mathcal{T}_{(t_{n+1})}$ have the same number of cells.

• c)

The cells of both simplex meshes are positively oriented with respect to the reference simplex

[TABLE]

The $d+1$ vertices of each simplex $K(t_{n})\in\mathcal{T}_{(t_{n})}$ are denoted by ${\bf v}_{1}^{n},\dots,{\bf v}_{d+1}^{n}$. We define for all $t\in\left[t_{n},t_{n+1}\right]$ and $\ell=1,...,d+1$ time-dependent straight lines

[TABLE]

These straight lines are for any $t\in\left[t_{n},t_{n+1}\right]$ the vertices of a time-dependent simplex cell given by

[TABLE]

where $\text{int}\left(\cdot\right)$ and $\text{conv}\left(\cdot\right)$ denote the interior and the convex hull of a set. In the following, the set of all time-dependent cells $K(t)$ is denoted by $\mathcal{T}_{\left(t\right)}$. Furthermore, for any cell $K\left(t\right)\in\mathcal{T}_{(t)}$ the diameter of the cell and the radius of the largest ball, contained in $K\left(t\right)$, are denoted by $h_{K\left(t\right)}$ as well as $\rho_{K\left(t\right)}$. Additionally, we define the following global length

[TABLE]

Henceforth, we assume:

1. (A1)

The domain $\overline{\Omega}$ is for all $t\in\left[0,T\right]$ exactly covered by the time-dependent cells (2.3) such that $\overline{\Omega}=\bigcup_{K\left(t\right)\in\mathcal{T}_{\left(t\right)}}\overline{K\left(t\right)}$. 2. (A2)

For all $t\in\left[0,T\right]$ and all cells $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$ $J_{K\left(t\right)}=\text{det}\left(\mathrm{{\bf A}}_{K\left(t\right)}\right)>0$. 3. (A3)

It exists constants $\kappa>0$ and $\tau>0$, independent of $h$, such that for all $t\in\left[0,T\right]$

[TABLE]

2.1.2 The ALE mapping and the grid velocity field

The time-dependent simplex cells (2.3) can be mapped to the time-independent reference simplex element (2.1) by the affine linear time-dependent mapping

[TABLE]

where the matrix $\mathrm{{\bf A}}_{K\left(t\right)}$ is given by

[TABLE]

We note that the matrix $\mathrm{{\bf A}}_{K\left(t\right)}$ is the Jacobian matrix of the mapping $\bm{\chi}_{K\left(t\right)}$ and the corresponding determinant is

[TABLE]

where $\left|K\left(t\right)\right|$ denotes the volume of the cell $K\left(t\right)$. In particular, $J_{K\left(t\right)}$ is independent of the spatial variables and belongs to $P^{d}\left(\left[t_{n-1},t_{n}\right]\right)$. It is worth to mention that in general for non-simplicial moving meshes $J_{K\left(t\right)}$ depends on spatial and temporal variables, since the shape of the elements can change when the corners move with different speed. In Figure 2.1 the two dimensional situation for a triangular element and a rectangular element is illustrated. The implementation of metric quantities which depend on spatial variables is not easy and requires caution (e.g. cf. Kopriva [16]).

Since the matrix ${\bf A}_{K\left(t\right)}$ is the Jacobian matrix of the mapping (2.5), we have the following metric transformations

[TABLE]

where $\bm{f}:\mathbb{R}\to\mathbb{R}^{d}$ is an arbitrary vector field with $\bm{f}^{*}=\bm{f}\circ\bm{\chi}_{K\left(t\right)}$, $u$ is a scalar function with $u^{*}=u\circ\bm{\chi}_{K\left(t\right)}$, $\mathbf{n}_{K(t)}$ is the normal of the cell $K\left(t\right)$ and $\mathbf{n}_{K_{\text{ref}}}$ is the reference normal. A proof of these metric transformations is given in Ciarlet [3, p. 461]. Moreover, the mapping (2.5) provides the grid velocity field in the point $\textbf{x}=\bm{\chi}_{K\left(t\right)}\left(\bm{\xi},t\right)$

[TABLE]

The definition (2.9) provides the following relation

[TABLE]

where $\omega_{j}$ are the coefficients of the grid velocity and $\left[\frac{d}{dt}\left(\mathrm{{\bf A}}_{K\left(t\right)}\right)\right]_{ji}$ are the coefficients of the matrix $\frac{d}{dt}\left(\mathrm{{\bf A}}_{K\left(t\right)}\right)$. We are also interested to find an identity for the time derivative of the determinant $J_{K\left(t\right)}$. Hence, we apply Jacobi’s formula (cf. Bellman [1]), using the equation relating the adjugate of $\mathrm{{\bf A}}_{K\left(t\right)}$ to the inverse $\mathrm{{\bf A}}_{K\left(t\right)}^{-1}$ and apply the identity (2.10). This results in the identity

[TABLE]

where $\text{tr}\left(\cdot\right)$ denotes the trace of a matrix, $\text{adj}\left(\cdot\right)$ denotes the adjoint of a matrix and $M_{K\left(t\right)}^{ij}$ is the $(i,j)$ minor of $\mathrm{{\bf A}}_{K\left(t\right)}$. Moreover, since the matrix $\mathrm{{\bf A}}_{K(t)}$ does not depend on spatial variables, we obtain

[TABLE]

Finally, we summarize some properties of the grid velocity. These properties will be used in the next sections.

Lemma 2.1**.**

The grid velocity $\bm{\omega}$ defined by (2.9) has the properties:

• (i)

For all $t\in\left[t_{n},t_{n+1}\right]$ the grid velocity belongs to the space $P^{1}\left(K_{\text{ref}},\mathbb{R}^{d}\right)$.

• (ii)

The grid velocity is time-independent for all points contained in the set $\partial{K_{\text{ref}}}$.

• (iii)

The divergence of the grid velocity satisfies

[TABLE]

Proof.

Since the matrix $\frac{d}{dt}\left(\mathrm{{\bf A}}_{K\left(t\right)}\right)$ is time independent, the property (i) follows from the definition of the grid velocity (2.9).

The property (ii) follows directly from the definitions of the time-dependent straight lines (2.2), the cells (2.3) and the grid velocity (2.9).

Finally, the property (iii) follows from the metric transformations (2.8) and the identity (2.12), since $\left[\frac{d}{dt}\left(\mathrm{{\bf A}}_{K\left(t\right)}\right)\right]_{ij}$ is time independent and the minors $M_{K\left(t\right)}^{ij}$, $i,j=1,\dots,d$, belong to the space $P^{d-1}\left(\left[t_{n-1},t_{n}\right]\right)$. ∎

2.2 The approximation space

We define the approximation space

[TABLE]

where $P^{k}\left(K_{\text{ref}}\right)$ denotes the space of polynomials in $K_{\text{ref}}$ of degree at most $k$. The functions from the space $\mathcal{V}_{h}(t)$ are discontinuous along the interface of two adjacent cells. Thus, we define for a function $v\in\mathcal{V}_{h}(t)$, an arbitrary cell $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$ and all $\nu=1,...,d+1$ the following limits

[TABLE]

where the vector $\mathbf{n}_{K\left(t\right)}^{\nu}$, $\nu=1,...,d+1$, is the outward normal of the cell $K\left(t\right)$ with respect to the simplex face $F_{K\left(t\right)}^{\nu}$. Then, the cell average and jump of the function $v$ along the simplex face $F_{K\left(t\right)}^{\nu}$ are defined by

[TABLE]

2.3 The semi-discrete ALE-DG method

For each cell $K(t)\in\mathcal{T}_{\left(t\right)}$, we approximate the solution $u$ of the problem (1.1) by the function

[TABLE]

where $r:=\frac{\left(k+d\right)!}{d!k!}$ and $\left\{\phi_{1}^{K(t)}\left({\bf x},t\right),...,\phi_{r}^{K(t)}\left({\bf x},t\right)\right\}$ is a basis of the space $\mathcal{V}_{h}\left(t\right)$ in the cell $K(t)$. The coefficients $u_{1}^{K(t)}\!\left(t\right)$,…,$u_{r}^{K(t)}\!\left(t\right)$ in (2.15) are the unknowns of the method. In order to determine these coefficients, we plug the function (2.15) in (1.1), multiply the equation by a test function $v\in\mathcal{V}_{h}(t)$ and use the change of variables theorem for integrals. This results in the equation

[TABLE]

where $v^{*}=v\circ\chi_{K\left(t\right)}$. The chain rule formula and the metric transformations (2.8) provide

[TABLE]

where $u_{h}^{*}=u_{h}\circ\chi_{K\left(t\right)}$. Next the identities (2.12) and (2.17) provide

[TABLE]

Then the equation (2.16) becomes

[TABLE]

where \tilde{\bm{g}}\left(\bm{\omega},u_{h}^{*}\right):=\mathrm{{\bf A}}_{K\left(t\right)}^{-1}\Big{(}\bm{g}\left(\bm{\omega},u_{h}^{*}\right)\Big{)} with

[TABLE]

At this point, we proceed similar as in the derivation of the standard DG method on a static mesh. First, we apply the integration by parts formula in the second integral in (2.19). Then, we replace the flux function $\tilde{\bm{g}}\left(\bm{\omega},u_{h}^{*,\text{int}_{K_{\text{ref}}}}\right)\cdot\mathbf{n}_{K_{\text{ref}}}$ in the surface integrals by a numerical flux function $\widehat{g}\left(\bm{\omega},u_{h}^{*,\text{int}_{K_{\text{ref}}}},u_{h}^{*,\text{ext}_{K_{\text{ref}}}},J_{K(t)}\tilde{\mathbf{n}}\left(t\right)\right)$ with $\tilde{\mathbf{n}}\left(t\right)=A_{K\left(t\right)}^{-T}\mathbf{n}_{K_{\text{ref}}}$. The numerical flux function needs to satisfy certain properties. These properties are discussed in the Section 2.4. Finally, on the reference cell, the semi-discrete ALE-DG method appears as the following problem:

Problem 1** (The semi-discrete ALE-DG method on the reference cell).**

Find a function $u_{h}\in\mathcal{V}_{h}(t)$, such that for all $v\in\mathcal{V}_{h}(t)$ and all cells $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$ holds

[TABLE]

where $\tilde{\mathbf{n}}\left(t\right)=A_{K\left(t\right)}^{-T}\mathbf{n}_{K_{\text{ref}}}$.

Since the test functions $v^{*}=v\circ\chi_{K\left(t\right)}$ are time independent on the reference cell $K_{\text{ref}}$, we obtain

[TABLE]

Therefore, on the physical domain, 1 is equivalent to:

Problem 2** (The semi-discrete ALE-DG method).**

Find a function $u_{h}\in\mathcal{V}_{h}(t)$, such that for all $v\in\mathcal{V}_{h}(t)$ and all cells $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$ holds

[TABLE]

2.4 The numerical flux function

The numerical flux in the ALE-DG method should satisfy:

1. (P1)

The flux is consistent with $\bm{g}\left(\bm{\omega},u\right)\cdot\mathbf{n}_{K(t)}$. 2. (P2)

The function $u\mapsto\widehat{g}\left(\bm{\omega},u,\cdot,\mathbf{n}_{K(t)}\right)$ is increasing and Lipschitz continuous. 3. (P3)

The function $u\mapsto\widehat{g}\left(\bm{\omega},\cdot,u,\mathbf{n}_{K(t)}\right)$ is decreasing and Lipschitz continuous. 4. (P4)

The flux is conservative such that

[TABLE]

Remark 1**.**

The properties ($P1$) - ($P4$) of the numerical flux supply for any cell $K(t)\in\mathcal{T}_{\left(t\right)}$ and all $v\in\left[\min\left(u_{1},u_{2}\right),\max\left(u_{1},u_{2}\right)\right]$

[TABLE]

The inequality (2.24) is the e-flux condition, which was introduced by Osher [28].

In general every numerical flux with these properties can be used in the ALE-DG method. A common example is the Lax-Friedrichs flux. This flux is given by

[TABLE]

We note that the functions $\widehat{g}_{\pm}\left(\bm{\omega},u,\mathbf{n}_{K(t)}\right)$ are increasing in the second argument.

3 Theoretical results for the semi-discrete method

In this section, we present some theoretical results for the semi-discrete ALE-DG method. We start with a proof for the $\mathrm{L}^{2}$-stability of the semi-discrete ALE-DG method. This proof requires techniques which were introduced by Jiang and Shu in [13] to proof a cell entropy inequality for the DG method. We note that for scalar conservation laws the function $\eta\left(u\right)=\frac{1}{2}u^{2}$ is an entropy. Thus, the $\mathrm{L}^{2}$-stability provides also entropy stability in the sense that the total entropy at a certain time point is bounded by the total entropy at initial time.

Proposition 3.1**.**

The solution $u_{h}$ of the semi-discrete ALE-DG method satisfies for any $t\in[0,T]$

[TABLE]

when the problem (1.1) is considered with periodic boundary conditions.

Proof.

Let $K(t)\in\mathcal{T}_{\left(t\right)}$ be an arbitrary cell. We use the ALE-DG solution $u_{h}^{*}=u_{h}\circ\chi_{K\left(t\right)}$ as test function in the equation (2.21) and obtain

[TABLE]

We obtain by the identity (2.12) and the change of variables theorem for integrals

[TABLE]

Next, we define the vector valued functions

[TABLE]

Then, we obtain by the metric transformations (2.8), the functions (3.4) and the change of variables theorem for integrals

[TABLE]

Furthermore, the divergence theorem gives

[TABLE]

where the vectors $\mathbf{n}_{K\left(t\right)}^{\nu}$, $\nu=1,\dots,d+1$, are the outward normals of the cell $K\left(t\right)$ with respect to the simplex faces $F_{K(t)}^{\nu}$, $\nu=1,\dots,d+1$. Likewise, the metric transformations (2.8) provide

[TABLE]

Next, we rearrange the equation (3.2) by applying the identities (3.3), (3.5), (3.6) and (3.7). This results in the equation

[TABLE]

Then, we sum the equation (3.8) over all cells $K(t)\in\mathcal{T}_{(t)}$ and obtain

[TABLE]

since we consider the problem (1.1) with periodic boundary conditions. The definition of the function $\bm{G}(\bm{\omega},u)$ in (3.4) yields $\partial_{u}\bm{G}(\bm{\omega},u)=\bm{g}(\bm{\omega},u)$. Thus, for any cell $K(t)\in\mathcal{T}_{(t)}$ the mean value theorem and the e-flux condition (2.24) provide for all $\nu=1,\dots,d+1$

[TABLE]

where $\theta_{K\left(t\right)}^{\nu}$ is a value between $\min\left(u_{h}^{\text{int}_{K(t)}},u_{h}^{\text{ext}_{K(t)}}\right)$ and $\max\left(u_{h}^{\text{int}_{K(t)}},u_{h}^{\text{ext}_{K(t)}}\right)$. Hence, the inequality (3.1) follows by integrating the equation (3.9) over the interval $[0,t]$. ∎

Furthermore, for sufficiently smooth solutions of the initial value problem (1.1), we have a suboptimal a priori error estimate in the sense of the $\mathrm{L}^{\infty}\left(0,T;\mathrm{L}^{2}\left(\Omega\right)\right)$-norm for the semi-discrete ALE-DG method.

Theorem 3.2**.**

Consider the initial value problem (1.1) with periodic boundary conditions and let $u\in\mathrm{W}^{1,\infty}\left(0,T;\mathrm{H}^{k+1}\left(\Omega\right)\right)$ be the exact solution, the flux function $\bm{f}\in\mathcal{C}^{3}\left(\mathbb{R},\mathbb{R}^{d}\right)$ and the grid velocity field $\bm{\omega}$ be bounded in $\Omega\times\left[0,T\right]$ and have bounded derivatives and $u_{h}$ be the solution of the semi-discrete ALE-DG method. The test function space (2.14) is given by piecewise polynomials of degree $k\geq\max\left\{1,\frac{d}{2}\right\}$. Furthermore, ($A1$) - ($A3$) are satisfied and the global length $h$ is given by (2.4). Then, it exists a constant $C$, which depends on the final time $T$ and is independent of $h$ and $u_{h}$, such that

[TABLE]

Theorem 3.2 can be proven with standard techniques from approximation theory (cf. Ciarlet [4]). In particular, a one dimensional proof is given in [14] and in [15] two dimensional error analysis for a semi-discrete ALE-DG method to solve the Hamilton-Jacobi equations is given. In addition, error analysis for the fully-discrete ALE-DG methods is given in [43]. Since, there are already these publications on error analysis for the ALE-DG method in the literature and the proof of Theorem 3.2 vary merely in technical details, we skip the proof in this paper.

Remark 2**.**

The proof of Theorem 3.2 requires the a priori assumption

[TABLE]

These a priori assumption is not necessary, if the problem (1.1) is considered with a linear flux function $\bm{f}\left(u\right)={\bf c}u$, ${\bf c}\in\mathbb{R}^{d}$. In the one-dimensional case, the a priori assumption

[TABLE]

can be applied, since for $d=1$ the inequality (3.13) and the inverse inequality [4, Theorem 3.2.6.] provide the a priori assumption (3.12). The assumption (3.13) was applied in [37, 35]. Moreover, the statement of theorem 3.2 can be also proven, for flux functions with less smoothness. Nevertheless, this requires a more restrictive a priori assumption and more restrictive bounds for the parameter $k$. This assumption was applied by Klingenberg et al. in [14].

4 The fully-discrete method

In this section, we discuss the time discretization of the ALE-DG method. In the first part of this section, we prove that the fully-discrete ALE-DG method satisfies the discrete geometric conservation law (D-GCL). Afterward, in two dimensions, we prove that the second and the third order accurate fully-discrete ALE-DG methods satisfy the maximum principle, when the bound preserving limiter developed by Zhang, Xia and Shu in [42] is applied.

4.1 Total-variation-diminishing Runge-Kutta (TVD-RK) methods

We apply the high order TVD-RK methods developed by Shu in [32] for the time discretization, which is also known as strong stability preserving Runge-Kutta (SSP-RK) methods [10, 11]. For a given ODE ${u}_{t}=q\left(u,t\right)$, a $s$-stage TVD-RK method can be written as

[TABLE]

where $t_{n+\gamma_{j}}:=t_{n}+\gamma_{j}\triangle t$. The coefficients of the $s$-stage TVD-RK method (4.1) need to satisfy

[TABLE]

for all $i=1,\dots,s$ and $j=0,\dots,s-1$. In the Section 4.3, we will investigate the commonly used second and third order TVD-RK methods from Shu [32]. The coefficients for these methods are given in the Tables 4.1 and 4.2.

[TABLE]

[TABLE]

4.2 The discrete geometric conservation law (D-GCL)

For the time discretization of the semi-discrete ALE-DG Problem 1, it is convenient to introduce the notation

[TABLE]

where $\tilde{\mathbf{n}}\left(t\right)=A_{K\left(t\right)}^{-T}\mathbf{n}_{K_{\text{ref}}}$. At this point, we note that the grid velocity satisfies the property (ii) in Lemma 2.1. Therefore, in order to avoid confusion with the time-dependency of the grid velocity, we will highlight the time variable if the grid velocity depends on time, otherwise the time variable will be omitted. For all $v\in\mathcal{V}_{h}(t)$ with $v^{*}=v\circ\chi_{K\left(t\right)}$ the equation (2.21) can be written as

[TABLE]

In accordance with Guillard and Farhat [12], we introduce the following definition.

Definition 4.1**.**

A fully-discrete moving mesh method for the initial value problem (1.1) satisfies a D-GCL, if for all $c\in\mathbb{R}$ and $n=0,\dots,\mathcal{N}-1$

[TABLE]

Lemma 4.2**.**

Let $c\in\mathbb{R}$ be a constant. Then for all $t\in\left[t_{n},t_{n+1}\right]$, $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$ and $v\in\mathcal{V}_{h}(t)$ with $v^{*}=v\circ\chi_{K\left(t\right)}$, it holds

[TABLE]

Proof.

Since $\bm{f}\left(c\right)$ contains merely constant coefficients and

[TABLE]

the integration by parts formula provides

[TABLE]

Furthermore, since $\tilde{\mathbf{n}}\left(t\right)=A_{K\left(t\right)}^{-T}\mathbf{n}_{K_{ref}}$, the property (P1) of the numerical flux provides

[TABLE]

Thus, we obtain the identity (4.5) by (4.6) and (4.7). ∎

Next, we assume that $u_{h}^{*}=c$ solves the semi-discrete ALE-DG method Problem 1. Then, we obtain by (4.4) and (4.5)

[TABLE]

The equation (4.8) and the ODE (2.12) are equivalent, since $c$ is an arbitrary constant, $v\in\mathcal{V}_{h}(t)$ with $v^{*}=v\circ\chi_{K\left(t\right)}$ is an arbitrary test function and the quantities

[TABLE]

are merely time-dependent. We note that the time evolution of the metric terms $J_{K\left(t\right)}$ needs to be respected in the time discretization of the semi-discrete ALE-DG method Problem 1. Therefore, we discretize the ODE (2.12) and (4.4) simultaneously by the same TVD-RK method. The stage solutions of the TVD-RK discretization for (2.12) will be used to update the metric terms in the TVD-RK discretization for (4.4).

The fully-discrete ALE-DG method: First, the ODE (2.12) is discretized by a s-stage TVD-RK method:

[TABLE]

where $K^{n+\gamma_{j}}:=K\left(t_{n}+\gamma_{j}\triangle t\right)$ and $\bm{\omega}^{n+\gamma_{j}}:=\bm{\omega}\left(t_{n}+\gamma_{j}\triangle t\right)$. The stage solutions $\left\{J_{K^{n,i}}\right\}_{i=0}^{s}$ are used to update the metric terms in the TVD-RK discretization of (4.4). The Runge-Kutta method needs to solve the ODE (2.12) exact such that

[TABLE]

where $\mathcal{T}_{(t_{n+1})}$ is the regular mesh of simplices which has been used in the Section 2.1 to construct the time-dependent cells (2.3). We note that in a $d$-dimensional space a TVD-RK method with order greater than or equal to $d$ is necessary to compute the metric term $J_{K^{n+1}}$ exactly, since the right hand side of the equation (2.12) belongs to the space $P^{d-1}\left(\left[t_{n},t_{n+1}\right]\right)$.

Problem 3** (The fully-discrete ALE-DG method on the reference cell).**

Find a function $u_{h}\in\mathcal{V}_{h}(t)$, such that for all $v\in\mathcal{V}_{h}(t)$ there holds

[TABLE]

where the metric terms $\left\{J_{K^{n,i}}\right\}_{i=0}^{s}$ are computed by (4.10).

We note that on a static mesh the cells $K(t)$ are time-independent. Thus, on a static mesh, the Problem 3 corresponds to the TVD-RK DG method developed by Cockburn and Shu in [6]. Next, we prove that the fully-discrete ALE-DG method satisfies the following statement.

Theorem 4.3**.**

Suppose a s-stage TVD-RK method with order greater than or equal to $d$ is used in (4.10) and Problem 3, and the solutions $\left\{J_{K^{n,i}}\right\}_{i=0}^{s}$ are used to compute the metric terms in Problem 3. Moreover, the solution at time level $t_{n}$ satisfies $u_{h}^{n,*}=c\in\mathbb{R}$. Then it is $u_{h}^{n,i,*}=c$ for all $i=0,\dots,s$.

Proof.

Let $i\in\left\{0,\dots,s\right\}$ be an arbitrary fixed index and $v\in\mathcal{V}_{h}(t)$ be an arbitrary test function. We are interested to investigate the $i$-th Runge-Kutta stage in the Problem 3. Hence, we can assume that $u_{h}^{n,j,*}=c$ for all $j=0,\dots,s-1$. Then, the equation (4.5) in Lemma 4.2 provides

[TABLE]

Next, we multiply the equation (4.10b) by $cv^{*}$ and integrate the result over the reference element $K_{\text{ref}}$. This provides the identity

[TABLE]

Since the metric terms $\left\{J_{K^{n,j}}\right\}_{j=0}^{i-1}$ in (4.12b) are computed by (4.10), the equation (4.14) can be plugged in the equation (4.13) and it follows

[TABLE]

Thus, it follows $u_{h}^{n,i,*}=c$, since the metric term $J_{K^{n,i}}$ is merely time-dependent and the test function $v\in\mathcal{V}_{h}(t)$ was chosen arbitrary. ∎

Theorem 4.3 provides that constant initial data is preserved in each Runge-Kutta stage. In particular, a direct consequence of the Theorem 4.3 is the following result.

Corollary 4.4**.**

Suppose a s-stage TVD-RK method with order greater than or equal to $d$ is used in (4.10) and Problem 3. Furthermore, the solutions $\left\{J_{K^{n,i}}\right\}_{i=0}^{s}$ given by (4.10) are used to compute the metric terms in Problem 3. Then the fully-discrete ALE-DG method satisfies the D-GCL in the sense of the Definition 4.1.

4.3 The discrete maximum principle

In this section, we investigate the fully-discrete ALE-DG method Problem 3 in two dimensions. First, we review the bound preserving limiter which was developed by Zhang, Xia and Shu in [42]. Then, we consider the second and third order accurate TVD-RK methods given in the Tables 4.1 and 4.2. We prove that for these methods the fully-discrete ALE-DG method satisfies the maximum principle when the bound preserving limiter is applied.

4.3.1 The bound preserving limiter

In the following, we briefly review the bound preserving limiter methodology for an arbitrary $s$-stage TVD-RK method (4.1). Let $u_{h}^{n}$ be the DG solution at time level $t=t_{n}$, $K\subseteq\Omega$ an arbitrary cell and $Q_{K}\subseteq K$ be a set of quadrature points. The corresponding quadrature formula needs to be exact for the integration of polynomials with degree $k$ on the cell $K$ and the set $Q_{K}$ needs to contain all quadrature points which are necessary to evaluate the surface integrals along the edges of $K$ in the DG method. A quadrature rule with these properties has been developed in [42] and will be presented in the Section 4.3.3. Then the cell average

[TABLE]

of the numerical solution $u_{h}$ can be computed exact by the quadrature formula. The limiter methodology is applied in two steps:

1. (L1)

It needs to be ensured that the $s$-stage TVD-RK method satisfies

[TABLE]

where $m$ as well as $M$ are given by (1.2), $u_{h}^{n,\ell}$ is the solution given by the $\ell$-th stage of the RK-DG method and $\overline{u}_{K}^{n,\ell}$ denotes the corresponding cell average with respect to the cell $K$. 2. (L2)

If (L1) is satisfied, the solution $u_{h}$ will be revised for all $\ell=1,...,s$ by

[TABLE]

[TABLE]

where $M_{K}:=\underset{{\bf x}\in Q_{K}}{\max}\,u_{h}^{n,\ell}\left({\bf x}\right)$, $m_{K}:=\underset{{\bf x}\in Q_{K}}{\min}\,u_{h}^{n,\ell}\left({\bf x}\right)$. Then it is ensured that $\widetilde{u}_{h}^{n,\ell}|_{K}\in\left[m,M\right]$, for all $\ell=1,...,s$.

The stability property (L1) is the maximum principle property for the cell-averages of the numerical solution $u_{h}$. We note that by an adjustment of the CFL constraint the stability property (L1) is satisfied for any high order TVD-RK method, if the forward Euler step satisfies (L1), since the TVD-RK methods are convex combinations of the forward Euler step (cf. Gottlieb and Shu [10]). However, the forward Euler method is merely first order accurate and thus this method is not an appropriate choice for the discretization of the ODE (2.12) and the semi-discrete ALE-DG method Problem 1 in two dimensions, since it does not provide that the equation (4.11) holds. Therefore, the convexity argument in Gottlieb and Shu [10] cannot be used to analyze a fully-discrete ALE-DG method which has a high order TVD-RK method as time integrator and uses the bound preserving limiter. For this reason, each high order TVD-RK method needs to be investigated separately.

4.3.2 High order time discretization methods

The step (L1) in the bound preserving limiter methodology is related to the cell averages of the numerical solution. Thus, a scheme satisfied by the cell averages of the ALE-DG solution $u_{h}$ needs to be investigated. In this section, we will present these schemes for the TVD-RK2 and the TVD-RK3 ALE-DG methods. Later, in the Section 4.3.4 and 4.3.5, these schemes will be used to prove that the TVD-RK2 and the TVD-RK3 ALE-DG methods satisfy indeed the maximum principle, when the bound preserving limiter (4.16) is applied.

The TVD-RK2 ALE-DG method: We note that by the identity (2.7) the metric term $J_{K(t)}$ can be also written as $2\left|K\left(t\right)\right|$. Thus, for the TVD-RK2 method in Table 4.1 the scheme (4.10) becomes

[TABLE]

Next, we apply the TVD-RK2 method and the test function $v^{*}=1$ in (4.12) and obtain

[TABLE]

where $\overline{u}_{K^{n}}^{n}$, $\overline{u}_{K^{n,1}}^{n,1}$ and $\overline{u}_{K^{n+1}}^{n+1}$ are the cell average values of the ALE-DG solution $u_{h}$ and the quantities $\left|K^{n}\right|$, $\left|K^{n,1}\right|$, $\left|K^{n+1}\right|$ are computed by (4.17). We note that $\left|K\left(t_{n+1}\right)\right|=\left|K^{n+1}\right|$, since in two space dimensions, the TVD-RK2 method in Table 4.1 solves the ODE (2.12) exactly.

The TVD-RK3 ALE-DG method: For the TVD-RK3 method in Table 4.2 the method (4.10) becomes

[TABLE]

where $t_{n+\frac{1}{2}}:=\frac{1}{2}\left(t_{n}+t_{n+1}\right)$ and $K^{n+\frac{1}{2}}=K\left(t_{n+\frac{1}{2}}\right)$. In two and three space dimensions, the third order TVD-RK3 method in Table 4.2 solves the ODE (2.12) exactly and thus $\left|K\left(t_{n+1}\right)\right|=\left|K^{n+1}\right|$. Next, we apply the test function $v^{*}=1$ in (4.12) and obtain the scheme

[TABLE]

where the values $\overline{u}_{K^{n}}^{n}$, $\overline{u}_{K^{n,1}}^{n,1}$, $\overline{u}_{K^{n,2}}^{n,2}$ and $\overline{u}_{K^{n+1}}^{n+1}$ are the cell average values of the ALE-DG solution $u_{h}$ and the quantities $\left|K^{n}\right|$, $\left|K^{n,1}\right|$, $\left|K^{n,2}\right|$, $\left|K^{n+1}\right|$ are computed by (4.19).

4.3.3 A quadrature rule to decompose the cell average value

In order to prove that the second and the third order fully-discrete ALE-DG methods satisfy the maximum principle, when the ALE-DG solution is revised by the limiter (4.16), we proceed similar to the derivation by Zhang, Xia and Shu in [42]. For this task, we need to apply a special quadrature formula developed by Zhang, Xia and Shu to decompose the cell average values of the ALE-DG solution. In the following this quadrature formula is briefly reviewed.

First of all, it should be noted that in the implementation of the ALE-DG method the edge integrals are approximated by a $k+1$-point $2k+1$ accurate Gauss quadrature formula. Hence, we obtain

[TABLE]

where $\tilde{\mathbf{n}}\left(t\right)={\bf A}_{K\left(t\right)}^{-T}\mathbf{n}_{K_{\text{ref}}}$, $\tilde{\mathbf{n}}_{F_{K_{\text{ref}}}^{\nu}}\left(t\right)={\bf A}_{K\left(t\right)}^{-T}\mathbf{n}_{F_{K_{\text{ref}}}^{\nu}}$, $F_{K_{\text{ref}}}^{\nu}$, $\nu=1,2,3$, are the edges of the reference cell $K_{\text{ref}}$, the corresponding normals and lengths of the edges are $\mathbf{n}_{F_{K_{\text{ref}}}^{\nu}}$ as well as $\ell_{F_{K_{\text{ref}}}^{\nu}}$. Moreover, we denote by $u_{\nu,\beta}^{*,\text{int}_{K_{\text{ref}}}}$ as well as $u_{\nu,\beta}^{*,\text{ext}_{K_{\text{ref}}}}$ the values of the mapped ALE-DG solution $u_{h}^{*}:=u_{h}\circ\bm{\chi}_{K\left(t\right)}$ evaluated in the $\beta$-th Gauss quadrature point on the edge $F_{K_{\text{ref}}}^{\nu}$. The corresponding quadrature weights for the interval $\left[-\frac{1}{2},\frac{1}{2}\right]$ are $\sigma_{\beta}$. Likewise, $\bm{\omega}_{\nu,\beta}$ are the values of the grid velocity $\bm{\omega}$ evaluated in the $\beta$-th Gauss quadrature point on the edge $F_{K_{\text{ref}}}^{\nu}$. We obtain for a constant $c\in\mathbb{R}$

[TABLE]

since the edge integrals are approximated by a $2k+1$ accurate Gauss quadrature formula and the grid velocity belongs to the space $P^{1}\left(K_{\text{ref}},\mathbb{R}^{d}\right)$ for all $t\in\left[t_{n},t_{n+1}\right]$ according to Lemma 2.1. Thus, it follows for a constant $c\in\mathbb{R}$

[TABLE]

since the numerical flux has the property (P1). Therefore, the statement of Theorem 4.3 stays true when we approximate the edge integrals in (4.3) by a $2k+1$ accurate Gauss quadrature formula.

The cell average values need to be decomposed by a quadrature formula, which includes the Gauss quadrature points for the edges $F_{K_{\text{ref}}}^{\nu}$, $\nu=1,2,3$. A quadrature formula with this property has been developed by Zhang et al. in [42]. Let us assume that $N$ is the smallest integer with $2N-3\geq k$. The $3(N-1)(k+1)$-point quadrature formula for triangular elements in [42], has the properties:

• •

The quadrature formula is exact for the integration of polynomials with degree $k$ on a triangular element.

• •

The quadrature points include the Gauss quadrature points for the edges $F_{K_{\text{ref}}}^{\nu}$, $\nu=1,2,3$.

• •

All the quadrature weights are positive and the weights for the Gauss quadrature points are given by

[TABLE]

where $\tilde{\sigma}$ corresponds to the 1-st and N-th Gauss-Lobatto quadrature weights for the interval $\left[-\frac{1}{2},\frac{1}{2}\right]$.

• •

The quadrature formula has $L:=3(N-2)(k+1)$ points in the interior of a triangular element.

This quadrature formula ensures that the cell average values can be written as follows

[TABLE]

We denote by $u_{\gamma}^{*,\text{int}_{K_{\text{ref}}}}$, $\gamma=1,...,L$, the values of the mapped ALE-DG solution $u_{h}^{*}:=u_{h}\circ\bm{\chi}_{K\left(t\right)}$ evaluated in the quadrature points which are lying in the interior of the reference cell $K_{\text{ref}}$. The corresponding quadrature weights are denoted by $\widetilde{\sigma}_{\gamma}$.

4.3.4 The maximum principle for the TVD-RK2 ALE-DG method

In this section, we prove that the TVD-RK2 ALE-DG method satisfies the maximum principle, when the ALE-DG solution is revised by the bound preserving limiter (4.16).

First of all, we show the property (L1). Therefore, similar to the process in [42], we decompose the stages of the scheme (4.18) in a sum of monotone increasing functions, which preserve constant states. Therefore, it is convenient to use vector notations. In particular, we define the set

[TABLE]

with $m$ and $M$ given by (1.2) and apply for any cell $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$ and all $\nu=1,2,3$ the vector notations

[TABLE]

Then, by applying the decomposition (4.25) of the cell average values, the scheme (4.18) can be written as

[TABLE]

where for all ${\bf a}\in\left[m,M\right]^{L}$ and ${\bf b},{\bf c}\in\mathbb{M}$

[TABLE]

In the Section 4.3.3, we mentioned that the result in Theorem 4.3 also holds when the edge integrals in (4.3) are approximated by a $2k+1$ accurate Gauss quadrature formula. Therefore for a vector $\left({\bf c},{\bf c},{\bf c}\right)\in\left[m,M\right]^{L}\cup\mathbb{M}$ with vector components given by the constant $c\in\left[m,M\right]$, we obtain by the Theorem 4.3

[TABLE]

[TABLE]

Henceforth, for the sake of simplicity, we will use the Lax-Friedrichs flux (2.25) for the analysis. However, it should be noted that the techniques which are presented in this section can be also applied to any other monotone flux. The use of the Lax-Friedrichs flux ensures to prove the following lemmas.

Lemma 4.5**.**

Let $\left|K^{n}\right|$, $\left|K^{n,1}\right|$, $\left|K^{n+1}\right|$ be given by (4.17) and ${\bf a}\in\left[m,M\right]^{L}$, ${\bf b},{\bf c}\in\mathbb{M}$. Then, under the CFL constraint

[TABLE]

with $\lambda^{n}$ given by (2.25c), it holds

[TABLE]

Proof.

First of all, we apply the Lax-Friedrichs flux (2.25) and rewrite the functions (4.28c), (4.28d) and (4.28e) as follows

[TABLE]

where $\tilde{\mathbf{n}}_{F_{K_{\text{ref}}}^{\nu}}\left(t\right):={\bf A}_{K\left(t\right)}^{-T}\mathbf{n}_{F_{K_{\text{ref}}}^{\nu}}$, $\nu=1,2,3$.

Next, we observe that for all $t\in\left[0,T\right]$, $\nu=1,2,3$ and $\beta=1,...,k+1$

[TABLE]

where $\lambda^{n}$ is given by (2.25c). Hence, we obtain by (4.17a) as well as the CFL constraint (4.31) for $\beta=1,...,k+1$

[TABLE]

and similar we get

[TABLE]

[TABLE]

In the following, we highlight by the symbolic notation $\uparrow$ that a function is increasing in the marked arguments. Likewise, we apply the notation $\bm{\uparrow}$ to highlight that a function with vector arguments increases in each vector component. Then, it follows by (4.33), (4.34) and (4.35)

[TABLE]

since $\widehat{g}_{\pm}\left(\bm{\omega},\uparrow,J_{K\left(t\right)}\tilde{\mathbf{n}}\left(t\right)\right)$ for all $t\in\left[t_{n},t_{n+1}\right]$ and all cells $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$. Therefore, we obtain

[TABLE]

since for all $\gamma=1,...,L$, it follows by (4.17a) and the CFL constraint (4.31)

[TABLE]

Finally, the equation (4.29) and (4.36) provide the inequality (4.32). ∎

Lemma 4.6**.**

Let $\left|K^{n}\right|$, $\left|K^{n,1}\right|$, $\left|K^{n+1}\right|$ be given by (4.17) and ${\bf a},\widetilde{{\bf a}}\in\left[m,M\right]^{L}$, ${\bf b},\widetilde{{\bf b}},\widetilde{{\bf c}}\in\mathbb{M}$. Then, under the CFL constraint

[TABLE]

with $\lambda^{n}$ given by (2.25c), we have

[TABLE]

Proof.

The equations (4.17a) and (4.17b) supply

[TABLE]

Thus, the same procedure as in the proof of lemma 4.6 provides

[TABLE]

by applying the identity (4.39) and the CFL constraint (4.37). Furthermore, by (4.17a) and (4.17b) follows

[TABLE]

Hence, we obtain for all $\gamma=1,...,L$

[TABLE]

by the identity (4.41) and the CFL constraint (4.37). In a similar way, it follows for all $\nu=1,2,3$ as well as $\beta=1,...,k+1$

[TABLE]

This ensures

[TABLE]

Therefore, the equation (4.30), (4.40) and (4.42) supply the inequality (4.38). ∎

We note that the assumption ($A3$) for the mesh parameter provides for all $t\in\left[0,T\right]$ and all cells $K\left(t\right)\in\mathcal{T}_{\left(t\right)}$

[TABLE]

since $\rho_{K\left(t\right)}$ denotes the radius of the largest ball contained in the cell $K\left(t\right)$. Therefore, the CFL constraints (4.31) and (4.37) can be generalized as

[TABLE]

Now, we apply the generalized CFL constraint and the previous lemmas to prove the maximum principle for the second order fully-discrete ALE-DG method.

Theorem 4.7**.**

Suppose ${\bf\widetilde{u}}^{n,*,\text{int}_{K_{\text{ref}}}}\in\left[m,M\right]^{L},{\bf u}^{n,*,\text{int}_{K_{\text{ref}}}},{\bf u}^{n,*,\text{ext}_{K_{\text{ref}}}}\in\mathbb{M}$. Furthermore, ($A1$) - ($A3$) and the CFL constraint (4.43) are satisfied. Then, the solution $u_{h}^{n+1}$ of the second order fully-discrete ALE-DG method revised by Zhang, Xia and Shu’s bound preserving limiter (4.16) belongs to the interval $\left[m,M\right]$.

Proof.

Let $K^{n}\in\mathcal{T}_{\left(t_{n}\right)}$ be an arbitrary cell. The main part of the proof follows in two steps.

Step 1. Since, ${\bf\widetilde{u}}^{n,*,\text{int}_{K_{\text{ref}}}}\in\left[m,M\right]^{L},{\bf u}^{n,*,\text{int}_{K_{\text{ref}}}},{\bf u}^{n,*,\text{ext}_{K_{\text{ref}}}}\in\mathbb{M}$ and the CFL constraint (4.43) is satisfied, it follows $\overline{u}_{K^{n,1}}^{n,1}\in\left[m,M\right]$ by the inequality (4.32) in lemma 4.5 and the equation (4.27a). Hence, the bound preserving limiter (4.16) ensures $u_{h}^{n,1}|_{K^{n,1}}\in\left[m,M\right]$.

Step 2. In the first step, it has been shown that $u_{h}^{n,1}|_{K^{n,1}}\in\left[m,M\right]$ when the bound preserving limiter (4.16) was applied. Hence, ${\bf\widetilde{u}}^{n,1*,\text{int}_{K_{\text{ref}}}}\in\left[m,M\right]^{L}$, ${\bf u}^{n,1,*,\text{int}_{K_{\text{ref}}}},{\bf u}^{n,1,*,\text{ext}_{K_{\text{ref}}}}\in\mathbb{M}$. Thus, the inequality (4.38) in lemma 4.6 and the equation (4.27b) supply $\overline{u}_{K^{n+1}}^{n+1}\in\left[m,M\right]$, since the CFL constraint (4.43) is satisfied. Finally, the bound preserving limiter (4.16) ensures that $u_{h}^{n+1}|_{K^{n+1}}\in\left[m,M\right]$.

In a similar way, we proceed for any other cell in $\mathcal{T}_{\left(t_{n+1}\right)}$. Thus, it follows $u_{h}^{n+1}\in\left[m,M\right]$. ∎

4.3.5 The maximum principle for the third order fully-discrete ALE-DG method

In this section, we show briefly how the result in theorem 4.7 can be extended to the third order TVD-RK3 ALE-DG method. We apply the decomposition (4.25) of the cell average values and the vector notations (4.26) to rewrite the scheme (4.20) as

[TABLE]

By the Theorem 4.3 we obtain the identities

[TABLE]

[TABLE]

[TABLE]

where $\left({\bf c},{\bf c},{\bf c}\right)\in\left[m,M\right]^{L}\cup\mathbb{M}$ with vector components given by the constant $c\in\left[m,M\right]$. Moreover, like in the previous section, we apply the Lax-Friedrichs flux (2.25). Then, by the same argumentation as in the proof of the lemma 4.5 and lemma 4.6, we obtain the following lemma.

Lemma 4.8**.**

Let $\left|K^{n}\right|$, $\left|K^{n,1}\right|$, $\left|K^{n,2}\right|$, $\left|K^{n+1}\right|$ be given by (4.19) and ${\bf a},{\bf\widetilde{a}}\in\left[m,M\right]^{L}$, ${\bf b},{\bf c},{\bf\widetilde{b}},{\bf\widetilde{c}}\in\mathbb{M}$. Then, under the CFL constraint (4.43), it holds

[TABLE]

[TABLE]

[TABLE]

This lemma provide the following analogue of the Theorem 4.7 for the second order fully-discrete ALE-DG method. The proof follows similar to the proof of Theorem 4.7. Hence, it is skipped.

Theorem 4.9**.**

Suppose ${\bf\widetilde{u}}^{n,*,\text{int}_{K_{\text{ref}}}}\in\left[m,M\right]^{L}$, ${\bf u}^{n,*,\text{int}_{K_{\text{ref}}}},{\bf u}^{n,*,\text{ext}_{K_{\text{ref}}}}\in\mathbb{M}$. Furthermore, ($A1$) - ($A3$) and the CFL constraint (4.43) are satisfied. Then, the solution $u_{h}^{n+1}$ of the third order fully-discrete ALE-DG method revised by Zhang, Xia and Shu’s bound preserving limiter (4.16) belongs to the interval $\left[m,M\right]$.

Remark 3**.**

It is also possible to apply other TVD-RK methods like the five stage fourth order method of Spiteri and Ruuth [33] as time integrator in a fully-discrete ALE-DG method. Then, it can be proven by the techniques presented in this section that these fully-discrete ALE-DG methods satisfy the maximum principle.

5 Numerical experiments

In this section, we demonstrate the performance of the ALE-DG method for conservation laws in two dimensions. In our simulation, the criss-triangular meshes are used. Furthermore, the third order TVD Runge-Kutta method is used in the first example (Example 5.1). In the other examples the five stage fourth order TVD Runge-Kutta method of Spiteri and Ruuth [33] is used for the time discretization. We observe that for this high order approximation in time our theoretical results hold numerically, too. In order to avoid complications with the stability of the explicit time integrator, we apply the suitable CFL condition dependent on equation (4.43).

To verify our theoretical results, we present numerical simulations for a linear advection equation and Burgers’ equation. Moreover, to highlight that the ALE-DG method can be also used for systems of conservation laws, we present a plain wave problem and a smooth vortex problem for the compressible Euler equations with a polytropic gas.

In all the numerical simulations, we consider two moving mesh scenarios. First a static uniform criss-triangular mesh with cell size $h_{0}$ is used. Next, a moving mesh with the grid point distribution

[TABLE]

is used. In equation (5.1) the points $(x_{j},y_{j})$ are the vertices of the triangular mesh and $t_{0}=\sqrt{10^{2}+5^{2}}$. The vertices $(x_{j}(0),y_{j}(0))$ at initial time are given by the same mesh as in the first moving mesh scenario. The grid point distribution (5.1) has also been used by Klingenberg [15] et al. and Persson et al. [29]. As an example, in Fig. 5.2 we draw a typical mesh at $t=0$ and the deformed one at $t=1$.

Example 5.1

(Linear advection equation) Here, we test the linear equation

[TABLE]

with the periodic boundary condition and the initial condition is taken to be $u_{0}(x,y)=1+0.5\sin(\pi(x+y))$. The exact solution is $u(x,y,t)=u_{0}(x-t,y-t)$ at time $t$.

The numerical solution on the static uniform grid is $u_{h}^{S}$ and the moving mesh solution is $u_{h}^{M}$. In Table 5.3, we show the $\mathrm{L}^{2}$-errors and the rates of convergence of the numerical solutions $u_{h}^{S}$ and $u_{h}^{M}$ for the advection equation at time $t=1$. In the computation, we used piecewise $P^{k}$ polynomial spaces with $k=1,2,3$ on static and moving triangular meshes with cell size $h_{0}$. We observe that both $u_{h}^{S}$ and $u_{h}^{M}$ have the optimal accuracy when $P^{k}$ polynomial spaces with $k=1,2,3$ are applied, but they do not satisfy the discrete maximum principle. This was expected, since the bound preserving limiter in Section 4.3.1 was not applied in the test case.

To verify the maximum principle, we run the same test case again, but this time with the bound preserving limiter. The results are given in Table 5.4, where $\widetilde{u}_{h}^{M}$ is the ALE-DG solution on the moving mesh with bound preserving limiter. We observe that the $\mathrm{L}^{2}$-errors and the rates of convergence for the moving mesh solution $\widetilde{u}_{h}^{M}$ are not affected by the use of the bound preserving limiter and thus the method still has optimal accuracy. Furthermore, we can see that the numerical solution $\widetilde{u}_{h}^{M}$ is limited in the same range $[0.5,1.5]$ as the initial data.

Example 5.2

(Burgers’ equation) Next, we investigate the Burgers’ equation

[TABLE]

In our simulation, the periodic boundary condition is used and the initial condition is also taken to be $u_{0}(x,y)=1+0.5\sin(\pi(x+y))$. We compute this example up to time $t=0.1$ before the shock front has been developed in the numerical solution.

In Table 5.5, the $\mathrm{L}^{2}$-errors and the rates of convergence for the numerical solutions $u_{h}^{S}$ and $u_{h}^{M}$, $\tilde{u}_{h}^{M}$ are presented. These functions are computed by the ALE-DG method with $P^{k}$, $k=1,2,3$, polynomial spaces. As in the previous example $u^{S}_{h}$ and $u^{M}_{h}$ are the numerical solutions of the ALE-DG method on the static uniform mesh and on the moving mesh. The solution $\tilde{u}_{h}^{M}$ is the moving mesh ALE-DG solution revised by the bound preserving limiter. We observe that the optimal accuracy is obtained for $u_{h}^{M},\widetilde{u}_{h}^{M}$ and $u_{h}^{S}$ in both moving mesh scenarios. Furthermore, maximum and minimum values of $\widetilde{u}_{h}^{M}$ are limited in the same range $[0.5,1.5]$ as the initial data when the bound preserving limiter is applied in the ALE-DG method on the moving grid.

To show that our proposed ALE-DG methods can handle the problem with shocks, we show the numerical solutions $u_{h}^{S}$ and $u_{h}^{M}$ of Burgers’ equation at time $t=0.45$ with piecewise $P^{1}$ polynomial approximation in Fig. 5.3. Here, we use the slope limiter developed by Cockburn et al. in [5]. From the results, it can be seen that the ALE-DG methods can capture the shocks well for the Burgers’ equation on both static and moving meshes.

Example 5.3

(Compressible Euler equations) We consider the two dimensional compressible Euler equations of gas dynamics for a polytropic gas

[TABLE]

with

[TABLE]

Here, $\rho$ is the density, $\mathbf{u}=\left(u,v\right)^{T}$ is the velocity field and $E$ is the total energy. Moreover, the adiabatic constant of air $\gamma=1.4$ is used, the pressure is given by $p=(\gamma-1)\left(E-\frac{1}{2}\rho\left|{\bf u}\right|^{2}\right)$ and ${\bf I}$ is the identity matrix. In our simulation, we test a plain wave problem and a smooth vortex problem.

First, we consider the plain wave problem and choose the domain related parameter in (5.4) as $x_{l}=y_{l}=0$ and $x_{r}=y_{r}=2$. The problem has the initial data

[TABLE]

and is investigated with the periodic boundary condition. The results in Table 5.6 show the $\mathrm{L}^{2}$-errors and the rates of convergence of the density $\rho_{h}$ given by the ALE-DG method with $P^{k}$, $k=1,2,3$, polynomial spaces. The numerical solutions of the ALE-DG method on the static uniform mesh and the moving mesh are $\rho^{S}_{h}$ and $\rho^{M}_{h}$. The numerical results show that we can obtain the optimal accuracy on both meshes.

Next, we consider the smooth vortex problem and choose the domain related parameter in (5.4) as $x_{l}=y_{l}=0$, $x_{r}=20$ and $y_{r}=15$. This problem was also presented by Persson et al. [29] and the initial condition is

[TABLE]

where $(\rho_{0},u_{0},v_{0},p_{0})^{T}=(1,1,1,1)^{T}$, $\theta=\arctan(0.5)$, $(x_{0},y_{0})=(5,5)$, $\epsilon=0.3$, $r_{0}=1.5$, $r=(1-(x-x_{0})^{2}-(y-y_{0})^{2})/r_{0}^{2}$ and $\alpha=\frac{(\gamma-1)\epsilon^{2}}{8\gamma\pi^{2}}$. We test the problem up to time $t=\sqrt{10^{2}+5^{2}}$ with the Direchlet boundary condition. The ALE-DG method with $P^{k},k=1,2,3$ approxiation is used to solve the problem on static uniform triangular meshes with cell size $h_{0}$ and on moving meshes with the grid point distribution (5.1) as before. The $\mathrm{L}^{2}$-errors and the rates of convergence for the numerical solutions of the density $\rho_{h}^{S},\rho_{h}^{M}$ and the pressure $p_{h}^{M}$, $p_{h}^{S}$ are shown in Table 5.7. We see that the numerical solutions are optimally accurate in both moving mesh scenarios.

Example 5.4

(Constant state preservation) The previous examples show that the ALE-DG method on moving meshes maintains the high order accuracy as the DG method on static meshes. The ability of the ALE-DG method to preserve constant states needs to be investigated, too. For this reason the linear advection equation (5.2) and the Burgers’ equation (5.3) are considered with the constant initial condition $u_{0}=1$. We solve these initial value problems with the ALE-DG method on moving triangular meshes with the grid point distribution (5.1). In Table 5.8 the results of the computations are listed and it can be seen that the ALE-DG method numerically satisfies the GCL. This result was expected, since we used a time discretization with an order greater than two and in the Section 4.2 it has been proven that a time discretization of this type is enough to ensure that the method preserves constant states.

To further show the D-GCL of ALE-DG methods, we adopted the meshes $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$ of $32$ cells in Fig. 5.4 with recursive refinement as the initial and final meshes. For the linear advection equation with the constant initial condition $u_{0}(x)=1$, we show the D-GCL errors at time $T=1.0$ in Table 5.9 by forward Euler, TVD-RK2 and TVD-RK3 methods respectively. Here, $P^{1}$ piecewise polynomial space is used in the ALE-DG method. We take time step size $\triangle t=\frac{h_{0}}{\max(|\bm{\omega}|)}$ with $\bm{\omega}=(\frac{(\mathbf{x_{1}}-\mathbf{x_{2}})}{T},\frac{(\mathbf{y_{1}}-\mathbf{y_{2}})}{T})$ and $(\mathbf{x_{1}},\mathbf{y_{1}}),(\mathbf{x_{2}},\mathbf{y_{2}})$ are vertices of meshes $\mathcal{T}_{1}$ and $\mathcal{T}_{2}$. The numerical results are consistent with the analysis on the D-GCL of ALE-DG methods.

6 Conclusions

In this paper, an ALE-DG method to solve conservation laws in several space dimensions on moving simplex meshes has been developed and analyzed. We began the paper with an analysis of the semi-discrete ALE-DG method and proved the $\mathrm{L}^{2}$-stability. Moreover, we presented a suboptimal a priori error estimate with respect to the $\mathrm{L}^{\infty}\left(0,T;\mathrm{L}^{2}\left(\Omega\right)\right)$-norm, where the suboptimality refers to the approximation properties of the discrete space.

Afterward, the fully-discrete ALE-DG method was investigated. In the context of Total-variation-diminishing Runge-Kutta methods, a relationship between the spatial dimension and the discrete geometric conservation law was elaborated. Furthermore, in two dimensions, second and third order fully-discrete ALE-DG methods were presented. We proved that these methods satisfy the maximum principle when the bound preserving limiter developed by Zhang, Xia and Shu [42] is applied. In a future work, it would be worthwhile to investigate if these methods are positive preserving when they are applied to the compressible Euler equations.

Beside our theoretical investigations, several numerical test examples for two moving mesh scenarios have been presented. These examples support our theoretical results and show that the ALE-DG method is numerically stable and uniformly high order accurate. In particular, the two test examples for the compressible Euler equations support the expectation that the ALE-DG method can be also applied to systems of conservation laws even when the development and analysis in this paper has been focused on scalar conservation laws in several space dimensions.

It should be mentioned that we did not use a moving mesh methodology in the numerical examples. The grid point distribution was specified for the calculations. The development of a suitable moving mesh methodology for our ALE-DG method is also a project for a future work.