Boundary treatment of implicit-explicit Runge-Kutta method for hyperbolic systems with source terms

TL;DR

This paper introduces a high-order finite difference boundary treatment for IMEX Runge-Kutta schemes solving hyperbolic systems with source terms, ensuring accuracy and stability on Cartesian meshes.

Contribution

It develops a novel boundary treatment combining RK schemes at the boundary with an inverse Lax-Wendroff procedure, applicable to arbitrary order IMEX schemes.

Findings

Achieves third-order accuracy in boundary treatment.

Demonstrates stability and accuracy on 1D and 2D examples.

Applicable to various PDEs with IMEX RK schemes.

Abstract

In this paper, we develop a high order finite difference boundary treatment method for the implicit-explicit (IMEX) Runge-Kutta (RK) schemes solving hyperbolic systems with possibly stiff source terms on a Cartesian mesh. The main challenge is how to obtain the solutions at ghost points resulting from the wide stencil of the interior high order scheme. We address this problem by combining the idea of using the RK schemes at the boundary and an inverse Lax-Wendroff procedure. The former preserves the accuracy of the RK schemes and the latter guarantees the stability. Our method is different from the widely used approach for the explicit RK schemes by imposing boundary conditions at intermediate stages, which could not be derived for the IMEX schemes. In addition, the intermediate boundary conditions are only available for explicit RK schemes up to third order while our method applies to arbitrary order IMEX and explicit RK schemes. Moreover, the present boundary treatment method may be adapted to IMEX RK schemes solving many other partial differential equations. For a specific third-order IMEX scheme, we demonstrate the good stability and third-order accuracy of our boundary treatment through both 1D examples and 2D reactive Euler equations.

Tables

Table 1. Table 6.1: Example 1: error table for the scalar nonlinear equation ( 6.1 ).

$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	7.62e-4		1.29e-3		4.22e-3
1/40	1.13e-4	2.75	2.17e-4	2.57	9.57e-4	2.14
1/80	1.41e-5	3.01	2.59e-5	3.07	1.07e-4	3.16
1/160	1.50e-6	3.23	2.42e-6	3.42	7.20e-6	3.89
1/320	1.57e-7	3.25	2.32e-7	3.38	5.53e-7	3.70

Table 2. Table 6.2: Example 2: error table for the linear system ( 6.3 ) with solution ( 6.4 ).

$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	9.77e-4		1.06e-3		2.37e-3
1/40	1.08e-4	3.03	1.25e-4	3.08	3.20e-4	2.89
1/80	1.46e-5	2.88	1.71e-5	2.87	6.59e-5	2.28
1/160	1.78e-6	3.03	1.95e-6	3.14	5.58e-6	3.56
1/320	2.35e-7	2.93	2.52e-7	2.95	7.24e-7	2.95

Table 3. Table 6.3: Example 2: Error table for the linear system ( 6.3 ) with solution ( 6.5 ).

			$ϵ=1$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	1.01e-3		1.18e-3		2.39e-3
1/40	1.20e-4	3.08	1.34e-4	3.14	3.17e-4	2.91
1/80	1.57e-5	2.93	1.78e-5	2.91	6.49e-5	2.29
1/160	1.90e-6	3.05	2.02e-6	3.14	5.52e-6	3.56
1/320	2.49e-7	2.93	2.62e-7	2.95	7.19e-7	2.94
			$ϵ={10}^{-10}$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	9.66e-4		1.16e-3		2.38e-3
1/40	1.16e-4	3.06	1.32e-4	3.13	3.22e-4	2.89
1/80	1.53e-5	2.93	1.77e-5	2.90	6.58e-5	2.29
1/160	1.84e-6	3.04	2.01e-6	3.14	5.72e-6	3.52
1/320	2.51e-7	2.88	2.68e-7	2.91	7.07e-7	3.02

Table 4. Table 6.4: Example 3: error table for the nonlinear system ( 6.8 ) with solution ( 6.4 ).

$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	9.25e-4		1.09e-3		2.39e-3
1/40	1.13e-4	3.03	1.29e-4	3.08	3.22e-4	2.89
1/80	1.52e-5	2.89	1.76e-5	2.87	6.61e-5	2.28
1/160	1.86e-6	3.03	2.01e-6	3.13	5.61e-6	3.56
1/320	2.44e-7	2.93	2.60e-7	2.95	7.28e-7	2.95

Table 5. Table 6.5: Example 3: error table for the nonlinear system ( 6.8 ) with solution ( 6.5 ).

			$ϵ=1$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	2.52e-4		2.99e-4		4.33e-4
1/40	4.28e-5	2.56	4.77e-5	2.65	6.74e-5	2.68
1/80	6.09e-6	2.81	6.58e-6	2.86	8.98e-6	2.91
1/160	7.86e-7	2.95	8.33e-7	2.98	1.10e-6	3.03
1/320	9.25e-8	3.09	9.66e-8	3.11	1.22e-7	3.18
			$ϵ={10}^{-10}$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	9.66e-4		1.16e-3		2.38e-3
1/40	1.16e-4	3.06	1.32e-4	3.13	3.22e-4	2.89
1/80	1.53e-5	2.93	1.77e-5	2.90	6.58e-5	2.29
1/160	1.84e-6	3.04	2.01e-6	3.14	5.72e-6	3.52
1/320	2.51e-7	2.88	2.68e-7	2.91	7.07e-7	3.02

Table 6. Table 6.6: Example 4: error table for 2D reactive euler equations with solution ( 6.9 ).

			$\eta =0.5$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	1.37e-2		1.64e-2		3.18e-2
1/40	4.12e-3	1.73	5.55e-3	1.56	1.29e-2	1.30
1/80	7.78e-4	2.40	1.27e-3	2.13	3.75e-3	1.78
1/160	7.35e-5	3.40	1.35e-4	3.23	5.25e-4	2.84
			$\eta =0.7$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	1.41e-2		1.64e-2		3.10e-2
1/40	4.10e-3	1.78	5.56e-3	1.56	1.27e-2	1.28
1/80	7.79e-4	2.40	1.27e-3	2.13	3.76e-3	1.76
1/160	7.34e-5	3.41	1.35e-4	3.23	5.32e-4	2.82
			$\eta =0.9$
$\Delta x$	$L^1$ error	order	$L^2$ error	order	$L^{\infty }$ error	order
1/20	1.44e-2		1.66e-2		3.08e-2
1/40	4.12e-3	1.81	5.58e-3	1.57	1.26e-2	1.28
1/80	7.70e-4	2.42	1.27e-3	2.14	3.69e-3	1.78
1/160	7.34e-5	3.39	1.35e-4	3.23	5.27e-4	2.81

