High Order Explicit Local Time-Stepping Methods For Hyperbolic Conservation Laws

TL;DR

This paper develops high order explicit local time-stepping methods for hyperbolic conservation laws, enabling variable time steps across regions while ensuring stability and conservation, demonstrated through numerical experiments.

Contribution

It introduces a novel predictor-corrector LTS framework based on SSP-RK schemes for hyperbolic conservation laws, achieving up to fourth-order accuracy.

Findings

Methods are proven to be conservative and stable for scalar laws.

Numerical experiments show high accuracy and efficiency.

Framework allows spatially variable time stepping with strong stability properties.

Abstract

In this paper we present and analyze a general framework for constructing high order explicit local time stepping (LTS) methods for hyperbolic conservation laws. In particular, we consider the model problem discretized by Runge-Kutta discontinuous Galerkin (RKDG) methods and design LTS algorithms based on strong stability preserving Runge-Kutta (SSP-RK) schemes, that allow spatially variable time step sizes to be used for time integrations in different regions. The proposed algorithms are of predictor-corrector type, in which the interface information along the time direction is first predicted based on the SSP-RK approximations and Taylor expansions, and then the fluxes over the region of interface are corrected to conserve mass exactly at each time step. Following the proposed framework, we detail the corresponding LTS schemes with accuracy up to the fourth order, and prove their conservation property and nonlinear stability for the scalar conservation laws. Numerical experiments are also presented to demonstrate excellent performance of the proposed LTS algorithms.

Tables

Table 1. Table 1 : [Linear advection with a smooth initial condition] L 1 superscript 𝐿 1 \displaystyle L^{1} relative errors at T = 2 𝑇 2 \displaystyle T=2 for different M 𝑀 \displaystyle M . The rates of convergence “CR” for fixed M 𝑀 \displaystyle M are shown in square brackets.

Linear advection
$\Delta x_{\text{coarse}}$	$M$	RK-DG2		RK-DG3		RK-DG4
		Rel. $L^1$ error	[CR]	Rel. $L^1$ error	[CR]	Rel. $L^1$ error	[CR]
$1/5$	$1$	5.70e-02	–	1.39e-03	–	4.43e-05	–
	$2$	3.59e-02	–	8.09e-04	–	2.40e-05	–
	$4$	3.13e-02	–	7.65e-04	–	2.32e-05	–
	$8$	3.04e-02	–	7.62e-04	–	2.33e-05	–
$1/10$	$1$	1.36e-02	[2.07]	1.66e-04	[3.07]	2.73e-06	[4.02]
	$2$	8.71e-03	[2.04]	9.76e-05	[3.05]	1.46e-06	[4.04]
	$4$	7.60e-03	[2.04]	9.20e-05	[3.06]	1.39e-06	[4.06]
	$8$	3.91e-03	[2.04]	9.17e-05	[3.06]	1.41e-06	[4.05]
$1/20$	$1$	3.30e-03	[2.04]	2.04e-05	[3.03]	1.71e-07	[4.00]
	$2$	2.14e-03	[2.03]	1.20e-05	[3.02]	9.07e-08	[4.01]
	$4$	1.87e-03	[2.02]	1.14e-05	[3.01]	8.66e-08	[4.01]
	$8$	1.81e-03	[2.03]	1.13e-05	[3.02]	8.83e-08	[4.00]
$1/40$	$1$	8.11e-04	[2.03]	2.54e-06	[3.01]	1.07e-08	[4.00]
	$2$	5.31e-04	[2.01]	1.49e-06	[3.01]	5.67e-09	[4.00]
	$4$	4.64e-04	[2.01]	1.41e-06	[3.02]	5.42e-09	[4.00]
	$8$	4.50e-04	[2.01]	1.41e-06	[3.00]	5.53e-09	[4.00]

Table 2. Table 2 : [Linear advection with a discontinuous initial condition] L 1 superscript 𝐿 1 \displaystyle L^{1} relative errors at T = 0.5 𝑇 0.5 \displaystyle T=0.5 and T = 1 𝑇 1 \displaystyle T=1 of RK-DG2 global time-stepping (GTS) and local time-stepping (LTS) schemes.

Linear advection, RK-DG2
$\Delta x_{\text{coarse}}$	$M$	At $T=0.5$			At $T=1$
		coarse GTS	fine GTS	LTS	coarse GTS	fine GTS	LTS
1/40	1	2.58e-02			2.53e-02
	2	3.46e-02	1.52e-02		3.46e-02	1.50e-02
	4	–	9.00e-03		–	8.99e-03
	8	–	5.40e-03		–	5.46e-03

Table 3. Table 3 : [Burger’s equation] L 1 superscript 𝐿 1 \displaystyle L^{1} relative errors at T = 0.3 𝑇 0.3 \displaystyle T=0.3 for different M 𝑀 \displaystyle M . The rates of convergence “CR” for fixed M 𝑀 \displaystyle M are shown in square brackets.

$\Delta x_{\text{coarse}}$	$M$	RK-DG2		RK-DG3		RK-DG4
		Rel. $L^1$ error	[CR]	Rel. $L^1$ error	[CR]	Rel. $L^1$ error	[CR]
$1/10$	$1$	6.49e-03	[1.97]	3.27e-04	[3.18]	2.26e-05	–
	$2$	3.97e-03	[2.07]	1.10e-04	[3.12]	6.75e-06	–
	$4$	3.37e-03	[2.09]	8.97e-05	[3.26]	5.97e-06	–
	$8$	3.22e-03	[2.09]	8.73e-05	[3.28]	5.93e-06	–
$1/20$	$1$	1.62e-03	[2.00]	4.10e-05	[3.00]	1.31e-06	[4.11]
	$2$	9.63e-04	[2.04]	1.45e-05	[2.92]	4.19e-07	[4.01]
	$4$	8.11e-04	[2.06]	1.17e-05	[2.94]	3.60e-07	[4.05]
	$8$	7.74e-04	[2.06]	1.14e-05	[2.94]	3.57e-07	[4.05]
$1/40$	$1$	4.03e-04	[2.01]	5.01e-06	[3.03]	8.68e-08	[3.92]
	$2$	2.38e-04	[2.02]	1.81e-06	[3.00]	2.63e-08	[3.99]
	$4$	1.99e-04	[2.03]	1.43e-06	[3.03]	2.30e-08	[3.97]
	$8$	1.89e-04	[2.03]	1.40e-06	[3.03]	2.28e-08	[3.97]
$1/80$	$1$	1.00e-04	[2.01]	6.18e-07	[3.02]	5.44e-09	[4.00]
	$2$	5.95e-05	[2.00]	2.25e-07	[3.01]	1.75e-09	[3.91]
	$4$	4.95e-05	[2.01]	1.78e-07	[3.01]	1.56e-09	[3.88]
	$8$	4.71e-05	[2.01]	1.74e-07	[3.01]	1.55e-09	[3.88]

Table 4. Table 4 : [Burger’s equation] L 1 superscript 𝐿 1 \displaystyle L^{1} relative errors in smooth regions at T = 1.1 𝑇 1.1 \displaystyle T=1.1 of different local time-stepping schemes.

At $T=1.1$, errors in smooth region $\Vert x-\text{shock}\Vert >=0.1$
$\Delta x_{\text{coarse}}$	$M$	2nd order LTS	3rd order LTS	4th order LTS
1/40	1	6.16e-05	4.25e-07	4.31e-09
	2	4.67e-05	1.73e-07	1.71e-09
	4	4.37e-05	1.53e-07	1.59e-09
	8	4.35e-05	1.52e-07	1.60e-09

Table 5. Table 5 : [Sod shock tube problem] L 1 superscript 𝐿 1 \displaystyle L^{1} relative errors at T = 2.0 𝑇 2.0 \displaystyle T=2.0 of the second order LTS algorithm for the Sod shock tube problem.

$\Delta x_{\text{coarse}}$	M	Density		Velocity		Pressure
		2 domains	3 domains	2 domains	3 domains	2 domains	3 domains
1/5	1	1.79e-02		3.88e-02		1.76e-02
	2	8.96e-03	8.28e-03	2.03e-02	1.75e-02	8.46e-03	7.22e-03
	4	6.35e-03	4.25e-03	1.39e-02	8.82e-03	6.82e-03	3.62e-03
	8	6.06e-03	2.22e-03	1.27e-02	4.46e-03	7.34e-03	1.81e-03
1/10	1	8.28e-03		1.78e-02		7.24e-03
	2	4.41e-03	4.23e-03	9.79e-03	8.83e-03	3.89e-03	3.61e-03
	4	2.88e-03	2.21e-03	6.11e-03	4.46e-03	2.76e-03	1.80e-03
	8	2.50e-03	1.16e-03	5.19e-03	2.30e-03	2.80e-03	9.03e-04

Table 6. Table 6 : [Lax shock tube problem] L 1 superscript 𝐿 1 \displaystyle L^{1} relative errors at T = 1.3 𝑇 1.3 \displaystyle T=1.3 of the second order LTS algorithm.

$\Delta x_{\text{coarse}}$	M	Density	Velocity	Pressure
1/5	1	4.74e-02	2.15e-02	1.44e-02
	2	2.94e-02	1.69e-02	1.07e-02
	4	2.11e-02	1.73e-02	9.12e-03
	8	1.18e-02	1.37e-02	7.63e-03
1/10	1	2.52e-02	1.20e-02	6.55e-03
	2	1.51e-02	9.07e-03	4.74e-03
	4	8.88e-03	7.13e-03	3.91e-03
	8	5.70e-03	7.41e-03	4.19e-03

Table 7. Table 7 : [ 30 ] Coefficients of the SSP-RK (5,4) scheme with 𝒞 = 1.652 𝒞 1.652 \displaystyle\mathcal{C}=1.652 .

