Journal of Scientific Computing https://doi.org/10.1007/s10915-018-0776-9

A Class of Low Dissipative Schemes for Solving Kinetic Equations Giacomo Dimarco1

· Cory Hauck2 · Raphaël Loubère3

Received: 18 May 2017 / Revised: 6 May 2018 / Accepted: 22 June 2018 © Springer Science+Business Media, LLC, part of Springer Nature 2018

Abstract We introduce an extension of the fast semi-Lagrangian scheme developed in J Comput Phys 255:680–698 (2013) in an effort to increase the spatial accuracy of the method. The basic idea of this extension is to modify the first-order accurate transport step of the original semi-Lagrangian scheme to allow for a general piecewise polynomial reconstruction of the distribution function. For each discrete velocity, we update the solution not at cell centers, but rather at the extreme points of the spatial reconstruction, the locations of which are different for each discrete velocity and change with time. Several approaches are discussed for evaluating the collision operator at these extreme points using only cell center values by making special assumption on the spatial variation of the collision operator. The result is a class of schemes that preserves the structure of the solution over very long times when compared to existing schemes in the literature. As a proof of concept, the new method is implemented in a one-dimensional setting, using piecewise linear reconstructions of the distribution function together with a related reconstruction of the collision operator. The method is derived both for the relatively simple Bhatnagar–Gross–Krook (BGK) operator as well as for the classical Boltzmann operator. Several numerical tests are used to assess the performance of the implementation, including comparisons with the original method in J Comput Phys 255:680–698 (2013) and with classical semi-Lagrangian methods of first and second order. In convergence tests, we observe uniform second-order accuracy across all flow regimes for the BGK operator and nearly second-order accuracy for the Boltzmann operator. In addition, we observe that the method outperforms the classical semi-Lagrangian approach, in particular when resolving fine solution structures in space.

This work was partly supported by Campus France under the PHC Gallileo 2015-2016 Number 32272UL. This manuscript has been authored, in part, by UT-Battelle, LLC, under Contract No. DE-AC0500OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for the United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan). The research of C.H. was sponsored by the Office of Advanced Scientific Computing Research and performed at the Oak Ridge National Laboratory, which is managed by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725. Extended author information available on the last page of the article

123

Journal of Scientific Computing

Keywords Kinetic equations · Boltzmann equation · Semi-Lagrangian schemes · High-order schemes · Rarefied gas dynamics Mathematics Subject Classification 65M25 · 82C40 · 76P05

1 Introduction Kinetic modeling is an important tool that is used to describe many different phenomena, ranging from classical fluid mechanics, plasmas and rarefied gases [2,6,21] to biological and socio-economic models [43]. The advantage of kinetic equations over classical fluid mechanical equations is the ability to take into account a broader range of physical regimes and phenomena. In particular, fluid models rely on the assumption of a local thermodynamic equilibrium state while kinetic models do not [7, Chapter 3]. Unfortunately, the additional fidelity provided by kinetic equations is accompanied by a variety of numerical challenges [21]. Typically, one would like to preserve conserved quantities at the numerical level (because they characterize states of local thermal equilibrium) and to handle the wide range of spatial and temporal scales that are encountered in practical applications. An even more fundamental problem is the size of the phase space: in a fluid description, the equations are defined over (x, t) ∈ R3 × R+ ; a kinetic equation, on the other hand, is defined over (x, v, t) ∈ R3 × R3 × R+ . This increase in dimension makes kinetic simulations of realistic, multi-dimensional problems very demanding, both in terms of computing time and memory resources. There are generally two approaches to simulating kinetic equations: deterministic methods such as finite volume, semi-Lagrangian, and spectral schemes [21,49] and probabilistic methods such as Direct Simulation Monte Carlo (DSMC) [1,5]. Both methodologies have strengths and weaknesses: while the first can obtain a high order of accuracy, the second is usually faster, especially for solving steady-state problems. However, probabilistic methods typically exhibit lower convergence rates and difficulties in describing non-stationary and slow-motion flows. Recent research efforts have focused on reducing the disadvantages of both deterministic and probabilistic approaches. A comprehensive review of current numerical techniques for deterministic discretization of kinetic equations can be found in [21]. Concerning the cost of evaluating collision operators, we recall the fast spectral approaches developed in [23,27, 42,44]. For tackling the multiple scales problem, we note the work in [19,20,33,36] and the references therein for a broad view of the numerical challenges and the techniques used to address them. For issues related to the excessive noise of Monte Carlo methods, we quote [5] for an overview on efficient and low variance Monte Carlo methods; and for applications of variance reduction techniques to kinetic equations, we recall the recent work in [34,35] and [12,15,18]. In this work, we focus on the construction of efficient numerical methods for the linear transport operator in a kinetic model of neutral particles. More specifically, we extend the semi-Lagrangian, Fast Kinetic (FK) Scheme developed in [16,17] to higher-order polynomial reconstructions of the distribution function with respect to the position variable x. This class of methods applies a semi-Lagrangian discretization [8–11,24,28,29,31,32,45–47,51] to the transport operator for each velocity component of a discrete velocity model (DVM) [4,38,40, 41]. Because it does not rely on remapping after each time step, the FK approach produces less numerical dissipation than classical semi-Lagrangian methods. Indeed, in the simple

123

Journal of Scientific Computing

case of a pure transport problem, the method is exact, whereas the dissipation introduced by standard semi-Lagrangian approaches will cause solutions to degenerate toward a constant value over very long time scales. The FK approach has been shown to be an efficient way of simulating kinetic equations when the transport operator is linear, enabling the simulation of full six-dimensional problems on a single-processor machine. However, the original implementation is based on a piecewise constant reconstruction of the distribution function for each velocity of the lattice. Thus solutions computed with this method are limited to first-order in space and time. In this paper, we extend the semi-Lagrangian approach of the original FK method to compute high-order solutions of the transport operator, and we couple these solutions via operator splitting to an approximation of the collision operator that, although first-order, is amenable to computation with fast spectral methods. After a general introduction, we consider the construction of the method both for the case of the BGK (Bhatnagar–Gross–Krook) operator [30] as well as for the general Boltzmann operator [6] by using piecewise-linear approximations of the distribution function to evaluate the transport operator. In a second part of the work, we observe numerically that this is enough to obtain second-order accuracy in space for all collisional regimes in the BGK case and nearly second-order in the Boltzmann case. The remainder of this article is organized as follows. In Sect. 2, we present the model and recall the Fast Kinetic Scheme for piecewise-constant reconstructions. In Sect. 3, we extend the method to piecewise-linear reconstructions. Results of several test problems are presented in Sect. 4 to assess the efficiency of the new method in comparison to more classical schemes. Conclusions and a discussion of future developments are given in Sect. 5.

2 Kinetic Equations and Fast Kinetic Scheme In this section, we briefly present the model and the numerical method designed in [16]. We refer the reader to papers cited below for a more complete description.

2.1 The Kinetic Model In the kinetic theory of rarefied gases, the state of a neutral particle system is characterized by a non-negative function f such that f (x, v, t) gives the density (with respect to the measure dvd x) of particles located at position x ∈ Ω ∈ Rdx and moving with velocity v ∈ Rdv at time t ∈ R+ . The evolution of f is assumed to be governed by a kinetic equation of the form ∂f + v · ∇x f = Q[ f ]. (1) ∂t Here the collision operator Q on the right-hand side of (1) describes the effects of particle interactions; its form depends on the details of the underlying microscopic dynamics. Independently of these details, Q satisfies the conservation property  Q[g](v)φ(v) dv = 0, ∀g ∈ Dom(Q), (2) R3

where Dom(Q) is the domain of definition of Q and the components of φ(v) = (1, v, |v|2 ) are the collision invariants. In addition, Q[g] = 0 if and only if g = M[g], ∀g ∈ Dom(Q),

(3)

123

Journal of Scientific Computing

where the Maxwelliandistribution M[g] = M[g](x, v, t) can be expressed in terms of the moments Ug (x, t) :=

R3

g(x, v, t)φ(v) dv: M[g](x, v, t) = E [Ug (x, t)](v),

(4)

with E specified below in (7). Integrating (1) against φ(v) yields a system of macroscopic conservation laws  ∂U f v f φ(v) dv = 0, (5) + ∇x · ∂t R3 which is not closed because U f does not provide enough information about f to evaluate the integral in (5). However, the fact that f ≈ M[ f ] near local thermal equilibrium motivates a fluid approximation for (5) in which f is replaced by E [U ], where U = (U0 , U1 , U2 ) satisfies the Euler equations ∂U + ∇x · F(U ) = 0, ∂t

 with F(U ) =

Rdv

(6)

v E [U ]φ(v) dv. In the case of point-wise particles with no internal struc-

ture, E [U ] is given by E [U ](v) =

  ρ |v − u|2 , exp (2π T )dv /2 2T

(7)

where the density ρ, velocity u, and temperature T are related to U by the (invertible) mapping  1 ρ|v|2 + dv ρT . (8) 2 A simple collision operator satisfying (2) and (3) is the Bathnagar–Gross–Krook (BGK) operator [30] U0 = ρ, U1 = ρu, U2 =

Q BG K [ f ] = ν(M[ f ] − f ),

(9)

where ν = ν(x, t) > 0 is the collision frequency. On the other hand, the classical Boltzmann operator [6] reads       QB( f ) = B (|v − v∗ |, ω) f (v  ) f v∗ − f (v) f (v∗ ) dv∗ dω, (10) Rdv

Sdv −1

where ω is a vector of the unit sphere Sdv −1 ⊂ Rdv . The velocities (v  , v∗ ) are related to the velocities (v, v∗ ) by 1 1 (11) (v + v∗ + |q|ω), v∗ = (v + v∗ − |q|ω) , 2 2 where q = v − v∗ is the relative velocity. The kernel B in (10) characterizes the details of binary interactions. It has the general form v =

B(|v − v∗ |, ω) = |v − v∗ |σ (|v − v∗ |, cos θ ),

(12)

where σ is the scattering cross-section, θ is the scattering angle between relative velocities and       v − v∗ · (v − v∗ ) v − v∗ · ω cos(θ ) = = . (13) |v − v∗ |2 |v − v∗ |

123

Journal of Scientific Computing

If the force between particles is given by an inverse power law with exponent k, then σ (|v − v∗ |, cos θ ) = bα (cos θ )|v − v∗ |α−1 ,

(14)

with α = (k − 5)/(k − 1). In the Maxwell pseudo-molecules model, k = 5 so that B is independent of the relative velocity: B(v, v∗ , ω) = b0 (cos θ ).

(15)

For numerical purposes, a widely used model is the variable hard sphere (VHS) model introduced by Bird [1]. In this model, bα (cos θ ) = Cα , where Cα is a positive constant, and hence σ (|v − v∗ |, cos θ ) = Cα |v − v∗ |α−1 .

(16)

In what follows, the formulation of the numerical scheme is first done in the general case of a collision operator which only satisfies conditions (2) and (3). Details are provided for the BGK model in Sect. 3.2 and the Boltzmann model in Sect. 3.3. For numerical simulations, we focus first on the BGK operator, which is often used as a proxy for the more complicated Boltzmann collision operator [6], with one physical dimension (dx = 1) and one velocity dimension (dv = 1). We then consider the Boltzmann model for Maxwell molecules with one physical dimension (dx = 1) and two velocity dimensions (dv = 2).

2.2 The Original Fast Kinetic Scheme The Fast Kinetic Scheme (FK) [16] belongs to the family of semi-Lagrangian schemes [10,11, 24], which are typically applied to a discrete velocity model (DVM) [4,38] approximation of the original kinetic equation. We briefly recall the basics of the scheme in the one-dimensional setting and refer the reader to [16] for details. We truncate the velocity space with some appropriate bounds and set a uniform grid with N points and spacing Δv. We then consider a vector whose components approximate the distribution function f at locations vk : f k (x, t) ≈ f (x, vk , t), k = 1, . . . , N .

(17)

The discrete velocity model is a set of N evolution equations for f k , 1 ≤ k ≤ N , of the form ∂t f k + vk · ∇x f k = Q k (f),

(18)

where f = ( f 1 , . . . , f N ) and Q k (f) is a suitable approximation of Q[ f ](vk ). Each equation in (18) is solved by a time splitting procedure. We recall here a first-order splitting approach in which the numerical solution of the left-hand side after a single step of the transport operator serves as the initial condition for the collision operator on the right-hand side. Transport stage −→ ∂t f k + vk · ∇x f k = 0,

(19a)

Collision stage −→ ∂t f k = Q k (f).

(19b)

A common issue with all discrete velocity models is that Q k does not satisfy a discrete version of (2). In other words,  φk Q k (f) Δv  = 0, (20) k

123

Journal of Scientific Computing

Vk+2

Vk+2 k

k

Vk+1

n+1

f k+2,j

