双曲守恒律方程的高精度ADER间断Galerkin方法

High Order ADER Discontinuous Galerkin Method for Hyperbolic Conservation Laws

摘要: 本文提出了一种全新的间断Galerkin (DG)方法，该方法使用单级ADER (任意时–空导数)方式进行时间离散。该方法利用微分变换步骤递归地将解的时–空展开系数通过低阶空间展开系数来表示，能够在空间和时间上达到任意高阶精度。与传统有限体积ADER格式相比较，该方法避免了在单元界面处求解广义Riemann问题。与多级Runge-Kutta DG (RKDG)方法相比较，由于不存在中间级，本方法需要较少的计算机内存。综上所述，所得到的方法是单步的、单级的、全离散的。最后，经典数值算例验证了该方法的良好性能：高精度、高分辨率、高效率。

Abstract: This article develops a new discontinuous Galerkin (DG) method using the one-stage ADER (Arbitrary DERivatives in time and space) approach for the temporal discretization. This current method employs the differential transformation procedure recursively to express the spatiotemporal expansion coefficients of the solution through the low order spatial expansion coefficients, which enables the method to easily achieve arbitrary high order accuracy in space and time. In comparison with the traditional ADER schemes, this method avoids solving the generalized Riemann problems at cell interfaces. Compared with the Runge-Kutta DG (RKDG) methods, the proposed method needs less computer memory storage due to no intermediate stages. In summary, the resulting method is one-step, one-stage, fully-discrete, and easily achieves arbitrary high order accuracy in space and time. Several examples illustrate the good performances of the present method: high order accuracy for smooth solutions, good resolution for discontinuous solutions and high efficiency.

1. 引言

在本文中，我们考虑一维双曲守恒律方程初值问题

$\{\begin{array}{l}{u}_t+f{\left(u\right)}_x=0,\\ u\left(x,0\right)=u_0\left(x\right),\end{array}x\in \left[a,b\right],$ (1)

的高精度数值方法。其中， $u=u\left(x,t\right)$ 为守恒变量， $f\left(u\right)=f\left(u\left(x,t\right)\right)$ 为物理流通量。该方程的特点是，即使初始条件充分光滑，随着时间的发展，解也会产生间断，因此对方法的要求非常高。间断Galerkin (DG)方法 [1] [2] [3] [4] 是求解守恒律方程的比较流行的一类方法。该方法的主要特点包括：高精度、方便处理复杂边界以及复杂区域。该方法是半离散方法，通常利用多级Runge-Kutta方法 [5] 进行时间离散，因此也被称作RKDG方法。但是，由于每一级均需要计算数值流通量和数值积分，所以RKDG方法计算成本较高。此外，由于多级的存在，该方法的计算机存储比较高。更进一步，频繁的内部信息交换阻碍了并行策略的有效实现。关于RKDG方法的简单回顾和最新研究进展，请参考文献 [6] [7]。

另外，Qiu、Dumbser、Shu [8] 提出一种借助于Lax-Wendroff时间离散方式的DG方法(LWDG) [8]。LWDG方法是单级的、显式的。此外，Dumbser和Munz提出了一种借助于单级ADER方法 [9] [10] [11] 进行时间离散的DG方法 [12] [13] [14]。随后，Dumber及其合作者针对该课题开展了一系列研究 [15] [16] [17] [18]。上述DG方法 [12] - [18] 的关键要素是在单元界面处构造高精度的数值流通量。单元界面处的未知量需要表示成关于时间的泰勒级数展开形式，其中时间导数通过反复使用控制方程本身并由空间导数来表示。这个步骤被称为Cauchy-Kowalewski步骤(在文献中也被称为Lax-Wendroff步骤)。此外，为了获得空间导数，还需要在单元界面处求解广义黎曼问题。然而，Cauchy-Kowalewski步骤是非常繁琐的，尤其对于高维问题。因此，替换或简化Cauchy-Kowalewski步骤是非常受欢迎的。值得注意的是，Dumbser和Munz [19] 提出了一种简化Cauchy-Kowalewski步骤的有效算法。最近，在DG方法的框架下，Dumber等人 [20]，以及汤华中等人 [21] 采用局部Galerkin预报指示子 [22] 来代替Cauchy-Kovalewski步骤。

2. 利用微分变换步骤构造ADER-DG方法

2.1. 相关符号

首先将空间计算区域 $\left[a,b\right]$ 离散为N个单元 $I_j=\left[x_{j-\frac{1}{2}},x_{j+\frac{1}{2}}\right],j=1,2,\cdots ,N$。其中 $a=x_{\frac{1}{2}}<x_{\frac{3}{2}}<\cdots <x_{N+\frac{1}{2}}=b$，并记单元中心、网格步长分别为 $x_j=\frac{1}{2}\left(x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}\right)$ 、 $\Delta x_j=x_{j+\frac{1}{2}}-x_{j-\frac{1}{2}}$。

