GPU parallel Monte Carlo light atom reaction sampling method

Abstract

The invention discloses a GPU (Graphics Processing Unit) parallel Monte Carlo light atom reaction sampling method, which comprises the following steps of: firstly, calculating a probability density function of scattering angle cosine based on a coherent scattering and Compton scattering differential section formula, and carrying out trapezoidal integration on the probability density function of the scattering angle cosine to obtain a probability distribution function of the scattering angle cosine; secondly, calculating probability distribution functions under all incident energy points to obtain a two-dimensional interpolation table of the probability distribution functions about incident energy and scattering angle cosine; and finally, sampling the scattering angle cosine in a Monte Carlo particle transport program, and calculating the energy of the secondary photons according to a sampling result. According to the GPU parallel Monte Carlo light atom reaction sampling method provided by the invention, the GPU parallel efficiency is greatly improved.

Description

Technical Field

[0001] The invention relates to the technical field of Monte Carlo particle transport calculation, and in particular to a GPU parallel Monte Carlo photoatom reaction sampling method. Background Art

[0002] The Monte Carlo particle transport method is a particle transport solution method based on probability theory. Because it can accurately simulate the physical process of particle transport in a medium, this method has been widely used in the past fifty years. The Monte Carlo particle transport method is random and requires the simulation of a large number of particle histories to achieve the desired statistical accuracy, which consumes a lot of computing resources. With the continuous development of GPUs and the natural suitability of the Monte Carlo method for parallelism, more and more people are developing Monte Carlo particle transport programs on GPUs.

[0003] However, the traditional Monte Carlo light atomic reaction sampling method is not suitable for GPU parallelism. The traditional Monte Carlo light atomic reaction sampling method has a large number of branch processes, and different threads may enter different code branches. However, threads on the same GPU core can only execute the same code, which makes the traditional method very inefficient in parallel on GPUs. Summary of the invention

[0004] In order to solve the problems existing in the above-mentioned prior art, the present invention proposes a GPU parallel Monte Carlo photoatom reaction sampling method, which is suitable for GPU parallelism and has a calculation accuracy comparable to that of traditional methods, solving the problem of low parallel efficiency of traditional sampling methods on GPUs.

[0005] To achieve the above object, the present invention adopts the following technical solutions:

[0006] A GPU parallel Monte Carlo photoatom reaction sampling method comprises the following steps:

[0007] Step 1: Calculate the probability density function of the scattering angle cosine by using the coherent scattering and Compton scattering differential cross section formulas, and discretize the probability density function into a two-dimensional interpolation table;

[0008] (1) Coherent scattering

[0009] The process of elastic scattering of photons by electrons bound in atomic nuclei is called coherent scattering. The coherent scattering differential cross section is calculated using the following formula:

[0010]

[0011] Where:

[0012] ——coherent scattering differential cross section / barn;

[0013] μ——cosine of scattering angle;

[0014] E——incident photon energy / MeV;

[0015] r e ——The classical radius of the electron, r e =2.81794×10 -13 cm

[0016] F(x,Z)——coherent scattering atomic form factor;

[0017] x——photon momentum transfer /

[0018] Z——atomic number.

[0019] The probability density function of coherent scattering is:

[0020]

[0021] Where:

[0022] P coh (E,μ)——Probability density function of coherent scattering;

[0023] σ coh (E)——Coherent scattering cross section / barn when the incident photon energy is E.

[0024] Among them, the photon momentum transfer amount x has the following relationship with the incident photon energy E and the cosine of the scattering angle μ:

[0025]

[0026] Where:

[0027] k is the ratio of incident photon energy to electron rest energy;

[0028] E——incident photon energy / MeV.

[0029] Combining formula (2) and formula (3), it can be seen that the probability density function of coherent scattering of the same nuclide is only related to the incident photon energy and the cosine of the scattering angle. Therefore, the incident photon energy grid (E1, E2, E3······E n ), scattering angle cosine grid (μ1,μ2,μ3······μ m ) is used to discretize the probability density function and use linear interpolation between grid points to obtain a two-dimensional interpolation table of the probability density function of coherent scattering.

[0030] (2) Compton scattering

[0031] The differential cross section for Compton scattering of photons with bound electrons is calculated using the following formula:

[0032]

[0033] Where:

[0034] ——Compton scattering differential cross section / barn;

[0035] k′——the ratio of the outgoing photon energy to the electron rest energy;

[0036] S(x,Z)——Compton scattering shape factor.

