Parallel feature line method calculation and storage optimization method for domestic heterogeneous architecture

Abstract

The invention provides a parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture, and belongs to the cross technical field of neutron physics and high-performance calculation. The method comprises the following steps: according to a reactor geometric construction numerical model, carrying out grid division to obtain flat source regions, and determining related data such as a filling material, a volume and a centroid of each flat source region; in the transportation calculation pretreatment stage, a 2D cyclic trajectory chain is generated on a radial plane, geometric information in the axial direction is searched for to generate an axial extrusion flat source area, 3D trajectory generation and related index information segmentation are executed and generated for an equipment end, and calculation information needed for executing a kernel function is copied to the equipment end; in the transport calculation stage, a kernel function is started at a host end, a device end sequentially executes flat source region source item updating, transport scanning calculation, standard flux updating and effective neutron multiplication factor and residual error calculation, and a PLA and TLA parallel algorithm is provided for transport scanning calculation. And copying the calculation result of the equipment end back to the host end, judging whether convergence is carried out or not, if not, returning to transport calculation, otherwise, ending iteration. Meanwhile, a transport scanning heterogeneous parallel performance model is provided, and an optimization scheme is given from the aspects of data transmission, storage overhead and the like according to the performance model.

Description

Technical Field

[0001] The present invention relates to the cross-technical field of neutron physics and high-performance computing, and in particular to a parallel characteristic line method computing and storage optimization method for a domestically produced heterogeneous architecture. Background Art

[0002] Neutron transport calculation is an important part of numerical reactor simulation. Studying the steady-state neutron transport problem, that is, solving the linear Boltzmann equation (BTE), can obtain the flux distribution of neutrons inside the reactor, which is of great significance for the initial reactor design and the prediction of neutron movement behavior inside the reactor.

[0003] The method of characteristics (MOC) is a common deterministic numerical method for solving the linear Boltzmann equation. Its solution idea is to transform the partial differential equation into an ordinary differential equation. The current new generation of reactor technology is becoming more and more complex, and the numerical simulation requires high fidelity, high precision and high efficiency. The direct three-dimensional method of characteristics (3D MOC) can realize the core pin-by-pin calculation and has high simulation accuracy.

[0004] With the refinement of simulation parameters and the scale of calculation examples, the direct three-dimensional characteristic line method for core transport solution shows the characteristics of huge storage and calculation. With the powerful computing power support, the supercomputing platform can meet the calculation and storage requirements of the program when calculating complex reactor problems, making the three-dimensional characteristic line neutron transport program have broad application prospects in engineering.

[0005] With the rapid development of domestic heterogeneous multi-core processors in recent years, heterogeneous node architecture has gradually become the dominant paradigm of high-performance computing platforms. Taking the CPU-DCU heterogeneous architecture as an example, DCU is a general-purpose coprocessor designed based on a large-scale parallel computing microstructure, and it cannot run calculations independently. However, it has strong parallel processing capabilities, high computing energy efficiency, and integrated on-chip high-bandwidth memory chips, with large storage bandwidth, and powerful floating-point computing capabilities. It is suitable for assisting the CPU in executing computationally intensive programs, and provides excellent parallel data processing capabilities in large-scale data calculations. It is widely used in high-performance computing. However, its unique hardware characteristics bring challenges to program transplantation and optimization, which are usually reflected in parallel algorithm design, reducing storage overhead, and optimizing program performance. Based on the above content, the present invention proposes a parallel characteristic line method calculation and storage optimization method for domestic heterogeneous architecture. Summary of the invention

[0006] The purpose of the present invention is to propose a parallel characteristic line method calculation and storage optimization method for a domestic heterogeneous architecture to solve the problems raised in the background technology; the present invention is aimed at a heterogeneous node architecture in which a CPU is equipped with a domestic DCU chip, and can use the domestic DCU chip to accelerate the three-dimensional characteristic line neutron transport calculation process while ensuring the calculation accuracy, improve the parallel computing efficiency, and reduce storage overhead.

[0007] In order to achieve the above object, the present invention adopts the following technical solutions:

[0008] A parallel characteristic line method calculation and storage optimization method for domestic heterogeneous architecture, including:

[0009] According to the numerical model of the reactor geometry, meshing is performed to obtain the flat source area and determine the filling material, volume, centroid and other related data of each flat source area;

[0010] In the pre-processing stage of transport calculation, a 2D loop trajectory chain is generated on the radial plane, the geometric information on the axial direction is searched to generate the axial extrusion flat source area, 3D trajectory and segmentation related index information are generated for the device side execution, and the calculation information required for executing the kernel function is copied to the device side;

[0011] In the transport calculation phase, the kernel function is started on the host side, and the device side sequentially performs source term update of the flat source area, transport scanning calculation, standard flux update, and calculation of effective neutron multiplication factor and residual;

[0012] The calculation results on the device side are copied back to the host side to determine whether they have converged. If not, return to the previous step. Otherwise, the iteration ends.

[0013] Furthermore, the solution of the source term update of the flat source region in each flat source region and each energy group is relatively independent. The total number of threads that can be run in parallel is the product of the number of flat source regions and the number of energy groups. The source term calculation of a single energy group in a single flat source region is assigned to a DCU thread (this parallel method can also be used to update the scalar flux). The corresponding kernel function needs to solve the total fission rate and the total scattering rate of the flat source region, and calculate the source term of the flat source region under a given energy group.