${\alpha }_{i\nu }$
1	0	0	0	0
0.261216512493821	0.738783487506179	0	0	0
0.623613752757655	0	0.376386247242345	0	0
0.444745181201454	0.120932584902288	0	0.434322233896258	0
0.213357715199957	0.209928473023448	0.063353148180384	0	0.513360663596212
${\beta }_{i\nu }$
0.605491839566400	0	0	0	0
0	0.447327372891397	0	0	0
0.000000844149769	0	0.227898801230261	0	0
0.002856233144485	0.073223693296006	0	0.262978568366434	0
0.002362549760441	0.127109977308333	0.038359814234063	0	0.310835692561898

Equations

$$\begin{array}[]{rll}\frac{\partial\boldsymbol{u}}{\partial t}+\sum_{i=1}^{d}\frac{\partial\boldsymbol{f}_{i}}{\partial x_{i}}(\boldsymbol{u})&=\boldsymbol{0},&\text{on}\;\mathbb{R}^{d}\times(0,T),\\ \boldsymbol{u}(\boldsymbol{x},0)&=\boldsymbol{u}_{0}(\boldsymbol{x}),&\text{in}\;\mathbb{R}^{d},\end{array}\vspace{-0.2cm}$$

$$\begin{array}[]{rcl}\boldsymbol{u}&\mapsto&\boldsymbol{f}_{i}(\boldsymbol{u}):=\left(f_{i1}(\boldsymbol{u}),\ldots,f_{im}(\boldsymbol{u})\right),\end{array}$$

$$\begin{array}[]{rll}\frac{\partial u}{\partial t}+\frac{\partial}{\partial x}f(u)&=0,&\text{on}\;\mathbb{R}\times(0,T),\vspace{4pt}\\ u(x,0)&=\;u_{0},&\text{in}\;\mathbb{R}.\end{array}$$

$$V_h=V_h^k=\left\{v\in L^1(\mathbb{R}):v{\mid }_{I_j}\in {\mathcal{P}}^k(I_j),\forall j\right\}\not\subset H^1(\mathbb{R}),$$

$$\begin{array}[]{l}\int_{I_{j}}\partial_{t}u_{h}(x,t)\,v_{h}(x)\,dx-\int_{I_{j}}f(u_{h}(x,t))\,\partial_{x}v_{h}(x)\,dx\vspace{0.1cm}\\ \hskip 42.67912pt+h(u_{h})_{j+\nicefrac{{1}}{{2}}}(t)\,v_{h}(x_{j+\nicefrac{{1}}{{2}}}^{-})-h(u_{h})_{j-\nicefrac{{1}}{{2}}}(t)\,v_{h}(x_{j-\nicefrac{{1}}{{2}}}^{+})=0,\vspace{0.2cm}\\ \int_{I_{j}}u_{h}(x,0)v_{h}(x)\,dx=\int_{I_{j}}u_{0}(x)v_{h}(x)\,dx.\end{array}$$

$$h(u{)}_{j+\text{nicefrac}12}(t)=h\left(u(x_{j+\text{nicefrac}12}^{-},t),u(x_{j+\text{nicefrac}12}^{+},t)\right).$$

$${\phi }_j^{(l)}:=P_l\left(\frac{2(x-x_j)}{{\Delta }_j}\right),\text{for }l=0,1,\dots ,k,$$

$$u_h(x,t)=\sum _{l=0}^k{u}_j^{(l)}(t){\phi }_j^{(l)}(x),\text{for}x\in I_j,$$

$$u_j^{(l)}(t):=\frac{2l+1}{{\Delta }_j}{\int }_{I_j}u(x,t){\phi }_j^{(l)}(x)dx,\text{for }l=0,1,\dots ,k.$$

$$\begin{array}[]{l}\left(\frac{1}{2l+1}\right)\frac{du_{j}^{(l)}(t)}{dt}-\frac{1}{\Delta_{j}}\int_{I_{j}}f(u_{h}(x,t))\,\partial_{x}\varphi_{j}^{(l)}\,dx\vspace{0.1cm}\\ \hskip 28.45274pt+\frac{1}{\Delta_{j}}\left[h(u_{h})_{j+\nicefrac{{1}}{{2}}}(t)-(-1)^{l}h(u_{h})_{j-\nicefrac{{1}}{{2}}}(t)\right]=0,\quad\forall\,l=0,1,\ldots,k,\end{array}$$

$$u_j^{(l)}(0)=\frac{2l+1}{{\Delta }_j}{\int }_{I_j}{u}_0(x){\phi }_j^{(l)}(x)dx.$$

$${\phi }_j^{(l)}(x_{j+\text{nicefrac}12}^{-})=P_l(1)=1,{\phi }_j^{(l)}(x_{j-\text{nicefrac}12}^{+})=P_l(-1)=(-1{)}^l.$$

$$u_{j+\text{nicefrac}12}^{-}(t)=\sum _{l=0}^k{u}_j^{(l)},u_{j-\text{nicefrac}12}^{+}(t)=\sum _{l=0}^k(-1{)}^l{u}_j^{(l)}.$$

$${\int }_{I_j}f(u_h){\partial }_x{\phi }_j^{(l)}dx={\int }_{I_j}f\left(u_{j-\text{nicefrac}12}^{+},(u_j^{(l)}{)}_{0\le l\le k},u_{j+\text{nicefrac}12}^{-}\right){\partial }_x{\phi }_j^{(l)}dx.\text{vspace}-0.2cm$$

$$\frac{d{\mathbf{U}}_h}{dt}={\mathbf{L}}_h({\mathbf{U}}_h),{\mathbf{U}}_h(0)={\mathbf{U}}_{h0},$$

$${\mathbf{L}}_h({\mathbf{U}}_h)={\left(L_{h,j}^{(l)}(u_{j-\text{nicefrac}12}^{\pm },{\mathbf{u}}_j,u_{j+\text{nicefrac}12}^{\pm })\right)}_{\forall j,l=0,1,\dots ,k},\text{vspace}-0.2cm$$

$$\begin{array}[]{ll}L_{h,j}^{(l)}(u_{j-\nicefrac{{1}}{{2}}}^{\pm},\boldsymbol{u}_{j},u_{j+\nicefrac{{1}}{{2}}}^{\pm})=&\frac{2l+1}{\Delta_{j}}\bigg{\{}\int_{I_{j}}f\left(u_{j-\nicefrac{{1}}{{2}}}^{+},\boldsymbol{u}_{j},u_{j+\nicefrac{{1}}{{2}}}^{-}\right)\,\partial_{x}\varphi_{j}^{(l)}\,dx\vspace{3pt}\\ &-\left[h(u_{j+\nicefrac{{1}}{{2}}}^{-},u_{j+\nicefrac{{1}}{{2}}}^{+})-(-1)^{l}h(u_{j-\nicefrac{{1}}{{2}}}^{-},u_{j-\nicefrac{{1}}{{2}}}^{+})\right]\bigg{\}},\end{array}$$

$${\mathbf{U}}_h^{n,(i)}=\sum _{\nu =0}^{i-1}{\alpha }_{i\nu }{\mathbf{U}}_h^{n,(\nu )}+{\beta }_{i\nu }\Delta t{\mathbf{L}}_h({\mathbf{U}}_h^{n,(\nu )}),\forall i=1,\dots ,s,$$

$$TV({\mathbf{u}}^{(l),n})=\sum _j\left|u_{j+1}^{(l),n}-u_j^{(l),n}\right|,\forall l=0,\dots ,k,\text{and}n=0,\dots ,N-1.$$

$$TV({\mathbf{u}}^{(l),n+1})\le TV({\mathbf{u}}^{(l),n}),\forall l=0,\dots ,k\text{and}n=0,\dots ,N-1,$$

$$TV({\mathbf{u}}^{(l),n+1})\le TV({\mathbf{u}}^{(l),0})+BT,\forall l=0,\dots ,k\text{and}n=0,\dots ,N-1,$$

