Projective integration for nonlinear BGK kinetic equations

TL;DR

This paper introduces a high-order, explicit projective integration scheme for nonlinear BGK kinetic equations that efficiently handles stiffness by combining small inner steps with an outer Runge-Kutta method, ensuring stability and accuracy.

Contribution

The paper develops an asymptotic-preserving, fully explicit projective integration method for nonlinear BGK equations, allowing larger time steps independent of stiffness.

Findings

Method achieves stability with outer time steps independent of stiffness

Numerical results confirm high-order accuracy in 1D and 2D

Inner steps effectively damp stiff components

Abstract

We present a high-order, fully explicit, asymptotic-preserving projective integration scheme for the nonlinear BGK equation. The method first takes a few small (inner) steps with a simple, explicit method (such as direct forward Euler) to damp out the stiff components of the solution. Then, the time derivative is estimated and used in an (outer) Runge-Kutta method of arbitrary order. Based on the spectrum of the linearized BGK operator, we deduce that, with an appropriate choice of inner step size, the time step restriction on the outer time step as well as the number of inner time steps is independent of the stiffness of the BGK source term. We illustrate the method with numerical results in one and two spatial dimensions.

Equations

$$\partial_{t}f^{\varepsilon}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f^{\varepsilon}=\frac{1}{\varepsilon}\mathcal{Q}\mathopen{}\mathclose{{}\left(f^{\varepsilon}}\right)(\mathbf{v}),$$

$$\partial_{t}f^{\varepsilon}+\mathbf{v}\cdot\nabla_{\mathbf{x}}f^{\varepsilon}=\frac{1}{\varepsilon}(\mathcal{M}_{\mathbf{v}}\mathopen{}\mathclose{{}\left(f^{\varepsilon}}\right)-f^{\varepsilon}),$$

$$\mathcal{M}_{\mathbf{v}}\mathopen{}\mathclose{{}\left(f^{\varepsilon}}\right)=\frac{\rho}{(2\pi T)^{D_{v}/2}}\exp{\mathopen{}\mathclose{{}\left(-\frac{|\mathbf{v}-\mathbf{\bar{v}}|^{2}}{2T}}\right)}:=\mathcal{M}_{\mathbf{v}}^{\rho,\mathbf{\bar{v}},T}.$$

$$\rho=\int_{\mathbb{R}^{D_{v}}}f^{\varepsilon}d\mathbf{v},\qquad\bar{v}^{d}=\frac{1}{\rho}\int_{\mathbb{R}^{D_{v}}}v^{d}f^{\varepsilon}d\mathbf{v},\qquad T=\frac{1}{D_{v}\rho}\int_{\mathbb{R}^{D_{v}}}\mathopen{}\mathclose{{}\left|\mathbf{v}-\mathbf{\bar{v}}}\right|^{2}f^{\varepsilon}d\mathbf{v},$$