[0014] Furthermore, each transport scan needs to calculate all the fluxes of the trajectory under each energy group, and derive a new flux to contribute to each FSR. A single transport scan involves five nested loops: azimuth, trajectory under azimuth, line segment on trajectory, energy group, and polar angle. The transport scan has the largest memory access and computational workload in the entire source iteration implementation. A large number of trajectories combined with angles, energy groups, etc. can provide excellent parallelism. Based on this, the SIMT thread parallelism provided by DCU can be used to perform large-scale parallel acceleration of transport scan calculations.

[0015] In one loop of the nested loop scanning, the corresponding kernel function needs to allocate computing tasks to the DCU thread according to different algorithms (PLA and TLA), complete the 3D trajectory generation and segmentation, and perform transport scanning along the trajectory to calculate the trajectory incident neutron angular flux, outgoing neutron angular flux and the difference between the two, and finally update the boundary flux according to the boundary conditions.

[0016] Furthermore, when calculating the angular flux variation of each 3D trajectory in the flat source area, the following two solutions are provided according to the different parallel granularity sizes:

[0017] The plane-level parallel algorithm (PLA) specifically refers to a DCU thread being responsible for trajectory generation segmentation, angular flux calculation and contribution value of a plane.

[0018] The Track Level Parallel Algorithm (TLA) specifically refers to a DCU thread being responsible for a 3D trajectory generation segment as well as angular flux calculation and contribution value.

[0019] In different parallel algorithm implementations, the logic of generating a single 3D trajectory remains unchanged, and the only change is the number.

[0020] Furthermore, the time performance of the five functions of the transport scan, namely, trajectory generation, segmentation, calculation of forward and reverse angular flux and contribution value of line segments, forward and reverse boundary flux of the transmission trajectory, and copying data between the host and the device, is evaluated to quantify the calculation and communication time, including:

[0021] For serial programs, the total serial scan time formula is:

[0022]

[0023] in, Refers to the time spent calculating the attribute values ​​of a single 3D trajectory. Refers to the time spent calculating the attribute values ​​of a single line segment. Refers to calculating the change in angular flux on a single energy group on a single line segment. Refers to the calculation of the contribution of the angular flux change on a single line segment and a single energy group to the standard flux in the flat source area. It is 0 because no floating point addition, subtraction, multiplication or other calculation operations are required when transmitting boundary flux. Refers to the total communication time of information related to a single 3D trajectory. Refers to the total communication time for accessing the line segment information on the projected 2D trajectory in the storage space, the axial extrusion flat source area, the flat source area centroid, and the material. Refers to the total communication time for accessing the line segment starting point angular flux, source term, material cross-section data and index calculation table lookup in the storage space when calculating the angular flux of a single energy group of a single line segment. Refers to the total communication time for accessing the weight information and the standard flux in the storage space, Refers to the communication time to access the boundary flux;

[0024] For parallel programs, if the total number of threads is equal to the number of Planes, for the PLA algorithm, the number of 3D trajectories processed by each thread is not exactly the same as the number of threads, so the total transport scanning time calculation formula is:

[0025]

[0026] If the total number of threads is equal to the number of 3D trajectories, for the TLA algorithm, the total number of 3D line segments on a single trajectory is n sdseg =1,n egrp Indicates the number of benchmark example energy groups; the number of 3D line segments processed by each thread in the TLA algorithm is not exactly the same, so the total transport scan time calculation formula is:

[0027]

[0028] Among them, T h2d represents the time it takes to copy data from the host to the device; λ∈(0,p) represents the factor that weakens the acceleration effect due to wavefront divergence; T comp represents the calculation time; p represents the total number of computing units used for floating-point calculations; T idle Indicates waiting time; T comm Indicates memory access communication; T h2d Indicates the time it takes for data to be copied from the host to the device; N thd Indicates the total number of threads; n egrp Indicates the number of benchmark example energy groups; t1′ comp Refers to the time it takes for a thread in a parallel program to calculate the attribute values ​​of a single 3D trajectory; t3′ comp Refers to the time taken by a thread in a parallel program to calculate the change in angular flux on a single energy group on a single line segment; t1′ comm Refers to the total communication time of the information related to a single 3D trajectory within a thread in a parallel program; t3′ comm Refers to the total communication time for a thread in a parallel program to access the line segment starting point angular flux, source term, material cross-section data and index calculation table lookup in the storage space when calculating the angular flux of a single energy group of a single line segment; ξ=0 means that the memory access time is 0; ξ=1 means that each thread accesses the memory completely serially;

[0029] Furthermore, for the T in the performance formula h2d 、T d2h, and proposed a data transmission optimization plan. When executing the transport scanning kernel function, it is necessary to copy the 2D trajectory, 2D line segment, 2D trajectory chain, axial flat source area information and 3D trajectory index information. Although the arrays storing this information do not occupy a large amount of memory, most of them are not one-dimensional arrays. It is necessary to call the hipMemcpyHtoD copy function multiple times, which causes the API startup overhead to be higher than the data transmission overhead itself. By changing the multi-dimensional data to one-dimensional data by adding the index array method on the host side, the number of times the copy function is opened is reduced, and the data transmission time is reduced.

[0030] Furthermore, a register compilation optimization scheme is proposed for ξ in the performance formula. The _launch_bounds compiler instruction provided by HIP is used to affect the register usage of the kernel function, ensure that the appropriate number of workgroups and wavefronts reside on the CU, and hide the overhead by switching execution units through hardware.