Vk+1

n

f k+2,j j

j n+1

Vk

f k+1,j

Vk

n

f k+1,j j

j

n

f k,j

n+1 f k,j

j

j

Fig. 1 Illustration of the transport scheme in one dimension for the first-order FK scheme. Left panel before transport stage, right panel after transport stage. Each discrete velocity (index k) drives its own transport equation with velocity vk . The representation of f is made by means of a piece-wise constant function. The shape of the entire function has not changed during the transport but the cell-centered values (bullets) may have

where φk = (1, vk , vk2 ) are the discretized collision invariants. As a consequence, the approximate moments  f (x, t) := U φk f k (x, t) Δv (21) k

are not conserved during the evaluation of a collision step. This problem arises in all numerical methods based on the discrete velocity models, and different strategies can be adopted to restore the correct macroscopic physical quantities [26,38,40,41]. Here we adopt a constrained optimization technique for the discretized distribution that, after each collision step, projects f into the space of discrete distributions with the correct moments. We do not enter into more details now but instead refer the reader to [26] for a description of the technique. The procedure is also used during initialization to ensure that the discrete distribution has the same moments as the initial kinetic distribution at time t = t 0 . We next introduce a uniform Cartesian grid in physical space made of M points with spacing Δx. Although we restrict ourselves here to uniform meshes in space, investigations are in progress to extend the present method to more complex geometries. We define a time discretization t n+1 = t n + Δt, n ≥ 0, where Δt is a time step that satisfies a CFL condition to be specified later. A single step of the FK method proceeds as follows. 0 be the approximate point-wise values at time t 0 of the distribution Transport stage Let f j,k f , i.e., f 0j = Pf in j , where P is the projection that results from the constrained optimization in = f (x , v , t 0 ). The idea is to solve the transport stage problem described above and f j,k j k (19a) continuously in space; see Fig. 1 for a sketch. To this aim, we define a piece-wise 0 0 for all x ∈ Ω , where Ω = (x constant function f k (x) = f j,k j j j−1/2 , x j+1/2 ) is a spatial 0 cell and Ω = j Ω j . Starting with data f k, j at time t 0 , the exact solution of (19a) is simply ∗,1

0

f k (x) = f k (x − vk Δt), for a.e. x ∈ Ω,

(22)

where the asterisk indicates that only the transport stage has been solved so far. In other words, 0 the entire function f k is advected with velocity vk for a time Δt. Applying this procedure to a generic time step n gives ∗,n+1

fk

123

n

(x) = f k (x − vk Δt),

for a.e. x ∈ Ω.

(23)

Journal of Scientific Computing n

The key point to observe in (23) is that the discontinuities of f k (x) are not aligned with the interfaces between adjacent cells. Instead, they depend entirely on the previous transport stage and thus they may be located anywhere in the physical space. In the absence of collisions, (23) gives the exact solution at each t n for initial data that is constant on each cell Ω j . Collision stage The effect of the collisions is to change the amplitude of f k (x). The idea is to solve the collision operator locally on the grid points and then extend these computed values to the full domain Ω. We solve the following ordinary differential equation ∂t f j,k = Q k (f j ),

where f j = ( f j,1 , . . . , f j,N ),

(24)

for all velocities vk , k = 1, . . . , N , and all grid points x j , j = 1, . . . , M. The initial data for solving this equation is furnished by the result of the transport stage obtained by (23) at points x j of the mesh at time t n+1 = t n + Δt. For example, if a forward Euler scheme is used, the approximation solution of (24) reads

 ∗,n+1 n+1 , (25) = f j,k + Δt Q k f ∗,n+1 f j,k j ∗,n+1

∗,n+1 where f j,k = fk (x j ). The equation in (25) furnishes a new value for f j at time t n+1 only in the cell centers of the spacial cells. However, one needs also the value of f j at all points of the space domain in order to perform the transport stage in the next time step. Therefore, we define a new piece-wise constant function Q k for each velocity of the lattice vk as

 n+1 ∗,n+1 ∗,n+1 , ∀x such that f k Q k (x) = Q k f ∗,n+1 (x) = f k (x j ). (26) j

Said differently, we make the fundamental assumption that for each k, the discontinuities of Q k coincide with those of f k . Using (26), one can rewrite the collision stage in terms of the spatially reconstructed functions as n+1

fk

∗,n+1

(x) = f k (x, t n + Δt) = f k

n+1

(x) + Δt Q k

(x).

(27)

This ends one time step of the FK scheme. As observed in [16], the transport scheme is stable for any Δt > 0. However, accuracy requirements dictate a constraint of the form   |vk | ≤ CFL = O(1). (28) Δt max k Δx Meanwhile, the time step constraint for the collision stage depends on the choice of the operator Q. To conclude this section, we note that the temporal accuracy of the scheme can be increased with high-order time splitting methods, while the spatial accuracy can be increased close to the fluid limit by employing high-order reconstruction of the fluid variables (see, for example, the details in [17] for a second-order implementation). In the next section, we introduce a procedure to improve the spatial accuracy for all regimes.

3 FK Schemes with Polynomial Reconstruction In this section, we extend the family of FK schemes to allow for piecewise polynomial reconstructions of the distribution function with respect to the spatial variable. The goal is to improve spatial accuracy while maintaining the efficiency of the original FK approach,

123

Journal of Scientific Computing

Vk+2

Vk+2 k

k n+1

n

f k+2,j

Vk+1

f k+2,j

Vk+1

j

Vk

j n f k+1,j

Vk

n+1

f k+1,j

j

j n+1

n

f k,j

f k,j j

j

Fig. 2 Illustration of the transport scheme in one dimension for the high-order FK scheme. Left panel before transport stage, right panel after transport stage. Each discrete velocity (index k) drives its own transport equation with velocity vk . The representation of f is made by means of piecewise polynomial function. The shape of the entire function has not changed during the transport but the cell-centered values (bullets) may have

in a way that is independent of the particular choice of collision operator. In Sect. 3.1, we present a general framework and then consider three possible approaches to reconstruct the collision operator in space. In Sect. 3.2, we discuss the details of the method for one of the resulting schemes (the one we consider most interesting) in the special case of piecewise linear reconstructions, and we implement it using the BGK operator. In Sect. 3.3, we extend the method to the Boltzmann operator, again with piecewise linear reconstructions. In Sect. 3.4, we preview work in progress to implement higher-order reconstructions.

3.1 General Approach As before, the starting point is a discrete velocity model for the original kinetic equation followed by an operator splitting of the collision and transport terms. Only a temporally first-order splitting is considered in this work; the extension to second-order Strang splitting appears straightforward, while the use of higher-order splittings is under study. We postpone such extensions to future work. Transport Stage We assume that at time t n we are given a piecewise polynomial reconstruction  f kn of a given order p for each velocity vk of the lattice. We assume further that the reconstruction is such that the support of each polynomial is an interval of length pΔx. Then, as for the piecewise constant case, we let ∗,n+1  (x) =  f kn (x − vk Δt), ∀x ∈ Ω, fk

(29)

be the exact solution of the transport stage (19a) after time Δt; see Fig. 2 for a sketch. In the absence of collisions, the next transport stage starts from where the previous one finished. This is in contrast to classical semi-Lagrangian schemes, which project the distribution function f back onto the piecewise polynomial space associated to the fixed mesh cells Ω j with j = 1, . . . , M. Collision Stage Given  f k∗,n+1 in (29), fix the lattice velocity vk , and let Ik be an interval of width pΔx such that  f k∗,n+1 is a polynomial of degree p when restricted to Ik . Let p+1 {yk,i }i=1 ⊂ Ik be the set of points at which we want to evaluate the collision operator. We call these control points. Using the forward Euler method, the collision stage at each yk,i is

123

Journal of Scientific Computing

n+1 (yk,i ), f kn+1 (yk,i ) =  f k∗,n+1 (yk,i ) + Δt Q k

(30)

n+1 (yk,i ) is a suitable approximation of the collision operator at yk,i . This approach where Q k differs substantially from standard semi-Lagrangian methods, which typically update the source terms at the cell centers of the fixed grid [49]. Here instead we allow for the freedom to choose the location of the control points for every velocity of the lattice vk in order to maintain local extrema. We stress that these points may and effectively, as detailed next, change with time during the time evolution of the problem. This particular choice permits to optimize the numerical method and reduce numerical diffusion. n+1 (yk,i ) is calculated. Three The only point that remains is to specify exactly how Q k possibilities are outlined below. The main reason for choosing one is to minimize the number of evaluations of the collision operator, since this is often the most expensive part of a kinetic simulation [22]. 1. Evaluate Q k at the cell centers and then interpolate Here the idea is to evaluate the collision operator at the cell centers {x j } M j=1 :   ∗,n+1 n+1 (x j ) = Q k  (x j ) , j = 1, . . . , M, (31) f Q k n+1 over Ω. and to use these values to construct a piecewise polynomial interpolant Q k The interpolant can then be evaluated directly at each yk,i . This approach requires the same number of collision operator evaluations M (the number of mesh points in physical space) as the original FK scheme. However, the interpolation procedure will likely induce the same loss of accuracy that one gets with a classical semi-Lagrangian discretization, since the interpolation points remain in the center of the cells. 2. Evaluate Q k directly at the control points In other words, set  ∗,n+1  n+1 (yk,i ) = Q k  Q (yk,i ) , i = 1, . . . , M. (32) f k Although a simple and natural choice, this approach is considerably more expensive when realistic collision operators are used. Indeed, fast methods for evaluating collision operators require a structured velocity grid, but the set of control points is in general different for different k. As a consequence, this approach may require as many as N × M evaluations of the collision operator, where the additional factor of N is due to the fact that the control points associated with different velocities do not overlap. As shown in [22], this number of evaluations is not practical for Boltzmann-like collision operators in multidimensional settings. 3. Evaluate Q k at the cell centers and apply the same reconstruction for the source and sink terms Specifically, set  ∗,n+1  n+1 (yk,i ) = Q k  (x(y ˆ k,i )) , (33) Q f k where x(y ˆ k,i ) = argmin x j |yk,i − x j |.

(34)

This choice is based on the fundamental assumption that the spatial variations in the source and sink terms of the collision operator should balance. It ensures that the collision operator needs to be evaluated only M times, and it is accurate in both the free-streaming and fluid limits, when the collision operator vanishes. In the following Section, we employ this method together with a piecewise linear reconstruction strategy for the distribution function. We observe in our numerical tests that this strategy yields a method that (1) is (nearly)second-order accurate across collisional regimes and (2) preserves the solution

123

Journal of Scientific Computing

Vk+2

Vk+2 k

k n+1

n

f k+2,j

Vk+1

f k+2,j

Vk+1

j

j

Vk

Vk

n f k+1,j

n+1

f k+1,j

j

j n

f k,j

n+1

f k,j

j

j

Fig. 3 Illustration of the transport scheme in 1D for the high order FK scheme with linear reconstruction. Left panel before transport stage, right panel after test transport stage. Each discrete velocity (index k) drives its own transport equation with velocity vk . The representation of f is a piecewise polynomial function. The shape of the entire function has not changed during the transport but the cell-centered values (bullets) may have

structure for very long times when compared to classical semi-Lagrangian approaches. However, as discussed later, the method may results in a loss of conservation and loss of positivity for which different remedies can be designed. Based on the discussion above, we proceed with the third option. After the collision operator is evaluated, Eq. (30) is used to update the values of the distribution function at each point yk,i . The new values of f kn+1 (yk,i ) are then used to define a polynomial interpolant  f kn+1 on Ik for the next transport stage. This completes one time step of the scheme.

3.2 Linear Polynomial Reconstruction Case for the BGK Model In this section, we focus on an implementation of the FK scheme that is based on a polynomial reconstruction of degree one. At the initial time we define, for each k and each interval (x j , x j+1 ), the piecewise linear reconstruction  f k0 (x) =

x j+1 − x 0 x − xj 0 f j,k + f j+1,k , Δx Δx

(35)

for all x ∈ (x j , x j+1 ). The transport stage has been detailed in Sect. 3.1; it can, however, be summarized by the sketch in Fig. 3 for piecewise linear reconstructions. For the sake of simplicity, we start by modeling the collisions with the BGK operator in (9). In this case, the collision stage can be evaluated exactly in time instead of using a forward Euler scheme. The result is 



 f kn+1 (x j ) = exp(−νΔt) f k∗,n+1 (x j ) + 1 − exp(−νΔt) Ek U ∗,n+1 , (36) j where Ek [U ∗,n+1 ] = E [U ∗,n+1 ](vk ) is a discretization of the Maxwellian distribution at the j j grid points that is uniquely defined by the moments = U ∗,n+1 j

 k

123

φk f k∗,n+1 (x j )Δv

(37)

Journal of Scientific Computing ambiguous definition

Vk

f

weighted average

