Neutron transport scanning operator optimization method based on built-in function and instruction level parallelism

Abstract

The invention provides a neutron transport scanning operator optimization method based on built-in function and instruction level parallelism, and the method comprises the steps: firstly reconstructing a data layout, declaring an array according to an energy group sequence, carrying out the memory alignment, and enabling a calculation process to adapt to an SIMD vectorization architecture; then implementing explicit vectorization calculation by adopting a built-in function, performing polar angle cycle expansion and instruction recombination, and fully mining instruction-level parallel potential; meanwhile, key index calculation is optimized into vectoring polynomial approximation. Based on the optimized data flow and the parallel computing mode, efficient updating of the angular flux and the standard flux is achieved, and convergence judgment and calculation are completed. According to the method, the calculation throughput rate and the hardware resource utilization rate of the scanning operator are remarkably improved, and the method can be coordinated with coarse-grained parallel strategies such as MPI / OpenMP and is suitable for a multi-architecture high-performance calculation environment.

Description

Technical Field

[0001] This invention relates to the fields of nuclear reactor and high-performance computing technology, and in particular to an optimization method for neutron transport scanning operators based on built-in functions and instruction-level parallelism. Background Technology

[0002] Reactor physics analysis is fundamental to nuclear reactor system design and safety assessment. Its core lies in obtaining core reactivity and fine power distribution across the entire reactor by solving neutron transport equations. The method of characteristics, due to its high computational accuracy, strong geometric adaptability, effective handling of anisotropic scattering, and multiple parallelizable dimensions, has become one of the mainstream methods for neutron transport calculation.

[0003] As the complexity of advanced nuclear reactors in terms of geometry, material distribution, and energy spectrum characteristics increases, engineering computational models are shifting from traditional few-group approximations to refined energy group discretization. This trend has led to a dramatic increase in the computational scale and memory bandwidth requirements of the Method of Characteristics (MOC) for neutron transport characteristics. Against this backdrop, the ray scanning operator, as the most time-consuming core component in MOC computation, directly determines the upper limit of the efficiency of the entire transport simulation. However, the implementation of existing MOC scanning operators faces a series of performance bottlenecks, making it difficult to meet the requirements of refined energy group computation: First, the calculation process involves multiple nested loops such as azimuth angle, ray, ray segment, and energy group, which have data dependencies and load imbalance problems, severely restricting the parallel efficiency of the loop level.

[0004] Secondly, the control flow and complex computational processes involved in flux integration and source term calculation hinder the compiler's automatic analysis and vectorization conversion, causing the computational process to fail to fully utilize the processor's SIMD vector units.

[0005] Third, the memory layout of core data such as cross-sections and throughput lacks targeted optimization, and no contiguous storage and address alignment declarations are made. This not only prevents efficient vector loading and storage operations but also causes a decrease in cache hit rate. At the same time, arrays that are not padded according to the vector length force the automatic vectorization process to degenerate into low-width or even scalar mode, severely limiting performance improvement.

[0006] Fourth, special functions such as exponent calculation are frequently called during the scanning process, and their inherent high computational latency and bandwidth requirements further exacerbate the overall performance bottleneck.

[0007] To address these challenges, the industry typically employs methods like MPI and OpenMP to perform coarse-grained parallel decomposition of regions or rays at the task level. While this achieves some scalability, it doesn't address the root causes of performance issues in the underlying computational microkernel. Its core computation still operates in scalar or weakly vectorized modes, failing to systematically utilize the wide vector registers, fused multiply-accumulate units, and instruction-level parallelism capabilities offered by modern processor architectures. Furthermore, compiler-automatic vectorization often fails when faced with data dependencies, complex branches, and non-contiguous memory accesses; even when successful, it often fails to achieve the desired speedup due to issues such as unaligned memory accesses, register allocation pressure, and execution port contention.

[0008] In summary, while existing technologies are mature in task-level parallelism, they generally lack an integrated performance optimization scheme that addresses the underlying computation of MOC scan operators, encompassing data layout, vectorization, instruction-level parallelism, cache optimization, portability, and numerical accuracy consistency. This makes it difficult to simultaneously achieve high throughput and stable, reproducible computational accuracy and efficiency when dealing with modern computing scenarios involving fine energy groups and massive parallelism. Summary of the Invention

