Importance Sampling Variance Reduction for the Fokker-Planck Rarefied Gas Particle Method

Models and methods that are able to accurately and efficiently predict the flows of low-speed rarefied gases are in high demand, due to the increasing ability to manufacture devices at micro and nano scales. One such model and method is a Fokker-Planck approximation to the Boltzmann equation, which can be solved numerically by a stochastic particle method. The stochastic nature of this method leads to noisy estimates of the thermodynamic quantities one wishes to sample when the signal is small in comparison to the thermal velocity of the gas. Recently, Gorji et al have proposed a method which is able to greatly reduce the variance of the estimators, by creating a correlated stochastic process which acts as a control variate for the noisy estimates. However, there are potential difficulties involved when the geometry of the problem is complex, as the method requires the density to be solved for independently. Importance sampling is a variance reduction technique that has already been shown to successfully reduce the noise in direct simulation Monte Carlo calculations. In this paper we propose an importance sampling method for the Fokker-Planck stochastic particle scheme. The method requires minimal change to the original algorithm, and dramatically reduces the variance of the estimates. We test the importance sampling scheme on a homogeneous relaxation, planar Couette flow and a lid-driven-cavity flow, and find that our method is able to greatly reduce the noise of estimated quantities. Significantly, we find that as the characteristic speed of the flow decreases, the variance of the noisy estimators becomes independent of the characteristic speed.


Introduction
Recent technological advances have resulted in manufacturing processes that have made possible the production of mechanical devices that operate on the scale of micro and nanometers [1]. Such technologies include lab-on-achip devices, micro-heat exchangers, gas chromatographers and micro-jet actuators for control in aerospace. At such small scales, the Navier-Stokes-Fourier (NSF) equations are no longer able to accurately model gas flows, due to the length scales of macroscopic gradients approaching the length of the molecules mean free path, λ. This results in the existence of a region known as the Knudsen layer near solid wall boundaries where the gas is prevented from relaxing to thermodynamic-equilibrium, invalidating the assumption that locally the gas is close to thermal equilibrium required for the NSF equations to be valid.
The Boltzmann equation is a mesoscopic model that is considered to provide the most accurate description of rarefied gases beyond Newton's laws. Before the advent of such small scale technologies, rarefied gas flows' largest application area was that of supersonic atmospheric flows, where the Mach number of flow, Ma > 1. Currently, the prevalent method for numerically approximating the solution to the Boltzmann equation in such regimes is a stochastic particle method called direct simulation Monte Carlo (DSMC) [2] [3]. Due to the stochastic nature of the method, DSMC becomes very inefficient for low-speed flows. Typically the Mach number, Ma 1 for flows within micro and nano technologies, and for a given level of statistical error, the computational costs of DSMC scale as Ma −2 [4]. This results in very long computation times for such calculations, and methods which are able to efficiently solve for low speed flows are highly desirable.
Currently, there are two methods that are able to greatly reduce the variance of the desired thermodynamic outputs of DSMC calculations. The first, low-variance DMSC (LVDSMC), works by adapting the DSMC collision routine to calculate the evolution of the deviation f d = f − f M from a Maxwellian distribution f M [5]. In low speed flows the deviation from equilibrium is small, allowing for a dramatic decrease in the variance of samples. An alternative method, variance reduced DSMC (VRDSMC), is able to work without significant changes to the DSMC algorithm [6]. The method relies on importance sampling, which allows the algorithm to sample from an equilibrium distribution where the thermodynamic variables are known aproiri, to create estimators with smaller variance.
More recently an alternative method to DSMC, where the Boltzmann collision operator is approximated by a Fokker-Planck operator, has been developed and shown to be more efficient than the basic DSMC algorithm [7] [8]. Like DSMC, it is solved stochastically using notional particles that represent a certain number of real particles in the gas to be simulated, and as such, the basic algorithm suffers from the same inhibitive scaling with the Mach number. Recently Gorji et al. [9] have proposed a method to reduce the variance of the Fokker-Planck solution algorithm that relies on creating a correlated equilibrium solution using the same set of random numbers that are used in the stochastic solution of the non-equilibrium process. The parallel correlated equilibrium process is in effect a control variate for the non-equilibrium process.
This purpose of this paper is to develop an importance sampling variance reduction scheme for the Fokker-Planck method and demonstrate its effectiveness in simple test cases. The paper is organised in the following way: in section 2 we introduce the Fokker-Planck model, and numerical stochastic particle scheme which we would like to create variance reduced estimators for. We then outline the general method that allows one to create variance reduced estimators by exploiting known information about how the macroscopic fields behave at equilibrium. In section 3 we describe the variance reduction scheme proposed by Gorji el al. [9], which creates a correlated equilibrium scheme. In section 4 we propose our importance sampling scheme, which we test in section 5 on a homogeneous relaxation, Couette flow and a lid-driven cavity flow. We then compare the importance sampling method against the results obtained by using a correlated equilibrium solution.