Equations

$${\partial }_{t}U+{\partial }_{x}F(U)=Q(U),$$

$$B(U(t,0),t)=0.$$

$$U^{(i)}=U^n-\Delta t\sum _{j=1}^{i-1}{\tilde{a}}_{ij}{\partial }_{x}F(U^{(j)})+\Delta t\sum _{j=1}^i{a}_{ij}Q(U^{(j)}),i=1,2,\dots ,s,$$

$$U^{n+1}=U^n-\Delta t\sum _{i=1}^s{\tilde{\omega }}_i{\partial }_{x}F(U^{(i)})+\Delta t\sum _{i=1}^s{\omega }_{i}Q(U^{(i)}).$$

$${\tilde{c}}_i=\sum _{j=1}^{i-1}{\tilde{a}}_{ij},c_i=\sum _{j=1}^i{a}_{ij}.$$

$$\frac{\Delta x}{2}=x_0<x_1<\cdots <x_N=1-\frac{\Delta x}{2}$$

$$U_j^n=\sum _{k=0}^2\frac{x_j^k}{k!}{U}^{n,(k)},j=-1,-2.$$

$$U_j^{(1)}=\sum _{k=0}^2\frac{x_j^k}{k!}{U}^{(1),(k)},j=-1,-2,$$

$$U^{(1)}(x_b)=U^n(x_b)-{\tilde{a}}_{10}\Delta t{\partial }_{x}F(U^n(x_b))+\Delta t{a}_{10}Q(U^{(1)}(x_b)).$$

$$\frac{\partial U^{(1)}}{\partial x}\big{|}_{x=x_{b}}=\partial_{x}U^{n}(x_{b})-\tilde{a}_{10}\Delta t\partial_{xx}F(U^{n}(x_{b}))+\Delta ta_{10}Q_{U}(U^{(1)}(x_{b}))\frac{\partial U^{(1)}}{\partial x}\big{|}_{x=x_{b}}.$$

$$A(U_{0}^{n})=\partial_{U}F(U)\big{|}_{U=U_{0}^{n}}$$

$$V_{j,m}=l_m{U}_j^n,m=1,2,\dots ,M,j=0,1,2.$$

$$V_{x_b,m}^{\ast (k)},k=0,1,2.$$

$$\begin{array}{rl}& l_m{U}^{n,(0)}=V_{x_b,m}^{\ast (0)},m=p+1,p+2,\dots ,M,\\ & B(U^{n,(0)},t_n)=0.\end{array}$$

$$B_U(U^{n,(0)},t_n){\partial }_t{U}^{n,(0)}+B_t(U^{n,(0)},t_n)=0,$$

$$B_U(U^{n,(0)},t_n){\partial }_t{U}^{n,(0)}=g(U^{n,(0)},t_n)$$

$$B_U(U^{n,(0)},t_n){\partial }_t{U}^{n,(0)}+B_U(U^{n,(0)},t_n)A(U^{n,(0)})U^{n,(1)}=B_U(U^{n,(0)},t_n)Q(U^{n,(0)}).$$

$$l_m{U}^{n,(1)}=V_{x_b,m}^{\ast (1)},m=p+1,p+2,\dots ,M,$$

$$B_U(U^{n,(0)},t_n)A(U^{n,(0)})U^{n,(1)}=B_U(U^{n,(0)},t_n)Q(U^{n,(0)})-g(U^{n,(0)},t_n).$$

$$U^{n,(2)}=L^{-1}{V}_{x_b,m}^{\ast (2)},$$

$$V^{*(k)}=\sum_{r=0}^{2}d_{r}\frac{\mbox{d}^{k}p_{r}(x)}{\mbox{d}x^{k}}\big{|}_{x=0},$$

$$V^{*(k)}=\sum_{r=0}^{2}\omega_{r}\frac{\mbox{d}^{k}p_{r}(x)}{\mbox{d}x^{k}}\big{|}_{x=0},$$

$${\omega }_r=\frac{{\alpha }_r}{{\sum }_{s=0}^2{\alpha }_s},{\alpha }_r=\frac{d_r}{(\epsilon +{\beta }_r{)}^2}.$$

$${\beta }_0=\Delta x^2,{\beta }_r=\sum _{l=1}^r{\int }_{x_{-1}}^{x_0}\Delta x^{2l-1}{\left(\frac{\text{mbox}{d}^l}{\text{mbox}d{x}^l}{p}_r(x)\right)}^2\text{mbox}dx,r=1,2.$$

$${\omega }_0=\mathcal{O}(\Delta x^2),{\omega }_1=\mathcal{O}(\Delta x),{\omega }_2=1-{\omega }_0-{\omega }_1$$

$${\partial }_{t}u+{\partial }_{x}u=0,0\le x\le 1,t>0,$$

$$u(t,0)=g(t)$$

$$u^{(1)}=u^n+\Delta t\mathcal{L}(u^n),$$

$$u^{(2)}=\frac{3}{4}{u}^n+\frac{1}{4}{u}^{(1)}+\frac{1}{4}\Delta t\mathcal{L}(u^{(1)}),$$

$$u^{n+1}=\frac{1}{3}{u}^n+\frac{2}{3}{u}^{(2)}+\frac{2}{3}\Delta t\mathcal{L}(u^{(2)}),$$

$$u_b^n=g(t_n).$$

$$u_b^{n,(1)}={\partial }_{x}u(t_n,x_b)=-{\partial }_{t}u(t_n,x_b)=-g^{\prime }(t_n).$$

$$u_b^{n,(2)}=g^{\prime \prime }(t_n),$$

Full text

Boundary treatment of implicit-explicit Runge-Kutta method for hyperbolic systems with source terms

Weifeng Zhao111Department of Applied Mathematics, University of Science and Technology Beijing, Beijing 100083, China. E-mail: wfzhao@ustb.edu.cn. Research supported by the National Natural Science Foundation of China (No. 11801030, No. 11861131004). , Juntao Huang222Department of Mathematics, Michigan State University, East Lansing, MI 48824, USA. E-mail: huangj75@msu.edu.

Abstract