假设现在的时间区间为 $\left[t^n,t^{n\text{+}1}\right]$，其中 $\Delta t=t^{n+1}-t^n$ 表示时间步长，然后记 ${\Omega }_j=I_j\times \left[t^n,t^{n+1}\right]$ 为时–空

控制单元。接下来，我们定义DG逼近空间如下：

$V_{\tau }^k=\left\{\varphi \left(x,t\right):{\varphi \left(x,t\right)|}_{{\Omega }_j}\in P^k\left({\Omega }_j\right),j=1,\cdots ,N\right\},$

其中 $P^k\left({\Omega }_j\right)$ 表示定义在时–空控制单元 ${\Omega }_j$ 上次数不超过k次的时–空多项式的集合。

2.2. ADER-DG方法的一般形式

首先在方程(1)两端同时乘一个任意空间测试函数 $\varphi \left(x\right)\in V_{\tau }^k$，并在时空控制单元 ${\Omega }_j$ 上使用分部积分，进而得到以下弱形式：

$\begin{array}{l}{\displaystyle {\int }_{I_j}u\left(x,t^{n+1}\right)\varphi \left(x\right)\text{d}x}-{\displaystyle {\int }_{I_j}u\left(x,t^n\right)\varphi \left(x\right)\text{d}x}+{\displaystyle {\int }_{t^n}^{t^{n+1}}f\left(x_{j+\frac{1}{2}},t\right)\text{d}t\cdot \varphi \left(x_{j+\frac{1}{2}}^{-}\right)}\\ -{\displaystyle {\int }_{t^n}^{t^{n+1}}f\left(x_{j-\frac{1}{2}},t\right)\text{d}t\cdot \varphi \left(x_{j-\frac{1}{2}}^{+}\right)}-{\displaystyle \underset{{\Omega }_j}{\iint }f\left(x,t\right){\varphi }_x\text{d}x\text{d}t}=0,j=1,\cdots ,N.\end{array}$

这里，二元函数 $f\left(x,t\right)=f\left(u\left(x,t\right)\right)$ 表示物理流通量。

假设未知解 $u\left(x,t\right)$ 在时–空控制单元 ${\Omega }_j$ 上的DG逼近为

$u_{\tau }\left(x,t\right)={\displaystyle \underset{k_t=0}{\overset{k}{\sum }}{\displaystyle \underset{k_x=0}{\overset{k-k_t}{\sum }}\tilde{u}\left(k_x,k_t\right){\left(x-x_j\right)}^{k_x}{\left(t-t^n\right)}^{k_t}}},$

这实际上是 $u\left(x,t\right)$ 的泰勒级数展开的截断形式。这里，我们的ADER-DG方法表述为：给定 $t^n$ 时刻的数值解 $u_{\tau }\left(x,t^n\right)\in V_{\tau }^K$，寻找下一时刻的数值解 $u_{\tau }\left(x,t^{n+1}\right)\in V_{\tau }^K$ 来满足以下等式：