[0031] Furthermore, for t' in the performance formula comm , a solution to reduce memory access overhead is proposed. The geometric boundary information used during trajectory generation and the weight information of the product group when calculating the angular flux contribution value are unchanged during the transport scan. Moreover, this part of information is small, and the threads in a wavefront almost access the same location. Putting this part of information in the constant memory of the DCU has lower memory access latency than global storage.

[0032] Furthermore, a pre-allocation space scheme is proposed. Each thread on the device side pre-allocates the space for the maximum number of segments that can be processed in the video memory. These spaces are partially redundant because the trajectory with the largest number of segments usually only occupies a small part. A direct mapping method is used to calculate the maximum number of segments in the trajectory that is offset by an integer multiple of the number of threads, that is, the maximum number of segments in the trajectory that this thread can process, rather than the maximum number of segments in the entire trajectory. The number of segments that each thread can process is summed up, the space pre-allocated in the video memory, and then a one-dimensional pointer array is used to store the starting position of each thread's allocated space according to the thread number.

[0033] Compared with the prior art, the present invention provides a parallel characteristic line method calculation and storage optimization method for a domestic heterogeneous architecture, which has the following beneficial effects:

[0034] The present invention designs transport scanning parallel algorithms PLA and TLA in view of the inherent parallelism characteristics of neutron transport programs when solving them. Starting from the performance analysis of the serial running time of transport scanning, a transport scanning heterogeneous parallel performance model is obtained to guide the optimization of heterogeneous algorithms. The present invention improves the utilization rate of hardware resources through register compilation optimization technology. At the same time, based on the merged memory access and shared storage technology, the memory access delay of the neutron transport scanning algorithm is reduced, and the variables that are reused in the algorithm and have low access transaction bandwidth efficiency are placed in the shared storage body. In order to reduce the phenomenon of high number of copy function calls but low data volume when data is transmitted between the host end and the device end, multidimensional data is encapsulated into one-dimensional data to reduce the number of copy function calls. The redundancy of the storage trajectory and segment space of PLA and TLA is analyzed, the data structure is adjusted, and the space utilization is optimized. BRIEF DESCRIPTION OF THE DRAWINGS

[0035] Figure 1 It is a flowchart of the parallel characteristic line method calculation and storage optimization method for domestic heterogeneous architecture proposed in Example 1 of the present invention;

[0036] Figure 2 This is a schematic diagram of the PLA algorithm thread tasks proposed in Example 1 of the present invention;

[0037] Figure 3 This is a schematic diagram of the TLA algorithm thread tasks proposed in Example 1 of the present invention;

[0038] Figure 4 This is a schematic diagram of a data transmission optimization solution proposed in Embodiment 1 of the present invention;

[0039] Figure 5 This is a schematic diagram of the domestic DCU multi-level storage system proposed in Example 1 of the present invention;

[0040] Figure 6 This is a schematic diagram of the track starting point merge access proposed in Embodiment 1 of the present invention;

[0041] Figure 7 Schematic diagram of the change of line segment storage space before and after the optimization proposed in Example 1 of the present invention. DETAILED DESCRIPTION

[0042] The technical solutions in the embodiments of the present invention will be described clearly and completely below in conjunction with the drawings in the embodiments of the present invention. Obviously, the described embodiments are only part of the embodiments of the present invention, rather than all the embodiments.

[0043] In order to better understand the present invention, the terms involved in the present invention are briefly explained as follows:

[0044] 1) Trajectory

[0045] A trajectory is a ray with start coordinates, end coordinates and direction (forward and backward), generated at a certain angle.

[0046] 2) Track chain

[0047] A trajectory chain refers to a ray loop formed by connecting trajectories end to end under corresponding boundary conditions. For example, a trajectory chain formed under total reflection boundary conditions is called a reflective trajectory chain.

[0048] 3) Pingyuan District

[0049] The commonly used approximation of the characteristic line method is called the flat source region approximation (or flat source approximation). It is assumed that the source term is constant in a discrete space unit. In this case, this discrete space unit is called a flat source region (FSRs), which means that the source term does not change from entering the flat source region to leaving it along the characteristic line.

[0050] 4) Axial extrusion geometry

[0051] Ray tracing is performed simultaneously on different radial planes. The planes that have completed ray tracing are superimposed together to obtain a plane containing all radial information, called an overlapping plane. The 2D geometric part on the overlapping plane is equivalent to the axial extrusion area, which has a corresponding unique identifier. It contains an axial grid and an FSRs array.

[0052] 5) HIP programming

[0053] The HIP parallel programming model is the upper layer of the Rocm ecosystem launched by AMD. It can provide a portable C / C++ heterogeneous parallel programming environment for domestic DCUs. The HIP programming model is very similar to CUDA in terms of interface and language specifications.

[0054] 6) SIMT

[0055] Both GPU and DCU are in single instruction multiple threads mode, that is, only the same instructions can be executed on each ALU.

[0056] 7) wavefront

[0057] The smallest unit of hardware scheduling in DCU is similar to the warp in AMD GPU, which contains a total of 64 threads. Threads in the same wavefront execute the same instructions on different data.

[0058] Based on the above content, the present invention describes the parallel characteristic line method calculation and storage optimization method of the domestic heterogeneous architecture proposed by the present invention in combination with the accompanying drawings and example data. The specific contents are as follows:

[0059] See also Figure 1 The embodiment of the present invention provides a parallel characteristic line method calculation and storage optimization method for domestic heterogeneous architecture, including:

[0060] According to the numerical model of the reactor geometry, meshing is performed to obtain the flat source area and determine the filling material, volume, centroid and other related data of each flat source area;

[0061] In the pre-processing stage of transport calculation, a 2D loop trajectory chain is generated on the radial plane, the geometric information on the axial direction is searched to generate the axial extrusion flat source area, 3D trajectory and segmentation related index information are generated for the device side execution, and the calculation information required for executing the kernel function is copied to the device side;

[0062] In the transport calculation phase, the kernel function is started on the host side, and the device side sequentially performs source term update of the flat source area, transport scanning calculation, standard flux update, and calculation of effective neutron multiplication factor and residual;

[0063] Copy the calculation results of the device side back to the host side to determine whether it has converged. If not, return to the previous step. Otherwise, the iteration ends.

[0064] The parallel characteristic line method calculation and storage optimization method of the domestic heterogeneous architecture proposed in the embodiment of the present invention, according to the numerical model of the reactor geometric structure, meshing is performed to obtain the flat source area and determine the filling material, volume, centroid and other related data of each flat source area; in the transport calculation pre-processing stage, a 2D loop trajectory chain is generated on the radial plane, the geometric information on the axial direction is searched to generate the axial extrusion flat source area, and the 3D trajectory and segmentation generation related index information are generated for the device end to execute, and the calculation information required for executing the kernel function is copied to the device end; in the transport calculation stage, the kernel function is started on the host end, and the device end sequentially performs the source term update of the flat source area, the transport scanning calculation, the update of the scalar flux and the calculation of the effective neutron increment factor and the residual, wherein the PLA and TLA parallel algorithms are proposed for the transport scanning calculation. The calculation results of the device end are copied back to the host end to determine whether they converge. If they do not converge, they return to the previous step, otherwise the iteration ends. In this way, the three-dimensional characteristic line neutron transport calculation process can be accelerated by using the domestic DCU chip under the premise of ensuring the calculation accuracy, and a transport scanning heterogeneous parallel performance model is proposed, and an optimization scheme is given from the aspects of data transmission and storage overhead according to the performance model.

[0065] In this embodiment, the heterogeneous parallel program consists of two parts, namely the host code running on the CPU side and the device code running on the DCU side, wherein the host side is responsible for cleaning and starting the kernel function running on the device side, and the device side is responsible for executing the kernel function to complete the computing tasks assigned by the host, that is, completing the neutron transport calculation according to the relevant information given by the host side transport calculation pre-processing process.

[0066] The main logic of the transport calculation is to assign an initial value of 1 to all source items on the grid during the first iterative solution. After calculating the scalar flux on all grids, use the scalar flux to calculate the source item Q under each FSR and energy group. After that, each source iteration calculation needs to calculate all equations in the iterative equation group. Calculate the change in angular flux of each trajectory in the flat source area, and then update the scalar flux with the change in angular flux passing through the same flat source area. Repeat the iteration until the final result meets the convergence rule. Figure 1 As shown in the figure, the source term update, transport scanning, standard flux update and effective value-added factor calculation in the flat source area are completed by three kernel functions: kernel 1, kernel 2 and kernel 3 respectively.

[0067] The source term update of the flat source area is solved relatively independently in each flat source area and each energy group. The total number of threads that can be run in parallel is the product of the number of flat source areas and the number of energy groups. The source term calculation of a single energy group in a single flat source area is assigned to a DCU thread. The execution steps of the corresponding kernel function kernel 1 are as follows:

[0068] A1, distributes computational load to threads.

[0069] A2, determine that the current thread is responsible for the flat source region i and energy group g, and then solve the total fission rate and the total scattering rate of this flat source region.

[0070] A3, based on the assumption of isotropic scattering, calculate the source term of the flat source region i under energy group g.

[0071] A4, arranges a new computing task for the current thread.

[0072] The parallel algorithm for scalar flux update is similar to that for source term update.

[0073] In this embodiment, each transport scan needs to calculate all fluxes of the trajectory under each energy group, and derive a new flux to contribute to each flat source region. A single transport scan involves five nested loops in azimuth, trajectory under azimuth, line segment on trajectory, energy group, and polar angle. The transport scan has the largest memory access and computational workload in the entire source iteration implementation. The boundary flux and starting flux accessed are both multiples of the number of 3D line segments, while the source term and scalar flux are multiples of the number of FSRs. In the actual benchmark test, the number of 3D line segments can reach 10 12 , and the number of flat source areas can reach 10 5 For serial computing programs, transport scanning accounts for more than 90% of the entire iteration time. In the analysis of transport scanning calculations, it is found that a large number of trajectories combined with angles, energy groups, etc. can provide excellent parallelism, which is essentially repeating a large number of identical operations on different data. Based on this, the SIMT thread parallel characteristics provided by DCU can be used to accelerate transport scanning calculations on a large scale in parallel.

[0074] In this embodiment, four layers of nested loops are required to calculate the angular flux, namely the number of 2D trajectories, the number of polar angles, the number of 3D trajectories under the polar angles, and the number of energy groups. Depending on the size of the parallel granularity, there are two parallel schemes to choose from, namely the Plane-level parallel algorithm (PLA) and the Track-level parallel algorithm (TLA).