[0037] Among them, k, k′, μ have the following relationship:

[0038]

[0039] Substituting formula (5) into formula (4), we can get the probability density function of Compton scattering:

[0040]

[0041] Where:

[0042] P inc (E,μ)——Probability density function of Compton scattering;

[0043] σ inc (E)——Compton scattering cross section when the incident photon energy is E / barn;

[0044] Among them, the photon momentum transfer amount x has the following relationship with the incident energy E and the cosine of the scattering angle μ:

[0045]

[0046] Similarly, the probability density function of Compton scattering of the same nuclide is only related to the incident photon energy and the cosine of the scattering angle. Therefore, the incident photon energy grid (E1, E2, E3... n ), scattering angle cosine grid (μ1,μ2,μ3······μ m ) is used to discretize the probability density function and use linear interpolation between grid points to obtain a two-dimensional interpolation table of the probability density function of Compton scattering.

[0047] Step 2: Integrate the two-dimensional interpolation table of the probability density function of coherent scattering and Compton scattering with respect to the cosine of the scattering angle to obtain a two-dimensional interpolation table of the probability distribution function with respect to the incident photon energy and the cosine of the scattering angle;

[0048] The probability distribution function of coherent scattering and Compton scattering is obtained by integrating the probability density function. The calculation formula is as follows using trapezoidal integration:

[0049]

[0050] Where:

[0051] C coh (E,μ j ), C inc (E,μ j )——Probability distribution function of coherent scattering and Compton scattering;

[0052] P coh (E,μ i ), P inc (E,μ i )——Coherent scattering, Compton scattering in μ i The probability density function at ;

[0053] P coh (E,μ i+1 ), P inc (E,μ i+1 )——Coherent scattering, Compton scattering in μ i+1 Probability density function at

[0054] i, j – counting identifiers.

[0055] Step 3: Sample the incident photon energy and probability distribution function, and interpolate the secondary photon scattering angle cosine according to the sampling results. Then calculate the secondary photon outgoing energy through the relationship between the secondary photon scattering angle cosine and the outgoing energy;

[0056] First, sample the incident photon energy and choose which probability distribution function to use under which incident photon energy. For photons with incident photon energy E, if E≤E1, use E1; if E≥E n Then use E n ; If E∈(E i ,E i+1 ),have Draw a random number ξ1, if ξ1≤f then use E i+1 Otherwise, use E i .

[0057] After sampling the incident photon energy, the scattering angle cosine is sampled. Assume that the incident photon energy is E i , draw a random number ξ2, and C(E i ,μ j )<ξ2C(E i ,μ m)<C(E i ,μ j+1 ), the cosine of the scattering angle can be obtained as follows:

[0058]

[0059] Where:

[0060] μ——cosine of scattering angle;

[0061]

[0062] P(E i ,μ j )——The incident energy is E i The cosine of the scattering angle is μ j The probability density function of ;

[0063] C(E i ,μ j )——The incident energy is E i The cosine of the scattering angle is μ j The probability distribution function when ;

[0064] μ m ——Grid point with maximum cosine of scattering angle.

[0065] After sampling the scattering angle cosine, the secondary photon outgoing energy needs to be calculated based on the relationship between the scattering angle cosine and the outgoing photon energy. For coherent scattering: the secondary photon outgoing energy is the same as the incident energy, E′=E. For Compton scattering, the secondary photon energy is:

[0066] Compared with the prior art, the present invention has the following advantages:

[0067] (1) The prior art has a large number of branch processes and the parallel efficiency on the GPU is very low; the present invention has no branch processes and can achieve very high parallel efficiency on the GPU.

[0068] (2) The prior art has low sampling efficiency when the photon energy is low and requires multiple cycles of calculation; the present invention only requires one calculation to complete the sampling. BRIEF DESCRIPTION OF THE DRAWINGS

[0069] Figure 1 is a flow chart of the method of the present invention.

[0070] Figure 2 (a), (b), (c) and (d) are the one-dimensional spherical shell problems. 56 Fe, 98 Mo, 182 W. 90Zr is the material, and the photon leakage spectra are calculated using the traditional method and the present method.

[0071] Figure 3 (a), (b), (c), (d), (e) and (f) are the photon leakage spectra at the detector calculated using the traditional method and this method in the RFNC benchmark with Al, Cu, Fe, Pb, H2O and NaCl as materials, respectively.

[0072] Figure 4 It is the parallel efficiency calculated by statistically calculating the time required to complete 100 million samples when using 1, 10, 50, 100, 500, and 1000 threads. DETAILED DESCRIPTION

