A high order bound preserving finite difference linear scheme for incompressible flows

TL;DR

This paper introduces a high-order finite difference scheme with a bound-preserving flux limiter for incompressible flows, offering high resolution and reduced dissipation compared to WENO schemes, suitable for oscillatory problems.

Contribution

A novel high-order linear scheme with a maximum-principle-preserving flux limiter for incompressible flows, combining efficiency, high resolution, and controlled oscillations.

Findings

Less dissipative than WENO schemes

High resolution due to Hermite reconstruction

Effective control of numerical oscillations

Abstract

We propose a high order finite difference linear scheme combined with a high order bound preserving maximum-principle-preserving (MPP) flux limiter to solve the incompressible flow system. For such problem with highly oscillatory structure but not strong shocks, our approach seems to be less dissipative and much less costly than a WENO type scheme, and has high resolution due to a Hermite reconstruction. Spurious numerical oscillations can be controlled by the MPP flux limiter. Numerical tests are performed for the Vlasov-Poisson system, the 2D guiding-center model and the incompressible Euler system. The comparison between the linear and WENO type schemes will demonstrate the good performance of our proposed approach.

Figures

Tables

Table 1. Table 3.1. Accuracy test for the 2D linear transport equation in Example 3.1 . ρ m ​ i ​ n subscript 𝜌 𝑚 𝑖 𝑛 \rho_{min} is the minimum numerical solution. “WO” stands for “without limiter”, “WL” stands for “with limiter”. T = 1 𝑇 1 T=1 .

	N	$L^1$ error	order	$L^{\infty }$ error	order	${\rho }_{min}$
WO	32	5.00e-04	–	1.29e-03	–	-0.001019
	64	1.90e-05	4.72	4.91e-05	4.72	-4.244e-05
	128	6.41e-07	4.89	1.68e-06	4.87	-1.242e-06
	256	2.05e-08	4.97	5.35e-08	4.97	-4.403e-08
WL	32	5.16e-04	–	1.33e-03	–	0.0003958
	64	1.90e-05	4.77	5.54e-05	4.59	5.522e-06
	128	6.41e-07	4.89	2.07e-06	4.74	9.832e-13
	256	2.05e-08	4.97	7.78e-08	4.73	1.712e-09

Table 2. Table 3.2. L 1 superscript 𝐿 1 L^{1} and L ∞ superscript 𝐿 L^{\infty} errors and orders for the VP system with initial condition ( 3.8 ) at T = 1 𝑇 1 T=1 . “WL” denotes the scheme with limiters, “WO” denotes the scheme without limiters. “ f m ​ i ​ n subscript 𝑓 𝑚 𝑖 𝑛 f_{min} ” is the minimum of the numerical solution. Mesh size N v = 2 ​ N x subscript 𝑁 𝑣 2 subscript 𝑁 𝑥 N_{v}=2N_{x} .

	$N_x$	$L^1$ error	order	$L^{\infty }$ error	order	$f_{min}$
WO	64	1.66e-05	–	0.0004144	–	-1.317e-06
	128	6.83e-07	4.60	1.913e-05	4.44	-2.197e-08
	256	2.54e-08	4.75	7.137e-07	4.75	-1.777e-09
	512	8.29e-10	4.94	2.342e-08	4.93	-5.535e-11
WL	64	1.68e-05	–	0.0004177	–	1.59e-32
	128	6.85e-07	4.62	1.915e-05	4.45	1.768e-34
	256	2.54e-08	4.75	7.136e-07	4.75	2.339e-35
	512	8.30e-10	4.94	2.342e-08	4.93	4.797e-36

Table 3. Table 3.3. Example 3.7 . ω m ​ i ​ n subscript 𝜔 𝑚 𝑖 𝑛 \omega_{min} and ω m ​ a ​ x subscript 𝜔 𝑚 𝑎 𝑥 \omega_{max} are the minimum and maximum numerical solutions respectively. “WO” stands for “without limiter”, “WL” stands for “with limiter”. T = 1 𝑇 1 T=1 .

	N	$L^1$ error	order	$L^{\infty }$ error	order	${\omega }_{max}$	${\omega }_{min}$
WO	32	2.34e-05	–	4.89e-05	–	2	-2
	64	8.86e-07	4.72	1.63e-06	4.91	2	-2
	128	2.93e-08	4.92	5.08e-08	5.00	2	-2
	256	9.43e-10	4.96	1.60e-09	4.99	2	-2
WL	32	2.34e-05	–	4.89e-05	–	2	-2
	64	8.86e-07	4.72	1.63e-06	4.91	2	-2
	128	2.93e-08	4.92	5.08e-08	5.00	2	-2
	256	9.43e-10	4.96	1.60e-09	4.99	2	-2

Equations

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\rho}{\partial t}\,+\,{\bf U}\cdot\nabla\rho\,=\,0,\\ \\ {\rm div}\,{\bf U}\,=\,0,\end{array}\right.$$

$$\frac{\partial f}{\partial t}+\mathbf{v}\cdot {\nabla }_{\mathbf{x}}f+\mathbf{E}\cdot {\nabla }_{\mathbf{v}}f=0,$$

$$\mathbf{E}(t,\mathbf{x})=-{\nabla }_{\mathbf{x}}\Phi (t,\mathbf{x}),-{\Delta }_{\mathbf{x}}\Phi (t,\mathbf{x})=\rho (t,\mathbf{x}).$$

$$\rho (t,\mathbf{x})={\int }_{{\mathbb{R}}^d}f(t,\mathbf{x},\mathbf{v})d\mathbf{v}.$$