$\begin{array}{l}{\displaystyle {\int }_{I_j}{u}_{\tau }\left(x,t^n\right)\varphi \left(x\right)\text{d}x}\\ ={\displaystyle {\int }_{I_j}{u}_{\tau }\left(x,t^n\right)\varphi \left(x\right)\text{d}x-{\hat{f}}_{j+\frac{1}{2}}\cdot \varphi \left(x_{j+\frac{1}{2}}^{-}\right)+{\hat{f}}_{j-\frac{1}{2}}\cdot \varphi \left(x_{j-\frac{1}{2}}^{+}\right)}+{\displaystyle \underset{{\Omega }_j}{\iint }{f}_{\tau }}\left(x,t\right){\varphi }_x\text{d}x\text{d}t,j=1,\cdots ,N,\end{array}$ (2)

这里 ${\hat{f}}_{j+\frac{1}{2}}$ 表示数值流通量，用来逼近单元界面处针对物理流通量在时间区间 $\left[t^n,t^{n\text{+}1}\right]$ 上的时间积分。

在本文中，我们采用简单高效的Lax-Friedrichs流通量：

${\hat{f}}_{j+\frac{1}{2}}=\frac{1}{2}{\displaystyle {\int }_{t^n}^{t^{n+1}}\left(f_{\tau }\left(x_{j\text{+}\frac{1}{2}}^{-},t\right)+f_{\tau }\left(x_{j+\frac{1}{2}}^{+},t\right)-\alpha \left(u_{\tau }\left(x_{j+\frac{1}{2}}^{+},t\right)-u_{\tau }\left(x_{j+\frac{1}{2}}^{-},t\right)\right)\right)\text{d}t}$ (3)

其中 $\alpha =\underset{u}{\max }\left|f^{\prime }\left(u\right)\right|$ 用来估计最大波速度，且最大值在整个空间计算区域上来获得。从公式(2)可以看出，如果想更新到时刻 $t^{n+1}$，我们需要构造单元界面处的数值流通量，计算空间积分，以及时–空积分。

2.3. 使用微分变换步骤计算时–空展开系数

传统ADER方法 [9] [10] [11] 一般采用Cauchy-Kowalewski步骤将解的时间导数通过其空间导数来表示。特别的，针对方程(1)的Cauchy-Kowalewski步骤如下：

$\begin{array}{l}{u}_t=-f^{\prime }\left(u\right)u_x,\\ u_{xt}=-f^{\prime }\left(u\right)u_{xx}-f^{″}\left(u\right){\left(u_x\right)}^2,\\ u_{tt}=-f^{″}\left(u\right)u_t{u}_x+f^{\prime }\left(u\right)f^{″}\left(u\right){\left(u_x\right)}^2+{\left(f^{\prime }\left(u\right)\right)}^2{u}_{xx},\\ ⋮\end{array}$

显然，对于高阶导数和高维问题，Cauchy-Kowalewski步骤将变得异常繁琐。

另外，文献 [23] 中提出了一种称为微分变换的替代步骤，用来简化Cauchy-Kowalewski步骤。在本文中，我们应用微分变换步骤来获得泰勒级数形式的数值解的时–空展开系数。同Cauchy-Kowalewski步骤相比较，微分变换步骤的实现非常简单有效 [24]。通常，时刻 $t^n$ 在单元 $I_j$ 中给定函数 $u\left(x,t\right)$，微分变换步骤用来获得如下系数：

$\tilde{u}\left(k_x,k_t\right)=\frac{1}{k_x!k_t!}{\frac{{\partial }^{k_x+k_t}u\left(x,t\right)}{\partial x^{k_x}\partial t^{k_t}}|}_{x=x_j,t=t^n},$