[0073] The present invention is further described in detail below with reference to the accompanying drawings and specific embodiments.

[0074] The present invention is a GPU parallel Monte Carlo photoatom reaction sampling method. Figure 1 As shown, the present invention comprises the following steps:

[0075] Step 1: Calculate the probability density function of the scattering angle cosine by using the coherent scattering and Compton scattering differential cross section formulas, and discretize the probability density function into a two-dimensional interpolation table;

[0076] (1) Coherent scattering

[0077] The process of elastic scattering of photons by electrons bound in atomic nuclei is called coherent scattering. The coherent scattering differential cross section is calculated using the following formula:

[0078]

[0079] Where:

[0080] ——coherent scattering differential cross section / barn;

[0081] μ——cosine of scattering angle;

[0082] E——incident photon energy / MeV;

[0083] r e ——The classical radius of the electron, r e =2.81794×10 -13 cm

[0084] F(x,Z)——coherent scattering atomic form factor;

[0085] x——photon momentum transfer /

[0086] Z——atomic number.

[0087] The probability density function of coherent scattering is:

[0088]

[0089] Where:

[0090] P coh (E,μ)——Probability density function of coherent scattering;

[0091] σ coh (E)——Coherent scattering cross section / barn when the incident photon energy is E.

[0092] Among them, the photon momentum transfer amount x has the following relationship with the incident photon energy E and the cosine of the scattering angle μ:

[0093]

[0094] Where:

[0095] k is the ratio of incident photon energy to electron rest energy;

[0096] E——incident photon energy / MeV.

[0097] Combining formula (2) and formula (3), it can be seen that the probability density function of coherent scattering of the same nuclide is only related to the incident photon energy and the cosine of the scattering angle. Therefore, appropriately selecting the incident photon energy grid (E1, E2, E3... n ), scattering angle cosine grid (μ1,μ2,μ3······μ m ) is used to discretize the probability density function and use linear interpolation between grid points to obtain a two-dimensional interpolation table of the probability density function of coherent scattering.

[0098] (2) Compton scattering

[0099] The differential cross section for Compton scattering of photons with bound electrons is calculated using the following formula:

[0100]

[0101] Where:

[0102] ——Compton scattering differential cross section / barn;

[0103] k′——the ratio of the outgoing photon energy to the electron rest energy;

[0104] S(x,Z)——Compton scattering shape factor.

[0105] Among them, k, k′, μ have the following relationship:

[0106]

[0107] Substituting formula (5) into formula (4), we can get the probability density function of Compton scattering:

[0108]

[0109] Where:

[0110] P inc (E,μ)——Probability density function of Compton scattering;

[0111] σ inc (E)——Compton scattering cross section when the incident photon energy is E / barn;

[0112] Among them, the photon momentum transfer amount x has the following relationship with the incident energy E and the cosine of the scattering angle μ:

[0113]

[0114] Similarly, the probability density function of Compton scattering of the same nuclide is only related to the incident photon energy and the cosine of the scattering angle. Therefore, the incident photon energy grid (E1, E2, E3... n ), scattering angle cosine grid (μ1,μ2,μ3······μ m ) is used to discretize the probability density function and use linear interpolation between grid points to obtain a two-dimensional interpolation table of the probability density function of Compton scattering.

[0115] Step 2: Integrate the two-dimensional interpolation table of the probability density function of coherent scattering and Compton scattering with respect to the cosine of the scattering angle to obtain a two-dimensional interpolation table of the probability distribution function with respect to the incident photon energy and the cosine of the scattering angle;

[0116] The probability distribution function of coherent scattering and Compton scattering is obtained by integrating the probability density function. The calculation formula is as follows using trapezoidal integration:

[0117]

[0118] Where:

[0119] C coh (E,μ j ), C inc (E,μ j )——Probability distribution function of coherent scattering and Compton scattering;

[0120] P coh (E,μ i ), P inc (E,μi )——Coherent scattering, Compton scattering in μ i The probability density function at ;

[0121] P coh (E,μ i+1 ), P inc (E,μ i+1 )——Coherent scattering, Compton scattering in μ i+1 Probability density function at

[0122] i, j – counting identifiers.

[0123] Step 3: Sample the incident photon energy and probability distribution function, and interpolate the secondary photon scattering angle cosine according to the sampling results. Then calculate the secondary photon outgoing energy through the relationship between the secondary photon scattering angle cosine and the outgoing energy;

