Research on GPU-accelerated algorithm in 3D finite difference neutron diffusion calculation method

Research on GPU-accelerated algorithm in 3D finite difference neutron diffusion calculation method 1 corresp first-author Corresponding author, q-xu09@mails.tsinghua.edu.cn Corresponding author, q-xu09@mails.tsinghua.edu.cn q-xu09@mails.tsinghua.edu.cn 1 1 1 Department of Engineering Physics, Tsinghua University, Beijing100084, China Department of Engineering Physics, Tsinghua University, Beijing100084, China In this paper, the adaptability of the neutron diffusion numerical algorithm on GPUs was studied, and a GPU-accelerated multi-group 3D neutron diffusion code based on finite difference method was developed. The IAEA 3D PWR benchmark problem was calculated in the numerical test. The results demonstrate both high efficiency and adequate accuracy of the GPU implementation for neutron diffusion equation. Neutron diffusion Finite difference Graphics Processing Unit (GPU) CUDA Acceleration I. INTRODUCTION 1 In the field of reactor physics, numerical solutions of 3-dimentional neutron diffusion equation are always required. Compared with the coarse mesh nodal techniques, the finite difference method is considered simpler and more precise, however, it costs unendurable computer time when analyzing a full-size reactor core. Since 2006, NVIDIA’s GPUs (Graphics Processing Units) has provided us with tremendous computational horsepower because of the release of CUDA [ 1 1 ]. In the field of nuclear reactor physics, the importance of the GPU+CPU heterogeneous platform has been growing gradually. Prayudhatama et al. [ 2 2 ] implemented a 1-D finite difference diffusion code on GPUs in 2010, and obtained up to 70× speedup compared to a corresponding CPU code. In 2011, Kodama et al. [ 3 3 ] ported the code SCOPE2 to GPUs, they got about 3 times speedup. In the same year, Gong et al. [ 4 4 ] exploited the parallelism of GPUs for the S n code Sweep3D, which was speeded up by about 2 to 8 times. In this work, a GPU-accelerated multi-group 3D neutron diffusion code based on finite difference method was implemented and optimized. The IAEA 3D PWR benchmark problem [ 5 5 ] was utilized to prove the high computational efficiency and accuracy of the GPU version code. The result in this work shows a bright future of GPU applications in nuclear reactor analysis. II. NEUTRON DIFFUSION EQUATION 1 According to the neutron diffusion theory, we have the multi-group neutron diffusion equation [ 6 6 ] as below, <math display="block" id="M1"><mrow><mo>−</mo><mo>∇</mo><mo>⋅</mo><msub><mi>D</mi><mi>g</mi></msub><mo>∇</mo><msub><mi>φ</mi><mi>g</mi></msub><mo>+</mo><msubsup><mi>Σ</mi><mi>t</mi><mi>g</mi></msubsup><mo stretchy="false">(</mo><mi>r</mi><mo stretchy="false">)</mo><msub><mi>φ</mi><mi>g</mi></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>Σ</mi><mi>s</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msubsup><msub><mi>φ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mo>+</mo><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mi>ν</mi><msubsup><mi>Σ</mi><mi>f</mi><msup><mi>g</mi><mo>′</mo></msup></msubsup><msub><mi>φ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mtext>,</mtext></mrow></math> where g is the energy group number, ranging from 1 to G , k eff is the effective multiplication factor, and φ g is the g th neutron flux. We solve this equation by the source iteration methodology [ 6 6 ], which includes inner and outer iterations. The inner iteration computes a group of linear algebra equations in the form of A X =B given the neutron scattering source and fission source. In the outer iteration, we use the neutron flux to update the neutron source and get ready for a next inner iteration. In this work, we focus on accelerating the inner iteration. Suppose the neutron source on the right hand of Eq. (1) is known, then after discretization on XYZ grid, we can get Eq. (2), which is a linear equation with a 7 diagonal positive definite matrix as the coefficient. <math display="block" id="M2"><mrow><mtable columnalign="left"><mtr columnalign="left"><mtd columnalign="left"><mrow><msub><mi>a</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mo>−</mo><mn>1,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub></mrow></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>b</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mo>−</mo><mn>1,</mn><mi>k</mi></mrow></msub><mo>+</mo><msub><mi>c</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi><mo>−</mo><mn>1</mn></mrow></msub><mo>+</mo><msub><mi>d</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>+</mo><msub><mi>e</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mo>+</mo><mn>1,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><mo>+</mo><msub><mi>f</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mo>+</mo><mn>1,</mn><mi>k</mi></mrow></msub><mo>+</mo><msub><mi>h</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><msub><mi>φ</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi><mo>+</mo><mn>1</mn></mrow></msub></mrow></mtd></mtr><mtr columnalign="left"><mtd columnalign="left"><mrow/></mtd><mtd columnalign="left"><mrow><mo>=</mo><msub><mi>s</mi><mrow><mi>i</mi><mn>,</mn><mi>j</mi><mn>,</mn><mi>k</mi></mrow></msub><mn>,</mn><mn>1</mn><mo>≤</mo><mi>i</mi><mo>≤</mo><mi>I</mi><mn>,</mn><mn>1</mn><mo>≤</mo><mi>j</mi><mo>≤</mo><mi>J</mi><mn>,</mn><mn>1</mn><mo>≤</mo><mi>k</mi><mo>≤</mo><mi>K</mi><mn>.</mn></mrow></mtd></mtr></mtable></mrow></math> Eq. (2) is a large-scale sparse matrix problem for full-size reactor analysis. In the perspective of numerical mathematics, Jacobi iteration is inefficient, thus, arithmetic techniques, such as CG, SOR, LSOR, ADI (Alternating Direction Implicit method), are needed for efficient calculation. III. GPU IMPLEMENTATION DETAILS 1 In order to test the computational capacity of GPUs, we do not resort to any mathematical skill; instead, the Jacobi iteration method is adopted for the inner iteration. The inner iteration is ported to a GTX TITAN GPU. In the inner iteration, the neutron flux is estimated according to the neutron source calculated with the neutron flux of the last source iteration. The outer iteration is still remained on CPUs to calculate the neutron effective multiplication factor according to the neutron fission source from GPUs. After an outer iteration, the effective multiplication factor is transferred from the CPU memory to the GPU memory to get the neutron source for the next source iteration. Fig. 1 Fig. 1 shows the tasks distribution and data transfer between GPUs and CPUs during one source iteration. A. Solving neutron flux 2 The neutron flux is solved via Jacobi inner iteration on GPUs. Because of the natural parallelism of the Jacobi iteration, it is of high possibility to implement this algorithm on GPUs with exciting speedups. According to CUDA, GPUs have two levels of parallelism, the first level is called grids of thread blocks, while the second level is blocks of threads. One thread block is designed to be mapped to an SM (Streaming Multiprocessor) on GPU chips, and one thread to be mapped to an SP (Streaming Processor) in an SM. In Jacobi iteration, the main part of computing tasks are production and addition operations at each flux point, which is shown by Eq. (2). To speed up such kind of iterations, the operations at each flux point ( i , j , k ) should be allocated to a specific GPU thread so that the computation tasks can be spread among the SPs on GPUs. Fig. 2 Fig. 2 demonstrates the mapping relationships between flux points and threads. As can be seen in Fig. 2 Fig. 2 , one flux point is mapped to one GPU thread, and each thread is responsible to update the flux at that point. For a large-scale 3D reactor model, there will be millions or even tens of millions of flux points needed to be updated using the surrounding old flux, however, the hardware resources of a GPU chip is limited to create as many threads as the flux points. To solve this problem, we update the neutron flux layer by layer as is illustrated in Fig. 2 Fig. 2 . After that, there will be enough computing resources for a GPU to accelerate the inner iteration procedure for each layer of flux points. B. Generating sources and data movements 2 When the neutron flux is solved after inner iterations, the fission source and the neutron source can be determined by the following equaptions: <math display="block" id="M3"><mrow><msub><mi>S</mi><mrow><mtext>fission</mtext></mrow></msub><mo>=</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><mi>ν</mi><msubsup><mi>Σ</mi><mi>f</mi><msup><mi>g</mi><mo>′</mo></msup></msubsup><msub><mi>φ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mtext>,</mtext></mrow></math> <math display="block" id="M4"><mrow><msub><mi>S</mi><mrow><mtext>neutro</mtext><mn>,</mn><mi>g</mi></mrow></msub><mo>=</mo><mfrac><mrow><msub><mi>χ</mi><mi>g</mi></msub></mrow><mrow><msub><mi>k</mi><mrow><mtext>eff</mtext></mrow></msub></mrow></mfrac><msub><mi>S</mi><mrow><mtext>fission</mtext></mrow></msub><mo>+</mo><mstyle displaystyle="true"><munderover><mo>∑</mo><mrow><msup><mi>g</mi><mo>′</mo></msup><mo>=</mo><mn>1</mn></mrow><mi>G</mi></munderover></mstyle><msubsup><mi>Σ</mi><mi>s</mi><mrow><msup><mi>g</mi><mo>′</mo></msup><mi>g</mi></mrow></msubsup><msub><mi>φ</mi><msup><mi>g</mi><mo>′</mo></msup></msub><mtext>,</mtext></mrow></math> where, S fission stands for the fission source and S neutron, g stands for the neutron source of energy group g , both of which are calculated by the neutron flux newly updated. In order to reduce data exchange between the CPU and GPU memories, these two sources are obtained on GPUs in parallel. As shown in Fig. 1 Fig. 1 , there are three data movements during one source iteration. The first data transfer happens after fission source was created, which moves fission source from device memory to host memory to calculate the effective multiplication factor k eff by accumulating the fission source of each flux point. The second data movement is for comparison between old and new neutron flux, during which the new neutron flux is transferred from device to host. The third one transfers k eff , a double type variable, back to device memory to get the neutron source. C. Data storage 2 For the fine grid finite difference method, large number of flux points lead to large memory space needs. Suppose there are N g energy groups, and N x , N y , N z flux points in the X, Y, Z direction respectively, then the memory space to store the eight coefficients (including the neutron source S i , j , k ) would be 4(bytes) × N g × N x × N y × N z × 8 bytes, and the memory for the neutron flux would be 4(bytes) × N g × N x × N y × N z bytes. When analyzing 3D full size reactors using GPUs, all the above data should be allocated to GPU memory, which is of limited volume. For GTX TITAN, the device memory is up to 6 GB under 64 bit operating systems. Under CUDA, a programmer is allowed to manage 7 different kinds of memory space, among which only the global memory and the texture memory are able to be utilized to store the coefficient data and the flux data. Because texture memory has a texture cache and higher bandwidth than global memory, it is advantageous to access data frequently from it. The only limitation of texture memory is that it is read-only. Thus, the coefficient data can be filled into texture memory, and the flux data be allocated to global memory. IV. PERFORMANCE TEST 1 In this section, we demonstrate the accuracy and efficiency of the GPU accelerated code. Besides, we also discuss a way of performance improvement by overclocking GPU processors. A. Experiment platform and benchmark problem 2 The accuracy of the GPU version diffusion code is tested by comparing the neutron flux computed by CITATION [ 7 7 ]. In order to prove the efficiency of the GPU code, we measure the performance of three diffusion codes listed in Table 1 Table 1 . 3DFD-CPU is a serial CPU version code which uses the Jacobi iteration method for inner iterations. 3DFD-GPU is obtained by accelerating the inner iteration part of 3DFD-CPU utilizing GPUs. HYPRE-8CORE [ 8 8 ] is a parallel diffusion code running on an 8-core CPU. The computing hardwares of these codes are also shown in Table 1 Table 1 . 1 1 Code 1 1 Computing platform 1 1 3DFD-CPU 1 1 CPU:Intel Quad Q9300 @ 2.50 GHz 1 1 DFD-GPU 1 1 GPU: NVIDIA Geforce GTX TITAN 1 1 HYPRE-8CORE 1 1 CPU:Intel Xeon E5520 @ 2.26 GHz The IAEA PWR benchmark problem, shown in Fig. 3 Fig. 3 , is used for the numerical experiment. This is an important benchmark problem widely used to test the performance of the neutron deterministic codes. The core is composed of 177 fuel assemblies, 9 of which are fully rodded and 4 of which are partially rodded. There are 64 reflector assemblies surrounding the core. The size of the assemblies is 20 cm×20 cm×340 cm, while the size of 1/4 core is 170 cm×170 cm×380 cm. B. Accuracy of the GPU code 2 To prove the accuracy of GPU computation, we compare the power distribution of 3DFD-CPU with that of CITATION. CITATION, developed by ORNL, is an industrial class code for solving the neutron diffusion equation. The comparison results are shown in Fig. 4 Fig. 4 . The convergence criterion is set so that the simulation comes to an end when the effective multiplication factor relative error is less than 1.0 × 10 -6 and the maximum point flux relative error is less than 1.0 × 10 -5 . The computing grid size used in Fig. 4 Fig. 4 is 2 cm, there are 1372750 spatial flux points. In Fig. 4 Fig. 4 , the relative error stands for the difference of the code result from the benchmark result. The power distribution of the GPU version code is close to that of CITATION. The accuracy comparison tells us that there is no need to worry about the accuracy and reliability of GPUs. C. Efficiency of the GPU code 2 We use the codes listed in Table 1 Table 1 to testify the computing power of GPUs. Firstly, 3DFD-GPU is compared with the 8-core CPU parallelized code HYPRE-8CORE, and then a comparison of computing time between 3DFD-GPU and 3DFD-CPU is made. According to Ref. [ 8 8 ], the author utilized the MPI-based parallelized linear algebra library HYPRE [ 9 9 ] to accelerate the diffusion code. Here we call the corresponding code in Ref. [ 8 8 ] as HYPRE-8CORE. HYPRE is a library developed by LLNL for solving large sparse linear systems of equations on massively parallel computers. On an 8-core tower server, the inner iteration part of the diffusion code is accelerated by the parallelized Conjugate Gradient algorithm. During simulation, the computing grid size is set to be 2.5 cm, and the convergence standard is that K eff relative error converges to 1.0×10 -5 and the maximum point flux relative error to 1.0 × 10 -4 . The computation speed comparison between 3DFD-GPU and HYPRE-8CORE is shown in Table 2 Table 2 . 1 1 Code 1 1 Computing time (s) 1 1 Speedup 1 1 HYPRE-8CORE 1 1 23.7 1 1 3.2 1 1 DFD-GPU 1 1 7.5 1 1 In Table 2 Table 2 , although the HYPRE-8CORE is accelerated by an 8-core server, 3DFD-GPU performs better. The performance comparison of 3DFD-GPU and 3DFD-CPU is shown in Fig. 5 Fig. 5 . We use six kinds of grid sizes, from 5 cm×5 cm×5 cm to 1 cm×1 cm×1 cm, to demonstrate the acceleration characteristic of GPUs for Jacobi iteration. Table 3 Table 3 lists the grid sizes and the corresponding grid numbers. The convergence criterion is that K eff relative error converges to 1.0× 10 -6 and the maximum point flux relative error to 1.0 × 10 -5 . 1 1 Grid Size (cm) 1 1 Grid Number 1 1 5 1 1 87856 1 1 10/3 1 1 296514 1 1 2.5 1 1 702848 1 1 2 1 1 1372750 1 1 1.25 1 1 5622784 1 1 1 1 1 10982000 Figure 5 Figure 5 shows the amazing accelerating power of GPUs compared with CPUs, especially when the grid size is set to 2 cm, a speedup factor of 86 was obtained. This phenomenon is caused by latency hiding, when the problem scales up and the amount of data increases, all cores on GPU are working at full capacity, then data transfer from the GPU memory by part of thread blocks can be operated while other blocks are executing the computational task. However, as can be seen from Fig. 5 Fig. 5 , it should be noticed that oversized data amount may decrease the speedups, because the communication overhead between the host and the device increases and the sequential part of the code may play an increasingly important role in the whole process. D. Performance improvement by overclocking 2 In order to get the same performance with lower energy consumption, NVIDIA decreased the base clock of GPUs of Kepler series, while increased the number of streaming processors in streaming multiprocessors (SMX). The base core clock of GTX TITAN is 837 MHz, which is lower than that of GTX 580 (Fermi architecture, 1544 MHz). We use the overclocking utility NVIDIA Inspector to set the core clock to be 1166 MHz and the memory clock to be 3334 MHz. Fig. 6 Fig. 6 shows the performance improvement after overclocking. In Fig. 6 Fig. 6 , the runtime and the speedup factor of 3DFD-GPU before and after overclocking are compared with each other, where the speedup factor is relative to the runtime of 3DFD-CPU. Through overclocking, the GPU acceleration effect is improved. The performance improvement depends on the scale of the analyzed problem, that is to say, more obvious performance enhancement can be obtained when the grid number increases. V. CONCLUSION 1 In this work, a GPU-accelerated multi-group 3D neutron diffusion code based on finite difference method was developed to speed up the finite difference methodology and examine the performance of GPUs. The IAEA 3D PWR benchmark problem is used as the problem model in the numerical experiment. By comparing the power distribution obtained from 3DFD-GPU and CITATION, we prove the accuracy of GPU computing. The performance advantage of GPUs is also demonstrated by comparing the runtime of 3DFD-GPU, 3DFD-CPU and HYPRE-8CORE. As to the future work, mathematical accelerating techniques, such as the Conjugate Gradient method and the Chebyshev extrapolation method, will be adopted to reduce the runtime of the GPU-based finite difference method to the same order of magnitude as the coarse mesh nodal methodology. References NVIDIA Corporation . CUDA C Programming Guide . 2012 , 3–4 and 71–75 . CUDA C Programming Guide Prayudhatama D , Waris A , Kurniasih N , et al . Proceedings of AIP Conference Proceedings , 2010 , 1244 : 121 - 126 . Kodama Y , Tatsumi M , Ohoka Y . Study on GPU Computing for SCOPE2 with CUDA . Proceedings of International Conference on Mathematics and Computational Methods Applied to Nuclear Science and Engineering (M&C2011) , Brazil , 2001 . Study on GPU Computing for SCOPE2 with CUDA Gong C , Liu J , Chi L , et al . GPU Accelerated Simulations of 3D Deterministic Particle Transport Using Discrete Ordinates Method . Journal of Computational Physics , 2011 , 230 ( 15 ): 6010 - 6022 . GPU Accelerated Simulations of 3D Deterministic Particle Transport Using Discrete Ordinates Method Argonne National Laboratory . Benchmark Problem Book. ANL-7416, Suppl.2 , 1977 , 277 - 280 . Benchmark Problem Book. ANL-7416, Suppl.2 Duderstadt J J and Hamilton L J , Nuclear Reactor Analysis . New York (USA) : John Wiley & Sons, Inc. , 1976 , 285 - 314 . Nuclear Reactor Analysis Fowler T B , Vondy D R , Cunningham G W . Nuclear Reactor Core Analysis Code CITATION, ORNL-TM-2496, Supplement 2 . ORNL , 1972 , 104 - 140 . Nuclear Reactor Core Analysis Code CITATION, ORNL-TM-2496, Supplement 2 Wu W B , Li Q , Wang K . Parallel Solution of 3D Neutron Diffusion Equation Based on HYPRE. Science and Technology on Reactor System Design Technology Laboratory Annual Report . Chengdu, China , 2010 , 35 - 40 (in Chinese) . Parallel Solution of 3D Neutron Diffusion Equation Based on HYPRE. Science and Technology on Reactor System Design Technology Laboratory Annual Report Lawrence Livermore National Laboratory . HYPRE User’s Manual (Version 2.7.0b) . 2011 , 1 - 6 . HYPRE User’s Manual (Version 2.7.0b) 10.13538/j.1001-8042/nst.25.010501 Supported by the 973 Program (No.2007CB209800) and National Natural Science Foundation of China (No.11105080) 2013-03-14 2013-09-20 2014-02-20 2014-02-20 2021-04-11 2014 Nuclear Science and Techniques 1 25