In this paper, we develop a high order finite difference boundary treatment method for the implicit-explicit (IMEX) Runge-Kutta (RK) schemes solving hyperbolic systems with possibly stiff source terms on a Cartesian mesh. The main challenge is how to obtain the solutions at ghost points resulting from the wide stencil of the interior high order scheme. We address this problem by combining the idea of using the RK schemes at the boundary and an inverse Lax-Wendroff procedure. The former preserves the accuracy of the RK schemes and the latter guarantees the stability. Our method is different from the widely used approach for the explicit RK schemes by imposing boundary conditions at intermediate stages [12], which could not be derived for the IMEX schemes. In addition, the intermediate boundary conditions are only available for explicit RK schemes up to third order while our method applies to arbitrary order IMEX and explicit RK schemes. Moreover, the present boundary treatment method may be adapted to IMEX RK schemes solving many other partial differential equations. For a specific third-order IMEX scheme, we demonstrate the good stability and third-order accuracy of our boundary treatment through both 1D examples and 2D reactive Euler equations.

Keywords: Boundary treatment; hyperbolic systems with source terms; implicit-explicit Runge-Kutta method; high order finite difference; inverse Lax-Wendroff procedure

1 Introduction

Implicit-explicit (IMEX) Runge-Kutta (RK) schemes [7] are a powerful tool for the time discretization of many partial differential equations (PDEs), such as convection-diffusion equations [28, 21, 29], hyperbolic systems with stiff sources [22, 23, 11, 8, 9, 10], and the Boltzmann equation [13, 16, 15, 19]. They achieve good stability by implicitly treating the stiff terms in the PDEs and high order accuracy by using several immediate stages. Though powerful and efficient, the IMEX RK schemes for the time-dependent PDEs are not without drawbacks. A long-standing problem is the lack of appropriate boundary treatment method [7], which is fundamental for practical applications of IMEX RK schemes.

For the explicit RK schemes solving hyperbolic equations, a popular boundary treatment method is to impose consistent boundary conditions for each intermediate stage [12]. However, these boundary conditions are derived only for explicit RK schemes up to third order. Though it is remedied for fourth order schemes in [1, 24], the methods therein only apply to one-dimensional (1D) scalar equations or systems with all characteristics flowing into the domain at the boundary. Moreover, the intermediate-stage boundary conditions are not available for IMEX RK schemes since the implicit part would result in nonlinear equations which cannot be analytically solved in general. On the other hand, boundary conditions for intermediate stages are derived and analyzed in [2, 3, 4, 5, 6] for implicit RK schemes solving linear scalar equations. At this point, we remark that another limitation of all the boundary treatments mentioned above [12, 1, 24, 2, 3, 4, 5, 6] is that they only deal with 1D case where the boundary is located at a grid point, restricting their applications for 2D problems on complex domains.

Perhaps due to the difficulties of extending the above boundary treatments of explicit and implicit schemes to the IMEX case, there are rare works on that of IMEX RK schemes. To our knowledge, the only literature on this topic is [29]. In that reference, a strategy of imposing intermediate boundary conditions is presented for one-dimensional (1D) linear convection-diffusion equations. There is a similar restriction that it does not apply to general nonlinear equations.

In this paper, we propose a boundary treatment method for the IMEX RK schemes solving hyperbolic systems with possibly stiff source terms. High order finite difference schemes on a Cartesian mesh are used for the space discretization. Due to the wide stencil of high order schemes, ghost points are needed near the boundary where the numerical stencil is partially outside of the computational domain. Then the main task of boundary treatment is to obtain the solutions at the ghost points. We address this problem by combining the idea of using the RK schemes at the boundary and an inverse Lax-Wendroff (ILW) procedure [25, 26]. The ILW procedure is employed in [25, 26] as a technique of boundary treatment for hyperbolic systems by using the PDEs to convert spatial derivatives of inflow boundaries to time derivatives, while those of the outflow boundaries are approximated with the solutions at interior points. Nevertheless, the boundary treatment in [25, 26] relies on the intermediate-stage boundary conditions [12], which are not available for IMEX schemes as aforementioned. To overcome this difficulty, we propose to use the RK schemes at the boundary directly, instead of imposing intermediate-stage boundary conditions. Combining this idea and the ILW procedure, we can obtain the solutions as well as their derivatives at the boundary and then the solutions of the ghosts points are computed by the Taylor expansion at the boundary. In this way, our boundary treatment not only preserves the accuracy of the IMEX RK schemes but also processes good stability. In addition, it allows the boundary to be arbitrarily located. While our method applies to general IMEX RK schemes, we demonstrate its effectiveness with a specific third-order scheme in [7] by fixing the CFL number as 0.8. The third-order accuracy of the method is verified for both 1D and 2D cases, and the numerical results of the 2D reactive Euler equations are comparable to those obtained by the third-order positivity-preserving schemes with reflective boundary conditions [17].

It is interesting to note that for 1D scalar equations we show that our method applied to third-order IMEX RK scheme can yield the widely used intermediate-stage boundary conditions in [12], indicating the correctness of our boundary treatment for IMEX RK schemes. Since the IMEX RK schemes include explicit schemes as special cases, the present boundary treatment also applies to explicit RK schemes. It does not rely on intermediate-stage boundary conditions and thus is expected to be valid for explicit RK schemes of order higher than three. This will be investigated in our next work. Moreover, our method may be adapted to IMEX RK schemes solving the other partial differential equations, *e.g. *convection-diffusion equations [28, 21, 29] and the Boltzmann equation [13, 16, 15, 19]. Many unresolved issues for the boundary treatment of RK schemes can be explored based on the idea of this work.

This paper is organized as follows. In Section 2, we introduce the IMEX RK schemes for hyperbolic systems with source terms. In Section 3, we use the 1D systems to illustrate our idea of boundary treatment. In Section 4, we show that for 1D scalar equations our boundary treatment can yield the widely used intermediate-stage boundary conditions. We extend our boundary treatment to the 2D reactive Euler equations in Section 5. Section 6 is devoted to numerical validations of our method. Finally, some conclusions and remarks are given in Section 7.

2 IMEX RK schemes for hyperbolic systems with source terms

Consider a one-dimensional nonlinear hyperbolic system with source term

[TABLE]

where $0\leq x\leq 1$ and $U,F(U),Q(U)\in\mathbb{R}^{M}$. In the domain $0\leq x\leq 1$, appropriate boundary conditions and initial data $U(0,x)$ are required to close the above system. We assume that the Jacobian matrix $F_{U}(U(t,0))$ always has $p$ positive eigenvalues and thus $p$ independent relations among incoming and outgoing modes are given at the left boundary $x=0$:

