Generation of the Maxwellian inflow distribution - Alejandro Garcia

16 Pages·2006·0.27 MB·English

JournalofComputationalPhysics217(2006)693–708
www.elsevier.com/locate/jcp
Generation of the Maxwellian inflow distribution
Alejandro L. Garcia a,*, Wolfgang Wagner b
aDepartmentofPhysics,SanJoseStateUniversity,ScienceBuilding,Room245,OneWashingtonSquare,SanJose,CA95192-0106,USA
bWeierstrassInstituteforAppliedAnalysisandStochastics,Mohrenstrasse39,10117Berlin,Germany
Received23February2005;receivedinrevisedform11January2006;accepted17January2006
Availableonline3March2006
Abstract
Thispaperpresentsseveralefficient,exactacceptance–rejectionmethodsforgeneratingtheMaxwellianinflowdistribu-
tion,thevelocitydistributionofgasmoleculescrossingaplane.Thenewmethodsaredemonstratedtobecomputationally
fasterand moreaccuratethanthe schemes commonly usedfor open boundaryconditions inparticle simulations.
(cid:2)2006Elsevier Inc. Allrights reserved.
Keywords: Randomnumbergeneration;DirectsimulationMonteCarlo;Openboundaryconditions;Acceptance–rejectionmethods
1. Introduction
Stochastic algorithms, commonly referred to as Monte Carlo methods, use random numbers generated
from a variety of distributions. Efficient generators have been developed for the most commonly used distri-
butions (e.g., uniform, Gaussian, and exponential) and general techniques (e.g., inversion) are available for
arbitrary distributions [1,2]. However, generating a complicated distribution by a generic method may not
be efficient or accurate, which is why specialized generators for specific applications are welcome.
This paper discusses the generation of random values z from the distribution
2ða(cid:2)zÞexpð(cid:2)z2Þ
p ðzÞ¼ p ; z<a; ð1:1Þ
a expð(cid:2)a2Þþa ffipffiffi½1þerfðaÞ(cid:3)
whereaissomereal-valuedparameteranderfdenotestheerrorfunction(cf.(A.12)).AsshowninSections2
and 3, this distribution arises when implementing the inflow boundary condition for particles crossing a sur-
face. Specifically, it is associated with the velocity distribution of particles, Maxwellian distributed in their
movingframeofreference,thatpassthroughaplane.Thisboundaryconditionisverycommoninmolecular
simulations of hydrodynamic flows in open systems [3,4]. With this association, we call p (z) the Maxwellian
a
inflow distribution. Note that for a=0 the distribution (1.1) reduces to the Rayleigh distribution, which is
trivial to generate [1].
* Correspondingauthor.Tel.:[phone];fax:[phone].
E-mailaddress:[email protected](A.L.Garcia).
0021-9991/$-seefrontmatter (cid:2)2006ElsevierInc.Allrightsreserved.
doi:10.1016/j.jcp.2006.01.025
694 A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708
Section 3 describes the general algorithm for generating the random velocities of particles for an inflow
boundarycondition.FortheMaxwellianinflowdistribution(1.1)therearethreemethodsforrandomvariate
generationthatareincommonuse.Thefirstisinversion(seeSection4),whichhasthedisadvantageofbeing
computationallyexpensive.Theothertwoareapproximateacceptance–rejectionschemes,discussedin[5]and
in Section 5.1, which are more efficient than inversion but not exact.
InSection5.2,wedevelopseveralexactacceptance–rejectionschemesthataremoreefficientthananyofthe
threemethodsincommonuse.ThecomputationalefficiencyofalltheschemesisdiscussedinSection6.Their
programming implementation in a practical example is summarized in Section 7. We conclude in Section 8
with some further remarks regarding applications of the present schemes and their extension to other
distributions.
2. Maxwellian inflow
This section establishes the mathematical formulation for inflow boundary conditions. Readers inter-
ested in a more physical introduction, presented in the context of a specific example, are directed to Sec-
tion 7.
The general inflow boundary condition for the Boltzmann equation [6] is
fðt;x;vÞðv;nðxÞÞ¼bðx;vÞ;
where tP0, x2oD and the velocity v2R3 is such that (v,n(x))>0. Here, oD denotes the boundary of the
spatialdomainD;n(x)istheunitinwardnormalvectoratx2oD;(v,n)isthescalar(dot)productofvectorsv
and n. The function b determines the inflow intensity (waiting time parameter)
1 Z Z
k¼ bðx;vÞdvrðdxÞ ð2:1Þ
g
oD ðv;nðxÞÞ>0
and the inflow law
1
bðx;vÞ; ð2:2Þ
gk
where r(dx) denotes the uniform surface measure (area) on oD and g is the weight of the incoming particles.
A case of special interest is Maxwellian inflow
(cid:3)M ðvÞðv;eÞ if x2C; ðv;eÞ>0;
.;V;T
bðx;vÞ¼
0 otherwise,
where
!
. kv(cid:2)Vk2
M ðvÞ¼ exp (cid:2)
.;V;T ð2pTÞ3=2 2T
is the Maxwellian distribution and
e¼nðxÞ 8x2C(cid:4)oD;
is the unit normal for some plane part C of the boundary. The inflow intensity (2.1) takes the form
1
k¼ rðCÞC; ð2:3Þ
g
where
Z
C ¼ M ðvÞðv;eÞdv. ð2:4Þ
.;V;T
ðv;eÞ>0
According to the inflow law (2.2), the position of the incoming particle is distributed uniformly on C, the in-
flowboundary,andtheparticleentersthedomainbyleavingthisboundaryatarandomtime[5,7,8].Itsveloc-
ity is generated according to the probability density
A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708 695
(
C(cid:2)1M ðvÞðv;eÞ if ðv;eÞ>0;
qðvÞ¼ .;V;T ð2:5Þ
0 otherwise.
3. General algorithm
We now turn to the question of how to generate a random variable n according to the probability density
(2.5).Employing a standard substitution, which forcompleteness ispresented inAppendixA, theproblem is
reduced to generating from the density
(
2mðaÞ(cid:2)1p(cid:2)1ðaþw Þexpð(cid:2)jjwjj2Þ if aþw >0;
~qðw ;w ;w Þ¼ 3 3
1 2 3
0 otherwise,
where
ðV;eÞ
a¼ p ð3:1Þ
ffiffiffiffiffiffi
2T
is the speed ratio and
p
mðxÞ¼expð(cid:2)x2Þþx ffipffiffi½1þerfðxÞ(cid:3); x2R. ð3:2Þ
The components of w are independent; the first two components are distributed according to the probability
density
1 exp(cid:4)(cid:2)w2(cid:2)w2(cid:5); w ;w 2R;
p 1 2 1 2
and the third is obtained as
w ¼(cid:2)z; ð3:3Þ
3
where z is distributed according to the Maxwell inflow distribution
2
p ðzÞ¼ ða(cid:2)zÞexpð(cid:2)z2Þ; z2ð(cid:2)1;aÞ. ð3:4Þ
a mðaÞ
The random variable n is obtained as (cf. (A.2), (3.3))
0 1
w(cid:5)
p 1
ffiffiffiffiffiffi
n¼V þ 2TQðeÞB w(cid:5) C, ð3:5Þ
@ 2 A
(cid:2)z(cid:5)
where Q is a coordinate rotation matrix (cf. (A.1)). The variables w(cid:5) and w(cid:5) are independent Gaussians with
1 2
zeromeanandvariance1/2.Inthecasea=0,oneobtainsfromtheinversetransformation(seeSection4)and
(A.10)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
z(cid:5) ¼(cid:2) (cid:2)logu; ð3:6Þ
whereuisuniformlydistributedon(0,1).Thegenerationofthevariablez*inthecasea6¼0,whichisthepoint
of this paper, will be discussed in the following sections.
4. Inverse transform
Onemethodtoobtainthecomponentz*tobeusedin(3.5)isbytheinversetransformmethod[1,2],thatis,
by solving numerically the equation
Z z
p ðxÞdx¼u; z<a; ð4:1Þ
a
(cid:2)1
696 A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708
where uisuniformly distributed on(0,1).Denotetheleft-hand side of Eq.(4.1)byF (z). According to(A.9)
a
(A.10), and (3.4), one obtains
1 p
F ðzÞ¼ fexpð(cid:2)z2Þþa ffipffiffi½1þerfðzÞ(cid:3)g.
a mðaÞ
Furthermore, F0ðxÞ¼p ðxÞ and
a a
2
p0ðzÞ¼ expð(cid:2)z2Þ½2z2(cid:2)2za(cid:2)1(cid:3).
a mðaÞ
Thus, the function F has an inflection point at
a
p
ffiffiffiffiffiffiffiffiffiffiffiffiffi
a(cid:2) a2þ2
zðaÞ¼ ; ð4:2Þ
2
and the function p takes its maximum there, that is,
a
maxp ðzÞ¼p ðzðaÞÞ. ð4:3Þ
a a
z<a
Note that z(a)<min(a,0).
Eq. (4.1) is solved by a Newton iteration, starting with the initial guess z =z(a). Calculate the error
0
E ¼F ðz Þ(cid:2)u; k ¼0;1;2;...;
k a k
and stop if it is small enough, specifically if |E |<E . Otherwise, calculate the new guess
k d
E mðaÞE
z ¼z (cid:2) k ¼z (cid:2) k
kþ1 k F0ðz Þ k 2ða(cid:2)z Þexpð(cid:2)z2Þ
a k k k
and continue the iteration.
The inversion method, as described above, has the disadvantage of being computationally expensive rela-
tivetoothergenerators,asdiscussedinSection6.AnalternativetoNewtoniterationistopre-computeF (z),
a
use(4.1)totabulateF(cid:2)1ðuÞandthengeneratezbyinterpolatedtablelookup.FortheMaxwellian inflowdis-
a
tribution there are several disadvantages to this approach: First, the infinite tail of the distribution must be
treated separately, typically by acceptance–rejection (see the following section). Second, the scheme is not
exactandalargetable,costingcomputermemory,withagoodinterpolationscheme,costingcomputertime,
are needed to maintain accuracy. Finally, the table of F(cid:2)1 must be recomputed for each different value of a.
a
Given the computational efficiency of the generators presented in this paper (see Section 6) it is unlikely that
inversion by table lookup, or other table methods [1], would be competitive.
5. Acceptance–rejection
This section discusses the use of the acceptance–rejection technique to generate the component z* to
be used in (3.5). This technique is based on selecting a suitable majorant (envelope) ^p such that (cf. (3.4))
p ðzÞ6^pðzÞ 8z<a. ð5:1Þ
a
The component z* is generated according to the probability density
1
Ra ^pðxÞdx^pðzÞ; z<a; ð5:2Þ
(cid:2)1
and is accepted with probability
p ðz(cid:5)Þ
a ð5:3Þ
^pðz(cid:5)Þ
so
1
Ra ^pðxÞdx ð5:4Þ
(cid:2)1
is the acceptance rate.
A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708 697
Theefficiencyoftheacceptance–rejectionmethodobviouslydependsonthechoiceoftheenvelopefunction
^p.Ideally,theacceptanceratewillbeclosetooneandtheprobabilitydensity(5.2)isadistributionthatcanbe
generated efficiently.
5.1. Approximate envelopes
This section presents two approximate acceptance–rejection methods that are in common use [5]. The
schemes are approximate because the envelopes do not satisfy (5.1) and thus the distributions they generate
only approximate p (z).
a
5.1.1. Box envelope
The first approximate method uses a rectangular box envelope (cf. (4.2), (4.3))
(cid:3)p ðzðaÞÞ; z ðaÞ<z<z ðaÞ;
^p ðzÞ¼ a < >
B
0 otherwise,
where z (a)<z (a)6a are some parameters to be specified later. Note that
< >
8
0; z<z ðaÞ;
Z z >< <
^p ðxÞdx¼ p ðzðaÞÞðz(cid:2)z ðaÞÞ; z ðaÞ<z<z ðaÞ;
B a < < >
(cid:2)1 >:p ðzðaÞÞðz ðaÞ(cid:2)z ðaÞÞ; z>z ðaÞ.
a > < >
One generates z* as
z(cid:5) ¼z ðaÞþðz ðaÞ(cid:2)z ðaÞÞu
< > <
which is accepted with probability (cf. (5.3))
p ðz(cid:5)Þ a(cid:2)z(cid:5)
a ¼ expðzðaÞ2(cid:2)ðz(cid:5)Þ2Þ.
p ðzðaÞÞ a(cid:2)zðaÞ
a
This method is approximate since the condition (cf. (5.1))
p ðzÞ6^p ðzÞ
a B
is not satisfied for z<z (a) and for z>z (a). The actual distribution generated by the box envelope is
< >
(cid:3)p ðzÞ=m ðaÞ; z ðaÞ<z<z ðaÞ;
~p ðzÞ¼ a B < >
a
0 otherwise,
where
Z z>
m ðaÞ¼ p ðzÞdz.
B a
z<
5.1.2. Reservoir envelope
The second approximate method uses the envelope
2
^p ðzÞ¼ ða(cid:2)z ðaÞÞexpð(cid:2)z2Þ;
R mðaÞ <
where z (a)<a is some parameter to be specified later. One generates z* according to the probability
<
density
1
p expð(cid:2)z2Þ ð5:5Þ
ffiffiffi
p
698 A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708
so that z* is a Gaussian with zero mean and variance 1/2. The generated value is rejected if z*>a, accepted
with probability (cf. (5.3))
p ðz(cid:5)Þ a(cid:2)z(cid:5)
a ¼ ð5:6Þ
^p ðz(cid:5)Þ a(cid:2)z ðaÞ
R <
if z*>z (a), and accepted with probability one otherwise.
<
The reservoir method is so named because it has the following physical interpretation: A particle is gener-
ated in a reservoir with position x and velocity v. The position is chosen uniformly from the interval ((cid:2)L,0)
andthevelocitychosenasv=a(cid:2)zwherezisdistributedas(5.5).Theparticleisacceptedifitmovespastthe
origin during a time interval s, that is,
xþvs>0 or ða(cid:2)zÞs>uL.
Taking
L
s¼
a(cid:2)z ðaÞ
<
gives us that z is accepted with probability (5.6).
Note that this method is approximate since the condition (cf. (5.1))
p ðzÞ6^p ðzÞ
a R
is not satisfied for z<z (a). The actual distribution generated by the reservoir envelope is
<
1 (cid:3)p ðzÞ; z ðaÞ<z<a;
~p ðzÞ¼ a <
a m ðaÞ ^p ðzÞ otherwise,
R R
where
Z a Z z<
m ðaÞ¼ p ðzÞdzþ ^p ðzÞdz.
R a R
z< (cid:2)1
5.1.3. Truncation errors
The error in these approximate methods is most easily seen from the absolute fractional error in the
moments, that is,
(cid:6)(cid:6)hzii(cid:2)h~zii(cid:6)(cid:6)
(cid:6) (cid:6);
(cid:6) hzii (cid:6)
where
Z a Z a
hzii¼ zip ðzÞdz and h~zii¼ zi~p ðzÞdz.
a a
(cid:2)1 (cid:2)1
This error is shown in Fig. 1 for the box and reservoir envelopes using typical values for the parameters (see
[5,7]) in these approximate envelopes. For the higher moments, the error is large (a few percent) when a<0.
This error can be reduced by suitable choice of the envelopes’ parameters but at the price of computational
efficiency (see Section 6).
5.2. Exact envelopes
This section presents the main results of this paper, namely, four exact envelope functions, two for each
caseofnegativeandpositivespeedratio(a<0anda>0).Welatershow(seeSection6)thattheseenvelopes
yieldefficientacceptance–rejectionschemesfortheMaxwellianinflowdistribution.Ofthefourenvelopespre-
sentedinthissection,thefirstandthethirdaremostefficientforlow-speedratioðjaj<1Þwhiletheothertwo
2
are efficient for a broad range of a.
A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708 699
100
10-1 <z3>
α = 1
1
α = 2
1
s
ent 10-2 <z2>
m
o
m
n
or i <z>
err
ac. 10-3
Fr
10-4
-5 -4 -3 -2 -1 0 1 2 3 4 5
a
10-1
10-2 <z3>
ents 10-3 <z2>
m
o
m
n
or i <z>
err
c. 10-4
a
Fr
10-5
10-6
-5 -4 -3 -2 -1 0 1 2 3 4 5
a
Fig. 1. Absolutefractional errorinthe first threemomentsasa functionof speedratio,a. Upperfigure isfor the boxenvelopewith
z (a)=min(a(cid:2)a, (cid:2)3), z (a)=min(a, 3) for a =1 (solid line) and a =2 (dashed line). Lower figure is the reservoir envelope with
< 1 > 1 1
z (a)=min(a(cid:2)1,(cid:2)3).
<
700 A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708
5.2.1. Envelope 1 (a<0)
For a<0, consider the envelope
2
^p ðzÞ¼ ð(cid:2)zÞexpð(cid:2)z2Þ
1 mðaÞ
and note that (cf. (A.10))
Z z 1
^p ðxÞdx¼ expð(cid:2)z2Þ. ð5:7Þ
1 mðaÞ
(cid:2)1
The acceptance rate (5.4) is
mðaÞ
expð(cid:2)a2Þ
and tends to 1 as a!0 (see Fig. 2). One obtains from (5.7)
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
z(cid:5) ¼(cid:2) a2(cid:2)logu ð5:8Þ
and the acceptance probability (5.3) is
a(cid:2)z(cid:5) (cid:3)0 if z(cid:5) !a;
!
(cid:2)z(cid:5) 1 if z(cid:5) !(cid:2)1.
5.2.2. Envelope 2 (a<0)
For a<0, consider the envelope (cf. (4.2), (4.3))
(cid:3)^p ðzÞ if z<bðaÞ;
1
^p ðzÞ¼
2
p ðzðaÞÞ if z2½bðaÞ;aÞ;
a
1
e 0.8
at
R
e 0.6
c
n
a
ept 0.4
c
c
A
0.2
0
-0.5 -0.4 -0.3 -0.2 -0.1 0 0.1 0.2 0.3 0.4 0.5
a
1
e0.8
at
R
e 0.6
c
n
a
ept0.4
c
c
A
0.2
0
-5 -4 -3 -2 -1 0 1 2 3 4 5
a
Fig.2. Acceptanceratesfor:Envelope1 (solidline, a<0);Envelope2(dashedline, a<0); Envelope3(solid line,a>0); Envelope4
(dashedline,a>0).
A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708 701
where the function b satisfies
bðaÞ6a.
The previous case is obtained for b(a)=a. One obtains
Z a 1
^p ðzÞdz¼ ½expð(cid:2)bðaÞ2Þþ2expð(cid:2)zðaÞ2Þða(cid:2)zðaÞÞ½a(cid:2)bðaÞ(cid:3)(cid:3)
2 mðaÞ
(cid:2)1
and the corresponding acceptance rate (5.4). Our particular choice is
bðaÞ¼a(cid:2)ð1(cid:2)aÞða(cid:2)zðaÞÞ;
which gives a favorable acceptance rate over a wide range of a (see Fig. 2).
For this piece-wise envelope, with probability
expð(cid:2)bðaÞ2Þ
;
expð(cid:2)bðaÞ2Þþ2expð(cid:2)zðaÞ2Þða(cid:2)zðaÞÞ½a(cid:2)bðaÞ(cid:3)
z* is generated according to the probability density
2expðbðaÞ2Þð(cid:2)zÞexpð(cid:2)z2Þ; z<bðaÞ;
so that (cf. (5.8))
qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
z(cid:5) ¼(cid:2) bðaÞ2(cid:2)logu;
and accepted with probability
a(cid:2)z(cid:5)
.
(cid:2)z(cid:5)
With probability
2expð(cid:2)zðaÞ2Þða(cid:2)zðaÞÞ½a(cid:2)bðaÞ(cid:3)
expð(cid:2)bðaÞ2Þþ2expð(cid:2)zðaÞ2Þða(cid:2)zðaÞÞ½a(cid:2)bðaÞ(cid:3)
z* is generated uniformly on [b(a),a] and accepted with probability
a(cid:2)z(cid:5)
expðzðaÞ2(cid:2)ðz(cid:5)Þ2Þ.
a(cid:2)zðaÞ
5.2.3. Envelope 3 (a>0)
For a>0, consider the envelope
(
p ðzÞ if z60;
^p ðzÞ¼ a
3 2mðaÞ(cid:2)1ða(cid:2)zÞ if z2ð0;aÞ;
and note that (cf. (A.9) and (A.10))
Z a 1 p
^p ðzÞdz¼ ða ffipffiffiþ1þa2Þ.
3 mðaÞ
(cid:2)1
The acceptance rate (5.4) is
mðaÞ
p ;
a ffipffiffiþ1þa2
which is shown in Fig. 2 to be close to unity for a<1.
With probability
p
ffiffiffi
a p
p
ffiffiffi
a pþ1þa2
702 A.L.Garcia,W.Wagner/JournalofComputationalPhysics217(2006)693–708
one generates z* according to the probability density
2
p expð(cid:2)z2Þ; z60;
ffiffiffi
p
so that z* is a one-sided Gaussian with parameters 0 and 1/2. With probability
1
p
a ffipffiffiþ1þa2
one generates z* according to the probability density
2ð(cid:2)zÞexpð(cid:2)z2Þ; z60;
so that z* is defined in (3.6). With probability
a2
p
ffiffiffi
a pþ1þa2
one generates z* according to the probability density
2
ða(cid:2)zÞ; z2ð0;aÞ;
a2
so that
p
z(cid:5) ¼að1(cid:2) ffiuffiffiÞ;
and accepts it with probability
expð(cid:2)ðz(cid:5)Þ2Þ.
5.2.4. Envelope 4 (a>0)
For a>0, consider the envelope
(
p ðzÞ if z60;
^p ðzÞ¼ a
4 2mðaÞ(cid:2)1aexpð(cid:2)z2Þ if z>0;
and note that (cf. (A.10))
Z a 1 p
ffiffiffi
^p ðzÞdz¼ ð2a pþ1Þ.
4 mðaÞ
(cid:2)1
The acceptance rate (5.4) is
mðaÞ
p ;
ffiffiffi
2a pþ1
which is shown in Fig. 2 to be favorable over a wide range of values, going to unity as a!1.
With probability
p
ffiffiffi
2a p
p
ffiffiffi
2a pþ1
one generates z* according to the probability density
1
p expð(cid:2)z2Þ; z2R;
ffiffiffi
p
so that z* is Gaussian with zero mean and variance 1/2. With probability
1
p
ffiffiffi
2a pþ1