A higher order moment preserving reduction scheme for the Stochastic Weighted Particle Method
aa r X i v : . [ m a t h . NA ] J a n A HIGHER ORDER MOMENT PRESERVING REDUCTIONSCHEME FOR THE STOCHASTIC WEIGHTED PARTICLEMETHOD ∗ SONAM LAMA † , JOHN ZWECK † , AND
MATTHEW GOECKNER ‡ Abstract.
The Stochastic Weighted Particle Method (SWPM) is a Monte Carlo techniquedeveloped by Rjasanow and Wagner that generalizes Bird’s Direct Simulation Monte Carlo (DSMC)method for solving the Boltzmann equation. To reduce computational cost due to the gradualincrease in the number of stochastic particles in the SWPM, Rjasanow and Wagner proposed severalparticle reduction schemes designed to preserve specified moments of the velocity distribution. Here,we introduce an improved particle reduction scheme that preserves all moments of the velocitydistribution up to the second order, as well as the raw and central heat flux both within each groupof particles to be reduced and for the entire system. Furthermore, we demonstrate that with the newreduction scheme the scalar fourth-order moment can be computed more accurately at a reducedcomputational cost.
Key words.
Boltzmann equation, stochastic weighted particle method, deterministic particlereduction, higher order moments
AMS subject classifications.
1. Introduction.
A fundamental problem in the computational modeling ofrarefied gases and plasmas is to determine the velocity probability density function(pdf) of each particle species. The evolution of these velocity pdfs is governed bythe Boltzmann equation, which models particle transport and collision processes[3, 4, 5]. Both deterministic and stochastic particle methods are used to solve theBoltzmann equation. Although deterministic methods avoid the uncertainties inher-ent in stochastic approaches, the cost of computing the Boltzmann collision opera-tor can still be prohibitively high, especially in the low probability tails of the pdf.However, recent theoretical advances in combination with increased computationalpower have led to the introduction of several promising deterministic spectral meth-ods [7, 8, 9, 10, 15, 16, 24]. For example, Gamba and Rjasanow recently proposeda Petrov-Gelerkin method whose computational efficiency is comparable to that ofstochastic methods [10]. Despite the recent reduction of their computational cost, de-terministic methods are not as flexible as stochastic methods for the modeling of thediverse range of collision, transport, and boundary surface phenomena, and particlegain and loss mechanisms that occur in experimental settings [11, 13].Accurate modeling of the low probability tails of the velocity distribution is also ofinterest to experimentalists. For example, reaction rates in plasmas are determined bythe overlap between the electron velocity pdf and the electron-impact cross sectionsof the various species. Therefore, accurate calculation of the low probability tailsof the electron velocity pdf is critical. If the plasma is in thermal equilibrium, theelectron velocity pdf can often be assumed to be Maxwellian. However, experimentalresults demonstrate that the Maxwellian assumption is often invalid [1, 6, 21, 22, 23],especially for pulsed plasmas where the velocity pdf may depend strongly on bothspatial position and on time [17]. Consequently, there is still a pressing need for ∗ Submitted to the editors DATE. † Department of Mathematics, The University of Texas at Dallas, Richardson, TX([email protected], [email protected]). ‡ Department of Physics, The University of Texas at Dallas, Richardson, TX ([email protected]). 1
S. LAMA, J. ZWECK, AND M. GOECKNER improved stochastic particle methods that have both greater computational efficiencyand higher accuracy, especially in the higher order moments and in the low probabilitytails of the distributions.Unlike deterministic methods, particle methods are not based on solving theBoltzmann equation directly. Rather they simulate a real system using stochasticparticles, each of which represents a group of physical particles that are in close prox-imity in phase space. Collisions between stochastic particles are designed so as toapproximate the collision processes modeled by the Boltzmann equation. Modernparticle methods are based on the Direct Simulation Monte Carlo (DSMC) methodwhich was developed by Bird [2]. A convergence proof for this method was given byWagner [26]. The DSMC method has many computational advantages over deter-ministic methods. However, the computational cost of accurately computing the lowprobability tails is still very high. To resolve the low probability tails with relativelylow computational cost, Rjasanow and Wagner introduced a generalization of theDSMC method which they called the Stochastic Weighted Particle Method (SWPM).One of the challenges for the SWPM is that the number of stochastic particlesgradually increases over the course of the simulation. To reduce the computationalcost, Rjasanow and Wagner proposed to use a particle reduction scheme in combi-nation with a clustering technique [18, 20]. With these methods, the particles arepartitioned into groups such that the particles are close together, and each group isreplaced by a small number of particles. A reduction scheme that does not requireclustering was proposed by Vikhansky and Kraft [25]. Their reduction scheme redis-tributes the statistical weights of the particles so as to conserve the mass, momentumand energy of the ensemble. They argue that the efficiency of a particle reductionscheme that relies on clustering primarily depends on the computational cost of theclustering algorithm. For the clustering algorithm used for the results in this paper,the computational cost scales linearly with the number of computational particles.In this context, it is important to note that the convergence theorem for the SWPMobtained by Rjasanow and Wagner [20] requires that the maximum diameter of thegroups of particles converges to zero as the initial number of computational parti-cles increases. As a result, there is a theoretical advantage to employing a clusteringtechnique.The reduction schemes proposed by Rjasanow and Wagner were designed to pre-serve a specified set of moments of the distribution. The particle reduction schemeof Rjasanow and Wagner that preserves the most moments is a deterministic reduc-tion scheme that preserves the total weight, momentum, energy and central heat fluxwithin each group [18, 20]. The total weight corresponds to the fraction of physicalparticles represented by the group. Although this reduction scheme preserves the cen-tral heat flux of each group, it does not preserve the raw heat flux, and consequentlyneither the raw nor the central heat flux are preserved for the entire system.In this paper, we improve upon the reduction scheme of Rjasanow and Wagnerby conserving all of the moments up to the second order (i.e. the full pressure andmomentum flux tensors), as well as both the raw and central heat flux, which are thirdorder moments. Conservation of all these moments within each group automaticallyguarantees that they are conserved for the entire system.We performed two series of simulation studies to evaluate the degree to whichour new deterministic particle reduction scheme improves upon that of Rjasanow andWagner’s deterministic reduction schemes. First, we present results which confirmthat the existing scheme of Rjasanow and Wagner does not conserve the raw heatflux within each group, while our new scheme conserves both raw and central heat
TOCHASTIC WEIGHTED PARTICLE METHOD
2. The stochastic weighted particle method.
In this section, we review thestochastic weighted particle method for the spatially homogeneous Boltzmann equa-tion. The stochastic weighted particle method is a particle method [19] that improvesupon Bird’s DSMC method [2] by decreasing the uncertainty in the computation ofrare events. In Bird’s method, each stochastic particle represents the same number ofphysical particles, and the number of stochastic particles is kept constant throughoutthe simulation. With the SWPM, the number of physical particles represented bya single stochastic particle varies over the course of the simulation. Each stochasticparticle represents a group of physical particles that are in close proximity in phasespace. Each stochastic particle is characterized by its velocity and weight. The weightquantifies the proportion of physical particles represented by the given stochastic par-ticle. The SWPM is based on a generalized version of the collision process used in theDSMC method in which only the physical particles corresponding to some portion ofthe weights of the colliding stochastic particles undergo collisions. For each stochasticcollision, this results in the creation of two new stochastic particles whose velocitiesare given by the post-collision velocities and whose weights quantify the proportion ofphysical particles involved in the collision process [18, 20]. The weights of the originalpair of colliding stochastic particles are reduced so as to keep the total weight of thesystem constant. As the number of stochastic particles increases due to collisions, thenumber of high velocity particles in the low probability tails of the velocity pdf alsoincreases. By periodically applying a clustering technique and a particle reductionscheme, the proportion of particles in the center of the distribution is reduced. Thecombined effect of these processes is to increase the fraction of stochastic particlesoccupying the low probability tails of the velocity pdf, which decreases the statisticaluncertainty in the tails.We consider the spatially homogeneous Boltzmann equation for a single species ofparticles with unit mass. This equation, which describes the evolution of the velocityprobability density function (pdf), f , due to collisions, is given by(2.1) ∂f∂t ( v , t ) = Z R Z S B ( v , w , Θ ) (cid:2) ( f ( v ′ , t ) f ( w ′ , t ) − f ( v , t ) f ( w , t ) (cid:3) d Θ d w , with an initial condition of the form(2.2) f ( v ) = f ( v , . Here, S denotes the unit sphere, B is the collision kernel, t is time, v and w are thepre-collision velocities, and v ′ and w ′ are the post-collision velocities. For simplicity, S. LAMA, J. ZWECK, AND M. GOECKNER for the results in this paper we consider isotropic Maxwell type interactions, for which(2.3) B ( v , w , Θ ) = 14 π . Assuming that the collisions are elastic, the post-collision velocities are given in termsof the pre-collision velocities and the direction vector, Θ , by(2.4) v ′ = 12 (cid:2) v + w − Θ | w − v | (cid:3) and w ′ = 12 (cid:2) v + w + Θ | w − v | (cid:3) . The state of the i th stochastic particle is given by ( v i , g i ), where v i and g i arethe velocity and weight, respectively. The state of the entire stochastic system is(2.5) z = { ( g , v ) , ( g , v ) , . . . , ( g m , v m ) } , where m is the current number of stochastic particles. To model a collision betweenthe stochastic particles indexed by i and j , we introduce the weight transfer function, γ coll ( z ; i, j ). This function encodes the proportion of physical particles representedby the stochastic particles indexed by i and j that undergo collisions when the stateof the system is z . The weight transfer function cannot exceed the minimum of theweights of the colliding particles,(2.6) 0 ≤ γ coll ( z ; i, j ) ≤ min( g i , g j ) . During a collision between the stochastic particles indexed by i and j , only thefraction of physical particles in the system represented by the weight γ coll ( z ; i, j ) un-dergo collisions. This process is modeled by adding one or two new stochastic particlesto the system. For the results in this paper, we use γ coll ( z ; i, j ) = min ( g i , g j ), whichalways results in two new stochastic particles. In this case, the state, [ J coll ( z ; i, j, Θ )] k ,of the k -th stochastic particle after a collision between particles i and j is given by[18, 20],(2.7) [ J coll ( z ; i, j, Θ )] k = ( v k , g k ) , if k ≤ m, k / ∈ { i, j } , ( v i , g i − γ coll ( z ; i, j )) , if k = i, ( v j , g j − γ coll ( z ; i, j )) , if k = j, ( v ′ i , γ coll ( z ; i, j )) , if k = m + 1 , ( v ′ j , γ coll ( z ; i, j )) , if k = m + 2 , resulting in a new system state,(2.8) z = { ( g , v ) , ( g , v ) , . . . , ( g m +1 , v m +1 ) , ( g m +2 , v m +2 ) } . After the collision, the fraction of physical particles corresponding to the weight γ coll ( z ; i, j ) are assigned the post-collision velocities, and the remaining fraction ofparticles is unchanged. The two new stochastic particles are indexed by m + 1 and m + 2. To keep the total weight constant this weight is subtracted from the weightsof the colliding stochastic particles, indexed by i and j . For elastic collisions givenby (2.4), this stochastic collision process conserves the total weight, momentum andenergy.To correctly model the evolution of the velocity pdf, we must relate the collisionfrequency for the stochastic system to that of the physical system. The total collisionfrequency in the physical system is given by(2.9) ν = Z R Z R Z S B ( v , w , Θ ) f ( v , t ) f ( w , t ) d Θ d w d v . TOCHASTIC WEIGHTED PARTICLE METHOD ν g i g j denote the frequency of collisions between the physical particles thatcorrespond to stochastic particles with states ( v i , g i ) and ( v j , g j ), then by (2.9) weobtain(2.10) ν g i g j = g i g j Z S B ( v i , v j , Θ ) d Θ . Furthermore, if we let e ν ij denote the frequency of collisions between particles i and j in the stochastic system, then by the definition of the weight transfer function, wehave that(2.11) e ν ij γ coll ( z ; i, j ) = ν g i g j . Therefore by (2.10), we obtain(2.12) e ν ij = g i g j γ coll ( z ; i, j ) Z S B ( v i , v j , Θ ) d Θ , and so, by (2.9) and (2.12), the total collision frequency in the stochastic system isgiven by(2.13) e ν ( z ) = 12 m X i =1 m X j =1 j = i g i g j γ coll ( z ; i, j ) Z Θ ∈ S B ( v i , v j , Θ ) d Θ . Using this frequency, we can obtain the waiting time between stochastic colli-sions. Since it is memoryless, this waiting time is a Poisson process which followsan exponential distribution. Therefore, the probability that a collision did not occurby time, t , is given by the survival function, P ( s > t ) = e − e ν ( z ) t . Since the survivalfunction has an uniform distribution on [0 , t = − ln( r ) / e ν ( z ), where r is a random number uniformly distributed on [0 , p ( z ; i, j ), of a collision between the stochastic particles i and j isgiven by the ratio of the frequency e ν ij of collisions between stochastic particles i and j given in (2.12) and the total collision frequency given in (2.13), that is,(2.14) p ( z ; k, l ) = g k g l γ coll ( z ; k,l ) Z Θ ∈ S B ( v k , v l , Θ ) d Θ m X i =1 m X j =1 j = i g i g j γ coll ( z ; i, j ) Z Θ ∈ S B ( v i , v j , Θ ) d Θ . Once a pair of colliding stochastic particles, k and l , has been randomly selected, thedirection vector Θ is chosen using the probability density function(2.15) η ( Θ ) = B ( v k , v l , Θ ) Z e Θ ∈ S B ( v k , v l , e Θ ) d e Θ . In the case of the constant collision kernel given by (2.3), the probability in (2.14)further simplifies to(2.16) p ( z ; k, l ) = g k g l γ coll ( z ; k,l ) m X i =1 m X j =1 j = i g i g j γ coll ( z ; i, j ) , S. LAMA, J. ZWECK, AND M. GOECKNER and η ( Θ ) = π . After the colliding pair of stochastic particles and the direction vectorhave been chosen, the velocities and weights of the colliding particles are updated using(2.7).The computation of the collision frequency in (2.13) and collision probability in(2.14) can be computationally expensive since in general these quantities need to beupdated after each collision. This issue also arises for the DSMC method when thecollision kernel is not constant. To overcome this computational issue, the techniqueof null collisions was developed by Koura [12] for the DSMC method. With thistechnique, an equal maximum collision frequency is assigned to all pairs of particles,which leads to an equal probability of collision for all pairs. Consequently the collidingpair can be selected at random from a uniform distribution. Once a pair is chosen,we decide whether the collision is an actual one or a null collision based on theprobability given by the ratio between the actual collision frequency and the assignedequal collision frequency. Rjasanow and Wagner generalized the technique of the nullcollisions to the SWPM [18, 20].
3. A reduction scheme conserving total weight, momentum, pressuretensor and heat flux.
As we explained in section 1, one of the challenges for theSWPM is that the number of stochastic particles gradually increases. For computa-tional feasibility, it is necessary to periodically reduce the number of particles. Thereare two steps in the reduction process. First, the stochastic particles need to be clus-tered into groups, and then each group of particles needs to be replaced by a smallnumber of particles.A number of clustering techniques have been proposed by Rjasanow, Wagnerand their collaborators [14, 18, 20]. One of these techniques is based on partitioningparticles into two groups with a cutting plane whose normal vector is in the directionof the eigenvector corresponding to the largest eigenvalue of the covariance matrix ofthe particles [18, 20]. This partitioning method is performed iteratively on each of thepartitioned groups using the group’s covariance matrix. The iteration continues untilthe product of the total weight and the standard deviation of the particle speeds withineach group is minimized, which results in a roughly uniform number of stochasticparticles in each group. We use this clustering method for the results in this paper.Rjasanow and Wagner also proposed several stochastic and deterministic particlereduction schemes to replace each group by a group with a small number of particles.These schemes are based on conserving a specific set of moments of the distributionwithin each group. The details of these reduction schemes can be found in [18, 20].We are interested in deterministic reduction schemes that conserve as many mo-ments as possible, so that the structure of the velocity pdf is preserved. In this paper,we propose a particle reduction scheme that conserves all the moments of the velocitypdf up to second order, given in Table 1, together with the raw and central heat flux,which are the most physically relevant third order moments. The raw heat flux vector, h , is computed relative to the origin, while the central heat flux, q , is relative to thedrift velocity, V . They are given by(3.1) h = 12 m X i =1 g i v i | v i | and q = 12 m X i =1 g i ( v i − V ) | v i − V | . In the following discussion, when we refer to third order moments we simply meanthe raw and central heat flux.The reduction scheme of Rjasanow and Wagner that preserves the most momentsand is closest to our scheme is the one that preserves the total weight, momentum,
TOCHASTIC WEIGHTED PARTICLE METHOD Table 1
Moments of the velocity pdf
Moment order Raw Moment Symbol Central Moment SymbolZero Total Weight ̺ — —First Momentum ̺ V — —Second Momentum Flux Tensor Π Pressure Tensor P energy, and central heat flux [18, 20]. With this scheme, although the central heat fluxis conserved within each group, the momentum flux tensor and the pressure tensorare not. Only the total energy, which is the trace of the momentum flux tensor, isconserved. As a consequence, the raw heat flux for each group is not conserved, andtherefore, the raw and central heat flux of the entire system are also not conserved.To conserve both the raw and central moments of a group, it is necessary andsufficient to conserve either of these moments and all of the lower order moments.Because of the additive property of raw moments, if a raw moment is conserved withineach group then it must also be conserved for the entire system. In particular, since thetotal weight and momentum are raw moments, conservation of ̺ and ̺ V within eachgroup ensures that these two moments are conserved for the entire system. Therefore,if we could conserve the total weight, momentum, pressure tensor and central heatflux within each group, then all of the raw and central moments up to the secondorder together with the raw and central heat flux would be conserved for the entiresystem.We formalize this idea as follows. Let m be the number of stochastic particles inthe system, and suppose that the particles have been partitioned into b n groups with m l stochastic particles in the l -th group, G l . Let g l,i and v l,i denote the weight andvelocity of the i -th particle in the l -th group. Then, the total weight, ̺ l , momentum, ̺ l V l , momentum flux tensor, Π l , and raw heat flux, h l , for the l -th group are givenby(3.2) ̺ l = m l X i =1 g l,i , ̺ l V l = m l X i =1 g l,i v l,i , Π l = m l X i =1 g l,i v l,i v Tl,i , h l = 12 m l X i =1 g l,i v l,i | v l,i | , where V l is the drift velocity of the l -th group. The pressure tensor, P l , and thecentral heat flux, q l , of the l -th group are given by(3.3) P l = m l X i =1 g l,i ( v l,i − V l ) ( v l,i − V l ) T and q l = 12 m l X i =1 g l,i (cid:0) v l,i − V l (cid:1) (cid:12)(cid:12) v l,i − V l (cid:12)(cid:12) . The energy, E l , and temperature, T l , are given by(3.4) E l = m l X i =1 g l,i | v l,i | , and 3 ̺ l T l = m l X i =1 g l,i | v l,i − V l | , where the quantity on the right hand side of the formula for T l is the trace of thepressure tensor. The energy is given in terms of the temperature by(3.5) E l = ̺ l | V l | + 3 ̺ l T l . S. LAMA, J. ZWECK, AND M. GOECKNER
The raw moments of the entire system are given by(3.6) ̺ = ˆ n X l =1 ̺ l , ̺ V = ˆ n X l =1 ̺ l V l , Π = ˆ n X l =1 Π l , and h = ˆ n X l =1 h l . Here ̺ is the total weight, V is the drift velocity, Π is the momentum flux tensor and h is the raw heat flux of the entire system.The relationship between the raw and central second order moments is given by(3.7) P l = Π l − ̺ l V l V Tl . Using this relationship, we observe that for a reduction scheme to preserve both of thesecond order moments, P l and Π l , it is sufficient to conserve the total weight, ̺ l , themomentum, ̺ l V l , and either P l or Π l . Since the raw moments are additive (see (3.6)),conservation of the total weight, momentum and momentum flux tensor within eachgroup leads to the conservation of the these moments for the entire system. Using (3.7)for the entire system, we conclude that the pressure tensor for the entire system is alsoconserved. Therefore, a reduction scheme that conserves the total weight, momentumand either of the second order moments for each group leads to the conservation ofthe moments up to second order for the entire system.Similarly, the relationship between the raw and central third order moment canbe determined using (3.2), (3.4), and (3.7) giving the equation,(3.8) q l = h l − P l V l − ̺ l V l (cid:12)(cid:12) V l (cid:12)(cid:12) − ̺ l T l V l , which relates the raw and central heat flux to each other via the lower order moments.As above, to conserve both the raw and central moments of a group up to third orderit is sufficient to conserve the total weight, momentum, pressure tensor, and either ofthe third order moments of the group. Similarly, using the additivity property of theraw moments, and the relationships between the moments given by (3.7) and (3.8)for the entire system, the pressure tensor and central heat flux of the system are alsoconserved together with all the raw moments. This verifies our claim that to conservethe raw and central moments of the system up to the third order during a reductionprocess, it is sufficient to conserve the total weight, momentum, pressure tensor andcentral heat flux of each group.Next, we present a novel particle reduction scheme that conserves the total weight,momentum, pressure tensor, and central heat flux in a group. First, we outline the ideabehind the conservation of these moments. Before describing this scheme, we brieflyrecall that for the reduction scheme that preserves the total weight and momentum of agroup, we simply replace all the stochastic particles in the group by a single stochasticparticle with the given weight and momentum [20]. The next higher order momentsare the momentum flux tensor and the pressure tensor. Since, the pressure tensoris a 3 × TOCHASTIC WEIGHTED PARTICLE METHOD
Theorem
Let G l be a group of stochastic particles. Suppose that the pres-sure tensor, P l , has k non-zero eigenvalues, λ , . . . , λ k , for some k ∈ { , , } , andan associated orthonormal set of eigenvectors, Θ , . . . , Θ k , with the direction of Θ i chosen so that b q l,i = Θ Ti q l > . Let e G l be the reduced group of k stochastic particleswhose weights and velocities, ( e v i , e g i ) , for i = 1 , . . . , k are given by (3.9) e v i = V l + γ i s k λ i ̺ l Θ i , e g i = ̺ l k
11 + γ i , e v i + k = V l − γ i s k λ i ̺ l Θ i , e g i + k = ̺ l k γ i γ i , for i = 1 , . . . , k where, (3.10) γ i = √ ̺ l b q l,i √ k λ i + s ̺ l b q l,i k λ i , for i = 1 , . . . , k. Then, e G l preserves the total weight, momentum, pressure tensor, and central heat fluxof G l , which leads to the preservation of all the moments up to the second order aswell as the raw and central heat flux of G l .Proof. We consider the case where the eigenvalues of the pressure tensor are allnonzero. We let the reduced group, e G l , consist of three pairs of particles, with eachpair of the form,(3.11) e v i = V l + α i Θ i , e v i +3 = V l − α i +3 Θ i , and e g i + e g i +3 = ̺ l , i ∈ { , , } , for some α i ∈ R and Θ i ∈ S . We derive the conditions on the unknown parameters, e g i , α i , and Θ i , so as to conserve the total weight, momentum, pressure tensor andheat flux.0 S. LAMA, J. ZWECK, AND M. GOECKNER
By construction, the total weight is conserved,(3.12) X i =1 e g i = ̺ l . Similarly, if we impose the condition(3.13) e g i α i = e g i +3 α i +3 , for i ∈ { , , } , we find that the momentum of the group is conserved, since(3.14) X i =1 e g i e v i = X i =1 e g i ( V l + α i Θ i ) + e g i +3 ( V l − α i +3 Θ i ) = ̺ l V l . Next, to conserve the pressure tensor, P l , we use the fact that it is a 3 × D = diag[ λ , λ , λ ] and an orthonormal matrix Q = [ Θ , Θ , Θ ] such that(3.15) D = Q T P l Q. That is, each { λ i , Θ i } is an eigenpair of the matrix P l . The condition, e P l = P l ,that the reduction scheme preserves the pressure tensor is therefore equivalent to thecondition(3.16) D = Q T h X i =1 e g i (cid:0)e v i − V l (cid:1)(cid:0)e v i − V l (cid:1) T i Q = X i =1 (cid:0)e g i α i + e g i +3 α i +3 (cid:1) (cid:0) Q T Θ i (cid:1) (cid:0) Q T Θ i (cid:1) T . Therefore, to conserve the pressure tensor, we require that(3.17) e g i α i + e g i +3 α i +3 = λ i , for i ∈ { , , } . In the basis of eigenvectors, the central heat flux, b q l , is given by(3.18) b q l = Q T q l . As in the statement of the theorem, we choose the direction of Θ i so that the i -thcomponent, b q l,i = Θ Ti q l , of b q l is positive. To conserve the central heat flux in thenew basis, we have(3.19) b q l = 12 X i =1 (cid:0)e g i α i − e g i +3 α i +3 (cid:1) Q T Θ i = 12 X i =1 (cid:0)e g i α i − e g i +3 α i +3 (cid:1) e i , and we obtain(3.20) b q l,i = 12 (cid:2)e g i α i − e g i +3 α i +3 (cid:3) . Next, to solve for α i and e g i , we apply a technique used by Rjasanow and Wagner[18, 20] . We introduce a new parameter, γ i , and express α i as(3.21) α i = γ i s λ i ̺ l . TOCHASTIC WEIGHTED PARTICLE METHOD λ i = e g i e g i +3 α i he g i +3 + e g i i . Substituting the expression for α i given by (3.21), and using (3.11) for the sum ofweights, we obtain(3.23) e g i e g i +3 γ i = 1 . Using this relationship in (3.11), we obtain(3.24) e g i = ̺ l γ i , α i = γ i s λ i ̺ l , and e g i +3 = ̺ l γ i γ i , α i +3 = 1 γ i s λ i ̺ l . To determine γ i , we substitute (3.24) into (3.20) to obtain(3.25) e g i α i − e g i +3 α i +3 = r ̺ l λ i γ i (cid:0) γ i − (cid:1) . Now by (3.20), b q l,i = (cid:2)e g i α i − e g i +3 α i +3 (cid:3) >
0. Therefore γ i > b q l,i = 12 r ̺ l λ i γ i (cid:0) γ i − (cid:1) = ⇒ γ i − √ ̺ l e q i √ λ i γ i − . Solving for the positive root, we obtain(3.27) γ i = √ ̺ l b q l,i √ λ i + s ̺ l b q l,i λ i . Therefore, the post reduction particles are given by (3.9) and (3.10) as required.If the pressure tensor has at least one zero eigenvalue, the moments can be con-served with fewer than six particles. The reason is that there is no need to introduceparticles whose heat flux is in the direction of the eigenvectors corresponding to thezero eigenvalues. In this situation, the result follows similarly to the calculationsabove.
4. Theoretical convergence of SWPM with the new reduction scheme.
In this section, we show that our new reduction scheme satisfies the assumptionsin Wagner’s convergence theorem for the SWPM [20, Thm. (3.22)]. This theoremprovides a collection of assumptions which guarantee that the sequence of empiricalmeasures of the Markov process produced by the SWPM converges to the weak so-lution of the Boltzmann equation as n → ∞ . These assumptions on the reductionscheme are given by [20, eq. (3.162)], and [20, eq. (3.164)]. According to Rjasanowand Wagner, assumption [20, eq. (3.162)] assures that the reduction is sufficientlyprecise, and [20, eq. (3.164)] restricts the increase in energy during reduction. Sincethe energy is conserved in our new reduction scheme, the second assumption relatedto the energy is satisfied. For assumption [20, eq. (3.162)], the arguments given byRjasanaw and Wagner for their deterministic reduction schemes also apply to our newdeterministic reduction scheme. Therefore, for this assumption to hold for our new2 S. LAMA, J. ZWECK, AND M. GOECKNER reduction scheme, it is sufficient to show that the inequality given by [20, eq. (3.273)]holds. This inequality states that(4.1) (cid:12)(cid:12)(cid:12)(cid:12) Z Z l Φ( e z l ) p red ( z l ; d e z l ) − Φ( z l ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ || ϕ || L (cid:20) m l X i =1 g i (cid:12)(cid:12) V l − v i (cid:12)(cid:12) + ̺ l p T l (cid:21) . Here,(4.2) z l = { ( g , v ) , ( g , v ) , . . . , ( g m l , v m l ) } is the state of a group G l prior to reduction, e z l is the post-reduction state, and p red ( z l ; d e z l ) is a measure that gives the probability that the post-reduction state liein the volume element, d e z l . The function Φ, which approximates the velocity pdf, isgiven by(4.3) Φ( z l ) = m l X i =1 g i ϕ ( v i ) , for the particles in the group G l . Here ϕ is an arbitrary test function. The norm forthe test function, || ϕ || L , is defined as(4.4) || ϕ || L = max (cid:26) || ϕ || ∞ , sup v = w ∈ R | ϕ ( v ) − ϕ ( w ) || v − w | (cid:27) . In the inequality (4.1), Z Z l Φ( e z l ) p red ( z l ; d e z l ) gives the expectation of Φ for the reducedsystem.For our new deterministic reduction scheme, in the spatially homogeneous case,for each group only one state is possible after reduction. Therefore, in the case whereall three eigenvalues of the pressure tensor are positive, p red ( z l ; d e z l ) = δ J red( zl ) ( d e z l ),where [ J red( z l ) ] i = ( e v i ( z l ) , e g i ( z l )), for i = 1 , . . . ,
6, is the post reduction state givenby Theorem 3.1. Therefore,(4.5) Z Z l Φ( e z l ) p red ( z l ; d e z l ) = Φ( J red ( e z l )) = X j =1 e g j ϕ ( e v j ) . Since X j =1 e g j = ̺ l , and applying the triangle inequality, we obtain(4.6) (cid:12)(cid:12)(cid:12)(cid:12) Z Z l Φ( e z l ) p red ( z l ; d e z l ) − Φ( z ) (cid:12)(cid:12)(cid:12)(cid:12) = (cid:12)(cid:12)(cid:12)(cid:12) X j =1 e g j ϕ ( e v j ) − m l X i =1 g i ϕ ( v i ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ X j =1 (cid:12)(cid:12)(cid:12)(cid:12) e g j ̺ l m l X i =1 g i ϕ ( e v j ) − e g j ̺ l m l X i =1 g i ϕ ( v i ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ X j =1 e g j ̺ l m l X i =1 g i (cid:12)(cid:12)(cid:12)(cid:12) ϕ ( e v j ) − ϕ ( v i ) (cid:12)(cid:12)(cid:12)(cid:12) ≤ || ϕ || L X j =1 e g j ̺ l m l X i =1 g i (cid:12)(cid:12)(cid:12)(cid:12)e v j − v i (cid:12)(cid:12)(cid:12)(cid:12) , TOCHASTIC WEIGHTED PARTICLE METHOD m l X i =1 g i = ̺ l , we obtain(4.7) X j =1 e g j ̺ l m l X i =1 g i (cid:12)(cid:12)(cid:12)(cid:12)e v j − v i (cid:12)(cid:12)(cid:12)(cid:12) ≤ m l X i =1 g i (cid:12)(cid:12) V l − v i (cid:12)(cid:12) + X j =1 e g j α j = m l X i =1 g i (cid:12)(cid:12) V l − v i (cid:12)(cid:12) + (cid:18) X j =1 e g j α j + 2 X j =1 X k>j e g j e g k α j α k (cid:19) ≤ m l X i =1 g i (cid:12)(cid:12) V l − v i (cid:12)(cid:12) + (cid:18) X j =1 e g j α j + X j =1 X k>j e g j e g k ( α j + α k ) (cid:19) = m l X i =1 g i (cid:12)(cid:12) V l − v i (cid:12)(cid:12) + (cid:18) ̺ l X j =1 e g j α j (cid:19) = m l X i =1 g i (cid:12)(cid:12) V l − v i (cid:12)(cid:12) + ̺ l p T l . Here, the final equality is obtained from (3.4) and (3.11). Therefore, the desired in-equality (4.1) holds for our new reduction scheme, and Wagner’s convergence theoremapplies in this context.
5. Numerical results.
In this section, we discuss our numerical results. Thealgorithm was implemented in C++ and all simulations were performed on a desktopmachine with a 3.6 GHz single processor. We verified that the total times for theparticle collisions and for the clustering and particle reductions both scale linearlywith the initial number of computational particles, m . The time taken to simulatethe clustering and particle reductions was approximately four times larger than thetime taken to simulate the particle collisions. However, as we will show in Table 3,for the results in Figures 1 to 3 below, the total computational time is only about 30seconds for N = 100 ensembles with m = 10 ,
240 particles per ensemble.First, to numerically verify the conclusions of Theorem 3.1, we study the sumover all the groups of the reduction errors for the raw and central heat flux. For thisstudy, we consider an initial Maxwellian distribution with temperature, T = 1, anddrift velocity, V = h , , i . We used a single ensemble to obtain these results, andthe initial number of computational particles is chosen to be m = 10,240. Once thenumber of computational particles reaches 4 m , we reduce it to e m ≈ m , which wasthe strategy that produced the largest errors for the deterministic reduction schemesof Rjasanow and Wagner [18]. We chose this strategy to demonstrate that our methodperforms well even under this condition.In Table 2, we compare three reduction schemes. All three schemes conservethe total weight and momentum, and in addition to these moments, the reductionschemes conserve the moments associated to their names. The first two schemes,energy conservation, and energy and central heat flux conservation (Ct. HF), arethe reduction schemes of Rjasanow and Wagner, and the third one is our reductionscheme which conserves the pressure tensor (PT) and central heat flux (Ct. HF). To4 S. LAMA, J. ZWECK, AND M. GOECKNER
Table 2
Maximum average relative errors in the central and raw heat flux for the three reduction schemes.
Reduction Scheme Central HeatFlux Error Raw Heat FluxErrorEnergy 1 0.01743Energy and Central Heat Flux (Ct. HF) 3.00844e-15 0.015072Pressure Tensor (PT) and Central Heat Flux (Ct. HF) 2.12406e-15 6.81888e-16 compare these schemes, we compute the relative 2-norm errors for each of the thirdorder moments of each group and take their average over all the groups, that is we let(5.1) E = 1 X l =1 || m After ,l − m Before ,l || || m Before ,l || . We obtained these average relative errors for the first ten reductions, and show themaximum of these errors in Table 2. The errors for the pressure tensor and centralheat flux scheme are smaller than 10 − , which is negligible. However, in the thirdcolumn of the table, we observe that for the energy and central heat flux scheme theraw heat flux error is about 2 × times larger than that for the pressure tensorand central heat flux scheme. These results support the theory in section 3 that theenergy and central heat flux scheme only conserves the central heat flux in each group,and does not conserve the raw heat flux, while the pressure tensor and central heatflux scheme conserves both. Furthermore, the energy scheme has the largest error forboth third order moments. Since this scheme replaces a group by two particles withequal weights and opposite velocities relative to the drift velocity, the central heatflux of the group after reduction is zero. This observation explains why the relativeerror in the central heat flux is 1 for the energy reduction scheme.In [20], Rjasanow and Wagner observed that the higher order moments of a dis-tribution are conserved statistically when averaged over a large number of ensembles,even if the reduction scheme only conserves the lower order moments. However, theyfound that the existing deterministic reduction schemes require a larger initial num-ber of computational particles for the convergence of the scalar fourth order momentthan for the lower moments. To examine this, for each reduction scheme we studiedthe convergence of (1,1)-component of the momentum flux tensor, Π , , the secondcomponent of the raw heat flux, h , and the scalar fourth order moment,(5.2) s = m X i =0 g i | v i | , as we increase m . For this study, we chose the initial condition to be a mixture ofMaxwellian distributions, since for this pdf there is an analytical formula for the givenmoments as a function of time, t [20]. The initial distribution is given by(5.3) f ( v ) = α M V ,T ( v ) + (1 − α ) M V ,T ( v ) , where M V ,T ( v ) and M V ,T ( v ) are Maxwellian distributions with drift velocities V and V , and temperatures T and T , respectively. We chose α = 0 . V = h− , , i , V = h , , i , and T = T = 1. We performed two sets of simulations inwhich we studied the short term (transient) behavior of s in the time interval [0 , TOCHASTIC WEIGHTED PARTICLE METHOD .
9% confidence interval as a function of time in the interval [0 , m is given by(5.4) E = | m anal − m || m anal | , where m = 1 N N X i =1 m i is the average of the simulated moments over the N ensembles.Similarly, the half-width of the relative confidence interval is given by(5.5) CI = z (1 − α ) | m anal | r σ N , where σ = P Ni =1 (cid:0) m i − m (cid:1) N − z (1 − α ) is the z -score for the confidence interval with α = 10 − . A statistical simulationcomputes a moment accurately if E < CI , that is if there is a high probability that m anal lies in the confidence interval centered at m , and that this confidence interval isrelatively narrow.In the first set of simulations, we used N = 500 ensembles and various initialnumbers of computational particles, m . In the left column of Figure 1, we showthe relative error, E , in (5.4) and confidence interval, CI , in (5.5) at time t = 3 forthe (1 , , , (top), the second component of theraw heat flux, h , (middle), and the scalar fourth-order moment, s , (bottom). Thepercentage relative error, for the energy, the energy and central heat flux, and thepressure tensor and central heat flux reduction schemes are shown using the symbolsin the legends. The half-width of the relative confidence intervals are shown using thecorresponding vertical lines. These quantities are plotted for the different values of m , which is displayed using a logarithmic scale. For each value of m , we have offsetthe results for the three reduction schemes from each other to aid comprehension.First, we observe that for each moment, the confidence intervals primarily dependon m rather than on the reduction scheme. Furthermore, with one slight exception,the errors for Π , and h are within the confidence intervals, even for a small numberof computational particles. For Π , , this result is to be expected since all threereductions schemes are designed to conserve momentum. However, as we saw inTable 2, the energy and energy and central heat flux reduction schemes do not preservethe raw heat flux. Therefore, the accuracy of the computation of h with these twoschemes is simply due to statistical averaging over the 500 ensembles. Significantly, inmost cases the errors for Π , and h are smaller with the pressure tensor and centralheat flux scheme than with the other two reduction schemes.The main advantage to be gained from using the new pressure tensor and heatflux reduction scheme can be seen in the results for the scalar fourth-order moment, s (see the bottom left panel of Figure 1). With our method, the error in s lies withinthe confidence interval for m ≥ m = 10,240.Therefore, the energy conservation, and energy and central heat flux conservationreduction schemes require more than 10 times the initial number of computationalparticles as the pressure tensor and heat flux conservation scheme to approximate thescalar fourth-order moment with the same degree of accuracy. As we see in Table 3,this requires at least seventeen times the computational time.6 S. LAMA, J. ZWECK, AND M. GOECKNER
256 512 1024 2048 1024000.250.50.751
Energy & Ct. HFPT & Ct. HF
256 512 1024 2048 1024000.250.50.7511.25
EnergyEnergy & Ct. HFPT & Ct. HF
256 512 1024 2048 1024000.511.52
EnergyEnergy & Ct. HFPT & Ct. HF
256 512 1024 2048 1024000.250.50.75
EnergyEnergy & Ct. HFPT & Ct. HF
256 512 1024 2048 102400246
EnergyEnergy & Ct. HFPT & Ct. HF
256 512 1024 2048 102400246
EnergyEnergy & Ct. HFPT & Ct. HF
Fig. 1 . Percentage relative error, E , in (5.4) at time t = 3 for selected moments, m , of thevelocity pdf. We compare the performance of the SWPM with the three different reduction schemesshown in the legends. We show results for the (1 , -component of the momentum, Π , , (top row),the second component of the raw heat flux, h , (middle row), and the scalar fourth-order moment, s , (bottom row). We plot the errors using symbols and the half confidence intervals with verticallines, so that E < CI when the symbol lies on the line. In the left column, we plot E as a functionof the number of particles, m , per ensemble for N = 500 ensembles. In the right column, we plot E as a function of m when N is chosen so that N × m = 1 , , . To further examine how accurately the three reductions schemes compute thescalar fourth-order moment, in Figure 2 we plot the evolution of s as a function oftime, together with 99 .
9% confidence intervals. The numerical results are shown withred-dashed lines and the true values are shown with solid blue lines. The results for
TOCHASTIC WEIGHTED PARTICLE METHOD True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
Fig. 2 . Evolution of the scalar fourth order moment, s , as a function of time together with . confidence intervals. In the different rows we show the results for different initial numbersof computational particles, m , per ensemble. We used N =
500 ensembles in all the panels. Weshow the results for the energy scheme (left column), energy and central heat flux scheme (middlecolumn), and pressure tensor and central heat flux scheme (right column). S. LAMA, J. ZWECK, AND M. GOECKNER
Table 3
Total computational time for the simulation results shown in Figure . These results wereobtained with the pressure tensor and central heat flux scheme. The computational times for thetwo reduction schemes of Rjasanow and Wagner were similar ( ≈ ± ). N = 500 N × m = 1 , , m t (sec) t (sec)256 2.61 20.85512 6.06 22.021,024 13.76 24.932,048 31.26 31.2610,240 205.56 39.05 the energy and the energy and central heat flux reduction schemes are shown in theleft and middle columns. The numerical results are visually close to the true valuesonly for m = 10 ,
240 (bottom left and middle panels). On the other hand, with thepressure tensor and central heat flux scheme (right column), the numerical resultsare reasonably accurate across the entire time range for m = 1 , s is significantly fasterfor the pressure tensor and central heat flux conservation scheme, than for the othertwo reduction schemes.To summarize our conclusions so far, the new reduction scheme provides im-proved accuracy at a significantly reduced computational cost. To provide additionalevidence for this conclusion, we performed a second set of simulations where we fixedthe total number of computational particles, m × N , to be 1,024,000. This valuewas kept constant to obtain approximately equal sized confidence intervals for all thesimulations. The results for these simulations are shown in the right column of Fig-ure 1 and in Figure 3. Comparing the errors and the half width of the confidenceintervals for Π , in Figure 1 (top right panel), we observe that our scheme is accurateeven with m = 256, while the other two schemes require a larger initial number ofcomputational particles. On the other hand, for h the results obtained with all threereduction schemes are acceptable for all the values of m . In the case of the fourthorder moment, for both reduction schemes of Rjasanow and Wagner, the error lieswithin the confidence interval only for m = 10,240. On the other hand, the errorsfor our scheme lies within the confidence intervals for m ≥ , m × N , and not on m or the choice of reductionscheme. Furthermore, in Table 3 we see that as m decreases the computational timedecreases. As a consequence, for each moment, for the same level of accuracy thecomputational time for our reduction scheme is significantly less than that for theother two reduction schemes. We also observe this phenomenon in Figure 3. Forexample, we observe the same degree of accuracy in the fourth order moment for ourreduction scheme with m = 1,024 and N = 1,000 (the right panel in the third row)as for the other two methods with m = 10,240 and N = 100 (the left and middlepanels in the last row). However, the computational time of 43 seconds for the tworeduction schemes of Rjasanow and Wagner is reduced by 42% to 25 seconds for ourscheme.
6. Conclusions.
We have confirmed that the reduction scheme of Rjasanowand Wagner that conserves total weight, momentum, energy and central heat flux ofa group does not conserve the raw heat flux in each group. Consequently, the raw andcentral heat flux of the entire system are not conserved. We resolved this problem by
TOCHASTIC WEIGHTED PARTICLE METHOD True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
True ValueEnergy
True ValueEnergy & Ct. HF
True ValuePT & Ct. HF
Fig. 3 . Evolution of the scalar fourth order moment, s, as a function of time, together with . confidence intervals, for different initial numbers of computational particles, m , per ensembleand different number of ensembles, N . For these results the total number of computational particles, m × N , was kept constant. We show the results for the energy scheme (left column), energy andcentral heat flux scheme (middle column), and pressure tensor and central heat flux scheme (rightcolumn). S. LAMA, J. ZWECK, AND M. GOECKNER devising a new reduction scheme that conserves the total weight, momentum, pressuretensor and heat flux within each group. Conservation of these moments within a groupresults in the conservation of all of the moments up to the second order, and bothraw and central heat flux among the third order moments of a group. This furtherleads to the preservation of these moments for the entire system.To examine the accuracy of our new reduction scheme, we performed simulationstudies to analyze the convergence of Π , , h , and the scalar fourth order momentfor the existing and new reduction schemes. The new reduction scheme leads tothe convergence of these moments, particularly the scalar fourth order moment, withsignificantly less computational cost compared to the existing deterministic reductionschemes. This shows that the preservation of additional moments in the new reductionscheme conserves the higher moments with better accuracy, and also minimizes thereduction error.Although the conservation of higher-order moments reduces the systematic errorintroduced by the reduction process, the clustering technique must also be carefullydesigned in order to accurately and efficiently compute the low-probability tails ofthe velocity pdf. Specifically, since the tails occupy a proportionately large volumeof phase space, we need to ensure that particles in the tails that are assigned to thesame group are sufficiently close together. In a forthcoming article we will use theproof of the convergence theorem for the SWPM [20] to develop such a clusteringalgorithm. In combination with the reduction scheme introduced in this paper, wewill demonstrate that this leads to a more efficient method for the computation oftail functionals. REFERENCES[1]
J. Allen , On the applicability of the Druyvesteyn method of measuring electron energy distri-butions , Journal of Physics D: Applied Physics, 11 (1978), p. L35.[2]
G. A. Bird and J. Brady , Molecular gas dynamics and the direct simulation of gas flows ,vol. 5, Clarendon press Oxford, 1994.[3]
J. A. Bittencourt , Fundamentals of plasma physics , Springer Science & Business Media,2013.[4]
L. Boltzmann , Weitere studien ¨uber das w¨armegleichgewicht unter gasmolek¨ulen , in KinetischeTheorie II, Springer, 1970, pp. 115–225.[5]
C. Cercignani , The Boltzmann equation , in The Boltzmann equation and its applications,Springer, 1988, pp. 40–103.[6]
J. V. DiCarlo and M. J. Kushner , Solving the spatially dependent Boltzmann’s equation forthe electron-velocity distribution using flux corrected transport , Journal of Applied Physics,66 (1989), pp. 5763–5774.[7]
I. Gamba, V. Panferov, and C. Villani , Upper Maxwellian bounds for the spatially homo-geneous Boltzmann equation , Archive for Rational Mechanics and Analysis, 194 (2009),pp. 253–282.[8]
I. M. Gamba and J. R. Haack , A conservative spectral method for the Boltzmann equation withanisotropic scattering and the grazing collisions limit , Journal of Computational Physics,270 (2014), pp. 40–57.[9]
I. M. Gamba, J. R. Haack, C. D. Hauck, and J. Hu , A fast spectral method for the Boltzmanncollision operator with general collision kernels , SIAM Journal on Scientific Computing,39 (2017), pp. B658–B674.[10]
I. M. Gamba and S. Rjasanow , Galerkin–Petrov approach for the Boltzmann equation , Jour-nal of Computational Physics, 366 (2018), pp. 341–365.[11]
J. Johannes, T. Bartel, G. A. Hebner, J. Woodworth, and D. J. Economou , Directsimulation Monte Carlo of inductively coupled plasma and comparison with experiments ,Journal of the Electrochemical Society, 144 (1997), pp. 2448–2455.[12]
K. Koura , Null-collision technique in the direct-simulation Monte Carlo method , The Physicsof fluids, 29 (1986), pp. 3509–3511.[13]
M. J. Kushner , Hybrid modelling of low temperature plasmas for fundamental investigations
TOCHASTIC WEIGHTED PARTICLE METHOD and equipment design , Journal of Physics D: Applied Physics, 42 (2009), p. 194013.[14] I. Matheis and W. Wagner , Convergence of the stochastic weighted particle method for theBoltzmann equation , SIAM Journal on Scientific Computing, 24 (2003), pp. 1589–1609.[15]
C. Mouhot and L. Pareschi , Fast algorithm for computing the Boltzmann collision operator ,Math. Comp. , 75 (2006), pp. 1833–1852.[16]
L. Pareschi and B. Perthame , A Fourier spectral method for homogeneous Boltzmann equa-tions , Transport Theory Statist., 25 (2002), pp. 369–382.[17]
J. Poulose, M. Goeckner, S. Shannon, D. Coumou, and L. Overzet , Driving frequencyfluctuations in pulsed capacitively coupled plasmas , The European Physical Journal D, 71(2017), p. 242.[18]
S. Rjasanow, T. Schreiber, and W. Wagner , Reduction of the number of particles in thestochastic weighted particle method for the Boltzmann equation , Journal of ComputationalPhysics, 145 (1998), pp. 382–405.[19]
S. Rjasanow and W. Wagner , A stochastic weighted particle method for the Boltzmannequation , Journal of Computational Physics, 124 (1996), pp. 243–253.[20]
S. Rjasanow and W. Wagner , Stochastic numerics for the Boltzmann equation , Springer,2005.[21]
T. Sheridan, M. Goeckner, and J. Goree , Electron velocity distribution functions in asputtering magnetron discharge for the E × B direction , Journal of Vacuum Science &Technology A: Vacuum, Surfaces, and Films, 16 (1998), pp. 2173–2176.[22]
C. Sozzi, E. De La Luna, D. Farina, J. Fessey, L. Figini, S. Garavaglia, G. Grossetti,S. Nowak, P. Platania, A. Simonetto, et al. , Measurements of electron velocity distri-bution function , in AIP Conference Proceedings, vol. 988, AIP, 2008, pp. 73–80.[23]
W. Tan , Langmuir probe measurement of electron temperature in a Druyvesteyn electronplasma , Journal of Physics D: Applied Physics, 6 (1973), p. 1206.[24]
M.-B. Tran , Nonlinear approximation theory for the homogeneous Boltzmann equation , arXivpreprint arXiv:1305.1667, (2013).[25]
A. Vikhansky and M. Kraft , Conservative method for the reduction of the number of particlesin the Monte Carlo simulation method for kinetic equations , Journal of ComputationalPhysics, 203 (2005), pp. 371–378.[26]