[TABLE]

Similarly, the boundary condition at $x=1$ are also imposed properly.

A general IMEX RK scheme for the hyperbolic system (2.1) is [11]

[TABLE]

Here the $s\times s$ matrices $\tilde{A}=(\tilde{a}_{ij})$ and $A=(a_{ij})$ satisfy $\tilde{a}_{ij}=0$ for $j\geq i$ and $a_{ij}=0$ for $j>i$ so that the scheme is explicit for the convection part and implicit for the source term. Along with the coefficient vectors $\tilde{w}=(\tilde{w}_{1},\tilde{w}_{2},\ldots,\tilde{w}_{s})^{T}$ and $w=(w_{1},w_{2},\ldots,w_{s})^{T}$, the scheme can be represented by a double Butcher tableau

[TABLE]

[TABLE]

where $\tilde{c}=(\tilde{c}_{1},\tilde{c}_{2},\ldots,\tilde{c}_{s})^{T}$, $c=(c_{1},c_{2},\ldots,c_{s})^{T}$ are defined as

[TABLE]

For the space discretization in (2.3), we use the finite difference WENO scheme on a Cartesian mesh [18], which needs a stencil of several points in each direction. As a consequence, ghost points are needed near the boundary $\partial\Omega$ where the numerical stencil is partially outside of $\Omega$. For explicit RK schemes, solutions at ghost points can be determined as in [25, 26] with the aid of an ILW procedure and the boundary condition of each stage constructed in [12]. However, this approach does not apply to IMEX schemes since the boundary conditions of intermediate stages could not be analytically derived for general nonlinear systems. Thus boundary treatment of IMEX RK schemes remains a challenge and we address this problem in the following.

3 Boundary treatment for 1D system

In this section, we use the one-dimensional system (2.1) to illustrate our idea of boundary treatment for the IMEX RK scheme (2.3), and then extend it to the two-dimensional case in Section 5. We employ the third-order IMEX RK scheme for concreteness, while our method applies to general IMEX schemes. Additionally, we use the third-order finite difference WENO scheme in space [20] and take $\Delta t=\mathcal{O}(\Delta x)$.

3.1 Computation of solutions at ghost points

We focus on the left boundary $x=0$ of the problem (2.1) and the method can be similarly applied to the right boundary. Following the notations in [26], we discretize the interval $(0,1)$ by a uniform mesh

[TABLE]

and set $x_{j},j=-1,-2$ as two ghost points near the left boundary $x=0$ (note that the boundary indeed can be arbitrarily located). Denote by $U_{j}^{n}$ the numerical solution of $U$ at position $x_{j}$ and time $t_{n}$. Assume that the interior solutions $U_{j},j=0,1,2,\ldots,N$ have been updated from time level $t_{n-1}$ to time level $t_{n}$. As in [26], we use a third-order Taylor approximation to construct the values at the ghost points $x_{j}$ for $j=-1,-2$, namely,

[TABLE]

Here $U^{n,(k)}$ denotes a $(3-k)$-th order approximation of the spatial derivative at the boundary point \frac{\partial^{k}U}{\partial x^{k}}\big{|}_{x=0,t=t_{n}}. With this formula, $U_{j}^{n}$ at the ghost points can be obtained once $U^{n,(k)},k=0,1,2$ are provided, which is well treated in [25, 26] (see also the next subsection for details). Then we can obtain $U_{j}^{(1)}$ at the interior points by the RK scheme (2.3a). The next task is to compute $U_{j}^{(1)},j=-1,-2$. We propose an approach to do this as follows.

Similar to the above procedure, we compute $U_{j}^{(1)}$ with the third-order Taylor expansion at the boundary point $x_{b}=0$:

[TABLE]

where $U^{(1),(k)}$ denotes a $(3-k)$-th order approximation of the spatial derivative \frac{\partial^{k}U^{(1)}}{\partial x^{k}}\big{|}_{x=0}. Then we turn to compute $U^{(1),(k)}$ for $k=0,1,2$. To this end, we apply the first stage of the RK solver (2.3a) for $U^{(1)}$ at the boundary point $x_{b}=0$:

[TABLE]

Notice that $U^{n}(x_{b})=U^{n,(0)}+\mathcal{O}(\Delta x^{3})$ and $\partial_{x}F(U^{n}(x_{b}))=F_{U}(U^{n}(x_{b}))\partial_{x}U^{n}(x_{b})=F_{U}(U^{n,(0)})U^{n,(1)}+\mathcal{O}(\Delta x^{2})$ are already known. Substituting these into (3.4) and solving the algebra equation, we can obtain an approximation of $U^{(1)}(x_{b})$ and denote it by $U^{(1),(0)}$. In addition, the error between $U^{(1),(0)}$ and $U^{(1)}(x_{b})$ defined by (3.4) is $\mathcal{O}(\Delta x^{3})$.

Furthermore, taking derivatives with respect to $x$ on both sides of (3.4) yields

[TABLE]

Here $\partial_{x}U^{n}(x_{b})=U^{n,(1)}+\mathcal{O}(\Delta x^{2})$ and $\partial_{xx}F(U^{n}(x_{b}))=F_{UU}(U^{n}(x_{b}))\partial_{x}U^{n}(x_{b})\partial_{x}U^{n}(x_{b})+F_{U}(U^{n}(x_{b}))\partial_{xx}U^{n}(x_{b})=F_{UU}(U^{n,(0)})U^{n,(1)}U^{n,(1)}+F_{U}(U^{n,(0)})U^{n,(2)}+\mathcal{O}(\Delta x)$. With these approximations, \frac{\partial U^{(1)}}{\partial x}\big{|}_{x=x_{b}} can be solved out from (3.5) and the resulting solution $U^{(1),(1)}$ is a second-order approximation of \frac{\partial U^{(1)}}{\partial x}\big{|}_{x=x_{b}}.

By taking higher order derivatives on both sides of (3.4), one can also compute \frac{\partial^{k}U^{(1)}}{\partial x^{k}}\big{|}_{x=x_{b}} for $k=2$. However, this procedure is quite complicated as it involves Jacobian of Jacobian. Here we simply approximate \frac{\partial^{k}U^{(1)}}{\partial x^{k}}\big{|}_{x=x_{b}} for $k=2$ by using the $(3-k)$-th order WENO type extrapolation with $U^{(1)}$ at interior points in [26], which will be introduced in subsection 3.3. In this way, we obtain $U^{(1),(k)}$ for $k=0,1,2$ with accuracy of order $(3-k)$. Then $U_{j}^{(1)}$ for $j=-1,-2$ can be computed by (3.3) with third-order accuracy. Having $U^{(1)}$ at the ghost points, we can then evolve from $U^{(1)}$ to $U^{(2)}$ using the interior difference scheme.

