- Email: [email protected]
The use of symmetry and other features of neutron paths for speeding up computations by the Monte Carlo method
THE USE OF SYMMETRY AND OTHER FEATURES OF NEUTRON PATHS FOR SPEEDING UP COMPUTATIONS BY THE MONTECARLO METHOD* I.G.DYAD’KIN and V. N. STARIKOV Ufa (Received
11 February
1967; reuised
10 April 1967)
SEVERAL means for improving the convergence for the Monte Carlo method are recommended in [l-61. The most important of them involves rational organization of the ,computer program. This includes: (1) achieving a solution of the problem as a whole, i.e. seeking some dependence or formula, and not computing separate points of the dependence in sequence and independently; (2) considering all the points of the dependence simultaneously correlated way, i.e. so that they have matched fluctuations;
and in a
(3) sharply reducing the amount of computation by exchanging computational information at adjacent points; ideally, the random path can be drawn at a single point of the dependence. Recommendations (2) and (3) are fairly easily realized in practice when solving Laplace’s equation and problems on gamma quanta propagation, but are extremely difficult in the case of slow neutrons. The reasons for this lie, firstly, in the variety of types of neutron elastic scattering by different nuclei, which makes it difficult to obtain a correlation of the neutron paths in media of different compositions; secondly, in the large number of collisions when a neutron is slowed from millions to a few parts of an electron-volt, which makes it impossible in practice to use a single path of given “value” for different media; and thirdly, in the low probability of neutron migration at great distances from point sources.
* Zh. uychfsl.
Mat. mat. Fiz. 8, 5, 1001-1012,
87
1968.
I. G. Dyad’kin
and V. N. Starikou
FIG.. Geometry
of the NGC problem:
Co = counters, 0 = polonium-beryllium F = lead/iron paths,
filter;
source;
A, B, C = conditional
drawn for simplicjty
in the zx plane
d,, 4, d, = pores of different
diamters.
In the present paper we propose a special algorithm for the draw of neutron directions, ensuring frequent coincidence of the directions of neutron motion after scattering in different media, even in the case of collisions with different nuclei and differences between the directions prior to scattering. To achieve high accuracy in the probability of neutron migration at great distances, we propose making use of the high degree of symmetry in individual parts of the path, even when there is no symmetry in the geometry of the problem. All these methods, discussed below, have enabled us to evaluate by computer for the first time the dependence of the readings of a neutron gamma counter (NGC) on the hydrogen content of strata and the pore diameter for a
Computations
89
by the Monte Carlo method
standard NGC (as used in all the industrial-geophysical
offices of the Union).
Since this problem will often be used to illustrate the proposed methods, we shall state it more fully. The NGC is mounted eccentrically in a cylindrical pore (see Fig.) filled with salt water. It consists of an iron tube one centimetre thick (its height can be assumed infinite), in which a lead/iron filter is mounted, for isolating the polonium-beryllium neutron source from the three SI-4G gamma quantum counters, located so that the distance between the source and mid-point of the counters is 60 cm. We want to find the dependence of the gamma quantum counter readings on the nuclear composition of the mineral stratum and the liquids filling the pore, the pore diameter and certain other parameters. The physical processes occurring in the NGC are: slowing of fast neutrons to energies of 0.1 eV, thermalization and diffusion of hot neutrons, radiative capture of neutrons with radiation of the corresponding spectrum of gamma quanta, and propagation of the latter in conditions of Compton, photo- and other effects with the possibility of their being registered in the indicator (the energy efficiency of this is known). Not all these processes are of equal importance, and some can be neglected, depending on the physical features of the specific problem of nuclear geophysics. 1. Correlation
between neutron different media
directions
in
Since, in accordance with what has been said, the problem will have to be solved in media with different nuclear compositions and different geometrical dimensions, we must first consider an algorithm whereby similar neutron paths will be obtained in different media. Let p and p1 be the neutron momentum before and after collision. If these quantities refer to different media, the medium number will be denoted as an argument, e.g. p(l) or p,(2). We take the neutron mass as unity, and the nuclear mass as M(i), where i is again the medium number. Obviously, b = pz / 2; R = p/p. If the usual draw algorithm [7l is used (centre of mass and two rotations), the correlation between R,(l) and Q,(2) is insufficient. Even when a set of random numbers is used, they diverge in the case of collisions with different nuclei, and what is more important, once the direction mismatch process has started, it becomes aggravated. With this algorithm, the probability R,(l) = Q,(2) when
1. G. Dyad’kin
90
and V. N. Starikm
M(l)# M(2)(andallthemorewhen Q(l) f Q(2) 1 becomes zero. We therefore used an algorithm of direct drawing Q, in the initial coordinate system, which ensures a high value for the above-mentioned coincidence probability. It is based on the following: the probability of the momentum lying in the interval tp,, p1 + dp,) after elastic scattering, symmetric in the centre of mase system, is
It can be written as the product of condi~on~
p~bability
densities:
where
We haveusedhere the rules for inversion with Dirac s-functions ‘[9f. In view of the fact that dW(Qi) --= d(QQi) we
(QQi + x)2(2xM - QQi) > 9 43r%~
t
have supW(&)=
w4nM2 + q2 .
The following ~gorit~ for the drawing Q, is now obvious (here and below the $3 with different subscripts denote random numbers uniformly distributed on (0, 11):
Computations
91
by the Monte Carlo method
(a) the Cartesian coordinates of the vector Q: in the initial coordinate system are drawn equiprobably on the unit sphere, for Qix’ -
i621,’ =
)‘[4fA (1-
pi)]
exp(2ni~2);
f&r =
281-
1,
v.9
or by Neumann’s method 171; (b) we check the inequality -(QQI’ + X&fY> XM’
p3 (M + l)2
(1.2)
M
or its equivalent M2(M -
1)2 >
l33[ (PC&‘)2 + M2 -
I] [(M + I)283 -
4MQQ,‘];
if it is satisfied, Q’, = Q,, and we pass to (19; otherwise, we return to (a) and repeat the procedure with different PI, &, /3,; (c)
we evaluate (1.3)
In particular, for hydrogen, (1.2) and (1.3) transform to QQ’, > p3, E, = G2Q1JzE, while for nuclei that are not particularly light, we have, when M >> 1,
M + 2QQ1’3 (M + 2) hi,
El
N
("$-;F1)'E.
We can easily compute the probability of the first of the Q’, satisfying (1.21, since it is equal to M’(M + l)-‘,
i.e. is equal to t/, even for hydrogen, while for the other
elements it is close to unity. We now find the probability x of Q,(l) and Q,(2) being identical when this algorithm is used in two differelf
media.
Starting from (1.2), we obtain for this probability: x =
SdQimin {~(Q(l),M(l),Qi);1L(Q(2),M(2),Qi))
s
daimax{~(Q(l),M(l),Bi);
$(Q(2),M(2),Qi)}’
(1.4)
1. G. Dyad’kin
92
and V. N. Starikov
where
M(i)
~(n(i),M(i),Pi)=,[M(i)+
+
CQ(i)Qi+ xaf((i)P
, ‘-
x&f(i)
The integrals in (1.4) are difficult to evaluate in the general case, but they can be completely investigated by taking the following cases. When Q(1) = Q(2), the dilemma is resolved by using one of the ti for the entire region of integration. In this case
x = mm
'
Q(l)= P(2),
Q,(l)=
(1.5)
Q2(2).
When Q(l) + Q(2): (a) if M(1) >> 1 and M(2) >> 1, we can estimate y by taking the minimum (maximum) possible value of the integrand outside the integral sign in the numerator (denominator):
(b) ifM(l)=l,M(2)>>1,
then x=0.25+0(1/M(2));
(c) if M(1) = M(2) = 1, the integrals can be evaluated by choosing the coordinate axes so that P (1) = (cos 8 = 0; 9 = 0) ; Q(2) = (cos8=O;cp=arccosC;1(~)Q(2)
a,
=
in which case .,c +min
{Q(~)QI+JQ(I)RIJ;
=
32(2)Qi+[Q(2)QiI)dQi
1-sin-
f
;
2
similar evaluation of the denominator of (1.4) gives cos (6/ 2), so that
Computations
x = tg
35 -
{
by the Monte Carlo method
ar~cos~(~)~~2) 4
),
M(l)
93
= ~(2)
= 1.
(47)
On the basis of (1.5)~(1.7) we can conclude that, independently of whether the directions of neutron motion are the same in different media initially, after scattering they will be the same: (1) with probability close to 1 if both nuclei are somewhat heavier than hydrogen; (2) with probability close to ‘/, if just one of the nuclei is a proton; (3) with probability varying between 0 and 1 if both nuclei are protons; assuming that the angle between R(l) and Q(2) is random, we have on the average
Thus the proposed algorithm not only gives a high probability of exact coincidence of the neutron directions in different media, but also operates in such a way that, if the directions are divergent, they tend to be brought into coincidence. Notice also that, as distinct from the algorithm of 171, the result is obtained directly in a fixed coordinate system, so that two rotations are ~necess~. This ~gorit~ has ensured good correlation between the paths in NCC problems. 2. Use of the symmetry of the neutron
paths
We now consider the difficulties due to the low probability of neutron migration at great distances and in the vicinity of indicators. The method proposed in the literator (e.g. 121) of drawing the most highly valued paths yields virtually nothing in problems with a complex geometry, where it is impossible to estimate with sufficient accuracy the required solution of the conjugate transport equation. Even a small divergence in its solution leads to sharp fluctuations of the corresponding statistical weights due to the large number of neutron collisions. But it is precisely the large number of collisions that suggests using the symmetry of different parts of the paths to obtaining a set of stopping points widely distribute t~oughout space (including in the vicinity of the indicator).
94
1. G. Dyad’kin
and V. N. Stnrikou
Even if the geometry of the problem is very asymmetric, there is still symmetry in the parts of the paths. We illustrate this idea by taking the NGC problem, where the geometry is typical for nuclear geophysics. Our figure is very lacking in symmetry. It maps into itself only when reflected in a plane passing through the pore axis and the inst~ment axis. Now take the trajectory A, the collision points on which are numbered in the direction of increasing time. (Of course, along an actual path there are on the average at least 25 collisions, even in water). Obviously, the last point A,,, where the neutron has been braked and captured, can be rotated about any of the axes Ai A j+l, where i = 9, 8, 7, 6, 4 (but not i = 5, 2, 11, since under these rotations the next piece in time of the path (A i+t A i+l.. .A J does not cross the boundaries of the.media. This fact alone enables us to obtain from A,, a whole series of points of radiative capture symmetric to it, including points which may be in the vicinity of the indicator. If we agree to rotate the terminal points of all the N drawn paths about one of the permissible axes so that the final point is rotated n times through the angle 2&n, we obtain nN path ends with identical weights l/n. In the case of a path of the type B or C, however, no axes of rotation may be available for the final point, or they may be unsuitable (close to the end or situated in the pore). In the case of these paths, we can rotate say the points B, or C, with a subsequent draw from all the symmetric points of the final pieces of paths up to radiative capture. If the same ~gorit~s are used, we can here in fact draw only some of the pieces, while the rest are obtained simply from symmetry. It is undesirable to rotate the end of the path B,, (or the point Clzl about the axes of symmetry lying in the pore (e.g. about B&,, or B,,B,,) in the NGC problem, because the pieces of Notice that the majority of the path in the pore contain only a few collisions. paths are of type A in this problem. A rotation can also be performed in every path about two axes: fist about one, at n points, then rotations of all these points about the other, again at n points. A total of n2 symmetric ends is thus obtained. We only have to bear in mind that the first rotation must be performed about the axis which is latest in time (for instance, first about AdSt,, then about A,A,, not vice versa). Then, by using the symmetry of individual pieces of the paths, we can obtain from one end of a path a set of other ends with very little expenditure of computer time; and what is specially import, the statistics weights are indent&l for all the ends of all the paths. The possibilities are even greater if the media consist of nuclei with strongly divergent masses (e.g. the case of porous limes with masses 1, 12, 16, and 401.
Computations
by the Monte Carlo method
95
From physical considerations, we can then neglect the asymmetry of the scatter in the laboratory system of coordinates and the energy losses at certain points of the path, where collision has occurred with heavier nuclei. For instance, let A i be one such point. Then obviously, A 1).can be transferred to another point of the sphere with centre A i and radius equal to the distance between A i and A,,, provided this sphere does not cut the boundary between the media. In this case, however, certain s~rnet~ ~~sfo~ations are always possible. For instance, A,, can be reflected in the plane passing through Ai and the pore axis (call this reflection $1, and also it can be reflected in the plane though A, perpendicular to the pore axis (call this a z reflection). These transformations, which are always possible for A-paths, are extremely convenient, since they enable us, when rotations are im~ssible or undesirable, to ‘reduce” such paths to the same statistical weights, i.e. they assist us to liquidate fluctuations of the statistical weights (concerning this, see Section 4). Notice that symmetry transformations in no way destroy the Monte Carlo stochastic scheme. Given a poor (but not incorrect) choice of the scheme of symmetry transformations, the only detriment is to increase the dispersion of the final results (which can be assessed by independent computations). As a rule, the symmetry scheme reducing the fluctuations of the neutron field close to the indicator can easily be selected on the basis of physical and geometric considerations. Its realization in the NGC problem will be described in Section 4. 3. The impossibility of using a single neutron path for different media It would be very attractive if a single neutron path could be used for media of different compositions, with suitable statistical weights ascribed to it (as in, for exsmple, electrostatic problems f61).Let us examine if this is possible. For simplicity, consider the stopping of neutrons without capture, due to elastic scattering, symmetric in the centre of mass system. The probability density of a neutron colliding at the point r, changing its momentum from p = ~$2to p1 = p&,, then travelling freely to the point r + sQ,, s > 0, is
(3.1)
This relation is normalized so that the integral with respect to dp,ds is equal to 1. Summation is over all K types of nucleus of mass Mi; Z(r) is the reciprocal of the complete macroscopic cross-section of neutron interaction with nuclei at the
96
I. D. Dyad’kin‘and
V. N. Starikov
point r, ci is the probability of collision with the i-th nucleus, 2
Ci =
1.
i-i
We shall assume that the nucleus com~sition is qu~itatively the same in the two differetit media (otherwise the same path cannot be used, since the presence of the &function would lead to a statistical weight 0 or M appearing in the case of.collisions with different nuclei). To estimate the fluctuations of the statistical weight of the path drawn in the first medium and assigned without variations of the second medium, we assume that the number of collisions is n, that the medium has no boundaries and that I and c are independent of the energy. Then, introducing the notation ci and 2 for the first medium, and ci + AC, I+ Al for the second (ci and Al, are not assumed small, though obviously, ACi-0)
i
, we obtain the following for the mean statistical
weight [, its mean
i=i
square Fz and its standard deviation d([’
fv&=
- 5’):
(3.2)
(3.3)
It is clear from (3.31 that it is actually better to draw the path in the medium with the longer path, this being the only possible way of avoiding infinite ~uctu~ions of the weights in the case Al > 1. #Yenow consider how advantageous it is to use a single path. Suppose that the machine time for computing the weights is negligible, and let the draw be made for Q media simultaneously. Then, is a time T is required for drawing N histories in the case of an indepedent draw,
Computations
Nq single draw was the other it is only
by the Monte Carlo method
97
paths will he drawn during this time. Hence, in the medium where the actually made, the error of the result will be less by a factor dq, and in media, by a factor dq/oE, than in the case of an independent draw. Hence advantageous to use a single path when OE<
1’9
(3.5)
From (3.4, uC can only be small if
In this case (3.5) gives
(3.6) Noting that n is of the order of 30-100, we now see that a single trajectory is only suitable for the_unrealistic case of a large number of media with very similar compositions, and then, by comparison with an independent draw, not with a correlated draw. But it can be, and often is, useful to use a single path up to the instance when the weight deviates from 1 by a certain amount, and thereafter branch the path, observing the correlation and symmetry methods described above. This in no way contradicts the Markov property of the process, since at any step of the path we preserve the ability either to draw it separately for all the media, or draw it for one, while introducing a statistical weight for any other. Or we might concervably adopt the device of using a single path for a large number of media, then “extinguishing” it by means of symmetry transformations in those media where the statistical weight becomes very large. Notice also tnat a single path can also almost be drawn in media of different qualitative compositions, provided they do not contain hydrogen. In this case we have to draw a specific type of nucleus in each medium separately, then, after obtaining the new direction of neutron movement in the first medium, multiply the statistical weight of the second medium by (see Section 1) M(2) / an,+-)l(M~(2)--_IPXQi12) M(1)\
szn,+~(@(1)-~QxQi12)
2 )I
~2(W--_(QxQ~12 iW(2)-+XQ,~~
-
(3*7)
I. G. Dyad’kin
98
and V. N. Starikov
Of course we must here assume a specific energy for each medium in accordance with (1.21, so that it can happen that the path has terminated in one medium, while in another its draw has to be completed up to the thermal energy. Since (3.7) is very nearly 1 and replaces (3.2), this method, using a specific nucleus and specific energy, is also suitable in the case of media of qualitatively similar composition. All this indicates that there is still no sound basis for using a single neutron path (though we can hope for some success from this approach). At the same time, the methods for Sections 1 and 2 work well. 4.
Features of NGC program organization; comparison of an experiment and check computation
On instructions from the Bashkirskii Geophysical Trust, we computed the dependence of NGC readings on the porosity of limestone (pores filled with water) and the pore diameter. The volumetric porosity was m = 1, 2, 3, 5, 10, 15, 20, 40 and 100%. The pore diameters were d = 5%, 7% and 11% in. The pore was filled with salt water. During the computations, account was taken of the following important physical phenomena, not accessible to analytic theory: the spectrum of the polonium-beryllium source, the neutron energy losses due to inelastic and elastic collisions, the spectrum of the radiative capture of gamma quanta, the Compton effect and photo-effect for propagation of gamma quanta, the spectral sensitivity of the SI-4G counters, the energy dependence of all the effective neutron and quantum cross-sections, and so on. No account was taken of the following: (1) neutron thermalization and diffusion; (A neutron was assumed to be captured where it was stopped. This is possible because, when the porosity varies, the stopping parameters vary by several hundred percent, while the diffusion lengths, which are one order less than the stopping lengths, have virtually no effect on the migration, so that the NGC dependence on the hydrogen content is not distorted). (2) gamma quanta of nonelastic the radiative quanta;
scatter,
since they are small compared with
(3) the effect of vapour formation during propagation of the gamma quanta, the cross-sections of which are small in limestone at gamma quantum energies up to 10 MeV. So much for the physical assumptions.
We now turn to the computer program
Co~put~~ons
by the Monte Carlo method
99
organization. 1. The computation was performed simultaneously for all m and d (30 variants), so that, for instance, we computed the neutron travel and the new coordinates in the m and d cycle, then the scatter in the m and d cycle, and so on. Many of the qu~tities required for all m and d (such as In ,6 for the travel, Q’ from (1.1) and so on) were computed once. 2. In order to retain a single path as long as possible for all 30 variants, [41, the coordinate axes for the different pore diameters were chosen so that the instrument position was the same for all d, while the pores themselves were in contact (Fig. 1). Hence, till the initial piece of path went beyond 5%“, it was the same for all 30 variants; if it were left after 5%“, but before 7?4”, it was the if it left after 7%“, but before same for 20 variants (while 10 were br~ch~); llW”, it was the same for 10 variants. 3. Correlation of the neutron path directions of Section 2.
was realized by the algorithm
4. Every neutron path of type A was subjected to symmetry transformations relative to the collision points with calcium or oxygen in accordance with what was said in Section 3, on the basis of the following considerations. Most of the paths leaving a pore’fspecially the longest in the case of low porosities) have a tendency to go a long way into the stratum and to migrate little in the z direction. The problem therefore amounts to extracting information from them about the stopping field in the zone round the pore and for z close to the indicator. This was done as follows. Call the points of collision with heavy nuclei in the stratum Dj, and the points of radiative capture A,. During the drawing of the path, we stored the D, for which the z coordinate was a maximum or minimum. These were the points with respect to which we reflected the point A, once. As a result, 4 points were obtained, reducing the fluctuations of the neutron field along the z axis. If the neutron path went beyond the pore wall by more than the stopping length in the stratum (for the energy at this point), the coordinates of the nearest point Di were stored, and for the rest of the path, we analyzed its departure from the line through D i parallel to the pore axis. It was decided on this basis whether a rotation of A, about the given line was possible; the rotation was through 21718(eight points). If no rotation was possible, single reflections were performed with respect to two points, co~es~nding to the rn~irn~ and minimum angular departure from the plane (4 points). The line of rotation was chosen so that the rotation (like the 4 reflection) commuted with the z reflection, thus simplifying the logic. In this way, either
100
I. G. Dyad’kin
and V. N. Starikov
TABLE
m 5%
Theoret. rel.error I + 15%
100
0.15
40 20 15 10
0.42 0.45 0.46 0.45
Pore & lab. data/ r;@rror 0
0.44 0.51 0.56 0.65
m %
Theoret. rel.error I f 15%
Pore & lab. data,+r;k;rror 0
5
0.69
0.78
3 2 1 0.5
0.83 0.86 1 0.80
0.85 0.92 1
32 or 16 symmetric points with weights 2-’ or 2-’ were obtained as a result of the -symmetry transformations. As regards the type B and C paths, since they are very rarely encountered in the present problem in the region of the indicator, it was decided not to complicate the program by performing symmetry transformations with them. They were left with the weight 1. 5. A gamma quantum path has to be drawn from each point of radiative capture. This was done as follows. For each m variant, gamma quantum paths were drawn in the homogeneous medium corresponding to the stratum, and we analyzed its departure from the line parallel to the pore axis and passing through the point of quantum birth. This path was then “fitted” against all the symmetric points of radiative capture; if it did not cross the pore during the fitting an answer was obtained directly about its being missed; if it did cut the pore, the answer was obtained after redrawing the piece cutting the pore. The fitting was carried out, not only for the initial path, but also for the three paths rotated through 2~14 from it about the above-mentioned line. In this way the weight of a quantum of radiative capture was ?4 (we do not include here the weight of the neutron path and the weight corresponding to the number of quanta per act of neutron capture). If the neutron path terminated at the pore, 4 rotated quantum paths were drawn from it in the true geometry in the usual way. The computations were performed on the Razdan-111 computer (speed approx. 10’ op/sec) of the Academy of Sciences of the Armenian SSR and the Erevan University Computer Centre, during 100 hours of machine time (as much time was spent on harmonizing the program, dealing with the working block of 10,000 cells, and checking the control variants). The program operation was monitored, in particular, by the method of printing out the neutron field. Experimental data were only obtained for a single pore with d = 5%“, bored out with complete removal of the core. The table gives a comparison of the theoretical curves with these data, obtained by M. S. Pernikov. It throws light on the suitability of the proposed
Computations
101
by the Monte Carlo method
mathematical model of the NGC. All this work was performed on the initiative of I. L. Dvorkin, who carried out all the geophysical analysis of the statement of the problem and of the results obtained (which are not quoted here). The authors thank the ByeleRussian assistance with the computations.
State University students for
REFERENCES
1.
KAHN,
H.
Use of different Monte Carlo sampling techniques. N.Y. 146-191, 1956.
Symposium
on Monte
Carlo Methods, 2.
GOERTZEL,
A. and KALOS,
Nucl. Energy,
Ser I, Phys.
M. Monte Carlo methods in transport problems. Math. 2, 315-395, 1958.
Prog.
3.
FROLOV, A. S. and CHENTSOV, N. N. On computation by the Monte Carlo method of definite integrals depending on a parameter. Zh. uychisl. Mat. mat. Fiz. 2, 4, 714-717, 1962.
4.
DYAD’KIN, I. G. LISENENKOV, A. T. and PONYATOV, G. I. On the convergence rate of the Monte Carlo method when solving problems of radioactive logging. Zh. vychisl. Mat. mat. Fiz. 5, 4, 763-768, 1965.
5.
DYAD’KIN and PONYATOV, G. I. On the theory and method of mathematical modelling of gamma quanta propagation. lzu. Akad. Nauk SSSR, Fiz., Zemli, 1,47-58, 1966.
6.
DYAD’KIN, I. G. and STARIKOV, V. N. The possibility of economizing on machine time when solving Laplace’s equation by the Monte Carlo method. Zh. uychisl. Mat. mat. Fiz. 5, 936-938, 1965.
7.
BUSLENKO, N. P. and GOLENKO, D. I. et al. Monte Carlo techniques statisticheskikh ispytaniil. Moscow, Fizmatgiz, 1962.
8.
AKHIEZER, voprosy
9.
IVANENKO, D. and SOKOLOV, A. Class ical field polya). Moscow, Gostekhizdat, 1951.
(Meted
A. and POMERANCHUK, I. Some aspects of nuclear theory (Nekotorye teorii yadral. Moscow-Leningrad, Gostekhizdat, 1959.
theory (Klassicheskaya
teoriya
Translated by D. E. Brown
11 February
1967; reuised
10 April 1967)
SEVERAL means for improving the convergence for the Monte Carlo method are recommended in [l-61. The most important of them involves rational organization of the ,computer program. This includes: (1) achieving a solution of the problem as a whole, i.e. seeking some dependence or formula, and not computing separate points of the dependence in sequence and independently; (2) considering all the points of the dependence simultaneously correlated way, i.e. so that they have matched fluctuations;
and in a
(3) sharply reducing the amount of computation by exchanging computational information at adjacent points; ideally, the random path can be drawn at a single point of the dependence. Recommendations (2) and (3) are fairly easily realized in practice when solving Laplace’s equation and problems on gamma quanta propagation, but are extremely difficult in the case of slow neutrons. The reasons for this lie, firstly, in the variety of types of neutron elastic scattering by different nuclei, which makes it difficult to obtain a correlation of the neutron paths in media of different compositions; secondly, in the large number of collisions when a neutron is slowed from millions to a few parts of an electron-volt, which makes it impossible in practice to use a single path of given “value” for different media; and thirdly, in the low probability of neutron migration at great distances from point sources.
* Zh. uychfsl.
Mat. mat. Fiz. 8, 5, 1001-1012,
87
1968.
I. G. Dyad’kin
and V. N. Starikou
FIG.. Geometry
of the NGC problem:
Co = counters, 0 = polonium-beryllium F = lead/iron paths,
filter;
source;
A, B, C = conditional
drawn for simplicjty
in the zx plane
d,, 4, d, = pores of different
diamters.
In the present paper we propose a special algorithm for the draw of neutron directions, ensuring frequent coincidence of the directions of neutron motion after scattering in different media, even in the case of collisions with different nuclei and differences between the directions prior to scattering. To achieve high accuracy in the probability of neutron migration at great distances, we propose making use of the high degree of symmetry in individual parts of the path, even when there is no symmetry in the geometry of the problem. All these methods, discussed below, have enabled us to evaluate by computer for the first time the dependence of the readings of a neutron gamma counter (NGC) on the hydrogen content of strata and the pore diameter for a
Computations
89
by the Monte Carlo method
standard NGC (as used in all the industrial-geophysical
offices of the Union).
Since this problem will often be used to illustrate the proposed methods, we shall state it more fully. The NGC is mounted eccentrically in a cylindrical pore (see Fig.) filled with salt water. It consists of an iron tube one centimetre thick (its height can be assumed infinite), in which a lead/iron filter is mounted, for isolating the polonium-beryllium neutron source from the three SI-4G gamma quantum counters, located so that the distance between the source and mid-point of the counters is 60 cm. We want to find the dependence of the gamma quantum counter readings on the nuclear composition of the mineral stratum and the liquids filling the pore, the pore diameter and certain other parameters. The physical processes occurring in the NGC are: slowing of fast neutrons to energies of 0.1 eV, thermalization and diffusion of hot neutrons, radiative capture of neutrons with radiation of the corresponding spectrum of gamma quanta, and propagation of the latter in conditions of Compton, photo- and other effects with the possibility of their being registered in the indicator (the energy efficiency of this is known). Not all these processes are of equal importance, and some can be neglected, depending on the physical features of the specific problem of nuclear geophysics. 1. Correlation
between neutron different media
directions
in
Since, in accordance with what has been said, the problem will have to be solved in media with different nuclear compositions and different geometrical dimensions, we must first consider an algorithm whereby similar neutron paths will be obtained in different media. Let p and p1 be the neutron momentum before and after collision. If these quantities refer to different media, the medium number will be denoted as an argument, e.g. p(l) or p,(2). We take the neutron mass as unity, and the nuclear mass as M(i), where i is again the medium number. Obviously, b = pz / 2; R = p/p. If the usual draw algorithm [7l is used (centre of mass and two rotations), the correlation between R,(l) and Q,(2) is insufficient. Even when a set of random numbers is used, they diverge in the case of collisions with different nuclei, and what is more important, once the direction mismatch process has started, it becomes aggravated. With this algorithm, the probability R,(l) = Q,(2) when
1. G. Dyad’kin
90
and V. N. Starikm
M(l)# M(2)(andallthemorewhen Q(l) f Q(2) 1 becomes zero. We therefore used an algorithm of direct drawing Q, in the initial coordinate system, which ensures a high value for the above-mentioned coincidence probability. It is based on the following: the probability of the momentum lying in the interval tp,, p1 + dp,) after elastic scattering, symmetric in the centre of mase system, is
It can be written as the product of condi~on~
p~bability
densities:
where
We haveusedhere the rules for inversion with Dirac s-functions ‘[9f. In view of the fact that dW(Qi) --= d(QQi) we
(QQi + x)2(2xM - QQi) > 9 43r%~
t
have supW(&)=
w4nM2 + q2 .
The following ~gorit~ for the drawing Q, is now obvious (here and below the $3 with different subscripts denote random numbers uniformly distributed on (0, 11):
Computations
91
by the Monte Carlo method
(a) the Cartesian coordinates of the vector Q: in the initial coordinate system are drawn equiprobably on the unit sphere, for Qix’ -
i621,’ =
)‘[4fA (1-
pi)]
exp(2ni~2);
f&r =
281-
1,
v.9
or by Neumann’s method 171; (b) we check the inequality -(QQI’ + X&fY> XM’
p3 (M + l)2
(1.2)
M
or its equivalent M2(M -
1)2 >
l33[ (PC&‘)2 + M2 -
I] [(M + I)283 -
4MQQ,‘];
if it is satisfied, Q’, = Q,, and we pass to (19; otherwise, we return to (a) and repeat the procedure with different PI, &, /3,; (c)
we evaluate (1.3)
In particular, for hydrogen, (1.2) and (1.3) transform to QQ’, > p3, E, = G2Q1JzE, while for nuclei that are not particularly light, we have, when M >> 1,
M + 2QQ1’3 (M + 2) hi,
El
N
("$-;F1)'E.
We can easily compute the probability of the first of the Q’, satisfying (1.21, since it is equal to M’(M + l)-‘,
i.e. is equal to t/, even for hydrogen, while for the other
elements it is close to unity. We now find the probability x of Q,(l) and Q,(2) being identical when this algorithm is used in two differelf
media.
Starting from (1.2), we obtain for this probability: x =
SdQimin {~(Q(l),M(l),Qi);1L(Q(2),M(2),Qi))
s
daimax{~(Q(l),M(l),Bi);
$(Q(2),M(2),Qi)}’
(1.4)
1. G. Dyad’kin
92
and V. N. Starikov
where
M(i)
~(n(i),M(i),Pi)=,[M(i)+
+
CQ(i)Qi+ xaf((i)P
, ‘-
x&f(i)
The integrals in (1.4) are difficult to evaluate in the general case, but they can be completely investigated by taking the following cases. When Q(1) = Q(2), the dilemma is resolved by using one of the ti for the entire region of integration. In this case
x = mm
'
Q(l)= P(2),
Q,(l)=
(1.5)
Q2(2).
When Q(l) + Q(2): (a) if M(1) >> 1 and M(2) >> 1, we can estimate y by taking the minimum (maximum) possible value of the integrand outside the integral sign in the numerator (denominator):
(b) ifM(l)=l,M(2)>>1,
then x=0.25+0(1/M(2));
(c) if M(1) = M(2) = 1, the integrals can be evaluated by choosing the coordinate axes so that P (1) = (cos 8 = 0; 9 = 0) ; Q(2) = (cos8=O;cp=arccosC;1(~)Q(2)
a,
=
in which case .,c +min
{Q(~)QI+JQ(I)RIJ;
=
32(2)Qi+[Q(2)QiI)dQi
1-sin-
f
;
2
similar evaluation of the denominator of (1.4) gives cos (6/ 2), so that
Computations
x = tg
35 -
{
by the Monte Carlo method
ar~cos~(~)~~2) 4
),
M(l)
93
= ~(2)
= 1.
(47)
On the basis of (1.5)~(1.7) we can conclude that, independently of whether the directions of neutron motion are the same in different media initially, after scattering they will be the same: (1) with probability close to 1 if both nuclei are somewhat heavier than hydrogen; (2) with probability close to ‘/, if just one of the nuclei is a proton; (3) with probability varying between 0 and 1 if both nuclei are protons; assuming that the angle between R(l) and Q(2) is random, we have on the average
Thus the proposed algorithm not only gives a high probability of exact coincidence of the neutron directions in different media, but also operates in such a way that, if the directions are divergent, they tend to be brought into coincidence. Notice also that, as distinct from the algorithm of 171, the result is obtained directly in a fixed coordinate system, so that two rotations are ~necess~. This ~gorit~ has ensured good correlation between the paths in NCC problems. 2. Use of the symmetry of the neutron
paths
We now consider the difficulties due to the low probability of neutron migration at great distances and in the vicinity of indicators. The method proposed in the literator (e.g. 121) of drawing the most highly valued paths yields virtually nothing in problems with a complex geometry, where it is impossible to estimate with sufficient accuracy the required solution of the conjugate transport equation. Even a small divergence in its solution leads to sharp fluctuations of the corresponding statistical weights due to the large number of neutron collisions. But it is precisely the large number of collisions that suggests using the symmetry of different parts of the paths to obtaining a set of stopping points widely distribute t~oughout space (including in the vicinity of the indicator).
94
1. G. Dyad’kin
and V. N. Stnrikou
Even if the geometry of the problem is very asymmetric, there is still symmetry in the parts of the paths. We illustrate this idea by taking the NGC problem, where the geometry is typical for nuclear geophysics. Our figure is very lacking in symmetry. It maps into itself only when reflected in a plane passing through the pore axis and the inst~ment axis. Now take the trajectory A, the collision points on which are numbered in the direction of increasing time. (Of course, along an actual path there are on the average at least 25 collisions, even in water). Obviously, the last point A,,, where the neutron has been braked and captured, can be rotated about any of the axes Ai A j+l, where i = 9, 8, 7, 6, 4 (but not i = 5, 2, 11, since under these rotations the next piece in time of the path (A i+t A i+l.. .A J does not cross the boundaries of the.media. This fact alone enables us to obtain from A,, a whole series of points of radiative capture symmetric to it, including points which may be in the vicinity of the indicator. If we agree to rotate the terminal points of all the N drawn paths about one of the permissible axes so that the final point is rotated n times through the angle 2&n, we obtain nN path ends with identical weights l/n. In the case of a path of the type B or C, however, no axes of rotation may be available for the final point, or they may be unsuitable (close to the end or situated in the pore). In the case of these paths, we can rotate say the points B, or C, with a subsequent draw from all the symmetric points of the final pieces of paths up to radiative capture. If the same ~gorit~s are used, we can here in fact draw only some of the pieces, while the rest are obtained simply from symmetry. It is undesirable to rotate the end of the path B,, (or the point Clzl about the axes of symmetry lying in the pore (e.g. about B&,, or B,,B,,) in the NGC problem, because the pieces of Notice that the majority of the path in the pore contain only a few collisions. paths are of type A in this problem. A rotation can also be performed in every path about two axes: fist about one, at n points, then rotations of all these points about the other, again at n points. A total of n2 symmetric ends is thus obtained. We only have to bear in mind that the first rotation must be performed about the axis which is latest in time (for instance, first about AdSt,, then about A,A,, not vice versa). Then, by using the symmetry of individual pieces of the paths, we can obtain from one end of a path a set of other ends with very little expenditure of computer time; and what is specially import, the statistics weights are indent&l for all the ends of all the paths. The possibilities are even greater if the media consist of nuclei with strongly divergent masses (e.g. the case of porous limes with masses 1, 12, 16, and 401.
Computations
by the Monte Carlo method
95
From physical considerations, we can then neglect the asymmetry of the scatter in the laboratory system of coordinates and the energy losses at certain points of the path, where collision has occurred with heavier nuclei. For instance, let A i be one such point. Then obviously, A 1).can be transferred to another point of the sphere with centre A i and radius equal to the distance between A i and A,,, provided this sphere does not cut the boundary between the media. In this case, however, certain s~rnet~ ~~sfo~ations are always possible. For instance, A,, can be reflected in the plane passing through Ai and the pore axis (call this reflection $1, and also it can be reflected in the plane though A, perpendicular to the pore axis (call this a z reflection). These transformations, which are always possible for A-paths, are extremely convenient, since they enable us, when rotations are im~ssible or undesirable, to ‘reduce” such paths to the same statistical weights, i.e. they assist us to liquidate fluctuations of the statistical weights (concerning this, see Section 4). Notice that symmetry transformations in no way destroy the Monte Carlo stochastic scheme. Given a poor (but not incorrect) choice of the scheme of symmetry transformations, the only detriment is to increase the dispersion of the final results (which can be assessed by independent computations). As a rule, the symmetry scheme reducing the fluctuations of the neutron field close to the indicator can easily be selected on the basis of physical and geometric considerations. Its realization in the NGC problem will be described in Section 4. 3. The impossibility of using a single neutron path for different media It would be very attractive if a single neutron path could be used for media of different compositions, with suitable statistical weights ascribed to it (as in, for exsmple, electrostatic problems f61).Let us examine if this is possible. For simplicity, consider the stopping of neutrons without capture, due to elastic scattering, symmetric in the centre of mass system. The probability density of a neutron colliding at the point r, changing its momentum from p = ~$2to p1 = p&,, then travelling freely to the point r + sQ,, s > 0, is
(3.1)
This relation is normalized so that the integral with respect to dp,ds is equal to 1. Summation is over all K types of nucleus of mass Mi; Z(r) is the reciprocal of the complete macroscopic cross-section of neutron interaction with nuclei at the
96
I. D. Dyad’kin‘and
V. N. Starikov
point r, ci is the probability of collision with the i-th nucleus, 2
Ci =
1.
i-i
We shall assume that the nucleus com~sition is qu~itatively the same in the two differetit media (otherwise the same path cannot be used, since the presence of the &function would lead to a statistical weight 0 or M appearing in the case of.collisions with different nuclei). To estimate the fluctuations of the statistical weight of the path drawn in the first medium and assigned without variations of the second medium, we assume that the number of collisions is n, that the medium has no boundaries and that I and c are independent of the energy. Then, introducing the notation ci and 2 for the first medium, and ci + AC, I+ Al for the second (ci and Al, are not assumed small, though obviously, ACi-0)
i
, we obtain the following for the mean statistical
weight [, its mean
i=i
square Fz and its standard deviation d([’
fv&=
- 5’):
(3.2)
(3.3)
It is clear from (3.31 that it is actually better to draw the path in the medium with the longer path, this being the only possible way of avoiding infinite ~uctu~ions of the weights in the case Al > 1. #Yenow consider how advantageous it is to use a single path. Suppose that the machine time for computing the weights is negligible, and let the draw be made for Q media simultaneously. Then, is a time T is required for drawing N histories in the case of an indepedent draw,
Computations
Nq single draw was the other it is only
by the Monte Carlo method
97
paths will he drawn during this time. Hence, in the medium where the actually made, the error of the result will be less by a factor dq, and in media, by a factor dq/oE, than in the case of an independent draw. Hence advantageous to use a single path when OE<
1’9
(3.5)
From (3.4, uC can only be small if
In this case (3.5) gives
(3.6) Noting that n is of the order of 30-100, we now see that a single trajectory is only suitable for the_unrealistic case of a large number of media with very similar compositions, and then, by comparison with an independent draw, not with a correlated draw. But it can be, and often is, useful to use a single path up to the instance when the weight deviates from 1 by a certain amount, and thereafter branch the path, observing the correlation and symmetry methods described above. This in no way contradicts the Markov property of the process, since at any step of the path we preserve the ability either to draw it separately for all the media, or draw it for one, while introducing a statistical weight for any other. Or we might concervably adopt the device of using a single path for a large number of media, then “extinguishing” it by means of symmetry transformations in those media where the statistical weight becomes very large. Notice also tnat a single path can also almost be drawn in media of different qualitative compositions, provided they do not contain hydrogen. In this case we have to draw a specific type of nucleus in each medium separately, then, after obtaining the new direction of neutron movement in the first medium, multiply the statistical weight of the second medium by (see Section 1) M(2) / an,+-)l(M~(2)--_IPXQi12) M(1)\
szn,+~(@(1)-~QxQi12)
2 )I
~2(W--_(QxQ~12 iW(2)-+XQ,~~
-
(3*7)
I. G. Dyad’kin
98
and V. N. Starikov
Of course we must here assume a specific energy for each medium in accordance with (1.21, so that it can happen that the path has terminated in one medium, while in another its draw has to be completed up to the thermal energy. Since (3.7) is very nearly 1 and replaces (3.2), this method, using a specific nucleus and specific energy, is also suitable in the case of media of qualitatively similar composition. All this indicates that there is still no sound basis for using a single neutron path (though we can hope for some success from this approach). At the same time, the methods for Sections 1 and 2 work well. 4.
Features of NGC program organization; comparison of an experiment and check computation
On instructions from the Bashkirskii Geophysical Trust, we computed the dependence of NGC readings on the porosity of limestone (pores filled with water) and the pore diameter. The volumetric porosity was m = 1, 2, 3, 5, 10, 15, 20, 40 and 100%. The pore diameters were d = 5%, 7% and 11% in. The pore was filled with salt water. During the computations, account was taken of the following important physical phenomena, not accessible to analytic theory: the spectrum of the polonium-beryllium source, the neutron energy losses due to inelastic and elastic collisions, the spectrum of the radiative capture of gamma quanta, the Compton effect and photo-effect for propagation of gamma quanta, the spectral sensitivity of the SI-4G counters, the energy dependence of all the effective neutron and quantum cross-sections, and so on. No account was taken of the following: (1) neutron thermalization and diffusion; (A neutron was assumed to be captured where it was stopped. This is possible because, when the porosity varies, the stopping parameters vary by several hundred percent, while the diffusion lengths, which are one order less than the stopping lengths, have virtually no effect on the migration, so that the NGC dependence on the hydrogen content is not distorted). (2) gamma quanta of nonelastic the radiative quanta;
scatter,
since they are small compared with
(3) the effect of vapour formation during propagation of the gamma quanta, the cross-sections of which are small in limestone at gamma quantum energies up to 10 MeV. So much for the physical assumptions.
We now turn to the computer program
Co~put~~ons
by the Monte Carlo method
99
organization. 1. The computation was performed simultaneously for all m and d (30 variants), so that, for instance, we computed the neutron travel and the new coordinates in the m and d cycle, then the scatter in the m and d cycle, and so on. Many of the qu~tities required for all m and d (such as In ,6 for the travel, Q’ from (1.1) and so on) were computed once. 2. In order to retain a single path as long as possible for all 30 variants, [41, the coordinate axes for the different pore diameters were chosen so that the instrument position was the same for all d, while the pores themselves were in contact (Fig. 1). Hence, till the initial piece of path went beyond 5%“, it was the same for all 30 variants; if it were left after 5%“, but before 7?4”, it was the if it left after 7%“, but before same for 20 variants (while 10 were br~ch~); llW”, it was the same for 10 variants. 3. Correlation of the neutron path directions of Section 2.
was realized by the algorithm
4. Every neutron path of type A was subjected to symmetry transformations relative to the collision points with calcium or oxygen in accordance with what was said in Section 3, on the basis of the following considerations. Most of the paths leaving a pore’fspecially the longest in the case of low porosities) have a tendency to go a long way into the stratum and to migrate little in the z direction. The problem therefore amounts to extracting information from them about the stopping field in the zone round the pore and for z close to the indicator. This was done as follows. Call the points of collision with heavy nuclei in the stratum Dj, and the points of radiative capture A,. During the drawing of the path, we stored the D, for which the z coordinate was a maximum or minimum. These were the points with respect to which we reflected the point A, once. As a result, 4 points were obtained, reducing the fluctuations of the neutron field along the z axis. If the neutron path went beyond the pore wall by more than the stopping length in the stratum (for the energy at this point), the coordinates of the nearest point Di were stored, and for the rest of the path, we analyzed its departure from the line through D i parallel to the pore axis. It was decided on this basis whether a rotation of A, about the given line was possible; the rotation was through 21718(eight points). If no rotation was possible, single reflections were performed with respect to two points, co~es~nding to the rn~irn~ and minimum angular departure from the plane (4 points). The line of rotation was chosen so that the rotation (like the 4 reflection) commuted with the z reflection, thus simplifying the logic. In this way, either
100
I. G. Dyad’kin
and V. N. Starikov
TABLE
m 5%
Theoret. rel.error I + 15%
100
0.15
40 20 15 10
0.42 0.45 0.46 0.45
Pore & lab. data/ r;@rror 0
0.44 0.51 0.56 0.65
m %
Theoret. rel.error I f 15%
Pore & lab. data,+r;k;rror 0
5
0.69
0.78
3 2 1 0.5
0.83 0.86 1 0.80
0.85 0.92 1
32 or 16 symmetric points with weights 2-’ or 2-’ were obtained as a result of the -symmetry transformations. As regards the type B and C paths, since they are very rarely encountered in the present problem in the region of the indicator, it was decided not to complicate the program by performing symmetry transformations with them. They were left with the weight 1. 5. A gamma quantum path has to be drawn from each point of radiative capture. This was done as follows. For each m variant, gamma quantum paths were drawn in the homogeneous medium corresponding to the stratum, and we analyzed its departure from the line parallel to the pore axis and passing through the point of quantum birth. This path was then “fitted” against all the symmetric points of radiative capture; if it did not cross the pore during the fitting an answer was obtained directly about its being missed; if it did cut the pore, the answer was obtained after redrawing the piece cutting the pore. The fitting was carried out, not only for the initial path, but also for the three paths rotated through 2~14 from it about the above-mentioned line. In this way the weight of a quantum of radiative capture was ?4 (we do not include here the weight of the neutron path and the weight corresponding to the number of quanta per act of neutron capture). If the neutron path terminated at the pore, 4 rotated quantum paths were drawn from it in the true geometry in the usual way. The computations were performed on the Razdan-111 computer (speed approx. 10’ op/sec) of the Academy of Sciences of the Armenian SSR and the Erevan University Computer Centre, during 100 hours of machine time (as much time was spent on harmonizing the program, dealing with the working block of 10,000 cells, and checking the control variants). The program operation was monitored, in particular, by the method of printing out the neutron field. Experimental data were only obtained for a single pore with d = 5%“, bored out with complete removal of the core. The table gives a comparison of the theoretical curves with these data, obtained by M. S. Pernikov. It throws light on the suitability of the proposed
Computations
101
by the Monte Carlo method
mathematical model of the NGC. All this work was performed on the initiative of I. L. Dvorkin, who carried out all the geophysical analysis of the statement of the problem and of the results obtained (which are not quoted here). The authors thank the ByeleRussian assistance with the computations.
State University students for
REFERENCES
1.
KAHN,
H.
Use of different Monte Carlo sampling techniques. N.Y. 146-191, 1956.
Symposium
on Monte
Carlo Methods, 2.
GOERTZEL,
A. and KALOS,
Nucl. Energy,
Ser I, Phys.
M. Monte Carlo methods in transport problems. Math. 2, 315-395, 1958.
Prog.
3.
FROLOV, A. S. and CHENTSOV, N. N. On computation by the Monte Carlo method of definite integrals depending on a parameter. Zh. uychisl. Mat. mat. Fiz. 2, 4, 714-717, 1962.
4.
DYAD’KIN, I. G. LISENENKOV, A. T. and PONYATOV, G. I. On the convergence rate of the Monte Carlo method when solving problems of radioactive logging. Zh. vychisl. Mat. mat. Fiz. 5, 4, 763-768, 1965.
5.
DYAD’KIN and PONYATOV, G. I. On the theory and method of mathematical modelling of gamma quanta propagation. lzu. Akad. Nauk SSSR, Fiz., Zemli, 1,47-58, 1966.
6.
DYAD’KIN, I. G. and STARIKOV, V. N. The possibility of economizing on machine time when solving Laplace’s equation by the Monte Carlo method. Zh. uychisl. Mat. mat. Fiz. 5, 936-938, 1965.
7.
BUSLENKO, N. P. and GOLENKO, D. I. et al. Monte Carlo techniques statisticheskikh ispytaniil. Moscow, Fizmatgiz, 1962.
8.
AKHIEZER, voprosy
9.
IVANENKO, D. and SOKOLOV, A. Class ical field polya). Moscow, Gostekhizdat, 1951.
(Meted
A. and POMERANCHUK, I. Some aspects of nuclear theory (Nekotorye teorii yadral. Moscow-Leningrad, Gostekhizdat, 1959.
theory (Klassicheskaya
teoriya
Translated by D. E. Brown