FULLY CONSERVATIVE SPECTRAL GALERKIN–PETROV METHOD FOR THE INHOMOGENEOUS BOLTZMANN EQUATION

. In this paper, we present an application of a Galerkin-Petrov method to the spatially one-dimensional Boltzmann equation. The three- dimensional velocity space is discretised by a spectral method. The space of the test functions is spanned by polynomials, which includes the collision invari- ants. This automatically insures the exact conservation of mass, momentum and energy. The resulting system of hyperbolic PDEs is solved with a ﬁnite volume method. We illustrate our method with two standard tests, namely the Fourier and the Sod shock tube problems. Our results are validated with the help of a stochastic particle method.


1.
Introduction. Numerically solving the Boltzmann equation has been a challenging problem for decades. Starting in the 1950s, first schemes had been proposed which were based on the Hilbert expansion or Monte Carlo integration of the collision operator, see [11] for a survey on modern deterministic schemes. Different authors also developed stochastic particle methods. The most popular among these is the Direct Simulation Monte Carlo (DSMC), which was introduced by Bird in the 1960s, see his monograph [3]. The DSMC, or modifications of it, are still in use today. Although particle methods can be applied to large range of problems, their usage can cause difficulties when considering slow flows, which require many repetitions in order to damp the stochastic fluctuations. That's why one is interested in competitive deterministic methods for the Boltzmann equation. Boltzmann or Boltzmann-like equations, The fastest schemes are based on the Fourier transform. This technique was first used by Pareschi and Perthame in [30], see also [31,32]. A first application to the inhomogeneous Boltzmann equation for oneand two-dimensional problems can be found in [15]. Stability and positivity for the homogeneous case were investigated in [13].
By the use of a Carlemann-like representation of the collision operator, this scheme was reformulated in [29] by Mouhot and Pareschi, now called the "Fast Spectral Scheme" (FSS) with quadratic complexity with respect to the number of degrees of freedom in case of the Maxwell molecules for the two-dimensional or the model of Hard Spheres for the three-dimensional collision operator. The computational cost was then reduced to almost linear complexity in [14] assuming the above restrictions on the collision kernel. In [38], Wu et al. applied FSS to more general collision kernels, including the Variable Hard Spheres model and approximations of the Lennard-Jones potential. This extended scheme was used for a systematic study of physically relevant collision kernels in [35], an analysis of micro- [37] and cavity flows in [36], as well as monoatomic mixtures [35] and polyatomic gases [39].
Recently, FSS was applied to the discretisation of the collision operator for threedimensional problems, both in space and velocity [10].
The Fourier representation of the collision operator in case of Maxwell pseudomolecules was used by Bobylev and Rjasanow [4], which was later extended in [5,6] and by Ibragimov and Rjasanow [24]. More recently, these results were extended to general cross-sections [19,20,17] by Gamba and Tharkabhushanam [19] and Gamba and Haack [17].
In general, Fourier-based methods only conserve mass and introduce a spectrally decaying error for smooth solutions on momentum and energy. Full conversation can be obtained by the use of Lagrange multipliers.
A discontinuous Galerkin discretisation based on subdivisions of a truncated velocity space was presented in [1], which was reduced to quadratic complexity in [2] using discrete Fourier transforms.
With an a priori decomposition of the velocity space into half-spaces, Ghiroldi and Gibelli [21] included discontinuous polynomial basis functions into their Galerkin scheme.
This work extends the ideas from [18] to the spatially inhomogeneous Boltzmann equation. Approaches similar to our method can be found in [26] and [16,23] for two dimensions, both in velocity and space.
In [16], Fonn, Grohs, and Hiptmair proposed a spectral scheme for the twodimensional spatially homogeneous Boltzmann equation. The particle density function is sought as a linear combination of basis functions consisting of polynomials, given in polar coordinates, multiplied with a Maxwellian weight. The homogeneous Boltzmann equation is discretised by a Galerkin ansatz with respect to a weighted inner product on L 2 (R 3 ). The physical quantities related to the collision invariants are not conserved in general but only for a special choice of the weight function. The numerical effort of computing the discretised collision operator is reduced by exploiting the simple representation of spherical harmonics and their rotations in two dimensions, namely Y l (ϕ) = exp(ilϕ). In [23], this ansatz was extended by Grohs, Hiptmair, and Pintarelli to the two-dimensional inhomogeneous Boltzmann equation. Strang splitting is used to decouple the advection and collision step. In order to fulfil the conservation laws, two methods are presented, one by adding a basis of the collision invariants to the test functions, which leads to an overdetermined system of ODEs. The other possibility is adding a Lagrange multiplier which corrects the coefficients of the basis functions at each time step. The advection equation is solved by least squares method on the phase space combined with a continuous first order FEM.
Kitzler and Schöberl [26] use a Galerkin-Petrov method in the velocity space and a Discontinuous Galerkin discretisation of the advection part for the twodimensional inhomogeneous Boltzmann equation. The test functions are spanned by tensor products of (continuous) polynomials in the velocity variable and piecewise continuous polynomials in the space variables. The basis functions are obtained by multiplying the test functions with a Maxwellian weight, where the bulk velocity and temperature may depend on the spatial variable. This automatically insures exact conservation laws. The translation invariance together with the simple parametrisation of S 1 is used to reduce the numerical effort of the collision step.
This paper is organised as follows. In Section 2, we introduce our notation for the Boltzmann equation and give a formulation of the one-dimensional problems which are discussed in this paper. In Section 3, we investigate the semi-discretised Boltzmann equation which arise after an integration over the velocity space. This also includes a small section on the Burnett functions which are used for the construction of our basis and test functions. Section 4 gives some details of the implementation, in particular numerical integration over R 3 and half-spaces as well as the formulation of boundary conditions and the discretisation of the hyperbolic system. In Section 5, we present the results of our scheme for two standard tests and different Knudsen numbers. Some concluding remarks can be found in Section 6.