$$\widetilde{m}(a_{1},\ldots,a_{\nu})=\left\{\begin{array}[]{ll}a_{1},&\text{if}\;|a_{1}|\leq C_{M}h^{2},\\ m(a_{1},\ldots,a_{\nu}),&\text{otherwise}.\end{array}\right.$$

$$m(a_{1},\ldots,a_{\nu})=s\min_{1\leq i\leq\nu}|a_{i}|,\;\;{\text{ with }}s=\left\{\begin{array}[]{ll}\text{sign}(a_{1}),&\text{if}\;\text{sign}(a_{1})=\ldots=\text{sign}(a_{\nu}),\\ 0,&\text{otherwise}.\end{array}\right.$$

$$\Lambda {\Pi }_h^1(u_h){|}_{I_j}=u_j^{(0)}+\tilde{m}(u_j^{(1)},u_{j+1}^{(0)}-u_j^{(0)},u_j^{(0)}-u_{j-1}^{(0)}){\phi }_j^{(1)}(x).$$

$$\begin{array}[]{ll}u_{j+\nicefrac{{1}}{{2}}}^{-\text{(mod)}}&=u_{j}^{(0)}+\widetilde{m}(u_{j+\nicefrac{{1}}{{2}}}^{-}-u_{j}^{(0)},u_{j+1}^{(0)}-u_{j}^{(0)},u_{j}^{(0)}-u_{j-1}^{(0)}),\vspace{3pt}\\ u_{j-\nicefrac{{1}}{{2}}}^{+\text{(mod)}}&=u_{j}^{(0)}-\widetilde{m}(u_{j}^{(0)}-u_{j-\nicefrac{{1}}{{2}}}^{+},u_{j+1}^{(0)}-u_{j}^{(0)},u_{j}^{(0)}-u_{j-1}^{(0)}),\end{array}\vspace{-0.1cm}$$

$$\Lambda\Pi^{k}_{h}(u_{h})|_{I_{j}}=\left\{\begin{array}[]{ll}u_{h}|_{I_{j}},&\text{if $\displaystyle u_{j+\nicefrac{{1}}{{2}}}^{-\text{(mod)}}=u_{j+\nicefrac{{1}}{{2}}}^{-}$ and $\displaystyle u_{j-\nicefrac{{1}}{{2}}}^{+\text{(mod)}}=u_{j-\nicefrac{{1}}{{2}}}^{+}$},\vspace{0.2cm}\\ \Lambda\Pi^{1}_{h}(u_{h})|_{I_{j}},&\text{otherwise}.\end{array}\right.$$

$$u_h^{\text{(mod)}}{|}_{I_j}:=\Lambda {\Pi }_h^k(u_h){|}_{I_j}=\sum _{l=0}^k{u}_j^{(l)\text{(mod)}}{\phi }_l,$$

$$\Lambda {\Pi }_h^k({\mathbf{U}}_h):={\mathbf{U}}_h^{\text{(mod)}}={\left[u_j^{(l)\text{(mod)}}\right]}_{\forall j,l=0,1,\dots ,k}.$$

$${\mathbf{U}}_h^{n,(i)\text{(mod)}}=\Lambda {\Pi }_h^k\left(\sum _{\nu =0}^{i-1}{\alpha }_{i\nu }{\mathbf{U}}_h^{n,(\nu )\text{(mod)}}+{\beta }_{i\nu }\Delta t{\mathbf{L}}_h\left({\mathbf{U}}_h^{n,(\nu )\text{(mod)}}\right)\right).$$

$$\left[t^n,t^{n+1}\right)=\bigcup _{p=0}^{M-1}\left[t^{n,p},t^{n,p+1}\right).\text{vspace}-0.1cm$$

Full text

1 Introduction

Numerical methods for hyperbolic conservation laws are a subject of great interest and importance as these laws are extensively used for modeling a wide range of physical phenomena such as gas dynamics, shallow water flow, advection of contaminants, traffic flows, etc. It is well known that these problems are often highly nonlinear and may have discontinuous solutions with sharp and moving fronts/shocks. To obtain accurate and stable numerical solutions to hyperbolic conservation laws, it is popular to use conservative high resolution methods in space together with explicit time stepping. Examples of such spatial discretization include the MUSCL (monotonic upwind scheme for conservation laws) [50], the ENO (essentially nonoscillatory) and WENO (weighted ENO) schemes [22, 23, 32, 27], and the RK-DG (Runge-Kutta discontinuous Galerkin) methods [6, 7, 8, 9]. Note that to guarantee numerical stability, the time step size needs to satisfy the CFL condition, which is determined by the spatial mesh size and wave speed. The use of local spatial refinements is efficient in resolving the sharp, moving fronts. However, as the CFL condition needs to hold everywhere, the step size for time integration would be controlled by the smallest cell size, or by the highest wave speed, which certainly increases the computational cost as a small time step size has to be used globally. Thus, to improve computational efficiency, the global CFL condition could be replaced by a local one so that the different time step sizes can be used in different regions: smaller time step sizes where the mesh is fine or the wave speed is high, and larger time step sizes where the mesh is coarse or the wave speed is low.

Explicit local time-stepping (LTS) algorithms have a long tradition. To the best of our knowledge, the first LTS algorithm for hyperbolic conservation laws was introduced in [36] for one-dimensional scalar case based on the forward Euler method in time. It is of predictor-corrector type and is first order accurate in both space and time. Extension to high resolution schemes with slope limiters for advection equations was presented in [12], and to second order in time for hyperbolic conservation laws in [13]. The numerical results on two-dimensional test problems confirm that these LTS schemes are very competitive to the global time-stepping with respect to the accuracy in time. The application of LTS schemes to the shallow water equations was investigated in [39] with a Godunov-type finite volume discretization in space and later in [49] using the RK-DG finite element methods. Note that the LTS scheme in [39] is only first order accurate in time, while the one in [49] is second order accurate in time on regions away from the LTS interface but its accuracy degrades to first order at the interface. The LTS scheme in [49] is based on the second-order strong stability preserving Runge-Kutta (SSP-RK) method, which is also known as a total variation diminishing (TVD) method introduced in [41, 45]. Higher order RK-based explicit LTS methods were introduced for conservation laws in [29, 1] and for wave propagation in [18]. In [14], a space-time fully adaptive multi-resolution method based on natural continuous extensions for RK methods was proposed, whose accuracy is of second order in both space and time. Other works related to LTS include the adaptive mesh refinement (AMR) method [2, 3], the multirate time-stepping method [11, 40] and the Implicit-Explicit (IMEX) based LTS methods [24, 19]. Among them, the AMR method involves the refinement in both space and time, i.e., small time step sizes are taken on the refined mesh and large time step sizes on the coarse mesh. It is different from our approach in the way that refined grids are placed over regions of the coarse grid and information is exchanged between the grids by means of injection and interpolation. The multirate time-stepping method allows different time step sizes in different regions but it requires buffer regions to accommodate the time scale transition between regions. An overview of LTS techniques over the last two decades can be found in [16].

In [25], inspired by the first order predictor-corrector scheme in [36], we have designed conservative second and third order explicit LTS algorithms, incorporating with SSP-RK, for the rotating shallow water equations. The model is discretized in space by a C-grid staggering finite volume method, namely the TRiSK scheme [48, 38], on orthogonal primal and dual meshes. Numerical results with parallel implementation show excellent performance of the LTS algorithms in terms of stability, accuracy, efficiency and scalability. In this work, we extend the approach to construct, in a systematic way, a framework of high order LTS algorithms for hyperbolic conservative laws. In order to derive high order LTS algorithms, the key idea is to find high order approximations on the interface at intermediate time levels to handle the coupling between coarse and fine time steppings. Our proposed schemes are also of predictor-corrector type: we derive the predictors based on Taylor series expansions of the solution at the current time level and the SSP-RK stepping algorithms at each intermediate time level. Our approach thus is different from the one proposed in [29, 1] where the predictors are based on RK time-stepping and interpolating polynomials. We present up to fourth order predictors within this framework, and show that the proposed LTS schemes preserve the accuracy in time over the entire domain. Concerning the corrector, it is designed to balance the fluxes from the regions with small time step sizes to the ones with large time step sizes. As high order SSP-RK methods consist of multiple stages, the fluxes at the same stage are accumulated over all the intermediate time levels to update the interface solution associated with that stage. As a consequence, the total mass is well conserved, though the corrector is no longer convex combinations of forward Euler steps as in the global SSP-RK methods. Nevertheless, we rigorously prove that the proposed LTS schemes for scalar conservation laws are total variation bounded (TVB). Such nonlinear stability is a crucial feature of any effective numerical method for hyperbolic conservation laws because it guarantees that the schemes can capture moving shocks without introducing nonphysical oscillations. Various numerical experiments are carried out to validate the accuracy, conservation and stability of our LTS schemes. Since time advancement of the simulations in the fine regions and in the coarse ones can be implemented in parallel (this will be discussed further in Section 3), the proposed LTS schemes preserve the natural parallelism of explicit stepping schemes.

Consider the initial value problem for hyperbolic conservation laws:

[TABLE]

where $\displaystyle\boldsymbol{u}(\boldsymbol{x}):=\left(u_{1}(\boldsymbol{x}),\ldots,u_{m}(\boldsymbol{x})\right)$ is an $\displaystyle m-$dimensional vector of unknowns and each flux function $\displaystyle\boldsymbol{f}_{i}:\mathbb{R}^{m}\longrightarrow\mathbb{R}^{m}$ defined by

[TABLE]

is vector-valued and is of $\displaystyle m$ components. Since we focus on the time discretization techniques in this paper, we shall only consider the one-dimensional case, $\displaystyle d=1$. In particular, our model problem is the following scalar hyperbolic conservation law, resulting from the system (1.1) with $\displaystyle d=m=1$:

[TABLE]

We shall construct and analyze high order Runge-Kutta discontinuous Galerkin algorithms with local time-stepping for (1.2). The proposed LTS algorithms can be straightforwardly extended to the case of one-dimensional systems of conservation laws ($\displaystyle m>1$), which will be presented in the numerical results, as well as to the higher dimensional problems ($\displaystyle d>1$).

The rest of this paper is structured as follows. In Section 2 we briefly introduce the RK-DG methods for scalar conservation laws (1.2). High order LTS algorithms are carefully derived in Section 3, and their conservation and stability properties are then proved in Section 4. Numerical results for various test cases are given in Section 5 to demonstrate the performance of the proposed LTS schemes. Additionally, coefficients of the SSP-RK methods used in the paper are given in Appendix A, and detailed derivation of the predictors for the proposed LTS schemes is presented in Appendix B.

2 Runge-Kutta discontinuous Galerkin methods

We first introduce the RK-DG methods and refer to [7] for a complete presentation of the methods. Within the framework of RK-DG, we first discretize equation (1.2) in space by the discontinuous Galerkin method, then integrate it in time by SSP-RK schemes, and finally apply a slope limiter to achieve stable and high order accurate numerical solutions.

2.1 Spatial discretization by the discontinuous Galerkin

Assume a partition of the real line $\displaystyle\mathbb{R}$ to have the $\displaystyle j$-th intervals as $\displaystyle I_{j}=\left(x_{j-\nicefrac{{1}}{{2}}},x_{j+\nicefrac{{1}}{{2}}}\right)$ and define $\displaystyle\Delta_{j}=x_{j+\nicefrac{{1}}{{2}}}-x_{j-\nicefrac{{1}}{{2}}}$ and $\displaystyle h=\max_{j}\Delta_{j}.$ Let $\displaystyle V_{h}$ be the finite dimensional space consisting of discontinuous, piecewise polynomial functions:

[TABLE]

where $\displaystyle\mathcal{P}^{k}(I_{j})$ is the space of polynomials of degree at most $\displaystyle k$ on $\displaystyle I_{j}$. Consider a weak formulation of (1.2) obtained from testing it by any function $\displaystyle v_{h}\in V_{h}$ over $\displaystyle I_{j}$:

For a.e. $\displaystyle t\in(0,T)$, find $\displaystyle u_{h}(t)\in V_{h}$ such that: $\displaystyle\forall\,j$ and $\displaystyle\forall\,v_{h}\in V_{h}$

[TABLE]

Note that we have replaced the nonlinear flux $\displaystyle f(u(x_{j+\nicefrac{{1}}{{2}}},t))$ in (2.1) by a Lipschitz, consistent, monotone numerical flux $\displaystyle h(u)_{j+\nicefrac{{1}}{{2}}}(t)$ which depends on the two values of $\displaystyle u$ at $\displaystyle x_{j+\nicefrac{{1}}{{2}}}$:

[TABLE]

The numerical flux $\displaystyle h(\cdot,\cdot)$ is required to satisfy the following properties: i) locally Lipschitz continuous; ii) consistent with the flux $\displaystyle f$, that is, $\displaystyle h(u,u)=f(u)$; and iii) nondecreasing in the first argument and nonincreasing in the second argument. Examples of such a flux include the Godunov flux, Engquist-Osher flux, Lax-Friedrichs flux and Roe flux.