[0009] To address the shortcomings of existing technologies, this invention provides an optimization method for neutron transport scanning operators based on built-in functions and instruction-level parallelism. This method solves the problems of low utilization of processor SIMD units and insufficient exploitation of instruction-level parallelism capabilities caused by discontinuous data layout, complex computational logic, and incomplete vectorization of computation in existing implementations. It significantly improves computational throughput and hardware resource utilization efficiency while maintaining numerical stability and physical correctness, and ensures its portability across different CPU architectures and vector instruction sets.

[0010] The technical solution of this invention is as follows: a neutron transport scanning operator optimization method based on built-in functions and instruction-level parallelism, which optimizes the low-level calculation of the ray scanning operator in the characteristic line method transport calculation, specifically including the following steps: S1) Start the program, the user inputs core modeling and sets boundary conditions, and generates a set of feature lines and line segments; S2) Declare the key arrays in transport calculations in energy group order and fill and align them with memory using vector length as the granularity to establish an energy group block aligned access mode. S3) Taking the innermost energy group iteration as the target, SIMD parallel computation is explicitly implemented using built-in functions, and the angular flux decay and source term update are uniformly rewritten in FMA form; S4) After cyclically expanding multiple polar angles, group similar instructions together to explicitly implement instruction-level parallelism (ILP) optimization; S5) Perform parallel calculations of exponential functions optimized by vectorization and polynomial approximation for multiple polar angles; S6) Write the outgoing angle flux of the current line segment into the incoming angle flux of the next line segment. After processing all line segments, perform global reduction of the source and flux, and update the fission source and effective multiplication coefficient. And determine the updated fission source and Check if convergence has occurred; if convergence has occurred, the program ends; otherwise, return to step S3) to proceed to the next source iteration.

[0011] Preferably, in step S2), the key array includes the total cross section, source term, angular flux, and standard flux.

[0012] Preferably, in step S2), the key arrays in the transport calculation are declared in energy group order and array filling and memory alignment are performed with vector length as the granularity to establish an energy group block access mode, as follows: Total cross-section array Source item array Angular flux array and standard flux array The storage order is uniformly set to follow the energy group. Continuous index; Declare and allocate a vector register with a length equal to the total number of energy groups, rounded up to the target platform length. An integer multiple array, while ensuring that the array's starting address is aligned with the target platform's vector register length, is used to establish a group block. Access mode.

[0013] Preferably, in step S3), the innermost energy group iteration is the target, and SIMD parallel computation is explicitly implemented using built-in functions, as follows: Traverse the energy group with the length of the target platform's vector register as the step size, processing one energy group block at a time. The change in angular flux is decomposed into two sequences of fused negative multiplication plus FNMADD instructions; The energy group block is completed in one go using vectorized built-in functions. Loading, computing, and storing operations.

[0014] Preferably, in step S4), after the polar angle loop is expanded, the instructions are arranged in the same type of instruction set. Vectorized energy group scanning is performed on each expanded polar angle using a rolling register. Software pipelining and parallelization are performed after the polar angle loop is expanded. Simultaneously, constants related to the polar angle but unrelated to the energy group in the loop are outsourced, and the constants that remain unchanged in the inner energy group loop are... Calculation and constants Move it outside the loop and use the SIMD broadcast instruction to broadcast each polar angle. and Load into the vector register for cyclic multiplexing.

[0015] Preferably, in step S5), a vectorized polynomial approximation method based on SIMD architecture is adopted to transform scalar exponent calculation into an efficient parallel computing mode. The exponential function calculation is decomposed into three parallel processing stages, namely: First, the input parameters are mapped to a finite interval through floating-point operations, and then multiple data are simultaneously transformed into intervals using vector instructions. Next, a parallel polynomial computation method based on fused multiply-accumulate instructions is used to perform batch approximation calculations on the mapped data; Finally, the polynomial calculation results are batch adjusted using processor vector scaling instructions to obtain the final exponential function output.

