[neutron-mc] Maxwell random deviate generator
George Apostolopoulos
gapost at ipta.demokritos.gr
Fri Oct 5 12:44:49 CEST 2007
Dear all,
I found today a much more elegant solution to this problem. Actually the
following three lines of code:
{
double u1,u2;
do u1 = rand01(); while (u1==0.0);
do u2 = rand01(); while (u2==0.0);
return sqrt(-log(u1*u2));
}
generate random velocities with the required Maxwell flux distribution
2*exp(-v^2)*v^3. (v is scaled to the thermal velocity v_T = 2kT/m)
What is most interesting for McStas is the following:
Generating velocities with the above routine is slower (about 0.6 times)
than the standard
McStas method of first generating velocities (or equivalently
wavelengths) uniformly and then calculating the particle weight. But:
the same accuracy
of Monte-Carlo results is reached with fewer particle histories when
sampling the true distribution. I get roughly a performance improvement
of 1.5 to 2 on my computer, depending on compiler.
As an example, the small attached program "maxwell2.c" calculates the
average <v^2> with both methods.
The output of "maxwell2.c" on my computer with gcc:
************************************************************************
********
Testing the Maxwell random deviate generator.
The value of <v^2> = 2 will be calculated by 2 methods:
1. Generating random v distributed uniformly in 0.1 < v <5 and
multiplying by the distribution weight w = 2*exp(-v^2)*v^3
2. Generating random v distributed according to 2*exp(-v^2)*v^3
1st method
N <v^2> delta<v^2>/<v^2> time(ms)
100 1.888390 1.463e-001 0.000
1000 2.157340 4.092e-002 0.000
10000 2.010312 1.352e-002 0.000
100000 2.004246 4.283e-003 15.000
1000000 1.998229 1.355e-003 172.000
10000000 1.999877 4.282e-004 1690.000
2nd method
N <v^2> delta<v^2>/<v^2> time(ms)
100 2.109436 6.675e-002 0.000
1000 2.007962 2.170e-002 0.000
10000 2.007131 7.057e-003 0.000
100000 2.000076 2.232e-003 16.000
1000000 1.998449 7.057e-004 235.000
10000000 1.999998 2.234e-004 2582.000
************************************************************************
*****
The pdf file is a short description of the method.
Best regards,
George
-----Original Message-----
From: neutron-mc-bounces at risoe.dk [mailto:neutron-mc-bounces at risoe.dk]
On Behalf Of George Apostolopoulos
Sent: Thursday, October 04, 2007 4:56 PM
To: neutron-mc at risoe.dk
Subject: [neutron-mc] Maxwell random deviate generator
Sorry, forgot to paste the c-code.
******************************************************
Dr. Georgios Apostolopoulos
Institute of Nuclear Technology & Radiation Protection
National Center for Scientific Research "Demokritos"
15310 Agia Paraskevi Attikis, Greece
tel +30-210-650-3731
fax +30-210-653-3431
email: gapost at ipta.demokritos.gr
******************************************************
-------------- next part --------------
A non-text attachment was scrubbed...
Name: maxwell2.pdf
Type: application/octet-stream
Size: 58804 bytes
Desc: maxwell2.pdf
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20071005/de44ced7/attachment.obj>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: maxwell2.c
Type: application/octet-stream
Size: 2638 bytes
Desc: maxwell2.c
URL: <http://mailman2.mcstas.org/pipermail/mcstas-users/attachments/20071005/de44ced7/attachment-0001.obj>
More information about the mcstas-users
mailing list