A local orthogonal basis of $\displaystyle V_{h}$ consists of functions $\displaystyle\varphi_{j}^{(l)}$ defined as, for any $\displaystyle j$,

[TABLE]

in which $\displaystyle P_{l}$ is the Legendre polynomial of degree $\displaystyle l$ and $\displaystyle x_{j}$ is the middle point of $\displaystyle I_{j}$. Consequently, the approximate solution $\displaystyle u_{h}$ is expressed uniquely as

[TABLE]

where the degrees of freedom $\displaystyle u_{j}^{(l)}(t)$ are determined by

[TABLE]

Note that $\displaystyle u_{j}^{(0)}$ is the cell average of $\displaystyle u$ in $\displaystyle I_{j}$. By taking $\displaystyle v_{h}=\varphi_{j}^{(l)}$ in (2.1), we obtain the following ODE for $\displaystyle u_{j}^{(l)}$ for any $\displaystyle j$:

[TABLE]

with the initial condition

[TABLE]

Note that in (2.3) we have used the following properties of Legendre polynomials:

[TABLE]

The numerical flux $\displaystyle h$ is computed by $\displaystyle h(u_{h})_{j+\nicefrac{{1}}{{2}}}(t)=h\left(u_{j+\nicefrac{{1}}{{2}}}^{-}(t),u_{j+\nicefrac{{1}}{{2}}}^{+}(t)\right),$ where $\displaystyle u_{j+\nicefrac{{1}}{{2}}}^{\pm}(t)=u_{h}(x_{j+\nicefrac{{1}}{{2}}}^{\pm},t)$ are defined by

[TABLE]

Approximating the integral in (2.3) by Gauss-Lobatto quadrature rules that involve the two endpoints of the interval yields (using the definition of $\displaystyle u_{h}$ in (2.2)):

[TABLE]

The system of ODEs (2.3) can be recast in an autonomous form as follows:

[TABLE]

where $\displaystyle\boldsymbol{U}_{h}=(\boldsymbol{u}_{j})_{\forall\,j}$ with $\displaystyle\boldsymbol{u}_{j}=(u_{j}^{(l)})_{l=0,\ldots,k}$, the right hand side

[TABLE]

with

[TABLE]

and the initial data \displaystyle\boldsymbol{U}_{h0}=\Big{[}\nicefrac{{(2l+1)}}{{\Delta_{j}}}\int_{I_{j}}u_{0}(x)\varphi_{j}^{(l)}(x)\,dx\Big{]}_{\forall\,j,\,l=0,1,\ldots,k}. Next, we solve (2.4) explicitly in time by the SSP-RK methods [44, 17].

2.2 Strong stability preserving Runge-Kutta time discretization

The SSP-RK methods have been proved to be effective for solving hyperbolic conservation laws with discontinuous solutions. Given a uniform partition of $\displaystyle(0,T)$, $\displaystyle 0=t_{0}<t_{1}<\ldots<t_{N-1}<t_{N}=T$, with the time step size $\displaystyle\textstyle\Delta t=\nicefrac{{T}}{{N}}$. The $\displaystyle s$-stage, $\displaystyle r$th-order SSP-RK methods, referred to as SSP-RK$\displaystyle(s,r)$, for solving the autonomous system (2.4) read as follows: for $\displaystyle n=0,\ldots,N-1$, compute

[TABLE]

where $\displaystyle\boldsymbol{U}_{h}^{n,(0)}=\boldsymbol{U}_{h}^{n}$ and set $\displaystyle\boldsymbol{U}_{h}^{n+1}=\boldsymbol{U}_{h}^{n,(s)}$. It is required that all the weights $\displaystyle\alpha_{i\nu},\beta_{i\nu}\geq 0$. To measure stability of RK-DG methods, we denote by $\displaystyle\boldsymbol{u}^{(l),n}=(u_{j}^{(l),n})_{\forall j}$ and define the total variation of numerical solutions by

[TABLE]

A numerical method is total variation diminishing (TVD) if

[TABLE]

and is total variation bounded (TVB) if

[TABLE]

for some constant $\displaystyle B$ independent of the time step size. The stability of the SSP-RK schemes is given by the following lemma.

Lemma 1** ([17]).**

If the forward Euler method $\displaystyle\;\boldsymbol{U}_{h}^{n+1}=\boldsymbol{U}_{h}^{n}+\Delta t\boldsymbol{L}_{h}(\boldsymbol{U}_{h}^{n})$ is TVD under the CFL condition $\displaystyle\Delta t\leq\Delta t_{FE}$, then the SSP-RK$\displaystyle(s,r)$ scheme (2.6) is TVD under the modified CFL condition: $\displaystyle\Delta t\leq\mathcal{C}\Delta t_{FE},$ where $\displaystyle\mathcal{C}:=\min_{i,\nu}\frac{\alpha_{i\nu}}{\beta_{i\nu}}$ is the SSP coefficient.

We present some commonly used SSP-RK schemes such as SSP-RK(2,2), SSP-RK(3,3) and SSP-RK(5,4) in detail in Appendix A.

2.3 TVB corrected slope limiter

In order to handle moving shocks while preserving high order accuracy in smooth regions, we follow [41] and define the TVB corrected minmod function $\displaystyle\widetilde{m}$:

[TABLE]

where $\displaystyle C_{M}>0$ is a constant and $\displaystyle m$ is the usual minmod function [20]:

[TABLE]

The corrected limiter leads to high order accuracy in any region where the solution is smooth, even at local extrema. The resulting scheme is no longer TVD, instead it is TVB. Next, we define the $\displaystyle(k+1)$th-order limiter $\displaystyle\Lambda\Pi^{k}_{h}$ as in [5]. When $\displaystyle k=1$, we have

[TABLE]

For $\displaystyle k>1$, we first compute

[TABLE]

then define

[TABLE]

We finally make the following notation

[TABLE]

and

[TABLE]

The complete RK-DG method with the TVB minmod limiter is given in Algorithm 1, in which $\displaystyle r=(k+1)$ to match the accuracy in space and in time, and $\displaystyle s\geq r$ is the number of stages in SSP-RK.

3 Local time stepping algorithms

In this section, we present high order LTS algorithms incorporated with the RK-DG methods for conservation laws. Given the solution $\displaystyle\boldsymbol{U}_{h}^{n\text{(mod)}}$ at $\displaystyle t^{n}$, possibly with moving shocks, we approximate the solution at $\displaystyle t^{n+1}$. To this end, we divide the domain into coarse and fine regions, and assume shocks only appear in the fine regions. This could be made possible by varying the LTS interfaces with time. Consequently, we can use spatially variable time steps: large step sizes in the coarse regions and small step sizes in the fine regions.

For simplicity of presentation, we decompose the domain into a coarse region $\displaystyle\Omega_{c}^{n}$ and a fine region $\displaystyle\Omega_{f}^{n}$. Extension to more complicated configurations with multiple subdomains is straightforward. Denoted by $\displaystyle x_{j_{0}^{n}+\nicefrac{{1}}{{2}}}$ the interface point at $\displaystyle t^{n}$, $\displaystyle\Omega_{c}^{n}=\left\{I_{j}:\,j\leq j_{0}^{n}\right\}$ the coarse region, and $\displaystyle\Omega_{f}^{n}=\left\{I_{j}:\,j\geq j_{0}^{n}+1\right\}$ the fine region. As depicted in Figure 1, we enforce a larger time step $\displaystyle\Delta t_{\text{coarse}}=\Delta t$ in $\displaystyle\Omega_{c}^{n}$ and a smaller time step $\displaystyle\Delta t_{\text{fine}}=\nicefrac{{\Delta t}}{{M}}$ in $\displaystyle\Omega_{f}^{n}$. We remark that the coarse time increment must be a union of fine time increments:

[TABLE]

To proceed in time in the fine region, one needs to find $\displaystyle(k+1)$th-order in time approximation of the flux at the interface at intermediate time levels $\displaystyle t^{n,p}$ for $\displaystyle p=1,\ldots,M-1$. This is obtained via a predictor based on $\displaystyle k$th-order Taylor expansions and the $\displaystyle(k+1)$th-order SSP-RK algorithm, assuming that the solution is smooth enough near the LTS interface. After advancing in the fine region to $\displaystyle t^{n+1}$, we will correct the flux at the interface in order to conserve mass exactly. The derivation of the predictors up to fourth order accuracy are presented in Appendix B. The proposed LTS algorithm of order $\displaystyle(k+1)$ consists of the following three steps:

Step 1: Predicting the interface values. We first compute the solution of the first $\displaystyle(s-1)$ stages of the SSP-RK$\displaystyle(s,k+1)$ scheme on the interface cell $\displaystyle I_{j_{0}^{n}}$ with a coarse time step:

[TABLE]