The Fokker-Planck collision operator
The Fokker Planck collision operator has appeared in several different contexts, originally derived for the distribution function of a Brownian particle in a fluid [10], but can also be derived from an expansion of a linear Boltzmann equation, when considering the evolution of density function for a particle in a heat bath [11]. It has been used to model electrons, dense liquids and more recently has received attention for its ability to model rarefied gas flows [12]. More recently it has been extended to describe flows of monatomic gas mixtures [13], diatomic molecules [14] and has been coupled to DSMC [15]. The equation for the one-particle distribution function f (x, v, t), over a state-space comprised of the position x ∈ R 3 , velocity v ∈ R 3 and time t ∈ R + takes the form: where τ is a relaxation time, c = v − u is the local relative molecular velocity, u is the mean velocity: T is the local temperature given by where R is the specific gas constant and ρ is the local density given by The collision operator A has the property of conserving mass, momentum and energy. That is where ψ = {1, v, v 2 } is the set of collisional invariants. The advantage of having a collision operator which can be written as a Fokker-Planck equation is that there exists an equivalent stochastic differential equation (SDE) representation for the dynamics of a random variable {X t , V t } whose distribution f evolves according to (2): where W t is a 3-dimensional Wiener process. An efficient scheme for evolving a collection of representative particles with positions and velocities {X j (t), V j (t)}, j = 1 . . . N, distributed according to the distribution f (x, v, t) in time was devised by Jenny et al. [7], and can be summarised as: where i = 1, 2, 3 indexes the dimension, and ∆t is the time-step, τ is the relaxation time and ξ are standard normal distributed random variables. The spatial domain is gridded into cells, and expectations of macroscopic quantities of interest are calculated during each timestep for each computational cell. The correct viscosity is obtained by choosing the relaxation time τ = 2µ/p, where µ is the viscosity and p is the pressure.

Variance reduction for Monte Carlo sampling
In this section we outline the general framework which allows one to reduce the variance of an estimator of a particular random variable, common to control variate and importance sampling schemes. Essentially, both methods exploit information about errors in estimates of known quantities. The general principal is as follows. Suppose we have a random variable X, and we wish to estimate E[X], let our estimate of E[X] be donated byX. Let Y be a different random variable with known expectation E[Y] with an estimator denoted byŶ. Then for any α ∈ R we can use the identity to create a new unbiased estimator for E[X], The variance of this estimator is and if we minimise this over possible choices of α, then minimiser α * is given by hence the variance for this choice of α is The only condition required for the variance of the estimator to be less than the variance of the original estimator is for Cov[X,Ŷ] > 0, and soX andŶ being dependent is a necessary condition. This is all supposing that we already know, α * which presupposes that we already know Cov[X,Ŷ]. In general this is not something that is not known a priori, but can be estimated throughout the simulation. In the next sections we will see how in practice it is possible to exploit this.

Parallel process variance reduction
We now briefly describe the method proposed by [9]. The objective of the method is to create a stochastic process Z t that is able to run in parallel to the original particle scheme, where crucially, the macroscopic fields are already known. If this is performed in a manner where the parallel process is correlated with the original stochastic process then the variance of the estimators can be reduced in the way described in the previous chapter. The coupling of the new stochastic process Z t to the original SDEs (8) is achieved in the following way: where the coefficients A and D are chosen to keep the marginal distribution of (X t , Z t ), which we denote f 0 (x, z, t), as a solution to a Fokker-Planck equation which when supplemented by appropriate boundary conditions, is solved by a Maxwellian The way to choose A and D to ensure that (23) is a solution of (22) is discussed in depth in [9]. The coupling of the processes Z t and V t , by the same Wiener process W t requires that when we use particles to generate numerical solutions, the same set of random numbers is used for both the equilibrium and non-equilibrium processes. This results in the correlation of estimations of expected quantities, allowing the kind of variance reduction outlined in the previous section. In this method, the parallel equilibrium process, with distribution f M , shares the same density ρ as the non-equilibrium process described which has f as its distribution function, and so the method cannot reduce the noise in the density calculation directly. Gorji et al propose that the density ρ is found from the continuity equation using conventional finite difference methods. The results they obtain for a homogeneous relaxation, Poisuille flow and lid-driven cavity flows show the method has the ability to substantially reduce noise. Because this method uses common random numbers to reduce the variance, we will refer to this method as a common random numbers (CRN) scheme.