Repeating the same procedure for each $U^{(i)},2\leq i\leq s$, we can compute the solution at the ghost points in the $i$-th intermediate stage. Then we can update the solution at all interior points in the $(i+1)$-th stage. Finally, we obtain $U^{n+1}$, i.e., the solution at the end of a complete RK cycle.

3.2 Computation of $U^{n,(k)}$ at the boundary

For the sake of completeness, we provide the method in [26] for computing $U^{n,(k)},k=0,1,2$, i.e., the $(3-k)$-th order approximation of $\frac{\partial^{k}U^{n}}{\partial x^{k}}$ at the boundary $x=0$.

3.2.1 $k=0$

We first do a local characteristic decomposition to determine the inflow and outflow boundary conditions as in [26]. Denote the Jacobian matrix of the flux evaluated at $x=x_{0}$ by

[TABLE]

and assume that it has $p$ positive eigenvalues $\lambda_{1},\lambda_{2},\ldots,\lambda_{p}$ and $(M-p)$ negative eigenvalues $\lambda_{p+1},\lambda_{p+2},\ldots,\lambda_{M}$ with $l_{1},l_{2},...l_{p}$ and $l_{p+1},l_{p+2},...l_{M}$ the corresponding left eigenvectors, respectively. Define by $V_{j,m}$ the $m$-th component of local characteristic variable at grid point $x_{j},j=0,1,2$, i.e.,

[TABLE]

We extrapolate the outgoing characteristic variable $V_{j,m},m=p+1,p+2,\ldots,M$ to the boundary $x_{b}$ with the WENO type extrapolation [26] (see also subsection 3.3), and denote the extrapolated $k$-th order derivative by

[TABLE]

With $V_{x_{b},m}^{*(0)}$ and the boundary condition (2.2), $U^{n,(0)}$ at the boundary can be solved out from the following equation

[TABLE]

3.2.2 $k=1$

Having $U^{n,(0)}$, we proceed to compute $U^{n,(1)}$ with the ILW procedure proposed in [25, 26]. To do this, we take derivative with respect to $t$ for the boundary condition (2.3)

[TABLE]

which can be written as

[TABLE]

with $g(U^{n,(0)},t_{n}):=-B_{t}(U^{n,(0)},t_{n})$. On the other hand, multiplying the equation (2.1) with $B_{U}(U^{n,(0)},t_{n})$ from the left yields

[TABLE]

With the above two equations and $V_{x_{b},m}^{*(1)}$ obtained by the WENO type extrapolation, $U^{n,(1)}$ can be solved out from

[TABLE]

This is the ILW procedure in [25, 26].

3.2.3 $k=2$

Following [26], we approximate $\frac{\partial^{2}U^{n}}{\partial x^{2}}|_{x=0}$ with the WENO type extrapolation. Specifically, since we have obtained the characteristic variables at grid points near the boundary in (3.6), the WENO type extrapolation is employed based on these characteristic variables to compute the second-order derivative $V_{x_{b},m}^{*(2)}$ for each $m=1,2,\cdots,M$. Then the approximation of $\frac{\partial^{2}U^{n}}{\partial x^{2}}|_{x=0}$ is given by

[TABLE]

where $L$ is local left eigenvector matrix composed of $l_{1},l_{2},\cdots,l_{M}$. Similarly, $\frac{\partial^{k}U^{n}}{\partial x^{k}}|_{x=0},k>2$ are also computed with the WENO type extrapolation if necessary for higher order boundary treatments.

3.3 WENO type extrapolation

Finally, we show how to obtain a $(3-k)$-th order approximation of $\frac{\partial^{k}V}{\partial x^{k}}|_{x=0}$, denoted by $V^{*(k)},k=0,1,2$, with the stencil $S_{3}=\{x_{0},x_{1},x_{2}\}$, which is introduced in [25]. When $V$ is smooth on $S_{3}$, we compute $V^{*(k)}$ as

[TABLE]

where $d_{0}=\Delta x^{2},d_{1}=\Delta x,d_{2}=1-\Delta x-\Delta x^{2}$ and $p_{r}(x)$ is the $r$-th order Lagrange polynomial constructed with $S_{3}$. When there is a discontinuity in the stencil $S_{3}$, the WENO type extrapolation is applied [25]:

[TABLE]

where $\omega_{r}$ are nonlinear weights depending on the value of $V_{j}$. Similar to the classical WENO scheme, the nonlinear weights are given by [25]

[TABLE]

Here $\varepsilon=10^{-6}$ and $\beta_{r}$ are the smoothness indicators determined by

[TABLE]

The explicit expressions of $\beta_{1}$ and $\beta_{2}$ are given in [25]. When $V$ is smooth on $S_{3}$, it is shown in [25] that

[TABLE]

so that the WENO type extrapolation (3.11) is $(3-k)$-th order accurate.

In the end of this section, we summarize our third-order boundary treatment for the hyperbolic system (2.1) as follows:

Step 1: Compute $U^{n,(k)}$, the $(3-k)$-th order approximation of $\frac{\partial^{k}U^{n}}{\partial x^{k}}$ at the boundary $x=0$ for $k=0,1,2$ as in subsection 3.2. Then impose $U^{n}_{j}$ at the ghost points $j=-1,-2$ with the Taylor expansion (3.2). This step is the same as that in [26] for explicit RK schemes. With $U^{n}$ at the ghost points, we can evolve from $U^{n}$ to $U^{(1)}$ with the interior difference scheme.

Step 2: Compute $U^{(1),(k)}$, the $(3-k)$-th order approximation of $\frac{\partial^{k}U^{(1)}}{\partial x^{k}}$ at the boundary $x=0$ for $k=0,1,2$ as in subsection 3.1. Specifically, $U^{(1),(0)}$ and $U^{(1),(1)}$ are solved out from the RK solver (3.4) and (3.5), respectively, and $U^{(1),(2)}$ is approximated with the WENO type extrapolation in subsection 3.3. Then impose $U^{(1)}_{j}$ at the ghost points $j=-1,-2$ with the Taylor expansion (3.3). Having $U^{(1)}$ at the ghost points, we can update $U^{(2)}$ at all interior points.