It is important to note that we compute $\displaystyle\boldsymbol{u}_{j_{0}^{n}}^{n,(i)}$ locally on $\displaystyle I_{j_{0}}$ by enforcing $\displaystyle u_{j_{0}^{n}-\nicefrac{{1}}{{2}}}^{n,(i),-}=u_{j_{0}^{n}-\nicefrac{{1}}{{2}}}^{n,(i),+}$ and $\displaystyle u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,(i),+}=u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,(i),-}$ in (2.5). This is obtained under the assumption that the solution near the LTS interface is continuous (for $\displaystyle k\leq 1$) or sufficiently smooth (for $\displaystyle k>1$). Thus, limiter is not necessary in this case and we have $\displaystyle\boldsymbol{u}_{j_{0}^{n}}^{n,(i)\text{(mod)}}=\boldsymbol{u}_{j_{0}^{n}}^{n,(i)}.$ We then use these values to predict the solution on the interface $\displaystyle x_{j_{0}^{n}+\nicefrac{{1}}{{2}}}$ at intermediate time levels $\displaystyle t^{n,p}$:

[TABLE]

where $\displaystyle u_{j_{0}^{n}}^{(l),n,p,(i)}$ are computed by the formulas in Appendix B. In particular:

For second order SSP-RK(2,2):

[TABLE]

for $\displaystyle l=0,1,$ where $\displaystyle\textstyle\theta_{p}=\frac{p}{M}$ and $\displaystyle\textstyle\eta_{p}=\frac{p+1}{M}$ for $\displaystyle p=0,1,\ldots,M-1$.

For third order SSP-RK(3,3):

[TABLE]

for $\displaystyle l=0,1,2,$ with $\displaystyle\theta_{p}$ and $\displaystyle\eta_{p}$ as above, and $\displaystyle\textstyle\widehat{\theta}_{p}=\frac{p^{2}}{M^{2}},\,\widehat{\eta}_{p}=\frac{p(p+2)}{M^{2}},\;\gamma_{p}=\frac{2p+1}{2M}$ and $\displaystyle\textstyle\widehat{\gamma}_{p}=\frac{2p^{2}+2p+1}{2M^{2}}$ for $\displaystyle p=0,1,\ldots,M-1$.

For fourth order SSP-RK(5,4): we approximate $\displaystyle u_{j_{0}^{n}}^{(l),n,p,(i)}$, for $\displaystyle p=0,1,\ldots,M-1$ as linear combinations of $\displaystyle u_{j_{0}^{n}}^{(l),n,(i)}$ for $\displaystyle l,i=0,\ldots,4,$ as presented in Appendix B.3.

Step 2: Advancing in the coarse and fine regions in parallel.

Step 2a). Advancing the coarse region excluding the interface cell: with the solution at the current time level, we advance solution to the next time level by running the SSP-RK with a coarse time step.

For all the cells $\displaystyle I_{j}$ with $\displaystyle j<j_{0}^{n}$, we perform:

For $\displaystyle i=1,\ldots,s,$

[TABLE] 2. 2.

Set $\displaystyle\boldsymbol{u}_{j}^{n+1}=\boldsymbol{u}_{j}^{n,(s)}$ for all $\displaystyle j<j_{0}^{n}$.

Step 2b). Advancing in the fine region: with the predicted values on the interface, we evaluate the interface flux $\displaystyle h(u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,-},u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,+})$ at the intermediate time levels, and consequently obtain the solution $\displaystyle u_{j}^{n,p}$ for all the cells $\displaystyle I_{j}$ with $\displaystyle j>j_{0}^{n}$ in the fine region. The TVB limiter is performed to obtain $\displaystyle u_{j}^{n,p\text{(mod)}}$ for $\displaystyle j>j_{0}^{n}$ and the predicted values are updated on the interface after limiting.

For all the cells $\displaystyle I_{j}$ with $\displaystyle j\geq j_{0}^{n}+1$, we perform:

For $\displaystyle p=0,\ldots,M-1$,

Set $\displaystyle u_{j}^{(l),n,p,(0)(mod)}=u_{j}^{(l),n,p(mod)}$, for $\displaystyle l=0,\ldots,k$. 2. 2.

For $\displaystyle i=1,\ldots,s$, we compute the solution at stage $\displaystyle i$:

[TABLE]

for $\displaystyle l=0,\ldots,k$. If $\displaystyle p<M-1$, limit the solution in the fine region

[TABLE]

and update the predicted interface value $\displaystyle u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,(i),-\text{(mod)}}$ (cf. (3.1)) after limiting:

[TABLE] 3. 3.

For all $\displaystyle j>j_{0}^{n}$, set:

[TABLE]

Step 3: Correcting the interface solution and limiting the global solution at $\displaystyle t^{n+1}$ locally. With the predicted interface value $\displaystyle u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,(\nu),-}$, we calculate the flux at the interface $\displaystyle x=x_{j_{0}^{n}+\nicefrac{{1}}{{2}}}$. Together with the flux at $\displaystyle x=x_{j_{0}^{n}-\nicefrac{{1}}{{2}}}$, which is frozen over $\displaystyle[t^{n},t^{n+1})$, we correct the solution of the interface cell $\displaystyle I_{j_{0}^{n}}$. Finally, a TVB limiter is applied, which can be implemented in parallel as [10], to limit the solution on $\displaystyle I_{j},\,\forall j$ in which only information on elements sharing edges with $\displaystyle I_{j}$ is necessary.

For $\displaystyle i=1,\ldots,s$, we compute the solution at stage $\displaystyle i$ at the interface:

[TABLE]

where $\displaystyle\widehat{u}_{j_{0}^{n}}^{(l),n,(0)}=u_{j_{0}^{n}}^{(l),n,(0)\text{(mod)}}$. 2. 2.

Set $\displaystyle u_{j_{0}^{n}}^{(l),n+1}=\widehat{u}_{j_{0}^{n}}^{(l),n,(s)}$ and perform the limiter: $\displaystyle\boldsymbol{U}_{h}^{n+1\text{(mod)}}=\Lambda\Pi_{h}^{k}(\boldsymbol{U}_{h}^{n+1}).$

4 Properties of LTS schemes

First, we notice that the proposed LTS schemes preserve the accuracy in time of the corresponding global SSP-RK methods due to the construction of the predictor and the corrector (see also Remark 4). In the following, we prove that the LTS schemes conserve mass exactly and importantly, they satisfy the TVB stability.

4.1 Conservation

Mass conservation of the proposed LTS schemes is obtained via the construction of the corrector. For simplicity, we assume that the solutions are obtained after performing the limiter defined in Subsection 2.3 and write $\displaystyle u_{j}^{(l),n}$ for $\displaystyle u_{j}^{(l),n\text{(mod)}}$.

Theorem 2**.**

The LTS schemes exhibit exact conservation of mass:

[TABLE]

Proof.

We only need to show that mass is conserved in the region of the LTS interface $\displaystyle x=x_{j_{0}^{n}+\nicefrac{{1}}{{2}}}$, $\displaystyle I_{j_{0}^{n}}\cup I_{j_{0}^{n}+1}$, under the assumption that no flux is imposed at $\displaystyle x_{j_{0}^{n}-\nicefrac{{1}}{{2}}}$ and $\displaystyle x_{j_{0}^{n}+\nicefrac{{3}}{{2}}}$:

[TABLE]

Next, we prove (4.1) for the second order LTS scheme based on SSP-RK$\displaystyle(2,2)$ (cf. Equations (A.1)). The proof for the third and fourth order LTS schemes can be done in a similar manner; in fact, the result holds for any high order LTS schemes with the corrector defined by (3.5).

For the fine cell $\displaystyle I_{j_{0}^{n}+1}$, the second order LTS algorithm reads:

[TABLE]

Thus, by recursion, we obtain:

[TABLE]

Taking $\displaystyle v_{h}=1$ in (2.1), using (4.2) and the definition of $\displaystyle L_{h,j}^{(l)}$ in (2.5), we have

[TABLE]

as no flux is imposed at $\displaystyle x_{j_{0}^{n}+\nicefrac{{3}}{{2}}}$.

For the interface cell $\displaystyle I_{j_{0}^{n}}$, the corrector (3.5) associated with SSP-RK$\displaystyle(2,2)$ is given by

[TABLE]

from which we deduce that

[TABLE]

As for the fine cell $\displaystyle j=j_{0}^{n}+1$, we choose $\displaystyle v_{h}=1$ in (2.1) and use (4.4) to obtain

[TABLE]

noting that no flux at $\displaystyle x_{j_{0}^{n}-\nicefrac{{1}}{{2}}}$ is assumed. Thus, the proof is completed by adding (4.3) and (4.5) together. ∎

4.2 Stability

Numerical methods for conservation laws need to satisfy certain nonlinear stability requirements in order to prevent spurious oscillations when the solution is discontinuous. In [36], the first order LTS scheme based on forward Euler is proved to be TVD with the predictor obtained by freezing the value at $\displaystyle t^{n}$:

[TABLE]

For higher order LTS schemes as proposed in Section 3, multiple stage time-stepping algorithms are employed and the predictors are obtained by taking linear combinations of the interface solution at different stages with a coarse time step size. Therefore, the proof of nonlinear stability for high order LTS schemes is not an obvious generalization from the first order one. Additionally, the corrector designed to conserve mass is not a convex combination of forward Euler steps as in the case of the global SSP-RK. As a consequence, the high order LTS schemes are not TVD anymore, instead they are TVB.

We next prove the stability of the second order LTS scheme by first showing that it is TVBM (total variation bounded in the means). The generalization to higher order LTS schemes can be done in a similar manner. We introduce some notation to be used in the proof. Denoted by $\displaystyle\Delta_{+}$ and $\displaystyle\Delta_{-}$ the forward and backward finite difference operators, respectively:

[TABLE]

Following [7], we decompose the interface values $\displaystyle u_{j+\nicefrac{{1}}{{2}}}^{\pm}$ as

[TABLE]

where $\displaystyle\overline{u}_{j}:=u_{j}^{(0)}$ is the mean value of $\displaystyle u$ on the cell $\displaystyle I_{j}$. As in [6], we denote:

[TABLE]