2.
One-dimensional problems for the Boltzmann equation. We consider the spatially one-dimensional Boltzmann equation which models the time evolution of the particle density function of a gas located in (0, 1). Suitable initial and boundary conditions will be discussed later.
The right-hand side of equation (1) is called the collision integral or collision operator. It is given by and σ denotes the differential cross-section which depends on the modelled gas interaction. Kn is called Knudsen number which is defined as the ratio of the mean free path λ MP and a characteristic length L of the system The collision kernel |u|σ(|u|, |u| −1 (u, e)) can usually be expressed as for some constants C λ , −3 < λ ≤ 1.
The weak formulation of the collision operator is particularly useful for our numerical scheme as in this form, the post-collision velocities v and w are transferred to the test function ψ.
As a direct consequence of Theorem 2.1 we get the following lemma.
Testing the Boltzmann equation with five functions, called collision invariants, we conclude the conservation of mass, momentum in all three directions and energy.
Proof. It is not hard to see that ψ 0 , . . . , ψ 4 fulfil The claim follows after an application of Theorem 2.1.
In order to get a well-defined problem for equation (1), initial and boundary conditions are imposed on the particle density function f . At the beginning, the particle density function is given by Maxwell distributions The boundaries can be opened, e.g. the shock tube, or closed in the case of the Fourier problem. For open boundaries, the particle density function has to fulfil inflow conditions If the boundary conditions should model walls, i.e. the boundaries are closed, diffusely reflecting conditions are imposed on f , : (w, n(x)) < 0} and n(0) = (1, 0, 0) and n(1) = (−1, 0, 0) denote the inner normals. This choice conserves the total mass of the gas.
3. Galerkin-Petrov method. The crucial part when solving the Boltzmann equation is the application of the collision operator. Every evaluation involves the computation of a five-fold integral. Furthermore, the operator is non-linear. However, the collision operator only acts on the velocity space and by the use of its weak form, we can exploit its bilinear structure. In this formulation, the collision operator is descretised by a sequence of symmetric matrices which are independent of time and space as well as chosen initial and boundary conditions. Once these matrices are computed and stored, they can be used for different initial and boundary values as well as for different time-space integration schemes.
This motivates the following ansatz for the approximation f (n) of the particle density function where Φ = ϕ 1 , . . . , ϕ n are called basis or ansatz functions, and is the vector of coefficients. Plugging this into the Boltzmann equation (1), we get the residual for t > 0, x ∈ (0, 1) and v ∈ R 3 . The coefficients (f j ) n j=1 are chosen such that for all t > 0 and x ∈ (0, 1), where the functions (ψ i ) n i=1 are called test functions. If ψ i = ϕ i for i = 1, . . . , n, this is the classical Galerkin method. However, choosing a set of test functions which includes the five collision invariants, the corresponding collision matrices vanish. We automatically conserve mass, momentum and energy.