$$\mathopen{}\mathclose{{}\left\{\begin{aligned} &\partial_{t}\rho+\operatorname{div}_{\mathbf{x}}(\rho\,\mathbf{\bar{v}})\,=\,0,\\ &\partial_{t}(\rho\,\mathbf{\bar{v}})+\operatorname{div}_{\mathbf{x}}\mathopen{}\mathclose{{}\left(\rho\,\mathbf{\bar{v}}\otimes\mathbf{\bar{v}}\,+\,\rho\,T\,{\rm\bf I}}\right)\,=\,\bm{0},\\ &\partial_{t}E+\operatorname{div}_{\mathbf{x}}\mathopen{}\mathclose{{}\left(\mathbf{\bar{v}}\mathopen{}\mathclose{{}\left(E+\rho\,T}\right)}\right)\,=\,0,\end{aligned}}\right.$$

$$\dot{\mathbf{f}}=\mathrm{D}_{t}(\mathbf{f}),\qquad\mathrm{D}_{t}(\mathbf{f})=-\mathrm{D}_{\bm{x},\bm{v}}(\mathbf{f})+\frac{1}{\varepsilon}\mathopen{}\mathclose{{}\left(\mathcal{M}_{\bm{v}}\mathopen{}\mathclose{{}\left(\mathbf{f}}\right)-\mathbf{f}}\right),$$

$${\mathbf{f}}^{k+1}=S_{\delta t}({\mathbf{f}}^k)={\mathbf{f}}^k+\delta t{\mathrm{D}}_t({\mathbf{f}}^k),k=0,1,\dots .$$

$${\mathbf{f}}^{n+1}={\mathbf{f}}^{n,K+1}+(\Delta t-(K+1)\delta t)\frac{{\mathbf{f}}^{n,K+1}-{\mathbf{f}}^{n,K}}{\delta t}.$$

$$s=1:$$

$$2\le s\le S:$$

$${\mathbf{f}}^{n+1}={\mathbf{f}}^{n,K+1}+(\Delta t-(K+1)\delta t)\sum _{s=1}^S{b}_s{\mathbf{k}}_s.$$

$$\sum _{s=1}^S{b}_s=1,\sum _{i=1}^{S-1}{a}_{s,i}=c_s,1\le s\le S.$$

$${\mathcal{M}}_{\text{lin}}(f^{\epsilon })(\mathbf{x},\mathbf{v},t)=\sum _{k=0}^{D_v+1}{\Psi }_k(\mathbf{v})({\Psi }_k,f^{\epsilon })(\mathbf{x},t),$$

$$(g,h)=\int_{\mathbb{R}^{D_{v}}}g(\mathbf{v})\overline{h(\mathbf{v})}\frac{1}{(2\pi)^{D_{v}/2}}\exp\mathopen{}\mathclose{{}\left(\frac{-\mathopen{}\mathclose{{}\left|\mathbf{v}}\right|^{2}}{2}}\right)d\mathbf{v}.$$

$$\big{(}\Psi_{0}(\mathbf{v}),\ldots,\Psi_{D_{v}+1}(\mathbf{v})\big{)}=\mathopen{}\mathclose{{}\left(1,v^{1},...,v^{D_{v}},\frac{\mathopen{}\mathclose{{}\left|\mathbf{v}}\right|^{2}-D_{v}}{2^{D_{v}/2}}}\right).$$

$${\partial }_t{f}^{\epsilon }+\mathbf{v}\cdot {\nabla }_{\mathbf{x}}{f}^{\epsilon }=-\frac{1}{\epsilon }(\mathcal{I}-{\Pi }_{\text{BGK}})f^{\epsilon },$$

$${\Pi }_{\text{BGK}}{f}^{\epsilon }=\sum _{k=0}^{D_v+1}{\Psi }_k(\mathbf{v})({\Psi }_k,f^{\epsilon }).$$

$$\big{(}\rho_{L},\bar{v}_{L},T_{L}\big{)}=(1,0,1),\qquad\quad\big{(}\rho_{R},\bar{v}_{R},T_{R}\big{)}=(0.125,0,0.25).$$

$$q^{d}=\frac{1}{2}\int_{\mathbb{R}^{D_{v}}}\mathopen{}\mathclose{{}\left|\mathbf{c}}\right|^{2}c^{d}f^{\varepsilon}d\mathbf{v},$$

$$\big{(}\rho_{L},\bar{v}^{x}_{L},\bar{v}^{y}_{L},T_{L}\big{)}=\mathopen{}\mathclose{{}\left(\frac{16}{7},\sqrt{\frac{5}{3}}\frac{7}{16},0,\frac{133}{64}}\right),\qquad\quad\big{(}\rho_{R},\mathbf{\bar{v}}_{R},T_{R}\big{)}=\mathopen{}\mathclose{{}\left(1,\bm{0},1}\right).$$

$$\rho(\mathbf{x},0)=1+1.5\exp\mathopen{}\mathclose{{}\left(-16\mathopen{}\mathclose{{}\left|\mathbf{x}-\mathbf{x}_{0}}\right|^{2}}\right).$$

Full text

MSC (2010): 82B40, 76P05, 65M08, 65L06.

0.1 Introduction

The Boltzmann equation constitutes the cornerstone of the kinetic theory of rarefied gases. In a dimensionless, scalar setting, it describes the evolution of the one-particle mass distribution function $f^{\varepsilon}(\mathbf{x},\mathbf{v},t)\in\mathbb{R}^{+}$ as:

[TABLE]

where $t\geq 0$ represents time, and $(\mathbf{x},\mathbf{v})\subset\mathbb{R}^{D_{x}\times D_{v}}$ are the $D_{x}$-dimensional particle positions and $D_{v}$-dimensional particle velocities. In equation (1), the dimensionless constant $\varepsilon>0$ determines the regime of the gas flow, for which we roughly identify the hydrodynamic regime $(\varepsilon\leq 10^{-4})$, the transitional regime $(\varepsilon\in[10^{-4},10^{-1}])$, and the kinetic regime $(\varepsilon\geq 10^{-1})$. Furthermore, the left hand side of (1) corresponds to a linear transport operator that comprises the convection of particles in space, whereas the right hand side contains the Boltzmann collision operator that entails velocity changes due to particle collisions. However, due to its high-dimensional and complicated structure, the Boltzmann collision operator is often replaced by simpler collision models that capture most essential features of the former. The most well-known such model is the BGK model Bhatnagar1954 , which models collisions as a linear relaxation towards thermodynamic equilibrium, and is given by:

[TABLE]

in which \mathcal{M}_{\mathbf{v}}\mathopen{}\mathclose{{}\left(f^{\varepsilon}}\right) denotes the local Maxwellian distribution, which, for a $D_{v}$-dimensional velocity space, is given by:

[TABLE]

The Maxwellian distribution contains the velocity moments of the distribution function $f^{\varepsilon}$, which are calculated as:

[TABLE]

where $\rho\in\mathbb{R}^{+}$, \mathbf{\bar{v}}=\mathopen{}\mathclose{{}\left(\bar{v}^{d}}\right)_{d=1}^{D_{v}}\in\mathbb{R}^{D_{v}} and $T\in\mathbb{R}^{+}$ are the density, macroscopic velocity and temperature, respectively, which all depend on space $\mathbf{x}$ and time $t$. Then, in the limit $\varepsilon\to 0$, the solution to equation (2) converges towards $\mathcal{M}_{\mathbf{v}}^{\rho,\mathbf{\bar{v}},T}$, whose moments in (4) are solution to the compressible Euler system:

[TABLE]

in which $E$ is the second moment of $f^{\varepsilon}$, namely its total energy.

In this paper, we construct a fully explicit, asymptotic-preserving, arbitrary order time integration method for the stiff equation (2). For a comprehensive review of numerical schemes for collisional kinetic equations such as equation (1), we refer to DimarcoPareschi15 . The asymptotic-preserving property Jin1999 implies that, in the limit when $\varepsilon$ tends to zero, an $\varepsilon$-independent time step constraint, of the form $\Delta t=O(\Delta x)$, can be used, in agreement with the classical hyperbolic CFL constraint for the limiting fluid equations (5). To achieve this, we will use a projective integration method, which was introduced in Gear2003projective and first applied to kinetic equations in lafitte2012 .

The remainder of this paper is structured as follows. We describe the projective integration method in more detail in section 0.2, after which we discuss (in section 0.3) the spectral properties of the linearized BGK operator, which are needed to ensure stability of the method. Some numerical experiments are done in section 0.4.

0.2 Projective integration

Projective integration Gear2003projective ; lafitte2012 combines a few small time steps with a naive (inner) timestepping method (here, a direct forward Euler discretization) with a much larger (projective, outer) time step. The idea is sketched in figure 1.

Inner integrators. We discretize equation (2) on a uniform, constant in time, periodic spatial mesh with spacing $\Delta x$, consisting of $I$ mesh points $x_{i}=i\Delta x$, ${1\leq i\leq I}$, with $I\Delta x=1$, and a uniform time mesh with time step $\delta t$ and discrete time instants $t^{k}=k\delta t$. Furthermore, we discretize velocity space by choosing $J$ discrete components denoted by $\mathbf{v}_{j}$. The numerical solution on this mesh is denoted by $f_{i,j}^{k}$, where we have dropped the superscript $\varepsilon$ on discretized quantities. We then obtain a semidiscrete system of ODEs of the form:

[TABLE]

where $\mathrm{D}_{\bm{x},\bm{v}}(\cdot)$ represents a suitable discretization of the convective derivative $\mathbf{v}\cdot\nabla_{\mathbf{x}}$ (for instance, using upwind differences), and $\mathbf{f}$ is a vector of size $I\cdot J$.

As inner integrator, we choose the (explicit) forward Euler method with time step $\delta t$, for which we will, later on, use the shorthand notation:

[TABLE]

Outer integrators. In system (6), the small parameter $\varepsilon$ leads to the classical time step restriction of the form $\delta t=O(\varepsilon)$ for the inner integrator. However, as $\varepsilon$ goes to [math], we obtain the limiting system (5) for which a standard finite volume/forward Euler method only needs to satisfy a stability restriction of the form $\Delta t\leq C\Delta x$, with $C$ a constant that depends on the specific choice of the scheme.

In lafitte2012 , it was proposed to use a projective integration method to accelerate such a brute-force integration; the idea, originating from Gear2003projective , is the following. Starting from a computed numerical solution $\mathbf{f}^{n}$ at time $t^{n}=n\Delta t$, one first takes $K+1$ inner steps of size $\delta t$ using (7), denoted as $\mathbf{f}^{n,k+1}$, in which the superscripts $(n,k)$ denote the numerical solution at time ${t^{n,k}=n\Delta t+k\delta t}$. The aim is to obtain a discrete derivative to be used in the outer step to compute $\mathbf{f}^{n+1}=\mathbf{f}^{n+1,0}$ via extrapolation in time:

[TABLE]

Higher-order projective Runge-Kutta (PRK) methods can be constructed by replacing each time derivative evaluation $\mathbf{k}_{s}$ in a classical Runge-Kutta method by $K+1$ steps of an inner integrator as follows:

[TABLE]

To ensure consistency, the Runge-Kutta matrix $\mathbf{a}=(a_{s,i})_{s,i=1}^{S}$, weights ${\mathbf{b}=(b_{s})_{s=1}^{S}}$, and nodes $\mathbf{c}=(c_{s})_{s=1}^{S}$ satisfy the conditions $0\leq b_{s}\leq 1$ and $0\leq c_{s}\leq 1,$ as well as:

[TABLE]

0.3 Spectral properties

To choose the method parameters (the size of the small and large time steps $\delta t$ and $\Delta t$, as well as the number $K$ of small steps), one needs to analyze the spectrum of the collision operator. In Melis2016 , this was done in the hyperbolic scaling for a system with a linear Maxwellian that serves as a relaxation of a nonlinear hyperbolic conservation law.

By linearizing the Maxwellian (3) around the global Maxwellian distribution $\mathcal{M}^{\rho^{\infty},\mathbf{\bar{v}}^{\infty},T^{\infty}}_{\mathbf{v}}=\mathcal{M}_{\mathbf{v}}^{1,0,1}$, it is shown in (Cercignani1988, , p.206) that the resulting linearized equilibrium can be written as:

[TABLE]

in which the scalar product is defined by:

[TABLE]

Furthermore, the orthonormal set of basis functions $\Psi_{k}(\mathbf{v})$ in (13) are obtained from a straightforward application of the Gram-Schmidt process to the $D_{v}+1$ collision invariants (1,\mathbf{v},\mathopen{}\mathclose{{}\left|\mathbf{v}}\right|^{2}), yielding:

[TABLE]

Using the linearized Maxwellian (13), the linearized version of the full BGK equation (2) reads:

[TABLE]

where $\mathcal{I}$ denotes the identity operator and $\Pi_{\text{BGK}}$ is the following rank-$(D_{v}+2)$ projection operator:

[TABLE]

This shows that the structure of the linearized Maxwellian (13) and the linearized BGK projection operator (17) are almost identical to those in Melis2016 . We can actually view these linear kinetic models as a special simplified case of the linearized BGK equation. Therefore, it is expected that the construction of stable, asymptotic-preserving projective integration methods for the full BGK equation (2) is practically identical to that in Melis2016 . In particular, the conclusion is that, when choosing $\delta t=\varepsilon$, one is able to choose $\Delta t=O(\Delta x)$ and $K$ independent of $\varepsilon$, resulting in a scheme with computational cost independent of $\varepsilon$.

0.4 Numerical experiments

BGK in 1D. As a first experiment, we focus on the nonlinear BGK equation (2) in 1D. We consider a Sod-like test case for $x\in[0,1]$ consisting of an initial centered Riemann problem with the following left and right state values:

[TABLE]

The initial distribution $f^{\varepsilon}(x,v,0)$ is then chosen as the Maxwellian (3) corresponding to the above initial macroscopic variables. We impose outflow boundary conditions and perform simulations for $t\in[0,0.15]$. As velocity space, we take the interval $[-8,8]$, which we discretize on a uniform grid using $J=80$ velocity nodes. In all simulations, space is discretized using the WENO3 spatial discretization with $\Delta x=0.01$. Below, we compare solutions for three gas flow regimes: $\varepsilon=10^{-1}$ (kinetic regime), $\varepsilon=10^{-2}$ (transitional regime) and $\varepsilon=10^{-5}$ (fluid regime).

In the kinetic $(\varepsilon=10^{-1})$ and transitional $(\varepsilon=10^{-2})$ regimes, we compute the numerical solution using the fourth order Runge-Kutta (RK4) time discretization with time step $\delta t=0.1\Delta x$. In the fluid regime $(\varepsilon=10^{-5})$, direct integration schemes such as RK4 become too expensive due to a severe time step restriction, which is required to ensure stability of the method. Exploiting that the spectrum of the linearized BGK equation is close to that of the linear kinetic models used in Melis2016 , see section 0.3, we construct a projective integration method to accelerate time integration in the fluid regime. As inner integrator, we select the forward Euler time discretization with $\delta t=\varepsilon$. As outer integrator, we choose the fourth-order projective Runge-Kutta (PRK4) method, using $K=2$ inner steps and an outer step of size $\Delta t=0.4\Delta x$.

The results are shown in figure 2, where we display the density $\rho$, macroscopic velocity $\bar{v}$ and temperature $T$ as given in (4) at $t=0.15$. In addition, we plot the heat flux $q$, which, in a general $D_{v}$-dimensional setting, is a vector \mathbf{q}=\mathopen{}\mathclose{{}\left(q^{d}}\right)_{d=1}^{D_{v}} with components given by:

[TABLE]

in which \mathbf{c}=\mathopen{}\mathclose{{}\left(c^{d}}\right)_{d=1}^{D_{v}}=\mathbf{v}-\mathbf{\bar{v}} is the peculiar velocity. The different regimes are shown by blue (kinetic), purple (transitional) and green (fluid) dots. The red line in each plot denotes the limiting $(\varepsilon\to 0)$ solution of each macroscopic variable, which all converge to the solution of the compressible Euler equations (5) with ideal gas law $P=\rho T$ and heat flux $q=0$. From this, we observe that the BGK solution is increasingly dissipative for increasing values of $\varepsilon$ since the rate with which $f^{\varepsilon}$ converges to its equilibrium \mathcal{M}_{v}\mathopen{}\mathclose{{}\left(f^{\varepsilon}}\right) becomes slower. In contrast, for sufficiently small $\varepsilon$, relaxation to thermodynamic equilibrium occurs practically instantaneous and the Euler equations (5) yield a valid description. Since this is a hyperbolic system, it allows for the development of sharp discontinuous and shock waves which are clearly seen in the numerical solution.

Shock-bubble interaction in 2D. Here, we consider the BGK equation in 2D and we investigate the interaction between a moving shock wave and a stationary smooth bubble, which was proposed in Torrilhon2006 , see also Cai2010 . This problem consists of a shock wave positioned at $x=-1$ in a spatial domain $\mathbf{x}=(x,y)\in[-2,3]\times[-1,1]$ traveling with Mach number $\mathit{Ma}=2$ into an equilibrium flow region. Over the shock wave, the following left $(x\leq-1)$ and right $(x>-1)$ state values are imposed Cai2010 :

[TABLE]

Due to this initial profile, the shock wave will propagate rightwards into the flow region at rest $(x>-1)$. Moreover, in this equilibrium region, a smooth Gaussian density bubble centered at $\mathbf{x}_{0}=(0.5,0)$ is placed, given by:

[TABLE]

Then, the initial distribution $f^{\varepsilon}(\mathbf{x},\mathbf{v},0)$ is chosen as the Maxwellian (3) corresponding to the initial macroscopic variables in (20)-(21). We impose outflow and periodic boundary conditions along the $x$- and $y$-directions, respectively, and we perform simulations for $t\in[0,0.8]$. As velocity space, we take the domain $[-10,10]^{2}$, which we discretize on a uniform grid using $J_{x}=J_{y}=30$. We discretize space using the WENO2 spatial discretization with $I_{x}=200$ and $I_{y}=25$. Furthermore, we consider a fluid regime by taking $\varepsilon=10^{-5}$.

We construct a PRK4 method with FE as inner integrator to speed up simulation in time. The inner time step is fixed as $\delta t=\varepsilon$ and we use $K=2$ inner steps in each outer integrator iteration. The outer time step is chosen as $\Delta t=0.4\Delta x$. To compare our results with those in Torrilhon2006 , where the smallest value of $\varepsilon$ is chosen as $\varepsilon=10^{-2}$, we regard the one-dimensional evolution of density and temperature along the axis $y=0$. For ${t\in\{0,0.2,0.4,0.6,0.8\}}$, we plot these intersections in figure 3. We conclude that we obtain the same solution structure at $t=0.8$ as in Torrilhon2006 . However, our results are sharper and less dissipative supposedly due to the particular small value of $\varepsilon$ ($10^{-5}$ versus $10^{-2}$). In contrast to Cai2010 , we nicely capture the swift changes in the temperature profile for $x\in[0.5,1]$ at $t=0.8$.