where

[TABLE]

Note that $\displaystyle h_{1}$ and $\displaystyle-h_{2}$ are nonnegative due to the monotonicity of $\displaystyle h(\cdot,\cdot)$. Then the flux associated with the mean value $\displaystyle\overline{u}_{j}$ (cf. Equation (2.5) with $\displaystyle l=0$) can be rewritten equivalently as

[TABLE]

Using the above notation, the second order LTS scheme as presented in Section 3 for the mean value $\displaystyle\overline{u}_{j}$ reads as follows: for $\displaystyle n=0,\ldots,N-1$,

Compute the predicted mean on the interface cell at the intermediate time levels from the solutions with a coarse time step: for $\displaystyle p=0,\ldots,M-1$,

[TABLE] 2. 2.

Advance in the coarse region, for all $\displaystyle j<j_{0}^{n}$:

[TABLE]

and in the fine region, for all $\displaystyle j>j_{0}^{n}$: for $\displaystyle p=0,\ldots,M-1$,

[TABLE]

Note that the interface $\displaystyle u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,(i),-}$, for $\displaystyle i=1,2,$ and $\displaystyle p=0,\ldots,M-1$ are computed by using the second order predictor (4.7) and (3.1). 3. 3.

Correcting the interface values for which the flux at $\displaystyle x=x_{j_{0}^{n}-\nicefrac{{1}}{{2}}}$ is frozen over $\displaystyle[t^{n},t^{n+1})$:

[TABLE]

The flux term in the right-hand side of (4.10) can be rewritten as

[TABLE]

where $\displaystyle u_{j_{0}^{n}-\nicefrac{{1}}{{2}}}^{n,p,+}$ is computed by the same predictor as $\displaystyle u_{j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,-}$. In addition, we write

[TABLE]

where

[TABLE]

Moreover, using the second order predictor (3.2), we deduce that

[TABLE]

Therefore, we can rewrite the correction (4.10) as follows:

[TABLE]

Similarly,

[TABLE]

in which $\displaystyle\mathfrak{h}_{1,j_{0}^{n}+\nicefrac{{1}}{{2}}}^{n,p,(1)}$ and $\displaystyle\mathfrak{h}_{2,j_{0}^{n}-\nicefrac{{1}}{{2}}}^{n,p,(1)}$ are defined in a similar way as in (4.11) but with the solutions of the first stage $\displaystyle u_{j_{0}^{n}\pm\nicefrac{{1}}{{2}}}^{n,p,(1),\pm}$ and $\displaystyle u_{j_{0}^{n}\pm\nicefrac{{1}}{{2}}}^{n,(1),\pm}$. The TVBM property of the second order LTS scheme is guaranteed by the following theorem.

Theorem 3** (TVBM).**

Assume that there exists some $\displaystyle\theta>0$ such that

[TABLE]

for $\displaystyle n=0,\ldots,N-1$, $\displaystyle p=0,\ldots,M-1$ and $\displaystyle i=0,1.$ In addition, if a local CFL condition is satisfied:

[TABLE]

where $\displaystyle h_{1}$ and $\displaystyle-h_{2}$ are the Lipschitz coefficients of $\displaystyle h(\cdot,\cdot)$ with respect to the first and second arguments respectively, and $\displaystyle\lambda_{j}^{n,p}$ is defined by

[TABLE]

Then the second order LTS scheme is TVBM.

Proof.

Following the techniques in [36], we first introduce some important facts that will be used later in the proof. From the monotonicity of $\displaystyle h(\cdot,\cdot)$ and (4.14), we deduce that, in (4.6):

[TABLE]

We may omit the superscripts for the ease of presentation. Given any two nonnegative numbers $\displaystyle\alpha,\beta$ and suppose $\displaystyle\lambda_{j}^{n,p}=\max(\alpha,\beta)$ that satisfies (4.15), we have

[TABLE]

and consequently,

[TABLE]

in which the functions must be evaluated at the same time level. Then, together with (4.16), we deduce that

[TABLE]

In the following, we compute the variation $\displaystyle|u_{j+1}-u_{j}|$ for all $\displaystyle j$. Particularly, we consider four cases:

i) If $\displaystyle j<(j_{0}^{n}-1)$: from (4.8), we find that

[TABLE]

Applying (4.17) with $\displaystyle\alpha=\frac{\Delta t}{\Delta x_{j+1}}$ and $\displaystyle\beta=\frac{\Delta t}{\Delta x_{j}}$, we deduce from (4.18) that

[TABLE]

From this we obtain

[TABLE]

or equivalently

[TABLE]

ii) If $\displaystyle j>j_{0}^{n}$: By the same argument applied to (4.9) with a fine time step, we find that

[TABLE]

Repeating this argument inductively, we obtain a similar bound as (4.19):

[TABLE]

iii) If $\displaystyle j=j_{0}^{n}$ (the interface cell)

We aim to show that (4.20) again holds for $\displaystyle j=j_{0}^{n}$, which is the main part of the proof. Using the formulation for the corrector (4.12)-(4.13), as well as the time-stepping scheme in the fine region, we obtain

[TABLE]

Since $\displaystyle u_{j_{0}^{n}-\nicefrac{{1}}{{2}}}^{n,(1),\pm}-u_{j_{0}^{n}-\nicefrac{{1}}{{2}}}^{n,\pm}=O(\Delta t)$, and by the CFL condition, we can bound

[TABLE]

We have

[TABLE]

Regarding the first two terms, let us write

[TABLE]

By definition of the second order predictor (3.2), the last term in (4.23) is given by

[TABLE]

On the other hand, the first and second terms in (4.23) can be computed by using the time-stepping in the fine region (4.9):

[TABLE]

Summing (4.25) over $\displaystyle p=0,\ldots,M-1$ yields

[TABLE]

Substituting (4.24) and (4.26) into (4.22) and then (4.21), and using (4.17) with $\displaystyle\textstyle\alpha=\frac{p}{M}\left(\frac{\Delta t}{\Delta x_{j_{0}^{n}+1}}\right)$ or $\displaystyle\textstyle\alpha=\frac{p+1}{M}\left(\frac{\Delta t}{\Delta x_{j_{0}^{n}+1}}\right)$ and $\displaystyle\textstyle\beta=\frac{\Delta t}{\Delta x_{j_{0}^{n}}}$, we obtain:

[TABLE]

Furthermore, by the definition of the second order predictor, we have $\displaystyle\textstyle\overline{u}_{j_{0}^{n}}^{n,p,(1)}\leavevmode\nobreak\ =\leavevmode\nobreak\ \overline{u}_{j_{0}^{n}}^{n,p}+\frac{1}{M}\left(\overline{u}_{j_{0}^{n}}^{n,(1)}-\overline{u}_{j_{0}^{n}}^{n}\right)=\overline{u}_{j_{0}^{n}}^{n,p}+O(\Delta t).$ This together with using (4.17) for $\displaystyle\textstyle\alpha\leavevmode\nobreak\ =\leavevmode\nobreak\ \frac{\Delta t}{M\Delta x_{j_{0}^{n}+1}}$ and $\displaystyle\beta=0$, we have that

[TABLE]

Plugging this into (4.27) yields:

[TABLE]

On the other hand, by the SSP-RK$\displaystyle(2,2)$ time-stepping in the fine cell $\displaystyle(j_{0}^{n}+1)$ and using (4.28), we deduce that

[TABLE]

in which we have used definition of the second order predictor to obtain

[TABLE]

Repeat the above argument inductively, we arrive at

[TABLE]

Consequently,

[TABLE]

where the last inequality is obtained by reversing the order of summation. Plug this into (4.29), we find that

[TABLE]

iv) If $\displaystyle j=j_{0}^{n}-1$: It remains to investigate the case of the interface cell and its neighbor in $\displaystyle\Omega_{c}^{n}$. Similarly to (4.21), we have

[TABLE]

By performing similar manipulations as for the case $\displaystyle j=j_{0}^{n}$, one arrives at

[TABLE]

Finally, we combine (4.19), (4.20), (4.30) and (4.31) and obtain:

[TABLE]

Hence, the second order LTS scheme is TVBM. ∎

The condition (4.14) is fulfilled if the solution is limited by the minmod function $\displaystyle m$ defined in (2.7) (see [7]). The scheme remains TVB when the modified minmod function $\displaystyle\widetilde{m}$ is used, which is achieved by Theorem 2.2 in [41] (see also [7, Lemma 2.3]). Finally, the TVB property of the means $\displaystyle\overline{u}_{j}$ can be passed to whole solution $\displaystyle u_{h}$ in the same manner as the RK-DG method [7, Propositition 2.11]. We remark that it is assumed that the solution near the time-dependent LTS interface is sufficiently smooth so that the condition (4.14) is satisfied in the region of the LTS interface without limiting. In practice, local time-stepping should be coupled with adaptive spatial meshing to achieve computational efficiency and accuracy when dealing with hyperbolic conservation laws.

5 Numerical experiments

We consider several standard test cases of one dimensional scalar conservation laws (Subsection 5.1) and system of conservation laws (Subsection 5.2). We aim to verify the accuracy, mass conservation and stability of the LTS schemes as predicted theoretically and compare with those by the global time-stepping (GTS) schemes. As a first step towards study the behavior of proposed schemes, we use a fixed LTS interface (i.e., $\displaystyle j_{0}^{n}=j_{0}$ for all $\displaystyle n$ ) in all the tests, instead of a time-varying LTS interface as discussed in Section 3, and leave the investigation on parallel performance of the proposed methods with space-time adaptive multiresolution meshes in two and three dimensions to future work. Note that for the prediction step, if the solution is discontinuous at the fixed LTS interface, then it is necessary to limit the solution on the interface with a coarse time step $\displaystyle\boldsymbol{u}_{j_{0}}^{n,(i)\text{(mod)}}$ for $\displaystyle i=1,\ldots,s-1$ before calculating the predicted interface values (3.1) at intermediate time levels.