TORSTEN KEßLER AND SERGEJ RJASANOW
For our method, the test functions are polynomials, including the collision invariants. The ansatz space is spanned by test functions multiplied with an Maxwellian weight. This choice of basis functions leads to consistent results in case of the steady-state and in the hydrodynamic limit, Kn → 0.
3.1. Burnett functions. Our basis and test functions are derived from the Burnett functions, which were introduced in kinetic theory by Burnett in 1935, see [7]. They were already used in [18] to construct basis and test functions. We adopt their main ideas with slight modifications regarding scaling factors.
Definition 3.1. Given k ∈ N, l ∈ N and m ∈ {−l, . . . , l}, the (normalised) Burnett function ψ k,l,m is defined in spherical coordinates, v = ρe, as where µ k,l is a normalisation constant, , L (l+ 1 2 ) k denotes the associated Laguerre polynomial of degree k and Y l,m the realvalued spherical harmonic, given by Here, Y m l denotes the usual complex-valued spherical harmonic.

Remark 1.
For v ∈ R 3 \ {0}, ρ ∈ [0, ∞) and e ∈ S 2 with v = ρe are uniquely given by If v = 0, ρ is equal to 0, but e can be chosen arbitrarily. However, the Burnett functions are well-defined. This can be seen by a case distinction of l ∈ N.
In the case of l = 0, Y 0,0 is a constant polynomial. Thus the value of ψ k,0,0 does not depend on e for all k ∈ N.
If l > 0, we have ρ l = 0 l = 0 and ψ k,l,m (0) = 0 for all k ∈ N and m ∈ {−l, . . . , l} independent of e. Proof. Let k ∈ N, l ∈ N and m ∈ {−l, . . . , l}. Then is a polynomial. The spherical harmonic Y l,m is the restriction of a homogeneous polynomial P l,m of degree l on S 2 . Hence it holds for all v = ρe ∈ R 3 ρ l Y l,m (e) = P l,m (v), which completes the proof.
Lemma 3.4. The collision invariants can be expressed as a linear combination of Burnett functions, namely