这里， $u\left(x,t\right)$ 是原始函数， $\tilde{u}\left(k_x,k_t\right)$ 是变换函数。实际上， $\tilde{u}\left(k_x,k_t\right)$ 表示函数 $u\left(x,t\right)$ 在点 $\left(x_j,t^n\right)$ 处的 $\left(k_x,k_t\right)$ 阶导数值。该步骤的逆运算为 $u\left(x,t\right)$ 的泰勒级数展开式：

$u\left(x,t\right)={\displaystyle \underset{k_x=0}{\overset{\infty }{\sum }}{\displaystyle \underset{k_t=0}{\overset{\infty }{\sum }}\tilde{u}\left(k_x,k_t\right)}}{\left(x-x_j\right)}^{k_x}{\left(t-t^n\right)}^{k_t}.$

下面列出了本文中用到的一些微分转换公式，具体见表1。

Table 1. Differential transformations of some functions encountered in this article

表1. 本文用到的有关微分变换公式

2.4. ADER-DG方法的详细步骤

基于微分变换步骤的ADER-DG方法的详细实现过程总结如下：

1) 从初始条件 $u\left(x,0\right)$ 出发，通过 $L^2$ 投影，得到 $\tilde{u}\left(k_x,0\right),k_x=0,1,\cdots ,k$，进而得到 $u_{\tau }\left(x,0\right)$ ；

2) 在时刻 $t=t^n$，在单元 $I_j$ 中，基于 $\tilde{u}\left(k_x,0\right)$ 和控制方程本身，并利用微分变换步骤，分别递归地得到

$\tilde{u}\left(k_x,k_t\right)$ 和 $\tilde{f}\left(k_x,k_t\right),k_t=0,1,\cdots ,k;k_x=0,1,\cdots ,k-k_t$，然后得到：

$u_{\tau }\left(x,t\right)={\displaystyle \underset{k_t=0}{\overset{k}{\sum }}{\displaystyle \underset{k_x=0}{\overset{k-k_t}{\sum }}\tilde{u}\left(k_x,k_t\right)}}{\left(x-x_j\right)}^{k_x}{\left(t-t^n\right)}^{k_t}\in {\Omega }_j,$

$f_{\tau }\left(x,t\right)={\displaystyle \underset{k_t=0}{\overset{k}{\sum }}{\displaystyle \underset{k_x=0}{\overset{k-k_t}{\sum }}\tilde{u}\left(k_x,k_t\right)}}{\left(x-x_j\right)}^{k_x}{\left(t-t^n\right)}^{k_t}\in {\Omega }_j.$

3) 通过公式(3)构造数值通量 ${\hat{f}}_{j+\frac{1}{2}}$。

4) 根据单步公式将数值解更新至下一个时间步 $t^{n+1}$，得到 $u_{\tau }\left(x,t^{n+1}\right)$。

5) 必要时针对解 $u_{\tau }\left(x,t^{n+1}\right)$ 使用斜率限制器。

注1：在使用微分变换步骤后，我们可以将数值解和流通量表示为时–空二元多项式。那么针对数值流通量和积分，我们可以采取精确计算，进而避免使用计算量较大的数值求积公式。

注2：实际上，微分变换步骤将控制方程本身转换为关于时–空泰勒级数的展开系数的递推关系。微分变换步骤使我们能够轻松地构造时–空上具有任意高阶精度的ADER-DG方法。此外，与Cauchy-Kowalewski步骤相比，微分变换步骤可以大大降低计算成本，并且应用范围更加广泛。

3. 数值结果

在本节中，我们将通过几个经典算例来验证该方法的良好性能。在整个计算过程中，我们分别采用时–空二次、四次多项式(即 $k=2$ 或者 $k=4$ )，并将CFL条件数分别取作0.18、0.1以保持数值稳定性。

3.1. 精度测试

我们首先通过如下拥有光滑解的初值问题