Step 3: Repeat the same procedure as in Step 2 for each $U^{(i)},2\leq i\leq s$ to compute the solution at the ghost points in the $i$-th intermediate stage. Then we can update the solution at all interior points in the $(i+1)$-th stage. Finally, we obtain $U^{n+1}$, i.e., the solution at the end of a complete RK cycle.

Remark 3.1**.**

Combining the idea of using the RK schemes of intermediate solutions at the boundary and the ILW procedure is the key of our boundary treatment. The former preserves the accuracy of the RK method and the latter guarantees the stability. This is different from the widely used approach for the explicit RK scheme by imposing boundary conditions for intermediate solutions [12], which are unavailable for the IMEX schemes. In addition, the intermediate boundary conditions are only available for explicit RK schemes up to third order while our method applies to general IMEX RK scheme.

Remark 3.2**.**

Though the idea of using the RK schemes of intermediate solutions has been stressed in [24], the method therein is only valid for 1D scalar equation or 1D system with all characteristics flowing into the domain at the boundary. Another limitation of the method in [24] is that the boundary has to be located at a grid point, which makes it difficult to be extended to 2D case.

Remark 3.3**.**

Since the IMEX RK schemes include explicit schemes as special cases, the present boundary treatment also applies to explicit RK schemes. It does not need to impose intermediate boundary conditions and thus is expected to be valid for explicit RK schemes of order higher than three. This will be investigated in our next work.

4 Analysis of the boundary treatment

In this section, we show that our boundary treatment applied to the third-order explicit RK scheme solving 1D scalar equations can yield the widely used boundary conditions for intermediate solutions in [12].

4.1 Linear equation

We consider a linear scalar equation

[TABLE]

and solve it with the third-order explicit strong-stability-preserving RK scheme [14]

[TABLE]

where $\mathcal{L}(u)=-\partial_{x}u$. We focus on the boundary treatment at the left boundary $x_{b}=0$.

According to the boundary treatment in the previous section, we first need to compute $u^{n}$ and its first- and second-order derivatives at the boundary point $x=x_{b}$, denoted by $u^{n}_{b}$ and $u^{n,(k)}_{b},k=1,2$, respectively. For the above scalar equation, $u^{n}_{b}$ is simply given by the boundary condition:

[TABLE]

Then according to subsection 3.2.2, the ILW procedure is employed to compute the first-order derivative $u^{n,(1)}_{b}$:

[TABLE]

Similarly, the second-order derivative $u^{n,(2)}_{b}$ computed by the ILW procedure is

[TABLE]

and by the third-order WENO extrapolation is

[TABLE]

Here the solution is assumed to be smooth near the boundary for (4.6).

Next we compute $u^{(1)}_{b}$ and its derivatives $u^{(1),(k)}_{b},k=1,2$. In our method, the RK scheme for $u^{(1)}$ is applied to determine its value at the boundary, i.e.,

[TABLE]

where (4.3) and (4.4) have been used. As shown in (3.5), we take derivatives with respect to $x$ on both sides of (4.2a) and evaluate at $x_{0}$ to obtain

[TABLE]

and thereby

[TABLE]

where (4.4) has been used. Applying the scheme (4.2b) for $u^{(2)}$ at the boundary $x_{0}$ and using (4.3), (4.7) and (4.8), we arrive at

[TABLE]

If we use (4.5) for $u^{n,(2)}_{b}$ obtained by the ILW procedure, then

[TABLE]

In this case, the relations (4.3), (4.7) and (4.10) are exactly the boundary conditions imposed at each stage in [12, 25].

If the WENO extrapolation is adopted for $u^{n,(2)}_{b}$, then

[TABLE]

which only differs from that in (4.10) with a third-order term $\mathcal{O}(\Delta t^{3})$ and thus maintains the third-order accuracy of the boundary treatment.

4.2 Nonlinear equation

We further show that the above analysis also applies to the nonlinear case. To this end, we use the third-order RK scheme (4.2) for a nonlinear scalar equation

[TABLE]

with $f^{\prime}(u(0,t))>0$, and apply our method at the left boundary $x_{b}=0$.

First, it is trivial that the boundary condition for $u^{n}$ is $u^{n}_{b}=g(t_{n})$. Using the ILW procedure, we can obtain [25]

[TABLE]

According to our boundary treatment, the RK scheme (4.2a) for $u^{(1)}$ is used at the boundary to determine $u^{(1)}_{b}$:

[TABLE]

Moreover, taking derivatives with respect to $x$ on both sides of the scheme (4.2a) and evaluating at the boundary point $x=x_{b}$, we have

[TABLE]

Thus $u^{(2)}_{b}$ determined by the scheme (4.2b) is

[TABLE]

Equations (4.13) and (4.15) are the boundary conditions for $u^{(1)}$ and $u^{(2)}$, respectively. In particular, (4.13) is the same with (4.7), and (4.15) degenerates to (4.10) when $f(u)$ is linear. Moreover, for nonlinear $f(u)$ we conduct the Taylor expansion

[TABLE]

and then (4.15) reduces to

[TABLE]

The above expression of $u^{(2)}_{b}$ only differs from (4.10) with $\mathcal{O}(\Delta t^{3})$ and thus the overall accuracy is still third order.

This analysis shows the correctness of our boundary treatment when it is applied to the explicit third-order RK scheme.

5 Boundary treatment for 2D reactive Euler equations

In this section, we extend the above idea of boundary treatment to the 2D case. For convenience of representation, we concentrate on the 2D reactive Euler equations

[TABLE]

where

[TABLE]

and $E=\frac{1}{2}\rho(u^{2}+v^{2})+\frac{p}{\gamma-1}+\rho qY$. Here $\rho$ is the density, $u$ and $v$ are the velocities in $x$ and $y$ directions, $E$ is the total energy, $p$ is the pressure, $Y$ is the reactant mass fraction, $q>0$ is the heat release of reaction and $\gamma$ is the specific heat ratio. The source term is assumed to be in an Arrhenius form

[TABLE]

where $T=p/\rho$ is the temperature, $\tilde{T}>0$ is the activation constant temperature and $\tilde{K}>0$ is a constant rate coefficient.

The IMEX RK scheme for the 2D reactive Euler equations (5.1) is

[TABLE]

