Giacomo Dimarco
email: giacomo.dimarco@unife.it

Raphaël Loubère

Jacek Narski

Thomas Rey

An efficient numerical method for solving the Boltzmann equation in multidimensions

Keywords: Boltzmann equation, kinetic equations, semi-Lagrangian schemes, spectral schemes, 3D/3D, GPU, CUDA, OpenMP, MPI. 2000 MSC: (82B40, 76P05, 65M70, 65M08, 65Y05, 65Y20

In this paper we deal with the extension of the Fast Kinetic Scheme (FKS) [J. Comput. Phys., Vol. 255, 2013, pp 680-698] originally constructed for solving the BGK equation, to the more challenging case of the Boltzmann equation. The scheme combines a robust and fast method for treating the transport part based on an innovative Lagrangian technique supplemented with fast spectral schemes to treat the collisional operator by means of an operator splitting approach. This approach along with several implementation features related to the parallelization of the algorithm permits to construct an efficient simulation tool which is numerically tested against exact and reference solutions on classical problems arising in rarefied gas dynamic. We present results up to the 3D×3D case for unsteady flows for the Variable Hard Sphere model which may serve as benchmark for future comparisons between different numerical methods for solving the multidimensional Boltzmann equation. For this reason, we also provide for each problem studied details on the computational cost and memory consumption as well as comparisons with the BGK model or the limit model of compressible Euler equations.

Introduction

Kinetic equations consider a representation of a gas as particles undergoing instantaneous collisions interspersed with ballistic motion [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF][START_REF] Dimarco | Numerical methods for kinetic equations[END_REF]. Nowadays, these models appear in a variety of sciences and applications such as astrophysics, aerospace and nuclear engineering, semiconductors, plasmas related to fusion processes as well as biology, medicine or social sciences. The common structure of such equations consists in a combination of a linear transport term with one or more interaction terms which furnishes the time evolution of the distribution of particles in the phase space. For its nature, the unknown distribution lives in a seven dimensional space, three dimensions for the physical space and three dimensions for the velocity space, plus the time. This makes the problem a real challenge from the numerical point of view, since the computational cost becomes immediately prohibitive for realistic multidimensional problems [START_REF] Dimarco | Numerical methods for kinetic equations[END_REF]. Aside from the curse of dimensionality problem, there are many other difficulties which are specific to kinetic equations. We recall two among the most important ones. The computational cost related to the evaluation of the collision operator involving multidimensional integrals which should be solved in each point of the physical space [START_REF] Filbet | Solving the Boltzmann equation in N log2 N[END_REF][START_REF] Pareschi | Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator[END_REF]. The second challenge is represented by the presence of multiple scales which requires the development of adapted numerical schemes to avoid the resolution of the stiff dynamics [START_REF] Dimarco | High order asymptotic-preserving schemes for the Boltzmann equation[END_REF][START_REF] Dimarco | Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations[END_REF][START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF][START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF][START_REF] Jin | Efficient asymptotic-preserving (ap) schemes for some multiscale kinetic equations[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Degond | Asymptotic-preserving schemes for fluid models of plasmas[END_REF] typically arising when dealing with problems with multiple regimes.

Historically, there exists two different approaches which are generally used to tackle kinetic equations from a numerical point of view: deterministic numerical methods such as finite volume, semi-Lagrangian and spectral schemes [START_REF] Dimarco | Numerical methods for kinetic equations[END_REF], and, probabilistic numerical methods such as Direct Simulation Monte Carlo (DSMC) schemes [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF][START_REF] Caflisch | Monte carlo and quasi-monte carlo methods[END_REF]. Both methodologies have strengths and weaknesses. While the first could normally reach high order of accuracy, the second are often faster, especially for solving steady problems but, typically, exhibit lower convergence rate and difficulties in describing non stationary and slow motion flows.

In this work, we deal with deterministic techniques. In particular, we focus on semi-Lagrangian approaches [START_REF] Crouseilles | A forward semi-lagrangian method for the numerical solution of the vlasov equation[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Güçlü | A high order cell-centered semi-Lagrangian scheme for multi-dimensional kinetic simulations of neutral gas flows[END_REF][START_REF] Shoucri | Numerical integration of the Vlasov equation[END_REF][START_REF] Filbet | Accurate numerical methods for the Boltzmann equation[END_REF][START_REF] Filbet | High order numerical methods for the space non-homogeneous Boltzmann equation[END_REF] for the transport part coupled with spectral methods [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF] for the interaction part. Our main goal is to tackle the challenges related to the high dimensionality of the equations and with the difficulties related to the approximation of the collision integral. More in details, we deal with the extension of our recent works [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF][START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF][START_REF] Dimarco | A multiscale fast semi-Lagrangian method for rarefied gas dynamics[END_REF][START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part III: High-performance-computing[END_REF] which were based on an innovative semi-Lagrangian technique for discretizing the transport part of a kinetic model (FKS method). This technique has been applied solely to the solution of a simple kinetic equation with a relaxation type collision operator, i.e. the BGK (Bhatnagar-Gross-Krook) operator [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF]. Here we extend it to the case of the full Boltzmann operator [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF][START_REF] Cercignani | The Boltzmann equation and its applications[END_REF] and we test his performances up to the six dimensional case. The FKS is based on the classical discrete velocity models (DVM) approach [START_REF] Bobylev | On approximation of the Boltzmann equation by discrete velocity models[END_REF][START_REF] Palczewski | A consistency result for a discrete-velocity model of the Boltzmann equation[END_REF][START_REF] Palczewski | Existence, stability, and convergence of solutions of discrete velocity models to the boltzmann equation[END_REF][START_REF] Mieussens | Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics[END_REF]. Successively, in order to overcome the problem of the excessive computational cost, we propose to use a Lagrangian technique which exactly solves the transport step on the entire domain, without reprojecting the solution on the grid at each time step. The FKS approach was shown to be an efficient way to solve linear transport equations, and, it has permitted the simulation of full six dimensions problems on a single processor machine [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF]. Unfortunately the solutions computed with this method are limited to a first order in space and time precision. Extension of the method to high order reconstruction is under consideration [START_REF] Dimarco | A new class of high order semi-Lagrangian schemes for rarefied gas dynamics[END_REF]. Concerning the discretization of the collision operator we rely on Fourier techniques. For the resolution of the Boltzmann integral, these techniques have been first introduced independently by L. Pareschi and B. Perthame in [START_REF] Pareschi | A fourier spectral method for homogeneous Boltzmann equations[END_REF] and by A. Bobylev and S. Rjasanow in [START_REF] Bobylev | Difference scheme for the Boltzmann equation based on the fast Fourier transform[END_REF]. Since then, this approach has been investigated by a many authors [START_REF] Filbet | Solving the Boltzmann equation in N log2 N[END_REF][START_REF] Pareschi | Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator[END_REF][START_REF] Gamba | Shock and boundary structure formation by spectral-Lagrangian methods for the inhomogeneous Boltzmann transport equation[END_REF][START_REF] Gamba | Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states[END_REF][START_REF] Gamba | A conservative spectral method for the Boltzmann equation with anisotropic scattering and the grazing collisions limit[END_REF][START_REF] Wu | Deterministic numerical solutions of the Boltzmann equation using the fast spectral method[END_REF][START_REF] Bobylev | Difference scheme for the Boltzmann equation based on the fast Fourier transform[END_REF][START_REF] Bobylev | Numerical solution of the Boltzmann equation using a fully conservative difference scheme based on the fast Fourier transform[END_REF][START_REF] Pareschi | Spectral methods for the non cut-off Boltzmann equation and numerical grazing collision limit[END_REF][START_REF] Filbet | Analysis of spectral methods for the homogeneous Boltzmann Equation[END_REF][START_REF] Wu | Influence of intermolecular potentials on rarefied gas flows: Fast spectral solutions of the Boltzmann equation[END_REF][START_REF] Alekssenko | Deterministic solution of the Boltzmann equation using a discontinuous Galerkin velocity discretization[END_REF][START_REF] Pareschi | On the stability of spectral methods for the homogeneous Boltzmann equation[END_REF][START_REF] Wu | A fast spectral method for the Boltzmann equation for monatomic gas mixtures[END_REF]. In this work we will make use in particular of the fast method described in [START_REF] Mouhot | Fast algorithms for computing the Boltzmann collision operator[END_REF][START_REF] Filbet | Solving the Boltzmann equation in N log2 N[END_REF] which has a complexity of the order of O(N dv log(N dv )) where N is the number of point in which the velocity space is discretized in one direction and d v the dimension of the velocity space. The method preserves mass, and approximates with spectral accuracy momentum and energy.

Combining opportunely the FKS method with the fast spectral approach we have developed a method for solving the Boltzmann equation up to the six dimensional case for unsteady flows. In order to obtain such result we have constructed a parallel version of our algorithm taking advantages of Graphical-Processor-Unit (GPU) under CUDA language. The results presented in this work show that we are nowadays ready and able to use kinetic equations to simulate realistic multidimensional flows. Up to our knowledge this is one of the first examples in literature of solution of the full multidimensional Boltzmann equation by means of deterministic schemes. Our main limitation to run extensive simulations remains at the present moment the lack of memory capacity due to the use of shared memory machines. However, in the present paper, we also report some preliminary results about the extension of our method to deal with distributed-memory computers but we postpone for a future research the detailed analysis of the Message Passive Interface (MPI) version of the algorithm as well as the analysis of its performances.

The article is organized as follows. In Section 2 we recall the Boltzmann equation. In Section 3 we present the Fast Kinetic Scheme and the fast spectral scheme. In Section 4 we detail the aspects related to the implementation of the resulting algorithm necessary to realize an efficient parallel numerical tool. Several tests starting from the 0D× 2D up to the 3D×3D case are studied in details in Section 5. These tests assess the validity of our approach as well as detail all the computational aspects. Conclusions and future works are finally exposed in Section 6.

The Boltzmann equation

In this section we briefly recall the Boltzmann equation and its main properties, we also recall some related models, i.e. the BGK model and the compressible Euler equations, which will be used for numerical comparisons. We refer the reader to [START_REF] Dimarco | Numerical methods for kinetic equations[END_REF][START_REF] Cercignani | The Boltzmann equation and its applications[END_REF] for an exhaustive description.

In the kinetic theory of rarefied gases, the non-negative function f (x, v, t) characterizes the state of the system and it defines the density of particles having velocity v ∈ R dv in position x ∈ R dx at time t ∈ R + , where d x is the physical dimension and d v the dimension of the velocity space. The time evolution of the particle system is obtained through the equation

∂f ∂t + v • ∇ x f = Q(f ). (1)
The operator Q(f ), on the right hand side in equation [START_REF] Cercignani | The Boltzmann equation and its applications[END_REF], describes the effects of particle interactions and its form depends on the details of the microscopic dynamic. Independently on the type of microscopic interactions considered, typically the operator is characterized by some conservation properties of the physical system. They are written as

R dv Q(f )φ(v) dv = 0, (2)
where φ(v) = (1, v, |v| 2 ) are commonly called the collision invariants. We denote by

U (x, t) = R dv f (x, v, t)φ(v) dv ∈ R 2+dv
the first three moments of the distribution function f , namely, U (x, t) = (ρ, ρu, E), where ρ is the density, ρu the momentum and E the energy. Integrating (1) against φ(v) yields a system of macroscopic conservation laws

∂ ∂t R dv f φ(v) dv + R dv v • ∇ x f φ(v) dv = 0. ( 3
)
The above moment system is not closed since the second term involves higher order moments of the distribution function f . However, using the additional property of the operator Q(f ) that the functions belonging to its kernel satisfy

Q(f ) = 0 iff f = M [f ], (4)
where the Maxwellian distributions

M [f ] = M [f ](x, v, t) can be expressed in terms of the set of moments U (x, t) by M [f ] = ρ (2πT ) dv/2 e -(v-u) 2 2T , 1 2 d v ρT = E - 1 2 ρ|u| 2 (5)
with T the temperature, one can get a closed system by replacing f with M [f ] in [START_REF] Filbet | Solving the Boltzmann equation in N log2 N[END_REF]. This corresponds to the set of compressible Euler equations which can be written as

∂U ∂t + ∇ x • F (U ) = 0, (6)
with

F (U ) = R dv M [f ]vφ(v) dv.
The simplest operator satisfying (2) and ( 4) is the relaxation operator [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-component systems[END_REF] Q

BGK (f ) = ν(M [f ] -f ), (7)
where ν = ν(x, t) > 0 defines the so-called collision frequency. Its values will be specified in the numerical test Section in order for the model to be as close as possible to the Boltzmann model described next. The classical Boltzmann operator reads

Q B (f ) = R dv S dv -1 B(|v -v * |, ω) (f (v )f (v * ) -f (v)f (v * )) dv * dω, ( 8
)
where ω is a vector of the unitary sphere S dv-1 ⊂ R dv . The post-collisional velocities (v , v * ) are given by the relations

v = 1 2 (v + v * + |q|ω), v * = 1 2 (v + v * + |q|ω), (9)
where q = v -v * is the relative velocity. The kernel B characterizes the details of the binary interactions, it has the form

B(|v -v * |, cos θ) = |v -v * |σ(|v -v * |, cos θ) ( 10
)
where the scattering cross-section σ, in the case of inverse k-th power forces between particles, can be written as

σ(|v -v * |, cos θ) = b α (cos θ)|v -v * | α-1 , ( 11
)
with α = (k -5)/(k -1). The special situation k = 5 gives the so-called Maxwell pseudo-molecules model with

B(v, v * , ω) = b 0 (cos θ). ( 12
)
For the Maxwell case the collision kernel is independent of the relative velocity. For numerical purposes, a widely used model is the variable hard sphere (VHS) model introduced by Bird [START_REF] Bird | Molecular gas dynamics and the direct simulation of gas flows[END_REF]. The model corresponds to b α (cos θ) = C α , where C α is a positive constant, and hence

σ(|v -v * |, cos θ) = C α |v -v * | α-1 . ( 13
)
In the numerical test Section we will consider the Maxwell molecules case when dealing with a velocity space of dimension d v = 2 and with the more realistic case of VHS molecules when dealing with the three dimensional case in velocity space:

d v = 3.
For comparison purposes we will also consider the simplified BGK model [START_REF] Jin | Asymptotic preserving (AP) schemes for multiscale kinetic and hyperbolic equations: a review[END_REF] and the compressible Euler system case (6).

The Fast Kinetic Scheme and the Fast Spectral Scheme

In this section we detail the numerical scheme. In the first part we discuss the Fast Kinetic Scheme (FKS) and in the second part the Fast Spectral method for the Boltzmann equation. The two solvers are connected by using splitting in time approaches as the ones described in [START_REF] Desvillettes | About the splitting algorithm for Boltzmann and B.G.K. equations[END_REF].

The Fast Kinetic Scheme (FKS)

The Fast Kinetic Scheme [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF][START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF] belongs to the family of so-called semi-Lagrangian schemes [START_REF] Crouseilles | A forward semi-lagrangian method for the numerical solution of the vlasov equation[END_REF][START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF][START_REF] Filbet | Conservative numerical schemes for the Vlasov equation[END_REF] which are typically applied to a Discrete Velocity Model (DVM) [START_REF] Bobylev | On approximation of the Boltzmann equation by discrete velocity models[END_REF][START_REF] Mieussens | Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics[END_REF] approximation of the original kinetic equation.

In order to introduce the scheme, let us truncate the velocity space by fixing some given bounds and set a cubic grid in velocity space of N points with ∆v the grid step which is taken equal in each direction. The continuous distribution function f is then replaced by a vector whose components are assumed to be approximations of the distribution function f at locations

v k : fk (x, t) ≈ f (x, v k , t). ( 14
)
The discrete velocity kinetic model consists then of a set of N evolution equations in the velocity space for fk , 1 ≤ k ≤ N , of the form

∂ t fk + v k • ∇ x fk = Q( fk ), (15)
where Q( fk ) is a suitable approximation of the collision operator Q(f ) at location k discussed in Section (3.3). Observe that, due to the truncation of the velocity space and to the finite number of points with which f is discretized, the moments of the discrete distribution function fk are such that

U (x, t) = k φ k fk (x, t) ∆v = U (x, t), ( 16
)
with

φ k = (1, v k , |v k | 2
) the discrete collision invariants, i.e. they are no longer those given by the continuous distribution f . This problem concerns all numerical methods based on the discrete velocity models and different strategies can be adopted to restore the correct macroscopic physical quantities [START_REF] Mieussens | Discrete velocity model and implicit scheme for the BGK equation of rarefied gas dynamics[END_REF][START_REF] Palczewski | A consistency result for a discrete-velocity model of the Boltzmann equation[END_REF][START_REF] Palczewski | Existence, stability, and convergence of solutions of discrete velocity models to the boltzmann equation[END_REF][START_REF] Gamba | Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states[END_REF]. At the same time, the spectral approach does not preserve the energy and the momentum of the system even if it approximates these quantities with spectral accuracy [START_REF] Pareschi | Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator[END_REF]. In order to solve this problem of loss of conservation we adopt an L 2 projection technique for the discretized distribution and the discretized collision operator which permits to project the discretized fk and Q( fk ) in the space of the distributions for which the moments are exactly the continuous macroscopic quantities we aim to preserve. We detail this approach in the following paragraph. For lightening the notations we omit the tilde in the rest of the paper, we only use it in the next paragraph in order to introduce the L 2 projection procedure. We will then suppose from now on that all the operations preserves the macroscopic moments if not differently stated. Let us now introduce a Cartesian uniform grid in the three dimensional physical space made of M points with ∆x a scalar which represents the grid step (the same in each direction) in the physical space. Further we define a time discretization t n+1 = t n + ∆t starting at t 0 , where ∆t is the time step defined by an opportune CFL condition discussed next. The time index n varies between 0 and N t so that the final time is t final = t nt = t 0 + N t ∆t. Each equation of system [START_REF] Crouseilles | Conservative semi-Lagrangian schemes for Vlasov equations[END_REF] is solved by a time splitting procedure. We recall here a first order splitting approach: first a transport step exactly solves the left-hand side, whereas a collision stage solves the right-hand side using the solution from the transport step as initial data:

f k,j n V k+1 f k+1,j n f k+2,j n f k+2,j n+1 f k+1,j n+1 f k,j n+1 V k j j k V k+2 j V k j j k V k+2 j V k+1
Transport stage -→ ∂ t f k + v k • ∇ x f k = 0, ( 17
)
Collisions stage -→ ∂ t f k = Q(f k ). ( 18
)
Transport step. Let f 0 j,k be the point-wise values at time t 0 of the distribution f , f 0 j,k = f (x j , v k , t 0 ). The idea behind the fast kinetic scheme is to solve the transport stage [START_REF] Shoucri | Numerical integration of the Vlasov equation[END_REF] continuously in space, see Fig. 1 for a sketch in the one dimensional setting. To this aim we define at the initial time the function f 0 k (x) as a piecewise continuous function for all x ∈ Ω j , where Ω j = [x j-1/2 ; x j+1/2 ] and Ω = j Ω j . Hence starting from data f 0 k,j at time index 0, the exact solution of ( 17) is simply

f * ,1 k = f 0 k (x -v k ∆t), ∀x ∈ Ω. ( 19
)
In other words, the entire function f 0 k is advected with velocity v k during ∆t unit of time and the * superscript indicates that only the transport step has been solved so far. The extension of this procedure to the generic time step n gives

f * ,n+1 k = f n k (x -v k ∆t), ∀x ∈ Ω, ( 20
)
where now, the key observation is that the discontinuities of the piecewise function f n k (x) do not lie on the interfaces of two different cells. Instead, the positions of the discontinuities depend entirely on the previous advection step and thus they may be located anywhere in the physical space. This means that if only the linear transport equation has to be solved, this approach gives the exact solution to the equation if the initial data is truly a piecewise constant function initially centered on the spacial mesh.

Collision step. The effect of the collisional step is to change the amplitude of f k (x). The idea is to solve the collision operator locally on the grid points, and, successively, extend these computed values to the full domain Ω. Thus we need to solve the following ordinary differential equation

∂ t f j,k = Q(f j,k ), (21)
where f j,k = f (x j , v k , t), for all velocities of the lattice k = 1, . . . , N and grid points j = 1, . . . , M . The initial data for solving this equation is furnished by the result of the transport step obtained by [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF] at points x j of the mesh at time t n+1 = t n + ∆t, i.e. f * ,n+1 k (x j ), for all k = 1, . . . , N , and j = 1, . . . , M . Then, the solution of ( 21), locally on the grid points, reads if, for simplicity, a forward Euler scheme is used as

f n+1 j,k = f * ,n+1 j,k + ∆t Q(f * ,n+1 j,k ), (22)
where f * ,n+1

j,k = f * ,n+1 k (x j ).
Many different time integrators can be employed to solve this equation. In particular special care is needed in the case in which the equation becomes stiff, refer to [START_REF] Dimarco | High order asymptotic-preserving schemes for the Boltzmann equation[END_REF][START_REF] Dimarco | Asymptotic preserving implicit-explicit Runge-Kutta methods for nonlinear kinetic equations[END_REF] for alternative strategies. Since the time integration of the collision term is not the issue considered in this paper, we considered the simplest possible scheme, but the FKS technique remains the same when other time integrators are employed. Equation [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF] furnishes a new value for the distribution f at time t n+1 only in the cell centers of the spacial cells for each velocity v k . However, one needs also the value of the distribution f in all points of the domain in order to perform the transport step at the next time step. Therefore, we define a new piecewise constant function Q k for each velocity of the lattice v k as

Q n+1 k (x) = Q(f * ,n+1 j,k ), ∀x such that f * ,n+1 k (x) = f * ,n+1 k (x j ). (23)
Said differently we make the fundamental assumption that the shape of Q(f k ) in space coincides with the one of f k . Thanks to the above choice one can rewrite the collision step in term of spatially reconstructed functions as

f n+1 k (x) = f k (x, t n + ∆t) = f * ,n+1 k (x) + ∆t Q n+1 k (x). ( 24
)
This ends one time step of the FKS scheme.

Concerning the transport part of the scheme, the time step ∆t is constrained by a CFL like condition of type ∆t max

k |v k | ∆x ≤ 1 = CFL. ( 25
)
The time step constraint for the collision step depends on the choice of the operator Q. Since in the numerical test section we used both a BGK and a Boltzmann operator, the time step constraint for the interaction part has been chosen as the minimum time step which gives stability in the solution of the ODE [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF] independently on the type of collision operator employed. As observed in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF] the transport scheme is stable for every choice of the time step, being the solution for a given fixed reconstruction performed exactly. Nonetheless the full scheme being based on a time splitting technique, the error is of the order of ∆t in the case of first order splitting or of order (∆t) q for a splitting of order q. This suggests to take the usual CFL condition for the transport part in order to maintain the time error small enough.

To conclude this Section let us observe that time accuracy can be increased by high order time splitting methods, while spacial accuracy can be increased close to the fluid limit to a nominally second-order accurate scheme by the use of piecewise linear reconstructions of state variables, see the details in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF]. Spacial accuracy for all regime can also be increased by using high order polynomial reconstruction for the distribution f . This work is in progress and results are discussed in [START_REF] Dimarco | A new class of high order semi-Lagrangian schemes for rarefied gas dynamics[END_REF].

Conservation of macroscopic quantities

In order to preserve mass, momentum and energy in the scheme, we employ the strategy proposed in [START_REF] Gamba | Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states[END_REF]. We consider one space cell, the same renormalization of f should be considered for all spatial cells. This step is performed at the beginning, i.e. t = 0. for the distribution f (x j , v k , t = 0) and after each collision step [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF] which causes loss of momentum and energy due to the spectral discretization. Let f = f 1 , f 2 , . . . , f N T be the distribution function vector at t = t n at the center of the cell and let f = (f 1 , f 2 , . . . , f N )

T be the unknown corrected distribution vector which fulfills the conservation of moments. Let

C (dv+2)×N =   (∆v) dv k v k (∆v) dv k |v k | 2 (∆v) dv k  
be a matrix of coefficients depending on the discretization parameters and U (dv+2)×1 = (ρ ρu E) T be the vector of conserved quantities, namely density, momentum and energy. Conservation can be imposed solving a constrained optimization formulation:

Given f ∈ R N , C ∈ R (dv+2)×N , and U ∈ R (dv+2)×1 , find f ∈ R N such that (26) f -f 2 2 is minimized subject to the constrain Cf = U.
The solution of this minimization problem can be analytically obtained by employing the Lagrange multiplier method. Let λ ∈ R dv+2 be the Lagrange multiplier vector. The corresponding scalar objective function to be minimized is given by

L(f, λ) = N k=1 | f k -f k | 2 + λ T (Cf -U ). (27)
Then, by nullifying the derivative of L(f, λ) with respect to f k we get

f = f + 1 2 C T λ ( 28
)
while the Lagrange multipliers are obtained by solving

CC T λ = 2(U -C f ). ( 29
)
In particular, the above expression says that the value of λ is uniquely determined by λ = 2(CC T ) -1 (U -C f ). Back substituting λ into (28) finally provides

f = f + C T (CC T ) -1 (U -C f ). ( 30
)
Note that matrix C T (CC T ) -1 is pre-computed and stored as being constant for each simulation. If a discretization of the Maxwellian distribution M [f ] is needed, as for example when using the BGK model, the same procedure should be applied to the function M (x j , v k , t n ) in order to assure conservation of the first three moments at each instant of time at which the Maxwellian function M (x j , v k , t n ) is invoked. In this case by defining with E[f ] the approximated Maxwellian and by Ẽ[f ] = M (x j , v k , t n ) the pointwise value of the equilibrium distribution, we get mimicking (30)

E[U ] = Ẽ[f ] + C T (CC T ) -1 (U -C Ẽ[f ]). ( 31
)
This ends the description of the procedure which permits to conserve macroscopic quantities.

Fast Spectral Scheme (FSS) to discretize the Boltzmann collision operator

The fast spectral discretization of the Boltzmann operator employed in this work is described in this section. To this aim, we focus again on a given cell x j at a given instant of time t n . The same computation is repeated for all cells x j , j = 1, .., M and times t n , n = 0, .., N t . Moreover, since the collision operator acts only on the velocity variable, to lighten the notation in this paragraph, only the dependency on the velocity variable v is considered for the distribution function f , i.e. f = f (v).

The first step to construct for our spectral discretization is to truncate the integration domain of the Boltzmann integral [START_REF] Jin | Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms[END_REF] as done for the distribution f . As a consequence, we suppose the distribution function f to have compact support on the ball B 0 (R) of radius R centered in the origin. Then, since one can prove that supp (Q(f )(v)) ⊂ B 0 ( √ 2R), in order to write a spectral approximation which avoid aliasing, it is sufficient that the distribution function f (v) is restricted on the cube [-T, T ] dv with T ≥ (2 + √ 2)R. Successively, one should assume f (v) = 0 on [-T, T ] dv \ B 0 (R) and extend f (v) to a periodic function on the set [-T, T ] dv . Let observe that the lower bound for T can be improved. For instance, the choice

T = (3 + √ 2)R/2
guarantees the absence of intersection between periods where f is different from zero. However, since in practice the support of f increases with time, we can just minimize the errors due to aliasing [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF] with spectral accuracy.

To further simplify the notation, let us take T = π and hence R = λπ with λ = 2/(3 + √ 2) in the following. We denote by Q R B (f ) the Boltzmann operator with cut-off. Hereafter, using one index to denote the d v -dimensional sums, we have that the approximate function f N can be represented as the truncated Fourier series by

f N (v) = N/2 k=-N/2 fk e ik•v , ( 32
) fk = 1 (2π) dv [-π,π] dv f (v)e -ik•v dv. ( 33
)
We then obtain a spectral quadrature of our collision operator by projecting ( 8) on the space of trigonometric polynomials of degree less or equal to N , i.e.

Qk =

[-π,π] dv Q R B (f N )e -ik•v dv, k = -N/2, . . . , N/2. ( 34
)
Finally, by substituting expression [START_REF] Bobylev | Difference scheme for the Boltzmann equation based on the fast Fourier transform[END_REF] in [START_REF] Gamba | Spectral-Lagrangian methods for collisional models of non-equilibrium statistical states[END_REF] one gets after some computations

Qk = N/2 l,m=-N/2 l+m=k fl fm β(l, m), k = -N, . . . , N, (35)
where

β(l, m) = B(l, m) -B(m, m) are given by B(l, m) = B0(2λπ) S dv -1 |q|σ(|q|, cos θ)e -i(l•q + +m•q -) dω dq. ( 36
) with q + = 1 2 (q + |q|ω), q -= 1 2 (q -|q|ω). (37)
Let us notice that the naive evaluation of [START_REF] Gamba | A conservative spectral method for the Boltzmann equation with anisotropic scattering and the grazing collisions limit[END_REF] requires O(n 2 ) operations, where n = N 3 . This causes the spectral method to be computationally very expensive, especially in dimension three. In order to reduce the number of operations needed to evaluate the collision integral, the main idea is to use another representation of ( 8), the so-called Carleman representation [START_REF] Carleman | Sur la théorie de l'équation intégrodifférentielle de Boltzmann[END_REF] which is obtained by using the following identity

1 2 S dv -1 F (|u|σ -u) dσ = 1 |u| d-2 R dv δ(2 x • u + |x| 2 ) F (x) dx.
This gives in our context for the Boltzmann integral

Q B (f ) = R dv R dv B(x, y)δ(x • y) [f (v + y) f (v + x) -f (v + x + y) f (v)] dx dy, ( 38
) with B(|x|, |y|) = 2 dv-1 σ |x| 2 + |y| 2 , |x| |x| 2 + |y| 2 (|x| 2 + |y| 2 ) -dv -2 2 . ( 39
)
This transformation yields the following new spectral quadrature formula

Qk = N/2 l,m=-N/2 l+m=k βF (l, m) fl fm , k = -N, ..., N (40)
where βF (l, m) = BF (l, m) -BF (m, m) are now given by

BF (l, m) = B0(R) B0(R) B(x, y) δ(x • y) e i(l•x+m•y) dx dy. ( 41
)
Now, in order to reduce the number of operation needed to evaluate [START_REF] Wu | Influence of intermolecular potentials on rarefied gas flows: Fast spectral solutions of the Boltzmann equation[END_REF], we look for a convolution structure.

The aim is to approximate each βF (l, m) by a sum βF (l, m)

A p=1 α p (l)α p (m),
where A represents the number of finite possible directions of collisions. This finally gives a sum of A discrete convolutions and, consequently, the algorithm can be computed in O(A N log 2 N ) operations by means of standard FFT technique [START_REF] Canuto | Spectral methods in fluid dynamics[END_REF].

In order to get this convolution form, we make the decoupling assumption

B(x, y) = a(|x|) b(|y|). ( 42
)
This assumption is satisfied if B is constant. This is the case of Maxwellian molecules in dimension two, and hard spheres in dimension three, the two cases treated in this paper. Indeed, using kernel [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF] in [START_REF] Filbet | Analysis of spectral methods for the homogeneous Boltzmann Equation[END_REF], one has

B(x, y) = 2 dv-1 C α (|x| 2 + |y| 2 ) -dv -α-2 2 , so that B is constant if d v = 2, α = 0 and d v = 3, α = 1.
We start by dealing with dimension 2 and B = 1, i.e. Maxwellian molecules. Here we write x and y in spherical coordinates x = ρe and y = ρ e to get

BF (l, m) = 1 4 S 1 S 1 δ(e • e ) R -R e iρ(l•e) dρ R -R
e iρ (m•e ) dρ de de .

Then, denoting φ 2 R (s) = R -R e iρs dρ, for s ∈ R, we have the explicit formula φ 2 R (s) = 2 R Sinc(Rs),
where Sinc(x) = sin(x)

x . This explicit formula is further plugged in the expression of BF (l, m) and using its parity property, this yields BF (l, m)

= π 0 φ 2 R (l • e θ ) φ 2 R (m • e θ+π/2 ) dθ.
Finally, a regular discretization of A equally spaced points, which is spectrally accurate because of the periodicity of the function, gives

BF (l, m) = π M A p=1 α p (l)α p (m), (43)
with

α p (l) = φ 2 R (l • e θp ), α p (m) = φ 2 R (m • e θp+π/2 ) ( 44
)
where θ p = πp/A. Now let us deal with dimension d = 3 with B satisfying the decoupling assumption [START_REF] Pareschi | On the stability of spectral methods for the homogeneous Boltzmann equation[END_REF]. First we switch to the spherical coordinates for BF (l, m):

BF (l, m) = 1 4 S 2 ×S 2 δ(e • e ) R -R ρa(ρ)e iρ(l•e) dρ R -R ρ b(ρ )e iρ (m•e ) dρ dede .
Then, integrating first e on the intersection of the unit sphere with the plane e ⊥ gives

BF (l, m) = 1 4 e∈S 2 φ 3 R,a (l • e) e ∈S 2 ∩e ⊥ φ 3 R,b (m • e ) de de
where

φ 3 R,a (s) = R -R ρ a(ρ) e iρs dρ.
This leads to the following decoupling formula

BF (l, m) = e∈S 2 + φ 3 R,a (l • e) ψ 3 R,b Π e ⊥ (m) de
where S 2 + denotes the half sphere and

ψ 3 R,b Π e ⊥ (m) = π 0 sin θ φ R,b |Π e ⊥ (m)| cos θ dθ.
Now, in the particular case where B = 1, i.e. the hard sphere model, we can explicitly compute the functions φ 3 R . These are

φ 3 R (s) = R 2 2 Sinc(Rs) -Sinc 2 (Rs/2) , ψ 3 R (s) = 2 R 2 Sinc 2 (Rs/2).
Taking a spherical parametrization (θ, ϕ) of e ∈ S 2 + and uniform grids of respective size A 1 and A 2 for θ and ϕ (again spectrally accurate because of the periodicity of the function) leads to the following quadrature formula for BF (l, m)

BF (l, m) = π 2 A 1 A 2 A1,A2 p,q=0 α p,q (l)α p,q (m) where α p,q (l) = φ 3 R,a l • e (θp,ϕq) , α p,q (m) = ψ 3 R,b Π e ⊥ (θp ,ϕq ) (m) , φ 3 R,a (s) = R -R ρ a(ρ) e iρs dρ, ψ 3 R,b (s) = π 0 sin θ φ 3 R,b (s cos θ) dθ,
and for all p and q (θ p , ϕ q ) = p π A 1 , q π A 2 .

Numerical implementation

We discuss in this part some aspects relative to the numerical implementation of the method.

Algorithm for the transport part

The method described in the previous sections can be resumed into two main actions: transport and interaction. In our implementation we have adopted a particle like interpretation of the FKS scheme which helps in reducing the computational effort and the memory requirement. In this interpretation, each point of the quadrature of the phase space is represented by a particle with a given mass, velocity and position. The transport phase causes the particle to move in the physical space, while the interaction phase causes the mass of each particle to change. In this setting, the distribution function f can be expressed as

f (x, v, t) = Np i=1 m i (t) δ(x -x i (t))δ(v -v i (t)), v i (t) = v k , ( 45
)
where x i (t) represents the particle position, v i (t) the particle velocity, m i the particle mass. Moreover, the particle velocity corresponds to the quadrature point chosen to discretize the velocity space. The effect of the transport is a simple shift of particles to their new positions according to

x i (t + ∆t) = x i (t) + v i (t)∆t. ( 46
)
The collision step acts only locally and changes the velocity distribution. The particle interpretation of this part of the scheme consists in changing the mass of each particle through the formula ( 22) which reads

m i (t + ∆t) = m i (t) + ∆t Q(v i ), (47)
where Q(v i ) corresponds to the approximation of the collision integral in the center of the cell evaluated at location v i by means of the fast spectral method described before. Let observe that thanks to the uniform grid in velocity and physical space, the number of particles is the same in each cell and remains constant in time. This permits to consider the motion of only a fixed subset of particles which belongs to only one chosen cell. The relative motion of all other particles in their cells being the same [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF]. This gives great computational advantage since the number of particle to track in time is greatly reduced. The only information which should be tracked remains the mass of each individual particle which is different for every velocity v i and position x i . This reinterpretation of the scheme permits to strongly reduce the computational cost related to the transport part as shown in the numerical test section.

Profiling and parallelization strategy

The Fast Kinetic scheme for the BGK equation was shown to be very efficient and extremely parallelizable (with a shown acceleration very close to the ideal scaling) in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part III: High-performance-computing[END_REF] on mild parallel architecture, namely using OpenMP on maximum 64 threads and on two Graphical Processing Units (GPU) using CUDA framework. These light parallel infrastructures are unfortunately not sufficient when Boltzmann operator is to be simulated. To this aim, let us briefly discuss the profiling of the collision algorithm. As we have already seen in the description of the method in the previous section, the resolution of Boltzmann operator requires several passages from the physical space to the Fourier space, which implies several calls to (inverse) Fourier transforms within nested loops which are due to the convolutive form of the integral. This interaction step covers about 98% of the total cost of a 2D×2D simulation on a serial machine, see table 1. The situation is alike in the case of 3D×3D simulations, and, thus, results are not reported. From table 1 we deduce that the collision step is the part of the code which demands to be dealt with great attention if any gain in performance is expected. More precisely, if we give a closer look to the computational cost related to the computation of the collision term Q(f ) in one fixed spacial cell, the major cost (beyond 90%) is spent by the Fourier transform routines. In our implementation, the collision is performed without any communication with the local neighborhood. In other words, once the N values of the distribution functions f j,k , k = 1, .., N are stored then Q(f ) is computed locally. Because N is not too large, typically the number of points chosen to discretize the velocity space goes from N = 16 3 to N = 64 3 , then the memory requirement is not extremely demanding. Finally, as seen in the previous Section the computation of Q(f ) is based on a loop over the number of angles θ p with p = 1, .., A in dimension two, and, a double loop over the angles θ p and

Main routines Cost CPU %

Sub-routines

Sub-cost CPU % φ q in dimension three with p = 1, ..A 1 and q = 1, .., A 2 with a typical number of chosen angles of eight. After a first Fourier transform of the distribution f , for each iterate of these loops, the coefficients of the quadrature formula are computed and then three calls to the inverse fast Fourier transform are done. The results of the FFT are gathered into temporary complex variables later arranged to form the solution of Q(f ). The costs related to the four main routines needed for the evaluation of the operator namely: Fourier transform, coefficients computation, inverse Fourier transform and sum of the obtained values to form Q(f ), are detailed in table 1.

We discuss now some aspects related to the parallelization strategy. The framework described is extremely well suited for parallelization. The strategy adopted is to divide the space domain in several subdomains because the interaction step, the most expensive one, does not demand any communication between spatial cells. Each cell can be therefore treated independently from the others. Then, the idea employed has been to use the OpenMP on shared memory systems and make each computational core responsible for a subset of a space mesh. This approach requires very little modification of a sequential code and gives strong scaling close to ideal as shown in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part III: High-performance-computing[END_REF] for the BGK case. A different possibility is to compute the collision kernel, and all the related Fast Fourier Transforms on a GPU. Even if, this involves substantial amount of slow communications between CPU and GPU, the computational complexity of this part is so large that the communication time does not play a fundamental role. This approach requires more programming effort, compared to the OpenMP parallelization strategy as several routines have to be completely rewritten. In particular, the code has to be adapted to the Single Instruction, Multiple Threads (SIMT) parallel programming mode imposed by the GPU architecture. That is to say, parts of the code that can be run in parallel on large number of threads have to be identified and rewritten as so called CUDA kernels -functions that are executed multiple times in parallel by different CUDA threads. Moreover, the GPU-CPU memory data transfer has to be taken into account. The payoff is however higher and, at least for computations which does not require too large memory storage, this is probably the best choice.

Unfortunately, when the number of points needed to perform a simulation grows in the case the Boltzmann operator is used, the above implementation strategies on shared memory systems are not sufficient any more. This is true particularly in the full three dimensional case. Luckily, the method described in this work can be easily adapted for its application in distributed memory systems. The idea adopted in this case is to distribute spatial degrees of freedom over computational nodes, keeping on every node a complete set of velocity space points. The computational domain is then decomposed into slices along the Z direction (see Fig. 2). The slices are successively distributed over different MPI processes. Every node will therefore compute the collision term for a part of computational domain. This can be done on each node using the OpenMP or GPU parallelization strategy, depending on the architecture at disposal. As the update of the primitive variables requires an exchange with neighboring spatial cells, some of the particle masses, which contain all the information related to the distribution function f in our implementation, need to be broadcast to other computational nodes. This data exchange, which is typically a bottleneck in the MPI computations, and, usually demands large efforts in order to be minimized the internodal communications, in the case of the Boltzmann equation, takes only a small part of the total runtime, even if, in our first MPI implementation, the domain decomposition is far from being optimal. We postpone to a future work the discussion related to a more efficient MPI parallelization strategy and all details related to this type of implementation.

Numerics

In this section, we validate the proposed SFKS (Spectral-FKS) method which couples the fast semi-Lagrangian and the fast spectral scheme. The testing methodology is divided into four parts.

Part 1. Sanity checks. We only consider the space homogeneous Boltzmann equation in two and three dimensions in velocity space. Two problems from [START_REF] Pareschi | Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator[END_REF] for which the analytical solution is known are simulated. The results show that our implementation of the fast spectral discretization of the Boltzmann operator is correct.

Part 2. Using the SFKS method in one dimension of physical space and two or three dimensions in velocity space, we show on Riemann like problems the differences between the BGK and the Boltzmann models.

The reported results justify the use of the more complex and costly collisional Boltzmann operator especially in situations far from the thermodynamical equilibrium. We also report the details of the computational costs as well as the numerical convergence of the method for an increasing number of points in the physical or velocity space.

Part 3. Using the SFKS method in the two dimensional case in space and velocity, we show that using the Boltzmann or the BGK operator leads to notably different results. We consider a regular vortex like problem and a reentry like problem with evolving angle of attack. Performances and profiling study of the SFKS approach supplemented with BGK or Boltzmann collisional operator are provided for refined grids in space and time. Details about the computational costs of the different part of the scheme, scalability with respect to the number of cells are also furnished in order to characterize the behavior of the method as precisely as possible.

Part 4. We test the SFKS method in the full three dimensional case in space and velocity for an unsteady test problem. The interaction of a flying object with the surrounding ambient gas is simulated. We compare the results obtained with the Boltzmann collisional operator with those given by a simpler # of points L 1 error L 2 error 8 2 6.3 × 10 -3 3.8 × 10 -3 16 2 5.7 × 10 -4 2.9 × 10 -4 32 2 6.5 × 10 -5 7.2 × 10 -5 BGK model. Moreover, we measure the cost of such a 3D×3D simulation in terms of CPU, cost per degree of freedom and scalability with respect to the number of cells. Some code profiling is also provided to measure the cost of the main components of the scheme.

Finally, the ability of the whole numerical code to run on several parallel environments is tested: CPU (OpenMP, MPI) and GPU (CUDA) types of parallelism frameworks are considered.

Part 1. Numerical results for the space homogeneous case

In this part, we validate our spectral discretization implementation of the Boltzmann collisional operator. To this aim, we re-employ the test cases considered in [START_REF] Pareschi | Numerical solution of the Boltzmann equation I: Spectrally accurate approximation of the collision operator[END_REF] and [START_REF] Filbet | Accurate numerical methods for the Boltzmann equation[END_REF] for the homogeneous two dimensional and three dimensional Boltzmann equation. For both situations, exact solutions are available and briefly recalled in the following.

Test 1.1. Convergence to equilibrium for the Maxwell molecules in dimension two.

We consider two dimensional in velocity Maxwellian molecules and the following space homogeneous initial condition

f (v, t = 0) = v 2 π exp(-v 2 ). (48)
The analytical solution (the so-called BKW one) is given for all times t by [START_REF] Bobylev | Exact solutions of the Boltzmann equation[END_REF][START_REF] Krook | Exact solutions of the Boltzmann equation[END_REF]:

f (v, t) = 1 2πS 2 exp(-v 2 /2S) 2S -1 + 1 -S 2S v 2 , ∀t > 0, (49)
with S = 1 -exp(-t/8)/2. The final time is set to t final = 10 and the time step is equal to ∆t = 0.02. The test is performed for N = 8 2 , 16 2 , 32 2 points and for eight discrete angles θ p . In order to keep the aliasing error smaller than the spectral error, the velocity domain is set to [-4, 4] 2 for N = 8 2 , [-6, 6] 2 for N = 16 2 and to [-9, 9] 2 for N = 32 2 points. This choice helps fighting back the aliasing error behavior since it increases with the number of points. Therefore, to have comparable aliasing errors for the three meshes, we enlarge the size of the velocity space when the number of points increases. We report the discrete L 1 and L 2 norms of the error for the distribution function f in Table 2. Moreover, in Figure 3, we report the distribution function f (v x , v y = 0, t) at different times (t = 1, 2 and 10) when N = 32. The solution is plotted versus the exact solution. The results clearly show the convergence of the method towards the exact solution. Last in Figure 3 right panel, we present the time evolution of the L 1 error for the three configurations. Consistently the errors decrease by about one order of magnitude. The last recorded results, at t final = 10 corresponds to the Figures reported in table 2.

Test 1.2. Space homogeneous comparison between the BGK and the Boltzmann model. Maxwellian molecules.

In Figure 4, we compare the convergence to equilibrium for the BGK and Boltzmann models on the test case described in the previous paragraph using N = 32 2 points on a domain [-9, 9] 2 . The L 1 norm of the difference between the two distribution functions f BGK (v, t) and f Boltz (v, t) as a function of time is shown. As expected the differences are increasing at early stages of the relaxation towards the equilibrium to reach a maximum difference around time t 1.66. Then, the two models slowly converge towards the same equilibrium solution. The red marked iterations in Figure 4 are further depicted in Figure 5, where we show the details of the difference in time between the two distribution functions. The same azimuthal scale is employed to ease the comparison.

Test 1.3. Convergence to equilibrium in dimension three. Hard sphere molecules.

The following initial condition is considered [START_REF] Filbet | Accurate numerical methods for the Boltzmann equation[END_REF] f

(v, t = 0) = 1 2(2πσ 2 ) 3/2 exp - |v -v 1 | 2 2σ 2 + exp - |v + v 1 | 2 2σ 2 , ( 50
)
where

σ 2 = 0.2 and v 1 is v 1 = (v x , v y , v z ) = (-1, -1, -0.25).
The final time is set to t final = 2, the time step is ∆t = 0.05, while the velocity domain is [-7; 7] 3 discretized with N = 32 3 points. In Figure 6, we report the initial data while in Figure 7, we report the relaxation to the equilibrium for the hard sphere model on a fixed plane passing through the points of coordinates ±(-1, -1, -0.25) and of normal (-0.1, -0.1, 1)1 .

In Figure 8, we report as before the comparison between the Boltzmann and the BGK model measuring the L 1 difference between the two distribution functions in time. Finally, in Figure 9 are presented the details of such differences between the distribution functions for the red marked iterations of figure 8 at v z = 0. As for the Maxwellian molecules in two dimensions, at the beginning of the simulation we observe the larger deviations between the two models while at the end they both converge to the same limit solution as expected.

The results of these homogeneous tests (in two and in three dimensions in velocity) are meant to verify and validate our implementation of the spectral discretization of the Boltzmann collisional operator. Moreover, the differences observed between the two models justify the necessity of the Boltzmann operator for space non homogeneous situations. We can now investigate the coupling of the Boltzmann operator with the Fast semi-Lagrangian method for several space non homogeneous cases.

Part 2. Numerical results for the one dimensional case in space.

In this part, we focus on solving the one dimensional in space Boltzmann and BGK equations considering a two or three dimensional dimensional velocity space setting. The purposes are twofold. First, we numerically demonstrate that the FKS and the spectral accurate Boltzmann kernel solver may be appropriately coupled and that they provide a valid kinetic solver for Boltzmann equations. Second, we show that BGK and Boltzmann models provide different results, justifying the use of a more sophisticated model. For these tests, a classical Riemann problem with Sod like initial data is considered

ρ = 1, u = 0, T = 2.5, if x ≤ L/2, ρ = 0.125, u = 0, T = 0.25, if x > L/2,
with Ω = [0, 2]. Initial local thermodynamics equilibrium is considered for all tests: f (x, v, t = 0) = M [f ](x, v, t = 0). The velocity space is set to L v = [-15, 15] for all cases. Dirichlet boundary conditions are set on the left/right boundaries of Ω. The two kinetic models are solved by employing a time rescaling factor in order to put in evidence the role of the collisions in the solutions. The rescaled equations reads

∂ t f + v • ∇ x f = 1 τ Q(f ), ( 51
)
where τ is the rescaled parameter (the frequency of relaxation) which plays the role of the non dimensional Knudsen number. Smaller is the relaxation frequency, faster is the relaxation of the distribution function towards the equilibrium state. However, the exact rate of convergence is dependent on the type of collision considered either BGK or Boltzmann (Maxwellian molecules or hard spheres).

Test 2.1. Numerical convergence of the Boltzmann equation. The two dimensional in velocity

Maxwellian molecules case. The two dimensional velocity space is first considered leading to a space/velocity mesh of the form M ×N with N = 64 2 and a varying number of space cells M . In Figure 10, we present the space convergence results for the density, the velocity and the temperature when the Boltzmann operator is solved. Successively refined (doubled) spatial meshes are employed, from 50 to 400 up to final time t final = 0.15. From these data we can observe that the simulation results seem to converge towards a given numerical solution in both cases τ = 10 -3 (left panels) and τ = 10 -4 (right panels). We can also observe that the increase in mesh resolution is profitable especially for smaller τ . This is the same behavior observed in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF]. In fact, the scheme precision decreases as the equilibrium state is approached, being virtually exact in non collisional or almost non collisional regimes. The loss of precision observed in fluid dynamic regimes can be recovered with a similar technique as the one proposed in [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF]. Here, however we do not consider this possibility. The CFL condition employed is the following

∆t ≤ min ∆x |v max | , τ ρ , ( 52
)
where the first term on the right hand side of the above equation comes from the will of keeping the error small enough in the splitting scheme, while the second term is due to the stability restriction in the solution of the space homogeneous problem when Maxwellian molecules are employed. In fact, in this case the loss part of the collision integral Q -(f ) can be estimated, giving

Q -(f ) = ρ(f )f . 0.

Test 2.2. Comparisons between the BGK model and the Boltzmann model. The two dimensional in velocity Maxwellian molecules case.

Here, the BGK and Boltzmann models are simulated for the same Sod-like problem. In order to have fairest as possible comparisons between the two models, we choose ν = ρ for the BGK model. This choice permits to have the same loss part for the two models since for Maxwellian molecules the loss part is close to ρf as stated in the previous paragraph. We fix M = 800 spatial cells and N = 64 2 cells in velocity space. This permits to consider almost converged results. In Figure 11, we present the results when τ = 10 -3 (left) and τ = 10 -4 (right) for the density (top), the velocity (middle top), the temperature (middle bottom) and the heat flux (bottom). The very first observation is the relative large differences between the two solutions. The Boltzmann solution systematically presents more dissipated waves than BGK solution. However, the main waves are located in the same positions for the two models. The difference in term of the macroscopic quantities can be of the order of 10% -25% of the solution. Systematically, the BGK model underestimates the values of the heat flux. In Figure 12 we report the absolute value of the difference between the two distribution functions f BGK (x, v, t) and f Boltz (x, v, t) at different locations x j and at final time t final as a function of the velocity variables v. The vertical scale is kept constant for all panels, only the color scale is adapted to the values reported.

Test 2.3. Numerical convergence of the Boltzmann equation. The three dimensional in velocity hard

sphere molecules case. Let us focus on the one dimension in space and three dimensions in velocity case. We consider the same Sod-like problem up to final time t final = 0.5 and two different collision frequencies τ = 10 -2 and τ = 2 10 -3 . The space/velocity mesh chosen is of the form M × N with N = 32 3 uniformly spread on a velocity domain [-16; 16] 3 and varying number of space cells M . In Figure 13, we present the space convergence results for the density, the velocity and the temperature for successively refined spatial meshes from 32 to 128. From these data we can observe that the simulation results seem to converge towards the same numerical solution in both cases τ = 10 -2 (left panels) and τ = 2 10 -3 (right panels). The CFL condition employed in this case is the following

∆t ≤ min ∆x |v max | , τ µ , ( 53
)
where the second term is due to the stability restriction in the solution of the space homogeneous problem when hard sphere molecules are employed. In this case, the loss part of the collision integral Q -(f ) can be only estimated, giving

Q -(f ) = L(f )f ≤ C α 4π (
2λπ) α and thus µ ≥ C α 4π(2λπ) α in order to ensure that the gain part Q + (f ) is positive and monotone. We have also performed comparisons between the BGK and the Boltzmann models in this one dimensional setting. The results are close to the ones obtained for the two dimensional Maxwell molecules, i.e. the BGK model tends to over-relax the distribution function to the equilibrium state and for this reason we do not report them. However, since the computational costs involved in the approximation of the Boltzmann integral in the three dimensional setting are much larger than those obtained with the simpler two dimensional model, in the next paragraph we analyze the performances of the scheme proposed in both cases.

Performances

In this part, we analyze the performances of our scheme in the one dimensional in space setting. In Table 3 we monitor the CPU time for the Boltzmann and BGK operators when a successively refined spatial mesh is considered for a fixed two and three dimensional velocity mesh. The data reported are relative to the same Sod-like problem considered in the previous paragraphs. We employ 64 2 velocity cells in the two dimensional case and 32 3 in the three dimensional case. The number of spatial cells is M = 100, 200 and 400 in the two dimensional case and M = 32, 64, and 128 in the three dimensional one. More precisely, we measure the CPU time T in seconds, the time per cycle as T cycle = T/N cycle , the time per cycle per spatial cell by T cell = T/N cycle /M and the time per cycle per degree of freedom (d.o.f)

T dof = T/N cycle /(M × N ).
The OpenMP parallel version of the scheme is used on a laptop having 8 threads (HP ZBook with Intel 8 Core i7-4940MX CPU @ 3.10GHz on Ubuntu 15.10 (64 bits)). These simulations have been run in parallel on non dedicated computer, as such the results are to be understood as rough estimates. By no mean we pretend that they cannot be improved. In the two dimensional case in velocity, the Boltzmann model results are about 15 times more expensive than the BGK model, while in three dimensions, it can reach about 250 times. Notice that, as expected, these ratio do increase when the collision frequency τ becomes smaller, for instance when τ = 2 × 10 -3 in three dimensions then the ratio in CPU time between the two models is of the order of 500. Let us notice that the number of cycles for hard sphere molecules is fixed for the reported simulations. This is due to the fact that the hard sphere model considered has a stability condition which is more restrictive than the CFL condition chosen for the transport part as opposite to the Maxwellian molecule model for the two dimensional case. For fair comparison, the number of cycles of the two models has been kept constant in the computation of the costs, even if, the space homogeneous BGK model does not present stability requirements. Next in table 4 we monitor the CPU time of the simulation when the number of cells is fixed in space to M = 100. The number of velocity cells increases from N = 8 2 to 64 2 in two and from N = 8 3 to 64 3 in three dimensions. The same Sod test case as in table 3 is simulated.

Because the CPU times have been obtained by parallel simulations, the analysis of such table must take into account the fact that the dimension of the mesh plays a role in the performance. For instance, in two dimensions, the smaller mesh (100 × 8 2 ) usually presents a larger CPU time per degree of freedom. This is due to the fact that a decent amount of cells is needed to observe the benefit of using a parallel machine.

For the Boltzmann simulations, in three dimensions, 8 velocity cells in each direction is not large enough to have a stable simulation. This is probably due to the large loss of energy caused by the spectral method which can be only partly cured by the L 2 projection technique detailed in Section 3.2. In fact, even if such renormalization permits to keep the correct energy after the collision step, it may transform the distribution function in an unphysical manner. This operation repeated multiple times give rise to instability in the spectral scheme. However, for completeness using the first completed iterations we report in the table an approximation of the CPU time that would be needed using this small amount of points. The Hard Sphere molecule simulations, as can be seen on the Figures reported, are extremely time consuming already in the one dimensional case in space. In two dimensions, the cost of the Boltzmann model is about 3 to 8 times (from the smallest to largest mesh) more expensive than the BGK model. In three dimensions, this ratio ranges between 40 and 700 times.

Sod like Riemann problem in 1D×2D τ = 10 -3 Model Velocity Cell # Cycle Time T/cycle T/cell T/d.o.f N Vel. M × N N cycle T (s) T cycle (s) T cell (s) T dof (s) BGK 64 2 [-

Part 3. Numerical results for the space two dimensional case.

In this part, we focus on solving the two dimensional in space and velocity Boltzmann and BGK equations. The purposes are twofold. First, we want to show the differences which arises between the two models. Second, we want to analyze the performances of the method in the two dimensional setting by monitoring the cost of such simulations and by performing a profiling of the scheme in terms of the main routines in order to highlight the eventual bottlenecks. This permits to understand in which part of the scheme one Sod like Riemann problem in 1D×2D τ = 10 The test case consists of an isentropic vortex in motion initially introduced for the collisional regime, i.e. the compressible Euler equations, in two dimensions in [START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]. This problem has an exact smooth solution expressed analytically in the fluid regime. The computational domain is Ω = [0, 10] 2 . The ambient flow is characterized by a gas with density, mean velocity and temperature respectively of ρ ∞ = 1.0, u x,∞ = 1.0, u y,∞ = 1.0, T ∞ = 1.0. A vortex is centered at (x v , y v ) = (5, 5) and supplemented at the initial time t = 0 with conditions u x (t = 0) = u x,∞ + δu, u y (t = 0) = u y,∞ + δv, T (t = 0) = T ∞ + δT with

δu x = -y β 2π exp 1 -r 2 2 , δu y = x β 2π exp 1 -r 2 2 , δT = - (γ -1)β 8γπ 2 exp 1 -r 2 ,
where r is the Euclidean distance in the two dimensional space, i.e. r = x 2 + y 2 , and x = x -x v , y = y -y v . The vortex strength depends on the value β fixed here to 5.0. The initial density is given by

ρ(t = 0) = ρ ∞ T (t = 0) T ∞ 1 γ-1 = 1 - (γ -1)β 8γπ 2 exp 1 -r 2 1 γ-1 . ( 55
)
Periodic boundary conditions are prescribed everywhere. At the final time chosen: t final = 10, the vortex is back to its original position and, in the collisional regime, the initial and final conditions are alike. In a rarefied regime, the exact solution is not known, but at least the cylindrical symmetry of the problem must be retrieved. M = 100 × 100 uniform spatial cells are considered on domain [0; 10] 2 with N = 32 × 32 uniform velocity cells on a velocity domain [-7.5; 7.5] 2 . The relaxation frequency is fixed to τ = 10 -1 . The initial data are presented in Figure 14. This problem is simulated using both the BGK and Boltzmann models with ν = ρ for the BGK case. The results are presented in Figure 15 where density, temperature and the velocity are plotted. The top line presents BGK model results while the bottom line presents Boltzmann ones. The same scale is used to ease the comparison. The velocity fluctuation (δu x , δu y ) is represented with the same scaling as to observe that the vortex is more dissipated by the Boltzmann model than by BGK one. In Figure 16 we replot the initial density along with the final BGK and Boltzmann results using the same color scale used for showing the initial data. From the Figures we can observe that both models reproduce a vortex at the correct final location but with different rates of dissipation. The BGK model furnishes results closer to the compressible Euler solution results which means over-relaxation.

We analyze now the performances of the scheme. In order to obtain the previous results both models used 1400 time steps to reach t final . This result is due to the choice of the time integrator: the same for both models. The Open MP version of the code is run on 8 cores and the total CPU times on a laptop HP ZBook Intel Core i7-4940MX CPU @ 3.10GHz×8 operated by Ubuntu 15.10 64bits, are of the order CPU BGK 102s, and CPU Boltz 2008s. We split the code in several distinct and conceptually important parts which are denoted as Transport, ToConservative, ToPrimitive, and Collision. The first routine implements the transport phase, the last routine implement the collision phase, while the two routines in the middle reconstruct the conserved macroscopic and primitive variables from the kinetic distribution. These two routines are necessary for defining the Maxwellian distribution for performing the collision step in the BGK model, while they only serve to show the results for the Boltzmann model. Moreover, since the macroscopic quantities change only due to the transport phase, one can think to associate the cost of these routines to the transport part. All routines are monitored during a simulation of the vortex problem using the same mesh reported above. The results are reported in BGK model to 97% for Boltzmann model. Indeed, while for the BGK model, the collision routine and the routine recomputing the conservative variables have some impact on the total CPU time, for the Boltzmann model the collision part is the only one participating to the global cost at least for the mesh used. We finally gather in Table 6 the results in terms of computational costs when the number of spatial cell in each direction doubles from M = 25 2 to M = 200 2 keeping the number of velocity cell fixed to N = 32 × 32 on velocity domain [-7.5; 7.5] 2 . We monitor the number of cycles, the CPU time T in second(s), and the time per cycle is obtained by T cycle = T/N cycle , the time per cycle per cell by T cell = T/N cycle /M and the time per cycle per degree of freedom T dof = T/N cycle /(M d × N dv ). Doubling the number of spatial cells leads to an increase in CPU time by a factor 7 to 8 for both models and the number of time cycles is exactly doubled. Boltzmann model demands on average 20 times more CPU resources than the BGK one. The time per cycle per cell T cell is rather constant (7.5 × 10 -6 for BGK and 1.4 × 10 -4 for Boltzmann (thanks to the almost linear complexity of the fast spectral solver), it becomes relatively easy to estimate the cost of further refined simulations keeping the velocity mesh fixed.

Test 3.2. Re-entry test in two dimensions with changing angle of attack in time.

This test is inspired from re-entry test cases described in [START_REF] Prasad | Flow simulation over re-entry bodies at supersonic and hypersonic speeds[END_REF][START_REF] Dimarco | A multiscale fast semi-Lagrangian method for rarefied gas dynamics[END_REF]. The computational domain is set to Ω = [0; 4] × [0; 4]. Within this domain we initiate three static objects, two small rectangles upfront ([x 0 ; x 1 ] × [y 0 ; y 1 ] and [x 0 ; x 1 ] × [y 0 ; y 1 ]) and one larger one behind ([x 0 ; x 1 ] × [(y 0 + y 0 )/2; (y 1 + y 1 )/2]) where x 0 = 1.5, x 1 = 1.7, x 0 = 1.8, x 1 = 2, y 0 = 1.7, y 1 = 1.95, y 0 = 2.05, y 1 = 2.3. The computational mesh in physical space is made of 800 × 800 square cells. The velocity space is [-10, 10] and is discretized with 32 points. We report the solution obtained in the fluid regime, i.e. τ = 0, obtained by projecting the distribution f over the equilibrium state M after each transport phase and the results obtained with τ = 10 -2 using both BGK and Boltzmann models. The initial density is set to ρ(t = 0) = 1, the velocity (u x , u y )(t = 0) = (3, 0) and the temperature to T (t = 0) = 1 everywhere. The final time is set to t final = 10. On the boundaries with the objects reflective boundary conditions are employed. Inflow boundary conditions are imposed to the west boundary whereas outflow boundary conditions are set elsewhere. The inflow boundary conditions are evolving in time and equal to

(u x , u y ) BC (t) =    (3, 0) if 0 ≤ t ≤ t 1 ( 9 -g 2 (t), g(t)) if t 1 ≤ t ≤ t 2 ( 3 √ 2 2 , 3 √ 2 2 ) if t 2 ≤ t ≤ t final (56)
where t 1 = 3/2, g(t) = t -t 1 and t 2 = 3 √ 2/2 + t 1 . Given these initial data, we expect a detached shock wave to occur upfront the objects and some complex wave pattern behind them. Moreover, setting the inflow boundary conditions to (56) splits the simulation into three stages. The first stage consists in the inflow boundary conditions facing the objects up to t 1 . For this stage the upfront detached shock and the complex flow structure behind the objects are formed but they are not yet steady. Next, for the second stage, the inflow boundary condition is changing its direction by smoothly increasing the y component of (u x , u y ) BC (t) up to t = t 2 . Note that this mimics a modification of the angle of attack of the objects with time. This change modifies the entire flow structure. Last, for the third stage the inflow boundary condition is fixed to (u x , u y ) BC = (u 0 , v 0 ) up to the final time. As such the flow reaches an almost steady state.

In Figure 17 we present the results for nine intermediate times when τ = 0 for the density. These results correspond to the fluid limit model that could be obtained when solving compressible Euler equations. Next, in Figures 18 and19, we compare the BGK results (top row) with the Boltzmann ones (middle row) for intermediate iterations 4500, 12000 and 20000 for τ = 10 -2 for respectively the density and the temperature. In the bottom row, it is shown the differences between the two models in terms of density and temperature. The color legend is the same for the first two rows for both Figures and it can be found in Figure 17, while for the bottom rows they are shown in the pictures.

Initial density BGK model results

Boltzmann model results In order to conclude this part, let us present some performance data related to those simulations. The total amount of CPU time needed to compute the N cycle = 26000 cycles for the BGK solution is 14.5h, while for Boltzmann model is 304h= 12.67d. The ratio is of the order 21 in favor of BGK consistently with the previous simulations. However, even if Boltzmann results demand a large amount of CPU time, we have seen that discrepancies do exist with respect to BGK model and in some cases, especially far from equilibrium, they cannot be ignored.

Part 4. Numerical results for the space three dimensional case.

In this last part, we present one numerical test in which we compare the relaxation model with the Boltzmann model in three space and velocity dimensions in a kinetic regime. Solving the full Boltzmann equation in three dimensions is extremely resource consuming, even if the fast spectral methods is used and consequently shared memory systems are not sufficient for this kind of problems. For this reason, only for this last case we adapted the method to distributed memory systems by employing MPI architecture as already stated and described in Section 4. The results reported for this situation are not to be intended as optimal since we are adapting the scheme to this kind of architecture and improvements in terms of efficiency are attended in the next future.

Test 4.1. Three dimensional re-entry test case.

The computational domain is set to Ω = [0, 2] 3 with a static cuboid placed in the center (see Fig. 20). The velocity space is [-10, 10] 3 and discretized with 32 3 points. The relaxation parameter is set to τ = 0.3. The initial density ρ is set to 1, the temperature T = 1 while the initial velocity is given by (u x , u y , u z ) = (2, 0, 0). The final time is set to t final = 0.6 leading to 379 time steps. The inflow boundary conditions are imposed on the left boundary (x = 0) while outflow boundary conditions on the remaining part of the boundary are imposed. Hard sphere molecules are considered for Boltzmann while for the BGK model the frequency ν is taken equal to µ = C α 4π(2λπ) α . For both models the CFL condition considered is consequently given by ∆t ≤ min ∆x

|v max | , τ µ . ( 57
)
The results are shown for the temperature and the density in Figure 20 while the discrepancies between the two in Figure 21. From the analysis of such results it clearly emerge a difference in the profiles of the macroscopic quantities between the two models. We now analyze the performances. This test case was run on the EOS supercomputer at CALMIP, Toulouse France (https://www.calmip.univ-toulouse.fr/). The supercomputer is equipped with 612 computational nodes, each of them containing two Intel R Ivybridge 2.8GHz 10 core CPUs and 64 GB of RAM. Each CPU is equipped with 25MB of cache memory. The code was compiled with gcc-5.3.0 and executed on 90 computational nodes. That is to say, on 1800 computational cores in parallel. In the case of the Boltzmann collision operator the Fast Fourier Transforms were computed by means of the fftw library, version 3.3.4. The total runtime was equal to t = 93713s (26h) for the Boltzmann model. This is equivalent to 46000 computational hours (3.25 years) on a sequential machine. On the other hand, the runtime for the BGK model was only t = 2174s (0.6h) with a ratio 43 in favor of the simpler relaxation model. The profiling data are summarized in Table 7. The communications and MPI synchronization take 38% of the computational time for BGK equation. On the other hand, for the Boltzmann model the time spent on communications and synchronization amounts merely to 3.4% of the total runtime. This is due to extreme computational complexity of the three dimensional Boltzmann collision kernel. As discussed in the previous paragraphs, the costs relatives to the routines ToConservative, ToPrimitive which are relative to the computation of the macroscopic variables from the distribution function can be imputed to the transport phase, since this is the only part of the scheme which modifies these quantities, since collisions preserves density, momentum and energy in the cell.

BGK-Boltzmann

Conclusion and perspectives

In this paper we have generalized the Fast Kinetic Scheme [START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part I: Basics on the BGK equation[END_REF][START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part II: The high order case[END_REF][START_REF] Dimarco | Towards an ultra efficient kinetic scheme. Part III: High-performance-computing[END_REF][START_REF] Dimarco | A new class of high order semi-Lagrangian schemes for rarefied gas dynamics[END_REF] to the challenging case of the Boltzmann collision integral. We have shown that it is possible to solve the full unsteady three dimensional Boltzmann equation with variable hard sphere kernel in a reasonable amount of time by using parallel architectures. Up to the author knowledge, the results reported in this paper represent one of the very first attempts of solving the seven dimensional Boltzmann model with deterministic numerical schemes. This has been made possible by combining a fast semi-Lagrangian approach for the transport part with a fast spectral method for the collision dynamic.

We have performed several numerical tests with the aim of detailing the behavior of the method in different situations in order to understand its strengths and weaknesses. A side scope of the paper has been to show the differences that arise between the Boltzmann model and the popular BGK relaxation model.

Differences have been observed to be large far from equilibrium situations.

In the future, we aim in working in the direction of additionally improving the fast spectral method since we have observed it to be one of the main bottleneck of the scheme by using for instance different grids during the transport and the collision phases. Another fundamental direction we aim to pursue is the optimization of the MPI parallelization which is necessary for considering more realistic applications as

Figure 1 :

 1 Figure 1: Illustration of the transport scheme for the FKS scheme. Left panel before transport step, right panel after transport step. Each discrete velocity (index k) drives its own transport equation with velocity v k . The representation of f is made by means of a piecewise constant function. The shape of the entire function has not changed during the transport but the cell-centered values (bullets) may have.

Figure 2 :

 2 Figure 2: Domain decomposition for the MPI parallelization. The domain slices are distributed over computational nodes. After the collision step the new information brought by the particles on the cells located on the subdomain boundary are communicated to the ghost cells of the neighboring nodes.

Figure 3 :

 3 Figure 3: Test 1.1. Sanity check -Left: f (vx, vy = 0, t) (symbols and straight line) as a function of velocity vx and time for t = 1, 2 and t final = 10 versus the exact solution (dashed line). Right: L 1 error as a function of time for 8, 16 and 32 points (or modes) in each direction. Maxwellian molecules in dimension two.

Figure 4 :

 4 Figure 4: Test 1.2. Differences between the distribution functions obtained with the BGK and the Boltzmann models. Maxwellian molecules in two dimensions. The distribution function values for the red marked iterations are depicted in Figure 5.

16

 16

Figure 5 :

 5 Figure 5: Test 1.2. Differences in the distribution functions obtained by using the BGK and the Boltzmann models at t = 0.2, t = 1, t = 2 and t = 6.4 from top-left to bottom-right. Maxwellian molecules in two dimensions.

Figure 6 :

 6 Figure 6: Test 1.3. Sanity check -Initial distribution function for f (vx, vy, 0) middle, f (vx, vy, 0.25) left and f (vx, vy, -0.25) right. Hard sphere molecules.

Figure 7 :

 7 Figure 7: Test 1.3. Sanity check for the hard sphere case at time t = 0, 0.25, 0.5 and t final = 2 for 32 points in each direction. The Figure shows the relaxation to the Maxwellian state for the distribution function f (vx, vy, vz, t) on a plane which passes from the points ±(-1, -1, -0.25) and of normal (-0.1, -0.1, 1).

Figure 8 :

 8 Figure 8: Test 1.3. Differences between the distribution functions obtained with the BGK and the Boltzmann models. Hard sphere molecules in three dimensions. The distribution function values for the red marked iterations are depicted in Figure 9.

Figure 9 :

 9 Figure 9: Test 1.3. Differences in the distribution functions obtained by using the BGK and the Boltzmann models at t = 0.05, t = 0.25, t = 0.5 and t = 1 from top-left to bottom-right for vz = 0. Hard sphere molecules in three dimensions.

Figure 10 :

 10 Figure 10: Test 2.1. One dimension in space and two dimension in velocity Boltzmann model with Maxwellian molecules for a Sod like test case at t final = 0.15. Mesh convergence results for τ = 10 -3 (left) and τ = 10 -4 (right). Density (top), velocity (middle) and temperature (bottom) are shown for M = 50, 100, 200, and 400 cells and N = 64 2 velocity cells.

Figure 11 :Figure 12 :

 1112 Figure 11: Test 2.2. One dimension in space and two dimension in velocity Boltzmann model with Maxwellian molecules (blue) and a BGK model (red) for a Sod like test case at t final = 0.15 for τ = 10 -3 (left) and τ = 10 -4 (right). Density (top), velocity (middle top), temperature (middle bottom) and heat flux (bottom) are shown for M = 800 cells and N = 64 2 velocity cells.

Figure 13 :

 13 Figure 13: Test 2.3: One dimension in space and three dimension in velocity for Sod like test case at t final = 0.5 for the hard sphere molecules case. Mesh convergence results for τ = 10 -2 (left) and τ = 2 10 -3 (right). Density (top), velocity (middle) and temperature (bottom) are shown for M = 32, 64 and 128 cells and N = 32 3 velocity cells.

Figure 15 :

 15 Figure 15: Test 3.1. Two dimensional vortex test case for τ = 10 -1 with M = 100 × 100 spatial cells and N = 32 2 velocity cells. Top/bottom line: results for BGK/Boltzmann models at t final = 10.0. Left to right: density, temperature, velocity vector (δux, δuy).

Figure 16 :

 16 Figure 16: Test 3.1. Two dimensional vortex test case for τ = 10 -1 with M = 100 × 100 spatial cells and N = 32 2 velocity cells. Left: initial density. Middle/right: results for BGK/Boltzmann model at t final = 10.0 using the same color scale of the initial density.

Figure 17 :

 17 Figure 17: Test 3.2. Two dimensional re-entry test case for τ = 0 with M = 800 × 800 spatial cells and N = 32 2 velocity cells. Density profile. Top-left to bottom-right iterations 500, 1400, 4500, 8000, 10000, 12000, 14000, 16000, 20000.

Figure 18 :

 18 Figure 18: Test 3.2. Two dimensional re-entry test case for τ = 10 -2 with M = 800 × 800 spatial cells and N = 32 2 velocity cells. BGK model (top row) Boltzmann model (middle row), difference between the two models (bottom row) at iterations 4500, 12000 and 20000. The color legend for the first two rows can be found in figure 17. Density profile.

Figure 19 :

 19 Figure 19: Test 3.2. Two dimensional re-entry test case for τ = 10 -2 with M = 800 × 800 spatial cells and N = 32 2 velocity cells. BGK model (top row) Boltzmann model (middle row), difference between the two models (bottom row) at iterations 4500, 12000 and 20000. The color legend for the first two rows can be found in figure 17. Temperature profile.

Figure 20 :

 20 Figure 20: Test 4.1. Three dimensional re-entry test case for τ = 0.3 with M = 90 × 90 × 90 spatial cells and N = 32 3 velocity cells. BGK model (left column), Boltzmann model (right column) at time t final = 0.6. Top row: temperature field with velocity streamlines, bottom row: isosurfaces of the density. Middle row: isosurfaces of the density difference between BGK and Boltzmann models.

Table 6 :

 6 Test 3.1. Two dimensional vortex problem simulations with spatial mesh variation. Monitoring of CPU time. Time per cycle is obtained by T cycle = T/N cycle , time per cycle per cell by T cell = T/N cycle /M and time per cycle per degree of freedom T dof = T/N cycle /(M d × N dv ).

Figure 21 :

 21 Figure 21: Test 4.1. Three dimensional re-entry test case for τ = 0.3 with M = 90 × 90 × 90 spatial cells and N = 32 3 velocity cells. Discrepancies between the BGK solution and the Boltzmann solution at time t final = 0.6. Density on the left and temperature on the right.

Table 1 :

 1 Schematic table of the average cost for each routine of the code on a cylindrical Sod like problem in 2D×2D. Simulation is run on a serial machine. The collision part demands almost all the CPU resources and, within this step, the calls to the (inverse) Fourier transform routine covers more than 85% of the total cost of one single step.

	Transport	2 %	-	2 %
	Update moments	0 %	-	0 %
			Compute F(f )	5 %
	Collision	≥ 98 %	θ p -Loop: compute coeff. θ p -Loop: compute F -1	5 % 87 %
			Update Q	3 %
	Total	100 %		100%

Table 2 :

 2 Test 1.1. Sanity check -L 1 and L 2 relative errors for the fast spectral method for different number of points at t final = 10. Maxwellian molecules in dimension two.

  10 -2 4.41 × 10 -5 1.08 × 10 -8 16.4 × 10 5 10 -1 6.33 × 10 -4 4.42 × 10 -8 8.2 × 10 5 400 × 64 2 443 113.495 2.56 × 10 -1 6.40 × 10 -4 1.56 × 10 -7 16.4 × 10 5

	100 × 64 2 4.1 × 10 5 200 × 64 2 8.2 × 10 5 400 × 64 2 1.77 × Boltzmann 111 0.935 8.42 × 10 -3 8.42 × 10 -5 2.06 × 10 -8 222 2.738 1.23 × 10 -2 6.17 × 10 -5 1.51 × 10 -8 15, 15] 436 7.822 64 2 100 × 64 2 111 8.174 7.36 × 10 -2 7.36 × 10 -4 5.64 × 10 -8 4.1 × 10 5 [-15, 15] 200 × 64 2 218 28.095 1.27 × Sod like Riemann problem in 1D×3D τ = 10 -2
	Model	Velocity N Vel.	Cell # M × N	Cycle N cycle	Time T (s)	T/cycle T cycle (s)	T/cell T cell (s)	T/d.o.f T dof (s)
	BGK	32	[-16, 16]	32 × 32 2 1.05 × 10 6 64 × 32 2 2.10 × 10 6 128 × 32 2 4.19 × 10 6	395 395 395	12.64 19.55 40.75	3.20 × 10 -2 1.00 × 10 -3 3.05 × 10 -8 4.95 × 10 -2 7.73 × 10 -4 2.36 × 10 -8 1.03 × 10 -1 8.06 × 10 -4 2.46 × 10 -8
	Boltzmann	32	[-16, 16]	32 × 32 2 1.05 × 10 6 64 × 32 2 2.10 × 10 6 128 × 32 2 4.19 × 10 6	395 395 395	2104.17 ∼ 35mn 4799.19 ∼ 1.3h 9114.56 ∼ 2.5h	5.33 12.10 23.10	1.66 × 10 -1 5.08 × 10 -6 1.90 × 10 -1 5.79 × 10 -6 1.80 × 10 -1 5.50 × 10 -6

Table 3 :

 3 One dimensional in space, two dimensional in velocity Maxwellian molecules and three dimensional in velocity hard sphere molecules simulations. Comparisons between the BGK and Boltzmann models for spatial mesh variation. Monitoring of CPU time. Time per cycle is obtained by T cycle = T/N cycle , time per cycle per cell by T cell = T/N cycle /M and time per cycle per degree of freedom T dof = T/N cycle /(M × N ).

	1										1
								Boltzmann 32 cells									Boltzmann 32 cells
								Boltzmann 64 cells									Boltzmann 64 cells
	0.9							Boltzmann 128 cells	0.9								Boltzmann 128 cells
	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.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	2										2
		Boltzmann 32 cells										Boltzmann 32 cells
		Boltzmann 64 cells										Boltzmann 64 cells
	1.8	Boltzmann 128 cells								1.8		Boltzmann 128 cells
	1.6										1.6
	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
	0										0
	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2
	6.5										6.5
		Boltzmann 32 cells										Boltzmann 32 cells
		Boltzmann 64 cells										Boltzmann 64 cells
		Boltzmann 128 cells										Boltzmann 128 cells
	6										6
	5.5										5.5
	5										5
	4.5										4.5
	4										4
	3.5										3.5
	3										3
	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2	0	0.2	0.4	0.6	0.8	1	1.2	1.4	1.6	1.8	2

Table 4 :

 4 One dimensional in space, two dimensional in velocity Maxwellian molecules and three dimensional in velocity hard sphere molecules simulations. Simulations are performed using the OpenMP version of the scheme run on 8 hreads. Comparisons between the BGK and Boltzmann models for velocity mesh variation. Monitoring of CPU time. Time per cycle is obtained by T cycle = T/N cycle , time per cycle per cell by T cell = T/N cycle /M and time per cycle per degree of freedomT dof = T/N cycle /(M × N ).should concentrate to improve the efficiency in the future. For all reported simulations Maxwellian molecules are considered for the Boltzmann model. The CFL condition employed is as for the one dimensional case and for all tests the following

	-3

Table 5 .

 5 The Boltzmann model is 20 times more expensive than BGK one. More in details, the relative cost of the collision steps jumps from 70% for Figure 14: Test 3.1. Two dimensional vortex test case for τ = 10 -1 with M = 100 × 100 spatial cells and N = 32 2 velocity cells. Initial density, temperature and velocity vector (δux, δuy).

																					10
																					8
	1										1
	0.99									1	0.99									1
	0.97 0.98									0.97 0.98 0.99	0.97 0.98									0.97 0.99 0.98	6
	0.96									0.96	0.96									0.96
										0.95										0.95
	0.95									0.94	0.95									0.94
	0.94									0.93	0.94									0.93	4
	0.93										0.93
																					2
										10										10
	0	2	4	6	8	10 0	2	4	6	8	0	2	4	6	8	10 0	2	4	6	8	0	0	2	4	6	8	10
										Cycle CPU Main routines Cost CPU		Cost
											(s)									vs total (s)	vs total (%)
															Transport			0.36			0.35%
								BGK		1400	102			ToConservative ToPrimitive Collision		30.26 0.19 70.83		29.77% 0.19% 69.69%
																=				101.64		100%
								Boltzmann		1400	2008			Transport ToConservative ToPrimitive Collision =		1.06 66.75 0.29 1939.63 2007.73		0.053% 3.325% 0.015% 96.61% 100%

Table 5 :

 5 Profiling of the average cost for each routine of the code on the test 3.1. (two dimensional vortex problem) simulated on a laptop using M = 100 2 and N = 32 2 points.

  2

							Density										Temperature								Velocity
																												10
	BGK results	0.96 0.98 1 1.02 1.04	0	2	4	6	8	10 0	2	4	6	8	10	0.95 0.96 0.97 0.98 0.99 1 1.01 1.02 1.03 1.04 1.05	0.96 0.98 1 1.02 1.04	0	2	4	6	8	10 0	2	4	6	8	10	1 1.005 1.01 1.015 1.02 1.025 1.03 1.035 1.04 1.045 1.05	8 2 4 6 0	0	2	4	6	8	10
	Boltzmann results	0.96 0.98 1 1.02 1.04	0	2	4	6	8	10 0	2	4	6	8	10	0.975 0.98 0.985 0.99 0.995 1 1.005 1.01 1.015 1.02	0.96 0.98 1 1.02 1.04	0	2	4	6	8	10 0	2	4	6	8	10	1.012 1.014 1.016 1.018 1.02 1.022 1.024 1.026	2 4 6 8 10 0	0	2	4	6	8	10

  × 10 -4 1.47 × 10 -7 = 64 × 10 4 50 2 × 32 2 701 245.03 0.350 1.40 × 10 -4 1.37 × 10 -7 = 256 × 10 4 100 2 × 32 2 1400 2008.56 1.435 1.44 × 10 -4 1.40 × 10 -7 = 1024 × 10 4 200 2 × 32 2 2800 15762 5.630 1.41 × 10 -4 1.37 × 10 -7 = 4096 × 10 4

				Vortex problem in two dimension
	Model	Velocity N Vel.	Cell # M d × N dv	Cycle N cycle	Time T (s)	T/cycle T cycle (s)	T/cell T cell (s)	T/d.o.f T dof (s)
				25 2 × 32 2 = 64 × 10 4	351	2.61	0.0035	5.60 × 10 -6 5.47 × 10 -9
	BGK	32 2	[-7.5, 7.5]	50 2 × 32 2 = 256 × 10 4 100 2 × 32 2 = 1024 × 10 4 200 2 × 32 2 = 4096 × 10 4	701 1400 2800	13.77 102.42 785.54	0.0196 0.0732 0.2806	7.84 × 10 -6 7.66 × 10 -9 7.32 × 10 -6 7.15 × 10 -9 7.02 × 10 -6 6.85 × 10 -9
				25 2 × 32 2	351	32.92	0.0938	1.50
	Boltzmann	32 2	[-7.5, 7.5]

Table 7 :

 7 Profiling of the average cost for each routine of the code on the test 4.1. (three dimensional re-entry test case) simulated using M = 90 3 and N = 32 3 points.

This plane passes through the center of the two spheres defined by[START_REF] Shu | Essentially non-oscillatory and weighted essentially non-oscillatory schemes for hyperbolic conservation laws[END_REF]

Acknowledgments

The authors would like to thanks Professor Francis Filbet from the University of Toulouse III, Professor Rémi Abgrall and Dott. Tulin Kaman from the institute of mathematics in Zürich, Switzerland for their suggestions and help. This work has been supported by the Galileo project G14 (Fast Asymptotic-Preserving and semi-Lagrangian schemes for High Performance Computing : applications to plasmas) from the Franco-Italian University and by the ANR project MOONRISE (MOdels, Oscillations and NumeRIcal SchEmes, 2015-2019). Thomas Rey was partially funded by Labex CEMPI (ANR-11-LABX-0007-01). This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation 2016-P1542 and the authors acknowledge the help from CALMIP.

well as the development of techniques for treating complex boundaries. Finally, the extension of the present scheme to plasmas is under study.