5.1 Scalar conservation laws

We first consider two model problems that obey the scalar conservation laws: the linear advection equation and Burgers’ equation. For problems with a smooth solution, we confirm the convergence order in time of our LTS algorithms. The effectiveness of LTS algorithms is demonstrated by comparing with the GTS schemes in terms of accuracy and CFL conditions.

Example 1: Linear problem

We solve the linear advection problem with a smooth initial condition

[TABLE]

in $\displaystyle-1\leq x\leq 1$ with periodic boundary conditions. The exact solution is given by $\displaystyle u(x,t)=\sin\pi(x-t)$. The spatial domain is divided into two subdomains, $\displaystyle\Omega_{1}=[-1,0]$ and $\displaystyle\Omega_{2}=[0,1]$. The mesh size and time step size in $\displaystyle\Omega_{i}$ are respectively $\displaystyle\Delta x_{i}$ and $\displaystyle\Delta t_{i}$, which are fine when $\displaystyle i=1$ and are coarse when $\displaystyle i=2$:

[TABLE]

for $\displaystyle M=1,2,4,8$, and for $\displaystyle k=1,2,3$ corresponding to second, third and fourth order LTS methods. The $\displaystyle L^{1}$ relative errors at $\displaystyle T=2$ of the three LTS algorithms are listed in Table 1. We observe that for all schemes, the errors decrease as $\displaystyle M$ increases; and the LTS schemes (with $\displaystyle M=2,4,8$) preserve the order of convergence as in the GTS case ($\displaystyle M=1$), regardless of how large $\displaystyle M$ is.

Now to show that the LTS algorithms are stable with a local CFL condition, we still consider the linear problem (5.1) but with a discontinuous initial condition

[TABLE]

We are interested in the behavior of the approximate solution near the discontinuity $\displaystyle x=-0.5$ at $\displaystyle T=0.5$ and $\displaystyle x=0$ (the LTS interface) at $\displaystyle T=1$. Again, the fine region is $\displaystyle\Omega_{1}=[-1,0]$ and the coarse one is $\displaystyle\Omega_{2}=[0,1]$. We use second order RK-DG method and consider three schemes as follows:

Coarse GTS scheme with a coarse global time step $\displaystyle\Delta t=\Delta t_{\text{coarse}}$. 2. 2.

Fine GTS scheme with a fine global time step $\displaystyle\Delta t=\nicefrac{{\Delta t_{\text{coarse}}}}{{M}}$. 3. 3.

LTS scheme with spatial variable time step $\displaystyle\Delta t_{1}=\nicefrac{{\Delta t_{\text{coarse}}}}{{M}}$ and $\displaystyle\Delta t_{2}=\Delta t_{\text{coarse}}$.

Note that the spatial mesh is refined in $\displaystyle\Omega_{1}$ by a factor of $\displaystyle M$. Table 2 shows the $\displaystyle L^{1}$ relative errors of the three schemes at $\displaystyle T=0.5$ and $\displaystyle T=1$ respectively. It is seen that the coarse GTS becomes unstable as the spatial mesh is refined due to the violation of the CFL condition, while the LTS scheme, with a valid local CFL condition, gives stable solution with the same accuracy as the fine GTS scheme.

Example 2: Burger’s equation

Next, we test the proposed algorithms on the Burgers’ equation with a smooth initial condition:

[TABLE]

in $\displaystyle-1\leq x\leq 1$. The exact solution of the problem is given by [23]:

[TABLE]

in which $\displaystyle v(x,t)$ is the solution of the Burgers’ equation with $\displaystyle v(x,0)=\sin\pi x$. We compute $\displaystyle v$ by Newton iterations to solve the characteristic relation:

[TABLE]

The solution $\displaystyle v$ in $\displaystyle(-1,0)$ is computed from $\displaystyle v$ in $\displaystyle(0,1)$ via: $\displaystyle v(-x,t)=-v(x,t)$. The solution of (5.2) is smooth up to $\displaystyle t=\nicefrac{{2}}{{\pi}}$ then it develops a moving shock. For details, see [23].

We divide the spatial domain into two zones and use the same discretization in space and in time as in Example 1. In Figure 2, we show the exact solution and the approximate solution by the fourth order LTS algorithm with $\displaystyle\Delta x_{\text{coarse}}=1/40$ and $\displaystyle M=4$. We see that LTS scheme clearly captures the shock with local refinement in space and in time. In Figure 3, mass evolution as a function of time of different LTS schemes with $\displaystyle\Delta x_{\text{coarse}}=1/40$ and $\displaystyle M=4$ is displayed. The LTS schemes conserve the mass in the region of the LTS interface, and thus in the whole domain. The relative $\displaystyle L^{1}$ errors at $\displaystyle T=0.3$ when the solution is still smooth are shown in Table 3. Again, the LTS schemes converge at the same order as the associated GTS schemes and the errors are improved as $\displaystyle M$ increases. At $\displaystyle T=1.1,$ the errors in the smooth regions ($\displaystyle 0.1$ away from the shock) are as the same magnitude as in the smooth case as displayed in Table 4.

5.2 Euler equations of gas dynamics

We next apply the proposed LTS algorithms to solve a system of one dimensional conservation laws. For the spatial discretization, we employ the DG methods for systems of equations presented in [8] with the local projection limiting in the characteristic fields. The time-stepping is still SSP-RK and thus it is straightforward to apply the proposed LTS algorithms for such a system. We consider the Euler equations of gas dynamics for a polytropic gas:

[TABLE]

with $\displaystyle P=(\gamma-1)\left(E-\nicefrac{{1}}{{2}}\rho q^{2}\right)$. Here $\displaystyle\rho,q,P$ and $\displaystyle E$ are the density, velocity, pressure and total energy, respectively; $\displaystyle m=\rho\,q$ is the momentum and $\displaystyle\gamma$ is the ratio of specific heats. In the following computation, we use $\displaystyle\gamma=1.4$ and present numerical results of applying the second order LTS algorithm to solve Riemann problems of Euler equations and the problem of interaction of blast waves. Note that for these test cases, there is no advantage of using higher order schemes as investigated in [8].

Example 3: Shock tube problem

Consider the Riemann problem

[TABLE]

with two sets of initial conditions:

• a)

The Sod problem [46]: $\displaystyle(\rho_{L},q_{L},P_{L})=(1,0,1)$ and $\displaystyle(\rho_{R},q_{R},P_{R})=(0.125,0,0.10)$;

• b)

The Lax problem [31]: $\displaystyle(\rho_{L},q_{L},P_{L})=(0.445,0.698,3.528)$ and $\displaystyle(\rho_{R},q_{R},P_{R})=(0.5,0,0.571)$.

The Sod problem has become a standard test problem of Euler equations with a monotone decreasing density profile. For this problem, we consider two settings of the decomposition into fine and coarse regions:

two subdomains with the coarse region $\displaystyle[-4.9,-0.5)$ and fine region $\displaystyle[-0.5,5.1)$; 2. 2.

three subdomains with the coarse regions $\displaystyle[-4.9,-2.9)$ and $\displaystyle[4,5.1]$, and the fine region $\displaystyle[-2.9,4)$.

The exact solution and approximation solution by the second order LTS algorithm at $\displaystyle T=2.0$ with $\displaystyle\Delta x_{\text{coarse}}=1/5$ and $\displaystyle M=4$ are shown in Figure 4 for the case of two subdomains and in Figure 5 for the case of three subdomains. The $\displaystyle L^{1}$ relative errors are displayed in Table 5 in which we observe that for both settings, the errors decrease as $\displaystyle M$ increases, especially for the three subdomain case, and the order of convergence is first order as the solution is discontinuous. The three subdomain setting, as expected, has a better performance since the fine region includes the contact discontinuities and the corners of rarefaction waves.

For the Lax problem, the density profile has a “built-up” intermediate state, thus we divide the domain into two subdomains with the coarse region $\displaystyle[-4.9,0)$ and the fine region $\displaystyle[0,5.1)$. The second order LTS approximate solution at $\displaystyle T=1.3$ with $\displaystyle\Delta x_{\text{coarse}}=1/5$ and $\displaystyle M=4$ is shown in Figure 6 together with the exact solution. We observe that the LTS scheme captures very well the “built-up” density profile with local refinement in space and in time. The $\displaystyle L^{1}$ relative errors are presented in Table 6 that match our expectation.

Example 4: Interaction of blast waves

We finally consider the problem of two interacting blast waves:

[TABLE]

with $\displaystyle(\rho_{L},q_{L},P_{L})=(1,0,10^{3})$, $\displaystyle(\rho_{M},q_{M},P_{M})=(1,0,10^{-2})$ and $\displaystyle(\rho_{R},q_{R},P_{R})=(1,0,100)$. Reflection boundary conditions are applied at $\displaystyle x=0$ and $\displaystyle x=1$. For details, see [51, 23].

We divide $\displaystyle\Omega=[0,1]$ into 3 subdomains $\displaystyle\Omega_{1}=[0,0.2)$, $\displaystyle\Omega_{2}=[0.2,0.9)$ and $\displaystyle\Omega_{3}=[0.9,1]$. The mesh and time step sizes in $\displaystyle\Omega_{2}$ are $\displaystyle M$ times smaller than those in $\displaystyle\Omega_{1}$ and $\displaystyle\Omega_{3}$. The solutions at $\displaystyle T=0.038$ with the second order global time-stepping ($\displaystyle M=1$) and local time-stepping ($\displaystyle M=2$) are shown in Figure 7. We see that the LTS scheme with a local refinement in space and time gives a much better resolution, especially for the density profile.

6 Conclusion