where the coefficients $\tilde{a}_{ij}$, $a_{ij}$, $\tilde{\omega}_{i}$ and $\omega_{i}$ are the same with those in (2.3). Assume the solutions of all the grid points inside $\Omega$ have been updated from time level $t_{n-1}$ to $t_{n}$. For a ghost point $P$, we find a point $P_{0}=(x_{0},y_{0})=\bm{x}_{0}$ on the boundary $\partial\Omega$ such that the normal $\bm{n}(\bm{x}_{0})$ at $P_{0}$ goes through $P$. Here $\bm{n}(\bm{x}_{0})$ points to the exterior of $\Omega$. We set up a local coordinate system at $P_{0}$ by

[TABLE]

where $\theta$ is the angle between the normal $\bm{n}(\bm{x}_{0})$ and the $x$-axis, and $\mathsf{T}$ is a rotational matrix. The $\hat{x}$-axis then points in the same direction as $\bm{n}(\bm{x}_{0})$ and the $\hat{y}$-axis points in the tangent direction. In this local coordinate system, the reactive Euler equations (5.1) become

[TABLE]

with

[TABLE]

The IMEX scheme in this coordinate is

[TABLE]

For a third-order boundary treatment, $\hat{U}^{n}$ at the ghost point $P$ is imposed by the Taylor expansion

[TABLE]

where $\Delta$ is the $\hat{x}$-coordinate of $P$ and $\hat{U}^{n,(k)}$ is a $(3-k)$-th order approximation of the normal derivative \frac{\partial^{k}\hat{U}}{\partial\hat{x}^{k}}\big{|}_{(x,y)=\bm{x}_{0},t=t_{n}}. With $\hat{U}^{n}$ at ghost points, $\hat{U}^{(1)}$ at interior points can be obtained with the scheme (5.5a). Since the computation of $\hat{U}^{n,(k)}$ is the same with that in [26] and its basic idea has been illustrated for the 1D case, we omit the details here.

The next step is to compute $\hat{U}^{(1)}$ at the ghost points. Similar to the 1D case, we use the Taylor expansion

[TABLE]

with $\hat{U}^{(1),(k)}$ a $(3-k)$-th order approximation of \frac{\partial^{k}\hat{U}^{(1)}}{\partial\hat{x}^{k}}\big{|}_{(x,y)=\bm{x}_{0}}, and turn to determine $\hat{U}^{(1),(k)}$, $k=0,1,2$. First, $\hat{U}^{(1),(0)}=\hat{U}^{(1)}(P_{0})$ can be computed with scheme (5.5a), i.e.,

[TABLE]

Furthermore, taking derivatives with respect to $\hat{x}$ and $\hat{y}$ on both sides of the above equation, we obtain

[TABLE]

Since $\hat{U}^{(1)}(P_{0})$ has already been obtained, and $\partial_{\hat{x}\hat{x}}F(\hat{U}^{n})(P_{0}),\partial_{\hat{x}\hat{y}}G(\hat{U}^{n})(P_{0}),\partial_{\hat{x}\hat{y}}F(\hat{U}^{n})(P_{0})$ and $\partial_{\hat{y}\hat{y}}G(\hat{U}^{n})(P_{0})$ can be approximated with $\hat{U}^{n}$ at interior grid points by the 2D WENO type extrapolation [26], then $\partial_{\hat{x}}\hat{U}^{(1)}(P_{0})$ and $\partial_{\hat{y}}\hat{U}^{(1)}(P_{0})$ can be solved out from the above equation. High order derivatives, $\partial_{\hat{x}\hat{x}}\hat{U}^{(1)}(P_{0})$, $\partial_{\hat{x}\hat{y}}\hat{U}^{(1)}(P_{0})$ and $\partial_{\hat{y}\hat{y}}\hat{U}^{(1)}(P_{0})$, are also computed by the 2D WENO type extrapolation [26]. Then $\hat{U}^{(1)}(P)$ can be computed with (5.7), and similarly $\hat{U}^{(k)}(P),k\geq 2$ can be obtained. The velocities at the ghost point can be transformed back to the original coordinate by

[TABLE]

while the scalar quantities remain unchanged in the two coordinates. In this way, both the solutions $U^{n}$ at time $t_{n}$ and intermediate solutions $U^{(i)}$ at ghost points are determined. Note that the 2D WENO type extrapolation has been been introduced detailedly in [26] and thus is omitted here.

6 Numerical validations

In this section we conduct several numerical experiments for both 1D and 2D problems to demonstrate the numerical stability and accuracy of our boundary treatment. In all the computations, we use the third-order finite difference WENO scheme [20] with the Lax-Friedrichs flux splitting to form the numerical fluxes [18] and the third-order IMEX RK scheme proposed in [7] with the double Butcher tableau

[TABLE]

[TABLE]

Additionally, we fix the CFL number to be $0.8$ in all the numerical examples. For all the 1D examples, the mesh for computation is taken as (3.1), namely, the boundary is located in the middle of two grid points. Different locations of the boundary are studied for 2D problems.

6.1 1D examples

Example 1. We first consider a 1D scalar hyperbolic equation

[TABLE]

in the domain $0\leq x\leq 1$ with analytical solution $u=\exp(t+x)$. Since $u>0$, the boundary condition only needs to be imposed at $x=0$. Specifically, the initial and boundary conditions are set as

[TABLE]

We compute the $L^{1},L^{2}$ and $L^{\infty}$ errors at the terminal time $t=1$, and the results are listed in Table 6.1. It can be seen that the convergence order is around three.

Example 2. Next we solve a 1D linear system

[TABLE]

in the domain $0\leq x\leq 1$. The Jacobian matrix has two eigenvalues $\pm 1$ and thus only one boundary condition should be imposed at $x=0$ and also $x=1$, respectively. We construct two sets of analytical solutions to (6.3):

[TABLE]

and

[TABLE]

There exists no boundary layer for the first solution (6.4) while the boundary layer occurs at $x=0$ for the second one (6.5) with small $\epsilon>0$. We test our method for the two cases by taking the initial data to match the analytical solution. The mesh is given as (3.2) and the solutions of $u$ at $t=1$ are used to evaluate the errors. In the second case, the three interior grid points near the left boundary are excluded in the computation of the errors due to the existence of the boundary layer.

For the solution (6.4), we specify $u$ and $v$ at the left and right boundaries, respectively, i.e.,

[TABLE]

Similarly, the boundary condition for (6.5) is given by

[TABLE]