Importance sampling
The method we propose is an importance sampling scheme. It differs from the importance sampling scheme utilised by VRDMSC [6], because the DSMC and Fokker-Planck method account for collisions in different ways. The principle that underpins the importance sampling remains the same however. Suppose we are interested in evaluating the expectation of g(v) where v is distributed according to the distribution function f , then given N independent samples {V 1 , . . . , V N } distributed according to f the following definition gives rise to the estimate: which we know from the Central Limit Theorem, has an error of order N −1/2 . We now define a weight function which is a measure of how likely one is to see a particle with velocity v, relative to how likely one is to observe this particle if it was distributed to a reference density f ref . This definition is well defined if the distribution f is absolutely continuous with respect to f ref , meaning that f ref (S ) = 0 whenever f (S ) = 0 for any subset S of the state-space. This definition can be viewed as a Radon Nikodym derivative. It can then be observed that the expectation of g(v) with respect to the reference distribution can be estimated using the original samples: This is significant as it allows one to sample from the reference distribution f ref , using the original set of samples from the distribution f . If the reference distribution is Maxwellian, f ref = f M where one knows the thermodynamic fields analytically, then one has the ability to create variance reduced estimators as described in section 2. In order to practically apply this, we need a method of evolving the weights and velocities For VRDSMC this is possible because it can be shown directly from the Boltzmann equation, that if two particles are chosen to collide with weights W i and W j then the post collision weights must be equal to 1 2 (W i + W j ). Because the Fokker-Planck dynamics have no explicit collisions, a different way to update the weights is needed.
Importance weights can be initialised exactly, because the initial velocities of the particles are distributed according to a prescribed initial distribution f 0 . As time is evolved during the calculation, the distribution of velocities will change and hence so must the weights attached to each particle. VRDSMC is able to do this by creating collision rules that ensure that post collision velocities are still able to sample from the same reference distribution. However, these rules are not relevant for the Fokker-Planck particle dynamics as there are no explicit collisions.

Weight update rule
Instead, let us suppose that a given particle updates its velocity from V t → V t+1 , where V t is distributed according to f t and V t+1 is distributed according to f t+1 , and that we know W t = W(V t ). In order to update the weight exactly, one would need to know f t+1 (V t+1 ), however this distribution function is unknown. A simple method to estimate the updated weight is to use the zeroth order Taylor expansions of the joint distributions of V t and V t+1 , resulting in the estimate: This has approximation immediately has some desirable properties. Firstly, the error of the approximation decays with ∆t. Also, it is possible to calculate this explicitly from the update rule V t → V t+1 given by equation (9). This conditional distribution will be a gaussian centred on V t plus the deterministic drift, with a temperature dependent variance. Further to this, it has the correct conditional expectation E[ W t+1 |W t ] = W t when the distribution is stationary. However, on its own it is not a suitable choice as if such a rule is repeated the variance of this approximation diverges, which is a common problem for this type of particle weight importance sampling method [6] [16]. This is a problem, because to reduce the variance of our estimators in a meaningful way, we require the weights to be close to unity. To avoid this problem we use the same kernel density estimator approach as used by the VRDSMC method [6]. Kernel density estimation (KDE) is a method that allows one to obtain an estimatef of a density function f from samples distributed according to that density function in the following way: where K r is a kernel function that integrates over the state-space to 1, and r is a smoothing parameter that controls the width of the kernel function. We use the same spherical kernels as [6]: which returns the reciprocal of the volume of a sphere of radius r if v i lies within the sphere of radius r centred on v, and otherwise returns a zero. If we combine this with (32), the update rule that is obtained is where S r (V i ) = V j : V j − V i < r is the set of samples whose members lie within the sphere of radius of r centred on V i . This KDE step has the effect of smoothing out the variation introduced by using a conditional probabilities to estimate a marginal probability, and making the scheme more stable. Increasing the smoothing parameter r results in an estimator with a smaller variance, however it also increases the bias of the estimation, so ideally r should be chosen to be as small as possible whilst maintaining an acceptable level of variation.

Boundary Conditions
We use the same boundary condition methodology as prescribed by the VRDSMC method, that is for diffusely reflecting fully accommodating walls, with temperature T wall and tangential velocity u wall . Supposing that the Maxwellian distribution at the boundary is given by f wall (v) = ρ wall P MB (v), where P MB is a gaussian probability density, and the boundary is the plane x = 0, then the no flux boundary condition is given by and similarly for the equilibrium solution The second integrals in the above equations are the particle fluxes, and can be estimated by counting the number of particles N in the cross through a wall of area ∆s in a time period ∆t by (1/∆s∆t)N in , and at equilibrium is estimated by (1/∆s∆t) N in i W i . Also, we can use the analytical properties of the Gaussian distribution to evaluate the first integrals: Therefore, a particle that changes velocity from V to V when colliding with a wall, changes its weight according to where typically, we choose the temperature of the equilibrium wall boundary condition to be equal to the temperature of non-equilibrium wall boundary boundary condition, i.e. T wall = T wall,eq .