[0075] The PLA parallel algorithm specifically refers to a DCU thread responsible for trajectory generation segmentation and angular flux calculation and contribution value of an axial plane. The corresponding DCU thread tasks are as follows: Figure 2 As shown in the figure, the three planes in the figure are managed by three DCU threads respectively. The trajectory segmentation of a plane is performed by nesting the trajectory stack according to the polar angle. When segmenting a trajectory stack, if z0(0) is defined as the ordinate of the intersection of the lowest position trajectory (numbered 0 in the plane) and the z axis, z i (0) represents the ordinate of the intersection of trajectory i and the z axis. Let the horizontal coordinate of a point on the trajectory be s, then the vertical coordinate of the point with horizontal coordinate s on the i-th trajectory is z i (s), which is calculated as:

[0076] z i (s)=z0(0)+iΔz+s*cotθ

[0077] The TLA parallel algorithm specifically refers to a DCU thread responsible for a 3D trajectory generation segment and angular flux calculation and contribution value. The corresponding DCU thread tasks are as follows: Figure 3 As shown in the figure, n 3D trajectories in an axial plane will be handled by n DCU threads. In the TLA parallel algorithm, the current thread tid is calculated through the grid-block-thread three-layer mapping relationship, and the thread number corresponds to the 3D trajectory number. By recording the index of the first trajectory number in the trajectory stack, the sequence number of the 3D trajectory in the trajectory stack is searched by binary search, and the 2D trajectory information and polar angle information projected by the 3D trajectory are obtained, and then the 3D trajectory and trajectory segment are generated.

[0078] In this embodiment, the time performance of each function of the transport scan is evaluated to quantify the calculation and communication time. For the serial program, assuming the benchmark example x, the number of flat source division loops N mg and the number of sectors N sec , the total number of energy groups N egrp , the total number of azimuth angles N a , radial plane spacing δ ρ , the total number of polar angles N p , axial plane spacing δ θ , the total number of flat source areas is N FSR , the total number of 2D trajectories N 2dtrk , the total number of 2D line segments N 2dseg, the total number of 3D trajectories N 3dtrk , the total number of 3D line segments N 3dseg The total time formula for transport scan serial calculation is:

[0079] T s =N 3dtrk t1+N 3dseg t2+2N 3dtrk N egrp t3+2N 3dseg N egrp t4+2N 3dseg N egrp t5

[0080] Among them, T s represents the total running time of the transport scan series, t1 represents the time to generate a single 3D trajectory, t2 represents the time to generate a single line segment, t3 represents the time to calculate the angular flux of a single energy group of a single line segment, t4 represents the contribution of the angular flux change of a single energy group of a single line segment to the flat source area, and t5 represents the boundary flux of a single energy group of a single trajectory.

[0081] According to the execution process of trajectory generation and segmentation, it can be concluded that the total number of 3D trajectories is related to the azimuth, polar angle and spacing, because the number of 3D line segments is related to the number of 3D trajectories and the number of flat source areas, and the number of flat source areas is related to the geometry of the benchmark problem and the parameters of the division ring and sector. Finally, the number of energy groups is related to the structure of the benchmark example, and the more complex the example, the larger the number. The above four positive correlations are represented by F1, F2, F3 and F4 respectively, then:

[0082] T s =F1(N a ,δ ρ ,N p ,δ θ )(t1+2F4(x)t5)+F2(F1(N a ,δ ρ ,N p ,δ θ ),F3(x,N rng ,N sec ))(t2+2F4(x)(t3+t4))

[0083] For each function, the time spent includes computing time and communication time. The computing time refers to the time to perform floating-point addition, subtraction, multiplication and other calculation operations on the chip, and the communication time refers to the time to access data from the storage space. Then:

[0084] t=t comp +t comm

[0085] Combining the above two equations, we can get the formula for the total serial time of transport scanning:

[0086]

[0087] in, Refers to the time spent calculating the attribute values ​​of a single 3D trajectory. Refers to the time spent calculating the attribute values ​​of a single line segment. Refers to calculating the change in angular flux on a single energy group on a single line segment. Refers to the calculation of the contribution of the angular flux change on a single line segment and a single energy group to the standard flux in the flat source area. It is 0 because no floating point addition, subtraction, multiplication or other calculation operations are required when transmitting boundary flux. Refers to the total communication time of information related to a single 3D trajectory. Refers to the total communication time for accessing the line segment information on the projected 2D trajectory in the storage space, the axial extrusion flat source area, the flat source area centroid, and the material. Refers to the total communication time for accessing the line segment starting point angular flux, source term, material cross-section data and index calculation table lookup in the storage space when calculating the angular flux of a single energy group of a single line segment. Refers to the total communication time for accessing the weight information and the standard flux in the storage space, Refers to the communication time to access the boundary flux.

[0088] For a parallel program, assume that the total number of computing units used for floating-point calculations is p, the total execution time of the transport scan is T, and the total number of threads is N. thd According to the DCU programming characteristics, the total execution time of a single node consists of three parts, including the time T for copying data from the host to the device. h2d , kernel function execution time on DCU T k The time T for copying data from the device to the host d2h . And T k Including calculation time T comp and memory access communication T comm Since DCU uses wavefront as the execution unit, some threads in the wavefront are idle due to the slowest thread, resulting in waiting time. Assume that the waiting time is T idle , when the cache hit rate of each thread task is the same, we can get T k Calculation formula:

[0089]

[0090] Among them, λ∈(0,p) indicates the factor that weakens the acceleration effect due to wavefront divergence, ζ=0 indicates that the memory access time is 0, and ζ=1 indicates that each thread accesses the memory completely serially. Since the wavefront thread memory access can be merged and the DCU hardware hides the memory access overhead by switching workgroup and wavefront, ξ∈(0,1). Theoretically, the task execution of a single thread is serial, the difference is that the task volume is different from serial. The serial and parallel execution time are related to their respective device indicators, that is, t comp and t ’comp They are related to the number of floating-point operations of the respective running devices, t comm and t ’comm The components of single-threaded parallel computing and memory access are similar to those of serial computing.

[0091] If the total number of threads is equal to the number of planes, for the PLA algorithm, the number of 3D trajectories processed by each thread is not exactly the same as the number of threads. The total transport scanning time calculation formula for the PLA algorithm is:

[0092]

[0093] If the total number of threads is equal to the number of 3D trajectories, for the TLA algorithm, the total number of 3D line segments on a single trajectory is n sdseg =1,n egrp Indicates the number of benchmark example energy groups. In the TLA algorithm, the number of 3D line segments processed by each thread is not exactly the same, so the total transport scan time calculation formula is:

[0094]

[0095] From the above PLA and TLA time performance formulas, it can be seen that the running time of the transport scan is coupled by multiple factors.

[0096] In this embodiment, for T in the performance formula h2d , T d2h, and proposes a data transmission optimization solution. Before executing the kernel function, the host-side data needs to be copied to the device side. In the CPU-DCU architecture, the host side and the device side are connected by a PCIe bus. When executing the transport scanning kernel function, it is necessary to copy 2D trajectories, 2D line segments, 2D trajectory chains, axial flat source area information, and 3D trajectory index information. Although the total amount of these data is not large, each type of information is composed of no less than two. For example, 2D trajectory information includes nx array, ny array, trajectory ID index array and deviceTrack2D array, and 2D trajectory chain information includes trajectory chain length array, trajectory chain number array, etc. Most of these arrays are not one-dimensional arrays. For one-dimensional arrays, only one call to the hipMemcpyHtoD copy function is needed to transfer information from the host side to the device side, but for high-dimensional arrays, the copy function needs to be called multiple times to transfer a small amount of data, resulting in the API startup overhead being higher than the data transmission overhead itself. Therefore, on the host side, multi-dimensional data is changed to one-dimensional data by adding an index array method to reduce the number of times the copy function is started. Figure 4 As shown, the original N a *5+N 2dtrk *4+3+2*N FSR +N eFSR The number of copy calls (in some cases up to 10,000) is reduced to 19, which improves the data transfer rate.

[0097] In this embodiment, a register compilation optimization scheme is proposed for ζ in the performance formula. When there is no register overflow, all local variables in the transport scan kernel function can be stored in registers. For the PLA parallel strategy, its transport scan kernel function contains azimuth and 2D trajectory loops, which requires the storage of a large number of local variables, and the register usage of a single thread is relatively high; while for the TLA parallel strategy, the register usage is relatively small. The DCU register resources are fixed, and the register usage of a single thread will affect the number of active wavefronts, reduce the effect of hiding the memory access overhead of wavefront and block switching, and greatly affect the parallel performance.

[0098] Use the _launch_bounds compiler directive provided by HIP to affect the registers of the kernel function. The compiler makes static decisions about register allocation based on two values: MAX_THREADS_PER_BLOCK and MIN_WARPS_PER_EU. The former ensures that the kernel is launched with a number of threads less than this value. Specifying it is less than 1024 effectively allows the compiler to use more resources than the default unconstrained. The latter instructs the compiler to minimize resource usage so that the required number of wavefronts can be active on the CU at the same time, mainly reducing the number of fixed resources (such as registers) used by the wavefront.

[0099] In this embodiment, for t' in the performance formula comm , and proposed a solution to reduce memory access overhead. In terms of program implementation, the data copied from the device side and the 3D trajectory and segment data requested on the device side are all stored in the video memory. Figure 5 The multi-level storage system of the domestic DCU is demonstrated, in which the memory access efficiency of constant storage and shared storage is much higher than that of global storage. The geometric boundary information used during trajectory generation is unchanged during the transport scan, and this part of information is small. The threads in a wavefront almost access the same location. This part of information can be placed in the constant memory of the DCU to reduce the thread's memory access latency. When calculating the angular flux contribution value, the weight array of the product group is also unchanged and small. The subscripts of this weight array are the azimuth and polar angle indexes. For the TLA algorithm, the threads in the same wavefront have almost the same azimuth and polar angle, that is, they access the same weight information, so this part of weight information can be stored in the DCU constant memory.