Basis and test functions.
Starting with the Burnett functions, we construct our basis and test functions. In order to approximate a broader range of temperatures, we rescale the Burnett functions by a characteristic temperature T > 0.
the test function with indices k, l, m.
the basis function with indices k, l, m.
Lemma 3.7. For all k, k ∈ N, l, l ∈ N and m ∈ {−l, . . . , l}, m ∈ {−l , . . . , l }, it holds If T = 1, we drop the superscript and simply write ϕ k,l,m and ψ k,l,m instead of ϕ T k,l,m and ψ T k,l,m . In the following, we will often write ϕ T j and ψ T i , meaning that we implicitly choose an enumeration of the basis and test functions.
3.3. Semi-discretised Boltzmann equation. The set of indices used for the basis and also for the test functions are chosen from a sequence of increasing sets which are indexed by two integers, indicating the highest used degree for radial and spherical polynomials, that is for K, L ∈ N the indices are given by This insures that for all non-trivial choices of K and L, the collision invariants are always elements of the set of test functions.
Plugging our basis and test functions into equation (6), we get a system of n = (K + 1)(L + 1) 2 partial differential equations. Rewritten in matrix form, these equations read as where M T is the mass matrix with entries D T is the so-called flow matrix whose entries are given by . . , n. We drop the superscripts in case of T = 1.
It will be shown in the following that the system of PDEs in equation (9) can be expressed as exposing the dependence of the matrices on the characteristic temperature T .
Mass matrix. The entries of the mass matrix M T can be rewritten as Thus the mass matrix is independent of T . The last line in the earlier equation shows that the mass matrix is the Gramian of a weighted L 2 -inner product on R 3 . Hence, it is symmetric and positive definite.
Flow matrix. The entries of the flow matrix D T can be simplified to It can be seen from the last line of the above equation that D T is symmetric.
Collision matrices. For the computation of the collision matrices, the weak form of the collision operator, as presented in Theorem 2.1, is used. Thus, for i, k, l = 1, . . . , n, Q T i [k, l] is given by where with v, w ∈ R 3 . As before, we scale the velocities with √ T by applying the substitution which has Jacobian T 3 , on the integral in equation (11). The collision velocities transform as for all v, w ∈ R 3 , e ∈ S 2 and therefore Combining this with equation (11) yields The substitution Consequently, the collision matrices are symmetric.
Lemma 3.8. There exist a basis of R n , denoted by (r i ) n i=1 , of generalised eigenvectors of D with respect to M and generalised eigenvalues ( The generalised spectrum of D is point symmetric with respect to 0, i.e. if λ is a generalised eigenvalue of D, then −λ is a generalised eigenvalue of D, too.
Proof. As D is symmetric and M is positive definite, the existence of the generalised eigensystem follows from standard theory. Let r ∈ R n be a generalised eigenvector of D with respect to M and λ ∈ R be the corresponding generalised eigenvalue. We definer ∈ R n byr The generalised spectrum of D with respect to M in case of K = 9, L = 9, i.e. n = 1000, is shown in Figure 1.
We can now decouple the left-hand side. Following the notation of Lemma 3.8, let R denote the matrix whose columns are the generalized eigenvectors (r i ) n i=1 and Λ = diag(λ 1 , . . . , λ n ). Equation (12) reads as Substituting f = Rg in equation (10), we get, after an application of equation (13), Initial and boundary conditions of the coefficient vector are obtained by minimising the L 2 -error which occurs when plugging the approximation f (n) into equations (3) and (4) or (5).
As the basis functions are orthonormal, we have .
Let us assume that there is an approximation of the outgoing flux σ out at the boundaries in case of diffuse reflection. If we impose inflow conditions, σ out = 1.
As the boundary conditions are only imposed on R 3 in , we minimise the L 2 -error on this half-space. Therefore, for all t > 0 and boundary points x, f (t, x) has to fulfil the linear system meaning that that this matrix is independent of T . M in is symmetric and positive definite. The approximation of the outgoing fluxes is discussed in Section 4.
4. Numerical realisation. The computation of the matrices from Section 3 involves three types of integrals, namely where H is a half-space, given by H = {v ∈ R 3 : (v, n) > 0} with n ∈ S 2 , and p is polynomial in three variables.

TORSTEN KEßLER AND SERGEJ RJASANOW
As the Burnett functions are given in spherical coordinates, the integrals over R 3 and H are rewritten as (16) The radial integrals are approximated by a Gauß quadrature generated by the weight function This quadrature rule is sometimes referred as Gauß-Maxwell quadrature [33]. The integrals over the sphere are approximated by the Lebedev quadrature [27]. Compared to a tensorised trapezoidal rule, the Lebedev quadrature is more efficient in terms of function evaluations in order to integrate polynomials up to a fixed degree exactly.
Although the Lebedev quadrature is suitable for the integration over S 2 , it fails if one integrates over half-spheres as demanded in (16) due to the discontinuity of the corresponding indicator function. In this case, we express p as a linear combination of spherical harmonics. The integral over the sphere is now reduced to integrals such as and β ∈ {0, 1}, which can be computed with the Funk-Hecke theorem [9]. The case β = 0 refers to projections an half-spaces, whereas integrals with β = 1 are needed for the computation of outflows.
Theorem 4.1 (Funk-Hecke). Let Y l be a spherical harmonic of degree l ∈ N, and g ∈ L 1 ([−1, 1]). For all n ∈ S 2 it holds: where P l is the Legendre polynomial of degree l.
An application of Theorem 4.1 yields for the earlier equation The one-dimensional integrals can be computed analytically, denotes the generalised binomial coefficient. By the use of these quadratures, one can compute the matrices from Section 3.3. Most of the time is spent on the computation of the collision matrices. That's why we will address their computation in further detail.
Computation of the collision matrices. In this section, we follow the approach from [18]. In contrast to [18], we use a Gauß-Maxwell quadrature for the computation of the radial integrals. Let us fix K ∈ N, the maximal degree of Laguerre polynomials, and L ∈ N, the maximal degree of spherical harmonics. We choose the indices of our basis and test functions as  (11) yields where v iv,jv = ρ iv e jv , w iw,jw = ρ iw e jw and (v iv,jv − w iw,jw , e j ) |v iv,jv − w iw,jw | × ψ i (v iv,jv,iw,jw (e j )) + ψ i (w iv,jv,iw,jw (e j )) − ψ i (v iv,jv ) − ψ i (w iw,jw )) ,
However, this form is not used for the computation of the matrix entries. Due to the product structure of the Burnett functions, there are some repetitive parts which have to be computed only once. Storing them reduces the computational effort.
Let us define the three arrays P GM , P L and P Q . The first array P GM is given by The approximation of the matrix entry of the collision matrix is now given by (P GM ) kv,lv,iv (P L ) lv,mv,jv (P GM ) kw,lw,iw (P L ) lw,mw,jw (P Q ) iv,jv,iw,jw,i The numerical effort for computing P GM and P L is sublinear in the number of basis functions and linear in the number of quadrature nodes. The computation time is negligible to the time spend on computing P Q . This is the crucial part of the numerical integration as its computation involves N 2 GM N 2 L evaluations of the spherical integrals, leading to N 2 GM N 3 L evaluations of the integrand which itself is rather expensive due to the computation of the collision velocities and the transformation from Cartesian to spherical coordinates in order to evaluate the test function ψ i . However, the computation of the entries of P Q can readily be parallelised. The same holds true if one computes the entries of Q i using the three precomputed arrays P GM , P L and P Q .
Numerical solution of the PDE system. For the solution of the PDE system in (14), we follow a finite volume approach. Each component of g is approximated a function which is constant on every cell of the spatial discretisation as depicted in Figure 2. The outgoing flow at the walls is approximated by for the left wall, Computational complexity. Although no splitting is performed in our scheme, we will discuss the computational complexity separately for the free flow and the collision phase. As before, let n denote the number of basis and test functions in velocity cells and let N x be the number of spatial cells. First, at each time step, appropriate boundary conditions are computed. This implies solving a linear system of dimension n × n. As only the right-hand side changes during the computation, but not the mass matrix of the corresponding half-space, this matrix can be precomputed for all boundary cells and can be stored in a factorised form, for instance using the Cholesky decomposition. This leads to a complexity of O(n 2 ) for each boundary cell.
For our one-dimensional problems, the number of boundary cells is 2, independent of the number of spatial cells. For two-or three-dimensional problems, the number of boundary cells depend on the number of spatial cells and so does the computation of boundary values.
The generalised eigenproblem in equation (12) is solved once at the beginning of the computation. Therefore computing g from f or vice versa is of complexity O(n 2 ). For the local flux reconstruction for the n uncoupled transport equations in (14), we get a complexity estimate of O(N x n). For the whole free flow phase we obtain a computational complexity of O(N x n 2 ). Evaluating the discretised collision operator evolves n bilinear forms, each of them having a computational complexity of O(n 2 ). This leads to a complexity of O(N x n 3 ) for the evaluation of the discretised collision operator in each spatial cell. The evaluation of the collision operator dominates the computation time, even for a small number of degrees of freedom.
One possibility to reduce the computational time without reducing the asymptotic complexity is the reformulation of the collision step, using the fact that the collision operator is defined locally and can therefore be evaluated for all spatial cells simultaneously. For the spatial cell with midpoint x j+ 1 2 , let f j denote the vector of coefficients of the basis function for this cell. A direct evaluation is done described in Algorithm 1, where C j ∈ R n denotes the discretised collision operator on cell j, see equation (14). This algorithm only performs matrix-vector and vector-vector operations Although it can easily be parallelised, it cannot fully benefit from modern CPU architectures. One way to reformulate parts of the algorithm one perform matrix-matrix operation instead of several matrix-vector operations.