[0124] First, sample the incident photon energy and choose which probability distribution function to use under which incident photon energy. For photons with incident photon energy E, if E≤E1, use E1; if E≥E n Then use E n ; If E∈(E i ,E i+1 ),have Draw a random number ξ1, if ξ1≤f then use E i+1 Otherwise, use E i .

[0125] After sampling the incident photon energy, the scattering angle cosine is sampled. Assume that the incident photon energy is E i , draw a random number ξ2, and C(E i ,μ j )<ξ2C(E i ,μ m )<C(E i ,μ j+1 ), the cosine of the scattering angle can be obtained as follows:

[0126]

[0127] Where:

[0128] μ——cosine of scattering angle;

[0129]

[0130] P(E i ,μ j )——The incident energy is E i The cosine of the scattering angle is μ j The probability density function of ;

[0131] C(E i ,μ j )——The incident energy is E i The cosine of the scattering angle is μ j The probability distribution function when ;

[0132] μ m ——Grid point with maximum cosine of scattering angle.

[0133] After sampling the scattering angle cosine, the secondary photon outgoing energy needs to be calculated based on the relationship between the scattering angle cosine and the outgoing photon energy. For coherent scattering: the secondary photon outgoing energy is the same as the incident energy, E′=E. For Compton scattering, the secondary photon energy is: Calculation accuracy verification

[0134] In order to verify the accuracy of the calculation of the GPU parallel Monte Carlo photoatom reaction sampling method proposed in the present invention, the one-dimensional spherical shell problem and the RFNC benchmark problem are calculated using the method of the present invention and the traditional method respectively, and the photon leakage spectra obtained by the two methods are compared.

[0135] One-dimensional spherical shell problem: The inner diameter of the spherical shell is 10 cm, the outer diameter is 19.95 cm, and an isotropic neutron source with an energy of 14.1 MeV is given at the center of the spherical shell. The material and density filled in the spherical shell are shown in Table 1. The calculation results are shown in Figure 2 As shown in (a), (b), (c) and (d), the results show that the calculation results of the method of the present invention are basically consistent with those of the traditional method, the deviation is mainly distributed around 1%, and the maximum deviation is -6.63%. The reason for the deviation is the error caused by discretizing the probability density function into a two-dimensional interpolation table, and the numerical error when the probability density function is numerically integrated to obtain the probability distribution function.

[0136] Table 1 Test materials and densities for one-dimensional spherical shell problems

[0137]

[0138] RFNC benchmark question: The design basis of the RFNC benchmark question is a series of experiments conducted by the Russian Federal Nuclear Center. This benchmark question places an isotropic neutron source with an energy of 14.04981MeV at the center of the sample, and the detector is located 850cm away from the neutron source. A metal rod 40cm long and 3cm in diameter is placed between the neutron source and the detector to reduce the background interference of photons. Spherical samples of six materials, Al, Cu, Fe, Pb, H2O, and NaCl, were selected to calculate the photon leakage spectrum at the detector. The density of each sample material is shown in Table 2. The calculation results are shown in Figure 3As shown in (a), (b), (c), (d), (e) and (f), it can be seen from the figure that the calculation results of the improved new method of the present invention are well consistent with those of the traditional method, with the maximum deviation being -5.129%. This shows that the calculation accuracy of the photoatom reaction Monte Carlo sampling method of the present invention is comparable to that of the traditional method.

[0139] Table 2 Main nuclides and mass density of each sample of RFNC benchmark

[0140]

[0141] Parallel efficiency verification

[0142] In order to prove that the method of the present invention solves the problem of low parallel efficiency of traditional methods on GPU, the parallel efficiency of the program developed based on the sampling method proposed in the present invention is tested on GPU. The hardware platform parameters of the program test are shown in Table 3.

[0143] Table 3 Computing platform hardware and software parameters

[0144]

[0145] The sampling process was simulated 100 million times. The photon energy was randomly generated in the energy range of 0.001 to 100 MeV using a random generator. The time required to complete 100 million samplings was calculated using 1, 10, 50, 100, 500, and 1000 threads. The parallel efficiency calculated is as follows: Figure 4 As shown in the figure, it can be seen that the method of the present invention has a very high parallel efficiency on the GPU, and still has a parallel efficiency of 85.6% when using 1000 threads for parallel computing.

Claims