[0016] The beneficial effects of this invention are as follows: 1. Based on MPI / OpenMP hybrid parallelism, this invention further improves the computing performance of hardware through systematic vectorization reconstruction from data layout to instruction flow; 2. This invention achieves stable and efficient SIMD vectorization in the energy group dimension in the scanning operator, and integrates instruction-level parallelism and targeted exponential function optimization, thereby fully stimulating the potential of wide vector units and multi-issue pipelines of modern processors at the hardware level. 3. This invention not only significantly improves computational throughput and parallel efficiency, but also has excellent cross-platform portability while ensuring numerical accuracy. It can be directly applied to modern multi-core, wide-vector, multi-issue, and superscalar high-performance computing environments with various architectures. 4. By reconstructing the data layout and setting the energy group dimension as the main loop dimension, this invention ensures that each vectorized loading, computation, and storage operation operates on continuously aligned memory addresses, thereby transforming the original unaligned memory access mode into a fully sequential and efficient pipeline, eliminating the performance loss caused by unaligned memory or step access. 5. This invention employs a vectorized polynomial approximation method based on SIMD architecture, transforming traditional scalar exponent calculation into an efficient parallel computing mode. Attached Figure Description

[0017] Figure 1 This is a flowchart illustrating the method of the present invention. Detailed Implementation

[0018] The specific embodiments of the present invention will be further described below with reference to the accompanying drawings: Example 1 like Figure 1 As shown, this embodiment provides a neutron transport scanning operator optimization method based on built-in functions and instruction-level parallelism. It optimizes the low-level computation of the ray scanning operator in the characteristic line method transport calculation, specifically including the following steps: S1) Start the program, the user inputs core modeling and sets boundary conditions, and generates a set of feature lines and line segments; S2) Declare the key arrays in transport calculations in energy group order and fill and align them with memory using vector length as the granularity to establish an energy group block aligned access mode. In this embodiment, the key array total cross-section array Source item array Angular flux array and standard flux array .

[0019] By using the total cross-section array Source item array Angular flux array and standard flux array The storage order is uniformly set to follow the energy group. Continuous index; Declare and assign a length to the energy group The total number is rounded up to the length of the target platform's vector register. An integer multiple array, while ensuring that the array's starting address is aligned with the target platform's vector register length, is used to establish a group block. Alignment access mode.

[0020] This embodiment ensures that each vectorized load, compute, and store operation operates on contiguously aligned memory addresses through the above operations, thereby transforming the original unaligned memory access mode into a fully sequential and efficient pipeline, eliminating the performance loss caused by unaligned memory or step access.

[0021] S3) Taking the innermost energy group iteration as the target, SIMD parallel computation is explicitly implemented using built-in functions, and the angular flux decay and source term update are uniformly rewritten in FMA form; In this embodiment, SIMD is a data-level parallelism technique that enables the simultaneous execution of the same operation on multiple data elements under the control of a single instruction. The built-in functions used in this embodiment can be directly mapped to underlying SIMD instructions, thereby achieving precise control over instruction-level behavior while maintaining the convenience of high-level language programming. In this embodiment, the built-in functions are a special type of function that the compiler can recognize. During the compilation process, they are directly converted into machine instructions or instruction sequences corresponding to the target architecture. Unlike regular function calls, built-in functions do not follow standard function calling conventions, and the compiler can perform deep optimizations on them, including inlining, instruction selection, and scheduling.

[0022] Specifically as follows: Traverse the energy group with the length of the target platform's vector register as the step size, processing one energy group block at a time. And the change in angular flux in equation (1) Decompose it into two sequences of fused negative multiplication plus FNMADD instructions, equations (5) and (6), and implement equation (3) using vector subtraction. Update and write back immediately; and complete the flux accumulation calculation of equation (4) with one FMADD and write back once; where: (1) ; (2) ; (3) ; (4) ; (5) ; (6) In the formula, The change in angular flux; Indicates the polar angle; Represents an energy group; The incident angular flux at the starting point of the ray segment; The length of the characteristic line segment; For the outward angle flux; Describe any energy group In the Standard flux on FSR in a flat source region ; It is a constant; Indicates the neutron optical path (optical length); This is a temporary parameter for the term in parentheses on the left side of equation (1); Indicates a neutron source; Represents the overall macroscopic cross-section; Indicates polar angle The cosecant function value; FMADD and FMADD respectively represent the fusion of negative multiplication, addition, and subtraction instructions and multiplication, addition, and subtraction instructions. FMADD ;⊙ represents element-wise multiplication.

[0023] Finally, the vectorized built-in function is used to complete the energy group block in one go. Loading, computing, and storing operations.