[0100] At the same time, the starting point of the 3D trajectory is placed in shared storage to reduce the algorithm memory access overhead. The starting point information of the 3D trajectory will be reused, but the memory access locations of the threads in the same wavefront are not adjacent in the video memory. When thread t needs to access the starting point of the i-th trajectory, the address value of this starting point is p. Because the 3D trajectory is stored continuously, and the compiler will perform memory alignment when compiling the 3D trajectory structure, the access location that thread t+1 needs to access the starting point of the i+1th trajectory is p+192 (assuming the machine is byte-addressed), that is, the threads in the same wavefront do not access the video memory continuously, resulting in only a very small amount of bandwidth being effective in the 512-byte memory transaction, and multiple requests for access to the video memory operation are required to meet the data required for the execution of the threads in the same wavefront, which prevents the hardware from fully merging the memory access operations. Therefore, in response to this phenomenon, shared memory is used to store the starting point information of the 3D trajectory. Since each iteration trajectory and line segment are newly generated, the trajectory starting point information can be considered as an intermediate value in the transport scanning process, and it is not necessary to write it back to the 3D trajectory information stored in the video memory. The trajectory starting point information required in each transport scan is replaced with the information stored in the shared memory. In the algorithm implementation, the shared memory is used for software data cache. Assume that the shared memory consists of 32 banks, each bank is 64 bytes wide. Figure 6 The optimized memory access merge process is shown. The threads of the same wavefront only need to access the shared memory twice to obtain the required data, which is more efficient than the non-merged memory access on the video memory.

[0101] In this embodiment, when executing the TLA algorithm, it is necessary to pre-allocate the maximum number of line segments that can be processed in the video memory for each thread on the device side. These spaces are partially redundant because the trajectory with the largest number of line segments often only occupies a small part. Taking the C5G7 benchmark example, when the azimuth angle is 16, the azimuth spacing is 0.05cm, the number of polar angles is 4, and the polar angle spacing is 0.1cm, 12047744 3D trajectories will be generated, and the maximum number of line segments in the trajectory is 146. When a line segment space of 146 is pre-allocated to each thread, 30% to 95% of the space is wasted, which does not meet the requirements for storage utilization. Therefore, a pre-allocation space scheme is proposed: the thread processes the trajectory according to the trajectory number, that is, the trajectory number processed by the thread is known before the kernel function is executed, and a direct mapping method can be used to calculate the maximum number of line segments in the trajectory offset by an integer multiple of the number of threads, that is, the maximum number of line segments in the trajectory that this thread can process, rather than the maximum number of line segments in the entire trajectory. The total number of line segments that each thread can process is allocated in advance in the video memory. Then, a one-dimensional pointer array is used to store the starting position of each thread's allocated space according to the thread number. This allocation method converts the two-dimensional address space into a one-dimensional space. Only one memory allocation function hipMalloc needs to be called, and space utilization is improved (PLA's idea of ​​reducing storage is similar to TLA). The changes in line segment storage space before and after optimization are as follows: Figure 7 shown.

[0102] In this embodiment, 3D C5G7 and Takeda full stack are used as examples to verify the correctness of the results, test the time performance and test the storage performance. The test shows that the parallel characteristic line calculation method for heterogeneous architecture constructed by the present invention can accelerate the transport scanning calculation and reduce the storage overhead while ensuring that the calculation results are consistent with the serial execution.

[0103] The parallel characteristic line method calculation and storage optimization method for domestic heterogeneous architecture provided in this embodiment also provides a considerable reference solution for the transplantation and optimization of parallel computing programs in other related physical fields. In addition, the performance model constructed in this patent quantifies the various factors that affect the performance of parallel programs, can expose bottlenecks, and accurately show where the code and system can be improved, thereby improving efficiency.

[0104] The above description is only a preferred specific implementation manner of the present invention, but the protection scope of the present invention is not limited thereto. Any technician familiar with the technical field can make equivalent replacements or changes according to the technical scheme and inventive concept of the present invention within the technical scope disclosed by the present invention, which should be covered by the protection scope of the present invention.

Claims

1. A parallel characteristic line method calculation and storage optimization method for domestic heterogeneous architecture, characterized in that: The method includes the following: S1. Construct a numerical model based on the reactor geometry, perform meshing to obtain the flat source area and determine the filling material, volume, and centroid data of each flat source area; S2, transport calculation pre-processing stage, generates a 2D loop trajectory chain on the radial plane, searches for axial geometric information to generate an axial extrusion flat source area, generates 3D trajectories and segmentation-related index information for device-side execution, and copies the calculation information required for executing the kernel function to the device side; S3, transport calculation stage, the kernel function is started on the host side, and the device side sequentially performs source term update of the flat source area, transport scanning calculation, standard flux update, and calculation of effective neutron multiplication factor and residual; S4. Copy the calculation results of the device side back to the host side to determine whether it has converged. If not, return to the previous step. Otherwise, the iteration ends.

2. According to claim 1, the parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture is characterized in that: The solution of the source term update of the flat source area described in S3 in each flat source area and each energy group is relatively independent, and the total number of parallel threads is the product of the number of flat source areas and the number of energy groups. The source term calculation of a single energy group on a single flat source area is assigned to a DCU thread, and its corresponding kernel function is the sum of the fission rate and the scattering rate of the flat source area, and the source term of the flat source area under a given energy group is calculated.

3. According to claim 2, the parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture is characterized in that: Each transport scan described in S3 calculates all fluxes of the trajectory under each energy group, and simultaneously derives a new flux contributing to each FSR; A single transport scan involves five layers of nested loop scans at azimuth, trajectory under azimuth, line segment on trajectory, energy group and polar angle; the transport scan calculation is based on the SIMT thread parallel characteristics provided by DCU to achieve large-scale parallel acceleration; In one cycle of the nested loop scanning, the corresponding kernel function allocates computing tasks to the DCU thread according to different algorithms, completes 3D trajectory generation and segmentation, and performs transport scanning along the trajectory, calculates the trajectory incident neutron angular flux, outgoing neutron angular flux and the difference between the two, and finally updates the boundary flux according to the boundary conditions.