The errors for the first case with solution (6.4) are given in Table 6.2, from which the desired third-order convergence can be observed. For the second case with solution (6.5), we take $\epsilon=1$ and $10^{-10}$ to test the convergence order. The results in Table 6.3 show that the present method of boundary treatment is third-order accurate for both choices of $\Delta x\ll\epsilon$ and $\Delta x\gg\epsilon$. Fig. 6.1 shows the solutions near the left boundary for different $\epsilon$ with $\Delta x=1/80$. It can be seen that our method behaves well even in the presence of boundary layer.

Example 3. We further construct the following nonlinear system

[TABLE]

in the domain $0\leq x\leq 1$ to test the accuracy of our method. The expressions (6.4) and (6.5) are also analytical solutions to (6.8) if cooperated with appropriate initial and boundary conditions. The Jacobian matrix of the flux is

[TABLE]

and it can be diagonalized as

[TABLE]

where $1+u+v>0$ and

[TABLE]

The boundary conditions in (6.6) and (6.7) are employed here.

The errors of the solution (6.4) with boundary condition (6.6) are given in Table 6.4 and the desired third-order accuracy is achieved again. Similar to the linear system, we take $\epsilon=1$ and $10^{-10}$ for the solution (6.5) with boundary condition (6.7) and the errors are shown in Table 6.5. It is clear that the convergence orders are about 3 for both choices of $\epsilon$. These results demonstrate the third-order accuracy of our method for nonlinear systems.

6.2 2D reactive Euler equations

All the following examples are devoted to validate our boundary treatment for the 2D reactive Euler equations (5.1). The parameters are fixed as in [27]:

[TABLE]

Example 4. We first verify the accuracy of the method through a problem with analytical solution

[TABLE]

in a square domain $[0,1]\times[0,1]$, subject to periodic boundary conditions at the left and right boundaries and $v=0$ at the top and bottom boundaries. The latter boundary condition is treated with our method. The mesh for computation is illustrated in Fig. 6.2. As shown in the figure, the distance between the bottom boundary and its adjacent gird line is $\eta\Delta x$ and that at the top boundary is $(1-\eta)\Delta x$. The left and right boundaries are placed in a similar way. We test the accuracy with $\eta=0.5,0.7$ and 0.9, and the results are listed in Table 6.6. It can be seen that the influence of $\eta$ on the errors is slight and the convergence order is about 3 for all the three choices of $\eta$. These demonstrate the third-order accuracy of our method for 2D problems.

Example 5. Next we test our method with the example in [27] containing discontinuity. The computational domain is $[0,2]\times[0,2]$ and the initial data is

[TABLE]

where $p$ and $Y$ are different from that in [27] to avoid possibly negative pressure or density. The boundary condition for the left and bottom boundaries are $u=0$ and $v=0$, respectively, and the solutions of the grid points outside the right and top boundaries are equal to those at the grid points inside the domain.

The mesh for computation is the same as that illustrated in Fig. 6.2 with a ratio $\eta$. We compare our results with those of the third-order positivity-preserving scheme in [17] with reflective boundary conditions (left and bottom) under the same mesh size $\Delta x=1/40$. The contour plots of $p$ and $\rho$ at $t=0.16$ are given in Fig. 6.3 and the solutions along the horizontal diagonal line $y=x$ are plotted in Fig. 6.4. Good agreement between the two methods can be observed. Additionally, we compare in Fig. 6.5 the solutions of different ratio $\eta$. It can be seen again that the influence of $\eta$ is small.

Example 6. We test the detonation diffraction in this example with the same configuration of boundaries in [27] (see Fig. 6.6 below). The spatial domain is $[0,5]\times[0,5]$ and the initial condition is

[TABLE]

which is constructed based on that in [27] to avoid negative pressure or density. The boundary condition is $\bm{n}\cdot\bm{u}=0$ with $\bm{n}$ being the normal direction of the boundary, except that at $x=0$, $(\rho,u,v,E,Y)=(11,6.18,0,970,1)$. The terminal time is $t=0.6$ and the mesh size is $\Delta x=1/20$. The ratio $\eta$ is fixed as 0.5 for all the boundaries.

The results of the present boundary treatment and the third-order positivity-preserving scheme in [17] with reflective boundary conditions are compared in Figs. 6.6 (contour plot) and 6.7 (solutions along $y=2.5$). It is clear that good agreements between the two methods are achieved again.

Example 7. Finally, we test our method through an example with multiple obstacles in [17], which is designed following [27]. In this example, the spatial domain is $[0,10]\times[0,10]$ and there are two obstacles with positions $[1,3]\times[0,3]$ and $[5,10]\times[0,5]$, respectively. The initial condition is set as (6.10) and the boundary condition is $\bm{n}\cdot\bm{u}=0$ for all the boundaries. We take $\Delta x=1/20$ and results of the present method and the third-order positivity-preserving scheme in [17] with reflective boundary conditions are shown in Fig. 6.8 at $t=1$ and Fig. 6.9 at $t=2$. It can be seen that the results of the two methods are very similar even with the very coarse mesh, which demonstrate the effectiveness of our boundary treatment for complex domains.

7 Conclusions and remarks

In this paper, we propose a high order finite difference boundary treatment method for the IMEX RK schemes solving hyperbolic systems with possibly stiff source terms on a Cartesian mesh. By combining the idea of using the RK schemes at the boundary and the ILW procedure, our method not only preserves the accuracy of the RK schemes but also processes good stability. Our method is different from the widely used approach for the explicit RK schemes by imposing boundary conditions for intermediate solutions [12], which could not be derived for the IMEX schemes. In addition, the intermediate boundary conditions are only available for explicit RK schemes up to third order while our method applies to arbitrary order IMEX and explicit RK schemes. Moreover, we show that our method applied to the third-order IMEX RK scheme solving 1D scalar equations can yield the widely used intermediate-stage boundary conditions in [12]. Finally, both 1D examples and 2D reactive Euler equations are used to demonstrate the good stability and third-order accuracy of our boundary treatment for a specific third-order IMEX scheme.

Since the IMEX RK schemes include explicit schemes as special cases, the present boundary treatment also applies to explicit RK schemes. It does not rely on intermediate-stage boundary conditions and thus is expected to be valid for explicit RK schemes of order higher than three. This will be investigated in our next work. Furthermore, our method may be adapted to IMEX RK schemes solving the other partial differential equations, *e.g. *convection-diffusion equations [28, 21, 29] and the Boltzmann equation [13, 16, 15, 19]. Many unresolved issues for the boundary treatment of RK schemes can be explored based on the idea of this work.