$\{\begin{array}{l}{u}_t+u_x=0,\hfill \\ u\left(x,0\right)=\sin \pi x,\hfill \end{array}x\in \left[0,2\right],$

来测试精度阶。该方程的微分变换的递推公式为

$\tilde{u}\left(k_x,k_t+1\right)=-\frac{k_x+1}{k_t+1}\tilde{u}\left(k_x+1,k_t\right).$

对于该算例，我仅仅取 $k=4$，并将时刻 $t=2$ 时的误差以及精度阶展示于表2。显然，我们获得了预期的五阶精度。

Table 2. Errors and order of accuracy

表2. 误差与精度阶

3.2. 线性对流方程

接下来考虑线性对流方程的初值问题：

$\{\begin{array}{ll}{u}_t+u_x=0,\hfill & x\in \left[-1,1\right],\hfill \\ u\left(x,0\right)=u_0\left(x\right),\hfill & 周期边界条件,\hfill \end{array}$

其中

$\begin{array}{l}{u}_0\left(x\right)=\{\begin{array}{ll}\frac{1}{6}\left(\vartheta \left(x,\beta ,z-\delta \right)+\vartheta \left(x,\beta ,z+\delta \right)+4\vartheta \left(x,\beta ,z\right)\right),\hfill & -0.8\le x\le -0.6,\hfill \\ 1,\hfill & -0.4\le x\le -0.2,\hfill \\ 1-\left|10\left(x-0.1\right)\right|,\hfill & 0\le x\le 0.2,\hfill \\ \frac{1}{6}\left(\Re \left(x,\alpha ,a-\delta \right)+\Re \left(x,\alpha ,a+\delta \right)+4\Re \left(x,\alpha ,a\right)\right),\hfill & 0.4\le x\le 0.6,\hfill \\ 0,\hfill & 其它,\hfill \end{array}\\ \vartheta \left(x,\beta ,z\right)={\text{e}}^{-\beta {\left(x-z\right)}^2},\\ \Re \left(x,\alpha ,a\right)=\sqrt{\max \left(1-{\alpha }^2{\left(x-a\right)}^2,0\right)}.\end{array}$

上面的常数分别取作

$a=0.5,z=-0.7,\delta =0.005,\alpha =10,\beta ={\log }_2/36{\delta }^2.$

该算例的初始条件包含非常复杂的形状。我们计算该算例至 $t=8$，三阶、五阶方法的数值结果见图1。

Figure 1. Solving the linear advection equation by third order (left) and fifth order (right) ADER-DG methods, t = 8

图1. 利用三阶(左)、五阶(右) ADER-DG方法求解线性对流方程t = 8

经过了长时间的计算，结果依然拥有非常好的分辨率。

3.3. 一维Burgers’方程

接下来，我们考虑一维非线性Burgers’方程：

$u_t+{\left(\frac{1}{2}{u}^2\right)}_x=0,$

这里， $f\left(u\right)=\frac{1}{2}{u}^2$ 表示物理流通量。该方程的微分变换的递推公式如下：

$\tilde{u}\left(k_x,k_t+1\right)=-\frac{k_x+1}{k_t+1}\tilde{f}\left(k_x+1,k_t\right),$

其中

$\tilde{f}\left(k_x,k_t\right)=\frac{1}{2}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\displaystyle \underset{s=0}{\overset{k_t}{\sum }}\tilde{u}\left(r,s\right)\tilde{u}}}\left(k_x-r,k_t-s\right).$

接下来，基于两种不同的初始条件，我们考虑以下两个算例：

算例1) $u\left(x,0\right)=\sin \left(\pi x\right),x\in \left[0,2\right]$ ；