TORSTEN KEßLER AND SERGEJ RJASANOW
Let F ∈ R n×Nx denote the matrix whose j-th column is given by the vector f j . The algorithm can be now reformulated with the help of F , see Algorithm 2. The new algorithm can now make use of the highly tuned matrix-matrix operations available through packages like BLAS in its different implementations. To validate our computations, we compare them with a stochastic simulation, carried out with the Direct Simulation Monte Carlo (DSMC).

Fourier problem.
We choose diffusely reflecting walls at the left and the right boundary, see (5). We seek for the steady state solution. Initially, the gas is in equilibrium state. The parameters of the Maxwellians in equation (3) are constant, ρ 0 = 1, V 0 = (0, 0, 0) and T 0 = 1.25. The left wall has temperature T l = 1.0, the right wall T r = 1.5. The characteristic temperature of the basis and test functions T is set to 1.25. This situation is depicted in Figure 3. We compare the results for three different Knudsen numbers, Kn = 1.0, 0.25 and 0.025. The final time t f is chosen such that on average, a particle has crossed the interval (0, 1) ten times. It can be calculated as t f ≈ 6.32. The number of cells N x for both the stochastic and the deterministic method is 512. The DSMC solutions were obtained by 20 000 repetitions and 1024 particles per cell in the beginning.
The time step size τ of the deterministic method is chosen as where CFL is the Courant number and (λ i ) n i=1 are the generalised eigenvalues of D with respect to M . In our tests, CFL was set to 0.5. The fluxes are approximated Figure 3. Sketch of the one-dimensional Fourier problem. We seek the particle density function along the axis labeled by x. T l , T r are the temperatures of the walls, T 0 is the initial temperature of the gas. by the fifth order WENO5 scheme [25]. The resulting system of ODEs is integrated by an explicit third order strong stability-preserving Runge-Kutta method [22]. The final densities and temperatures for the stochastic method and several choices of K and L, c.f. Table 1, are shown in Figure 7 for Kn = 1.0, in Figure 9 for Kn = 0.25 and in Figure 11 for Kn = 0.025, respectively. They are also plotted near the walls in Figures 8, 10 and 12 for the Knudsen numbers 1.0, 0.25 and 0.025, respectively. Plots of the final particle densities restricted to a straight line for the Knudsen numbers 1.0, 0.25 and 0.025 can be found in Figures 13a, 13b and 13c, respectively. In Figure 6, a contour plot of the particle density function for Kn = 0.25 restricted to the plane {v ∈ R 3 : v 3 = 0} is shown.
The lower the Knudsen number, the better our method performs compared to the stochastic solution. For Knudsen numbers above 1.0, free flow is dominating over the collision interaction. This leads to almost discontinuous solutions of the Boltzmann equation, especially close to the walls, where only the inflowing particles are affected by the boundary conditions. Due to the lack of sufficient interaction for Kn = 1.0, these streams of different temperatures are not well mixed. Therefore the particle density functions have different asymptotic behaviour towards ±∞. These functions cannot be approximated well with our basis functions, see Figure 13a. That's why the highest errors in Figures 7a and 7b appear near the walls.
Lowering the Knudsen number leads to a higher importance of the collisions. The collision operator smoothes the particle density function. Moreover, the higher impact of the collisions also reduces the distance in which the particles preserve their temperatures given by the boundary conditions. The different temperatures of particles emerging from the walls and the rest of the gas are mixed. This leads to more symmetric particle density functions which can be approximated well by our basis function, c.f. Figure 13c.
Therefore, our method shows a very good agreement with the stochastic simulation in case of Kn = 0.025.