In this work, high order explicit local time-stepping algorithms have been proposed and analyzed for hyperbolic conservation laws. The approaches are of predictor-corrector type, and algorithms of up to fourth order accuracy are constructed in a general setting of Runge-Kutta discontinuous Galerkin methods with the modified minmod limiter. With our LTS schemes, different time-step sizes can be used based on a local CFL condition instead of the more restrictive global CFL condition. Thus, they outperform the global time-stepping for simulations on multi-resolution meshes or of multiple scales. In addition, we rigorously prove the conservation property and nonlinear stability of these schemes. Numerical results confirm their accuracy and efficiency. Future work includes the coupling of adaptive multi-resolution meshes with our local time-stepping to carry out simulations in parallel and further investigations on the numerical performance of the proposed approaches for large scale simulations on modern supercomputer systems.

Appendices

Appendix A SSP-RK time stepping schemes

We first present the SSP-RK(2,2) and SSP-RK(3,3) schemes for solution of the system (2.4), which are optimal in the sense that the number of stages equals the order of accuracy and the coefficients $\displaystyle\beta_{i,\nu}$ are nonnegative [17]. In both cases, the schemes possess the SSP coefficient $\displaystyle\mathcal{C}=1$.

Second order SSP-RK(2,2): $\displaystyle\alpha_{10}=\beta_{10}=1$, $\displaystyle\alpha_{20}=\alpha_{21}=\beta_{21}=\nicefrac{{1}}{{2}}$ and $\displaystyle\beta_{20}=0$, which is equivalent to the Heun’s method:

[TABLE]

Third order SSP-RK(3,3): $\displaystyle\alpha_{10}=\beta_{10}=1$, $\displaystyle\alpha_{20}=\nicefrac{{3}}{{4}},\,\alpha_{21}=\beta_{21}=\nicefrac{{1}}{{4}}$, $\displaystyle\alpha_{30}=\nicefrac{{1}}{{3}},\,\alpha_{32}=\beta_{32}=\nicefrac{{2}}{{3}}$ and $\displaystyle\beta_{20}=\beta_{30}=\alpha_{31}=\beta_{31}=0$, or explicitly:

[TABLE]

For higher order schemes $\displaystyle r\geq 4$, we can not avoid negative $\displaystyle\beta_{i\nu}$ without using additional stages. We shall use the SSP-RK(5,4) proposed in [30] with the SSP coefficient $\displaystyle\mathcal{C}=1.652$, and the coefficients of the scheme are listed in Table 7.

Appendix B Derivation of the predictors

We derive the predictors up to fourth order accuracy $\displaystyle(k=1,2,3)$ in the DG-RK setting. To simplify the notation, in the following, we consider an arbitrary equation in the system (2.4) and write it in the form:

[TABLE]

where $\displaystyle\boldsymbol{w}=\left(w_{j}(t)\right)_{\forall\,j}$ ($\displaystyle w_{j}$ represents $\displaystyle u_{j}^{(l)}$ and we have dropped the superscript $\displaystyle(l)$ for simplicity) and $\displaystyle L_{j}$ is a multivariable, real-valued function.

Given the time partition with coarse and fine time steps as defined in Section 3 and assume that the solution $\displaystyle\boldsymbol{w}^{n}$ at $\displaystyle t^{n}$ is known, we shall construct the approximation of $\displaystyle\boldsymbol{w}$ at the interface $\displaystyle x_{j+\nicefrac{{1}}{{2}}}$ at the intermediate time levels $\displaystyle t^{n,p}$ for $\displaystyle p=1,\ldots,M-1$. Performing Taylor expansion of $\displaystyle w_{j}$ at $\displaystyle t^{n}$ yields:

[TABLE]

The time derivatives of $\displaystyle w_{j}$ up to order $\displaystyle k$ are approximated by the SSP-RK solution of the first $\displaystyle(s-1)$ stages with a coarse time step, $\displaystyle w_{j}^{n,(i)}$ for $\displaystyle i=1,\ldots,s-1$. These approximations are detailed in the following for second, third and fourth order schemes respectively.

B.1 Predictor for the SSP-RK(2,2) scheme

We obtain the approximation for $\displaystyle w_{j}^{n,p}=w_{j}^{n,p,(0)}$ by truncating (B.2) to the second term:

[TABLE]

where $\displaystyle\Delta t_{p}=p\nicefrac{{\Delta t}}{{M}}$. To compute the first time derivative, we first notice that $\displaystyle\partial_{t}w_{j}^{n,(0)}=L_{j}(\boldsymbol{w}^{n,(0)})$ according to (B.1). Thus, by using the solution at stage 1 of SSP-RK(2,2) with a coarse time step

[TABLE]

we deduce that

[TABLE]

Substituting this into (B.3) yields:

[TABLE]

We also need to predict $\displaystyle w_{j}^{n,p,(1)}$, the solution at stage 1 at intermediate time levels, for $\displaystyle p=0,\ldots,M-1$. By definition of the SSP-RK(2,2) scheme, we have:

[TABLE]

As $\displaystyle L_{j}(\boldsymbol{w}^{n,p,(0)})=L_{j}(\boldsymbol{w}^{n,(0)})+O(\Delta t)$, we deduce from (B.6) that

[TABLE]

or equivalently via (B.5)

[TABLE]

It is clear from (B.5) and (B.6) that the proposed predictor gives second order accurate approximations of the solutions at intermediate time levels $\displaystyle t^{n,p}$, for $\displaystyle p=0,\ldots,M-1$.

B.2 Predictor for the SSP-RK(3,3) scheme

As in the second order case, we approximate $\displaystyle w_{j}^{n,p,(0)}$ by truncating (B.2), but now to the third term:

[TABLE]

The time derivatives are computed from the solutions at stage 1 and stage 2 of SSP-RK(3,3) with a coarse time step:

[TABLE]

The first time derivative can be obtained as in (B.4), for the second time derivative, by the chain rule we deduce from (B.1) that

[TABLE]

and we will compute the right-hand side by using (B.9). In particular, by performing first order Taylor expansion, we obtain:

[TABLE]

Substituting this into (B.9) yields

[TABLE]

from which we deduce that

[TABLE]

where the last equality is obtained by (B.4). Inserting this and (B.4) into (B.7), we arrive at

[TABLE]

for $\displaystyle p=0,1,\ldots,M-1$. Similarly, we can approximate $\displaystyle w_{j}^{n,p,(1)}$ and $\displaystyle w_{j}^{n,p,(2)}$, the solutions at stage 1 and stage 2 of SSP-RK(3,3) as follows:

[TABLE]

and

[TABLE]

B.3 Predictor for the SSP-RK(5,4) scheme

Again, we approximate $\displaystyle w_{j}^{n,p,(0)}$ by truncating (B.2), now to the fourth term:

[TABLE]

for $\displaystyle p=0,1,\ldots,M-1$. As for the second and third order cases, we approximate the time derivatives

[TABLE]

by using the solution the first four stages of SSP-RK(5,4) with a coarse time step:

[TABLE]

Denote by

[TABLE]

with

[TABLE]

From (B.13) we have

[TABLE]

Next, we approximate the flux by Taylor expansion with $\displaystyle O(\Delta t^{3})$ truncated error:

[TABLE]

in which $\displaystyle\Delta t^{n,(1)}$ is defined in (B.17) and the last equality is obtained by substituting the derivatives in time (B.11)-(B.12). Similarly,

[TABLE]

and

[TABLE]

For the fourth order SSP-RK scheme, the number of stages is larger than the order of the scheme. Consequently, we can compute different approximations of the time derivatives $\displaystyle\partial_{tt}w_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}w_{j}^{n,(0)}$ using either equations (B.19)-(B.20) or (B.20)-(B.21). In particular, if we substitute the equations (B.19)-(B.20) into (B.14)-(B.15), we obtain the following system for $\displaystyle\partial_{tt}w_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}w_{j}^{n,(0)}$:

[TABLE]

where

[TABLE]

By solving (B.22)-(B.23), we can compute fourth order approximations of $\displaystyle\partial_{tt}w_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}w_{j}^{n,(0)}$, denoted by $\displaystyle\partial_{tt}\overline{w}_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}\overline{w}_{j}^{n,(0)}$, as linear combinations of the solutions at different stages with a coarse time step $\displaystyle w_{j}^{n,(0)},\,w_{j}^{n,(1)},\,w_{j}^{n,(2)}$ and $\displaystyle w_{j}^{n,(3)}$:

[TABLE]

Similarly, we can substitute the equations (B.20)-(B.21) into (B.15)-(B.16) to obtain alternative approximations of $\displaystyle\partial_{tt}w_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}w_{j}^{n,(0)}$, denoted by $\displaystyle\partial_{tt}\overline{\overline{w}}_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}\overline{\overline{w}}_{j}^{n,(0)}$, as linear combinations of the solutions $\displaystyle w_{j}^{n,(0)},\,w_{j}^{n,(1)},\,w_{j}^{n,(3)}$ and $\displaystyle w_{j}^{n,(4)}$:

[TABLE]

To take into account values at all four stages, we choose the average of (B.24) and (B.25) as approximations of $\displaystyle\partial_{tt}w_{j}^{n,(0)}$ and $\displaystyle\partial_{ttt}w_{j}^{n,(0)}$ respectively and insert them into (B.10) to obtain a fourth order approximation of $\displaystyle w_{j}^{n,p,(0)}$.

Next, we approximate $\displaystyle w_{j}^{n,p,(i)},\,i=1,2,3,4$ based on the definition of the SSP-RK(5,4) scheme:

[TABLE]

Denote by

[TABLE]

with $\displaystyle\gamma^{(i)}$ defined in (B.18). We approximate $\displaystyle L_{j}(\boldsymbol{w}^{n,p,(i)}),\,i=0,1,2,3,$ as in (B.19)-(B.21) and use (B.24) and (B.25) to approximate the time derivatives:

[TABLE]

Using this we can compute $\displaystyle w_{j}^{n,p,(i)},$ $\displaystyle i=1,2,3,4$ with fourth order accuracy in time.

Remark 4**.**

By construction, the predictors for SSP-RK$\displaystyle(2,2)$, SSP-RK$\displaystyle(3,3)$ and SSP-RK$\displaystyle(5,4)$ are respectively second, third and fourth order accurate in time.