Homogeneous Relaxation to Equilibrium
We will demonstrate the effectiveness of this method first with a homogeneous relaxation to equilibrium, i.e. when f (t, x, v) = f (t, v) has no spatial component. We start from an initial distribution of particles which will relax towards the Maxwellian distribution f M (v, 0, (4/3)c 2 0 ). In figures (1a)-(1b) we show how the variance reduced estimator performs against the standard estimator, when estimating |v 1 | using 100 particles, with and without the KDE stabilisation procedure. In both cases, the variance of the new estimator is smaller than the standard estimator, but the estimator with stabilisation from the KDE reduces the variance even further.

Couette Flow
To test the particle weight variance reduction, we have applied the scheme to sample from a steady-state planar Couette flow, and compare to results obtained using a common random number scheme. A Couette flow is a flow where the fluid is bounded by two parallel walls moving in opposite directions within their planes, with velocity ±U wall . For Knudsen numbers Kn = 0.05, 0.5, 1.0 respectively, Figures (2), (3), (4) show the variance reduced and standard Monte Carlo estimators of the steady-state flow velocity field parallel to the wall, v 2 (x 1 ), (left) as well as the temperature profile across the channel T (x 1 ), for a Couette flow with wall velocity v wall = 0.01c 0 , Kn = 0.5, 20 cells and 100 particles per cell. All the results show a significant improvement of performance over the unweighted standard Monte Carlo estimator.
Next we compare the VRFP importance sampling scheme to the CRN correlated equilibrium scheme. Because we are interested in the noise of the estimate of the velocity profile across the channel, and the speed of the flow is small, we make the simplifying assumption that the steady-state density across the channel is constant. This allows us to choose the coefficients A = z/τ and D = √ 2RT wall /τ so that the correlated equilibrium process is distributed according to the global Maxwellian, Figure 5 compares the noise-to-signal ratio of the CRN scheme, VRFP scheme and standard Monte Carlo estimator against signal strength for the samples taken from steady state Couette flow estimator. The results show that for the CRN scheme there is a reduction in the noise-to-signal ratio by a factor of 10 over all tested levels of signal, corresponding to a speed up of over 100 times. However, it suffers the same scaling properties with signal size as the standard Monte Carlo estimator. In contrast, the importance weighted variance reduced estimator has a noise-to-signal ratio that is independent of the signal size as the signal size decreases. This results in an unbounded speed-up over the standard Monte-Carlo estimator as the signal size decreases to zero. Because of the independence of the signal strength on the noise-to-signal ratio, there is a signal strength where for larger signal strengths, the CRN scheme outperforms the particle weight scheme, and for Couette flow we estimate this to be at a Mach number greater than 0.1.

Lid-Driven Cavity
To further demonstrate the effectiveness of the method, we apply it to a lid-driven cavity flow, where the fluid is bounded in two dimensions by a square box in the (x, y) plane, with translational symmetry in the z axis. Three of the bounding walls are stationary, and one of the bounding walls moves within its plane at constant velocity U wall , giving rise to a circulatory flow within the cavity.   Figures (6a)-(6b) show the velocity and non-dimensional temperature field (T/T 0 − 1) of the steady state flow, with a lid velocity of U wall = 0.001c 0 for the standard Monte Carlo and variance reduced sampling schemes. The results have been averaged over 5000 time-steps, and 10 independent ensembles on a 50 × 50 grid, with an average of 30 particles per cell. The standard Monte Carlo scheme is not able to pick up the signal, whereas we see clearly that the importance sampling scheme is able to recover the signal. In Figure 7 we compute the streamlines of the variance reduced flow, and the sheer stress π 12 .

Conclusion
In this paper we have developed an importance sampling method for the Fokker-Planck rarefied gas model, that assigns weights to each stochastic particle allowing one to sample from an equilibrium distribution. We have demonstrated its effectiveness in reducing the variance of estimates of thermodynamics quantities for low Mach number flows over a range of Knudsen numbers. Significantly, the level of noise in the estimators becomes independent of the Mach number for low-speed flows. We believe it to be a versatile and robust method, and because it doesn't alter the basic algorithm of the particle solution scheme, it is able to be used in conjunction with other variance reduction schemes such as the CRN method.