5.2.
Conservation of total mass. In its semi-discrete form, see equations (10) and (14), our scheme conserves mass, momentum and energy. However, in its full-discretised formulation, several errors are introduced, for instance the approximation error of initial and boundary particle density functions in velocity space, or the discretisation error in time and space. These errors may also be coupled.
In case of diffusely reflecting walls, the outflows at (intermediate) time steps are approximated by the outflow of the particle density functions next to the walls from the last (intermediate) time step. We will inspect the error in mass conservation for different discretisation parameters and Knudsen numbers.
We consider a gas in the interval [0, 1] with two diffusely reflecting walls with temperature T w = 1. Depending on the choice of discretisation parameters, Maxwellians are represented exactly. For a non-zero initial L 2 -error, we choose an initial condition far from equilibrium, namely a mixture of two Maxwellians, see equation (3), In the following, we use α = 1 2 , ρ 1 = ρ 2 = 1, In particular, we have a constant density equals 1 and therefore a total mass M 0 = 1. The relative L 2 -error for different degrees of freedom and characteristic temperature T = 1 can be found in Table 2. We observe the expected spectral convergence.
With a fixed end time t f = 5 and CFL = 0.5, we run our simulation with different Knudsen numbers, which affect the smoothness of the solution, and different choices of basis functions and number of spatial cells. At each time step, the total mass is computed by summation of the local densities, Note, that the midpoint rule is exact in this case, as the coefficients and therefore the macroscopic quantities are piecewise constant. Figure 4 shows the typical course of the total mass. At the beginning, the total mass is larger than 1 because of the error of the L 2 -projection. It then decreases until the gas is sufficiently close to the equilibrium distribution. In this case, the outflow relaxes too and becomes (nearly) constant. From this point, the total mass is perfectly conserved, albeit with a non-vanishing error compared to the analytic result.
The Tables 3, 4, 5 show the relative L ∞ -error of the total mass for different discretisation parameters, for the Knudsen numbers Kn = 1, Kn = 0.1 and Kn = 0.01, respectively. For Kn = 1, the spectral accuracy is lost. In this regime the particle density functions shows discontinuities in low order derivatives. For a fixed number of spatial cells, only a decrease of linear order in mass conservation can be observed. Furthermore, the approximation error in velocity space dominates the error in space or time.
Lowering the Knudsen number leads to a smoother particle density function. For the choices K = 1, L = 2 and K = 2, L = 4 and K = 3, L = 6, the L 2 -error of the initial condition dominates, compare Table 2. For K = 4, L = 8, the smoothness of the solution of the Boltzmann equation limits the error, nearly independent of the spatial discretisation.
In case of Kn = 0.01, the solution of the Boltzmann equations seems to be smooth enough, such that the total mass is conserved up to the spectral error introduced by the approximation of the initial condition. 5.3. Shock tube problem. As a third example, we consider a Sod-like shock tube problem [34]. The parameters of the initial Maxwellians are piecewise constant, see Figure 5, We show the results for three different Knudsen numbers, Kn = 0.1, 0.01, 0.001 and Table 2. Relative L 2 -error for the mixture from equation (18) with parameters given in equation (19). For the DSMC solution, we chose 512 cells with a typical number of 1024 particles per cell. For the Galerkin-Petrov method, the index set for our basis functions is given by I 3,3 . For the characteristic temperature we choose T = 1. Increasing the radial or the spherical degree hardly changes the results. We use 512 cells in space. The system of ODEs is integrated with a adaptive Runge-Kutta-Fehlberg Figure 5. Sketch of the initial situation of the shock tube problem. Two areas of same bulk velocities and temperatures but different densities are separated by diaphragm (dashed line), which is removed at t = 0.
method [12] to address the stiffness which is caused by the small Knudsen numbers in this example. Before discussing the particle density functions, we first compare the computed macroscopic quantities. Density, bulk velocity, temperature and pressure at the end time t f are shown in Figures 14, 15 and 16 for the Knudsen numbers 0.1, 0.01 and 0.001, respectively.
The physical quantities from the Galerkin-Petrov method are in good agreement with those from the DSMC method for all three Knudsen numbers.
Comparing the results from the Boltzmann equation with the exact solution of the Euler equations, one can clearly see the differences between them become smaller when lowering the Knudsen number, which is conform to the Hilbert expansion. The best agreement is observed in the rarefaction wave and shock discontinuity. Due to a non-vanishing physical (and numerical) viscosity, the edges of the solutions of the Boltzmann equation are smeared out.
However near the contact discontinuity, the solutions of the Boltzmann equation and the Euler equations differ even for small Knudsen numbers. This can be seen best by studying the bump in the bulk velocity, which is located at the contact discontinuity. Generated at t = 0 by the removal of the diaphragm, it travels with the speed of the contact discontinuity, which seems to be independent of the Knudsen number.
Our numerical experiments indicate that, in contrast to the solution of the Euler equations, the bulk velocity is not constant around the contact discontinuity, but slightly lower left from the bump and larger at the right. This means that the two discontinuities are travelling with different speeds.
As already mentioned in the discussion of the macroscopic quantities, three parts of the shock dynamic are of interest, namely the rarefaction wave, the shock and the contact discontinuity. A good agreement on the macroscopic scale does not necessarily imply small errors in the particle density function, see the section on the Fourier problem.
We fix three spatial positions, x = 0.4, 0.5 and 0.6 and plot the marginal particle density functions when the gradient of the density at one of these points is steepest, see In case of Kn = 0.1, the particle density functions obtained from DSMC are discontinuous ( Figure 18 for x = 0.6) or have discontinuous derivatives ( Figure  17 at x = 0.5). Therefore, spectral convergence is lost and adding basis and test functions of higher degree will not hardly decrease the error.
With Kn = 0.01, the approximation becomes better, although in the shock front, the particle density function still has discontinuous derivatives in the beginning, see Figure 20 at x = 0.5.
Comparing the particle density functions for Kn = 0.001, the descrisation error of the transport operator now dominates the error of the particle density functions. Especially at the jumps of the macroscopic quantities, the error both introduced by the WENO reconstruction for the Galerkin-Petrov method and the Lagrangian method for DSMC is larger than the approximation error for smooth functions in the velocity space.
6. Conclusion. The Galerkin-Petrov method presented in this paper is fully conservative independent of the time stepping method. Most of the computation time is spent on the collision matrices. However, their computation can be easily parallelised. Furthermore, having computed them for a specific choice of basis and test functions, they can be reused for different problems, as they are independent of the initial and boundary conditions, the characteristic temperature or the spatial and temporal discretisation. This leads to a very fast deterministic method, as a computation of the collision operator is reduced to the evaluation of several quadratic forms. Depending on the problem, one can obtain very good results with around 100 basis functions.      (17). Figures 10a and 10b show the density at the time t f near the left and the right wall, respectively. Figures 10c and  10d show the temperature at the time t f near the left and the right wall, respectively.   (17). Figure 11a shows the density at the time t f , whereas in Figure 11b, the temperature at the time t f is shown.  (17). Figures 12a and 12b show the density at the time t f near the left and the right wall, respectively. Figures 12c and  12d show the temperature at the time t f near the left and the right wall, respectively. x x = 0.6 K = 3, L = 3 DSMC Figure 17. Particle density functions for Kn = 0.1 shortly after the diaphragm is removed. x x = 0.6 K = 3, L = 3 DSMC Figure 18. Particle density functions for Kn = 0.1 when the shock discontinuity reaches the third evaluation point. x x = 0.6 K = 3, L = 3 DSMC Figure 19. Particle density functions for Kn = 0.1 when the contact discontinuity reaches the third evaluation point at time t f . x x x = 0.6 K = 3, L = 3 DSMC Figure 21. Particle density functions for Kn = 0.01 when the shock discontinuity reaches the third evaluation point. x x = 0.6 K = 3, L = 3 DSMC x x = 0.6 K = 3, L = 3 DSMC x x = 0.5 x x = 0.6 K = 3, L = 3 DSMC Figure 24. Particle density functions for Kn = 0.001 when the shock discontinuity reaches the third evaluation point. x x = 0.6 K = 3, L = 3 DSMC Figure 25. Particle density functions for Kn = 0.001 when the contact discontinuity reaches the third evaluation point at time t f .