1. A GPU parallel Monte Carlo photoatom reaction sampling method, characterized by: The following steps are involved: Step 1: Calculate the probability density function of the scattering angle cosine by using the coherent scattering and Compton scattering differential cross section formulas, and discretize the probability density function into a two-dimensional interpolation table; (1) Coherent scattering The process of elastic scattering of photons by electrons bound in atomic nuclei is called coherent scattering. The coherent scattering differential cross section is calculated using the following formula: Where: μ——cosine of scattering angle; E——incident photon energy / MeV; r e ——The classical radius of the electron, r e =2.81794×10 -13 cm F(x,Z)——coherent scattering atomic form factor; x——photon momentum transfer / Z——atomic number; The probability density function of coherent scattering is: Where: P coh (E,μ)——Probability density function of coherent scattering; σ coh (E)——coherent scattering cross section when the incident photon energy is E / barn; Among them, the photon momentum transfer amount x has the following relationship with the incident photon energy E and the cosine of the scattering angle μ: Where: k is the ratio of incident photon energy to electron rest energy; E——incident photon energy / MeV; Combining formula (2) and formula (3), it can be seen that the probability density function of coherent scattering of the same nuclide is only related to the incident photon energy and the cosine of the scattering angle; therefore, the incident photon energy grid (E1, E2, E3······E n ), scattering angle cosine grid (μ1,μ2,μ3······μ m ) discretizes the probability density function and uses linear interpolation between grid points to obtain a two-dimensional interpolation table of the probability density function of coherent scattering; (2) Compton scattering The differential cross section for Compton scattering of photons with bound electrons is calculated using the following formula: Where: ——Compton scattering differential cross section / barn; k′——the ratio of the outgoing photon energy to the electron rest energy; S(x,Z)——Compton scattering shape factor; Among them, k, k′, μ have the following relationship: Substituting formula (5) into formula (4), the probability density function of Compton scattering is: Where: Among them, the photon momentum transfer amount x has the following relationship with the incident energy E and the cosine of the scattering angle μ: Similarly, the probability density function of Compton scattering of the same nuclide is only related to the incident photon energy and the cosine of the scattering angle; therefore, the incident photon energy grid (E1, E2, E3······E n ), scattering angle cosine grid (μ1,μ2,μ3······μ m ) discretizes the probability density function and uses linear interpolation between grid points to obtain a two-dimensional interpolation table of the probability density function of Compton scattering; Step 2: Integrate the two-dimensional interpolation table of the probability density function of coherent scattering and Compton scattering with respect to the cosine of the scattering angle to obtain a two-dimensional interpolation table of the probability distribution function with respect to the incident photon energy and the cosine of the scattering angle; The probability distribution function of coherent scattering and Compton scattering is obtained by integrating the probability density function. The calculation formula is as follows using trapezoidal integration: Where: C coh (E,μ j ), C inc (E,μ j )——Probability distribution function of coherent scattering and Compton scattering; P coh (E,μ i ), P inc (E,μ i )——Coherent scattering, Compton scattering in μ i The probability density function at ; P coh (E,μ i+1 ), P inc (E,μ i+1 )——Coherent scattering, Compton scattering in μ i+1 Probability density function at i, j – counting identifiers; Step 3: Sample the incident photon energy and probability distribution function, interpolate the sampling results to obtain the cosine of the secondary photon scattering angle, and then calculate the secondary photon outgoing energy through the relationship between the cosine of the secondary photon scattering angle and the outgoing energy; First, sample the incident photon energy and choose the probability distribution function under which incident photon energy to use; for photons with incident photon energy E, if E≤E1, use E1; if E≥E n Then use E n ; If E∈(E i ,E i+1 ),have Draw a random number ξ1, if ξ1≤f then use E i+1 Otherwise, use E i ; After sampling the incident photon energy, the scattering angle cosine is sampled; assuming that the incident photon energy is E i , draw a random number ξ2, and C(E i ,μ j )<ξ2C(E i ,μ m )<C(E i ,μ j+1 ), the cosine of the scattering angle is given by: Where: μ——cosine of scattering angle; P(E i ,μ j )——The incident energy is E i The cosine of the scattering angle is μ j The probability density function of ; C(E i ,μ j )——The incident energy is E i The cosine of the scattering angle is μ j The probability distribution function when ; μ m ——The maximum grid point of the cosine of the scattering angle; After completing the sampling of the scattering angle cosine, the secondary photon outgoing energy needs to be calculated through the relationship between the scattering angle cosine and the outgoing photon energy; for coherent scattering: the secondary photon outgoing energy is the same as the incident energy E′=E; for Compton scattering, the secondary photon energy is:

Images