[0024] S4) Perform cyclic expansion on multiple polar angles and explicitly implement instruction-level parallelism (ILP) optimization by grouping similar instructions; the details are as follows: S41) Use the vector loading instruction to sequentially load the relevant data of the current energy group block into the vector register; S42) Read the FSR standard flux corresponding to the current energy group block through a single vector loading operation. To the vector register; S43), for different polar angles The memory access, computation, and storage instructions are reorganized, and similar instructions are grouped together to explicitly implement ILP; S44) Using a vector merging and write-back method, the scalar flux of multiple polar angles is... The update results are written back to memory all at once; Among them, the standard flux The update expression is: ; (4) In the formula, Describe any energy group In the scalar flux on each FSR ; It is a constant.

[0025] In this embodiment, before entering the innermost energy group cycle, the polar angle cycle is expanded and the instructions are reorganized and arranged to carry out static ILP optimization. The goal is to allow more independent instructions to enter the out-of-order executor and execution port simultaneously through source code-level transformation without changing data dependencies and numerical semantics, thereby masking long latency and improving throughput per unit time. After polar-angle cycle expansion, software pipeline parallelization is implemented. Simultaneously, by hoisting out the cycle invariants, the invariants in the inner energy group cycle are... Calculation and constants Move it outside the loop and use the SIMD broadcast instruction to broadcast each polar angle. and Load the vector register for cyclic reuse of the power group, reducing redundant loading and scalar computation.

[0026] S5) Initiate parallel calculations of exponential functions optimized by vectorization and polynomial approximation for multiple polar angles to mask their latency, thereby stabilizing wide vector and FMA throughput and reducing memory access overhead. In this embodiment, a vectorized polynomial approximation method based on SIMD architecture is adopted to transform scalar exponential calculation into an efficient parallel computing mode. In this embodiment, the exponential function calculation is decomposed into three parallel processing stages, namely: First, the input parameters are mapped to a finite interval through floating-point operations, and then multiple data are simultaneously transformed into intervals using vector instructions. Next, a parallel polynomial computation method based on fused multiply-accumulate instructions is used to perform batch approximation calculations on the mapped data; Finally, the polynomial calculation results are batch adjusted using processor vector scaling instructions to obtain the final exponential function output.

[0027] S6) Write the outgoing angle flux of the current line segment into the incoming angle flux of the next line segment. After processing all line segments, perform global reduction of the source and flux, and update the fission source and effective multiplication coefficient. And determine the updated fission source and Check if convergence has occurred; if convergence has occurred, the program ends; otherwise, return to step S3) to proceed to the next source iteration.

[0028] In this embodiment, the effective multiplication coefficient is updated using the neutron balance method (i.e., the ratio of total fission rate to the sum of total absorption rate and leakage rate). Utilizing the new effective multiplication coefficient and the normalized standard flux reconstructs the fission source; Subsequently, the root-mean-square relative residual of the reconstructed fission source and the new effective multiplication coefficient were calculated. The absolute change, if the residual is less than the set threshold and the new effective multiplication coefficient If the change is less than 1 pcm, then convergence is determined and the program ends; otherwise, return to step S3) and use the updated source term for the next source iteration.

[0029] The embodiments and descriptions above are merely illustrative of the principles and preferred embodiments of the present invention. Various changes and modifications may be made to the present invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed.

Claims

1. A method for neutron transport scan operator optimization based on built-in functions and instruction level parallelism, characterized in that, The method comprises the following steps: S1), starting the program, the user inputs the core modeling and sets the boundary conditions, generates the characteristic line and the line segment set; S2), the key array in the transport calculation is declared in the energy group order and the array filling and memory alignment are performed with the vector length as the granularity, the energy group block access mode is established; S3), the innermost energy group iteration is taken as the target, the SIMD parallel calculation is explicitly implemented by using the built-in function, the angular flux attenuation and the source term update are uniformly rewritten into the FMA form; S4), after the loop unrolling of multiple polar angles, the same type of instructions are grouped, and the instruction level parallel ILP optimization is explicitly implemented; S5), the exponential function calculation optimized by vectorization and polynomial approximation is initiated in parallel for multiple polar angles; S6), write the exit angle flux of the current line segment into the next line segment incident angle flux, after completing the processing of all line segments, perform global scaling of the source and flux, and update the fission source and effective multiplication factor ; and determining whether the updated fission source and converge; if so, the procedure ends, otherwise, returning to step S3) for the next source iteration.

2. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 1, wherein: In step S2), the key array in the transport calculation is declared in the energy group order and the array filling and memory alignment are performed with the vector length as the granularity, the energy group block access mode is established, and the specific process is as follows: The storage order of the total cross section array , the source term array , the angular flux array , and the score flux array is set to be continuous along the energy group index. A vector register length of the target platform is declared and assigned to the length of the array An integer multiple of the array length is calculated while ensuring that the first address of the array is aligned to the vector register length of the target platform, and a vector group block is established Access pattern.

3. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 2, wherein: In step S3), the innermost energy group iteration is taken as the target, the SIMD parallel calculation is explicitly implemented by using the built-in function, and the specific process is as follows: Traverse the capability groups in steps of the vector register length of the target platform, processing one capability group block at a time and decomposing the angular flux change into a sequence of two fused negative multiply add (FNMADD) instructions The energy group block is completed in one go using vectorized built-in functions. Loading, computing, and storing operations.

4. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 3, wherein: In step S3), the angular flux change is decomposed into two sequences of fusion negative multiplication addition FNMADD instructions, and the specific process is as follows: The angle flux variation amount of formula (1) is decomposed into two sequences of fusion negative multiply-add FNMADD instructions of formula (5) and formula (6), and formula (3) is implemented by vector subtraction The update of formula (4) is immediately written back; And one FMADD is used to complete the flux accumulation calculation of formula (4) and write back once; Wherein: ；（1） ； （2） ； （3） ； （4） ； （5） ； （6） In the formula, is the angular flux change amount; denotes the polar angle; denotes the energy group; is the incident angular flux at the start point of the ray segment; is the characteristic line segment length; is the exit angular flux; denotes any energy group On the FSR of the first flat source region, the tag flux ; is the characteristic line segment weight; denotes the neutron optical path; is the temporary parameter of the left bracket item in formula (1); denotes the neutron source; denotes the macroscopic total cross section; denotes the cotangent function value of the polar angle ; , FMADD respectively denote fusion negative multiplication and addition, multiplication and subtraction instructions; ; FMADD ; is an element-by-element multiplication.

5. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 4, wherein: In step S4), multiple polar angles are looped out, and the same type of instructions are grouped to explicitly implement the instruction level parallel ILP optimization; The specific process is as follows: S41), the vector loading instruction is used to load the related data of the current energy group block into the vector register in sequence; S42), reading the FSR tag flux corresponding to the current energy bin block by a single vector load operation to a vector register; S43), reorganize the memory, calculation and storage instructions with different polar angles , and arrange the same kind of instruction set to explicitly realize ILP; S44) Using a vector merging and write-back method, the scalar flux of multiple polar angles is... The update results are written back to memory all at once.

6. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 5, wherein: In step S4), software pipelining is performed after the polar angle loop is unrolled, while the invariant computations and constants are moved outside the loop and loaded into vector registers for each polar angle using SIMD broadcast instructions, for reuse in the energy group loop. and are loaded into vector registers for each polar angle using SIMD broadcast instructions, for reuse in the energy group loop.

7. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 6, wherein: In step S5), the vectorization polynomial approximation method based on the SIMD architecture is used to convert the scalar exponential calculation into an efficient parallel calculation mode; Specifically: The exponential function calculation is divided into three parallel processing stages, which are: Firstly, the input parameters are mapped to a limited interval through floating point operation, and the interval transformation is performed on multiple data at the same time by using the vector instruction; Then, the parallel polynomial calculation method based on the fusion multiplication addition instruction is used to perform batch approximation calculation on the mapped data; Finally, the polynomial calculation result is adjusted in batches by the processor vector scaling instruction to obtain the final exponential function output result.

8. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 1, wherein: In step S6, the effective multiplication factor is updated using the neutron balance method The fission source is reconstructed using the new effective multiplication factor and the normalized tally flux.

9. The built-in function and instruction level parallel based neutron transport scan operator optimization method of claim 8, wherein: In step S6), the root mean square relative residual of the reconstructed fission source and the new effective multiplication coefficient are calculated. The absolute change, if the residual is less than the set threshold and the new effective multiplication coefficient If the change is less than 1 pcm, then convergence is determined and the program ends; otherwise, return to step S3) and use the updated source term for the next source iteration.

Images