SAME SLOPE

SAME SLOPE

M[f]

xj−1/2xj xj+1/2 linear reconstr. −

M[f]

−

+ +

min/max chosen values

fj

M[ fj ]

+

xj−1/2

xj+1/2

xj−1/2xj xj+1/2

Fig. 4 Illustration of the reconstruction technique employed in the R-FK scheme—from the distribution function data one reconstructs an upwind piece-wise linear polynomial (green on the left panel). From the distribution functions known at the cell centers we evaluate the Maxwellian also at cell centers (blue cross on the right panel). Last, on the top panel, we employ the same slope from the distribution function which passes through the Maxwellian cell center value. An ambiguous definition at each node is clarified by choosing the less extreme values in case the slopes are of different signs, a weighted averaged value otherwise, see the black/red bullets (Color figure online)

via the consistent relation 

  Δv = U (x j , t ∗,n+1 ). φk Ek U ∗,n+1 j

(38)

k

Equation (36) furnishes the new value of the distribution f at time t n+1 only at the centers of the spatial cells for each velocity vk ; see the green bullets in Fig. 4. We still need to identify, for each k, the points at which the collision operator is evaluated. These control points should be defined for each velocity of the lattice vk and for each interval Ik upon which the reconstruction is a linear polynomial. A natural strategy is to align them with cell centers at t = 0, i.e., yk,i (t = 0) = x j for some j, and then let them move with the advective flow in each transport step. Updating the solution at these points means that the distribution function f is always updated at local extrema. In general, the locations of these extrema are different for each vk : at the initial time, they are all located at cell centers; however, during each transport stage, they are shifted by vk Δt in space and therefore become misaligned (see Fig. 3). Once the kinetic distribution is known, evaluating the BGK collision operator amounts to evaluating the Maxwellian distribution. Thus the general strategy presented in Sect. 3.1 (strategy 3) amounts to setting the slopes in the reconstruction of Maxwellian equal to the slopes in the reconstruction of the kinetic distribution. However, in order to enforce continuity of the reconstruction at yk,i and to avoid the creation of extrema that lead to spurious oscillations, an additional limiter is applied. The entire procedure is depicted in Fig. 4. More specifically, for each evaluation point yk,i , there exists a j such that x j ≤ yk,i < x j+1 . Given this j, we set   +,n+1 + ak, j+1 (yk,i − x j+1 ) Ek (yk,i ) = Ek U ∗,n+1 (39) j+1 and −,n+1

Ek

  + ak, j (yk,i − x j ), (yk,i ) = Ek U ∗,n+1 j

(40)

123

Journal of Scientific Computing

where ak, j and ak, j+1 are the slopes of f k to the left and right of yk,i , respectively. In general, +,n+1 Ek  = Ek−,n+1 ; refer to the top picture in Fig. 4. Thus we set ⎧ 

x −y −,n+1 ⎪ Ek (yk,i ) j+1Δx k,i ⎪ ⎪ ⎪ 

⎪ ⎪ ⎨ + E +,n+1 (yk,i ) yk,i −x j if ak, j ak, j+1 > 0, k Δx 

Ekn+1 (yk,i ) = (41) +,n+1 −,n+1 ⎪ (yk,i ), Ek (yk,i ) if ak, j ≤ 0, ak, j+1 ≥ 0, ⎪ ⎪ max Ek ⎪  ⎪ ⎪ ⎩ min E +,n+1 (yk,i ), E −,n+1 (yk,i ) if ak, j ≥ 0, ak, j+1 ≤ 0. k k Finally, the new distribution f kn+1 at any extreme point yk,i is defined by



f kn+1 (yk,i ) = exp(−νΔt) f k∗,n+1 (yk,i ) + 1 − exp(−νΔt) Ekn+1 (yk,i ).

(42)

These points are used to define a new piecewise linear interpolant: for each k and each interval (yk,i , yk,i+1 ), we set yk,i+1 − x n+1 x − yk,i n+1  f kn+1 (x) = f k (yk,i ) + f k (yk,i+1 ) (43) Δx Δx for all x ∈ (yk,i , yk,i+1 ). This reconstruction then serves as the initial condition for the transport stage in the next time step. We conclude this section with two remarks: – An alternative to (41) for defining Ekn+1 (yk,i ) is to (i) reconstruct the macroscopic quantity U (yk,i ) from the knowledge of the distribution function values in these points and then (ii) substitute U (yk,i ) into (7) and evaluate the Maxwellian. However, this strategy will be quite expensive in the case of the BGK operator since it will require the computation of N × M Maxwellian distribution and it will not be applicable to more general operators, such as the Boltzmann operator, that cannot be defined solely in terms of the macroscopic quantities. – While the reconstruction in (41) does not guarantee conservation of macroscopic quantities, numerical tests suggest that the loss in conservation is quite small. When strict conservation is necessary, one can recover the correct moments by applying the constrained optimization technique that was designed in [26] and mentioned briefly in Sect. 2.2 above to the distribution function f˜kn+1 .

3.3 Linear Polynomial Reconstruction Case for the Boltzmann Model In this section, we focus on the construction of the FK scheme in the case of the Boltzmann model and, as in the previous section, we make use of a polynomial reconstruction of degree one. Thus, at the initial time we define, for each k and each interval (x j , x j+1 ), the piecewise linear reconstruction x j+1 − x 0 x − xj 0  f k0 (x) = f j,k + f j+1,k , (44) Δx Δx for all x ∈ (x j , x j+1 ). The transport stage, depicted in Fig. 3, is identical to the BGK case. The collision stage is similar and requires two pieces: (i) an algorithm to evaluate the Boltzmann operator at the cell centers and (ii) a method that approximation the operator at the control points based on the values at the cell centers. Once the collision operator is evaluated at the the control points, we employ a forward Euler scheme to integrate the equation. The method here described can be extended to other time integrators to increase the temporal order or to handle stiff dynamics [19–21].

123

Journal of Scientific Computing

3.3.1 Evaluation of the Boltzmann Operator We use a fast spectral approach [3,22,25,26,39] to approximate the Boltzmann operator. While the algorithm is presented here for completeness, the reader is referred to the references for further explanation and details. We focus on a single cell x j at a given instant of time t n in the description of the fast spectral method. The same computation is repeated for all cells x j , j = 1, . . . , M and times t n , n = 0, . . . , Nt . The algorithm is based on the Carleman representation of Q B [22,25]   ˜ QB( f ) = B(z, y)δ(z · y) Rdv

Rdv

[ f (v + y) f (v + z) − f (v + z + y) f (v)] dz dy, where ˜ B(|z|, |y|) = 2dv −1 σ



|z|2

+ |y|2 ,





|z| |z|2

(45)

+ |y|2

(|z|2 + |y|2 )−

dv −2 2

.

(46)

and, in order to lighten the notation, we have suppress the x and t dependence in all the variables, writing for example f = f (v). The same holds true for the other variables involved in the approximation of the collision operator, only the velocity dependence is kept. The first step is to truncate the integration domain of the Boltzmann integral (10). If f has compact support on the ball B0 (R) of radius R centered in the origin, then to avoid √ aliasing [39], we restrict the velocity domain to the cube [−T , T ]dv , where T ≥ (2 + 2)R. We then set f (v) = 0 on [−T , T ]dv \ B0 (R) and extend it to a periodic function on the set [−T , T ]dv . To simplify the notation, we may assume (upon rescaling of v) that T = π so √ that R = λπ with λ = 2/(3 + 2). We denote by Q BR ( f ) the Boltzmann operator with the cut-off introduced. We then approximate f by a truncated Fourier series by f N (v) =

N /2 

fˆk eik·v ,

k=−N /2

fˆk =

1 (2π)dv

 [−π,π ]dv

f (v)e−ik·v dv.

(47)

where the summation is understood in a multi-index sense. We then obtain a spectral quadrature of our collision operator by projecting (10) with the Carleman representation for Q B , onto the space of trigonometric polynomials of degree less or equal to N . The result is ∂t fˆk = Qˆ k , k = −N /2, . . . , N /2.

(48)

where, using the Carleman representation,  Q BR ( f N )e−ik·v dv Qˆ k = [−π,π ]dv

=

N /2 

βˆ F (l, m) fˆl fˆm , k = −N /2, . . . , N /2,

(49)

l,m=−N /2 l+m=k

with βˆ F (l, m) = Bˆ F (l, m) − Bˆ F (m, m) and   ˜ Bˆ F (l, m) = B(x, y) δ(x · y) ei(l·x+m·y) d x d y. B0 (R)

B0 (R)

(50)

123

Journal of Scientific Computing

From the spectral representation of the collision operator, one can separate the so-called gain from the loss part. This gives N /2 

Qˆ k =

βˆ F (l, m) fˆl fˆm = Pˆk − Lˆ k , k = −N /2, . . . , N /2

(51)

l,m=−N /2 l+m=k

with N /2 

Pˆk =

Bˆ F (l, m) fˆl fˆm ,

Lˆ k =

N /2 

Bˆ F (m, m) fˆl fˆm .

l,m=−N /2

l,m=−N /2

l+m=k

l+m=k

(52)

Finally, to reduce the number of operations in (49) we look for a convolution structure of the approximated operator (49). In other words, we search to approximate βˆ F (l, m) by a sum βˆ F (l, m) 

A 

α p (l)α p (m),

(53)

p=1

where A represents the number of finite possible directions of collisions. The formula in (53) reduces the evaluation of the collision integral to a sum of A discrete convolutions. In total, the sum requires O(A N 3 log2 N 3 ) operations when using standard FFT techniques [25,39]. ˜ However, in order to obtain the convolution structure, B(x, y) must be decomposable, i.e., ˜ B(x, y) = a(|x|) b(|y|).

(54)

Such a condition is satisfied when B˜ is constant. For the kernel (12), ˜ B(x, y) = 2dv −1 Cα (|x|2 + |y|2 )−

dv −α−2 2

,

(55)

which is constant if dv = 2 and α = 0 or if dv = 3 and α = 1. We consider here the case dv = 2 and α = 0, in which case B˜ = 1. By writing z and y in spherical coordinates: z = ρe and y = ρ  e , one gets  R   R    1  iρ(l·e) iρ  (m·e )  ˆ δ(e · e ) e dρ e dρ de de . (56) B F (l, m) = 4 S1 S1 −R −R R Then, denoting φ 2R (s) = −R eiρs dρ, for s ∈ R, we have the explicit formula φ 2R (s) = 2 R Sinc(Rs),

(57)

ˆ where Sinc(x) = sin(x) x . This explicit formula is further plugged in the expression of B F (l, m) and using its parity property, this yields  π Bˆ F (l, m) = φ 2R (l · eθ ) φ 2R (m · eθ +π/2 ) dθ. (58) 0

Finally, a regular discretization of A equally spaced points, which is spectrally accurate because of the periodicity of the function, gives A π  2 φ R (l · eθ p ) φ 2R (m · eθ p +π/2 ), Bˆ F (l, m) = M p=1

where θ p = π p/A.

123

(59)

Journal of Scientific Computing

3.3.2 Reconstruction The algorithm in the previous section is used to find a numerical approximation of Q B at each cell center x j and time step t n . To extend the approximation from the cell centers to the control points, we recall the assumption made in Sect. 3.1. That is the spatial variation of the source and sink (52) terms balance and therefore the slopes of the numerical representation are locally zero ((33)). For the BGK case, such an assumption was equivalent to setting the local slopes of the Maxwellian equal to that of f . However, to ensure a continuous profile for the following transport stage, an addition limiter was imposed on the Maxwellian distribution value ((41)). For the Boltzmann case, a similar limiting procedure can be applied directly to f in the control points. Let

 f k+,n+1 (yk,i ) = f k∗,n+1 (yk,i ) + Δt Pk+,n+1 (yk,i ) − L +,n+1 (y ) , (60) k,i k (yk,i ) are computed the solution of the Boltzmann equation where Pk+,n+1 (yk,i ) and L +,n+1 k using the solution of the Boltzmann operator on the center of the cell which lies on the right of the control point yk,i :



 n+1 n+1 Pk+,n+1 (yk,i ) − L +,n+1 (y ) = P (x ) − L (x ) = Q(x j+1 ). (61) k,i j+1 j+1 k k k Let also f k−,n+1 (yk,i ) = f k∗,n+1 (yk,i ) + Δt



 Pk−,n+1 (yk,i ) − L −,n+1 (yk,i ) , k

(62)

the solution of the Boltzmann equation where Pk−,n+1 (yk,i ) and L −,n+1 (yk,i ) are computed k using the solution of the Boltzmann operator on the center of the cell which lies on the left of the control point yk,i :



 Pk−,n+1 (yk,i ) − L −,n+1 (yk,i ) = Pkn+1 (x j ) − L n+1 (63) k k (x j ) = Q(x j ). Then, the new value of the distribution function in this control point is computed in a similar fashion as the case of the BGK model as 