$$\left\{\begin{array}[]{l}\displaystyle\frac{\partial\rho}{\partial t}\,+\,{\bf U}\cdot\nabla\rho\,=\,0,\\ \\ -\Delta\,\Phi\,=\,\rho,\end{array}\right.$$

$${\partial }_t\rho +{\mathrm{div}}_{\mathbf{x}}\left(\mathbf{U}\rho \right)=0,$$

$$\begin{array}[]{l}x_{min}\,=\,x_{0}<x_{1}<\cdots<x_{N_{x}-1}<x_{N_{x}}\,=\,x_{max},\\ \\ y_{min}\,=\,y_{0}<y_{1}<\cdots<y_{N_{y}-1}<y_{N_{y}}\,=\,y_{max},\end{array}$$

$${\rho }_{i,j}^{n+1}={\rho }_{i,j}^n-\Delta t\left(\frac{{\hat{H}}_{i+\frac{1}{2},j}-{\hat{H}}_{i-\frac{1}{2},j}}{\Delta x}+\frac{{\hat{G}}_{i,j+\frac{1}{2}}-{\hat{G}}_{i,j-\frac{1}{2}}}{\Delta y}\right),$$

$${\partial }_t\rho +{\partial }_x(U\rho )=0,$$

$${\rho }_i^{n+1}={\rho }_i^n-\frac{\Delta t}{\Delta x}\left({\hat{h}}_{i+\frac{1}{2}}-{\hat{h}}_{i-\frac{1}{2}}\right),$$

$$\hat{h}_{i+\frac{1}{2}}\,=\,\left\{\begin{array}[]{l}\displaystyle h^{-}_{i+\frac{1}{2}},\qquad\text{ if }\frac{U^{n}_{i}+U^{n}_{i+1}}{2}>0,\\ \\ h^{+}_{i+\frac{1}{2}},\qquad\text{ otherwise }.\end{array}\right.$$

$$h_{i+\frac{1}{2}}^{-}=\frac{1}{27}\left(-8h_{i-1}+19h_i+19h_{i+1}+3G_{i-3/2}^{\prime }-6G_{i+3/2}^{\prime }\right),$$

$$G^{\prime}_{i+\frac{1}{2}}\,=\,\frac{1}{60}\Big{[}(h_{i+3}-h_{i-2})-8(h_{i+2}-h_{i-1})+37(h_{i+1}-h_{i})\Big{]},$$

$$h_{i+\frac{1}{2}}^{-}={\omega }_l{h}_l(x_{i+\frac{1}{2}})+{\omega }_c{h}_c(x_{i+\frac{1}{2}})+{\omega }_r{h}_r(x_{i+\frac{1}{2}}).$$

$$\left\{\begin{array}[]{l}\displaystyle h_{l}(x_{i+\frac{1}{2}})=-2h_{i-1}+2h_{i}+G^{\prime}_{i-\frac{3}{2}},\\ \\ \displaystyle h_{c}(x_{i+\frac{1}{2}})=\frac{-h_{i-1}+5h_{i}+2h_{i+1}}{6},\\ \\ \displaystyle h_{r}(x_{i+\frac{1}{2}})=\frac{h_{i}+5h_{i+1}-2G^{\prime}_{i+\frac{3}{2}}}{4},\end{array}\right.$$

$${\omega }_k=\frac{{\alpha }_k}{{\alpha }_l+{\alpha }_c+{\alpha }_r},{\alpha }_k=\frac{c_k}{(ϵ+{\beta }_k{)}^2},k=l,c,r.$$

$${\beta }_l$$

$$=\frac{13}{16}{s}_1^2+\frac{3}{16}(s_1-4s_2{)}^2,\text{ with }{s}_1=h_{i-1}-h_i,s_2=-3h_{i-1}+h_0+G_{i-\frac{3}{2}}^{\prime },$$

$${\beta }_c$$

$$=\frac{1}{4}{s}_1^2+\frac{13}{12}{s}_2^2,\text{ with }{s}_1=h_{i+1}-h_{i-1},s_2=h_{i+1}-2h_0+h_{i-1},$$

$${\beta }_r$$

$$=\frac{13}{16}{s}_1^2+\frac{3}{16}(s_1-4s_2{)}^2,\text{ with }{s}_1=h_{i+1}-h_i,s_2=-3h_{i+1}+h_0+G_{i+\frac{3}{2}}^{\prime }.$$

$${\rho }^{n+1}={\rho }^n+\Delta tL({\rho }^n),$$

$$\left\{\begin{array}[]{l}\rho^{(1)}\,=\,\rho^{n}+\frac{1}{2}\,\Delta t\,L(\rho^{n}),\\ \\ \rho^{(2)}\,=\,\rho^{n}+\frac{1}{2}\,\Delta t\,L(\rho^{(1)}),\\ \\ \rho^{(3)}\,=\,\rho^{n}+\Delta t\,L(\rho^{(2)}),\\ \\ \rho^{n+1}\,=\,\rho^{n}+\frac{1}{6}\,\Delta t\,\left(L(\rho^{n})+2L(\rho^{(1)})+2L(\rho^{(2)})+L(\rho^{(3)})\right),\end{array}\right.$$

$${\rho }_{i,j}^{n+1}={\rho }_{i,j}^n-\Delta t\left(\frac{{\hat{H}}_{i+\frac{1}{2},j}^{rk}-{\hat{H}}_{i-\frac{1}{2},j}^{rk}}{\Delta x}+\frac{{\hat{G}}_{i,j+\frac{1}{2}}^{rk}-{\hat{G}}_{i,j-\frac{1}{2}}^{rk}}{\Delta y}\right),$$

$${\hat{H}}^{rk}=\frac{1}{6}\left({\hat{H}}^{(0)}+2{\hat{H}}^{(1)}+2{\hat{H}}^{(2)}+{\hat{H}}^{(3)}\right),$$

$$\{\begin{array}{l}{\partial }_t\rho +{\mathrm{div}}_{\mathbf{x}}(\mathbf{U}\rho )=0,\\ {\mathrm{div}}_{\mathbf{x}}\mathbf{U}=0,\end{array}$$

$$\{\begin{array}{l}{\rho }_{i,j}^{n+1}={\rho }_{i,j}^n-{\lambda }_x\left(U_{i+\frac{1}{2},j}^{-}{\rho }_{i,j}^n+U_{i+\frac{1}{2},j}^{+}{\rho }_{i+1,j}^n-U_{i-\frac{1}{2},j}^{-}{\rho }_{i-1,j}^n-U_{i-\frac{1}{2},j}^{+}{\rho }_{i,j}^n\right)\\ -{\lambda }_y\left(U_{i,j+\frac{1}{2}}^{-}{\rho }_{i,j}^n+U_{i,j+\frac{1}{2}}^{+}{\rho }_{i,j+1}^n-U_{i,j-\frac{1}{2}}^{-}{\rho }_{i,j-1}^n-U_{i,j-\frac{1}{2}}^{+}{\rho }_{i,j}^n\right),\\ \frac{1}{\Delta x}\left(U_{i+\frac{1}{2},j}^{-}+U_{i+\frac{1}{2},j}^{+}-U_{i-\frac{1}{2},j}^{-}+U_{i-\frac{1}{2},j}^{+}\right)+\frac{1}{\Delta y}\left(U_{i,j+\frac{1}{2}}^{-}+U_{i,j+\frac{1}{2}}^{+}-U_{i,j-\frac{1}{2}}^{-}-U_{i,j-\frac{1}{2}}^{+}\right)=0,\end{array}$$

$$\Delta t\le \frac{1}{2\max |\mathbf{U}|}\frac{\Delta x\Delta y}{\Delta x+\Delta y}$$

$$U_{i+\frac{1}{2},j}^{-}\ge 0,U_{i+\frac{1}{2},j}^{+}<0,U_{i,j+\frac{1}{2}}^{-}\ge 0,U_{i,j+\frac{1}{2}}^{+}<0,$$

$${\rho }_{i,j}^{n+1}={\alpha }_{i,j}{\rho }_{i,j}^n+{\alpha }_{i+1,j}{\rho }_{i+1,j}^n+{\alpha }_{i-1,j}{\rho }_{i-1,j}^n+{\alpha }_{i,j+1}{\rho }_{i,j+1}^n+{\alpha }_{i,j-1}{\rho }_{i,j-1}^n,$$

$$\{\begin{array}{l}{\alpha }_{i,j}=1-{\lambda }_x\left(U_{i+\frac{1}{2},j}^{-}-U_{i-\frac{1}{2},j}^{+}\right)-{\lambda }_y\left(U_{i,j+\frac{1}{2}}^{-}-U_{i,j-\frac{1}{2}}^{+}\right),\\ {\alpha }_{i+1,j}=-{\lambda }_x{U}_{i+\frac{1}{2},j}^{+},{\alpha }_{i-1,j}={\lambda }_x{U}_{i-\frac{1}{2},j}^{-},\\ {\alpha }_{i,j+1}=-{\lambda }_y{U}_{i,j+\frac{1}{2}}^{+},{\alpha }_{i,j-1}={\lambda }_y{U}_{i,j-\frac{1}{2}}^{-}.\end{array}$$

$${\alpha }_{i,j}+{\alpha }_{i+1,j}+{\alpha }_{i-1,j}+{\alpha }_{i,j+1}+{\alpha }_{i,j-1}=1,$$

Full text

Keywords. High order linear scheme; MPP flux limiter; Incompressible transport equations.

1. 1 Introduction

2. 2 Finite difference scheme with Hermite reconstruction

3. 2.1 Conservative finite difference scheme

4. 2.2 Hermite linear reconstruction

5. 2.3 Hermite WENO reconstruction

6. 2.4 High order Runge-Kutta time discretization

7. 2.5 Bound preserving MPP flux limiter

8. 3 Numerical examples

9. 3.1 Linear transport equations

10. 3.2 Vlasov-Poisson system

11. 3.3 Kelvin-Helmholtz instability

12. 3.4 Incompressible Euler equations

13. 4 Conclusion

1. Introduction

In this paper, we are interested in the numerical approximation of incompressible transport equations as

[TABLE]

where ${\bf U}$ represents the advection field and $\rho$ is a nonnegative density.

The typical example of application is the well known Vlasov-Poisson (VP) system arising in collisionless plasma physics. It describes the time evolution of particles under the effects of self-consistent electrostatic field and reads

[TABLE]

$f\,:=\,f(t,{\bf x},{\bf v})$ is the distribution function in the phase space $({\bf x},{\bf v})\in\mathbb{R}^{d}\times\mathbb{R}^{d},\,\,d=1,2,3$. ${\bf E}\,:=\,{\bf E}(t,{\bf x})$ is the electric field, which can be determined by the Poisson’s equation from an electric potential function $\Phi(t,{\bf x})$

[TABLE]

The charge density $\rho(t,{\bf x})$ is defined as

[TABLE]

Another example is the two dimensional guiding-center model, which describes the evolution of the charge density $\rho$ in a highly magnetized plasma in the transverse plane of a tokomak, is given by

[TABLE]

where ${\bf U}\,=\,(-\Phi_{y},\Phi_{x})$ is a divergence free velocity. The two dimensional guiding center model can also be referred as an asymptotic model of the VP system by averaging in the velocity phase space, for details, see [22]. We notice that the guiding-center model (1.4) is in the same form as the two dimensional incompressible Euler equations in the vorticity stream function formulation, which describes the evolution of vortices in fluid hydrodynamics.

For all the models we mention above, they all have an transport equation coupled with a Poisson’s equation for the advection velocity, and moreover, the advection velocity is divergence free. In the following, we refer them as incompressible flow models.

Many numerical schemes have been proposed for solving these models, especially recently, high order schemes are very attractive due to their high resolutions for such problems with rich solution structures. For example, deterministic methods, there are finite difference, finite volume and finite element Eulerian methods [8, 37, 32, 40, 6, 14, 15, 41, 13, 39, 38], semi-Lagrangian methods [3, 34, 5, 25, 33, 18, 10, 26, 27, 28, 9, 24, 4, 31], and discontinuous Galerkin finite element methods [17, 42, 23, 7, 11, 19, 2], also see many other references therein. However, due to the highly oscillatory structure of such problems, linear type schemes for these problems would show significant spurious numerical oscillations, which might get worse with increased orders. Weighted essentially nonoscillotry (WENO) reconstruction, which was originally developed in the presence of both shocks and small fine structures for fluid hydrodynamics of hyperbolic conservation laws, see for example [29], is frequently adopted in most of the finite difference and finite volume Eulerian or semi-Lagrangian methods mentioned above. For WENO type schemes, we can often observe the good performance of the WENO reconstruction on suppressing spurious oscillations [25, 33, 37, 32, 40], which can also be seen from our examples in the numerical section. However, for incompressible flow problems, their solutions are highly oscillatory but without strong discontinuities. We might expect excessive usage of the WENO reconstruction, which is computationally more expensive than a pure linear type scheme and too dissipative for certain classes of problems [20], especially for high dimensional problems with long time simulations. Although hybrid approaches of coupling a linear scheme and a WENO scheme can save some cost from the WENO reconstruction, for example the very recent work [12] and references therein, they still focus on problems with strong discontinuities and WENO may not be avoided.

In this paper, we propose to solve the incompressible flow problems with a high order linear scheme without WENO reconstruction. Linear schemes have the following several good properties: (1) less dissipative: for example, it preserves the $L^{2}$ norm (also energy and entropy for the VP system) better than the WENO type scheme; (2) less costly and easier implementation: without WENO reconstruction, it saves a lot computational cost and can be easily implemented, especially when extended to high dimensional problems; (3) higher resolution if with a Hermite reconstruction. We adopt the scheme in [37] but with a Hermite linear reconstruction. In order to control the spurious numerical oscillations due to a linear type scheme, we seek to combine it with a newly developed high order bound preserving maximum-principle-preserving (MPP) flux limiter. The MPP flux limiter was first proposed by Xu et. al. [35, 21], and then improved by Xiong et. al. [32, 33]. It can be seen as a very weak limiter, which just pulls the numerical overshootings and undershootings back to its physical range, without excessive dissipating the solution within the range. Moreover, there is no further time step restriction on this MPP flux limiter from the original scheme for linear stabilities. We refer to [36] for a review of recent works on this bound preserving high order flux limiter. The coupling of the linear scheme and the MPP flux limiter keeps the original high order accuracy while maintaining a large CFL number. Due to the highly oscillatory but non discontinuous solutions of incompressible flow problems, the MPP flux limiter can serve as a necessary auxiliary tool for the linear scheme. The extra work from applying the MPP flux limiter at each final stage of a multi-stage Runge-Kutta time discretization, is much less than the WENO reconstruction. Numerical experiments, especially the bump-on-tail instability of long time simulation and the Kelvin-Helmoholtz instability of the 2D guiding-center model are as benchmark tests, will be performed to demonstrate the good performance of our proposed approach.

The rest of the paper is organized as follows. In Section 2, we will describe the conservative finite difference scheme with both linear and WENO reconstructions, for the completeness of comparison in the numerical section. The bound preserving MPP flux limiter will also be briefly reviewed. In Section 3, numerical experiments including the VP system, the Kelvin-Helmholtz instability of the 2D guiding center model and the incompressible Euler system will be studied. Conclusions are made in Section 4.

2. Finite difference scheme with Hermite reconstruction

In this section, we will describe our scheme to solve the VP system (1.2)-(1.3), the 2D guiding-center model (1.4) as well as the incompressible Euler system. Here we just consider $d=1$ for the VP system, in common we will have a 2D transport equation with a divergence free velocity, coupled with a Poisson’s type equation, which is 1D for the VP system, and 2D for the other two. In the following, we will take the 2D guiding-center model (1.4) as an example to describe our schemes. The other two models can be applied similarly.

We propose a finite difference scheme with a Hermite linear reconstruction for solving the 2D conservative transport equation. We adopt the scheme developed by Filbet and Yang [37], which has a Hermite WENO reconstruction. In the following, we will briefly recall the 2D conservative finite difference scheme and describe both the Hermite linear and WENO reconstructions. We note that the smooth indicators for the Hermite WENO reconstruction are modified as compared to [37]. The Poisson’s equation for the electric potential function $\Phi$ will be solved by fast Fourier transform (FFT) for periodic boundary condition on an interval in one-dimension (1D) or periodic boundary conditions on a rectangular domain in two-dimension (2D), which will be omitted here. We refer to [37] for more details.

2.1. Conservative finite difference scheme

We consider the 2D transport equation in a conservative form

[TABLE]

with ${\bf U}={\bf U}(t,{\bf x})$ such that ${\rm div}_{\bf x}{\bf U}=0$ and ${\bf x}=(x,y)$. For simplicity, we assume a uniform discretization of the computational domain $[x_{min},x_{max}]\times[y_{min},y_{max}]$ with $N_{x}\times N_{y}$ grid points

[TABLE]

where the mesh sizes are $\Delta x=x_{i+1}-x_{i}$ and $\Delta y=y_{j+1}-y_{j}$ for $0\leq i\leq N_{x},\,\,0\leq j\leq N_{y}$. A conservative finite difference scheme with Euler forward time discretization is defined as follows:

[TABLE]

where the time step is $\Delta t=t^{n+1}-t^{n}$. $\rho^{n}_{i,j}$ is the numerical value at time level $t^{n}$ on the grid point $(x_{i},y_{j})$. $\hat{H}_{{i+\frac{1}{2}},j}$, $\hat{G}_{i,{j+\frac{1}{2}}}$ are the numerical fluxes in the $x$ and $y$ directions respectively.

2.2. Hermite linear reconstruction

For a finite difference scheme (2.1), the numerical fluxes $\hat{H}_{{i+\frac{1}{2}},j}$ and $\hat{G}_{i,{j+\frac{1}{2}}}$ are reconstructed dimension by dimension. Here we will take a 1D transport equation to illustrate on how to obtain a flux by a Hermite linear reconstruction. $\hat{H}_{{i+\frac{1}{2}},j}$ and $\hat{G}_{i,{j+\frac{1}{2}}}$ are obtained in this way along each of its own direction. The Hermite linear reconstruction is what we propose in this paper. A corresponding Hermite WENO reconstruction will be described in the next subsection, which is used for the comparison in the numerical section.

Let us consider a prototype 1D conservative transport equation

[TABLE]

with velocity $U=U(t,x)$ and a uniform discretization with mesh size $\Delta x=x_{k+1}-x_{k}$ for $0\leq k<N_{x}$. A conservative finite difference scheme for (2.2) can be written as

[TABLE]

where $\rho^{n}_{i}$ is the numerical point value at time level $t^{n}$ on the grid point $x_{i}$. $\hat{h}_{i+\frac{1}{2}}$ can be chosen as an upwind numerical flux, which is

[TABLE]

$h^{\pm}_{i+\frac{1}{2}}$ are fluxes reconstructed from $\{h^{n}_{i}=U^{n}_{i}\rho^{n}_{i}\}_{i}$ by a Hermite linear reconstruction, from the left and right sides of $x_{i+\frac{1}{2}}$ respectively, where $U^{n}_{i}$ is the numerical velocity approximating $U(t^{n},x_{i})$. For simplicity, we drop the superscript $n$ for $h^{n}_{i}$ and $h^{-}_{{i+\frac{1}{2}}}$ is simply reconstructed by

[TABLE]

where the derivative of the primitive function $G^{\prime}_{i+\frac{1}{2}}$ is given by a 6th order central difference approximation

[TABLE]

$h^{+}_{{i+\frac{1}{2}}}$ can be obtained in mirror symmetric with respect to $x_{i+\frac{1}{2}}$.

2.3. Hermite WENO reconstruction

Now we will describe a fifth order Hermite WENO reconstruction to compute $h^{-}_{i+\frac{1}{2}}$ corresponding to (2.3). The procedure is outlined as follows. $h^{+}_{i+\frac{1}{2}}$ can also be obtained in mirror symmetric with respect to $x_{i+\frac{1}{2}}$. We drop the superscript $n$ for $h^{n}_{i}$ and we have

[TABLE]

The three polynomials $h_{l}(x)$, $h_{c}(x)$ and $h_{r}(x)$ evaluating at $x_{i+\frac{1}{2}}$ are

[TABLE]

whereas the weights $\omega_{l}$, $\omega_{c}$ and $\omega_{r}$ are the nonlinear WENO weights and determined according to the smoothness indicators

[TABLE]

In order to match (2.3), the linear coefficients are $c_{l}=1/9$ and $c_{c}=c_{r}=4/9$. The small parameter $\epsilon=10^{-6}$ is to avoid the denominator to be [math].

To evaluate the smooth indicators $\beta_{l}$, $\beta_{c}$ and $\beta_{r}$, we measure them on the cell $[x_{i-\frac{1}{2}},x_{i+\frac{1}{2}}]$ instead of $[x_{i},x_{i+1}]$ as in [37]. In this way, the smooth indicators are symmetric with respect to $x_{i}$, as we can see below:

[TABLE]

2.4. High order Runge-Kutta time discretization

The first order Euler forward time discretization (2.1) can be generalized to high order Runge-Kutta (RK) time discretization [30]. For example, if we write (2.1) in the following form

[TABLE]

here the subscripts for $\rho$ are dropped and the operator $L$ denotes the spatial discretization, then a 4th order RK time discretization can be written as

[TABLE]

with a CFL number $\frac{2}{3}$ for linear stability [30]. The last stage of (2.6) can be written in the same form as (2.1), which is

[TABLE]

where the numerical flux $\hat{H}^{rk}_{{i+\frac{1}{2}},j}$ is accumulated from all formal stages, that is

[TABLE]

where $\hat{H}^{(l)}$ is the numerical flux at the corresponding $l$-th stage, and the subscript $({i+\frac{1}{2}},j)$ is dropped for clarity. $\hat{G}^{rk}_{i,{j+\frac{1}{2}}}$ is defined similarly.

2.5. Bound preserving MPP flux limiter

For a high order bound preserving limiter, we adopt the parametrized maximum-principle-preserving (MPP) flux limiter developed in [32, 33]. The MPP flux limiter is only applied at the final stage of (2.6), that is (2.7). Due to the same form of (2.1) and (2.7), the MPP flux limiter is defined as a convex combination of a first order monotone flux and the high order flux in (2.7). In the following, we first describe a first order monotone MPP scheme for a general incompressible flow system, then we recall a specific first order monotone MPP scheme developed in [32] for the incompressible flow system with the Poisson’s equation.

For a general incompressible flow system with divergence free condition

[TABLE]

a first order scheme can be defined as

[TABLE]

where $\lambda_{x}=\Delta t/\Delta x$ and $\lambda_{y}=\Delta t/\Delta y$. If the scheme (2.9) is a consistent discretization to (2.8), we have the following statement:

Proposition 2.1**.**

For a first order consistent scheme (2.9) solving the incompressible system (2.8), it is monotone and MPP, if the time step is small enough, e.g.,

[TABLE]

and the coefficients have signs that

[TABLE]

for all $i,j$.

Proof.

We can rewrite the first equation in the scheme (2.9) to be

[TABLE]

with

[TABLE]

It is easy to check that $\alpha_{i,j},\alpha_{i\pm 1,j},\alpha_{i,j\pm 1}$ are all positive if the two conditions (2.10) and (2.11) are satisfied. Moreover,

[TABLE]

due to the enforced discrete divergence free condition in (2.9). Since (2.12) is a convex combination, so it is monotone and MPP. ∎

A specific first order monotone MPP scheme depends on how to choose $U^{\pm}_{{i+\frac{1}{2}},j}$ and $U^{\pm}_{i,{j+\frac{1}{2}}}$. By using a potential function $\Phi$, where ${\bf U}=(-\Phi_{y},\Phi_{x})$, a first order monontone MPP scheme with Lax-Friedrichs flux defined in [32] is to take:

[TABLE]

where

[TABLE]

In this paper, the potential function $\Phi$ is obtained by solving the Poisson’s equation $\Delta\Phi=\rho$ with FFT.

Now to describe the MPP flux limiter, we write the first order monotone MPP scheme in a conservative flux difference form

[TABLE]

with

[TABLE]

$\hat{h}_{{i+\frac{1}{2}},j}$ and $\hat{g}_{i,{j+\frac{1}{2}}}$ are first order monotone numerical fluxes. Similarly, (2.7) is

[TABLE]

In order to ensure maximum principle, we are looking for type of limiters

[TABLE]

such that

[TABLE]

(2.17) and (2.18) form coupled inequalities for the limiting parameters $\theta_{i+1/2,j},\theta_{i,j+1/2}$. We will need to find a group of numbers $\Lambda_{L,i,j},\Lambda_{R,i,j},\Lambda_{D,i,j},\Lambda_{U,i,j}$, such that (2.18) is satisfied with

[TABLE]

This is achieved in a decoupled way for the minimum and maximum parts. For the maximum value case, we let

[TABLE]

and denote

[TABLE]

then (2.18) can be rewritten as

[TABLE]

We shall now focus on decoupling the inequalities (2.21): for the single node $(i,j)$,

1. (1)

Identify positive values out of the four locally defined numbers $F_{i-1/2,j}$, $F_{i+1/2,j}$, $F_{i,j-1/2}$, $F_{i,j+1/2}$; 2. (2)

Corresponding to those positive values, collectively, the limiting parameters can be defined. For example, if $F_{i+1/2,j},F_{i-1/2,j}>0$ and $F_{i,j-1/2},F_{i,j+1/2}\leq 0$, then

[TABLE]

that is, all high order fluxes which possibly contribute (beyond that of the first order fluxes, which is not good) to the overshooting or undershooting of the updated value shall be limited by the same scaling. Similarly we can find $\Lambda^{M}_{i,j\pm 1/2}$ and also similarly for the minimum value case of $\Lambda^{m}_{i\pm 1/2,j}$ and $\Lambda^{m}_{i,j\pm 1/2}$. Therefore the range of the limiting parameters satisfying MPP for a single node $(i,j)$ is defined by

[TABLE]

Considering the limiters from neighboring nodes, finally we let

[TABLE]

Substituting (2.24) into (2.17), our final scheme is (2.16) by replacing the fluxes from (2.17).

3. Numerical examples

In this section, we will apply the 5th order finite difference scheme with Hermite linear reconstruction (2.3) and (2.4), coupled with the bound preserving MPP flux limiter, which we denote as “HLinear5 MPP”, to test its good performance. The Hermite linear scheme without bound preserving MPP flux limiter is denoted as “HLinear5”, while those with Hermite WENO reconstruction (2.5) are denoted as “HWENO5 MPP” and “HWENO5” for with and without MPP flux limiter correspondingly. All the schemes without Hermite reconstruction, for example the one in [32], will be denoted as “WENO5 MPP” and “WENO5”, and the corresponding linear schemes as “Linear5 MPP” and “Linear5”, respectively. They will also be used for comparison in the following. The time step is taken as

[TABLE]

where $\alpha_{x}=\|U_{x}\|_{L^{\infty}}$, $\alpha_{y}=\|U_{y}\|_{L^{\infty}}$ and we take the CFL number to be $\text{CFL}=0.6$ in all the following tests.

We would emphasize that to apply the MPP limiter, the minimum and maximum values of the bounds are taken from the initial data, which can be explicitly determined for all the tests here.

3.1. Linear transport equations

Example 3.1**.**

We first take a 2D linear transport equation $\rho_{t}+\rho_{x}+\rho_{y}=0$ with initial condition

[TABLE]

on the domain $[0,\,2\pi]^{2}$ with periodic boundary conditions, to test the accuracy of the “HLinear5 MPP” and “HLinear5” schemes. The exact solution is

[TABLE]

The solution is within $[0,\,2]$. From Table 3.1, we can clearly observe undershootings without the MPP limiter, while with limiter the lower bound is well preserved (as is for the upper bound) and 5th order accuracy is maintained. We would mention that for all other schemes we considered in this section, the accuracies are all tested and perform similarly. We omit them to save space.

Then we take a 1D linear transport equation $\rho_{t}+\rho_{x}=0$ with highly oscillatory solutions. Here the initial condition is chosen to be

[TABLE]

on the domain $[0,\,2\pi]$ with periodic boundary condition, and the exact solution is

[TABLE]

In Fig. 3.1, we show the results for both linear and WENO type schemes, without and with MPP limiter respectively. From the figures, we can observe that the “HLinear5” scheme has better resolution than the “Linear5” scheme, due to a Hermite reconstruction for such a highly oscillatory solutions, while linear schemes are better than WENO type schemes due to less dissipation. We can also see that the resolution of these schemes do not significantly affect by the MPP limiter, which is used to control the solutions within the global physical range $[-1,\,1]$. So in the following, we will focus on the “HLinear5” and “HLinear5 MPP” schemes, and compare them to their corresponding WENO type approaches, which are “HWENO5” and “HWENO5 MPP” respectively. We mainly would like to show that a 5th order linear scheme with a Hermite reconstruction and an MPP limiter, that is “HLinear5 MPP”, is a very good scheme for highly oscillatory solution without discontinuities.

3.2. Vlasov-Poisson system

In this subsection, we consider the VP system (1.2) and (1.3) with 1D in ${\bf x}$ and 1D in ${\bf v}$, and periodic boundary condition in both directions. We take ${\bf x}=(x,v)$ and use $N_{v}$ instead of $N_{y}$ for the special meaning of the argument $v$ here. The cut-off domain in the $v$ direction is taken to be $[-2\pi,2\pi]$ if without specifications. The VP system preserves several norms, which should remain constant in time:

1. (1)

$L^{p}$ norm, $1\leq p<\infty$:

[TABLE] 2. (2)

Energy:

[TABLE]

where $E(x,t)$ is the electric field. 3. (3)

Entropy:

[TABLE]

We will track the relative deviations of these quantities numerically to measure the quality of our numerical scheme. The relative deviation is defined to be the deviation away from the corresponding initial value divided by the magnitude of the initial value. The mass conservation over time $\int_{v}\int_{x}f(x,v,t)dxdv$ is obvious for conservative schemes, which is the same as the $L^{1}$ norm when $f$ is positive, so only the $L^{1}$ norm is shown below. We would also note that for the VP system, although the minimum part is known as positivity preserving, however, due to a cut-off domain in the $v$ direction, the minimum value might be close to but above [math]. We clearly indicate the minimum and maximum values from the initial data for the tests below.

Example 3.2**.**

(Accuracy test) We first consider the VP system with the following initial condition

[TABLE]

and periodic boundary conditions on the computational domain $[0,4\pi]\times[-4\pi,4\pi]$, where $k=0.5$, to test the accuracy of the schemes for this system. From the initial data, the solution should be within the range $[0,\,\frac{1}{\sqrt{2\pi}}]$.

In Table 3.2, we show the $L^{1}$ and $L^{\infty}$ errors and orders for the “HLinear5” and “HLinear5 MPP” scheme respectively. For this example, $5$th order accuracy can also be observed as mesh refinement. Without limiter, the solution of the distribution function does not preserve positivity, while with limiter, all values are above [math]. Here we measure the errors on two consecutive mesh sizes by double refinement, since the exact solution is not explicitly available.

Example 3.3**.**

(Strong Landau damping) We then consider the strong Landau damping with the initial condition:

[TABLE]

where $\alpha=0.5$ and $k=0.5$. The length of the domain in the x-direction is $L=\frac{2\pi}{k}$, which is similar in the following two examples. For this problem, from the initial data, the solution should be within the range $[\frac{1}{\sqrt{2\pi}}(1-\alpha)\exp(-\frac{(2\pi)^{2}}{2}),\,\frac{1}{\sqrt{2\pi}}(1+\alpha)]$. In Fig. 3.2, we plot the surface of the distribution function $f$ at $T=50$ in the range of $[0,0.6]$. The mesh grid is $256\times 256$. We can observe that without limiter, the solution can be negative (white spots), while negative values are eliminated by the scheme with limiter, similarly for the following examples. Then in Fig. 3.3, we show the time evolution of the electric field in $L^{2}$ norm and $L^{\infty}$ norm, the relative derivation of the discrete $L^{1}$ norm, $L^{2}$ norm, kinetic energy and entropy. From this figure, we can see that the electric field for all schemes are almost the same. The linear type scheme can preserve the $L^{2}$ norm and entropy better than the WENO type scheme, where the MPP limiter does not significantly affect these two quantities. However, for the $L^{1}$ norm and energy, the linear scheme without limiter performs much worse than the WENO type scheme, but after with limiter, it is the best. Especially for the energy, “HLinear5 MPP” is much better than all other schemes. This example has clearly show the good performance of our approach.

Example 3.4**.**

(Symmetric two stream instability) We now consider the symmetric two stream instability with the initial condition:

[TABLE]

where $\alpha=0.05$, $u=0.99$, $v_{th}=0.3$ and $k=\frac{2}{13}$. Similarly, from the initial data, the minimum value of the solution is taking $v=2\pi$ and $\cos(kx)=-1$ in (3.10), while the maximum value is taking $v=u$ or $v=-u$, and $\cos(kx)=1$. We plot the numerical solution at $T=70$ in Fig. 3.4 for the “HLinear5” and “HLinear5 MPP” schemes respectively. The mesh grid is $256\times 256$. We can also clearly observe that without limiter, the solution becomes negative, which, however, can be eliminated by the scheme with limiter. We show the time evolution of the electric field in $L^{2}$ norm and $L^{\infty}$ norm, the relative derivation of the discrete $L^{1}$ norm, $L^{2}$ norm, kinetic energy and entropy in Fig. 3.5. Similar results as the last example are obtained.

Example 3.5**.**

(Bump-on-tail instability) We consider the bump-on-tail instability problem [1, 33] with the initial condition as

[TABLE]

where the bump-on-tail distribution is

[TABLE]

and $n_{p}=0.9$, $n_{b}=0.2$, $v_{b}=4.5$, $v_{t}=0.5$, $\alpha=0.04$, $k=0.3$. The computational domain is $[0,\frac{2\pi}{k}]\times[-3\pi,3\pi]$. We take the cut-off domain in $v$ a little larger for the consideration of the shifting $v_{b}$ in the velocity. For this problem, the minimum value of the initial data is taking $v=-3\pi$ and $\cos(kx)=-1$ in (3.11) and (3.12), while the maximum value is taking $v=0$ and $\cos(kx)=1$. We run a long time simulation up to $T=1000$, which is a good test to show the effectiveness of the linear scheme in preserving some important quantities and the saving of computational cost as compared to the WENO type scheme. We first show the surface of the distribution function at $T=500$ in Fig. 3.6, we can clearly observe without limiter, the solution has undershootings, while the scheme with limiter produces very good result. We show the time evolution of the electric field in $L^{2}$ norm and $L^{\infty}$ norm in Fig. 3.7, and compare the linear scheme with the WENO type scheme. The $L^{\infty}$ norm of the electric filed $E_{max}$ matches those in [33, 1]. However, we can see that the linear scheme preserves the electric field better than the WENO type scheme after some large time, e.g., $T=400$, while the results are almost the same between with and without limiter. We then show the relative derivation of the discrete $L^{1}$ norm, $L^{2}$ norm, kinetic energy and entropy in Fig. 3.8. The MPP limiter almost does not affect the $L^{2}$ norm and entropy, however, again we can see it improves the $L^{1}$ norm and the energy a lot from eliminating negative numerical values, which indicates the great performance of the MPP limiter on the linear scheme. Besides, the linear scheme preserves the $L^{2}$ norm, entropy and energy much better than the WENO type scheme, especially the energy. Moreover, for the computational cost, the linear scheme is at least about 3-4 times faster than the WENO type scheme. This example with long time simulation clearly demonstrates that the linear scheme is better than the WENO type scheme, which is much less dissipative and computationally much more efficient.

3.3. Kelvin-Helmholtz instability

Example 3.6**.**

The Kelvin-Helmholtz instability comes from the 2D guiding-center model (1.4) [16] with the initial data

[TABLE]

and periodic boundary conditions on the domain $[0,4\pi]\times[0,2\pi]$. We let $k=0.5$, which will create a Kelvin-Helmholtz instability. The exact solution should be within the range $[-1.015,1.015]$.

For this example, we show the solution at $T=40$ with mesh grid $256\times 256$ for “HLinear5” and “HLinear5 MPP” in Fig. 3.9 respectively. For the figures drawing in the physical range $[-1.015,1.015]$, we can observe white spots for “HLinear5” which is without MPP limiter, while “HLinear5 MPP” with MPP limiter preserves the bounds very well. In Fig. 3.10, we compare the time evolution of the relative $L^{2}$ norm and the numerical minimum values for both linear and WENO type schemes without and with limiters respectively. First, we can find that the linear scheme is less dissipative than the WENO one, as the $L^{2}$ norm preserves better for the linear scheme. The schemes with and without MPP limiter preserves the $L^{2}$ norm almost the same, so that there is no significant affection from the MPP limiter. For the minimum numerical values, the linear scheme without MPP limiter has shown large undershootings, the WENO type scheme without MPP limiter performs even worse for this example (which is also similar for the classic WENO scheme “WENO5”, we omit the figure here), but both schemes with MPP limiter preserve the lower bound very well. Similarly for the upper bound.

From this example, we can see that for problems with highly oscillatory but without discontinuous solutions, a high order linear scheme with the MPP limiter can result in a very good scheme, which is less dissipative and without significant spurious oscillations. The original idea to apply WENO reconstruction to suppress numerical oscillations seems to be even fail for this example. We can take this example as a benchmark test to support our main idea in this paper.

3.4. Incompressible Euler equations

Example 3.7**.**

(Accuracy test) We first consider the incompressible Euler system on the domain $[0,2\pi]\times[0,2\pi]$ with an initial condition $\omega_{0}(x,y)=-2\sin(x)\sin(y)$. The exact solution will stay stationary with $\omega(x,y,t)=-2\sin(x)\sin(y)$, which is in the range of $[-2,2]$. For this example, from Table 3.3, we can see that the numerical solution with and without the limiter all stay in the right range, which shows that the limiter does not destroy the high order accuracy.

Example 3.8**.**

(Vortex patch). In this example, we consider the incompressible Euler equations for the vortex patch problem with the initial condition given by

[TABLE]

We show the surface of $\omega$ at $T=10$ in Fig. 3.11. The mesh grid is $256\times 256$. We can observe the good performance of the MPP flux limiter on this problem. We also show the time evolution of the relative difference of $L^{2}$ norm as compared to the initial data, and the minimum numerical solution on the bottom of Fig. 3.11. Here we can still see that the linear scheme preserves the $L^{2}$ norm better than the WENO type scheme and the MPP limiter can eliminate the oscillations around two extreme values $1$ and $-1$ from the schemes without limiter. However, we would mention that for this example with discontinuous initial data, the MPP limiter cannot remove the small numerical oscillations around [math] (figures are omitted here). In this case, the WENO will be needed to completely remove oscillations for middle discontinuities.

4. Conclusion

In this paper, we have proposed a linear scheme combined with a high order bound-preserving MPP flux limiter, for solving incompressible flow problems. As compared to WENO type schemes, our approach is less dissipative and much less costly, so that is much more efficient for high dimensional problems with long time simulations. Applications to the Vlasov-Poisson system, 2D guiding-center model in plasma physics, as well as the incompressible Euler equations in fluid hydrodynamics have demonstrated the good performance of our proposed approach.

Acknowledgement

T. Xiong would like to thank Dr. F. Filbet for many helpful discussions. T. Xiong acknowledges support from the Marie Skłodowska-Curie Individual Fellowships H2020-MSCA-IF-2014 of the European Commission, under the project HNSKMAP 654175. This work has also been supported by the Fundamental Research Funds for the Central Universities No. 20720160009, NSFC grant 11601455, NSAF grant U1630247 and NSF grant of Fujian Province 2016J05022.