4. According to claim 3, the parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture is characterized in that: When calculating the angular flux change of each 3D trajectory in the flat source area, the following two solutions are provided according to the different parallel granularity sizes: ①Plane-level parallel algorithm, denoted as PLA, specifically refers to a DCU thread responsible for trajectory generation segmentation, angular flux calculation and contribution value of a plane; Track-level parallel algorithm, denoted as TLA, specifically refers to a DCU thread responsible for a 3D track generation segment and angular flux calculation and contribution value; In different parallel algorithm implementations, the logic of generating a single 3D trajectory remains unchanged, and the only change is the number.

5. According to claim 4, the parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture is characterized in that: The time performance of the five functions of the transport scan, including trajectory generation, segmentation, calculation of forward and reverse angular flux and contribution value of line segments, forward and reverse boundary flux of the transmission trajectory, and copying data between the host and device ends, is evaluated to quantify the calculation and communication time, including: For serial programs, the total serial scan time formula is: in, Refers to the time spent calculating the attribute values ​​of a single 3D trajectory. Refers to the time spent calculating the attribute values ​​of a single line segment. Refers to calculating the change in angular flux on a single energy group on a single line segment. Refers to the calculation of the contribution of the angular flux change on a single line segment and a single energy group to the standard flux in the flat source area. It is 0 because no floating point addition, subtraction, multiplication or other calculation operations are required when transmitting boundary flux; Refers to the total communication time of information related to a single 3D trajectory. Refers to the total communication time for accessing the line segment information on the projected 2D trajectory in the storage space, the axial extrusion flat source area, the flat source area centroid, and the material. Refers to the total communication time for accessing the line segment starting point angular flux, source term, material cross-section data and index calculation table lookup in the storage space when calculating the angular flux of a single energy group of a single line segment. Refers to the total communication time for accessing the weight information and the standard flux in the storage space, Refers to the communication time to access the boundary flux; For parallel programs, if the total number of threads is equal to the number of Planes, for the PLA algorithm, the number of 3D trajectories processed by each thread is not exactly the same as the number of threads, so the total transport scanning time calculation formula is: If the total number of threads is equal to the number of 3D trajectories, for the TLA algorithm, the total number of 3D line segments on a single trajectory is n sdseg =1,n egrp Indicates the number of benchmark example energy groups; the number of 3D line segments processed by each thread in the TLA algorithm is not exactly the same, so the total transport scan time calculation formula is: Among them, T h2d represents the time it takes to copy data from the host to the device; λ∈(0,p) represents the factor that weakens the acceleration effect due to wavefront divergence; T comp represents the calculation time; p represents the total number of computing units used for floating-point calculations; T idle Indicates waiting time; T comm Indicates memory access communication; T h2d Indicates the time it takes for data to be copied from the host to the device; N thd Indicates the total number of threads; n egrp Indicates the number of benchmark example groups; Refers to the time it takes for a thread in a parallel program to calculate the attribute values ​​of a single 3D trajectory; Refers to the time it takes for a thread in a parallel program to calculate the change in angular flux on a single energy group on a single line segment; Refers to the total communication time of information related to a single 3D trajectory within a thread in a parallel program; It refers to the total communication time for a thread in a parallel program to access the segment starting point angular flux, source term, material cross-section data and index calculation table lookup in the storage space when calculating the angular flux of a single energy group of a single segment; ξ = 0 means that the memory access time is 0; ξ = 1 means that each thread accesses the memory completely serially.

6. According to claim 5, the parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture is characterized in that: Regarding T in the performance formula h2d , T d2h , a data transmission optimization plan is proposed, which specifically includes the following contents: when executing the transport scanning kernel function, the multidimensional data is changed into one-dimensional data on the host side by adding the index array method, reducing the number of times the hipMemcpyHtoD copy function is opened and reducing the data transmission time; the multidimensional data includes 2D trajectory, 2D line segment, 2D trajectory chain, axial flat source area information and 3D trajectory index information.

7. According to claim 6, the parallel characteristic line method calculation and storage optimization method of a domestic heterogeneous architecture is characterized in that: Aiming at ξ in the performance formula, a register compilation optimization scheme is proposed, which includes the following contents: using the _launch_bounds compiler instruction provided by HIP to affect the register usage of kernel functions, ensuring that the appropriate number of workgroups and wavefronts reside on the CU, and hiding the overhead by switching execution units by hardware.

8. A method for parallel characteristic line calculation and storage optimization of a domestic heterogeneous architecture according to claim 7, characterized in that: For the t in the performance formula ’comm , a scheme to reduce memory access overhead is proposed, which specifically includes the following contents: ,Put the geometric boundary information used during trajectory generation and the weight information of the ,product group when calculating the angular flux contribution value into the constant memory of ,DCU to reduce the memory access latency.

9. A method for parallel characteristic line calculation and storage optimization of a domestic heterogeneous architecture according to claim 1, characterized in that: Aiming at the TLA algorithm, a pre-allocation space scheme is proposed, which specifically includes the following contents: using direct mapping method to calculate the maximum number of line segments in the trajectory with an offset of an integer multiple of the number of threads; comprehensively considering the number of line segments that each thread can process, pre-allocate space in the video memory, and then use a one-dimensional pointer array to store the starting position of each thread allocation space according to the thread number.

Images