⎧ x j+1 −yk,i −,n+1 ⎪ f (y ) k,i ⎪ k Δx ⎪ 

⎪ ⎨ yk,i −x j +,n+1 + f if ak, j ak, j+1 > 0, (y ) n+1 k,i k Δx (64) f k (yk,i ) = ⎪ +,n+1 −,n+1 ⎪ (yk,i ), f k (yk,i )) if ak, j ≤ 0, ak, j+1 ≥ 0, ⎪ max( f k ⎪ ⎩ min( f k+,n+1 (yk,i ), f k−,n+1 (yk,i )) if ak, j ≥ 0, ak, j+1 ≤ 0. These updated points are used to define a new piecewise linear interpolant: for each k and each interval (yk,i , yk,i+1 ), we set yk,i+1 − x n+1 x − yk,i n+1  f kn+1 (x) = f k (yk,i ) + f k (yk,i+1 ) Δx Δx

(65)

for all x ∈ (yk,i , yk,i+1 ). This reconstruction then serves as the initial condition for the transport stage in the next time step. Since, the approach proposed does not guarantee the conservation of the macroscopic quantities, one possible solution which permits to restore the exact conservation consists in projecting the reconstructed distribution function  f kn+1 (x) into the space of functions which preserve the macroscopic moments and which are close in the L 2 sense to the original function. This step could be performed by employing a constrained optimization technique for the discretized distribution that projects f into the space of discrete

123

Journal of Scientific Computing

distributions with the correct moments [26]. However, we do not employ such a technique in this work and postpone to future investigations a possible extension of the proposed method to handle perfect conservation.

3.4 About High Order Polynomial Reconstructions In this section, we briefly discuss extensions of the method to polynomial reconstructions of order two and higher. The transport step and evaluation of the collision operator at cell centers are essentially the same; we refer to the previous subsections for those details. The main differences when going to higher order come from the reconstruction procedure. Below we outline the main issues to be addressed, leaving specifics of the method as well as implementation and testing for future work. The major considerations are as follows: – Selection of control points For the piecewise linear reconstruction, the control points are aligned with the cell centers x j at t = 0 and are then displaced by vk Δt at each time step. It is easy to see that these points are local extrema and thus a natural choice. For higher order polynomials, this is no longer the case. Thus, while it is still possible to use this set of points, it is not the only choice. One may instead adjust the control point to line up with local extrema. – Interpolation For the piecewise linear reconstruction, there is also a natural choice of interpolation. However for higher polynomials, one must decide among several choice. For example, in the case of parabolic reconstruction, each control point will belong to three different three-point stencils, each of which determines a different parabolic fit, and a strategy must be devised to choose among them. One strategy for doing so is to choose a reconstruction which minimizes the total variation at the cell centers. Moreover, additional limiting may be incorporated to ensure that the local extrema are advected with the flow, so that the control points can be chosen in the same way as is done in the piecewise linear reconstruction. – Enforcing continuity As in the case of piecewise linear reconstruction, the piecewise constant representation of the of the collision operator will induce discontinuities at the control points. Thus a limiting strategy similar to (64) must be employed at the control points in order to enforce continuity there.

4 Numerical Experiments In this section, we perform several numerical experiments for the FK scheme using piecewise linear reconstructions (referred to as R-FK: reconstructed FK) applied to the kinetic equation (1) with both the BGK collision operator and the Boltzmann collision operator. We compare the R-FK scheme with the original FK scheme (referred to simply as FK), which uses piecewise constants to approximate the transport operator, and also with two standard semi-Lagrangian (SL) schemes: a first-order scheme (referred to as ‘SL-Upwind’) that uses piece-wise linear reconstructions and a second-order scheme (referred to as ‘SL-MUSCL’) that employs a piece-wise quadratic reconstruction with a Van Leer flux limiter [37,49]. In terms of accuracy, we expect the R-FK scheme to perform better than either of the first-order schemes and at least as well as the SL-MUSCL scheme. In terms of computational cost, we expect the R-FK scheme to be less expensive than the SL-MUSCL scheme. Indeed, in an ideal situation, the cost of R-FK should be close to that of FK, the latter of which has already been shown to be a cheaper alternative to the classical SL strategy [16].

123

Journal of Scientific Computing

We consider first the BGK model. We test convergence on a problem with a smooth solution and observe second-order convergence for a wide range of collision frequencies. We also compare the four semi-Lagrangian schemes on two benchmark problems: a Riemann problem with Sod-like initial data and a problem with highly oscillatory initial data. The second test is performed in a strongly kinetic regime and therefore requires a scheme with very little dissipation in order to maintain accurate results. We perform a mesh convergence study, in physical space and we perform a cost analysis with space and velocity refinements to further elucidate this point. We also demonstrate, for the oscillatory problem, that R-FK is the most efficient of the four schemes; that is, it requires the least computational cost to achieve the same error. This result portends well for more realistic computations in two and three dimensions. The final test is a planar Couette flow, which requires non-trivial boundary conditions. We simulate this problem with several different values of the collision frequency and compare the results of the R-FK scheme with those of the second-order SL method. We also measure the numerical convergence of the R-FK scheme and observe that, in the case, the convergence rate oscillates, but on average is second-order or higher. In the second part, we discuss numerical results for the Boltzmann equation. We compare the R-FK method with piecewise linear reconstructions with the second order MUSCL-SL method for a Riemann type problem. As done for the BGK model, we also measure the numerical convergence rate on a smooth solution. Results show that also for this model second-order convergence is attained by the proposed method.

4.1 Convergence Test Problem for the BGK Model To investigate the convergence of R-FK, we consider a problem with a smooth solution, defined on a spatial domain Ω = [0, 1] and truncated velocity domain L v = [− 15, 15]. The boundary conditions are periodic in x; and the initial condition is specified by a Maxwellian distribution (7) with density, mean velocity and temperature profiles that are smoothly varying in space: ρ(x, t = 0) = 1 +

1 1 sin(2π x), u(x, t = 0) = 0, T (x, t = 0) = 5 + sin(2π x). 2 2 (66)

The simulation is run to a final time tfinal = 0.025. The meshes are uniform in both physical and velocity spaces: M () = 100 ×(2 ) cells in physical space with  = 0, 1, 2, 3, 4, 5, 6 and N = 50 points in velocity space. In Fig. 5, we plot the initial data and the final solution for ν = 101 , 102 , 104 . Density, mean velocity and temperature are displayed for the numerical solution with M (0) = 100 cells. Due to our initialization, the cell centers of mesh  are also cell centers of mesh  + 1, for  = 0, 1, 2, 3, 4, 5. Consequently we can estimate the numerical order of convergence using three consecutive meshes with indices ,  + 1,  + 2 and  = 0, 1, 2, 3, 4 [37]: ⎛

pρ()

⎞ | ⎜ ⎟ j=1 ⎟, = log2 ⎜ ⎝  M () (+1) ⎠ (+2) |ρ2 j − ρ4 j |  M ()

j=1

()

⎛

(+1)

|ρ j − ρ2 j

()

pT

⎞ | ⎜ ⎟ j=1 ⎟, = log2 ⎜ ⎝  M () (+1) ⎠ (+2) |T2 j − T4 j |  M ()

()

|T j

(+1)

− T2 j

j=1

(67)

123

Journal of Scientific Computing

1.2

5.2

1

5

0.6

4.9 4.8

0.4

4.7

-0.2

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x

1.6

Density Velocity Temperature

1.4 1.2

0.2

-0.4

4.4 1

-0.6

5.2

1.6

5.05

0.4

5

0.2 -0.2

4.9

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x Density Velocity Temperature

1.4 1.2

4.95

0

4.95

0 -0.2

5.1

0.6

5

0.4

4.5

0.8

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x

5.25 5.2

5.1

0.8 0.6

5.05

0.4

5

0.2

4.95

0

4.9 4.85

-0.4 4.85 1

4.85 1

5.15

1

-0.2

4.9

-0.4 -0.6

0.6

5.15

1

5.1 5.05

0.8

4.6

0

5.15

-0.6

Temperature

5.1

Density Velocity Temperature

Temperature

5.3

0.8

0.2

Density, Velocity

1.4

Temperature

1

1.6

5.4

Temperature

Density, Velocity

1.2

5.5

Density, Velocity

Density Velocity Temperature

1.4

Density, Velocity

1.6

0

0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 x

4.8 1

Fig. 5 BGK Model, convergence test case, results for M (0) = 100 cells — Top panel: initial data (left), final solution for ν = 10 (right). Bottom panel: final solution for ν = 102 (left), final solution for ν = 104 (right)

for the density and temperature, respectively. We also measure the convergence rate of the distribution function f , using the same meshes and the following formula: ⎛

 M ()

N

⎜ () p f = log2 ⎜ ⎝

k=1 N

k=1

() | f j,k j=1

 M () j=1

−

(+1)

(+1) f 2 j,k | (+2)

| f 2 j,k − f 4 j,k |

⎞ ⎟ ⎟ ,  = 0, 1, 2, 3, 4. ⎠

(68)

In order to properly measure the error due to the space discretization, we make the time discretization error very small by fixing the time step Δt to the value needed for the finest mesh; that is, Δt = (Δx)(6) / max(|v|)  10−5 , where (Δx)(6) = 1/6400. Moreover, in order to observe the convergence rate for different regimes, we consider six different collision frequencies: ν = 10−2 , 10−1 , 100 , 101 , 102 and 104 . Results of the convergence tests are gathered in Tables 1, 2 and 3. We observe that for all quantities of interest and for all six collisional frequencies, the nominal orders of accuracy are close to two.

4.2 Riemann Problem for the BGK Model We next consider a classical Riemann problem with Sod-like initial data [48]. The spatial domain is again Ω = [0, 1] and the truncated velocity domain is L v = [−20, 20]. The initial condition is generated by a Maxwellian distribution (7) with piecewise constant fluid variables:

123

Journal of Scientific Computing Table 1 Convergence rate of the density in the BGK model for a smooth solution computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (67) and for six collisional frequencies Meshes M () , M (+1) , M (+2)

()

Order of convergence, pρ ν = 10−2

ν = 10−1

ν = 100

ν = 101

ν = 102

ν = 104

100, 200, 400

1.927

1.927

1.927

1.928

1.930

2.395

200, 400, 800

2.043

2.043

2.043

2.041

2.050

2.311

400, 800, 1600

2.028

2.028

2.028

2.028

2.027

2.157

800, 1600, 3200

2.002

2.002

2.002

2.002

2.002

2.089

1600, 3200, 6400

1.998

1.998

1.998

1.998

1.997

2.048

Table 2 Convergence rate of the temperature in the BGK model for a smooth solution computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (67) and for six collisional frequencies Meshes M () , M (+1) , M (+2)

()

Order of convergence, pT ν = 10−2 ν = 10−1 ν = 100

ν = 101

ν = 102

ν = 104

100, 200, 400

0.942

0.943

0.954

1.078

1.915

2.996

200, 400, 800

3.074

3.073

3.059

2.874

2.103

2.698

400, 800, 1600

1.805

1.805

1.809

1.795

1.906

2.326

800, 1600, 3200

2.183

2.183

2.178

2.141

2.073

2.269

1600, 3200, 6400

1.974

1.974

1.972

1.974

1.984

2.154

Table 3 Convergence rate of the kinetic distribution in the BGK model for the smooth solution test problem. The rate is computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (68) and for six collisional frequencies Meshes M () , M (+1) , M (+2)

()

Order of convergence, p f ν = 10−2

ν = 10−1

ν = 100

ν = 101

ν = 102

ν = 104

100, 200, 400

1.950

1.950

1.950

1.951

1.955

2.467

200, 400, 800

2.021

2.021

2.021

2.018

2.024

2.315

400, 800, 1600

2.024

2.024

2.024

2.024

2.018

2.156

800, 1600, 3200

2.005

2.005

2.005

2.005

2.007

2.094

1600, 3200, 6400

1.998

1.998

1.998

1.996

1.997

2.042

ρ(x, t = 0) = 1, u(x, t = 0) = 0, T (x, t = 0) = 5,

if x ≤ 1/2,

ρ(x, t = 0) = 0.125, u(x, t = 0) = 0, T (x, t = 0) = 4,

if x > 1/2.

The same initial conditions are used to run three simulations corresponding to three different collision frequencies: ν = 102 , 103 and 104 . In each case, the simulation is run to a final time tfinal = 0.07. The meshes are uniform in both the physical and velocity spaces, with N = 50 and M = 300, respectively. The CFL condition is the same for every test so that Δt = 0.5Δx/ max(|v|), where max(|v|) = 20. Simulation results for the density, mean velocity, and temperature are shown in Figs. 6, 7, and 8, respectively, for collision frequencies ν = 104 (top panels), ν = 103 (middle panels),

123

Journal of Scientific Computing 1

0.65

SL-Upwind SL-MUSCL FKS R-FKS

0.9 0.8

SL-Upwind SL-MUSCL FKS R-FKS

0.6 0.55 0.5

0.7

0.45

0.6

0.4 0.5

0.35

0.4

0.3

0.3

0.25

0.2

0.2

0.1

0.15

0

0.1

0.2

0.3

0.4

0.5

0.6

1

0.7

0.8

0.9

1

0.9 0.8

0.54

0.56

0.58

0.6

0.62

0.65

SL-Upwind SL-MUSCL FKS R-FKS

0.64

0.66

SL-Upwind SL-MUSCL FKS R-FKS

0.6 0.55 0.5

0.7

0.45

0.6

0.4 0.5

0.35

0.4

0.3

0.3

0.25

0.2

0.2

0.1

0.15

0

0.1

0.2

0.3

0.4

0.5

0.6

1

0.7

0.8

0.9

1

0.9 0.8

0.54

0.56

0.58

0.6

0.62

0.55

SL-Upwind SL-MUSCL FKS R-FKS

0.64

0.66

SL-Upwind SL-MUSCL FKS R-FKS

0.5 0.45

0.7 0.6

0.4

0.5

0.35

0.4

0.3

0.3 0.25

0.2 0.1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

0.54

0.56

0.58

0.6

0.62

0.64

0.66

Fig. 6 Riemann problem for the BGK model—solution at tfinal = 0.07 for the density. Left full solution, right zoom close to the region in which the solution develops a contact discontinuity in the limit of infinite collisions. Collision frequency ν = 104 (top), 103 (middle), 102 (bottom)

and ν = 102 (bottom panels). On the left column, we plot the solution over the entire domain, while on the right side, we focus on interesting parts of the solution profiles. The density profile in Fig. 6 is magnified around the location of the contact discontinuity that appears in the fluid dynamic limit (x  0.59). The mean velocity profile in Fig. 7 is magnified near the location of the shock wave that also forms in the fluid limit (x  0.86). For the temperature profile in Fig. 8, we zoom around the location of the contact discontinuity when ν = 104 and around the region of the shock wave for ν = 103 and ν = 102 .1 The results obtained by the four schemes are reported on each plot: FK in red, R-FK in blue, SL-Upwind in magenta and SL-MUSCL in green. 1 We consider different locations to highlight areas where we observe the largest differences in the methods.

123

Journal of Scientific Computing 1.4

1.4

1.2

1.2

1

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2 SL-Upwind SL-MUSCL FKS R-FKS

0 -0.2

SL-Upwind SL-MUSCL FKS R-FKS

0

0.1

0.2

0.3

0.4

0.5

0 0.6

0.7

0.8

0.9

1

1.6 1.4

-0.2 0.83 1.4

0.84

0.85

0.86

0.87

1.2

0.9

1

1 0.8

0.8

0.6

0.6

0.4

0.4

0.2

SL-Upwind SL-MUSCL FKS R-FKS

0 0

0.1

0.2

0.3

0.4

0.5

0.2 0.6

0.7

0.8

0.9

1

1.8

0 0.83 1

1.6

0.9

1.4

0.8

1.2

0.7

1

0.6

0.8

0.5

0.6

0.4

0.4

0

0.1

0.2

0.3

0.4

0.5

0.84

0.85

0.86

0.87

0.88

0.89

0.9

SL-Upwind SL-MUSCL FKS R-FKS

0.3

SL-Upwind SL-MUSCL FKS R-FKS

0.2 0

0.89

SL-Upwind SL-MUSCL FKS R-FKS

1.2

-0.2

0.88

0.2 0.6

0.7

0.8

0.9

1

0.1 0.83

0.84

0.85

0.86

0.87

0.88

0.89

0.9

Fig. 7 Riemann problem for the BGK model—solution at tfinal = 0.07 for the mean velocity. Left full solution, right zoom close to the region in which the solution develops a shock wave in the limit of infinite collisions. Collision frequency ν = 104 (top), 103 (middle), 102 (bottom)

The results show that SL-MUSCL and R-FK are more accurate than SL-Upwind and FK in the highly collisional case ν = 104 . From the magnifications on the right side of Figs. 6, 7 and 8 (top panels), it is clear that the second-order schemes require fewer points to resolve discontinuities than do the first-order schemes. We also observe that for the density, the R-FK scheme resolves the contact discontinuity better than the SL-MUSCL scheme, while for the temperature and mean velocity, the two solutions are almost superimposed. Similar behavior is observed for the shock wave in the fluid dynamic limit: the solution is resolved better (less points on the shock wave are present) by the R-FK scheme than by the other three. On the other hand, the complexity of R-FK is comparable to SL-Upwind since in both cases first-order polynomials are employed. As ν decreases, the wave fronts begin to smooth out and the differences between the schemes becomes less pronounced. However, for ν = 103 we still observe a difference

123

Journal of Scientific Computing 9

9

SL-Upwind SL-MUSCL FKS R-FKS

8 7

7

6

6

5

5

4

4

3

3

2

0

0.1

0.2

0.3

9

0.4

0.5

0.6

0.7

0.8

0.9

1

2

7

6

6

5

5

4

4

3

3 0

0.1

0.2

0.3

8

0.4

0.5

0.6

0.7

0.8

0.9

1

0.58

0.6

0.62

0.64

0.66

0.68

0.7

0.64

0.66

0.68

0.7

0.62 0.64 0.66 0.68

0.7

SL-Upwind SL-MUSCL FKS R-FKS

2

0.54

7.5

SL-Upwind SL-MUSCL FKS R-FKS

7

0.56

8

7

2

0.54

9

SL-Upwind SL-MUSCL FKS R-FKS

8

SL-Upwind SL-MUSCL FKS R-FKS

8

7 6.5

0.56

0.58

0.6

0.62

SL-Upwind SL-MUSCL FKS R-FKS

6

6

5.5 5

5 4.5

4

4 3.5

3

3 2

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

2.5

0.54 0.56 0.58

0.6

Fig. 8 Riemann problem for the BGK model—solution at tfinal = 0.07 for the temperature. Left full solution, right zoom close to the region in which the solution develops a contact discontinuity (top) or a shock wave (middle and bottom) in the limit of infinite collisions. Collision frequency ν = 104 (top), 103 (middle), 102 (bottom)

between the four solutions; in particular, the semi-Lagrangian schemes are more diffusive than the FK schemes.2 In addition, we observe that the R-FK scheme produces slightly less diffusive results than the FK scheme does. Finally, when ν = 102 , kinetic effects are strong and the four solutions are almost superimposed. At this point, it is difficult to appreciate the differences between the methods, although the SL-Upwind profiles are still a bit less sharp than the rest.

2 A mesh convergence test, not reported, confirms that the wave locations for all schemes converge to the

same location when the number of points increases.

123

Journal of Scientific Computing Initial velocity

1

Initial velocity

1

0.5

0.5

0

0

-0.5

-0.5

-1

-1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

0.25

0.3

0.35

0.4

Fig. 9 Initial mean velocity field for the highly oscillating kinetic test case—δ = 0.02 and 600 mesh cells are used—left: full view on Ω, right: zoom on [0.2, 0.4]

4.3 Highly Oscillating Kinetic Problem for the BGK Model The test problem in this section is designed to illustrate the benefit that characterizes the FK method: reduced numerical dissipation when compared to classical semi-Lagrangian schemes in a kinetic regime. The spatial domain is Ω = [0, 1] and the velocity domain is truncated to L v = [− 25, 25]. Initially the gas is in a Maxwellian state (cf. (7)) with uniform density and temperature: ρ(x, t = 0) = 1, T (x, t = 0) = 5,

(69)

but a discontinuous, highly oscillatory mean velocity field with period δ  L = 1. Given M , we set the spatial mesh {xi }i=1 ⎧ ⎨ +1 if x ∈ (x2k − δ/2, x2k + δ/2) ∩ (0.25; 0.75), (70) u(x, t = 0) = −1 if x ∈ (x2k+1 − δ/2, x2k+1 + δ/2) ∩ (0.25; 0.75), ⎩ 0 else for all integers k ∈ [0, M/2]. A depiction of the initial mean velocity field is provided in Fig. 9. We set the collision frequency to ν = 102 so that kinetic effects are non-negligible. In this situation, the mean velocity field induces waves that emanate from each discontinuity. These waves eventually interact, creating secondary waves that lead to further interactions and additional waves. If a scheme introduces too much dissipation, then the structure of these waves is lost. To assess the performance of each scheme, we perform a mesh convergence study and an efficiency study.

4.3.1 Mesh Convergence Study The final time is fixed at tfinal = 0.025, and the period of the initial velocity data is set to δ = 0.02. Then an increasing sequence of spatial meshes: M = 600, 1200, 2400 and 4800 is used for each of the four schemes, while the velocity domain and the number of velocity cells are kept fixed. In Fig. 10, we present the density, mean velocity and temperature profiles for the four schemes when 600 cells are used. The results show that the schemes perform comparably well in capturing the large-scale structures in the solution. To investigate the small-scale structures in the solution, we zoom in on each of the macroscopic variables. Density, mean velocity, and temperature profiles are given in Figs. 11, 12, and 13, respectively. From these figures, we observe that most of the small-scale structures in

123

Journal of Scientific Computing 1.04

SL-Upwind SL-MUSCL FKS R-FKS

1.03 1.02 1.01 1 0.99 0.98 0.97

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.15

0.8

0.9

1

SL-Upwind SL-MUSCL FKS R-FKS

0.1 0.05 0 -0.05 -0.1 -0.15

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

6.2

0.8

0.9

1

SL-Upwind SL-MUSCL FKS R-FKS

6 5.8 5.6 5.4 5.2 5

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

Fig. 10 Oscillating problem for the BGK model with 600 spatial cells—solution at tfinal = 0.025 for the density (left), mean velocity (middle) and temperature (right) with ν = 102 and the SL-Upwind, SL-MUSCL, FK and R-FK schemes

the SL-Upwind solution have been damped out by numerical dissipation. A similar, but less pronounced effect can be observed in the solution generated by SL-MUSCL which, due to the embedded limiter, reduces to first-order at local extrema (other limiters that clip extrema will have a similar effect). As the mesh is refined, SL-Upwind begins to capture the physical oscillations induced by the initial condition. However, the stair-case profile is never fully recovered, even at the finest mesh. On the other hand, both FK schemes are able to maintain the small-scale structures with relatively coarse meshes.

4.3.2 Cost and Efficiency Study In the one-dimensional case, all of the schemes are relatively inexpensive to run. Nonetheless, it is important to assess the cost of each one, with an eye toward multi-dimensional situations. We do this in two ways: first we measure the CPU time at fixed mesh, regardless of the

123

Journal of Scientific Computing 1.032

1.032

SL-Upwind SL-MUSCL FKS R-FKS

1.03 1.028

1.028

1.026

1.026

1.024

1.024

1.022

1.022

1.02

1.02

1.018

1.018

1.016 0.77

0.78

0.79

0.8

0.81

0.82

1.032

0.83

0.84

0.85

1.016 0.77

1.028

1.026

1.026

1.024

1.024

1.022

1.022

1.02

1.02

1.018

1.018 0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.79

0.8

0.81

0.82

0.85

1.016 0.77

0.83

0.84

0.85

SL-Upwind SL-MUSCL FKS R-FKS

1.03

1.028

1.016 0.77

0.78

1.032

SL-Upwind SL-MUSCL FKS R-FKS

1.03

SL-Upwind SL-MUSCL FKS R-FKS

1.03

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

Fig. 11 Oscillating problem for the BGK model—zooms on density for 600, 1200, 2400 and 4800 cells (from top-left to bottom-right)—solution at tfinal = 0.025 with ν = 102 for the SL-Upwind, SL-MUSCL, FK and R-FK schemes

accuracy obtained; then we measure the CPU time needed to achieve a prescribed level of accuracy. Because an exact solution does not exist for this problem, we use a reference solution to decide when a given scheme has attained an acceptable level of accuracy. Cost at fixed meshes We present in Table 4 the cost of the simulations from the mesh convergence study in physical space. For each simulation, we report the number of time steps Ncycle and the CPU time T in second(s). We also compute the time per cycle, Tcycle = T/Ncycle , and time per cycle per physical cell, Tcell = T/Ncycle /Nc . Lastly, we calculate the ratio of CPU time for each scheme to the CPU time of SL-Upwind. From this table, we see that for a fixed mesh (i) FK is slightly faster than SL-Upwind; (ii) SL-MUSCL is about 1.4 times slower than SL-Upwind and FK; (iii) R-FK is about 2.3 times slower than SL-Upwind and FK; and (iv) R-FK is about 1.6 times slower than SL-MUSCL. In Table 5, we report the results when the spacial mesh is fixed to Nc = 2400 cells while the number of cells in the velocity space is refined: Nv = 50, 100, 200 and 400. The other simulation parameters are kept the same. We report the value of the CPU time T, and the number of time steps Ncycle and the time needed to update one velocity cell: Tvel.cell = T/Ncycle /Nv . The costs to update one velocity cell for the FK scheme and the SL-Upwind methods are comparable ( roughly 3 × 10−4 s), while the SL-MUSCL and R-FK are roughly twice as expensive (7 × 10−4 s). The same conclusions of the space mesh refinement can be drawn: the costs of the FK and SL-Upwind are close, while R-FK and SL-MUSCL are two times more expensive. In both cases, the FK class of methods perform better in terms of precision than standard semi-Lagrangian methods. Cost at fixed accuracy While there is no exact solution to the problem with oscillating initial data, we have observed that the SL-MUSCL and the FK family convergence to the same solution as the physical mesh is refined. Based on this converged solution, we estimate that

123

Journal of Scientific Computing 0.13

SL-Upwind SL-MUSCL FKS R-FKS

0.125 0.12 0.115

0.13

0.12 0.115

0.11

0.11

0.105

0.105

0.1

0.1

0.095

0.095

0.09

0.09

0.085

0.085

0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84

0.13

0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84

0.13

SL-Upwind SL-MUSCL FKS R-FKS

0.125

SL-Upwind SL-MUSCL FKS R-FKS

0.125

0.12

0.12

0.115

0.115

0.11

0.11

0.105

0.105

0.1

0.1

0.095

0.095

0.09

SL-Upwind SL-MUSCL FKS R-FKS

0.125

0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84

0.09

0.66 0.68 0.7 0.72 0.74 0.76 0.78 0.8 0.82 0.84

Fig. 12 Oscillating problem for the BGK model—zooms on velocity for 600, 1200, 2400 and 4800 cells (from top-left to bottom-right)—solution at tfinal = 0.025 with ν = 102 for the SL-Upwind, SL-MUSCL, FK and R-FK schemes 5.78

SL-Upwind SL-MUSCL FKS R-FKS

5.76 5.74 5.72

5.78

5.74 5.72

5.7

5.7

5.68

5.68

5.66

5.66

5.64

5.64

5.62

5.62

5.6 0.28 0.29

0.3

5.78

0.31 0.32 0.33 0.34 0.35 0.36 0.37 SL-Upwind SL-MUSCL FKS R-FKS

5.76 5.74 5.72

5.6 0.28 0.29

0.3

5.78

0.31 0.32 0.33 0.34 0.35 0.36 0.37 SL-Upwind SL-MUSCL FKS R-FKS

5.76 5.74 5.72

5.7

5.7

5.68

5.68

5.66

5.66

5.64

5.64

5.62

5.62

5.6 0.28 0.29

SL-Upwind SL-MUSCL FKS R-FKS

5.76

0.3

0.31 0.32 0.33 0.34 0.35 0.36 0.37

5.6 0.28 0.29

0.3

0.31 0.32 0.33 0.34 0.35 0.36 0.37

Fig. 13 Oscillating problem for the BGK model—zooms on temperature for 600, 1200, 2400 and 4800 cells (from top-left to bottom-right)—solution at tfinal = 0.025 with ν = 102 for the SL-Upwind, SL-MUSCL, FK and R-FK schemes

123

50

50

50

50

SL-Upwind

SL-MUSCL

FK

R-FK

220 441 882 1764 220 441 882 1764 220 441 882 1764 220 441 882 1764

600 × 50 = 3 × 104

1200 × 50 = 6 × 104

2400 × 50 = 12 × 104

4800 × 50 = 24 × 104

600 × 50 = 3 × 104

1200 × 50 = 6 × 104

2400 × 50 = 12 × 104

4800 × 50 = 24 × 104

600 × 50 = 3 × 104

1200 × 50 = 6 × 104

2400 × 50 = 12 × 104

4800 × 50 = 24 × 104

600 × 50 = 3 × 104

1200 × 50 = 6 × 104

2400 × 50 = 12 × 104

4800 × 50 = 24 × 104

[−25, 25]

[−25, 25]

[−25, 25]

[−25, 25]

Cycle Ncycle

Cell # Nc × Nv

Vel

121.552 s

29.612 s

7.404 s

1.852 s

42.124 s

10.408 s

2.636 s

0.644 s

81.716 s

16.960 s

4.596 s

1.132 s

58.380 s

13.368 s

3.076 s

0.752 s

Time T (s)

Time per cycle is obtained by Tcycle = T/Ncycle and time per cycle per cell by Tcell = T/Ncycle /Nc

Nv

Method

1.0 1.0 1.51 1.49

6.18 × 10−6 8.58 × 10−6 8.68 × 10−6

Average →

1.27 1.40 1.42 0.57 0.57

9.65 × 10−6 8.73 × 10−6 4.88 × 10−6 4.98 × 10−6

0.04463 Average →

0.61 0.52 0.57 2.88 2.81

4.97 × 10−6 6.18 × 10−6 1.40 × 10−5 1.40 × 10−5

0.02388 Average →

2.85 2.89 2.86

1.44 × 10−5 1.41 × 10−5

0.0689 Average →

0.03356

1.40 × 10−5

0.01379

0.00842

0.01180

4.92 × 10−6

0.00598

0.00293

0.01923

8.01 × 10−6

0.01042

0.00515

0.01516

1.0

1.0

5.81 × 10−6 6.89 × 10−6

1.0

5.70 × 10−6

0.03310

Time ratio vs SL-Upwind

T/cell Tcell (s)

6.32 × 10−6

0.00698

0.00342

T /cycle Tcycle (s)

Table 4 Oscillating test case for the BGK model, results for 600, 1200, 2400 and 4800 cells for the SL-Upwind, SL-MUSCL, FK and R-FK schemes

Journal of Scientific Computing

123

123

2400

2400

2400

2400

SL-Upwind

SL-MUSCL

FKS

R-FKS

882 891 895 897 882 891 895 897 882 891 895 897 882 891 895 897

2400 × 50 = 12 × 104

2400 × 100 = 24 × 104

2400 × 200 = 48 × 104

2400 × 400 = 96 × 104

2400 × 50 = 3 × 104

2400 × 100 = 6 × 104

2400 × 200 = 12 × 104

2400 × 400 = 24 × 104

2400 × 50 = 12 × 104

2400 × 100 = 24 × 104

2400 × 200 = 48 × 104

2400 × 400 = 96 × 104

2400 × 50 = 12 × 104

2400 × 100 = 24 × 104

2400 × 200 = 48 × 104

2400 × 400 = 96 × 104

[−25, 25]

[−25, 25]

[−25, 25]

[−25, 25]

Cycle Ncycle

Cell # Nc × Nv

Vel

246.85 s

123.39 s

61.40 s

29.61 s

88.02 s

43.86 s

21.49 s

10.41 s

481.38 s

121.65 s

46.22 s

16.96 s

108.61 s

55.25 s

28.64 s

13.368 s

Time T (s)

1.0 1.0 1.27

3.09 × 10−4 3.85 × 10−4

Average →

0.79 0.81 0.78 2.22 2.14

3.09 × 10−4 6.71 × 10−4 6.89 × 10−4

0.09813 Average →

2.23 2.27 2.22

6.88 × 10−4 6.84 × 10−4

0.27520 Average →

0.13787

6.89 × 10−4

0.06891

0.03357

0.04901

3.03 × 10−4

0.75

3.21 × 10−4 3.09 × 10−4

0.02412

2.38 0.78

Average → 0.01180

3.03 × 10−4

4.43

1.34 × 10−3

0.53666

7.31 × 10−4

1.61 2.20

6.80 × 10−4

0.13592

0.05187

5.19 × 10−4

0.01923

0.06173

1.0

1.0

3.21 × 10−4 3.03 × 10−4

1.0

3.03 × 10−4

0.12108

Time ratio versus SL-Upwind

T/vel.cell Tvel.cell (s)

3.09 × 10−4

0.03214

0.01516

T/cycle Tcycle (s)

Time per cycle is obtained by Tcycle = T/Ncycle and time per cycle per velocity cell by Tvel.cell = T/Ncycle /Nv

Nc

Method

Table 5 Oscillating test case for the BGK model, results for Nc = 2400 spacial cells and Nv = 50, 100, 200 or 400 cells in velocity space for the SL-Upwind, SL-MUSCL, FKS and R-FK schemes

Journal of Scientific Computing

Journal of Scientific Computing Table 6 Oscillating test case for the BGK model—estimation of the CPU time T needed to capture the plateaus of estimated characteristics length δl

N δl

SL-Upwind

SL-MUSCL

FK

R-FK

150 4.8 × 10−3

15 4.8 × 10−3

2 4.8 × 10−3

2 4.8 × 10−3

Nc = 2N /δl

62500

6250

833

833

Ncycle = (0.3675) Nc

22969

2297

306

306

Tcell (s)

6.18 × 10−6

8.73 × 10−6

6.18 × 10−6

1.41 × 10−5

T = Nc × Ncycle × Tcell

8871.52s

125.34s

1.58s

3.60s

(2h28mn)

(2mn5s)

The schemes diffuse discontinuities over N cells. The required number of cells is Nc = 2N /δl, the number of time steps Ncycle = (0.3675) Nc , where 0.3675 = 1764/4800 the ratio time steps/Nc observed in Table 4, T = Nc × Ncycle × Tcell . Tcell is chosen as the average of the times reported in Table 4

the characteristic distance between two discontinuities is roughly Δl = 4.8 × 10−3 . (This corresponds to 23 cells of size Δx = 1/4800.) We use this distance to estimate the time needed for each scheme to resolve the plateaus between discontinuities. We compute this estimate as follows, using the SL-MUSCL scheme as an illustrative example: – From the SL-MUSCL results, we approximate the transition layer between two plateaus in the SL-MUSCL scheme to be roughly N = 15 cells wide. Assuming that we desire a transition layer that is half as wide as the plateau and recalling that spatial domain has length L = 1, we estimate that the SL-MUSCL scheme requires about Nc ≥

L N = 6250 Δl/2

(71)

cells to resolve a plateau. – According to Table 4, the ratio between the number of time steps Ncycle and the number of spatial cells Nc is 0.3675 (1764/4800 = 882/2400 = 441/1200  220/600). Hence the number of required time steps is Ncycle = (0.3675) Nc . – Finally, given the average time Tcell  8.73 × 10−6 s from Table 4 needed to update one cell during one time step, we can estimate the CPU time needed by the SL-MUSCL scheme via equation T = Nc × Ncycle × Tcell  125.34s = 2min 5s.

(72)

In Table 6 we report values of T for all four schemes, along with the other parameters used to calculate these values. As expected, numerical dissipation of the SL-Upwind and SL-MUSCL schemes drastically affects their efficiency. Specifically, they would need approximately 2.5 h and 2 min, respectively, to resolve the oscillations in this test problem. The FK family, on the other hand, requires only 2–4 s. We conclude this section by stating that the computational cost analysis is a difficult task to accomplish and the numbers reported in the tables are are only meant to provide a rough idea of the amount of resources that can already be saved in 1D by the FK family of schemes. More significant gains are expected when 2D×2D and, more importantly, 3D×3D simulations are performed.

123

Journal of Scientific Computing

4.4 Planar Couette Flow for the BGK Model In this part, we discuss a Couette type flow. The physical setting considered is one dimension in physical space and two dimensions in velocity space, i.e., f (x, v, t) ∈ R × R2 × R+ . The spatial domain is Ω = [0, 1] and the truncated velocity domain is L v = [− 8, 8]2 . The initial condition is generated by a Maxwellian distribution (7) with constant fluid variables: ρ(x, t = 0) = 1, u x (x, t = 0) = 0, u y (x, t = 0) = 0, T (x, t = 0) = 1. The same initial conditions are used to run three simulations corresponding to three different collision frequencies: ν = 50, ν = 100 and ν = 500. In each case, the simulation is run until a stationary solution is reached, in practice when t = tfinal = 15. The meshes are uniform in both the physical space, with M = 80 points, and velocity space, with N = 322 points. The CFL condition is the same for every test so that Δt = 0.5 Δx/ max(|v|), where max(|v|) = 8. In the Couette problem, the right boundary is moving at constant speed u y (0, t) = 1 and the wall has a constant temperature T (0, t) = 1. On the other hand, the left boundary mimics a plate at rest with constant temperature T (1, t) = 1. These conditions cause mean velocity of the gas in the y-direction to evolve until it reaches an monotonically increasing, equilibrium profile. Simulation results for the density, mean velocity in the y-direction, and temperature are shown in Fig. 14 for the three different collision frequencies. On the left column, we plot the solution using the R-FK method while on the right columns we report the solution obtained with the second order SL-MUSCL method. Results of the two different schemes match closely for the regimes tested. It is important to stress that the first order FK method is not able to reproduce correctly the solution to the Couette problem while the second order R-FK one is. This is due to the low order reconstruction of the distribution function in space which does not allow to catch enough accurately the solution at the boundaries. We finally measure the numerical convergence rate of the R-FK method in this particular setting, i.e., for the Couette flow situation. For performing such study, we run the same simulation as above with uniform meshes in both the velocity and physical spaces with N = 322 points in velocity space and M () = 20 × (2 ) cells in space with  = 0, 1, 2, 3, 4 for ν = 50, ν = 100 and ν = 500. We estimate the numerical order of convergence using three consecutive meshes with indexes ,  + 1,  + 2 and  = 0, 1, 2 [37] using Eq. (67). In order to properly measure the error due to the space discretization, we make the time discretization error very small by fixing the time step Δt to the value needed for the finest mesh; that is, Δt = 0.5(Δx)(6) / max(|v|), where (Δx)(6) = 1/640. Results of the convergence tests are gathered in Tables 7 and 8. We observe that for the two quantities of interest and for the three collisional frequencies, the nominal orders of accuracy are oscillating but they remain close to two. This is an interesting result since, as before stated, the first order FKS is not able to reproduce a Couette flow correctly due to its low accuracy. On the contrary the R-FKS is able to correctly approximate such a solution with an high order of accuracy.

4.5 Convergence Test for the Boltzmann Model To investigate the convergence of R-FK in the case of the Boltzmann operator, we consider a problem with a smooth solution similar to the one proposed for the BGK case. The spatial domain is Ω = [0, 1], the truncated velocity domain is L v = [− 8, 8]2 , the boundary

123

Journal of Scientific Computing 1.04

1.04 =50 =100 =500

=50 =100 =500

1.03

1.03

1.02

1.02

1.01

1.01

1

1

0.99

0.99

0.98

0.98 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.06

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.5

0.6

0.7

0.8

0.9

1

1.06 =50 =100 =500

=50 =100 =500

1.05

1.05

1.04

1.04

1.03

1.03

1.02

1.02

1.01

1.01

1

1

0.99

0.99 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

0

0.1

1 =50 =100 =500

0.9

=50 =100 =500

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0

0.1

0.2

0.3

0.4

Fig. 14 Couette flow, BGK model—solution at tfinal = 15 for the density (top panel), the temperature (middle panel) and the mean velocity in the y-direction (bottom panel). Collision frequency from ν = 50 to ν = 500. Left panels solution of the R-FK method, right panels solution of the second order MUSCL-SL method

conditions are periodic in x and the initial condition is specified by a Maxwellian distribution (7) with density, mean velocity and temperature profiles given by ρ(x, t = 0) = 1.5 +

1 1 sin(2π x), u(x, t = 0) = 0, T (x, t = 0) = . (73) 2 ρ(x, t = 0)

123

Journal of Scientific Computing Table 7 Convergence rate of the density for the Couette flow problem using the BGK model, computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (67) and for three collisional frequencies

Table 8 Convergence rate for the temperature for the Couette flow problem using the BGK model, computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (67) and for three collisional frequencies

Meshes M () , M (+1) , M (+2)

()

Order of convergence, pρ ν = 50

ν = 100

ν = 500

20,40,80

3.4316

3.0228

3.2420

40,80,160

1.2931

1.6771

1.4660

80,160,320

2.5326

2.5517

2.5232

160,320,640

4.1590

3.6960

4.1534

Meshes M () , M (+1) , M (+2)

Order of convergence, pT

()

ν = 50

ν = 100

ν = 500

20, 40, 80

2.7106

2.9785

2.5577

40, 80, 160

1.5706

1.8956

1.7148

80, 160, 320

6.2579

3.4912

4.1968

160, 320, 640

3.0973

3.2253

3.2724

Table 9 Convergence rate of the density in the Boltzmann for a smooth solution computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (67) and for five different Knudsen numbers. Boltzmann model Meshes M () , M (+1) , M (+2)

()

Order of convergence, pρ ε=1

ε = 10−1

ε = 10−2

ε = 10−3

ε = 10−4

20,40,80

1.1034

1.1017

1.0334

2.1021

1.8825

40,80,160

2.1920

2.1734

2.3264

2.4082

1.9527

80,160,320

1.9004

1.9401

1.7475

1.7128

2.5474

160,320,640

2.0066

1.9024

1.6922

2.0138

1.9557

The simulation is run to a final time tfinal = 0.01. The meshes are uniform in both the physical and velocity spaces with M () = 20 × (2 ),  = 0, 1, 2, 3, 4, 5, cells in physical space, N = 322 points in velocity space while A = 16 angles are employed. The CFL condition is fixed in such a way that Δt = 0.5Δx (6) / max(|v|) with Δx (6) = 1/640. Due to our initialization, the cell centers of mesh  are also cell centers of mesh  + 1, for  = 0, 1, 2, 3, 4. Consequently we estimate the numerical order of convergence using three consecutive meshes with indices ,  + 1,  + 2 and  = 0, 1, 2, 3 [37] using (67) for the density and temperature variables. In our study, we consider five different Knudsen numbers: ε = 1, 10−1 , 10−2 , 10−3 , 10−4 . Results of the convergence tests are gathered in Tables 9 and 10. We observe that for all quantities of interest the nominal orders of accuracy are close to two for all values of the Knudsen number except for ε = 10−2 . For the value of  the convergence order is around 1.8 for the finest mesh. Further numerical test are required to determined is full second-order convergence is achieved in this case.

123

Journal of Scientific Computing Table 10 Convergence rate of the temperature in the Boltzmann for a smooth solution computed from three successively refined meshes M () , M (+1) , M (+2) via Eq. (67) and for five different Knudsen numbers. Boltzmann model Meshes M () , M (+1) , M (+2)

()

Order of convergence, pT ε=1 ε = 10−1

ε = 10−2

ε = 10−3

ε = 10−4

20, 40, 80

1.4448

1.4399

1.4459

2.5499

2.5133

40, 80, 160

2.1462

2.1543

2.2073

1.9304

1.8278

80, 160, 320

1.9049

1.8850

1.8107

1.9088

2.3540

160, 320, 640

2.0130

1.9638

1.7671

2.6166

1.8599

4.6 Riemann Problem for the Boltzmann Equation We now consider a classical Riemann problem with Sod-like initial data [48] for the Boltzmann model. The setting is similar of the case of the BGK case: the spatial domain is Ω = [0, 1] and the truncated velocity domain is L v = [− 8, 8]2 . The Boltzmann equation is rescaled in such a way that a parameter ε appears in front of the collision operator. We then have ∂f 1 + v · ∇x f = Q[ f ]. ∂t ε

(74)

This parameter is called the Knudsen number and it is proportional to the mean free path between two consecutive collisions. This rescaling is done to make it easier the analysis of the scheme for the different regimes, as 1/ε plays a role analogous the the collision frequency ν in the BGK model. The initial condition is generated by a Maxwellian distribution (7) with piecewise constant fluid variables: ρ(x, t = 0) = 1, u(x, t = 0) = 0, T (x, t = 0) = 1,

if x ≤ 1/2,

ρ(x, t = 0) = 0.125, u(x, t = 0) = 0, T (x, t = 0) = 0.8,

if x > 1/2.

The same initial conditions are used to run three simulations corresponding to three different Knudsen numbers: ε = 10−3 , ε = 10−4 and ε = 5 × 10−5 . In each case, the simulation is run to a final time tfinal = 0.15. The choices done for the Knudsen number value reflect the fact that our aim is to explore the intermediate regimes. In fact, in standard kinetic theory studies, Knudsen number grater than 10−3 corresponds to the slip regime while a gas can be considered in free molecular regime only for Knudsen greater than 10. For smaller Knudsen numbers the effects of the collision term cannot be generally neglected. Consequently, when the Knudsen number is smaller than ε = 5 × 10−5 , the solution is very close to equilibrium, the distribution function is close to the Maxwellian distribution and results are analog to the BGK case. On the other hand, for Knudsen numbers larger than ε = 10−3 , the effects of the collision operator are very low and solutions are close to the free transport regime, where the scheme is, by construction, exact. The meshes are chosen to be uniform in both physical and velocity spaces, with M = 200 Δx and N = 322 , respectively. The corresponding CFL condition is Δt = min(0.5 max(|v|) , ρε ), where max(|v|) = 8. Concerning the collision operator discretization, we choose a number of discretization angles A = 8. This is enough to guarantee a good accuracy of the results in many situations [21,25]. The number of modes in Fourier space is the same as the number of points in velocity space, i.e. Nˆ = 322 .

123

Journal of Scientific Computing

Although the problem setting is like that of the BGK case, we expect to see differences in the solution. The main difference is that the velocity space is now two-dimensional, and the equation of state that relates the energy to the temperature in the fluid limit depends on dv . As a result the wave speeds in the fluid limit are different. Another difference is that the relaxation toward local thermal equilibrium (which eventually occurs in both models) occurs at different rates in the two models. As a result the physical diffusion is not the same. Specifically, we expect the Boltzmann model to relax at a slower rate so that the density, mean velocity and temperature profiles are less steep than in the BGK case. Simulation results for the density, mean velocity, and temperature are shown in Figs. 15, 16, and 17, respectively, for Knudsen numbers ε = 5 × 10−5 (top panels), ε = 10−4 (middle panels) and ε = 10−3 (bottom panels). As for the BGK case, on the left column, we plot the solution over the entire domain, while on the right side, we focus on the regions close to the waves (the magnifications are done close to the location where the shock wave forms in the fluid dynamic limit.) In these regions, it is easy to observe the differences between the schemes. R-FK results are displayed in red, while SL-MUSCL results are in blue. The results show that the R-FK is more accurate than SL-MUSCL method for all regime tested.

5 Conclusion and Perspectives In this paper we have presented a new class of semi-Lagrangian methods, based on polynomial reconstructions of the distribution function, that extends the Fast Kinetic Scheme developed in [16,17]. The basic strategy of the method is to update the extreme points of the distribution function, rather than the cells centers, in order to reduce the effects of numerical dissipation. In the first part of the paper, the method is presented in a general setting and several different approaches to compute the collision operator are proposed. Given the extreme cost of evaluating the collision operator, the best choice appears to compute the collision operator only once time per cell and to make the shape of the solution in space not dependent on the collision operator. This choice, which motivated by the observation that the collision operator acts only in velocity space, is later supported by numerical results. In the second part of the paper, details are given for the specific case of a piece-wise linear polynomial reconstruction for the BGK collision operator as well as for the Boltzmann operator. For the BGK operator case, the strategy for the collision operator amounts to the assumption that the distribution function and its associated Maxwellian share the same local variations in space for each velocity of the lattice. However, an additional limiter is needed in order to treat the discontinuities of the Maxwellian function at interfaces. For the Boltzmann case, the idea is similar, and we assume that the spatial variation of the distribution function is not affected by collisions, but only transport. In the third part of the paper, the new R-FK scheme is tested on four different problems for the BGK case. The first is a convergence test that exhibits second-order convergence across a wide range of collision frequencies. The second test, a Sod-like Riemann problem, shows that in the fluid regime, the results of the R-FK scheme are comparable to that of a standard second-order, semi-Lagrangian approach. For moderately collisional regimes, the new scheme performs better than the original FK scheme and the second-order semiLagrangian scheme. In a full kinetic regime, the differences between the FK and semiLagrangian schemes become negligible. The third test uses highly oscillatory initial data in a kinetic regime. The initial data generates many discontinuities that are damped by classical

123

Journal of Scientific Computing 1

0.25 R-FKS SL-MUSCL

R-FKS SL-MUSCL

0.9

0.8

0.7

0.2

0.6

0.5

0.4

0.15

0.3

0.2

0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

0.1 0.75

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

0.25 R-FKS SL-MUSCL

R-FKS SL-MUSCL

0.9

0.8

0.7

0.2

0.6

0.5

0.4

0.15

0.3

0.2

0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1

0.1 0.75

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

0.25 R-FKS SL-MUSCL

R-FKS SL-MUSCL

0.9

0.8

0.7

0.2

0.6

0.5

0.4

0.15

0.3

0.2

0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.1 0.75

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

Fig. 15 Riemann problem for the Boltzmann model—solution at tfinal = 0.15 for the density. Left full solution, right zoom close to the region in which the solution develops a contact discontinuity in the limit of infinite collisions. Knudsen numbers ε = 5 × 10−5 (top), ε = 10−4 (middle), ε = 10−3 (bottom)

123

Journal of Scientific Computing 0.8

0.8

R-FKS SL-MUSCL

R-FKS SL-MUSCL

0.7

0.7

0.6

0.6

0.5

0.5 0.4

0.4 0.3

0.3 0.2

0.2 0.1

0.1

0

0 0.75

-0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

0.8

0.8

R-FKS SL-MUSCL

R-FKS SL-MUSCL

0.7

0.7

0.6

0.6

0.5

0.5 0.4

0.4 0.3

0.3 0.2

0.2 0.1

0.1

0

0 0.75

-0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.9

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

0.8 R-FKS SL-MUSCL

0.8

R-FKS SL-MUSCL

0.7 0.7

0.6 0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2 0.1

0.1 0

-0.1 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0 0.75

0.8

0.85

Fig. 16 Riemann problem for the Boltzmann model—solution at tfinal = 0.15 for the velocity. Left full solution, right zoom close to the region in which the solution develops a shock wave in the limit of infinite collisions. Knudsen numbers ε = 5 × 10−5 (top), ε = 10−4 (middle), ε = 10−3 (bottom)

semi-Lagrangian schemes, but are maintained by the R-FK scheme due to its natural antidissipative formulation. The results obtained by the FK schemes are shown in this test to be superior in terms of accuracy and computational cost. As a final test, a problem with nontrivial boundary conditions, the Couette flow, is considered. For this problem, comparisons with a

123

Journal of Scientific Computing 1.5

1.5 R-FKS SL-MUSCL

R-FKS SL-MUSCL

1.4

1.4

1.3 1.3

1.2 1.2

1.1 1

1.1

0.9

1

0.8 0.9

0.7 0.8

0.6 0.5 0

0.2

0.4

0.6

0.8

1

1.5

0.7 0.75

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

1.5 R-FKS SL-MUSCL

1.4

R-FKS SL-MUSCL

1.4

1.3 1.3

1.2

1.1

1.2

1

1.1

0.9 1

0.8 0.9

0.7 0.8

0.6

0.5 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

1.4

0.7 0.75

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

1.5 R-FKS SL-MUSCL

R-FKS SL-MUSCL

1.3 1.4

1.2 1.3

1.1 1.2

1 1.1

0.9 1

0.8 0.9

0.7

0.8

0.6

0.5 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

0.7 0.75

0.76

0.77

0.78

0.79

0.8

0.81

0.82

0.83

0.84

0.85

Fig. 17 Riemann problem for the Boltzmann model—solution at tfinal = 0.15 for the temperature. Left full solution, right zoom close to the region in which the solution develops a contact discontinuity (top) or a shock wave (middle and bottom) in the limit of infinite collisions. Knudsen numbers ε = 5 × 10−5 (top), ε = 10−4 (middle), ε = 10−3 (bottom)

semi-Lagrangian method are performed and the numerical accuracy is measured. On average, a numerical second-order convergence rate is observed. A second battery of tests has been performed for the case of the Boltzmann model. The first test is a problem with a smooth solution, where near second-order convergence is observed

123

Journal of Scientific Computing

across collisional regimes. The second test is a Sod-like Riemann problem. In this test, as in the BGK case, the R-FK method provides better accuracy than standard semi-Lagrangian methods in all regimes tested. In the future, we will investigate extensions of this approach to the multi-dimensional case. While a similar reconstruction technique may apply, an analysis should be carried out to determine how to handle lines and surfaces of discontinuity in two and three dimensions. Finally, we intend to perform an analytical study of the convergence error of the proposed scheme in order to show that the numerically observed second order of convergence in space can be recovered theoretically.

References 1. Bird, G.A.: Molecular Gas Dynamics and Direct Simulation of Gas Flows. Clarendon Press, Oxford (1994) 2. Birsdall, C.K., Langdon, A.B.: Plasma Physics Via Computer Simulation. Series in Plasma Physics. Institute of Physics (IOP), London (2004) 3. Bobylev, A.V., Rjasanow, S.: Difference scheme for the Boltzmann equation based on the fast Fourier transform. Eur. J. Mech. B. Fluids 16(2), 293–306 (1997) 4. Bobylev, A.V., Palczewski, A., Schneider, J.: On approximation of the Boltzmann equation by discrete velocity models. C. R. Acad. Sci. Paris Ser. I. Math 320, 639–644 (1995) 5. Caflisch, R.E.: Monte Carlo and Quasi-Monte Carlo methods. Acta Numer. 7, 1–49 (1998) 6. Cercignani, C.: The Boltzmann Equation and Its Applications. Springer, New York (1988) 7. Cercignani, C., Illner, R., Pulvirenti, Mario: The Mathematical Theory of Dilute Gases, vol. 106. Springer Science & Business Media, New York (2013) 8. Chacón, L., del-Castillo-Negrete, D., Hauck, C.D.: An asymptotic-preserving semi-Lagrangian algorithm for the time-dependent anisotropic heat transport equation. J. Comp. Phys. 272, 719–746 (2014) 9. Cheng, C.Z., Knorr, G.: The integration of the Vlasov equation in configuration space. J. Comput. Phys. 22, 330–351 (1976) 10. Crouseilles, N., Respaud, T., Sonnendrucker, E.: A forward semi-Lagrangian method for the numerical solution of the Vlasov equation. Comp. Phys. Commun. 180(10), 1730–1745 (2009) 11. Crouseilles, N., Mehrenberger, M., Sonnendrucker, E.: Conservative semi-Lagrangian schemes for Vlasov equations. J. Comp. Phys. 229, 1927–1953 (2010) 12. Degond, P., Dimarco, G., Pareschi, L.: The moment guided Monte Carlo method. Int. J. Numer. Methods Fluids 67, 189–213 (2011) 13. Desvillettes, L., Mischler, S.: About the splitting algorithm for Boltzmann and BGK equations. Math. Mod. Methods Appl. Sci. 6, 1079–1101 (1996) 14. Dimarco, G., Hauck, C., Loubère, R.: Multidimensional high order FK schemes for the Boltzmann equation (in progress) 15. Dimarco, G.: The hybrid moment guided Monte Carlo method for the Boltzmann equation. Kin. Rel. Models 6, 291–315 (2013) 16. Dimarco, G., Loubère, R.: Towards an ultra efficient kinetic scheme. Part I: basics on the BGK equation. J. Comput. Phys. 255, 680–698 (2013) 17. Dimarco, G., Loubère, R.: Towards an ultra efficient kinetic scheme. Part II: the high order case. J. Comput. Phys. 255, 699–719 (2013) 18. Dimarco, G., Pareschi, L.: A fluid solver independent hybrid method for multiscale kinetic equations. SIAM J. Sci. Comput. 32, 603–634 (2010) 19. Dimarco, G., Pareschi, L.: Asymptotic preserving implicit-explicit Runge–Kutta methods for non linear kinetic equations. SIAM J. Numer. Anal. 49, 2057–2077 (2011) 20. Dimarco, G., Pareschi, L.: Exponential Runge–Kutta methods for stiff kinetic equations. SIAM J. Numer. Anal. 51, 1064–1087 (2013) 21. Dimarco, G., Pareschi, L.: Numerical methods for kinetic equations. Acta Numer. 23, 369–520 (2014) 22. Dimarco, G., Loubère, R., Narski, J., Rey, T.: An efficient numerical method for solving the Boltzmann equation in multidimensions. J. Comp. Phys. 353, 46–81 (2018) 23. Filbet, F., Russo, G.: High order numerical methods for the space non-homogeneous Boltzmann equation. J. Comput. Phys. 186, 457–480 (2003)

123

Journal of Scientific Computing 24. Filbet, F., Sonnendrücker, E., Bertrand, P.: Conservative numerical schemes for the Vlasov equation. J. Comput. Phys. 172, 166–187 (2001) 25. Filbet, F., Mouhot, C., Pareschi, L.: Solving the Boltzmann equation in N log2 N. SIAM J. Sci. Comput. 28(3), 1029–1053 (2007) 26. Gamba, I.M., Tharkabhushaman, S.H.: Spectral-Lagrangian based methods applied to computation of non-equilibrium statistical states. J. Comput. Phys. 228, 2012–2036 (2009) 27. Gamba, I.M., Haack, J.R., Hauck, C.D., Hu, J.: A fast spectral method for the Boltzmann collision operator with general collision kernels. SIAM J. Sci. Comput. 39, B658–B674 (2017) 28. Groppi, M., Russo, G., Stracquadanio, G.: High order semi-Lagrangian methods for the BGK equation. Commun. Math. Sci. 14, 389414 (2016) 29. Groppi, M., Russo, G., Stracquadanio, G.: Boundary conditions for semi-Lagrangian methods for the BGK model. Commun. Appl. Ind. Math 7(3), 138164 (2016) 30. Gross, E.P., Bathnagar, P.L., Krook, M.: A model for collision processes in gases. I. small amplitude processes in charged and neutral one-component systems. Phys. Rev. 94, 511–525 (1954) 31. Güçlü, Y., Hitchon, W.N.G.: A high order cell-centered semi-Lagrangian scheme for multi-dimensional kinetic simulations of neutral gas flows. J. Comput. Phys. 231, 3289–3316 (2012) 32. Güçlü, Y., Christlieb, A.J., Hitchon, W.N.G.: Arbitrarily high order convected scheme solution of the Vlasov–Poisson system. J. Comput. Phys. 270, 711–752 (2014) 33. Hauck, C.D., McClarren, R.G.: A collision-based hybrid method for time-dependent, linear, kinetic transport equations. Multiscale Model. Simul. 11, 1197–1227 (2013) 34. Homolle, T., Hadjiconstantinou, N.: A low-variance deviational simulation Monte Carlo for the Boltzmann equation. J. Comput. Phys. 226, 2341–2358 (2007) 35. Homolle, T., Hadjiconstantinou, N.: Low-variance deviational simulation Monte Carlo. Phys. Fluids 19, 041701 (2007) 36. Jin, S.: Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations. SIAM J. Sci. Comput. 21, 441454 (1999) 37. LeVeque, R.J.: Numerical Methods for Conservation Laws, Lectures in Mathematics. Birkhauser Verlag, Basel (1992) 38. Mieussens, L.: Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamic. Math. Mod. Methods App. Sci. 10, 1121–1149 (2000) 39. Mouhot, C., Pareschi, L.: Fast algorithms for computing the Boltzmann collision operator. Math. Comp. 75, 1833–1852 (2006) 40. Palczewski, A., Schneider, J.: Existence, stability, and convergence of solutions of discrete velocity models to the Boltzmann equation. J. Stat. Phys. 91, 307–326 (1998) 41. Palczewski, A., Schneider, J., Bobylev, A.V.: A consistency result for a discrete-velocity model of the Boltzmann equation. SIAM J. Numer. Anal. 34, 1865–1883 (1997) 42. Pareschi, L., Russo, G.: Numerical solution of the Boltzmann equation I: spectrally accurate approximation of the collision operator. SIAM J. Numer. Anal. 37, 12171245 (2000) 43. Pareschi, L., Toscani, G.: Interacting Multi-agent Systems. Kinetic Equations and Monte Carlo Methods. Oxford University Press, New York (2013) 44. Pareschi, L., Russo, G., Toscani, G.: Fast spectral methods for the Fokker PlanckLandau collision operator. J. Comput. Phys. 165, 216236 (2000) 45. Qiu, J.-M., Christlieb, A.: A Conservative high order semi-Lagrangian WENO method for the Vlasov equation. J. Comput. Phys. 229, 1130–1149 (2010) 46. Qiu, J.-M., Shu, C.-W.: Conservative semi-Lagrangian finite difference WENO formulations with applications to the Vlasov equation. Commun. Comput. Phys. 10, 979–1000 (2011) 47. Shoucri, M., Knorr, G.: Numerical integration of the Vlasov equation. J. Comput. Phys. 14, 8492 (1974) 48. Sod, G.: A survey of several finite difference methods for systems of nonlinear hyperbolic conservation laws. J. Comput. Phys. 27, 1–31 (1978) 49. Sonnendrücker, E., Roche, J., Bertrand, P., Ghizzo, A.: The semi-Lagrangian method for the numerical resolution of the Vlasov equation. J. Comput. Phys. 149, 201220 (1999) 50. Strang, G.: On the construction and the comparison of difference schemes. SIAM J. Numer. Anal. 5, 506–517 (1968) 51. Xiong, T., Qiu, J.-M., Xu, Z., Christlieb, A.: High order maximum principle preserving semi-Lagrangian finite difference WENO schemes for the Vlasov equation. J. Comput. Phys. 73, 618639 (2014)

123

Journal of Scientific Computing

Affiliations Giacomo Dimarco1

B

· Cory Hauck2 · Raphaël Loubère3

Giacomo Dimarco [email protected] Cory Hauck [email protected] Raphaël Loubère [email protected]

1

Department of Mathematics and Computer Science, University of Ferrara, Ferrara, Italy

2

Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, TN, USA

3

CNRS and Institut de Mathématiques de Bordeaux UMR 5251, Université de Bordeaux, 33 405 Talence, France

123