算例2) $u\left(x,0\right)=\{\begin{array}{ll}-0.5\hfill & x\in \left[0,0.5\right],\hfill \\ 1\hfill & x\in \left(0.5,1\right],\hfill \\ 0\hfill & x\in \left(1,1.5\right],\hfill \end{array}x\in \left[0,1.5\right]$。

我们基于200个单元的网格分别计算这两个算例，数值结果见图2。数值结果能够非常好地拟合精确解，且具有较高的分辨率。

Figure 2. Solving the one-dimensional Bugers’ equation by third order (upper) and fifth order (lower) ADER-DG methods. Left: example 1), $t=1.5/\pi$ ; Right: example 2), $t=0.5$

图2. 利用三阶(上)、五阶(下) ADER-DG方法求解一维Bugers’方程。左：算例1)， $t=1.5/\pi$ ；右：算例2)， $t=0.5$

3.4. 二维Burgers’方程

最后，我们求解二维Burgers’方程：

$\{\begin{array}{l}{u}_t+{\left(\frac{u^2}{2}\right)}_x+{\left(\frac{u^2}{2}\right)}_y=0,\\ u\left(x,0\right)=0.5+\sin \left(\frac{\pi }{2}\left(x+y\right)\right),\end{array}\left[x,y\right]\in \left[0,4\right]\times \left[0,4\right].$

这里， $f\left(u\right)=g\left(u\right)=\frac{1}{2}{u}^2$ 分别表示x方向和y方向的物理流通量。另外，该方程的微分变换步骤的递推公式如下：

$\tilde{u}\left(k_x,k_y,k_t+1\right)=-\frac{k_x+1}{k_t+1}\tilde{f}\left(k_x+1,k_y,k_t\right)-\frac{k_y+1}{k_t+1}\tilde{g}\left(k_x,k_y+1,k_t\right)，$

其中

$\tilde{f}\left(k_x,k_y,k_t\right)=\tilde{g}\left(k_x,k_y,k_t\right)=\frac{1}{2}{\displaystyle \underset{r=0}{\overset{k_x}{\sum }}{\displaystyle \underset{s=0}{\overset{k_y}{\sum }}{\displaystyle \underset{h=0}{\overset{k_t}{\sum }}\tilde{u}}}}\left(r,s,h\right)\tilde{u}\left(k_x-r,k_y-s,k_t-h\right).$

我们基于包含100 × 100一致单元的网格将该算例计算至 $t=1.5/\pi$，数值结果见图3。显然，数值结果与精确解吻合的非常好，并且保持了较高的分辨率。

Figure 3. Solving the two-dimensional Bugers’ equation by third order (upper) and fifth order (lower) ADER-DG methods,. Left: the cutting-plot along; Right: three-dimensional bird’s eye view

图3. 利用三阶(上)、五阶(下) ADER-DG方法求解二维Burgers’方程，。左：沿着的截面图；右：三维鸟瞰图

最后，我们还针对算例3.2~算例3.4，关于三阶ADER-DG方法与三阶RKDG方法就CPU时间进行了比较，见表3。显然，ADER-DG方法能够节约非常可观的计算成本，拥有更高的效率，从而更加适应实际问题的需求。

Table 3. Comparison about the CPU time

表3. CPU时间的比较

4. 结论

在本文中，基于微分变换步骤，我们提出了一种全新的DG方法用于求解标量双曲守恒律方程，该方法使用单级ADER方式进行时间离散。与传统的多级RKDG方法相比，该方法是单步的、单级的，并且是全离散的。此外，微分转换步骤使得该方法能够比较容易地在时–空上达到任意高阶精度。关于将时–空导数转换为由空间导数来表示的操作，与Cauchy-Kowalewski步骤相比较，微分变换步骤更加简洁高效。而且，数值结果表明该方法针对光滑问题保持预期的高阶精度，针对间断解拥有较高的分辨率。与同阶RKDG方法相比较，该方法能节约较为可观的CPU时间。将来，我们将这种方法扩展到方程组以及二维问题。

致谢

本研究得到国家自然科学基金项目(11771228)的资助。

NOTES

^*通讯作者。